Pyゼミ1.07 CT画像にDICOM Color LUTを使ってヒートマップ画像などを重ねてFusion画像を作ってみる
今回は何らかの結果（ヒートマップや推論マップなど）とその元画像を重ねるプログラムを作成します。
何らかの結果にはPydicomが提供するとColor LUTを利用してカラー画像に変換します。
Ｐyゼミ1.06ではImageJのColor LUTを使いましたが、ここではPydicomが提供するColor LUTを使います。
Color LUTを適用した画像をDICOM画像に重ねてFusion画像を作成してみます。
医療の世界ではFusionと呼んでいますが、Computer Visonの世界ではSuperimposeと呼んでいます。
具体的には下に示すような画像を作成します（グレイのCT画像に模擬的な結果のカラー画像を重ねています）。
ポイント
- 何らかの結果を疑似的に2次元ガウス分布を用いて作成する。
- pydiocmのLUT適応機能を使ってグレイ画像をカラー画像に変換する。
- グレイ画像のDICOM画像を読み込み、RGBのモノクロ画像（3チャンネル）を作成する。
- OpneCVのaddweighted()メソッドを使ってCT画像と重ね合わせる。
プログラムはここから
サンプルのDICOM画像はここから
まずは必要なライブラリについて
以下のライブラリをインポートします。
import cv2, math
import numpy as np
import pydicom
from pydicom.pixel_data_handlers.util import apply_color_lut
pydicomが提供するLUTを設定する
pydicomは8つのLUTを提供しています。
DICOM画像にLUTを適応するためにはapply_color_lut( )メソッドを使います。
引数にはpydicomが提供するLUTの名称を与える必要があります。
ここでは、このLUTの名称を予めリストlutarrに格納しておきましょう。
lutarr =
["PET", #0
"HOT_IRON", #1
"HOT_METAL_BLUE", #2
"PET_20_STEP", #3
"SPRING", #4
"SUMMER", #5
"FALL", #6
"WINTER"] #7
LUTの名称をindexの値で変数lutnmに格納しておくと、数値を変えるだけでLUTを変更することができて便利です。
ガウス分布を作成する関数gaussianを作成する
この関数は引数にガウス分布の中心の座標、cx、cyと標準偏差sgm、そし生成する画像のサイズsizeを与えます。
はじめに２次元のガウス分布を用いてグレイ画像を作成するためにgaussian()関数を作成します。
# Gaussian
distribution
def
gaussian(cx, cy, sgm, size):
mapimg = np.zeros((size, size))
for h in range(size):
ph =
1/math.sqrt(2.0*math.pi*sgm)*math.exp(-math.pow((h - cy), 2)/2/sgm)
for w in range(size):
pw =
1/math.sqrt(2.0*math.pi*sgm)*math.exp(-math.pow((w-cx), 2)/2/sgm)
mapimg[h,w] = ph*pw
vmax = np.max(mapimg)
print(">>max gaussinan=",
vmax)
mapimg = mapimg / vmax
return mapimg
ガウス分布の中心座標cxとcy、ガウス分布の広がりsgmと生成画像サイズsizeを指定します。
ここでは、以下の値を与えて、gaussian関数を呼びだし、2次元ガウス分布をpmapに代入します。2次元ガウス分布は最大値が1.0になるように正規化されています。
# Set gaussian
condition
size = 512 # map size
x, y = 350,
256 # coordination of Gaussian center
sgm = 500 # Sigma of Gaussian distribution
pmap =
gaussian(x,y,sgm, size)
pmapの最大値は1.0なので、8ビットの最大値255になるように変換して画像に保存します。
pmap *=
255 # convert [0.0~1.0] to
[0.0~255.0]
print(">>Max
=", np.max(pmap))
cv2.imwrite("mapOrg.png",
pmap)
保存されたmapOrg.pngは以下のようになります。
この画像にPydicomが提供するColor LUTを適用してみましょう。
cimg =
apply_color_lut(pmap, palette=lutnm)
cimg =
cimg[ : , : , : : -1]
cv2.imwrite("mapLut.png",
cimg)
カラー画像保存する際、BGRをRGBに変換するので cimg = cimg[ : , : , : : -1] とします。
結果の画像mapLut.pngは次のように表示されます。
グレイスケールを画像に入れる
画素値と色の関係がよくわからないので、グレイスケールを左に入れてみましょう。
このグレイスケールはColor LUTを適用するとカラースケールになります。
# Add color
scale
for y in
range(256):
pmap[350-y, :15] = y
cimg =
apply_color_lut(pmap, palette=lutnm)
cimg =
cimg[ : , : , : : -1]
cv2.imwrite("mapLut2.png",
cimg)
このスケールはpmapの縦350の位置から上に向かって画素値0から255を代入します。
スケールの幅は15画素です（:15で指定しています)。
結果mapLut2.pngは以下のように表示されます。
ＤICOM画像を読み込む
それではDICOM画像Brain01を読み込んでみます。
# Read Dicom
file
dcmfnm =
"../dcmdir1/Brain01"
ds =
pydicom.dcmread(dcmfnm)
print(">>ds
=",ds)
wc = ds.WindowCenter
ww = ds.WindowWidth
img =
ds.pixel_array
max = wc +
ww/2
min = wc -
ww/2
img = (img -
min)/(max - min) * 255
img[img >
255] = 255
img[img <
0] = 0
img =
img.astype('uint8') # float32 → uint8
print(">>
image shape=",img.shape, img.dtype)
# Save png
image and read png image
cv2.imwrite("./brain01.png",
img)
DICOMのタグ情報を読み込みウインドウセンタwcとウインドウ幅wwを使ってウインドニング処理を行います。
また、デフォルトfloat32をuint8の符号なし8ビットに変換しているます。
この画像のshapeとタイプを表示すると次のように表示されます。
>1> image
shape= (512, 512) uint8
shapeは(512,512)と表示され、512×512のグレイ画像であることを示しています。
Fusionを行うためには2つの画像がカラー画像であることが必要です。
グレイ画像をカラー画像に変換する３つの方法
方法1：保存した画像をimread()メソッドで読み込む。
CV2のimreadのデフォルトはグレイ画像もカラー画像として読み込まれます。
img2 = cv2.imread("./brain01.png")
print(">2> image2 shape=",
img2.shape, img2.dtype)
結果は以下のように、読み込まれた画像img2が3チャンネルである、つまりRGBのカラー画像であることを示しています。
>2>
image2 shape= (512, 512, 3) uint8
方法2：cv2のmerge()メソッドを使います。
# Merge 3 imags
img3 =
cv2.merge((img, img, img))
print(">3>
image3 shape=", img3.shape, img3.dtype)
RGBに同じimgをマージして3チャンネルの画像を作成しています。
結果は次のように、確かに512×512画素の3チャンネルの画像img3ができています。
>3>
image3 shape= (512, 512, 3) uint8
方法3：numpyのzerosメソッドを使って3チャンネルの画像を作成する。
# Assing image
to 3 chanel image
img4 = np.zeros((512,
512, 3), dtype=np.uint8)
img4[ : , : , 0] = img
img4[ : , : , 1] = img
img4[ : , : , 2] = img
print(">4>
image4 shape=", img4.shape, img4.dtype)
あらかじめ作成したimg4の0、1、2の各チャンネルに画像imgを代入します。
実行結果は3チャンネルの画像img4ができています。
>4>
image4 shape= (512, 512, 3) uint8
画像にカラースケールの背景を入れる
3チャンネルに変換した画像img2にカラースケールの発色をよくするため画素値128で埋めます。
これを入れないとカラースケールが暗くなります。
# For color
scale
img2[350-256 : 350,
: 15] = 128
２つの画像を重ね合わせる
それではaddWeighted()メソッドを使って２つの画像img2とcimgを重ね合わせます。
# Create fusion
image with addWeighted()
rimg =
cv2.addWeighted(src1=img2, alpha=1.0, src2=cimg, beta=0.5, gamma=0)
cv2.imwrite("fusion.png",
rimg)
結果画像fusion.pngは次のように表示されます。
addWeighted()メソッドの引数alphaとbetaは２つの画像の重ねるときの重みになります。
主たる画像、ここではBrain
CT画像のimg2は1.0とし、ガウス分布画像のcimgは0.5と重みを小さくして重ね合わせています。
bataの値については場合により調整してみましょう。
Post Script:
実はaddWeighted()メソッドを使った画像の重ね合わせは結構はまりました。
大前提は2つの画像はカラー画像であることです。
しかし、DICOM画像をpydicomで読み込んだときは１チャンネルの画像、その画像をいったん保存し、imread()メソッドで読み込むとデフォルトでカラー画像として読み込まれます。
したがって、ある時はaddWeighted()で上手く重なり、あるときは上手くいかない。
この原因を明らかにするのにはまってしまいました。気を付けましょう。
