提问

以一张png做演示，不依赖gdal和rasterio第三方库，用python一步一步演示遥感影像几何校正的过程，写成函数，帮我学习几何校正，帮助理解几何校正的本质

以下是它的思考：

嗯，用户想让我用Python一步步演示遥感影像的几何校正过程，而且不依赖GDAL和rasterio这两个库。他们希望写成函数，帮助他们理解几何校正的本质。好的，首先我得理清楚几何校正的基本步骤，然后考虑如何用纯Python实现这些步骤，而不用那些专业的库。

首先，几何校正的主要步骤通常包括确定控制点、构建变换模型、然后对图像进行重采样。用户可能想了解的是多项式变换，比如一阶或二阶多项式，然后用最小二乘法来求解系数。接着是反向映射，用双线性插值等方法填充校正后的图像。不过，使用纯Python处理图像的话，可能需要用到PIL库来读取和保存图像，以及处理像素数据。

那用户的需求是不用GDAL和rasterio，所以需要找到替代的库。PIL或者Pillow应该可以处理基本的图像操作。但处理地理空间数据的话，可能需要自己解析坐标，不过用户可能只是需要一个演示，所以可能控制点是由用户提供的，而不是从地理文件中读取的。比如，用户自己输入控制点的坐标对。

接下来，构建多项式变换模型。比如，一阶多项式有六个系数，需要至少三个控制点；二阶的话需要更多。然后用最小二乘法来解这个方程组。这里可能需要用到numpy的线性代数求解功能，比如numpy.linalg.lstsq。不过用户可能希望不用太多第三方库，但numpy应该是可以接受的，毕竟纯Python处理矩阵运算会非常麻烦。

然后是图像的重采样部分。这一步需要遍历校正后的图像的每个像素，找到在原图中的对应位置，然后通过插值得到像素值。反向映射是比较常用的方法，比如对于每个校正后的像素点，用变换矩阵的逆来找到原图中的位置，然后用双线性插值计算像素值。这里需要处理的是，原图中的坐标可能超出图像范围，这时候需要处理边缘情况，比如填充0或者进行镜像等，但简单的处理可能是填充黑色。

接下来，我需要考虑如何将这些步骤拆分成函数。比如，一个函数用来计算变换系数，另一个函数用来应用变换并进行重采样。控制点的输入可能需要用户以列表形式提供，每个控制点包含原图和目标图的坐标。

可能的问题点：如何确保变换矩阵的正确性，尤其是在处理多项式的时候。比如，一阶变换的模型是x' = a0 + a1x + a2y，y' = b0 + b1x + b2y。这需要为x和y分别构建系数矩阵，然后求解。或者，是否应该将x和y分开处理，或者一起处理？通常，每个坐标方向的变换是独立的，所以可能需要分别求解a和b的系数。

例如，对于控制点，我们有多个点，每个点有原图的坐标（x, y）和目标坐标（x', y'）。对于一阶多项式，每个坐标的方程是x' = a0 + a1x + a2y。每个控制点可以构成一个方程，所有的控制点组成一个方程组，用最小二乘法求解a0, a1, a2。同样地，对y'方向也是一样，得到b0, b1, b2。所以，可能需要分别求解x和y的系数。

那在代码中，我需要为x和y分别构建系数矩阵。例如，对于一阶多项式，每个控制点的系数矩阵的行是[1, x, y]，对应的目标值是x'。同样，对于y方向，目标值是y'。然后使用numpy.linalg.lstsq来解这个方程组。

接下来是反向变换的问题。因为通常几何校正是将校正后的图像每个像素映射回原图的位置，所以需要逆变换。不过用户提供的控制点可能是从原图到目标图的变换，所以在进行校正时，可能需要逆变换的参数。或者，假设用户提供的控制点是原图到目标图的对应关系，那么当我们要生成目标图像时，每个目标图像的像素点需要找到对应的原图位置，这时候可能需要使用逆变换。不过这可能比较复杂，特别是当变换矩阵不是可逆的时候。或者，是否应该直接使用原变换的正向变换，将原图的点映射到目标图？这可能会产生问题，因为目标图的像素可能不是整数，或者存在空洞。因此，反向映射更为常见，即遍历每个目标图像的像素，通过逆变换找到在原图中的位置，然后进行插值。

