给定三点计算仿射变换

Question

我有两张图像，并使用 sift 找到了三个相似的 2D 点。我需要计算图像之间的仿射变换。不幸的是，我错过了讲座，而且那里的信息对我来说有点密集。计算这个 2x3 矩阵的一般方法是什么？

我有一个 2x3 矩阵 [x1 y1;x2 y2;x3 y3] 中的点矩阵，但我从那里迷路了。 感谢您的帮助。

Answer 1

通常，二维点的仿射变换表示为

x' = A*x

其中x 是原始二维位置的三向量[x; y; 1]，x' 是变换点。仿射矩阵A是

A = [a11 a12 a13;
     a21 a22 a23;
       0   0   1]

当x 和A已知并且您希望恢复x' 时，此表单很有用。

但是，您可以用不同的方式表达这种关系。 让

X = [xi yi 1  0  0  0;
      0  0 0 xi yi  1 ]

而a 是一个列向量

a = [a11; a12; a13; a21; a22; a23]

然后

X*a = [xi'; yi']

适用于所有对应点对x_i, x_i'。

当您知道点对之间的对应关系并希望恢复 A 的参数时，这种替代形式非常有用。
将所有点堆叠在一个大矩阵X（每个点两行）中，您将拥有 2*n-by-6 矩阵 X 乘以 6 未知向量 a 等于 2*n-by -1 堆叠对应点的列向量（记为x_prime）：

X*a = x_prime

求解a：

a = X \ x_prime

以最小二乘的方式恢复a 的参数。

祝你好运，别逃课了！

Answer 2

很抱歉没有使用 Matlab，但我只使用了一点 Python。我认为这段代码可以帮助你（抱歉代码风格不好——我是数学家，不是程序员）

import numpy as np
# input data
ins = [[1, 1], [2, 3], [3, 2]]  # <- points
out = [[0, 2], [1, 2], [-2, -1]] # <- mapped to
# calculations
l = len(ins)
B = np.vstack([np.transpose(ins), np.ones(l)])
D = 1.0 / np.linalg.det(B)
entry = lambda r,d: np.linalg.det(np.delete(np.vstack([r, B]), (d+1), axis=0))
M = [[(-1)**i * D * entry(R, i) for i in range(l)] for R in np.transpose(out)]
A, t = np.hsplit(np.array(M), [l-1])
t = np.transpose(t)[0]
# output
print("Affine transformation matrix:\n", A)
print("Affine transformation translation vector:\n", t)
# unittests
print("TESTING:")
for p, P in zip(np.array(ins), np.array(out)):
  image_p = np.dot(A, p) + t
  result = "[OK]" if np.allclose(image_p, P) else "[ERROR]"
  print(p, " mapped to: ", image_p, " ; expected: ", P, result)

这段代码演示了如何将仿射变换恢复为矩阵和向量，并测试初始点是否映射到它们应该映射到的位置。您可以使用Google colab 测试此代码，因此您无需安装任何东西。也许，你可以把它翻译成 Matlab。

关于此代码背后的理论：它基于“Beginner's guide to mapping simplexes affinely”中提出的方程，矩阵恢复在“规范符号的恢复”一节中描述。同一作者发表了“Workbook on mapping simplexes affinely”，其中包含许多此类实际示例。

Answer 3

为了在 OpenCV 中实现，您可以使用 cv2.getAffineTransform

getAffineTransform() [1/2]

Mat cv::getAffineTransform  (   const Point2f   src[],
    const Point2f   dst[] 
)   

Python:

cv.getAffineTransform(  src, dst    ) ->    retval

Answer 4

这适用于任何想要在 C 中执行此操作的人。该算法基于此 MathStackExchange post。作者展示了如何使用 Gauss-Jordan 消元法求仿射变换系数的公式，

/*
*Function: Affine Solver
*Role: Finds Affine Transforming mapping (X,Y) to (X',Y')
*Input: double array[A,B,C,D,E,F], 
*   int array[X-Coordinates], int array[Y-Coordinates],
*   int array[X'-Coordinates],int array[Y'-Coordinates]
*Output:void - Fills double array[A,B,C,D,E,F]
*/

