ITK系列30_水平集（快速步进）算法对脑部PNG图像进行二维分割

ITK系列30_⽔平集（快速步进）算法对脑部PNG图像进⾏⼆

维分割

⽔平集分割

⽔平集是跟踪轮廓和表⾯运动的⼀种数字化⽅法。不直接对轮廓进⾏操作，⽽是将轮廓设置成⼀个⾼维函数的零⽔平集，这个⾼维函数叫做⽔平集函数： Ψ(X,t) 。然后⽔平集函数运动成为⼀个微分⽅程。在任何时候，通过从输出中提取零⽔平集来得到运动的轮廓。使⽤⽔平集的主要优点是可以对任何复杂的结构进⾏模式化和拓扑变换，⽐如暗中操作融合和分离。通过使⽤基于图像的诸如亮度均值、梯度和边缘之类的特征的微分⽅程的解答，⽔平集就可以⽤来对图像进⾏分割。在⼀个典型的⽅法中，⽤户对⼀个轮廓进⾏初始化并运动这个轮廓直到它符合图像中的⼀个解剖结构。在这个⽂献中，发表了许多这些基本概念的执⾏和变量。赛思 Sethian 已经对这⼀领域做了⼀个综述。

快速步进分割

当微分⽅程控制⽔平集运动时有⼀个⾮常简单的形式，可以使⽤⼀个快速运动算法叫做快速步进。接下来的例⼦阐述了

itk::FastMarchingImageFilter 的⽤法。这个滤波器对⼀个简单⽔平集运动问题执⾏⼀个快速⾏进⽅法。

在这个例⼦中，微分⽅程中使⽤的速率系数是通过⽤户以⼀种图像的形式来提供的。通常作为⼀个梯度值的函数来计算这个图像。在⽂献中有⼏个映射是很常见的，⽐如负指数 exp( − x) 和倒数 1/(1+x) 。在当前的例⼦中我们决定使⽤⼀种 S 函数Sigmoid 函数，因为它可以提供⼀个控制参数的好⽅法，这个参数可以定制成⼀个很好的速率图像形状。这个映射按以下⽅式进⾏：当图像梯度很⾼时前⾯的传播速率很慢⽽在梯度很低的区域传播速率很快。这种⽅式将导致轮廓不断传播直到它到达图像中的解剖结构边缘并在那些边

缘前将速率降低下来。 FastMarchingImageFilter 的输出是⼀个时间交叉图，对每⼀个像素，它还表达了到达这个像素位置之前所花费的时间。然后在输出图像中⼀个门限的应⽤等同于在轮廓运动过程中的⼀个特殊时间上对轮廓取⼀个快照。轮廓横渡到⼀个特殊解剖结构边缘被认为将花费更长的时间。这将导致在接近结构边缘的时间交叉像素值发⽣很⼤的变化。通过调⽤⼀个时间范围⽤这个滤波器来执⾏分割，在这个时间范围内图像空间中⼀个区域将在⼀段长时间内包含轮廓。如图 9-15 所⽰显⽰了包含在⼀个分割任务中的FastMarchingImageFilter 的应⽤的主要成员。它包括⼀个使⽤ itk::CurvatureAnisotropicDiffusionImageFilter 进⾏平滑的最初阶段。平滑后的图像作为输⼊传递给 itk::GradientMagnitudeRecursiveGaussianImageFilter ，然后再给

itk::SigmoidImageFilter 。最后，将 FastMarchingImageFilter 的输出传递给 itk::BinaryThresholdImageFilter 以便产⽣⼀个binary mask 来表达分割对象。

接下来例⼦中的代码阐述了⼀个⽤来执⾏快速步进分割的管道的典型设置。⾸先，使⽤⼀个边缘保护滤波器对输⼊图像进⾏平滑；然后计算它的梯度值并传递给⼀个 sigmoid 滤波器。 sigmoid 滤波器的结果是图像的潜能，将⽤来影响微分⽅程的速率系数。

实例30 快速步进算法对脑部PNG图像进⾏⼆维分割

//包含⽤来从输⼊图像中去除噪声头⽂件

#include "itkCurvatureAnisotropicDiffusionImageFilter.h"

//这两个滤波器连在⼀起将产⽣调节描述⽔平集运动的微分⽅程中的速率系数的图像潜能。

#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"

#include "itkSigmoidImageFilter.h"

#include "itkFastMarchingImageFilter.h"

//从 FastMarchingImageFilter 得到的时间交叉图将使⽤ BinaryThresholdImageFilter 来进⾏门限处理

#include "itkBinaryThresholdImageFilter.h"

//使⽤ itk::ImageFileReader 和 itk::ImageFileWriter 来对图像进⾏读写

#include "itkImageFileReader.h"

#include "itkImageFileWriter.h"

#include "itkRescaleIntensityImageFilter.h"

