中直方图的绘制与保存_ITK_医学图像配准2_灰度直方图

/*绘制联合直方图，可以较精确的观察配准时内部的变化
图像的评价值。这部分描述了在优化器响应时间的参考区间内被连续的联合直方图监测的交 互信息的评价值的变化机制。



联合直方图的原理：
    1、设置了一个动态的调整范围函数：为后续图像的变换不断的调整
    2、创建联合直方图类：
        定义了很多东西
        联合直方图的写入
            针对迭代；针对名字；将其写入
    3、设置监控类：每次响应后，调用优化器的复写器去写联合直方图
    4、主函数：
        1、设置基本参数：参考图像、浮动图像、图像像素点
        2、基本操作：转换、优化、插值、配准并将其实例化链接
        3、对联合直方图进行参数设置
        4、迭代；观察
        5、输入图像，归一化，高斯滤波，平移优化
        6、输出初始直方图；迭代；输出每次迭代的参数
        7、获取配准后的联合直方图
        8、重采样滤波
        9、输出



*/


#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkNormalizeImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"

#include <iomanip>


#include "itkHistogramToEntropyImageFilter.h"
#include "itkMutualInformationHistogramImageToImageMetric.h"


#include "itkCommand.h"

#include <stdio.h>

template< class TInput >
class RescaleDynamicRangeFunctor//动态调整范围函数
{
public:
  typedef unsigned char OutputPixelType;
  RescaleDynamicRangeFunctor() {};
  ~RescaleDynamicRangeFunctor() {};
  inline OutputPixelType operator()( const TInput &A )
    {
    if( (A > 0.0) )
      {
      if( -(30.0 * std::log(A)) > 255 )
        {
        return static_cast<OutputPixelType>( 255 );
        }
      else
        {
        return itk::Math::Round<OutputPixelType>( -(30.0 * std::log(A)) );
        }
      }
    else
      {
      return static_cast<OutputPixelType>(255);
      }
    }
};






/*
   联合直方图的原理：
        忽略书上的灰度直方图概念，参考图像与浮动图像相同位置的点亮与不亮表示在直方图上

*/

//创建类HistogramWriter：作用是把联合直方图用ITK支持的格式保存下来 
//每一次优化器响应时，会对其进行调用。复写器会把联合直方图保存到xxx.mhd的文件中
//其中xxx是次数
namespace
{
class HistogramWriter//写联合直方图
{
public:
  typedef float InternalPixelType;
  itkStaticConstMacro( Dimension, unsigned int, 2);

  typedef itk::Image< InternalPixelType, Dimension > InternalImageType;

  typedef itk::MutualInformationHistogramImageToImageMetric<//互信息
                                        InternalImageType,
                                        InternalImageType >    MetricType;
  // Software Guide : EndCodeSnippet

  typedef MetricType::Pointer   MetricPointer;

  // Software Guide : BeginCodeSnippet
  typedef MetricType::HistogramType   HistogramType;

  typedef itk::HistogramToEntropyImageFilter< HistogramType, InternalImageType>
                                HistogramToEntropyImageFilterType;

  typedef HistogramToEntropyImageFilterType::Pointer
                                HistogramToImageFilterPointer;

  typedef HistogramToEntropyImageFilterType::OutputImageType OutputImageType;

  typedef itk::ImageFileWriter< OutputImageType > HistogramFileWriterType;
  typedef HistogramFileWriterType::Pointer        HistogramFileWriterPointer;
  // Software Guide : EndCodeSnippet

  typedef HistogramToEntropyImageFilterType::OutputPixelType OutputPixelType;
  
  




  
  
  //直方图的写入：
  HistogramWriter()://链接
    m_Metric(0)
    {
    this->m_Filter = HistogramToEntropyImageFilterType::New();
    this->m_HistogramFileWriter = HistogramFileWriterType::New();//输出
    this->m_HistogramFileWriter->SetInput( this->m_Filter->GetOutput() );//复写器的输入


    }

  ~HistogramWriter() { };

  void SetMetric( MetricPointer metric )
    {
    this->m_Metric = metric;
    }

  MetricPointer GetMetric() const
    {
    return this->m_Metric;
    }




