您现在的位置是:首页 > 正文

【ChatGPT+Python】Landsat卫星图像黑边去云及旋转校正

2024-01-30 20:49:47阅读 0

引言

下图是一张Landsat图像的示例(右图)。我们可以明显地看到四周的黑边和倾斜的角度。这是由于卫星传感器成像导致的。一般情况下,我们是不需要去除黑边和选择的,因为这样做之后投影信息和位置信息就不正确了。但对于做深度学习图像处理任务的同学来说,可能不需要考虑此类约束。因此,出于节省存储和内存的考虑,借助ChatGPT直接帮助我们生成了黑边去云及旋转校正的Python代码,最终效果如左图。

在这里插入图片描述

代码

Prompt: 现有一RGB彩色图像,图像内仅有一个旋转的矩形物体,其周围的像素值全是0,要求裁剪出其中的矩形物体,并将其旋转校正,注意矩形物体的边界也可能有值为0的像素,输出结果保存为彩色图像。使用Python实现。

  • 安装依赖
pip install opencv-python numpy
  • 算法流程
  1. 读取图像
  2. 转换为灰度图像,然后通过阈值化操作将所有非零像素转换为白色(255),所有零像素保持为黑色(0)。
  3. 使用OpenCV的findContours方法找到最大的轮廓,这应该是你的矩形物体。
  4. 使用minAreaRect获取轮廓的最小外接矩形(包括旋转角度)。
  5. 通过旋转原始图像将矩形校正。
  6. 裁剪并保存结果图像。
  • 黑边去除及旋转校正的代码如下:
import cv2
import numpy as np

# 1. 读取图像
img = cv2.imread('input.jpg')

# 2. 转换为灰度图像
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 3. 找到最大的轮廓
_, thresh = cv2.threshold(gray, 1, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
contours, _ = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
contour = max(contours, key=cv2.contourArea)

# 4. 获取最小外接矩形
rect = cv2.minAreaRect(contour)
box = cv2.boxPoints(rect)
box = np.int0(box)

# 5. 旋转原始图像以校正矩形
width = int(rect[1][0])
height = int(rect[1][1])
src_pts = box.astype("float32")
dst_pts = np.array([[0, height-1],
                    [0, 0],
                    [width-1, 0],
                    [width-1, height-1]], dtype="float32")
M = cv2.getPerspectiveTransform(src_pts, dst_pts)
warped = cv2.warpPerspective(img, M, (width, height))

# 6. 显示中间过程(用于可解释性研究)
cv2.drawContours(img, [box], 0, (0, 0, 255), 2)
cv2.imshow('edge', img)
cv2.waitKey(0)
cv2.destroyAllWindows()

# 7. 裁剪并保存结果图像
output = warped[0:height, 0:width]
cv2.imwrite('output.jpg', output)

以上代码首先读取图像,并通过二值化将图像分割成背景和物体。然后找到物体的外接矩形,使用仿射变换将物体旋转到水平位置,最后对旋转后的图像进行裁剪。

但请注意,对于非常复杂或噪声很大的图像,这种方法可能不适用。你可能需要对图像进行更多的预处理,比如滤波、形态学操作等等。

实验

表1:边缘(轮廓)检测方法采用不同函数的对比结果

在这里插入图片描述

拓展

上述代码不仅对于Landsat卫星图像有效,对于普通的图像(周围是黑色填充,内接一个矩形物体)也同样有效。下面给出一个演示结果。

在这里插入图片描述

此外,上述代码中寻找最大轮廓除使用 cv2.findContours() 函数外,也可以使用如下方式替代,只不过代码的执行时间会变长,但是鲁棒性会更好,即抗噪声等干扰的能力强。

import cv2
import numpy as np

# 1. 读取图像
img = cv2.imread('input.jpg')

# 2. 转换为灰度图像
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# 3. 找到最大的轮廓
# 使用Canny边缘检测算法检测图像中的边缘
edges = cv2.Canny(gray, 50, 150, apertureSize=3)

# 使用霍夫变换检测图像中的直线
lines = cv2.HoughLines(edges, 1, np.pi/180, 200)

# 获取检测到的直线中最长的一条
longest_line = None
max_length = 0

for line in lines:
    for rho, theta in line:
        a = np.cos(theta)
        b = np.sin(theta)
        x0 = a*rho
        y0 = b*rho
        x1 = int(x0 + 1000*(-b))
        y1 = int(y0 + 1000*(a))
        x2 = int(x0 - 1000*(-b))
        y2 = int(y0 - 1000*(a))
        length = np.sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1))
        if length > max_length:
            max_length = length
            longest_line = line

# 获取直线的端点坐标
for rho, theta in longest_line:
    a = np.cos(theta)
    b = np.sin(theta)
    x0 = a*rho
    y0 = b*rho
    x1 = int(x0 + 1000*(-b))
    y1 = int(y0 + 1000*(a))
    x2 = int(x0 - 1000*(-b))
    y2 = int(y0 - 1000*(a))

# 计算直线的角度
angle = np.arctan2(y2-y1, x2-x1) * 180 / np.pi

