Pyゼミ1.08 PET画像を読む

医療画像の画素値を定量評価に使う

PET画像など画素値を定量評価に使う場合があります。

以前に紹介した「Pyゼミ1.02　DICOMからJPEGなどへの変換」では、CT画像をウインドニングして0から255の値に変換しています。

したがって、本来のCT値を読むことはできません。

今回はPET画像を対象にPETの本来の画素値を取り出す方法を紹介します。

※PET, Positron Emission Tomography（陽電子放出断層撮影）

プログラミング　スタート

プログラムはここから

PET画像はここから

初めに準備です。    

importするライブラリは以下の通りです。

本記事ではpydicomのapply_modality_lutの機能を紹介しますのでimportしておきます。

ファイル名を変数dcmfnmにセットしておきます。

import pydicom

import numpy as np 

import cv2

from pydicom.pixel_data_handlers.util import apply_modality_lut

dcmfnm = "./an_pet.50.dcm"

PET画像を読み16ビットの画素値を取り出す

PETの1画素に割り当てらえているビット数は16ビット（Bits Allocated)です。

そして有効なデータとして15ビット（Bits Stored)を使って格納されています。

16ビットの深さでリスケール変換して取り出すメソッドreadimageWindowing16メソッドを作成します。

def readimageWindowing16(filepath):
    print("Image Widowing16")
    ds = pydicom.dcmread(filepath)
    #print("ds=\n",ds)
    img = ds.pixel_array
    print("\timg shape=",img.shape, type(img))
    print("\tOrg img max=",np.max(img),"min=",np.min(img))
    
    ri  = ds.RescaleIntercept
    rs  = ds.RescaleSlope
    
    img = img*rs + ri   # Convert pixel value to PET value
    img = img.astype(np.uint16) # Unsigned 16bit
    print("\trsIntersept =",ri, "\n\trsSlope     =",rs)
    
    return img

img16  = readimageWindowing16(dcmfnm)

print("\nWindowing 16bit img: max =", np.max(img16), "min=",np.min(img16), "type=",img16.dtype)

cv2.imwrite("win16.png", img16)

このメソッドは、引数filepathにDICOMファイルのpathを与えると、PETの正しい画素値に変換された画像データが返ります。

dcmreadメソッドを使い、DICOMのデータセットをdsに読み込んでいます。

そして画像データpixel_arrayをimgに代入しています。

正しいPETの画素値に変換するためにリスケール切片とリスケール傾斜をそれぞれ変数riとrsに代入します。

PETの画素値への変換は、読み込んだ画像imgにリスケール傾斜rsを掛けてリスケール切片riを加えるだけです。

最後にこのメソッドを呼び出し、PETの正しい画素値に変換された画像をimg16に代入しています。

また、変換後の最大と最小を表示して、PNGファイル（win16.png)に保存しています。

実行してみましょう

実行するとコンソールに以下のように表示されます。

Image Widowing16

        img shape= (128, 128) <class 'numpy.ndarray'>

        Org img max= 32767 min= 0

        rsIntersept = 0 

        rsSlope     = 1.728227

Converted Windowing 16bit img: max = 56628 min= 0 type= uint16

画像サイズは128×128です。

最大画素値は32,767です。これは2進数で見ると15ビットすべて1になります。

つまり15ビット長の最大値です。

最小画素値は0です。

変換後の最大画素値（PETの画素値）は56,628となっています。

これは画素値にリスケール傾斜を掛けているために値が1.728227倍になっています。

画素値を格納するビット数が15ビットの場合、最大値は32,767に決定されます。

しかし、実際には15ビットで扱える最大値を超える場合があります。このとき、リスケール傾斜で調整を行っています。

結果画像を見てみましょう。

PETファントムで収集した再構成画像です。きれいに表示されています。

img16の配列の各要素にはPETの各値が格納されています。

この値を使って定量評価ができますね。

apply_modality_lutを使ってみる

pydicomのpixel_data_handlers.utilライブラリの中にapply_modality_lutメソッドがあります。

これは、読み込んだ画像に対して簡単にモダリティのLUTやリスケール操作を行うことができます。

それではreadimageModalityLutメソッドを作成して、このメソッドを呼び出してます。

変換された画像imgModの最大値と最小値を表示して、画像を保存します。

