CT 폐 분할, HU 임계값으로 충분할까? 공개 CT 6건으로 lungmask U-net과 직접 비교

지난 실습에서는 공개 흉부 CT를 HU 볼륨으로 만들고 조직별 HU를 직접 쟀다. 폐는 대략 −850 HU 근처였다. 그러면 바로 다음 질문이 나온다. "HU가 낮은 곳만 고르면 폐 마스크가 되지 않을까?" 인터넷의 폐 분할 튜토리얼 대부분이 이렇게 한다. 그런데 기준값은 −150, −320, −400, −500 HU로 글마다 다르고, 딥러닝 모델과 비교해 얼마나 맞는지 잰 글은 찾기 어렵다.

그래서 직접 비교했다. 공개 CT 6건에 같은 임계값 방법을 −300부터 −800 HU까지 돌리고, 공개 딥러닝 폐 분할 모델인 lungmask(U-net, R231)의 결과와 Dice, 폐 부피, 처리 시간을 비교했다. 이 글의 숫자는 모두 아래 코드를 직접 돌려서 얻은 값이다. 환경: Apple M4, 16 GB RAM, macOS 26.6.2, Python 3.12, CPU만 사용.

요약하면 이렇다.

• −400 HU가 가장 잘 맞았다. Dice 중앙값 0.984(0.974–0.988), 폐 부피 차이 −14~+64 mL. −800 HU는 한 건에서 Dice 0.729, 부피 −1.42 L로 무너졌다.
• 불일치는 대부분 정해진 곳에서 나온다. 기관·주기관지(건당 35–70 mL), 폐 가장자리 2픽셀 안쪽(31–62 mL), 그리고 슬라이스별 구멍 채우기가 폐로 잘못 채운 연부조직(3건에서 5–16 mL)이다.
• 속도는 약 60배 차이 났다. 임계값 0.7–1.4초, U-net(CPU) 44–90초.

1. 데이터: LIDC-IDRI 공개 CT 6건

지난 글과 같이 TCIA의 LIDC-IDRI 컬렉션을 쓴다. 이번에는 환자 0001–0006의 CT series 6개를 받았다. 슬라이스 간격은 1.25 mm 2건, 2.5 mm 4건이고 픽셀 간격은 0.625–0.82 mm다. 로그인 없이 받을 수 있다.

mkdir -p data && cd data
while read p uid; do
  mkdir -p $p && curl -s -o $p.zip \
    "https://services.cancerimagingarchive.net/nbia-api/services/v1/getImage?SeriesInstanceUID=$uid" \
    && unzip -q $p.zip -d $p && rm $p.zip
done <<'EOF'
0001 1.3.6.1.4.1.14519.5.2.1.6279.6001.179049373636438705059720603192
0002 1.3.6.1.4.1.14519.5.2.1.6279.6001.619372068417051974713149104919
0003 1.3.6.1.4.1.14519.5.2.1.6279.6001.170706757615202213033480003264
0004 1.3.6.1.4.1.14519.5.2.1.6279.6001.323541312620128092852212458228
0005 1.3.6.1.4.1.14519.5.2.1.6279.6001.129007566048223160327836686225
0006 1.3.6.1.4.1.14519.5.2.1.6279.6001.132817748896065918417924920957
EOF

6건을 풀면 525 MB다. LIDC-IDRI는 CC BY 3.0 라이선스이므로 쓸 때 데이터셋과 논문을 인용해야 한다(글 끝 참고문헌). 또 TCIA 이용 조건에 따라 데이터로 참여자를 식별하려 해서는 안 된다.

pip install lungmask SimpleITK pydicom scipy matplotlib   # lungmask가 PyTorch를 함께 설치한다

이번에는 pydicom 대신 SimpleITK로 읽는다. ImageSeriesReader가 슬라이스를 위치 순서로 정렬하고 RescaleSlope/Intercept도 적용해 주기 때문에, 지난 글에서 손으로 처리한 함정 두 개(파일 순서, 정수 연산)가 한 번에 해결된다. lungmask도 SimpleITK 이미지를 입력으로 받는다.

2. 임계값 방법: 세 단계

lung_threshold.py:

import numpy as np
import SimpleITK as sitk
from scipy import ndimage


def read_ct(folder):
    """DICOM series -> (SimpleITK 이미지, HU 배열 (z, y, x), 복셀 크기 mm (z, y, x))."""
    reader = sitk.ImageSeriesReader()
    reader.SetFileNames(reader.GetGDCMSeriesFileNames(str(folder)))  # 위치 기준 정렬
    img = reader.Execute()  # RescaleSlope/Intercept가 적용된 HU
    hu = sitk.GetArrayFromImage(img).astype(np.int16)
    return img, hu, img.GetSpacing()[::-1]