所以在代码中，可能需要先确定变换模型，然后计算逆变换。或者，用户提供的控制点是从目标图到原图的对应关系？或者需要明确控制点的方向？

这里可能需要用户提供的控制点对是原图坐标到目标图坐标的映射。例如，原图中的点(x_src, y_src)对应目标图中的点(x_dst, y_dst)。然后，要构建从目标到原图的变换，因为在校正的时候，我们需要对于每个目标图像的位置，找到对应的原图位置，并采样。或者，可能用户提供的控制点是目标到原图的？这需要明确。

或者，可能用户提供的控制点对是原图到参考（目标）坐标的，即原图的点对应到目标坐标的点。因此，要构建一个从目标坐标到原图坐标的变换，以便将目标图像的每个像素映射回原图，然后进行采样。例如，当我们校正时，生成的目标图像中的每个点(x', y')需要找到对应的原图中的(x, y)，这可能需要逆变换，或者构建另一个方向的变换模型。

这里可能需要重新考虑控制点的使用方式。比如，假设我们有原图的坐标点（比如，存在几何畸变的图像）和它们对应的正确的地理坐标（或目标坐标）。几何校正的目标是将原图校正到目标坐标系中。因此，需要建立一个从目标坐标到原图坐标的映射，因为对于每个目标坐标的像素，需要知道在原图中的哪个位置获取像素值。因此，可能控制点对应该是目标坐标（x_dst, y_dst）与原图中的坐标（x_src, y_src）之间的对应关系。因此，我们需要建立一个从目标坐标到原图坐标的变换模型。这可能与用户提供的控制点的顺序有关，可能需要用户以这样的方式提供：即每个控制点包含目标坐标和对应的原图坐标。例如，控制点可能是这样的： [(x_dst1, y_dst1, x_src1, y_src1), ...]。或者可能反过来？

这可能会影响如何构建变换模型。例如，如果用户提供的控制点是原图的坐标到目标坐标的映射，即原图中的点(x_src, y_src)应该被映射到目标坐标(x_dst, y_dst)。那么，我们需要构建一个变换模型，将目标坐标映射回原图坐标。这可能比较复杂，因为变换可能是不可逆的。或者，用户可能需要用目标到原图的控制点对，这样可以直接构建从目标坐标到原图坐标的变换，这样就不需要求逆。

例如，假设用户提供的每个控制点是（dst_x, dst_y, src_x, src_y），这样我们直接可以用dst坐标作为自变量，src坐标作为因变量，建立模型。这样，在生成目标图像时，每个目标像素的位置（dst_x, dst_y）可以通过模型计算出对应的src_x和src_y，从而在原图中采样。这样更直接，不需要求逆变换。

所以，用户可能需要将控制点以这种方式提供。因此，在构建多项式模型时，自变量是目标坐标的x和y，因变量是原图的x和y。比如，原图的x坐标可以表示为多项式函数，参数是目标坐标的x和y。同样，原图的y坐标也是如此。

因此，构建两个多项式模型：一个用于预测原图的x坐标（基于目标坐标的x和y），另一个用于预测原图的y坐标。这样，当处理每个目标像素的坐标时，带入多项式模型，得到对应的原图坐标，然后进行插值。

这可能更合理，因此用户提供的控制点应该是目标坐标到原图坐标的对应关系。这可能更直接，因为几何校正的目的是将原图校正到目标坐标系，所以对于每个目标像素的位置，我们需要找到在原图中的对应位置。

这样的话，控制点的结构应该是：对于每个控制点，给出其在目标图像中的位置（比如，正确的坐标）和其在原图中的位置（可能畸变的坐标）。例如，可能用户知道某些地物在目标坐标系中的正确位置，以及这些地物在原图中的位置，从而建立对应关系。

现在，构建模型的话，比如，使用二阶多项式模型，原图的x和y坐标可以表示为：

src_x = a0 + a1*dst_x + a2*dst_y + a3*dst_x^2 + a4*dst_x*dst_y + a5*dst_y^2

