(2013)Comparison of Thresholding Techniques for SVD Coefficients in CT Perfusion Image Analysis

0. 요약

SVD(Singular Value Decomposition) 기반의 deconvolution 방식은 CT Perfusion image 해석에서 가장 널리 사용되는 알고리즘이다. 이 알고리즘은 noise의 영향을 줄이기 위해 SVD 계수를 thresholding하는 과정이 사용된다. 이 때 threshold 값으로 고정된 값을 사용하거나 미리 정해진 Oscillation Index에 따른 threshold 값이 사용된다. 본 논문에서는 두 thresholding 방식의 accuracy를 비교하기 위한 몬테 칼로 모의 실험 방식을 제안했다. 또한 CTP image 해석에서 사용하는 평활화 과정이 알고리즘의 정확도에 미치는 영향을 측정하기 위해 이 실험 방식을 확장하였다.

1. 서론

Perfusion data를 얻기 위해서 CT나 MR을 이용하는 경우 contrast를 주입한 후 30 ~ 80초 동안 1 ~ 2초 간격으로 취득된 30 ~ 40장의 slice images를 얻고 각 voxel 별로 contrast의 농도 변화를 추적하여 CBF(Cerebral Blood Flow), CBV(Cerebral Blood Volume), MTT(Mean Transit Time) 등의 Perfusion metric에 대한 map을 생성한다.

Perfusion 분석에서 가장 중요한 것은 cerebral artery과 tissue에서 얻은 두 개의 함수를 deconvolution하는 과정인데, SVD 기반의 알고리즘이 가장 효과적인 방식으로 알려져 있다. 최근에는 block circulant SVD 알고리즘이 자주 사용되고 있다.

Deconvolution은 inverse 문제를 푸는 것 이기 때문에 불안정한 연산이므로 noise의 영향을 많이 받는데 정확성을 높이기 위해 SVD coefficient에 대한 thresholding 과정이 필요하다.

Block circulant SVD 알고리즘에서는 coefficient의 thresholding을 위해 Oscillation Index를 사용하는 방식과 Fixed threshold를 사용하는 방식이 많이 사용되고 있다. 두 가지 방식을 비교하기 위해 본 논문에서는 실제와 유사한 조건으로 생성된 농도 함수를 이용한 몬테 칼로 모의 실험을 통해 두 thresholding 방식을 비교하였다.

CT에서는 contrast 농도 함수를 얻기 위해 현재 image와 배경 image의 차이를 구하는 image subtraction 알고리즘이 사용된다. 따라서 noise가 심하게 발생한다.

Noise를 제거하기 위해서는 filter를 이용해 평활화를 수행하게 된다.

본 논문에서는 평활화 과정이 SVD coefficient의 thresholding에 어떤 영향을 미치는지 알아본다.

2. 뇌관류 영상 해석 기법

1) 이론적 배경

CT Perfusion image는 contrast를 주사한 다음 1 ∼2초 간격으로 동일한 CT image를 30∼40장 정도 취득 하여 얻어진다.뇌의 각 부분에서의 perfusion parameter를 구하기 위해 cerebral artery에서의 농도 변화를 나타내는 AIF(Arterial Input Function) $C_a(t)$와 perfusion parameter를 구하고자 하는 voxel에서의 contrast 농도 변화를 나타내는 TF(Tissue Function) $C_v(t)$가 사용된다.

위 그림을 보면 1로 표시한 점이 cerebral artery 부분이고 2, 3번이 tissue이다.

TF $C_v(t)$는 AIF $C_a(t)$와 residue function $R(t)$와의 convolution의 결과로 식을 세운다.

$R(t)$는 해당 voxel에 남아있는 contrast의 잔량 비율을 나타내는데, contrast가 주입된 시점을 $t=0$인 순간으로 정의하면, 이 때 contrast 전체가 남아있어야 하므로 $R(0)=1$이다. 또한 시간이 지나 contrast가 모두 빠져나가면 0이된다.

위 세 함수는 다음의 관계를 가진다.

$C_v(t)=F_v(C_a(t)\otimes R(t))=F_v\int_0^tC_a(\tau)R(t-\tau)d\tau$

여기서 $F_v$는 CBF(Cerebral Blood Flow)라고 정의되는 parameter이다. CBF는 단위 시간동안 일정 크기의 조직을 통해 흐르는 혈 액의 양을 나타내는데, 보통 100g의 조직을 통해 1분 동안 흐르는 혈액의 양을 ml로 나타내므로 ml/100g/min, 또는 조직의 양을 부피로 나타낼 때 ml/100ml/min 의 단위를 가진다. 아래 그림은 CBF map 예시이다.