  //写入联合直方图：针对迭代
  void WriteHistogramFile( unsigned int iterationNumber )
  {
    std::string outputFileBase = "JointHistogram";
    std::ostringstream outputFilename;
    outputFilename << outputFileBase
                   << "."
                   << std::setfill('0') << std::setw(3) << iterationNumber
                   << "."
                   << "mhd";
    m_HistogramFileWriter->SetFileName( outputFilename.str() );
    this->m_Filter->SetInput( this->GetMetric()->GetHistogram() );
    this->m_Filter->Modified();

    try
      {
      m_Filter->Update();
      }
    catch( itk::ExceptionObject & err )
      {
      std::cerr << "ERROR: ExceptionObject caught !" << std::endl;
      std::cerr << err << std::endl;
      }

    try
      {
      m_HistogramFileWriter->Update();
      }
    catch( itk::ExceptionObject & excp )
      {
      std::cerr << "Exception thrown " << excp << std::endl;
      }

    std::cout << "Joint Histogram file: ";
    std::cout << outputFilename.str() << " written" << std::endl;

    }




//写入联合联合直方图：针对文件名
  void WriteHistogramFile( std::string &outputFilename  )
    {
    // Software Guide : EndCodeSnippet


    // Software Guide : BeginCodeSnippet
    this->m_Filter->SetInput( this->GetMetric()->GetHistogram() );
    this->m_Filter->Modified();


    typedef itk::Image< unsigned char, Dimension > RescaledOutputImageType;

    typedef RescaleDynamicRangeFunctor<
                              OutputPixelType
                                      > RescaleDynamicRangeFunctorType;

    typedef itk::UnaryFunctorImageFilter<
                                OutputImageType,
                                RescaledOutputImageType,
                                RescaleDynamicRangeFunctorType
                                      > RescaleDynamicRangeFilterType;

    RescaleDynamicRangeFilterType::Pointer rescaler =
                                RescaleDynamicRangeFilterType::New();

    rescaler->SetInput( m_Filter->GetOutput() );

    typedef itk::ImageFileWriter< RescaledOutputImageType > RescaledWriterType;

    RescaledWriterType::Pointer rescaledWriter =
                                RescaledWriterType::New();

    rescaledWriter->SetInput( rescaler->GetOutput() );

    rescaledWriter->SetFileName( outputFilename );
      m_Filter->Update();
      rescaledWriter->Update();
    std::cout << "Joint Histogram file: " << outputFilename <<
        " written" << std::endl;
    }




private:
  MetricPointer                   m_Metric;
  HistogramToImageFilterPointer   m_Filter;
  HistogramFileWriterPointer      m_HistogramFileWriter;
  // Software Guide : EndCodeSnippet
  std::string   m_OutputFile;
};





// 监控类：每次响应后，调用优化器的复写器去写联合直方图
class CommandIterationUpdate : public itk::Command
{
public:
  typedef  CommandIterationUpdate   Self;
  typedef  itk::Command             Superclass;
  typedef  itk::SmartPointer<Self>  Pointer;
  itkSimpleNewMacro( Self );

protected:
  CommandIterationUpdate()
    {
    m_WriteHistogramsAfterEveryIteration = false;
    }

public:

  typedef   itk::RegularStepGradientDescentOptimizer     OptimizerType;
  typedef   const OptimizerType *                        OptimizerPointer;

  void Execute(itk::Object *caller, const itk::EventObject & event) ITK_OVERRIDE
    {
    Execute( (const itk::Object *)caller, event);
    }

