CT HU 값, 직접 재 보기: pydicom으로 공개 CT를 3D HU 볼륨으로 만드는 실습

의료 AI 개발자가 처음 알아야 할 것: DICOM, HU, CT Normalization에서는 DICOM, HU, windowing, normalization이 어떻게 이어지는지 개념 중심으로 정리했다. 그 글을 찾아 들어오는 검색어를 보면 질문은 더 구체적이다. "CT HU 값이 얼마인가", "HU 뜻", "Hounsfield unit 정의".

이번 글은 그 질문에 손으로 답해 보는 실습편이다. 누구나 받을 수 있는 공개 흉부 CT 한 건을 내려받아, DICOM 폴더를 3D HU 볼륨으로 만들고, 공기·폐·지방·근육·조영된 혈관의 HU를 직접 잰다. 그 과정에서 실제 데이터에서 나온 함정 세 가지도 숫자로 확인한다. 이 글의 코드는 모두 아래 환경에서 실제로 실행했고, 본문의 숫자는 그 실행 결과다.

• Python 3.12, pydicom 3.0.2, NumPy 2.5.4, SciPy 1.18.1, matplotlib
• 데이터: LIDC-IDRI-0001 흉부 CT 1개 series (TCIA 공개, CC BY 3.0)

HU 한 줄 정의

HU(Hounsfield Unit)는 각 voxel의 X선 감쇠계수 μ를 물과 공기에 맞춰 다시 눈금을 매긴 값이다.

HU = 1000 × (μ − μ_water) / (μ_water − μ_air)

정의에 따라 물은 0 HU, 공기는 −1000 HU다. 나머지 조직의 값은 정의가 아니라 측정값이어서, 문헌마다 범위가 조금씩 다르고 촬영 조건과 조영제에 따라 달라진다. 그래서 표를 외우기보다 직접 재 보는 편이 이해가 빠르다.

1. 공개 CT 받기

미국 국립암연구소(NCI)가 지원하는 공개 의료영상 저장소 The Cancer Imaging Archive(TCIA)의 LIDC-IDRI 컬렉션에서 첫 번째 환자의 CT series를 쓴다. 공개 API라 로그인 없이 받을 수 있다.

# 환자 LIDC-IDRI-0001의 series 목록 (JSON)
curl -s "https://services.cancerimagingarchive.net/nbia-api/services/v1/getSeries?Collection=LIDC-IDRI&PatientID=LIDC-IDRI-0001&format=json"

# 그중 CT series (133장) 내려받기: zip 약 35 MB, 풀면 67 MB
curl -s -o series.zip "https://services.cancerimagingarchive.net/nbia-api/services/v1/getImage?SeriesInstanceUID=1.3.6.1.4.1.14519.5.2.1.6279.6001.179049373636438705059720603192"
mkdir -p lidc && unzip -q series.zip -d lidc

LIDC-IDRI는 CC BY 3.0 라이선스이며 사용 시 데이터셋과 논문을 인용해야 한다(글 끝 참고문헌). 또 TCIA 이용 조건상 데이터로 참여자를 식별하려 해서는 안 된다.

2. 픽셀보다 헤더를 먼저 본다

첫 장의 헤더에서 HU 계산과 공간 정보에 필요한 값만 뽑으면 다음과 같다.

3. 함정 1: 파일 이름 순서는 슬라이스 순서가 아니다

파일 이름으로 정렬한 처음 다섯 장의 InstanceNumber와 z 위치는 이렇다.

이름순으로 쌓으면 장기가 뒤죽박죽 섞인 볼륨이 나오고, 에러는 나지 않는다. 그래서 ImageOrientationPatient로 슬라이스 면의 법선 n을 구하고, 각 슬라이스의 ImagePositionPatient를 n에 투영한 값으로 정렬한다. 정렬 후 133장의 위치는 −340 mm부터 −10 mm까지였고, 132개 간격이 모두 정확히 2.5 mm였다. 이 series에서는 InstanceNumber가 z와 정확히 반대 방향으로 붙어 있었다(가장 아래 슬라이스가 133번). 그래서 InstanceNumber로 정렬해도 "순서"는 맞지만 위아래가 뒤집힌다.

