CT 영상에서 골다공증 판별 방법의 성능 향상

A Performance Enhancement of Osteoporosis Classification in CT images

Abstract

Classification methods based on dual energy X-ray absorptiometry, ultrasonic waves, and quantitative computed tomography have been proposed. Also, a classification method based on machine learning with bone mineral density and structural indicators extracted from the CT images has been proposed. We propose a method which enhances the performance of existing classification method based on bone mineral density and structural indicators by extending structural indicators and using principal component analysis. Experimental result shows that the proposed method in this paper improves the correctness of osteoporosis classification 2.8% with extended structural indicators only and 4.8% with both extended structural indicators and principal component analysis. In addition, this paper proposes a method of automatic phantom analysis needed to convert the CT values to BMD values. While existing method requires manual operation to mark the bone region within the phantom, the proposed method detects the bone region automatically by detecting circles in the CT image. The proposed method and the existing method gave the same conversion formula for converting CT value to bone mineral density.

1. 서 론

고령화 사회가 되어 감에 따라 골다공증은 노인들의 중요한 건강 문제로 받아들여지고 있다. 이로 인하여 골다공증의 진단을 위한 방법들로 이중에너지 X선 흡수 계측법(Dual energy X-ray Absorptiometry; DXA), 초음파를 이용한 방법, 정량적 컴퓨터 단층촬영(quantitative computed tomography)을 이용한 방법들이 개발되었다.

이중에너지 X선 흡수 계측법은 요추와 대퇴골 2군데의 2차원 영상에서 단위 면적당 무게로 표현되는 골밀도를 계산한다[1-3]. 골밀도 값으로부터 T 값을 추출하는데, T 값은 20~30대 동일 성별 정상인의 골밀도와 차이를 표준편차에 의해 수치화 값이다. 세계 보건 기구에서는 T 값이 -2.0이하이면 골연화증으로 -2.5 이하이면 골다공증으로 진단하도록 권고하고 있다[4]. 이중에너지 X선 흡수 계측법이 2차원 영상을 이용하는 것에 반해, 정량적 컴퓨터단층촬영(quantitative computed tomography) 방법은 3차원 영상에서 골밀도를 측정한다[5-7]. 이 방법에서는 주로 요추(lumbar spine)의 해면뼈 영역에 대하여 단위 부피당 무게로 표현되는 골밀도를 계산한다. 이 방법은 주위 조직의 영향을 받지 않고 골밀도를 검사할 수 있으므로 이중에너지 X선 흡수 계측법에 비하여 보다 정확할 수 있다. 초음파를 이용한 방법도 사용되고 있는데, 이 방법은 장비의 크기가 작아 공간을 적게 차지하여 간편하게 이용할 수 있으며 방사선 노출이 없는 장점이 있다[8-10].

위의 방법들이 골밀도를 이용하여 골다공증을 진단하고 있는데, CT 영상에서 추출된 골밀도와 함께 다른 구조적 지표자를 함께 사용함으로써 골다공증을 진단하는 방법이 참고문헌 [11]에서 제안되었다. 이 방법에서는 요골 원위부(distal radius)를 촬영하여 골밀도와 함께 6개의 구조적 지표자들을 검출하였으며, 기계학습 방법을 이용하여 골다공증을 진단하였다.

본 논문에서는 지표자의 개수를 16개로 확장하고 주성분 분석을 통하여 데이터를 가공하는 방법을 이용하여 골다공증의 인식률을 향상시키는 방법을 제안한다. 참고문헌 [11]에서는 SAD(Sum of Absolute Difference), SVM(Support Vector Machine), ANN(Aritifical Neural Network) 세 가지 방법을 이용하여 골다공증을 진단하는 실험을 했는데, 골밀도와 지표자를 함께 사용하였을 경우에 평균적으로 SAD, SVM, ANN 각각 88.0%, 86.5%, 79.9%의 인식률을 보였다. 본 논문에서는 실험 결과 지표자를 16개로 확장함으로써 평균적으로 SAD, SVM, ANN 각각 89.2%, 87.4%, 84.7%의 인식률을 얻을 수 있었다. 또한 주성분 분석을 이용함으로써 SAD, SVM, ANN 각각 90.9%, 89.5%, 86.0%의 인식률을 얻을 수 있었다.

