使用Python和交叉相关进行图像配准
我有两张显示完全相同内容的图片:2D高斯形状的斑点。我把这两张16位的png文件叫做“left.png”和“right.png”。不过因为它们是通过稍微不同的光学设置得到的,所以对应的斑点(实际上是一样的)出现在稍微不同的位置。这意味着右边的图片看起来有点拉伸、扭曲,或者说是以一种非线性的方式变形。因此,我想要找到从左边到右边的变换关系。
具体来说,我想为左边每个像素的x和y坐标找到一个函数,这个函数能给我一个位移向量的组成部分,指向右边对应的像素。
在之前的尝试中,我试图找到对应斑点的位置,以获得相对距离deltaX和deltaY。然后我把这些距离拟合到泰勒展开的二阶,得到了每个左边像素(x,y)指向右边对应像素(x',y')的位移向量的x和y分量。
为了得到一个更通用的结果,我想使用归一化互相关。为此,我将左边每个像素的值与右边对应像素的值相乘,并对这些乘积求和。我想要的变换应该连接那些能使这个和最大的像素。所以当和最大时,我就知道我乘的是对应的像素。
我真的尝试了很多,但没能成功。我的问题是,有没有人有想法或者做过类似的事情。
import numpy as np
import Image
left = np.array(Image.open('left.png'))
right = np.array(Image.open('right.png'))
# for normalization (http://en.wikipedia.org/wiki/Cross-correlation#Normalized_cross-correlation)
left = (left - left.mean()) / left.std()
right = (right - right.mean()) / right.std()
如果我可以让这个问题更清楚,请告诉我。我还需要了解如何使用latex发布问题。
非常感谢大家的意见。
[left.png] https://i.stack.imgur.com/oSTER.png [right.png] https://i.stack.imgur.com/Njahj.png
我担心,在大多数情况下,16位图像在我使用的系统上看起来只是黑色 :( 但当然里面是有数据的。
更新 1
我试图澄清我的问题。我在寻找一个向量场,这个向量场的位移向量指向从每个像素在left.png到对应的像素在right.png。我的问题是,我不确定我有哪些限制。
其中向量r(x和y分量)指向left.png中的一个像素,而向量r-prime(x-prime和y-prime分量)指向right.png中的对应像素。对于每个r都有一个位移向量。
我之前做的是,手动找到向量场d的分量,并将它们拟合到二次多项式:
所以我拟合了:
这样说有道理吗?是否可以通过互相关得到所有的delta-x(x,y)和delta-y(x,y)?如果通过位移向量将对应的像素连接在一起,互相关应该是最大的,对吗?
更新 2
我想到的算法如下:
- 变形right.png
- 获取互相关的值
- 进一步变形right.png
- 获取互相关的值并与之前的值进行比较
- 如果更大,说明变形效果好;如果不大,重新变形并尝试其他方法
- 在最大化互相关值后,知道变形的情况 :)
关于变形:能否先沿x和y方向移动以最大化互相关,然后在第二步中进行x和y方向的拉伸或压缩,第三步再进行二次的x和y方向变形,并重复这个过程?我在使用整数坐标时真的遇到问题。你认为我需要对图片进行插值,以获得连续的分布吗?我还需要再考虑一下 :( 感谢大家的参与 :)
3 个回答
你可以看看 bunwarpj,这个工具已经能做到你想要的功能。虽然它不是用Python写的,但我在类似的情况下用过它。你可以导出一个普通文本格式的样条变换,如果你想的话,可以使用这个功能。
我觉得交叉相关在这里可能帮不上忙,因为它只能给你整个图像的一个最佳位移。这里有三个我会考虑的替代方案:
对点的子集进行交叉相关。比如说,取右上角的三个点,通过交叉相关找到最佳的x-y位移。这可以给你左上角的粗略变换。对尽可能多的点群重复这个过程,以获得合理的变换图。然后用你的泰勒展开式来拟合,可能会得到比较接近的结果。不过,要让交叉相关有效,点之间的位移差必须小于点的大小,否则你就无法让一个点群中的所有点同时重叠。在这种情况下,第二个选项可能更合适。
如果位移相对较小(我认为这是第一个选项的条件),那么我们可以假设左图中的某个点在右图中最近的点就是对应的点。因此,对于左图中的每个点,我们找到右图中最近的点,并将其作为该位置的位移。通过40多个分布良好的位移向量,我们可以通过拟合你的泰勒展开式来获得实际位移的合理近似。
这可能是最慢的方法,但如果你有较大的位移(而第二个选项因此不适用),这可能是最稳健的:使用类似进化算法的方法来寻找位移。应用一个随机变换,计算剩余误差(你可能需要将其定义为原始图像和变换图像中点之间最小距离的总和),并根据这些结果改进你的变换。如果你的位移相当大,你可能需要进行非常广泛的搜索,因为你可能会在你的搜索空间中遇到很多局部最小值。
我会尝试第二个选项,因为你的位移似乎可能足够小,可以轻松地将左图中的点与右图中的点对应起来。
更新
我猜测你的光学设备引入了非线性失真,并且有两个独立的光束路径(每个路径使用不同的滤镜?)会使这两幅图像之间的关系更加非线性。PiQuer建议的仿射变换可能是一个合理的方法,但可能永远无法完全覆盖实际的失真。
我认为你用低阶泰勒多项式进行拟合的方法是可行的。这在我所有类似条件的应用中都有效。最高的阶数可能应该是xy^2和x^2y;超过这个阶数你可能就察觉不到了。
另外,你也许可以先对每幅图像的失真进行校准,然后再进行实验。这样你就不依赖于点的分布,而可以使用高分辨率的参考图像来获得最佳的变换描述。
上面提到的第二个选项仍然是我建议的让两幅图像重叠的方法。这可以完全自动化,我不太明白你所说的更一般的结果是什么意思。
更新2
你提到在两幅图像中匹配点时遇到了困难。如果是这样,我认为你的迭代交叉相关方法可能也不太稳健。你的点非常小,因此只有当两幅图像之间的差异很小时,它们才能重叠。
原则上,你提出的解决方案没有问题,但它是否有效很大程度上取决于你的变形大小和优化算法的稳健性。如果你一开始重叠很少,那么可能很难找到一个好的优化起点。然而,如果一开始就有足够的重叠,那么你应该能够先找到每个点的变形,但在评论中你表示这并不奏效。
也许你可以尝试一种混合方案:找到点群的交叉相关,以获取优化的起点,然后使用你在更新中描述的过程来调整变形。因此:
- 对于一个NxN像素的片段,找到左图和右图之间的位移
- 对大约16个这样的片段重复这个过程
- 使用这16个点计算变形的近似值
- 将其作为你优化方法的起点
OpenCV(还有它的Python绑定)有一个叫做StarDetector的类,它实现了这个算法。
另外,你也可以看看OpenCV的SIFT类,SIFT代表的是尺度不变特征变换。
更新
关于你的评论,我明白“正确”的变换会最大化图像之间的互相关,但我不太明白你是如何选择要最大化的变换集合的。也许如果你知道三个匹配点的坐标(可以通过一些启发式方法或手动选择),并且如果你期望有仿射变换,你可以使用类似cv2.getAffineTransform的东西来获得一个好的初始变换,这样可以帮助你进行最大化过程。然后你可以使用一些小的附加变换来形成一个可以进行最大化的集合。但我觉得这个方法有点像是在重新发明SIFT已经能处理的事情。
要实际变换你的测试图像,你可以使用cv2.warpAffine,它也可以处理边界值(例如用0填充)。要计算互相关,你可以使用scipy.signal.correlate2d。
更新
你的最新更新确实让我澄清了一些观点。但我认为位移的向量场并不是最自然的寻找方式,这也是误解的来源。我更倾向于寻找一个全局变换T,应用于左侧图像的任意点(x,y)后,得到右侧的(x',y')=T(x,y),而且T对每个像素的形式是相同的。例如,这可以是位移、旋转、缩放,可能还有一些透视变换的组合。我不能说找到这样的变换是否现实,这取决于你的设置,但如果场景在两侧是物理上相同的,我认为期待某种仿射变换是合理的。这就是我建议使用cv2.getAffineTransform的原因。当然,从这样的T计算你的位移向量场是很简单的,因为这只是T(x,y)-(x,y)。
这样做的一个大优点是,你的变换自由度非常少,而我认为在位移向量场中有2N个自由度,其中N是亮点的数量。
如果确实是仿射变换,我建议使用类似这样的算法:
- 在左侧识别三个明亮且孤立的点
- 为这三个点中的每一个定义一个边界框,以便你可以希望在右侧图像中找到对应的点
- 找到对应点的坐标,例如使用某种相关性方法,像在cv2.matchTemplate中实现的,或者仅仅在边界框内找到最亮的点。
- 一旦你有了三个匹配的坐标对,使用cv2.getAffineTransform计算将一组坐标转换为另一组的仿射变换。
- 将这个仿射变换应用于左侧图像,作为检查,如果你找到的变换是正确的,你可以计算整体的归一化互相关是否超过某个阈值,或者在你将一幅图像相对于另一幅图像移动时是否显著下降。
- 如果你需要,还可以从你的变换T简单地计算位移向量场。
更新
看起来cv2.getAffineTransform期望输入数据类型为'float32',这有点麻烦。假设源坐标是(sxi,syi)
,目标坐标是(dxi,dyi)
,其中i=0,1,2
,那么你需要的是
src = np.array( ((sx0,sy0),(sx1,sy1),(sx2,sy2)), dtype='float32' )
dst = np.array( ((dx0,dy0),(dx1,dy1),(dx2,dy2)), dtype='float32' )
result = cv2.getAffineTransform(src,dst)