#include "itkConnectedThresholdImageFilter.h"//连接门限头文件
#include "itkImage.h"
#include "itkCastImageFilter.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
//图像中存在的噪声将大大降低滤波器生长大面积区域的能力。当面对噪声图像时，通常
//是使用一个边缘保留平滑滤波器。
int main( int argc, char *argv[])
{/*if( argc < 7 ){std::cerr << "Missing Parameters " << std::endl;std::cerr << "Usage: " << argv[0];std::cerr << " inputImage  outputImage seedX seedY lowerThreshold upperThreshold" << std::endl;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;typedef itk::CastImageFilter< InternalImageType, OutputImageType >CastingFilterType;CastingFilterType::Pointer caster = CastingFilterType::New();//图像读取与图像写类型定义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_huizhi.png" );//使用图像类型作为模板参数来对平滑滤波器进行实例化typedef itk::CurvatureFlowImageFilter< InternalImageType, InternalImageType >CurvatureFlowImageFilterType;//调用 New() 方式来创建滤波器并将接指向 itk::SmartPointer//平滑滤波器实例化对象smoothingCurvatureFlowImageFilterType::Pointer smoothing =CurvatureFlowImageFilterType::New();//声明区域生长滤波器的类型，本例中使用 ConnectedThresholdImageFiltertypedef itk::ConnectedThresholdImageFilter< InternalImageType,InternalImageType > ConnectedFilterType;//使用 New( ) 方式构造这种类的一个滤波器//连接门限滤波器实例化对象connectedThresholdConnectedFilterType::Pointer connectedThreshold = ConnectedFilterType::New();//读取图像进行平滑滤波smoothing->SetInput( reader->GetOutput() );//滤波后进行连接门限connectedThreshold->SetInput( smoothing->GetOutput() );/*由于只有一小部分图像文件格式支持浮点型数据类型，所以使用 cast filter 将浮点型数据类型转换成整型*/caster->SetInput( connectedThreshold->GetOutput() );//处理后输出数据到writerwriter->SetInput( caster->GetOutput() );/*CurvatureFlowImageFilter（平滑滤波器）需要定义两个参数。下面是一个二维图像的常见值。当然它们也需要根据输入图像存在的噪声的数量进行适当的调整*/smoothing->SetNumberOfIterations( 5 );smoothing->SetTimeStep( 0.125 );/*ConnectedThresholdImageFilter（连通门限图像滤波）有两个主要的参数（lowerThreshold和upperThreshold）需要定义，它们分别是为了确定是否包含在区域中的亮度值而制定的标准的上门限和下门限。这两个值设定得太接近势必会降低区域生长的机动性，而设定得太远必将整个图像都卷入区域中*/const InternalPixelType lowerThreshold = atof( "180" );const InternalPixelType upperThreshold = atof( "210" );connectedThreshold->SetLower(  lowerThreshold  );connectedThreshold->SetUpper(  upperThreshold  );/*这个滤波器的输出是一个二值图像，这个二值图像除了分割出的区域外到处都是零值像素。区域中的亮度值是由 SetReplaceValue() 方式来选择的*/connectedThreshold->SetReplaceValue( 255 );InternalImageType::IndexType  index;/*这个算法的实例化需要用户提供一个种子点index。将这个点选在被分割的解剖学结构的典型区域是很便捷的。种子是以一种 itk::Index 的形式传递给 SetSeed() 方式的*/index[0] = atoi( "107" );index[1] = atoi( "69" );connectedThreshold->SetSeed( index );/*writer 上的 Updata() 方法引发了管道的运行。通常在出现错误和抛出异议时, 从一个try / catch 模块调用 updata ：*/try{writer->Update();}catch( itk::ExceptionObject & excep ){std::cerr << "Exception caught !" << std::endl;std::cerr << excep << std::endl;}return EXIT_SUCCESS;
}

#include "itkOtsuThresholdImageFilter.h"//Otsu分割头文件
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"int main( int argc, char * argv[] )
{/*if( argc < 5 ){std::cerr << "Usage: " << argv[0];std::cerr << " inputImageFile outputImageFile ";std::cerr << " insideValue    outsideValue   "  << std::endl;return EXIT_FAILURE;}*///确定作为输入和输出图像的像素类型typedef  unsigned char  InputPixelType;typedef  unsigned char  OutputPixelType;const unsigned int      Dimension = 2;//使用输入输出图像的像素类型和维来定义它们的图像类型typedef itk::Image< InputPixelType,  Dimension >   InputImageType;typedef itk::Image< OutputPixelType, Dimension >   OutputImageType;//使用上面定义的输入输出图像的类型对滤波器类型进行实例化typedef itk::OtsuThresholdImageFilter<InputImageType, OutputImageType >  FilterType;//对一个 itk::ImageFileReader 类也进行实例化以便从文件中读取图像数据typedef itk::ImageFileReader< InputImageType >  ReaderType;typedef itk::ImageFileWriter< OutputImageType >  WriterType;//调用 New( ) 方式创建滤波器和 reader 并将结果指向 itk::Smartpointers ReaderType::Pointer reader = ReaderType::New();FilterType::Pointer filter = FilterType::New();WriterType::Pointer writer = WriterType::New();//由 reader 得到的图像作为输入传递给 OtsuThresholdImageFilterreader->SetFileName("BrainProtonDensitySlice.png");filter->SetInput(reader->GetOutput());writer->SetInput( filter->GetOutput() );//outsideValue:给定阈值的上下限范围之外的像素的亮度值//insideValue:给定阈值的上下限范围之内的像素的亮度值const OutputPixelType outsideValue = atoi( "0" );//范围之外设置为纯黑const OutputPixelType insideValue  = atoi( "255" );//范围之内设置为纯白/*用 SetOutsideValue() 方式定义指向那些亮度值在给定阈值的上下限范围之外的像素的亮度值，用 SetInsideValue() 方式定义指向那些亮度值在给定阈值的上下限范围之内的像素的亮度值*/filter->SetOutsideValue( outsideValue );filter->SetInsideValue(  insideValue  );try{filter->Update();}catch( itk::ExceptionObject & excp ){std::cerr << "Exception thrown " << excp << std::endl;}//由滤波器内部计算得到的阈值。为了做到这些我们调用 GetThreshold 方式int threshold = filter->GetThreshold();std::cout << "Threshold = " << threshold << std::endl;writer->SetFileName( "BrainProtonDensitySlice_OTSU2.png" );try{writer->Update();}catch( itk::ExceptionObject & excp ){std::cerr << "Exception thrown " << excp << std::endl;}return EXIT_SUCCESS;
}