인식률 향상 방법에 더불어 본 논문에서는 골밀도 계산을 위한 팬텀 분석 작업에 있어서, 참고문헌[11]에서 사용하였던 수작업을 자동으로 수행하는 방법을 제안한다. 참고문헌 [11]과 본 논문에서 사용한 CT는 정량적 CT가 아니므로 골밀도 값이 제공되지 않는다. CT 영상의 값을 골밀도 값으로 변환하기 위해 팬텀을 이용하는데, 팬텀의 해면뼈에 해당하는 부분을 분석하는 것이 필요하다. 팬텀에는 6개의 뼈 영역이 있는데, 각 영역을 분석하기 위해 참고문헌 [11]에서는 각 영역을 수작업으로 지정하였다. 본 논문에서는 영상에서 원을 검출하는 방법을 이용하여 자동으로 6개의 영역을 찾는 방법을 제안한다. 제안된 방법으로 구한 각 영역의 평균 CT 값은 수작업으로 영역을 지정하여 구한 평균 CT 값과 99.4% 이상 일치하였다.

 

2. 골밀도 추출

본 논문에서 사용한 CT는 정량적 CT가 아니므로 복셀의 값을 HU(Housefield Unit) 단위로 생성한다. HU 값을 골밀도 값으로 변환하기 위해서 QRM사의 팬텀을 사용하였다[12]. 이 팬텀은 팔뚝 관절에 대해서 HU 값을 골밀도 값으로 변환하는데 사용될 수 있도록 만들어진 팬텀으로서, Fig. 1의 단면도에 나타나 있는 바와 같이 6개의 영역(S1~S6)으로 구성된다. S1, S2, S3는 대, 중, 소 크기의 뼈를 모방한 것으로서, 각각 28mm, 21mm, 14mm 크기의 지름을 갖는 원기둥 형태의 구조물이다. 이들의 칼슘 수산화인회석(hydroxyapatite; HA) 밀도는 200, 100, 50 mg HA/cm3 로 만들어져 있다. 이들의 가장자리 부분은 치밀뼈에 해당하는, 밀도가 800 mg HA/cm3인 물체로 채워져 있는데, 두께가 S1과 S3는 1.2mm이고 S2는 0.6mm이다. S4, S5, S6은 각각 10mm, 14mm, 10mm 크기의 지름을 갖는 원기둥 형태의 구조물이다. 내부는 물과 같은 밀도를 가지는 레진으로 채워져 있으며, 가장자리는 두께가 각각 2.5mm, 1mm, 2.5mm인 치밀뼈에 해당하는 물질로 채워져 있다.

Fig. 1.Cross-sectional diagram of the phantom.

팬텀을 CT로 촬영한 단면도 영상이 Fig. 2에 나타나 있다. Fig. 2(a)의 왼쪽에는 영역 S1이 오른쪽에는 영역 S5가 나타나 있고 Fig. 2(b)의 왼쪽에는 영역 S3이 오른쪽에는 영역 S5가 나타나 있다. Fig. 2(c)의 왼쪽에는 영역 S3이 오른쪽에는 영역 S5가 나타나 있고 Fig. 2(d)의 왼쪽에는 영역 S4가 오른쪽에는 영역 S6이 나타나 있다.

Fig. 2.Slice image of the phantom (a) Section S1 and S4 (b) Section S2 and S4 (c) Section S3 and S4 (b) Section S4 and S6.

HU 값과 골밀도 값과의 관계를 얻기 위해서 S1~S6 각각의 해면뼈 영역과 치밀뼈 영역의 평균 HU 값을 구한다. 이를 위해 먼저 각 영역을 분할하고 평균 HU 값을 계산한다. 참고문헌 [11]에서는 각 영역을 수작업으로 지정한 다음에, 평균 HU 값을 계산하였다. 본 논문에서는 번거로운 수작업을 없애고 자동으로 여섯 개 영역을 분할하여 평균 HU 값을 계산하는 방법을 제안한다. 단면도에서 각 영역은 원 형태이므로, 본 논문에서는 원을 검출하는 방법을 사용하였다. 원을 검출하기 위해 먼저 단면도 영상을 이진화하였다. 이진화를 위한 임계값은 치밀뼈에 해당하는 영역이 객체가 되고 나머지가 배경이 되도록 하는 값으로 사람이 CT 영상을 분석하여 적절한 값을 지정해 주었다. Fig. 3에는 이진화 영상이 나타나 있다. Fig. 3(a)의 영상은 S1 영역 중간에 있는 슬라이스 영상을 이진화한 것이고 Fig. 3(b)의 영상은 S1 영역과 S2 영역이 만나는 경계 부분에 해당하는 슬라이스 영상을 이진화한 것이다. S1 영역과 S2 영역이 만나는 경계 부분에서는 두 영역이 완벽하게 분리되지 않고 그림과 같이 두 영역이 겹쳐서 나타날 수도 있고, 이때에는 원의 형태가 불완전하게 나타날 수가 있다. 이와 같이 두 영역이 겹쳐서 나타나는 부분은 영역의 평균 HU 값을 계산하는 데에서 제외시켜야 한다.

