CTのグレー画像をLUTを使ってカラー画像にしてみる。
keyword: DICOM, ImageJ, LUT(Lookup table)
CTやMR画像などはグレー画像で色がありません。
これに色を付けるには１つの画素値を赤（R)，緑（G)と青（B)の３つの色に分ける必要があります。
グレーの画素値（0から255）をRGBに分解する対応表をLookup table（LUT）と呼びます。
ImageJのLUTでCT画像に色を付ける
画像処理ツールのImageJにはたくさんのLUTが揃っています。
 |
| ImageJのImage/Lookup TablesにたくさんのLUTを見ることができます。 |
これを使ってImageJでグレー画像に疑似カラーの色を付けることができます。
オリジナルのDICOM画像を表示します。
 |
| ImageJによるDICOM画像表示 |
この画像に「16 Color」のLUTを使ってカラー化してみます。
 |
| LUT処理後のDICOM画像 |
ImageJのLUTはImageJがインストールされたディレクトリImageの下のlutsのディレクトリの中に .lut ファイルとして存在します。
ImageJのLUTファイルを読み込む
プログラムreadImageJLut.pyは初めにLUTを入れるリストのためにnumpyとLUTを表示するためにmatplotlibをインポートします。
LUTファイル名を変数lutfnmに代入します。
LUTがファイルをOpenしてdtにLUTのデータを”rb”で，バイナリデータとして読み込みます。
f = open(lutfnm, "rb")
dt = f.read()
f.close()
dtにはRGBの値が1次元のバイナリデータとして格納されているため整数型のnumpyのリストに変換します。
for文を使ってdtの要素を一つずつ取り出しリストrgbに入れます。
rgb = np.array( [ i for i in dt ] )
リストrgbにはRの成分256個，Gの成分56個，Bの成分256個が1次元に格納されているので，reshapeメソッドを使って，3ｘ256の2次元リストに変換しています。
これをR,G,Bの各成分のリストr，ｇ，ｂに代入します。
r, g, b = np.reshape( rgb, (3,256) )
これで各画素値に対応するRGB成分がリストr,g,bに格納されました。
グレーの画素値がvのとき，各RGB成分は r[v], g[v], b[v]として取り出すことができます。
readImageJLut.py
import numpy as np
from matplotlib import pyplot as plt
lutfnm = "D:\Dropbox\Temporary\ImageJ\ij152-win-java8\ImageJ\luts\\16_Colors.lut"
f = open(lutfnm, "rb")
dt = f.read()
f.close()
print(">1>要素数：", len(dt), ",型:",type(dt))
rgb = np.array([i for i in dt]) # [r0・・・ｒ２５５ g0・・・g255 b0・・・b255]
r, g, b = np.reshape(rgb,(3,256)) # 赤，緑と青の成分をリストr,gとbに分解
print(">2>",len(r),len(g),len(b))
# 画素値（Gray）とRGB成分の表示
print("LUT\nGray\tRed\tGreen\tBlue")
for i in range(256):
print(i,"\t", r[i],"\t", g[i],"\t", b[i])
# LUT値のグラフ表示
x = np.linspace(0,255,256) # 0〜255を256分割する
plt.plot(x, r, color='red')
plt.plot(x, g, color='green')
plt.plot(x, b, color='blue')
plt.show()
プログラムを実行してみよう
プログラムを実行するとはじめに次の2行が表示されます。
>1>要素数： 768 ,型: <class 'bytes'> >2> 256 256 256
要素数は768と表示されています。これは256ｘRGBになります。
そしてImageJのLUTはバイナリデータでファイルに格納されているためバイト型になっています。
次にバイト型が整数型のリストに変換され，LUTが表示されます。
LUT Gray Red Green Blue 0 0 0 0 1 0 0 0 2 0 0 0 3 0 0 0 4 0 0 0 5 0 0 0 6 0 0 0 7 0 0 0 8 0 0 0 9 0 0 0 10 0 0 0 11 0 0 0 12 0 0 0 13 0 0 0 14 0 0 0 15 0 0 0 16 1 1 171 17 1 1 171 18 1 1 171 19 1 1 171 ・・・・・・・・・ 236 222 180 222 237 222 180 222 238 222 180 222 239 222 180 222 240 255 255 255 241 255 255 255 242 255 255 255 243 255 255 255 244 255 255 255 245 255 255 255 246 255 255 255 247 255 255 255 248 255 255 255 249 255 255 255 250 255 255 255 251 255 255 255 252 255 255 255 253 255 255 255 254 255 255 255 255 255 255 255
これは左の列のGrayの値の画素値がRGB成分に分解されるときの値を示しています。
この表がLUTの実体です。
最後にLUTの各値をプロットしたグラフが現れます。
 |