#include "itkOtsuMultipleThresholdsCalculator.h"//包含头文件#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkScalarImageToHistogramGenerator.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkNumericTraits.h"#include <iomanip>
#include <stdio.h>
//这个例子阐述了如何使用 itk::OtsuMultipleThresholdsCalculator
int main( int argc, char * argv[] )
{/*if( argc < 5 ){std::cerr << "Usage: " << argv[0];std::cerr << " inputImageFile outputImageFileBase ";std::cerr << "  outputImageFileExtension numberOfThresholdsToCalculate "  << std::endl;return EXIT_FAILURE;}*///Convenience typedefstypedef  unsigned short  InputPixelType;typedef  unsigned char   OutputPixelType;typedef itk::Image< InputPixelType,  2 >   InputImageType;typedef itk::Image< OutputPixelType, 2 >   OutputImageType;/*OtsuMultipleThresholdsCalculator 以最大化类之间的变异的原则来计算一个给定的直方图的门限。我们使用 ScalarImageToHistogramGenerator 来得到直方图*/typedef itk::Statistics::ScalarImageToHistogramGenerator<InputImageType > ScalarImageToHistogramGeneratorType;typedef ScalarImageToHistogramGeneratorType::HistogramType HistogramType;typedef itk::OtsuMultipleThresholdsCalculator< HistogramType > CalculatorType;typedef itk::ImageFileReader< InputImageType >  ReaderType;typedef itk::ImageFileWriter< OutputImageType > WriterType;//一旦计算得到门限值，我们将使用 BinaryThresholdImageFilter 来对图像进行分割typedef itk::BinaryThresholdImageFilter<InputImageType, OutputImageType >  FilterType;ScalarImageToHistogramGeneratorType::Pointer scalarImageToHistogramGenerator= ScalarImageToHistogramGeneratorType::New();CalculatorType::Pointer calculator = CalculatorType::New();FilterType::Pointer filter = FilterType::New();ReaderType::Pointer reader = ReaderType::New();WriterType::Pointer writer = WriterType::New();scalarImageToHistogramGenerator->SetNumberOfBins( 128 );//设置Otsu分割阈值数目calculator->SetNumberOfThresholds( atoi("3") );const OutputPixelType outsideValue = 0;const OutputPixelType insideValue = 255;filter->SetOutsideValue( outsideValue );filter->SetInsideValue(  insideValue  );//Connect Pipelinereader->SetFileName( "BrainProtonDensitySlice.png" );//管道如下所示scalarImageToHistogramGenerator->SetInput( reader->GetOutput() );calculator->SetInputHistogram(scalarImageToHistogramGenerator->GetOutput() );filter->SetInput( reader->GetOutput() );writer->SetInput( filter->GetOutput() );//读取报错异常try{reader->Update();}catch( itk::ExceptionObject & excp ){std::cerr << "Exception thrown while reading image" << excp << std::endl;}scalarImageToHistogramGenerator->Compute();//calculator报错异常try{calculator->Compute();}catch( itk::ExceptionObject & excp ){std::cerr << "Exception thrown " << excp << std::endl;}//使用 GetOutput 方式可以得到门限值const CalculatorType::OutputType &thresholdVector = calculator->GetOutput();std::string outputFileBase = "BrainProtonDensitySlice_OtsuMu";InputPixelType lowerThreshold = itk::NumericTraits<InputPixelType>::min();InputPixelType upperThreshold;typedef CalculatorType::OutputType::const_iterator ThresholdItType;for( ThresholdItType itNum = thresholdVector.begin();itNum != thresholdVector.end();++itNum ){std::cout << "OtsuThreshold["<< (int)(itNum - thresholdVector.begin())<< "] = "<< static_cast<itk::NumericTraits<CalculatorType::MeasurementType>::PrintType>(*itNum)<< std::endl;upperThreshold = static_cast<InputPixelType>(*itNum);filter->SetLowerThreshold( lowerThreshold );filter->SetUpperThreshold( upperThreshold );lowerThreshold = upperThreshold;std::ostringstream outputFilename;outputFilename << outputFileBase<< std::setfill('0') << std::setw(3) << (itNum - thresholdVector.begin())<< "."<< "png";writer->SetFileName( outputFilename.str() );try{writer->Update();}catch( itk::ExceptionObject & excp ){std::cerr << "Exception thrown " << excp << std::endl;}}upperThreshold = itk::NumericTraits<InputPixelType>::max();filter->SetLowerThreshold( lowerThreshold );filter->SetUpperThreshold( upperThreshold );std::ostringstream outputFilename2;outputFilename2 << outputFileBase<< std::setfill('0') << std::setw(3) << thresholdVector.size()<< "."<< "jpg";writer->SetFileName( outputFilename2.str() );try{writer->Update();}catch( itk::ExceptionObject & excp ){std::cerr << "Exception thrown " << excp << std::endl;}return EXIT_SUCCESS;
}

