CTを対象としたMetal Artifact Reduction（MAR）の雰囲気を味わう記録

ビビッときた

これは一度体験しておかないといけないと思ったから、記録します。

データ

gammex 467 : https://ro-journal.biomedcentral.com/articles/10.1186/1748-717X-9-114#additional-information
abd : Case courtesy of Dinesh Brand, Radiopaedia.org, rID: 67631

免責

利用させていただくデータは16 bit 画像でないので、HU値で処理できません。なので、16 bit画像なら、もっとこう(HU値を使ったしきい値処理など)できてより良い検討ができるんじゃないかなってことを知るために、まずは8 bitで試してみた初歩的な内容になります。
また、MAR界隈で利用されているSinogramへの画素補間や、しきい値処理は、もっと複雑な工程があるはずです。例えば、Contourで関心領域を分けたり、処理の過程で途中でしきい値処理を挟んだりです。こういった工程は今回は無視しています。
（どなたか詳しい方、教えていただけましたら幸いです。）

環境

GoogleDrive+Colaboratory

コード

まずはDriveとつなげます。（画像がDriveに入っているという前提でいきます。）

from google.colab import drive
drive.mount('/content/drive')

利用するライブラリとパッケージです。

import os
import numpy as np
import matplotlib.pyplot as plt

from skimage import io,img_as_float
from skimage.transform import radon, iradon
from scipy.interpolate import interp1d

from sklearn.impute import KNNImputer
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer

Driveをお使いの方（かた）は、、、

os.chdir('画像データが保存されているフォルダへ')

データ

データはこのようになります。

# 論文から8 bit画像を画面キャプチャーでお借りました。
gammex 467 : https://ro-journal.biomedcentral.com/articles/10.1186/1748-717X-9-114#additional-information
# リンクからダウンロードしました。
adb : Case courtesy of Dinesh Brand, Radiopaedia.org, rID: 67631

これらを読み込みます。

gam_path = 'Gammex467.tif'
abd_path = 'ct-metal-artifact-reduction-single-energy-metal-artifact-reduction-semar-algorithm.jpg'
semar_path = 'ct-metal-artifact-reduction-single-energy-metal-artifact-reduction-semar-algorithm_SEMAR.jpg'

画像を読み込みます。