def fill_holes_2d(mask, spacing=None, max_hole_cm2=None):
    """슬라이스마다 구멍 채우기. max_hole_cm2를 주면 그보다 큰 구멍(종격동 등)은 남긴다."""
    out = mask.copy()
    for z in range(mask.shape[0]):
        holes = ndimage.binary_fill_holes(mask[z]) & ~mask[z]
        if max_hole_cm2 is not None:
            lab, n = ndimage.label(holes)
            area = np.bincount(lab.ravel(), minlength=n + 1) * spacing[1] * spacing[2] / 100
            area[0] = np.inf
            holes = np.isin(lab, np.flatnonzero(area <= max_hole_cm2))
        out[z] |= holes
    return out


def lung_mask(hu, threshold=-400, spacing=None, max_hole_cm2=None, min_ratio=0.1):
    """threshold HU 미만의 '공기'에서 몸 밖 공기를 지우고 큰 덩어리만 남긴다."""
    air = hu < threshold
    labels, n = ndimage.label(air)  # 3D 6-연결
    # 1) 어느 슬라이스에서든 영상 가장자리에 닿는 성분 = 몸 밖 공기 (FOV 밖 -2048 포함)
    edge = np.unique(np.concatenate([labels[:, 0, :].ravel(), labels[:, -1, :].ravel(),
                                     labels[:, :, 0].ravel(), labels[:, :, -1].ravel()]))
    sizes = np.bincount(labels.ravel(), minlength=n + 1)
    sizes[0] = 0
    sizes[edge] = 0
    # 2) 가장 큰 성분 + 그 10% 이상인 성분 (좌우 폐가 따로 잡히는 경우)
    keep = np.flatnonzero(sizes >= min_ratio * sizes.max())
    mask = np.isin(labels, keep)
    # 3) 폐 안의 혈관·결절은 HU가 높아 빠져 있으니 슬라이스마다 구멍을 채운다
    return fill_holes_2d(mask, spacing, max_hole_cm2)

1. 몸 밖 공기 지우기. 임계값보다 낮은 복셀을 3D로 연결해 덩어리로 나누고, 영상 가장자리에 닿는 덩어리는 버린다. 몸 밖 공기와, 지난 글에서 본 촬영 범위(FOV) 밖 −2048 HU 영역이 여기서 함께 빠진다.
2. 큰 덩어리만 남기기. 장 속 가스처럼 작은 공기 덩어리를 버린다. 좌우 폐는 기관으로 이어져 한 덩어리가 되는 경우가 많지만, 따로 잡혀도 둘 다 남도록 "가장 큰 것의 10% 이상"을 기준으로 했다.
3. 구멍 채우기. 폐 안의 혈관은 HU가 높아서 마스크에 구멍으로 남는다. 슬라이스마다 binary_fill_holes로 채운다. 아래에서 보겠지만 문제가 생기는 곳도 여기다.

3. 비교 기준: lungmask U-net

lungmask는 다양한 질환이 있는 흉부 CT로 학습한 공개 U-net이다(Hofmanninger 외, 2020). 사용법은 세 줄이다. 처음 실행할 때 가중치 119 MB를 내려받는다.

from lungmask import LMInferer

inferer = LMInferer(modelname="R231", force_cpu=True)
seg = inferer.apply(img)       # (z, y, x), 0 = 배경, 1 = 오른쪽 폐, 2 = 왼쪽 폐
ref = seg > 0

주의할 점: LIDC-IDRI에는 사람이 그린 폐 마스크가 없다. 그래서 이 글의 Dice는 "정확도"가 아니라 U-net과 얼마나 같은가를 나타낸다. U-net도 틀릴 수 있다. 그래서 숫자만 보지 않고, 어디서 서로 다른지를 그림으로 직접 확인했다.

비교 스크립트 compare.py는 6건 각각에 대해 U-net을 한 번 돌리고, 임계값 6개(−300 ~ −800 HU)로 마스크를 만들어 Dice, 폐 부피, "임계값에만 있는 부피"와 "U-net에만 있는 부피"(mL)를 기록한다. 임계값 방법의 시간은 3번 잰 값의 중앙값이다. 전체 코드는 글 끝에 있다.

4. 결과 1: 임계값은 몇 HU가 좋을까