#include "itkNeighborhoodConnectedImageFilter.h"
#include "itkImage.h"
#include "itkCastImageFilter.h"
//使用 itk::CurvatureFlowImageFilter 在保护边缘时平滑图像
#include "itkCurvatureFlowImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"int main( int argc, char *argv[] )
{/*if( argc < 7 ){std::cerr << "Missing Parameters." << std::endl;std::cerr << "Usage: " << argv[0]<< " inputImage outputImage"<< " seedX seedY"<< " lowerThreshold upperThreshold" << std::endl;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;typedef itk::CastImageFilter< InternalImageType, OutputImageType >CastingFilterType;CastingFilterType::Pointer caster = CastingFilterType::New();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( "Neighborhood_huizhi2.png");typedef itk::CurvatureFlowImageFilter<InternalImageType, InternalImageType>CurvatureFlowImageFilterType;//调用 New( ) 方式来定义滤波器并将结果指向一个 itk::SmartPointer CurvatureFlowImageFilterType::Pointer smoothing =CurvatureFlowImageFilterType::New();//现在我们声明区域生长滤波器的类型。这种情况下是 NeighborhoodConnectedImageFiltertypedef itk::NeighborhoodConnectedImageFilter<InternalImageType,InternalImageType > ConnectedFilterType;//使用 New( ) 方式对这个类的一个文件进行结构化ConnectedFilterType::Pointer neighborhoodConnected= ConnectedFilterType::New();/*现在到了连接一个简单的线性管道的时候。在管道的开头添加一个 reader 文件并在末尾添加一个 cast filter 和 writer 。由于只有仅仅一小部分图像文件格式支持浮点型数据类型，所以使用 cast filter 将浮点型数据类型转换成整型*/smoothing->SetInput( reader->GetOutput() );neighborhoodConnected->SetInput( smoothing->GetOutput() );caster->SetInput( neighborhoodConnected->GetOutput() );writer->SetInput( caster->GetOutput() );//滤波/*CurvatureFlowImageFilter 需要定义两个参数。下面是一个二维图像的常见值。当然它们也需要根据输入图像存在的噪声的数量进行适当的调整*/smoothing->SetNumberOfIterations( 5 );smoothing->SetTimeStep( 0.125 );const InternalPixelType lowerThreshold = atof( "180" );const InternalPixelType upperThreshold = atof( "210" );/*NeighborhoodConnectedImageFilter 有两个主要的参数需要设定，它们分别是为了确定是否包含在区域中的亮度值而制定的标准的上门限和下门限。这两个值设定得太接近势必会降低区域生长的机动性，而设定得太远必将整个图像都卷入区域中*/neighborhoodConnected->SetLower( lowerThreshold );neighborhoodConnected->SetUpper( upperThreshold );/*这里我们增加定义邻域大小的一个至关重要的参数，用来决定一个像素是否在这个区域中。邻域越大，这个滤波器在输入图像中滤掉噪声的稳定性就越强，但是计算需要的时间也将更长。这里我们选择每个维上两个长度r（2r+1)的一个滤波器，这将产生一个 5×5 像素的邻域*/InternalImageType::SizeType radius;radius[0] = 2;   // two pixels along Xradius[1] = 2;   // two pixels along YneighborhoodConnected->SetRadius( radius );/*由于在 ConnectedThresholdImageFilter 中，现在我们就必须提供在区域中能被输出像素所接受的亮度值而且至少是一个种子点来定义最初的区域*/InternalImageType::IndexType index;index[0] = atoi( "107" );index[1] = atoi( "69");neighborhoodConnected->SetSeed( index );neighborhoodConnected->SetReplaceValue( 255 );/*Writer 上的 Updata() 方法引发了管道的运行。通常在出现错误和抛出异议时, 从一个try / catch 模块调用 updata*/try{writer->Update();}catch( itk::ExceptionObject & excep ){std::cerr << "Exception caught !" << std::endl;std::cerr << excep << std::endl;}return EXIT_SUCCESS;
}

