[Insight-users] Ups Sorry the code

javier silva bravo javier_silva_bravo at hotmail.com
Fri May 20 11:39:43 EDT 2005


An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20050520/adee4e6c/attachment.htm
-------------- next part --------------
#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
//#include "itkLinearInterpolateImageFunction.h"
#include "itkBSplineInterpolateImageFunction.h"
#include "itkImage.h"
#include "itkTimeProbesCollectorBase.h"
#include "itkBSplineDeformableTransform.h"
//#include <itkElasticBodySplineKernelTransform.h>
#include "itkLBFGSBOptimizer.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
//#include "itkSquaredDifferenceImageFilter.h"
#include "itkCommand.h"
#include <fstream>

class CommandIterationUpdate : public itk::Command
{
public:
  typedef  CommandIterationUpdate   Self;
  typedef  itk::Command             Superclass;
  typedef itk::SmartPointer<Self>  Pointer;
  itkNewMacro( Self );
protected:
  CommandIterationUpdate() {};
  std::ofstream m_file;

public:
  typedef itk::LBFGSBOptimizer     OptimizerType;
  typedef   const OptimizerType   *    OptimizerPointer;

  void SetFileName(const char* str)
    {
	  try{
    m_file.open(str, std::ios::trunc);
	  }
	  catch(itk::ExceptionObject & err)
	  {
    std::cerr << "ExceptionObject caught !" << std::endl;
	std::cerr << "Error al tratar de leer el archivo de los valores" << 
std::endl;
	std::cerr << err << std::endl;
     };
    }


  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)
    {
      OptimizerPointer optimizer =
        dynamic_cast< OptimizerPointer >( object );
      if( typeid( event ) != typeid( itk::IterationEvent ) )
        {
        return;
        }
      std::cout << optimizer->GetCurrentIteration() << "          ";
      std::cout << optimizer->GetValue() << "    ";
      std::cout << optimizer->GetInfinityNormOfProjectedGradient() << 
std::endl;

	  m_file << optimizer->GetCurrentIteration() << "   ";
      m_file << optimizer->GetValue() << "    ";
      m_file << optimizer->GetInfinityNormOfProjectedGradient() <<std::endl;

    }
};


int main( int argc, char *argv[] )
{
  /*if( argc < 4 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile  movingImageFile outputImagefile  ";
    std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";
    std::cerr << " [deformationField] ";
    return 1;
    }*/
	if( argc < 3 )
    {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0];
    std::cerr << " fixedImageFile  movingImageFile outputImagefile  ";
    std::cerr << " [ValuesField] ";
    return 1;
    }


  std::cout << "Inicio del programa " << std::endl;
  const    unsigned int    ImageDimension = 2;
  typedef  float           PixelType;

  typedef itk::Image< PixelType, ImageDimension >  FixedImageType;
  typedef itk::Image< PixelType, ImageDimension >  MovingImageType;


  const unsigned int SpaceDimension = ImageDimension;
  const unsigned int SplineOrder = 3;
  typedef double CoordinateRepType;

  typedef itk::BSplineDeformableTransform<
                            CoordinateRepType,
                            SpaceDimension,
                            SplineOrder >     TransformType;


/* typedef itk::ElasticBodySplineKernelTransform<
                            CoordinateRepType,
                            SpaceDimension>     TransformType;*/

  typedef itk::LBFGSBOptimizer       OptimizerType;


  typedef itk::MeanSquaresImageToImageMetric<
                                    FixedImageType,
                                    MovingImageType >    MetricType;

/* typedef itk:: LinearInterpolateImageFunction<
                                    MovingImageType,
                                    double          >    InterpolatorType;*/
  typedef itk::BSplineInterpolateImageFunction<
	                                 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();

  interpolator->SetSplineOrder(SplineOrder);

  metric->SetInterpolator(  interpolator  );

  registration->SetMetric(        metric        );
  registration->SetOptimizer(     optimizer     );
  registration->SetInterpolator(  interpolator  );



  TransformType::Pointer  transform = TransformType::New();
  registration->SetTransform( transform );
  metric->SetTransform( transform );

  typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
  typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;

  FixedImageReaderType::Pointer  fixedImageReader  = 
FixedImageReaderType::New();
  MovingImageReaderType::Pointer movingImageReader = 
MovingImageReaderType::New();

  fixedImageReader->SetFileName(  argv[1] );
  movingImageReader->SetFileName( argv[2] );

  FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
  MovingImageType::ConstPointer movingImage = 
