Pyゼミ1.02 DICOMからJPEGなどへの変換

DICOM画像をJPEGやPNGなどへの一般的な画像に変換しよう。

keyword: DICOM, OpenCV，JPEG, PNG，ImageJ

DICOM画像は専用のうDICOMビューワがないと表示することができません。
さいわいImageJのような高機能なフリーのビューワソフトも存在しますが，今回は一般的なJPEG画像など扱いやすい形に変換するプログラムを作成します。

プログラムはDICOM画像をウィンドニング処理（表示に最適な濃度変換処理）を加えてJPEGに変換します。
JPEG画像は非可逆な圧縮であり，圧縮時に画質を指定して保存することができます。
画質とファイルサイズの関係，画質低下に伴うエラーの状態を確認してみます。
さらに可逆圧縮（完全に元画像に復元できる圧縮方法）であるPNGファイルについても検討してみます。

OpenCVの画像ファイルの入出力関数を使う

OpneCV（Open Source Computer Vision Library）は画像処理をサポートするライブラリで，様々な関数が用意されて非常に便利なライブラリです。
このライブラリの中に画像の入出力もあるので，ここでは画像の書き込みにOpenCVのメソッドを使うことにします。
OpenCVの他のメソッドについては別途Pyゼミで紹介します。

OpenCVのインストールとバージョンの確認

OpenCVをインストールします。

$ pip␣install␣opencv-python⏎   

pythonを起動して，opneCVをインポートしてバージョンを確認してみます。

>>> import cv2
>>> cv2.__version__
'3.1.0'

OpenCVの3.1がインストールされているのが確認できました。

DICOM画像をJPEGに画像に変換して保存する

DICOM画像（この画像は実験にのみ使用してください）を読み込んでJPEGファイルに保存するプログラムを入力して実行してみましょう。

dicom2jpeg.py

import pydicom
import cv2
from matplotlib import pyplot as plt
    
def readDicom2jpeg(dcmfnm, savefnm):       #DICOM画像を読込，ウィンドニングしてJPEGに保存
    ds  = pydicom.read_file(dcmfnm)        #DICOM画像を読み込む
    wc  = ds.WindowCenter                　#ウィンドウセンター値を代入
    ww  = ds.WindowWidth                 　#ウィンドウ幅を代入       
    img = ds.pixel_array                 　#画素値を代入        
        
    #表示画素値の最大と最小を計算する
    max = wc + ww/2                     
    min = wc - ww/2      
    print("wc=",wc,"ww=",ww,"→ max=",max," min=",min)  
   
    #ウインドニング処理
    img = 255*(img - min)/(max - min)   　 #最大と最小画素値を0から255に変換        
    img[img > 255] = 255                  #255より大きい画素値は255に変換
    img[img < 0] = 0     　　　　　　　　 #JPEG画像として保存
    cv2.imwrite(savefnm, img, [cv2.IMWRITE_JPEG_QUALITY, 100])

if __name__ == '__main__':    
    readDicom2jpeg("../dcmdir1/Brain01","brain.jpg")   

このプログラムを実行するとカレントディレクトリにbrain.jpgが作成されます。

$ python␣dicom2jpeg.py⏎

生成されたbrain.jpgファイル

プログラムの解説

このプログラムは，プログラム内で使用するライブラリのインポート部とDICOM画像を読み込んでJPEG画像に保存するreadDicom2jpeg関数とプログラムの実行時に最初に呼び出されるmain関数から構成されています。

インポート部はDICOM画像を読み込むためのpydicomとJPEG画像を保存するためのOpenCVのインポートです。
後の実験のためにpyplotもインポートします。

readDicom2jpeg関数は，読み込むDICOMファイル名dcmfnmとJPEG画像を保存するファイル名savefnmを引数にとります。

はじめにDICOM画像を読み込み変数dsに代入します。

このdsからウィンドウセンタ，ウィンドウ幅と画像データをそれぞれ変数wc, ww, imgに代入しています。
imgにはCT画像の画素値が16ビットで格納されています。

次にウィンドニング処理を行います。
ウィンドニング処理はCT値の上限（max）を255に下限（min）を0になるように変換します。

maxとminはウィンドウセンターwcとウィンドウ幅wwから計算します。

    　　　max = wc + ww/2                   
    　　　min = wc - ww/2 　　　　

minとmaxの間のCT値を0から255の8bitの値に次式で変換します。

　　　　img = 255*(img - min)/(max - min)

従来（Javaなど）であれば，二重のfor文で各画素を変換するとろ，pythonは上のようにたった1行で処理してくれます。
Python素晴らしい。

変換後のimgの画素値が0未満の値（負の値）は0に，255より大きい値は255に変換します。