#include "itkConfidenceConnectedImageFilter.h"//包含置信连接类
//图像中存在的噪声会降低这个滤波器生长大面积区域的能力。当面对噪声图像时，通常
//是使用一个边缘保留平滑滤波器
#include "itkCastImageFilter.h"//滤波器
#include "itkCurvatureFlowImageFilter.h"//边缘保留平滑滤波器
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"//ConfidenceConnected Image Filter 使用的标准是基于当前区域的简单统计上的。首先，算法计算包含在区域中的所
//有像素亮度的平均值和标准差。用户提供一个因子用来乘以标准差并定义一个平均值的范
//围。相邻像素中亮度值在这个范围内的将包含到这个区域中。当没有更多的像素符合这个标
//准时，算法将结束它的第一次迭代。使用包含在区域内的所有像素再次计算亮度值的平均值
//和标准差。这个平均值和标准差定义一个新的亮度范围，用来查看当前区域的邻域并评价它
//们的亮度是否落在这个范围内。重复这个迭代过程直到没有新的像素再加进来或者已经达到
//了迭代器的最大数目。下面的式子阐述了这个滤波器的使用范围：
//I(X) ∈[m− fσ, m + fσ]
//其中 m 和 σ 分别是区域亮度的平均值和标准差， f 是用户提供的一个因子， I() 是图像，
//X 是区域中特殊相邻像素的位置int main( int argc, char *argv[] )
{//if( argc < 5 )//  {//  std::cerr << "Missing Parameters " << std::endl;//  std::cerr << "Usage: " << argv[0];//  std::cerr << " inputImage  outputImage seedX seedY " << std::endl;//  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;typedef itk::CastImageFilter< InternalImageType, OutputImageType >CastingFilterType;CastingFilterType::Pointer caster = CastingFilterType::New();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_connect2.png" );//使用图像类型作为模板参数来对平滑滤波器进行实例化typedef itk::CurvatureFlowImageFilter< InternalImageType, InternalImageType >CurvatureFlowImageFilterType;//调用 New( ) 方式来创建滤波器并将接指向 itk::SmartPointerCurvatureFlowImageFilterType::Pointer smoothing =CurvatureFlowImageFilterType::New();/*现在我们声明区域生长滤波器的类型，本例中使用 ConnectedThresholdImageFilter*/typedef itk::ConfidenceConnectedImageFilter<InternalImageType, InternalImageType> ConnectedFilterType;//然后我们使用 New( ) 方式构造这种类的一个滤波器ConnectedFilterType::Pointer confidenceConnected = ConnectedFilterType::New();/*现在到了创建一个简单的线性管道的时候。在管道的开头添加一个 reader 件头并在末尾添加一个 cast filter 和 writer 。由于只有仅仅一小部分图像文件格式支持浮点型数据类型，所以使用 cast filter 将浮点型数据类型转换成整型*/smoothing->SetInput( reader->GetOutput() );//滤波confidenceConnected->SetInput( smoothing->GetOutput() );//区域增长置信连接caster->SetInput( confidenceConnected->GetOutput() );//类型转换，浮点型数据类型转换成整型writer->SetInput( caster->GetOutput() );//滤波/*CurvatureFlowImageFilter 需要定义两个参数。下面是一个二维图像的常见值。当然它们也需要根据输入图像存在的噪声的数量进行适当的调整*/smoothing->SetNumberOfIterations( 5 );smoothing->SetTimeStep( 0.125 );//区域增长-置信连接/*ConfidenceConnectedImageFilter 需要定义两个参数。第一个，定义亮度范围大小的因子f 。乘法因子太小将限制当前区域中很多有相似亮度的像素；乘法因子太大将放松接受条件并导致区域过度生长。值太大就会使得这些区域生长到邻近区域，而这些邻近区域实际上可能属于独立的解剖结构*/confidenceConnected->SetMultiplier( 2.5 );//定义亮度范围大小的因子f/*迭代器的数目是基于被分割的解剖学结构亮度的均匀性之上指定的。均匀性高的区域仅需要一对迭代器。带有 ramp 效果的区域，如带有不同一领域的 MRI 图像，就需要更多的迭代器。实际上，精心选择乘法因子似乎比迭代器的数目更加重要。然而，务必紧记这个算法毫无疑问需要集中于一个稳定的区域。在这个算法上使用过多的迭代器将很有可能使这个区域吞没整个图像*/confidenceConnected->SetNumberOfIterations( 5 );//迭代器数目（迭代次数）/*这个滤波器的输出是一个二值图像，这个二值图像除了分割出的区域外到处都是零值像素。区域中的亮度值是由 SetReplaceValue() 方式来选择的*/confidenceConnected->SetReplaceValue( 255 );/*这个算法的实例化需要用户提供一个种子点。将这个点选在被分割的解剖学结构的典型区域是很便捷的。种子区域周围的一个小的邻域将用来计算包含在标准里的最初的平均值和标准差。种子是以一种 itk::Index 的形式传递给 SetSeed() 方式的*/InternalImageType::IndexType  index;//设置种子点坐标位置index[0] = atoi( "107" );index[1] = atoi( "69" );confidenceConnected->SetSeed( index );/*种子区域周围邻域最初的大小是使用 SetInitialNeighborhoodRadius() 方式来定义的。这个邻域将定义成一个拥有 2r + 1 个像素的 N（二维图像为2，三维为3） 维矩形区域，其中 r 是作为最初邻域范围的传递值*/confidenceConnected->SetInitialNeighborhoodRadius( 2 );//设置最初邻域范围的传递值r  即5像素矩阵/*writer 上的 Updata() 方法引发了管道的运行。通常在出现错误和抛出异议时, 从一个try / catch 模块调用 updata*/try{writer->Update();}catch( itk::ExceptionObject & excep ){std::cerr << "Exception caught !" << std::endl;std::cerr << excep << std::endl;}return EXIT_SUCCESS;
}

#include "itkConfidenceConnectedImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"int main( int argc, char *argv[] )
{//if( argc < 3 )//  {//  std::cerr << "Missing Parameters " << std::endl;//  std::cerr << "Usage: " << argv[0];//  std::cerr << " inputImage  outputImage " << std::endl;//  return EXIT_FAILURE;//  }typedef   float           InternalPixelType;const     unsigned int    Dimension = 3;typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;typedef unsigned char                            OutputPixelType;typedef itk::Image< OutputPixelType, Dimension > OutputImageType;typedef itk::CastImageFilter< InternalImageType, OutputImageType >CastingFilterType;CastingFilterType::Pointer caster = CastingFilterType::New();typedef  itk::ImageFileReader< InternalImageType > ReaderType;typedef  itk::ImageFileWriter<  OutputImageType  > WriterType;ReaderType::Pointer reader = ReaderType::New();WriterType::Pointer writer = WriterType::New();reader->SetFileName( "BrainProtonDensity3Slices.mha" );writer->SetFileName( "BrainProtonDensity3Slices_3.mha" );typedef itk::CurvatureFlowImageFilter< InternalImageType, InternalImageType >CurvatureFlowImageFilterType;CurvatureFlowImageFilterType::Pointer smoothing =CurvatureFlowImageFilterType::New();typedef itk::ConfidenceConnectedImageFilter<InternalImageType, InternalImageType>ConnectedFilterType;ConnectedFilterType::Pointer confidenceConnected = ConnectedFilterType::New();smoothing->SetInput( reader->GetOutput() );confidenceConnected->SetInput( smoothing->GetOutput() );caster->SetInput( confidenceConnected->GetOutput() );writer->SetInput( caster->GetOutput() );smoothing->SetNumberOfIterations( 2 );//smoothing->SetTimeStep(0.05);//每步迭代时间confidenceConnected->SetMultiplier( 2.5 );//；设置乘法因子f 2.5  可调confidenceConnected->SetNumberOfIterations( 5 );//设置迭代次数为5（迭代器数目）confidenceConnected->SetInitialNeighborhoodRadius( 2 );//设置领域范围为2confidenceConnected->SetReplaceValue( 255 );//设置种子点1InternalImageType::IndexType index1;index1[0] = 63;//X 轴index1[1] = 67;//Y 轴index1[2] = 1;//Z 轴 （第几个切片）confidenceConnected->AddSeed(index1);/* InternalImageType::IndexType index1;index1[0] = 118;index1[1] = 133;index1[2] = 92;confidenceConnected->AddSeed( index1 );*/try{writer->Update();}catch( itk::ExceptionObject & excep ){std::cerr << "Exception caught !" << std::endl;std::cerr << excep << std::endl;return EXIT_FAILURE;}return EXIT_SUCCESS;
}