CBF 이외의 parameter로 CBV(Cerebral Blood Volume)와 MTT(Mean Transit Time)가 있는데 다음과 같이 정의된다.

$CBV=\frac{\int C_v(t)dt}{\int C_a(t)}$

$MTT=\frac{CBV}{F_v}$

CBV는 해당 복셀에서 혈액이 차지하는 양을 나타내는 것으로 일반적으로 ml/100g, 또는 조직의 무게 대신 부 피를 사용하는 경우 ml/100ml의 단위를 가지게 된다. MTT는 해당 복셀 영역을 혈액이 통과하는 평균시간을 초 단위로 나타내는데, CBF와 CBV를 알면 할 수 있다.

2) SVD를 이용한 파라미터 추정

CBF는 일반적으로 $C_v(t)$와 $C_a(t)$간의 deconvolution을 통해 $F_vR(t)(\equiv k(t))$를 구한 다음, $R(t)$의 최대값은 1이라는 성질을 이용하여 $F_v=\max[k(t)]$에 의해 구해진다.

$C_v(t)$와 $C_a(t)$의 길이를 $N$이라 하면 다음과 같이 discretized 형식으로 나타낼 수 있다.

$C_v(t_j)=\Delta t\ \cdot\ F_v\sum\limits_{i=0}^jC_a(t_i)R(t_j-t_i),\ (j=0,1,...,N-1)$

여기서 $\Delta t$는 CT image를 얻는 시간 간격을 나타낸다. 위 식을 행렬로 표현하면

$\begin{bmatrix} C_v(t_0) \\ C_v(t_1) \\ \vdots \\ C_v(t_{N-1}) \end{bmatrix}=\Delta t\begin{bmatrix} C_a(t_0) & 0 & \cdots & 0 \\ C_a(t_1) & C_a(t_0) & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ C_a(t_{N-1}) & C_a(t_{N-2}) & \cdots & C_a(t_0) \end{bmatrix}\times\begin{bmatrix} R(t_0) \\ R(t_1) \\ \vdots \\ R(t_{N-1}) \end{bmatrix}\ \cdot\ F_v$

또는 간략하게 $C=Ar$로 표현할 수 있다. 여기서 $A$는 $\Delta t$와 $C_a(t_i)$ 관련 element들을 포함하는 $N\times N$ 크기의 행렬, $r$은 $R(t_i)$에 $F_v$를 곱한 $N\times1$ 행렬, $C$는 $C_v(t_i)$ 관련 $N\times1$ 행렬을 의미한다.

행렬 $A$의 rank는 보통 $N$보다 작은데, $r$을 구하기 위해 SVD 방식을 이용한다. 행렬 $A$는 $A=USV^T$로 분할할 수 있는데 이 때 $A$의 역행렬은 $A^{-1}=VWU^T$가 된다.

마지막으로 $r=A^{-1}C=VWU^TC$에 의해 $r$을 구하고 최대값을 취함으로 CBF인 $F_v$를 결정하게 된다.

3) 블록 순환형 행렬을 이용한 SVD 방식

CT Perfusion 연구가 더욱 진행되면서 tissue에서 voxel의 위치에 따라 contrast가 도달하는 시간이 달라지는 현상이 관찰되었다. 지연 시간의 차이에 의해 TF에 dispersion이 발생되면 측정된 perfusion parameter에 오차가 발생하게 되는데, 이러한 오차를 줄이기 위한 방법으로 행렬 $A$를 block circulant matrix로 정의하는 방식이 제안되었다. 우선 zero-padding 기법으로 AIF, TF의 길이를 $L(L\geq2N)$로 늘린다음, 행렬 $A$를 $L\times L$ 크기의 block circulant matrix $\widetilde{A}$로 대체한다. 행렬 $A$의 element를 $a_{ij}$라 할 때, $\widetilde{A}$의 element $\tilde{a}_{ij}$는 다음과 같다.

$\tilde{a}_{ij}= \begin{cases} a_{ij}, &j\leq i \\ a_{L+i-j,0}, &j>i \end{cases}$

마찬가지로 $\widetilde{A}$도 SVD를 적용하여 $\widetilde{A}=\widetilde{U}\widetilde{S}\widetilde{V}^T$로 분해할 수 있다. Zero-padding된 TF의 벡터를 $\widetilde{C}$라 하면, residue function은 $\widetilde{r}=\widetilde{V}\widetilde{W}\widetilde{U}^T$로 주어진다.

