ITK中的水平集算法-快速步进分割

本文不在理论部分纠缠太多，大部分都是ITK中算法的代码实现。并且代码模块参考于参考文献，非原创。

ITK中通过itkFastMarchingImageFliter实现快速步进，方法原型是点与点之间的路径搜索，在搜索的过程中添加判断条件，进行点的取舍，从而构成了图像分割模型。FastMarching与Dijkstra类似，不同点在于Dijkstra用欧氏距离进行节点之间的更新，而FastMarching利用由Eikonal方程化简得到的近似偏微分方程进行更新。

快速步进水平集分割模型分三步：

1.使用边缘保护滤波器对输入图像平滑 CurvatureNDAnisotropicDiffusionFunction

2.计算梯度值后将其传给sigmoid滤波器 GradientMagnitudeRecursiveGaussianImageFilter  SigmoidImageFilter

3.将图像用FastMarchingImageFliter进行分割处理 FastMarchingImageFliter

其中，itkSigmoidImageFilter ，类itk::SigmoidImageFilter对输入图像的灰度值进行变换的规则如上式所示，用户调用类itk::SigmoidImageFilter里的方法SetOutputMaximum(), SetOutputMinimum()设置输出图像的灰度值最大值最小值，也就是式中的Max和Min。I表示输入图像的灰度值，I’表示输出图像的灰度值，α表示输入图像的灰度范围，β表示输入图像灰度范围的中心位置(这两个值有点类似窗宽和窗位)。它们的值的设置由方法：SetAlpha(), SetBeta()来完成。一般认为图像的灰度范围为[-3α,3α]之间，α和β的值可正可负。

也就是：对于输入图像的灰度值，如果其像素灰度值低于β-3α或者高于β+3α时，则这些灰度值将被映射到用户设置的Minimum和Maximum的值，如果输入图像的灰度值是在[β-3α, β+3α]之间时，则就会调用以上公式计算出输出图像的灰度值。

注意这部分显示的时候，不能直接用cast转换，否则纯黑，应该用RescaleIntensityImageFilter将图像灰阶转换到0-255才行。

itkCurvatureAnisotropicDiffusionImageFilter  // 实例化异向分散滤波器，对图像进行平滑操作

typedef   itk::CurvatureAnisotropicDiffusionImageFilter< InternalImageType, 

InternalImageType >  SmoothingFilterType;// 创建滤波器，对图像边缘进行保护和平滑操作

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

smoothing->SetInput( image_temp);

smoothing->SetTimeStep( 0.125 );//步长

smoothing->SetNumberOfIterations(  5 );//迭代次数

smoothing->SetConductanceParameter( 9.0 );// 电导系数

利用改进的曲率扩散方程的图像（MCDE）

itkAnisotropicDiffusionFunction and itkCurvatureNDAnisotropicDiffusionFunction.

输入和输出滤波器必须是一个标量itk:：数值类型的图像（float或double）

此过滤器的默认的时间步长被设置为最大的理论上稳定的值：0.5 / 2^n，其中n是图像的维数。

对于二维图像，有效的时间步骤是低于0.1250。

