광자 수송을 위한 몬테카를로 방법

몬테카를로 방법으로 광자 전파를 모델링하는 것은 광자 전송을 시뮬레이션하기 위한 유연하지만 엄격한 접근법이다. 이 방법에서 광자 운송의 국지적 규칙은 광자간 상호작용 현장 사이의 광자 이동의 단계적 크기와 산란 사건이 발생할 때 광자의 궤적에서 편향의 각도를 설명하는 확률 분포로 표현된다. 이는 미분방정식을 이용한 광자의 움직임을 설명하는 복사전달식(RTE)에 의한 광자수송을 분석적으로 모델링하는 것과 같다. 그러나, RTE의 폐쇄형 형태 솔루션은 종종 가능하지 않다. 일부 기하형상의 경우, 확산 근사치를 사용하여 RTE를 단순화할 수 있다. 비록 이것이 많은 부정확성, 특히 선원과 경계 근처를 야기하지만 말이다. 이와는 대조적으로 몬테카를로 시뮬레이션은 추적된 광자의 수를 증가시킴으로써 임의로 정확해질 수 있다. 예를 들어, 초기 탄도 광자 흐름과 후기 확산 전파를 모두 반무한 매체에서 발생한 연필 빔 사건을 몬테카를로 시뮬레이션한 영화를 보라.

몬테카를로 방법은 반드시 통계적 방법이기 때문에 정밀도를 얻기 위해서는 상당한 계산 시간이 필요하다. 또한 몬테카를로 시뮬레이션은 원하는 공간 및 시간 분해능을 통해 동시에 여러 물리적 양을 추적할 수 있다. 이러한 유연성은 몬테카를로 모델링을 강력한 도구로 만든다. 따라서 계산상 비효율적이지만, 몬테카를로 방법은 많은 생물의학 용도에 대한 광자 전송의 시뮬레이션 측정을 위한 표준으로 간주되는 경우가 많다.

몬테카를로 방법의 생물학적 응용

바이오메디컬 영상화

생물 조직의 광학적 특성은 생물의학 이미징에 대한 접근법을 제공한다. 혈액과 멜라닌에서 흡수되고 신경세포와 암세포핵에서 산란되는 등 내생적인 대조가 많다. 또한 형광 탐침은 많은 다른 조직을 대상으로 할 수 있다. 현미경 검사 기법(콘코칼, 투포톤, 광학 일관 단층촬영법 포함)은 이러한 성질을 높은 공간 분해능으로 영상화할 수 있는 능력이 있지만, 탄도 광자에 의존하기 때문에 깊이 침투가 몇 밀리미터로 제한된다. 광자가 산재해 있는 조직으로 더 깊이 영상화하려면 그러한 환경에서 광자의 많은 수의 통계적 행동에 대한 더 깊은 이해가 필요하다. 몬테카를로 방법은 조직 깊은 곳의 광학적 특성을 재구성하기 위해 다른 기법에 의해 사용되어 온 유연한 프레임워크를 제공한다. 이 중 몇 가지 기술에 대한 간략한 소개가 여기에 제시되어 있다.

• 광음향 단층 촬영 PAT에서는 확산 레이저 광선이 흡수되어 국소 온도 상승을 일으킨다. 이 국소 온도 변화는 초음파 변환기를 통해 감지되는 열가소성 확장을 통해 초음파를 생성한다. 실제로 다양한 설정 매개변수(즉, 광파장, 변환기 숫자 구멍)가 다양하며, 그 결과 몬테카를로 모델링은 실험 방법에 앞서 조직 반응을 예측하는 데 귀중한 도구가 된다.
• 확산광학단층촬영(DOT)은 근적외선 광원과 검출기를 배열해 생물조직의 광학적 특성을 측정하는 영상기술이다. 옥시와 디옥시헤모글로빈(기능성 신경영상 또는 암 검출용)에 의한 흡수, 형광 탐침의 농도 등 다양한 대비가 가능하다. 영상을 재구성하기 위해서는 빛이 주어진 소스에서 주어진 검출기로 이동하는 방식과 측정이 광학적 특성(전방 모델이라고 함)의 분포와 변화에 따라 어떻게 달라지는지 알아야 한다. 생물조직의 산란성이 높아 그러한 경로가 복잡하고 민감도 기능이 분산된다. 전방 모델은 몬테카를로 방법을 사용하여 종종 생성된다.

