KR101471646B1 - Ｇｐ－ｇｐｕ를 이용한 3ｄ 의료 영상 정합의 병렬처리방법 - Google Patents

Abstract

본 발명은 같은 환자에 대하여 서로 다른 시간에 획득한 의료영상으로부터 질환을 모니터링하기 위해 필요한 의료영상의 정합기법에 관한 것으로, 본 발명에 따르면, 종래의 분할기반의 기법과 화소값기반 기법의 장점을 취합하여, 3차원 의료영상 데이터를 보다 빠르고 정확하게 처리할 수 있도록, GP-GPU 기반 병렬처리를 적용하는 것을 특징으로 하는 GP-GPU를 이용한 3D 의료영상 정합의 병렬처리방법이 제공된다.

Description

ＧＰ－ＧＰＵ를 이용한 3Ｄ 의료 영상 정합의 병렬처리방법{Parallel processing of 3D medical image registration by GP-GPU}

본 발명은, 3차원 의료영상의 정합 처리를 고속으로 처리하는 방법에 관한 것으로, 더 상세하게는, GP-GPU(General-Purpose computation on Graphics Processing Units)를 이용한 3D 의료 영상 정합의 병렬처리방법에 관한 것이다.

종래, 예를 들면, 관절염(Osteoarthritis)과 같은 질환을 치료하기 위해서는, 일반적으로, 도 7에 나타낸 바와 같이, 질환이 있는 관절의 위치를 CT(Computed Tomography), MRI(Magnetic Resonance Imaging), PET(Positron Emission Tomography) 및 SPECT(Single Photon Emission Computed Tomography) 등과 같은 스캐너(scanner)를 사용하여 의료영상을 2차원 영상 또는 이미지 슬라이드로 촬영하고, 이러한 영상들을 이용하여 질환을 모니터링(monitoring) 하거나 치료를 행하게 된다.

또한, 질환을 모니터링 하기 위해서는, 동일한 환자에 대하여 적어도 두 세트의 의료영상을 이용해야 하며, 여기서, 두 세트의 의료영상이란, 예를 들면, 환자가 처음으로 입원할 때 찍었던 영상과, 6개월이나 1년 후에 찍은 영상과 같이, 서로 다른 시간에 각각 얻어진 영상을 의미하고, 이러한 두 영상을 각각 입력(source) 영상과 목표(target) 영상이라 한다.

아울러, 질환이 있는 위치의 영상을 획득할 때, 환자가 움직이거나 또는 질환이 발생한 경우, 의사들이 2차원 영상을 한 장씩 눈으로 보고 판단하기는 용이하지 못하므로, 이를 위해, 영상처리기법 중 하나인 3차원 영상 정합기법이 필요하게 된다.

즉, 일반적으로, 영상 정합기법은, 하나의 대상을 다른 시간이나 관점에서 촬영한 영상이나, 서로 다른 좌표계에서 얻어진 영상을 같은 좌표계로 다시 매칭(matching) 해주는 기법을 말한다.

또한, 의료영상에 있어서, 영상 정합은, 서로 다른 스캐너에서 획득한 두 개의 영상을 통합할 때나, 동일한 스캐너를 사용하여 서로 다른 시간에 얻은 두 개의 영상을 통합할 때 중요한 역할을 한다.

여기서, 예를 들면, 머리 CT 영상들을 매칭하는 정합기법과 무릎 CT 영상들을 매칭하는 정합기법은, 원하는 결과가 상이하므로 그 정합기법이 서로 다를 수 있다.

또한, 기존에 개발된 정합기법들은, 정합 방식에 따라 크게 특징기반(Feature-based) 기법, 분할기반(Segmentation-based) 기법 및 화소값 기반(Intensity-based) 기법으로 분류된다.

먼저, 특징기반 기법은, 입력 영상과 목표 영상에서 해당된 특징들을 검출하고, 그 특징들 또는 랜드마크(Landmark)들을 이용하여 정합을 수행하며, 여기서, 영상처리에서의 특징이란, 한 영상을 대표하는 모든 것들을 그 영상의 특징이라고 말할 수 있다.

즉, 영상에서 통속적인 특징들은, 영상의 점(points), 코너(corners), 선(Lines), 경계(boundary), 구분선(contour) 등이며, 인체의 부위에 따라 사용되는 특징도 다르다.

아울러, 이러한 특징들을 찾을 때에는, 전문의가 수동으로 찾을 수도 있고, 또는, 예를 들면, SIFT, SURF 방법들과 같이, 영상에서 특징을 검출하는 자동 영상처리 기반 알고리즘들을 사용하여 자동으로 찾도록 할 수도 있다.

따라서 특징기반의 기법은, 특징만 매칭하기 때문에 정합 수행시간이 오래 걸리지 않지만 정확성이 떨어지는 단점이 있고, 반면, 전문가가 수동으로 특징을 찾으면 정확성은 높일 수 있으나 수행시간이 오래 걸리게 된다.

또한, 분할기반의 정합기법은, 특징기반 정합기법과 방식이 유사하며, 즉, 분할기반 기법은, 특징기반 기법과 마찬가지로 두 영상(입력 영상과 목표 영상)에서 변화가 없는 부분을 찾고, 그 부분에 대하여 매칭을 수행한다.

더 상세하게는, 예를 들면, 무릎 영상에서 뼈 부분을 분할할 때, 수동으로 분할하는 방식과 자동으로 분할하는 방식이 있고, 영상처리 분할 알고리즘 중 ACM(Active Contour model)과 Mean Shift 알고리즘들은 의료영상처리 분야에 자주 사용되고 있는 알고리즘들이며, 분할기반 기법도 특징기반 기법과 마찬가지로 영상의 분할한 영역만 매칭하기 때문에 정합시간이 짧다.

또한, 상기한 특징기반의 기법 및 분할기반의 기법과는 별개로, 화소값 기반의 정합기법은, 영상에서 특징을 검출하거나 분할영역을 찾지 않고 영상의 화소값을 그대로 이용하는 기법이다.

아울러, 이러한 화소값을 이용한 기법은, 전체 이미지(whole image)를 사용하는 경우가 있고, 서브 이미지(sub-image)만 사용하는 경우가 있다.

먼저, 서브 이미지로 매칭하는 방식을 템플릿 매칭(template matching) 방식이라 하며, 의료영상에서는 서브 이미지로 매칭하는 방식보다는 전체 이미지를 사용하여 매칭하는 방식이 많이 이용되고, 이는, 이러한 방식이 임상 응용(clinical applications)을 위해 신뢰할 수 있는 결과를 제공하기 때문이다.

또한, 두 영상(입력 영상과 목표 영상)의 화소값을 비교할 때, 유사성 측정 메트릭(similarity measure metric)을 이용하여 계산을 행하며, 원본 이미지에 따라 사용되는 메트릭도 다르다.

예를 들면, MRI 영상과 CT 영상을 정합할 때 적용하는 메트릭과, CT 영상과 CT영상을 정합할 때 적용하는 메트릭은 상이하며, 일반적인 메트릭들은 MI(mutual information), NMI(normalized mutual information), NCC(normalized cross correlation), SSD(sum of squared differences) 및 SAD(sum of absolute differences) 등이 사용된다.

여기서, MI 및 NMI는, 예를 들면, MRI/CT와 같이, 서로 다른 스캐너로 획득한 두 영상을 정합할 때 사용하고 NCC, SSD, SAD는, 예를 들면, CT/CT와 같이 같은 스캐너로 획득한 두 영상을 정합할 때 적용된다.

따라서 상기한 바와 같이, 화소값을 이용한 기법은 전체 이미지의 화소값에 메트릭을 사용하여 비교하기 때문에 높은 정확도를 보장하지만, 수행시간이 오래 걸린다.

