ITKのWatershedフィルターによる領域分割(2)

前回は簡単にWatershed法の原理とITKの2種類のWatershedフィルターについてお話した。今回はそのフィルターのひとつであるitk::WatershedImageFilterの使い方について説明したい。

文章を読みやすくするために、以下より”itk::WatershedImageFilter”を ”WatershedImageFilter” と記述していく。

WatershedImageFilterの特徴は、以下の５項目にまとめることができる。

(1) ある地形に雨が降ったとき、等高線（標高）の高いところから低いところに雨滴が滴り落ち集水域（Watershed）を形成するという地球物理学的なモデルを基にしたアルゴリズムが実装されている。

(2) 入力画像の画素値（輝度値、濃度値）を高さとしたとき、位置\({\bf x}\)の値は高さ関数\(h({\bf x})\)で表すことができ、最急降下法によって\(h({\bf x})\)の局所的な最小値に移動する\({\bf x}\)の集合を集水域と定義する。

(3) 初期セグメンテーション領域をノードとしてマージツリーを作成する。

(4) 初期セグメンテーション領域は指定したフラッドレベル（Flood level）までマージされる。

(5) 無視できる局所的最小値のしきい値を指定することができる。

(1)~(3)については前回までにお話したので、今回は(4)と(5)について説明し、WatershedImageFilterを実際に画像に適用するためのソースコードを紹介する。

フラッドレベルを設定することで、初期セグメンテーション領域をどの程度マージ（統合）するかを決定することができる。フラッドレベル値の設定はWatershedImageFilter::SetLevel関数によって行う。パラメーターは入力画像の最大最小値間における任意のパーセンテージであり、0を設定した場合はマージされず初期セグメンテーションの状態のままとなる。1.0であれば全ての画素は1つの領域に統合される。

しきい値（Threshold）

しきい値を設定することによって、この値より小さい最小値を無視することができる。しきい値はWatershedImageFilter::SetThreshold関数を用いて設定することができる。このとき、フラッドレベルと同様にパーセンテージを指定する必要がある。

セグメンテーション処理の概要

ここからは実際にWatershedImageFilterを適用する方法について説明する。画像を入力しセグメンテーション画像を得るための処理のステップ（パイプライン）は以下となる。

なぜStep2での勾配強度画像（勾配の大きさを画素の値とする画像）を求めるかというと、エッジにあたる、他の画素値に比べ相対的に高い値を持つ画素がオブジェクトの境界と一致する高さ画像を作成するためである。ここでのオブジェクトとは、抽出したい対象物のこと。つまり勾配強度画像にWatershed法を適用すると、オブジェクトの境界が理論的には流水域の境界(Watershed line)と一致するのである。（実際には前述したフラッドレベルやしきい値による調整が必要であるが。）Step3で今回のテーマであるWatershedフィルターWatershedImageFilterを使用する。

プログラムの実装

それでは、実際にプログラムを実装し画像を出力していきたいと思う。まずは、前述のStep1、Step2の部分についてお見せする。メモリの初期化や変数の宣言は一部省略しており、このままコピー＆ペーストしてもエラーとなるのでお気を付け頂きたい。

元画像はRGB画像だが、グレースケール化した画像を入力画像(inputImage)とする。入力画像からの勾配強度画像を得るために、平滑化も行ってくれるGradientMagnitudeRecursiveGaussianImageFilterを用いる。

// 勾配フィルター処理
// 勾配強度画像作成
typedef itk::GradientMagnitudeRecursiveGaussianImageFilter  FloatGradientFilterType;
FloatGradientFilterType::Pointer floatGradientMagnitudeFilter = FloatGradientFilterType::New();
floatGradientMagnitudeFilter->SetInput(importImage->GetOutput());
floatGradientMagnitudeFilter->SetSigma(2.0);
floatGradientMagnitudeFilter->Update();

次にStep3とStep4について説明する。Watershedフィルターの入力は勾配強度画像であり、出力はラベリング画像となる。ラベリング画像での各領域の画素は0や1、2など、領域を表すためのラベル値を持ち、同じ領域に属する画素はすべて同じラベル値を持つ。このラベリング画像は視認性がよくないので、カラーマッピングを適用する。以下のコードがStep3とStep4の実装部分であり、前述したフラッドレベルは0.5（watershedImageFilter->SetLevel(0.5)）、しきい値は0.005（watershedImageFilter->SetThreshold(0.005)）としている。

// Watershedフィルター処理
typedef itk::WatershedImageFilter WatershedImageFilterType;
WatershedImageFilterType::Pointer watershedImageFilter = WatershedImageFilterType::New();
watershedImageFilter->SetInput(floatGradientMagnitudeFilter->GetOutput());
watershedImageFilter->SetThreshold(0.005);
watershedImageFilter->SetLevel(0.5);
watershedImageFilter->Update();

// カラーマッピング処理
typedef itk::Functor::ScalarToRGBPixelFunctor< unsigned long > ColorMapFunctorType;
typedef itk::RGBPixel RGBPixelType;
typedef itk::Image RGBImageType;
typedef itk::UnaryFunctorImageFilter UnaryFunctorImageFilterType;
UnaryFunctorImageFilterType::Pointer colorMapImageFilter = UnaryFunctorImageFilterType::New();
colorMapImageFilter->SetInput(watershedImageFilter->GetOutput());
colorMapImageFilter->Update();

ちなみに、フラッドレベル値やしきい値を0に設定した場合（watershedImageFilter->SetLevel(0) 、 watershedImageFilter->SetThreshold(0)）は初期セグメンテーション領域が出力される。一般的に初期セグメンテーション領域はオーバーセグメント（過分割）となり、オブジェクトの領域を正確に得ることはできない。つまりWatershedImageFilterにおいて、ユーザーによるフラッドレベル値としきい値の設定は非常に重要であるといえる。

WatershedImagFilterの問題点

冒頭でWatersheImageFilterの特徴として、

(1)　ある地形に雨が降ったとき、等高線（標高）の高いところから低いところに雨滴が滴り落ち集水域（Watershed）を形成するという地球物理学的なモデルを基にしたアルゴリズムが実装されている。

と説明したが、雨滴が滴り落ちず、その場に留まるフラットゾーンはどの集水域に所属するのか、という問題が出てくる。WatershedImageFilterでは、フラットゾーンに隣接する最も低い最小値をもつ集水域に分類される。つまりフラットゾーン上の画素はその集水域と同じラベル値をを持つことになる。しかしながらこの方法で正確にオブジェクト領域がセグメンテーションされるとは限らない。

まとめ

WaterhsedImageFilterでは、雨滴が高いところから低いところへと滴り落ち集水域を形成するという地球物理学的モデルがそのままアルゴリズム化されている。しかしながら、初期セグメンテーション領域の統合が高コストであることや、フラットゾーンの分類問題などが発生してしまう。

次回は、効率化されたアルゴリズムを持ち、さらに集水域の境界線を得ることができるMorphologicalWatershedImageFilterについて紹介したい。

参考文献とリンク：

WatershedImageFilter

GradientMagnitudeRecursiveGaussianImageFilter

“The watershed transform in ITK – discussion and new developments”

バックナンバー：

第一回：ITKのWatershedフィルターによる領域分割(1)

ITKのWatershedフィルターによる領域分割(2)

目次