3. CT 영상 해석에서의 SVD 계수 임계화 방식

일반적으로 행렬 $A$의 rank는 $N$보다 작으므로 $r$은 noise의 영향으로 인해 실제값과 다르다.

SVD를 적용 시 noise의 영향을 줄이는 방법은 $\widetilde{S}$에서 $\widetilde{W}$를 구하기 위해 다음과 같은 thresholding 과정을 수행하는 것이다.

$\widetilde{W}_{ii}= \begin{cases} 1/\widetilde{S}_{ii}, &\text{if}\ \widetilde{S}\geq T_{SVD} \\ 0, &\text{otherwise} \end{cases}$

여기서 $T_{SVD}$는 threshold 값으로 이 값이 커질수록 residue function의 $r$의 변동이 줄어들게 된다. $T_{SVD}$는 일반적으로 행렬 $\widetilde{S}$의 element 중 최대값인 $\widetilde{S}_{11}$에 비례하는 값으로 선택한다.

$T_{SVD}=P_{SVD}S_{11}\ (0<P_{SVD}<1)$

일반적으로 CT image에서는 voxel의 noise가 심할수록 큰 $P_{SVD}$를 사용해야 한다.

1) 고정 경계치를 이용한 임계화 (FT 방식)

가장 간단한 thersholding 방식은 전체 image에 고정된 $P_{SVD}$를 사용하는 것이다. Image의 모든 voxel에 동일한 $P_{SVD}$를 사용하여 $\widetilde{W}$를 구하고 $k(t)$를 구한다. 연산 시간이 최소화된다는 장점이 있지만 정확도가 떨어진다는 단점을 가진다.

2) 진동지수(OI)에 의한 임계화 (OI 방식)

Voxel에 따라 noise level이 변하는 경우에 $T_{SVD}$를 변화 시키는 방법으로 OI(Oscillation Index) 방식이 있다. OI는 추정된 $k(t)$의 변화가 얼마나 심한지 나타내는 계수로 다음과 같다.

$OI=\frac{1}{L}\frac{1}{k_{\max}}(\sum\limits_{n=1}^{L-2}|k(n-1)-2k(n)+k(n+1)|)$

OI 방식은 각 voxel에서 OI가 일정하도록 threshold 값을 정한다. 우선 충분히 작은 $T_{SVD}$를 적용하여 $\widetilde{W}$를 구하고 $k(t)$를 구한 다음, OI를 계산한다. 이 값이 원하는 OI 보다 작아질 때 까지 $T_{SVD}$를 조금씩 증가시키면서 이 과정을 반복한다.

4. 임계화 방식의 성능 비교 방법

1) 모의 실험 방법

CT Perfusion image의 경우 검증된 data 값은 존재하지 않다. 이러한 어려움을 극복하기 위해 몬테 칼로 모의 실험 방식을 이용한 성능 평가를 수행한다. 모의 실험을 위해 우선 AIF를 생성하는데, 이 함수는 일반적으로 변형 감마 함수 형태를 가지므로 다음의 모형 함수를 사용한다.

$C_a(t)= \begin{cases} 0, &t<t_0 \\ C_0(t-t_0)^\alpha e^{-(t-t_0)/\beta}, &t\geq t_0 \end{cases}$

여기서 $t_0$는 주입된 contrast가 cerebral artery 도착하는 시간을 나타내고, $\alpha$와 $\beta$는 파형의 세부 형태를 표현하는 parameter이다. 다른 연구들과 유사하게 $t_0=5s,\ \alpha=3.0,\ \beta=1.5s$ 등으로 설정하였다.

TF는 AIF와 residue function의 convolution에 의해 생성한다. 이 대 TF를 얻으려면 residue function을 정의해야 하는데 다음과 같은 지수 함수 형태를 사용하였다.

$R(t)=e^{-t/r}$

TF를 생성하려면 우선 perfusion parameter인 CBF와 MTT를 지정해야 한다. 일반적으로 tissue는 CBF의 범위가 10-60 ml/100 ml/min, MTT는 2-12 s 정도가 된다. 지정된 $\tau$ 값과 함수 유형에 의해 $R(t)$가 결정되면 이 함수를 1초 간 격으로 표본화한 다음, convolution을 적용하여 TF $C_v(t)$를 생성한다.

