Mathematical morphology に基づくバイオイメージからの構造情報の抽出と解析
Extraction of Structural Features of Biological Images and Its Evaluationby Mathematical Morphology
木森義隆ａ，ｂ，ｃ，諸根信弘ｂ，片山栄作ｃYoshitaka Kimori, Nobuhiro Morone and Eisaku Katayamaａ（財）医療機器センターｂ

国立精神・神経センター神経研究所ｃ

東京大学医科学研究所

要旨

Mathematical morphology は， 画 像 中 の 物 体 形 状 を 操 作 し 計 測 す る 理 論 的 枠 組 み で あ る． 本 稿 で は， バ イ オ イ メ ー ジ に 対 し，mathematical morphology に基づく画像処理を適用することにより，その構造情報を抽出し，定量的に計測する手法について述べる．

キーワード：画像処理，mathematical morphology，電子顕微鏡

処理や，周波数選択型の処理にはないユニークな特性である．
1.

はじめにまた，全ての演算は，Minkowsky 和と差という単純な 2 項

光学または電子顕微鏡を用いて in vitro 条件，あるいは組演算を基本とするため，計算機での実装が容易である．バイ

織・細胞中における生体高分子の分布や構造を捉え，それらオイメージングの分野においてモルフォロジはまだ公知ではの機能や動態を探るバイオイメージング技術は，現代のライないが，光学顕微鏡像3），電子顕微鏡像2） さらに DNA マイ

フサイエンスにとって必須の手段である．しかし，画像デークロアレイの画像4）などに対して適用した例がある．いずれ

タの解釈には，観察者の主観が介在し，客観的とは言い難いも，対象物のセグメンテーションや特徴的な構造の抽出に用面もある．これは，計算機による画像処理・認識技術に未整いられ，効果的な処理方法が実現されている．

備な部分が多いことに起因する．観察者が像を解釈する場合，本稿では，様々なバイオイメージから生物学的に有意な情経験や知識に基づいて“かたち”を読み解くプロセスを経る．

報を発掘・整理して観察者に提示し，新たな知見につなげる

しかし，従来の画像処理手法ではそれを定量的に取り扱う手ための一連のワークフローを計算機上に構築するモルフォロ立てがなく，また，互いに関連性のない様々なフィルタを寄ジカル画像処理理論について概説する．

せ集めて実行することが多く，解析に対する論理性があいま2.

いであった．
1）

Minkowski 集合演算とモルフォロジの基本演算

これを踏まえ，我々は，Mathematical morphology （数理

モルフォロジは，1960 年代の後半，鉱石の顕微鏡画像中形態学，以降モルフォロジ）に基づく非線形的な画像解析手に分布する鉱物の幾何学的特性と物理的特性の関連を解析す法を開発し，画像の定量的，客観的な評価を目指して従来とるための手段として G. Matheron と J. Serra により考案され2）は異なる観点から新たな試みを行った ．

た．理論の整備と計算機環境の発展に伴い，1980 年代になっモルフォロジは，集合論を基盤に，様々な画像処理を一貫て，テクスチャ解析や様々な画像処理に適用されるようにした理論体系のもとで表現し実行する数理体系である．物体なった．モルフォロジは，対象画像に構造要素と呼ばれる基の“丸い”，“線維状”などのかたちの表現形式を直接取り扱底関数を作用させる集合論的操作により，画像中の構造を単うことができるため，形態に対し一定の論理的根拠を備えた純な“かたち”の集合として処理できることに特徴がある．

定量的な解釈が可能となる．これは，輝度分布に基づく画像

これを活かし，画像中の形状に関わる特徴抽出や識別に適しｂ〒 187–8502 東京都小平市小川東町 4–1–1
TEL: 042–341–2712; FAX: 042–346–1749
E-mail: kimori@ncnp.go.jp
2008 年 11 月 17 日受付

講座

た手法として用いられる．モルフォロジは一般に N 次元空間における集合論として展開されるが，本稿では，バイオイメージへの適用を鑑み，2 次元モルフォロジ演算について説明する．

Mathematical morphology に基づくバイオイメージからの構造情報の抽出と解析

41【著作権者：社団法人 日本顕微鏡学会】

