1

我需要找到从两个不同相机捕获的同一场景的两个图像中提取的两组独立特征之间的匹配。我正在使用 HumanEvaI 数据集，所以我有相机的校准矩阵（在本例中为 BW1 和 BW2）。

我不能使用简单相关、SIFT 或 SURF 之类的方法来解决问题，因为相机彼此相距很远并且还旋转。所以图像之间的差异很大，也有遮挡。

我一直专注于根据我已经拥有的校准信息而能够构建的地面实况点匹配来寻找捕获的图像之间的 Homography。一旦我有了这个单应性，我将使用完美匹配（匈牙利算法）来找到最佳对应关系。这里单应性的重要性在于我必须计算点之间的距离。

到目前为止一切似乎都很好，我的问题是我无法找到一个好的单应性。我已经尝试过 RANSAC 方法、具有 sampson 距离的黄金标准方法（这是一种非线性优化方法），并且主要来自 Richard Hartley 的《计算机视觉中的多视图几何》第二版一书中的所有内容。

到目前为止，我已经在 matlab 中实现了所有内容。

还有另一种方法可以做到这一点吗？我在正确的道路上吗？如果是这样，我可能做错了什么？

4

3 回答 3

2

你可以试试这两种方法：

  1. 一种用于非刚性配准的新点匹配算法（使用薄板样条） - 相对较慢。
  2. 点集注册：相干点漂移（我觉得更快更准确）。

就第二种方法而言，我觉得它在存在异常值的情况下给出了非常好的配准结果，它速度快并且能够恢复复杂的变换。但是第一种方法在注册领域也是众所周知的方法，您也可以尝试一下。

尝试理解算法的核心，然后继续阅读在线可用的代码。

  1. 此处的薄板样条- 下载 TPS-RPM 演示。
  2. 点集注册：相干点漂移在这里
于 2013-02-11T04:47:54.980 回答
2

你可能会觉得我的解决方案很有趣。它是 Coherent Point Drift 算法的纯 numpy实现。

这是一个例子：

from functools import partial
from scipy.io import loadmat
import matplotlib.pyplot as plt
import numpy as np
import time

class RigidRegistration(object):
    def __init__(self, X, Y, R=None, t=None, s=None, sigma2=None, maxIterations=100, tolerance=0.001, w=0):
    if X.shape[1] != Y.shape[1]:
        raise 'Both point clouds must have the same number of dimensions!'

    self.X             = X
    self.Y             = Y
    (self.N, self.D)   = self.X.shape
    (self.M, _)        = self.Y.shape
    self.R             = np.eye(self.D) if R is None else R
    self.t             = np.atleast_2d(np.zeros((1, self.D))) if t is None else t
    self.s             = 1 if s is None else s
    self.sigma2        = sigma2
    self.iteration     = 0
    self.maxIterations = maxIterations
    self.tolerance     = tolerance
    self.w             = w
    self.q             = 0
    self.err           = 0

def register(self, callback):
    self.initialize()

    while self.iteration < self.maxIterations and self.err > self.tolerance:
        self.iterate()
        callback(X=self.X, Y=self.Y)

    return self.Y, self.s, self.R, self.t

def iterate(self):
    self.EStep()
    self.MStep()
    self.iteration = self.iteration + 1

def MStep(self):
    self.updateTransform()
    self.transformPointCloud()
    self.updateVariance()

def updateTransform(self):
    muX = np.divide(np.sum(np.dot(self.P, self.X), axis=0), self.Np)
    muY = np.divide(np.sum(np.dot(np.transpose(self.P), self.Y), axis=0), self.Np)

    self.XX = self.X - np.tile(muX, (self.N, 1))
    YY      = self.Y - np.tile(muY, (self.M, 1))

    self.A = np.dot(np.transpose(self.XX), np.transpose(self.P))
    self.A = np.dot(self.A, YY)

    U, _, V = np.linalg.svd(self.A, full_matrices=True)
    C = np.ones((self.D, ))
    C[self.D-1] = np.linalg.det(np.dot(U, V))

    self.R = np.dot(np.dot(U, np.diag(C)), V)

    self.YPY = np.dot(np.transpose(self.P1), np.sum(np.multiply(YY, YY), axis=1))

    self.s = np.trace(np.dot(np.transpose(self.A), self.R)) / self.YPY

    self.t = np.transpose(muX) - self.s * np.dot(self.R, np.transpose(muY))

