1.問題描述：
給點兩個d維空間中的點集合
常用的轉換關系可能包括以下幾種：
（1）等距變換，需要求R和T，等距變換前后長度，面積，線之前的角度不變，常見于坐標轉換、物體移動，姿態變化等，自由度6（3+3）
（2）相似變換，需要求R和T及縮放系數s， 常用于非剛體移動，或不同目標匹配對齊，目標縮放等，自由度7（6+1）
（3）仿射變換（正交投影），平移變換T 和非均勻變換（A）的復合，A是可逆矩陣，不要求是正交矩陣，仿射變換對影像旋轉+平移+縮放+切變（shear，也叫傾斜變換），相比前兩種，變換影像的形狀發生了改變，但是原圖中平行線仍然保持平行，自由度12（9+3）
(4)射影變換（透視變換）
射影變換的不變數：重合關系、長度的交比，自由度：15（右下角v=1）
射影變換是對影像的旋轉+平移+縮放+切邊+射影，平行線變換后可能不再平行，而是交于一點，如平視圖轉換為鳥瞰圖，
關于二維空間（影像）的變化可參考：基礎知識：二維常見變換
2.1 等距變換求解
其求解數學運算式如下：
wi 表示每個點對之前的權重，
對于該問題共有6個自由度，至少需要兩組點對，對矩陣求最小二乘解，該方法我們暫不討論，
主流的做法是使用去中心化點集協方差矩陣SVD求解，步驟如下：
（1）構建上述問題模型：
（2）去中心化及協方差矩陣SVD分解求
計算中心值：
點集去中心化，去除轉移矩陣的影響：
計算協方差矩陣：
SVD分解，求旋轉矩陣
（3）計算轉移矩陣
具體推導程序可參見：
利用SVD求得兩個對應點集合的旋轉矩陣R和轉移矩陣t的數學推導
Least-Squares Rigid Motion Using SVD
python 代碼參考 計算兩個對應點集之間的旋轉矩陣R和轉移矩陣T
如下：
from numpy import *
from math import sqrt
def estimate_quilong_transform_3D(A, B):
assert len(A) == len(B)
N = A.shape[0];
mu_A = mean(A, axis=0)
mu_B = mean(B, axis=0)
AA = A - tile(mu_A, (N, 1))
BB = B - tile(mu_B, (N, 1))
H = transpose(AA) * BB
U, S, Vt = linalg.svd(H)
R = Vt.T * U.T
if linalg.det(R) < 0:
print "Reflection detected"
Vt[2, :] *= -1
R = Vt.T * U.T
t = -R * mu_A.T + mu_B.T
return R, t
2.2 相似變換求解
在等距變換的基礎上，點集P、Q的比例不一致，典型的應用為：3D模型融合，比如將給一個人頭3D模型加上眼鏡3d模型，兩模型尺寸不一致，模型匹配時我們在人頭和眼鏡上各選取4個關鍵點，求解相似變換矩陣，
點集合P、Q的比例s不影響旋轉矩陣R的求解，但是會影響T的結果，這是因為放大縮小是基于某一個參考坐標進行的，默認情況下對集合P的坐標乘以s進行縮放，等效于參考坐標為集合P的原點，那么集合P的中心位置也等效于乘了s，對應到T有平移分量，更合理的做法是將集合P的中心移到坐標原點，然后乘以縮放系數，
相似變換求解的難點通常在于點集P、Q的選取較為隨意，無嚴格定位和對齊，因此估計出來的變化矩陣誤差較大，
我們依然使用去中心化點集協方差矩陣SVD求解相似變換矩陣，
（1）估計出縮放系數s
首先計算點集合P、Q中所有點連成的線，線段個數為：n*(n-1)/2，然后計算所有線段的Ld階長度，然后對兩集合中所有線段長度相除得到n*(n-1)/2縮放系數，去除過短的線段和例外的縮放系數，剩下的取平均即可得到縮放系數的估計s，
（2）旋轉矩陣R與轉移矩陣T
旋轉矩陣R求解程序與2.1相同，
不同于2.1的是，在2.1中計算轉移矩陣時加入s
代碼如下：
def get_all_side_length(points):
all_dis=[]
for i in range(len(points)-1):
for j in range(i+1,len(points))：
all_dis.append(points[i]-points[j])
all_dis=np.array(all_dis)
return np.linalg.norm(all_dis,axis=1)
def get_scale(A,B):
dis_A=get_all_side_length(np.array(A))
dis_B=get_all_side_length(np.array(B))
scale=np.abs(dis_B/dis_A）
mask=np.abs(scale)>0.0001
scale_sort=np.sort(scale[mask].reshape(-1))
d_n=len(scale_sort)
s_mean=scale_sort[int(d_n/4):int(d_n*3/4)].mean() #only use medium data
return s_mean
def estimate_similarity_transform_3D(A, B):
## 求解R
assert len(A) == len(B)
N = A.shape[0];
mu_A = mean(A, axis=0)
mu_B = mean(B, axis=0)
AA = A - tile(mu_A, (N, 1))
BB = B - tile(mu_B, (N, 1))
H = transpose(AA) * BB
U, S, Vt = linalg.svd(H)
R = Vt.T * U.T
if linalg.det(R) < 0:
print "Reflection detected"
Vt[2, :] *= -1
R = Vt.T * U.T
s_mean=get_scale(A,B)
t = -s_mean * R * mu_A.T + mu_B.T
return R, t ,s_mean
2.3 仿射變換求解
仿射變換A的求解是典型的最小二乘求解方法，我們使用np.linalg.lstsq()方法，
自由度12，至少需要4對三維點（不共線）
代碼：
def estimate_affine_transform_3D(A, B):
'''
A:[n,3] B[n,3]
return:
P_Affine:(3,4) the third row is [0,0,0,1]
'''
A_homo=np.hstack((A,np.ones([A.shape[0],1]))) #nx4
P=np.linalg.lstsq(A_homo,Y)[0].T #Affine matrix 3 x 4
return P
def P2sRt(P)
''' decompositing camera matrix P
P:(3,4) Affine Camera Matrix
Returns:
s:
R:(3,3) rotation matrix
t:(3,) translation
'''
t=P[:,3]
R1=P[0:1,:3]
R2=P[1:2,:3]
s=(np.linalg.norm(R1) + np.linalg.norm(R2))/2.0
r1=R1/np.linalg.norm(R1)
r2=R2/np.linalg.norm(R2)
r3=np.cross(r1,r2) #叉乘，三維空間兩向量決定的平面法向量
R=np.concatenate((r1,r2,r3),0)
return s,R,t
2.4 射影變換求解
與仿射變換求解類似
轉載請註明出處，本文鏈接：https://www.uj5u.com/qukuanlian/290934.html
標籤：區塊鏈