src_y = b0 + b1*dst_x + b2*dst_y + b3*dst_x^2 + b4*dst_x*dst_y + b5*dst_y^2

对于每个控制点，我们可以建立两个方程：一个是关于src_x的，另一个是src_y的。如果有足够的控制点，可以用最小二乘法求解这些系数。

然后，在应用变换时，对于目标图像中的每个像素点（dst_x, dst_y），计算对应的src_x和src_y，然后从原图中采样（使用双线性插值）得到该位置的像素值，赋给目标图像。

现在，如何将这些步骤转化为代码？

首先，读取原图像，可以使用Pillow库，将图像转换为numpy数组。然后，定义控制点列表，每个元素包含目标坐标和原图坐标。例如，control_points的结构可能是[(dst_x1, dst_y1, src_x1, src_y1), ...]。

然后，构建多项式矩阵。例如，对于二阶多项式，每个控制点对应一个行向量，包含各项的值。比如，对于dst_x和dst_y，生成的行向量是[1, dst_x, dst_y, dst_x^2, dst_x*dst_y, dst_y^2]。将这些行向量组成矩阵A，然后解方程组A * coeffs_x = src_x_vector，同理coeffs_y = ...。使用numpy.linalg.lstsq来解。

得到coeffs_x和coeffs_y后，就可以对每个目标坐标点计算对应的原图坐标。

然后，生成目标图像，遍历每个像素的坐标（dst_x, dst_y），计算对应的src_x和src_y。然后，根据这些坐标在原图中进行插值。

双线性插值的实现需要处理原图中的四个邻近像素，如果坐标超出原图范围，则填充0或者进行其他处理。

现在，将这个过程写成函数。可能需要以下步骤：

1. 读取原图像，转换为numpy数组。

2. 定义控制点。

3. 根据控制点构建多项式模型，求解系数。

4. 创建目标图像数组，尺寸可能需要用户指定，或者根据原图尺寸和变换自动计算。

5. 遍历目标图像的每个像素，计算对应的原图坐标，进行插值。

6. 将结果保存为PNG图像。

可能的函数结构：

- 函数geometric_correction，参数为原图像路径，控制点列表，多项式阶数，目标图像尺寸，输出路径等。

- 辅助函数，如计算多项式系数，双线性插值等。

需要注意的问题：

- 控制点的数量必须足够。比如，一阶多项式需要至少3个点，二阶需要至少6个点，否则无法解唯一解，或者需要使用最小二乘法当点数多于所需时。

- 目标图像的尺寸可能需要用户指定，或者自动计算为原图的大小，或者根据控制点的范围确定。

- 处理图像时，原图的坐标可能是浮点数，需要插值。

- 颜色通道的处理，比如RGB图像需要分别处理每个通道。

例如，原图是三维数组（高度，宽度，通道），所以插值时需要对每个通道进行操作。

现在，编写代码的大致步骤：

首先，导入必要的库：PIL的Image，numpy，以及math。

然后，定义函数：

def geometric_correction(src_image_path, control_points, order=2, output_size=None, output_path='corrected.png'):

# 步骤1：读取原图像

src_img = Image.open(src_image_path)

src_array = np.array(src_img)

height_src, width_src = src_img.size  # 注意PIL的size是（宽，高）

channels = src_array.shape[2] if len(src_array.shape) ==3 else 1

# 步骤2：处理控制点，构建多项式方程组

# 控制点格式应为 [(dst_x, dst_y, src_x, src_y), ...]

# 提取目标坐标和原坐标

dst_points = np.array([(pt[0], pt[1]) for pt in control_points])

src_x = np.array([pt[2] for pt in control_points])

src_y = np.array([pt[3] for pt in control_points])

# 构建多项式矩阵A

# 根据阶数生成对应的项

A = []

for (x, y) in dst_points:

terms = []

if order >= 1:

# 一阶项：1, x, y

terms = [1, x, y]

if order >=2:

# 二阶项：x², xy, y²

terms += [x**2, x*y, y**2]

if order >=3:

# 三阶及以上可能需要更多项，这里假设只处理1和2阶

pass # 暂时不支持

A.append(terms)

A = np.array(A)

# 解最小二乘问题

coeffs_x, residuals, rank, singular_values = np.linalg.lstsq(A, src_x, rcond=None)

coeffs_y, residuals, rank, singular_values = np.linalg.lstsq(A, src_y, rcond=None)

# 步骤3：确定目标图像尺寸

if output_size is None:

# 可能需要根据控制点的dst坐标的最大最小值来确定

dst_x_coords = [pt[0] for pt in control_points]

dst_y_coords = [pt[1] for pt in control_points]

min_x, max_x = min(dst_x_coords), max(dst_x_coords)

min_y, max_y = min(dst_y_coords), max(dst_y_coords)

output_width = int(max_x - min_x + 1)

output_height = int(max_y - min_y + 1)

# 或者直接使用原图的尺寸？

# 这里假设用户指定或者使用原图尺寸

# 可能需要更合理的处理，比如根据变换后的原图范围

# 但为了简单，可能直接使用原图尺寸？

output_width, output_height = src_img.size

else:

output_width, output_height = output_size

# 创建目标图像数组

if channels ==1:

corrected_array = np.zeros((output_height, output_width), dtype=np.uint8)

else:

corrected_array = np.zeros((output_height, output_width, channels), dtype=np.uint8)

# 步骤4：遍历目标图像的每个像素，计算对应的原图坐标，并插值

for dst_y in range(output_height):

for dst_x in range(output_width):

# 计算多项式项

terms = []

x = dst_x

y = dst_y

if order >=1:

terms = [1, x, y]

if order >=2:

terms += [x**2, x*y, y**2]

terms = np.array(terms)

# 计算对应的原图坐标

src_x = np.dot(terms, coeffs_x)

src_y = np.dot(terms, coeffs_y)

# 进行双线性插值

# 需要处理src_x和src_y是否在原图范围内

x0 = int(np.floor(src_x))

y0 = int(np.floor(src_y))

x1 = x0 +1

y1 = y0 +1

# 边界检查

if x0 <0 or y0 <0 or x1 >= width_src or y1 >= height_src:

# 超出范围，设置为0或透明

pixel_value = 0

else:

# 双线性插值

dx = src_x - x0

dy = src_y - y0

for c in range(channels) if channels !=1 else [0]:

# 获取四个邻域的像素值

top_left = src_array[y0, x0, c] if channels !=1 else src_array[y0, x0]

top_right = src_array[y0, x1, c] if channels !=1 else src_array[y0, x1]

bottom_left = src_array[y1, x0, c] if channels !=1 else src_array[y1, x0]

bottom_right = src_array[y1, x1, c] if channels !=1 else src_array[y1, x1]

# 插值计算

top = top_left * (1 - dx) + top_right * dx

bottom = bottom_left * (1 - dx) + bottom_right * dx

pixel_value = top * (1 - dy) + bottom * dy

pixel_value = int(np.clip(pixel_value, 0, 255))

if channels ==1:

corrected_array[dst_y, dst_x] = pixel_value

else:

corrected_array[dst_y, dst_x, c] = pixel_value

# 保存图像

corrected_img = Image.fromarray(corrected_array)

corrected_img.save(output_path)

但这段代码可能存在问题：

1. PIL的Image.open读取的图像，转换为numpy数组后的形状是（高度，宽度，通道），而PIL的size是（宽度，高度）。所以在处理的时候需要将坐标转换正确。例如，原图的宽度是src_img.size[0]，高度是src_img.size[1]。在numpy数组中，是src_array[y, x]。

所以在代码中，src_array的形状应该是（高度，宽度， 通道）或者（高度， 宽度）对于灰度图。例如，当原图是RGB时，src_array的形状是 (height, width, 3)。

所以在遍历的时候，原图的坐标是src_x和src_y，对应到数组中的位置是 [y, x]。因此，在双线性插值时，需要注意x和y的顺序是否与数组的索引对应。