방사선요법

방사선 치료의 목표는 일반적으로 이온화 방사 형태의 에너지를 암 조직에 전달하면서 주변의 정상 조직을 보호하는 것이다. 몬테카를로 모델링은 선형 가속기의 콜리메이션 업스트림에서 산란뿐만 아니라 환자 조직에서 산란으로 인해 환자가 경험하게 될 말초 선량을 결정하기 위해 방사선 치료에 일반적으로 사용된다.

광역학요법

광역학 치료(PDT)에서는 화학요법제를 활성화하기 위해 빛을 사용한다. PDT의 특성상 적절한 수준의 빛이 전달되어 항암화학요법을 활성화할 수 있도록 조직 내 산란과 흡수를 모델링하는 몬테카를로 방법을 사용하는 것이 유용하다.

산란 매체에서의 광자 수송 구현

여기에 제시된 것은 광자 몬테 카를로 방법의 동종 무한 매체에 대한 모델이다. 그러나 이 모델은 다중 레이어드 미디어를 위해 쉽게 확장된다. 비균형 매체인 경우 경계를 고려해야 한다. 반무한 매체(광자가 상단 경계를 벗어나면 분실된 것으로 간주되는 경우)에 대해서는 특별히 고려해야 한다. 자세한 내용은 페이지 하단의 링크를 참조하십시오. 무한히 작은 포인트 소스(공간과 시간의 디락 델타 함수로 분석적으로 표현됨)를 이용하여 문제를 해결한다. 임의의 소스 기하학에 대한 반응은 그린의 함수(또는 충분한 공간 대칭성이 존재하는 경우 콘볼루션)의 방법을 사용하여 구성할 수 있다. 필요한 파라미터는 흡수계수, 산란계수, 산란상함수 등이다. (경계가 고려되면 각 매체에 대한 굴절률도 제공되어야 한다.) 광학 경로 길이를 사용하여 광자 비행의 총 경과 시간을 추적함으로써 시간 분해 반응을 찾는다. 임의의 시간 프로파일이 있는 출처에 대한 반응은 시간 내 콘볼루션을 통해 모델링될 수 있다.

단순화된 모델에서는 다음과 같은 분산 감소 기법을 사용하여 계산 시간을 단축한다. 광자를 개별적으로 전파하는 대신 특정 무게의 광자 패킷(일반적으로 단일성으로 초기화)을 생성한다. 광자가 탁한 매질에서 상호작용하면서 흡수에 따른 중량을 퇴적시키고 나머지 중량은 매질의 다른 부분으로 흩어지게 된다. 특정 애플리케이션의 관심도에 따라 임의의 수의 변수가 도중에 기록될 수 있다. 각 광자 패킷은 그것이 종료, 반사 또는 전송될 때까지 반복적으로 다음과 같은 번호의 단계를 거친다. 그 과정은 오른쪽의 도식도에 도해되어 있다. 시뮬레이션된 측정 결과가 원하는 신호 대 잡음 비를 가질 때까지 원하는 수의 광자 패킷을 시작하고 모델링할 수 있다. 몬테카를로 모델링은 난수를 포함하는 통계적 프로세스인 만큼, 우리는 많은 계산에 대한 의사 난수로서 전체에서 변수 ξ을 사용할 것이다.

1단계: 광자 패킷 시작

