Deep Learning for undersampled MRI reconstruction

Introduction

  자기 공명 영상(Magnetic Impedance Imaging)은 강한 자기장 내의 인체에 전파를 전사해서 반향되는 자장을 측정하여 영상을 얻는 방법이다.  구체적으로, MRI 기계에서는 frequency와 phase 두 가지 방향으로 데이터를 수집한다.  frequency encoding 방향으로는 빛의 속도에 가깝게 수집되지만 phase encoding 방향으로는 물리적인 이유로 인해 상대적으로 느리게 측정된다. 이렇게 측정된 데이터를 수학적인 변환(Fourier transform)을 이용하여 변환하면 원하는 영상을 얻을 수 있다. 하지만,  기존의 MRI 방식은 이미지를 얻기 위해 필요한 데이터 수집 시간이 다른 영상 처리 방식보다 느리다는 단점이 있고 이러한 단점을 극복하기 위해 phase encoding 방향의 측정을 적게 하면서 본래의 MRI영상을 얻고자 하는 노력을 여러 사람들이 하고 있다. 이러한 영상처리 개선 방식을 총칭하여 undersampled MRI reconstruction이라고 한다.  우리는 deep learning을 이용하여 undersampled MRI reconstruction problem를 해결 하고자 한다.

 

Method

  Undersampled MRI reconstruction problem을 수학적으로 표현하기 위해 다음과 같은 변수들을 정의하자.

  • Fully-sampled k-space data is denoted by x_{\mathrm{full}}.
  • Subsampled k-space data is denoted by x.
  • Subsampling function from x_{\mathrm{full}} to x is denoted by \mathcal{S}.
  • Reconstruction image from fully sampled data is denoted by y.
  • Direct reconstruction image from subsampled data is denoted by y_{_{\mathcal{S}}}

위 정의와 함께 Undersampled MRI reconstruction problem의 수학적인 구조는 아래와 같은 그림으로 간략히 표현 할 수 있다.

1

 

다음과 같은 forward problem를 고려하자.

S \circ \mathcal{F}^{-1} (y) = x

forward problem은 데이터를 어떻게 subsampling을 하는가에 따라서 의존한다. 그러므로, subsampling function \mathcal{S}가 forward problem를 구체화 시키고 구체화된 문제에 해당하는 inverse problem이 우리가 풀고자하는 undersampled MRI reconstruction problem(finding f : x \rightarrow y)에 상응한다.

 

1) Subsampling Strategy

Inverse problem을 살펴보기 전에 Forward problem의 well-definedness에 대해서 살펴보아야 한다. 주어진 forward problem은 subsampling function \mathcal{S}에 의해서 결정이 되기에  적절한 \mathcal{S}가 주어진 경우에 inverse problem을 풀 수 있을 것이다.  아래의 예시가 어떤 \mathcal{S}가 적절한지에 대한 직관적인 이해를 도움을 준다. 예를 들어,  f를 전통적인 방식의 복원 변환인 Inverse Fourier transform에 절댓값을 취하는 함수라고 가정하고 \mathcal{S}_0=100% fully sampling, \mathcal{S}_1=50% uniform subsampling, \mathcal{S}_2=50% uniform subsampling added some low frequency에 대한 복원 결과는 Figure 1과 같다.

2

Figure 1. 각각의 이미지는 (a) \mathcal{S}_1 (b) \mathcal{S}_0 (c) \mathcal{S}_2에 의해서 복원된 이미지

Figure 1는 뇌암이 있는 환자의 MRI 영상을 이용하여 얻은 이미지 들이다. (b)에 해당하는 영상은 100% sampling을 통해서 완벽하게 복원한 실제 환자의 MRI 영상이고 (a)와 (c) 이미지는 subsampling으로 얻은 영상들이다. (a)와 (c)에서 확인 할 수 있는 두 subsampling의 큰 차이점은 환자의 뇌에 있는 암세포의 위치가 (a)인 경우에는 어느 곳에 암세포가 있는지를 결정 할 수 없지만 (c)의 경우에는 위치를 결정 할 수 있다.  수식적으로 forward problem의 separability로 해석 할 수 있다.

we say h has separability if h(x) \neq h(y) for any x,y satisfying x \neq y.

그러므로, Figure 1에서 얻은 결과로 해석해보면 uniform sampling에 약간의 low frequency를 더해준 경우에 data의 separability를 보장해줄 것이라고 생각할 수 있다. 그래서, 실제 implementation에서는 25% uniform sampling에 4% low frequency를 더해준 29%의 subsampling을 \mathcal{S}로 사용하였다.

 

2) Reconstruction map f

우리는 undersampled MRI reconstruction f를 딥러닝 기술을 이용하여 구현하려고 하였다.  f를 구성하는 구조는 U-net을 이용한 심층학습 과정과 측정한 데이터 보존을 위한 k-공간 수정 과정 두 가지로 나눌 수 있다. 구체적으로, f는 다음과 같이 정의 할 수 있다.

f = |\mathcal{F}^{-1}| \circ f_{cor} \circ \mathcal{F}\circ f_d \circ |\mathcal{F}|^{-1} 

여기에서 f 는 trained u-net, f_{cor}k-space correction function 이고 \mathcal{F}는 Fourier transform이다. 아래의 알고리즘 도식도를 보면  함수 f의 구조를 직관적으로 쉽게 이해 할 수 있다.

3.JPG

Experiment

우리가 설계한 f를 가지고 얻은 복원 결과는 다음과 같다.

4.JPG

 

written by Chang Min Hyun

Update data : 11 September 2017

# 더욱 자세한 이론, 결과 해석은 논문을 참고하면 얻을 수 있습니다.