例如，假设原图的某个点的坐标为 (src_x, src_y)，在数组中对应的位置是 [y, x]。例如，src_x是列数，对应到数组的第二个轴，src_y是行数，对应到数组的第一个轴。这可能容易出错。

例如，如果原图是宽度为W，高度为H，那么数组的索引是 [y][x]，其中x的范围是0到W-1，y的范围是0到H-1。所以，当计算得到原图的坐标为 (src_x, src_y)，需要转换为数组的索引为 (int(src_y), int(src_x))。或者，可能坐标的x对应数组的列（第二个轴），y对应数组的行（第一个轴）。

这可能是一个容易出错的地方。例如，在代码中，当用户提供的控制点中的src_x和src_y是原图中的坐标，比如原图的左上角是 (0,0)，向右为x增加，向下为y增加。那么在数组中，对应的是 [y][x]。所以在处理时，比如，原图的坐标是 (x, y)，对应的数组索引是 [y][x]。例如，src_x的范围是 0到width_src-1，src_y的范围是0到height_src-1。

所以，在双线性插值时，原图的坐标需要转换为数组的索引，即：

x0 = floor(src_x)

y0 = floor(src_y)

x1 = x0 + 1

y1 = y0 + 1

然后，检查这些x0, x1是否在 [0, width_src -1] 范围内，y0, y1是否在 [0, height_src -1] 范围内。如果越界，则处理为0或者跳过。

但是，在代码中，原图的宽度和高度是从src_img.size得到的，而src_img.size是 (width, height)。所以在数组中的形状是 (height, width, channels)。因此，在访问数组时，必须确保y0和y1不超过height-1，x0和x1不超过width-1。

例如，假设原图的宽度是width_src，高度是height_src。数组的shape是 (height_src, width_src, ...)。所以在检查时，应确保x0 >=0，x1 < width_src，y0 >=0，y1 < height_src。或者，在双线性插值的时候，如果坐标超出范围，就直接返回0。

另一个问题，当处理RGB图像时，每个像素有三个通道，所以在双线性插值时，必须对每个通道进行插值。例如，在代码中，可能需要循环每个通道，或者使用向量化操作。

此外，在计算多项式项的时候，代码中的dst_x和dst_y是否是正确的目标坐标？比如，假设目标图像的坐标从0开始到output_width-1和output_height-1。例如，当output_size设置为原图的大小时，那么遍历的dst_x的范围是0到width_src-1，dst_y是0到height_src-1。这可能不正确，因为控制点中的目标坐标可能是一个不同的范围。例如，用户可能希望将原图校正到一个更大的坐标系中，或者不同的位置。因此，这里可能需要更灵活的处理，比如允许用户指定输出图像的范围，或者根据控制点的目标坐标自动计算。

这可能比较复杂，但为了示例，可能让用户指定output_size，或者默认使用原图的尺寸。例如，在示例中，用户可能希望校正后的图像与原图尺寸相同，或者根据控制点的范围调整。

但为了简化问题，可能在示例中将output_size设置为原图的尺寸，或者根据控制点的最大坐标自动计算。例如，如果控制点中的目标坐标是相对于校正后的图像的左上角，那么可能需要确定校正后的图像的范围。

这可能超出示例的范围，所以可能假设output_size由用户提供，或者在函数中默认为原图的尺寸。

另一个问题是，当多项式阶数较高时，比如二阶，需要更多的控制点。如果用户提供的控制点不足，numpy.linalg.lstsq会给出一个最小二乘解，但可能不准确。所以函数需要确保控制点的数量足够，比如对于二阶，至少需要6个点。否则，解可能不可靠。

现在，关于双线性插值的实现，代码中的循环是否正确？比如，对于每个通道，单独处理。或者，是否可以利用numpy的向量化操作来加速？

在示例代码中，为了清晰，可能会逐个像素处理，但在实际中这会非常慢，特别是对大图像。不过，用户可能只是需要一个演示，所以可以接受。

此外，双线性插值的计算是否正确？例如，dx应该是src_x - x0，而dy是src_y - y0。然后，top是左右之间的线性插值，bottom同理，最后是上下之间的线性插值。这部分公式是否正确？

是的，正确的双线性插值公式应该是这样。

