An Accelerated Approach to Dose Distribution Calculation in Inverse Treatment Planning for Brachytherapy

근접 치료에서 역방향 치료 계획의 선량분포 계산 가속화 방법

Abstract

With the recent development of static and dynamic modulated brachytherapy methods in brachytherapy, which use radiation shielding to modulate the dose distribution to deliver the dose, the amount of parameters and data required for dose calculation in inverse treatment planning and treatment plan optimization algorithms suitable for new directional beam intensity modulated brachytherapy is increasing. Although intensity-modulated brachytherapy enables accurate dose delivery of radiation, the increased amount of parameters and data increases the elapsed time required for dose calculation. In this study, a GPU-based CUDA-accelerated dose calculation algorithm was constructed to reduce the increase in dose calculation elapsed time. The acceleration of the calculation process was achieved by parallelizing the calculation of the system matrix of the volume of interest and the dose calculation. The developed algorithms were all performed in the same computing environment with an Intel (3.7 GHz, 6-core) CPU and a single NVIDIA GTX 1080ti graphics card, and the dose calculation time was evaluated by measuring only the dose calculation time, excluding the additional time required for loading data from disk and preprocessing operations. The results showed that the accelerated algorithm reduced the dose calculation time by about 30 times compared to the CPU-only calculation. The accelerated dose calculation algorithm can be expected to speed up treatment planning when new treatment plans need to be created to account for daily variations in applicator movement, such as in adaptive radiotherapy, or when dose calculation needs to account for changing parameters, such as in dynamically modulated brachytherapy.

최근 근접 치료에서 방사선 차폐막을 사용하여 선량 분포를 변조하여 선량을 전달하는 정적 및 동적 변조 근접 치료 방법이 개발됨에 따라 새로운 방향성 빔 세기 변조 근접 치료에 적합한 역방향 치료 계획 및 치료 계획 최적화 알고리즘에서 선량 계산에 필요한 파라미터 및 데이터의 양이 증가하고 있다. 세기 변조 근접 치료는 방사선의 정확한 선량 전달이 가능하지만, 파라미터와 데이터의 양이 증가하기 때문에 선량 계산에 필요한 경과 시간이 증가한다. 본 연구에서는 선량 계산 경과 시간의 증가를 줄이기 위해 그래픽 카드 기반의 CUDA 가속 선량 계산 알고리즘을 구축하였다. 계산 과정의 가속화 방법은 관심 체적의 시스템 행렬 계산 및 선량 계산의 병렬화를 이용하여 진행하였다. 개발된 알고리즘은 모두 인텔(3.7GHz, 6코어) CPU와 단일 NVIDIA GTX 1080ti 그래픽 카드가 장착된 동일한 컴퓨팅 환경에서 수행하였으며, 선량 계산 시간은 디스크에서 데이터를 불러오고 전처리를 위한 작업 등의 추가 적으로 필요한 시간은 제외하고 선량 계산 시간만 측정하여 평가하였다. 그 결과 가속화된 알고리즘은 CPU로만 계산할 때보다 선량 계산 시간이 약 30배 단축된 것으로 나타났다. 가속화된 선량 계산 알고리즘은 적응방사선치료와 같이 매일 변화되는 어플리케이터의 움직임을 고려하여 새로운 치료 계획을 수립해야 하는 경우나 동적 변조 근접 치료와 같이 선량 계산에 변화되는 파라미터를 고려해야 하는 경우 치료 계획 수립 속도를 높일 수 있을 것으로 판단된다.

Ⅰ. INTRODUCTION

최신 방사선 치료는 종양의 범위를 보다 구체적으로 정확하게 설정하고 정상 조직의 방사선 조사를 최소화하여 치료하는 입체세기조절회전 방사선 치료, 영상유도 방사선 치료, 정위 방사선 치료와 같은 외부 방사선 치료 방식의 발전에도 불구하고, 근접 치료(Brachytherapy; BT)는 여전히 많은 자궁경부암 환자에게 필수적인 치료 방법이다^[1-3]. 근접치료는 자궁경부에 있는 종양에 많은 방사선을 조사함과 동시에 종양과 인접한 건강한 조직에 대한 선량 감소를 제공할 수 있어, 외부 방사선 치료(External Beam Radiation Therapy; EBRT) 방식보다 정상조직에 전달되는 선량을 최소화하고 종양 중심으로 선량을 전달하여 치료범위의 선량분포를 향상시킬 수 있다^[4,5]. 하지만 이러한 근접 치료에서 사용되는 방사선원은 등방성 선량분포를 가지고 있고, 불규칙한 모양의 종양이나 위험 장기(Organ at Risk; OAR)에 근접한 종양의 경우 등방성 선량 분포를 이용하면 표적 체적(Target volume)에 원하는 선량이 전달되지 않을 수 있다^[6,7].