#include "itkIsolatedConnectedImageFilter.h"
#include "itkImage.h"
#include "itkCastImageFilter.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"int main( int argc, char *argv[] )
{/*if( argc < 7 ){std::cerr << "Missing Parameters " << std::endl;std::cerr << "Usage: " << argv[0];std::cerr << " inputImage  outputImage seedX1 seedY1";std::cerr << " lowerThreshold seedX2 seedY2" << std::endl;return EXIT_FAILURE;}*///我们使用一个像素类型和一个特殊维来定义图像的类型：typedef   float           InternalPixelType;const     unsigned int    Dimension = 2;typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;//下面几行是对 IsolatedConnectedImageFilter 进行实例化的代码typedef unsigned char                            OutputPixelType;typedef itk::Image< OutputPixelType, Dimension > OutputImageType;typedef itk::CastImageFilter< InternalImageType, OutputImageType >CastingFilterType;CastingFilterType::Pointer caster = CastingFilterType::New();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( " Isolated_baizhi.png" );typedef itk::CurvatureFlowImageFilter< InternalImageType, InternalImageType >CurvatureFlowImageFilterType;CurvatureFlowImageFilterType::Pointer smoothing =CurvatureFlowImageFilterType::New();typedef itk::IsolatedConnectedImageFilter<InternalImageType,InternalImageType> ConnectedFilterType;//使用 New( ) 方式对这个类的一个文件进行结构化ConnectedFilterType::Pointer isolatedConnected = ConnectedFilterType::New();//现在连接管道smoothing->SetInput( reader->GetOutput() );isolatedConnected->SetInput( smoothing->GetOutput() );caster->SetInput( isolatedConnected->GetOutput() );writer->SetInput( caster->GetOutput() );/*IsolatedConnectedImageFilter 期望用户指定一个门限和两个种子。在这个例子中，我们从命令行得到它们*/smoothing->SetNumberOfIterations( 5 );smoothing->SetTimeStep( 0.125 );InternalImageType::IndexType  indexSeed1;//白质种子点indexSeed1[0] = atoi( "61" );indexSeed1[1] = atoi( "140" );//下门限值（要分割的白质的下门限值）const InternalPixelType lowerThreshold = atof( "150" );InternalImageType::IndexType  indexSeed2;//灰质种子点indexSeed2[0] = atoi( "63" );indexSeed2[1] = atoi( "43" );/*由于在 ConnectedThresholdImageFilter 中，现在我们就必须指定在区域中能被输出像素所接受的亮度值以及至少一个种子点来定义最初的区域*/isolatedConnected->SetLower(  lowerThreshold  );isolatedConnected->AddSeed1( indexSeed1 );isolatedConnected->AddSeed2( indexSeed2 );isolatedConnected->SetReplaceValue( 255 );/*writer 上的 Updata() 方法触发管道的运行。通常在出现错误和抛出异议时, 从一个 try / catch模块调用 updata*/try{writer->Update();}catch( itk::ExceptionObject & excep ){std::cerr << "Exception caught !" << std::endl;std::cerr << excep << std::endl;}/*这个亮度值允许我们对两个区域进行分割，使用 GetIsolatedValue() 方式可以对区域进行恢复*/std::cout << "Isolated Value Found = ";std::cout << isolatedConnected->GetIsolatedValue()  << std::endl;return EXIT_SUCCESS;
}

#include "itkImage.h"
#include "itkCastImageFilter.h"
#include "itkCurvatureFlowImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"int main(int argc, char* argv[])
{//定义浮点像素类型 InternalPixelTypetypedef   float           InternalPixelType;const     unsigned int    Dimension = 2;//定义浮点、维数2的图像类型 InternalImageTypetypedef itk::Image< InternalPixelType, Dimension >  InternalImageType;//定义无符号字符 输出像素类型 OutputPixelTypetypedef unsigned char                            OutputPixelType;//定义无符号字符、维数2 输出图像类型 OutputImageTypetypedef itk::Image< OutputPixelType, Dimension > OutputImageType;//定义滤波输入图像类型InternalImageType   滤波输出图像类型OutputImageTypetypedef itk::CastImageFilter< InternalImageType, OutputImageType >CastingFilterType;//实例化滤波器对象casterCastingFilterType::Pointer caster = CastingFilterType::New();//定义图像读取类型ReaderTypetypedef  itk::ImageFileReader< InternalImageType > ReaderType;//定义图像写类型WriterTypetypedef  itk::ImageFileWriter<  OutputImageType  > WriterType;ReaderType::Pointer reader = ReaderType::New();WriterType::Pointer writer = WriterType::New();reader->SetFileName("BrainProtonDensitySlice.png");writer->SetFileName(" lvbo_CastingFilter.png");//实例化滤波对象smoothingtypedef itk::CurvatureFlowImageFilter< InternalImageType, InternalImageType >CurvatureFlowImageFilterType;CurvatureFlowImageFilterType::Pointer smoothing =CurvatureFlowImageFilterType::New();现在连接管道smoothing->SetInput(reader->GetOutput());//输出滤波后的效果图caster->SetInput(smoothing->GetOutput());writer->SetInput(caster->GetOutput());/*IsolatedConnectedImageFilter 期望用户指定一个门限和两个种子。在这个例子中，我们从命令行得到它们*/smoothing->SetNumberOfIterations(10);smoothing->SetTimeStep(0.125);try{writer->Update();}catch (itk::ExceptionObject & excep){std::cerr << "Exception caught !" << std::endl;std::cerr << excep << std::endl;}return EXIT_SUCCESS;
}