static void PrintCommandLineUsage( const int argc, const char * const argv[] )

{

std::cerr << "Missing Parameters " << std::endl;

std::cerr << "Usage: " << argv[0];

std::cerr << " inputImage outputImage seedX seedY";

std::cerr << " Sigma SigmoidAlpha SigmoidBeta TimeThreshold StoppingValue";

std::cerr << " smoothingOutputImage gradientMagnitudeOutputImage sigmoidOutputImage" << std::endl;

for (int qq=0; qq< argc; ++qq)

{

std::cout << "argv[" << qq << "] = " << argv[qq] << std::endl;

}

int main( int argc, char *argv[] )

{

/*if (argc != 13)

{

PrintCommandLineUsage(argc, argv);

return EXIT_FAILURE;

}*/

/*现在我们使⽤⼀个像素类型和⼀个特殊维来定义图像类型。由于需要⽤到平滑滤波器，

所以在这种情况下对像素使⽤浮点型数据类型*/

typedef float InternalPixelType;

const unsigned int Dimension = 2;

typedef itk::Image< InternalPixelType, Dimension > InternalImageType;

//另⼀⽅⾯，输出图像被声明为⼆值的

typedef unsigned char OutputPixelType;

typedef itk::Image< OutputPixelType, Dimension > OutputImageType;

//下⾯使⽤图像内在的类型和输出图像类型来对 BinaryThresholdImageFilter 的类型进⾏实

//例化：

typedef itk::BinaryThresholdImageFilter< InternalImageType,

OutputImageType > ThresholdingFilterType;

ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();

//设置阈值

const InternalPixelType timeThreshold = atof("100");

/*传递给 BinaryThresholdImageFilter 的上门限将定义那个从时间交叉图得来的时间快照。药膳论文

在⼀个理想的应⽤中，⽤户拥有使⽤视觉反馈交互式的选择这个门限的能⼒。在这⾥，由于

是⼀个最⼩限度的例⼦，从命令⾏来得到这个值*/

thresholder->SetLowerThreshold( 0.0 );

thresholder->SetUpperThreshold( timeThreshold );

thresholder->SetOutsideValue( 0 );

制定班级公约

thresholder->SetInsideValue( 255 );

//下⾯⼏⾏我们对读写器类型进⾏实例化

typedef itk::ImageFileReader< InternalImageType > ReaderType;

typedef itk::ImageFileWriter< OutputImageType > WriterType;

ReaderType::Pointer reader = ReaderType::New();

WriterType::Pointer writer = WriterType::New();

//输⼊图像

reader->SetFileName("BrainProtonDensitySlice.png");

//输出图像宝钢电子商务平台

writer->SetFileName("BrainProtonDensitySlice_FastMarching.png");

typedef itk::RescaleIntensityImageFilter<

InternalImageType,

OutputImageType > CastFilterType;

//使⽤图像内在的类型对 CurvatureAnisotropicDiffusionImageFilter 类型进⾏实例化

typedef itk::CurvatureAnisotropicDiffusionImageFilter<

InternalImageType,

InternalImageType > SmoothingFilterType;

//然后，通过调⽤ New( ) ⽅式来创建滤波器并将结果指向⼀个 itk::SmartPointer

SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();

/*使⽤图像内在的类型来对 GradientMagnitudeRecursiveGaussianImageFilter 和

SigmoidImageFilter 的类型进⾏实例化*/

typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<

InternalImageType,

InternalImageType > GradientFilterType;

typedef itk::SigmoidImageFilter<

InternalImageType,

InternalImageType > SigmoidFilterType;

//使⽤ New( ) ⽅式来对相应的滤波器进⾏实例化

GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();

SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();

/*使⽤ SetOutputMinimum() 和 SetOutputMaximum() ⽅式来定义 SigmoidImageFilter 输出的

最⼩值和最⼤值。在我们的例⼦中，我们希望这两个值分别是 0.0 和 1.0 ，以便得到⼀个最佳

的速率图像来给 MarchingImageFilter 。在 6.3.2 ⼩节详细地表达了 SigmoidImageFilter 的使⽤⽅

法。*/

sigmoid->SetOutputMinimum( 0.0 );

sigmoid->SetOutputMaximum( 1.0 );

//现在我们声明 FastMarchingImageFilter 的类型：

typedef itk::FastMarchingImageFilter< InternalImageType,

InternalImageType > FastMarchingFilterType;

/然后我们使⽤ New( ) ⽅式结构化这个类的⼀个滤波器

FastMarchingFilterType::Pointer fastMarching

= FastMarchingFilterType::New();

//现在我们使⽤下⾯的代码在⼀个管道中连接这些滤波器

smoothing->SetInput( reader->GetOutput() );

gradientMagnitude->SetInput( smoothing->GetOutput() );

sigmoid->SetInput( gradientMagnitude->GetOutput() );

fastMarching->SetInput( sigmoid->GetOutput() );

thresholder->SetInput( fastMarching->GetOutput() );

writer->SetInput( thresholder->GetOutput() );

/*CurvatureAnisotropicDiffusionImageFilter 类需要定义⼀对参数。下⾯是⼆维图像中的典

型值。然⽽它们可能需要根据输⼊图像中存在的噪声的数量进⾏适当的调整。在 6.7.3 ⼩节中

讨论了这个滤波器*/

smoothing->SetTimeStep( 0.125 );

smoothing->SetNumberOfIterations( 5 );

smoothing->SetConductanceParameter( 9.0 );

//设置Sigma(σ )

const double sigma = atof("1.0");

/*执⾏ GradientMagnitudeRecursiveGaussianImageFilter 就等同于通过⼀个派⽣的操作符紧

跟⼀个⾼斯核的⼀个回旋。这个⾼斯的 sigma(σ) 可以⽤来控制图像边缘影响的范围。在 6.4.2

⼩节中讨论了这个滤波器*/

gradientMagnitude->SetSigma( sigma );

/设置SigmoidAlpha(α) SigmoidBeta(β)

const double alpha = atof("-0.5");

const double beta = atof("3.0");

/*SigmoidImageFilter 类需要两个参数来定义⽤于 sigmoid 变量的线性变换。使⽤ SetAlpha()和 SetBeta() ⽅式来传递这些参数。在这个例⼦的背景中，参数⽤来强化速率图像中的⾼低值

区域之间的区别。在理想情况下，解剖结构中的均匀区域内的速率应该是 1.0 ⽽在结构边缘

附近应该迅速地衰减到 0.0 。下⾯是寻这个值的启发。在梯度量值图像中，我们将沿被分

割的解剖结构地轮廓的最⼩值称为 K1 ，将结构中间的梯度强度的均值称为 K2 。这两个值表

达了我们想要映射到速率图像中间距为[0:1] 的动态范围。我们期望 sigmoid 将 K1 映射为 0.0 ，

K2 为 1.0 。由于 K1 的值要⽐ K2 更⼤⽽我们期望将它们分别映射为 0.0 和 1.0 ，所以我们对 alpha(α)选择⼀个负值以便 sigmoid 函数也做⼀个反转亮度映射。这个映射将产⽣⼀个速率图像，其

中⽔平集将在均匀区域快速⾏进⽽在轮廓上停⽌。 Beta(β) 的参考值是(K1 + K2) / 2 ⽽α的参考值是(K

2 - K1) / 6 必须是⼀个负数。在我们的简单例⼦中，这些值是⽤户从命令⾏提供的。⽤户

可以通过留⼼梯度量值图像来估计这些值*/

sigmoid->SetAlpha( alpha );

sigmoid->SetBeta( beta );

/*FastMarchingImageFilter 需要⽤户提供⼀个轮廓扩张的种⼦点。⽤户实际上不仅要提供⼀个种⼦点⽽是⼀个种⼦点集。⼀个好的种⼦点集将增⼤分割⼀个复杂对象⽽不丢失数据的

可能性。使⽤多种⼦也会减少访问整个对象所需要的时间同时减少前⾯访问的区域边缘的漏

洞的风险。例如，当分割⼀个拉长的对象时，在其中⼀端设置⼀个种⼦将导致传播到另⼀端

会需要很长⼀段时间。沿着轴线设置⼏个种⼦⽆疑将是确定整个对象尽早被捕获的最佳策

略。⽔平集的⼀个重要的性质就是它们不需要任何额外辅助来融合⼏个起点的⾃然能⼒。使

⽤多种⼦正是利⽤了这个性质。

这些种⼦储存在⼀个容器中。在 FastMarchingImageFilter 的特性中将这个容器的类型定

义为 NodeContainer ：*/

typedef FastMarchingFilterType::NodeContainer NodeContainer;

typedef FastMarchingFilterType::NodeType NodeType;

NodeContainer::Pointer seeds = NodeContainer::New();

InternalImageType::IndexType seedPosition;

//设置种⼦点位置

seedPosition[0] = atoi("99");

seedPosition[1] = atoi("114");

//节点作为堆栈变量来创建，并使⽤⼀个值和⼀个 itk::Index 位置来进⾏初始化

NodeType node;

const double seedValue = 0.0;

node.SetValue( seedValue );

node.SetIndex( seedPosition );

//节点列表被初始化，然后使⽤ InsertElement( ) 来插⼊每个节点：

seeds->Initialize();

seeds->InsertElement( 0, node );

//现在使⽤ SetTrialPoints( ) ⽅式将种⼦节点集传递给 FastMarchingImageFilter ：

fastMarching->SetTrialPoints( seeds );

/*FastMarchingImageFilter 需要⽤户指定作为输出产⽣的图像的⼤⼩。使⽤ SetOutputSize()来完成。注意：这⾥的⼤⼩是从平滑滤波器的输出图像得到的。这个图像的⼤⼩仅仅在直接或间接地调⽤了这个滤波器的 Updata() 之后才是有效的。*/

//writer 上的 Updata( ) ⽅法引发了管道的运⾏。通常在出现错误和抛出异议时 , 从⼀个

//smoothingOutputImage

//使⽤边缘保护平滑滤波器平滑滤波后的图像

try

{

CastFilterType::Pointer caster1 = CastFilterType::New();

WriterType::Pointer writer1 = WriterType::New();

caster1->SetInput( smoothing->GetOutput() );

writer1->SetInput( caster1->GetOutput() );

writer1->SetFileName("smoothingOutputImage.png");

caster1->SetOutputMinimum( 0 );

caster1->SetOutputMaximum( 255 );

writer1->Update();

}

catch( itk::ExceptionObject & err )

{

std::cerr << "ExceptionObject caught !" << std::endl;

std::cerr << err << std::endl;

return EXIT_FAILURE;

}

// gradientMagnitudeOutputImage

//滤波后图像的梯度强度图像

try

{

CastFilterType::Pointer caster2 = CastFilterType::New();

WriterType::Pointer writer2 = WriterType::New();

caster2->SetInput( gradientMagnitude->GetOutput() );

writer2->SetInput( caster2->GetOutput() );

writer2->SetFileName("gradientMagnitudeOutputImage.png");

caster2->SetOutputMinimum( 0 );

caster2->SetOutputMaximum( 255 );

writer2->Update();

}

catch( itk::ExceptionObject & err )

{

std::cerr << "ExceptionObject caught !" << std::endl;

std::cerr << err << std::endl;

return EXIT_FAILURE;

}

// sigmoidOutputImage

//梯度强度的 sigmoid 函数图像

try

{

CastFilterType::Pointer caster3 = CastFilterType::New();

WriterType::Pointer writer3 = WriterType::New();

caster3->SetInput( sigmoid->GetOutput() );

writer3->SetInput( caster3->GetOutput() );

// FastMarchingImageFilter 分割处理产⽣的结果图像

writer3->SetFileName("sigmoidOutputImage.png");

caster3->SetOutputMinimum( 0 );

caster3->SetOutputMaximum( 255 );

writer3->Update();

}

catch( itk::ExceptionObject & err )

{

std::cerr << "ExceptionObject caught !" << std::endl;

std::cerr << err << std::endl;

return EXIT_FAILURE;

}

fastMarching->SetOutputSize(

够快物流网reader->GetOutput()->GetBufferedRegion().GetSize() );

//设置stoppingTime，⽐阈值（门限值⾼⼀点）

const double stoppingTime = atof("103");

/*由于表⽰轮廓的起点将从头到尾不断传播，⼀旦到达⼀个特定的时间就希望停⽌这⼀过程。这就允许我们在假定相关区域已经计算过的假设下节省计算时间。使⽤ SetStopping Value() ⽅式来定义停⽌过程的值。原则上，这个停⽌值应该⽐门限值⾼⼀点*/

fastMarching->SetStoppingValue( stoppingTime );

// try / catch 模块调⽤ updata 。

try

{

writer->Update();

}

catch( itk::ExceptionObject & excep )

{

std::cerr << "Exception caught !" << std::endl;

std::cerr << excep << std::endl;

return EXIT_FAILURE;

}

try

被禁止的爱

{

CastFilterType::Pointer caster4 = CastFilterType::New();

WriterType::Pointer writer4 = WriterType::New();

caster4->SetInput( fastMarching->GetOutput() );

writer4->SetInput( caster4->GetOutput() );

writer4->SetFileName("FastMarchingFilterOutput4.png");

caster4->SetOutputMinimum( 0 );

caster4->SetOutputMaximum( 255 );

writer4->Update();

}

catch( itk::ExceptionObject & err )

{

std::cerr << "ExceptionObject caught !" << std::endl;

std::cerr << err << std::endl;

return EXIT_FAILURE;

}

typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;

InternalWriterType::Pointer mapWriter = InternalWriterType::New();

mapWriter->SetInput( fastMarching->GetOutput() );

mapWriter->SetFileName("FastMarchingFilterOutput4.mha");

mapWriter->Update();

金属型铸造机InternalWriterType::Pointer speedWriter = InternalWriterType::New();

speedWriter->SetInput( sigmoid->GetOutput() );

本文发布于:2024-09-23 00:41:36，感谢您对本站的认可！

本文链接：https://www.17tex.com/xueshu/307856.html

上一篇：裂纹检测算法

下一篇：数字图像处理复习题

标签：图像时间轮廓速率梯度

留言与评论（共有 0 条评论）