이렇게 생성된 $C_a(t)$와 $C_v(t)$에 주어진 SNR에 맞추어 gaussian noise를 추가하면 분석하고자 하는 입력 신호가 생성된다. 두 신호에 gaussian noise를 추가하였는데, noise의 정도를 신호대 noise비 SNR을 이용하여 나타낸다. TF의 최대치가 $C_{max}$, gaussian noise의 표준 편차가 $\sigma_n$ 일 때 TF의 $SNR_C$를 다음과 같이 정의한다.

$SNR_C=\frac{C_max}{\sigma_n}$

위 그림은 생성된 TF $C_v(t)$의 샘플을 보여주는데, 여기서는 noise가 없을 때 와 중간 수준(SNR=8), 그리고 심한 수준(SNR=2)의 noise가 추가된 신호를 나타내고 있다.

생성된 입력 신호를 평활화한 다음, 두 신호에 대해 SVD 방식을 적용하여 perfusion parameter $\hat{F}$와 $\hat{\tau}$를 추정하고, 이 값들을 원래 값과의 오차를 구해 각 평활화 방식의 성능을 평가한다. 본 실험에서는 CBF를 성능 평가 지표로 사용하였는데, 측정된 CBF의 오차는 다음과 같이 정의한다.

$E_F=\frac{|F-\hat{F_v|}}{F_v}$

2) 영상 평활화 효과의 생성

모의 실험 방식에서는 각 voxel의 TF를 독립적으로 처리하여 parameter를 측정한 다. 그러나 실제 CT perfusion에서는 일반적으로 여러 voxel에 걸치도록 설정하는 경우가 많다. 또한 noise의 영향을 줄이기 위해 입력 image에 대해 평활화 filter를 적용하게 된다.

본 연구에서는 perfusion parameter가 일정하게 유지되는 영역에 서로 독립적이고 균일한 분포의 noise가 추가되었다고 가정한 모의 실험을 수행하였다. 이 상황에서 $M\times M$크 기의 box filtering이 수행된다면, AIF와 TF의 신호 성분은 변하지 않고, noise 성분만 filtering되는 효과 가 나타나게 된다. 따라서 box filter에 의한 평활화 효과를 내기 위해 다음과 같은 방법을 사용하였다.

마찬가지로 noise가 없는 신호 $C_a(t)$와 $C_v(t)$를 생성한다. 다음 단계에서 주어진 SNR에 맞추어 $2M^2$개의 gaussian noise를 생성하고, 이들 noise을 $M^2$개씩으로 구성된 두 개의 그룹으로 나 눈 다음, 각 그룹의 평균을 $C_a(t)$와 $C_v(t)$에 더한다.

5. 실험 결과

1) 모의 실험 결과

신호의 SNR에 따른 최적의 threshold 값의 변화 그래프이다. 각 SNR에서 일 정한 CBF와 2∼12 사이의 임의의 MTT를 가지는 1,000 개의 AIF와 TF를 생성한다. 이와 같이 생성된 함수들에 대해 SVD를 수행하고 $P_{SVD}$를 0 ∼0.5 사이에서 0.005 간격으로 변화 시키면서 CBF를 측 정한다. 각각의 threshold 값을 이용하여 구한 CBF에 대해 오차 $E_F$를 측정하여, 이를 최소화하는 threshold value $\hat{P}_{SVD}$를 구한다. 예상하는 바와 같이 SNR이 증가할수록 threshold value가 낮아지는 것을 볼 수 있다.

각 threshold 값에서 SNR과 오차와의 관계를 나타내는 그래프이다. FT 방식의 경우 신호의 SNR이 10 이상 되는 경우 threshold 값이 작을수록 결과의 오차가 작지만, 낮은 SNR에서 작은 threshold 값(0.05)는 큰 오차를 초래하는 것을 볼 수 있다. OI 방식의 경우 threshold 값 0.2 정도에서 최선의 결과를 산출하고, 이보다 커지면 오차가 다시 증가하는 현상이 관찰되어, 모든 실험에서 0.1과 0.2의 threshold 값을 적용하였다.

입력 영상을 3x3와 5x5 크기의 box filter를 이용한 평활화 과정이 결과의 정확도에 어떤 영향을 미치는지 알아본다. 왼쪽 그래프가 3x3, 오른쪽 그래프가 5x5이다. 앞서 설명한 바와 같이 9개와 25개의 잡음 신 호의 평균을 샘플 신호에 더한 다음, FT와 OI 방식의 thresholding 과정을 적용하였다.