# 4. 获取最小外接矩形
h, w = img.shape[:2]
rect = cv2.minAreaRect(np.array([(x, y) for x in range(w) for y in range(h) if edges[y, x] > 0]))
box = cv2.boxPoints(rect)
box = np.int0(box)

# 5. 旋转原始图像以校正矩形
width = int(rect[1][0])
height = int(rect[1][1])
src_pts = box.astype("float32")
dst_pts = np.array([[0, height-1],
                    [0, 0],
                    [width-1, 0],
                    [width-1, height-1]], dtype="float32")
M = cv2.getPerspectiveTransform(src_pts, dst_pts)
warped = cv2.warpPerspective(img, M, (width, height))

# 6. 显示中间过程(用于可解释性研究)
cv2.drawContours(img, [box], 0, (0, 0, 255), 2)
cv2.imshow('edge', img)
cv2.waitKey(0)
cv2.destroyAllWindows()

# 7. 裁剪并保存结果图像
output = warped[0:height, 0:width]
cv2.imwrite('output.jpg', output)

注意

Landsat的遥感影像四个角有黑色区域,这是正常的。表明那个黑色区的地方没有数据,在实际应用中也不需要去掉四个角的黑色区域。一般我们使用shp矢量量面来对遥感影像进行剪裁,提取出我们所需要研究或者显示的区域就可以啦。但是如果你真的想去掉黑色区域的话,你可以使用重分类,把黑色的区域变成白色,这样和背景就一致了,在发布服务的时候设为白色透明就可以了。还一个问题,Landsat8遥感影像是具有准确的投影坐标系的。倾斜是因为在卫星扫描的时候就是这个样子,拿到的影像是由准确的地理坐标的,不能够把它给旋转正了。旋转正了的话,它的投影信息和位置信息就不正确了,如果感觉不好看,可以直接剪裁出来正方形,或者把相邻的影像进行拼接,然后再剪裁。

参考

https://www.zhihu.com/question/497775240/answer/2216988184

网站文章

  • Struts的三种传参方式

    Struts的三种传参方式

    1、通过属性传参数(Attr) UserAction.java package com.bebig.struts2.user.action; import com.opensymphony.xwork2.ActionSupport; public class UserAction extends ActionSupport { private String name;

    2024-01-30 20:49:39
  • 【OpenCV C++&Python】(一)图像读取、显示和保存

    【OpenCV C++&Python】(一)图像读取、显示和保存

    文章目录OpenCV简介Mat图像存储方式显式创建Mat对象输出格式图像读取、显示和保存C++Python OpenCV简介 OpenCV(开源计算机视觉库)是一个开源库,是基于C/C++开发的: O...

    2024-01-30 20:49:33
  • 模型评估指标

    模型评估指标

    模型评估指标一、回归(Regression)算法指标1. Mean Absolute Error 平均绝对误差2. Mean Squared Error 均方误差3. Root Mean Square...

    2024-01-30 20:48:58
  • 【闲聊杂谈】HTTPS原理详解

    【闲聊杂谈】HTTPS原理详解

    HTTP虽然使用极为广泛, 但是却存在不小的安全缺陷, 主要是其数据的明文传送和消息完整性检测的缺乏, 而这两点恰好是网络支付, 网络交易等新兴应用中安全方面最需要关注的。关于 HTTP的明文数据传输...

    2024-01-30 20:48:49
  • org.hibernate.MappingException: Unknown entity:

    最近学习JEECG框架,使用代码自动生成功能并导入的过程后出现以下问题:页面能正常访问,但是首先前台页面出现NULL,相继后台打印出【org.jeecgframework.core.common.exception.MyExceptionHandler]java.lang.NullPointerException】点击确定后进行增删改查操作,编辑内容新增,提交后前台显示Unknown en

    2024-01-30 20:48:41
  • 管理oracle控制文件

    每一个oracle数据库都有一个控制文件。控制文件是一个小型的二进制文件,可以记录数据库的物理结构,包含以下的内容:数据库名称、相关数据文件和联机重做日志文件的名称和位置、数据库创建的时标、当前日志的序号、检验点信息。 无论何时打开数据库,控制文件必须能够由oracle数据库服务器写入内容。没有控制文件,数据库就不能装载,且很难恢复。oracle数据库控制文件在数据库创建的同时创

    2024-01-30 20:48:32
  • javascript 【2018.11.29】

    <html><head> <title></title></head><body> <script type="text/javascript"> for(i=0;i<10;i++) { if(i==3) break;

    2024-01-30 20:48:04
  • TCP/IP 三次握手

    TCP/IP 三次握手

    TCP/IP 三次握手

    2024-01-30 20:47:58
  • Java中高位转低位溢出的计算过程

    System.out.println((byte) 129);System.out.println((byte) -129);System.out.println("~b2: " + ~10);结果是:-127127~b2: -11计算机中是以补码进行计算正数的反码补码都是原码,如:10原码: 1010反码: 1010补码:1010负数 -10原码 10000000...

    2024-01-30 20:47:51
  • PAT - 乙级 1040 有几个PAT

    1040. 有几个PAT(25)时间限制120 ms内存限制65536 kB代码长度限制8000 B判题程序Standard作者CAO, Peng字符串APPAPT中包含了两个单词“PAT”,其中第一...

    2024-01-30 20:47:22