registration/clitkConvertBSplineDeformableTransformToVFGenericFilter.cxx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://www.centreleonberard.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ===========================================================================**/
  18 #ifndef clitkConvertBSplineDeformableTransformToVFGenericFilter_cxx
  19 #define clitkConvertBSplineDeformableTransformToVFGenericFilter_cxx
  20
  21 /* =================================================
  22  * @file   clitkConvertBSplineDeformableTransformToVFGenericFilter.cxx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "clitkConvertBSplineDeformableTransformToVFGenericFilter.h"
  31
  32
  33 namespace clitk
  34 {
  35
  36
  37   //-----------------------------------------------------------
  38   // Constructor
  39   //-----------------------------------------------------------
  40   ConvertBSplineDeformableTransformToVFGenericFilter::ConvertBSplineDeformableTransformToVFGenericFilter()
  41   {
  42     m_Verbose=false;
  43     m_InputFileName="";
  44   }
  45
  46   //-------------------------------------------------------------------
  47   // Update with the number of dimensions
  48   //-------------------------------------------------------------------
  49   template<>
  50   void
  51   ConvertBSplineDeformableTransformToVFGenericFilter::UpdateWithDim<3>(std::string PixelType, int Components)
  52   {
  53     // Components
  54     if (Components !=3)
  55       {
  56         std::cerr<<"Number of components is "<<Components<<"! Only 3 components is supported."<<std::endl;
  57         return;
  58       }
  59     if (PixelType != "double")
  60       {
  61         std::cerr<<"PixelType is  "<<PixelType<<"! Only double coefficient images are supported."<<std::endl;
  62         std::cerr<<"Reading image as double..."<<std::endl;
  63       }
  64
  65     // ImageTypes
  66     const unsigned int Dimension=3;
  67     typedef itk::Vector<double, Dimension> InputPixelType;
  68     typedef itk::Vector<float, Dimension> OutputPixelType;
  69     typedef itk::Image<InputPixelType, Dimension> InputImageType;
  70     typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
  71
  72     // Read the input
  73     typedef itk::ImageFileReader<InputImageType> InputReaderType;
  74     InputReaderType::Pointer reader = InputReaderType::New();
  75     reader->SetFileName( m_InputFileName);
  76     reader->Update();
  77     InputImageType::Pointer input= reader->GetOutput();
  78
  79
  80     // -----------------------------------------------
  81     // Filter
  82     // -----------------------------------------------
  83     typedef itk::TransformToDeformationFieldSource<OutputImageType, double> ConvertorType;
  84     ConvertorType::Pointer filter= ConvertorType::New();
  85
  86     //Output image info
  87     if (m_ArgsInfo.like_given)
  88       {
  89         typedef itk::ImageFileReader<OutputImageType> ReaderType;
  90         ReaderType::Pointer reader2=ReaderType::New();
  91         reader2->SetFileName(m_ArgsInfo.like_arg);
  92         reader2->Update();
  93
  94         OutputImageType::Pointer image=reader2->GetOutput();
  95         filter->SetOutputParametersFromImage(image);
  96       }
  97     else
  98       {
  99         unsigned int i=0;
 100         if(m_ArgsInfo.origin_given)
 101           {
 102             OutputImageType::PointType origin;
 103             for(i=0;i<Dimension;i++)
 104               origin[i]=m_ArgsInfo.origin_arg[i];
 105             filter->SetOutputOrigin(origin);
 106           }
 107         if (m_ArgsInfo.spacing_given)
 108           {
 109             OutputImageType::SpacingType spacing;
 110             for(i=0;i<Dimension;i++)
 111               spacing[i]=m_ArgsInfo.spacing_arg[i];
 112             filter->SetOutputSpacing(spacing);
 113           }
 114         if (m_ArgsInfo.spacing_given)
 115           {
 116             OutputImageType::SizeType size;
 117             for(i=0;i<Dimension;i++)
 118               size[i]=m_ArgsInfo.size_arg[i];
 119             filter->SetOutputSize(size);
 120           }
 121       }
 122
 123     if (m_Verbose)
 124       {
 125         std::cout<< "Setting output origin to "<<filter->GetOutputOrigin()<<"..."<<std::endl;
 126         std::cout<< "Setting output spacing to "<<filter->GetOutputSpacing()<<"..."<<std::endl;
 127         std::cout<< "Setting output size to "<<filter->GetOutputSize()<<"..."<<std::endl;
 128       }
 129
 130
 131     // -----------------------------------------------
 132     // Transform
 133     // -----------------------------------------------
 134     typedef clitk::BSplineDeformableTransform< double, Dimension, Dimension> TransformType;
 135     TransformType::Pointer transform=TransformType::New();
 136
 137     // Spline orders:  Default is cubic splines
 138     InputImageType::RegionType::SizeType splineOrders ;
 139     splineOrders.Fill(3);
 140     if (m_ArgsInfo.order_given)
 141       for(unsigned int i=0; i<Dimension;i++)
 142         splineOrders[i]=m_ArgsInfo.order_arg[i];
 143     if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 144
 145     // Mask
 146     typedef itk::ImageMaskSpatialObject<  Dimension >   MaskType;
 147     MaskType::Pointer  spatialObjectMask=NULL;
 148     if (m_ArgsInfo.mask_given)
 149       {
 150         typedef itk::Image< unsigned char, Dimension >   ImageMaskType;
 151         typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 152         MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 153         maskReader->SetFileName(m_ArgsInfo.mask_arg);
 154
 155         try
 156           {
 157             maskReader->Update();
 158           }
 159         catch( itk::ExceptionObject & err )
 160           {
 161             std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 162             std::cerr << err << std::endl;
 163             return;
 164           }
 165         if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 166
 167         // Set the image to the spatialObject
 168         spatialObjectMask = MaskType::New();
 169         spatialObjectMask->SetImage( maskReader->GetOutput() );
 170       }
 171
 172
 173     // Samplingfactors
 174     InputImageType::SizeType samplingFactors;
 175     for (unsigned int i=0; i< Dimension; i++)
 176       {
 177         samplingFactors[i]= (int) ( input->GetSpacing()[i]/ filter->GetOutputSpacing()[i]);
 178         if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 179       }
 180
 181
 182     // Set
 183     transform->SetSplineOrders(splineOrders);
 184     transform->SetMask(spatialObjectMask);
 185     transform->SetLUTSamplingFactors(samplingFactors);
 186     transform->SetCoefficientImage(input);
 187     filter->SetTransform(transform);
 188
 189
 190     // -----------------------------------------------
 191     // Update
 192     // -----------------------------------------------
 193     if (m_Verbose)std::cout<< "Converting the BSpline transform..."<<std::endl;
 194     try
 195       {
 196         filter->Update();
 197       }
 198     catch (itk::ExceptionObject)
 199       {
 200         std::cerr<<"Error: Exception thrown during execution convertion filter!"<<std::endl;
 201       }
 202
 203     OutputImageType::Pointer output=filter->GetOutput();
 204
 205
 206     // -----------------------------------------------
 207     // Output
 208     // -----------------------------------------------
 209     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 210     WriterType::Pointer writer = WriterType::New();
 211     writer->SetFileName(m_ArgsInfo.output_arg);
 212     writer->SetInput(output);
 213     writer->Update();
 214
 215   }
 216
 217
 218   //-------------------------------------------------------------------
 219   // Update with the number of dimensions
 220   //-------------------------------------------------------------------
 221   template<>
 222   void
 223   ConvertBSplineDeformableTransformToVFGenericFilter::UpdateWithDim<4>(std::string PixelType, int Components)
 224   {
 225     // Components
 226     if (Components !=3)
 227       {
 228         std::cerr<<"Number of components is "<<Components<<"! Only 3 components is supported."<<std::endl;
 229         return;
 230       }
 231     if (PixelType != "double")
 232       {
 233         std::cerr<<"PixelType is  "<<PixelType<<"! Only double coefficient images are supported."<<std::endl;
 234         std::cerr<<"Reading image as double..."<<std::endl;
 235       }
 236
 237     // ImageTypes
 238     const unsigned int Dimension=4;
 239     const unsigned int SpaceDimension=3;
 240     typedef itk::Vector<double, SpaceDimension> InputPixelType;
 241     typedef itk::Vector<float, SpaceDimension> OutputPixelType;
 242     typedef itk::Image<InputPixelType, Dimension> InputImageType;
 243     typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
 244
 245     // Read the input
 246     typedef itk::ImageFileReader<InputImageType> InputReaderType;
 247     InputReaderType::Pointer reader = InputReaderType::New();
 248     reader->SetFileName( m_InputFileName);
 249     reader->Update();
 250     InputImageType::Pointer input= reader->GetOutput();
 251
 252
 253     // -----------------------------------------------
 254     // Filter
 255     // -----------------------------------------------
 256     typedef clitk::TransformToDeformationFieldSource<OutputImageType, double> ConvertorType;
 257     ConvertorType::Pointer filter= ConvertorType::New();
 258
 259     //Output image info
 260     if (m_ArgsInfo.like_given)
 261       {
 262         typedef itk::ImageFileReader<OutputImageType> ReaderType;
 263         ReaderType::Pointer reader2=ReaderType::New();
 264         reader2->SetFileName(m_ArgsInfo.like_arg);
 265         reader2->Update();
 266
 267         OutputImageType::Pointer image=reader2->GetOutput();
 268         filter->SetOutputParametersFromImage(image);
 269       }
 270     else
 271       {
 272         unsigned int i=0;
 273         if(m_ArgsInfo.origin_given)
 274           {
 275             OutputImageType::PointType origin;
 276             for(i=0;i<Dimension;i++)
 277               origin[i]=m_ArgsInfo.origin_arg[i];
 278             filter->SetOutputOrigin(origin);
 279           }
 280         if (m_ArgsInfo.spacing_given)
 281           {
 282             OutputImageType::SpacingType spacing;
 283             for(i=0;i<Dimension;i++)
 284               spacing[i]=m_ArgsInfo.spacing_arg[i];
 285             filter->SetOutputSpacing(spacing);
 286           }
 287         if (m_ArgsInfo.spacing_given)
 288           {
 289             OutputImageType::SizeType size;
 290             for(i=0;i<Dimension;i++)
 291               size[i]=m_ArgsInfo.size_arg[i];
 292             filter->SetOutputSize(size);
 293           }
 294       }
 295     //Output image info
 296     if (m_Verbose)
 297       {
 298         std::cout<< "Setting output origin to "<<filter->GetOutputOrigin()<<"..."<<std::endl;
 299         std::cout<< "Setting output spacing to "<<filter->GetOutputSpacing()<<"..."<<std::endl;
 300         std::cout<< "Setting output size to "<<filter->GetOutputSize()<<"..."<<std::endl;
 301       }
 302
 303
 304     // -----------------------------------------------
 305     // Transform
 306     // -----------------------------------------------
 307     typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform< double, Dimension, Dimension > TransformType;
 308     TransformType::Pointer transform=TransformType::New();
 309     transform->SetTransformShape(m_ArgsInfo.shape_arg);
 310
 311     // Spline orders:  Default is cubic splines
 312     InputImageType::RegionType::SizeType splineOrders ;
 313     splineOrders.Fill(3);
 314     if (m_ArgsInfo.order_given)
 315       for(unsigned int i=0; i<Dimension;i++)
 316         splineOrders[i]=m_ArgsInfo.order_arg[i];
 317     if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 318
 319     // Mask
 320     typedef itk::ImageMaskSpatialObject<  Dimension >   MaskType;
 321     MaskType::Pointer  spatialObjectMask=NULL;
 322     if (m_ArgsInfo.mask_given)
 323       {
 324         typedef itk::Image< unsigned char, Dimension >   ImageMaskType;
 325         typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 326         MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 327         maskReader->SetFileName(m_ArgsInfo.mask_arg);
 328
 329         try
 330           {
 331             maskReader->Update();
 332           }
 333         catch( itk::ExceptionObject & err )
 334           {
 335             std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 336             std::cerr << err << std::endl;
 337             return;
 338           }
 339         if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 340
 341         // Set the image to the spatialObject
 342         spatialObjectMask = MaskType::New();
 343         spatialObjectMask->SetImage( maskReader->GetOutput() );
 344       }
 345
 346
 347     // Samplingfactors
 348     InputImageType::SizeType samplingFactors;
 349     for (unsigned int i=0; i< Dimension; i++)
 350       {
 351         samplingFactors[i]= (int) ( input->GetSpacing()[i]/ filter->GetOutputSpacing()[i]);
 352         if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 353       }
 354     if( !(m_ArgsInfo.shape_arg%2) )samplingFactors[Dimension-1]=5;
 355
 356     // Set
 357     transform->SetSplineOrders(splineOrders);
 358     transform->SetMask(spatialObjectMask);
 359     transform->SetLUTSamplingFactors(samplingFactors);
 360     transform->SetCoefficientImage(input);
 361     filter->SetTransform(transform);
 362
 363
 364     // -----------------------------------------------
 365     // Update
 366     // -----------------------------------------------
 367     if (m_Verbose)std::cout<< "Converting the BSpline transform..."<<std::endl;
 368     try
 369       {
 370         filter->Update();
 371       }
 372     catch (itk::ExceptionObject)
 373       {
 374         std::cerr<<"Error: Exception thrown during execution convertion filter!"<<std::endl;
 375       }
 376
 377     OutputImageType::Pointer output=filter->GetOutput();
 378
 379
 380     // -----------------------------------------------
 381     // Output
 382     // -----------------------------------------------
 383     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 384     WriterType::Pointer writer = WriterType::New();
 385     writer->SetFileName(m_ArgsInfo.output_arg);
 386     writer->SetInput(output);
 387     writer->Update();
 388
 389   }
 390
 391
 392
 393   //-----------------------------------------------------------
 394   // Update
 395   //-----------------------------------------------------------
 396   void ConvertBSplineDeformableTransformToVFGenericFilter::Update()
 397   {
 398
 399     // Read the Dimension and PixelType
 400     int Dimension, Components;
 401     std::string PixelType;
 402     ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components);
 403
 404     // Call UpdateWithDim
 405     //if(Dimension==2) UpdateWithDim<2>(PixelType, Components);
 406     if(Dimension==3) UpdateWithDim<3>(PixelType, Components);
 407     else if (Dimension==4) UpdateWithDim<4>(PixelType, Components);
 408     else
 409       {
 410         std::cout<<"Error, Only for 3 Dimensions!!!"<<std::endl ;
 411         return;
 412       }
 413   }
 414
 415 } //end clitk
 416
 417 #endif  //#define clitkConvertBSplineDeformableTransformToVFGenericFilter_cxx