引言
下图是一张Landsat图像的示例（右图）。我们可以明显地看到四周的黑边和倾斜的角度。这是由于卫星传感器成像导致的。一般情况下，我们是不需要去除黑边和选择的，因为这样做之后投影信息和位置信息就不正确了。但对于做深度学习图像处理任务的同学来说，可能不需要考虑此类约束。因此，出于节省存储和内存的考虑，借助ChatGPT直接帮助我们生成了黑边去云及旋转校正的Python代码，最终效果如左图。
代码
Prompt: 现有一RGB彩色图像，图像内仅有一个旋转的矩形物体，其周围的像素值全是0，要求裁剪出其中的矩形物体，并将其旋转校正，注意矩形物体的边界也可能有值为0的像素，输出结果保存为彩色图像。使用Python实现。
- 安装依赖
pip install opencv-python numpy
- 算法流程
- 读取图像
- 转换为灰度图像，然后通过阈值化操作将所有非零像素转换为白色（255），所有零像素保持为黑色（0）。
- 使用OpenCV的findContours方法找到最大的轮廓，这应该是你的矩形物体。
- 使用minAreaRect获取轮廓的最小外接矩形（包括旋转角度）。
- 通过旋转原始图像将矩形校正。
- 裁剪并保存结果图像。
- 黑边去除及旋转校正的代码如下：
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
