如何從兩個不同相機拍攝的同一場景的兩幅圖像中提取兩個獨立特徵集之間的匹配？

Question

我需要找到從兩個不同相機捕獲的同一場景的兩幅圖像中提取的兩組獨立特徵之間的匹配。 我使用的是HumanEvaI數據集，因此我有相機的校準矩陣（在這種情況下爲BW1和BW2）。

我不能使用像簡單相關，SIFT或SURF這樣的方法來解決這個問題，因爲這些相機距離彼此很遠並且也是旋轉的。所以圖像之間的差異很大，並且也存在遮擋。

我一直專注於根據基於我已經擁有的校準信息已經能夠構建的地面實況點匹配的捕獲圖像找到一個Homography。 一旦我有了這個單應性，我將使用完美的匹配（匈牙利算法）來找到最好的對應關係。這裏單應性的重要性在於我必須計算點之間的距離。

到目前爲止，一切似乎都很好，我的問題是我一直無法找到一個好的單應性。我已經嘗試了RANSAC方法，具有桑普森距離的金標準方法（這是一種非線性優化方法），並且主要來自理查德哈特利編着的「計算機視覺中的多視圖幾何」第二版。

到目前爲止，我已經在matlab中實現了一切。

是否有另一種方法可以做到這一點？我在正確的道路上？如果是的話，我可能會做錯什麼？

Answer 1

我想你可能會發現有用的另一種方法是描述here。
該方法嘗試將局部模型擬合到一組點。當存在多個不同的局部模型時，其全局優化方法允許其優於RANSAC。
我也相信他們有可用的代碼。

Answer 2

可以嘗試這兩種方法：

1. 用於非剛性配準一個新的點匹配算法（使用薄板樣條） - 相對較慢。
2. 點集註冊：相干點漂移（我覺得更快，更準確）。

就第二種方法而言，我覺得它在異常值的情況下給出了非常好的配準結果，它很快並且能夠恢復複雜的轉換。但第一種方法在註冊領域也是一個衆所周知的方法，您也可以嘗試。

試着理解算法的核心，然後轉到在線提供的代碼。

1. 薄板樣條here - 下載TPS-RPM演示。
2. 點集註冊：相干點漂移here

Answer 3

您可能會發現我的解決方案很有趣。它是相干點漂移算法的純粹的numpy implementation。

下面是一個例子：

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()