#include "itkVectorConfidenceConnectedImageFilter.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRGBPixel.h"int main( int argc, char *argv[] )
{/*if( argc < 7 ){std::cerr << "Missing Parameters " << std::endl;std::cerr << "Usage: " << argv[0]<< " inputImage  outputImage"<< " seedX seedY"<< " multiplier iterations" << std::endl;return EXIT_FAILURE;}*//*现在我们使用一个像素类型和一个特殊的维来定义图像类型。由于需要使用平滑滤波器，所以这种情况下对像素使用浮点型数据类型*/const unsigned int Dimension = 2;typedef unsigned char                           PixelComponentType;typedef itk::RGBPixel< PixelComponentType >     InputPixelType;typedef itk::Image< InputPixelType, Dimension > InputImageType;typedef unsigned char                            OutputPixelType;typedef itk::Image< OutputPixelType, Dimension > OutputImageType;typedef itk::ImageFileReader< InputImageType >  ReaderType;typedef itk::ImageFileWriter< OutputImageType > WriterType;ReaderType::Pointer reader = ReaderType::New();WriterType::Pointer writer = WriterType::New();reader->SetFileName( "VisibleWomanEyeSlice.png" );writer->SetFileName( "VisibleWomanEyeSlice_Vitreo.png" );/*声明区域生长滤波器的类型，本例中使用 itk::VectorConfidenceConnectedImageFilter*/ typedef itk::VectorConfidenceConnectedImageFilter< InputImageType,OutputImageType > ConnectedFilterType;/*然后我们使用 New() 方式对一个滤波器对象结构化*/ConnectedFilterType::Pointer confidenceConnected= ConnectedFilterType::New();//然后我们创建一个简单的线性处理管confidenceConnected->SetInput( reader->GetOutput() );writer->SetInput( confidenceConnected->GetOutput() );//乘法因数fconst double multiplier = atof("3");//VectorConfidenceConnectedImageFilter 需要定义两个参数。第一个，定义亮度范围大小//的因子 f 。乘法因子太小将限制当前区域中很多有相似亮度的像素；乘法因子太大将放松接//受条件并将导致区域过度生长。值太大就会使得这些区域生长到邻近区域，而这些邻近区域//实际上可能属于独立的解剖结构。confidenceConnected->SetMultiplier( multiplier );/*迭代器的数目是基于被分割的解剖学结构亮度的均匀性之上指定的。均匀性高的区域仅需要一对迭代器。带有 ramp 效果的区域，如带有不同一领域的 MRI 图像，就需要更多的迭代器。实际上，精心选择乘法因子似乎比迭代器的数目更加重要。然而，务必记住这个算法需要集中于一个稳定的区域。在这个算法上使用过多的迭代器将很有可能是这个区域吞没整个图像*///迭代器数目const unsigned int iterations = atoi("1");confidenceConnected->SetNumberOfIterations( iterations );/*这个滤波器的输出是一个二值图像，这个二值图像除了分割出的区域外到处都是零值像素。区域中的亮度值是由 SetReplaceValue() 方式来选择的*/confidenceConnected->SetReplaceValue( 255 );InputImageType::IndexType  index;/*这个算法的实例化需要用户提供一个种子点。将这个点选在被分割的解剖学结构的典型区域是很便捷的。种子区域周围的一个小的邻域将用来计算包含在标准里的最初的平均值和标准差。种子以一种 itk::Index 的形式传递给 SetSeed() 方式*///种子位置index[0] = atoi("66");index[1] = atoi("66");confidenceConnected->SetSeed( index );/*种子区域周围邻域最初的大小是使用 SetInitialNeighborhoodRadius() 方式来定义的。这个邻域将定义成一个拥有 2r + 1 个像素的 N 维矩形区域，其中 r 是作为最初邻域范围的传递值*/confidenceConnected->SetInitialNeighborhoodRadius( 3 );/*writer 上的 Updata() 方法引发了管道的运行。通常在出现错误和抛出异议时, 从一个try / catch 模块调用 updata*/try{writer->Update();}catch( itk::ExceptionObject & excep ){std::cerr << "Exception caught !" << std::endl;std::cerr << excep << std::endl;}/*使用 GetMean() 方式和 GetCoaraiance() 方式可以查看最后一个迭代器最终的均值向量和协方差矩阵的值*/typedef ConnectedFilterType::MeanVectorType       MeanVectorType;typedef ConnectedFilterType::CovarianceMatrixType CovarianceMatrixType;const MeanVectorType & mean = confidenceConnected->GetMean();const CovarianceMatrixType & covariance= confidenceConnected->GetCovariance();std::cout << "Mean vector = "       << mean       << std::endl;std::cout << "Covariance matrix = " << covariance << std::endl;return EXIT_SUCCESS;
}

