Re - ImageJで学ぶ!: 2015-09-20

2015年9月26日土曜日

第26回 連続画像処理で学ぶ!

前回まで、主に一枚の単一画像の処理や解析について説明してきました。
ImageJは、1つのウインドウに2枚以上の複数画像を表示し、処理や解析などを行うことが可能です。
これらの複数画像を一度に処理するには、画像群を束にして処理する(stack:スタック)処理を利用します。
lmageJは、このスタックに対しても、入出力から動画像フォーマット変換処理や三次元画像変換など多彩な機能を提供してくれています。今回はImageJの連続画像処理について、入門的な内容について述べます。

ImageJによるスタック処理機能


スタックを形成するひとつひとつの画像をスライス、およびフレームと言います。スタックにおけるすべてのスライスは、同じマトリックスサイズで、濃淡値(ビット)も同じでなければなりません。ImageJのメニュー〈File/lmport/lmage Sequence>はスタックとしてフォルダ内の画像を開くことができます。連続するスライスのDICOMファイルもスライス番号ごとに整理されて表示されます。画像は、表示された画面ウィンドウ下部のスクロールバーでめくっていく(スライス送り)ことができます。


(目的のシリーズが保存されたフォルダを選択)
(任意のスライスをスクロール表示)

また、ImageJにおけるフィルタ処理などは、スタックのすべてのスライスに対して実行されます。
次の図は、モンタージュを表示した例です。スタック処理は〈lmage/Stacks/_〉の下に基本処理が用意されています。


画像のリスライスも可能です。


この例では、CT画像のスライス間隔をDICOMヘッダ情報からキャッチして、(lmage/Show Info_でヘッダ情報を表示できる)、〈lmage/ Stack/Reslice>を選び、元画像のスライス間隔以下の値をOutput Z Spacingの欄に入力することでMPR像の奥行きのスライス間隔を決定します。Avoid interpolationを選択すると、画像間の重なりのないスタックが作成できますが、元画像のスライス厚よりも薄いスライス厚で再構成する場合には、チェックを外しましょう(画像が増えた分、引き伸ばされた画像になってしまうので)。

最後に、元画像のリスライス方向を決めて実行すると、MPRが完成します。この例では、「Left」から切り出す設定にしているので、サジタル画像が作成できていることがわかります。

スタックを利用した基本三次元処理機能


スタック画像を表示した後、〈lmage/ Stacks>を選択すると、いくつかの三次元画像変換処理が実装されています。
次の図はOrthogonal ViewsのMPR例です。任意の基準線が同期します。





次に、最大値投影法(maximum intensity projection)、最小値投影法、Ray-Sum法、平均値投影法、標準偏差投影法および中央値投影法を説明します。
これらの処理はあらかじめMPRスタックをつくってから、〈lmage/Stack/Z-Projection>を選択します。不要なスライスは、投影像のノイズになるため再構成に含めないように設定できます。

次の図は、リスライスで作成したPET画像のサジタル画像を最大値、平均値、最小値および標準偏差値を投影した画像例です。画像作成時に、関心領域を含む画像のみを指定しています。


(MIP)

(AVG)

(MIN)

(STD)

この中で、市販されている3Dワークステーションには、あまり標準搭載されていない機能として、標準偏差値投影法があります。これは、最大値投影法で目的とする関心領域がわかりずらい場合に、その補助的な画像として、関心領域とその周囲組織との全体像を把握するのに良好な画像を提供します。

また、ImageJの標準機能では、一方向のMPRだけではなく、360度方向に回転させ、動画表示やその動画の保存が可能です。次の図は、くlmage/Stack/3D-Project...〉より、〈Brightest Point(最大値投影法)〉を選択し、Y軸を中心に10 °ずつ回転させてアニメーションを作成した例です。スタックを作成後、Save as...で"avi"を選択しました。



ImageJの複数画像の計測処理機能について


