tools/clitkVFMerge.cxx

   1 #ifndef CLITKVFMERGE_CXX
   2 #define CLITKVFMERGE_CXX
   3
   4 /**
   5  * @file   clitkVFMerge.cxx
   6  * @author Jef Vandemeulebroucke <jefvdmb@gmail.com>
   7  * @date   June 15  10:14:53 2007
   8  *
   9  * @brief  Read in one VF (ex mhd, vf)  and write to another.  Transforming from mm to voxels needed for the vf format is implemented in clitkVfImageIO.cxx .
  10  *
  11  *
  12  */
  13
  14 // clitk include
  15 #include "clitkVFMerge_ggo.h"
  16 #include "clitkIO.h"
  17 #include "clitkIOCommon.h"
  18
  19 // itk include
  20 //#include "itkReadTransform.h"
  21 #include "itkImageFileWriter.h"
  22 #include <iostream>
  23 #include "itkImageFileReader.h"
  24 //#include "itkRawImageIO.h"
  25 //#include "macro.h"
  26
  27 int main( int argc, char *argv[] )
  28 {
  29
  30   // Init command line
  31   GGO(clitkVFMerge, args_info);
  32   CLITK_INIT;
  33
  34 const unsigned int SpaceDimension = 3;
  35 const unsigned int ModelDimension = 4;
  36 typedef itk::Vector< float, SpaceDimension > Displacement;
  37 typedef itk::Image< Displacement, SpaceDimension > DeformationFieldType;
  38 typedef itk::Image< Displacement, ModelDimension > DynamicDeformationFieldType;
  39 typedef itk::ImageFileReader< DeformationFieldType > DeformationFieldReaderType;
  40 typedef itk::ImageFileWriter< DynamicDeformationFieldType > DynamicDeformationFieldWriterType;
  41
  42
  43 //declare the dynamic
  44  DynamicDeformationFieldType::Pointer  dynamicDeformationField=DynamicDeformationFieldType::New();
  45
  46
  47  //declare their iterators
  48  typedef itk::ImageRegionIterator< DynamicDeformationFieldType> DynamicDeformationFieldIteratorType;
  49  DynamicDeformationFieldIteratorType *dynamicIteratorPointer= new DynamicDeformationFieldIteratorType;
  50
  51  for (unsigned int i=0 ; i< args_info.inputs_num ; i ++)
  52    {
  53      //read in the deformation field i
  54      DeformationFieldReaderType::Pointer deformationFieldReader = DeformationFieldReaderType::New();
  55      deformationFieldReader->SetFileName( args_info.inputs[i]);
  56      if (args_info.verbose_flag) std::cout<<"Reading VF number "<< i+1 << std::endl;
  57      deformationFieldReader->Update();
  58      DeformationFieldType::Pointer currentDeformationField = deformationFieldReader->GetOutput();
  59
  60      //create an iterator for the current deformation field
  61      typedef itk::ImageRegionIterator<DeformationFieldType> FieldIteratorType;
  62      FieldIteratorType fieldIterator(currentDeformationField, currentDeformationField->GetLargestPossibleRegion());
  63
  64      //Allocate memory for the dynamic components
  65      if (i==0)
  66        {
  67          DynamicDeformationFieldType::RegionType dynamicDeformationFieldRegion;
  68          DynamicDeformationFieldType::RegionType::SizeType dynamicDeformationFieldSize;
  69          DeformationFieldType::RegionType::SizeType deformationFieldSize;
  70          deformationFieldSize= currentDeformationField->GetLargestPossibleRegion().GetSize();
  71          dynamicDeformationFieldSize[0]=deformationFieldSize[0];
  72          dynamicDeformationFieldSize[1]=deformationFieldSize[1];
  73          dynamicDeformationFieldSize[2]=deformationFieldSize[2];
  74          dynamicDeformationFieldSize[3]=args_info.inputs_num;
  75          dynamicDeformationFieldRegion.SetSize(dynamicDeformationFieldSize);
  76          DynamicDeformationFieldType::IndexType start;
  77          start.Fill(0);
  78          dynamicDeformationFieldRegion.SetIndex(start);
  79          dynamicDeformationField->SetRegions(dynamicDeformationFieldRegion);
  80          dynamicDeformationField->Allocate();
  81
  82
  83             //Set the spacing
  84              DeformationFieldType::SpacingType spacing= currentDeformationField->GetSpacing();
  85              DynamicDeformationFieldType::SpacingType dynamicSpacing;
  86             dynamicSpacing[0]=spacing[0];
  87             dynamicSpacing[1]=spacing[1];
  88             dynamicSpacing[2]=spacing[2];
  89             dynamicSpacing[3]=args_info.spacing_arg; //JV par exemple
  90             dynamicDeformationField->SetSpacing(dynamicSpacing);
  91         DynamicDeformationFieldType::PointType origin;
  92         origin[0]=args_info.xorigin_arg;
  93         origin[1]=args_info.yorigin_arg;
  94         origin[2]=args_info.zorigin_arg;
  95         origin[3]=0; //temporal origin is always 0
  96         dynamicDeformationField->SetOrigin(origin);
  97
  98             // Creating iterators for the currentDeformationField and the DynamicDeformationField
  99             DynamicDeformationFieldIteratorType *dynamicIterator= new DynamicDeformationFieldIteratorType(dynamicDeformationField, dynamicDeformationField->GetLargestPossibleRegion());
 100             dynamicIteratorPointer=dynamicIterator;
 101             dynamicIteratorPointer->GoToBegin();
 102        }
 103      if (args_info.verbose_flag) std::cout<<"Merging VF number "<< i+1 << std::endl;
 104      //Copy the current component of the input into dynamicDeformationFieldComponent
 105      fieldIterator.GoToBegin();
 106      while(!fieldIterator.IsAtEnd())
 107        {
 108          dynamicIteratorPointer->Set(fieldIterator.Get());
 109          ++fieldIterator;
 110          ++(*dynamicIteratorPointer);
 111        }
 112    }
 113
 114
 115 //Write the vector field
 116  DynamicDeformationFieldWriterType::Pointer writer = DynamicDeformationFieldWriterType::New();
 117  writer->SetInput( dynamicDeformationField );
 118  writer->SetFileName( args_info.output_arg );
 119  if (args_info.verbose_flag) std::cout<<"Writing the dynamic VF"<< std::endl;
 120
 121
 122  try {
 123
 124    writer->Update( );
 125  }
 126  catch( itk::ExceptionObject & excp ) {
 127    std::cerr << "Problem writing the output file" << std::endl;
 128    std::cerr << excp << std::endl;
 129    return EXIT_FAILURE;
 130  }
 131  return EXIT_SUCCESS;
 132 }
 133 #endif
 134
 135