상기한 바와 같이, 기존에 제시된 종래의 영상 정합기법들은 각자 장단점을 가지고 있으며, 즉, 특징기반의 기법과 분할기반의 기법은 정합시간이 빠르나 매칭하기 전에 특징이나 분할 영역을 구해야 하기 때문에 정합 정확도를 보장하지 못하며, 반면, 화소값 기반의 방식은 영상의 화소값을 그대로 계산하기 때문에 정확도를 보장할 수 있으나 픽셀값을 하나씩 비교하므로 수행시간이 많이 걸린다는 단점이 있다.

따라서 의료영상 처리분야에 있어서 상기한 바와 같은 종래기술의 문제점을 해결하기 위하여는, 정확도를 보장하면서 동시에 수행시간도 짧은 새로운 정합기법의 개발이 요구되나, 아직까지 그러한 요구를 모두 만족시키는 의료영상의 정합기법은 제시된 바 없었다.

아울러, 의료 분야에 있어서, 고가의 연산처리 비용(computational cost)으로 인해 3D 영상을 적용하는 경우는 매우 제한적이었으나, 최근, 멀티코어 CPU 및 그래픽 전용의 처리장치(GPU)의 사용이 증가하면서, 병렬 처리에 의해 연산 시간의 문제는 점차 해소되고 있다.

더 상세하게는, GPU를 이용하여 영상정합을 가속하는 기술에 있어서, 종래, 예를 들면, NVIDIA사의 GPU를 프로그래밍하기 위해 C/C++과 유사한 프로그래밍 환경인 CUDA(Compute Unified Device Architecture)를 이용하는 GPU 병렬처리 방법이 있으며, 이러한 CUDA에 기반하여 다양한 영상처리의 가속방법이 제안되어 왔다.

그러나 상기한 바와 같은 CUDA를 이용한 알고리즘들은, NVIDIA사의 하드웨어에서만 동작한다는 문제점이 있었고, 따라서 이러한 문제점을 해결하기 위하여는,특정 회사의 특정 제품이 아닌, 범용으로 어느 하드웨어에서나 동작 가능한 알고리즘을 제공하는 것이 바람직하다.

즉, 예를 들면, OpenCL과 같이, 특정한 하드웨어를 위해 개발된 것이 아닌 보다 일반적인 병렬처리 구조를 이용하여 GPU에 의해 3D 의료영상의 정합처리를 더욱 가속할 수 있는 알고리즘이나 방법을 제공하는 것이 바람직하나, 아직까지 그러한 요구를 모두 만족시키는 알고리즘이나 방법 또한 제시된 바 없었다.

본 발명은 상기한 바와 같은 종래기술의 문제점을 해결하고자 하는 것으로, 따라서 본 발명의 목적은, 서로 다른 시간에 획득한 같은 환자의 의료 영상에서 질환을 모니터링하기 위해 필요한 의료영상의 정합 기법에 있어서, 종래의 분할기반의 기법과 화소값기반 기법의 장점을 종합하여, 보다 빠르고 정확한 의료영상의 정합방법을 제공하고자 하는 것이다.

또한, 본 발명의 다른 목적은, 3차원 데이터를 빠르게 처리할 수 있도록 GP-GPU(General-Purpose computation on Graphics Processing Units) 기반 병렬처리를 적용하여, 종래의 의료영상 정합기법에 비하여 보다 빠르고 정확한 3D 의료영상 정합의 병렬처리방법을 제공하고자 하는 것이다.

상기한 바와 같은 목적을 달성하기 위해, 본 발명에 따르면, 서로 다른 시간에 획득한 동일한 환자의 2개의 영상을 정합하는 처리를 고속으로 처리하기 위한 영상 정합의 병렬처리방법에 있어서, GP-GPU(General-Purpose computation on Graphics Processing Units)를 이용하여 상기 2개의 영상을 정합하는 처리가 CPU와 GPU에 의해 병렬처리되는 것을 특징으로 하는 영상 정합의 병렬처리방법이 제공된다.

여기서, 상기 방법은, 서로 다른 시간에 획득된 동일한 환자의 2개의 영상을 각각 입력받는 입력 단계; 상기 입력받는 단계에서 입력된 영상들을 스케일링(scaling) 하는 스케일링 단계; 상기 스케일링된 영상들로부터 초기 변환 파라미터(initial transformation parameter)를 추출하는 초기화 단계(initialization step); 상기 초기화 단계에서 추출된 상기 초기 변환 파라미터들을 최적화하여 최종 파라미터를 구하는 최적화 단계(optimization step); 및 상기 최적화 단계에서 최적화된 상기 최종 파라미터를 이용하여 상기 입력된 영상들을 정합하고 표시장치를 통하여 표시하는 시각화 단계(visualization step)를 포함하여 구성되며, 상기 입력 단계, 상기 스케일링 단계, 상기 초기화 단계 및 상기 시각화 단계는 CPU에 의해 수행되고, 상기 최적화 단계는 GPU에 의해 상기 CPU의 처리와 함께 병렬처리되는 것을 특징으로 한다.

또한, 상기 CPU는 C 언어로 작성된 프로그램에 의해 제어되고, 상기 GPU는 OpenCL 또는 CUDA를 이용하여 제어되도록 구성되는 것을 특징으로 한다.

아울러, 상기 스케일링 단계는, 상기 입력된 영상들의 복셀 사이즈(voxel size)에 따른 스케일링 파라미터를 이용하여 상기 입력된 영상들의 스케일링이 수행되는 것을 특징으로 한다.

더욱이, 상기 복셀 사이즈에 대한 정보는, 입력 영상의 DICOM(Digital Imaging and Communication in Medicine) 이미지 헤더 파일로부터 얻어지며, 상기 복셀 사이즈는, x-y 평면에서 픽셀 스페이싱(pixel spacing)을 통해 추출되고, z-평면에서는 입력 DICOM 영상의 헤더 파일 정보의 슬라이스(slice)로부터 얻어지고, 제 1 입력영상 V₁의 복셀 사이즈가 제 2 입력영상 V₂의 복셀 사이즈보다 작은 것으로 가정하면, 상기 스케일링 파라미터는 이하의 수학식을 이용하여 계산되며,

(여기서, V_1i는 볼륨 V₁의 복셀 사이즈이고, V_2i는 볼륨 V₂의 복셀 사이즈이다.)

스케일링 파라미터가 계산된 후, 상기 V₁은 스케일 다운되고, 상기 V₂는 기준 볼륨(reference volume)(V_r)으로 정의되며, 상기 V₁은 부동 볼륨(float volume)(V_f)으로 설정되는 것을 특징으로 한다.

또한, 상기 초기화 단계에서, 상기 초기 변환 파라미터는, 3개의 회전 파라미터와 3개의 변환 파라미터를 포함하며, 상기 기준 볼륨 V_r과 상기 부동 볼륨 V_f 사이의 상대 위치(related position)는, 변환 파라미터 P = {t_x, t_y, t_z, α, β, γ}(여기서, t_x, t_y, t_z는 변환량(translation quanta)이고 α, β, γ는 각각 기준 볼륨에 대하여 3D 축에 따른 부동 볼륨의 회전각)의 집합에 의해 정의되는 것을 특징으로 한다.

