1 //包含用来从输入图像中去除噪声头文件
2 #include "itkCurvatureAnisotropicDiffusionImageFilter.h"
3 //这两个滤波器连在一起将产生调节描述水平集运动的微分方程中的速率系数的图像潜能。
4 #include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
5 #include "itkSigmoidImageFilter.h"
6
7 #include "itkFastMarchingImageFilter.h"
8 //从 FastMarchingImageFilter 得到的时间交叉图将使用 BinaryThresholdImageFilter 来进行门限处理
9 #include "itkBinaryThresholdImageFilter.h"
10 //使用 itk::ImageFileReader 和 itk::ImageFileWriter 来对图像进行读写
11 #include "itkImageFileReader.h"
12 #include "itkImageFileWriter.h"
13 #include "itkRescaleIntensityImageFilter.h"
14
15 static void PrintCommandLineUsage( const int argc, const char * const argv[] )
16 {
17 std::cerr << "Missing Parameters " << std::endl;
18 std::cerr << "Usage: " << argv[0];
19 std::cerr << " inputImage outputImage seedX seedY";
20 std::cerr << " Sigma SigmoidAlpha SigmoidBeta TimeThreshold StoppingValue";
21 std::cerr << " smoothingOutputImage gradientMagnitudeOutputImage sigmoidOutputImage" << std::endl;
22
23 for (int qq=0; qq< argc; ++qq)
24 {
25 std::cout << "argv[" << qq << "] = " << argv[qq] << std::endl;
26 }
27 }
28
29 int main( int argc, char *argv[] )
30 {
31 /*if (argc != 13)
32 {
33 PrintCommandLineUsage(argc, argv);
34 return EXIT_FAILURE;
35 }*/
36
37 /*现在我们使用一个像素类型和一个特殊维来定义图像类型。由于需要用到平滑滤波器，
38 所以在这种情况下对像素使用浮点型数据类型*/
39 typedef float InternalPixelType;
40 const unsigned int Dimension = 2;
41 typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
42 //另一方面，输出图像被声明为二值的
43 typedef unsigned char OutputPixelType;
44 typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
45 //下面使用图像内在的类型和输出图像类型来对 BinaryThresholdImageFilter 的类型进行实
46 //例化：
47 typedef itk::BinaryThresholdImageFilter< InternalImageType,
48 OutputImageType > ThresholdingFilterType;
49 ThresholdingFilterType::Pointer thresholder = ThresholdingFilterType::New();
50 //设置阈值
51 const InternalPixelType timeThreshold = atof("100");
52 /*传递给 BinaryThresholdImageFilter 的上门限将定义那个从时间交叉图得来的时间快照。
53 在一个理想的应用中，用户拥有使用视觉反馈交互式的选择这个门限的能力。在这里，由于
54 是一个最小限度的例子，从命令行来得到这个值*/
55 thresholder->SetLowerThreshold( 0.0 );
56 thresholder->SetUpperThreshold( timeThreshold );
57
58 thresholder->SetOutsideValue( 0 );
59 thresholder->SetInsideValue( 255 );
60 //下面几行我们对读写器类型进行实例化
61 typedef itk::ImageFileReader< InternalImageType > ReaderType;
62 typedef itk::ImageFileWriter< OutputImageType > WriterType;
63
64 ReaderType::Pointer reader = ReaderType::New();
65 WriterType::Pointer writer = WriterType::New();
66 //输入图像
67 reader->SetFileName("BrainProtonDensitySlice.png");
68 //输出图像
69 writer->SetFileName("BrainProtonDensitySlice_FastMarching.png");
70
71 typedef itk::RescaleIntensityImageFilter<
72 InternalImageType,
73 OutputImageType > CastFilterType;
74 //使用图像内在的类型对 CurvatureAnisotropicDiffusionImageFilter 类型进行实例化
75 typedef itk::CurvatureAnisotropicDiffusionImageFilter<
76 InternalImageType,
77 InternalImageType > SmoothingFilterType;
78 //然后，通过调用 New( ) 方式来创建滤波器并将结果指向一个 itk::SmartPointer
79 SmoothingFilterType::Pointer smoothing = SmoothingFilterType::New();
80 /*使用图像内在的类型来对 GradientMagnitudeRecursiveGaussianImageFilter 和
81 SigmoidImageFilter 的类型进行实例化*/
82 typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
83 InternalImageType,
84 InternalImageType > GradientFilterType;
85 typedef itk::SigmoidImageFilter<
86 InternalImageType,
87 InternalImageType > SigmoidFilterType;
88 //使用 New( ) 方式来对相应的滤波器进行实例化
89 GradientFilterType::Pointer gradientMagnitude = GradientFilterType::New();
90 SigmoidFilterType::Pointer sigmoid = SigmoidFilterType::New();
91 /*使用 SetOutputMinimum() 和 SetOutputMaximum() 方式来定义 SigmoidImageFilter 输出的
92 最小值和最大值。在我们的例子中，我们希望这两个值分别是 0.0 和 1.0 ，以便得到一个最佳
93 的速率图像来给 MarchingImageFilter 。在 6.3.2 小节详细地表达了 SigmoidImageFilter 的使用方
94 法。*/
95 sigmoid->SetOutputMinimum( 0.0 );
96 sigmoid->SetOutputMaximum( 1.0 );
97 //现在我们声明 FastMarchingImageFilter 的类型：
98 typedef itk::FastMarchingImageFilter< InternalImageType,
99 InternalImageType > FastMarchingFilterType;
100 //然后我们使用 New( ) 方式结构化这个类的一个滤波器
101 FastMarchingFilterType::Pointer fastMarching
102 = FastMarchingFilterType::New();
103 //现在我们使用下面的代码在一个管道中连接这些滤波器
104 smoothing->SetInput( reader->GetOutput() );
105 gradientMagnitude->SetInput( smoothing->GetOutput() );
106 sigmoid->SetInput( gradientMagnitude->GetOutput() );
107 fastMarching->SetInput( sigmoid->GetOutput() );
108 thresholder->SetInput( fastMarching->GetOutput() );
109 writer->SetInput( thresholder->GetOutput() );
110 /*CurvatureAnisotropicDiffusionImageFilter 类需要定义一对参数。下面是二维图像中的典
111 型值。然而它们可能需要根据输入图像中存在的噪声的数量进行适当的调整。在 6.7.3 小节中
112 讨论了这个滤波器*/
113 smoothing->SetTimeStep( 0.125 );
114 smoothing->SetNumberOfIterations( 5 );
115 smoothing->SetConductanceParameter( 9.0 );
116 //设置Sigma(σ )
117 const double sigma = atof("1.0");
118 /*执行 GradientMagnitudeRecursiveGaussianImageFilter 就等同于通过一个派生的操作符紧
119 跟一个高斯核的一个回旋。这个高斯的 sigma(σ) 可以用来控制图像边缘影响的范围。在 6.4.2
120 小节中讨论了这个滤波器*/
121 gradientMagnitude->SetSigma( sigma );
122 //设置SigmoidAlpha(α) SigmoidBeta(β)
123 const double alpha = atof("-0.5");
124 const double beta = atof("3.0");
125 /*SigmoidImageFilter 类需要两个参数来定义用于 sigmoid 变量的线性变换。使用 SetAlpha()
126 和 SetBeta() 方式来传递这些参数。在这个例子的背景中，参数用来强化速率图像中的高低值
127 区域之间的区别。在理想情况下，解剖结构中的均匀区域内的速率应该是 1.0 而在结构边缘
128 附近应该迅速地衰减到 0.0 。下面是寻找这个值的启发。在梯度量值图像中，我们将沿被分
129 割的解剖结构地轮廓的最小值称为 K1 ，将结构中间的梯度强度的均值称为 K2 。这两个值表
130 达了我们想要映射到速率图像中间距为[0:1] 的动态范围。我们期望 sigmoid 将 K1 映射为 0.0 ，
131 K2 为 1.0 。由于 K1 的值要比 K2 更大而我们期望将它们分别映射为 0.0 和 1.0 ，所以我们对 alpha(α)
132 选择一个负值以便 sigmoid 函数也做一个反转亮度映射。这个映射将产生一个速率图像，其
133 中水平集将在均匀区域快速行进而在轮廓上停止。 Beta(β) 的参考值是(K1 + K2) / 2 而 α 的参考值
134 是(K2 - K1) / 6 必须是一个负数。在我们的简单例子中，这些值是用户从命令行提供的。用户
135 可以通过留心梯度量值图像来估计这些值*/
136 sigmoid->SetAlpha( alpha );
137 sigmoid->SetBeta( beta );
138 /*FastMarchingImageFilter 需要用户提供一个轮廓扩张的种子点。用户实际上不仅要提供
139 一个种子点而是一个种子点集。一个好的种子点集将增大分割一个复杂对象而不丢失数据的
140 可能性。使用多种子也会减少访问整个对象所需要的时间同时减少前面访问的区域边缘的漏
141 洞的风险。例如，当分割一个拉长的对象时，在其中一端设置一个种子将导致传播到另一端
142 会需要很长一段时间。沿着轴线设置几个种子无疑将是确定整个对象尽早被捕获的最佳策
143 略。水平集的一个重要的性质就是它们不需要任何额外辅助来融合几个起点的自然能力。使
144 用多种子正是利用了这个性质。
145 这些种子储存在一个容器中。在 FastMarchingImageFilter 的特性中将这个容器的类型定
146 义为 NodeContainer ：*/
147 typedef FastMarchingFilterType::NodeContainer NodeContainer;
148 typedef FastMarchingFilterType::NodeType NodeType;
149 NodeContainer::Pointer seeds = NodeContainer::New();
150
151 InternalImageType::IndexType seedPosition;
152 //设置种子点位置
153 seedPosition[0] = atoi("99");
154 seedPosition[1] = atoi("114");
155 //节点作为堆栈变量来创建，并使用一个值和一个 itk::Index 位置来进行初始化
156 NodeType node;
157 const double seedValue = 0.0;
158
159 node.SetValue( seedValue );
160 node.SetIndex( seedPosition );
161 //节点列表被初始化，然后使用 InsertElement( ) 来插入每个节点：
162 seeds->Initialize();
163 seeds->InsertElement( 0, node );
164 //现在使用 SetTrialPoints( ) 方式将种子节点集传递给 FastMarchingImageFilter ：
165 fastMarching->SetTrialPoints( seeds );
166 /*FastMarchingImageFilter 需要用户指定作为输出产生的图像的大小。使用 SetOutputSize()
167 来完成。注意：这里的大小是从平滑滤波器的输出图像得到的。这个图像的大小仅仅在直接
168 或间接地调用了这个滤波器的 Updata() 之后才是有效的。*/
169 //writer 上的 Updata( ) 方法引发了管道的运行。通常在出现错误和抛出异议时 , 从一个
170 //smoothingOutputImage
171 //使用边缘保护平滑滤波器平滑滤波后的图像
172 try
173 {
174 CastFilterType::Pointer caster1 = CastFilterType::New();
175 WriterType::Pointer writer1 = WriterType::New();
176 caster1->SetInput( smoothing->GetOutput() );
177 writer1->SetInput( caster1->GetOutput() );
178 //
179 writer1->SetFileName("smoothingOutputImage.png");
180 caster1->SetOutputMinimum( 0 );
181 caster1->SetOutputMaximum( 255 );
182 writer1->Update();
183 }
184 catch( itk::ExceptionObject & err )
185 {
186 std::cerr << "ExceptionObject caught !" << std::endl;
187 std::cerr << err << std::endl;
188 return EXIT_FAILURE;
189 }
190 // gradientMagnitudeOutputImage
191 //滤波后图像的梯度强度图像
192 try
193 {
194 CastFilterType::Pointer caster2 = CastFilterType::New();
195 WriterType::Pointer writer2 = WriterType::New();
196 caster2->SetInput( gradientMagnitude->GetOutput() );
197 writer2->SetInput( caster2->GetOutput() );
198 writer2->SetFileName("gradientMagnitudeOutputImage.png");
199 caster2->SetOutputMinimum( 0 );
200 caster2->SetOutputMaximum( 255 );
201 writer2->Update();
202 }
203 catch( itk::ExceptionObject & err )
204 {
205 std::cerr << "ExceptionObject caught !" << std::endl;
206 std::cerr << err << std::endl;
207 return EXIT_FAILURE;
208 }
209 // sigmoidOutputImage
210 //梯度强度的 sigmoid 函数图像
211 try
212 {
213 CastFilterType::Pointer caster3 = CastFilterType::New();
214 WriterType::Pointer writer3 = WriterType::New();
215 caster3->SetInput( sigmoid->GetOutput() );
216 writer3->SetInput( caster3->GetOutput() );
217 // FastMarchingImageFilter 分割处理产生的结果图像
218 writer3->SetFileName("sigmoidOutputImage.png");
219 caster3->SetOutputMinimum( 0 );
220 caster3->SetOutputMaximum( 255 );
221 writer3->Update();
222 }
223 catch( itk::ExceptionObject & err )
224 {
225 std::cerr << "ExceptionObject caught !" << std::endl;
226 std::cerr << err << std::endl;
227 return EXIT_FAILURE;
228 }
229
230
231 fastMarching->SetOutputSize(
232 reader->GetOutput()->GetBufferedRegion().GetSize() );
233 //设置stoppingTime，比阈值（门限值高一点）
234 const double stoppingTime = atof("103");
235 /*由于表示轮廓的起点将从头到尾不断传播，一旦到达一个特定的时间就希望停止这一过
236 程。这就允许我们在假定相关区域已经计算过的假设下节省计算时间。使用 SetStopping
237 Value() 方式来定义停止过程的值。原则上，这个停止值应该比门限值高一点*/
238 fastMarching->SetStoppingValue( stoppingTime );
239 // try / catch 模块调用 updata 。
240 try
241 {
242 writer->Update();
243 }
244 catch( itk::ExceptionObject & excep )
245 {
246 std::cerr << "Exception caught !" << std::endl;
247 std::cerr << excep << std::endl;
248 return EXIT_FAILURE;
249 }
250
251 try
252 {
253 CastFilterType::Pointer caster4 = CastFilterType::New();
254 WriterType::Pointer writer4 = WriterType::New();
255 caster4->SetInput( fastMarching->GetOutput() );
256 writer4->SetInput( caster4->GetOutput() );
257 writer4->SetFileName("FastMarchingFilterOutput4.png");
258 caster4->SetOutputMinimum( 0 );
259 caster4->SetOutputMaximum( 255 );
260 writer4->Update();
261 }
262 catch( itk::ExceptionObject & err )
263 {
264 std::cerr << "ExceptionObject caught !" << std::endl;
265 std::cerr << err << std::endl;
266 return EXIT_FAILURE;
267 }
268
269 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
270
271 InternalWriterType::Pointer mapWriter = InternalWriterType::New();
272 mapWriter->SetInput( fastMarching->GetOutput() );
273 mapWriter->SetFileName("FastMarchingFilterOutput4.mha");
274 mapWriter->Update();
275
276 InternalWriterType::Pointer speedWriter = InternalWriterType::New();
277 speedWriter->SetInput( sigmoid->GetOutput() );
278 speedWriter->SetFileName("FastMarchingFilterOutput3.mha");
279 speedWriter->Update();
280
281 InternalWriterType::Pointer gradientWriter = InternalWriterType::New();
282 gradientWriter->SetInput( gradientMagnitude->GetOutput() );
283 gradientWriter->SetFileName("FastMarchingFilterOutput2.mha");
284 gradientWriter->Update();
285
286 return EXIT_SUCCESS;
287 }
使用的终止值都为 100 。如图 9-16 所示表达了图 9-15 所示中阐述的流程的中间输出。从左到右分别是：各向异性扩散滤波器的输出，平滑后图像的梯度强度和最后作为 FastMarchingImageFilter 的速率图像使用的梯度强度的 sigmoid 。
基于 FastMarchingImageFilter 分割处理产生的结果：
注意：图像中的灰质部分并没有被完全分割，这就阐述了在对被分割的解剖结构并没有占据图像的延伸区域进行分割时水平集方法的弱点。这在当结构的宽度和梯度滤波器所产生的衰减带的大小做比较时显得尤为真实。对这个限制可行的一个工作区是沿着拉长的对象使用多种子分布，然而，白质相对于灰质分割来说是一个价值不高的工作，可能需要一个比这个简单例子中使用的更精细的方法。