Fig. 3.Binarized image (a) a slice within section 1 (b) a slice between section 1 and section 2.

이진화를 한 다음에는 객체 검출을 용이하게 하기 위해, 객체 내부에 있을 수 있는 빈 공간을 채우는 작업을 수행하였다. 내부 빈 공간을 채운 다음에는 침식연산을 적용하고 그다음에 팽창 연산을 적용하여 잡음을 제거하였다. Fig. 4에는 Fig. 3의 두 슬라이스 영상에 대해 내부를 채운 영상, 침식 연산을 적용한 영상, 그리고 팽창 연산을 적용한 영상이 나타나 있다. Fig. 4(a)의 S1 영역과 S5 영역 모두 내부가 채워졌지만 Fig. 4(b)에서는 S1 영역의 경우 원의 경계선이 증간에 끊겨서 열려져 있으므로 내부가 채워지지 않았다. Fig. 4(a), (b) 영상의 왼쪽 아래 구석 부분에 점들이 나타나 있는데 이는 CT 받침대의 일부 영역이다. 이 점들이 침식 연산과 팽창 연산에 의하여 제거된 것을 Fig. 4(c)-(f)에서 볼 수 있다. 또한 Fig. 4(b)에 나타나 있는 S1 영역 내부의 점들도 침식 연산과 팽창 연산에 의해 제거된 것을 볼 수 있다.

Fig. 4.Filling and noise removing (a), (b) hole filled image (c), (d) eroded image (e), (f) dilated image.

잡음을 제거한 다음에는 객체의 윤곽선을 검출하였다. 윤곽선을 추출하는 방법으로는 여러 가지가 있는데, 본 논문에서는 물체의 가장자리에 있는 한 픽셀로부터 시작하여 가장자리에 있는 픽셀을 계속 따라가다가 시작점에 도착하면 추적을 멈추는 방법을 사용하였다. Fig. 5에는 Fig. 4(e)와 (f)의 영상에 대하여 윤곽선을 검출한 결과가 나타나 있다.

Fig. 5.Extracted contour image.

각 단면도 영상에는 검출하고자 하는 객체로 S1~S6 영역 중 2개의 영역이 포함되어 있다. 이진화 영상에는 이들 2개 영역의 객체 이외에도 다른 객체들이 존재 할 수 있다. 또한 S1~S6 영역이 서로 겹치는 부분에서는 영역이 완벽하게 분리되지 않고 두 영역이 겹쳐서 나타날 수도 있다. 따라서 S1~S6 영역에 대한 온전한 형태의 객체만 검출하기 위해, 객체의 면적과 원형지수를 계산하였다.

윤곽선 내부의 면적은 다각형 면적을 계산하는 공식을 이용하여 구하였다. N-각형의 각 정점의 좌표가 연결되어 있는 순서대로 (x0, y0), (x1, y1), (x2, y2), ⋯, (xN-1, yN-1)이라 할 때에 면적은 식 1와 같이 구할 수 있다. 여기에서 xN = x0, yN = y0이다.

원형 지수는 객체의 모양이 얼마나 원에 가까운지를 나타내는 값으로서 이상적인 원의 경우 원형지수 값이 1이다. 이 값이 작을수록 더 둥글지 않다는 것을 나타낸다. 원형 지수 계산식은 식 2와 같다.

원형지수가 0.7보다 크고 윤곽선 내부면적이 영역 S4 면적의 90% 이상인 객체만 검출하였다. 원형지수의 임계값을 1에 가까운 값으로 설정하지 않은 이유는, 이상적인 원과는 달리 이진화 영상에서 원의 윤곽선이 픽셀로 구성됨으로 인해 매끄럽지 못하고 울퉁불퉁하게 되면서 윤곽선의 길이가 더 길어져서 원형지수 값이 작아지기 때문이다. Fig. 6에는 원형지수와 내부면적 조건을 만족하는 윤곽선을 검출한 결과가 나타나 있다.