아울러, 상기 초기화 단계는, 상기 기준 볼륨 및 상기 부동 볼륨이 모두 2진화되는(binarized) 단계; 2진화된 상기 기준 볼륨 및 상기 부동 볼륨의 픽셀의 좌표로부터 3D 벡터를 형성되고, 형성된 상기 3D 벡터로부터 중심(centroiod) 및 관성행렬(inertia matrix)이 산출되는 단계; 각각의 상기 관성행렬의 고유벡터로부터 각 볼륨의 회전각이 산출되는 단계; 상기 기준 볼륨의 회전각(x, y, z)으로부터 상기 부동 볼륨의 회전각(x, y, z)을 감산함으로써(subtracting) 3개의 초기 회전 파라미터가 산출되는 단계; 및 상기 기준 볼륨의 중심(x, y, z)으로부터 상기 부동 볼륨의 중심(x, y, z)을 감산함으로써 3개의 초기 변환 파라미터가 산출되는 단계를 포함하는 것을 특징으로 한다.

여기서, 상기 2진화되는 단계는, B(x, y, z)를 3D 볼륨 V(x, y, z)의 초기 2진화 볼륨이라 하면, 이하의 수학식을 이용하여 상기 기준 볼륨 및 상기 부동 볼륨이 2진화되는 것을 특징으로 한다.

(여기서, x, y, z는 이미지의 복셀의 좌표(coordinates)이고, τ는 2진화 볼륨을 정의하는 임계값(threshold value)이다.)

또한, 상기 관성행렬은, 이하의 수학식으로 정의되며, 상기 관성행렬의 고유벡터(eigenvectors)에 의해 관성 주축(principal axes)이 정의되는 것을 특징으로 한다.

(여기서,

는 객체 모멘트(object moments) 이고, 함수 f(x, y, z)는 복셀 데이터의 이미지 내용(image content)을 나타내며, x_c, y_c, z_c는 상기 객체의 중심을 나타내고, 상기 관성행렬로부터 계산된 3개의 고유벡터는 상기 객체의 관성 주축을 나타낸다.)

아울러, 상기 고유벡터의 행렬 형태는 이하의 수학식으로 나타내지며,

회전행렬 R에 대하여 상기 행렬 E를 푸는 것에 의해(E = R), 상기 회전각 α, β, γ가 이하의 수학식으로 계산되고,

상기 회전행렬은, 이하의 수학식으로 나타내지는 것을 특징으로 한다.

R = R_γ×R_β×R_α

(여기서, α, β, γ는 각각 3D 축에 대한 오일러 각도(Euler angles) 이며,

이다.)

더욱이, 상기 최적화 단계는, 상기 부동 볼륨의 각 복셀에 대하여, 강체 변환 행렬(rigid body transfomation matrix) 및 3차-선형 보간법(tri-linear interpolation) 연산(operate)을 이용하여 상기 부동 볼륨을 변환하는 단계; 변환된 상기 부동 볼륨과 상기 기준 볼륨 사이의 유사성 스코어(similarity score)를 측정하는 단계; 및 상기 유사성 스코어가 최대가 될 때까지 모든 복셀에 대하여 상기한 단계들을 반복하는 단계를 포함하는 것을 특징으로 한다.

또한, 상기 최적화 단계는, 두 입력 영상 사이의 관계를 3개의 변환 파라미터 및 3개의 회전 파라미터에 의해 정의되는 강체변환(rigid body transformation) 행렬 M이라 하면, 상기 강체변환행렬 M은 이하의 수학식으로 나타내지는 것을 특징으로 한다.

M = T(t)R

(여기서, 상기 T(t) 및 상기 (R)은 동일 좌표계(homogeneous coordinates)에서의 변환벡터와 회전행렬이다.)

아울러, 상기 유사성 스코어의 측정은, 상기 최적화 단계에서 두 볼륨의 유사성의 정도(degree)를 정량화(quantify) 하기 위해 이하의 NCC(normalized cross-correlation) 함수를 이용하여 수행되는 것을 특징으로 한다.

(여기서, n은 복셀의 총 수(total number), i와 j는, x가 V_r(x_i) 및 V_f(x_j)의 x-, y-, z- 좌표를 나타내고 V_r과 V_f가 평균 화소값(mean intensity value)을 나타낼 때, 포인트 x_i와 x_j에서의 화소값을 나타내는 복셀 인덱스(voxel index) 이다.)

더욱이, 상기 방법은, 상기 최적화 단계에서, 상기 부동 볼륨을 변환하는 단계 및 상기 유사성 스코어를 측정하는 단계의 처리가 상기 GPU에 의해 병렬처리되는 것을 특징으로 한다.

상기한 바와 같이, 본 발명에 따르면, 종래의 분할기반의 기법과 화소값기반 기법의 장점만을 취합함으로써, 종래의 의료영상 정합기법들보다 더욱 빠르고 정확한 의료영상의 정합방법을 제공할 수 있다.

또한, 본 발명에 따르면, GP-GPU(General-Purpose computation on Graphics Processing Units) 기반 병렬처리를 적용함으로써, 종래의 분할기반의 기법과 화소값기반 기법의 장점만을 취합하여 3차원 의료영상 데이터를 보다 빠르고 정확하게 처리할 수 있는 GP-GPU를 이용한 3D 의료영상 정합의 병렬처리방법을 제공할 수 있다.

도 1은 2개의 영상을 정합하기 위한 영상정합 처리방법의 전체적인 처리 흐름을 개략적으로 나타내는 플로차트이다.
도 2는 본 발명의 실시예에 따른 GP-GPU를 이용한 3D 의료영상 정합의 병렬처리방법의 전체적인 구성을 개략적으로 나타내는 도면이다.
도 3은 OpenCL의 플랫폼 모델의 구성을 개략적으로 나타내는 도면이다.
도 4는 OpenCL의 NDRange 구성을 개략적으로 나타내는 도면이다.
도 5는 OpenCL의 메모리 모델의 구성을 개략적으로 나타내는 도면으로, OpenCL에 의해 정의되는 메모리 계층 구조(memory hierarchy)의 다이어그램을 나타내는 도면이다.
도 6은 3D 변환을 위한 OpenCL 구현의 설계방법의 개념을 개략적으로 나타내는 도면이다.
도 7은 영상 정합에 이용되는 환자의 무릎 이미지를 나타내는 도면이다.

이하, 첨부된 도면을 참조하여 본 발명에 따른 GP-GPU를 이용한 3D 의료영상 정합의 병렬처리방법의 상세한 내용에 대하여 설명한다.

여기서, 이하에 설명하는 내용은 본 발명을 실시하기 위한 실시예일 뿐이며, 본 발명은 이하에 설명하는 실시예의 내용으로만 한정되는 것은 아니라는 사실에 유념해야 한다.

즉, 본 발명은, 후술하는 바와 같이, 서로 다른 시간에 획득한 같은 환자의 영상에서 질환을 모니터링하기 각각의 영상을 정합하는 방법에 있어서, 종래의 분할기반의 기법과 화소값기반 기법의 장점만을 취합하여 빠르고 정확한 의료영상의 정합방법을 제공하고자 하는 것이다.

이를 위해, 본 발명에 따르면, 후술하는 바와 같이, 먼저, 분할기반의 기법을 이용하여 초기 변환 파라미터를 고려하고, 고려한 초기 변환 파라미터 최적화하기 위해서 화소값기반 기법을 적용하여 자동적으로 최종 파라미터를 구하며, 이때, 화소값기반 최적화는 수행시간이 오래 걸리므로, 3차원 데이터를 빠르게 처리할 수 있도록 GP-GPU 기반의 병렬처리기법을 적용하여 처리속도를 높이는 것을 특징으로 하는 GP-GPU를 이용한 3D 의료영상 정합의 병렬처리방법이 제공된다.

계속해서, 첨부된 도면을 참조하여, 상기한 바와 같은 본 발명에 따른 GP-GPU를 이용한 3D 의료영상 정합의 병렬처리방법의 구체적인 실시예에 대하여 상세히 설명한다.

여기서, 이하의 실시예에서는, CT 스캐너를 이용하여 서로 다른 시간에 얻어진 3D 무릎 이미지의 영상을 정합하는 예를 통하여 본 발명을 설명하나, 본 발명은 이러한 경우로만 한정되는 것은 아니다.

