Bias Field Correction에 대하여

의료 영상 분석에서 정확한 intensity(강도) 표현은 매우 중요합니다.

하지만 MRI 영상에는 종종 불균일한 밝기 분포, Bias Field (또는 Intensity Non-uniformity)라는 것이 존재하는데, 이런 요소는 제거해주지 않으면 후처리 과정에 방해가 됩니다.

 

이번 포스팅에서는 이 bias field가 무엇인지, 왜 발생하는지, 그리고 어떻게 bias field correction을 수행하는지에 대해 알아보겠습니다.

---

Bias Field란?

Bias Field란 실제 조직의 물리적 신호가 아닌 기계적, 물리적 원인에 의해 덧붙여진 저주파 성분을 말합니다.

 

MRI 영상을 보다보면, 어느 부위만 전반적으로 어두워져서 실제로 같은 조직인데도 intensity가 다른 시리즈를 볼 수 있습니다.

 

이러한 영상은 당연히 실제 조직의 intensity가 왜곡되었기 때문에 segmentation, registration, intensity-based analysis 등과 같은 작업에서 오류 유발할 수 있는 요소가 될 것입니다.

---

Bias Field 발생 원인

1. Radio Frequency Excitation Field (B₁⁺ field)

• MRI는 인체에 고주파(RF)를 쏴서 수소 원자를 자극함
• 이 자극 강도가 몸 전체에 고르게 퍼지지 않으면, 같은 조직도 다른 밝기로 보임

→ MRI가 신호를 보내는 쪽(RF 송신 필드)의 불균일성

 

2. Magnetic Field Inhomogeneity (B₀ 불균일성)

• MRI의 기본 자기장(예: 1.5T, 3T)이 이론적으로는 균일해야 하지만, 실제로는 기계 구조, 체형, 주변 금속 등에 따라 공간마다 강도가 달라짐
  → 특히 고자기장(3T 이상)에서 더 심하게 나타남

→ 기본 자기장(B₀)이 공간마다 다름

 

3. Non-uniform Reception Coil Sensitivity (수신 코일 민감도 불균일)

• 인체 주변에 여러 개의 수신 코일(phased-array coil)이 배치됨
• 코일에 가까운 부위는 신호가 강하게, 먼 부위는 약하게 수신됨

→ MRI 신호를 받아들이는 코일의 감도 차이

---

MRI 영상의 수학적 모델

Bias field correction 원리에 대해 알아보기 전에, 먼저 MRI 영상의 수학적 모델에 대해 짚고 넘어가야 합니다.

$$I(x) =B(x)*T(x)+n(x)$$

$I(x)$	획득된 영상 (bias 포함)
$T(x)$	True tissue intensity
$B(x)$	bias field (저주파 밝기 왜곡, multiplicative)
$n(x)$	Noise

 

Bias field correction은 간단히 말하면 $I(x)$에 추정된 bias field $B(x)$를 나눠줘서 실제 조직 intensity $T(x)$를 복원하는 작업입니다.

---

Bias Field Correction 원리

1. Input Image

먼저 우리가 보정하고 싶은 MRI 영상을 불러옵니다.

---

2. Mask Generation

이미지 전체가 아닌, 뇌처럼 실제로 의미 있는 조직 부분만 보정하고 싶기 때문에 마스크를 생성해야합니다.

• 사용자가 직접 주는 방식
• Otsu 알고리즘을 통해 foreground(조직)와 background(공기 등)를 나눠서 생성하는 방식

---

3. Shrink Image (optional)

계산 속도를 빠르게 하기 위해 이미지 해상도를 줄입니다. (예: 1/2 크기) 

실험해보니, 결과는 아래와 같았습니다. (26 slices)

• shrink 1 : 76.49 sec
• shrink 2 : 22.15 sec
• shrink 4 : 10.21 sec

⇒ 퀄리티는 당연하게도 shrink를 많이 할 수록 떨어졌구요.

---

4. Multi-resolution Pyramid

이미지를 해상도에 따라 여러 단계로 축소된 버전으로 만들어, coarse-to-fine 방식 (작은 해상도부터 시작해서 고해상도로 올라가는 방식)으로 점진적으로 보정을 수행합니다.

• 낮은 해상도에서는 큰 구조나 패턴을 빠르게 감지
• 높은 해상도에서는 세부 디테일을 보정

coarse level에서 찾은 bias field를 다음 level의 초기값으로 넘겨줌으로써 효율적으로 계산을 수행하게 됩니다.

---

5. B-spline Bias Field Estimation (BFGS 최적화)

Bias Field는 부드럽게 변화하는 함수이기 때문에, B-spline 곡선으로 모델링하는 것이 적합합니다.

$$logB(x)=∑c_k⋅ϕ_k(x)$$