　　　　img[img > 255] = 255            

　　　　img[img < 0] =  0
この処理も従来（Javaなど）であれば，2重のfor文で各画素を変換していきますが，pythonでは上の2行で変換を行います。
しかも，for文で書くよりも高速に処理します。
Python素晴らしい。 

OpneCvのimwriteメソッドを使って，ウインドニング処理した画像imgをJPEG画像としてファイル名savefnmで保存します。

JPEG画像は劣化する

JPEG画像を保存するときに画質をcv2.IMWRITE_JPEG_QUALITYで指定することができます。
この値（Quality Factor, 以下QF）が100のとき，最高画質で保存します。
QFを小さくすると画質は低下しますが，圧縮したファイルサイズは小さくなります。

 それでは以下のプログラムをreadDicom2jpeg関数の最後に（JPEG保存のimwriteメソッドの下に）追加して実験してみましょう。
JPEGとは別に可逆圧縮（画質の劣化がない）のPNGでも保存してみましょう（ファイル名の拡張子を.pngにするだけ）。

    cv2.imwrite("QF70" + savefnm, img, [cv2.IMWRITE_JPEG_QUALITY, 70])
    cv2.imwrite("QF50" + savefnm, img, [cv2.IMWRITE_JPEG_QUALITY, 50])
    cv2.imwrite("QF30" + savefnm, img, [cv2.IMWRITE_JPEG_QUALITY, 30])
    cv2.imwrite("Brain.png", img)

作成されたJPEG画像のファイルサイズを見るとQFが小さくなると圧縮率も高くなり，ファイルサイズは小さくなっています。

　　　　　brain.jpg　　　   100.7kB
　　　　　QF70brain.jpg　　26.7kB
　　　　　QF50brain.jpg　　20.7kB
　　　　　QF30brain.jpg　　15.7kB

画像を拡大表示してみると

QFが70程度では画像の劣化は目立ちません（下図中）。
QFが30ではブロックノイズと呼ばれる8ｘ8画素の大きさのブロックが目立つようになります（下図右）。
JPEGは8ｘ8画素単位でDCT（discrete cosine transform，離散コサイン変換）を行い，高周波数成分を除くことで圧縮を行っています。
したがって，QFを小さくするとブッロク間の境界が目立つことになります（ブロックノイズ）。

オリジナルの画像との誤差を見てみよう

元のDICOMを読み込んだ画像データimgとJPEG圧縮した画像データimg2との差分により２つの画像の間の誤差を見てみましょう。

readDicom2jpeg関数の最後に下のプログラムを追加して実行してみよう。
imreadメソッドの引数の 0 は画像読み込みのフラグで，0はgray画像（1チャンネルの画像）として読み込むことを示します。
この引数はディフォルトでは3で，カラー画像として読み込みます。

    max = 5
    img2 = cv2.imread("QF30" + savefnm, 0)
    dif = abs(img - img2)
    plt.imshow(dif, cmap = 'gray', vmax = max, vmin = 0)
    plt.show()
    
    img2 = cv2.imread(savefnm, 0)
    dif = abs(img - img2)
    plt.imshow(dif, cmap = 'gray', vmax = max, vmin = 0)
    plt.show()
    
    img2 = cv2.imread("Brain.png", 0)
    dif = abs(img - img2)
    plt.imshow(dif, cmap = 'gray', vmax = max, vmin = 0)
    plt.show()

実行すると元のDICOM画像とQFを変えて保存した画像との差分画像（誤差画像）が表示されます。

QFが30との差分画像は明らかにブロックノイズが認められ，ノイズが多いのがわかります（下図左）。

QFを100にした最高品質のJPEG画像でも僅かにノイズがあるのがわかります（下図中）。

PNG画像は可逆圧縮であるため基本的に元画像との間に際はないはずだが，JPEGのQF100と同程度のノイズが認められます（下図右）。
これはOpenCVの内部的なアルゴリズムが原因かもしれません。

参考　ImageJによる差分画像

参考までにImageJで同様の実験結果を示します。
ImageJでDICOM画像を読み込み，8bitに変換した後にQFを変えたJPEG画像との差分（subtraction）を行いました。

QFが30のとき，差分画像は濃度変化の大きおエッジ部分に誤差が認められます。
QFが100のとき，わずかに誤差があるのが認められます。
PNG画像の差分画像には全く誤差が認められません。つまりPNGは完全に元に戻る画像圧縮フォーマットということです。

おわりに

今回はDICOM画像をウィンドニング処理を行い一般的な画像フォーマットに変換する方法を学びました。
この処理はこの後の様々なプログラムで再利用します。