Multi-Resolution 多分辨策略示例

/* 
*多分辨策略,或称为多分辨尺度, 是图像配准中广泛应用的一种方法 
*用于加快计算速度, 改进配准精度, 以及配准的健壮性 
*基本思想为: 首先使用一种粗糙的尺度对少量的图像像素进行处理,然后在下一层使用一种精确的 
*尺度, 并用上一层的结果对其参数进行初始化. 迭代该过程, 直到达到最精确的尺度. 
*这种由粗到细, 在大尺度上看整体, 在小尺度上看细节的方法能够极大程度地提高配准成功率 
*本例演示了使用:itk::MultiResolutionImageRegistrationMethod 解决一个简单的图像配准问题 
*采用多分辨策略进行图像配准, 需要两个图像金字塔组件, 用来创建向下采样的图像序列. 
*/

ITK 提供了一个兼容的多分辨策略配准框架,该框架有另外两个组件：固定图像及浮动图像金字塔，用于对固定图像及浮动图像向下采样。图像金字塔根据用户定义的收缩因子对图像进行平滑及二次采样处理。如下图所示：

多分辨策略通过由粗到精的方式解决配准问题，如下图所示。首先对固定及浮动图像金字塔的第一层进行处理，使用一种粗尺度方法，然后用其结果初始化下一层较精细的变换参数，循环，直至达到最精确的结果。

 

[cpp:showcolumns] view plaincopy

·········10········20········30········40········50········60········70········80········90········100·······110·······120·······130·······140·······150