#include <iostream>
#include "itkVectorGradientAnisotropicDiffusionImageFilter.h"
#include "itkVectorGradientMagnitudeImageFilter.h"
#include "itkWatershedImageFilter.h"#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkVectorCastImageFilter.h"
#include "itkScalarToRGBPixelFunctor.h"int main( int argc, char *argv[] )
{/*if (argc < 8 ){std::cerr << "Missing Parameters " << std::endl;std::cerr << "Usage: " << argv[0];std::cerr << " inputImage outputImage conductanceTerm diffusionIterations lowerThreshold outputScaleLevel gradientMode " << std::endl;return EXIT_FAILURE;}*//*现在我们声明图像和像素类型用来势力化滤波器。所有的滤波器都需要实数类型的像素类型来正常工作。预处理阶段直接使用向量类型的数据而分割时使用浮点型标量数据。使用itk::VectorCastImageFilter 来将图像从 RGB 像素类型转换为数字向量类型*/typedef itk::RGBPixel< unsigned char >       RGBPixelType;typedef itk::Image< RGBPixelType, 2 >        RGBImageType;typedef itk::Vector< float, 3 >              VectorPixelType;typedef itk::Image< VectorPixelType, 2 >     VectorImageType;typedef itk::Image< itk::IdentifierType, 2 > LabeledImageType;typedef itk::Image< float, 2 >               ScalarImageType;//使用上面创建的类型来声明各种图像域处理滤波器，并最终将它们用在处理过程中typedef itk::ImageFileReader< RGBImageType >   FileReaderType;typedef itk::VectorCastImageFilter< RGBImageType, VectorImageType >CastFilterType;typedef itk::VectorGradientAnisotropicDiffusionImageFilter<VectorImageType, VectorImageType >DiffusionFilterType;typedef itk::VectorGradientMagnitudeImageFilter< VectorImageType >GradientMagnitudeFilterType;typedef itk::WatershedImageFilter< ScalarImageType >WatershedFilterType;typedef itk::ImageFileWriter<RGBImageType> FileWriterType;FileReaderType::Pointer reader = FileReaderType::New();reader->SetFileName("VisibleWomanEyeSlice.png");CastFilterType::Pointer caster = CastFilterType::New();/*接下来我们实例化这些滤波器并设置它们的参数。第一步在图像预处理管道中使用一个各向异性扩散滤波器来扩散输入彩色图像。对于这种滤波器类， CFL 条件要求对二维图像的time step 不能超过 0.25 ，对三维不能超过 0.125 。迭代器的数量和 conductance term 将从命令行得到。参见 6.7.3 小节将得到更多关于 ITK 各向异性扩散滤波器的信息*/DiffusionFilterType::Pointer diffusion = DiffusionFilterType::New();//设置迭代次数diffusion->SetNumberOfIterations( atoi("10") );//设置参数 conductancediffusion->SetConductanceParameter( atof("2.0") );diffusion->SetTimeStep(0.125);//对向量类型图像 ITK 梯度大小滤波器可以随意地选择几个参数，这里我们仅允许那些主//要成分分析的可用和不可用GradientMagnitudeFilterType::Pointergradient = GradientMagnitudeFilterType::New();//设置主要组成部分gradient->SetUsePrincipleComponents(atoi("on"));/*最后我们设置分水岭滤波器。它有两个参数。水平 Level 控制分水岭深度，而门限Threshold 控制输入的最低门限值。这两个参数都作为输入图像中最大深度的一个百分比(0.0 - 1.0) 来设置的*/WatershedFilterType::Pointer watershed = WatershedFilterType::New();//设置levelwatershed->SetLevel( atof("0.2") );/*分水岭滤波器的输出是一个无符号长整型 lable 图像，其中 lable 表示在一个特殊分割区域中一个像素的成员关系。这种形式对视觉是不现实的，所以为了这个例子的目的，我们将把它转换为 RGB 像素。 RGB 图像具有可以被保存为 png 文件和使用标准视图软件观看的优点。itk::Functor::ScalarToRGBPixelFunctor 类是可以将一个标量值转换成一个 itk::RGBPixel 的一个特殊功能对象。将这个算符写进 itk::UnaryFunctorImageFilter 来创建一个图像滤波器，用来把标量图像转换成 RGB 图像*///设置阈值watershed->SetThreshold( atof("0.01") );typedef itk::Functor::ScalarToRGBPixelFunctor<unsigned long>ColorMapFunctorType;typedef itk::UnaryFunctorImageFilter<LabeledImageType,RGBImageType, ColorMapFunctorType> ColorMapFilterType;ColorMapFilterType::Pointer colormapper = ColorMapFilterType::New();FileWriterType::Pointer writer = FileWriterType::New();//保存的文件名writer->SetFileName("VisibleWomanEyeSlice_FSL1.png");//这种滤波器连接到一个单一的管道，在每个结尾都有 readers 和 writerscaster->SetInput(reader->GetOutput());diffusion->SetInput(caster->GetOutput());gradient->SetInput(diffusion->GetOutput());watershed->SetInput(gradient->GetOutput());colormapper->SetInput(watershed->GetOutput());writer->SetInput(colormapper->GetOutput());try{writer->Update();}catch (itk::ExceptionObject &e){std::cerr << e << std::endl;return EXIT_FAILURE;}return EXIT_SUCCESS;
}

//包含用来从输入图像中去除噪声头文件
#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::SmartPointerSmoothingFilterType::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 而 α 的参考值是(K2 - 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");/*由于表示轮廓的起点将从头到尾不断传播，一旦到达一个特定的时间就希望停止这一过程。这就允许我们在假定相关区域已经计算过的假设下节省计算时间。使用 SetStoppingValue() 方式来定义停止过程的值。原则上，这个停止值应该比门限值高一点*/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() );speedWriter->SetFileName("FastMarchingFilterOutput3.mha");speedWriter->Update();InternalWriterType::Pointer gradientWriter = InternalWriterType::New();gradientWriter->SetInput( gradientMagnitude->GetOutput() );gradientWriter->SetFileName("FastMarchingFilterOutput2.mha");gradientWriter->Update();return EXIT_SUCCESS;
}