ImageJは複数画像を一挙に計測する機能を多く有しています。
次の図は、スタックした画像群のうちの1画像にROIを設定し、そのROIと同じ位置の平均値をすべての画像について自動計算しプロットする“Z-axis Profle”機能(lmage/Stack/Plot Z-axis Profile_)を施工した例です。


(188枚目にROIを設置)
各スライスにおけるROIの平均信号値)

(その詳細:Result機能)

グラフはZ軸方向のスライスがX軸、平均値がY軸としてプロットされています。プロットしたグラフに、近似曲線のカーブフィッティングを施工したい場合は、そのグラフを表示したまま〈Analyze/Tools/CurveFitting...〉を選択します。
CurveFitterダイアログが表示されるので、プロットの結果(Resultテーブル)を一度SaveAsで保存し、エクセルファイルにしてから、任意の値をこの中にコピー&ペーストし、相当する近似関数をプルダウンメニューから選択してFitボタンでフィッティング、もし、グラフ画像のみ欲しい場合はさらにApplyボタンを押します。

(ResultテーブルのSet measurementを調整して、不要な項目を非表示にすれば、Resultテーブルから直接コピペもできます。)

(エクセルから貼り付け:スライス番号対最大値)

(Fit実行:線形近似)

(Fit実行:ガウシアン分布)


この際、R2値やそのほかの計算のための詳細パラメータは、logダイアログ内にすべてリスト表示されます。近似曲線は、その 近似曲線R2値が 1 に近い場合に最も信頼性が高くなります。上の例では、線形近似よりも、ガウシアン分布の方が適していることがわかります。ちょうど、山になっている部分が、先ほど示した関心領域のあるスライス位置です。
このテクニックは、よく造影剤のダイナミック解析や臨床における日常の画質管理などで使用される計測処理であるので覚えておくと便利です。

また、メニューのProcess下にある機能のフィルタや2値化など、スタックに対する処理の際、スタック処理をするか、スタック内に表示した1画像だけに処理をかけるかは、その都度、ポップアップでメッセージが表示されます。

以上、今回はImageJによる複数画像に対する基本的かつ入門的な処理方法を紹介しました。
次回はImageJの動画フォーマットの取り扱いとそれに関連するプラグインの説明を行う予定です!

Reference
  • 山本修司:ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(22・3) 2007, p104-107」

2015年9月25日金曜日

第25回 2値化の応用 医用画像の特徴量解析(フラクタル解析)で学ぶ!

今回は、2値画像の医用への応用の1例として、フラクタル解析を紹介します。
画像の特徴量解析は医用画像上に表示されるさまざまな症例の形状や濃度情報の特徴を定量評価します。例えば、がんの良性・悪性の判断材料や、時系列で形状や濃淡値が変化していくような病気の定量解析に有用です。今回は2値画像を用いたフラクタル解析の中で最も簡単で単純なアルゴリズムであるボックスカウンティング(box-counting)について、 lmageJを用いて解説していきます 

フラクタルとは


フラクタルは、自己相似性を持つ図形や構造などの総称で、1975年に数学者のべノア・マンデルブロ(Benoit B. Mandelbrot)氏によって作られた造語です
この「自己相似性」という用語は、図形の一部を拡大すると、ほかの部分または全体と形が一致する性質のことで、その中でも、図形のどの部分をとっても細部を拡大すると全体と同じ形になる、一様な自己相似性をフラクタル(fractal)と言います。

自己相似性を持つ図形には、曲線をつくるものと、特定の図形内で分割していくものとに分けられます。曲線をつくるものには、コッホ曲線・コッホ雪片、分割していくものにはシェルピンスキーのギャスケットなどがあります。

本稿では、構造の複雑さを示すことを目的として、前者のフラクタル図形の最も標準的な曲線である コッホ(koch)曲線を用いて自己相似性を説明します。

コッホ曲線