Fig. 6.Extracted circles.

발견된 객체가 S1~S6 영역 중 어느 영역에 속하는지는 일단은 객체의 면적에 따라 구분하였다. 그런데 영역 S3과 S5 그리고 S4와 S6이 각각 면적이 같은데, 이때에는 객체의 위치에 따라 구분하였다. 영역 S1, S3, S4는 영역의 중심 위치가 같으므로 영역의 중심 위치를 비교하여 영역 S3과 S5 그리고 S4와 S6를 구분하였다.

모든 슬라이스 영상에 대해 두 개의 영역에 대한 윤곽선을 구한 다음에는, 슬라이스 영상을 처음부터 마지막까지 스캔하면서 중심 좌표가 비슷하고 면적이 비슷한 윤곽선이 나타나기 시작하는 슬라이스와 마지막으로 나타나는 슬라이스를 검출하여 각 영역의 시작 위치와 끝 위치를 검출한다.

각 영역에 대하여 시작 슬라이스와 마지막 슬라이스를 검출한 다음에는 참고문헌 [11]의 방법에 의해 치밀뼈 영역과 해면뼈 영역을 추출한다. 해면뼈 영역을 구하기 위해 먼저, 이진화를 통하여 치밀뼈 내부를 채운 영상을 구하고 치밀뼈 영상을 반전시킨 영상을 구한다. 치밀뼈 내부를 채운 영상과 치밀뼈 영역을 반전한 영상의 교집합을 구하면 해면뼈 영역을 얻게 된다. S1~S6에 대하여 이렇게 구한 해면뼈 영역과 치밀뼈 영역에 대한 영상이 Fig. 7에 나타나 있다.

Fig. 7.3D image of the extracted cortical and trabecular bone region (a) Section 1 (b) Section 2 (c) Section 3 (d) Section 4 (e) Section 5 (f) Section 6.

해면뼈 영역과 치밀뼈 영역을 구한 다음에는 각 영역에 대한 평균 HU 값을 구한다. Table 1에는 S1~S6의 각 부분에 대해 해면뼈 영역과 치밀뼈 영역의 평균 HU 값이 나타나 있다. 표에 나타나 있는 치밀뼈에 대한 HU 값을 살펴보면 치밀뼈의 두께가 큰 경우에 HU 값도 증가하는 것을 볼 수 있다. 해면뼈에 해당하는 영역을 살펴보면 S4, S5, S6 영역에 같은 물질이 채워져 있음에도 불구하고 S4와 S6 영역의 HU 값은 서로 비슷하고 S5 영역의 HU 값이 작게 나왔는데, 이것은 바깥쪽의 치밀뼈의 두께가 다른 것에 기인한 것으로 해석된다.

Table 1.Estimated HU Value for the Six Sections of the Phantom

이렇게 구한 영역별 평균 HU 값은 참고문헌 [11]에서 구한 값과 99.4% 이상 일치하였다. 평균 HU 값에 약간의 차이가 나는 것은 영역의 슬라이스 범위를 수작업으로 지정한 경우와 자동으로 추츨한 경우에 있어서 슬라이스 범위가 약간 차이가 나기 때문이다. Table 1을 살펴보면 S4, S5, S6 영역의 밀도는 동일한데, HU 값이 S4와 S6은 비슷하고 S5는 차이가 나는 것을 볼 수 있다. 영역 S5의 HU 값이 작은 것은 치밀뼈의 두께가 S4와 S6보다 얇기 때문인 것으로 분석된다. 영역 S1, S2, S3의 경우에도 밀도가 증가하는 정도와 HU 값이 증가하는 정도가 일정하지 않음을 알 수 있다. 이와 같이 치밀뼈의 두께에 영향을 받으므로, 치밀뼈 두께가 같은 S1과 S3을 기준으로 HU 값과 밀도 값 사이의 관계식을 구하였다. 관계식은 2차원 평면에서 두 점을 지나는 직선의 방정식을 이용하여 구하였으며, HU 값이 X이고 밀도값이 Y라면 X와 Y는 식 3의 관계를 갖는다. 본 논문의 방법으로 구한 영역 S1과 S3의 평균 HU 값은 참고문헌 [11]과 동일한 값을 가지므로 식 (3)도 참고문헌 [11]과 동일하다.

 

3. 구조적 지표자 추출

