Medical Image Registration - Warping / Interpolation

Gyuha Park·2021년 8월 26일
Computer VisionImage RegistrationLinear Algebra
0

Medical Image Analysis

목록 보기
17/21
post-thumbnail

1. Warping

1) Forward Warping

[abecdf001][xy1]=[x′y′1]\begin{bmatrix} a & b & e \\ c & d & f \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} x\\y\\1 \end{bmatrix}=\begin{bmatrix} x'\\y' \\ 1 \end{bmatrix}

[abecdf001][111]=[1.11.31]→[111]\begin{bmatrix} a & b & e \\ c & d & f \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} 1\\1\\1 \end{bmatrix}=\begin{bmatrix} 1.1\\1.3 \\ 1 \end{bmatrix}\rightarrow\begin{bmatrix} 1\\1 \\ 1 \end{bmatrix}

[abecdf001][211]=[2.31.11]→[211]\begin{bmatrix} a & b & e \\ c & d & f \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} 2\\1\\1 \end{bmatrix}=\begin{bmatrix} 2.3\\1.1 \\ 1 \end{bmatrix}\rightarrow\begin{bmatrix} 2\\1 \\ 1 \end{bmatrix}

[abecdf001][311]=[3.61.41]→[411]\begin{bmatrix} a & b & e \\ c & d & f \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} 3\\1\\1 \end{bmatrix}=\begin{bmatrix} 3.6\\1.4 \\ 1 \end{bmatrix}\rightarrow\begin{bmatrix} 4\\1 \\ 1 \end{bmatrix}

위 식과 같이 moving image의 (1,1,1), (2,1,1), (3,1,1)(1,1,1),\ (2,1,1),\ (3,1,1) 좌표에 transformation matrix를 곱해 새로운 좌표를 얻은 후 반올림하면 (1,1,1), (2,1,1), (4,1,1)(1,1,1),\ (2,1,1),\ (4,1,1)이 나온다. 구한 좌표에 해당하는 moving image의 intensity를 입력하면 된다.

이 때, (3,1,1)(3,1,1)이 비게 되는 문제가 발생하게 된다. 이것이 Forward warping의 단점이다.

2) Backward Warping

[abecdf001]−1[x′y′1]=[xy1]\begin{bmatrix} a & b & e \\ c & d & f \\ 0 & 0 & 1 \end{bmatrix}^{-1}\begin{bmatrix} x'\\y'\\1 \end{bmatrix}=\begin{bmatrix} x\\y \\ 1 \end{bmatrix}

[abecdf001]−1[111]=[2.11.31]→[211]\begin{bmatrix} a & b & e \\ c & d & f \\ 0 & 0 & 1 \end{bmatrix}^{-1}\begin{bmatrix} 1\\1\\1 \end{bmatrix}=\begin{bmatrix} 2.1\\1.3 \\ 1 \end{bmatrix}\rightarrow\begin{bmatrix} 2\\1 \\ 1 \end{bmatrix}

Backward warping은 반대로 transformed image의 좌표에 transformation matrix의 역행렬을 곱함으로 얻은 moving image 좌표의 intensity를 찾아가는 방법이다. 일반적으로 가장 많이 사용하는 방법이다.

2. Interpolation

1) Bilinear Interpolation

f(x,y1)=x2−xx2−x1Q11+x−x1x2−x1Q21f(x,y_1)=\cfrac{x_2-x}{x_2-x_1}Q_{11}+\cfrac{x-x_1}{x_2-x_1}Q_{21}

f(x,y2)=x2−xx2−x1Q12+x−x1x2−x1Q22f(x,y_2)=\cfrac{x_2-x}{x_2-x_1}Q_{12}+\cfrac{x-x_1}{x_2-x_1}Q_{22}

f(x,y)=y2−xy2−y1f(x,y1)+y−y1y2−y1f(x,y2)f(x,y)=\cfrac{y_2-x}{y_2-y_1}f(x,y_1)+\cfrac{y-y_1}{y_2-y_1}f(x,y_2)

=y2−xy2−y1(x2−xx2−x1Q11+x−x1x2−x1Q21)+=\cfrac{y_2-x}{y_2-y_1}\left( \cfrac{x_2-x}{x_2-x_1}Q_{11}+\cfrac{x-x_1}{x_2-x_1}Q_{21} \right)+

    y−y1y2−y1(x2−xx2−x1Q12+x−x1x2−x1Q22)\ \ \ \ \cfrac{y-y_1}{y_2-y_1}\left( \cfrac{x_2-x}{x_2-x_1}Q_{12}+\cfrac{x-x_1}{x_2-x_1}Q_{22} \right)

위 식에서 f(x,y)f(x,y)는 좌표 (x,y)(x,y)의 intensity 값을 나타낸다.