즉, 본 발명의 실시예에 따른 의료영상 정합의 병렬처리방법은, 먼저, 입력 영상인 무릎뼈 영상의 2개의(binary) 뼈 구조의 관성 주축(principal axes)으로부터 산출되는 초기 변환 파라미터(initial transformation parameters)를 추출한다.

다음으로, 다운힐 심플렉스 옵티마이저(downhill simplex optimizer)의 유사 메트릭(similarity metric)을 최소화함으로써(minimizing), 상기한 초기 변환 파라미터를 최적화한다.

여기서, 정합될 이미지가 동일한 형태(modality)로부터 얻어지므로, 유사메트릭으로서 NCC(Normalized-cross correltion)가 적용된다.

또한, 본 발명에 따르면, 종래의 EMP-MI와 달리, 정확한 정합 결과를 제공하기 위해 뼈 영역의 전체 3D 볼륨을 사용하며, 아울러, 스케일링 팩터(scaling factor)는 정합을 위해 계산되는 파라미터를 감소하도록 입력 볼륨의 복셀 사이즈(voxel size)에 근거하여 계산된다.

더욱이, 알고리즘 중 가장 많은 시간이 소비되는 부분인 보간법(interpolation)을 포함하는 변환(transformation) 및 NCC 스코어(score) 계산은, NDIVIA사의 Tesla 이외에, AMD사의 APU(Accelerated Processing Unit), 인텔사의 MIC 아키텍쳐(Many Integrated Core architecture)와 같은 다중 CPU 및 다중 GPU를 동시에 활용 가능한 OpenCL 구현(implementation)을 통하여 병렬처리된다(parallelized).

상기한 변환 파라미터는, 최적화 단계에 이어지는 기준(reference) 및 변환된 볼륨 사이의 복셀 유사성 측정(voxel similarity measurement)을 계산함으로써 얻어진다.

최적화 파라미터는 유사성이 최소 또는 최대가 될 때까지 조사되고, 유사성 측정의 역할은 2개의 이미지가 얼마나 일치하는가를(matched) 나타내는 값을 반환하는(return) 것이다.

최적화 파라미터 행렬(optimal parameter matrix)을 정의하기 위한 비용함수(costfunction)는 다음과 같다.

[수학식 1]

여기서, S는 유사 메트릭(similarity metric)이고, V_r 및 V_f는 각각 기준 볼륨(reference volume) 및 부동 볼륨(float volume)이며, M은 그들 사이의 기하학적 변환행렬(geometric transform matrix)이다.

또한, 변환행렬 M의 형식은 정합 도메인(registration domain)에 따르고, S는 정합 또는 이미지의 특징에 포함되는 형태(modality)의 수에 따라 선택되며, 반복 회수를 줄이기 위해, 초기 파라미터는 관성 주축 기반 접근(principal axes based approach)에 의해 추정된다.

더 상세하게는, 도 1을 참조하면, 도 1은 2개의 영상을 정합하기 위한 처리방법의 전체적인 알고리즘을 나타내는 플로차트이다.

즉, 도 1에 나타낸 바와 같이, 2개의 영상을 정합하기 위한 처리방법은, 초기화 단계 및 최적화 단계의 주요한 두 단계로 구분될 수 있고, 초기화 단계 전에, 이후 단계에서 계산되어야 할 파라미터들의 수를 감소하기 위해 복셀 사이즈에 따른 스케일링이 고려되며, 스케일링 파라미터는 입력 볼륨의 복셀 사이즈에 따라 고려된다.

또한, 초기화 단계는, 바이너리 볼륨 추출(binary volume extraction), 중심 계산(centroid calculation) 및 관성 주축 계산(principal axes calculation)의 세 개의 단계를 포함하여 구성된다.

아울러, 각각의 파라미터들은, 후술하는 바와 같이 하여 두 번째 단계인 최적화 단계에서 최적화된다.

즉, 두 입력 볼륨 사이의 관계를 세 개의 변환 파라미터 및 세 개의 회전 파라미터의 6개의 파라미터에 의해 정의되는 강체변환(rigid body transformation) 행렬로 가정하면, 강체변환행렬 M은, 동일 좌표계(homogeneous coordinates)에서 변환벡터 T(t)와 회전행렬(R)의 연속(concatenation)이라 할 수 있다.

따라서 M은 다음과 같은 형태를 가진다.

[수학식 2]

M = T(t)R

여기서, 회전행렬을 기술하는 일반적인 방법은, 다음의 일련의 행렬들의 분해형(decomposed form) 이다.

[수학식 3]

여기서, α, β, γ는 각각 3D 축에 대한 오일러 각도(Euler angles) 이다.

따라서 회전행렬은, 다음과 같이 나타낼 수 있다.

[수학식 4]

R = R_γ×R_β×R_α

스케일링 팩터 s_x, s_y 및 s_z는, 회전 및 변환 파라미터를 계산하기 전에 별도로 고려되며, 본 발명에 따르면, 스케일링 파라미터는 입력 볼륨의 복셀 사이즈에 따라 고려된다.

복셀 사이즈에 대한 정보는, 입력 이미지의 DICOM(Digital Imaging and Communication in Medicine) 이미지 헤더 파일로부터 얻어지며, x-y 평면에서, 복셀 사이즈는 픽셀 스페이싱(pixel spacing)을 통해 추출되고, z-복셀 사이즈는 입력 DICOM 이미지의 헤더파일 정보의 슬라이스(slice)로부터 얻어진다.

V₁의 복셀 사이즈가 V₂의 복셀 사이즈보다 작은 것으로 가정하면, 스케일링 파라미터는 다음과 같이 계산된다.

[수학식 5]

여기서, V_1i는 볼륨 V₁의 복셀 사이즈이고, V_2i는 볼륨 V₂의 복셀 사이즈이다.

스케일링 파라미터가 계산된 후, V₁은 스케일 다운되고, V₂는 기준 볼륨(reference volume) V_r로 정의되며, 스케일된 볼륨 V₁은 부동 볼륨(float volume) V_f으로 설정된다.

정합 처리 전에 스케일링 파라미터가 최초에 계산되므로, 회전에 대한 3개와 변환에 대한 3개의 6개의 파라미터가 초기화 및 최적화 단계 모두에서 고려되며, 따라서 기준 볼륨 V_r과 부동 볼륨 V_f 사이의 상대 위치(related position)는 변환 파라미터 P = {t_x, t_y, t_z, α, β, γ}의 집합에 의해 정의된다(여기서, t_x, t_y, t_z는 변환량(translation quanta)이고 α, β, γ는 각각 기준 볼륨에 대하여 3D 축에 따른 부동 볼륨의 회전각이다).

계속해서, 기준 볼륨에 대하여 부동 볼륨을 최적으로 정렬하는 최적화 파라미터 집합 P를 구하는 방법에 대하여 설명한다.

6개의 초기 변환 파라미터를 구하기 위해, 정합될 기준 및 부동 CT 볼륨은 모두 2진화되며(binarized), 2진화 볼륨은, 관성 주축(principal axes)을 추출하기 위해 3D의 기하학적 형상(geometric shape)을 나타내기 위한 특징으로서 사용된다.

B(x, y, z)를 3D 볼륨 V(x, y, z)의 초기 2진화 볼륨이라 하면, 다음과 같이 나타낼 수 있다.

[수학식 6]

여기서, x, y, z는 이미지의 복셀의 좌표(coordinates)이고, τ는 2진화 볼륨을 정의하는 임계값(threshold value)이다.