ct_volume.py:

import sys
from collections import defaultdict
from pathlib import Path

import numpy as np
import pydicom


def load_series(folder):
    """폴더 안의 DICOM을 SeriesInstanceUID별로 묶는다 (픽셀은 아직 읽지 않음)."""
    series = defaultdict(list)
    for path in Path(folder).rglob("*.dcm"):
        ds = pydicom.dcmread(path, stop_before_pixels=True)
        if ds.get("Modality") == "CT" and "LOCALIZER" not in ds.get("ImageType", []):
            series[ds.SeriesInstanceUID].append((path, ds))
    return series


def sort_slices(items):
    """파일명·InstanceNumber가 아니라 슬라이스 법선 위의 위치로 정렬한다."""
    iop = np.array(items[0][1].ImageOrientationPatient, dtype=float)
    normal = np.cross(iop[:3], iop[3:])
    pos = [float(np.dot(normal, np.array(ds.ImagePositionPatient, dtype=float))) for _, ds in items]
    order = np.argsort(pos)
    return [items[i] for i in order], np.array(pos)[order]


def to_hu(ds):
    """Stored value -> HU. 슬라이스마다 slope/intercept를 따로 읽고, 반드시 float로 계산한다."""
    slope = float(ds.get("RescaleSlope", 1))
    intercept = float(ds.get("RescaleIntercept", 0))
    return ds.pixel_array.astype(np.float32) * slope + intercept


def load_volume(items):
    items, pos = sort_slices(items)
    gaps = np.diff(pos)
    if not np.allclose(gaps, gaps[0], atol=0.01):
        raise ValueError(f"슬라이스 간격이 일정하지 않음: {np.unique(gaps.round(3))}")
    vol = np.stack([to_hu(pydicom.dcmread(p)) for p, _ in items])  # (z, y, x)
    row_mm, col_mm = (float(v) for v in items[0][1].PixelSpacing)
    return vol, (abs(float(gaps[0])), row_mm, col_mm), items


if __name__ == "__main__":
    series = load_series(sys.argv[1])
    uid, items = max(series.items(), key=lambda kv: len(kv[1]))
    vol, spacing, items = load_volume(items)
    print("shape (z,y,x):", vol.shape, "spacing mm (z,y,x):", spacing)
    print("HU min/max:", vol.min(), vol.max())

$ python ct_volume.py lidc
shape (z,y,x): (133, 512, 512) spacing mm (z,y,x): (2.5, 0.703125, 0.703125)
HU min/max: -2048.0 3071.0

정렬된 133장의 픽셀을 읽어 HU 볼륨으로 쌓는 데 노트북(Apple M4)에서 약 0.1초 걸렸다. 간격 검사는 일부러 엄격하게 두었다. 슬라이스가 빠졌거나 중복된 series를 조용히 통과시키는 것보다 멈추는 편이 낫다.

4. 함정 2: 정수 그대로 계산하면 값이 감긴다

이 series의 저장값은 부호 있는 int16이라 문제가 없었다. 하지만 PixelRepresentation이 0인 데이터는 pydicom이 uint16 배열을 준다. 이때 intercept를 정수 연산으로 더하면 NumPy는 경고 없이 값을 감는다.

import numpy as np

stored = np.array([24, 3000], dtype=np.uint16)
print(stored - 1024)            # [64536  1976]  <- 24 - 1024 가 64536 이 된다
print(stored.astype(np.float32) * 1.0 - 1024)   # [-1000.  1976.]

위에서 −1000 HU(공기)여야 할 값이 64536이 됐다. to_hu()가 먼저 float32로 바꾸는 이유다. pydicom의 pydicom.pixels.apply_rescale()을 써도 된다. 이 데이터에서 직접 계산한 값과 apply_rescale() 결과는 완전히 같았다(최대 차이 0, 결과 dtype은 float64).

5. 함정 3: 볼륨의 21%는 몸도 공기도 아니다

위 출력의 최솟값 −2048 HU는 공기(−1000)보다도 낮다. 볼륨 전체에서 가장 많은 값을 세어 보면 이렇다.