此外，当x0或y0超出原图范围时，如何处理？在示例代码中，如果x0 <0或y0 <0，或者x1 >= width_src或y1 >= height_src，那么该像素会被设为0。这可能是一种简单的处理方式，但在实际中可能需要其他方法，如镜像、重复边缘像素等。但为了简单，示例可能采用填0的方式。

现在，测试这个函数是否可行。例如，假设用户有一个简单的例子，比如原图是某种几何畸变的图像，用户提供几个控制点，然后运行函数生成校正后的图像。

例如，假设原图是一个被水平错切的图像，用户提供四个控制点，校正后的图像应该是一个矩形。控制点可能是在目标图像中的四个角点，对应的原图中的错切后的坐标。

但如何生成这样的示例呢？可能用户需要手动定义控制点对。

例如，假设原图被向右错切，每个点的x坐标增加了y坐标的某个比例。例如，原图中的点 (x, y) 在目标图像中应为 (x - 0.5*y, y)。这可能是一个一阶变换。用户提供的控制点可以是四个角点：

例如，原图的四个角点可能被错切成：

目标坐标（0,0）对应原图（0,0）

目标坐标（width,0）对应原图（width, 0)

目标坐标（0,height）对应原图（0.5*height, height)

目标坐标（width,height）对应原图(width + 0.5*height, height)

这样的话，控制点可能帮助计算出一阶变换的系数。

但用户可能需要更多的控制点，比如至少三个点，因为一阶多项式有6个系数（每个坐标方向3个）。或者，可能这里的目标到原图的变换模型是一阶，每个坐标方向需要三个参数，因此需要至少三个控制点。

例如，假设用户提供的控制点足够，函数将正确计算变换参数，并应用校正。

现在，编写一个示例的使用代码：

假设原图是一个带有错切的图像，用户希望校正它。例如，原图的左上角在目标图像的（0,0），原图的右上角在目标图像的（width,0），原图的左下角在目标图像的（0,height），而原图中的左下角实际位于（0, height）可能被错切成（delta, height）。因此，用户提供的控制点可能包括：

control_points = [

(0, 0, 0, 0),          # 目标(0,0) 对应原图(0,0)

(width, 0, width, 0),  # 目标(width,0) 对应原图(width,0)

(0, height, delta, height)  # 目标(0,height)对应原图(delta, height)

]

这里，假设delta是某个值，比如原图在y方向有水平错切。例如，假设delta = 50，当height=100时。

此时，使用一阶多项式，构建变换模型，其中原图的x坐标由目标x和y的线性组合得到，例如x_src = a0 + a1*dst_x + a2*dst_y。而根据这三个控制点：

对于点 (0,0)，x_src=0 = a0 + a1*0 + a2*0 → a0=0.

对于点 (width,0)，x_src=width = a0 + a1*width + a2*0 → width = a0 + a1*width. 因为a0=0，所以a1=1.

对于点 (0,height)，x_src=delta = a0 + a1*0 + a2*height → delta = a2*height → a2=delta/height.

所以，原图的x_src = dst_x + (delta/height)*dst_y.

这样，当应用这个变换时，校正后的图像每个点的dst_x和dst_y将被转换为原图的x_src = dst_x + (delta/height)*dst_y，而y_src = dst_y.

这样，在目标图像中的每个点（dst_x,dst_y）对应的原图坐标是（dst_x + (delta/height)*dst_y, dst_y）。这样，校正后的图像将消除错切，因为原图中的错切被逆变换了。

这可能是一个合理的测试案例。

在代码中，用户需要将这些控制点传入函数，然后生成校正后的图像。

综上，这个函数的大体结构应该能工作，但需要注意坐标系的转换，以及PIL和numpy数组的索引顺序是否正确。

最后，用户可能希望这个函数能帮助他们理解几何校正的本质，即通过控制点建立坐标映射关系，然后利用插值填充校正后的图像。这个过程的核心是坐标变换模型的构建和重采样，这正是函数所要演示的。

代码

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
# @Time : 2025/2/12 

'''

'''
import numpy as np
from PIL import Image
import math

