OpenCV中轮廓之间的仿射变换

Question

我有一系列的海底图像，这些图像是从胶片上扫描得来的，需要进行注册处理。

from pylab import *
import cv2
import urllib

urllib.urlretrieve('http://geoport.whoi.edu/images/frame014.png','frame014.png');
urllib.urlretrieve('http://geoport.whoi.edu/images/frame015.png','frame015.png');

gray1=cv2.imread('frame014.png',0)
gray2=cv2.imread('frame015.png',0)
figure(figsize=(14,6))
subplot(121);imshow(gray1,cmap=cm.gray);
subplot(122);imshow(gray2,cmap=cm.gray);

我希望利用每张图片左侧的黑色区域进行配准，因为该区域位于相机内部且应在时间上固定。因此，我只需要计算黑色区域之间的仿射变换。

我通过阈值处理和寻找最大轮廓来确定这些区域:

def find_biggest_contour(gray,threshold=40):
    # threshold a grayscale image 
    ret,thresh = cv2.threshold(gray,threshold,255,1)
    # find the contours
    contours,h = cv2.findContours(thresh,mode=cv2.RETR_LIST,method=cv2.CHAIN_APPROX_NONE)
    # measure the perimeter
    perim = [cv2.arcLength(cnt,True) for cnt in contours]
    # find contour with largest perimeter
    i=perim.index(max(perim))
    return contours[i]

c1=find_biggest_contour(gray1)
c2=find_biggest_contour(gray2)

x1=c1[:,0,0];y1=c1[:,0,1]
x2=c2[:,0,0];y2=c2[:,0,1]

figure(figsize=(8,8))
imshow(gray1,cmap=cm.gray, alpha=0.5);plot(x1,y1,'b-')
imshow(gray2,cmap=cm.gray, alpha=0.5);plot(x2,y2,'g-')
axis([0,1500,1000,0]);

蓝色线条是第一帧中最长的轮廓，绿色线条则是第二帧中最长的轮廓。

如何最好地确定蓝色和绿色轮廓之间的旋转和偏移？

我只想在台阶周围的某个区域使用轮廓的右侧，类似于箭头之间的区域。

当然，如果有更好的方法来对齐这些图像，我很乐意听取建议。我已经尝试了标准的特征匹配方法，但效果不太好。

- Rich Signell

---

你只需要绿色和蓝色轮廓之间的向量，还是需要向量、旋转角度和缩放比例？ - cyriel

仅偏移和旋转（无缩放） - Rich Signell

Answer 1

根据Shambool建议的方法，我得出了以下结果。我使用了Ramer-Douglas-Peucker算法来简化感兴趣区域内的轮廓，并确定了两个转折点。我原本想使用这两个转折点来获取我的三个未知数（x偏移、y偏移和旋转角度），但第二个转折点太靠右了，因为RDP在这个区域简化了较平滑的曲线。所以我改用了指向第一个转折点的线段的角度。将图像1和图像2之间的这个角度差分给我旋转角度。我对这个解决方案还不完全满意。它对这两张图片的效果已经足够好了，但我不确定它是否能在整个图像序列上表现良好。我们拭目以待。

最好将轮廓拟合到黑色边框的已知形状上。

# select region of interest from largest contour 
ind1=where((x1>190.) & (y1>200.) & (y1<900.))[0]
ind2=where((x2>190.) & (y2>200.) & (y2<900.))[0]
figure(figsize=(10,10))
imshow(gray1,cmap=cm.gray, alpha=0.5);plot(x1[ind1],y1[ind1],'b-')
imshow(gray2,cmap=cm.gray, alpha=0.5);plot(x2[ind2],y2[ind2],'g-')
axis([0,1500,1000,0])

def angle(x1,y1):
    #  Returns angle of each segment along an (x,y) track
    return array([math.atan2(y,x) for (y,x) in zip(diff(y1),diff(x1))])

