19

我有对象让我们在模型图像上说。我想计算模型图像上的对象和目标图像上的对象之间的变换(位移、缩放、旋转)。我想假设对象可以被视为 2D,因此应该只计算 2D 转换。

首先,我想以手动辅助的方式进行。用户在模型图像上选择基点,然后在目标图像上选择目标点。点数应由用户定义(但不少于至少 2-3 个点)。当点给出不同的信息时,应该对转换进行平均,例如,可以从中计算匹配的质量。

所以问题是关于计算两组点的转换,但是因为我想在图像上做这件事,所以我添加了图像处理标签。

特别受欢迎的是带有一些代码或伪代码的参考和建议。

两点是很简单的问题,只需要线的旋转,比例和位移,但是如何用更多的点来做,并且平均它并计算一些质量因素。

目前的解决方案是:

void transformFnc(std::vector<PointF> basePoints, std::vector<PointF> targetPoints,
                  PointF& offset, double rotation, double scale)
{
    std::vector<Line> basePointsLines;
    std::vector<Line> targetPointsLines;
    assert(basePoints.size() == targetPoints.size());
    int pointsNumber = basePoints.size();

    for(int i = 0; i < pointsNumber; i++)
    {
         for(int j = i + 1; j < pointsNumber; j++)
         {
             basePointsLines.push_back(Line(basePoints[i], basePoints[j]));
             targetPointsLines.push_back(Line(targetPoints[i], targetPoints[j]));
         }
    }
    std::vector<double> scalesVector;
    std::vector<double> rotationsVector;
    double baseCenterX = 0, baseCenterY = 0, targetCenterX = 0, targetCenterY = 0;
    for(std::vector<Line>::iterator it = basePointsLines.begin(), i = targetPointsLines.begin();
        it != basePointsLines.end(), i != targetPointsLines.end(); it++, i++)
    {
        scalesVector.push_back((*i).length()/(*it).length());
        baseCenterX += (*it).pointAt(0.5).x(); 
        baseCenterY += (*it).pointAt(0.5).y();
        targetCenterX += (*i).pointAt(0.5).x();
        targetCenterY += (*i).pointAt(0.5).y();
        double rotation;
        rotation = (*i).angleTo((*it));
        rotationsVector.push_back(rotation);
    }
    baseCenterX = baseCenterX / pointsNumber;
    baseCenterY = baseCenterY / pointsNumber;
    targetCenterX = targetCenterX / pointsNumber;
    targetCenterY = targetCenterY / pointsNumber;

    offset = PointF(targetCenterX - baseCenterX, targetCenterY - baseCenterY);
    scale = sum(scalesVector) / scalesVector.size();
    rotation = sum(rotationsVector) / rotationsVector.size();
}

我能在这段代码中找到的唯一优化是从比例和旋转中消除那些与其他值相差太大的值。

我正在寻找解决方案命题的代码或伪代码。它也可以是一些代码的引用。

到目前为止,我知道的答案是:

  • 可以使用RANSAC算法
  • 我需要寻找一些最小二乘意义上的仿射变换计算算法
4

5 回答 5

27

首先用一个 3x3 仿射变换矩阵将问题概括为一个简单的仿射变换:即

[M11 M12 M13]
[M21 M22 M23]
[M31 M32 M33]

由于我们已经知道第三行总是 [0 0 1] 我们可以简单地忽略它。

现在我们可以将问题描述为以下矩阵方程

[xp0]     [x0 y0 1  0  0  0 ]
[yp0]     [0  0  0  x0 y0 1 ]     [M11]
[xp1]     [x1 y1 1  0  0  0 ]     [M12]
[yp1]  =  [0  0  0  x1 y1 1 ]  *  [M13]
[xp2]     [x2 y2 1  0  0  0 ]     [M21]
[yp2]     [0  0  0  x2 y2 1 ]     [M22]
[xp3]     [x3 y3 1  0  0  0 ]     [M23]
[yp3]     [0  0  0  x3 y3 1 ]

其中 xp 和 yp 是投影坐标,x 和 y 是原始坐标。

让我们称之为

proj = M * trans

然后,我们可以通过以下方式计算适合转换的最小二乘法

trans = pinv(M) * proj

其中 pinv 是伪逆。

这给了我们一个仿射变换,它最适合在最小二乘意义上给出的点。

现在显然这也会产生剪切、坐标翻转以及您不想要的非均匀缩放,因此我们需要以某种方式限制仿射变换以避免剪切。事实证明这很容易,我们可以使用单个向量来描述旋转(向量的方向)和缩放(向量的大小),另一个向量将与其正交。这将自由度减少了两个。

M21 = -M12
M22 = M11

所以减少到

[xp0]     [x0  y0 1 0]
[yp0]     [y0 -x0 0 1]
[xp1]     [x1  y1 1 0]     [M11]
[yp1]  =  [y1 -x1 0 1]  *  [M12]
[xp2]     [x2  y2 1 0]     [M13]
[yp2]     [y2 -x2 0 1]     [M23]
[xp3]     [x3  y3 1 0]
[yp3]     [y3 -x3 0 1]

求解上述矩阵方程后,由 M12 和 M11 计算 M21 和 M22。

于 2013-01-11T03:27:35.010 回答
6

我会尝试迭代最近点算法。

在这里,您可以找到具有缩放功能的实现。(SICP)

另一个有用的链接

于 2013-01-11T07:52:40.600 回答
2

为简单起见,假设您的输入x1,...,xn和输出y1,...,yn是复数。

  1. 您可以通过计算 来计算平均位移avg(y) - avg(x),一旦完成,您可以减去两者的平均值xy使它们以 0 为中心。

  2. 您现在想要找到旋转和比例。您可以将两者都表示为单个复数z,因此x*z应尽可能接近y. Butx*z是: 的(实)坐标的 R 线性函数,因此可以使用经典的线性代数来求解在最小二乘意义上尽可能接近的问题。(zx,zy)zzx*zy

综上所述,这为您提供了最小二乘意义上的最佳变换。

于 2012-07-27T13:47:57.397 回答
1

您的变换是仿射变换,可以写成 3*3 矩阵。所以你的问题基本上是计算从一组点到其他点的最小均方误差仿射变换。

这个问题在普通的计算几何文献中很简单地解决了。一本不错的经典书籍是:http ://www.robots.ox.ac.uk/~vgg/hzbook/hzbook1.html (没有广告,它只是大多数人的参考书)。您将找到有关 2D 和 3D 几何的所有信息。快速搜索“仿射变换 LMSE”等词也会为您提供信息和代码。

此外,您还可以使用其他类型的算法,稳健的算法,如 RANSAC。根据您的应用程序,朝那个方向走得更远可能会很有趣。

于 2012-07-27T13:38:18.483 回答
0

Matlab中更简单明了的代码可以给你转换。

And more complex C++ code(using VXL lib) with python and matlab wrapper included.

Or you can use some modificated ICP(iterative closest point) algorithm that is robust to noise.

于 2013-04-29T05:58:42.447 回答