movingImageReader->GetOutput();

  registration->SetFixedImage(  fixedImage   );
  registration->SetMovingImage(   movingImage   );

  metric->SetFixedImage(  fixedImage   );
  metric->SetMovingImage(   movingImage   );
  interpolator->SetInputImage(movingImage);
  try{
  fixedImageReader->Update();
  movingImageReader->Update();
  }
  catch(itk::ExceptionObject & err)
  {
	  std::cerr<<"Error al leer Archivo"<<std::endl;
	  std::cerr<<err<<std::endl;
	  return -1;
  }
  FixedImageType::RegionType fixedRegion = 
fixedImage->GetLargestPossibleRegion();//fixedImage->GetLargestPossibleRegion();->GetBufferedRegion();

  registration->SetFixedImageRegion( fixedRegion );
  metric->SetFixedImageRegion( fixedRegion );

  typedef TransformType::RegionType RegionType;
  RegionType bsplineRegion;
  RegionType::SizeType   gridSizeOnImage;
  RegionType::SizeType   gridBorderSize;
  RegionType::SizeType   totalGridSize;
  /*typedef TransformType::IndexType IndexType;
  IndexType index;*/
  gridSizeOnImage.Fill( 5 );
  gridBorderSize.Fill( 3 );    // Border for spline order = 3 ( 1 lower, 2 
upper )
  totalGridSize = gridSizeOnImage + gridBorderSize;

  /*index[0]=0;
  index[1]=3;
  bsplineRegion.SetIndex(index);*/
  bsplineRegion.SetSize( totalGridSize );

  typedef TransformType::SpacingType SpacingType;
  SpacingType spacing = fixedImage->GetSpacing();

  typedef TransformType::OriginType OriginType;
  OriginType origin = fixedImage->GetOrigin();

  FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
  std::cout<<"Tamaño de la region fija"<<"   ";
  std::cout<<fixedRegion.GetSize()<<std::endl;
  for(unsigned int r=0; r<ImageDimension; r++)
    {
    spacing[r] *= floor( static_cast<double>(fixedImageSize[r] - 1)  /
                  static_cast<double>(gridSizeOnImage[r] - 1) );
	std::cout<<"spacing[r]"<<"  ";
	std::cout<<spacing[r]<<std::endl;
	origin[r]  -=  spacing[r];
    }

  transform->SetGridSpacing( spacing );
  transform->SetGridOrigin( origin );
  transform->SetGridRegion( bsplineRegion );
  std::cout<<"spacing"<<"  ";
  std::cout<<spacing<<std::endl;
  std::cout<<"origin"<<"   ";
  std::cout<<origin<<std::endl;
  std::cout<<"BSplineRegion"<<"  ";
  std::cout<<bsplineRegion<<std::endl;

  typedef TransformType::ParametersType     ParametersType;

  const unsigned int numberOfParameters =
               transform->GetNumberOfParameters();

  ParametersType parameters( numberOfParameters );

  parameters.Fill( 0.0 );

  transform->SetParameters( parameters );

  registration->SetInitialTransformParameters( transform->GetParameters() );
  metric->SetTransformParameters( transform->GetParameters() );
  metric->SetComputeGradient(metric->GetComputeGradient());

  std::cout << metric->GetComputeGradient() << std::endl;

  OptimizerType::BoundSelectionType boundSelect( 
transform->GetNumberOfParameters() );
  OptimizerType::BoundValueType upperBound( 
transform->GetNumberOfParameters() );
  OptimizerType::BoundValueType lowerBound( 
transform->GetNumberOfParameters() );

  boundSelect.Fill( 0 );
  upperBound.Fill( 0.0 );
  lowerBound.Fill( 0.0 );

  optimizer->SetBoundSelection( boundSelect );
  optimizer->SetUpperBound( upperBound );
  optimizer->SetLowerBound( lowerBound );
  optimizer->SetCostFunction(metric);
  optimizer->SetCostFunctionConvergenceFactor( 1e+7 );
  optimizer->SetProjectedGradientTolerance( 0.02 );//1.0 0.155555
  optimizer->SetMaximumNumberOfIterations( 50 );
  optimizer->SetMaximumNumberOfEvaluations( 500 );
  optimizer->SetMaximumNumberOfCorrections( 5 );
  std::cout << " SetProjectedGradientTolerance 0.02  " << std::endl;

  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
  try{
  observer->SetFileName(argv[4]);
  }
  catch(itk::ExceptionObject & err)
  {
	  std::cerr <<" Error al leer Archivo "<<std::endl;
	  std::cerr <<err<<std::endl;
  }
  optimizer->AddObserver( itk::IterationEvent(), observer );



  itk::TimeProbesCollectorBase collector;

  std::cout << std::endl << "Starting Registration" << std::endl;
  std::cout << "Iteracion  Valor   InfinityNormOfProjectedGradient  " << 