우리 모델에서, 우리는 굴절률과 일치하지 않는 매체를 입력하는 것과 관련된 초기 규격 반사율을 무시하고 있다. 이것을 염두에 두고, 우리는 단지 초기 방향뿐만 아니라 광자 패킷의 초기 위치를 설정할 필요가 있다. 글로벌 좌표계를 이용하면 편리하다. 우리는 위치를 결정하기 위해 세 개의 데카르트 좌표를 사용할 것이고, 전파 방향을 결정하기 위해 세 개의 방향 코사인을 사용할 것이다. 초기 시작 조건은 용도에 따라 달라지지만, 원점에서 초기화된 펜슬빔의 경우 초기 위치와 방향 코사인을 다음과 같이 설정할 수 있다(각 패킷의 초기 방향을 무작위화하여 등방성 선원을 쉽게 모델링할 수 있다).

${\displaystyle {\begin}x&=0\\{\text{position: }y&=0\\\\\\mu_{x}=0\\\\\direction 코사인: }}\mu _{y}=0\mu _{z}=1\end}}}}}}}정렬}}}}}}}}}}}}}}}}}}}}}}}}}}}}}ndata.$

2단계: 단계 크기 선택 및 광자 패킷 이동

단계 크기 s는 광자 패킷이 상호작용 사이트 간에 이동하는 거리를 의미한다. 스텝 사이즈 선택 방법은 다양하다. 다음은 동종 모델에 사용하는 광자 스텝 크기 선택(역분포 방법과 Beer-Lambert 법칙을 사용하여 도출)의 기본 형식이다.

${\displaystyle s=-{\frac {\ln \xi }{\mu _{t}}}}$

여기서 ${\displaystyle \xi}$는 난수이고 μ t {\$displaystyle {\$mu_${$는 총 상호작용 계수(즉, 흡수 및 산란 계수의 합계)이다$$.

스텝 크기를 선택하면 광자 패킷이 방향 코사인에 의해 정의된 방향으로 거리 s에 의해 전파된다. 이는 단순히 좌표를 다음과 같이 업데이트하기만 하면 쉽게 이루어질 수 있다.

${\displaystyle {\signed}x&\왼쪽 화살표 x+\mu _{x}s\\y&\왼쪽 화살표 y+\mu _{y}s\z&\leftarrow z+\mu _{z}s\ended}}}}}}$

3단계: 흡수 및 산란

광자 중량의 일부는 각 상호작용 현장에서 흡수된다. 중량의 이 부분은 다음과 같이 결정된다.

${\displaystyle \Delta W={\frac {\mu _{a}{\mu _{t}}}W}$

여기서 $\mu {\displaystyle {}_{}}$ ${\displaystyle a_{}}$ ${\$는 흡수$$ 계수다.

그런 다음 흡수 분포를 특정 연구에 관심이 있는 경우 가중치를 배열로 기록할 수 있다. 광자 패킷의 중량은 다음과 같이 업데이트되어야 한다.

$W-\Delta W\,}$

흡수에 이어 광자 패킷이 흩어진다. 광자 산란 각도의 코사인 가중 평균은 -1과 1 사이의 값을 갖는 산란 음이소트로피(g)로 알려져 있다. 광학 음이소트로피가 0이면 일반적으로 산란이 등방성임을 나타낸다. g가 값 1에 근접하면 이는 산란 현상이 주로 전방 방향임을 나타낸다. 광자 패킷(따라서 광자 방향 코사인)의 새로운 방향을 결정하기 위해서는 산란 위상 함수를 알아야 한다. 종종 헤니-그린스타인 위상 기능이 사용된다. 그 후 산란각인 θ을 다음과 같은 공식을 이용하여 결정한다.

${\displaystyle \cos \theta ={\begin{cases}{\frac {1}{2g}}\left[1+g^{2}-\left({\frac {1-g^{2}}{1-g+2g\xi }}\right)^{2}\right]&{\text{ if }}g\neq 0\\1-2\xi &{\text{ if }}g=0\end{cases}}}$

그리고 극각 φ은 일반적으로 0과 $2{\displaystyle \mathrm{}}$ ${\displaystyle \pi }$ 사이에 균일하게 분포되어 있는 것으로 가정한다$([\displaystyle 2\pi$}). 이러한 가정에 근거하여 다음과 같이 설정할 수 있다$.$

${\displaystyle \varphi =2\pi \xi {\frac {}{}}}$

이러한 각도와 원래의 방향을 바탕으로 새로운 방향 코사인을 찾을 수 있다. 새로운 전파 방향은 다음과 같이 글로벌 좌표계로 나타낼 수 있다.

${\displaystyle{\begin{정렬}\mu '_{)}&, ={\frac{\sin \theta(\mu_{x}\mu_{z}\cos\varphi -\mu_{y}\sin \varphi)}{\sqrt{1-\mu_{z}^{2}}}}+\mu _{)}\cos, ={\frac{\sin \theta(\mu_{y}\mu_{z}\cos\varphi +\mu_{x}\sin \varphi)}{\sqrt{1-\mu_{z}^{2}}}}+\mu _{y}\cos\\\mu '_{z}& \theta\\\mu '_{y}& \theta, =-{\sqrt{1-\mu_{z}^{2}}}\sin. \theta \cos \varphi +\mu _{z}\cos \theta \\\\ended}}}$

특별한 경우

${\displaystyle {\displaysty}\mu _{z}=1\ended}}}$

사용하다

${\displaystyle {\signified}\mu '_{x}&=\sin \cos \varphi \\\mu '_{y}&=\sin \theta \theta \\cos \theeta \\\ended{x}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}?$

또는

${\displaystyle {\displaysty}\mu _{z}=-1\ended}}}$