1. #if defined(_MSC_VER)  
2. #pragma warning ( disable : 4786 )  
3. #endif  
4.   
5. #include "itkMultiResolutionImageRegistrationMethod.h"  
6. #include "itkTranslationTransform.h"  
7. #include "itkMattesMutualInformationImageToImageMetric.h"  
8. #include "itkLinearInterpolateImageFunction.h"  
9. #include "itkRegularStepGradientDescentOptimizer.h"  
10. #include "itkMultiResolutionPyramidImageFilter.h"  
11. #include "itkImage.h"  
12.   
13. #include "itkImageFileReader.h"  
14. #include "itkImageFileWriter.h"  
15.   
16. #include "itkResampleImageFilter.h"  
17. #include "itkCastImageFilter.h"  
18. #include "itkCheckerBoardImageFilter.h"  
19.   
20. #include "itkCommand.h"  
21.   
22. //  
23. //典型的情况, 用户可能会在不同的层间调整组件的参数设置, 或者使用不同的组件  
24. //这可以通过使用 ITK 实现的 Command/Observer 设计模式实现  
25. //在每一层开始配准前, MultiResolutionImageRegistrationMethod 会调用一个 IterationEvent  
26. //配准组件的改变, 可以通过实现响应该事件的 itk::Command 来实现  
27. //新接口的其模板参数为 multi-resolution registration type  
28. //本例通过改变优化器在不同层间的参数来演示该机制  
29. //  
30. template <typename TRegistration>  
31. class RegistrationInterfaceCommand : public itk::Command   
32. {  
33. public:  
34.     typedef  RegistrationInterfaceCommand   Self;  
35.     typedef  itk::Command                   Superclass;  
36.     typedef  itk::SmartPointer<Self>        Pointer;  
37.     itkNewMacro( Self );  
38. protected:  
39.     RegistrationInterfaceCommand() {};  
40.   
41. public:  
42.     typedef   TRegistration                              RegistrationType;  
43.     typedef   RegistrationType *                         RegistrationPointer;  
44.     typedef   itk::RegularStepGradientDescentOptimizer   OptimizerType;  
45.     typedef   OptimizerType *                            OptimizerPointer;  
46.   
47.     void Execute(itk::Object * object, const itk::EventObject & event)  
48.     {  
49.         if( !(itk::IterationEvent().CheckEvent( &event )) )  
50.         {  
51.             return;  
52.         }  
53.   
54.         RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );  
55.   
56.         //first resolution level 设置 step length, minimum step length 成较大值.  
57.         //然后在随后的每一层, 减小 minimum step length 10 倍, 以使优化器集中在较小的关键区域  
58.         //maximum step length 设置为当前步长, 这样, 优化器下一层配准前重新初始化时  
59.         //step length 将使用上一层的最后值.   
60.         //这将保证 the continuity of the path taken by the optimizer through the parameter space.  
61.         OptimizerPointer optimizer =   
62.                 dynamic_cast< OptimizerPointer >( registration->GetOptimizer() );  
63.   
64.         std::cout << "-------------------------------------" << std::endl;  
65.         std::cout << "MultiResolution Level : "  
66.                   << registration->GetCurrentLevel()  << std::endl;  
67.         std::cout << std::endl;  
68.   
69.         if ( registration->GetCurrentLevel() == 0 )  
70.         {  
71.             //第一层处理时设置的参数  
72.             optimizer->SetMaximumStepLength( 16.00 );    
73.             optimizer->SetMinimumStepLength( 0.01 );  
74.         }  
75.         else  
76.         {  
77.             //以后每进行到下一层, 则根据定义的收缩因子使用上一层的参数进行初始化  
78.             //这里分别缩小 4 倍, 10 倍.  
79.             optimizer->SetMaximumStepLength( optimizer->GetMaximumStepLength() / 4.0 );  
80.             optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() / 10.0 );  
81.         }  
82.     }  
83.   
84.     void Execute(const itk::Object * , const itk::EventObject & )  
85.     {  
86.         return;   
87.     }  
88. };  
89.   
90. class CommandIterationUpdate : public itk::Command   
91. {  
92. public:  
93.     typedef  CommandIterationUpdate   Self;  
94.     typedef  itk::Command             Superclass;  
95.     typedef  itk::SmartPointer<Self>  Pointer;  
96.     itkNewMacro( Self );  
97. protected:  
98.     CommandIterationUpdate() {};  
99. public:  
100.     typedef   itk::RegularStepGradientDescentOptimizer     OptimizerType;  
101.     typedef   const OptimizerType                         *OptimizerPointer;  
102.   
103.     void Execute(itk::Object *caller, const itk::EventObject & event)  
104.     {  
105.         Execute( (const itk::Object *)caller, event);  
106.     }  
107.   
108.     void Execute(const itk::Object * object, const itk::EventObject & event)  
109.     {  
110.         OptimizerPointer optimizer =   
111.                                 dynamic_cast< OptimizerPointer >( object );  
112.         if( !(itk::IterationEvent().CheckEvent( &event )) )  
113.         {  
114.             return;  
115.         }  
116.         std::cout << optimizer->GetCurrentIteration() << "   ";  
117.         std::cout << optimizer->GetValue() << "   ";  
118.         std::cout << optimizer->GetCurrentPosition() << std::endl;  
119.     }  
120. };  
121.   
122.   
123. int main( int argc, char *argv[] )  
124. {  
125.     if( argc < 4 )  
126.     {  
127.         std::cerr << "Missing Parameters " << std::endl;  
128.         std::cerr << "Usage: " << argv[0];  
129.         std::cerr << " fixedImageFile  movingImageFile ";  
130.         std::cerr << " outputImagefile [backgroundGrayLevel]";  
131.         std::cerr << " [checkerBoardBefore] [checkerBoardAfter]";  
132.         std::cerr << " [useExplicitPDFderivatives ] " << std::endl;  
133.         std::cerr << " [numberOfBins] [numberOfSamples ] " << std::endl;  
134.         return EXIT_FAILURE;  
135.     }  
136.   
137.   
138.     //1.定义图像类型  
139.     const    unsigned int    Dimension = 2;  
140.     typedef  unsigned short  PixelType;  
141.   
142.     typedef itk::Image< PixelType, Dimension >  FixedImageType;  
143.     typedef itk::Image< PixelType, Dimension >  MovingImageType;  
144.   
145.     //由于图像金字塔在计算降低采样率的图像时是递归的,所以输出图像要求 real pixel types:  
146.     typedef   float     InternalPixelType;  
147.     typedef itk::Image< InternalPixelType, Dimension > InternalImageType;  
148.   
149.   
150.     //2.定义配准框架的基本组件: 变换,优化,插值,测度, 及用于连接各组件的配准方法组件  
151.     //对于多分辨策略, 这里需要另外两个组件, 固定图像及浮动图像金字塔组件  
152.     typedef itk::TranslationTransform< double, Dimension > TransformType;   
153.     typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;   
154.     typedef itk::LinearInterpolateImageFunction<   
155.                             InternalImageType, double > InterpolatorType;   
156.     typedef itk::MattesMutualInformationImageToImageMetric<   
157.                        InternalImageType, InternalImageType > MetricType;   
158.   
159.     //多分辨策略配准方法组件,用于连接各组件  
160.     typedef itk::MultiResolutionImageRegistrationMethod<   
161.                 InternalImageType, InternalImageType >  RegistrationType;  
162.   
163.     //在多分辨策略框架中(multi-resolution framework)  
164.     //使用类 itk::MultiResolutionPyramidImageFilter  
165.     //创建一个降低采样率的图像(downsampled image)金字塔.   
166.     //每幅 downsampled image 的大小由用户以收缩因子的形式指定.   
167.     typedef itk::MultiResolutionPyramidImageFilter<  
168.             InternalImageType,InternalImageType >  FixedImagePyramidType;  
169.     typedef itk::MultiResolutionPyramidImageFilter<  
170.             InternalImageType,InternalImageType > MovingImagePyramidType;  
171.   
172.     TransformType::Pointer      transform     = TransformType::New();  
173.     OptimizerType::Pointer      optimizer     = OptimizerType::New();  
174.     InterpolatorType::Pointer   interpolator  = InterpolatorType::New();  
175.     RegistrationType::Pointer   registration  = RegistrationType::New();  
176.     MetricType::Pointer         metric        = MetricType::New();  
177.   
178.     FixedImagePyramidType::Pointer fixedImagePyramid =   
179.                                             FixedImagePyramidType::New();  
180.     MovingImagePyramidType::Pointer movingImagePyramid =  
181.                                            MovingImagePyramidType::New();  
182.   
183.     //3.使用配准方法组件,将 优化,变换,插值,测度 四个基本组件连接起来  
184.     //对于多分辨策略, 还需要连接 固定图像与浮动图像金字塔组件  
185.     registration->SetOptimizer(     optimizer     );  
186.     registration->SetTransform(     transform     );  
187.     registration->SetInterpolator(  interpolator  );  
188.     registration->SetMetric( metric  );  
189.     registration->SetFixedImagePyramid( fixedImagePyramid );  
190.     registration->SetMovingImagePyramid( movingImagePyramid );  
191.   
192.     typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;  
193.     typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;  
194.   
195.     FixedImageReaderType::Pointer  fixedImageReader    
196.                                             = FixedImageReaderType::New();  
197.     MovingImageReaderType::Pointer movingImageReader   
198.                                             = MovingImageReaderType::New();  
199.   
200.     fixedImageReader->SetFileName(  argv[1] );  
201.     movingImageReader->SetFileName( argv[2] );  
202.   
203.     //4.设置待配准图像及配准区域, 并对图像类型进行适当转换处理  
204.     //对输入图像进行额外的处理  
205.     //固定,浮动图像均来自文件,在将图像连接至配准过程前  
206.     //需要将它们转化为内部图像类型,使用 itk::CastImageFilter  
207.     typedef itk::CastImageFilter<   
208.                     FixedImageType, InternalImageType > FixedCastFilterType;  
209.     typedef itk::CastImageFilter<   
210.                   MovingImageType, InternalImageType > MovingCastFilterType;  
211.   
212.     FixedCastFilterType::Pointer fixedCaster   = FixedCastFilterType::New();  
213.     MovingCastFilterType::Pointer movingCaster = MovingCastFilterType::New();  
214.   
215.     fixedCaster->SetInput(  fixedImageReader->GetOutput() );  
216.     movingCaster->SetInput( movingImageReader->GetOutput() );  
217.   
218.     registration->SetFixedImage(    fixedCaster->GetOutput()    );  
219.     registration->SetMovingImage(   movingCaster->GetOutput()   );  
220.   
221.     fixedCaster->Update();  
222.     registration->SetFixedImageRegion(   
223.                         fixedCaster->GetOutput()->GetBufferedRegion() );  
224.   
225.   
226.     //5.设置各组件的参数，变换、优化,测度等组件的参数   
227.     typedef RegistrationType::ParametersType ParametersType;  
228.     ParametersType initialParameters( transform->GetNumberOfParameters() );  
229.   
230.     initialParameters[0] = 0.0;  // Initial offset in mm along X  
231.     initialParameters[1] = 0.0;  // Initial offset in mm along Y  
232.     registration->SetInitialTransformParameters( initialParameters );  
233.   
234.     //Matte 互信息测度需要选择两个参数:  
235.     //1. number of bins: 用来计算熵值  
236.     //2. number of spatial samples: 用来计算密度估计(density estimates)  
237.     //典型的情况 50 histogram bins 足够了,注意, bins 值对优化器行为有着动态的影响  
238.     //spatial samples 值依赖图像的内容, 如果图像比较光滑且没有包含大量的细节  
239.     //则可以使用总像素的 1% 作为其值. 否则, 可以使用更多的比例, 如 20%  
240.     //本例设置如下:  
241.     metric->SetNumberOfHistogramBins( 128 );  
242.     metric->SetNumberOfSpatialSamples( 50000 );  
243.   
244.     if( argc > 8 )  
245.     {  
246.         metric->SetNumberOfHistogramBins( atoi( argv[8] ) );  
247.     }  
248.   
249.     if( argc > 9 )  
250.     {  
251.         metric->SetNumberOfSpatialSamples( atoi( argv[9] ) );  
252.     }  
253.   
254.     //  
255.     metric->ReinitializeSeed( 76926294 );  
256.   
257.     if( argc > 7 )  
258.     {  
259.         metric->SetUseExplicitPDFDerivatives( atoi( argv[7] ) );  
260.     }  
261.       
262.     //设置迭代次数限制及松驰因子  
263.     //关于松驰因子, 需参考 regular step gradient descent 优化方法  
264.     //当该优化器遇到参数在参数空间运动的方向改变时, 便减小 step length  
265.     //缩小的比例通过松驰因子控制(relaxation factor). 该值默认为 0.5  
266.     optimizer->SetNumberOfIterations( 200 );  
267.     optimizer->SetRelaxationFactor( 0.9 );  
268.   
269.     //6.实例化一 Command/Observer 对象, 监视配准过程的执行, 并触发配准过程的执行.    
270.     //多分辨策略还需要实例化另一 Command/Observer 对象, 用于改变层间组件及参数  
271.     CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();  
272.     optimizer->AddObserver( itk::IterationEvent(), observer );  
273.   
274.     typedef RegistrationInterfaceCommand<RegistrationType> CommandType;  
275.     CommandType::Pointer command = CommandType::New();  
276.     registration->AddObserver( itk::IterationEvent(), command );  
277.   
278.     //对于多分辨策略, 需要设置其层数  
279.     //这里设置 multi-resolution levels 为 3, 并触发配准过程的执行.  
280.     registration->SetNumberOfLevels( 3 );  
281.   
282.     try   
283.     {   
284.         registration->StartRegistration();   
285.     }   
286.     catch( itk::ExceptionObject & err )   
287.     {   
288.         std::cout << "ExceptionObject caught !" << std::endl;   
289.         std::cout << err << std::endl;   
290.         return EXIT_FAILURE;  
291.     }   
292.   
293.     //获取配准完成后的各个参数的最终值  
294.     ParametersType finalParameters = registration->GetLastTransformParameters();  
295.     double TranslationAlongX = finalParameters[0];  
296.     double TranslationAlongY = finalParameters[1];  
297.     unsigned int numberOfIterations = optimizer->GetCurrentIteration();  
298.     double bestValue = optimizer->GetValue();  
299.   
300.     std::cout << "Result = " << std::endl;  
301.     std::cout << " Translation X = " << TranslationAlongX  << std::endl;  
302.     std::cout << " Translation Y = " << TranslationAlongY  << std::endl;  
303.     std::cout << " Iterations    = " << numberOfIterations << std::endl;  
304.     std::cout << " Metric value  = " << bestValue          << std::endl;  
305.   
306.     //7.重采样, 将变换后的浮动图像映射到固定图像空间中,保存配准结果  
307.     typedef itk::ResampleImageFilter<   
308.                     MovingImageType, FixedImageType > ResampleFilterType;  
309.   
310.     TransformType::Pointer finalTransform = TransformType::New();  
311.     finalTransform->SetParameters( finalParameters );  
312.     finalTransform->SetFixedParameters( transform->GetFixedParameters() );  
313.   
314.     ResampleFilterType::Pointer resample = ResampleFilterType::New();  
315.     resample->SetTransform( finalTransform );  
316.     resample->SetInput( movingImageReader->GetOutput() );  
317.   
318.     FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();  
319.   
320.     PixelType backgroundGrayLevel = 100;  
321.     if( argc > 4 )  
322.     {  
323.         backgroundGrayLevel = atoi( argv[4] );  
324.     }  
325.   
326.     resample->SetSize(    fixedImage->GetLargestPossibleRegion().GetSize() );  
327.     resample->SetOutputOrigin(  fixedImage->GetOrigin() );  
328.     resample->SetOutputSpacing( fixedImage->GetSpacing() );  
329.     resample->SetOutputDirection( fixedImage->GetDirection() );  
330.     resample->SetDefaultPixelValue( backgroundGrayLevel );  
331.   
332.     //转换输出图像类型 CastImageFilter  
333.     typedef  unsigned char  OutputPixelType;  
334.     typedef itk::Image< OutputPixelType, Dimension > OutputImageType;  
335.     typedef itk::CastImageFilter<   
336.                     FixedImageType, OutputImageType > CastFilterType;  
337.     typedef itk::ImageFileWriter< OutputImageType >  WriterType;  
338.   
339.     WriterType::Pointer      writer =  WriterType::New();  
340.     CastFilterType::Pointer  caster =  CastFilterType::New();  
341.   
342.     writer->SetFileName( argv[3] );  
343.     caster->SetInput( resample->GetOutput() );  
344.     writer->SetInput( caster->GetOutput()   );  
345.     writer->Update();  
346.   
347.     //生成配准前后, 浮动图像与固定图像的对比图, 使用一棋盘状图显示差异  
348.     typedef itk::CheckerBoardImageFilter< FixedImageType > CheckerBoardFilterType;  
349.     CheckerBoardFilterType::Pointer checker = CheckerBoardFilterType::New();  
350.   
351.     checker->SetInput1( fixedImage );  
352.     checker->SetInput2( resample->GetOutput() );  
353.   
354.     caster->SetInput( checker->GetOutput() );  
355.     writer->SetInput( caster->GetOutput()   );  
356.   
357.     resample->SetDefaultPixelValue( 0 );  
358.   
359.     //配准前浮动图像与固定图像间的差异  
360.     TransformType::Pointer identityTransform = TransformType::New();  
361.     identityTransform->SetIdentity();  
362.     resample->SetTransform( identityTransform );  
363.   
364.     if( argc > 5 )  
365.     {  
366.         writer->SetFileName( argv[5] );  
367.         writer->Update();  
368.     }  
369.   
370.     //配准后  
371.     resample->SetTransform( finalTransform );  
372.     if( argc > 6 )  
373.     {  
374.         writer->SetFileName( argv[6] );  
375.         writer->Update();  
376.     }  
377.   
378.     return EXIT_SUCCESS;  
379. }  

 

