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 #if ITK_VERSION_MAJOR >= 4
  84     typedef itk::TransformToDisplacementFieldSource<OutputImageType, double> ConvertorType;
  85 #else
  86     typedef itk::TransformToDeformationFieldSource<OutputImageType, double> ConvertorType;
  87 #endif
  88     ConvertorType::Pointer filter= ConvertorType::New();
  89
  90     //Output image info
  91     if (m_ArgsInfo.like_given)
  92       {
  93         typedef itk::ImageFileReader<OutputImageType> ReaderType;
  94         ReaderType::Pointer reader2=ReaderType::New();
  95         reader2->SetFileName(m_ArgsInfo.like_arg);
  96         reader2->Update();
  97
  98         OutputImageType::Pointer image=reader2->GetOutput();
  99         filter->SetOutputParametersFromImage(image);
 100       }
 101     else
 102       {
 103         unsigned int i=0;
 104         if(m_ArgsInfo.origin_given)
 105           {
 106             OutputImageType::PointType origin;
 107             for(i=0;i<Dimension;i++)
 108               origin[i]=m_ArgsInfo.origin_arg[i];
 109             filter->SetOutputOrigin(origin);
 110           }
 111         if (m_ArgsInfo.spacing_given)
 112           {
 113             OutputImageType::SpacingType spacing;
 114             for(i=0;i<Dimension;i++)
 115               spacing[i]=m_ArgsInfo.spacing_arg[i];
 116             filter->SetOutputSpacing(spacing);
 117           }
 118         if (m_ArgsInfo.spacing_given)
 119           {
 120             OutputImageType::SizeType size;
 121             for(i=0;i<Dimension;i++)
 122               size[i]=m_ArgsInfo.size_arg[i];
 123             filter->SetOutputSize(size);
 124           }
 125       }
 126
 127     if (m_Verbose)
 128       {
 129         std::cout<< "Setting output origin to "<<filter->GetOutputOrigin()<<"..."<<std::endl;
 130         std::cout<< "Setting output spacing to "<<filter->GetOutputSpacing()<<"..."<<std::endl;
 131         std::cout<< "Setting output size to "<<filter->GetOutputSize()<<"..."<<std::endl;
 132       }
 133
 134
 135     // -----------------------------------------------
 136     // Transform
 137     // -----------------------------------------------
 138     typedef clitk::BSplineDeformableTransform< double, Dimension, Dimension> TransformType;
 139     TransformType::Pointer transform=TransformType::New();
 140
 141     // Spline orders:  Default is cubic splines
 142     InputImageType::RegionType::SizeType splineOrders ;
 143     splineOrders.Fill(3);
 144     if (m_ArgsInfo.order_given)
 145       for(unsigned int i=0; i<Dimension;i++)
 146         splineOrders[i]=m_ArgsInfo.order_arg[i];
 147     if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 148
 149     // Mask
 150     typedef itk::ImageMaskSpatialObject<  Dimension >   MaskType;
 151     MaskType::Pointer  spatialObjectMask=NULL;
 152     if (m_ArgsInfo.mask_given)
 153       {
 154         typedef itk::Image< unsigned char, Dimension >   ImageMaskType;
 155         typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 156         MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 157         maskReader->SetFileName(m_ArgsInfo.mask_arg);
 158
 159         try
 160           {
 161             maskReader->Update();
 162           }
 163         catch( itk::ExceptionObject & err )
 164           {
 165             std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 166             std::cerr << err << std::endl;
 167             return;
 168           }
 169         if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 170
 171         // Set the image to the spatialObject
 172         spatialObjectMask = MaskType::New();
 173         spatialObjectMask->SetImage( maskReader->GetOutput() );
 174       }
 175
 176
 177     // Samplingfactors
 178     InputImageType::SizeType samplingFactors;
 179     for (unsigned int i=0; i< Dimension; i++)
 180       {
 181         samplingFactors[i]= (int) ( input->GetSpacing()[i]/ filter->GetOutputSpacing()[i]);
 182         if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 183       }
 184
 185
 186     // Set
 187     transform->SetSplineOrders(splineOrders);
 188     transform->SetMask(spatialObjectMask);
 189     transform->SetLUTSamplingFactors(samplingFactors);
 190     transform->SetCoefficientImage(input);
 191     filter->SetTransform(transform);
 192
 193
 194     // -----------------------------------------------
 195     // Update
 196     // -----------------------------------------------
 197     if (m_Verbose)std::cout<< "Converting the BSpline transform..."<<std::endl;
 198     try
 199       {
 200         filter->Update();
 201       }
 202     catch (itk::ExceptionObject)
 203       {
 204         std::cerr<<"Error: Exception thrown during execution convertion filter!"<<std::endl;
 205       }
 206
 207     OutputImageType::Pointer output=filter->GetOutput();
 208
 209
 210     // -----------------------------------------------
 211     // Output
 212     // -----------------------------------------------
 213     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 214     WriterType::Pointer writer = WriterType::New();
 215     writer->SetFileName(m_ArgsInfo.output_arg);
 216     writer->SetInput(output);
 217     writer->Update();
 218
 219   }
 220
 221
 222   //-------------------------------------------------------------------
 223   // Update with the number of dimensions
 224   //-------------------------------------------------------------------
 225   template<>
 226   void
 227   ConvertBSplineDeformableTransformToVFGenericFilter::UpdateWithDim<4>(std::string PixelType, int Components)
 228   {
 229     // Components
 230     if (Components !=3)
 231       {
 232         std::cerr<<"Number of components is "<<Components<<"! Only 3 components is supported."<<std::endl;
 233         return;
 234       }
 235     if (PixelType != "double")
 236       {
 237         std::cerr<<"PixelType is  "<<PixelType<<"! Only double coefficient images are supported."<<std::endl;
 238         std::cerr<<"Reading image as double..."<<std::endl;
 239       }
 240
 241     // ImageTypes
 242     const unsigned int Dimension=4;
 243     const unsigned int SpaceDimension=3;
 244     typedef itk::Vector<double, SpaceDimension> InputPixelType;
 245     typedef itk::Vector<float, SpaceDimension> OutputPixelType;
 246     typedef itk::Image<InputPixelType, Dimension> InputImageType;
 247     typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
 248
 249     // Read the input
 250     typedef itk::ImageFileReader<InputImageType> InputReaderType;
 251     InputReaderType::Pointer reader = InputReaderType::New();
 252     reader->SetFileName( m_InputFileName);
 253     reader->Update();
 254     InputImageType::Pointer input= reader->GetOutput();
 255
 256
 257     // -----------------------------------------------
 258     // Filter
 259     // -----------------------------------------------
 260     typedef clitk::TransformToDeformationFieldSource<OutputImageType, double> ConvertorType;
 261     ConvertorType::Pointer filter= ConvertorType::New();
 262
 263     //Output image info
 264     if (m_ArgsInfo.like_given)
 265       {
 266         typedef itk::ImageFileReader<OutputImageType> ReaderType;
 267         ReaderType::Pointer reader2=ReaderType::New();
 268         reader2->SetFileName(m_ArgsInfo.like_arg);
 269         reader2->Update();
 270
 271         OutputImageType::Pointer image=reader2->GetOutput();
 272         filter->SetOutputParametersFromImage(image);
 273       }
 274     else
 275       {
 276         unsigned int i=0;
 277         if(m_ArgsInfo.origin_given)
 278           {
 279             OutputImageType::PointType origin;
 280             for(i=0;i<Dimension;i++)
 281               origin[i]=m_ArgsInfo.origin_arg[i];
 282             filter->SetOutputOrigin(origin);
 283           }
 284         if (m_ArgsInfo.spacing_given)
 285           {
 286             OutputImageType::SpacingType spacing;
 287             for(i=0;i<Dimension;i++)
 288               spacing[i]=m_ArgsInfo.spacing_arg[i];
 289             filter->SetOutputSpacing(spacing);
 290           }
 291         if (m_ArgsInfo.spacing_given)
 292           {
 293             OutputImageType::SizeType size;
 294             for(i=0;i<Dimension;i++)
 295               size[i]=m_ArgsInfo.size_arg[i];
 296             filter->SetOutputSize(size);
 297           }
 298       }
 299     //Output image info
 300     if (m_Verbose)
 301       {
 302         std::cout<< "Setting output origin to "<<filter->GetOutputOrigin()<<"..."<<std::endl;
 303         std::cout<< "Setting output spacing to "<<filter->GetOutputSpacing()<<"..."<<std::endl;
 304         std::cout<< "Setting output size to "<<filter->GetOutputSize()<<"..."<<std::endl;
 305       }
 306
 307
 308     // -----------------------------------------------
 309     // Transform
 310     // -----------------------------------------------
 311     typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform< double, Dimension, Dimension > TransformType;
 312     TransformType::Pointer transform=TransformType::New();
 313     transform->SetTransformShape(m_ArgsInfo.shape_arg);
 314
 315     // Spline orders:  Default is cubic splines
 316     InputImageType::RegionType::SizeType splineOrders ;
 317     splineOrders.Fill(3);
 318     if (m_ArgsInfo.order_given)
 319       for(unsigned int i=0; i<Dimension;i++)
 320         splineOrders[i]=m_ArgsInfo.order_arg[i];
 321     if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 322
 323     // Mask
 324     typedef itk::ImageMaskSpatialObject<  Dimension >   MaskType;
 325     MaskType::Pointer  spatialObjectMask=NULL;
 326     if (m_ArgsInfo.mask_given)
 327       {
 328         typedef itk::Image< unsigned char, Dimension >   ImageMaskType;
 329         typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 330         MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 331         maskReader->SetFileName(m_ArgsInfo.mask_arg);
 332
 333         try
 334           {
 335             maskReader->Update();
 336           }
 337         catch( itk::ExceptionObject & err )
 338           {
 339             std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 340             std::cerr << err << std::endl;
 341             return;
 342           }
 343         if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 344
 345         // Set the image to the spatialObject
 346         spatialObjectMask = MaskType::New();
 347         spatialObjectMask->SetImage( maskReader->GetOutput() );
 348       }
 349
 350
 351     // Samplingfactors
 352     InputImageType::SizeType samplingFactors;
 353     for (unsigned int i=0; i< Dimension; i++)
 354       {
 355         samplingFactors[i]= (int) ( input->GetSpacing()[i]/ filter->GetOutputSpacing()[i]);
 356         if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 357       }
 358     if( !(m_ArgsInfo.shape_arg%2) )samplingFactors[Dimension-1]=5;
 359
 360     // Set
 361     transform->SetSplineOrders(splineOrders);
 362     transform->SetMask(spatialObjectMask);
 363     transform->SetLUTSamplingFactors(samplingFactors);
 364     transform->SetCoefficientImage(input);
 365     filter->SetTransform(transform);
 366
 367
 368     // -----------------------------------------------
 369     // Update
 370     // -----------------------------------------------
 371     if (m_Verbose)std::cout<< "Converting the BSpline transform..."<<std::endl;
 372     try
 373       {
 374         filter->Update();
 375       }
 376     catch (itk::ExceptionObject)
 377       {
 378         std::cerr<<"Error: Exception thrown during execution convertion filter!"<<std::endl;
 379       }
 380
 381     OutputImageType::Pointer output=filter->GetOutput();
 382
 383
 384     // -----------------------------------------------
 385     // Output
 386     // -----------------------------------------------
 387     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 388     WriterType::Pointer writer = WriterType::New();
 389     writer->SetFileName(m_ArgsInfo.output_arg);
 390     writer->SetInput(output);
 391     writer->Update();
 392
 393   }
 394
 395
 396
 397   //-----------------------------------------------------------
 398   // Update
 399   //-----------------------------------------------------------
 400   void ConvertBSplineDeformableTransformToVFGenericFilter::Update()
 401   {
 402
 403     // Read the Dimension and PixelType
 404     int Dimension, Components;
 405     std::string PixelType;
 406     ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components);
 407
 408     // Call UpdateWithDim
 409     //if(Dimension==2) UpdateWithDim<2>(PixelType, Components);
 410     if(Dimension==3) UpdateWithDim<3>(PixelType, Components);
 411     else if (Dimension==4) UpdateWithDim<4>(PixelType, Components);
 412     else
 413       {
 414         std::cout<<"Error, Only for 3 Dimensions!!!"<<std::endl ;
 415         return;
 416       }
 417   }
 418
 419 } //end clitk
 420
 421 #endif  //#define clitkConvertBSplineDeformableTransformToVFGenericFilter_cxx