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 = 2;
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( "BrainProtonDensitySlice.png" );
51 //保存图像
52 writer->SetFileName( "naoshi_Threshold_LevelSet.png" );
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("81");
100 seedPosition[1] = atoi("112");
101 //表面到种子点的最初距离
102 const double initialDistance = atof("5");
103
104 NodeType node;
105
106 const double seedValue = - initialDistance;
107
108 node.SetValue( seedValue );
109 node.SetIndex( seedPosition );
110
111 seeds->Initialize();
112 seeds->InsertElement( 0, node );
113
114 fastMarching->SetTrialPoints( seeds );
115
116 fastMarching->SetSpeedConstant( 1.0 );
117 /*调用writer上的Updata()方法引发了管道的运行。通常在出现错误和抛出异议时，从一
118 个try / catch模块调用updata：*/
119 try
120 {
121 reader->Update();
122 const InternalImageType * inputImage = reader->GetOutput();
123 fastMarching->SetOutputRegion( inputImage->GetBufferedRegion() );
124 fastMarching->SetOutputSpacing( inputImage->GetSpacing() );
125 fastMarching->SetOutputOrigin( inputImage->GetOrigin() );
126 fastMarching->SetOutputDirection( inputImage->GetDirection() );
127 writer->Update();
128 }
129 catch( itk::ExceptionObject & excep )
130 {
131 std::cerr << "Exception caught !" << std::endl;
132 std::cerr << excep << std::endl;
133 return EXIT_FAILURE;
134 }
135
136 std::cout << std::endl;
137 std::cout << "Max. no. iterations: " << thresholdSegmentation->GetNumberOfIterations() << std::endl;
138 std::cout << "Max. RMS error: " << thresholdSegmentation->GetMaximumRMSError() << std::endl;
139 std::cout << std::endl;
140 std::cout << "No. elpased iterations: " << thresholdSegmentation->GetElapsedIterations() << std::endl;
141 std::cout << "RMS change: " << thresholdSegmentation->GetRMSChange() << std::endl;
142
143 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
144 //fastMarching快速步进输出水平集
145 InternalWriterType::Pointer mapWriter = InternalWriterType::New();
146 mapWriter->SetInput( fastMarching->GetOutput() );
147 mapWriter->SetFileName("fastMarchingImage.mha");
148 mapWriter->Update();
149 //阈值水平集分割的速度（传播系数P）图像
150 InternalWriterType::Pointer speedWriter = InternalWriterType::New();
151 speedWriter->SetInput( thresholdSegmentation->GetSpeedImage() );
152 speedWriter->SetFileName("speedTermImage.mha");
153 speedWriter->Update();
154
155 return EXIT_SUCCESS;
156 }