해면뼈의 구조적 지표자는 BS(Volume Surface), BV(Bone Volume), TS(Total Surface), TV(Total Volume), BS/BV, BV/TV, BS/TV, Tb.Th(Trabecular Thickness), Tb.Sp(Trabecular Separation), Tb.N(Trabecular Number), Tb.pf(Trabecular Pattern Factor), SMI(Structural Model Index), Ct.Th(Cortical Thickness), Ct.Ar(Cortical Area), D(Displacement)의 15가지 지표자를 추출하였다. 이 지표자들 중 BV/TV, Tb.Th, Tb.Sp, Tb.N은 참고문헌[13]에서 설명한 방법으로 추출하였다. 치밀뼈에 대한 지표자는 Ct라는 이름으로 시작되며, 나머지 지표자는 모두 해면뼈에 대한 지표자이다. BS는 해면뼈의 표면적 넓이이고, BV는 해면뼈가 실제로 차지하는 공간의 부피이다. TS와 TV는 해면뼈 사이의 빈 공간을 포함하여 해면뼈가 존재하는 전체 공간의 표면적과 부피이다. BS, BV, TS, TV 계산을 위하여 해면뼈가 차지하는 전체 영역과 실제로 뼈가 차지하는 영역을 3차원 메쉬로 나타낸 다음에 메쉬가 차지하는 표면적과 부피를 계산하였다. 메쉬가 차지하는 표면적과 부피는 VTK(Visulaization Toolkit)의 라이브러리 함수 중 3차원 메쉬의 표면적과 부피를 계산하는 함수를 사용하여 구했다[14].

지표자 Tb.Th는 해면뼈의 평균 두께를 나타내는데, 해면뼈의 두께는 일정하지 않고 위치에 따라 변화하므로 뼈의 각 위치에서의 두께를 구하여 평균값을 구한다. 해면뼈의 두께를 구하기 위해서는 3차원상에서 해면뼈 영역에 구를 정합하여 두께를 구하는 방법[12]을 이용하였다. 지표자 Tb.Sp는 해면뼈 사이의 평균 거리를 나타내는데, Tb.Sp는 해변뼈가 차지하는 전체 영역에서 빈 공간의 점들에 대해 구 정합을 적용하여 구한다. 지표자 Tb.N은 단위 길이를 통과하는데 얼마나 많은 해면뼈를 거치는지를 나타내는데, 식 4와 같이 간략화된 방법을 이용하여 구한다.

지표자 Tb.pf는 해면뼈의 상대적인 볼록함 또는 오목함을 나타내는 지표자이다. 볼록함의 정도가 크다는 것은 뼈의 연결이 중간에 끊어진 곳이 많아서 뼈가 약하다는 것을 의미하고 오목함의 정도가 크다는 것은 뼈의 연결이 잘되어 있어 뼈가 강하다는 것을 의미한다. 이 지표자는 해면뼈 영역을 팽창시키기 전 후의 표면적과 부피 구하여 식 5와 같이 구하였다[15].

여기에서 S1과 V1은 원래 해면 뼈의 표면적과 부피이고 S2와 V2는 뼈 영역을 팽창시킨 후의 표면적과 부피이다. 지표자 SMI는 막대 형대의 뼈와 평판 형태의 뼈 중에 어느 쪽이 더 우세한 지를 나타내는 지표자이다. 해면뼈는 평판 형태의 뼈가 우세하나 뼈의 소실에 따라 막대 형태의 뼈가 더 우세해지게 된다. 이 지표자도 해면뼈 영역을 팽창시키기 전 후의 표면적과 부피 구하여 식 6과 같이 구하였다[16].

여기에서 S와 V는 원래 해면 뼈의 표면적과 부피이고 S‘는 뼈 영역을 팽창으로 인한 표면적의 변화량이다. 지표자 Ct.Th는 치밀뼈 영역의 두께로서 분할한 치밀뼈 영역에 구정합을 적용하여 구한다. 지표자 Ct.Ar은 치밀뼈의 단면도 상에서의 면적의 합이다. 이 지표자는 치밀뼈 영상의 각 슬라이스 별로 치밀뼈 영역의 면적을 구한 다음 이를 합하여 구한다.