  void Execute(const itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE
    {
    OptimizerPointer optimizer = static_cast< OptimizerPointer >( object );
    if( ! itk::IterationEvent().CheckEvent( &event ) || optimizer == ITK_NULLPTR )
      {
      return;
      }
    std::cout << optimizer->GetCurrentIteration() << "   ";
    std::cout << optimizer->GetValue() << "   ";
    std::cout << optimizer->GetCurrentPosition() << std::endl;
    //针对联合直方图的监控
    if( optimizer->GetCurrentIteration() == 0 )
      {

      m_JointHistogramWriter.WriteHistogramFile( m_InitialHistogramFile );//复写器
      
      }
    if( m_WriteHistogramsAfterEveryIteration )
      {
      m_JointHistogramWriter.WriteHistogramFile(
                                           optimizer->GetCurrentIteration() );//复写器
      }
    }

  void SetWriteHistogramsAfterEveryIteration( bool value )
    {
    m_WriteHistogramsAfterEveryIteration = value;
    }

  void SetInitialHistogramFile( const char * filename )
    {
    m_InitialHistogramFile = filename;
    }

  HistogramWriter m_JointHistogramWriter;

private:
  bool              m_WriteHistogramsAfterEveryIteration;//是否迭代
  std::string       m_InitialHistogramFile;

};

} 







int main(  )
{
  

  typedef  unsigned char  PixelType;

  const unsigned int Dimension = 2;

  typedef itk::Image< PixelType, Dimension >         FixedImageType;//参考图像
  typedef itk::Image< PixelType, Dimension >         MovingImageType;//浮动图像
  typedef   float                                    InternalPixelType;//图像种类
  typedef itk::Image< InternalPixelType, Dimension > InternalImageType;//像素种类




  //定义
  typedef itk::TranslationTransform< double, Dimension > TransformType;//转换
  typedef itk::RegularStepGradientDescentOptimizer       OptimizerType;//优化
  typedef itk::LinearInterpolateImageFunction<
                                    InternalImageType,
                                    double             > InterpolatorType;//插值
  typedef itk::ImageRegistrationMethod<
                                    InternalImageType,
                                    InternalImageType >  RegistrationType;//配准

  typedef itk::MutualInformationHistogramImageToImageMetric<
                                          InternalImageType,
                                          InternalImageType >    MetricType;//度量


  //实例化
  TransformType::Pointer      transform     = TransformType::New();//转换
  OptimizerType::Pointer      optimizer     = OptimizerType::New();//优化
  InterpolatorType::Pointer   interpolator  = InterpolatorType::New();//插值
  RegistrationType::Pointer   registration  = RegistrationType::New();//配准
  MetricType::Pointer         metric        = MetricType::New();//度量



  //安放
  registration->SetOptimizer(     optimizer     );
  registration->SetTransform(     transform     );
  registration->SetInterpolator(  interpolator  );




  //输出图像设置
  unsigned int numberOfHistogramBins = atoi( "256" );//联合直方图的个数；这个函数是用来将字符串转化为数字的
  MetricType::HistogramType::SizeType histogramSize;
  histogramSize.SetSize(2);
  histogramSize[0] = numberOfHistogramBins;
  histogramSize[1] = numberOfHistogramBins;
  metric->SetHistogramSize( histogramSize );//放入；渲染？



  //参数设置
  const unsigned int numberOfParameters = transform->GetNumberOfParameters();//参数
  typedef MetricType::ScalesType ScalesType;
  ScalesType scales( numberOfParameters );
  scales.Fill( 1.0 );
  metric->SetDerivativeStepLengthScales(scales);




  //迭代放入
  //观察者迭代
  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
 //度量
  observer->m_JointHistogramWriter.SetMetric( metric );
  registration->SetMetric( metric  );



  //输入图像
  typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
  typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;

  FixedImageReaderType::Pointer  fixedImageReader  = FixedImageReaderType::New();
  MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
  fixedImageReader->SetFileName("E:ProgramITKSourceInsightToolkit-4.13.2ExamplesDataBrainProtonDensitySliceBorder20.png");//输入图像文件
  movingImageReader->SetFileName("E:ProgramITKSourceInsightToolkit-4.13.2ExamplesDataBrainProtonDensitySliceShifted13x17y.png");

  //归一化：为了让亮度更加的标准
  typedef itk::NormalizeImageFilter<
                                FixedImageType,
                                InternalImageType
                                        > FixedNormalizeFilterType;

  typedef itk::NormalizeImageFilter<
                                MovingImageType,
                                InternalImageType
                                              > MovingNormalizeFilterType;

  FixedNormalizeFilterType::Pointer fixedNormalizer =
                                            FixedNormalizeFilterType::New();

  MovingNormalizeFilterType::Pointer movingNormalizer =
                                            MovingNormalizeFilterType::New();
  typedef itk::DiscreteGaussianImageFilter<
                                      InternalImageType,
                                      InternalImageType
                                                    > GaussianFilterType;


  //高斯滤波
  GaussianFilterType::Pointer fixedSmoother  = GaussianFilterType::New();
  GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();

  fixedSmoother->SetVariance( 2.0 );
  movingSmoother->SetVariance( 2.0 );
  fixedNormalizer->SetInput(  fixedImageReader->GetOutput() );
  movingNormalizer->SetInput( movingImageReader->GetOutput() );

  fixedSmoother->SetInput( fixedNormalizer->GetOutput() );
  movingSmoother->SetInput( movingNormalizer->GetOutput() );

  registration->SetFixedImage(    fixedSmoother->GetOutput()    );
  registration->SetMovingImage(   movingSmoother->GetOutput()   );
  fixedNormalizer->Update();


  registration->SetFixedImageRegion(
       fixedNormalizer->GetOutput()->GetBufferedRegion() );




  //获得每次平移的位置；优化参数
  typedef RegistrationType::ParametersType ParametersType;
  ParametersType initialParameters( transform->GetNumberOfParameters() );

  initialParameters[0] = 0.0;  // Initial offset in mm along X
  initialParameters[1] = 0.0;  // Initial offset in mm along Y

  registration->SetInitialTransformParameters( initialParameters );


  optimizer->SetMaximumStepLength( 4.00 );
  optimizer->SetMinimumStepLength( 0.01 );
  optimizer->SetRelaxationFactor(  0.90 );
  optimizer->SetNumberOfIterations( 10 );
  optimizer->MaximizeOn();
  optimizer->AddObserver( itk::IterationEvent(), observer );

  //初始化直方图文件
  observer->SetInitialHistogramFile( "E://2.png" );//一开始两个图像的直方图
  //迭代
  observer->SetWriteHistogramsAfterEveryIteration( true );
   



  //迭代步骤


    registration->Update();
    std::cout << "Optimizer stop condition: "
              << registration->GetOptimizer()->GetStopConditionDescription()
              << std::endl;



//获得每次迭代的参数然后输出
  ParametersType finalParameters = registration->GetLastTransformParameters();
  double TranslationAlongX = finalParameters[0];
  double TranslationAlongY = finalParameters[1];
  unsigned int numberOfIterations = optimizer->GetCurrentIteration();
  double bestValue = optimizer->GetValue();
  std::cout << "Result = " << std::endl;
  std::cout << " Translation X = " << TranslationAlongX  << std::endl;
  std::cout << " Translation Y = " << TranslationAlongY  << std::endl;
  std::cout << " Iterations    = " << numberOfIterations << std::endl;
  std::cout << " Metric value  = " << bestValue          << std::endl;




  //获取配准后的联合直方图
  std::string histogramAfter("E://3.png");
  observer->m_JointHistogramWriter.WriteHistogramFile( histogramAfter );//复写



  //重采样滤波器
  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() );

  FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();

  resample->SetSize(    fixedImage->GetLargestPossibleRegion().GetSize() );
  resample->SetOutputOrigin(  fixedImage->GetOrigin() );
  resample->SetOutputSpacing( fixedImage->GetSpacing() );
  resample->SetOutputDirection( fixedImage->GetDirection() );
  resample->SetDefaultPixelValue( 100 );





  //正常的输出
  typedef  unsigned char  OutputPixelType;

  typedef itk::Image< OutputPixelType, Dimension > OutputImageType;

  typedef itk::CastImageFilter<
                        FixedImageType,
                        OutputImageType > CastFilterType;

  typedef itk::ImageFileWriter< OutputImageType >  WriterType;
  WriterType::Pointer      writer =  WriterType::New();
  CastFilterType::Pointer  caster =  CastFilterType::New();
  writer->SetFileName( "E://4.png" );//输出
  caster->SetInput( resample->GetOutput() );
  writer->SetInput( caster->GetOutput()   );
  writer->Update();
  return EXIT_SUCCESS;

}

结果：

配准后的图像，原始直方图，迭代10次的直方图，迭代20次的直方图，迭代100次的直方图