void ITKTestClass::ImageRegistration1(std::string fixedImageFile, std::string movingImageFile, std::string outputImagefile, std::string differenceImageAfter, std::string differenceImageBefore, double initialStepLength) { 
    itk::PNGImageIOFactory::RegisterOneFactory(); itk::TIFFImageIOFactory::RegisterOneFactory(); itk::BMPImageIOFactory::RegisterOneFactory(); itk::GDCMImageIOFactory::RegisterOneFactory(); itk::JPEGImageIOFactory::RegisterOneFactory(); //1.定义待配准图像类型: 维数, 像素类型 const unsigned int Dimension = 2; typedef unsigned char PixelType; typedef itk::Image< PixelType, Dimension > FixedImageType; typedef itk::Image< PixelType, Dimension > MovingImageType; //2.定义配准框架的基本组件: 变换, 优化, 测度, 插值, 配准组件 //用于 2D 图像的一个刚体变换,该变换的唯一参数为: 空间坐标的类型 typedef itk::CenteredRigid2DTransform< double > TransformType; typedef itk::RegularStepGradientDescentOptimizer OptimizerType; typedef itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType > MetricType; typedef itk::LinearInterpolateImageFunction< MovingImageType, double > InterpolatorType; //配准组件: 该组件用于连接其它各组件 typedef itk::ImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType; MetricType::Pointer metric = MetricType::New(); OptimizerType::Pointer optimizer = OptimizerType::New(); InterpolatorType::Pointer interpolator = InterpolatorType::New(); RegistrationType::Pointer registration = RegistrationType::New(); //3.使用配准组件将：变换, 优化, 测度, 插值 四个基本组件连接至一起 registration->SetMetric(metric); registration->SetOptimizer(optimizer); registration->SetInterpolator(interpolator); TransformType::Pointer transform = TransformType::New(); registration->SetTransform(transform); //4.设置待配准图像及变换区域 //定义待配准两幅输入图像的 reader typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType; typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType; FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New(); MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New(); fixedImageReader->SetFileName(fixedImageFile); movingImageReader->SetFileName(movingImageFile); //本例, 固定图像与浮动图像都来自文件 //需要 itk::ImageRegistrationMethod 从 file readers 的输出获取输入 registration->SetFixedImage(fixedImageReader->GetOutput()); registration->SetMovingImage(movingImageReader->GetOutput()); //更新 readers,确保在初始化变换时图像参数(size,origin,spacing)有效. fixedImageReader->Update(); //可以将配准限制在固定图像的一特定区域,通过将其作为测度(metric)计算的输入. //该区域通过方法SetFixedImageRegion() 定义. 通过使用这个特征, 可以避免图像 //中某些不需要的的对象影响配准输出,或者减少计算时间,本例中使用整幅图像.  //此区域通过 BufferedRegion 标识. registration->SetFixedImageRegion( fixedImageReader->GetOutput()->GetBufferedRegion()); //更新固定,浮动图像 reader, 使其 origin, size, spacing 有效, fixedImageReader->Update(); movingImageReader->Update(); //5.设置各组件的参数，主要是变换、优化组件的参数 //本例使用的为变换为: 先进行一个旋转变换,再进行一平移变换. typedef FixedImageType::SpacingType SpacingType; typedef FixedImageType::PointType OriginType; typedef FixedImageType::RegionType RegionType; typedef FixedImageType::SizeType SizeType; FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput(); const SpacingType fixedSpacing = fixedImage->GetSpacing(); const OriginType fixedOrigin = fixedImage->GetOrigin(); const RegionType fixedRegion = fixedImage->GetLargestPossibleRegion(); const SizeType fixedSize = fixedRegion.GetSize(); //本例使用固定图像的中心点做为旋转中心 //使用固定图与浮动图中心之间距离向量作为平移量 //首先计算固定图像中心点 TransformType::InputPointType centerFixed; centerFixed[0] = fixedOrigin[0] + fixedSpacing[0] * fixedSize[0] / 2.0; centerFixed[1] = fixedOrigin[1] + fixedSpacing[1] * fixedSize[1] / 2.0; MovingImageType::Pointer movingImage = movingImageReader->GetOutput(); const SpacingType movingSpacing = movingImage->GetSpacing(); const OriginType movingOrigin = movingImage->GetOrigin(); const RegionType movingRegion = movingImage->GetLargestPossibleRegion(); const SizeType movingSize = movingRegion.GetSize(); //然后计算浮动图像中心点 TransformType::InputPointType centerMoving; centerMoving[0] = movingOrigin[0] + movingSpacing[0] * movingSize[0] / 2.0; centerMoving[1] = movingOrigin[1] + movingSpacing[1] * movingSize[1] / 2.0; //设置旋转中心 transform->SetCenter(centerFixed); //设置平移量 transform->SetTranslation(centerMoving - centerFixed); //设置初始旋转角度 transform->SetAngle(0.0); //然后将当前的变换参数作为初始参数, 传递给配准过程 registration->SetInitialTransformParameters(transform->GetParameters()); //设置优化组件参数. //注意: 旋转和变换的"单位刻度" 是不同的,旋转用弧度度量, 平移以毫米度量 typedef OptimizerType::ScalesType OptimizerScalesType; OptimizerScalesType optimizerScales(transform->GetNumberOfParameters()); const double translationScale = 1.0 / 1000.0; optimizerScales[0] = 1.0; optimizerScales[1] = translationScale; optimizerScales[2] = translationScale; optimizerScales[3] = translationScale; optimizerScales[4] = translationScale; optimizer->SetScales(optimizerScales); //本例使用的优化器是 itk::RegularStepGradientDescentOptimizer,  //梯度下降法的一种,下面定义优化参数:  //relaxation fator, initial step length, minimal step length, number of iterations optimizer->SetRelaxationFactor(0.6); //松驰因子 optimizer->SetMaximumStepLength(initialStepLength); //最大步长 //最小步长和迭代最大次数, 用来限制优化过程迭代的执行 optimizer->SetMinimumStepLength(0.001); //最小步长 optimizer->SetNumberOfIterations(200); //最大迭代次数,以防止无限次迭代 /* *应该注意: 配准过程是一个优化的过程, 优化过程每迭代一次, *就与相似性测度进行一次比较 *当相似性测度值达到我们设定的预期值,或者达到设定的迭代上限时 *就停止迭代,得到最终结果. */ //6.实例化一 Command/Observer 对象, 监视配准过程的执行, 并触发配准过程的执行. //实例化一个 Command/Observer 对象, 将其注册为 optimizer 的观察者 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New(); optimizer->AddObserver(itk::IterationEvent(), observer); //取得配准执行完毕的最终参数 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters(); const double finalAngle = finalParameters[0]; const double finalRotationCenterX = finalParameters[1]; const double finalRotationCenterY = finalParameters[2]; const double finalTranslationX = finalParameters[3]; const double finalTranslationY = finalParameters[4]; const unsigned int numberOfIterations = optimizer->GetCurrentIteration(); const double bestValue = optimizer->GetValue(); const double finalAngleInDegrees = finalAngle * 180.0 / vnl_math::pi; std::cout << "Result = " << std::endl; std::cout << " Angle (radians) = " << finalAngle << std::endl; std::cout << " Angle (degrees) = " << finalAngleInDegrees << std::endl; std::cout << " Center X = " << finalRotationCenterX << std::endl; std::cout << " Center Y = " << finalRotationCenterY << std::endl; std::cout << " Translation X = " << finalTranslationX << std::endl; std::cout << " Translation Y = " << finalTranslationY << std::endl; std::cout << " Iterations = " << numberOfIterations << std::endl; std::cout << " Metric value = " << bestValue << std::endl; //7.重采样, 将变换后的浮动图像映射到固定图像空间中,保存配准结果 //一般情况,最后一步为: 利用最终的变换结果将浮动图像映射到固定图像空间 //可通过 itk::ResampleImageFilter 完成 typedef itk::ResampleImageFilter< MovingImageType, FixedImageType > ResampleFilterType; TransformType::Pointer finalTransform = TransformType::New(); finalTransform->SetParameters(finalParameters); finalTransform->SetFixedParameters(transform->GetFixedParameters()); //设置重采样过滤器的相应参数 ResampleFilterType::Pointer resample = ResampleFilterType::New(); resample->SetTransform(finalTransform); //将浮动图连接至重采样过滤器的输入端 resample->SetInput(movingImageReader->GetOutput()); resample->SetSize(fixedImage->GetLargestPossibleRegion().GetSize()); resample->SetOutputOrigin(fixedImage->GetOrigin()); resample->SetOutputSpacing(fixedImage->GetSpacing()); resample->SetOutputDirection(fixedImage->GetDirection()); //设置默认灰度值,用来 "突出" 显示映射后在浮动图像之外的区域 resample->SetDefaultPixelValue(100); //保存配准后图像 typedef itk::ImageFileWriter< FixedImageType > WriterFixedType; WriterFixedType::Pointer writer = WriterFixedType::New(); writer->SetFileName(outputImagefile); writer->SetInput(resample->GetOutput()); //重采样过滤器的输出 try { 
    writer->Update(); } catch (itk::ExceptionObject & excp) { 
    std::cerr << "ExceptionObject while writing the resampled image !" << std::endl; std::cerr << excp << std::endl; } //通过 itk::SubtractImageFilter 比较"经过变换的浮动图像"与"固定图像"之间的差异 typedef itk::Image< float, Dimension > DifferenceImageType; typedef itk::SubtractImageFilter< FixedImageType, FixedImageType, DifferenceImageType > DifferenceFilterType; DifferenceFilterType::Pointer difference = DifferenceFilterType::New(); typedef unsigned char OutputPixelType; typedef itk::Image< OutputPixelType, Dimension > OutputImageType; //由于两幅图像之间的差异可能对应非常小的亮度值 //使用 itk::RescaleIntensityImageFilter 重新调节亮度值,使其更加明显 typedef itk::RescaleIntensityImageFilter< DifferenceImageType, OutputImageType > RescalerType; RescalerType::Pointer intensityRescaler = RescalerType::New(); intensityRescaler->SetOutputMinimum(0); intensityRescaler->SetOutputMaximum(255); difference->SetInput1(fixedImageReader->GetOutput()); difference->SetInput2(resample->GetOutput()); resample->SetDefaultPixelValue(1); intensityRescaler->SetInput(difference->GetOutput()); typedef itk::ImageFileWriter< OutputImageType > WriterType; WriterType::Pointer writer2 = WriterType::New(); writer2->SetInput(intensityRescaler->GetOutput()); try { 
    //保存配准后, 浮动图与固定图之间差异的对比 writer2->SetFileName(differenceImageAfter); writer2->Update(); //保存配准前, 浮动图与固定图之间差异的对比, 这里仍然需要重采样 //因为浮动图与固定图不一定有相同的 origin, size, spacing.  //可以过一个恒等变换(Identity Transform)完成 TransformType::Pointer identityTransform = TransformType::New(); identityTransform->SetIdentity(); resample->SetTransform(identityTransform); writer2->SetFileName(differenceImageBefore); writer2->Update(); } catch (itk::ExceptionObject & excp) { 
    std::cerr << "Error while writing difference images" << std::endl; std::cerr << excp << std::endl; } }