3D 라벨링 처리는 초기 임계값을 추출한 후에 이어진다. 2진화 이미지의 연결된 성분(connected components)에 대한 라벨링은 2진화 이미지를 심볼릭 이미지로 변환하여 각각의 연결된 성분에 고유한(unique) 라벨을 할당하도록 한다.

라벨링 후, 이미지 필링(filling) 및 이로전(erosion)과 같은 형태학적 연산(morphological operations)이 이미지의 구멍(hole)을 채우고 이미지의 원하지 않는 작은 영역을 제거하기 위해 적용된다.

6개의 파라미터의 계산은, 2진 객체(binary objects)의 관성 주축을 이용하여 달성되며, 3D 객체의 관성 주축은 3×3 관성행렬(inertia matrix)의 고유벡터(eigenvectors)에 의해 정의된다.

여기서, 상기한 관성행렬은 다음과 같이 정의된다.

[수학식 7]

여기서,

는 객체 모멘트(object moments) 이고, 함수 f(x, y, z)는 복셀 데이터의 이미지 내용(image content)을 나타내며, x_c, y_c, z_c는 객체의 중심을 나타내며, 관성행렬로부터 계산된 3개의 고유벡터는 직접적으로 객체의 관성 주축을 나타낸다.

고유벡터의 행렬 형태는 다음과 같이 나타낼 수 있다.

[수학식 8]

회전행렬 R에 대하여 행렬 E를 푸는 것으로(즉, E = R), 회전각 α, β, γ가 다음과 같이 계산될 수 있다.

[수학식 9]

또한, 초기 파라미터를 추출하는 단계는, 먼저, 2진화된 볼륨의 픽셀의 좌표로부터 3차원 벡터가 형성되면, 이어서, 3D 좌표 벡터로부터 중심 및 관성행렬이 산출되며, 최종적으로, 각각의 관성행렬의 고유벡터로부터 각 볼륨의 회전각이 산출된다.

기준 볼륨의 회전각 x, y, z로부터 부동 볼륨의 회전각 x, y, z를 감산함으로써(subtracting) 3개의 초기 회전각이 산출되며, 마찬가지로, 기준 볼륨의 중심 x, y, z로부터 부동 볼륨의 중심 x, y, z를 감산함으로써 3개의 초기 변환 파라미터가 산출된다.

계속해서, 일반적으로 널리 이용되는 다운힐 심플렉스 방법(downhill simplex method)을 적용하여 파라미터를 최적화하는 방법의 상세한 내용에 대하여 설명한다.

즉, 상기한 방법은, 최적화 처리를 N+1 포인트로 초기화하고, N차원(N-dimensional) 파라미터 공간의 초기 심플렉스(simplex)를 정의한다.

여기서, 상기한 심플렉스는 반복적으로 변형되며(deformed), 즉, 상기한 심플렉스는, 비용함수(cost function) S의 최소값을 향하여 심플렉스의 정점들(vertices)이 천이하도록(shift), 반복되는 단계에서 반사(reflection), 확장(expansion) 또는 단축(contraction)에 의해 형성된다.

또한, 다운힐 심플렉스 방법은, 예를 들면, "Convergence properties of the Nelder-Mead simplex method in low dimensions", J. C. Lagarias, J. A. Reeds, M. H. Wright and P. E. Wright, SIAM Journal on Optimization, 9(1999) 112-147에 개시된 바와 같은 종래의 알고리즘에 근거하여 구현 가능하며, 일반적으로, 심플렉스는 제로값(zero value)으로 초기화되고 유사성 스코어(similarity score)가 최소값이 될 때까지 최적화된다.

본 발명의 실시예에서는, 정합 시간을 알고리즘적으로 개선하기 위해 기존의 초기화 단계에 근거하여 얻어진 6개의 변환 파라미터로 상기한 심플렉스가 초기화되는 것으로 하여 본 발명을 설명한다.

먼저, 변환 과정에 대하여 설명하면, 최적화의 반복시에 강체 변환(rigid transformation)이 존재하게 되면, 처리는 다음과 같이 수행될 수 있다.

(1) 기준 볼륨과 동일한 기하학적 형상(geometry)을 가지는 빈 이미지 볼륨(empty image volume)(임시 변환 볼륨(temporary transform volume))을 생성.

(2) 빈 이미지 볼륨의 각 복셀에 대하여,

1) 수학식 M = T(t)R에 기재된 바와 같은 변환 행렬을 이용하여, 부동 볼륨의 대응하는 복셀 포인트를 계산.

2) 정규 그리드(regular grid) 상에 계산된 대응하는 복셀의 3-선형 삽입(tri-linear interpolation)을 연산(operate).

3) 삽입된 복셀을 빈 이미지 볼륨에 설정(set).

(3) 모든 복셀에 대하여 상기한 단계들을 반복.

다음으로, 유사성 측정(similarity measure)에 대하여 설명하면, 다음과 같다.

즉, 정합의 인증(validation)은, 두 볼륨 사이의 상관 계수(correlation coefficient)에 의해 평가되며, 따라서 최적화 단계에서 볼륨의 유사성의 정도(degree)를 정량화(quantify) 하기 위해 NCC(normalized cross-correlation)가 적용된다.

여기서, NCC 함수는 다음과 같이 정의된다.

[수학식 10]

여기서, n은 복셀의 총 수(total number), i와 j는 복셀 인덱스(voxel index), x가 V_r(x_i) 및 V_f(x_j)의 x-, y-, z- 좌표를 나타내고, V_r과 V_f가 평균 화소값(mean intensity value)을 나타낼 때, 포인트 x_i와 x_j에서의 화소값 이다.

아울러, 최종 변환 파라미터는 V_r과 V_f 사이에서 NCC 스코어가 최대일 때 최적인 것으로 간주된다.

계속해서, 도 2를 참조하여, GPU 구현에 대하여 설명한다.

즉, 도 2를 참조하면, 도 2는 본 발명의 실시예에 따른 GP-GPU를 이용한 3D 의료영상 정합의 병렬처리방법의 전체적인 구성을 개략적으로 나타내는 도면이다.

더 상세하게는, 도 2에 나타낸 바와 같이, 본 발명의 실시예에 따른 GP-GPU를 이용한 3D 의료영상 정합의 병렬처리방법은, 먼저, 두 개의 영상을 입력(source) 영상과 목표(target) 영상으로 각각 입력받고, 각 영상의 스케일링 파라미터를 찾는다.

이어서, 파라미터의 초기화를 위해, 도 2에 나타낸 바와 같이, 각 영상의 바이너리 볼륨을 추출하고, 중심(centroid) 및 관성 주축(principal axes)을 계산하며, 두 영상과 계산된 초기 파라미터를 "MEX" 인터페이스를 통해 CPU 및 GPU로 각각 보내고, CPU와 GPU는 각각 C 언어 및 OpenCL을 이용하여 최적화 작업을 병렬 처리한다.

상기한 바와 같이, 파라미터의 최적화 과정은 CPU와 GPU에 의해 각각 병렬처리되며, 즉, CPU는 C 언어로 작성된 프로그램에 의해 제어되고, GPU는 OpenCL을 이용하여 제어된다.

또한, 도 2에 있어서, 3차원 강체변환(rigid transformation)은 입력 영상들로부터 얻어지는 좌표계로 정의되는 임의의 점 P(x,y,z)를 절대좌표계에서 정의되는 점 P'(x,y,z)로 변환하는 것을 의미하며, 다음과 같이 나타낼 수 있다.

[수학식 11]

P' = RP+t

여기서, R은 3×3 회전행렬(rotation matrix)이고, t는 3×1 이동벡터(translation vector) 이며, R은 x, y, z 좌표계의 세 회전행렬을 곱하여 얻은 행렬이다.

또한, 회전행렬은 다음과 같이 나타낼 수 있다.

[수학식 12]