def simplify(x,y, tolerance=40, min_angle = 60.*pi/180.): 
    """
    Use the Ramer-Douglas-Peucker algorithm to simplify the path
    http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
    Python implementation: https://github.com/sebleier/RDP/
    """
    from RDP import rdp   
    points=vstack((x,y)).T
    simplified = array(rdp(points.tolist(), tolerance))
    sx, sy = simplified.T

    theta=abs(diff(angle(sx,sy)))
    # Select the index of the points with the greatest theta
    # Large theta is associated with greatest change in direction.
    idx = where(theta>min_angle)[0]+1
    return sx,sy,idx

sx1,sy1,i1 = simplify(x1[ind1],y1[ind1])
sx2,sy2,i2 = simplify(x2[ind2],y2[ind2])
fig = plt.figure(figsize=(10,6))
ax =fig.add_subplot(111)

ax.plot(x1, y1, 'b-', x2, y2, 'g-',label='original path')
ax.plot(sx1, sy1, 'ko-', sx2, sy2, 'ko-',lw=2, label='simplified path')
ax.plot(sx1[i1], sy1[i1], 'ro', sx2[i2], sy2[i2], 'ro', 
    markersize = 10, label='turning points')
ax.invert_yaxis()
plt.legend(loc='best')

# determine x,y offset between 1st turning points, and 
# angle from difference in slopes of line segments approaching 1st turning point
xoff = sx2[i2[0]] - sx1[i1[0]]
yoff = sy2[i2[0]] - sy1[i1[0]]
iseg1 = [i1[0]-1, i1[0]]
iseg2 = [i2[0]-1, i2[0]]
ang1 = angle(sx1[iseg1], sy1[iseg1])
ang2 = angle(sx2[iseg2], sy2[iseg2])
ang = -(ang2[0] - ang1[0])
print xoff, yoff, ang*180.*pi

-28 14 5.07775871644

# 2x3 affine matrix M
M=array([cos(ang),sin(ang),xoff,-sin(ang),cos(ang),yoff]).reshape(2,3)
print M

[[  9.99959685e-01   8.97932821e-03  -2.80000000e+01]
 [ -8.97932821e-03   9.99959685e-01   1.40000000e+01]]

# warp 2nd image into coordinate frame of 1st
Minv = cv2.invertAffineTransform(M)
gray2b = cv2.warpAffine(gray2,Minv,shape(gray2.T))

figure(figsize=(10,10))
imshow(gray1,cmap=cm.gray, alpha=0.5);plot(x1[ind1],y1[ind1],'b-')
imshow(gray2b,cmap=cm.gray, alpha=0.5);
axis([0,1500,1000,0]);
title('image1 and transformed image2 overlain with 50% transparency');

- Rich Signell

---

这里的Stackoverflow礼仪是什么？我应该选择Shambool的建议，因为他启发了我最终使用的解决方案吗？还是我应该选择我的解决方案，因为它包含代码，并且这是我实际使用的内容？ - Rich Signell

Answer 2

好问题。

一种方法是将轮廓表示为2D点云，然后进行配准。 更简单明了的Matlab代码可以给你仿射变换。

还有更复杂的C++代码(使用VXL库)，包括Python和Matlab包装器。 或者您可以使用一些修改版的迭代最近点(ICP)算法，该算法对噪声具有鲁棒性，并且可以处理仿射变换。

另一种方法是使用使用像素值的某种配准方法。 Matlab代码(我认为它使用某种最小化+交叉相关度量) 也许在医学影像中使用了一些光流配准(或其他类型的配准)。

您还可以使用SIFT(SURF)等点特征。

您可以在 FIJI(ImageJ) 中快速尝试， 也可以查看这个链接。

1. 打开两个图像
2. 插件->特征提取->sift(或其他)
3. 将期望转换设置为仿射
4. 查看在ImageJ日志中估计的变换模型[3,3]同态矩阵。 如果好用，那么您可以使用OpenCV在Python中实现它，或者使用带有ImageJ的Jython。

如果您贴出原始图像并描述所有条件(似乎图像在帧之间发生了变化)，效果会更好。

- mrgloom

---

我已经发布了原始图像。这些图像确实会随着海底条件的变化而改变。这就是为什么我想使用黑色区域来注册图像，因为它不应该发生变化。 - Rich Signell

Answer 3

您可以用它们各自的椭圆来表示这些轮廓。这些椭圆以轮廓的质心为中心，并朝向主密度轴。您可以比较质心和方向角。

