《itk实用demo》-FastMarching

FastMarching

#include "itkGradientMagnitudeImageFilter.h"
#include "itkThresholdSegmentationLevelSetImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkMaskImageFilter.h"
class Seeds
{
public:
    int _x;
    int _y;
    int _z;
    int value;
};
void FastMarching(ImageType::Pointer sourceimage, list<Seeds> list_seedPoint, float timeThreshold, float sigma, float alpha, float beta, double stoppingTime)
{

    typedef itk::RescaleIntensityImageFilter< ImageType, ImageType >  CastFilterType_f; //先将像素值类型转换为float型，以便进行平滑处理
    //typedef itk::RescaleIntensityImageFilter< InternalImageType, OutputImageType >  CastFilterType;
    CastFilterType_f::Pointer caster_f = CastFilterType_f::New();
    caster_f->SetOutputMinimum( -2048.0 );
    caster_f->SetOutputMaximum( 2048.0 );

    typedef itk::BinaryThresholdImageFilter<ImageType, ImageType> ThresholdingFilterType;
    ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
    thresholder->SetLowerThreshold(0.0);
    thresholder->SetUpperThreshold(timeThreshold);   //UpperThreshold相当于拍照的那一时刻
    thresholder->SetOutsideValue(  0  );
    thresholder->SetInsideValue(  255 );   //将时间为0至UpperThreshold的像素都置为白色

    typedef itk::CurvatureAnisotropicDiffusionImageFilter< ImageType, ImageType >  SmoothingFilterType;  //平滑
    SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
    typedef  itk::GradientMagnitudeRecursiveGaussianImageFilter< ImageType, ImageType >  GradientFilterType;  //求图像梯度
    typedef  itk::SigmoidImageFilter<ImageType, ImageType >  SigmoidFilterType;  //求传播速率
    GradientFilterType::Pointer  gradientMagnitude = GradientFilterType::New();
    SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();

    sigmoid->SetOutputMinimum( 0.0 );
    sigmoid->SetOutputMaximum( 1.0 );  //根据图像的梯度得出传播速率，速率从0.0至1.0
    sigmoid->SetAlpha( alpha );
    sigmoid->SetBeta( beta );

    typedef  itk::FastMarchingImageFilter< ImageType, ImageType >   FastMarchingFilterType;
    FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();
    //CastFilterType::Pointer caster = CastFilterType::New();
    //caster->SetOutputMinimum( 0 );
    //caster->SetOutputMaximum( 255 );
    //平滑系数
    smoothing->SetTimeStep( 0.125 );
    smoothing->SetNumberOfIterations(  5 );
    smoothing->SetConductanceParameter( 9.0 );

    gradientMagnitude->SetSigma( sigma );  //设置求梯度图像时用到的参数

    typedef FastMarchingFilterType::NodeContainer  NodeContainer;
    typedef FastMarchingFilterType::NodeType  NodeType;
    NodeContainer::Pointer seeds = NodeContainer::New();
    ImageType::IndexType  seedPosition;
    seeds->Initialize();
    NodeType node;
    const double seedValue = 0.0;  //种子点的值为0说明种子点的到达时间为0，然后行进过程中，其中像素点的值不断进行累加

    list<Seeds>::iterator i,iend;
    iend = Seeds.end();
    int j;
    for(i=Seeds.begin(),j=0; i!=iend; i++,j++)  //遍历list容器
    {
        seedPosition[0] = (*i).index_x;
        seedPosition[1] = (*i).index_y;
        node.SetValue( seedValue );
        node.SetIndex( seedPosition );
        seeds->InsertElement( j, node );
    }

    fastMarching->SetTrialPoints(  seeds  );   //设置一系列种子点，然后根据种子点进行传播
    fastMarching->SetOutputSize( caster_f->GetOutput()->GetBufferedRegion().GetSize() );
    fastMarching->SetStoppingValue(  stoppingTime  );

    //在管道中连接这些滤波器
    caster_f->SetInput( sourceimage );
    smoothing->SetInput( caster_f->GetOutput() );
    gradientMagnitude->SetInput( smoothing->GetOutput() );
    sigmoid->SetInput( gradientMagnitude->GetOutput() );
    fastMarching->SetInput( sigmoid->GetOutput() );
    thresholder->SetInput( fastMarching->GetOutput() );


    typedef itk::MaskImageFilter< ImageType, ImageType, ImageType >  MaskFilterType;
    MaskFilterType::Pointer maskFilter = MaskFilterType::New();
    maskFilter->SetInput1( sourceimage );
    maskFilter->SetInput2( thresholder->GetOutput());    

}