R = R_γ×R_β×R_α

여기서,

이고, t 벡터는 t=[t _x t _y t _z ] ^T 로 표시하며, 스케일링 변환 행렬은 아래와 같다.

[수학식 13]

도 2에 나타낸 바와 같이, 강체변환 파라미터 P={tx, ty, tz, α, β, γ}를 찾기 전에, 이하의 식을 이용하여, 입력된 두 이미지의 복셀(voxel) 정보를 이용하여 3개의 스케일링 파라미터를 계산한다.

[수학식 14]

여기서 v_1i는 v_2i보다 작으며, v_1i는 첫 번째 입력된 이미지의 복셀정보이고, v_2i 는 두 번째 입력된 이미지의 복셀정보이다.

상기한 바와 같이 두 이미지의 크기를 맞추고 나서 강체변환의 6개의 파라미터를 고려하기 위해 다음과 같은 다섯 개의 단계을 수행한다.

1. 입력된 완본 이미지 두 개에서 뼈 부분을 분할하는 단계

2. 분할한 바이너리 이미지의 각 중심을 찾는 단계

3. 바이너리 이미지의 관성 행렬(inertia matrix)을 통해 고유벡터 (eigenvectors)를 고려하고, 각 고유벡터에서 각각 회전 파라미터를 계산하는 단계

4. 첫 번째 입력된 이미지의 회전 파라미터에서 두 번째 입력된 이미지의 회전 파라미터를 감산하고, 두 이미지의 회전 파라미터를 고려하는 단계

5. 두 이미지의 중심을 감산하고 변환 파라미터를 고려하는 단계

상기한 바와 같이 하여 6개의 초기 파라미터를 관성 주축 기반으로 찾고 나서, 그 파라미터들을 화소값 기반으로 최적화한다.

여기서, 최적화의 방법으로서는, 상기한 바와 같이 다운힐 심플렉스(Downhill simplex) 방법을 적용하고, 유사성 측정 메트릭은 NCC(normalized cross correlation )를 적용한다.

또한, 수행시간이 오래 걸리는 변환 계산과 메트릭 계산은 그래픽 카드의 GPU를 이용하여 처리함으로써, CPU와 GPU의 병렬처리 기법을 사용하여 처리속도를 높일 수 있다.

다음으로, 상기한 바와 같은 OpenCL을 이용한 병렬처리의 구현방법에 대하여 상세히 설명한다.

즉, 초기화 단계에 의해 정합 시간을 알고리즘적으로 감소시킨다 하더라도, 삽입 및 유사성 스코어 계산을 포함하는 대량의 반복적인 최적화 처리로 인해 정합 처리의 속도 문제는 여전히 남아 있게 된다.

더 상세하게는, 반복을 위한 NCC 스코어를 얻기 위해, 가산(addition), 감산(subtraction), 승산(multiplication), 제산(division), 평균 제곱(mean square), 제곱근 연산(square root operation)과 같은 다수의 수학적 연산이 각 픽셀에 적용되고, 주어진 변환행렬에 의한 하나의 좌표계로부터 다른 좌표계로 볼륨의 변환은 볼륨 내의 모든 복셀에 걸쳐 행해져야(traverse) 하며, 각 복셀의 좌표는 변환행렬과 곱해지고, 화소값(intensity value)은 삽입 후에 새로운 볼륨으로 보내진다.

이러한 연산은 복잡하지는 않으나, 볼륨 내의 모든 복셀을 거쳐야 하고, 일반적으로 3D 의료영상은 큰 사이즈의 볼륨을 가지므로, 직렬 프로그래밍에서는 매우 시간이 걸린다.

NCC 스코어 계산과 변환의 문제는, SIMD 문제로서 나타내질 수 있고, OpenCL은 이러한 종류의 문제를 다루기 위해 적용가능한 병렬 프로그래밍 모델을 제공한다.

또한, OpenCL의 상세한 내용에 대하여는, 예를 들면, "OpenCL Programming Guide", A. Munshi, B. Gaster, T. G. Mattson, J. Fung and D. Ginsburg, Addison-Wesley Professional, 2011.에 게시된 내용을 참조할 수 있으며, 즉, OpenCL은 다양한 형태의 연산장치의 사용에 관한 이종 연산모델(heterogeneous computing model)이며, GPU 및 다른 프로세서들이 CPU와 직렬로(tandem) 동작하도록 설계된 개방적이고 자유로운 API(Application Programming Interface) 이다.

OpenCL 구조는, 플랫폼 모델(platform model), 실행 모델(execution model), 프로그래밍 모델(programming model) 및 메모리 모델(memory model)의 4가지 모델을 이용하여 기술될 수 있다.

더 상세하게는, 도 3을 참조하면, 도 3은 OpenCL의 플랫폼 모델의 구성을 개략적으로 나타내는 도면이다.

즉, 도 3에 나타낸 바와 같이, 플랫폼 모델은, 호스트 시스템(CPU)과 다양한 OpenCL 구성요소 사이의 전체적인 관계를 기술하는 것이며, 여기서, 호스트는 하나 이상의 다른 연산장치와 연결될 수 있다.

또한, 각각의 연산장치는 연산 유닛(코어)의 집합(collection)이며, 각각의 연산 유닛은 처리소자(processing elements)(PEs)로 구성되고, 처리소자는 프로세서의 적절한 사용과 SIMD(Single Instruction, Multiple Data) 또는 SPMD(Single Programm, Multiple Data)와 같은 코드의 실행을 수행한다.

다음으로, OpenCL의 실행 모델은, 장치의 커널 실행(kernel execution) 및 호스트의 순차 프로그램(sequential program) 또는 호스트 프로그램(host program) 실행의 두 부분으로 구성된다.

여기서, 커널은 컴퓨터 장치에서 실행되는 실행가능한 코드의 기본 단위(basic unit)이며 C의 함수와 유사하고, 각각의 커널 실행은 작업 아이템(work-item)이라 불리며 작업 아이템의 그룹은 작업 그룹(work-groups)으로 나타낸다.

아울러, 실행을 위해 커널이 선택되면, "NDRange"라 불리는 인덱스 공간(index space)이 정의되며, 이때, NDRange는 N이 1, 2 또는 3일 때 N차원 인덱스 공간을 의미한다.

즉, 도 4를 참조하면, 도 4는 OpenCL의 NDRange 구성을 개략적으로 나타내는 도면이다.

도 4에 나타낸 바와 같이, 작업 그룹 크기의 특정은 또한 NDRange 구성(NDRange configuration)이라고도 불리며, 동기화(synchronization)는 오직 동일한 작업 그룹 내의 작업 아이템 사이에서만 허용되고, 다른 작업 그룹의 작업 아이템과는 동기화가 불가능하다.

또한, 실행 모델의 호스트 프로그램은 디바이스 콘텍스트(device context)를 정의하는 호스트 시스템에서 동작하고, 명령 큐(command queue)를 이용하는 커널 실행을 정렬(queue) 한다.

아울러, 프로그래밍 모델은 작업 기반(task-based) 또는 데이터 기반(data-based)이 될 수 있으며, 작업 기반 병렬 컴퓨팅은 작업 그룹이 다른 모든 작업 그룹에 대하여 독립적으로 실행되는 반면, 데이터 기반 병렬 컴퓨팅은 장치 커널의 다중 인스턴스(multiple instance)에 의해 병렬 처리된다.

계속해서, 도 5를 참조하여, OpenCL의 메모리 모델의 구성에 대하여 설명한다.

도 5에 나타낸 바와 같이, OpenCL은 프라이빗 메모리(private memory)로부터 글로벌 메모리(global memory)까지 메모리 레인징(memory ranging)을 통하여 메모리 레벨을 정의하며, 즉, OpenCL은 프라이빗(private), 로컬(local), 상수(constant) 및 글로벌(global) 메모리의 4개의 메모리 공간(memory space)을 정의한다.