std::endl;

  try
    {
    collector.Start( "Registration" );
    registration->StartRegistration();
    collector.Stop( "Registration" );
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return -1;
    }
      collector.Report();

  OptimizerType::ParametersType finalParameters =
                    registration->GetLastTransformParameters();

  transform->SetParameters( finalParameters );

  typedef itk::ResampleImageFilter<
                            MovingImageType,
                            FixedImageType>    ResampleFilterType;

  ResampleFilterType::Pointer resample = ResampleFilterType::New();

  resample->SetTransform( transform );
  resample->SetInput( movingImageReader->GetOutput() );
  resample->SetInterpolator(interpolator);
  resample->SetSize(    fixedImage->GetLargestPossibleRegion().GetSize() );
  resample->SetOutputOrigin(  fixedImage->GetOrigin() );
  resample->SetOutputSpacing( fixedImage->GetSpacing() );
  resample->SetDefaultPixelValue( 100 );

  typedef  unsigned char  OutputPixelType;

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

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

  typedef itk::ImageFileWriter< OutputImageType >  WriterType;


  WriterType::Pointer      writer =  WriterType::New();
  CastFilterType::Pointer  caster =  CastFilterType::New();


  writer->SetFileName( argv[3] );


  caster->SetInput( resample->GetOutput() );
  writer->SetInput( caster->GetOutput()   );

  std::cout << "Creando Imagen De Salida ..." << std::endl;
  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & err )
    {
    std::cerr << "ExceptionObject caught !" << std::endl;
    std::cerr << err << std::endl;
    return -1;
    }

/*

  typedef itk::SquaredDifferenceImageFilter<
                                  FixedImageType,
                                  FixedImageType,
                                  OutputImageType > DifferenceFilterType;

  DifferenceFilterType::Pointer difference = DifferenceFilterType::New();

  WriterType::Pointer writer2 = WriterType::New();
  writer2->SetInput( difference->GetOutput() );


  // Compute the difference image between the
  // fixed and resampled moving image.
  if( argc >= 5 )
    {
    difference->SetInput1( fixedImageReader->GetOutput() );
    difference->SetInput2( resample->GetOutput() );
    writer2->SetFileName( argv[4] );
    try
      {
      writer2->Update();
      }
    catch( itk::ExceptionObject & err )
      {
      std::cerr << "ExceptionObject caught !" << std::endl;
      std::cerr << err << std::endl;
      return -1;
      }
    }


  // Compute the difference image between the
  // fixed and moving image before registration.
  if( argc >= 6 )
    {
    writer2->SetFileName( argv[5] );
    difference->SetInput1( fixedImageReader->GetOutput() );
    difference->SetInput2( movingImageReader->GetOutput() );
    try
      {
      writer2->Update();
      }
    catch( itk::ExceptionObject & err )
      {
      std::cerr << "ExceptionObject caught !" << std::endl;
      std::cerr << err << std::endl;
      return -1;
      }
    }



  // Generate the explicit deformation field resulting from
  // the registration.
  if( argc >= 7 )
    {

    typedef itk::Vector< float, ImageDimension >  VectorType;
    typedef itk::Image< VectorType, ImageDimension >  DeformationFieldType;

    DeformationFieldType::Pointer field = DeformationFieldType::New();
    field->SetRegions( fixedRegion );
    field->SetOrigin( fixedImage->GetOrigin() );
    field->SetSpacing( fixedImage->GetSpacing() );
    field->Allocate();

    typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;
    FieldIterator fi( field, fixedRegion );

    fi.GoToBegin();

    TransformType::InputPointType  fixedPoint;
    TransformType::OutputPointType movingPoint;
    DeformationFieldType::IndexType index;

    VectorType displacement;

    while( ! fi.IsAtEnd() )
      {
      index = fi.GetIndex();
      field->TransformIndexToPhysicalPoint( index, fixedPoint );
      movingPoint = transform->TransformPoint( fixedPoint );
      displacement[0] = movingPoint[0] - fixedPoint[0];
      displacement[1] = movingPoint[1] - fixedPoint[1];
      fi.Set( displacement );
      ++fi;
      }



    typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
    FieldWriterType::Pointer fieldWriter = FieldWriterType::New();

    fieldWriter->SetInput( field );

    fieldWriter->SetFileName( argv[6] );
    try
      {
      fieldWriter->Update();
      }
    catch( itk::ExceptionObject & excp )
      {
      std::cerr << "Exception thrown " << std::endl;
      std::cerr << excp << std::endl;
      return EXIT_FAILURE;
      }
    }*/
  std::cout << "El programa termino satisfactoriamente" << std::endl;
  return 0;
}






More information about the Insight-users mailing list