コッホ曲線は具体的には以下のような手順で生成できます。 

  1. 適当な長さの線分を引く。(これをイニシエータと呼びます)。
  2. その線分を3等分し、中央部のみジェネレータ(三角形の底辺なしの図形) に置き換える。
  3. 4本の線分ができているので、それぞれの線分を3等分し、それぞれの中央をジェネレータに置き換える。
  4. 16本の線分ができているので,それ をまたジェネレータに置きi換える。 
  5. 以降、この操作を繰り返すと、上述のような無限段階のコッホ曲線になる。 


このように、線分を3等分し、分割した2点を頂点とする正三角形の作図を無限に繰り返すことによって得られる図形がコッホ曲線です。

ここで、この図形がどのくらい複雑になったかを示すフラクタル次元というものを算出します。正確には、より細かなスケールへと拡大するにつれあるフラクタルがどれだけ完全に空間を満たしているように見えるかを示す統計的な量がフラクタル次元です。

ラクタル次元Dは、 元の図形のサイズ1/nの相似比、かつ、できた同じ図形の個数をmとすると、次の式で表せます。

D = lognm・・・式(1)

したがって、上記のコッホ曲線の場合、相似比は、線分の長さが1回の操作ごと1/3となっており、その個数が常に4つずつ増えていますので、D=log4/log3= 1.26次元であることがわかります。線は一次元、面は二次元なので、コッホ曲線は少し複雑な次元の線です。

このように、フラクタル次元は図形の複雑さを数値で表していると言えます。また、正方形の次元を考える場合、普通は直感で二次元とわかりますが、フラクタル次元で説明すると、例えば、正方形の各辺を2分すると4つの相似図形になることを考えると、式(1)から、D=log4なので二次元となります。同じように、立方体の次元を計算すると、三次元であることが式 (1)から計算できます。

しかし、フラクタル次元のほとんどは、無理数(分子・分母ともに整数の分数ではない数) になります。フラクタル次元は、数学的に厳密に定義するのは非常に難しいと言われていますが、医学を含め、経済学や自然科学など広い分野で応用されているようです。 

生体のフラクタルについて


人間の体は、ガス交換や栄養物質の血液中の交換などを効率良く行わなくてはならない構造になっています。そのため、限りある体積に最大限の表面積を包含できる組織構造が必要で、その代表的なものとして、肺構造、毛細血管の分岐構造や腸絨毛構造はフラクタル構造(細部の形状が一様に連なっている)であることが知られています。

ImageJによる医用画像のフラクタル次元の測定方法


医用画像のフラクタル次元を測るために、ImageJにも搭載されているシンプルな方法「ボックスカウンティング法」を紹介します。この方法は2値化画像を処理します。ボックスカウンティング法は、解析対象の画面上にボックスを設定し、そのボックスに含まれる信号値を持つ部分の数を数え(N(x))、その大きさ(一辺の長さx)を変化させていきます。そして得られたカウントを両対数 [log(X) v.s. log(N(x))]グラフとし、直線近似してその直線の傾きDをフラクタル次元とします。

このImageJのボックスカウンティング機能はソースコードが公開されています。(FractalBoxCounter.java) “ijmeasure.CurveFitter”クラスでlog換したボックスサイズ“sizes” “ boxCountSum”の値を引数として代入して初期化し、“ij.measure.CurveFitter. doFit(int fitType)”の引数として “STRAIGHT LINE”を代入して直線近似するアルゴリズムになっています。

さらに、傾きを求めるために、上述の”CurveFitter”クラスの“getParams();”メソッドで傾きDを抽出して、フラクタル次元としています。 

それでは,ImageJを用いて実際にフラクタル解析をやってみましょう。
画像はCTから得られた肺がん(扁平上皮癌)例です。


手順は2つあります。共通の手順は、CT画像を読み込んで、肺がん部分を切り出し、8-bitグレイスケールに変換しておくところまで同じです。
(Process>Binaryの機能も活用していきます!)

1つ目の手順は、ImageJデフォルトのオートで2値化処理をかけ、ホールをうめる「Fill Holes」処理を行い、関心領域を黒塗り状態にして、そのままOutlines処理をかけ、最後に反転処理を行います。