더 상세하게는, 도 5를 참조하면, 도 5는 OpenCL의 메모리 모델의 구성을 개략적으로 나타내는 도면으로, OpenCL에 의해 정의되는 메모리 계층 구조(memory hierarchy)의 다이어그램을 나타내고 있다.

도 5에 있어서, 프라이빗 메모리는 단일 계산유닛(single compute unit)에 의해서만 사용될 수 있는 메모리이고, 싱글 코어 CPU의 레지스터와 유사하다.

또한, 로컬 메모리는 작업 그룹 내의 작업 아이템에 의해 사용될 수 있는 메모리이고, 상수 메모리는 커널의 실행 동안 장치 내의 모든 계산유닛에 의해 읽기 전용으로(read only access) 상수 데이터(constant data)를 저장하기 위한 메모리이며, 글로벌 메모리는 장치 내의 모든 계산 유닛에 의해 사용될 수 있다.

여기서, 상기한 메모리들을 충분히 활용하기 위하여는, 각 종류의 메모리의 특성이 알고리즘의 특징 및 데이터 구조를 만족해야 한다.

계속해서 OpenCL의 구현방법의 상세한 내용에 대하여 설명한다.

GPU가 데이터를 처리하기 전에, 해야 할 중요한 작업은 CPU로부터 GPU로 데이터의 할당(allocation) 이다.

먼저, 입력 볼륨의 3차원 데이터가 효율적인 계산을 위해 1차원 배열로 형성되고, 볼륨 데이터의 사이즈가 크고 또한 모든 작업 아이템에 의해 엑세스 가능해야 하므로, 그러한 데이터는 호스트로부터 글로벌 메모리에 복사된다(copied).

또한, 변환 행렬은, 국지적으로(locally)만 필요하고 부동형의 4×4 사이즈를 가지므로 로컬 메모리에 복사되고, 볼륨의 폭(width), 높이(height) 및 깊이(depth) 정보와 같은 다른 상수 변수(constant vatiables)는 프라이빗 메모리에 복사된다.

아울러, NCC 스코어의 계산은, 주로 각 복셀의 합(summation)에 대한 루프(loop)를 포함하고, GPU 메모리의 코어의 각 복셀값의 합을 위해, 예를 들면, "http://www.nvidia.com/content/cudazone/download/OpenCL/NVIDIA_OpenCL_ProgrammingOverview.pdf"에 제시된 바와 같은, 병렬 감소 방법(parallel reduction stratege)이 사용된다.

또한, 작업 아이템의 합은 작업 그룹 내에서 완료되며, 결과는 작업그룹 IDs에 저장되고 호스트에 전달되며, 최종 합산은 루프의 수가 작업 그룹의 사이즈와 같아질 때까지 반복되어 호스트 측에서 완료된다.

상기한 NCC 함수의 식에 나타낸 바와 같이, 해당 식은 각 볼륨의 평균 및 제곱 평균의 계산을 포함하고 있으며, 따라서 합산을 위한 감소 커널(reduction kenel)과 분모(denominator) 및 분자(nominator)의 분할 계산(partial calculation)을 위한 평균 제곱 커널(mean square kernel)의 2개의 커널이 NCC 스코어를 계산하기 위해 분할된다.

이어서, 최종 계산은, 결과의 합산을 수행해야 하는 호스트 프로그램 측에서 완료된다.

또한, 상기한 바와 같은 병렬 변환을 위한 OpenCL 구현의 설계는 도 6에 나타낸 바와 같다.

즉, 도 6을 참조하면, 도 6은 3D 변환을 위한 OpenCL 구형의 설계방법의 개념을 개략적으로 나타내는 도면이다.

도 6에 나타낸 바와 같이, 병렬화되어야 할 픽셀의 총 수는, 할당되어야 할 글로벌 메모리의 사이즈가 폭×높이×깊이의 총계(amount)와 같도록, 이미지의 폭, 높이 및 깊이의 곱과 같다.

작업 그룹 사이즈는, 이미지 볼륨의 높이와 동일한 사이즈로 정의되고, 512의 로컬 사이즈가 변환 및 NCC 스코어 계산의 병렬화(parallelization)를 위해 사용된다.

여기서, 병렬화 처리의 단계는 다음과 같다.

1. 호스트로부터 장치로 메모리 복사(볼륨 데이터, 변환 파리미터, 상수값)

2. 최적화기에서 각각 반복

(1) 변환 커널 런칭

1) 글로벌 인덱스를 복셀 좌표로 변환

2) 변환 행렬에 복셀 좌표를 곱함

3) 3선형 삽입 후 부동 볼륨의 이전 좌표로부터 임시 볼륨의 새 좌표로 화소값을 복사

(2) 감소 커널 런칭

4) 임시 볼륨의 화소값을 합산하고, 각 작업 그룹에 대한 결과를 각각의 작업 그룹 IDs에 저장

5) 장치로부터 호스트로 결과를 복사하고 평균값의 계산은 호스트측에서 완료

(3) 평균 제곱 커널 런칭

6) 기준 및 변환 볼륨의 평균 제곱을 계산

7) 각 볼륨에 대한 평균 제곱의 합을 계산하기 위해 감소를 계산

8) 장치로부터 호스트로 관련된 값을 복사

(4) 호스트에서 NCC 스코어 계산을 완료, 여기서, 작업 그룹, 제곱, 제곱근 및 제산의 추가적인 합을 계산해야 함.

3. NCC 스코어를 최적화기에 반환

여기서, 상기한 (2)-4) 단계에서 임시 볼륨의 화소값을 합산하는 방법 또한, "http://www.nvidia.com/content/cudazone/download/OpenCL/NVIDIA_OpenCL_ProgrammingOverview.pdf"에 게시된 바와 같은 감소 방법을 이용할 수 있다.

따라서 상기한 바와 같이, 본 발명의 실시예에 따르면, 3차원 관성주축 기반 방법과 화소값 기반 방법의 장점만을 결합하여 보다 빠르고 정확하며 효과적인 3차원 영상의 정합방법을 제공할 수 있다.

이상, 상기한 바와 같은 본 발명의 실시예를 통하여 본 발명에 따른 GP-GPU를 이용한 3D 의료영상 정합의 병렬처리방법의 상세한 내용에 대하여 설명하였으나, 본 발명은 상기한 실시예에 기재된 내용으로만 한정되는 것은 아니며, 따라서 본 발명은, 본 발명이 속하는 기술분야에서 통상의 지식을 가진 자에 의해 설계상의 필요 및 기타 다양한 요인에 따라 여러 가지 수정, 변경, 결합 및 대체 등이 가능한 것임은 당연한 일이라 하겠다.

Claims (14)

