第一个ITK程序
1、CMakeLists.txt
1 # This is the root ITK CMakeLists file.
2 cmake_minimum_required(VERSION 3.10)
3
4 # This project is designed to be built outside the Insight source tree.
5 project(ITK_demo)
6
7 # Find VTK
8 set(ITK_DIR D:/ProgramFiles/ITK-5.1.1/lib/cmake/ITK-5.1) # add wmz
9 find_package(ITK REQUIRED)
10 include_directories(${ITK_INCLUDE_DIRS})
11 message("ITK dir = ${ITK_INCLUDE_DIRS}")
12 message("ITK lib = ${ITK_LIBRARIES}")
13
14 include(${ITK_USE_FILE})
15 #aux_source_directory(src SRC_LIST)
16 set(SRC_LIST
17 ./itk_demo.cpp)
18
19 add_executable(itk_demo ${SRC_LIST} )
20 target_link_libraries(itk_demo ${ITK_LIBRARIES})
关于 include(${ITK_USE_FILE}) 的说明可以在 UseITK.cmake 中找到：
1 # -------------
2 #
3 # This file is not part of the ITK API. It exists purely as an
4 # implementation detail. This CMake module may change from version to
5 # version without notice, or even be removed.
6 #
7 # We mean it.
8 #
9
10 # This file sets up include directories, link directories, IO settings and
11 # compiler settings for a project to use ITK. It should not be
12 # included directly, but rather through the ITK_USE_FILE setting
13 # obtained from ITKConfig.cmake.
2、 测试数据
测试数据下载路径：https://github.com/InsightSoftwareConsortium/ITK/tree/master/Examples/Data
其实编译ITK时的目录下就有需要的测试数据，比如我的ITK-5.1.1目录下。
\ITK-5.1.1\Examples\Data
3、代码
作为第一个示例程序本来应该写一个很简单的像HelloWorld的程序，但是一些比较简单的官网的程序 要么依赖VTK，要么版本高于ITK5.1.1.
所以就找了一个比较长的程序，是一个配准的程序。
代码来自：https://github.com/InsightSoftwareConsortium/ITK/blob/master/Examples/RegistrationITKv4/MultiResImageRegistration1.cxx
我找了一个其他人做过的中文注释版
1 #include "itkImageRegistrationMethodv4.h"
2 #include "itkTranslationTransform.h"
3 #include "itkMeanSquaresImageToImageMetricv4.h"
4 #include "itkRegularStepGradientDescentOptimizerv4.h"
5 #include "itkImageFileReader.h"
6 #include "itkImageFileWriter.h"
7 #include "itkPNGImageIOFactory.h"
8 #include "itkResampleImageFilter.h"
9 #include "itkCastImageFilter.h"
10 #include "itkRescaleIntensityImageFilter.h"
11 #include "itkSubtractImageFilter.h"
12 /*****************************************************************************************************************
13 * 本例子是一个图像配准的Demo
14 * 0、创建了一个Command对象，用于监控配准的过程，被后面的对象调用
15 * 1、首先要定义像素的维度以及像素类型：进进而链接参考图像以及浮动图像
16 * 2、定义框架的基本组件：
17 * 确定变换种类：TransformType：二维变换
18 * 确定优化方法：OptimizerType：梯度下降
19 * 确定相似度度量：MetricType：链接两个图像：浮动图像以及参考图像
20 * 3、创建图像组件，并且通过创建上述框架，进而进行设置（链接）
21 * 4、设置插值方法：LinearInterpolateImageFunction并且链接在一起
22 * 5、6：通过ImageFileReader方法进行读取，链接到 registration并更新
23 * 7、针对前面的TransformType进行实例化：平移变换用于配准SetInitialTransformParameters：用于设置初始值
24 * 8、针对优化方法的设置：OptimizerType：前面在创建的时候已经设置了其梯度下降方法，此步骤用于对其微调：初始步长，收敛公差，最大迭代次数
25 * 9、通过RegistrationParameterScalesFromPhysicalShift：将每一个配准要素链接到配准方法中执行，
26 * 10、实例化Common对象，监控配准过程的执行，触发配准过程--迭代
27 * 11、通过update函数触发配准的执行
28 * 12、配准结果定义空间变换的参数序列：其结果由GetLastTransformParameters( )获得并且输出
29 * X、Y的变换：TranslationAlongX；TranslationAlongY
30 * 迭代次数：numberOfIterations
31 * 最后的结果：bestValue
32 * 通过CompositeTransform：AddTransform将转换添加到堆栈的背面，并且拥有可优化的参数。
33 * 也就是说：添加堆栈，副本？？
34 * 13、14、15、16、ResampleFilterType方法：
35 * 用变换参数将两幅图像进行叠加比较，并设置重采样滤波器：输入两幅图像
36 * 输出的是一个变换
37 * 对滤波器进行相关参数的设置：大小、原点、间距、位置
38 * 并通过CastFilterType：setInput：weiter进行相关的输出
39 * 此时：这个图象就是配准结束后的图像
40 * 17、通过itk::SubtractImageFilter对两幅图像进行比较：
41 * fixedImageReader；resampler
42 * 18、对图像进行处理 itk::RescaleIntensityImageFilter：调节一下亮度；并进行输出
43 * 19、一致性转发计算参考图像与正在移动图像之间的不同，输出图片5
44 ******************************************************************************************************************/
45 /*CommandIterationUpdate 类：
46 继承Command，监视配准过程的执行。每调用一次，输出相应参数
47 object类指向事件的观察者
48 Execute方法，类似cellbake，回转
49 observer方法：
50
51 */
52 class CommandIterationUpdate : public itk::Command
53 {
54 public:
55 typedef CommandIterationUpdate Self;
56 typedef itk::Command Superclass;
57 typedef itk::SmartPointer<Self> Pointer;
58 itkNewMacro(Self);//宏，包装了所有的new（）所有代码
59 protected:
60 CommandIterationUpdate() {};
61
62 public:
63
64 typedef itk::RegularStepGradientDescentOptimizerv4<double> OptimizerType;
65 typedef const OptimizerType* OptimizerPointer;
66
67 void Execute(itk::Object* caller, const itk::EventObject& event) ITK_OVERRIDE
68 {
69 Execute((const itk::Object*)caller, event);
70 }
71 //Object表示激活事件的对象，event表示被激活的事件
72 void Execute(const itk::Object* object, const itk::EventObject& event) ITK_OVERRIDE
73 {
74 OptimizerPointer optimizer = static_cast<OptimizerPointer>(object);
75 //checkEvent表示是否观察的对象
76 if (!itk::IterationEvent().CheckEvent(&event))
77 {
78 return;
79 }
80
81 std::cout << optimizer->GetCurrentIteration() << " = ";
82 std::cout << optimizer->GetValue() << " : ";
83 std::cout << optimizer->GetCurrentPosition() << std::endl;
84 }
85
86 };
87
88
89 int main()
90 {
91
92
93 //1、定义图像的维度以及像素执行
94 const unsigned int Dimension = 2;//定义维度
95 typedef float PixelType;//图像像素类型
96 typedef itk::Image< PixelType, Dimension > FixedImageType;//输入数据的类型：参考图像
97 typedef itk::Image< PixelType, Dimension > MovingImageType;//浮动图像
98
99
100
101 //2、定义配准框架的基本组件：变换、优化、测度配准组件
102 //用于2D图像的一个刚性配准，变换的唯一参数是：空间坐标类型
103 //配准
104 typedef itk::TranslationTransform< double, Dimension > TransformType;//把参考图像的空间映射到待配准图像的映射
105 //优化
106 typedef itk::RegularStepGradientDescentOptimizerv4<double> OptimizerType;//优化算法：牛顿梯度下降法
107 //度量
108 typedef itk::MeanSquaresImageToImageMetricv4<//相似度测量：均方根
109 FixedImageType,
110 MovingImageType > MetricType;
111
112
113
114 //3、该组件用用于连接其他组件
115 typedef itk::ImageRegistrationMethodv4<
116 FixedImageType,
117 MovingImageType,
118 TransformType > RegistrationType;
119
120 MetricType::Pointer metric = MetricType::New();
121 OptimizerType::Pointer optimizer = OptimizerType::New();
122 RegistrationType::Pointer registration = RegistrationType::New();
123 //连接组件：变换、优化组件
124 registration->SetMetric(metric);
125 registration->SetOptimizer(optimizer);
126
127 //4、插值方法
128 typedef itk::LinearInterpolateImageFunction<//选择校对机类型，校对机会对配准图像在非网格位置的程度进行评估
129 FixedImageType,
130 double > FixedLinearInterpolatorType;
131 typedef itk::LinearInterpolateImageFunction<
132 MovingImageType,
133 double > MovingLinearInterpolatorType;
134 FixedLinearInterpolatorType::Pointer fixedInterpolator =//每一个配准要素需要其new创建
135 FixedLinearInterpolatorType::New();
136 MovingLinearInterpolatorType::Pointer movingInterpolator =
137 MovingLinearInterpolatorType::New();
138 metric->SetFixedInterpolator(fixedInterpolator);
139 metric->SetMovingInterpolator(movingInterpolator);
140
141 //5、设置待配准图像以及变换区域
142 typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
143 typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
144 FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
145 MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
146
147
148 //6、读图像！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！！
149 itk::PNGImageIOFactory::RegisterOneFactory();
150 fixedImageReader->SetFileName("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\data\\BrainProtonDensitySliceBorder20.png");//输入图像文件
151 movingImageReader->SetFileName("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\data\\BrainProtonDensitySliceShifted13x17y.png");
152 //因为图像是从文件读取的，所以下面方法是用于获取图像数据的
153 //需要 itk::ImageRegistrationMethod 从 file readers 的输出获取输入
154 registration->SetFixedImage(fixedImageReader->GetOutput());
155 registration->SetMovingImage(movingImageReader->GetOutput());
156 //更新reader，确保其有效
157 // fixedImageReader->Update();
158 // movingImageReader->Update();
159
160
161 //7、平移变换用于配准
162 TransformType::Pointer movingInitialTransform = TransformType::New();
163 TransformType::ParametersType initialParameters(
164 movingInitialTransform->GetNumberOfParameters());
165 initialParameters[0] = 0.0; // Initial offset in mm along X
166 initialParameters[1] = 0.0; // Initial offset in mm along Y
167 movingInitialTransform->SetParameters(initialParameters);
168 registration->SetMovingInitialTransform(movingInitialTransform);
169
170 //8、准备执行配准方法：对优化器参数进行微调
171 TransformType::Pointer identityTransform = TransformType::New();
172 identityTransform->SetIdentity();
173
174 registration->SetFixedInitialTransform(identityTransform);
175 //初始振幅的长度用SetMaximumStepLength( ) 定义
176 //建立迭代的次数需要谨慎。最大数用SetNumberOfIterations()定义：
177 optimizer->SetLearningRate(4);
178 optimizer->SetMinimumStepLength(0.001);//优化器的收敛公差
179 optimizer->SetRelaxationFactor(0.5);
180
181 //9、将每一个配准要素连接到配准方法执行中
182 bool useEstimator = false;
183
184 //useEstimator = atoi(argv[6]) != 0;
185
186
187 if (useEstimator)
188 {
189
190
191 typedef itk::RegistrationParameterScalesFromPhysicalShift<MetricType> ScalesEstimatorType;
192 ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
193 scalesEstimator->SetMetric(metric);
194 scalesEstimator->SetTransformForward(true);
195 optimizer->SetScalesEstimator(scalesEstimator);
196 optimizer->SetDoEstimateLearningRateOnce(true);
197 }
198 optimizer->SetNumberOfIterations(200);//最大迭代次数
199
200
201 //10、实例化commend对象，监视配准过程的执行，并处触发配准过程
202 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
203 optimizer->AddObserver(itk::IterationEvent(), observer);
204
205
206 const unsigned int numberOfLevels = 1;
207
208 RegistrationType::ShrinkFactorsArrayType shrinkFactorsPerLevel;
209 shrinkFactorsPerLevel.SetSize(1);
210 shrinkFactorsPerLevel[0] = 1;
211
212 RegistrationType::SmoothingSigmasArrayType smoothingSigmasPerLevel;
213 smoothingSigmasPerLevel.SetSize(1);
214 smoothingSigmasPerLevel[0] = 0;
215
216 registration->SetNumberOfLevels(numberOfLevels);
217 registration->SetSmoothingSigmasPerLevel(smoothingSigmasPerLevel);
218 registration->SetShrinkFactorsPerLevel(shrinkFactorsPerLevel);
219
220
221
222 //11、通过调用Update函数触发配准执行
223 try
224 {
225 registration->Update();
226 std::cout << "Optimizer stop condition: "
227 << registration->GetOptimizer()->GetStopConditionDescription()
228 << std::endl;
229 }
230 catch (itk::ExceptionObject& err)
231 {
232 std::cerr << "ExceptionObject caught !" << std::endl;
233 std::cerr << err << std::endl;
234 return EXIT_FAILURE;
235 }
236
237 //12、配准结果是一系列定义空间变换的参数序列，结果由get获得
238 TransformType::ConstPointer transform = registration->GetTransform();
239 TransformType::ParametersType finalParameters = transform->GetParameters();
240 const double TranslationAlongX = finalParameters[0];//队列中每个元素对应着沿着一个空间维度的平移
241 const double TranslationAlongY = finalParameters[1];
242
243 //优化器能够询问抵达收敛的迭代的实际次数并通过GetCurrentIteration()返回出来
244 const unsigned int numberOfIterations = optimizer->GetCurrentIteration();//迭代次数
245 //最终参数集合的图像量规值通过优化器的GetValue();
246 const double bestValue = optimizer->GetValue();//最优化的度量
247
248 //将上述输出
249 std::cout << "Result = " << std::endl;
250 std::cout << " Translation X = " << TranslationAlongX << std::endl;//输出移动X的值
251 std::cout << " Translation Y = " << TranslationAlongY << std::endl;//输出移动Y的值
252 std::cout << " Iterations = " << numberOfIterations << std::endl;//输出迭代次数
253 std::cout << " Metric value = " << bestValue << std::endl;//输出优化的度量
254
255 typedef itk::CompositeTransform<
256 double,
257 Dimension > CompositeTransformType;
258 CompositeTransformType::Pointer outputCompositeTransform =
259 CompositeTransformType::New();
260 outputCompositeTransform->AddTransform(movingInitialTransform);
261 outputCompositeTransform->AddTransform(
262 registration->GetModifiableTransform());
263
264
265 //13、用变换结果将待配准图映射到参考图像中
266 typedef itk::ResampleImageFilter<
267 MovingImageType,
268 FixedImageType > ResampleFilterType;
269
270 //14、创建一个重采样滤波器，输入待配准图像
271 ResampleFilterType::Pointer resampler = ResampleFilterType::New();
272 resampler->SetInput(movingImageReader->GetOutput());
273 //配准函数生成的变换也作为重采样滤波器的输入被传递
274 resampler->SetTransform(outputCompositeTransform);
275
276 //15、ResampleImageFilter要求指定额外的参数，特别是输出图像的间 距、原点和大小
277 FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
278 resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());//尺寸
279 resampler->SetOutputOrigin(fixedImage->GetOrigin());//原点
280 resampler->SetOutputSpacing(fixedImage->GetSpacing());//间距
281 resampler->SetOutputDirection(fixedImage->GetDirection());//位置
282 resampler->SetDefaultPixelValue(100);
283
284 //16、滤波器的输出被传递给一个在文件中存储图像的writer
285 typedef unsigned char OutputPixelType;
286
287 typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
288
289 typedef itk::CastImageFilter<//转化重采样的像素类型到最终的writer类型
290 FixedImageType,
291 OutputImageType > CastFilterType;
292 typedef itk::ImageFileWriter< OutputImageType > WriterType;
293 //调用new函数创建新的滤波器
294 WriterType::Pointer writer = WriterType::New();
295 CastFilterType::Pointer caster = CastFilterType::New();
296 writer->SetFileName("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\output\\RegistrationITKv4Moving13x17yInputType.png");//写到文件夹位置
297 caster->SetInput(resampler->GetOutput());
298 writer->SetInput(caster->GetOutput());
299 writer->Update();//触发更新
300
301 //17、参照图像和被变换的待配准图像很容易用itk::SubtractImageFilter比较
302 //pixel-wise滤波器 计算两幅输入的同源像素的不同：
303 typedef itk::SubtractImageFilter<
304 FixedImageType,
305 FixedImageType,
306 FixedImageType > DifferenceFilterType;
307
308 DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
309
310 difference->SetInput1(fixedImageReader->GetOutput());//不同
311 difference->SetInput2(resampler->GetOutput());
312
313 //18、两幅图像的不同也许比较暗，我们用下面方法对其进行调节亮度，使之更加的明显
314 typedef itk::RescaleIntensityImageFilter<
315 FixedImageType,
316 OutputImageType > RescalerType;
317
318 RescalerType::Pointer intensityRescaler = RescalerType::New();
319
320 intensityRescaler->SetInput(difference->GetOutput());
321 intensityRescaler->SetOutputMinimum(0);
322 intensityRescaler->SetOutputMaximum(255);
323
324 resampler->SetDefaultPixelValue(1);
325
326
327 //输出到另外一个位置（调亮）
328 WriterType::Pointer writer2 = WriterType::New();
329 writer2->SetInput(intensityRescaler->GetOutput());
330
331 writer2->SetFileName("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\output\\Moving13x17yInputType.png");
332 writer2->Update();
333
334 //设置了一致性转换，计算参考图像的不同
335
336 resampler->SetTransform(identityTransform);
337
338 writer2->SetFileName("E:\\documents\\vs2019\\itk_demo\\build\\RelWithDebInfo\\output\\DifferenceBeforeRegistration.png");
339 writer2->Update();
340
341 return EXIT_SUCCESS;
342
343 }
4、结果
输入图像
输出图像