def transformPointCloud(self, Y=None):
    if not Y:
        self.Y = self.s * np.dot(self.Y, np.transpose(self.R)) + np.tile(np.transpose(self.t), (self.M, 1))
        return
    else:
        return self.s * np.dot(Y, np.transpose(self.R)) + np.tile(np.transpose(self.t), (self.M, 1))

def updateVariance(self):
    qprev = self.q

    trAR     = np.trace(np.dot(self.A, np.transpose(self.R)))
    xPx      = np.dot(np.transpose(self.Pt1), np.sum(np.multiply(self.XX, self.XX), axis =1))
    self.q   = (xPx - 2 * self.s * trAR + self.s * self.s * self.YPY) / (2 * self.sigma2) + self.D * self.Np/2 * np.log(self.sigma2)
    self.err = np.abs(self.q - qprev)

    self.sigma2 = (xPx - self.s * trAR) / (self.Np * self.D)

    if self.sigma2 <= 0:
        self.sigma2 = self.tolerance / 10

def initialize(self):
    self.Y = self.s * np.dot(self.Y, np.transpose(self.R)) + np.repeat(self.t, self.M, axis=0)
    if not self.sigma2:
        XX = np.reshape(self.X, (1, self.N, self.D))
        YY = np.reshape(self.Y, (self.M, 1, self.D))
        XX = np.tile(XX, (self.M, 1, 1))
        YY = np.tile(YY, (1, self.N, 1))
        diff = XX - YY
        err  = np.multiply(diff, diff)
        self.sigma2 = np.sum(err) / (self.D * self.M * self.N)

    self.err  = self.tolerance + 1
    self.q    = -self.err - self.N * self.D/2 * np.log(self.sigma2)

def EStep(self):
    P = np.zeros((self.M, self.N))

    for i in range(0, self.M):
        diff     = self.X - np.tile(self.Y[i, :], (self.N, 1))
        diff    = np.multiply(diff, diff)
        P[i, :] = P[i, :] + np.sum(diff, axis=1)

    c = (2 * np.pi * self.sigma2) ** (self.D / 2)
    c = c * self.w / (1 - self.w)
    c = c * self.M / self.N

    P = np.exp(-P / (2 * self.sigma2))
    den = np.sum(P, axis=0)
    den = np.tile(den, (self.M, 1))
    den[den==0] = np.finfo(float).eps

    self.P   = np.divide(P, den)
    self.Pt1 = np.sum(self.P, axis=0)
    self.P1  = np.sum(self.P, axis=1)
    self.Np = np.sum(self.P1)

def visualize(X, Y, ax):
    plt.cla()
    ax.scatter(X[:,0] ,  X[:,1], color='red')
    ax.scatter(Y[:,0] ,  Y[:,1], color='blue')
    plt.draw()
    plt.pause(0.001)

def main():
    fish = loadmat('./data/fish.mat')
    X = fish['X'] # number-of-points x number-of-dimensions array of fixed points
    Y = fish['Y'] # number-of-points x number-of-dimensions array of moving points

    fig = plt.figure()
    fig.add_axes([0, 0, 1, 1])
    callback = partial(visualize, ax=fig.axes[0])

    reg = RigidRegistration(X, Y)
    reg.register(callback)
    plt.show()

if __name__ == '__main__':
main()
于 2016-10-31T21:53:00.453 回答
0

此处描述了我认为您可能会发现有用的另一种方法。
这种方法试图将局部模型拟合到一组点。当存在多个不同的局部模型时，其全局优化方法使其性能优于 RANSAC。
我也相信他们有可用的代码。

于 2013-02-11T07:38:09.487 回答