의학물리 : 제 3 권제 호 0 음향광학단층촬영 (Acousto-Optical Tomography) 의수학적모델과수치해석적시뮬레이션 이화여자대학교의과대학 * 의과학연구소, 방사선종양학과 남혜원 * ㆍ허장용 * ㆍ김소영 * ㆍ이레나 본연구는최근의공학분야에서중요한영역으로대두되고있는광학과초음파의장점을결합한영상방법인 AOT (Acousto-Optical Tomography) 의수학적모델을제시하였다. AOT는광학필드를초음파기둥에서변화시켜서초음파기둥의위치정보를이용하여영상을재구성하는방법이다. AOT의수학적모델은두단계로나뉠수있다. 첫번째단계에서는광학필드의복원을하고, 두번째단계에서는획득한광학필드를기반으로확산방정식의역문제를풀어흡수함수μ (absorption coefficient) 를산출한다. 본연구에서는두번째단계에해당하는역문제의해를구하기위하여수치해석적인최소화문제로변환하고, 수치적팬텀을이용하여시뮬레이션하였다. 전통적인기울기하강방법을이용하여역문제시뮬레이션의결과를보였다. 전변동정규화기반의최소화문제를제안하여기울기하강방법의결과에서보인번짐효과를개선하였다. 중심단어 : 음향광학, 확산방정식 서 의료영상분야는 9세기후반의 x-ray를시작으로하여많은발전이이루어지고있다. 현재는 CT, PET, MRI, 초음파영상등여러가지다른의료영상방법이대두되고있다. 초음파영상과광학영상은인체에무해하고비전리, 그리고다른의료영상에비해비용효과가크다는장점을가지고있다. 그러나초음파영상은각조직마다낮은대조도때문에의료진단에사용하기힘든단점을가지고있다. 반면광학영상은높은대조도를갖지만낮은해상도와깊이제한에서오는단점때문에의료진단에는한계를가지고있다. 따라서최근들어새로운영상방법으로초음파와광학을결합하는영상분야가연구되고있다. 그중대표적인두가지는광음향이라고불리는 PAT (photo-acoustic tomography) 와 UMOT (ultrasound-modulated optical tomography) 라고도불리는 AOT (acousto-optical tomography) 이 본연구는보건복지부보건의료연구개발사업의지원에의하여이루어진것임 (A098). 이논문은 0 년 월 3 일접수, 0 년 3 월 5 일채택되었음. 책임저자 : 이레나, (58-70) 서울시양천구목동 9- 이화여자대학교의과대학부속목동병원방사선종양학과 Tel: 0)650-5337, Fax: 0)654-0363 E-mail: renalee@ewha.ac.kr 론 다. ) PAT와 AOT의개발동기는 Table 을보면잘설명될수있다. ) 순수한광학영상은생물학적조직 (biological tissue) 에서강한빛의산란으로선명하지않은영상을보여주고있다. 광학영상은 OCT (optical coherence tomography) 와같이제한된깊이를갖거나, 또는 DOT (diffuse optical tomography) 와같이제한된해상도를갖게된다. 광전자와비교할때, 초음파는생물학적조직에서 3배정도낮은산란을갖는다. 그결과, 초음파영상은높은해상도를갖는다. 그러나, 초음파영상은생물학적조직의특성을검출을기반으로영상이이루어지므로질병의조기관측에는어려움을갖는다즉, 산소포화도, 헤모글로빈의밀집등은감지하기가힘들다. AOT에서영상의대비는광학의성질에기반을하고, 높은해상도는광학의산란이나발산에영향을받지않는다. ) AOT의실험장치는 Fig. a에나타나있다. ) Fig. b에서보이는것과같이 AOT에서는레이저빔이생물학적조직을통과할때일부의산란된광입자는초음파기둥에의해초음파파장으로인해파장이변한다. 이변화한파장은초음파주파수를가지므로초음파기둥을지나는광입자를변환광입자 (tagged photon) 라고부른다. 초음파주파를갖는변환광입자와초음파기둥의위치정보를기반으로영상을재구성한다. 본논문에서는 AOT의수학적모델링을세우고그것의 - 4 -
의학물리 : 제 3 권제 호 0 Table. Comparison table: properties of ultrasound and optical tomography. Properties Modalities OCT DOT US (Ultrasound) AOT/PAT Contrast Resolution Imaging depth Speckle artifacts Scattering coefficient Good Excellent ( 0μm) Poor ( mm) Excellent Poor ( 5 mm) Excellent ( 5 cm) None Poor for early cancers Excellent & scalable ( 50μm) Good & scalable ( 3 cm) Weak Excellent (=DOT) Excellent (=US) Good & scalable (=US) None Motivation for acousto-optical tomography (AOT). For comparison purposes, OCT (Optical coherence tomography) and DOT (diffuse optical tomography) as the pure optical imaging modalities and US (ultrasonography) at 6 MHz as the pure ultrasound imaging modality. Fig.. (a) AOT Experimental setup. (b) Sketch of AOT. 수치해석적인시뮬레이션을하고자한다. 재료및방법. 배경지식 ) 확산방정식 : 생물학적조직에서입자의이동은 Boltzmann 수송방정식으로나타낼수있다. 3,4) 만일흡수계수 (absorption coefficient) μ a 와산란계수 (scattering coefficient) μ s 가 μ a/μ s 의관계를만족하면 Boltzmann 방정식은아래식 () 과같은확산방정식으로근사시킬수있다. 5,6) ( + μaϕ( = S( 3( μ s + μa ) () 확산방정식 () 은광학영상에서역수송문제들을해결하는데많이사용되고있다. ) 초음파기둥과광입자와의상호작용 : 초음파와광학영상의결합을알려면생체조직내에서광학필드와초음파빔의상호작용을알아야한다. 그것은다음과같은상관작용함수로나타난다. ) k G( τ ) 3 3l I( ) P( ) A ( )( cosϖ sτ ) exp( iδφ' ( τ ) d () - 43 -
남혜원외 3 인 : 제목음향광학단층촬영 (Acousto-Optical Tomography) 의수학적모델과수치해석적시뮬레이션 n( α ) ' exp( i Δφ' ( τ ) = < exp( iq j Δr j ) > (3) j= 위방정식에서 r은위치를나타내는벡터이고, k는광파동벡터, ωs는주파수, A( 은이동의크기를의미한다. 거리값 A( 은초음파압력에따라달라진다. 또한 I( 은광학부분의크기로, 앞의식 () 에서 φ( 과같은값이다. 상 ' 수 l은평균자유행정이동의평균거리값이고, Δ r j ) 는 j 번째산란에서변화량, qj는산란벡터, α는관측시간 t와 t+τ 사이의궤도를의미한다. 기호 <, >는시간으로의평균을의미한다. (3) 번식은초음파의부재중관측될수있는값이다. 시간 τ로평균을하면, 수식 () 는상수 c를포함하는아래식 (4) 와같이된다. G( = c P( ) A ( ) ϕ ( ) d (4) 함수 G 은위치 r에서의시간으로평균한상관작용함수이다. E를레이저관측하는 CCD 카메라의표면이라고할때, 초음파기둥의평행이동에따라변화하지않으므로임의의 x에서의측정치는식 (5) 로나타낼수있다. G( x) = G ( dr = c P( ) A ( x) ϕ( ) d dr E E = c P( ) A ( x) ϕ( ) d. AOT 의수학적모델 AOT의수학적모델은두단계로나뉠수있다. 첫번째는측정된데이터 G(x) 로부터광학필드 φ를재구성하는것이고 ( 식 (5)), 두번째단계에는 단계에서복원한 φ로부터흡수율함수를복원하는것이다 ( 식 ()). ) 광학필드 φ 복원 : 식 (5) 에서와같이 G(x) 가주어졌을때함수 φ를복원한다. 함수 T를다음과같은회선연산함수로정의한다. (5) 상수 c는위의식 (5) 에서와같이주어져있다. 주어진측정치 g o(x) 로부터광학필드 φ를복원하는것이첫번째단계이다. 이동거리를나타내는함수 A는초음파파장방정식 (acoustic wave equation) 의해다. 즉, g는연속함수로써 L () 에속한함수이다. 그러면함수 T: L () L () 는선형유계컴팩트연산자이다. 또한 T의수반연산자 T* 가존재한다. 위의식 (6) 을바로푸는대신정규화를포함한아래와같은정규연산방정식을풀도록한다. ( T * T + εi ) u = T * go 관측된측정치 g o(x) 로구해진광학필드는다음과같다. ( T * T + I ) T * g P = ε ϕ ) 흡수함수의복원 : 첫번째단계에서구한광학함수 ϕ 와확산방정식 () 을이용하여흡수계수함수 μ a 를구한다. 일반적으로건강한조직과그렇지않은조직은광학흡수력이다르다. 이것을이용하여영상을복원하고자한다. 방정식 () 에서 ϕ 가주어졌을때 μ a 를구하는것은역문제 (inverse problem) 로구성할수있다. 편미분방정식의해를구하기위하여경계조건을정의해야한다. 영상을구성하고자하는도메인을 로놓고, 그경계치는 Fig. 에서와같이빛이들어가는영역 Γ 0 와빛이나오는부분인 Γ 로나뉠수있다. 흡수계수함수를구하는역문제는다음과같이정의될수있다. 함수 γ(x) 를다음과같이정의하자. = 3( μ + μ s o ) Tu( x) = = g( x y) u( y) dy ( g u) = g o ( x) ( x) (6) 여기서함수 u( x) = P( y) ϕ( y) 이고 g( z) = ca ( z), 함수 A와 Fig.. Sketch of the boundary. - 44 -
의학물리 : 제 3 권제 호 0 경계조건을포함한확산방정식은 γ + μϕ = 0 in ϕ = f on Γ (7) 0 n ϕ + αϕ = 0 on Γ n 경계조건에나타난함수 n 는입력유량을의미한다. 경계 Γ0에서의유량은함수 f로주어진다. 여기서함수 f는레이저소스에의해결정된다. 레이저빛이나가는경계인 Γ에서산출유량은입력유량에비례한다. 위의역문제를최소화문제로접근한다. 주어진광학필드 φ와경계치방정식 (boundary condition) 을통하여흡수계수함수 μ를구하고자한다. 이것을최소화문제 (minimization problem) 로나타내보면아래식과같다. ϕ min J = F( μ) ϕ μ X L ( ) (8) 여기서함수 F는 X에서 Sobolev 공간 H () 로가는함수로써 F( μ ) = ϕ 를만족한다. 집합 X는 L () 의부분집합으로아래와같이정의하였다. X = { μ L ( ) : 0 < μ μ( x) μ <, μ is Lipschitz continuous with Lip c} 결과본연구에서최소화문제 (8) 을기울기하강 (gradient descent) 방법을이용하여수치적으로시뮬레이션하였다. 첫번째시뮬레이션에서는정규화없이기울기하강방법을사용하였고두번째방법에서는전변동정규화 (total variation regularization) 방법을첨가하여시뮬레이션하였다.. 기울기하강방법 (Gradient Descent Method) 수치해석적으로해를구하는가장일반적인방법인기울기하강방법을적용하기위하여최소화문제 J(μ) 의도함수를구해야한다. L norm은 H 상에서미분가능하므로함수 J(μ) 는 Fréchet 미분가능하다. Fréchet 도함수를구하여보면, * DJ ( μ ) = DF ( μ ) ( F ( μ ) ϕ ) 로나타내어진다. 여기서 * 는수반연산자를의미한다. DF( μ ) = ψ 라고놓으면 ψ 는아래의방정식의해로나타난다. γ ψ + μψ = sϕ + sd in ψ = f on Γ0 n ψ + αψ = 0 on Γ n 함수 g = ϕ ϕ 라고놓고, 함수 u를다음편미분방정식의해라고정의하면, γ u + μu = g in u = f on Γ0 n u + αu = 0 on Γ n J(μ) 의 Fréchet 도함수는 Fig. 3. Result of the gradient descent method. (a) Actual solution. (b) Reconstruction result of (a) with gradient descent method. - 45 -
남혜원외 3 인 : 제목음향광학단층촬영 (Acousto-Optical Tomography) 의수학적모델과수치해석적시뮬레이션 DJ ( δ μ ) = ψg =< δ, ( uϕ + Dγ u) > μ 가되며, 다시정리하자면 DJ = ( uϕ + Dγ u) 로나타내어진다. 참 고문헌 7의정리 8.8에의해, 함수 u와 φ는 W, () 에속한함수이고, 그러므로함수 u, φ, u, φ는 W, () 에속한다. 공간 W, () 로부터 L 4 () 로의컴팩트임베딩이존재한다는정리 (Theorem 7.6 5) 에따라 J의 Fréchet 도함수 DJ(μ) 는 L 4 () 에속한함수가된다. Fig. 3b는팬텀흡수함수 μ를 Fig. 3a와같이놓고기울기하강을이용하여복원한흡수함수 μ이다. 팬텀에서의최대은 로놓았는데시뮬레이션을통하여구해진최대치는.6 정도이고, 번짐으로인하여직사각형인팬텀이둥근형태로복원되는결과를볼수있다.. 전변동최소정규화 Fig. 3b에서보여진번짐효과를없애는이미지제고 (image enhancement) 기술중 Dobson과 Santosa 6,7) 가제안한비선형정규화를적용하는방법이있다. 이기술은전변동을최소화하는정규화를넣음으로인해좀더번짐을없애면서복원을가능하게한다. 이기술을적용하면, 위의식 (8) 은다음과같이바뀐다. min J ( μ X μ) = F( μ) ϕ + β μ L ( ) (9) Fig. 4. Reconstruction result of Fig. 3(a) with gradient descent method with total variation regularization. 였다. AOT는초음파와광학영상을결합한새로운영상단층촬영기법으로초음파의단점인낮은대조와광학영상의단점인낮은해상도와낮은깊이를보안한방법이다. 본연구에서는간단한수학적모델과그것의수치해석적방법을통하여영상복원방법을제시하였다. 좀더정확한시뮬레이션을위하여향후리서치는이영상방법이잡음에견고한지여부와현재 차원에서의시뮬레이션에서 3차원시뮬레이션으로확장할예정이다. 식 (9) 에서의 μ는분포개념에서의경사를의미한다. 즉, μ는벡터값을갖는라돈척도 (Radon measure) 이고식 (9) 에서의 μ 는 μ의전변동을의미한다. 양의상수 β의크기로부터정규화의정도를결정할수있다. Fig. 3a 의팬텀으로전변동정규화를포함하여시뮬레이션한결과는 Fig. 4와같다. Fig. 3b의결과와비교해보면팬텀의네모난모양은찾았으나최대치가.로써전병동정규화를넣지않았을경우보다 0% 더낮아졌고, 실제팬텀의크기보다두배이상더크게복원되는결과를얻었다. 고찰및결론 본연구에서는초음파와레이저를결합한새로운의학영상방법에대한수학적모델과그것의수치해석적시뮬레이션을하였다. 수치적시뮬레이션은일반적인방법인기울기하강방법과전변동정규화를첨가한방법을사용하 참고문헌. Wang LV: Ultrasound-mediated biophotonic imaging: A review of acousto-optical tomography and photo-acoustic tomography. Disease Markers 9:3-38 (003, 004). Kempe M, Larionov M, Zaslavsky D, Genack AZ: Acoustooptic tomography with multiply scattered light. J Opt Soc Am A 4:5-58 (997) 3. Dorn O: A transport-back transport method for optical tomography. Inverse Problems 4:07-30 (998) 4. Dorn O: Scattering and absorption transport sensitivity functions for optical tomography. Optics Express 7:49-506 (000) 5. Gilbarg D, Trudinger NS: Elliptic Partial Differential Equation of Second Orde Springer-Verlag, Berlin, Heidelberg (998) 6. Dobson DC: Recovery of blocky images in electrical impedance tomography. Inverse Problems in Medical Imaging and Nondestructive Testing, Springer-Verlag, (997), pp. 43-64 7. Dobson DC, Santosa F: An image-enhancement technique for electrical impedance tomography. Inverse Problems 0: 37-334 (994) - 46 -
Mathematical Model for Acousto-Optical Tomography and Its Numerical Simulation Haewon Nam*, Jangyong Hur*, Soyoung Kim*, Rena Lee 의학물리 : 제 3 권제 호 0 Departments of *Medical Science, Radiation Oncology, School of Medicine, Ewha Womans University, Seoul, Korea In this pape Acousto-Optical tomography is modeled by a linear integral equation and an inverse problem involving a diffusion equation in n-spatial dimensions. We make two-step mathematical model. First, we solve a linear integral equation. Assuming the optical energy fluence rate has been recovered from the previous equation, the absorption coefficient μ is then reconstructed by solving an inverse problem. Numerical experiments are presented for the case n=. The traditional gradient descent method is used for the numerical simulations. The result of the gradient descent method produces the blurring effect. To get rid of the blurring effect, we suggest the total variation regularization for the minimization problem. Key Words: Acousto-optical tomography, Diffusion equation - 47 -