itk::SymmetricForcesDemonsRegistrationFilter

对称Demons方法测试心得：

这里有个问题，二维的时候不会出错；

但是三维的时候会报错说两个图像集的phsical space 不一致！

这个错误在release下运行不了，但是在debug下出错的时候直接点击继续是不会有问题的，很奇怪，我把两个的图像集在输入之前原点手动设置为同一个，还是会报错！

Description: itk::ERROR: itk::ERROR: DiffeomorphicDemonsRegistrationFilter(0000019575A09930): Inputs do not occupy the same physical space!
InputImage Origin: [-1.6968750e+02, -1.5796870e+02, -3.3950000e+02], InputImageMovingImage Origin: [-1.6968750e+02, -1.5796870e+02, -3.2950000e+02]
        Tolerance: 5.0000000e-06
 

而且在实验itk::DiffeomorphicDemonsRegistrationFilter（微分同胚配准）的时候也遇到了这个问题。。。先记下来，后面找到原因再来更新吧

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkSymmetricForcesDemonsRegistrationFilter.h"
#include "itkHistogramMatchingImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkWarpImageFilter.h"
// Software Guide : EndCodeSnippet
//  The following section of code implements a Command observer
//  that will monitor the evolution of the registration process.
//
class CommandIterationUpdate : public itk::Command
{
public:
    typedef  CommandIterationUpdate                     Self;
    typedef  itk::Command                               Superclass;
    typedef  itk::SmartPointer<CommandIterationUpdate>  Pointer;
    itkNewMacro(CommandIterationUpdate);
protected:
    CommandIterationUpdate() {};
    typedef itk::Image< float, 3 >            InternalImageType;
    typedef itk::Vector< float, 3 >           VectorPixelType;
    typedef itk::Image<  VectorPixelType, 3 > DisplacementFieldType;
    typedef itk::SymmetricForcesDemonsRegistrationFilter<
        InternalImageType,
        InternalImageType,
        DisplacementFieldType>   RegistrationFilterType;
public:
    void Execute(itk::Object *caller, const itk::EventObject & event)
    {
        Execute((const itk::Object *)caller, event);
    }
    void Execute(const itk::Object * object, const itk::EventObject & event)
    {
        const RegistrationFilterType * filter =
            dynamic_cast<const RegistrationFilterType *>(object);
        if (!(itk::IterationEvent().CheckEvent(&event)))
        {
            return;
        }
        std::cout << filter->GetMetric() << std::endl;
    }
};
int main(int argc, char *argv[])
{

    // Software Guide : BeginLatex
    //
    // Second, we declare the types of the images.
    //
    // Software Guide : EndLatex
    // Software Guide : BeginCodeSnippet
    const unsigned int Dimension = 3;
        typedef unsigned short PixelType;
        typedef itk::Image< PixelType, Dimension > ImageType;
        typedef itk::ImageFileReader< ImageType > ImageReaderType;
    
        std::string pathfix = "D:/fix";
        std::string pathmov = "D:/mov";
        ImageType::Pointer fixedImage, movingImage;
        if (Dimension == 2) {
            ImageReaderType::Pointer fixedImageReader = ImageReaderType::New();
            ImageReaderType::Pointer movingImageReader = ImageReaderType::New();
    
            fixedImageReader->SetFileName(pathfix);
            movingImageReader->SetFileName(pathmov);
            fixedImage = fixedImageReader->GetOutput();
            movingImage = movingImageReader->GetOutput();
        }
        else
        {
            fixedImage = load_dicom(pathfix);
            movingImage = load_dicom(pathmov);
        }
    

    
    
    // Image file readers are set up in a similar fashion to previous examples.
    // To support the re-mapping of the moving image intensity, we declare an
    // internal image type with a floating point pixel type and cast the input
    // images to the internal image type.
    //
    // Software Guide : EndLatex
    // Software Guide : BeginCodeSnippet
    typedef float                                      InternalPixelType;
    typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
    typedef itk::CastImageFilter< ImageType,
        InternalImageType >  FixedImageCasterType;
    typedef itk::CastImageFilter< ImageType,
        InternalImageType >  MovingImageCasterType;
    FixedImageCasterType::Pointer fixedImageCaster = FixedImageCasterType::New();
    MovingImageCasterType::Pointer movingImageCaster
        = MovingImageCasterType::New();
    fixedImageCaster->SetInput(fixedImage);
    movingImageCaster->SetInput(movingImage);
   
    typedef itk::HistogramMatchingImageFilter<
        InternalImageType,
        InternalImageType >   MatchingFilterType;
    MatchingFilterType::Pointer matcher = MatchingFilterType::New();
   
    matcher->SetInput(movingImageCaster->GetOutput());
    matcher->SetReferenceImage(fixedImageCaster->GetOutput());
  
    matcher->SetNumberOfHistogramLevels(1024);
    matcher->SetNumberOfMatchPoints(7);
   
    matcher->ThresholdAtMeanIntensityOn();
   
    typedef itk::Vector< float, Dimension >                VectorPixelType;
    typedef itk::Image<  VectorPixelType, Dimension >      DisplacementFieldType;
    typedef itk::SymmetricForcesDemonsRegistrationFilter<
        InternalImageType,
        InternalImageType,
        DisplacementFieldType> RegistrationFilterType;
    RegistrationFilterType::Pointer filter = RegistrationFilterType::New();
   
    CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
    filter->AddObserver(itk::IterationEvent(), observer);
   
    filter->SetFixedImage(fixedImageCaster->GetOutput());
    filter->SetMovingImage(matcher->GetOutput());
   
    filter->SetNumberOfIterations(50);
    filter->SetStandardDeviations(1.0);
   
    filter->Update();
   
    typedef itk::WarpImageFilter<
        ImageType,
        ImageType,
        DisplacementFieldType  >     WarperType;
    typedef itk::LinearInterpolateImageFunction<
        ImageType,
        double          >  InterpolatorType;
    WarperType::Pointer warper = WarperType::New();
    InterpolatorType::Pointer interpolator = InterpolatorType::New();
  
    warper->SetInput(movingImage);
    warper->SetInterpolator(interpolator);
    warper->SetOutputSpacing(fixedImage->GetSpacing());
    warper->SetOutputOrigin(fixedImage->GetOrigin());
    warper->SetOutputDirection(fixedImage->GetDirection());
 
    warper->SetDisplacementField(filter->GetOutput());
    // Software Guide : EndCodeSnippet
    // Write warped image out to file
    typedef  unsigned char                           OutputPixelType;
    typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
    typedef itk::CastImageFilter<
        ImageType,
        OutputImageType >          CastFilterType;
    typedef itk::ImageFileWriter< OutputImageType >  WriterType;
    WriterType::Pointer      writer = WriterType::New();
    CastFilterType::Pointer  caster = CastFilterType::New();
    writer->SetFileName("0.mha");
    caster->SetInput(warper->GetOutput());
    writer->SetInput(caster->GetOutput());
    writer->Update();
   
    if (argc > 4) // if a fourth line argument has been provided...
    {
        // Software Guide : BeginCodeSnippet
        typedef itk::ImageFileWriter< DisplacementFieldType > FieldWriterType;
        FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
        fieldWriter->SetFileName(argv[4]);
        fieldWriter->SetInput(filter->GetOutput());
        fieldWriter->Update();
        // Software Guide : EndCodeSnippet
        // Software Guide : BeginLatex
        //
        // Note that the file format used for writing the deformation field must be
        // capable of representing multiple components per pixel. This is the case
        // for the MetaImage and VTK file formats for example.
        //
        // Software Guide : EndLatex
    }
    return EXIT_SUCCESS;
}