2つ目の手順は、2値化する前に、エッジ検出を行い、それを反転させてから、スケルトン処理を行います。

いずれの操作でも、最後に関心領域外の余計な領域に残ったものをマニュアルで削除しなければなりません。もしくは、Snakeなどの自動辺縁抽出を利用するといいかもしれません。

これらの画像を作成したら、Analyze>Tools>Fractal box count...で、背景が黒塗りの場合はBlack backgroundにチェックをして実行しましょう。

今回は、2つ目の手順で試した結果をご紹介します。
比較対象として、円形ROIの辺縁抽出結果をフラクタル解析した結果と比べてみます。
やはり、構造が複雑な形状の方が、フラクタル次元が高い値を示していることがわかります。

今回の解析では、肺がんのスピキュラと呼ばれる部分の形状に合わせてフラクタル図形が計算され、フラクタル図形を含むボックス(ピクセル)がカウントされ、前述のアルゴリズムに従ってフラクタル次元が算出されています。このように、フラクタル次元は構造の辺縁の複雑さを客観的に示すことができます。もし解析対象(例えばスピキュラ)の構造の複雑さが疾病の悪性度と相関することが知られている場合は、構造の複雑さを客観的に示す指標であるフラクタル次元が活用できますね。
フラクタル次元は、応用がいろいろとできそうなイメージングバイオマーカの1つです。

(切り出した領域)
(関心領域の辺縁を抽出)



(フラクタル解析 D = 1.24)
(円形ROIの辺縁を抽出)

(フラクタル解析 D = 1.12)


次回は3D画像処理の元となる複数の画像群(スタック)の画像処理について触れていきます!

Reference
  • 山本修司:ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(222) 2007, p110-112

2015年9月24日木曜日

第24回 2値画像の形状解析で学ぶ!

医用画像や顕微鏡から得られた細胞画像より、注目する領域を抽出し、その特徴量を計算することはよく行われます。
この際、最も困難なことは、画像の中からターゲットとする特定の臓器、組織、病巣(腫瘍や特徴的な広がりを持つ炎症など)や細胞などをコンピュータによって自動抽出することです。この抽出方法(セグメンテーション)については、多くのアルゴリズムが論文などで紹介ますが、lmageJのセグメンテーションのためのプラグインも多く紹介されています。

lmageJでの実践を交えて説明していきます。

画像の形状特徴量の計測


画像の形状特徴量の計測方法として、ImageJの代表的な機能を紹介します。あらかじめ、画像の大きさなどが既知の場合、ImageJのAnalyze>Set Scaleで画素値と実測スケールが正しいかチェックしましょう。次に、2値画像の特徴量として計算可能なものを、Analyze>SetMeasurementsのパラメータの中から選んでチェックボックスをオンにします。図に、胚細胞のImageJサンプルデータ(Rat_Hippocampal_Neuron)の細胞形状の特徴量計測例を示します。




 (Set measurements)
(Measure:result)

このようなさまざまな指標によって、細胞の形状特徴量を測定した結果を定量的に観察し、することができます
以下、チェックした特徴量の一部を簡単に解説します。

①面積(Area)
画像の特徴量の面積は、画素数を数えることにより求められます。2値化した画像の対象(黒または白画素)の場合は、ヒストグラムから各濃淡値のピクセルカウントを確認できます。

②周囲長(Perimeter)
周囲長は、まず境界追跡アルゴリズム(Snakeなど)を用いて境界を求め、画素の長さを加算していきます。このとき、斜め方向に隣接する画素をルート2倍にするなど、誤差補正が必要です。
2値画像での2点間距離は、一般的な数学で使用されるユークリッド距離や、4近傍 (上下左右)処理によって2点間で接する画素を追っていく方法(4近傍距離または市街地距離と言う)、および8近傍処理によって2点間で最も接する面積が多い画素を追っていく方法(8近傍距離またはチェス盤距離と言う)で求められます。
ImageJでは、Bounding Rectangle (限界長方形:選択枠の最小長方形)という指標もあり、出力結果では、BX、BY、Width、Heightという値が結果テーブルに表示されます。BX、BYは、長方形の左上の角の座標を示しています。