최근에는 원자번호가 높은 물질로 차폐된 어플리케이터를 사용하여 근접 치료의 방사선원으로 부터 전달되는 선량의 강도를 조절하는 기술이 연구 되고 있고^[8], 선량 전달 중에 차폐체가 선원 또는 주변 조직에 비해 움직이지 않는 정적(static)^[9-11] 그리고 치료 중 선원 또는 주변 조직에 대해 차폐체의 방향이 바뀌는 기술인 동적(dynamic) 변조 근접 치료로 나눌 수 있다^[12-14].

정적 및 동적 차폐체를 이용한 선량의 강도 조절기술은 기존의 선량 계산 방법에 필요한 선원 위치및 선원이 머무르는 시간(dwell time)에 더하여 방사선 방출각도라는 변수를 추가하여 선량을 계산하여야 하기 때문에, 선원에 의해 전달되는 선량의 정확한 계산을 위해 세가지 변수를 최적화할 수 있는 역방향 치료계획(Inverse Treatment Planning)이 필요하다. 특히, 정적 변조 근접 치료에서 역방향 치료계획은 선원의 위치, 위험 장기(OAR)의 위치, 표적 체적의 위치 및 차폐체의 위치를 고려하여 선량분포를 계산해야 하며, 동적 변조 근접 치료에서는 시간에 따라 변하는 차폐체에 의한 선량분포를 추가적으로 고려해야만 한다. 따라서, 근접 치료의 역방향 치료계획에서 정확한 선량분포를 얻기 위하여 계산량이 많아지며 계산속도가 느려질 수 밖에 없다. 이에 본 연구에서는 근접 치료에서 선량 강도 조절 기술의 적용으로 증가된 선량 최적화 계산량을 GPU (Graphics Processing Unit)기반의 CUDA(Compute Unified Device Architecture) 가속 프로그래밍을 적용하여 선량 계산 속도를 향상시키는 연구를 진행하였다.

Ⅱ. MATERIAL AND METHODS

1. 역방향 치료계획

선량 계산에 사용된 영상 데이터는 CT 장비(Somatom Definition Flash, Siemens Healthcare, Forchheim, Germany)를 사용하여 획득하였고, 복셀크기는 0.5 × 0.5 × 4.0 mm³이었다. 원형 타입으로 구성된 CT/MR (Computed Tomography/Magnetic Resonance Imaging) 호환 애플리케이터 세트 (텐덤-원형 타입 애플리케이터, Elekta AB, 스웨덴 스톡홀름)를 사용하여 치료계획을 작성하였으며, 어플리케이터는 원형 부분의 각도는 30도, 선원 위치 간격은 2.5 mm, 첫 번째 소스에서 애플리케이터 팁까지의 거리는 7 mm로 정의되어 있다. Fig. 1과 같이 선원, 표적 체적 및 위험 장기 내부의 복셀모델링에 대한 역방향 치료계획에서의 선량으로 Eq. (1)을 이용하여 계산하였다.

A_kx = b_k       (1)

Fig. 1. Object and source modeling of the contour point-lists extracted from the DICOM-RT structure files, with resampled points on the three-dimensional model for accurate dose calculation based on dwell time of each source.

여기서, A_k 는 선량계산이 필요한 관심 체적의 시스템 행렬, x는 선원의 체류 시간으로, 선원의 체류시간이 결정되면 역방향 치료 계획에서의 관심 체적에 전달되는 선량분포 b_k가 결정된다. 역방향 치료계획에서 A_k 의 시스템 행렬은 필수적이다^[15].

2. 선량분포 계산 가속화 방법

선량분포 계산은 고도의 연산 과정이 필요하며, 연산 과정의 효율을 높이기 위해 계산 과정에서 처리할 데이터의 양과 알고리즘의 복잡도에 따라 병렬적으로 수행해야 한다. 본 연구에서는 계산 효율을 높이기 위해 Fig. 2와 같이 선원 위치에 따른 표적 체적 및 위험 장기(OAR) 내부의 복셀에 대한 선량 행렬만 계산하였다. 즉, 컴퓨팅 가속을 위해 DICOM-RT 구조의 윤곽이 입력되어 GPU 텍스처 메모리로 전송하여, 가져온 윤곽 목록에서 포인트 샘플링을 수행하여 실제 선량분포 계산 과정에 필요한 관심 부분의 체적 행렬을 생성하였다.