1. 삭제
2. 서로 다른 시간에 획득한 동일한 환자의 2개의 영상을 CPU(Central Processing Unit)와 GPU(Graphics Processing Unit)에 의해 병렬처리하여 정합하는 영상 정합의 병렬처리방법에 있어서,
   서로 다른 시간에 획득된 동일한 환자의 2개의 영상을 각각 입력받는 입력 단계;
   상기 입력받는 단계에서 입력된 영상들의 복셀 사이즈(voxel size)에 대한 정보를 입력 영상들의 DICOM(Digital Imaging and Communication in Medicine) 이미지 헤더 파일로부터 얻고, 상기 복셀 사이즈를 x-y 평면에서 픽셀 스페이싱(pixel spacing)을 통해 추출하며, z-평면에서는 입력 DICOM 영상의 헤더 파일 정보의 슬라이스(slice)로부터 얻고, 제 1 입력영상 V₁의 복셀 사이즈가 제 2 입력영상 V₂의 복셀 사이즈보다 작을 경우, 이하의 수학식을 기초로 산출된 스케일링 파라미터를 이용하여 스케일링(scaling) 하는 스케일링 단계;
   상기 V₁을 스케일 다운하고, 상기 V₂을 기준 볼륨(reference volume)(V_r)으로 설정하며, 상기 V₁을 부동 볼륨(float volume)(V_f)으로 설정하는 단계;
   상기 스케일링된 영상들로부터 초기 변환 파라미터(initial transformation parameter)를 추출하는 초기화 단계(initialization step);
   상기 초기화 단계에서 추출된 상기 초기 변환 파라미터들을 상기 부동 볼륨 및 상기 기준 볼륨 사이의 유사성 스코어가 최대가 되도록 하는 최종 파라미터를 산출하는 최적화 단계(optimization step); 및
   상기 최적화 단계에서 산출된 상기 최종 파라미터를 이용하여 상기 입력된 영상들을 정합하고 표시장치를 통하여 표시하는 시각화 단계(visualization step)를 포함하는 것을 특징으로 하는 영상 정합의 병렬처리방법:
   

   

   
   여기서, V_1i는 볼륨 V₁의 복셀 사이즈이고, V_2i는 볼륨 V₂의 복셀 사이즈이다.
   
3. 삭제
4. 삭제
5. 제 2항에 있어서,
   상기 CPU는 C 언어로 작성된 프로그램에 의해 제어되고, 상기 GPU는 OpenCL 또는 CUDA(Compute Unified Device Architecture)를 이용하여 제어되도록 구성되는 것을 특징으로 하는 영상 정합의 병렬처리방법.
   
6. 제 2항에 있어서,
   상기 초기화 단계에서, 상기 초기 변환 파라미터는, 3개의 회전 파라미터와 3개의 변환 파라미터를 포함하며,
   상기 기준 볼륨 V_r과 상기 부동 볼륨 V_f 사이의 상대 위치(related position)는, 변환 파라미터 P = {t_x, t_y, t_z, α, β, γ}(여기서, t_x, t_y, t_z는 변환량(translation quanta)이고 α, β, γ는 각각 기준 볼륨에 대하여 3D 축에 따른 부동 볼륨의 회전각)의 집합에 의해 정의되는 것을 특징으로 하는 영상 정합의 병렬처리방법.
   
7. 제 6항에 있어서,
   상기 초기화 단계는,
   상기 기준 볼륨 및 상기 부동 볼륨이 모두 2진화되는(binarized) 단계;
   2진화된 상기 기준 볼륨 및 상기 부동 볼륨의 픽셀의 좌표로부터 3D 벡터를 형성되고, 형성된 상기 3D 벡터로부터 중심(centroiod) 및 관성행렬(inertia matrix)이 산출되는 단계;
   각각의 상기 관성행렬의 고유벡터로부터 각 볼륨의 회전각이 산출되는 단계;
   상기 기준 볼륨의 회전각(x, y, z)으로부터 상기 부동 볼륨의 회전각(x, y, z)을 감산함으로써(subtracting) 3개의 초기 회전 파라미터가 산출되는 단계; 및
   상기 기준 볼륨의 중심(x, y, z)으로부터 상기 부동 볼륨의 중심(x, y, z)을 감산함으로써 3개의 초기 변환 파라미터가 산출되는 단계를 포함하는 것을 특징으로 하는 영상 정합의 병렬처리방법.
   
8. 제 7항에 있어서,
   상기 2진화되는 단계는,
   B(x, y, z)를 3D 볼륨 V(x, y, z)의 초기 2진화 볼륨이라 하면, 이하의 수학식을 이용하여 상기 기준 볼륨 및 상기 부동 볼륨이 2진화되는 것을 특징으로 하는 영상 정합의 병렬처리 방법:
   
   

   

   
   
   여기서, x, y, z는 이미지의 복셀의 좌표(coordinates)이고, τ는 2진화 볼륨을 정의하는 임계값(threshold value)이다.
   
9. 제 8항에 있어서,
   상기 관성행렬은, 이하의 수학식으로 정의되며,
   상기 관성행렬의 고유벡터(eigenvectors)에 의해 관성 주축(principal axes)이 정의되는 것을 특징으로 하는 영상 정합의 병렬처리방법:
   
   

   

   
   
   여기서,

   

   는 객체 모멘트(object moments) 이고, 함수 f(x, y, z)는 복셀 데이터의 이미지 내용(image content)을 나타내며, x_c, y_c, z_c는 객체의 중심을 나타내고, 상기 관성행렬로부터 계산된 3개의 고유벡터는 상기 객체의 관성 주축을 나타낸다.
   
10. 제 9항에 있어서,
    상기 고유벡터의 행렬 형태는 이하의 수학식으로 나타내지며,
    
    

    

    
    
    회전행렬 R에 대하여 상기 행렬 E의 연산에 의해(E = R), 상기 회전각 α, β, γ가 이하의 수학식으로 계산되고,
    
    

    

    
    
    상기 회전행렬은, 이하의 수학식으로 나타내지는 것을 특징으로 하는 영상 정합의 병렬처리방법:
    
    R = R_γ×R_β×R_α
    
    여기서, α, β, γ는 각각 3D 축에 대한 오일러 각도(Euler angles) 이며,
    

    

    
    이다.
    
11. 제 10항에 있어서,
    상기 최적화 단계는,
    상기 부동 볼륨의 각 복셀에 대하여, 강체 변환 행렬(rigid body transfomation matrix) 및 3차-선형 보간법(tri-linear interpolation) 연산(operate)을 이용하여 상기 부동 볼륨을 변환하는 단계;
    변환된 상기 부동 볼륨과 상기 기준 볼륨 사이의 유사성 스코어(similarity score)를 측정하는 단계; 및
    상기 유사성 스코어가 최대가 될 때까지 모든 복셀에 대하여 상기 부동 볼륨을 변환하는 단계 및 상기 유사성 스코어를 측정하는 단계를 반복하는 단계를 포함하는 것을 특징으로 하는 영상 정합의 병렬처리방법.
    
12. 제 11항에 있어서,
    상기 최적화 단계는,
    두 입력 영상 사이의 관계를 3개의 변환 파라미터 및 3개의 회전 파라미터에 의해 정의되는 강체변환(rigid body transformation) 행렬 M이라 하면, 상기 강체변환행렬 M은 이하의 수학식으로 나타내지는 것을 특징으로 하는 영상 정합의 병렬처리방법:
    
    M = T(t)R
    
    여기서, 상기 T(t) 및 상기 (R)은 동일 좌표계(homogeneous coordinates)에서의 변환벡터와 회전행렬이다.
    
13. 제 12항에 있어서,
    상기 유사성 스코어의 측정은, 상기 최적화 단계에서 두 볼륨의 유사성의 정도(degree)를 정량화(quantify) 하기 위해 이하의 NCC(normalized cross-correlation) 함수를 이용하여 수행되는 것을 특징으로 하는 영상 정합의 병렬처리방법:
    
    

    

    
    
    여기서, i와 j는, x가 V_r(x_i) 및 V_f(x_j)의 x-, y-, z- 좌표를 나타내고 V_r과 V_f가 평균 화소값(mean intensity value)을 나타낼 때, 포인트 x_i와 x_j에서의 화소값을 나타내는 복셀 인덱스(voxel index) 이다.
    
14. 제 13항에 있어서,
    상기 최적화 단계에서, 상기 부동 볼륨을 변환하는 단계 및 상기 유사성 스코어를 측정하는 단계의 처리가 상기 GPU에 의해 병렬처리되는 것을 특징으로 하는 영상 정합의 병렬처리방법.
    

Images