1 #include "itkImage.h"
2 #include "itkThresholdSegmentationLevelSetImageFilter.h"
3
4 #include "itkFastMarchingImageFilter.h"
5 #include "itkBinaryThresholdImageFilter.h"
6 #include "itkImageFileReader.h"
7 #include "itkImageFileWriter.h"
8 #include "itkZeroCrossingImageFilter.h"
9
10
11 int main( int argc, char *argv[] )
12 {
13 /*if( argc < 8 )
14 {
15 std::cerr << "Missing Parameters " << std::endl;
16 std::cerr << "Usage: " << argv[0];
17 std::cerr << " inputImage outputImage";
18 std::cerr << " seedX seedY InitialDistance";
19 std::cerr << " LowerThreshold";
20 std::cerr << " UpperThreshold";
21 std::cerr << " [CurvatureScaling == 1.0]";
22 std::cerr << std::endl;
23 return EXIT_FAILURE;
24 }*/
25 /*现在我们使用一个特殊的像素类型和维来定义图像类型。在这种情况下我们将使用二维
26 浮点型图像*/
27 typedef float InternalPixelType;
28 const unsigned int Dimension = 3;
29 typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
30 //接下来的几行使用New()方式实例化一个ThresholdSegmentationLevelSetImageFilter：
31 typedef unsigned char OutputPixelType;
32 typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
33 typedef itk::BinaryThresholdImageFilter<InternalImageType, OutputImageType>
34 ThresholdingFilterType;
35
36 ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
37
38 thresholder->SetLowerThreshold( -1000.0 );
39 thresholder->SetUpperThreshold( 0.0 );
40
41 thresholder->SetOutsideValue( 0 );
42 thresholder->SetInsideValue( 255 );
43
44 typedef itk::ImageFileReader< InternalImageType > ReaderType;
45 typedef itk::ImageFileWriter< OutputImageType > WriterType;
46
47 ReaderType::Pointer reader = ReaderType::New();
48 WriterType::Pointer writer = WriterType::New();
49 //输入图像
50 reader->SetFileName( "BrainProtonDensity3Slices.mha" );
51 //保存图像
52 writer->SetFileName( "naoshi_Threshold_LevelSet.mha" );
53
54 typedef itk::FastMarchingImageFilter< InternalImageType, InternalImageType >
55 FastMarchingFilterType;
56
57 FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
58
59 typedef itk::ThresholdSegmentationLevelSetImageFilter< InternalImageType,
60 InternalImageType > ThresholdSegmentationLevelSetImageFilterType;
61 ThresholdSegmentationLevelSetImageFilterType::Pointer thresholdSegmentation =
62 ThresholdSegmentationLevelSetImageFilterType::New();
63 /*对于ThresholdSegmentationLevelSetImageFilter，缩放比例参数用来平衡从等式（9 - 3）
64 中的传播(膨胀)和曲率(表面平滑)系数的影响。这个滤波器中不使用水平对流系数。使用
65 SetPropagationScaling()和SetCurvatureScaling()方式来设置这些系数。在这个例子中这两个
66 系数都设置为1.0*/
67 thresholdSegmentation->SetPropagationScaling( 1.0 );
68 /*if ( argc > 8 )
69 {
70 thresholdSegmentation->SetCurvatureScaling( atof(argv[8]) );
71 }
72 else
73 {*/
74 thresholdSegmentation->SetCurvatureScaling( 1.0 );
75 /* }*/
76
77 thresholdSegmentation->SetMaximumRMSError( 0.02 );
78 thresholdSegmentation->SetNumberOfIterations( 1200 );//设置最大迭代次数
79 /* 收敛标准MaximumRMSError和MaximumIterations设置的和前面例子中的一样。现在我
80 们设置上下门限U和L，以及在初始化模型时使用的等值面值*/
81 //上门限值（U）
82 thresholdSegmentation->SetUpperThreshold( ::atof("250") );
83 //下门限值（L）
84 thresholdSegmentation->SetLowerThreshold( ::atof("210") );
85 thresholdSegmentation->SetIsoSurfaceValue(0.0);
86
87 thresholdSegmentation->SetInput( fastMarching->GetOutput() );
88 thresholdSegmentation->SetFeatureImage( reader->GetOutput() );
89 thresholder->SetInput( thresholdSegmentation->GetOutput() );
90 writer->SetInput( thresholder->GetOutput() );
91
92 typedef FastMarchingFilterType::NodeContainer NodeContainer;
93 typedef FastMarchingFilterType::NodeType NodeType;
94
95 NodeContainer::Pointer seeds = NodeContainer::New();
96
97 InternalImageType::IndexType seedPosition;
98 //设置种子点位置
99 seedPosition[0] = atoi("84");
100 seedPosition[1] = atoi("126");
101 seedPosition[2] = atoi("2");
102 //表面到种子点的最初距离
103 const double initialDistance = atof("5");//5
104
105 NodeType node;
106
107 const double seedValue = - initialDistance;
108
109 node.SetValue( seedValue );
110 node.SetIndex( seedPosition );
111
112 seeds->Initialize();
113 seeds->InsertElement( 0, node );
114
115 fastMarching->SetTrialPoints( seeds );
116
117 fastMarching->SetSpeedConstant( 1.0 );
118 /*调用writer上的Updata()方法引发了管道的运行。通常在出现错误和抛出异议时，从一
119 个try / catch模块调用updata：*/
120 try
121 {
122 reader->Update();
123 const InternalImageType * inputImage = reader->GetOutput();
124 fastMarching->SetOutputRegion( inputImage->GetBufferedRegion() );
125 fastMarching->SetOutputSpacing( inputImage->GetSpacing() );
126 fastMarching->SetOutputOrigin( inputImage->GetOrigin() );
127 fastMarching->SetOutputDirection( inputImage->GetDirection() );
128 writer->Update();
129 }
130 catch( itk::ExceptionObject & excep )
131 {
132 std::cerr << "Exception caught !" << std::endl;
133 std::cerr << excep << std::endl;
134 return EXIT_FAILURE;
135 }
136
137 std::cout << std::endl;
138 std::cout << "Max. no. iterations: " << thresholdSegmentation->GetNumberOfIterations() << std::endl;
139 std::cout << "Max. RMS error: " << thresholdSegmentation->GetMaximumRMSError() << std::endl;
140 std::cout << std::endl;
141 std::cout << "No. elpased iterations: " << thresholdSegmentation->GetElapsedIterations() << std::endl;
142 std::cout << "RMS change: " << thresholdSegmentation->GetRMSChange() << std::endl;
143
144 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
145 //fastMarching快速步进输出水平集
146 InternalWriterType::Pointer mapWriter = InternalWriterType::New();
147 mapWriter->SetInput( fastMarching->GetOutput() );
148 mapWriter->SetFileName("fastMarchingImage.mha");
149 mapWriter->Update();
150 //阈值水平集分割的速度（传播系数P）图像
151 InternalWriterType::Pointer speedWriter = InternalWriterType::New();
152 speedWriter->SetInput( thresholdSegmentation->GetSpeedImage() );
153 speedWriter->SetFileName("speedTermImage.mha");
154 speedWriter->Update();
155
156 return EXIT_SUCCESS;
157 }