결과를 검토하면 공통적으로 다음과 같은 현상들이 관찰된다.

• FT 방식의 경우 SNR이 작은 경우에는 threshold 값이 높을수록 유리하지만, SNR이 어느 이상(5∼10 정 도) 커지면 threshold 값이 낮을수록 CBF 오차가 줄어든다.
• SNR이 어느 이상 (10 정도) 커지면 FT 방식의 오차는 비교적 일정하게 유지되지만, OI 방식의 경우 SNR에 따라 오차가 계속 줄어드는 현상을 보인다.
• SNR이 일정 수준 (5∼15 정도) 이상에서는 OI 방식 의 오차가 FT에 비해 작다.
• FT 방식의 경우 보다 정확한 결과를 얻기 위해 신호의 SNR과 평활화 filter의 크기에 따라 다른 threshold 값을 사용하여야 하지만, OI의 경우 고정된 threshold 값을 적용할 수 있다.
• SNR이 어느 이상 레벨(평활화 정도에 따라 5∼15) 에서 OI 방식이 FT 방식에 비해 정확도가 높다.
• OI 방식을 적용하더라도 측정된 CBF는 10∼15% 정도의 오차범위를 보인다.

2) 실제 CT 영상을 이용한 실험 결과

실험에 사용된 영상은 2초 간격으로 취득된 36장의 CT image로 구성된다. Noise의 영향을 줄이기 위해 각각의 영상을 5x5 box filter로 평활화하였다. 또한 각 voxel에서의 농도 함수값은 다시 5x5 영역에서의 평균값으로 얻었다. 결과적으로 농도 함수를 얻기 위해 5x5 평균 filter를 두 번 적용한 것과 같다.

본 연구에서는 slice image에서 뼈 부분을 제외한 영역에서 CBF와 MTT를 구한 다음, 그 값이 일정한 범위 내에 들어오는 voxel들에 대해서만 CBF 값을 기록하였다. 이 perfusion 해석 절차를 전체 voxel에 적용 하여 CBF 값들을 얻은 다음, 다시 5x5 크기의 평활화 필터를 적용하여 최종 CBF map을 구했다.

위 그림은 threshold 값을 0.12로 설정하여 얻은 OI 방식의 결과를 보인 것이다. CBF 값은 0∼150 ml/100ml/min 범위에 분포하고 있는데, 위 그림에서는 각 voxel에서의 CBF 값을 pseudo-color로 나타낸 것이다. 이 결과에서는 tissue에서 혈관과 백질(white matter), 회색질(gray matter) 등이 잘 구분되어 임상적으로 올바른 결과라고 판단된다.

위 그림은 FT와 OI 방식의 결과를 비교해보기 위해 같은 image에 대해 4개의 서로 다른 threshold 값을 사용하여 FT 방식을 적용해 보았다. 각각의 voxel에서 OI와 FT 방식으로 구한 CBF 값을 점으로 나타냈는데, 두 결과의 차이가 적을수록 점이 대각선에 가까이 오게 된다.

FT와 OI 결과와의 상관계수 를 구했는데, 그 결과가 위 표에 나와 있다. 표에서는 FT 가 0.08과 0.12 일 때 두 thresholding 방식의 상관관계가 매우 높은 것을 볼 수 있다. 결과적으로 실제 영상에서 FT로 얻은 결과는 threshold 값에 따라 OI와 매우 유사한 결과를 얻을 수 있지만, 부분 적으로 일치하지 않는 점들이 있는 것을 알 수 있다.

6. 결론

본 논문에서는 SVD 기반 CT Perfusion image 해석에서 계수의 thresholding 방식들의 정확성을 비교하였다. FT와 OI 방식의 정확도를 측정하기 위 한 몬테 칼로 모의 실험 방법을 제안하였고, 이 방식을 확장하여 image의 평활화가 perfusion parameter의 정확도에 미치는 영향을 측정하는 실험 기법을 제안하였다. 모의 시험을 수행한 결과, OI 방식이 FT 방 식에 비해 전반적으로 정확한 결과를 만들어준다는 점을 확인하였다. 또한 FT 방식에서는 보다 정확한 결과를 얻기 위해 신호의 SNR에 따라 다른 threshold 값을 적용 해야 하는 반면, OI 방식에서는 고정된 threshold 값을 사용 할 수 있음을 보게 되었다. 그러나 FT 방식은 계산량 측면에서는 OI 방식에 비해 우수하다.