# 8 bit image (本当は16-bitで検討したいのですが、データがないので)
'''
gam = io.imread(gam_path) # numpy.ndarray
gam = img_as_float(gam)
abd = io.imread(abd_path)
abd = img_as_float(abd)
semar = io.imread(semar_path)
semar = img_as_float(semar)

正規化関数です。
今回は簡略化して、MinMaxスケールでいきます。

def normalize(image):
    im = image.copy()
    im = (im - np.min(im)) / (np.max(im) - np.min(im))
    return np.clip(im, 0.0, 1.0)

しきい値処理の関数です。
しきい値未満は0.0、以上はそのままにします。

def clip(image, threshold):
    im = image.copy()
    clipped = np.where(im<threshold, 0.0, im)
    return clipped

画像とサイノグラムを表示する関数です。

def show_image(image, sinogram):
    plt.subplot(1,2,1)
    plt.imshow(image, cmap='gray')
    plt.subplot(1,2,2)
    plt.imshow(sinogram, cmap='gray')
    plt.show()

サイノグラムの画素補間のための関数です。
(本当は、縦横だけでなく、斜めとかからも補間できるとよいのだろうか)

def interpolate_nan_columns(arr, kind='linear'):
    '''
    arr : adarray
    kind: linear, nearest, cubic
    '''
    for j in range(arr.shape[1]):
        y = arr[:, j]
        nans = np.isnan(y)
        if nans.all():
            continue  # 全ての値がNaNの場合、補間しない
        x = np.arange(len(y))
        f = interp1d(x[~nans], y[~nans], kind=kind, fill_value='extrapolate')
        arr[nans, j] = f(x[nans])
    return arr

def interpolate_nan_rows(arr, kind='linear'):
    for i in range(arr.shape[0]):
        y = arr[i, :]
        nans = np.isnan(y)
        if nans.all():
            continue  # 全ての値がNaNの場合、補間しない
        x = np.arange(len(y))
        f = interp1d(x[~nans], y[~nans], kind=kind, fill_value='extrapolate')
        arr[i, nans] = f(x[nans])
    return arr

画素補間でなく、欠損データの補間として扱うこともできます。
例えば、kNNやMICEなどですね。
ただ、今回はあまりうまく行かなかったので、紹介のみです。

def knn_impute(arr, nn=5):
    imputer = KNNImputer(n_neighbors=nn)
    return imputer.fit_transform(arr)

MARの関数です。
私の気分で書いていますので、アルゴリズムは、実際と違うところが含まれると思います。

def calc_mar(image, threshold, filter='ramp', interp='cubic'):
    '''
    filter
        ramp        高解像度が必要な場合
        shepp-Logan	バランス重視の再構成
        cosine      ノイズの多いデータ
        hamming	    ノイズと解像度のバランス
        hann	    非常にノイズの多いデータ
    interp
        ‘linear’, ‘nearest’, and ‘cubic’
    '''
    # preprocess
    # 正規化（MinMaxスケール）
    im_norm = normalize(image)
    # メタル画像の取得
    clipped = clip(im_norm, threshold)
    '''
    calc sinogram0 from org image
    '''
    # ラドン変換 180度の投影
    theta = np.linspace(0.0, 180.0, max(image.shape), endpoint=False)
    # org sinogram
    sinogram0 = radon(im_norm, theta=theta, circle=False)
    # metal sinogram
    sinogram1 = radon(clipped, theta=theta, circle=False)

    # original
    show_image(im_norm, sinogram0)
    print()
    # clip
    show_image(clipped, sinogram1)
    print()

    '''
    calculate no metal sinogram
    '''
    sinogram2 = sinogram0 - sinogram1
    # 改善の余地あり
    # 殆どがバックグラウンドな信号なので、簡易で中央値未満をNaNに。
    th = np.percentile(sinogram2, 50)
    # 0.0はバックグラウンド
    # 小さい値はアーチファクト領域など
    sinogram2[(sinogram2 > 0.0) & (sinogram2 < th)] = np.nan

    # interpolation
    sinogram2 = interpolate_nan_columns(sinogram2, kind=interp)
    sinogram2 = interpolate_nan_rows(sinogram2, kind=interp)

    # impute by knn
    # sinogram2 = knn_impute(sinogram2)

    # impute by imputer
    # imp_mean = IterativeImputer(max_iter=5, min_value=th, random_state=0)
    # imp_mean.fit(sinogram2)
    # sinogram2 = imp_mean.transform(sinogram2)

    '''
    reconstruct
    '''
    # filter : ramp, shepp-logan, cosine, hamming, hann
    # interpolation: ‘linear’, ‘nearest’, and ‘cubic’
    recon_nometal = iradon(sinogram2, theta=theta, filter_name=filter, interpolation='cubic', circle=False)

    # center crop (if you use circle=True, no need following cutout.)
    reconstructed_center = np.array(recon_nometal.shape) // 2
    original_center = np.array(image.shape) // 2
    start = reconstructed_center - original_center
    end = start + image.shape
    recon_nometal = recon_nometal[start[0]:end[0], start[1]:end[1]]

    # no metal
    show_image(recon_nometal, sinogram2)

    '''
    reduction
    '''
    sinogram3 = sinogram0 + sinogram2
    sinogram3 = normalize(sinogram3)
    recon_suppress = iradon(sinogram3, theta=theta, filter_name=filter, interpolation='cubic', circle=False)
    recon_suppress = recon_suppress[start[0]:end[0], start[1]:end[1]]
    recon_suppress = normalize(recon_suppress)

    show_image(recon_suppress, sinogram3)

    '''
    最後にメタル画像を加算
    '''
    mar_im = recon_suppress+(clipped)/2 # すべて加算すると強すぎるかもなので

    return mar_im

結果

gammex 467

res = calc_mar(gam, 0.9)
show_image(gam, res)

元画像

高信号領域のみの画像

差分画像

合成画像

MAR画像（左：ORG、右：処理結果）

あまり、アーチファクト消えていない（泣）

abd1

res1 = calc_mar(abd, 0.9)
show_image(gam, res1)

元画像

高信号領域のみの画像

差分画像

合成画像

MAR画像（左：ORG、右：処理結果）

サイノグラムの中心部の穴を補間でうまいことなくせたら、もっとうまくいくのだろう。

abd2(SEMAR済み)

元画像

高信号領域のみの画像

差分画像

合成画像

MAR画像（左：ORG、右：処理結果）

なんか、「おっ」ってなるけど、黒くなっただけ感あります。

追記

サイノグラムに、inpaint_biharmonic()を使うと、いい場合がある。
しかし、やっぱり、歯とか、人工股関節とかの複雑なメタルアーチファクトだと、マスクを広げすると、ぼかしすぎる感じがする。気持ち細めにマスクをつくると、いい感じ。

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
from skimage.transform import radon, iradon, rescale
from skimage.data import shepp_logan_phantom
from skimage.restoration import inpaint_biharmonic # 追加

# --- 1. 準備: 補間関数の定義 ---
def interpolate_sinogram_by_kind(sinogram, mask, kind='cubic'):
    # (上記の interpolate_sinogram_by_kind 関数をここにコピー)
    inpainted_sinogram = np.copy(sinogram)
    num_angles, num_detectors = sinogram.shape
    detector_indices = np.arange(num_detectors)
    min_points_for_kind = {'linear': 2, 'quadratic': 3, 'cubic': 4}
    required_points = min_points_for_kind.get(kind, min_points_for_kind['cubic'])

    for i in range(num_angles):
        row_data = sinogram[i, :]
        row_mask = mask[i, :]
        valid_indices = detector_indices[~row_mask]
        valid_values = row_data[~row_mask]
        nan_in_valid_values = np.isnan(valid_values)
        if np.any(nan_in_valid_values):
            valid_indices = valid_indices[~nan_in_valid_values]
            valid_values = valid_values[~nan_in_valid_values]

        current_kind_effective = kind
        if len(valid_indices) < required_points:
            if len(valid_indices) >= min_points_for_kind['linear']:
                current_kind_effective = 'linear'
            elif len(valid_indices) == 1:
                inpainted_sinogram[i, row_mask] = valid_values[0]
                continue
            else:
                inpainted_sinogram[i, row_mask] = 0
                continue
        
        final_required_points = min_points_for_kind.get(current_kind_effective, 2)
        if len(valid_indices) < final_required_points:
             if len(valid_indices) == 1: inpainted_sinogram[i, row_mask] = valid_values[0]
             else: inpainted_sinogram[i, row_mask] = 0
             continue

        fill_val = 0.0
        if len(valid_values) > 0: # 境界値補外のための値
            fill_val = (valid_values[0], valid_values[-1])

        interp_func = interp1d(
            valid_indices, valid_values, kind=current_kind_effective,
            bounds_error=False, fill_value=fill_val
        )
        indices_to_interpolate = detector_indices[row_mask]
        if len(indices_to_interpolate) > 0:
            inpainted_values = interp_func(indices_to_interpolate)
            inpainted_sinogram[i, indices_to_interpolate] = inpainted_values
    return inpainted_sinogram


# --- 2. ファントムデータの準備 (前回と同様) ---
image_size = 128
theta = np.linspace(0., 180., image_size, endpoint=False)
phantom_true_original_hr = shepp_logan_phantom()
phantom_true_original = rescale(phantom_true_original_hr, scale=image_size / phantom_true_original_hr.shape[0], mode='reflect', anti_aliasing=False)
metal_value = 2.0
phantom_with_metal = np.copy(phantom_true_original)
metal_size_x = 5; metal_size_y = 10
metal_pos_x_start = image_size // 2 + 15; metal_pos_y_start = image_size // 2 - 20
phantom_with_metal[metal_pos_y_start:metal_pos_y_start + metal_size_y, metal_pos_x_start:metal_pos_x_start + metal_size_x] = metal_value

# --- 3. 金属トレースマスクの作成 (前回と同様) ---
metal_threshold = 1.5
metal_segmentation_mask_image = phantom_with_metal >= metal_threshold
sinogram_metal_trace = radon(metal_segmentation_mask_image.astype(float), theta=theta, circle=True)
sinogram_metal_mask = sinogram_metal_trace > 1e-3

# --- 4. サイノグラムの生成と各種補間 ---
sinogram_corrupted = radon(phantom_with_metal, theta=theta, circle=True)
sinogram_true_original = radon(phantom_true_original, theta=theta, circle=True)

# 4a. 線形補間
inpainted_sinogram_linear = interpolate_sinogram_by_kind(sinogram_corrupted, sinogram_metal_mask, kind='linear')

# 4b. 3次スプライン補間 (平滑化なし)
inpainted_sinogram_cubic = interpolate_sinogram_by_kind(sinogram_corrupted, sinogram_metal_mask, kind='cubic')

# 4c. biharmonic inpainting
# マスク部分はinpaint_biharmonicによって上書きされる
inpainted_sinogram_biharmonic = inpaint_biharmonic(
    sinogram_corrupted,
    sinogram_metal_mask,
    multichannel=False # For grayscale
)


# --- 5. 画像の再構成 ---
filter_name = 'ramp'
reconstructed_true_original = iradon(sinogram_true_original, theta=theta, circle=True, filter_name=filter_name)
reconstructed_corrupted_no_mar = iradon(sinogram_corrupted, theta=theta, circle=True, filter_name=filter_name)
reconstructed_linear_mar = iradon(inpainted_sinogram_linear, theta=theta, circle=True, filter_name=filter_name)
reconstructed_cubic_mar = iradon(inpainted_sinogram_cubic, theta=theta, circle=True, filter_name=filter_name)
reconstructed_biharmonic_mar = iradon(inpainted_sinogram_biharmonic, theta=theta, circle=True, filter_name=filter_name)

# --- 6. 結果の比較表示 ---

# サイノグラムの表示
fig_sino, axes_sino = plt.subplots(1, 5, figsize=(25, 5), sharex=True, sharey=True) # 5つのサイノグラムを表示
s_vmin = np.min(sinogram_true_original); s_vmax = np.max(sinogram_true_original)

axes_sino[0].imshow(sinogram_corrupted, cmap=plt.cm.Greys_r, aspect='auto', vmin=s_vmin, vmax=s_vmax)
alpha_mask_sino = np.zeros((*sinogram_metal_mask.shape, 4)); alpha_mask_sino[sinogram_metal_mask] = [1,0,0,0.3]
axes_sino[0].imshow(alpha_mask_sino, aspect='auto')
axes_sino[0].set_title('Corrupted (Mask Red)')

axes_sino[1].imshow(inpainted_sinogram_linear, cmap=plt.cm.Greys_r, aspect='auto', vmin=s_vmin, vmax=s_vmax)
axes_sino[1].set_title('Linear Inpainting')

axes_sino[2].imshow(inpainted_sinogram_cubic, cmap=plt.cm.Greys_r, aspect='auto', vmin=s_vmin, vmax=s_vmax)
axes_sino[2].set_title('Cubic Spline Inpainting')

axes_sino[3].imshow(inpainted_sinogram_biharmonic, cmap=plt.cm.Greys_r, aspect='auto', vmin=s_vmin, vmax=s_vmax)
axes_sino[3].set_title('Biharmonic Inpainting')

axes_sino[4].imshow(sinogram_true_original, cmap=plt.cm.Greys_r, aspect='auto', vmin=s_vmin, vmax=s_vmax)
axes_sino[4].set_title('Ideal Sinogram')

for ax in axes_sino: ax.set_xlabel('Projection Angle'); ax.set_ylabel('Detector Position')
fig_sino.tight_layout()
plt.show()

# 再構成画像の表示
fig_recon, axes_recon = plt.subplots(1, 5, figsize=(25, 5), sharex=True, sharey=True) # 5つの再構成画像
r_vmin = np.min(reconstructed_true_original) * 0.9; r_vmax = np.max(reconstructed_true_original) * 1.1

axes_recon[0].imshow(reconstructed_corrupted_no_mar, cmap=plt.cm.Greys_r, vmin=r_vmin, vmax=r_vmax); axes_recon[0].set_title('No MAR')
axes_recon[1].imshow(reconstructed_linear_mar, cmap=plt.cm.Greys_r, vmin=r_vmin, vmax=r_vmax); axes_recon[1].set_title('Linear MAR')
axes_recon[2].imshow(reconstructed_cubic_mar, cmap=plt.cm.Greys_r, vmin=r_vmin, vmax=r_vmax); axes_recon[2].set_title('Cubic Spline MAR')
axes_recon[3].imshow(reconstructed_biharmonic_mar, cmap=plt.cm.Greys_r, vmin=r_vmin, vmax=r_vmax); axes_recon[3].set_title('Biharmonic MAR')
axes_recon[4].imshow(reconstructed_true_original, cmap=plt.cm.Greys_r, vmin=r_vmin, vmax=r_vmax); axes_recon[4].set_title('Ideal (No Metal)')

for ax in axes_recon: ax.axis('off')
fig_recon.tight_layout()
plt.show()

# 特定の行（角度）のプロットで詳細比較
angle_to_plot = sinogram_corrupted.shape[0] // 2
if np.sum(sinogram_metal_mask[angle_to_plot,:]) == 0 and sinogram_corrupted.shape[0] > 20:
    for i_angle_search in range(sinogram_corrupted.shape[0]):
        if np.sum(sinogram_metal_mask[i_angle_search,:]) > 0 : angle_to_plot = i_angle_search; break

plt.figure(figsize=(14, 7))
plt.plot(sinogram_corrupted[angle_to_plot, :], 'o', markersize=3, label='Corrupted Data (Input)', alpha=0.5)
plt.plot(inpainted_sinogram_linear[angle_to_plot, :], 'r.-', label='Linear Inpainted', alpha=0.7)
plt.plot(inpainted_sinogram_cubic[angle_to_plot, :], 'm.--', label='Cubic Spline Inpainted', alpha=0.7)
plt.plot(inpainted_sinogram_biharmonic[angle_to_plot, :], 'c-', label='Biharmonic Inpainted', linewidth=2)
plt.plot(sinogram_true_original[angle_to_plot, :], 'k:', label='Ideal (No Metal)', linewidth=2, alpha=0.6)

masked_indices_row = np.where(sinogram_metal_mask[angle_to_plot, :])[0]
y_min_plot, y_max_plot = plt.ylim()
if len(masked_indices_row)>0:
    for idx in masked_indices_row:
        plt.fill_betweenx([y_min_plot, y_max_plot], idx-0.5, idx+0.5, color='red', alpha=0.05, label='Masked Region' if idx == masked_indices_row[0] else "")

plt.title(f'Sinogram Data at Angle {angle_to_plot} - Comparison of Inpainting Methods')
plt.xlabel('Detector Bin'); plt.ylabel('Intensity')
plt.legend(fontsize=9); plt.grid(True); plt.ylim(y_min_plot, y_max_plot)
plt.show()

References

• Ballhausen, H., Reiner, M., Ganswindt, U. et al. Post-processing sets of tilted CT volumes as a method for metal artifact reduction. Radiat Oncol 9, 114 (2014). https://doi.org/10.1186/1748-717X-9-114
• Brand D, CT metal artifact reduction - single energy metal artifact reduction (SEMAR) algorithm. Case study, Radiopaedia.org (Accessed on 26 Jun 2024) https://doi.org/10.53347/rID-67631
• Ichikawa K, Kawashima H, Takata T. An image-based metal artifact reduction technique utilizing forward projection in computed tomography. Radiol Phys Technol. 2024 Jun;17(2):402-411. doi: 10.1007/s12194-024-00790-1. Epub 2024 Mar 28. PMID: 38546970; PMCID: PMC11128408.

終わり

免責で書きましたが、ここでは、8-bit画像で試した結果のみを記録しました。本来、CTは16-bitで、かつ、HU値で操作できるので、16-bitでトライしたら、また違う結果になると思います。機会があればトライしたいです。

共同研究など、お声がけお待ちしております。

Stay visionary