2) Bicubic Interpolation

주변 16개의 좌표가 p[i+1,j+1]=f(i,j)p[i+1,j+1]=f(i,j)로 주어지면 아래 그림과 같이 나타낼 수 있다.

Interpolated surface는 다음과 같다.
f(x,y)=∑i=03∑j=03aijxiyjf(x,y)=\sum\limits_{i=0}^3\sum\limits_{j=0}^3a_{ij}x^iy^j
이때, 16개의 aija_{ij}를 찾아야 한다. 이를 위해선 16개의 조건이 필요한데 16개의 주변 좌표를 이용해 얻을 수 있다.

  • 꼭짓점의 값 (4개)
    f(0,0)=a00f(0,0)=a_{00}
    f(1,0)=a00+a10+a20+a30f(1,0)=a_{00}+a_{10}+a_{20}+a_{30}
    f(0,1)=a00+a01+a02+a03f(0,1)=a_{00}+a_{01}+a_{02}+a_{03} f(1,1)=a00+a10+a20+a30+a01+a11+a21+a31+a02+a12+a22+a32+a03+a13+a23+a33f(1,1)=a_{00}+a_{10}+a_{20}+a_{30}+a_{01}+a_{11}+a_{21}+a_{31}+a_{02}+a_{12}+a_{22}+a_{32}+a_{03}+a_{13}+a_{23}+a_{33}

  • 꼭짓점의 미분계수 (8개)
    fx(0,0)=12(f(1,0)−f(−1,0))=a10f_x(0,0)=\cfrac{1}{2}(f(1,0)-f(-1,0))=a_{10}
    fx(1,0)=a10+2a20+3a30f_x(1,0)=a_{10}+2a_{20}+3a_{30}
    fx(0,1)=a10+a11+a12+a13f_x(0,1)=a_{10}+a_{11}+a_{12}+a_{13} fx(1,1)=a10+2a20+3a30+a11+2a21+3a31+a12+2a22+3a32+a13+2a23+3a33f_x(1,1)=a_{10}+2a_{20}+3a_{30}+a_{11}+2a_{21}+3a_{31}+a_{12}+2a_{22}+3a_{32}+a_{13}+2a_{23}+3a_{33}
    fy(0,0)=12(f(0,1)−f(0,−1))=a01f_y(0,0)=\cfrac{1}{2}(f(0,1)-f(0,-1))=a_{01}
    fy(1,0)=a01+a11+a21+a31f_y(1,0)=a_{01}+a_{11}+a_{21}+a_{31}
    fy(0,1)=a01+2a02+3a03f_y(0,1)=a_{01}+2a_{02}+3a_{03}
    fy(1,1)=a01+a11+a21+a31+2a02+2a12+2a22+2a32+3a03+3a13+3a23+3a33f_y(1,1)=a_{01}+a_{11}+a_{21}+a_{31}+2a_{02}+2a_{12}+2a_{22}+2a_{32}+3a_{03}+3a_{13}+3a_{23}+3a_{33}

  • 꼭짓점의 교차 미분계수 (4개)
    fxy(0,0)=14(f(1,1)−f(1,−1)−f(−1,1)+f(−1,−1))=a11f_{xy}(0,0)=\cfrac{1}{4}(f(1,1)-f(1,-1)-f(-1,1)+f(-1,-1))=a_{11}
    fxy(1,0)=a11+2a21+3a31f_{xy}(1,0)=a_{11}+2a_{21}+3a_{31}
    fxy(0,1)=a11+2a12+3a13f_{xy}(0,1)=a_{11}+2a_{12}+3a_{13}
    fxy(1,1)=a11+2a21+3a31+2a12+4a22+6a32+3a13+6a23+9a33f_{xy}(1,1)=a_{11}+2a_{21}+3a_{31}+2a_{12}+4a_{22}+6a_{32}+3a_{13}+6a_{23}+9a_{33}

위 식 들은 미지수 16개를 갖는 16개의 연립 방정식이므로 해를 구할 수 있다.
v=[a00,a10,a20,a30,a01,a11,a21,a31,a02,a12,a22,a32,a03,a13,a23,a33]T\text{v}=[a_{00},a_{10},a_{20},a_{30},a_{01},a_{11},a_{21},a_{31},a_{02},a_{12},a_{22},a_{32},a_{03},a_{13},a_{23},a_{33}]^T
f=[f(0,0),f(1,0),f(0,1),f(1,1),fx(0,0),fx(1,0),fx(0,1),fx(1,1),fy(0,0),fy(1,0),fy(0,1),fy(1,1),fxy(0,0),fxy(1,0),fxy(0,1),fxy(1,1)]\text{f}=[f(0,0),f(1,0),f(0,1),f(1,1),f_x(0,0),f_x(1,0),f_x(0,1),f_x(1,1),\\f_y(0,0),f_y(1,0),f_y(0,1),f_y(1,1),f_{xy}(0,0),f_{xy}(1,0),f_{xy}(0,1),f_{xy}(1,1)]

