DICOM画像を用いてopenCVの領域抽出をやってみた① | 診療放射線技師がPythonをはじめました。

3回前の記事で、openCVの領域抽出を汎用画像を用いて行いました。今回は、その領域抽出をDICOM画像を用いてやってみたいと思います。

Contents

1 プログラムの流れ
2 前回のコード
3 コードを組んでいく

プログラムの流れ

プログラム作成の流れを組んでみたいと思います。

まずは、前回組んだ

領域抽出の設定、閾値設定を確認できるプログラムを組んでみたをベースに機能拡張していきたいと思います。

画像表示部分を一つ追加し、そこに領域抽出した画像を表示します。

また、openCVの領域抽出は4種類ありますので、それら抽出方法を選べるように4つのボタンを配置します。

また、領域抽出した領域の面積および、外周の長さを表示するようにします。

最後に、プログラムの組み方が悪いのか、私のPCのスペックが低いのか、スケールを動かした時の動作が遅くなるので、それぞれのスケールの値を一定数増加、減少することができるボタンを配置したいと思います。

画像領域を一つ追加
４種類の領域抽出ボタンを追加
抽出結果を表示
スケールの増減ボタンを配置

最終的には以下のようになります。

前回のコード

前回のコードは以下になります

import tkinter as tk
import cv2
import matplotlib.pyplot as plt
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
import numpy as np
import pydicom
import fileselect as fs