Fig. 2. Schematic diagram of the dose delivered to one voxel by each source.

그리고, 실제 선량분포 계산은 병렬적으로 GPU 텍스처 메모리에서 각 선원 위치와 관심 부분의 체적 형렬을 Eq. (1)과 AAPM TG-43 선량 권고 및 GET ESTRO 작업 그룹의 지침에 따라 계산 하여, 역방향 치료계획을 작성하였다^[7,16]. 전반적으로 병렬 연산속도의 향상뿐만 아니라 그래픽 메모리와 연산 장치간의 데이터 전송 향상에도 중점을 두었다.

3. 정량적 평가

선량분포 계산 가속 속도를 측정하기 위해 CPU 및 GPU 프로그래밍을 사용하여 각 선량 계산 방식에 소요되는 계산 시간을 비교하였다. 두 계산 알고리즘은 모두 인텔(3.7 GHz, 6코어) CPU와 단일 GPU (NVIDIA GTX 1080ti) 카드가 장착된 동일한 컴퓨팅 환경에서 수행되었다.

선량 계산 시간은 디스크에서 데이터를 로드하고 사전작업을 위한 작업 등의 추가적으로 필요한 시간은 제외하고 선량 계산 시간만 측정하였다.

Ⅲ. RESULT

1. 선량 계산 정확도

근접 치료에서 많이 사용되는 외경이 1.0 mm인 원통형 스테인리스강 캡슐에 밀봉되어 있는 3.5 mm 192-Ir 선원을 이용하였다. 선원은 특정 한 지점에 위치시켰으며, 표적 체적 및 위험 장기의 16개 무작위 위치를 선정하여 선량 계산 정확도를 확인하기 위하여, Table 1에 임상에서 사용되는 상업용 선량분포 계산 프로그램과 본 연구에서 개발한 선량 계산 알고리즘을 이용하여 계산된 선량의 결과를 나타내었다. 관심 위치에서 무작위로 선택한 지점들의 계산된 선량 차이의 평균은 – 0.0983%로 높은 계산 정확도를 확인하였다.

Table 1. Comparison between doses calculated by the commercial RTP system and computed by the in-house program utilizing a point source

2. 역방향 치료 계획

Fig. 3은 실제 임상에서 자궁경부암 근접 치료 시 권고되어 사용하고 있는 선량을 기준 (5 ~ 6 분할 치료, 총 30 Gy 선량 전달)을 이용하여 일반적으로 사용되어지는 근접 치료 계획(Nominal plan)과 Eq. (1)을 이용하여 역방향 치료 계획(Inverse plan)을 이용한 선원의 체류 시간을 산출하여 나타낸 결과이다.

Fig. 3. The corresponding dwell times for a (a) nominal plan and (b) inverse plan method at each of the 15 dwell position.

선원이 머무를 수 있는 각 위치들은 사용된 어플리케이터의 특성에 따라, 15개의 선원 위치(dwell position)를 기반으로 선원의 체류 시간을 산출하여 나타내었다.

각 위치들에서 선원의 체류시간이 선량 계산 방법에 따라 다르게 결정된 것을 확인 할 수 있었다. Fig. 3에서 두 방법에 의해 결정된 선원이 머무르는 시간의 총합은 각각 Fig.3 - (a)에서 218.38초, Fig.3 - (b)에서 215.51초로 나타났으며, 2.87초로 미세하게 역방향 치료 계획에서 선원이 머무르는 총 시간의 감소를 확인할 수 있었다. Table 2는 표적 체적 및 위험 장기들의 계산된 선량 차이를 나타낸 결과이며, 선량은 Table 3의 GET ESTRO 작업 그룹의 지침인 EMBRACE 권고에 따라서 계산 하였다^[7]. 각 관심 위치들에서 계산된 선량 차이의 전체 평균은 -2.88%로 역방향 치료계획에서 선원의 총 체류 시간을 감소시키면서 각 관심 위치의 선량은 EMBRACE 권고를 만족하도록 치료계획을 작성할 수 있음을 확인하였다.

Table 2. Comparison between doses calculated by the commercial RTP system and computed by the in-house program

Table 3. The dose constraints in accordance with EMBRACE recommendation and the total EQD2 (EBRT: 50.40 Gy + BT) of each target and OAR

3. 선량 계산 가속화

