モチベーション

医療画像関連の仕事をしているので、周波数処理とかで、ちょっと拘っているものだったりとか、試してみたくなる。ここでは、ハイブリッド処理を（ほんの気持ち）体験した記録を残す。

従来の周波数強調処理

古典的な画像フィルターや、Unsharpマスク処理、エッジ強調処理などに加えて、マルチ周波数処理、MUSICA（Multi-Scale Image Contrast Amplification）処理、UNIQUE（Unified Image Quality Enhancement）処理、周波数処理とダイナミックレンジ圧縮を組み合わせたハイブリッド処理などがある。

MUSICA、UNIQUEやハイブリッド処理は、処理のアルゴリズムの中に、周波数処理やエッジ抽出、コントラスト調整などの複数の工程が含まれる。理論上、処理対象とする周波数帯域を選択的に操作することができる。

ハイブリッド処理の初歩

プロセスは次の通り。

ボケの加減が異なる画像を多数作成する
隣り合うボケ画像同士の差分画像を得る
ボケ画像の差分の総和を得て、これを高周波成分画像とする
オリジナル画像からボケ画像の総和を差し引くことで、低周波成分画像を作る
高周波成分画像に周波数強調処理を行う（ここでは簡単にラプラシアンフィルターを用いる）
低周波成分画像にダイナミックレンジ圧縮を行う（ここでは簡単にガンマ補正を用いる）
これらを足しあわせる

コード

# 利用ライブラリ（DICOM使うときはpydicomを）
import os
import numpy as np
import matplotlib.pyplot as plt
import cv2

# サンプル画像
# https://pixabay.com/ja/photos/x%E7%B7%9A%E6%92%AE%E5%BD%B1-%E8%A8%BA%E6%96%AD-%E8%A7%A3%E5%89%96%E5%AD%A6-%E6%80%AA%E6%88%91-3057768/

# 正規化関数
def norm(img):
    max = img.max()
    min = img.min()
    return (img - min)/(max - min)

# 画像読み込み
img = cv2.imread(file)
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) # 8-bit
img = np.array(img, dtype=np.float32)
img = norm(img)

def getBlurImage(img, ksize_array):
    # img : 入力画像
    # ksize_array : カーネルサイズ配列
    results = np.zeros((len(ksize_array), img.shape[0], img.shape[1]))
    for i, ksize in enumerate(ksize_array):
        results[i] += cv2.GaussianBlur(img.copy(), (ksize, ksize), 0)
    return results

def getSubtraction(org, blurs):
    # org : オリジナル画像
    # blurs : ガウシアンぼかし画像配列
    results = np.zeros(blurs.shape)
    for i , blur in enumerate(blurs):
        if i == 0:
            results[i] += org - blur
            continue
        results[i] += blurs[i-1] - blur
    return results

def sum_all(subtractions):
    sum_im = np.zeros((subtractions.shape[1], subtractions.shape[2]))
    for subtraction in subtractions:
        sum_im += subtraction
    return sum_im

def frequencyEnhancement_2(sum_im):
    # ラプラシアンフィルタを作成
    laplacian = cv2.Laplacian(sum_im, cv2.CV_64F)
    # 絶対値を取り、データ型を戻す
    laplacian_abs = np.absolute(laplacian)
    # 高周波成分を強調
    sharpened_image = cv2.addWeighted(sum_im, 1.5, laplacian_abs, -0.5, 0)
    return norm(sharpened_image)

def dynamicRangeControl_2(low_freq_img, gamma):
    # ガンマ補正を適用するためのルックアップテーブルを作成
    im = norm(low_freq_img)
    im = np.array(im*255, dtype=np.uint8)
    inv_gamma = 1.0 / gamma
    table = np.array([((i / 255.0) ** inv_gamma) * 255 for i in range(256)], dtype="uint8")
    # 画像にガンマ補正を適用
    return norm(cv2.LUT(im, table))

def showImages(images, row, col):
    for i, image in enumerate(images):
        plt.subplot(row, col, i+1)
        plt.axis("off")
        plt.imshow(image, cmap="gray")
    plt.show()

# 100のアンシャープマスク画像（ガウシアンブラーリングで代用）でテスト
# 奇数を100個生成
# カーネルが小さいほど高周波成分が残る
odd_numbers = np.arange(1, 200, 2)
blurs = getBlurImage(img, odd_numbers)
showImages(blurs, 10, 10)

subtractions = getSubtraction(img, blurs)
showImages(subtractions, 10, 10)

sum_im = sum_all(subtractions)
plt.imshow(sum_im, cmap="gray")
plt.show()

sum_im = sum_all(subtractions)
plt.imshow(sum_im, cmap="gray")
plt.show()

gamma = 1.2 # 大きいほど暗い部分が明るくなる
low_freq_img = img - sum_im
low_freq_img = dynamicRangeControl_2(low_freq_img, gamma)
plt.imshow(low_freq_img, cmap="gray")
plt.show()

# 低周波成分と高周波成分を足しわせる
result = high_freq_img + low_freq_img
result = norm(result)

plt.subplot(1, 2, 1)
plt.axis("off")
plt.title("original")
plt.imshow(img, cmap="gray")
plt.subplot(1, 2, 2)
plt.axis("off")
plt.title("result")
plt.imshow(result, cmap="gray")
plt.show()

References

Yoshida, Makoto & Saeki, Yusuke & Miyai, Masahiro & Shirono, Hiroki & Kato, Katsuya. (2024). Examination of Optimal Frequency Processing in X-ray Images for Occult Femoral Neck Fractures大腿骨頚部不顕性骨折に対する単純X線画像を用いた最適な周波数処理の検討. Japanese Journal of Radiological Technology. 80. 10.6009/jjrt.2024-1447.

最後に

古典的な手法は自分のアイディアの源泉になると思っている。
従来の方法に新しい方法を組み合わせれば、新規性はひねり出せる。

Stay visionary