• $logB(x)$ : bias field (로그 공간에서 다룸)
• $c_k$ : control points (coefficients) → 추정해야 하는 값들
• $ϕ_k(x)$ : B-spline basis functions

즉, 전체 bias field는 basis function들의 weighted sum이 된다는 것입니다.

 

그런 다음 BFGS라는 수학적 기법으로 가장 적절한 bias field를 탐색합니다.

이것은 일종의 최적화 문제인데요, 조직의 intensity가 일관되게 보이도록 만드는 것(= 같은 조직이면 영상 어디에 있든 비슷한 밝기가 되도록)이 목적인 최적화 문제라고 볼 수 있습니다.

 

최적화 방식은 아래와 같습니다.

1. $\log I(x) - \sum c_k \cdot \phi_k(x)$로 residual image 계산
2. 이 residual image의 histogram이 sharp한 peak를 가지도록 계수 $c_k$ 최적화

objective funcion은 residual histogram의 sharpness(entropy ↓ 또는 kurtosis ↑)가 될 것입니다.

---

6. Log Bias Field → Exp

추정된 bias field는 로그값 형태 $logB(X)$이므로, 지수 함수(exp)를 씌워서 실제 bias field $B(x)$를 구해줍니다.

---

7. Bias Field Correction

$$\text{Corrected Image} = \frac{\text{Input Image}}{\text{Bias Field}}$$

⇒ 각 픽셀에서 bias를 제거함으로써 원래 조직의 밝기를 복원합니다.

---

코드 예시

import os
import SimpleITK as sitk
import time


def apply_n4_bias_correction(
    input_path,
    output_path,
    log_bias_output_path,
    shrink_factor=2,
    mask_path=None,
    num_iterations=50,
    num_fitting_levels=4,
):
    start_time = time.time()

    # 1. 입력 이미지 (float32)
    input_image = sitk.ReadImage(input_path, sitk.sitkFloat32)
    image = input_image

    # 2. 마스크: 지정 or Otsu 자동 생성
    if mask_path and os.path.exists(mask_path):
        mask_image = sitk.ReadImage(mask_path, sitk.sitkUInt8)
    else:
        mask_image = sitk.OtsuThreshold(input_image, 0, 1, 200)

    # 3. Shrink
    if shrink_factor > 1:
        image = sitk.Shrink(input_image, [shrink_factor] * input_image.GetDimension())
        mask_image = sitk.Shrink(mask_image, [shrink_factor] * input_image.GetDimension())

    # 4. 필터 생성 및 파라미터 설정
    corrector = sitk.N4BiasFieldCorrectionImageFilter()
    corrector.SetMaximumNumberOfIterations([num_iterations] * num_fitting_levels)

    # 5. 보정 실행
    corrector.Execute(image, mask_image)

    # 6. log bias field 얻기
    log_bias_field = corrector.GetLogBiasFieldAsImage(input_image)
    corrected_image_full = input_image / sitk.Exp(log_bias_field)

    # 7. 결과 저장
    sitk.WriteImage(corrected_image_full, output_path)
    # sitk.WriteImage(log_bias_field, log_bias_output_path)

    print(f"Saved corrected image to: {output_path}")
    # print(f"Saved log bias field to: {log_bias_output_path}")

    elapsed = time.time() - start_time
    print(f"Bias field correction took {elapsed:.2f} seconds")

    return corrected_image_full, log_bias_field


if __name__ == "__main__":
    # 경로 설정
    input_file = "/home/ykseo/main/test_code/image_preprocessing/brain/bias_field_correction/input/T2_NECT_preprocessed_T2.nii.gz"
    mask_file = "/home/ykseo/main/test_code/image_preprocessing/brain/bias_field_correction/input/T2_NECT_preprocessed_T2_mask.nii.gz"
    output_file = "/home/ykseo/main/test_code/image_preprocessing/brain/bias_field_correction/output/T2_NECT_preprocessed_T2_corrected.nii.gz"
    log_bias_output_file = "/home/ykseo/main/test_code/image_preprocessing/brain/bias_field_correction/output/T2_NECT_preprocessed_T2_log_bias_field.nii.gz"

    apply_n4_bias_correction(
        input_path=input_file,
        output_path=output_file,
        log_bias_output_path=log_bias_output_file,
        shrink_factor=2,
        mask_path=mask_file,
        num_iterations=100,
        num_fitting_levels=4,
    )

파라미터	설명	값 예시	높이면	낮추면
shrink_factor	이미지 축소 배율	1~4	속도↑, 정밀도↓	속도↓, 정밀도↑
num_iterations	단계당 반복 횟수	30~100	정밀도↑, 시간↑	빠르지만 부정확
num_fitting_levels	해상도 단계 수	3~5	global→local 정밀 보정	coarse bias만 제거

---

결과 예시

왼쪽: 원본, 오른쪽: 보정

왼쪽: 원본, 오른쪽: 보정