#include "itkGeodesicActiveContourLevelSetImageFilter.h"#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
#include "itkSigmoidImageFilter.h"
#include "itkFastMarchingImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
//程序第一阶段是使用 itk::CurvatureAnisotropicDiffusionImageFilter 来进行平
//滑。将平滑后的图像作为输入传递给 itk::GradientMagnitudeRecursiveGaussianImageFilter ，并
//且再传递给 itk::SigmoidImageFilter 以便产生潜在图像的边缘。一系列用户提供的种子传递给
//FastMarchingImageFilter 以便计算距离映射。从这个映射中减去一个常数以便得到一个水平
//集，其中零水平集表示最初的轮廓。这个水平集作为输入传递给 GeodesicActiveContourLevel
//SetImageFilter 。
//最后，将 GeodesicActiveContourLevelSetImageFilter 产生的水平集传递给一个 Binary
//ThresholdImageFilter 以便创建一个二值模板来表达分割对象。int main( int argc, char *argv[] )
{//现在我们使用一个特殊的像素类型和一个维来定义图像类型。由于需要用到平滑滤波//    器，这种情况下对像素使用浮点型数据类型typedef   float           InternalPixelType;const     unsigned int    Dimension = 2;typedef itk::Image< InternalPixelType, Dimension >  InternalImageType;//接下来几行对的类型进行实例化。并使用 New( ) 方式创建这个类型的一个对象typedef unsigned char                            OutputPixelType;typedef itk::Image< OutputPixelType, Dimension > OutputImageType;typedef itk::BinaryThresholdImageFilter<InternalImageType,OutputImageType    >       ThresholdingFilterType;ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();thresholder->SetLowerThreshold( -1000.0 );thresholder->SetUpperThreshold(     0.0 );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_G_ActiveContour.png");typedef itk::RescaleIntensityImageFilter<InternalImageType,OutputImageType >   CastFilterType;typedef   itk::CurvatureAnisotropicDiffusionImageFilter<InternalImageType,InternalImageType >  SmoothingFilterType;SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();typedef   itk::GradientMagnitudeRecursiveGaussianImageFilter<InternalImageType,InternalImageType >  GradientFilterType;typedef   itk::SigmoidImageFilter<InternalImageType,InternalImageType >  SigmoidFilterType;GradientFilterType::Pointer  gradientMagnitude = GradientFilterType::New();SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();sigmoid->SetOutputMinimum(  0.0  );sigmoid->SetOutputMaximum(  1.0  );typedef  itk::FastMarchingImageFilter<InternalImageType,InternalImageType >    FastMarchingFilterType;FastMarchingFilterType::Pointer  fastMarching = FastMarchingFilterType::New();typedef  itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType,InternalImageType >    GeodesicActiveContourFilterType;GeodesicActiveContourFilterType::Pointer geodesicActiveContour =GeodesicActiveContourFilterType::New();//设置扩展系数const double propagationScaling = atof("10.0" );/*对于 GeodesicActiveContourLevelSetImageFilter ，缩放比例参数用来对传播(膨胀) 系数、曲率(平滑) 系数和水平对流系数之间进行交替换位。这些参数使用 SetPropagationScaling() 、SetCurvatureScaling() 和 SetAdvectionScaling() 方式来设置。在这个例子中，我们将设置曲率和水平对流缩放比例作为一个参数并把传播缩放比例设置为一个命令行变量*///传播 ( 膨胀 ) 系数geodesicActiveContour->SetPropagationScaling( propagationScaling );//曲率(平滑) 系数geodesicActiveContour->SetCurvatureScaling( 1.0 );//水平对流系数geodesicActiveContour->SetAdvectionScaling( 1.0 );geodesicActiveContour->SetMaximumRMSError( 0.02 );geodesicActiveContour->SetNumberOfIterations( 800 );//现在使用接下来的几行代码将这些滤波器连接到图 9-21 所示的一个流程：smoothing->SetInput( reader->GetOutput() );gradientMagnitude->SetInput( smoothing->GetOutput() );sigmoid->SetInput( gradientMagnitude->GetOutput() );geodesicActiveContour->SetInput(  fastMarching->GetOutput() );geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );thresholder->SetInput( geodesicActiveContour->GetOutput() );writer->SetInput( thresholder->GetOutput() );smoothing->SetTimeStep( 0.125 );smoothing->SetNumberOfIterations(  5 );smoothing->SetConductanceParameter( 9.0 );//设置σconst double sigma = atof("0.5");gradientMagnitude->SetSigma(  sigma  );//设置αconst double alpha =  atof("-0.3");//设置β const double beta  =  atof("2.0");sigmoid->SetAlpha( alpha );sigmoid->SetBeta(  beta  );typedef FastMarchingFilterType::NodeContainer  NodeContainer;typedef FastMarchingFilterType::NodeType       NodeType;NodeContainer::Pointer seeds = NodeContainer::New();InternalImageType::IndexType  seedPosition;//设置种子点seedPosition[0] = atoi("40");seedPosition[1] = atoi("90");//设置间距const double initialDistance = atof("5.0");NodeType node;const double seedValue = - initialDistance;node.SetValue( seedValue );node.SetIndex( seedPosition );seeds->Initialize();seeds->InsertElement( 0, node );fastMarching->SetTrialPoints(  seeds  );fastMarching->SetSpeedConstant( 1.0 );CastFilterType::Pointer caster1 = CastFilterType::New();CastFilterType::Pointer caster2 = CastFilterType::New();CastFilterType::Pointer caster3 = CastFilterType::New();CastFilterType::Pointer caster4 = CastFilterType::New();WriterType::Pointer writer1 = WriterType::New();WriterType::Pointer writer2 = WriterType::New();WriterType::Pointer writer3 = WriterType::New();WriterType::Pointer writer4 = WriterType::New();caster1->SetInput( smoothing->GetOutput() );writer1->SetInput( caster1->GetOutput() );writer1->SetFileName("GeodesicActiveContourImageFilterOutput1.png");caster1->SetOutputMinimum(   0 );caster1->SetOutputMaximum( 255 );writer1->Update();caster2->SetInput( gradientMagnitude->GetOutput() );writer2->SetInput( caster2->GetOutput() );writer2->SetFileName("GeodesicActiveContourImageFilterOutput2.png");caster2->SetOutputMinimum(   0 );caster2->SetOutputMaximum( 255 );writer2->Update();caster3->SetInput( sigmoid->GetOutput() );writer3->SetInput( caster3->GetOutput() );writer3->SetFileName("GeodesicActiveContourImageFilterOutput3.png");caster3->SetOutputMinimum(   0 );caster3->SetOutputMaximum( 255 );writer3->Update();caster4->SetInput( fastMarching->GetOutput() );writer4->SetInput( caster4->GetOutput() );writer4->SetFileName("GeodesicActiveContourImageFilterOutput4.png");caster4->SetOutputMinimum(   0 );caster4->SetOutputMaximum( 255 );fastMarching->SetOutputSize(reader->GetOutput()->GetBufferedRegion().GetSize() );//Writer 上的 Updata( ) 方法触发流程的运行。通常在出现错误和抛出异议时，从一个//try / catch 模块调用 updata ：try{writer->Update();}catch( itk::ExceptionObject & excep ){std::cerr << "Exception caught !" << std::endl;std::cerr << excep << std::endl;return EXIT_FAILURE;}std::cout << std::endl;std::cout << "Max. no. iterations: " << geodesicActiveContour->GetNumberOfIterations() << std::endl;std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;std::cout << std::endl;std::cout << "No. elpased iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;writer4->Update();typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;InternalWriterType::Pointer mapWriter = InternalWriterType::New();mapWriter->SetInput( fastMarching->GetOutput() );mapWriter->SetFileName("GeodesicActiveContourImageFilterOutput4.mha");mapWriter->Update();InternalWriterType::Pointer speedWriter = InternalWriterType::New();speedWriter->SetInput( sigmoid->GetOutput() );speedWriter->SetFileName("GeodesicActiveContourImageFilterOutput3.mha");speedWriter->Update();InternalWriterType::Pointer gradientWriter = InternalWriterType::New();gradientWriter->SetInput( gradientMagnitude->GetOutput() );gradientWriter->SetFileName("GeodesicActiveContourImageFilterOutput2.mha");gradientWriter->Update();return EXIT_SUCCESS;
}