对于3D图像，有效的时间步骤是低于0.0625。

    // 实例化内部图像类型    typedef   unsigned char           InternalPixelType;    const     unsigned int    Dimension = 2;    typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;    // 实例化外部图像类型    typedef   unsigned char                            OutputPixelType;    typedef itk::Image< OutputPixelType, Dimension > OutputImageType;    //-------------------------------------------    typedef   float           floatPixelType;    typedef itk::Image< floatPixelType, Dimension >  floatImageType;    typedef itk::CastImageFilter<InternalImageType, floatImageType> CastFilterTypefloat;    CastFilterTypefloat::Pointer castFilter = CastFilterTypefloat::New();    castFilter->SetInput(sliceimage_current);    castFilter->Update();    // 实例化异向分散滤波器，对图像进行平滑操作    typedef   itk::CurvatureAnisotropicDiffusionImageFilter<        floatImageType,        floatImageType >  SmoothingFilterType;    // 创建滤波器，对图像边缘进行保护和平滑操作    SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();    smoothing->SetInput( castFilter->GetOutput());    smoothing->SetTimeStep( 0.125 );    smoothing->SetNumberOfIterations(  5 );    smoothing->SetConductanceParameter( 9.0 );// 电导系数    // 实例化高斯梯度放大递归滤波器    typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter<        floatImageType,        floatImageType >  GradientFilterType;    // 创建梯度强度滤波器，对平滑后的图像像素生成梯度    GradientFilterType::Pointer  gradientMagnitude = GradientFilterType::New();    gradientMagnitude->SetInput( smoothing->GetOutput() );    const double sigma =0.5;    gradientMagnitude->SetSigma(  sigma  );// 设置σ    // 实例化Sigmoid滤波器，用于计算前向传播的速率系数    typedef   itk::SigmoidImageFilter<floatImageType,floatImageType >  SigmoidFilterType;    // 创建sigmoid函数，用于计算前向传播的速率系数    SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();    sigmoid->SetInput( gradientMagnitude->GetOutput() );    sigmoid->SetOutputMinimum(  0.0  );    sigmoid->SetOutputMaximum(  1.0  );    const double alpha =  -3;    const double beta  =  3.0;    sigmoid->SetAlpha( alpha );    sigmoid->SetBeta(  beta  );    //-----------------------------------------------------------    // 实例化快速步进滤波器    typedef  itk::FastMarchingImageFilter< floatImageType,floatImageType > FastMarchingFilterType;    // 快速步进滤波器    FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();    fastMarching->SetInput( sigmoid->GetOutput() );    // 设置轮廓扩张的种子点    typedef FastMarchingFilterType::NodeContainer           NodeContainer;    typedef FastMarchingFilterType::NodeType                NodeType;    NodeContainer::Pointer seeds = NodeContainer::New();    InternalImageType::IndexType  seedPosition;// 种子的坐标    seeds->Initialize();                       // 初始化，并插入种子    for (int i=0;i<centerline.size();i++)    {        //double *point=points->GetPoint(i);        InternalImageType::IndexType  seedPosition;// 种子的坐标        seedPosition[0] = centerline.at(i).index_x;        seedPosition[1] = centerline.at(i).index_y;        NodeType node;        // 种子值        node.SetValue( 0 );        node.SetIndex( seedPosition );        seeds->InsertElement( i, node );    }    fastMarching->SetTrialPoints(  seeds  );    fastMarching->SetOutputSize(sliceimage_current->GetBufferedRegion().GetSize() );// 设置输出图像尺寸    /*const double stoppingTime = 30;*/    fastMarching->SetStoppingValue(  fastmaichingparameters.stoppingTime  );    fastMarching->SetNormalizationFactor(1);    // 实例化二值门限图像滤波器，对fastMarching输出的时间交叉图进行像素处理    typedef itk::BinaryThresholdImageFilter< floatImageType,        floatImageType> ThresholdingFilterType;    ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();    thresholder->SetInput( fastMarching->GetOutput() );    // 设置上下门限值    const InternalPixelType  timeThreshold = 100;    thresholder->SetLowerThreshold(           0.0  );    thresholder->SetUpperThreshold( timeThreshold  );    // 设置门限内、外像素值（亮度）    thresholder->SetOutsideValue(  0  );    thresholder->SetInsideValue(  255 );    //-----------------------------------------------------------    // 像素亮度强度滤波器，将输入像素投射到输出图像的范围内    typedef itk::RescaleIntensityImageFilter<        InternalImageType,        OutputImageType >   CastFilterType;    CastFilterType::Pointer caster1 = CastFilterType::New();    caster1->SetInput( thresholder->GetOutput() );    caster1->SetOutputMinimum(   0 );    caster1->SetOutputMaximum( 255 );

参考文献：

《三维图像编程实验》

http://www.cnblogs.com/shushen/archive/2016/04/12/5381753.html

http://blog.163.com/jacky_ling0/blog/static/137392571201092183013782/