CT는 원형 영역만 재구성하는데, 512 × 512 정사각형 배열의 네 귀퉁이는 어딘가에 값을 넣어야 한다. 이 장비는 거기에 −2048 HU를 넣었다. 그런데 헤더의 PixelPaddingValue는 63536이다. 16비트 부호 있는 정수로 읽으면 −2000이고, 실제 볼륨에 −2000은 하나도 없었다. 헤더의 채움값을 믿지 말고 히스토그램을 확인해야 한다.

이 값을 그대로 두면 생기는 일:

• 통계가 틀어진다. 슬라이스 80의 평균 HU는 전체로 계산하면 −861.1, FOV 안만 계산하면 −540.3이다.
• min-max 정규화가 망가진다. 볼륨 최솟값이 −1000이 아니라 −2048이 되어, 실제 조직이 쓰는 범위가 좁아진다. 정규화 전에 고정 HU 범위로 clipping해야 하는 이유 중 하나다.
• 리샘플링하면 경계로 번진다. 선형 보간은 FOV 경계에서 −2048과 몸 밖 공기를 섞는다. 아래 리샘플링 코드는 그래서 먼저 FOV 밖을 −1024로 바꾼다.

6. 조직별 HU를 직접 재 보기

정렬된 볼륨에서 두 슬라이스를 골라 원형 ROI(관심 영역)를 놓고 평균과 표준편차를 쟀다. 슬라이스 80은 대동맥궁 아래 흉부, 슬라이스 40은 간과 비장이 보이는 상복부다. 이 CT는 조영제를 쓴 검사다(대동맥이 하얗게 보인다).

import numpy as np

def roi_stats(vol, k, x, y, r):
    """슬라이스 k에서 (x, y) 중심, 반지름 r 픽셀 원 안의 HU 평균·표준편차."""
    yy, xx = np.mgrid[: vol.shape[1], : vol.shape[2]]
    v = vol[k][(xx - x) ** 2 + (yy - y) ** 2 <= r * r]
    return round(float(v.mean()), 1), round(float(v.std()), 1), int(v.size)

print(roi_stats(vol, 80, 440, 330, 4))   # 피하지방 -> (-110.3, 6.5, 49)

참고 범위는 Wikipedia "Hounsfield scale" 표의 값이고, 그 표가 인용한 원 문헌은 각 항목마다 다르다(참고문헌 4). 읽는 법:

• 정의와 맞는 값: 공기 −1010, 지방 −110, 근육 45 HU는 기대 범위 안이다. 공기가 −1000보다 10 낮은 것도 이상하지 않다. 측정에는 잡음과 장비 보정 오차가 있다.
• 조영제가 바꾼 값: 대동맥 436 HU는 혈액이 아니라 요오드 조영제의 값이다. 간 72 HU도 표의 60 ± 6보다 높고, 비장 173 HU는 표의 "조영 CT의 연부조직" 범위에 들어간다. 같은 장기라도 조영 여부와 촬영 시점(phase)에 따라 HU가 크게 달라지므로, 학습 데이터의 phase가 섞여 있는지 먼저 확인해야 한다.
• 참고 범위와 다른 값: 폐는 −877~−904 HU로 표의 −700~−600보다 낮았다. 폐 HU는 공기를 얼마나 들이마셨는지에 따라 크게 달라진다. 척추체 116 HU도 표의 해면골 범위보다 낮다. 해면골 HU는 나이와 골밀도, 측정 위치에 따라 달라지므로, 이 숫자 하나로 골밀도를 판단해서는 안 된다.
• ROI 위치가 숫자를 바꾼다: 근육 ROI를 처음에는 척추 반대편에 놓았는데 8 HU(SD 23)가 나왔다. 같은 슬라이스의 등 근육 위에서 반지름 4픽셀 ROI를 여섯 곳에 놓아 보니 평균이 −16 HU부터 45 HU까지 흩어졌다. 근육 사이 지방이 섞이면(partial volume) 평균이 내려간다. "근육은 40 HU"처럼 표의 숫자를 그대로 threshold로 쓰면 이런 차이를 놓친다.

7. 같은 HU, 다른 window

앞 글에서 설명한 windowing을 같은 HU 배열에 적용하면 다음과 같다. 이 series 헤더에 저장된 window는 폐 window(−600 / 1600)였다.

