Medical Image Registration - Non-rigid registration via ICP

1. Iterative Closest Point

위 그림에서 moving image에서 빨간색 점 선, fixed image에서 파란색 점 선의 feature를 뽑았다고 가정하자.

첫 번째로 moving image의 feature의 중심 좌표를 구한 후 fixed image의 feature의 중심 좌표로 translation한다.

두 번째로 translation 후 moving image의 feature와 fixed image의 feature의 서로 가장 가까운 closest point를 찾아 matching되는 점들 사이의 rotation matrix를 pseudo inverse를 이용해 구한다.

이와 같은 과정을 closest point 간의 거리가 일정 수치 이내로 들어올 때 까지 계속 반복하는 것이 ICP이다.

2. Non-rigid ICP Registration

Non-rigid registration은 rigid registration과 같이 모든 moving image의 points를 fixed image의 point로 이동할 하나의 transformation을 구하는 것이 아니라, 각 point $v$마다 transformation $X$를 구하는 것이 목적이다.

Non-rigid ICP registration은 아래 세가지 식 들을 최소화 함으로 최적의 $X$를 반복해서 구한다.

• Distance term

  $\bar{E}_d(X):=\sum\limits_{v_i\in V}w_i\text{dist}^2(u_i,X_iv_i)$

  Fixed image의 point $T$와 transformed moving image의 point $X_iv_i$의 거리를 나타낸다. $w_i$는 각각의 point마다 주는 가중치이다.

  위 식을 변형하면 다음과 같다.

  $W:=\text{diag}(w_1,\ldots,w_n)$

  $I_n=n\times n \ \ \text{identity matrix}$

  $\otimes =\text{Kronecker product}$

  $\sum\limits_{v_i\in V}w_i\text{dist}^2(u_i,X_iv_i)=\left|\left|(W\otimes I_3)\left(\begin{bmatrix}X_1& & \\ &\ddots& \\ & &X_n\end{bmatrix}\begin{bmatrix}v_1\\\vdots\\v_n\end{bmatrix}-\begin{bmatrix}u_1\\\vdots\\u_n\end{bmatrix}\right)\right|\right|^2$

  $D:=\begin{bmatrix}v_1^T& & & \\ &v_2^T\\ & &\ddots& \\ & & &v_n^T\end{bmatrix}$

  $U:=[u_1,\ldots,u_n]^T$

  $\bar{E}_d(X)=||W(DX-U)||_F^2$

• Stiffness term

  $\bar{E}_s(X):=\sum\limits_{\{i,j\}\in \epsilon}||(X_i-X_j)G||_F^2$

  $G:=\text{diag}(1,1,1,\gamma)$

  $M=\text{node-arc incidence matrix}$

  $E_s(X)=||(M\otimes G)X||_F^2$

  서로 인접한 points의 transformation $X$의 차이가 최소가 되도록 하여 stiffness 효과를 주는 것이다. $G$는 affine transformation과 translation term 사이의 비율을 조정하는 역할을 한다.

• Landmark term

  $\bar{E}_l(X):=\sum\limits_{(v_i,l)\in\mathcal{L}}||X_iv_i-l||^2$

  $E_l(X)=||D_LX-U_L||_F^2$

  Landmark와 그에 해당되는 $X_iv_i$의 거리 차를 나타낸다.

• Complete cost function

  전체 식 들을 조합하면 다음과 같다.

  $\bar{E}(X):=\bar{E}_d(X)+\alpha E_s(X)+\beta E_l(X)$

  $\alpha$값은 iteration 마다 점점 줄여나간다.

  그리고 아래와 같이 $AX-B$ 형태로 변경이 가능하다.

  $\bar{E}(X)=\left|\left|\begin{bmatrix}\alpha M\otimes G\\WD\\\beta D_L\end{bmatrix}X-\begin{bmatrix}0\\WU\\U_L\end{bmatrix}\right|\right|_F^2=||AX-B||_F^2$

  즉, pseudo inverse로 최적의 $X$를 구할 수 있게 된다.

  $X=(A^TA)^{-1}A^TB$