U-net 기준 폐 부피는 3.44–5.87 L였다. 그러니 −400 HU에서 부피 차이는 2% 이내다.

−800 HU가 무너진 이유도 숫자로 확인했다. 폐 HU가 −850 근처라고 해서 폐의 모든 복셀이 −800보다 낮은 것은 아니다. U-net 마스크 안에서 −800 HU보다 높은 복셀의 비율은 16.5–25.7%였고, 0004번은 48.8%였다(폐 HU 중앙값 −802, 다른 5건은 −844 ~ −886). 같은 환자라도 숨을 얼마나 들이쉬었는지, 어느 부위가 아래쪽에 깔렸는지에 따라 폐 밀도가 달라진다. 그래서 폐 "평균"에 가까운 임계값은 위험하다. −400 HU 위에 있는 폐 복셀은 2.2–5.2%뿐이었다. 이 여유가 −400을 안정적으로 만든다.

5. 결과 2: 불일치는 어디서 나오나

Dice 0.98은 높아 보이지만, 폐처럼 큰 장기에서는 Dice가 너그럽다. −400 HU에서도 건당 70–146 mL가 "임계값에만" 있었다. 그래서 그 부피를 3D 덩어리 단위로 쪼개 보았다.

• 기관·주기관지: 34.6–70.3 mL. 6건 모두에서 "임계값에만 있는" 덩어리 중 가장 큰 것은 영상 가운데를 위아래로 지나는 관이었다(평균 −768 ~ −883 HU). 기관 안은 공기라서 임계값으로는 폐와 구별할 수 없다. U-net은 기관을 폐에서 뺀다. 폐 부피나 폐 안의 병변 비율을 잴 때는 이 차이가 그대로 오차가 된다.
• 폐 가장자리: 31.4–62.1 mL. 나머지 대부분은 U-net 마스크 경계에서 2픽셀 안쪽에 있었다. 흉막 경계의 부분용적(partial volume) 때문에 생기는, 어느 쪽이 맞다고 하기 어려운 차이다. 반대로 "U-net에만 있는" 복셀은 평균 −207 ~ −274 HU로, 경계와 폐 안쪽의 밀도 높은 부분이었다.
• 구멍 채우기가 만든 오류: 4.8–16.2 mL (3건). 0001, 0003, 0004번에서는 평균 −15 ~ +20 HU, 즉 연부조직 덩어리가 폐로 채워졌다. 위 그림 오른쪽(0003번 폐 바닥)을 보면 횡격막 아래의 위(胃) 윗부분이 이 슬라이스에서 폐에 사방으로 둘러싸여 있다. 그래서 2D 구멍 채우기가 이것을 "혈관"과 똑같이 채워 버렸다. 바로 아래 슬라이스들을 보면 조영제가 찬 위라는 것을 확인할 수 있다.

작은 수정: 큰 구멍은 채우지 않기

혈관 단면은 작고 위나 심장 단면은 크다. 그래서 max_hole_cm2=10으로 10 cm²보다 큰 구멍은 채우지 않게 했다. 결과는 이렇다.

Dice는 거의 그대로지만, 엉뚱한 장기 16 mL가 폐로 들어가는 일은 막았다. 다만 10 cm²는 이 6건에서 정한 값이다. 큰 결절이나 종괴가 폐 안에 있으면 그것도 구멍으로 남는다. 결절 분석이 목적이라면 이 값을 그대로 쓰지 말고 자기 데이터에서 다시 확인해야 한다.

기관을 빼는 것은 더 어렵다. 우리는 맨 위 슬라이스의 기관에서 시작해 −950 HU 미만으로 영역을 키우는 간단한 방법을 시도했다. 하지만 6건 중 의미 있게 기관을 잡은 것은 1건뿐이었다(0005, "임계값에만 있는 부피" 74.6 → 43.0 mL). 기관 제거까지 필요하다면, 이 부분은 학습된 모델이 확실히 낫다.

6. 결과 3: 시간

Apple M4 CPU 기준으로 57–64배 차이였고, 둘 다 슬라이스 수에 거의 비례했다. lungmask를 CPU에서 돌리면 512×512 슬라이스 하나에 약 0.33초가 걸린다. GPU를 쓰면 훨씬 빨라지겠지만, 이번에는 GPU에서 재지 않았다.

7. 언제 무엇을 쓸까