def level():
    global ww_low, ww_high, w_level, w_width

    window_level, window_width = level_sc.get(), w_width_sc.get()

    ww_low = window_level – (window_width // 2)
    ww_high = window_level + (window_width // 2)

def window(self):
    global ww_low, ww_high, img_arr, ax1

    level()

    ax1.imshow(img_arr, cmap=’bone’, vmin=ww_low, vmax=ww_high)
    fig.canvas.draw()

def thresho(self):
    global img_arr, ax2, fig

    ret, thresh = cv2.threshold(img_arr, int(self), 255, cv2.THRESH_BINARY)
    ax2.imshow(thresh, cmap=’bone’)

    fig.canvas.draw()


def main():
    global ww_low, ww_high, img_arr, thre_sc, ax2, fig,\
        img_arr, ax1, level_sc, w_width_sc, ww_low, ww_high, w_level, w_width


    filename = fs.single_fileselect()
    dcm = pydicom.dcmread(filename)
    img_arr = np.array(dcm.pixel_array)
    w_level = int(dcm&#91;0x0028,0x1050].value)
    w_width = int(dcm&#91;0x0028, 0x1051].value)

    fig = plt.figure(figsize=(12, 6))
    ax1 = fig.add_subplot(1, 2, 1)
    ax2 = fig.add_subplot(1, 2, 2)

    ax1.axes.xaxis.set_visible(False), ax1.axes.yaxis.set_visible(False)
    ax2.axes.xaxis.set_visible(False), ax2.axes.yaxis.set_visible(False)


    root = tk.Tk()
    root.geometry(“1410×600”)

    Canvas = FigureCanvasTkAgg(fig, master=root)
    Canvas.get_tk_widget().grid(row=0, column=0, rowspan=10)


    var_scale_level = tk.IntVar()
    level_sc = tk.Scale(root,
        label=’Window Level’,
        variable=var_scale_level,
        orient=tk.HORIZONTAL,
        length=200,
        from_= np.min(img_arr),
        to=np.max(img_arr),
        command=window)
    level_sc.set(w_level)
    level_sc.grid(row=1, column=1)

    var_scale_w_width = tk.IntVar()
    w_width_sc = tk.Scale(root,
        label = ‘Window Width’,
        variable=var_scale_w_width,
        orient=tk.HORIZONTAL,
        length=200,
        from_=0,
        to=(np.max(img_arr) – np.min(img_arr))//2,
        command=window)
    w_width_sc .set(w_width)
    w_width_sc.grid(row=3, column=1)


    thre = 0

    var_thre = tk.IntVar()
    thre_sc = tk.Scale(root,
        label=’Threshold’,
        variable=var_thre,
        orient=tk.HORIZONTAL,
        length=200,
        from_=np.min(img_arr),
        to=np.max(img_arr),
        command=thresho)
    thre_sc.set(thre)
    thre_sc.grid(row=5, column=1)

    window(1)
    thresho(1)

    root.mainloop()

if __name__ == “__main__”:
    main()

100

101

102

103

104

105

106

import tkinter as tk

import cv2

import matplotlib.pyplot as plt

from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg

import numpy as np

import pydicom

import fileselect as fs

def level():

global ww_low, ww_high, w_level, w_width

window_level, window_width = level_sc.get(), w_width_sc.get()

ww_low = window_level – (window_width // 2)

ww_high = window_level + (window_width // 2)

def window(self):

global ww_low, ww_high, img_arr, ax1

level()

ax1.imshow(img_arr, cmap=‘bone’, vmin=ww_low, vmax=ww_high)

fig.canvas.draw()

def thresho(self):

global img_arr, ax2, fig

ret, thresh = cv2.threshold(img_arr, int(self), 255, cv2.THRESH_BINARY)

ax2.imshow(thresh, cmap=‘bone’)

fig.canvas.draw()

def main():

global ww_low, ww_high, img_arr, thre_sc, ax2, fig,\

img_arr, ax1, level_sc, w_width_sc, ww_low, ww_high, w_level, w_width

filename = fs.single_fileselect()

dcm = pydicom.dcmread(filename)

img_arr = np.array(dcm.pixel_array)

w_level = int(dcm[0x0028,0x1050].value)

w_width = int(dcm[0x0028, 0x1051].value)

fig = plt.figure(figsize=(12, 6))

ax1 = fig.add_subplot(1, 2, 1)

ax2 = fig.add_subplot(1, 2, 2)

ax1.axes.xaxis.set_visible(False), ax1.axes.yaxis.set_visible(False)

ax2.axes.xaxis.set_visible(False), ax2.axes.yaxis.set_visible(False)

root = tk.Tk()

root.geometry(“1410×600”)

Canvas = FigureCanvasTkAgg(fig, master=root)

Canvas.get_tk_widget().grid(row=0, column=0, rowspan=10)

var_scale_level = tk.IntVar()

level_sc = tk.Scale(root,

label=‘Window Level’,

variable=var_scale_level,

orient=tk.HORIZONTAL,

length=200,

from_= np.min(img_arr),

to=np.max(img_arr),

command=window)

level_sc.set(w_level)

level_sc.grid(row=1, column=1)

var_scale_w_width = tk.IntVar()

w_width_sc = tk.Scale(root,

label = ‘Window Width’,

variable=var_scale_w_width,

orient=tk.HORIZONTAL,

length=200,

from_=0,

to=(np.max(img_arr) – np.min(img_arr))//2,

command=window)

w_width_sc .set(w_width)

w_width_sc.grid(row=3, column=1)

thre = 0

var_thre = tk.IntVar()

thre_sc = tk.Scale(root,

label=‘Threshold’,

variable=var_thre,

orient=tk.HORIZONTAL,

length=200,

from_=np.min(img_arr),

to=np.max(img_arr),

command=thresho)

thre_sc.set(thre)

thre_sc.grid(row=5, column=1)

window(1)

thresho(1)

root.mainloop()

if __name__ == “__main__”:

main()

コードを組んでいく

それでは、流れに沿ってコードを組んでいきたいと思います

画像領域を一つ追加

ここは単純に上記コードの49行目と52行目にax3として追加していきます。

追加コードは以下

ax3 = fig.add_subplot(1, 3, 3)

ax3.axes.xaxis.set_visible(False), ax3.axes.yaxis.set_visible(False)

ax3の追加に伴い47行目、48行目は真ん中の数字が”２”から”３”に変更になります。

ax1 = fig.add_subplot(1, 3, 1)

ax2 = fig.add_subplot(1, 3, 2)

４種類の領域抽出ボタンを追加

tkinterのボタンの設置は

ボタンの名前 = tkinter.Button(配置するフィールド, text=”ボタン上に表示する文字”, command=ボタンを押した時の動作, width = ボタンの幅)

で設定し、

ボタンの名前.grid(row=8, column=1)

で表示場所の指定をします。

今回はopenCVの領域抽出方法である”EXTERNAL”、”LIST”、”CCOMP”、”TREE”の4つのボタンを配置します。

コードは以下となります。

ext_but = tk.Button(root, text=”EXTERNAL”, command=external, width = 10)
ext_but.grid(row=8, column=1)

list_but = tk.Button(root, text=”LIST”, command=list, width = 10)
list_but.grid(row=8, column=2)

ccomp_but = tk.Button(root, text=”CCOMP”, command=ccomp, width = 10)
ccomp_but.grid(row=9, column=1)

TREE_but = tk.Button(root, text=”TREE”, command=tree, width = 10)
TREE_but.grid(row=9, column=2)

ext_but = tk.Button(root, text=“EXTERNAL”, command=external, width = 10)

ext_but.grid(row=8, column=1)

list_but = tk.Button(root, text=“LIST”, command=list, width = 10)

list_but.grid(row=8, column=2)

ccomp_but = tk.Button(root, text=“CCOMP”, command=ccomp, width = 10)

ccomp_but.grid(row=9, column=1)

TREE_but = tk.Button(root, text=“TREE”, command=tree, width = 10)

TREE_but.grid(row=9, column=2)

上記コードを

window(1)

の前、大体100行目辺りになると思います。に記載しましょう。

抽出結果を表示

抽出結果を一番右側の画像表示領域に表示するコードに入ります。

初めに、領域を抽出するのはDICOM画像の配列でもできますが8ビット階調に変換する必要があります。

結果を画像上に重ね合わせるのは8ビット階調の汎用画像（jpgや、png画像）でなければならないという事。

なので、まずは4つのボタンが押された時に、まず汎用画像を作成する事を行わなければなりません。

手間は増えますが、階調を合わせた一番左の画像に重ね合わせたいので一番左の画像を一度汎用画像として保存してそれを使いたいと思います。

汎用画像に保存

一番左側の画像条件で画像を一度汎用画像に保存します。DICOM画像は16ビット画像なので8ビット画像に変換する必要があります。

とりあえず、画像配列をコピーしウインドウ幅の最低値以下のピクセル値を最低の値に、ウインドウ幅の最高値以上のピクセル値を最大の値に変えます。

その後、ウインドウ幅を0～256の階調に割り振ります。

そして、img_8bit.pngという名前で、プログラムがあるフォルダに保存します。

これらコードを関数としておきましょう。今回はcontという名前で作成しました。


def cont():
    img_unit8 = copy.deepcopy(img_arr)　#img_unit8としてコピー

    np.clip(img_unit8, ww_low, ww_high, out= img_unit8)
          #コピーした配列の最低値以下をウインドウ幅最低値に、
          　#　最高値以上も同様に）

    img_unit8 -= img_unit8.min()
          #配列全体をウインドウ幅最低値を引くことで0からに変更）

    np.floor_divide(img_unit8, (img_unit8.max() + 1) / 256,
                out = img_unit8, casting=’unsafe’)
          #ウインドウ幅を256分割する。

    cv2.imwrite(‘./img_8bit.png’,img_unit8)

def cont():

img_unit8 = copy.deepcopy(img_arr)　#img_unit8としてコピー

np.clip(img_unit8, ww_low, ww_high, out= img_unit8)

#コピーした配列の最低値以下をウインドウ幅最低値に、

　#　最高値以上も同様に）

img_unit8 -= img_unit8.min()

#配列全体をウインドウ幅最低値を引くことで0からに変更）

np.floor_divide(img_unit8, (img_unit8.max() + 1) / 256,

out = img_unit8, casting=‘unsafe’)

#ウインドウ幅を256分割する。

cv2.imwrite(‘./img_8bit.png’,img_unit8)

領域抽出のプログラム

領域抽出は、閾値設定した配列　”thresh” で行います。(以前のコード29行目で定義）以前のプログラムで閾値設定した配列は ”thresh” という変数名で持っています。

しかし、先ほども記載しましたが領域抽出は8ビット階調の配列でなければいけないので16ビットの配列を8ビット階調に変換します。

thresh_8bit = thresh.astype(‘u1’)

1	thresh_8bit = thresh.astype(‘u1’)

ここで、”u1”とは符号なしの8ビット整数型を示します。

今、階調変換した配列を用いて領域抽出を行います。

ret, contours, hierarchy = cv2.findContours(thresh_8bit, 0, cv2.CHAIN_APPROX_SIMPLE)

1	ret, contours, hierarchy = cv2.findContours(thresh_8bit, 0, cv2.CHAIN_APPROX_SIMPLE)

最新バージョンのfindContoursの戻り値は3つから2つに変わり、retの戻り値は廃止になった模様です。もし、エラーが出た場合、retを削除してから試してください。

領域抽出の結果はcontoursに入っています。それらを先ほど保存した汎用画像の上に表示していきます。

表示には

cv2.drawContours(img_png, contours, -1, (255, 255, 0), 1)

で表示できますので、これをfor文で表示していきます。drowContoursについて詳しくはチュートリアルをご覧ください

以上の工程を領域抽出方法4個分の関数として作成していきます。

例えば、 cv2.RETR_LISTの場合

def list():
    cont()
    #ping画像の作成関数
    img_png = cv2.imread(‘img_8bit.png’)  #作成した画像の読み込み

    thresh_8bit = thresh.astype(‘u1′)
 #閾値配列を8ビット化

    ret, contours, hierarchy = cv2.findContours(thresh_8bit, 1,cv2.CHAIN_APPROX_SIMPLE)

    for contour in contours:
        cv2.drawContours(img_png, contours, -1, (255, 255, 0), 1)
    ax3.imshow(img_png, cmap=’bone’)

    fig.canvas.draw()

def list():

cont()

#ping画像の作成関数

img_png = cv2.imread(‘img_8bit.png’) #作成した画像の読み込み

thresh_8bit = thresh.astype(‘u1’)

#閾値配列を8ビット化

ret, contours, hierarchy = cv2.findContours(thresh_8bit, 1,cv2.CHAIN_APPROX_SIMPLE)

for contour in contours:

cv2.drawContours(img_png, contours, –1, (255, 255, 0), 1)

ax3.imshow(img_png, cmap=‘bone’)

fig.canvas.draw()

ちなみに、findContoursの第2引数は領域抽出の方法なのですが、4つの種類が用意されていて以下の様になっています。

抽出方法	番号	内容
cv2.RETR_EXTERNAL	0	輪郭のうち、最も外側のみ抽出
cv2.RETR_LIST	1	全ての輪郭を抽出
cv2.RETR_CCOMP	2	全ての輪郭を２レベルの階層に分けて出力
cv2.RETR_TREE	3	抽出した内側の輪郭も抽出

抽出方法の”cv2.RETR_LIST”などの文字列でも指定できますが、番号でも指定できます。

4つの抽出方法の関数は以下のようになります。


def external():
    img_png = cv2.imread(‘img_8bit.png’)
    thresh_8bit = thresh.astype(‘u1′)
    ret, contours, hierarchy = cv2.findContours(thresh_8bit, 0,   cv2.CHAIN_APPROX_SIMPLE)

    for contour in contours:
        cv2.drawContours(img_png, contours, -1, (255, 255, 0), 1)

    ax3.imshow(img_png, cmap=’bone’)
    fig.canvas.draw()

def list():
    img_png = cv2.imread(‘img_8bit.png’)
    thresh_8bit = thresh.astype(‘u1′)
    ret, contours, hierarchy = cv2.findContours(thresh_8bit, 1,  cv2.CHAIN_APPROX_SIMPLE)

    for contour in contours:
        cv2.drawContours(img_png, contours, -1, (255, 255, 0), 1)

    ax3.imshow(img_png, cmap=’bone’)
    fig.canvas.draw()

def ccomp():
    img_png = cv2.imread(‘img_8bit.png’)
    thresh_8bit = thresh.astype(‘u1′)
    ret, contours, hierarchy = cv2.findContours(thresh_8bit, 2, cv2.CHAIN_APPROX_SIMPLE)

    for contour in contours:
        cv2.drawContours(img_png, contours, -1, (255, 255, 0), 1)

    ax3.imshow(img_png, cmap=’bone’)
    fig.canvas.draw()

def tree():
    img_png = cv2.imread(‘img_8bit.png’)
    thresh_8bit = thresh.astype(‘u1′)
    ret, contours, hierarchy = cv2.findContours(thresh_8bit, 3,   cv2.CHAIN_APPROX_SIMPLE)

    for contour in contours:
        cv2.drawContours(img_png, contours, -1, (255, 255, 0), 1)

    ax3.imshow(img_png, cmap=’bone’)
    fig.canvas.draw()

def external():

img_png = cv2.imread(‘img_8bit.png’)

thresh_8bit = thresh.astype(‘u1’)

ret, contours, hierarchy = cv2.findContours(thresh_8bit, 0, cv2.CHAIN_APPROX_SIMPLE)

for contour in contours:

cv2.drawContours(img_png, contours, –1, (255, 255, 0), 1)

ax3.imshow(img_png, cmap=‘bone’)

fig.canvas.draw()

def list():

img_png = cv2.imread(‘img_8bit.png’)

thresh_8bit = thresh.astype(‘u1’)

ret, contours, hierarchy = cv2.findContours(thresh_8bit, 1, cv2.CHAIN_APPROX_SIMPLE)

for contour in contours:

cv2.drawContours(img_png, contours, –1, (255, 255, 0), 1)

ax3.imshow(img_png, cmap=‘bone’)

fig.canvas.draw()

def ccomp():

img_png = cv2.imread(‘img_8bit.png’)

thresh_8bit = thresh.astype(‘u1’)

ret, contours, hierarchy = cv2.findContours(thresh_8bit, 2, cv2.CHAIN_APPROX_SIMPLE)

for contour in contours:

cv2.drawContours(img_png, contours, –1, (255, 255, 0), 1)

ax3.imshow(img_png, cmap=‘bone’)

fig.canvas.draw()

def tree():

img_png = cv2.imread(‘img_8bit.png’)

thresh_8bit = thresh.astype(‘u1’)

ret, contours, hierarchy = cv2.findContours(thresh_8bit, 3, cv2.CHAIN_APPROX_SIMPLE)

for contour in contours:

cv2.drawContours(img_png, contours, –1, (255, 255, 0), 1)

ax3.imshow(img_png, cmap=‘bone’)

fig.canvas.draw()

ちょっと、長くなってしまいました。

今回はここまでとします。

次回、４つのボタンから最後に作成した関数に繋げて、領域抽出を行いその結果を表示するところまでやっていきたいと思います。

お疲れ様でした。