A=(1000000000000000111100000000000010001000100010001111111111111111010000000000000001230000000000000100010001000100012301230123012300001000000000000000111100000000000010002000300000001111222233330000010000000000000001230000000000000100020003000000012302460369)\text{A}=\begin{pmatrix}1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\1&1&1&1&0&0&0&0&0&0&0&0&0&0&0&0\\1&0&0&0&1&0&0&0&1&0&0&0&1&0&0&0\\1&1&1&1&1&1&1&1&1&1&1&1&1&1&1&1\\0&1&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\0&1&2&3&0&0&0&0&0&0&0&0&0&0&0&0\\0&1&0&0&0&1&0&0&0&1&0&0&0&1&0&0\\0&1&2&3&0&1&2&3&0&1&2&3&0&1&2&3\\0&0&0&0&1&0&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&1&1&1&1&0&0&0&0&0&0&0&0\\0&0&0&0&1&0&0&0&2&0&0&0&3&0&0&0\\0&0&0&0&1&1&1&1&2&2&2&2&3&3&3&3\\0&0&0&0&0&1&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&0&1&2&3&0&0&0&0&0&0&0&0\\0&0&0&0&0&1&0&0&0&2&0&0&0&3&0&0\\0&0&0&0&0&1&2&3&0&2&4&6&0&3&6&9\end{pmatrix}

A−1=(10000000000000000000100000000000−3300−2100000000002−2001100000000000000000010000000000000000000100000000000−3300−2−100000000002−2001100−30300000−20−1000000000−30300000−20−109−9−9963−6−36−63−34221−666−6−3−333−44−22−2−2−1−120−20000010100000000020−2000000010−666−6−4−242−33−33−2−1−2−14−4−4422−2−22−22−21111)\text{A}^{-1}=\begin{pmatrix}1&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&1&0&0&0&0&0&0&0&0&0&0&0\\-3&3&0&0&-2&1&0&0&0&0&0&0&0&0&0&0\\2&-2&0&0&1&1&0&0&0&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&1&0&0&0&0&0&0&0\\0&0&0&0&0&0&0&0&0&0&0&0&1&0&0&0\\0&0&0&0&0&0&0&0&-3&3&0&0&-2&-1&0&0\\0&0&0&0&0&0&0&0&2&-2&0&0&1&1&0&0\\-3&0&3&0&0&0&0&0&-2&0&-1&0&0&0&0&0\\0&0&0&0&-3&0&3&0&0&0&0&0&-2&0&-1&0\\9&-9&-9&9&6&3&-6&-3&6&-6&3&-3&4&2&2&1\\-6&6&6&-6&-3&-3&3&3&-4&4&-2&2&-2&-2&-1&-1\\2&0&-2&0&0&0&0&0&1&0&1&0&0&0&0&0\\0&0&0&0&2&0&-2&0&0&0&0&0&0&0&1&0\\-6&6&6&-6&-4&-2&4&2&-3&3&-3&3&-2&-1&-2&-1\\4&-4&-4&4&2&2&-2&-2&2&-2&2&-2&1&1&1&1\end{pmatrix}

f=Av\text{f}=\text{A}\text{v}
v=A−1f\text{v}=\text{A}^{-1}\text{f}

구한 결과는 다음과 같다.
a00=p11a_{00}=p_{11}
a01=−12p10+12p12a_{01}=-\cfrac{1}{2}p_{10}+\cfrac{1}{2}p_{12}
⋮\vdots
a33=a_{33}= 14p00−34p01+34p02−14p03−34p10+94p11−94p12+34p13+34p20−94p21+94p22−34p23−14p30+34p31−34p32+14p33\cfrac{1}{4}p_{00}-\cfrac{3}{4}p_{01}+\cfrac{3}{4}p_{02}-\cfrac{1}{4}p_{03}-\cfrac{3}{4}p_{10}+\cfrac{9}{4}p_{11}-\cfrac{9}{4}p_{12}+\cfrac{3}{4}p_{13}+\cfrac{3}{4}p_{20}-\cfrac{9}{4}p_{21}+\cfrac{9}{4}p_{22}-\cfrac{3}{4}p_{23}-\cfrac{1}{4}p_{30}+\cfrac{3}{4}p_{31}-\cfrac{3}{4}p_{32}+\cfrac{1}{4}p_{33}

profile
Gyuha Park
Medical Imaging & AI
이전 포스트

Medical Image Registration - Transformation matrix in 3D

다음 포스트

Medical Image Registration - Similarity measure

0개의 댓글