def readimageModalityLut(filepath):

    print("\nConvert with Modalyty_Lut")

    ds = pydicom.dcmread(filepath)   

    pxarry = ds.pixel_array

 

    img = apply_modality_lut(pxarry, ds)

    img = img.astype(np.uint16)

    

    return img

imgMod = readimageModalityLut(dcmfnm) 

print("\nConverted ModalityLut img: max =", np.max(imgMod), "min=",np.min(imgMod), "type=",imgMod.dtype)    

cv2.imwrite("petModLut.png", imgMod)

実行してみましょう

実行すると次のようにコンソール表示されます。

Convert with Modalyty_Lut

Converted ModalityLut img: max = 56628 min= 0 type= uint16

変換後の最大値は56,628となっています。

先のreadimageWindowing16メソッドの結果と同じです。

しかし、プログラムはリスケール傾斜や切片を読むことなく画素データpxarryをapply_modality_lutメソッドで直接PETの画素値に変換しています。

画像も以下のようにきれいに表示されています。

8ビットでウインドニングしてしまう

16ビットの画素値を取り出さなくてはいけないデータを8ビットでウインドニングして出力するとどうなるでしょう。

まず、8ビットでウインドニングするメソッドを作成します。

def readimageWindowing8(filepath):

    print("Image Widowing 8bit")

    ds = pydicom.dcmread(filepath)

    #print("ds=\n",ds)

    img = ds.pixel_array

    print("\timg shape=",img.shape, type(img))

    print("\tOrg img max=",np.max(img),"min=",np.min(img))

    

    wc  = ds.WindowCenter

    ww  = ds.WindowWidth

    ri  = ds.RescaleIntercept

    rs  = ds.RescaleSlope

    

    img = img*rs + ri   # Convert pixel value to PET value

    print("\twc=",wc,type(wc), "\n\tww=", ww, type(ww))

    print("\trsIntersept =",ri, "\n\trsSlope     =",rs)

    

    maxval = wc + ww//2

    minval = wc - ww//2

    img8 = (img - minval)/(maxval - minval)*255

    img8 = np.clip(img8, 0, 255).astype(np.uint8)

       

    return img8   

img8   = readimageWindowing8(dcmfnm)  

print("\nConverted Windowing  8bit img: max =", np.max(img8), "min=",np.min(img8), "type=",img8.dtype)

cv2.imwrite("win8.png", img8)

readimageWindowing8メソッドの中で画像imgをリスケール傾斜と切片を使って変化するところまでは同じです。

その後、ウインドウセンタwcとウインドウ幅wwを用いて表示する最大値maxvalと最小値minvalを決めています。

そして画像のウインドニング処理を行っています。

実行してみると次のようにコンソールに表示されます。

Image Widowing 8bit

        img shape= (128, 128) <class 'numpy.ndarray'>

        Org img max= 32767 min= 0

        wc= 6351.233 <class 'pydicom.valuerep.DSfloat'> 

        ww= 12702.47 <class 'pydicom.valuerep.DSfloat'>

        rsIntersept = 0 

        rsSlope     = 1.728227

Converted Windowing  8bit img: max = 255 min= 0 type= uint8

ウインドウセンタwcとウインドウ幅がそれぞれ6351.233と12702.47となり、これから計算される最大値maxvalと最小値minvalは12702.233と0.2330となります。最大値は変換後の画素の持つ最大値56,628に遠く及びません。

出力したPNG画像は白くつぶれた画像になってしまいます。 

おわりに

ある先生からPETの定量解析をしたいという要望があり本記事を紹介しました。

以前紹介した8ビット画像で出力する方法ですと最大値は255に固定されてしまいます。

また、ウインドセンタとウインドウ幅の適切な値がDICOMタグに入っていない場合出力した画像自体おかしな表示になります。

ImageJでも表示しましたがもちろん同じ結果でした。

この原因として考えられるのは、最初に表示した画像でウインドウセンタとウインド幅を調整した後、それ以降の画像は調整されず、過去の調整結果が引き継がれたのでしょう。

ファントム実験のデータなのでそういうこともありそう。

最大値が56,628、最小値が0の画像を表示するのであれば、適切なウインドウセンタが28,314でウインドウ幅が56,628となるはずです。

ブログに上げようと思い、掲載まで約1年かかってしまいました。

ほんとうに年月の経つのは早いですね。