Reference

[1] E.J. Cand`es, J. Romberg and T. Tao,“Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information,”  IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489–509, 2006.

[2] D.L. Donoho, ”Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52, pp. 1289–1306, 2006.

[3] D.L. Donoho, ”For most large underdetermined systems of linear equations the minimal 1-norm solution is also the sparsest solution,” Communications on pure and applied mathematics, vol. 59, pp. 797–829,  2004.

[4] Google, ”TensorFlow: Large-scale machine learning on heterogeneous systems,” URL http://tensorflow.org/, 2015.

[5] K. Hammernik, T. Klatzer, E. Kobler, M.P. Recht, D.K. Sodickson, T. Pock and F. Knoll, ”Learning ad Variational Network for Reconstruction of Accelerated MRI Data,” arXiv preprint, arXiv:1704.00447, 2017.

[6] T. Tieleman and G. Hinton, ”Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude,” COURSERA: Neural Networks for Machine Learning, 2012.

[7] D.J. Larkman and R.G. Nunes, ”Parallel magnetic resonance imaging,” Phys. Med. Biol., vol. 52, pp. R15–R55, 2007.

[8] P.C. Lauterbur, ”Image Formation by Induced Local Interactions: Examples of Employing Nuclear Magnetic Resonance,” Nature, vol. 242, pp. 190–191, 1973.

[9] D. Lee, J. Yoo, J.C. Ye, ”Deep artifact learning for compressed sensing and parallel MRI,” arXiv preprint, arXiv:1703.01120, 2017.

[10] M. Lustig, D.L. Donoho and J.M. Pauly, ”Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging,” Magnetic Resonance in Medicine, Vol. 58, pp. 1182–1195, 2007.

[11] H. Nyquist, ”Certain topics in telegraph transmission theory,” Trans. AIEE, vol. 47, pp. 617–644, 1928.

[12] K.P. Pruessmann, M. Weiger, M.B. Scheidegger and P. Boesiger, ”SENSE: sensitivity encoding for fast MRI,” Magn. Reson. Med., vol. 42, pp. 952–962, 1999.

[13] O. Ronneberger, P. Fischer, and T. Brox, ”U-net: Convolutional networks for biomedical image segmentation,” in Int. Conf. on Medical Image Computing and Computer-Assisted Intervention, Springer, pp. 234–241, 2015.

[14] J. K. Seo and E. J. Woo, ”Nonlinear inverse problems in imaging,” Chichester, U.K.: John Wiley & Sons, 2013.

[15] J. K. Seo, E. J. Woo, U. Katscher, and Y. Wang, ”Electro-Magnetic Tissue Properties MRI,” Imperial College Press, 2014.

[16] D.K. Sodickson and W.J. Manning, ”Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays,” Magn. Reson. Med., vol. 38, pp. 591–603, 1997.

[17] Z. Wang, A. C. Bovik, H.R. Sheikh, E.P. Simoncelli, ”Image Quality Assessment: From Error Visibility to Structural Similarity,” IEEE Trans. on Image Processing Vol. 13, no. 4, pp. 600-612, 2004.

[18] S. Wang, Z. Su, L. Ying, X. Peng, S. Zhu, F. Liang, D. Feng and D. Liang, ”Accelerating magnetic resonance imaging via deep learning,” in Proc. IEEE 13th Int. Conf. Biomedical Imaging, pp. 514–517, 2016.

Fourier Transform

1. 1-Dimension Discrete Fourier Transform

1.1. 1D Discrete Fourier Transform

Why do we need to do Discrete Fourier Transform(DFT)? A Discrete Fourier Transform has been proposed since we cannot apply Continuous Fourier Transform in the real world. Specifically, a continuous fourier transform requires us to integrate from negative infinite to the positive infinite, which is impossible. Therefore, countable samples should be set in the domain so as to apply discrete fourier transform. This process converts a wave or an image from time domain to continuous domain. (By doing this process, we can get some benefits: saving data, extract the information we want, etc)

continunous: \displaystyle X(F) = \int_{-\infty}^{\infty} x(t) e^{-i2\pi Ft} dt

discrete : \displaystyle X_k = \sum_{n=0}^{N-1}{x_n \cdot e^{- {i2\pi kn}/N}}

where \displaystyle F \cong k/N and \displaystyle (t \cong n)

F: Frequency  \\  k: k^{th} sample  \\  N : Sample Frequency  \\  t: time

그림3

Let’s take a look at a simple example of sine wave. This sine wave is y = sin(2 \pi x) , whose frequency is 1 Hz. In order to see the simplest version of Discrete Fourier Transfrom, let’s assume that we are dividing 1 period sine wave into 4 samples. x_0 = 0, x_1 = 1, x_2 = 0, x_3 = -1
Now, we want to transform these four points into new points by Discrete Fourier Transform,

\displaystyle X_k = \sum_{n=0}^{N-1}{x_n \cdot e^{- {i2\pi kn}/N}}

Specifially, take a look at x_1 be transformed into X_1

x1

Similarly, computations would be made on X_0, X_2, X_3, and this can be expressed visually by matrix form.

그림2.jpg

Then, as we can see the below, the time domain – input sequence has been transformed into the frequency-domain – magnitude response.

그림4.jpg

But, here we can come up with a question:
Why did 3Hz frequency mapping has showed up although there were only 1Hz wave?
⇒ The answer can be found from the “sampling order”
To simply illustrate this, below images could be a good example.

그림as1

Originally, while we do Fourier Integral or Fourier Transform, the domain is from -\infty to \infty. As we assume that we support the sine wave into [-0.5, 0.5], the sample numbering should be done in the order of [x_0, x_1, x_2, x_3] as above. However, actual Discrete Fourier Transform is sampled in the order of [x_2, x_3, x_0, x_1]. Thus, we can expect that the DFT results will be symmetric to the continuous version of Fourier Transform or Fourier Integral.그림13

This way of calculation concludes the shift of the ramification. As we can see below, the half of the outcome on the negative frequency domain is shifted to the positive frequency domain. (Since this always happens to Fourier Transform, scientists usually use “fft shift” function to visualize it)

그림1231

Since the image is always symmetric with respect to the half of the Sample Frequency. It is called a Nyquist Limit.

Nyquist Limit = Sample Frequency/2

Likewise, sin(2\pi x) + sin(4\pi x) wave could be transformed into the frequency domain, resulting frequency 1Hz, and 2Hz.

그림7

그림8그림9.jpg

Thus, the actual frequency domain until Nyquist limit is called the “Nyquist Domain” and the copied frequency is called the “Periodic Mirrored Mapping”

There is an another question:
What happens if we have more sample frequency?
Let’s look at the below discrete transform with 16 samples.
그림rasjkrj2

The result can be compared to the first graph with 4 samples.
Since the above samples has 4 times more samples, the intensity is 4 times stronger. And the frequency range has become wider.
Moreover, the Nyquist limit is moved from 2Hz(4/2=2) to 8Hz(16/2).
Therefore, the copied frequency appeared on the 15Hz.
(In medical imaging, lower frequencies are more crucial than higher frequencies. Hence, fourier shifting is essential to handle the information easily.)

1.2. The General Form of Discrete Fourier Transform

Generalizing the samples with an integer N yields the following relations.

그림22313

1 Dimension Discrete Fourier Transform can be generalized as a matrix form by substituting W = e^{-{ i2\pi}/N}.

1.3. Inverse Discrete Fourier Transform

How can we implement Dicrete Fourier Transform of the frequency signals?
Surprisingly, by applying Inverse Discrete Fourier Transform, we can perfectly restore the information we have transformed. The relations can be generalized from the below example.

Since F_N's complex conjugate is \displaystyle \overline{F_N} = \sum_{n=0}^{N-1} x_n \cdot \overline{w^{nk}} , the simplest multiplication of F_2 and \overline{F_2} is as below.

그림2333

Did you find the simple relation of F_2 and \overline{F_2} ?
If we generalize that to the integer N, (we skip the proof of it using mathematical induction) we can elicit the below relation.

그림2abc.jpg

Finally, we get the Inverse Discrete Fourier Transform.

idft.jpg

In summary, calculating complex conjugate of F_N can be simply done by just changing the sign of the Discrete Fourier Transform. And after that, dividing by N, sample numbers, results in the form of Inverse Fourier Transform matrix.

1.4. 1D Fast Fourier Transform

Since the discrete transform’s cost is too expensive(N^2), An alternative method is developed by Yates, Frank which is extremely faster than that of Discrete Fourier Transform.^{[3]} The following example can be helpful to clarify this.

(Example Matrix with 4 samples)

그림2afadfasdfas

The strategy of Fast Fourier Transform is seperating even terms and odd terms to lower the computing cost.
The first step is changing the odd column and even column to gather same parities.
The second step is using the characteristic of complex plane. Since W = e^{ - {2\pi i}/N} , W^4 = e^{ {- 2\pi i /4} \cdot 4} = e^{-2\pi i} = e^{0} = 1 = W^{0}, the higher terms of W can be simplified to lower terms.
The third step is simplifying the matrix to lower the cost of computation.

Computating cost is changed from N^2 to Nlog_{2} N.
The more samples we have, the stronger the FFT is.

The Generalized Fast Fourier Transform can be illustrated as below diagram.

그림222.jpg

To implement FFT, we need to seperate the even terms and odd terms first. And then adding and subtracting each pair of even parities and odd parities will outcome the beautiful ramifications.

2.1. 2D Discrete Fourier Transform & Fast Fourier Transform

The actual goal of conducting Fourier Transform is applying it to 2 Dimension images such as MRI images or CT images. It sounds hard, but the algorithm is actually simple if we understood 1 dimension algorithm.

그림2adasdasda.jpg

It is nothing but just doing 1D DFT/FFT 2 times to rows and to columns.
But one thing we need to notice is that as we can see from the above algorithm, the lower frequency is centered at the each corners (as we remind 1D Fourier Transform results are biased), especially left top corner.
So, by fourier shifting( the built in function is fftshift2(img) in matlab), the lower frequencies are shifted to the center of the image and the higher frequencies are located to the edge.

(The simple example of 2D DFT & FFT)

그림2exa

As we transform the sine waves, the DFT/FFT catches the 1 hz frequency each. The x-component is of the row frequency, and the y-component is of the column frequency.
If we apply Inverse Fourier Transform, it perfectly recovers the original image.

(The example of phantom(512 * 512))그림2azxcxcxc.jpg

Here we can think of another question in the current research topic:
Can we recover the image even though we don’t have 100% samples from frequency domain?
We come up with this question because when we first get images from CT scan, the results are handed in the form of 2 dimension frequency images. And then recover them to the real image.
However, it takes almost 1 hour to scan a human body without any movement. If the woman moves a little bit, or sneeze, the image can be completely distorted. This is a huge loss in monetary, and time aspect.

If we Inverse Fourier Transform the 50% sampled frequency image, the result is weird as above. This effect can be understood by Poisson Summation Formula. (We do not handle this in this page)
This problem has been solved recently by using deep learning technique. If you are interested, I recommend you to refer the research of Changmin Hyun, Department of Computational Science and Engineering, Yonsei University.

 

Writer: Keunsik Jung, Youngbin Cho, Hyunchul Cho

 


Reference

  1. Discrete Fourier Transform – Simple Step by Step. Simon Xu. https://www.youtube.com/watch?v=mkGsMWi_j4Q&t=328s
  2. Kreyszig, E. (2006). Advanced engineering mathematics (10th ed.). Hoboken, NJ: Wiley.
  3. Yates, Frank (1937). “The design and analysis of factorial experiments”. Technical Communication no. 35 of the Commonwealth Bureau of Soils.
  4. Fast Fourier Transform. (n.d.)l. in Wikipedia. Retrieved August 30, 2017,  from https://ko.wikipedia.org/wiki/%EA%B3%A0%EC%86%8D_%ED%91%B8%EB%A6%AC%EC%97%90_%EB%B3%80%ED%99%98

전기임피던스 단층촬영 기술을 이용한 심폐기능 영상화

Introduction

연간 5천만명의 환자들이 인공호흡기의 도움을 받아 호흡하고 있다. 인공호흡기는 주입할 공기의 양을, 압력 또는 부피 파라미터를 이용해 조절한다. 하지만, 이는 환자 개개인의 폐 용적량을 고려하지 않은 값으로, 과도한 양의 공기가 주입될 수도 있다. 과도하게 주입된 공기는 폐포의 손상을 유발하고 이로 인한 제2차 감염으로 사망 할 수도 있다. 인공호흡기에 의한 폐포 손상을 줄이고 적절한 공기 주입량을 정하기 위해서, 공기 주입에 따른 폐 용적량의 변화를 영상화할 필요가 있다. 이러한 실시간 심폐기능 영상화는 전기 임피던스 토모그래피(Electrial impedance tomography, 이하 EIT) 기술을 사용하면 가능하다.

draeger_EITdevice
그림 1. Draeger사의 상용화된 EIT기기 PulmoVista500. 환자의 흉부에 임피던스 측정을 위해 전극이 부착되어 있다.

EIT는 비침습적이며 방사선을 사용하지 않아 인체에 무해한 영상 기법으로 생체전기(bioimpedance)로부터 인체 내부의 전기 도전율을 영상화하는 기술이다. 숨을 들이쉴때, 폐의 용적이 증가하고 공기가 폐 안으로 유입된다. 이때, 공기는 전기가 통하지 않으므로 폐의 도전율은 감소한다. 즉, 인체 내부의 도전율 분포 변화로부터 폐의 크기 변화를 알 수 있는 것이다. EIT는 환자의 흉부에 여러 개의 전극을 부착하여, 소량의 전류를 인가하고 전압을 측정한다. EIT를 이용한 심폐기능영상은 측정된 전압 데이타로부터 도전율 분포를 영상화 함으로써, 폐의 용적 변화를 시각화하는 것을 목표로 한다.

imlkh1

Measurement

전형적으로 EIT는 환자의 흉부의 다음과 같이 16 또는 32개의 전극을 붙여서 임피던스를 측정한다. 그림 1과 같이 부착된 전극으로부터 위, 아래 방향으로 거리가 멀어짐에 따라 임피던스의 영향도 줄어들게 된다 [1]. 이런 임피던스의 영향을 고려하여 이차원 영상화 기술을 개발하게 되며 또한 두 줄 이상의 전극 layer를 사용하여 3차원 영상화도 가능하다.

impedance_contribution
그림 2. 임피던스 영향  [1], 부착된 전극으로부터 위, 아래방향으로 거리가 멀어질수록 임피던싀 영향이 줄어든다.
이 글에서는 16개의 전극을 이용한 이차원 영상을 예로들어 설명하겠다. 영상화하고자 하는 흉부 단면을 \Omega라고 표기하고, 그것의 경계 곡선을 \partial\Omega라 하자. 경계에 16개의 전극 \{\mathcal{E}_1,\mathcal{E}_2,\cdots,\mathcal{E}_{16}\}을 등간격으로 부착한다. 상용화 기기에 탑재되어 있는 인접한 전류 주입 패턴 (Adjacent current injection pattern)으로 임피던스를 측정하자. 즉, 인전한 전극 쌍 (\mathcal{E}_{j},\mathcal{E}_{j+1})으로 전극을 주입하고 나머지 인접한 전극 쌍 (\mathcal{E}_i,\mathcal{E}_{i+1})들에서 전압을 측정한다(그림3 참고). 이때, \mathcal{E}_{16+1}\mathcal{E}_1로 간주한다. 시간 t에서, j번째 전극쌍에서 주입된 전류에 대한, i번째 전극쌍에서 측정된 전압을 V_{j,i}(t)로 표기하자. 각 j (=1,2,\cdots,16)에 대해서, 세 전압 데이타 V_{j,j-1},V_{j,j},V_{j,j+1}는 살갗과 전극의 접촉 임피던스(skin-electrode contact impedance)의 영향을 많이 받아 사용하지 않는다 [2, 3, 4] (여기서, V_{1,1-1}:=V_{1,16} 그리고 V_{16,16+1}:=V_{16,1}). 따라서, 각 주입 전류당 13개의 전압 데이터를 얻고, 16개의 다른 주입 전류들로 부터 총 208=16\times 13의 전압 데이터를 얻어 한장의 이미지 복원을 위해 사용된다. 이 전압들을 모아 다음과 같이 열벡터 \mathbf{V}로 표기하자.

\displaystyle \mathbf{V}=\left[V_{1,3},\cdots,V_{1,15},V_{2,4},\cdots,V_{2,16},\cdots,V_{16,15}\right]^T

사실, Reciprocity law [4]에 의해서 V_{j,i}=V_{i,j}이고, 208의 절반인 104개만이 독립적인 정보를 갖는다.

injection_pattern_cropped
그림 3. 인접한 전류 인가 및 전압 측정 방식 (Adjacent current injection pattern)

Forward Problem

부착된 전극을 통해서, \omega의 각진동수(angular frequency)를 갖는 교류 전류를 I mA만큼 주입한다고 하자. 그리고, 시간 t, 위치 \mathbf{r}\in \Omega에서의 복소수 도전율을 \gamma(\mathbf{r},t)=\sigma(\mathbf{r},t)+\sqrt{-1}\omega\varepsilon(\mathbf{r},t)라고 표기하자. 여기서, \sigma는 도전율(electrical conductivity), \varepsilon는 유전율(electrical permittivity)이다. j번째 전극쌍 (\mathcal{E}_j,\mathcal{E}_{j+1})에서 전류가 주입되면, 인체 내부에 전기장이 생성된다. 이때, 생성되는 전기장은 time-harmonic 전기포텐셜 u_j^\gamma로 표현될 수 있고, 다음의 편미분방정식을 만족한다 [4].

\label{eq:govern} \left\{\begin{array}{rl} \nabla\cdot(\gamma\nabla u_j^\gamma)=0&~~\mbox{in}~~{\Omega}\\ (\gamma\nabla u_j^\gamma)\cdot\mathbf{n}=0&~~\mbox{on}~~\partial{\Omega}\setminus \cup_i^{16}\mathcal{E}_i\\ \int_{\mathcal{E}_i}\gamma \nabla u_j^\gamma\cdot\mathbf{n}=0&~~\mbox{for}~~i\in\{1,2,\ldots,16\}\setminus\{j,j+1\}\\ u_j^\gamma+z_{i}(\gamma\nabla u_j^\gamma\cdot\mathbf{n})&=U_i^j ~~\mbox{on}~~\mathcal{E}_i~~\mbox{for}~~ i=1,2,\ldots,16\\ \int_{\mathcal{E}_j}\gamma\nabla u_j^\gamma\cdot \mathbf{n}\,ds=I&= \int_{\mathcal{E}_{j+1}}\gamma\nabla u_j^\gamma\cdot \mathbf{n} \,ds \end{array}\right.

이때, \mathbf{n}\partial\Omega에 대한 바깥 방향의 단위 수직 벡터이다. z_ii번째에서의 살갗과 전극 접촉 임피던스이다. U_i^jj번째 주입 전류에 대한, 전극 \mathcal{E}_i에서의 전기포테셜 값을 나타낸다. i번째 전극 쌍 (\mathcal{E}_i,\mathcal{E}_{i+1})에서 측정된 전압 V_{j,i}은 다음과 같이 표현된다.

V_{j,i}=U^j_i-U^j_{i+1}

전극 접촉 임퍼던스  z_i 때문에, 엄밀히 말하면  u^\gamma_j|_{\mathcal{E}_i}\neq U^j_i이다. 하지만,

전류를 주입하는데 사용되는 전극을, 전압 측정하는데에 동시에 사용하지 않는다면, z_i(\gamma \nabla u_j^\gamma\cdot\mathbf{n})\approx 0가 되고, u^\gamma_j|_{\mathcal{E}_i}\approx U^j_i임을 가정할 수 있다. 부분적분 또는 Green’s formula에 의해 데이타 V_{j,i}와 도전율 분포 \gamma는 근사적으로 다음의 관계를 갖는다.

V_{i,j}(t)=V^{j,i}(t)=\int_{{\Omega}}\gamma(\mathbf{r},t)\nabla u_i^\gamma(\mathbf{r})\cdot\nabla u_j^\gamma(\mathbf{r})\,d\mathbf{r}

여기서, \mathbf{r}=(x,y)\in\Omega, d\mathbf{r}는 면적원소(area element)를 뜻한다.

심폐영상에서의 역문제(Inverse problem)는 시간차(time-difference) 데이터 \frac{d}{dt}\mathbf{V}로 부터 시간차 영상 \frac{\partial}{\partial t}\gamma를 복원하는 것을 목표로 한다.

Image Reconstruction

가장 널리 쓰이는 영상 복원 방법 중의 하나인 선형 영상복원 방법(the linearized sensitivity method)에 대해 설명하겠다. \frac{d}{dt}\mathbf{V}\frac{\partial}{\partial t}\gamma의 관계를 조사하기 위해, \eref{eq:eq}의 양변을 시간에 대해 미분하면 다음식을 얻는다.

\label{eq:nonlinear} \dot{V}_{i,j}(t)=-\int_{{\Omega}}\dot{\gamma}(\mathbf{r},t)\nabla u_i^\gamma(\mathbf{r})\cdot\nabla u_j^\gamma(\mathbf{r})\,d\mathbf{r}

이때, \dot{V}_{i,j}\dot{\gamma}은 각각 V_{i,j}\gamma의 시간 $t$에 대한 미분을 뜻한다. \eref{eq:nonlinear}에서 u_j^\gamma\gamma에 의존하기 때문에, 데이터 \dot{V}_{i,j}\dot{\gamma}에 비선형적으로 의존한다. u_j^\gamma를 다른 포텐셜 함수로 대체 함으로써, 다음과 같이 선형화 시킨다.

\label{eq:linear} \dot{V}_{i,j}(t)\approx-\int_{{\Omega}}\dot{\gamma}(\mathbf{r},t)\nabla u_i(\mathbf{r})\cdot\nabla u_j(\mathbf{r})\,d\mathbf{r}

여기서, u_jj번째 전류에 대응되는 포텐셜로 다음의 편미분 방정식의 해이다.

\label{eq:shunt} \left\{\begin{array}{rl} \nabla\cdot\nabla u_j=0&~~\mbox{in}~~{\Omega}\\ \nabla u_j\cdot\mathbf{n}=0&~~\mbox{on}~~\partial{\Omega}\setminus \cup_i^{16}\mathcal{E}_i\\ \int_{\mathcal{E}_i} \nabla u_j\cdot\mathbf{n}=0&~~\mbox{for}~~i\in\{1,\ldots,16\}\setminus\{j,j+1\}\\ \mathbf{n}\times\nabla u_j=0& ~~\mbox{on}~~\mathcal{E}_i~~\mbox{for}~~ i=1,\ldots,16\\ \int_{\mathcal{E}_j}\nabla u_j\cdot \mathbf{n}\,ds=I&= \int_{\mathcal{E}_{j+1}}\nabla u_j\cdot \mathbf{n} \,ds\\ \sum_{i=1}^{16} u_j|_{\mathcal{E}_i} =0 & \end{array}\right.

선형식 방정식 \eref{eq:linear}를 단순하게 표현하면,

\label{eq:linear_vector} \dot{\mathbf{V}}(t)\approx\int_{\Omega}\dot{\gamma}(\mathbf{r},t)\mathbf{s}(\mathbf{r})\,d\mathbf{r},

여기서, \dot{\mathbf{V}}(t)는 시간 t에서 전압 데이터에 대한 시간차 벡터이다, 그리고 \mathbf{s}(\mathbf{r})는 다음과 같이 주어진다.

\label{sense} \hspace{-0.0cm}\mathbf{s}(\mathbf{r}):=\left[ s^{1,3}(\mathbf{r}),s^{1,4}(\mathbf{r}),\ldots, s^{16,14}(\mathbf{r})\right]^T~~~\mbox{with}~~~s^{i,j}(\mathbf{r})=-\nabla u_i(\mathbf{r})\cdot\nabla u_j(\mathbf{r}).

선형 방정식 \eref{eq:linear_vector}을 수치적으로 풀기 위해, 도메인 \Omega을 유한개로 이산화 시킨다 (\Delta_k, k=1,2,\cdots, n_{elem}):

\Omega=\bigcup_{k=1}^{n_{elem}} \Delta_k

이산화된 \dot{\gamma}(t)를 다음의 벡터 \dot{\boldsymbol{\gamma}}(t)로 표현하면, 선형방정식은 다음 행렬 방정식으로 표현된다:

\label{eq:lienarized}\dot{\mathbf{V}}(t)~\approx~\mathbb{S}\, \dot{\boldsymbol{\gamma}}(t)

여기서, \mathbb{S}는 센시티비티 행렬(sensitivity matrix) 또는 야코비안 행렬(Jacobian matrix)로 불리우고 다음과 같이 주어진다.

\label{sense_discrete} \mathbb{S}=\left[ \begin{array}{ccc} &| & \\ \cdots&\mathbf{S}_k&\cdots\\ &| & \end{array} \right]\quad\mbox{with}~~~ \mathbf{S}_k:=\int_{\Delta_k} \mathbf{s}(\mathbf{r})\,d\mathbf{r}

센시티비티 행렬은 영상 도메인 \Omega과 전극 위치, 전류 인가 패턴에 의존하며, 행과 열의 갯수는 각각 208 과 n_{elem}이다. 센시티비티의 k번째 열 \mathbf{S}_kk번째 픽셀 \Delta_k에 대응하여, 그 픽셀에서 \gamma가 단위만큼 변할때의 측정된 경계 전압의 변화를 뜻한다.

선형 시스템 \eref{eq:lienarized}을 풀면 \dot{\boldsymbol{\gamma}}를 복원할 수 있지만, 센시티비티 행렬 \mathbb{S}가 불량하기 때문에, 정규화(regularization) 과정이 필요하다. 다음의 정규화된 최소제곱법(the regularized least-squares data-fitting method)을 흔히 사용한다:

\label{eq:min} \left( \mathbb{S}^T\mathbb{S}+\lambda \mathcal{R}^T\mathcal{R}\right)^{-1} \mathbb{S}^T \dot{\mathbf{V}}(t)

여기서, \mathcal{R}은 정규화 오퍼레이터로, 단위 행렬(\mathcal{R}=\mathbb{I})인 경우 Tikhonov 정규화를 뜻한다. \lambda는 정규화 파라미터로 정규화 정도를 결정한다. \lambda가 클 경우 데이터 오차에대해 안정적이지만, 너무 클 경우 영상이 희미해지고 정확도가 떨어진다.

다음은 실제 사람의 흉부에 측정된 데이터에 영상 복원 알고리즘을 적용하여 얻어진 심폐영상들이다. 여기서, 측정 데이터는 경희대의 IIRC그룹 [5] 에서 제작된 기기를 사용하였고, 영상 복원 방법은 Fidelity-embedded regularization (FER) method [6]이 사용되었다.

실시간 EIT 실험 동영상을 보려면 다음 링크를 방문하면 된다. Movie1, Movie2

exp_sciospec
그림 4. EIT를 이용한 실시간 심폐기능 영상
[1] E. Teschner, I. Michael, and S. Leonhardt, “Electrical Impedance Tomography: The Realisation of Regional Ventilation Monitoring (2nd Edition)”, Tech. rep. Draeger GmbH, 2015.
[2] V. Kolehmainen, M. Vauhkonen, P. A. Karjalainen, J. P. Kaipio, “Assessment of errors in static electrical impedance tomography with adjacent and trigonometric current patterns”, Physiol. Meas., vol. 18, no. 4, pp. 289-303, 1997.
[3] W. R. B. Lionheart, “EIT reconstruction algorithms: pitfalls, challenges and recent developments”, Physiol. Meas., vol. 25, no. 1, pp. 125-142, 2004.
[4] J. K. Seo and E. J. Woo, Nonlinear Inverse Problems in Imaging, Wiley, 2013.
[5] 경희대학교 임피던스 영상신기술연구센터 Impedance Imaging Research Center(IIRC), http://iirc.khu.ac.kr.
[6] K. Lee, E. J. Woo and J. K. Seo, “A Fidelity-embedded Regularization Method for Robust Electrical Impedance Tomography”, arXiv:1703.01821 [math.AP], 2017.

Machine-learning-based nonlinear decomposition of CT images for metal artifact reduction

CT_metal_artifact
Metal-Artifacted CT Image

노년층 인구가 증가함에 따라, 인체 내 인공 보철물 소유자가 증가하고 있다. CT영상에서 인체 내 금속성 보철물이 포함되면 metal artifacts가 생성되어 영상을 왜곡시킨다. 이런 영상 왜곡의 근본적인 이유는 금속 물질 때문에 CT sinogram이 아래 그림과 같이 라돈 변화의 Range에 속하지 않기 때문이다. (즉, X-ray 데이터가 비선형적으로 변화되기 때문이다.) 이 현상을 보완하기 위해, Convolution Neural Network (CNN) 기반의 machine learning 기술을 이용 하고자 한다.

radon_range_space
Beam-Hardening Structure
goal
Our Goal

우리는 CNN중 U-net을 사용하여 CT 영상에서 생성되는 metal-artifact를 학습시켰다.  처음에 아래 그림과 같은 dataset 구성으로 학습되기를 시도했다. 그리고  이 학습된 네트워크를 통해, 생성된 artifact를 원래 영상에서 빼줌으로써 artifact가 없는 CT 영상을 얻고자 했다. 하지만 이 dataset (training & label)은 임의의 background (금속 밖의 영역)에 대해 학습이 잘 되어야 하기때문에 네트워크가 잘 수렴하지 않았다. 게다가 실제로 얻기 힘든 dataset 이므로(ground-truth image를 얻기 힘들기 때문에),  좋은 dataset구성이 아니었다.

dl_mar_try1
Metal-Artifact Decomposition using U-net 1

이 문제를 해결하기위해,  ex-vivo metal (생체 밖의 금속)에 U-net을 적용하였다. 이 방식은 복잡한 background에 대한 영향을 줄일수 있을 뿐만 아니라, 금속 데이터만을 필요로 하기 때문에, (즉, 금속을 포함한 환자의 CT 영상과 ground-truth CT영상이 필요 없다.) data 습득면에서도 U-net 1 보다 상대적으로 큰 장점을 가진다. 학습을 위해, 금속에 의해 생성되는 대략적인 numerical metal artifact를 형성하여,  input data로 사용하고, 그 금속이 생성하는 실제 artifact CT image를 output data로 사용한다.  이 학습된 ntwork를 이용하여, 우리는 다음과 같은 작업을 진행한다. 먼저 metal을 포함하는 실제 CT 영상에서 metal을 분할하고 metal artifact의 대략적인 형태를 만든 뒤에 학습된 네트워크를 통해 실제 artifact를 생성한다 (아래 그림 참조). 그리고 그 artifact를 원래 영상에서 빼주어 artifact가 없는 CT영상을 얻었다.

dl_mar_try2
Metal-Artifact Decomposition using U-net 2

아래 링크를 통해, 세부적인 내용과 실험결과를 확인 할 수 있다.

관련 논문 : Machine-learning-based nonlinear decomposition of CT images for metal artifact reduction

Ubuntu 16.04 LTS에 tensorflow gpu버전 설치하기

ubuntu 16.04 LTS

Anaconda를 이용하여 Tensorflow 를 설치하는법

GPU setting

1.  NVIDIA-Driver

Screenshot_from_2017-01-16_20-47-13

All Settings – Software & Updates – Additional Drivers – Unsing NVIDIA binary driver 를 선택후 Apply Changes를 누른다.

(드라이버 업데이트를 하는동안 튕기면1.  NVIDIA-Driver, 무한로그인에 걸리게 되므로 다른작업을 하지않는것을 권장함.)

그러면 엔비디아 드라이버가 설치된다.

재부팅후 터미널에 nvidia-smi를 치면 그래픽카드 정보가 출력된다.

Screenshot_from_2017-01-16_22-21-15.jpg

 

2. cuda toolkit

전반적인 설치는 텐서플로우 홈페이지를 참조하였다.

https://www.tensorflow.org/get_started/os_setup

먼저 엔비디아 사이트에서 cuda toolkit을 받는다.

https://developer.nvidia.com/cuda-downloads

Screenshot_from_2017-01-16_20-59-46.jpg

Linux – x86_64 – Ubuntu – 16.04 – deb(network)를 선택후 Download를 누른다.

그다음 터미널을 열어준다.

cd Downloads/ (파일을 받은곳으로 이동)

sudo dpkg -i cuda-repo-ubuntu1604_8.0.44-1_amd64.deb

sudo apt-get update

sudo apt-get install cuda

를 입력하면 cuda가 설치된다.

bashrc file 수정

cd /home/username/

gedit .bashrc

그리고

export CUDA_HOME=/usr/local/cuda-8.0

export PATH=/usr/local/cuda-8.0/bin${PATH:+:${PATH}}

export LD_LIBRARY_PATH=/usr/local/cuda-8.0/lib64${LD_LIBRARY_PATH:+:${LD_LIBRARY_PATH}}

이를 맨 아래 추가후 저장하고 종료한다.

그리고 source ~/.bashrc를 입력

Screenshot_from_2017-01-16_21-16-22.jpg

터미널에 nvcc –version를 입력하면 위와같이 cuda가 설치된것을 확인할수 있다.

3. Install cuDNN 5.1

https://developer.nvidia.com/cudnn

Screenshot_from_2017-01-16_21-07-41.jpg

링크에 들어가서 Downloads를 누르면 로그인 화면이 나온다.

이부분은 NVIDIA에 회원가입을 해야 할수있기때문에, 회원가입후 로그인을 하면 다운로드 화면이 나오게된다.

cuDNN v5.1을 선택하고, 그 밑으로 나오는 cuDNN v5.1 Library for Linux를 다운로드 한다.

그다음 터미널을 열어서

cd Downloads/

tar xvzf cudnn-8.0-linux-x64-v5.1.tgz (또는 직접 압축풀기 하여도됨)

sudo mv cuda/include/cudnn.h /usr/local/cuda-8.0/include

sudo mv cuda/lib64/libcudnn* /usr/local/cuda-8.0/lib64

sudo chmod a+r /usr/local/cuda-8.0/include/cudnn.h /usr/local/cuda-8.0/lib64/libcudnn*

위의 명령어를 입력하면

다운받은 cuda폴더에 있던 내용물이 /usr/local/cuda-8.0폴더에 복사됨.

본격적으로 tensorflow를 설치하기전에

pip을 이용해서 tensorflow를 설치할것이기 때문에 pip부터 설치해준다.

pip설치

sudo apt-get install python-pip

sudo apt-get install python3-pip

pip install –upgrade pip

pip3 install –upgrade pip을 이용해서 업데이트

(그외 텐서플로우에 필요한 다른 라이브러리들)

sudo apt-get install python-numpy python-dev python-wheel

sudo apt-get instal  python3-numpy python3-dev python3-wheel

sudo apt-get install python-matplotlib
sudo apt-get install python3-matplotlib

4. Install Anaconda

https://www.tensorflow.org/get_started/os_setup#anaconda_installation

Tensorflow 홈페이지에서 Anaconda installation항목을 참조하였다.

https://www.continuum.io/downloads

위 링크를 통하여 아나콘다를 다운받는 홈페이지로 이동한다.

Screenshot_from_2017-01-16_13-06-22.jpg

자신이 사용하는 python version에 맞게 아나콘다 설치파일을 받는다.

그다음 터미널을 열어서

cd Downloads/

bash Anaconda3-4.2.0-Linux-x86_64.sh

를 입력하면 아나콘다가 설치된다.(파이썬 3.5의 경우)

설치 마지막에

Screenshot_from_2017-01-16_13-15-28.jpg

bashrc파일에 PATH를 추가할꺼냐고 물어보는데 yes를 입력한다.

그리고 설치가 완료되면 source .bashrc를 입력해주고

터미널에 python을 입력하면

Screenshot_from_2017-01-16_13-16-49.jpg

이와같이 Anaconda 환경에서 python이 실행되는것을 확인할수 있다.

그리고 conda update conda를 입력하여 아나콘다를 업데이트 해준다.

 참고)————–———————–———————–———————–———————–———————-

bashrc에 PATH가 추가가 안된경우(위와같이 python을 입력했는데 Anaconda가 안나올때)

파이썬과 아나콘다를 연결해주기 위해서 bashrc파일을 수정해줘야한다.

bash파일에 추가하기 위해서

cd를 입력해서 홈폴더를 들어간후

gedit .bashrc를 입력후

맨 아랫줄에

export PATH=/home/kangcheol/anaconda3/bin:$PATH

를 추가하고 save후 끈다음 그상태에서

source .bashrc를 입력해서 바뀐것을 적용

Screenshot_from_2017-01-16_13-56-13.jpg

—————————————————————————————————————————————–

 

5. Create Tensorflow Enviroment

터미널에

$ conda create -n tensorflow python=3.5

를 입력해서 콘다 환경을 만든다

Screenshot_from_2017-01-16_13-20-12.jpg

설치가 완료되면 마지막에 이와같이 뜬다.

Screenshot_from_2017-01-16_13-20-12 (1).jpg

source activate tensorflow를 입력하여 환경을 활성화 할수있다.

6. Install Tensorflow

tensorflow 가 activate된 상태에서

1) 온라인에서 받으면서 설치

export TF_BINARY_URL=https://storage.googleapis.com/tensorflow/linux/gpu/tensorflow_gpu-0.12.1-cp35-cp35m-linux_x86_64.whl

를 입력하고

그후에

pip3 install –upgrade $TF_BINARY_URL

or

pip3 install –ignore-installed –upgrade $TF_BINARY_URL

를 입력하여 텐서플로우를 설치

2) 바이너리 파일을 받은후 설치 (좀더 안정적)
위 링크로 들어가면 whl파일을 받을수있다.
그리고 터미널에 들어가서 아래와같은 명령어를 입력한다.
cd Downloads/
pip install tensorflow_gpu-0.12.1-cp35-cp35m-linux_x86_64.whl
pip3 install tensorflow_gpu-0.12.1-cp35-cp35m-linux_x86_64.whl

마지막으로 tensorflow가 잘 작동하는지 확인할수 있다.

Screenshot_from_2017-01-16_21-58-01