CMakelists.txt:

CMAKE_MINIMUM_REQUIRED(VERSION 2.4)

PROJECT(MultiResImageRegistration1)

FIND_PACKAGE(ITK)

IF(ITK_FOUND) 
    INCLUDE(${ITK_USE_FILE})

ELSE(ITK_FOUND) 
    MESSAGE(FATAL_ERROR 
            "ITK not found. Please set ITK_DIR.") 
ENDIF(ITK_FOUND)

ADD_EXECUTABLE(MultiResImageRegistration1 MultiResImageRegistration1.cxx )

TARGET_LINK_LIBRARIES(MultiResImageRegistration1 ITKIO ITKNumerics)

测试图像如下：固定图像，浮动图像

  

结果如下：配准后图像、配准前浮动图像与固定图像差异，配准后的浮动图像与固定图像的差异

  

  
程序中我们设置了最大层数为 3，配准过程处理了 3 层后成功并停止，最终参数为：

>MultiResImageRegistration1.exe fixed.png moving.png out.png 100 befor.png after.png

------------------------------------- 
MultiResolution Level : 0

0   -0.90675   [15.6488, 3.33401] 
。。。。 
。。。。 
73   -2.2638   [13.1505, 17.1896] 
74   -2.2659   [13.1597, 17.1936] 
------------------------------------- 
MultiResolution Level : 1