지표자 D는 변위값인데 변위란 뼈를 일정 하중으로 눌렀을 때에, 뼈의 윗면이 얼마만큼 아래 방향으로 눌려서 내려갔는지를 의미한다. 변위값은 영상으로부터 유한 요소 모델링을 이용하여 뼈를 모델링하고 유한 요소 분석을 이용하여 뼈에 압력을 가하고 뼈의 변화를 시뮬레이션하는 방법을 이용하여 구하였다. 다른 지표자들은 해면뼈와 치밀뼈를 분할하여 구하였는데 변위값 계산에서는 해면뼈와 치밀뼈 모두 포함된 뼈에 대해 시뮬레이션을 수행하였다. 본 논문에서는 유한 요소 모델로 Fig. 8(a)와 같은 8개의 노드를 가지는 육면체를 사용하였다. Fig. 8(b)에는 분석 대상 뼈 영상이 나타나 있다. Fig. 8(b)의 분석 대상 뼈 영역을 유한 요소 모델로 변환한 결과가 Fig. 8(c)에 나타나 있다. Fig. 8(c)의 왼쪽에는 전체 유한 요소 모델이 나타나 있고 오른쪽에는 일부를 확대한 모습이 나타나 있다. 그림에서 맨 윗면 사각형의 가장자리 선의 색이 다르게 표시되어 있는 것은, 이 유한 요소 모델의 윗면에 하중 조건을 부여했기 때문이다. 유한 요소 모델의 윗면 노드들에는 1000N의 힘이 주어 졌고 맨 아랫면 노드들은 위치를 고정하여 움직이지 않도록 하였다. 유한 요소 분석을 위해서 유한 요소 모델의 재질을 지정해야 하는데, 본 논문에서는 탄성계수를 10GPa로 설정하였으며 가로-세로 변형 비율인 포아송 비율을 0.3으로 설정하였다. 이러한 탄성계수와 포아송 비율은 뼈에 대한 유한 요소 분석에서 널리 사용하는 값이다[17].

Fig. 8.Finite element modeling (a) hexahedron model (b) bone image (c) generated finite element model.

유한 요소 분석은 ANSYS 프로그램을 이용하여 수행하였다. Fig. 9에는 시뮬레이션 결과를 도식화한 그림이 나타나 있다. 이 그림에서는 노드들의 변위 값에 따라 색을 달리하여 표시하였다. 노드들 마다 변위가 다르므로 지표자 D 값은 윗 면에 있는 노드들의 변위의 평균값으로 설정하였다.

Fig. 9.Finite element analysis result.

 

4. 실험 결과

본 논문에서는 정상군과 골다공증군 각각 18명의 손목뼈 영상을 CT로 촬영하여, 각각의 영상에서 분석영역을 지정한 다음에 해면뼈 영역과 치밀뼈 영역의 골밀도를 측정하였고 15가지의 구조적 지표자를 추출하였다. Table 2에는 골밀도와 지표자들의 평균, 표준편차, 두 군의 값에 대한 T-test의 P 값이 나타나 있다. P 값이 0.05보다 작으면 지표 값이 두 군 사이에 서로 구별될 수 있을 만큼 서로 차이가 난다는 것을 의미한다. 지표자 중 TS와 TV의 경우에만 P 값이 0.5보다 커서 골다공증군과 정상군에서 이 두 지표자 값이 거의 비슷한 것으로 나타났다. TS와 TV는 해면뼈의 빈 공간을 포함하는 전체 영역의 표면적과 부피이므로 정상군과 골다공증군 사이에 비슷한 값을 가지는 것이 당연하다.

Table 2.Extracted BMD and Structural Parameter Values

Tb.BMD는 해면뼈 영역에 대한 골밀도 값이고 Ct.BMD는 치밀뼈 영역에 대한 골밀도 값인데, Tb,BMD 값이 Ct.BMD 값에 비하여 정상군과 골다공증군 사이에 더 큰 차이가 있는 것은 골다공증이 해면뼈에 주로 발생하는 것을 알 수 있다. BS, BV, BV/TV, BS/TV 값들 모두 골다공증군에서 골 소실로 인하여 값이 작아진 것을 볼 수 있다. 또한 골다공증 군에서 골 소실로 인하여 해면뼈의 두께는 얇아지고 해면뼈 사이의 거리는 커지므로 Tb.Th 값은 작아지고 Tb.Sp 값은 커진 것을 알 수 있다. 골다공증 군에서 치밀뼈의 두께와 면적은 감소하기는 했지만 해면뼈에 비하면 그리 크지 않은 것을 알 수 있다.

