マイブログ リスト

医療言語処理講座

2022年10月19日水曜日

Pyゼミ0.06 備忘録 Ubuntu20.04にCUDA11.4、TensorFlowをインストールする 

はじめに

Ubuntuを20.04にして、NVIDIAのドライバーをインストールするとCUDA12.0がインストールされます(2023年2月20日)。

問題なくインストールされnvidia-smiを正常に起動されます。

しかし、TensorflowがGPUを認識しません。

そこで、今回CUDAを11.4にバージョンダウンして利用する方法を投稿します。


1)Ubuntu20.04インストール後 ソフトウェアアップデート

$ sudo apt update
$ sudo apt upgrade

$ sudo reboot


2)CUDA&Driverインストール

下記URLにてCUDA Toolkit Archiveのページへ
https://developer.nvidia.com/cuda-toolkit-archive

バージョン11.4の下記を選択する
CUDA Toolkit Documentation v11.4.4 

Select Target Platformは以下を選択
  • OS:Linux
  • Architecture:x86_64
  • Distribution:Ubuntu
  • Version:20.04
  • Installer Type:dev(local)

注意:cudaのインストール時に明示的に11.4をcuda-11-4と指定すること


$ wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2004/x86_64/cuda-ubuntu2004.pin
$ sudo mv cuda-ubuntu2004.pin /etc/apt/preferences.d/cuda-repository-pin-600
$ wget https://developer.download.nvidia.com/compute/cuda/11.4.4/local_installers/cuda-repo-ubuntu2004-11-4-local_11.4.4-470.82.01-1_amd64.deb
$ sudo dpkg -i cuda-repo-ubuntu2004-11-4-local_11.4.4-470.82.01-1_amd64.deb
$ sudo apt-key add /var/cuda-repo-ubuntu2004-11-4-local/7fa2af80.pub
$ sudo apt-get update
$ sudo apt-get -y install cuda-11-4


$ sudo reboot



確認1

$ nvidia-smi

Sun Mar 26 13:51:22 2023       
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 470.82.01    Driver Version: 470.82.01    CUDA Version: 11.4     |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|                               |                      |               MIG M. |
|===============================+======================+======================|
|   0  NVIDIA GeForce ...  On   | 00000000:01:00.0  On |                  N/A |
| 38%   33C    P0    28W / 120W |    188MiB /  6069MiB |      0%      Default |
|                               |                      |                  N/A |
+-------------------------------+----------------------+----------------------+
                                                                               
+-----------------------------------------------------------------------------+
| Processes:                                                                  |
|  GPU   GI   CI        PID   Type   Process name                  GPU Memory |
|        ID   ID                                                   Usage      |
|=============================================================================|
|    0   N/A  N/A      1084      G   /usr/lib/xorg/Xorg                 35MiB |
|    0   N/A  N/A      1569      G   /usr/lib/xorg/Xorg                 58MiB |
|    0   N/A  N/A      1698      G   /usr/bin/gnome-shell               84MiB |
+-----------------------------------------------------------------------------+


確認2

/usr/localにcudaの関連ディレクトリがあるか確認する

$ sudo find /usr/local | grep cuda

/usr/local/cuda-11.4
/usr/local/cuda-11.4/compute-sanitizer
   ・・・

   



3)CUDAパス設定

CUDA のライブラリにパスを通します

$ gedit .bashrc

次の3行を .bashrc ファイルの最後に追加

#CUDA
export PATH="/usr/local/cuda/bin:$PATH"
export LD_LIBRARY_PATH="/usr/local/cuda/lib64:$LD_LIBRARY_PATH"


確認

$ source .bashrc
$ nvcc -V
nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2021 NVIDIA Corporation
Built on Mon_Oct_11_21:27:02_PDT_2021
Cuda compilation tools, release 11.4, V11.4.152
Build cuda_11.4.r11.4/compiler.30521435_0



4)cuDNNのインストール

下記のURLからNVIDIA Developer Program Membership Requiredへ
https://developer.nvidia.com/rdp/cudnn-download

ログイン(またはJoin nowにて新規登録)

cuDNN Downloadのダウンロードページへ

今回はcuDNN v8.6をインストールするため、最新のボタンの下の下記を選択
            Archived cuDNN Releases 

リストから下記を選択しファイルをダウンロード
            Download cuDNN v8.6.0 (October 3rd, 2022), for CUDA 11.x
            Local Installer for Ubuntu20.04 x86_64 (Deb)

ダウンロードしたファイル名:cudnn-local-repo-ubuntu2004-8.6.0.163_1.0-1_amd64.deb

ファイルをhomeディレクトリに移動し、以下を実行

$ sudo dpkg -i cudnn-local-repo-ubuntu2004-8.6.0.163_1.0-1_amd64.deb

$ sudo cp /var/cudnn-local-repo-ubuntu2004-8.6.0.163/cudnn-local-B0FE0A41-keyring.gpg /usr/share/keyrings/
$ sudo apt update
$ sudo apt -y install libcudnn8 libcudnn8-dev libcudnn8-samples