③楕円指標(Fit Elipse)
ImageJでは,メニューのAnalyze>Set ScaleのFit Elipseに相当します。セグメンテーション後の境界枠に適した楕円を計算し、Major(長軸)、Minor(短軸)、Angleが指標として出力されます。MajorやMinorは、最も適する楕円の第1軸、第2軸で、Angleは画像のx軸と平行な線と、楕円径の第1軸と交わる角度を示します。
Set ScaleダイアログにあるPixel Aspect Ratioの値が1.0で有効となり、それ以外のときは使用できません。

④円形度(CircularityまたはRoundness)
対象がどれだけ円に近いかを表す尺度のことで、面積をS、周囲長をLとしたとき、4πS/L2乗で計算されます。対象が円である場合は1.0で、複雑になるにしたがって数値が下がります。

⑤2点間最大距離(Feret’s Diameter)
ImageJでは、“フェレーの直径(Feret’s Diameter)”という名称で実装されています。選択枠内の任意の2点間の最大距離を計測します。


(Pharmacopeial Forum:Volume No.30(6) Page 2212, fig 1.)

2値画像のモーメント特徴


一般的に、モーメントとは物体を回転させる力の大きさを表しますが、これを画像の統計量として捉えて、特徴量抽出に利用することもできます。
2値画像のモーメント特徴とは、画素の位置の重みづけをして合計した数値です。モーメントは、多値画像においても、ImageJを用いて画素値の平均値の周りの1次モーメントから4次モーメントまで特徴量を簡単に計測できます。
2値画像のモーメントは、ImageJでは標準では実装されていない特徴量ですが、重要な指標ですので簡単に説明します。

画像のモーメントは、一般的に次式で定義される(p+q)次のモーメントM(p,q)を用いて計算されます。




ここで、fijは画素値であり、対象の画素は1で、対象より外側のバックグランドは0です。
重心座標(I,J)を決定する場合、I = M(1,0) / M(0,0)、J = M(0, 1) / M(0,0)がよく用いられます。主軸の方向を示すtanθは、対象画像が伸びている方向を表し、次式によって計算されます。