골밀도와 지표자를 이용하여 골다공증군과 정상군을 판별하는 실험을 수행하였다. 부류를 판별하는 방법으로 여러 가지 방법이 사용되는데, 본 논문에서는 SAD(Sum of Absoulte Difference), SVM(Support Vector Machine), 인공 신경망(ANN: Artificial Neural Network) 세 가지 방법을 사용하였다. 정상군과 골다공증군 각각 18명의 영상 중에 일부는 학습에 사용하고 나머지는 판별에 사용하는 방식으로 판별 실험을 수행하였으며, 학습과 판별에 사용되는 영상을 돌려가며 사용하였다. 예를 들어 17개 영상을 학습에 사용하고 1개 영상을 판별에 사용한다고 할 때, 1번을 제외한 나머지 영상으로 학습을 하고 1번으로 판별하는 경우, 2번을 제외한 영상으로 학습하고 2번으로 판별하는 경우, 마지막으로 18번을 제외한 나머지 영상으로 학습을 하고 1번으로 판별하는 경우 까지 총 18 가지의 경우에 대하여 판별 실험을 수행하였다. 학습에 사용되는 영상의 수도 17개부터 9개까지 변화시켜 가면서 실험을 수행하였다.

SAD를 이용한 방법에서는 정상군과 골다공증군에 대하여 각각 학습에 이용되는 영상에서 추출된 골밀도와 지표자들의 평균값을 구한다. 그리고 판별에 사용되는 영상에서 추출된 골밀도와 지표자과 평균 값 사이의 차이의 제곱 값의 합을 정성군과 골다공증군에 대하여 구한 다음에 이 값이 작은 군에 속하는 것으로 판별하게 된다. SVM과 인공신경망은 패턴 분류에 널리 사용되는 도구인데, 본 논문에서는 SVM 구현 중에서 참고문헌 [18]의 라이브러리를 사용하였고 인공신경망 구현으로는 참고문헌 [19]의 라이브러리를 사용하였다.

골밀도만을 이용하여 판별 실험을 수행하였고 골밀도와 구조적 지표자들을 함께 사용하여 판별실험을 수행하였다. Table 3에는 Tb.BMD 만을 이용하여 분류 실험을 한 결과가 나타나 있다. 학습에 사용된 데이터 수를 정상군과 골다공증군 각각에 대하여 9개부터 17개까지 변화를 주면서 실험을 수행하였다. Tb.BMD 만을 사용한 경우에는 인공신경망, SVM, SAD 순으로 정확도가 향상되었다. 또한, 대체적으로 학습에 사용된 데이터가 많을수록 분류 정확도가 향상되는 것을 알 수 있다.

Table 3.Classification results for the Trabecular BMD

해면뼈의 골밀도 값과 더불어 구조적 지표자를 사용하여 분류 실험을 수행하였다. Table 2에서 해면뼈의 골밀도 값을 제외한 나머지 지표자 값을 P 값에 따라 정렬하여, P 값이 작은 것부터 함께 사용하였으며, 추가 지표자의 수를 변화시켜 가면서 분류 실험을 하였다. P 값이 작은 지표자로부터 나열하면, Tb.pf, Tb.Sp, BV/TV, D, Tb.N, SMI, BV, Ct.BMD, BS/TV, BS, Tb.Th, BS/BV, Ct.Th, Ct.Ar, TV, TS 순이었다. Table 4에는 Tb.BMD 만 사용한 경우와 더불어, Tb.BMD와 한 개의 지표자를 함께 사용한 경우부터 열다섯 개의 지표자를 모두 함께 사용한 경우까지의 분류 실험 결과가 나타나 있다. 여기에 나타나 있는 분류 정확도 비율은 학습 데이터의 수를 9개부터 17개까지 변화시켜가면서 분류 실험한 값을 평균한 값이다. 지표자를 여러개 사용할 때, 지표자들 사이의 값의 범위가 서로 다르면 문제가 발생하게 된다. Table 2의 지표자들의 평균값을 살펴보면 1보다 작은 값이 있는 반면에 수천단위의 값도 있다. 이들 값을 그대로 사용하면 큰 값을 가지는 지표자가 판별에 더 큰 비중을 차지하게되는 문제가 발생한다. 따라서 본 논문에서는 골밀도와 지표자 값들을 -1에서 1 사이의 값을 가지도록 정규화한 다음에 SAD, SVM, ANN 방법을 적용하였다. Tb.BMD 값만을 사용한 경우보다 지표자들을 함께 사용함으로써 분류의 정확도가 대체적으로 향상되는 것을 알 수 있다. SAD는 9개, SVM은 10개, ANN은 12개의 지표자를 함께 사용하는 경우에 최대의 정확도를 보였다. 함께 사용한 지표자의 수가 증가함에 따라 정확도가 향상되다가 나중에는 오히려 정확도가 떨어지는 것으로 나타났는데, 이는 나중에 추가된 지표자일수로 P 값이 큰 것으로서 골다공증군과 정상군 사이의 변별력이 떨어지는 지표자이기 때문인 것으로 분석된다. 최대 정확도는 SVM을 사용한 경우에 91.0%인 것으로 나타났다. 이와 같은 실험을 통하여 골밀도 하나의 값을 사용하는 경우보다 다른 지표자들을 함께 사용함으로써 분류의 정확도가 향상됨을 알 수 있었다.