0   -2.18173   [10.1725, 14.5334] 
1   -1.29464   [13.4157, 16.0959] 
2   -1.91879   [11.6505, 18.8127] 
。。。。。 

102   -2.25089   [13.0004, 16.9996] 
103   -2.25085   [12.9996, 17.0003] 
-------------------------------------

------------------------------------- 
MultiResolution Level : 2

0   -2.80063   [13.9996, 17.0093] 
1   -1.45265   [13.0998, 16.992] 
2   -2.56817   [12.2098, 17.1264] 
3   -1.61491   [13.0139, 17.0287] 
4   -2.7529   [12.6293, 16.4095] 
5   -1.71975   [12.943, 16.9857] 
6   -2.66095   [13.5739, 17.1659] 
7   -1.80306   [13.0141, 16.978] 
8   -2.7611   [12.6566, 17.448] 
9   -1.84529   [12.9994, 17.0419] 
。。。。。

。。。。。 
109   -2.80063   [12.9998, 17.0004] 
110   -2.80062   [12.9998, 17.0003] 
111   -2.80062   [12.9998, 17.0003] 
Result = 
Translation X = 12.9998 
Translation Y = 17.0003 
Iterations    = 113 
Metric value  = -2.80062

可以看到，浮动图像平移了大约（13, 17)。