def bilinear_interpolation(src_array, x, y):
    """双线性插值实现"""
    x0, y0 = math.floor(x), math.floor(y)
    x1, y1 = x0 + 1, y0 + 1

    # 边界检查
    height, width = src_array.shape[:2]
    if not (0 <= x0 < width and 0 <= x1 < width and
            0 <= y0 < height and 0 <= y1 < height):
        return 0

    # 计算权重
    dx, dy = x - x0, y - y0

    # 获取四个邻域像素值
    if len(src_array.shape) == 3:  # 彩色图像
        top_left = src_array[y0, x0, :]
        top_right = src_array[y0, x1, :]
        bottom_left = src_array[y1, x0, :]
        bottom_right = src_array[y1, x1, :]
    else:  # 灰度图像
        top_left = src_array[y0, x0]
        top_right = src_array[y0, x1]
        bottom_left = src_array[y1, x0]
        bottom_right = src_array[y1, x1]

    # 插值计算
    top = top_left * (1 - dx) + top_right * dx
    bottom = bottom_left * (1 - dx) + bottom_right * dx
    return top * (1 - dy) + bottom * dy


def polynomial_coefficients(control_points, order):
    """计算多项式变换系数"""
    dst_points = np.array([[pt[0], pt[1]] for pt in control_points])
    src_x = np.array([pt[2] for pt in control_points])
    src_y = np.array([pt[3] for pt in control_points])

    # 构建多项式矩阵
    A = []
    for (x, y) in dst_points:
        terms = [1]
        if order >= 1:  # 一阶项
            terms += [x, y]
        if order >= 2:  # 二阶项
            terms += [x ** 2, x * y, y ** 2]
        A.append(terms)

    A = np.array(A)
    coeffs_x = np.linalg.lstsq(A, src_x, rcond=None)[0]
    coeffs_y = np.linalg.lstsq(A, src_y, rcond=None)[0]
    return coeffs_x, coeffs_y


def geometric_correction(src_path, control_points, order=1, output_size=None):
    """
    几何校正主函数
    :param src_path: 原始图像路径
    :param control_points: 控制点列表，格式[(dst_x, dst_y, src_x, src_y)]
    :param order: 多项式阶数（1或2）
    :param output_size: 输出图像尺寸（宽，高）
    """
    # 读取原始影像
    src_img = Image.open(src_path)
    src_array = np.array(src_img)
    src_height, src_width = src_array.shape[:2]

    # 计算变换系数
    coeffs_x, coeffs_y = polynomial_coefficients(control_points, order)

    # 确定输出尺寸
    if output_size is None:
        output_width, output_height = src_img.size
    else:
        output_width, output_height = output_size

    # 创建输出图像
    if len(src_array.shape) == 3:
        corrected = np.zeros((output_height, output_width, 3), dtype=np.uint8)
    else:
        corrected = np.zeros((output_height, output_width), dtype=np.uint8)

    # 遍历输出图像的每个像素
    for y in range(output_height):
        for x in range(output_width):
            # 计算对应的原图坐标
            terms = [1]
            if order >= 1:
                terms += [x, y]
            if order >= 2:
                terms += [x ** 2, x * y, y ** 2]

            src_x = np.dot(terms, coeffs_x)
            src_y = np.dot(terms, coeffs_y)

            # 执行插值
            src_array = src_array[:,:,:3]
            if 0 <= src_x < src_width and 0 <= src_y < src_height:
                corrected[y, x] = bilinear_interpolation(src_array, src_x, src_y)

    return Image.fromarray(corrected)


# 示例使用
if __name__ == "__main__":
    # 定义控制点示例（目标坐标X，目标坐标Y，原图坐标X，原图坐标Y）
    control_points = [
        (0, 0, 0, 50),  # 左上角
        (300, 0, 250, 30),  # 右上角
        (0, 300, 40, 250),  # 左下角
        (600, 500, 560, 570)  # 右下角
    ]

    # 执行几何校正
    corrected_img = geometric_correction(
        "11.png",
        control_points,
        order=1,
        output_size=(1500, 1500)
    )
    corrected_img.save("corrected.png")

测试