| 16_Color.lutのカラー成分のグラフ |
Grayが15以下ではRGBすべて0なので黒で表示され，240以上ではRGBすべて255なので白で表示されます。
Grayが16以上239以下ではRGBそれぞれの割合で混合され様々な色を出力することができます。
CT画像をLUTを使ってカラー画像にする
プログラムimageLUT.pyはDICOM画像を読み込むのためのpydicomモジュールと画像保存のためのcv2モジュールをインポートします。
プログラムの前半はLUTファイルを読み込みrgb各成分をリストr，g，bに読み込みます。
今回はLUTに”BRGBCMYW.lut”を使いました。
DICOM画像を読み込み変数dに代入します。
ウィンドウニング情報を変数wwとwcに，画像の高さと幅をrowとcolに，画像情報をimgに代入します。
DICOM画像のウィンドニング処理した後に，カラースケールを画像中に表示するために255～0の帯を画像に入れます。
for y in range(256): #カラースケース
for x in range(50): #スケールの幅
img[350-y, x] = y #画像中にスケールを入れる。
1行目のfor y in range(256):は，range(256)により変数yに0から255の値が入ります。
yは画素値でもあり，カラースケールの帯の縦方向の表示位置でもあります。
次のfor x in range(50):は，スケールの帯の幅を50にします。
最後に img[350-y, x] = y は，画像img中の350-ｙ,ｘの位置に画素値yを代入します。
カラースケールの入った画像をグレー画像のまま”original.png”のファイル名で保存します。
ここで350は画像の上縁から350画素の位置からカラースケールが始まることを示します。
次にLUTを使って画像imgをRGBの各成分に分けるために，画素値がすべて0のrimg，gimg，bimgを作成します。
rimg = np.zeros((row, col))
次に二重のfor文を使って，x，yの画素の値img[y,x]をリストr，g，bを使って各成分の値に変換して画像rimg，gimg，bimgに代入します。
rimg[y, x] = r[int(img[y, x])]
Opencvのmergeメソッド使ってRGB成分の画像をマージしてカラー画像cimgを作成します。
cimg = cv2.merge((rimg, gimg, bimg))
pyplotを使ってコンソールに表示します。
 |
| pyplotでコンソールに表示された疑似カラー画像 |
OpneCVの処理系はカラー画像をBGRとの順で管理しているため他のPILなどと扱いが異なるため注意が必要です。
BGRのままで画像を保存する次のようになります。
 |
| BGR順の画像で正しく色が表示されていません。 |
OpneCVで画像を保存するとき，その直前でRGBをBGRに変換して保存します。
BGR→RBG変換にはいくつかの方法があります。
ひとつはmergeするときに順番をBGRとしてカラー画像を作成する方法です。
cimg = cv2.merge((bimg, gimg, rimg))
先のmergeとはrimgとbimgが反対になっています。
もう一つの方法はリストの技を使う方法です。
cimg = cimg[:, :, ::-1]
既に存在するRGBG画像をBGR画像に変換するのに便利です。
imageLUT.py
import numpy as np
from matplotlib import pyplot as plt
import pydicom
import cv2
lutfnm = "ImageJインストールディレクトリ\ImageJ\luts\BRGBCMYW.lut"
#lutfnm = "ImageJインストールディレクトリ\ImageJ\luts\\16_Colors.lut"
#lutfnm = "ImageJインストールディレクトリ\ImageJ\luts\Red_Hot.lut"
f = open(lutfnm, "rb")
dt = f.read()
f.close()
print(">1>", len(dt), type(dt))
rgb = np.array([i for i in dt]) # [r0・・・ｒ２５５ g0・・・g255 b0・・・b255]
r, g, b = np.reshape(rgb,(3,256)) # [[r0・・・ｒ２５５] [g0・・・g255] [b0・・・b255]]
x = np.linspace(0,255,256) # 0〜255を256分割する
plt.plot(x, r, color='red')
plt.plot(x, g, color='green')
plt.plot(x, b, color='blue')
plt.show()
# DICOM画像を読み込む
d = pydicom.dcmread('../dcmdir1/Brain01')
wc = d.WindowCenter #ウィンドウセンターの取得
ww = d.WindowWidth #ウィンドウ幅の取得
row= d.Rows #画像の行数の取得
col= d.Columns #画像の列数を取得
img = d.pixel_array #画像を取得
#DICOM画像のウィンドニング
max = wc + ww / 2 #表示最大画素値
min = wc - ww / 2 #表示最小画素値
img = 255*(img-min)/(max-min)
img[img > 255] =255
img[img < 0] = 0
# カラースケールを画像に入れる
for y in range(256): #カラースケース
for x in range(50): #スケールの幅
img[350-y, x] = y #画像中にスケールを入れる。
cv2.imwrite("original.png", img)
#LUTを使ってR,G,Bのチャンネルを作成する。
rimg = np.zeros((row, col)) #Rチャンネルの初期化
gimg = np.zeros((row, col)) #Gチャンネルの初期化
bimg = np.zeros((row, col)) #Bチャンネルの初期化
for y in range(row):
for x in range(col):
#LUTで画像を各チャンネルに変換
rimg[y, x] = r[int(img[y, x])]
gimg[y, x] = g[int(img[y, x])]
bimg[y, x] = b[int(img[y, x])]
cimg = cv2.merge((rimg, gimg, bimg)) # RGB
plt.imshow(cimg)
plt.show()
#OpenCV は画像のチャンネル順を BGR として扱うためBGRの順ででマージする。
#cimg = cv2.merge((bimg, gimg, rimg)) # BGR
cimg = cimg[:, :, ::-1] # BGR->RGB変換
cv2.imwrite("temp.png", cimg)
プログラムを実行してみる
LUTには”BRGBCMYW.lut”を用いています。
実行するとLUTのグラフとLUT処理した画像がコンソールに現れます。
カレントディレクトリにカラースケールを入れたグレー画像original.pngとLUT処理された画像temp.pngが保存されます。
 |
| LUT処理前のグレー画像 左にカラースケール（まだグレーですが）が入っています。 |
 |
| LUT（BRGBCMYW.lut）処理後の画像 左のカラースケールに色がついている |
ImageJはたくさんのLUTを提供しています。是非いろいろ試してみてください。
３つほどカラースケールを下に示しました。
目的に応じて使い分けるとよいでしょう。
おわりに
今回はLUTの処理について，ImageJのLUTファイルを読み込むところから説明をしました。
様々な処理結果を見やすくすにはLUTを使った疑似カラー表示が適している場合が多いと思います。
LUTによる表示は実験結果の効果的な表現に便利です。