void AffineSolver(double *AtoF, int *X, int *Y, int *XP, int *YP)
{
    AtoF[0] = (double)(XP[1]*Y[0] - XP[2]*Y[0] - XP[0]*Y[1] + XP[2]*Y[1] + XP[0]*Y[2] -XP[1]*Y[2]) / 
          (double)(X[1]*Y[0] - X[2]*Y[0] - X[0]*Y[1] + X[2]*Y[1] + X[0]*Y[2] - X[1]*Y[2]);
          
    AtoF[1] = (double)(XP[1]*X[0] - XP[2]*X[0] - XP[0]*X[1] + XP[2]*X[1] + XP[0]*X[2] -XP[1]*X[2]) / 
          (double)(-X[1]*Y[0] + X[2]*Y[0] + X[0]*Y[1] - X[2]*Y[1] - X[0]*Y[2] +X[1]*Y[2]);
          
    AtoF[2] = (double)(YP[1]*Y[0] - YP[2]*Y[0] - YP[0]*Y[1] + YP[2]*Y[1] + YP[0]*Y[2] -YP[1]*Y[2]) / 
          (double)(X[1]*Y[0] - X[2]*Y[0] - X[0]*Y[1] + X[2]*Y[1] + X[0]*Y[2] - X[1]*Y[2]);
          
    AtoF[3] = (double)(YP[1]*X[0] - YP[2]*X[0] - YP[0]*X[1] + YP[2]*X[1] + YP[0]*X[2] -YP[1]*X[2]) / 
          (double)(-X[1]*Y[0] + X[2]*Y[0] + X[0]*Y[1] - X[2]*Y[1] - X[0]*Y[2] +X[1]*Y[2]);
          
    AtoF[4] = (double)(XP[2]*X[1]*Y[0] - XP[1]*X[2]*Y[0]-XP[2]*X[0]*Y[1] + XP[0]*X[2]*Y[1]+
           XP[1]*X[0]*Y[2] - XP[0]*X[1]*Y[2]) / 
          (double)(X[1]*Y[0] - X[2]*Y[0] - X[0]*Y[1] + X[2]*Y[1] + X[0]*Y[2] - X[1]*Y[2]);
          
    AtoF[5] = (double)(YP[2]*X[1]*Y[0] - YP[1]*X[2]*Y[0]-YP[2]*X[0]*Y[1] + YP[0]*X[2]*Y[1] + YP[1]*X[0]*Y[2] -            YP[0]*X[1]*Y[2]) / 
          (double)(X[1]*Y[0] - X[2]*Y[0] - X[0]*Y[1] + X[2]*Y[1] + X[0]*Y[2] - X[1]*Y[2]);
}

/*
*Function: PrintMatrix
*Role: Prints 2*3 matrix as //a b e
                //c d f
*Input: double array[ABCDEF]
*Output: voids
*/

void PrintMatrix(double *AtoF)
{
    printf("a = %f ",AtoF[0]);
    printf("b = %f ",AtoF[1]);
    printf("e = %f\n",AtoF[4]);
    printf("c = %f ",AtoF[2]);
    printf("d = %f ",AtoF[3]);
    printf("f = %f ",AtoF[5]);
}

int main()
{
/*Test*/
/*Find transform mapping (0,10),(0,0),(10,0) to (0,5)(0,0)(5,0)*/
/*Expected Output*/
//a = 0.500000 b = 0.000000 e = -0.000000
//c = -0.000000 d = 0.500000 f = -0.000000
/*Test*/
    double *AtoF = calloc(6, sizeof(double));
    int  X[] = { 0, 0,10};
    int  Y[] = {10, 0, 0};
    int XP[] = { 0, 0, 5};
    int YP[] = { 5, 0, 0};
    AffineSolver(AtoF,X,Y,XP,YP);
    PrintMatrix(AtoF);
    free(AtoF);
    return 0;
}