사용하다

${\displaystyle {\signified}\mu '_{x}&=\sin \cos \varphi \\\\mu '_{y}&=-\sin \theta \cos \theta \\comparted}}}}}}}}$

C-코드:

/*********************** Indicatrix ********************* *New direction cosines after scattering by angle theta, fi. * mux new=(sin(theta)*(mux*muz*cos(fi)-muy*sin(fi)))/sqrt(1-muz^2)+mux*cos(theta) * muy new=(sin(theta)*(muy*muz*cos(fi)+mux*sin(fi)))/sqrt(1-muz^2)+muy*cos(theta) * muz new= - sqrt(1-muz^2)*sin(theta)*cos(fi)+muz*cos(theta) *-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- */ void Indicatrix (double muxs, double muys, double muzs, double mutheta, double fi, double *muxd, double *muyd, double *muzd) {  double costheta = mutheta;  double sintheta = sqrt(1.0-costheta*costheta); // sin(theta)  double sinfi = sin(fi);  double cosfi = cos(fi);  if (muzs == 1.0) {    *muxd = sintheta*cosfi;    *muyd = sintheta*sinfi;    *muzd = costheta;  } elseif (muzs == -1.0) {    *muxd = sintheta*cosfi;    *muyd = -sintheta*sinfi;    *muzd = -costheta;  } else {    double denom = sqrt(1.0-muzs*muzs);    double muzcosfi = muzs*cosfi;    *muxd = sintheta*(muxs*muzcosfi-muys*sinfi)/denom + muxs*costheta;    *muyd = sintheta*(muys*muzcosfi+muxs*sinfi)/denom + muys*costheta;    *muzd = -csn*sintheta*cosfi + muzs*costheta; }} 

4단계: 광자 종단

만약 광자 패킷이 많은 상호작용을 경험했다면, 대부분의 애플리케이션에서 패킷에 남아 있는 무게는 거의 중요하지 않다. 그 결과, 충분히 작은 무게의 광자 패킷을 종단하는 수단을 결정할 필요가 있다. 간단한 방법으로는 임계값을 사용할 수 있으며, 광자 패킷의 무게가 임계값보다 낮으면 패킷은 죽은 것으로 간주된다. 앞서 언급한 방법은 에너지를 절약하지 않기 때문에 한계가 있다. 총 에너지를 일정하게 유지하기 위해, 러시아 룰렛 기법은 광자에 대해 특정 중량 한계치 이하에 채택되는 경우가 많다. 이 기술은 광자의 생존 여부를 결정하기 위해 룰렛 상수 m을 사용한다. 광자 패킷은 m에서 생존할 수 있는 한 번의 기회를 가지며, 이 경우 W가 초기 중량인 mW의 새로운 중량을 받게 된다(이 새로운 중량은 평균적으로 에너지를 절약한다). 다른 모든 경우 광자 중량은 0으로 설정되고 광자는 종료된다. 이는 다음과 같이 수학적으로 표현된다.

${\displaystyle W={\begin{case}mW&\xi \leq 1/m\\0&\xi >1/m\end{case}}}$

광자 전송의 그래픽 처리 장치(GPU) 및 고속 몬테카를로 시뮬레이션

탁한 매체에서의 광자 이동의 몬테카를로 시뮬레이션은 많은 수의 광자가 독립적으로 전파되지만 동일한 규칙과 다른 무작위 번호 순서에 따라 전파되는 매우 병렬 가능한 문제다. 이러한 특수 유형의 몬테카를로 시뮬레이션의 병렬 특성은 그래픽 처리 장치(GPU)에서 실행하기에 매우 적합함을 렌더링한다. 프로그램 가능한 GPU의 출시는 그러한 개발을 시작했으며, 2008년부터는 광자 이주의 고속 몬테카를로 시뮬레이션에 GPU를 사용하는 것에 관한 몇 가지 보고가 있었다.^[1][2][3][4]

이 기본 접근법 자체는 함께 연결된 여러 GPU를 사용함으로써 병렬화될 수 있다. 저자의 홈페이지(GPU 클러스터 기반 다층 탁도 매체에서 광수송의 몬테카를로 시뮬레이션)에서 다운로드할 수 있는 「GPU 클러스터 MCML」이 대표적이다.http://bmp.hust.edu.cn/GPU_Cluster/GPU_Cluster_MCML.HTM

참고 항목

• 생물조직에서 광자수송을 위한 복사전달방정식과 확산이론
• 몬테카를로법
• 산란용 매체에서 광학 광폭반응을 위한 콘볼루션
• 몬테카를로 전자 운송 방법

기타 몬테카를로 자원 링크

• 미국 워싱턴대학교 광학영상연구소 루이(MCML)
• 오리건 메디컬 레이저 센터
• 룬드 대학의 광자 이주 몬테 카를로 연구, 스웨덴 GPU 가속 몬테 카를로 시뮬레이션 및 확장 가능한 몬테 카를로 연구. 다운로드할 소스 코드를 여십시오.
• 탁한 산란 매질에서 광수송을 위한 구름 기반 몬테카를로 연구 및 비상업 활동에 자유롭게 사용할 수 있는 도구.
• 몬테카를로 시뮬레이션(C++ 소스 코드 포함)의 예로서 조직 내 광 전송

참조

• Wang, L-H; Wu Hsin-I (2007). Biomedical Optics: Principles and Imaging. Wiley.
• L.-H. Wang; S. L. Jacques; L.-Q. Zheng (1995). "MCML—Monte Carlo modeling of light transport in multi-layered tissues". Computer Methods and Programs in Biomedicine. 47 (2): 131–146. doi:10.1016/0169-2607(95)01640-F.
• L.-H. Wang; S. L. Jacques; L.-Q. Zheng (1997). "Conv—convolution for responses to a finite diameter photon beam incident on multi-layered tissues" (PDF). Computer Methods and Programs in Biomedicine. 54 (3): 141–150. doi:10.1016/S0169-2607(97)00021-7.
• S. L. Jacques; L.-H. Wang (1995). "Monte Carlo modeling of light transport in tissues" (PDF). In A. J. Welch; M. J. C. van Gemert (eds.). Optical Thermal Response of Laser Irradiated Tissue. New York: Plenum Press. pp. 73–100.
• L.-H. Wang; S. L. Jacques (1994). "Optimized radial and angular positions in Monte Carlo modeling" (PDF). Medical Physics. 21 (7): 1081–1083. Bibcode:1994MedPh..21.1081W. doi:10.1118/1.597351. PMID 7968840.

인라인 참조