• 대량 CT의 빠른 전처리(폐 영역 crop, 학습 데이터 bounding box, 품질 점검): −400 HU 임계값으로 충분하다. U-net과의 차이는 대부분 기관과 가장자리였다.
• 폐 부피·폐 안 병변 비율처럼 숫자를 보고하는 경우: 기관이 35–70 mL씩 섞여 들어가므로 U-net 같은 학습 모델을 쓰거나, 기관 제거 단계를 따로 검증해야 한다.
• 임계값을 낮추고 싶을 때: −600 HU 아래로 내리면 폐 밀도가 높은 스캔에서 폐가 크게 빠진다(−800 HU에서 한 건은 1.4 L). 폐 "평균"이 아니라 폐 HU 분포의 위쪽 꼬리를 보고 정해야 한다.

체크리스트

• 슬라이스 정렬과 slope/intercept 적용을 확인한다(SimpleITK ImageSeriesReader가 둘 다 해 준다).
• 몸 밖 공기는 "가장자리에 닿는 3D 성분"으로 지운다. FOV 밖 −2048 영역도 함께 빠진다.
• 임계값은 −400 HU 근처에서 시작하고, 자기 데이터의 폐 HU 분포로 검증한다.
• 2D 구멍 채우기는 크기 제한을 두고, 결과를 그림으로 확인한다.
• 기관·주기관지가 포함되는지 확인한다. 부피를 보고한다면 반드시 빼야 한다.
• Dice만 보지 말고 mL 단위의 불일치와 그 위치를 같이 본다.

비교 스크립트

"""usage: python compare.py DATA_DIR OUT_DIR    (DATA_DIR/<case>/ 아래에 DICOM series)"""
import json, statistics, sys, time
from pathlib import Path

import numpy as np
from lungmask import LMInferer
from scipy import ndimage

from lung_threshold import lung_mask, read_ct

THRESHOLDS = [-300, -400, -500, -600, -700, -800]


def dice(a, b):
    return 2 * np.logical_and(a, b).sum() / (a.sum() + b.sum())


def main(data_dir, out_dir):
    out = Path(out_dir)
    out.mkdir(parents=True, exist_ok=True)
    inferer = LMInferer(modelname="R231", force_cpu=True)
    rows = []
    for case in sorted(p for p in Path(data_dir).iterdir() if p.is_dir()):
        img, hu, sp = read_ct(case)
        ml = float(np.prod(sp)) / 1000  # 복셀 1개의 부피 (mL)
        t0 = time.perf_counter()
        ref = inferer.apply(img) > 0
        unet_s = time.perf_counter() - t0
        row = {"case": case.name, "unet_s": round(unet_s, 1),
               "unet_L": round(ref.sum() * ml / 1000, 3), "by_threshold": {}}
        for t in THRESHOLDS:
            times = []
            for _ in range(3):
                t0 = time.perf_counter()
                m = lung_mask(hu, t)
                times.append(time.perf_counter() - t0)
            row["by_threshold"][t] = {
                "dice": round(float(dice(m, ref)), 4), "L": round(m.sum() * ml / 1000, 3),
                "extra_mL": round((m & ~ref).sum() * ml, 1), "missed_mL": round((ref & ~m).sum() * ml, 1),
                "s": round(statistics.median(times), 2)}
        # -400 HU에서 '임계값에만 있는' 부피를 3D 덩어리로 나눠 본다
        m = lung_mask(hu, -400)
        lab, n = ndimage.label(m & ~ref)
        sizes = np.bincount(lab.ravel())[1:]
        row["largest_extra_mL"] = [round(s * ml, 1) for s in sorted(sizes)[::-1][:2]]
        rows.append(row)
        print(json.dumps(row))
    json.dump(rows, open(out / "compare.json", "w"), indent=1)


if __name__ == "__main__":
    main(*sys.argv[1:3])

위 코드는 글에 싣기 위해 줄인 것이다. 실제로 돌린 버전(U-net 결과 캐시, 덩어리별 위치·평균 HU, 가장자리 2픽셀 분석, 구멍 크기 제한 비교 포함)과 그림 코드, 결과 JSON은 함께 공개한 파일에 있다.

참고문헌

1. S. G. Armato 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. LIDC-IDRI 데이터셋, The Cancer Imaging Archive. doi:10.7937/K9/TCIA.2015.LO9QL9SX (CC BY 3.0)
3. K. Clark 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. J. Hofmanninger et al. Automatic lung segmentation in routine imaging is primarily a data diversity problem, not a methodology problem. European Radiology Experimental 4, 50, 2020. doi:10.1186/s41747-020-00173-2 · lungmask (GitHub)
5. SimpleITK 문서