1 #include "itkGeodesicActiveContourLevelSetImageFilter.h"
2
3 #include "itkCurvatureAnisotropicDiffusionImageFilter.h"
4 #include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
5 #include "itkSigmoidImageFilter.h"
6 #include "itkFastMarchingImageFilter.h"
7 #include "itkRescaleIntensityImageFilter.h"
8 #include "itkBinaryThresholdImageFilter.h"
9 #include "itkImageFileReader.h"
10 #include "itkImageFileWriter.h"
11 //程序第一阶段是使用 itk::CurvatureAnisotropicDiffusionImageFilter 来进行平
12 //滑。将平滑后的图像作为输入传递给 itk::GradientMagnitudeRecursiveGaussianImageFilter ，并
13 //且再传递给 itk::SigmoidImageFilter 以便产生潜在图像的边缘。一系列用户提供的种子传递给
14 //FastMarchingImageFilter 以便计算距离映射。从这个映射中减去一个常数以便得到一个水平
15 //集，其中零水平集表示最初的轮廓。这个水平集作为输入传递给 GeodesicActiveContourLevel
16 //SetImageFilter 。
17 //最后，将 GeodesicActiveContourLevelSetImageFilter 产生的水平集传递给一个 Binary
18 //ThresholdImageFilter 以便创建一个二值模板来表达分割对象。
19
20 int main( int argc, char *argv[] )
21 {
22 //现在我们使用一个特殊的像素类型和一个维来定义图像类型。由于需要用到平滑滤波
23 // 器，这种情况下对像素使用浮点型数据类型
24 typedef float InternalPixelType;
25 const unsigned int Dimension = 2;
26 typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
27 //接下来几行对的类型进行实例化。并使用 New( ) 方式创建这个类型的一个对象
28 typedef unsigned char OutputPixelType;
29 typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
30 typedef itk::BinaryThresholdImageFilter<
31 InternalImageType,
32 OutputImageType > ThresholdingFilterType;
33
34 ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
35
36 thresholder->SetLowerThreshold( -1000.0 );
37 thresholder->SetUpperThreshold( 0.0 );
38
39 thresholder->SetOutsideValue( 0 );
40 thresholder->SetInsideValue( 255 );
41
42
43 typedef itk::ImageFileReader< InternalImageType > ReaderType;
44 typedef itk::ImageFileWriter< OutputImageType > WriterType;
45
46 ReaderType::Pointer reader = ReaderType::New();
47 WriterType::Pointer writer = WriterType::New();
48 //设置读取文件路径
49 reader->SetFileName("BrainProtonDensitySlice.png");
50 //设置保存文件路径
51 writer->SetFileName("BrainProtonDensitySlice_G_ActiveContour.png");
52
53 typedef itk::RescaleIntensityImageFilter<
54 InternalImageType,
55 OutputImageType > CastFilterType;
56
57 typedef itk::CurvatureAnisotropicDiffusionImageFilter<
58 InternalImageType,
59 InternalImageType > SmoothingFilterType;
60
61 SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
62
63
64 typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
65 InternalImageType,
66 InternalImageType > GradientFilterType;
67 typedef itk::SigmoidImageFilter<
68 InternalImageType,
69 InternalImageType > SigmoidFilterType;
70
71 GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
72
73 SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
74
75 sigmoid->SetOutputMinimum( 0.0 );
76 sigmoid->SetOutputMaximum( 1.0 );
77
78
79 typedef itk::FastMarchingImageFilter<
80 InternalImageType,
81 InternalImageType > FastMarchingFilterType;
82
83
84 FastMarchingFilterType::Pointer fastMarching = FastMarchingFilterType::New();
85
86 typedef itk::GeodesicActiveContourLevelSetImageFilter< InternalImageType,
87 InternalImageType > GeodesicActiveContourFilterType;
88 GeodesicActiveContourFilterType::Pointer geodesicActiveContour =
89 GeodesicActiveContourFilterType::New();
90 //设置扩展系数
91 const double propagationScaling = atof("10.0" );
92 /*对于 GeodesicActiveContourLevelSetImageFilter ，缩放比例参数用来对传播(膨胀) 系数、
93 曲率(平滑) 系数和水平对流系数之间进行交替换位。这些参数使用 SetPropagationScaling() 、
94 SetCurvatureScaling() 和 SetAdvectionScaling() 方式来设置。在这个例子中，我们将设置曲率
95 和水平对流缩放比例作为一个参数并把传播缩放比例设置为一个命令行变量*/
96 //传播 ( 膨胀 ) 系数
97 geodesicActiveContour->SetPropagationScaling( propagationScaling );
98 //曲率(平滑) 系数
99 geodesicActiveContour->SetCurvatureScaling( 1.0 );
100 //水平对流系数
101 geodesicActiveContour->SetAdvectionScaling( 1.0 );
102
103 geodesicActiveContour->SetMaximumRMSError( 0.02 );
104 geodesicActiveContour->SetNumberOfIterations( 800 );
105 //现在使用接下来的几行代码将这些滤波器连接到图 9-21 所示的一个流程：
106 smoothing->SetInput( reader->GetOutput() );
107 gradientMagnitude->SetInput( smoothing->GetOutput() );
108 sigmoid->SetInput( gradientMagnitude->GetOutput() );
109
110 geodesicActiveContour->SetInput( fastMarching->GetOutput() );
111 geodesicActiveContour->SetFeatureImage( sigmoid->GetOutput() );
112
113 thresholder->SetInput( geodesicActiveContour->GetOutput() );
114 writer->SetInput( thresholder->GetOutput() );
115
116 smoothing->SetTimeStep( 0.125 );
117 smoothing->SetNumberOfIterations( 5 );
118 smoothing->SetConductanceParameter( 9.0 );
119 //设置σ
120 const double sigma = atof("0.5");
121 gradientMagnitude->SetSigma( sigma );
122
123
124 //设置α
125 const double alpha = atof("-0.3");
126 //设置β
127 const double beta = atof("2.0");
128
129 sigmoid->SetAlpha( alpha );
130 sigmoid->SetBeta( beta );
131
132
133
134 typedef FastMarchingFilterType::NodeContainer NodeContainer;
135 typedef FastMarchingFilterType::NodeType NodeType;
136
137 NodeContainer::Pointer seeds = NodeContainer::New();
138
139 InternalImageType::IndexType seedPosition;
140 //设置种子点
141 seedPosition[0] = atoi("40");
142 seedPosition[1] = atoi("90");
143
144 //设置间距
145 const double initialDistance = atof("5.0");
146
147 NodeType node;
148
149 const double seedValue = - initialDistance;
150
151 node.SetValue( seedValue );
152 node.SetIndex( seedPosition );
153
154
155
156 seeds->Initialize();
157 seeds->InsertElement( 0, node );
158
159
160 fastMarching->SetTrialPoints( seeds );
161
162 fastMarching->SetSpeedConstant( 1.0 );
163
164
165 CastFilterType::Pointer caster1 = CastFilterType::New();
166 CastFilterType::Pointer caster2 = CastFilterType::New();
167 CastFilterType::Pointer caster3 = CastFilterType::New();
168 CastFilterType::Pointer caster4 = CastFilterType::New();
169
170 WriterType::Pointer writer1 = WriterType::New();
171 WriterType::Pointer writer2 = WriterType::New();
172 WriterType::Pointer writer3 = WriterType::New();
173 WriterType::Pointer writer4 = WriterType::New();
174
175 caster1->SetInput( smoothing->GetOutput() );
176 writer1->SetInput( caster1->GetOutput() );
177 writer1->SetFileName("GeodesicActiveContourImageFilterOutput1.png");
178 caster1->SetOutputMinimum( 0 );
179 caster1->SetOutputMaximum( 255 );
180 writer1->Update();
181
182 caster2->SetInput( gradientMagnitude->GetOutput() );
183 writer2->SetInput( caster2->GetOutput() );
184 writer2->SetFileName("GeodesicActiveContourImageFilterOutput2.png");
185 caster2->SetOutputMinimum( 0 );
186 caster2->SetOutputMaximum( 255 );
187 writer2->Update();
188
189 caster3->SetInput( sigmoid->GetOutput() );
190 writer3->SetInput( caster3->GetOutput() );
191 writer3->SetFileName("GeodesicActiveContourImageFilterOutput3.png");
192 caster3->SetOutputMinimum( 0 );
193 caster3->SetOutputMaximum( 255 );
194 writer3->Update();
195
196 caster4->SetInput( fastMarching->GetOutput() );
197 writer4->SetInput( caster4->GetOutput() );
198 writer4->SetFileName("GeodesicActiveContourImageFilterOutput4.png");
199 caster4->SetOutputMinimum( 0 );
200 caster4->SetOutputMaximum( 255 );
201
202 fastMarching->SetOutputSize(
203 reader->GetOutput()->GetBufferedRegion().GetSize() );
204
205 //Writer 上的 Updata( ) 方法触发流程的运行。通常在出现错误和抛出异议时，从一个
206 //try / catch 模块调用 updata ：
207 try
208 {
209 writer->Update();
210 }
211 catch( itk::ExceptionObject & excep )
212 {
213 std::cerr << "Exception caught !" << std::endl;
214 std::cerr << excep << std::endl;
215 return EXIT_FAILURE;
216 }
217
218 std::cout << std::endl;
219 std::cout << "Max. no. iterations: " << geodesicActiveContour->GetNumberOfIterations() << std::endl;
220 std::cout << "Max. RMS error: " << geodesicActiveContour->GetMaximumRMSError() << std::endl;
221 std::cout << std::endl;
222 std::cout << "No. elpased iterations: " << geodesicActiveContour->GetElapsedIterations() << std::endl;
223 std::cout << "RMS change: " << geodesicActiveContour->GetRMSChange() << std::endl;
224
225 writer4->Update();
226
227
228 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
229
230 InternalWriterType::Pointer mapWriter = InternalWriterType::New();
231 mapWriter->SetInput( fastMarching->GetOutput() );
232 mapWriter->SetFileName("GeodesicActiveContourImageFilterOutput4.mha");
233 mapWriter->Update();
234
235 InternalWriterType::Pointer speedWriter = InternalWriterType::New();
236 speedWriter->SetInput( sigmoid->GetOutput() );
237 speedWriter->SetFileName("GeodesicActiveContourImageFilterOutput3.mha");
238 speedWriter->Update();
239
240 InternalWriterType::Pointer gradientWriter = InternalWriterType::New();
241 gradientWriter->SetInput( gradientMagnitude->GetOutput() );
242 gradientWriter->SetFileName("GeodesicActiveContourImageFilterOutput2.mha");
243 gradientWriter->Update();
244
245 return EXIT_SUCCESS;
246 }
左脑室、右脑室、白质和灰质分别得到的输出窗口：
注意：分割白质部分需要一个相关的更大的传播缩放比例。有两个原因：白质边界的低对比和组织结构的复杂形状。不幸的是这些缩放比例参数的最优值仅仅通过实验来得到。在一个真实的应用中，我们可以可以想象一个通过用户管理轮廓运动的交互式机制从而调节这些参数。