確認


$ dpkg -l | grep cudnn

ii  cudnn-local-repo-ubuntu2004-8.6.0.163      1.0-1                               amd64        cudnn-local repository configuration files
ii  libcudnn8                                  8.6.0.163-1+cuda11.8                amd64        cuDNN runtime libraries
ii  libcudnn8-dev                              8.6.0.163-1+cuda11.8                amd64        cuDNN development libraries and headers
ii  libcudnn8-samples                          8.6.0.163-1+cuda11.8                amd64        cuDNN samples



5)Anacondaのインストール

Python 3.9をインストールするため、昨年のファイルを使用

https://www.anaconda.com/distribution/

ダウントード:Anaconda3-2023.03-Linux-x86_64.sh


$ bash Anaconda3-2023.03-Linux-x86_64.sh


確認

$ source .bashrc

$  conda -V
conda 4.12.0


$ python
[GCC 7.5.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> 
Python 3.9.12 (main, Apr  5 2022, 06:56:58) 


6)TnesorFlow と Kerasのインストール

Tensorflowの2.8.2をインストール。

次を実行することでGPU対応のTensorflowとKerasがインストールされる

$ pip install tensorflow==2.8.2


GPU確認

python起動後

>>> from tensorflow.python.client import device_lib
>>> device_lib.list_local_devices()
2023-03-26 14:06:15.580531: I tensorflow/core/platform/cpu_feature_guard.cc:151] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
2023-03-26 14:06:15.612622: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:15.706883: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:15.708182: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:16.083901: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:16.084200: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:16.084453: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:936] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
2023-03-26 14:06:16.084689: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1525] Created device /device:GPU:0 with 5239 MB memory:  -> device: 0, name: NVIDIA GeForce GTX 1060 6GB, pci bus id: 0000:01:00.0, compute capability: 6.1
[name: "/device:CPU:0"
device_type: "CPU"
memory_limit: 268435456
locality {
}
incarnation: 3029067233870127766
xla_global_id: -1
, name: "/device:GPU:0"
device_type: "GPU"
memory_limit: 5493751808
locality {
  bus_id: 1
  links {
  }
}
incarnation: 7369121378162548002
physical_device_desc: "device: 0, name: NVIDIA GeForce GTX 1060 6GB, pci bus id: 0000:01:00.0, compute capability: 6.1"
xla_global_id: 416903419
]
>>> 



7)OpenCVのインストール

$ pip install opencv-python

Collecting opencv-python
  Downloading opencv_python-4.6.0.66-cp36-abi3-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (60.9 MB)
     |████████████████████████████████| 60.9 MB 10.5 MB/s 
Requirement already satisfied: numpy>=1.14.5 in ./anaconda3/lib/python3.9/site-packages (from opencv-python) (1.21.5)
Installing collected packages: opencv-python
Successfully installed opencv-python-4.6.0.66


確認

>>> import cv2
>>> cv2.__version__
'4.7.0'



8)Pydicomインストール

$ pip install  pydicom

Collecting pydicom
  Downloading pydicom-2.3.0-py3-none-any.whl (2.0 MB)
     |████████████████████████████████| 2.0 MB 3.1 MB/s 
Installing collected packages: pydicom
Successfully installed pydicom-2.3.1



9)DCMTKのインストール

$ sudo apt install dcmtk


確認

$ dcmdjpeg --version

(以下の内容が表示されればOK)

$dcmtk: dcmdjpeg v3.6.4 2018-11-29 $
dcmdjpeg: Decode JPEG-compressed DICOM file
Host type: Debian
Character encoding: UTF-8
External libraries used:
- ZLIB, Version 1.2.11
- IJG, Version 6b  27-Mar-1998 (modified)





10)ImageJのインストール

ImageJのページからLinux版をダウンロード

https://imagej.nih.gov/ij/

ダウンロード:ij153-linux64-java8.zip

Homeで展開、起動してアイコンをお気に入りに登録する


GUI画面が小さい場合

Edit/Options/AppearanceのGUI scaleを1.5にしてみる。




11)その他の設定

○メニューに空のドキュメントを追加する

# フォルダーが日本語のままの場合

$ touch ~/テンプレート/空のドキュメント


○エディタ

・行番号の表示

・Tab 4文字空白で設定

・フォントサイズ 16


-----------------------------------

 おまけ)Dropboxのインストール

$ cd ~ && wget -O - "https://www.dropbox.com/download?plat=lnx.x86_64" | tar xzf -


起動

$ ~/.dropbox-dist/dropboxd



2022年5月6日金曜日

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年かかってしまいました。
ほんとうに年月の経つのは早いですね。


2022年5月2日月曜日

Pyゼミ1.0 番外編 DICOM RDSRを読む

 Pydicomを使ってRDSRを読んで必要な情報を取り出してみる