def window(img, center, width):
    lo, hi = center - width / 2, center + width / 2
    return ((np.clip(img, lo, hi) - lo) / (hi - lo) * 255).astype(np.uint8)

window는 표시와 모델 입력 범위를 고르는 일일 뿐, HU 자체를 바꾸지 않는다. 그래서 측정과 통계는 항상 window 전의 HU 배열에서 한다.

8. 1 mm 등방성 볼륨으로 리샘플링

이 볼륨의 voxel은 0.70 × 0.70 × 2.5 mm로 z 방향만 3.6배 길다. 3D 모델에 넣거나 길이를 재려면 보통 등방성으로 맞춘다.

from scipy import ndimage

outside = vol <= -2000                                  # FOV 밖 채움값
vol_clean = np.where(outside, -1024, vol)              # 공기 값으로 바꾼 뒤 보간
iso = ndimage.zoom(vol_clean, zoom=spacing, order=1)   # spacing = (2.5, 0.703125, 0.703125)
print(iso.shape)                                       # (332, 360, 360)

zoom 배율은 "현재 간격 ÷ 목표 간격"이므로 목표가 1 mm면 현재 간격 그대로가 배율이 된다. 결과는 (133, 512, 512)에서 (332, 360, 360)이 됐고, 선형 보간(order=1)으로 0.86초, float32 기준 172 MB였다. 같은 배열 크기라도 물리적 범위가 달라질 수 있다는 앞 글의 설명이 여기서 숫자로 보인다. 원본 512 픽셀은 360 mm를, z 133장은 330 mm를 덮는다. 라벨 마스크를 같이 리샘플링할 때는 order=0(최근접)을 써야 클래스 값이 섞이지 않는다.

체크리스트

• SeriesInstanceUID로 묶고, localizer는 뺀다.
• 파일명이나 InstanceNumber가 아니라 ImagePositionPatient · 법선으로 정렬하고, 간격이 일정한지 검사한다.
• slope/intercept는 슬라이스마다 읽고, float로 바꾼 뒤 계산한다.
• 히스토그램으로 FOV 밖 채움값을 찾는다. 헤더의 PixelPaddingValue만 믿지 않는다.
• 조영 여부와 phase를 확인한 뒤 HU를 해석한다.
• 측정·통계는 window 전 HU로 하고, 정규화 전에 고정 범위로 clipping한다.
• 리샘플링 전에 채움값을 정리하고, 마스크는 최근접 보간을 쓴다.

위 코드 블록만으로 본문의 표와 숫자를 재현할 수 있다. ROI 좌표는 표의 번호 순서대로 (슬라이스, x, y, 반지름) = (80, 256, 40, 10), (80, 150, 300, 10), (80, 360, 300, 10), (80, 440, 330, 4), (80, 280, 405, 4), (40, 150, 220, 12), (40, 390, 300, 10), (80, 268, 328, 8), (80, 238, 356, 6)이다. 슬라이스 번호는 z 오름차순으로 정렬한 볼륨 기준이다.

참고문헌

1. Armato SG III, et al. The Lung Image Database Consortium (LIDC) and Image Database Resource Initiative (IDRI): A completed reference database of lung nodules on CT scans. Medical Physics 38(2):915–931, 2011. doi:10.1118/1.3528204
2. Armato SG III, et al. Data From LIDC-IDRI [Data set]. The Cancer Imaging Archive, 2015. doi:10.7937/K9/TCIA.2015.LO9QL9SX
3. Clark K, et al. The Cancer Imaging Archive (TCIA): Maintaining and Operating a Public Information Repository. Journal of Digital Imaging 26(6):1045–1057, 2013. doi:10.1007/s10278-013-9622-7
4. Wikipedia, Hounsfield scale (HU 정의식과 조직별 참고 범위 표, 각 항목의 원 문헌은 해당 페이지 참고).
5. DICOM PS3.3, C.8.2.1 CT Image Module (Rescale Slope/Intercept), C.7.6.2 Image Plane Module (ImagePositionPatient, ImageOrientationPatient, PixelSpacing).
6. pydicom 문서, pydicom.pixels.apply_rescale.