1）填充轮廓=>使用thickness=CV_FILLED绘制轮廓

2）查找矩=>cvMoments()

3）然后使用它们。

质心：{ x，y } = {M10 / M00，M01 / M00 }

方向（theta）： 编辑：我为您的情况定制了遗留（enteringblobdetection.cpp）的示例代码。

            /* Image moments */
            double      M00,X,Y,XX,YY,XY;
            CvMoments   m;
            CvRect      r = ((CvContour*)cnt)->rect;
            CvMat       mat;
            cvMoments( cvGetSubRect(pImgFG,&mat,r), &m, 0 );
            M00 = cvGetSpatialMoment( &m, 0, 0 );
            X = cvGetSpatialMoment( &m, 1, 0 )/M00;
            Y = cvGetSpatialMoment( &m, 0, 1 )/M00;
            XX = (cvGetSpatialMoment( &m, 2, 0 )/M00) - X*X;
            YY = (cvGetSpatialMoment( &m, 0, 2 )/M00) - Y*Y;  
            XY = (cvGetSpatialMoment( &m, 1, 1 )/M00) - X*Y; 

            /* Contour description */
            CvPoint myCentroid(r.x+(float)X,r.y+(float)Y);
            double myTheta =  atan( 2*XY/(XX-YY) );

另外，请查看使用OpenCV 2.0示例的this。

- LovaBill

---

这是一个很酷的想法，但我认为在这种情况下不会起作用，因为轮廓的顶部和左侧的那些东西会影响到矩形的计算，而这些东西我不想使用。 - Rich Signell

你可以裁剪对象。 - LovaBill

我也考虑过这个方法，但是在我看来，为了使这种方法给出正确的答案，需要使用旋转和平移的边界框对对象进行裁剪，这就需要已经知道答案了。不是吗？ - Rich Signell

Answer 4

如果您不想找到两幅图像之间的单应性，而是要找到仿射变换，则需要三个未知数：旋转角度（R）以及在x和y方向上的位移（X，Y）。因此，需要至少两个点（每个点有两个已知值）才能找到未知数。两个图像之间应该匹配两个点或两条线，每个点或线都有两个已知值：截距和斜率。如果采用点匹配方法，则两个点之间的距离越远，找到的变换对噪声的鲁棒性就越强（如果您记得误差传播规则，这很简单）。

在两点匹配法中：

1. 在第一幅图像I1中找到两个点（A和B），并在第二幅图像I2中找到它们对应的点（A'和B'）。 2. 找到A和B之间的中点C，以及A'和B'之间的中点C'。 3. 差值C和C'（C-C'）给出了图像之间的平移（X和Y）。 4. 使用C-A和C'-A'的点积可以找到旋转角度（R）。

为了检测稳健点，我会找到你已经找到的柜台侧面上二阶导数（海森矩阵）绝对值最高的点，并尝试匹配它们。由于你提到这是视频镜头，你可以轻松地假设每两帧之间的变换很小，以拒绝离群值。

- fireant

---

这似乎是最有希望的。我正在努力获取那两个点。 - Rich Signell

我今天太忙了，无法帮助编码。一旦您有候选匹配点，由于图像之间的简单转换，有几种简单的方法可以消除不匹配/异常值，因此一种方法是，而不是努力寻找最佳匹配，您可以找到一组候选匹配，然后将它们细化，直到获得最佳的两个匹配。明天可能我能够协助编码，以防您在那时还没有找到一个稳健的解决方案。 - fireant

我大部分都按照您的建议并发布了代码。但是我对于确定我的三个未知数的方式仍不完全满意。如果您有更好的想法，请告诉我。 - Rich Signell

你能上传视频到YouTube或其他地方，这样我们就可以看一下吗？ - fireant

你想要原始图像的视频，还是使用上面的代码来注册它们？ - Rich Signell

正如我所料，这种方法在图像序列上效果不佳。最大的轮廓并不总是被捕捉到，而且拐点也不总是准确的。我想我要尝试另一种完全不同的方法，试图匹配黑色边框右侧的曲线线条。 - Rich Signell