Table 4.Classification results for the additional parameters

위와 같이 골밀도와 지표자 값을 직접 사용하는 방법에 더불어, 주성분 분석을 수행하여 값을 변환하여 분류 실험을 수행하였다. 주성분 분석은 데이터로부터 중요한 성분 값을 추출하는 방법으로서 패턴 분류에서 입력 데이터의 차원을 축소하는데 널리 사용되는 방법이다. 본 논문에서 사용한 데이터의 값은 총 16개이기 때문에 차원을 축소할 필요는 없지만, 주성분 분석을 통하여 추출한 값을 사용함으로써 분류의 정확도에 영향을 받을 수 있기 때문에 주성분 분석을 활용해 보았다. Table 5에는 주성분 분석을 적용했을 경우의 분류 실험 결과가 나타나 있다. 골밀도와 지표자 값을 포함해서 총 16개의 데이터 값에 대하여 주성분 분석을 수행하였다. 그리고 성분값을 중요도에 따라 정렬하고 제일 중요한 성분값부터 시작하여, 사용되는 성분값 개수를 증가시켜 가면서 SAD, SVM, ANN을 적용하였다. 주성분값의 개수가 적을 때에 정확도가 높은 것으로 나타났으며, 주성분값의 개수에 따라 정확도가 약간 감소하기는 했지만 Table 4의 결과에 비하면 편차가 적은 편이었다. Table 4와 비교하여 보면 평균적으로는 주성분 분석을 활용함으로써 분류의 정확도가 향상되는 것을 알 수 있다.

Table 5.Classification results with principal component analysis

Table 4와 5를 살펴보면 참고문헌 [11]의 방법은 SAD, SVM, ANN 각각 88.0%, 86.5%, 79.9%의 인식률을 보였고, 본 논문에서는 실험 결과 지표자를 16개로 확장함으로써 평균적으로 SAD, SVM, ANN 각각 89.2%, 87.4%, 84.7%의 인식률을 얻어 평균적으로 2.8%의 인식률 향상을 수 있었다. 또한 지표자 확장과 주성분 분석을 함께 이용함으로써 SAD, SVM, ANN 각각 90.9%, 89.5%, 86.0%의 인식률을 얻어 참고문헌 [11]에 비해 평균적으로 4.8%의 인식률 향상을 얻을 수 있었다.

 

5. 결 론

본 논문에서는 구조적 지표자를 확장하고 주성분 분석기법을 이용하여 기존의 골다공증 판별의 성능을 향상시키는 방법을 제안하였다. 기존의 방법에서 사용한 6개의 지표자에 추가로 10개의 지표자를 CT 영상으로부터 추출하였다. 특히, 유한 요소 분석을 이용하여 뼈에 힘이 가해졌을 때에 발생할 수 있는 변위를 측정하여 뼈의 강도를 나타내는 지표도 추가하였다. 실험 결과 본 논문에서 제안된 방법은 기존의 방법에 비하여 골다공증 판별 정확도를 지표자 확장을 이용해서 2.8% 향상시켰고 지표자 확장과 주성분분석을 함께 사용해서는 4.8% 향상시켰다. 또한 본 논문은 CT 값을 골밀도 값으로 변환하기 위해 필요한 팬텀 분석을 자동을 하는 방법을 제안하였다. 기존의 방법에서는 팬텀 내의 뼈 영역을 수작업으로 지정해주어야 했는데, 본 논문의 방법에서는 CT 영상에서 원을 자동으로 검출함으로써 뼈 영역을 검출한다. 제안된 방법은 기존 방법의 수작업 과정을 자동화하면서도 같은 분석 결과를 생성하였다.