定义二维网格的算法

Question

假设一个网格由一组网格参数定义：它的原点 (x0, y0)，一边和 x 轴之间的角度，增量和 - 请看下图

网格上有已知坐标的散点，但并不完全落在网格交点上。有没有一种算法可以找到一组网格参数来定义网格，使点最适合网格交点？

假设已知坐标为： (2, 5.464), (3.732, 6.464), (5.464, 7.464) (3, 3.732), (4.732, 4.732), (6.464, 5.732) (4, 2), (5.732, 3), (7.464, 4)。我希望算法找到原点 (4, 2)，角度 30 度，增量都是 2.

Answer 1

你可以通过找到一个矩阵来解决这个问题，该矩阵将点从位置 (0, 0), (0, 1), ... (2, 2) 转换到给定的点上。

虽然网格只有5个自由度（原点位置+角度+尺度），但使用2x3矩阵A定义变换更容易，因为在这种情况下问题可以线性化。

设索引为(x0, y0)的点在网格上变换为点(x0', y0')，例如(0, 0) -> (2, 5.464) 设a_ij 是矩阵 A 的系数。然后这对点产生 2 个方程：

a_00 * x0 + a_01 * y0 + a_02 = x0'
a_10 * x0 + a_11 * y0 + a_12 = y0'

未知数是a_ij，所以这些方程可以写成

形式

a_00 * x0 + a_01 * y0 + a_02 + a_10 * 0 + a_11 * 0 + a_12 * 0 = x0'
a_00 * 0 + a_01 * 0 + a_02 * 0 + a_10 * x0 + a_11 * y0 + a_12 = y0'

或矩阵形式

K0 * (a_00, a_01, a_02, a_10, a_11, a_12)^T = (x0', y0')^T

哪里

K0 = (
    x0, y0, 1,  0,  0,  0
    0,  0,  0,  x0, y0, 1
)

每对点的这些等式可以组合成一个等式

K * (a_00, a_01, a_02, a_10, a_11, a_12)^T = (x0', y0', x1', y1', ..., xn', yn')^T

或K * a = b 其中

K = (
    x0, y0, 1,  0,  0,  0
    0,  0,  0,  x0, y0, 1
    x1, y1, 1,  0,  0,  0
    0,  0,  0,  x1, y1, 1
    ...
    xn, yn, 1,  0,  0,  0
    0,  0,  0,  xn, yn, 1
)

和(xi, yi), (xi', yi')是一对对应点

这可以作为非齐次线性方程组求解。在这种情况下，解决方案将最小化从每个点到最近的网格交点的距离的平方和。假设点从具有正态分布噪声的网格交叉点移动，也可以考虑此变换以最大化整体可能性。

a = (K^T * K)^-1 * K^T * b

如果有线性代数库，这个算法很容易实现。以下是 Python:

中的示例

import numpy as np

n_points = 9
aligned_points = [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)]
grid_points = [(2, 5.464), (3.732, 6.464), (5.464, 7.464), (3, 3.732), (4.732, 4.732), (6.464, 5.732), (4, 2), (5.732, 3), (7.464, 4)]

K = np.zeros((n_points * 2, 6))
b = np.zeros(n_points * 2)
for i in range(n_points):
    K[i * 2, 0] = aligned_points[i, 0]
    K[i * 2, 1] = aligned_points[i, 1]
    K[i * 2, 2] = 1
    K[i * 2 + 1, 3] = aligned_points[i, 0]
    K[i * 2 + 1, 4] = aligned_points[i, 1]
    K[i * 2 + 1, 5] = 1
    
    b[i * 2] = grid_points[i, 0]
    b[i * 2 + 1] = grid_points[i, 1]
    
# operator '@' is matrix multiplication
a = np.linalg.inv(np.transpose(K) @ K) @ np.transpose(K) @ b
A = a.reshape(2, 3)
print(A)

[[ 1.     1.732  2.   ]
 [-1.732  1.     5.464]]

然后可以从这个矩阵中提取参数：

theta = math.degrees(math.atan2(A[1, 0], A[0, 0]))
scale_x = math.sqrt(A[1, 0] ** 2 + A[0, 0] ** 2)
scale_y = math.sqrt(A[1, 1] ** 2 + A[0, 1] ** 2)
origin_x = A[0, 2]
origin_y = A[1, 2]

theta = -59.99927221917264
scale_x = 1.99995599951599
scale_y = 1.9999559995159895
origin_x = 1.9999999999999993
origin_y = 5.464

不过还有一个小问题：矩阵A对应的是仿射变换。这意味着不能保证网格轴是垂直的。如果这是一个问题，那么可以修改矩阵的前两列，使变换保留角度。

Answer 2

更新： 我修正了错误并解决了符号歧义，所以现在这个算法产生了预期的结果。然而，应该测试它是否所有情况都被正确处理。

这是解决这个问题的另一种尝试。这个想法是将变换分解为非均匀缩放矩阵和旋转矩阵A = R * S，然后在给定r1^2 + r2^2 = 1 的限制条件下求解这些矩阵的系数sx, sy, r1, r2。此处描述了最小化问题：How to find a transformation (non-uniform scaling and similarity) that maps one set of points to another?

def shift_points(points):
    n_points = len(points)
    shift = tuple(sum(coords) / n_points for coords in zip(*points))
    shifted_points = [(point[0] - shift[0], point[1] - shift[1]) for point in points]
    return shifted_points, shift

n_points = 9
aligned_points = [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)]
grid_points = [(2, 5.464), (3.732, 6.464), (5.464, 7.464), (3, 3.732), (4.732, 4.732), (6.464, 5.732), (4, 2), (5.732, 3), (7.464, 4)]

aligned_points, aligned_shift = shift_points(aligned_points)
grid_points, grid_shift = shift_points(grid_points)
c1, c2 = 0, 0
b11, b12, b21, b22 = 0, 0, 0, 0
for i in range(n_points):
    c1 += aligned_points[i][0] ** 2
    c2 += aligned_points[i][0] ** 2
    b11 -= 2 * aligned_points[i][0] * grid_points[i][0]
    b12 -= 2 * aligned_points[i][1] * grid_points[i][0]
    b21 -= 2 * aligned_points[i][0] * grid_points[i][1]
    b22 -= 2 * aligned_points[i][1] * grid_points[i][1]
k = (b11 ** 2 * c2 + b22 ** 2 * c1 - b21 ** 2 * c2 - b12 ** 2 * c1) / \
    (b21 * b11 * c2 - b12 * b22 * c1)
# r1_sqr and r2_sqr might need to be swapped
r1_sqr = 2 / (k ** 2 + 4 + k * math.sqrt(k ** 2 + 4))
r2_sqr = 2 / (k ** 2 + 4 - k * math.sqrt(k ** 2 + 4))
for sign1, sign2 in [(1, 1), (-1, 1), (1, -1), (-1, -1)]:
    r1 = sign1 * math.sqrt(r1_sqr)
    r2 = sign2 * math.sqrt(r2_sqr)
    scale_x = -b11 / (2 * c1) * r1 - b21 / (2 * c1) * r2
    scale_y = b12 / (2 * c2) * r2 - b22 / (2 * c2) * r1
    if scale_x >= 0 and scale_y >= 0:
        break
theta = math.degrees(math.atan2(r2, r1))

r1_sqr和r2_sqr的选择可能有歧义。可以从aligned_shift和grid_shift估计原点，但我还没有实现。

theta = -59.99927221917264
scale_x = 1.9999559995159895
scale_y = 1.9999559995159895

定义二维网格的算法

Algorithm to define a 2d grid

optimization

inversion