Fig. 4의 선량 계산 가속화 알고리즘은 Python 기반 CUDA 언어로 작성되었으며, 인텔(3.7GHz, 6코어) CPU와 단일 GPU (NVIDIA GTX 1080ti) 카드가 장착된 컴퓨터에서 실행되어졌다. Table 4는 CPU와 GPU를 사용한 선량 계산 시간 비교 결과 나타낸다. 역방향 치료 계획에서 선량 계산 시간의 경우 GPU를 사용할 경우 CPU를 사용한 수행보다 약 30.086배 만큼 가속화 비율을 보였다.

Table 4. Total elapsed time for producing the dose calculation in inverse plan

Fig. 4. Pseudo-programming of parallel computation in the sampling and dose calculation operations.

Ⅳ. DISCUSSION

전통적인 방사선치료 방법 중 등방성의 선량분포를 가지는 근접 치료의 경우 관심 체적 부위에 원하는 비대칭 형태의 선량을 정확하게 전달하기 위해 다양한 방법이 고안되어 지고 있다. 본 연구에서는 근접 치료에서 새롭게 고안되어지고 있는 비대칭 형태의 선량 전달 방법을 효과적으로 적용시킬 수 있는 역방향 치료 계획을 작성하여 확인하였다. 또한, 근접 치료의 역방향 치료 계획에서 선량분포를 계산함에 있어 속도를 개선하는 연구는 지금까지 수행된 바가 없다. 따라서, 변수의 증가로 많아진 계산량으로 인하여 발생된 선량 계산 속도의 저하를 CUDA 가속 프로그래밍을 적용하여 선량 계산 속도를 향상시키는 연구를 진행하였으며, 선량 계산에 필요한 경과 시간을 39.443초에서 1.311초로 약 30배 만큼 단축시킬 수 있어 만족스러운 결과를 보여주었다. 본 연구에서는 Eq. (1)에서 볼 수 있듯이 선량 계산 시 관심 체적의 시스템 행렬, 선원이 위치할 수 있는 곳에서의 선원의 체류 시간을 최적화하여 관심 체적에 전달되는 선량분포를 계산해야만 하기 때문에, 선량 계산량이 증가되어 계산 속도가 현저히 느려질 수 밖에 없다. 따라서, 과심 체적의 시스템 행렬과 선원의 체류시간을 통한 선량 최적화 계산의 병렬화를 진행하였으며, CUDA 가속 프래그래밍을 기반으로 선량계산 알고리즘의 추가적인 가속화를 진행 할 수 있었다. 본 연구에서 단일 GPU 카드를 사용하여 가속화 알고리즘을 적용하였지만, 고성능 GPU 및 많은 수의 GPU 카드를 병렬화 하여 선량 계산을 진행하면 추가적인 계산 속도의 가속화 결과를 기대할 수 있을 것이다. 또한, 표적 체적 및 위험 장기 내부의 복셀모델링에 대한 역방향 치료계획에서의 선량을 계산하기위해 사용된 Eq. (1)에 표적 체적 및 위험 장기에 대한 제약조건을 추가하면 선량 계산의 정확도를 향상시켜 최적화할 수 있을 것으로 판단되며, 향후 연구에서 고성능 GPU 카드의 병렬화를 통한 선량 계산속도 향상 및 선량 계산 알고리즘의 고도화를 기반으로 정확도의 향상을 진행하고자 한다.

Ⅴ. CONCLUSION

근접 치료의 역방향 치료 계획에서의 선량 계산을 GPU 기반의 CUDA 가속 프로그래밍을 적용하여 선량 계산 속도를 가속화 하였다. 본 연구에서는 인텔(3.7GHz, 6코어) CPU와 단일 GPU (NVIDIA GTX 1080ti) 카드를 사용하여 약 30배의 가속화 결과를 얻을 수 있었다. 특히, 주변 조직에 대해 차폐체의 방향이 바뀌는 기술인 동적 (dynamic) 변조 근접 치료에서 선량 계산 시, 추가적으로 “회전”에 관련된 인자를 추가하여 최적화해야 하기 때문에 GPU 기반의 CUDA 가속 프로그래밍 적용을 통한 가속화를 기대할 수 있었다. 추가적으로, 근접 치료에서 어플리케이터의 고정을 완벽하지 않을 경우 기존의 치료 계획과 다른 선량이 전달될 수 있기 때문에 적응방사선치료(adaptive radiation therapy)가 고려되어야 하며, 선량 계산 가속화 알고리즘을 기반으로 한 실시간 방사선치료계획의 수정으로 빠른 치료 계획의 작성이 가능하여, 근접 치료에서 선량 전달의 정확성을 향상시킬 수 있을 것이다.