前回（第三回）はリージョングローイング法を使って「人体の外側の人体組織ではない部分」を抽出し、「肺野内の空気領域」を得る方法について説明しました。またリージョングローイングはITKを使って実装しました。おなじみの第一回のおさらいとなりますが、次の方法で肺のセグメンテーションを行います。
(1)胸部CT画像入力
(2)K-means法で人体組織とそうでない部分を分ける
(3)人体組織でない部分（外側）をリージョングローイング法で抽出する
(4)肺野内の空気領域を得る
(5)クロージングで肺内部の空洞を埋める
(6)セグメンテーション完了！
前回は(3)と(4)について説明しました。今回は「(5)クロージングで肺内部の空洞を埋める」処理について説明していきたいと思います。前回と同様、この処理でもITKを使用します！
このクロージングとは、簡単に説明すると、アルゴリズムを使って対象物体内部にある穴を埋めていく処理のことです。以下で詳しく説明していきます。図1に示すように、肺内部にはたくさん空洞がありますが、これらのほとんどが血管です（肺がんであることもあります）。前の処理では空気領域を検出することが目的であったので、これらは現時点では検出されていないのです（白い部分が検出されている領域です）。肺野領域をすべて検出するには、血管も検出される必要があり、今回の目的がまさにそれなのです。では、参りましょう！
穴を埋めるだけクァ。
目次
クロージングとはClosingのことです。クロージングはモルフォロジー演算の一種です。画像処理におけるモルフォロジー演算（変換）では、カーネルと呼ばれるフィルタを画像上で走査し、対象領域において「膨張」や「収縮」とよばれる演算などを行います。膨張を行った後に収縮を行う処理は「クロージング」と呼ばれています。文章で説明するとわかりにくにので、まずは図で説明したいと思います。
因みに、英語で収縮はErosion、膨張はDilationというクァ。
図2において、カーネルを元画像で走査（スライド）させていき、カーネル内の領域に含まれる画素の値が全て1であれば1を出力し、そうでなければ0を出力します。したがって、収縮を行った後の画像は浸食されたような形状となります。
膨張処理ではカーネル内の領域に値が1である画素がひとつでも存在すれば、1を出力します。したがって、図3のように元の物体に比べると膨張しているのがわかりますね。
クロージングはこの膨張の後に収縮を行う処理のことで、穴を埋めるために膨張処理をして、元の外縁を得るために収縮処理を行うようなイメージです。
それでは、ITKを使ってクロージング処理を実装していきましょう。ITKは一見難しそうに見えますが、慣れてしまえば非常にわかりやすく使いやすいシロモノです。まずはITKではおなじみ、変数型を決定して実体化します。
// イメージインポートフィルター typedef itk::ImportImageFilter< short, 3 > ImportImageFilterType; // クロージング処理構造体 typedef itk::BinaryBallStructuringElement< short, 3 > StructuringElementType; // 収縮フィルター typedef itk::BinaryErodeImageFilter< ImageType3D, ImageType3D, StructuringElementType > ErodeFilterType; // 膨張フィルター typedef itk::BinaryDilateImageFilter< ImageType3D, ImageType3D, StructuringElementType > DilateFilterType;
BinaryBallStructuringElementとはカーネルです。short型の3次元としているので、ボールのようなカーネルですね。また、ImportImageFilterはITKにおける画像データテンプレートクラスと考えて下さい。ImageType3Dについては、前回（第四回）説明していますので、そちらをご覧ください。
次に、あくまでもvtkImageDataを主たる画像データクラスとしているので、前提として以下の変換を行います。
vtkImageData ⇔ ImportImageFilter
この変換を行うために、以下の処理を行います。かなり地道な処理です。また、メモリの初期化や変数の宣言は一部省略しており、このままコピペしてもエラーとなりますのでお気を付けください。
// 画像サイズの取得、設定 // inImg, outImgはVTKにおける画像データ // vtkImageData* inImg, outImg inImg->GetDimensions( dims ); ImportImageFilterType::SizeType size; size[0] = dims[0]; // xサイズ size[1] = dims[1]; // yサイズ size[2] = dims[2]; // zサイズ ImportImageFilterType::IndexType start; start.Fill( 0 ); ImportImageFilterType::RegionType region; region.SetIndex( start ); region.SetSize( size ); // 画像原点の取得、設定 inImg->GetOrigin( org ); ImportImageFilterType::OriginType origin; origin[0] = org[0]; origin[1] = org[1]; origin[2] = org[2]; // ピクセルスペーシングの設定 inImg->GetSpacing( space ); ImportImageFilterType::SpacingType spacing; spacing[0] = space[0]; spacing[1] = space[1]; spacing[2] = space[2]; // 画像データの大きさ numberOfPixels = size[2] * size[1] * size[0]; sPtr_vin = static_cast< short* >( inImg->GetScalarPointer() ); sPtr_vout = static_cast< short* >( outImg->GetScalarPointer() ); // itk画像データの作成 ImportImageFilterType::Pointer importImage = ImportImageFilterType::New(); importImage->SetRegion( region ); importImage->SetOrigin( origin ); importImage->SetSpacing( spacing ); // vtkデータからitkデータへのコピー localBuffer = new short[ numberOfPixels ]; idx_start = 0*dims[1]*dims[0] + 0*dims[0] + 0; idx_end = (dims[2] -1)*dims[1]*dims[0] + (dims[1]-1)*dims[0] + (dims[0]-1); memcpy( localBuffer, &sPtr_vin[idx_start], sizeof(short)*numberOfPixels ); // メモリは自動的に解放 const bool importImageFilterWillOwnTheBuffer = true; importImage->SetImportPointer( localBuffer, numberOfPixels, importImageFilterWillOwnTheBuffer );
vtkImageDataから画像の原点やサイズ、ピクセルスペーシングそして画像データなどを取得し、ImportImageFilterにそれらををコピーします。
ここからいよいよクロージング処理です。
// クロージング処理
// 膨張count回→収縮count回
// 構造体の決定
StructuringElementType structuringElement;
structuringElement.SetRadius( /*お好きなサイズ*/ );
structuringElement.CreateStructuringElement();
// 膨張処理
DilateFilterType::Pointer* bdArr;
bdArr = new DilateFilterType::Pointer[count];
for( num = 0; num < count; num++ ){
DilateFilterType::Pointer bd;
bdArr[num] = DilateFilterType::New();
bd = bdArr[num];
bd->SetKernel( structuringElement );
bd->SetDilateValue( 1 );
bd->SetBackgroundValue( 0 );
if( !num ){
bd->SetInput( importImage->GetOutput() );
}
else{
bd->SetInput( bdArr[num-1]->GetOutput() );
}
bd->Update();
}
// 収縮処理
ErodeFilterType::Pointer* beArr;
beArr = new ErodeFilterType::Pointer[count];
for( num = 0; num < count; num++ ){
ErodeFilterType::Pointer be;
beArr[num] = ErodeFilterType::New();
be = beArr[num];
be->SetKernel(structuringElement);
be->SetErodeValue( 1 );
be->SetBackgroundValue( 0 );
if( !num ){
be->SetInput( bdArr[count-1]->GetOutput() );
}
else{
be->SetInput( beArr[num-1]->GetOutput());
}
be->Update();
}
StructuringElementType::SetRadiusではカーネルの半径を設定することができます。上記の実装では膨張をcount回行った後、収縮をcount回行っています。また、膨張ではSetDilateValue、収縮ではSetErodeValueに前景領域の値を設定します。今回は二値化画像（1:前景、0:背景)を用いるので両方とも1となります。
これでクロージングは終了です。
最後に、ITKでの画像データをVTKでの画像データにコピーして終わり！
short* sPtr_itk = beArr[num-1]->GetOutput()->GetBufferPointer(); short* sPtr_vout = static_cast< short* >( outImg->GetScalarPointer() ); // itkデータからvtkへのコピー unsingned long idx_start = 0; memcpy( &sPtr_vout[idx_start], sPtr_itk, sizeof(short)*numberOfPixels ); outImg->Update();
こちらがクロージング処理の結果です。これでセグメンテーションの全処理が終了です。
VTKによって結果を三次元画像、スライス画像の両方で表示しています。
いかがでしたでしょうか。今回のポイントは以下となります。
- クロージング処理によって穴を埋めることができる
- クロージングは膨張(Dilation)処理と収縮(Erosion)処理を行うこと
- カーネルはitk::BinaryBallStructuringElementクラス、膨張処理はBinaryDilateImageFilterクラス、収縮処理はBinaryErodeImageFilterクラスによって定義される
以上で肺のセグメンテーションは完了となります。
第一回から第四回にかけて、肺のセグメンテーションを例にしてITKの使い方を説明しました。変数の定義、クラスの実体化など、若干ややこしい印象を受けるかもしれませんが、慣れると非常に便利です。なぜなら、「フィルターの各種設定を行うだけで、自動的に目的の画像を作成」してくれるのですから。モルフォロジー変換といった簡単なアルゴリズムでも、自分で一から実装するのはかなり大変です。二次元ならまだしも、CT画像のような三次元になるとかなりの労力となります。また、思わぬバグも出てくるかもしれません。ITKを使えば、そのような労力や時間的コストを省くことができます。余った時間を新しいアルゴリズムやセオリーを開発するのに割り当てることもできますね。
「ITKを使った肺のセグメンテーション」の連載はこの回で終了です！
クァクァ。終わりクァ。