Motion correction in optical coherence tomography volumes on a per A-scan basis using orthogonal scan patterns
1. Introduction
신체 조직을 따로 떼어내지 않고 직접적인 접촉 없이 마이크로미터 단위로 내부 구조를 볼 수 있다는 점에서 OCT (Optical Coherence Tomography)는 중요한 도구로 자리잡아 가고 있다. Spectical, Swept source OCT의 개발로 속도가 빨라진 OCT는 이제 3D 볼륨을 얻는데 몇초 정도가 소요된다. 하지만 이러한 속도 향상에도 불구하고 머리나 눈의 움직임으로 인한 모션 아티팩트 문제를 뛰어넘기는 아직 힘들어보인다. 이러한 모션 아티팩트의 원인으로, axial 방향으로 움직임을 일으키는 심장 박동이나 호흡이 있다. 또한 saccades라 불리는 눈의 drift나 시점 이동 등 trasversal 방향의 영향을 주는 요인도 존재한다. 이러한 움직임이 없을 경우에만 완벽한 3D 볼륨을 얻을 수가 있으나 현실적으로 이는 어려운 일이며, 이러한 이유로 얻은 3D 볼륨은 같은 지점을 여러번 스캔하여도 실제로는 다른 부분의 결과를 얻게 되거나, 연속적으로 얻어져야 할 위치들 사이에 빈 공간이 생기기도 한다.
이전의 연구들에서는 다른 하드웨어와의 연동을 통하여 해결하거나 cross-corrleation을 이용하는 등의 방법들이 제안되었다. 이 논문에서는 다른 하드웨어를 사용하지 않고 같은 영역에서 수행된 2개 이상의 3D볼륨 데이터를 이용하여 스캔 중에 있었던 움직임을 추정해보고자 한다. 이 때, 각 3D 볼륨 데이터들은 서로 uncorrelated 하도록 스캔의 방향이 달라야 한다.
2. Method
2.1 Scan process modeling
3D 볼륨은 한 지점의 A-scan 결과를 시간적으로 연속 측정한 결과를 모아 생성하게 된다. 따라서 시간 t에서 한 지점 $ p $에서의 A-scan 결과를 다음과 같은 함수로 표현해보자.
$$AScan(p, t) : R^3 \times R \rightarrow R^d$$
지점 $ p = (x, y, z)^T $ 는 스캔 장치를 기준으로 하는 좌표계 상에서의 스캔 위치를, t는 Ascan이 얻어진 시각 t를 나타낸다. 이 AScan 함수는 axial 방향의 인텐시티 값을 갖는 d 차원의 벡터를 결과로 하는데, 이 d는 OCT 시스템의 깊이 방향 해상도가 된다.
OCT는 스캔 지점 p를 망막을 가로지게 이동시키면서 스캔을 하게 되는데, 대상 물체가 크게 기울여지지 않는한 둘의 좌표계는 상대적으로 서로 일치하다고 볼 수 있다. 이런 과정으로 수식으로 나타내어보면, 각 위치 $ p_i $를 각각 해당되는 $ t_i $ 에서 스캔($ i = i ... N_{ASCANS} $) 하여 3D 볼륨을 얻게 된다. 각 $ p_i $는 움직임 보정 알고리즘에 따른 위치 이동 패턴을 띈다.
보통 3D 볼륨을 얻는 스캔 패턴은 x 방향의 B-scan 스캔을 각 y 좌표마다 반복하는 방식으로 얻기 때문에 x 방향으로는 빠른 스캔이 가능하나, y 방향으로는 상대적으로 느린 특정을 가진다. 이러한 방식에서의 빠른 스캔 방향을 $ d_{fast} $, 느린 스캔 방향을 $ d_{slow} $ 이라고 하고, $ d_{fast} $가 X축인 스캔을 XFAST, Y축인 스캔을 YFAST라고 이야기해보자. 이전의 움직임 보정을 위한 방법들은 복잡한 스캔 패턴을 가지는 방법들도 더러 있었으나 이 논문에서는 간단하게 XFAST 패턴과 YFAST 패턴을 반복하여 2개 이상의 볼륨을 얻는 패턴을 이용하도록 하자.
2.2 Modeling motion
3D 볼륨을 얻기 위한 모든 AScan을 얻는 동안 눈이 움직이면서 모션 아티팩트가 생기게 된다. 이러한 눈의 움직임을 나타내기 위하여 OCT 스캐너 기준의 좌표계를 이용하여 나타내어 보자.
$$Disp(t) : R \rightarrow R^3$$
이 함수는 특정 시각 t 에서 의도되었던 스캔 위치를 기준으로 실제 스캔된 위치가 얼마나 떨어져 있는지를 나타내는 3차원 벡터이다. 따라서 스캔 시작 시각을 0 으로 생각한다면 아래와 같이 나타낼 수 있다.
$$AScan(p, t) = AScan(p - Disp(t), 0)$$
이 함수 Disp는 스캐너와 물체 사이의 상대적인 위치 이동을 나타내기 때문에, 만약 물체의 이동이 없다면 Disp의 값들은 이정한 값으로 유지될 것이다.
2.3 Resampling
움직임을 보정하였을 때의 실제 의도된 스캔 위치 $ p_i $의 AScan 정보 $ AScan(p_i, 0) $ 를 계산하기 위하여 interpolation 함수
$$Interp(p_i) = ASCan(p_i, t_i) = AScan(p_i - Disp(t_i), 0)$$
을 정의해보자. 샘플 공간은 grid형태를 띄고 있으므로 이 위치 $ p_i $ 는 연속적인 interpolation을 이용하여 가상 AScan 결과를 생성해 낸다. 우측의 Disp를 제거하여 식을 아래와 같이 변형시킬 수 있다.
$$Interp(p_i + Disp(t_i)) = AScan(p_i, 0) + e_{interp}$$
여기서 $ e_{interp} $ 는 interpolation한 결과의 오차이나, 물체의 3D 볼륨이 빈 공간이 없이 꽉 차 있고, 각 위치에 따른 랜덤 노이즈 또한 비슷하므로 이는 무시할 만한 수준이다. 따라서 오직 움직임만이 높은 오차 e를 야기한다고 생각할 수 있다.
2.4 Objective function
이제 2개 이상의 3D 볼륨을 서로 수직인 $ d_{fast} $ 를 갖도록 연속적으로 얻도록 하자. 최종 목표는 어떤 알지 못하는 볼륨에 대하여 각 스캔의 $ Disp(t) $ 를 추정하는 것이다. 이를 위하여 전체 움직임을 아우르는 Objective function을 정의하도록 하자.
$$F(Disp(t)) = \sum_{v1=1}^{N_v} \sum_{v2=v1+1}^{N_v} \theta(v1, v2) || R_{v1, v2} || ^ 2 + \alpha \int_{t} || ( \frac{\delta Disp(t)}{\delta t})^T \begin{pmatrix} 1 && 0 && 0 \\ 0 && 1 && 0 \\ 0 && 0 && \sqrt{\alpha_z} \end{pmatrix} ||^2 \delta t$$
식은 전체 얻어진 3D 볼륨간 쌍 v1, v2에 대하여 계산된 값들의 합으로 이루어지는데, v1과 v2이 순서가 뒤바뀌어도 같은 값을 갖기 때문이 중복된 쌍을 배제하기 위한 고려가 되어있다.
첫번째 항의 R은 Disp(t)에 대하여 보정된 볼륨간의 차이를 나타내는 것으로 이는 간단하게 두 볼륨의 인텐시티 값의 L2-norm으로 계산된다.
두번째 항은 눈의 움직임 특성상 그 값이 급격하게 변하지 않고 어느정도 연속성을 띔을 고려한 것으로써 시간에 따른 그 움직임 정도를 regularizer로 이용한 것이다.
두 항 사이의 중요도를 $ \alpha $로 조절하였다.
이 함수를 최적화, 즉 보정된 3D 볼륨간의 차이를 가장 적게 하면서도, 그 움직임의 변화가 적은 움직임 Disp(t)을 찾아낸다면 그것이 움직임에 대한 프로파일이라 할 수 있고, 이를 통해 볼륨을 재구성한다면 움직임이 보정된 볼륨을 구할 수 있게 된다.
2.5 Optimization and multi-resolution technique
이렇게 식을 구성하였지만, 올바른 Disp(t)가 존재하는 공간은 너무 크다. 예를 들어 400*400 크기의 볼륨 2개를 최적화할 경우, 2*3*400*400개의 변수에 대한 계산을 해내야 한다. 따라서 Gaussian reduction을 이용한 resolution pyramid를 구성한 뒤, 낮은 해상도부터 해상도를 높여가며 반복적으로 계산을 수행한다.
각 단계에서의 최적화 방법은, L-BFGS (limited memory Broyden-Fletcher-Goldfarb-Shanno) 방법을 사용하였다.
2.6 Pre-processing
최적화를 수행하기 전에 노이즈에 대한 영향을 줄이기 위하여, Axial 방향으로만 Median 필터링을 적용하고, 그 크기를 반으로 줄이는 downsampling을 해 준다. 또한 픽셀 값은 평균이 0, 분산이 1이 되도록 정규화한다.
2.7 Volume merging
최적화된 결과를 바탕으로 가중치합, interpolation 등의 방법을 적용하여 최종 볼륨을 만들어낸다.
3 Evalution and discussion
3.1 Evalution
$ \alpha $ 와 $ \alpha_z $ 의 변화에 따른 성능 변화에 대한 실험도 수행하였다. SNR의 상승 또한 보였으며, 시각적으로 좋은 성능을 보였음을 보였다.