ImageJのプラグインでは、ジョンズ・ ホプキンス大学でからMomentMacroJマクロとして紹介されているので、試してみてはいかがでしょうか。(http://www.hopkinsmedicine.org/fae/mmacro.html).

※テキストをダウンロードして、ImageJのプラグイン>マクロからインストールでダウンロードしたテキストファイルを選択します。これで、マクロのリストに加わるので、もう一度プラグインのマクロから選択して起動します。

(補足)
参考記事に別のプラグインの紹介があるので、こちらでもご紹介いたします。
ImageJを用いた拡散テンソル画像処理プラグイン
https://sites.duke.edu/dblab/jdti/#comment-34

このプラグインは教育や研究用に配布されています。
入手するには、上記ページからリクエスト(コメント)を送信して承諾してもらう必要があります。

テストデータも入手できますので、勉強には最適ですね。

さて、次回は特徴量解析の少し専門的な技術をご紹介します!

Reference
  • 「ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(22・1) 2007, p94-95」

2015年9月23日水曜日

第23回 2値画像処理の基本で学ぶ!

前回は2値化の手法について代表例を紹介しました。
画像を2値化することによって元画像のさまざまな特徴の抽出や特徴量の計算もしくは計測が可能になります。今回も引き続き、2値化についてふれつつ、画像処理方法や解析・計測するための基本的な方法について説明します。

2値化の連結性


2値画像は画素値を2つの値(例えば、1と0)で表現するため、同じ画素値の画素が群集をつくったりする様子などを観察でき、この状態を定義するために「連結」という概念が用いられます。

一般的に、2値画像上の任意の画素(i,j)に対して、その周辺を埋める(i+x, j+y) における適当な整数の対(xとy)によって配置される画素の集合を“近傍”と言います。

二次元のデジタル画像では、その任意の画素(i,j)の1画素分の上下左右を最近傍とした4近傍と、対角方向の画素も含める8近傍が使用されます。
次の図に示すとおり、任意画素(i,j)の4近傍は、画素(i+1,j)、(j,j+1)、(i-1,j)および (i, j-1)であり、8近傍では、これらに(i+1,j+1)、(i-1,j+1)、(i+1,j-1)、および(i-1, j-1)が加わります。

(参考記事より引用)

また、8近傍を上図のx0からx8のように反時計方向に回るような画素の位置で表記することもあります。
この2値画像中で互いに連結している画素の集合を1つのクラスにまとめると、2値画像の、例えば、いくつかの画素値0を持つクラスと画素値1を持つクラスができ、この各々のクラスを“連結成分”と言います。

連結成分の特徴量


2値画像の中で、画素値1の連結を「4連結」でクラス分けする場合と「8連結」でクラス分けする場合で連結性は異なります。画素値1について4連結で連結成分を抽出した場合、抽出されなかった画素値0の連結成分を観察すると、8連結で考えなければ矛盾が生じるような連結構成になっています。また、その逆で、画素値1について8連結で連結成分を抽出すれば、画素値0の連結成分は、4連結で説明しなければ矛盾が生じるような連結構成です。

画素値0の連結成分について、その周辺の画素が1つも連結しない場合は、これを“孔(ホール)”と言います。

画素値1の連結成分が孔を含まないときに、これを単連結成分、孔を1つでも含むとき、多重連結成分と言います。

次の図に、これらの説明図を示します。図中aは,4連結の場合で、成分数が5つ、孔はなし、図中b は8連結の場合で、成分数が3、孔の数が3です。また、成分数から孔の数を引いた数を“オイラー数(示性数)”と言い、2値画像の特徴量計算などに使用されます。
例えば。孔がない粒子は、オイラー数が成分数と等しいため、成分数の推定が可能となります。

(参考記事より引用)

連結数と特徴点について


2値画像において、ある画素の値を変更しても画像全体の連結性が保持されることを“消去可能”と言います。これを調べる際、“連結数”を用いると簡単に画素の消去可能性がわかります。

連結数は注目する画素の8近傍において、連続する画素を1つの成分とみなし、その成分数を数えることによって求められます。
この場合、4連結と8連結はいずれも0〜4の値をとります。次の図に、連結数とその連結数から分類される特徴点を示します。図のように、3×3マトリックスの中央の画素に注目し、その周辺の8近傍の画素が連結している成分数が連結数であり、連結数0が内部点、連結数1が端点、連結数2が連結点、連結数3が分岐点、連結数4が交差点となります。

(参考記事より引用)

収縮と膨張処理


注目画素の近傍を一回り剥ぎ取るような処理を“収縮(エロージョン)”と言います。この処理は、成分の形状によっては複数の成分に分割されるために、前述したような消去可能な処理とは異なります。また、逆に、注目画素に近傍画素を一回り加えるような処理を“膨張(ダイレージョン)”と言います。

これらの処理は数学的にはミンコフスキーの和と差として定義されています。

この処理は、孔がある場合,それを埋める効果があります。また、この処理も、非連結成分を連結して1つの成分に統合することがあるため、連結性を保持しない処理です。

ダイレーションやエロージョンは構造要素よりも小さな凹凸を画像から取り除く効果があります。しかし、実際にそれらが単独で用いられることは少ないです。これは、ダイレーションを行えば図形は拡大し、エロージョンを行うと縮小し、いずれの処理も処理の前後で画像のサイズが基本的に変わってしまうからです。

このため、両者を組み合わせ、大体の大きさが変わらないような処理が多くの場合に用いられます。それがオープニングとクロージングです。

同じ回数だけ膨張して収縮する処理を“クロージング”と言い、この組み合わせ処理によって小さな孔を除くことができます。またその逆で、同じ回数だけ収縮して膨張する処理を“オープニング”と言います。
どちらの処理でも、画像の小さな孤立成分を除くことができます。

以下、ImageJによる各処理結果を示します。
メニューのProcess/Binaryのプルダウンメニューに、2値化画像の各画像処理(Erode, Dilate, Open, Close)が並んでいます。

元画像

2値化
(図中の黄色枠を切り出す)

狭窄部のみ切り出し

エロージョン

ダイレージョン

オープニング

クロージング

収縮ではノイズ除去効果、膨張では穴埋め効果が得られているのが観察できます。実際には、オープニングとクロージングにおける収縮と膨張の回数は画像の質を考慮して設定しなければなりません。
オープニング処理は、外側に突き出した突起が削られて滑らかになり、図形の辺縁を滑らかにする一種の平滑化処理です。だだし、内側からの平滑化ですので、入江のような部分はそのまま残ります。また、狭い部分は図形の分離が起こります。孤立した小さな領域は消滅することから雑音除去の目的で利用することもできます。
クロージング処理はオープニング処理同様に平滑化された図形が得られます。入江のような部分や小さな穴がふさがれる効果があります。

次回も引き続き、2値画像処理について話を進めていきます。

参考文献:
  • 山本修司:ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(21・10) 2006, p75-77」
  • 田村秀行:2値画像の連結性と距離.2値画像処理,コンピュータ画像処理,東京,オーム社, 143一147,2002.
  • 小畑秀文:モルフォロジー,東京,コロナ社,p43,1996.


2015年9月21日月曜日

第22回 画像の2値化についてで学ぶ!

この記事は参考記事を援用して、筆者の考えも交えつつ、記述しています。

前回まで扱ってきた画像は、すべて8ビット〜32ビット階調などの複数の濃淡値を持っていました。今回は、画素値が0か1しかとらない2値画像化の手法(2値化処理)を説明します。


2値画像処理は、日常臨床で使用されるソフトウエアや三次元画像処理ワークステーションなど、いろいろな場面で適用される処理です。
CT画像から、体脂肪測定、骨や肺の空気などの領域抽出(セグメンテーション)、形状の特徴量の解析(フラクタル次元の計測、フーリエ記述子、カオスなど)などがその例です。

2値化の方法


画像の濃淡を0か1などの2値に分ける作業を”2値化処理”といいます。この作業のための境界値を閾値(threshold)と呼び、例えば、その値以上の画素値を1、それ未満の画素値を0へと変換することにより、2値化処理ができます。

(元画像)

二値化画像

(二値化設定:Image>Adjust>threshold)

閾値を決定するための最も単純な方法は、画像のヒストグラムで上限と下限の閾値を決める方法です。

p-タイル法


p-タイル法とは、パターン領域の面積比率(パターンデータ数/全データ数)を指定して閾値を決定する方法で、パターン領域の面積比率があらかじめ推定できなければなりません。

そのため、対象画像に対して事前の画素情報が必要です。医用画像処理では情報が未知の画像が多く、この方法の適用が困難な場合が多いです。図面や文書画像など、切り出すべき対象図形の面積が推定できる場合に有効な方法です。

(この例では黄色枠の左端の縦線位置の濃淡値が閾値になります。)


モード法


モード法は対象図形と背景の濃度分布の濃度差が大きい場合にヒストグラムが双峰型の分布となることを利用し、その山の部分と谷の部分が明確なときに、その谷底を閾値に決定する方法です。この方法はノイズの多い画像や複雑な画像ではそのヒストグラムから谷底を見つけるのが困難であり、あらかじめノイズ軽減処理などの前処理を加えた後に使用される場合が多いようです。


(参考記事より引用)


判別分析法


判別分析法とは、濃度ヒストグラムから2つのクラスを分割するために、その閾値を判別分析する方法です。

2つのクラスの平均値の分散(クラス間分散)と、各クラスの分散(分離度もしくは判別比)が最大になるような閾値を決定します。
ImageJでは、”Otsu Thresholding”として判別分析法を実装しています。


2値化処理による臨床応用例(CT画像からの脂肪面積測定)


簡単な臨床応用例として、ImageJによるCT画像を用いた体脂肪面積の測定の方法を紹介します。


方法は非常に簡単で、図に示すように、臍部の断層CT像から、脂肪領域の面積を計算するだけです。(輪郭抽出の際に用いる”領域抽出法”は次回以降で紹介予定ですが、ここでは、Spline Snakeプラグインや、Segmentation Assistantツールが便利です)

一通りやってみます。

最初に、面積、体積の算出にはピクセルサイズを使いますから、ピクセルサイズをImage>Propertiesで確認しておきます。

(X:約0.97mm, Y:約0.97mm, Z:2.0mm)
(DICOM画像は自動的にピクセルキャリブレーションされます)

まず、余計な領域であるバックグラウンドを削除します。ROIを体表に設定します。
設定したら、Edit>Clear Outsideします。

(ポリゴンROIを設定)

(Edit>Clear Outside)

もう1つROIを追加して内臓脂肪および皮下脂肪の領域を分けられるようにします。



(ROIの記録は"ROI Manager"を使います。ROIを微調整したら、Updateを忘れずに。)

次に、Threshold機能を使って、脂肪とそれ以外の画素値を2値化します。
今回はROIを2つ設定しています。全体に閾値処理をしたいので、体表のROIを選択した状態でThresholdします。大雑把ですが、筋骨格と脂肪が分離できます。



(Thresholdを設定してApply)

脂肪と筋骨格がどのピクセル値に割り当てられたか確認しておきます。

(キャプチャなのでマウスカーソルが見えていませんが、脂肪(黒い領域)がピクセル値「0」(ImageJ上のvalue=0)であることがわかります。それ以外(白い領域)は筋骨格でピクセル値255です。体表の白い部分はノイズとして捉えて下さい。)

ここで、それぞれのROIから、0のピクセル数と255のピクセル数をカウントします。
カウントはAnalyze>Histogramで確認します。
ヒストグラムを表示する際は、ヒストグラムにしたいROIをアクティブにしてから実行します。
ヒストグラムが表示できたら、カラーバーのところで対応するピクセル値にマウスを持っていき、ヒストグラムウィンドウ右下のカウントを確認します。
例えば、体表のROIをアクティブにしてヒストグラムを確認してみると、ピクセル値0:45725、ピクセル値255:28027であることがわかります。

(ピクセル値0:45725)

(ピクセル値255:28027)

念の為、これらのカウントを合計して確かめておきます。
体表ROI内側のピクセル数はCount:73752です。そして、画像から45725+28027=73752であることがわかりますから、きっちり2値化できていることが確認できます。

同様に、筋骨格を囲ったROIでもカウントを確認します。総数は52487、ピクセル値0:26314、ピクセル値255:26173です(画像省略)。



最後に、ピクセルカウント数から面積、体積を算出します。

筋骨格

26173 * 0.97mm * 0.97mm = 約24626mm^2(平方ミリメートル)
26173 * 0.97mm * 0.97mm * 2.0mm = 約49252mm^3(立方ミリメートル)

内臓脂肪

26314 * 0.97mm * 0.97mm = 約24758mm^2(平方ミリメートル)

皮下脂肪

(45725-26314) * 0.97mm * 0.97mm = 約18263mm^2(平方ミリメートル)


(ImageJは、DICOM画像を解析する場合、自動でDICOMヘッダーからピクセルサイズを調整します。)

次回も引き続き2値画像処理について述べます!


Visionary Imaging Services, Inc.


  • メディカルイメージングテクノロジーサポート
  • 臨床研究、臨床試験サポート
  • 医用画像処理解析ソフトウェア開発
  • コンタクト:お問い合わせはこちらまで

References

  • 参考記事:「山本修司:ImageJで学ぶ実践医用・バイオ画像処理.INNERVISION(21・10) 2006, p75-77