RDSRはRadiation Dose Structured Reportの略で、CTや血管造影などX線を用いた検査でどのくらX線量(Radiation dose)を照射したのか、構造化された書式(Structured report)にて出力したものです。

RDSRはDICOM規格で規定されています。詳細は規格書を読みましょう。

今回は時間がなく規格書を理解せず、帰納的にSRの解析に取り組みました。
解析結果が正しいかどうか、医療現場で使用しているソフトウェアの結果で検証しています。
検証した範囲では問題なさそうですが、利用は自己責任でお願いします。

プログラムはこちらから
RDSRファイルはこちらから

(この記事は2022年4月24日に日本放射線技術学会 北海道支部 第78回春季学術大会 医療情報セミナーで発表したものです)


PydicomはDICOM画像だけでなくSRも読むことができます。RDSRファイルを読み込んで表示するだけであれば数行で書くことができます。


import pydicom
fnm = "./CT_rdsr.SR"
ds = pydicom.dcmread(fnm)   
print(ds)


このたった4行でRDSRのDICOMファイル"rdsr.sr"の内容をpydicomを用いて変数dsに読み込み、それをプリントすることでRDSRの内容を確認することができる。


DICOM RDSRの内容を見てみる


DICOMのファイルはデータセット呼ばれるデータエレメントの集合でできています。
データエレメントは、タグと呼ばれる16進数のグループ番号とエレメント番号から始まり、キーワード、値表現(VR)、値の長さ(VL)と値(VF, Value field)で構成されています。

print(ds)で出力された内容を見てみしょう。
表示された各行はデータエレメントを表していて次のように表示されています。

                (グループ番号, エレメント番号) キーワード VR:値

(0002,0002)と(0008,0016)のSOP Class UIDが”X-Ray Radiation Dose SR Storage"と記述されています。
また、(0008,0060)のModalityは”SR”とあり、このファイルがStructured reportであることがわかります。

Dataset.file_meta -------------------------------
(0002, 0000) File Meta Information Group Length  UL: 196
(0002, 0001) File Meta Information Version       OB: b'\x00\x01'
(0002, 0002) Media Storage SOP Class UID         UI: X-Ray Radiation Dose SR Storage
(0002, 0003) Media Storage SOP Instance UID      UI: 1.3.12.2.1107.5.1.4.65939.30000022022322414137700000249
(0002, 0010) Transfer Syntax UID                 UI: Explicit VR Little Endian
(0002, 0012) Implementation Class UID            UI: 1.3.12.2.1107.5.4.7
(0002, 0013) Implementation Version Name         SH: 'SIEMENSXLEOVB21C'
-------------------------------------------------
(0008, 0005) Specific Character Set              CS: 'ISO_IR 100'
(0008, 0016) SOP Class UID                       UI: X-Ray Radiation Dose SR Storage
(0008, 0018) SOP Instance UID                    UI: 1.3.12.2.1107.5.1.4.65939.30000022022322414137700000249
(0008, 0020) Study Date                          DA: '20220224'
(0008, 0021) Series Date                         DA: '20220224'
(0008, 0023) Content Date                        DA: '20220224'
(0008, 0030) Study Time                          TM: '074907.240000'
(0008, 0031) Series Time                         TM: '075458.275000'
(0008, 0033) Content Time                        TM: '075458.315000'
(0008, 0050) Accession Number                    SH: '6010016637053'
(0008, 0060) Modality                            CS: 'SR'



上記以外の情報は通常のDICOM画像とそう変わりありません。