2.1.

2 値のモルフォロジ演算

ここでは，まず，2 値のモルフォロジについて説明する．
2 値画像は，黒い画素（値 0）の背景中に，白い画素（値 1）が集まってできた物体で構成されるものとする．この物体の存在領域が関心領域である．このとき，物体は「白い画素の集合」と捉えられる．個々の白い画素は集合の要素（元）であり，画像の原点からのベクトルとして表現されるので，物体の形状をベクトルの集合として捉えなおすことができる．
本稿では，集合（画像）は，アルファベットの大文字で，要素は小文字で表す．
モルフォロジの演算は以下に示す Minkowski 和と差が基
本 に な り， そ れ ら の 組 み 合 わ せ で 演 算 が 構 成 さ れ る．
Minkowski 和とは，2 つの集合 X，B があるとき，それぞれの集合の要素 x，b の全ての組み合わせからできる集合のことである．集合 X を b だけ平行移動して得られる集合を (X(X)bと書くと，Minkowski 和はその和集合であり，以下のように

図 1 原画像 X の構造要素 B による dilation, erosion, openingおよび closing 演算．点線は原画像の輪郭を，灰色の領域は演算後の画像を表す．

定義できる．
(X)bX⊕B = ∪ (Xb∈B（1）

一方，Minkowski 差は，集合 ((XX)b の積集合として以下のように与えられる．
(X)bX⊖B = ∩(Xb∈B( ⊖Bs)⊕B
XB = (X（6）

一方，closing を XB とすると，( ⊕Bs)⊖B
XB = (X（2）

ここで，記号⊕，⊖はそれぞれ，Mikowski 和，差の演算子（7）と定義される．これらは，dilation と erosion の演算順序が異なる．
図 1 に opening と closing の例を示す．opening は境界を滑を意味する．
また，集合 B を原点を中心として 180° 回転させて得られらかにし，円盤 B が入り込めないような地峡形状や尖ったる集合を，集合 B の対称集合と呼ぶ．ここでは，Bs と表し，岬形状を除去する平滑化作用をもつ．一方 closing は，B が以下のように定義する．

入り込めないような間隙や穴を埋め，一方で突起領域をそのBs = {−b:b ∈ B}

（3）

まま残す作用をもつ．
2.2.

多値のモルフォロジ演算

集合 B をその対称集合 Bs で置き換えたときの Minkowski

多値（濃淡）画像は，画素の座標（画像の空間領域）を x,和と差をそれぞれ，dilation, erosion と呼ぶ．これらはモルフォ

y の 2 つの次元で表し，画素値 f(, ) を次元 z で表す 3 次元
f(x,yロジの基本演算であり，dilation は，“ずらし重ね”，erosionEuclid 空間上の集合と考える．多値画像に対するモルフォロは“掻き取り”と訳される．dilation は次のように定義される．
(X)− bX⊕Bs = ∪(Xb∈B（4）b∈Bを満たす領域であり，これは 2 値集合として表現することができるため，濃淡画像が 2 値の集合へ変換される．

一方，erosion は次のように定義される．
(X)− bX⊖Bs = ∩(X

ジ演算は，陰影（Umbra）という概念を用いれば，2 値画像,)の場合の演算と平行な議論が可能となる．陰影とは，z ≤ f(f(x,y（5）これらの演算の適用例を図 1 に示す．2 値画像を表す閉集関数 f の関数 g による Minkowski 和は，次のようになる．
((f⊕g)(x( ) = max {f{f(x − z) + g(z( )}x−z− ∈F,z∈G（8）合 X に，2 次元の円盤状の閉集合 B（構造要素）を作用させここで，F,F G はそれぞれ，関数 f,f g の定義域である．Minkowskiる．この場合，B は原点に対し対称であるため，dilation と差は，次のようになる．

erosion の結果は，それぞれ，Minkowski 和と差の結果に等しくなる．X に B をずらし重ねた結果，dilation は B の半径分だけ外側に膨張させる．一方，erosion は B の半径の大きさ だ け X を 掻 き 取 り， そ の 分 だ け 収 縮 さ せ る． ま た，
dilation と erosion を 組 み 合 わ せ た 演 算 と し て opening とclosing がある．集合 X の構造要素 B による opening を XB と表すと，以下のように定義できる．
42((f⊖g)(x( ) = min {f{f(x − z) − g(z( )}x−z− ∈F,z∈G（9）関数 g(x( ) の原点に対する対称関数 gs(x( ) を定義すると以下のようになる．
( ) = g(− x)gs(x（10）2 値のモルフォロジ演算と同様，この対称関数を用いた
顕微鏡 Vol. 44, No. 1（2009）【著作権者：社団法人 日本顕微鏡学会】

Minkowski 和と差をそれぞれ dilation と erosion と呼び，多値画像に対する opening, closing もこれらを用いて定義される．
Opening: fg = (f( ⊖ggs)⊕gg（11）( ⊕ggs)⊖gg
f g = (f（12）Closing:

ロジ演算を用いた平滑化フィルタを提案している6）．その具体例として opening と closing を組み合わせた，open-closing（OC）演算，close-opening（CO）演算があり，これらの演算 を さ ら に 線 形 結 合 し た LOCO（Linear combination of OCand CO）フィルタ7）が知られる．しかし，これらのフィルタには，処理結果に前述した構造要素の形状由来のアーチファクトが生じる場合がある．それを改善するため，OC，CO 演

3. Image rotation mathematical morphology に基づく新規
フィルタ

算を式（13，14）で書き直し（それぞれ hOhC, hChO とする），それらを線形結合した IRMM による新たな平滑化フィルタ

生物試料の形態は微細かつ複雑な様相を呈し，画像取得の段階で不可逆的なノイズが重畳される．それを軽減し原画像に含まれる情報を有効に利用するためには，微細構造の特性（SF）を定義する．
F( ChO)/2SF = ((hOhC)/2 + (h（15）を良好に保存しつつ，背景ノイズや，画素値の変調にロバス

図 2 はノイズ除去を目的として，原画像（図 2a）に対し，トに対応できる画像処理法が望ましい．しかし，前章で示しメディアンフィルタ，LOCO フィルタおよび SF による平滑た通常のモルフォロジ演算では構造要素の作用方向が一定の化処理の結果を比較した例である．図 2b は，メディアンフィため，対象画像中に含まれる様々な向きの構造物に対し，均ルタでの結果である．ウィンドウサイズは，9 × 9 画素とした．

質な処理が困難である2）．そのまま実施すれば，解析対象のノイズの除去とともに，解析対象の輪郭がなまっている．図

変形，構造要素の形状に依存したアーチファクトなどが生じ

2c は，LOCO フィルタでノイズを除去した結果である．やる場合がある．

はり 9 × 9 画素の正方形構造要素を用いたが，この場合も輪そこで我々は，モルフォロジの特性を活かしつつバイオイ郭形状が変形している．さらに，構造要素の形状が重畳されメージに対応できる画像処理を目指して新たなモルフォロジアーチファクトが生じた．対象物は，6 回の回転対称性をも演 算（Image Rotation Mathematical Morphology（IRMM） 演

つため，この輪郭形状を保持するには，それぞれの突起の方2,5）算）を開発した．これは対象画像をあらゆる方向に回転さ向に沿った平滑化処理が必要である．SF による処理では，せて処理することにより，任意の方向を向いた構造の処理を9 × 9 画素の正方形構造要素を用い，原画像を 22.5 度ずつ回可能にするものである．

転させ，そのつど平滑化処理したものを統合し出力した

いま，対象画像を f，構造要素をfg とする．f をその中心座（図 2d）．その結果，ノイズは除去され，対象物の輪郭形状

標を回転の中心として角度 180·i/i/N 度（i = 0, 1, …, N−N 1）ごが保たれている．図 2c と比較すると，構造要素の処理方向とに回転させ，そのつど演算を行う．N は 180 度を等分したの多方向性が結果の向上をもたらしている．

数である．結果は，各角度ごとの処理画像 N 枚のセットを統合して出力する．
この opening および closing 演算は以下の最大（max），最小（min）演算として定義される．
Opening: hO(x( , y) = max[h[ 0(x( , y), h1(x( , y),…, hN−( , y)]（13）N 1(xClosing:N 1( , y) = min[h[ 0(x( , y), h1(x( , y),…, hN−( , y)]（14）(xhC(x

ここで，hi は画像 fi の g による opening，hi は closing を意味する．
3.1.

IRMM による平滑化フィルタ

画像処理において平滑化は重要な役割を担う．単に画像のノイズ除去の目的のみではなく，背景と対象物体の分割（segmentation）や特徴パターン抽出の前処理として使われる．平滑化フィルタには，対象物の領域が十分に均質な濃度値で表現されるとともに，そのエッジや比較的小さな構造的な特徴を保存することが要求される．これまで，線形や非線形の様々な平滑化フィルタが提案されているが，この要求を満たすロバストな処理は困難であった．
Maragos らは，従来のフィルタでは実現できない，幾何学的な解析を行うことが大きな優位性をもつとして，モルフォ講座

図 2 様々な平滑化処理の比較．（a）原画像（311 × 311 画素），（b）メディアンフィルタ，（c）LOCO フィルタ，（d）IRMM による平滑化フィルタ（SF）の結果．
F

Mathematical morphology に基づくバイオイメージからの構造情報の抽出と解析

43【著作権者：社団法人 日本顕微鏡学会】

3.2.

IRMM による特徴抽出フィルタでも孤立した輝点として観察できるが，その分子の位置の決さらに，top-hat 変換と呼ばれる演算も IRMM で定義した．
TH(H(f) = f − hO

（16）

top-hat 変換は，原画像から，それを opening した結果を差定には個々の輝点の像を識別して抽出する必要がある．ここでは，輝点の抽出過程を自動化するための手法を説明し，抽出した輝点を既存の解析法に適用可能なかたちに変換する手法を述べる．

し引く演算である．opening（hO）は，一種の平滑化処理で実際の画像に先立ち，モデル画像に本手法を適用した例をあり，構造要素に比べて小さな凸部は削られる．これは，示す．図 4a は，不均一な背景画像（231 × 231 画素）上に輝low-pass filter であり，この出力を元の画像から引き去れば点のモデル像（半径 7 画素のガウス分布）を配置したもので

high-pass filter（opening によって削られた部分が残る）となある．図中の点線は，画素値のプロファイルを計測した位置る．構造要素のスケールを変化させるとその通過帯域が変化であり，以降の各処理による画素値の変化を図 4e に示す．

するという特性を活かし，画像に含まれる特定の形状の抽出

図 4b は，a にガウスノイズ（σ = 40）を加えコントラストをに用いられる．

低下させたものである．原画像との S/N 比は，12.36［dB］

図 3 は，様々な大きさの構造要素で TH による処理を行っであり，画素値のプロファイルを見ても，輝点と背景はほとた例である．原画像（図 3a）は，305 × 305 画素の画像であり，んど区別が付かない．この画像に TH を用いた結果を図 4c

画素値を高さとした複雑な地形を表現している．図 3b では，に示す．ただし，抽出された領域を 1，他を 0 とする 2 値化

構造要素を 5 × 5 画素の正方形とすることにより，5 画素以下を行っている．構造要素は，15 × 1 画素の線分を用いた．の幅をもつ凸状の領域を抽出した．図 3c では，25 画素以下，top-hat 変換では 1 階調でも差があればその領域を抽出でき

図 3d では，45 画素以下の幅をもつ凸状の領域を抽出した．

る．抽出したそれぞれの輝点領域の重心位置を決定し，それ構造要素をパラメータとして TH を用いれば，このようにを頂点に半径 7 画素のガウス分布を発生させて輝点を復元

原画像から特定の形状情報を選択的に抽出することが可能で（規格化）した（図 4d）．頂点位置の誤差は，画像劣化の影ある．実際の処理では，顕微鏡などの分解能を考慮しながら響により 0.88 ± 0.06（平均値±標準偏差）画素であった．

抽出すべき領域の大きさを実測し，それに基づいて構造要素

従来の輝点解析手法は，輝点の点像分布関数（point spreadの大きさを決定する．

function）をガウス関数と仮定したテンプレートマッチングを基本とする手法が多い8）．従って，この規格化および頂点

4. 光学顕微鏡で可視化された粒子の自動解析を可能にする画像処理

位置の定量化は，既存のシステムに適用する際に重要になる．
次に，これを踏まえ培養細胞中のカベオリンに蛍光抗体を細胞内部に存在する分子に蛍光タグを融合して蛍光観察す

付け，光学顕微鏡で可視化した像を対象として処理したることにより，その局在や挙動を可視化できる．蛍光 1 分子（ 図 5）． モ デ ル 画 像 と 同 様 に 輝 点 を 抽 出， 規 格 化 し た．
図 5a，b の左側は，それぞれ別の細胞の原画像（515 × 515画素）を，そこから輝点を抽出した結果をその右側に示す．
輝点は，非常に込み合って存在しているため，大きな幅を持つ構造要素では，近接する複数の輝点をそれぞれ独立したものとして抽出できない．そこで，構造要素は幅 1 画素の線状

図 3 IRMM による特徴抽出フィルタ（TH）（a）
H による構造抽出．
原画像（305 × 305 画素）から，5 画素以下（b），25 画素以下（c），45 画素以下（d）の凸状の構造を抽出した．
44

図 4 輝点モデルの抽出．
（a）原画像（231 × 231 画素）．画素
値のプロファイルを測定した位置を点線で示す．（b）劣化画像．
（c）（b）に IRMM による特徴抽出フィルタ（TH）を施し，輝H点領域を抽出した後，2 値化した．
（d）輝点の復元（規格化）．
（e）
各処理後の画素値プロファイル．

顕微鏡 Vol. 44, No. 1（2009）【著作権者：社団法人 日本顕微鏡学会】

図 5 光学顕微鏡画像中の輝点の抽出．（a），（b）の左側は，カベオリンを可視化した蛍光顕微鏡像（原画像，515 × 515 画素）．
それぞれ，右は抽出された輝点の像．
（c）は，（b）の左図の矩形領域の拡大．
（d）は，（b）の右図の矩形領域の拡大．それぞれの右側に，画素値を高さに変換した図を表す．とした．実際には，あらかじめ輝点のサイズを実測し，7 × 1のミスなどが低減することが期待でき，より効果的であると画素とした．これにより，直径が 7 画素以下の輝点を抽出で考えられる．

き，また近接している場合は，それらが 1 画素でも離れていれば，それぞれを独立したものとして抽出できる．最終的に輝点は，半径 3 画素のガウス分布として規格化した．図 5b

5. アクトミオシン系の解析：蛋白質複合体の電子顕微鏡画像への応用の矩形領域を拡大したものが図 5c，d である．図 5c が原画像，筆者らは，精製された蛋白質の電子顕微鏡像の中から，対図 5d が抽出された輝点の画像である．それぞれの右側に輝象蛋白質像をセグメンテーションし，あらかじめ作成された度を高さとして示す．背景の不均一な輝度分布の中から，輝“像の見え方”のデータベース（対象蛋白質の結晶構造から点のみが抽出されていることがわかる．既存の輝点トラッキ構築された 3D モデルをあらゆる角度で 2D に投影した画像

ングシステムなどへの入力に図 5d を用いれば，マッチングの辞書）とマッチングを行うという画像認識システムを構築

図 6 蛋白質複合体への応用．
（a）原画像．337 × 337 画素．背景を含むアクチンとミオシンの複合体の電子顕微鏡像．
（b）
平滑化処理．（c）2 値化処理．（d）複合体の領域の抽出．（e）複合体領域のセグメンテーション．（f）ミオシンの領域の抽出．（g）
構造特徴（サブドメイン）を強調したミオシン像（左）と対応するモデル像（右）．主要なサブドメインを図中に示す．それ
ぞれ，U50（L50），upper（lower）50-kD subdomain，N は，NH2-terminal barrel，ELC（RLC）は essential light chain（regulatorylight chain）．

講座

Mathematical morphology に基づくバイオイメージからの構造情報の抽出と解析

45【著作権者：社団法人 日本顕微鏡学会】

した2）．ここでは，機能中の分子モーター蛋白質複合体（ア能であることがわかる．これは，周波数選択型の線形フィル

クチン，ミオシン）を急速凍結ディープエッチレプリカ法でタでは実現困難な処理であり，また，このような汎用的特徴捉えた画像を対象に解析を行った例を示す．まず，原画像は，ロバストなフィルタリングアルゴリズムの創出を支えて（図 6a，337 × 337 画素）から，アクトミオシンの複合体をいる．

抽出し，次に，複合体からのミオシン分子像のみをセグメンテーションする．このため，まず，原画像を平滑化処理し6.

（図 6b，SF による．構造要素は直径 21 画素の円盤），その後，本稿では，モルフォロジの基礎演算から，バイオイメージ

9）

自動閾値処理 により 2 値化した（図 6c）．

おわりにに適すよう拡張した IRMM 演算について説明した．そして，その後，最大の面積をもつ領域（すなわち複合体の領域）光学顕微鏡の輝点抽出，電子顕微鏡の蛋白質複合体に適用しの み を 残 し， あ と は 削 除 し た（ 図 6d）． そ の 画 像 にた例を示した．数学的基盤に裏打ちされた本画像処理理論でWatershed 変換10）を施すことにより，複合体の領域を複数のは，フィルタの設計，処理結果の解釈などに対して客観的な断片にセグメンテーションした（図 6e）．各セグメントは，議論が可能となる．解析対象のかたちを定量的に取り扱い，．この中それぞれ異なる濃淡値で表している（ラベル処理）

ノイズに埋もれた情報を顕に提示する本演算手法は，今後まから，ミオシンの領域を選択し，その部分をマスクとして，すます必要になると思われる．

原画像からミオシン像を抽出した（図 6f）．この画像を，ミオシン分子のさまざまな状態の結晶構造と比較し，最もよく合致する 3D 構造とその投影角の推定を行う．ただ，今回解析したのはアクチンの滑り運動を起こしつつある反応中間体の一瞬の構造であり，その結晶構造は得られていない．そこ
で，既存の結晶構造をサブドメイン単位に分割し，それらの構成部品を，電子顕微鏡像にマッチングさせ，全体として 3次元的に合致するモデルを再構築した（未発表データ）．モ
デル像との照合には，ミオシンのレプリカ像そのものは用いず，分子内サブドメインで構成される構造パターン抽出した
特徴画像を用いた（図 6g 左側）．構造パターンは，レプリカ像の空間分解能を考慮し，1 辺が 2 nm の正方形を構造要素として TH によって抽出した．図 6g の右側は，最もよく合致する原子モデルの像である．原子モデル像との対応により，ミオシンのサブドメイン配置の同定まで可能となった．
構造パターンを抽出した TH は，4 章で輝点を抽出したフィルタと同じである．様々な大きさや形状の異なる抽出対象物
でも，構造要素を変更することにより同じ手続きで抽出が可

46

文献1）Serra, J.: Image Analysis and Mathematical Morphology, Academic
Press, London,（1982）
2）Kimori, Y., Oguchi, Y., Ichise, N, et al.: Ultramicroscopy, 107，25–39（2007）3）Beare, R.: IEEE Trans. Pattern Analysis and Machine Intelligence,
28，1063–1074（2006）
4）Angulo, J. and Serra, J.: Bioinformatics, 19，553–562（2003）
5）「回転式モルフォロジ―フィルター及びそれを用いた画像処理法」出願番号：2006-280441
6）Maragos, P. and Schafer, R.W.: IEEE Trans. Acoust. Speech Signal
Process., 35，1153–1184（1987）
7）Schulze, M.A. and Pearce, J.A.: Proc. of IEEE. Int. Conf. on
Acoustics Speech and Signal Processing, 1993, Vol. 5, p. 57
8）Cheezum, M.K., Walker, W.F. and Guilford, W.H.: Biophys. J., 81，2378–2388（2001）9）Otsu, N.: IEEE Trans. Syst. Man Cybern., 9，62–66（1979）
10）Mayer, F. and Vachier, C.: Math. Morphology, proc. 6th Int. Symp.,
2002, p. 69

顕微鏡 Vol. 44, No. 1（2009）【著作権者：社団法人 日本顕微鏡学会】