しかし、もう少し下のほうを見てみると、(0040,0275)のRequest Attributes Sequence以下で字下げ(インデント)が見られ、
さらに(0040,0008) Scheduled Protocol Code Sequence以下でインデントが発生しています。
これはVRが"SQ"に設定されているためです。
”SQ”は「ある意味のある集合のシーケンス(連続した項目)を含み、このシーケンスは入れ子状態をとることができる」ものです。
詳しくはDICOM Part5(http://www.jira-net.or.jp/dicom/dicom_data_02_01.html)を参照してください。

(0029, 1008) [syngo Report Type]                 CS: 'CTDOSEREPORT'
(0029, 1009) [syngo Report]                      LO: '1.0'
(0029, 1015) [SR Variant]                        US: 6
(0029, 1017) [SC SOP Instance UID]               UI: 1.3.12.2.1107.5.1.4.65939.30000022022322414137700000247
(0040, 0275)  Request Attributes Sequence  1 item(s) ---- 
   (0040, 0008)  Scheduled Protocol Code Sequence  1 item(s) ---- 
      (0008, 0100) Code Value                          SH: '3003100550020000'
      (0008, 0102) Coding Scheme Designator            SH: '99'
      (0008, 0103) Coding Scheme Version               SH: '6010016637053'
      (0008, 0104) Code Meaning                        LO: 'AG?F?~???ETAE?F?S?g'
      ---------
   (0040, 0009) Scheduled Procedure Step ID         SH: '6010016637053'
   (0040, 1001) Requested Procedure ID              SH: '6010016637053'
   ---------


この"SQ"をさらに追いかけてみます。
すると、(0040, a730) Content Sequence が20項目あると表示されています。
つまり、これ以下に何らかの意味ある項目が20項目あることを示しています。

(0040, a050) Continuity Of Content               CS: 'SEPARATE'
(0040, a372)  Performed Procedure Code Sequence  0 item(s) ---- 
(0040, a491) Completion Flag                     CS: 'COMPLETE'
(0040, a493) Verification Flag                   CS: 'UNVERIFIED'
(0040, a504)  Content Template Sequence  1 item(s) ---- 
   (0008, 0105) Mapping Resource                    CS: 'DCMR'
   (0040, db00) Template Identifier                 CS: '10011'
   ---------
(0040, a730)  Content Sequence  20 item(s) ---- 
   (0040, a010) Relationship Type                   CS: 'HAS CONCEPT MOD'
   (0040, a040) Value Type                          CS: 'CODE'
   (0040, a043)  Concept Name Code Sequence  1 item(s) ---- 
      (0008, 0100) Code Value                          SH: '121058'
      (0008, 0102) Coding Scheme Designator            SH: 'DCM'
      (0008, 0104) Code Meaning                        LO: 'Procedure reported'
      ---------
   (0040, a168)  Concept Code Sequence  1 item(s) ---- 
      (0008, 0100) Code Value                          SH: 'P5-08000'
      (0008, 0102) Coding Scheme Designator            SH: 'SRT'
      (0008, 0104) Code Meaning                        LO: 'Computed Tomography X-Ray'
      ---------
   (0040, a730)  Content Sequence  1 item(s) ---- 
      (0040, a010) Relationship Type                   CS: 'HAS CONCEPT MOD'
      (0040, a040) Value Type                          CS: 'CODE'
      (0040, a043)  Concept Name Code Sequence  1 item(s) ---- 
         (0008, 0100) Code Value                          SH: 'G-C0E8'
         (0008, 0102) Coding Scheme Designator            SH: 'SRT'
         (0008, 0104) Code Meaning                        LO: 'Has Intent'
         ---------

さらに、下を見てみると次のように表示されています。
(0008, 0104)  Code Meaningには'CT Accumulated Dose Data'とあり、
その下に
(0008, 0104)  Code Meaningには'Total Number of Irradiation Events'とあり、
またその下に
(0008, 0104)  Code Meaningには'events'とあり、
次に
(0040, a30a) Numeric Valueには"7.0"とあります。

さらに、その下には
(0008, 0104) Code Meaning には'CT Dose Length Product Total'とあり、
その下に
(0008, 0104) Code Meaningには'mGy.cm'と単位があり、
次に
(0040, a30a) Numeric Valueには"377.94"と数値があります。

   (0040, a010) Relationship Type                   CS: 'CONTAINS'
   (0040, a040) Value Type                          CS: 'CONTAINER'
   (0040, a043)  Concept Name Code Sequence  1 item(s) ---- 
      (0008, 0100) Code Value                          SH: '113811'
      (0008, 0102) Coding Scheme Designator            SH: 'DCM'
      (0008, 0104) Code Meaning                        LO: 'CT Accumulated Dose Data'
      ---------
   (0040, a050) Continuity Of Content               CS: 'SEPARATE'
   (0040, a730)  Content Sequence  2 item(s) ---- 
      (0040, a010) Relationship Type                   CS: 'CONTAINS'
      (0040, a040) Value Type                          CS: 'NUM'
      (0040, a043)  Concept Name Code Sequence  1 item(s) ---- 
         (0008, 0100) Code Value                          SH: '113812'
         (0008, 0102) Coding Scheme Designator            SH: 'DCM'
         (0008, 0104) Code Meaning                        LO: 'Total Number of Irradiation Events'
         ---------
      (0040, a300)  Measured Value Sequence  1 item(s) ---- 
         (0040, 08ea)  Measurement Units Code Sequence  1 item(s) ---- 
            (0008, 0100) Code Value                          SH: '{events}'
            (0008, 0102) Coding Scheme Designator            SH: 'UCUM'
            (0008, 0103) Coding Scheme Version               SH: '1.4'
            (0008, 0104) Code Meaning                        LO: 'events'
            ---------
         (0040, a30a) Numeric Value                       DS: "7.0"
         ---------
      ---------
      (0040, a010) Relationship Type                   CS: 'CONTAINS'
      (0040, a040) Value Type                          CS: 'NUM'
      (0040, a043)  Concept Name Code Sequence  1 item(s) ---- 
         (0008, 0100) Code Value                          SH: '113813'
         (0008, 0102) Coding Scheme Designator            SH: 'DCM'
         (0008, 0104) Code Meaning                        LO: 'CT Dose Length Product Total'
         ---------
      (0040, a300)  Measured Value Sequence  1 item(s) ---- 
         (0040, 08ea)  Measurement Units Code Sequence  1 item(s) ---- 
            (0008, 0100) Code Value                          SH: 'mGy.cm'
            (0008, 0102) Coding Scheme Designator            SH: 'UCUM'
            (0008, 0103) Coding Scheme Version               SH: '1.4'
            (0008, 0104) Code Meaning                        LO: 'mGy.cm'
            ---------
         (0040, a30a) Numeric Value                       DS: "377.94"
         ---------
      ---------
   ---------

このように、Code Meaningの値やNumeric Valueの値を抽出することで必要な情報を得ることができるのではないかと考えました。

少し解析のフローをまとめると

  1. データエレメントのKeywordが“Content Sequence”が見つかったら解析を始める
  2. データエレメントのVRが“SQ“のとき、KeywordとValueを取得する
  3. データエレメントのkeywordが“Code Meaning”のとき、Valueの値を出力する
  4. 取得したValue(データエレメント)について2.の解析を繰り返す。
つまり、VRに"SQ"が見つかったらその下の階の項目を取り出し、その中のVRに"SQ"が見つかったら、さらにその下の。。。と繰り返すことで各階のCode Meaningの情報を取り出すこと考えました。

しかし、階層が3階層より下が上手く取り出すことができませんでした。
これはPydicomの仕様(バグ?)かわかりません。
そこで、タグ、キーワード、VRと値を取り出すメソッド(getInfoDCM)を作成して対応しました。



RDSRを読むPythonプログラム


RDSRを読み込むための準備を最初に書きます。
インポートするのPydicomだけです。
RDSRのファイル名を変数fnmに代入します。
rstxt, dstxtとsmryはSRの処理結果を格納する変数です。
nLvは"SQ"による階層の深さ、インデントの深さを代入する変数です。
prevLは前の"SQ"の行(データエレメント)の深さを格納する変数です。
inkeywordsには、抽出すべき値のキーワードを格納しています。


import pydicom

fnm = "./CT_rdsr.SR"    # RDSR DICOM file name
rstxt = ""      # For result text using getInfoDCM method
dstxt = ""      # For DICOM dataset
smry  = ""      # For target text
nLv   = 0       # Depth Level
prevL = 0       # For depth level of previous SQ line

# For extraction keyword
inkeywords = "CodeMeaning,TextValue,NumericValue,UID,DateTime"  


データエレメントの深さ求めるメソッドとデータエレメントの各要素求めるメソッドについて


1.データエレメントの深さを求めるgetDepthメソッド

 インデントは各行(データエレメント)の前の空白の数を数えて返します。
 引数のdataはひとつのデータエレメントです。

def getDepth(data):
    cnt = 0         
    for i in range(len(data)):
        if data[i] == ' ':
            cnt += 1
        else:
            break
    return cnt


2.データエレメントの各要素求めるgetInfoDCMメソッド
 このメソッドはデータエレメント中のタグ、キーワード、VRと値を抽出して返します。
 引数のdataには (grp,elm) Key word VR: 'value'  のようにひとつのデータエレメントが格納されています。

Pydicomでは"SQ"のシーケンスの終わりに'---------'が表示されます。
メソッドの最初にこの処理を記述します。
与えれたdataの内容が'---------'であった場合、
        tagとVRには何も値を入れず、
        キーワードには、”NoKeyword”と
        値には、'---------'をそのまま入れて返します。


def getInfoDCM(data):
    # ex) data is "(grp,elm) Key word VR: 'value' "
    data = data.strip(" ")  # Remove blanks on both sides of string
    if data == '---------':
        tag = ["",""]
        kwd = "NoKeyword"
        vr  = ""
        val = '---------'
        return tag, kwd, vr, val
        
    # Tag
    tag = data.split(")")[0]    # "(grp, elm"
    tag = tag.replace("(", "").replace(" ","")  # "grp, elm"
    tag = tag.split(",")    # [group, elemnt]
    
    # vr
    vr = ""
    if "item(s)" in data:   # ex) (0040, a043)  Concept Name Code Sequence  1 item(s) ----
        vr = "SQ"
    else:
        vr = data.split(":")[0][-2:]    # ex) "(grp,elm) Key word VR"->"VR"
        
    # value
    val = ""
    if ":" in data:
        val = data.split(":")[-1]       # [・・・ VR:value] -> value
    elif "item(s)" in data:
        val = ' '.join(data.split(" ")[-3:-1])  # -> n item(s) 
    val = val.replace("\'","").replace("\"","") # remove " and '
    
    # keyword
    kwd  = ""
    stat = data.find(")")+1 # Start position in keyword
    end  = data.find(vr)    # End position in keyword
    kwd  = data[stat:end]   # Extract keyword (remove (g,e))
    if "item(s)" in kwd:    # Case of VR="SQ"
        kwd = kwd.split("item(s)")[0]   # [kwd n,----] -> (g,e) kwd n 
        kwd = kwd.strip(" ").split(" ")[0:-1] # Extract the string before "n items"
        kwd = ''.join(kwd).strip(" ")  # Joint keyword and remove blanks      
    else:
        kwd = kwd.replace(" ","")   # Remove blank in keyword
    
    return tag, kwd, vr, val


はじめにTagの処理を見ます。
        (grp,elm) Key word VR: 'value'
グループ番号とエレメント番号は()の中にありますのでこれを頼りに2つの番号を抽出し、リストtagに代入します。

次に、VRです。
        (grp,elm) Key word VR: 'value'
        (grp,elm) Key word 99 item(s) ----
変数dataの中に"item(s)"があればVRは”SQ”です。
"SQ"以外の場合(else)、コロン(:)で変数dataの値を2つ分割し、前半[0]の最後の2文字[-2:]を抽出し変数vrに代入します。

次に、値です。
        (grp,elm) Key word VR: 'value'
        (grp,elm) Key word 99 item(s) ----
変数data中にコロン(:)があれば2つに分割し後半[1]を値として変数valに代入します。
コロンがなく、”item(s)”があれば、変数dataを空白で分割し”99 items(s)"のみ抽出しvalに代入します。
最後にvalの値の前後のシングル/ダブルクォーテーション( ’ / " )を削除します。

最後にキーワードです。
        (grp,elm) Key word VR: 'value'
        (grp,elm) Key word 99 item(s) ----
データエレメントの”)”の後ろの文字(stat番目)からVRの前(end)の文字列を抽出data[stat:end]して変数kwdに代入します。
しかし、この処理では、VRが"SQ"のときキーワードに"n item(s)"が存在します。
これを除く処理を行いました。
また、抽出したキーワード中の空白は除くようにしました。つまり”Study Date”は”StudyDate”になります。
これはPydicomのキーワードの表現に揃えるためです。



メインのプログラムを読む


はじめにpydicomを使ってRDSRのファイルを読み込みデータセットを変数dsに代入します。
SOPClassUIDなど必要な基本情報を表示します。
smryは必要な抽出項目の要約を代入する変数です。

変数dsは'pydicom.dataset.FileDataset'の型を取りますが、これを文字列型に変換して変数txdsに代入します。
txdsを改行(\n)で分割して行単位に分割してリストlndsに代入します。
読み込んだ情報などを表示していったん処理が停止(pause)します。
Enterキーで処理が進みます。

for文でlndsから1行(1データエレメント)ずつ取り出して処理します。

valueに'ContentSequence'が見つかったら、smryflagをTrueにして必要情報を抽出するように処理します。

変数rstxtには行番号、深さのレベル数とともに抽出した結果を追加代入(+=)しています。
一方、dstxtにはデータセットに行番号とデータエレメントを追加代入しています。
この2つの出力ファイルを比較することにより、間違いなく抽出できているか確認できます。
 
変数smryには注目するキーワードのデータエレメントの値を追加代入しています。

# ===== MAIN =====        
ds = pydicom.dcmread(fnm)   # Read DICOM SR

#print("ds=\n",ds)
print("SOP Class UID:", ds.SOPClassUID)
print("Modality     :", ds.Modality)
print("Manufacturer :", ds.Manufacturer)
print("Study Date   :", ds.StudyDate)
print("Pt Name      :", ds.PatientName)
print("Pt ID        :", ds.PatientID)
print("Modality     :", ds.Modality)
print("Content Sequence:", len(ds.ContentSequence))

input("pause. Type ENTER key")

smry  = "Study Date," + str(ds.StudyDate) + "\n"
smry += "Study Time," + str(ds.StudyTime) + "\n"

txds = str(ds)      # Convert pydicom.dataset to String

print(txds)
print("\nData type :",type(ds),"\n\t\t-->",type(txds))
lnds = txds.split("\n")     # Split text to lines
print("len(ds)=",len(lnds)) 
input("pause. Type ENTER key")

for i, ln in enumerate(lnds):
    tag, keyword, vr, value = getInfoDCM(ln)
    nLv = getDepth(ln)//3

    if value in 'ContentSequence':  # Start anlysis from this line
        smryflag = True

    sIndent = " "*(nLv*3)           # Indent
    rstxt += str(i)+",L:"+str(nLv)+sIndent+"("+str(tag)+"),"+keyword+","+vr+","+value+"\n"    
    dstxt += str(i)+","+ln + "\n"
    
    print("value = >",value,"<")
        
    # CR flag
    if prevL > nLv or nLv <= 1:
        crflag = True
    else:
        crflag = False
    
    if smryflag and keyword in inkeywords:
        #smry += "["+str(i)+"L"+str(+nLv)+"]"+str(value)+", "
        smry += str(value)+", "

    if vr == "SQ" and crflag:   # VR is SQ, crflag is True
        smry += "\n"            # Then add CR in smry

    if vr == "SQ":              # VR is SQ,
        prevL = nLv             # Then assign now level to prevL

smry += rmCRline(smry)

with open("rdsrResult.txt","w") as f:
    f.write(rstxt)
with open("rdsrDataset.txt","w") as f:
    f.write(dstxt)
with open("rdsrSummary.csv","w") as f:
    f.write(smry)



さて、ここで問題です。変数smryに値を追加していくだけですと、抽出したい情報がたった1行で出力されてしまいます。
そこで、意味ある情報の塊の境界で改行を入れる必要があります。

改行するかどうかを変数crflagを用いて制御しています。
    if prevL > nLv or nLv <= 1:
        crflag = True
    else:
        crflag = False

これは、1つ前の"SQ"のデータエレメントの階層レベルprevLが現在の階層レベルnLvより大きい(つまり現在の階層は前の階層より上がるとき)、または現在の階層が1以下のとき(つまり階層が0または1の深さにある)ときはcrflagをTrueにしています。
それ以外はcrflagをFalseに設定します。

そして、VRの値とcrflagの状態を見て変数smryに改行を入れています。
    if vr == "SQ" and crflag:   # VR is SQ, crflag is True
        smry += "\n"            # Then add CR in smry


このままでは改行が至る所に入ってしまうため、無駄な改行を削除するためrmCRlineメソッドを作成しました。
    smry += rmCRline(smry)



実行してみる


プログラム実行してみましょう
初めに基本的な情報を表示して停止します。


SOP Class UID: 1.2.840.10008.5.1.4.1.1.88.67
Modality     : SR
Manufacturer : SIEMENS
Study Date   : 20220224
Pt Name      : NO^Name
Pt ID        : 12345678
Modality     : SR
Content Sequence: 20

pause. Type ENTER key


ここで、「ENTER」キーを打つと、処理が次に進みます。

                ・・・・・・・
      (0008, 0104) Code Meaning                        LO: 'Source of Dose Information'
      ---------
   (0040, a168)  Concept Code Sequence  1 item(s) ---- 
      (0008, 0100) Code Value                          SH: '113856'
      (0008, 0102) Coding Scheme Designator            SH: 'DCM'
      (0008, 0104) Code Meaning                        LO: 'Automated Data Collection'
      ---------
   ---------

Data type : <class 'pydicom.dataset.FileDataset'
                --> <class 'str'>
len(ds)= 3159

pause. Type ENTER key


読み込んだデータセットのデータの型が'pydicom.dataset.FileDataset'から'str'に変換されました。
「ENTER」キーを打つと、処理が次に進みます。

・・・・・・・
value = >  CONTAINS <
value = >  CODE <
value = > 1 item(s) <
value = >  113854 <
value = >  DCM <
value = >  Source of Dose Information <
value = > --------- <
value = > 1 item(s) <
value = >  113856 <
value = >  DCM <
value = >  Automated Data Collection <
value = > --------- <
value = > --------- <


最後に3つのファイルが出力されます。
rdsrDataset.txt        データエレメントに行番後を付加して出力します。
rdsrResult.txt          getInfoDCMメソッドで抽出した結果を行番号と一緒に出力します。
rdsrSummary.txt     必要な情報を行単位で出力します。


summayの結果を見てみる


rdsrSummary.csvファイルの出力結果は次のようになります。

最初に基本的な情報が表示されます。


以下ではCTのX線量に関する情報が単位とともに表示されています。


























おわりに


今回、2022年4月24日に日本放射線技術学会 北海道支部 第78回春季学術大会 医療情報セミナーで発表したものを掲載しました。

DICOM SRの読み込みはPydicomのサンプルrtplan.srで試していましたが、今回RDSRは初めての経験でした。
読み込むまでは良かったのですが、階層が6階層まで下りて行く必要があり、Pydicomの標準機能ではうまく情報を取り込むことができませんでした。
CTのほかにAG(血管造影)のRDSRも解析しましたが、同じ結果でした。

本来はDICOM規格をしっかり理解したうえでプログラミングすべきでしたが、短期間に結果を出す必要があり、力ずくの対応でした。
したがってまだ気づかぬバグなどが潜んでいるかもしれません。プログラムの使用には注意をお願いします。


2021年6月13日日曜日

 

Pyゼミ1.07 CT画像にDICOM Color LUTを使ってヒートマップ画像などを重ねてFusion画像を作ってみる

 

今回は何らかの結果(ヒートマップや推論マップなど)とその元画像を重ねるプログラムを作成します。

何らかの結果にはPydicomが提供するとColor LUTを利用してカラー画像に変換します。

Pyゼミ1.06ではImageJのColor LUTを使いましたが、ここではPydicomが提供するColor LUTを使います。

Color LUTを適用した画像をDICOM画像に重ねてFusion画像を作成してみます。


医療の世界ではFusionと呼んでいますが、Computer Visonの世界ではSuperimposeと呼んでいます。

 

具体的には下に示すような画像を作成します(グレイのCT画像に模擬的な結果のカラー画像を重ねています)。

     

fusion.png

 

ポイント

  • 何らかの結果を疑似的に2次元ガウス分布を用いて作成する。
  • pydiocmLUT適応機能を使ってグレイ画像をカラー画像に変換する。
  • グレイ画像のDICOM画像を読み込み、RGBのモノクロ画像(3チャンネル)を作成する。
  • OpneCVaddweighted()メソッドを使ってCT画像と重ね合わせる。


プログラムはここから

サンプルのDICOM画像はここから


まずは必要なライブラリについて

以下のライブラリをインポートします。


import cv2, math

import numpy as np

import pydicom

from pydicom.pixel_data_handlers.util import apply_color_lut



pydicomが提供するLUTを設定する

pydicom8つの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

 lutnm = lutarr[0]           # Set LUT Name with index of lutarr


LUTの名称をindexの値で変数lutnmに格納しておくと、数値を変えるだけでLUTを変更することができて便利です。


 

ガウス分布を作成する関数gaussianを作成する

 

この関数は引数にガウス分布の中心の座標、cxcyと標準偏差sgm、そし生成する画像のサイズsizeを与えます。

はじめに2次元のガウス分布を用いてグレイ画像を作成するために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


ガウス分布の中心座標cxcy、ガウス分布の広がり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は以下のようになります。

mapOrg.png

 

この画像にPydicomが提供するColor LUTを適用してみましょう。

 

cimg = apply_color_lut(pmap, palette=lutnm)

cimg = cimg[ : , : , : : -1]

cv2.imwrite("mapLut.png", cimg)

 

カラー画像保存する際、BGRRGBに変換するので cimg = cimg[ : , : , : : -1] とします。

 

結果の画像mapLut.pngは次のように表示されます。

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は以下のように表示されます。


mapLut2.png

 

DICOM画像を読み込む

それでは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を使ってウインドニング処理を行います。

また、デフォルトfloat32uint8の符号なし8ビットに変換しているます。

この画像のshapeとタイプを表示すると次のように表示されます。

 

>1> image shape= (512, 512) uint8

 

shape(512,512)と表示され、512×512のグレイ画像であることを示しています。

Fusionを行うためには2つの画像がカラー画像であることが必要です。

 

グレイ画像をカラー画像に変換する3つの方法

 1チャンネルのグレイ画像を3チャンネルのカラー画像に変換する方法を紹介します。


方法1:保存した画像をimread()メソッドで読み込む。

CV2imreadのデフォルトはグレイ画像もカラー画像として読み込まれます。

 

img2 = cv2.imread("./brain01.png")

print(">2> image2 shape=", img2.shape, img2.dtype)

 

結果は以下のように、読み込まれた画像img23チャンネルである、つまりRGBのカラー画像であることを示しています。

 

>2> image2 shape= (512, 512, 3) uint8


明示的にimread()グレイ画像として読み込むときには引数にcv2.IMREAD_GRAYSCALEを使って、v2.imread(filename, cv2.IMREAD_GRAYSCALE)として読み込むか、cv2.IMREAD_GRAYSCALEの代わりに を使ってもOKです。

 

方法2cv2merge()メソッドを使います。

merge()メソッドを使って同じ3枚の画像imgを統合します。

 

# 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


逆に複数のチャンネルをもつ画像を分割するにはsplit()メソッドを使って、r, g, b = cv2.split(img) のように記述します。

  

方法3numpyzerosメソッドを使って3チャンネルの画像を作成する。

 zeros()メソッドを使い、512×512の3チャンネル、符号なし8ビットの画像img4をあらかじめ作成します。


# 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の012の各チャンネルに画像imgを代入します。

実行結果は3チャンネルの画像img4ができています。


>4> image4 shape= (512, 512, 3) uint8

 

画像にカラースケールの背景を入れる

3チャンネルに変換した画像img2にカラースケールの発色をよくするため画素値128で埋めます。

これを入れないとカラースケールが暗くなります。

 

# For color scale

img2[350-256 : 350, : 15] = 128

 


2つの画像を重ね合わせる

それではaddWeighted()メソッドを使って2つの画像img2cimgを重ね合わせます。

 

# 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は次のように表示されます。

fusion.png

 

addWeighted()メソッドの引数alphabetaは2つの画像の重ねるときの重みになります。

主たる画像、ここではBrain CT画像のimg21.0とし、ガウス分布画像のcimg0.5と重みを小さくして重ね合わせています。

bataの値については場合により調整してみましょう。

 

 

Post Script:

実はaddWeighted()メソッドを使った画像の重ね合わせは結構はまりました。

大前提は2つの画像はカラー画像であることです。

しかし、DICOM画像をpydicomで読み込んだときは1チャンネルの画像、その画像をいったん保存し、imread()メソッドで読み込むとデフォルトでカラー画像として読み込まれます。

したがって、ある時はaddWeighted()で上手く重なり、あるときは上手くいかない。

この原因を明らかにするのにはまってしまいました。気を付けましょう。