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 // clitk include
  31 #include "clitkIO.h"
  32 #include "clitkCommon.h"
  33 #include "clitkImageCommon.h"
  34 #include "clitkConvertBSplineDeformableTransformToVF_ggo.h"
  35 #include "clitkBSplineDeformableTransform.h"
  36 #include "clitkTransformToDeformationFieldSource.h"
  37 #include "clitkShapedBLUTSpatioTemporalDeformableTransform.h"
  38 #include "itkImageMaskSpatialObject.h"
  39
  40 #include "clitkConvertBSplineDeformableTransformToVFGenericFilter.h"
  41 #include "clitkVectorImageToImageFilter.h"
  42 #if ITK_VERSION_MAJOR >= 4
  43 #include "itkTransformToDisplacementFieldSource.h"
  44 #else
  45 #include "itkTransformToDeformationFieldSource.h"
  46 #endif
  47 #include "itkBSplineDeformableTransform.h"
  48
  49
  50 namespace clitk
  51 {
  52
  53
  54   //-----------------------------------------------------------
  55   // Constructor
  56   //-----------------------------------------------------------
  57   ConvertBSplineDeformableTransformToVFGenericFilter::ConvertBSplineDeformableTransformToVFGenericFilter()
  58   {
  59     m_Verbose=false;
  60     m_InputFileName="";
  61   }
  62
  63   //-------------------------------------------------------------------
  64   // Update with the number of dimensions
  65   //-------------------------------------------------------------------
  66   template<>
  67   void
  68   ConvertBSplineDeformableTransformToVFGenericFilter::UpdateWithDim<3>(std::string PixelType, int Components)
  69 {
  70     // Components
  71     if (Components !=3)
  72     {
  73         std::cerr<<"Number of components is "<<Components<<"! Only 3 components is supported."<<std::endl;
  74         return;
  75     }
  76     if (PixelType != "double")
  77     {
  78         std::cerr<<"PixelType is  "<<PixelType<<"! Only double coefficient images are supported."<<std::endl;
  79         std::cerr<<"Reading image as double..."<<std::endl;
  80     }
  81
  82     // ImageTypes
  83     const unsigned int Dimension=3;
  84     typedef itk::Vector<double, Dimension> InputPixelType;
  85     typedef itk::Vector<float, Dimension> OutputPixelType;
  86     typedef itk::Image<InputPixelType, Dimension> InputImageType;
  87     typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
  88
  89     // Read the input
  90     typedef itk::ImageFileReader<InputImageType> InputReaderType;
  91     InputReaderType::Pointer reader = InputReaderType::New();
  92     reader->SetFileName( m_InputFileName);
  93     reader->Update();
  94     InputImageType::Pointer input= reader->GetOutput();
  95
  96
  97     // -----------------------------------------------
  98     // Filter
  99     // -----------------------------------------------
 100 #if ITK_VERSION_MAJOR >= 4
 101     typedef itk::TransformToDisplacementFieldSource<OutputImageType, double> ConvertorType;
 102 #else
 103     typedef itk::TransformToDeformationFieldSource<OutputImageType, double> ConvertorType;
 104 #endif
 105     ConvertorType::Pointer filter= ConvertorType::New();
 106
 107     //Output image info
 108     if (m_ArgsInfo.like_given)
 109     {
 110 /*        typedef itk::ImageFileReader<OutputImageType> ReaderType;
 111         ReaderType::Pointer reader2=ReaderType::New();
 112         reader2->SetFileName(m_ArgsInfo.like_arg);
 113         reader2->Update();
 114
 115         OutputImageType::Pointer image=reader2->GetOutput();
 116         filter->SetOutputParametersFromImage(image);*/
 117
 118         typedef itk::ImageIOBase ImageIOType;
 119         typename ImageIOType::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(m_ArgsInfo.like_arg, itk::ImageIOFactory::ReadMode);
 120         imageIO->SetFileName(m_ArgsInfo.like_arg);
 121         imageIO->ReadImageInformation();
 122
 123         typename ConvertorType::SizeType output_size;
 124         typename ConvertorType::SpacingType output_spacing;
 125         typename ConvertorType::OriginType output_origin;
 126         typename ConvertorType::DirectionType output_direction;
 127         for (unsigned int i = 0; i < Dimension; i++) {
 128           output_size[i] = imageIO->GetDimensions(i);
 129           output_spacing[i] = imageIO->GetSpacing(i);
 130           output_origin[i] = imageIO->GetOrigin(i);
 131           for (unsigned int j = 0; j < Dimension; j++)
 132             output_direction[i][j] = imageIO->GetDirection(i)[j];
 133         }
 134
 135         filter->SetOutputOrigin(output_origin);
 136         filter->SetOutputSpacing(output_spacing);
 137         filter->SetOutputSize(output_size);
 138         filter->SetOutputDirection(output_direction);
 139     }
 140     else
 141     {
 142         unsigned int i=0;
 143         if (m_ArgsInfo.origin_given)
 144         {
 145             OutputImageType::PointType origin;
 146             for (i=0;i<Dimension;i++)
 147                 origin[i]=m_ArgsInfo.origin_arg[i];
 148             filter->SetOutputOrigin(origin);
 149         }
 150         if (m_ArgsInfo.spacing_given)
 151         {
 152             OutputImageType::SpacingType spacing;
 153             for (i=0;i<Dimension;i++)
 154                 spacing[i]=m_ArgsInfo.spacing_arg[i];
 155             filter->SetOutputSpacing(spacing);
 156         }
 157         if (m_ArgsInfo.spacing_given)
 158         {
 159             OutputImageType::SizeType size;
 160             for (i=0;i<Dimension;i++)
 161                 size[i]=m_ArgsInfo.size_arg[i];
 162             filter->SetOutputSize(size);
 163         }
 164     }
 165
 166     if (m_Verbose)
 167     {
 168         std::cout<< "Setting output origin to "<<filter->GetOutputOrigin()<<"..."<<std::endl;
 169         std::cout<< "Setting output spacing to "<<filter->GetOutputSpacing()<<"..."<<std::endl;
 170         std::cout<< "Setting output size to "<<filter->GetOutputSize()<<"..."<<std::endl;
 171     }
 172
 173
 174     // -----------------------------------------------
 175     // Transform
 176     // -----------------------------------------------
 177     typedef clitk::BSplineDeformableTransform< double, Dimension, Dimension> BLUTTransformType;
 178     BLUTTransformType::Pointer blut_transform=BLUTTransformType::New();
 179
 180     typedef itk::BSplineDeformableTransform< double, Dimension, Dimension> ITKTransformType;
 181     ITKTransformType::Pointer itk_transform=ITKTransformType::New();
 182
 183     typedef itk::Transform< double, Dimension, Dimension> GenericTransformType;
 184     typename GenericTransformType::Pointer transform;
 185
 186     // Mask
 187     typedef itk::ImageMaskSpatialObject<  Dimension >   MaskType;
 188     MaskType::Pointer  spatialObjectMask=NULL;
 189     if (m_ArgsInfo.mask_given)
 190     {
 191         typedef itk::Image< unsigned char, Dimension >   ImageMaskType;
 192         typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 193         MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 194         maskReader->SetFileName(m_ArgsInfo.mask_arg);
 195
 196         try
 197         {
 198             maskReader->Update();
 199         }
 200         catch ( itk::ExceptionObject & err )
 201         {
 202             std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 203             std::cerr << err << std::endl;
 204             return;
 205         }
 206         if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 207
 208         // Set the image to the spatialObject
 209         spatialObjectMask = MaskType::New();
 210         spatialObjectMask->SetImage( maskReader->GetOutput() );
 211         blut_transform->SetMask(spatialObjectMask);
 212     }
 213
 214     if (m_ArgsInfo.type_arg != 0 ) { // using BLUT
 215       // Spline orders:  Default is cubic splines
 216         if (m_Verbose) {
 217           std::cout << "Using clitk::BLUT." << std::endl;
 218           std::cout << "Setting spline orders and sampling factors." << std::endl;
 219         }
 220       InputImageType::RegionType::SizeType splineOrders ;
 221       splineOrders.Fill(3);
 222       if (m_ArgsInfo.order_given)
 223           for (unsigned int i=0; i<Dimension;i++)
 224               splineOrders[i]=m_ArgsInfo.order_arg[i];
 225       if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 226
 227       // Samplingfactors
 228       InputImageType::SizeType samplingFactors;
 229       for (unsigned int i=0; i< Dimension; i++)
 230       {
 231           samplingFactors[i]= (int) ( input->GetSpacing()[i]/ filter->GetOutputSpacing()[i]);
 232           if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 233       }
 234       blut_transform->SetSplineOrders(splineOrders);
 235       blut_transform->SetLUTSamplingFactors(samplingFactors);
 236       blut_transform->SetCoefficientImage(input);
 237
 238       transform = blut_transform;
 239     }
 240     else { // using ITK transform
 241         if (m_Verbose) {
 242           std::cout << "Using itk::BSpline" << std::endl;
 243           std::cout << "Extracting components from input coefficient image and creating one coefficient image per-component" << std::endl;
 244         }
 245
 246         typedef  double PixelType;
 247         typedef itk::Vector<PixelType, Dimension> CoefficientType;
 248         typedef itk::Image<CoefficientType, Dimension> CoefficientImageType;
 249         typedef clitk::VectorImageToImageFilter<CoefficientImageType, typename ITKTransformType::ImageType> FilterType;
 250         typename FilterType::Pointer component_filter[Dimension];
 251
 252         typename ITKTransformType::ImagePointer coefficient_images[Dimension];
 253         for (unsigned int i=0; i < Dimension; i++) {
 254             component_filter[i] = FilterType::New();
 255             component_filter[i]->SetInput(input);
 256             component_filter[i]->SetComponentIndex(i);
 257             component_filter[i]->Update();
 258             coefficient_images[i] = component_filter[i]->GetOutput();
 259         }
 260         itk_transform->SetCoefficientImage(coefficient_images);
 261
 262         transform = itk_transform;
 263     }
 264
 265     filter->SetTransform(transform);
 266
 267
 268     // -----------------------------------------------
 269     // Update
 270     // -----------------------------------------------
 271     if (m_Verbose)std::cout<< "Converting the BSpline transform..."<<std::endl;
 272     try
 273     {
 274         filter->Update();
 275     }
 276     catch (itk::ExceptionObject)
 277     {
 278         std::cerr<<"Error: Exception thrown during execution convertion filter!"<<std::endl;
 279     }
 280
 281     OutputImageType::Pointer output=filter->GetOutput();
 282
 283
 284     // -----------------------------------------------
 285     // Output
 286     // -----------------------------------------------
 287     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 288     WriterType::Pointer writer = WriterType::New();
 289     writer->SetFileName(m_ArgsInfo.output_arg);
 290     writer->SetInput(output);
 291     writer->Update();
 292
 293   }
 294
 295 /*
 296   //-------------------------------------------------------------------
 297   // Update with the number of dimensions
 298   //-------------------------------------------------------------------
 299   template<>
 300 void
 301 ConvertBSplineDeformableTransformToVFGenericFilter::UpdateWithDim<4>(std::string PixelType, int Components)
 302 {
 303     // Components
 304     if (Components !=3)
 305     {
 306         std::cerr<<"Number of components is "<<Components<<"! Only 3 components is supported."<<std::endl;
 307         return;
 308     }
 309     if (PixelType != "double")
 310     {
 311         std::cerr<<"PixelType is  "<<PixelType<<"! Only double coefficient images are supported."<<std::endl;
 312         std::cerr<<"Reading image as double..."<<std::endl;
 313     }
 314
 315     // ImageTypes
 316     const unsigned int Dimension=4;
 317     const unsigned int SpaceDimension=3;
 318     typedef itk::Vector<double, SpaceDimension> InputPixelType;
 319     typedef itk::Vector<float, SpaceDimension> OutputPixelType;
 320     typedef itk::Image<InputPixelType, Dimension> InputImageType;
 321     typedef itk::Image<OutputPixelType, Dimension> OutputImageType;
 322
 323     // Read the input
 324     typedef itk::ImageFileReader<InputImageType> InputReaderType;
 325     InputReaderType::Pointer reader = InputReaderType::New();
 326     reader->SetFileName( m_InputFileName);
 327     reader->Update();
 328     InputImageType::Pointer input= reader->GetOutput();
 329
 330
 331     // -----------------------------------------------
 332     // Filter
 333     // -----------------------------------------------
 334     typedef clitk::TransformToDeformationFieldSource<OutputImageType, double> ConvertorType;
 335     ConvertorType::Pointer filter= ConvertorType::New();
 336
 337     //Output image info
 338     if (m_ArgsInfo.like_given)
 339     {
 340         typedef itk::ImageFileReader<OutputImageType> ReaderType;
 341         ReaderType::Pointer reader2=ReaderType::New();
 342         reader2->SetFileName(m_ArgsInfo.like_arg);
 343         reader2->Update();
 344
 345         OutputImageType::Pointer image=reader2->GetOutput();
 346         filter->SetOutputParametersFromImage(image);
 347     }
 348     else
 349     {
 350         unsigned int i=0;
 351         if (m_ArgsInfo.origin_given)
 352         {
 353             OutputImageType::PointType origin;
 354             for (i=0;i<Dimension;i++)
 355                 origin[i]=m_ArgsInfo.origin_arg[i];
 356             filter->SetOutputOrigin(origin);
 357         }
 358         if (m_ArgsInfo.spacing_given)
 359         {
 360             OutputImageType::SpacingType spacing;
 361             for (i=0;i<Dimension;i++)
 362                 spacing[i]=m_ArgsInfo.spacing_arg[i];
 363             filter->SetOutputSpacing(spacing);
 364         }
 365         if (m_ArgsInfo.spacing_given)
 366         {
 367             OutputImageType::SizeType size;
 368             for (i=0;i<Dimension;i++)
 369                 size[i]=m_ArgsInfo.size_arg[i];
 370             filter->SetOutputSize(size);
 371         }
 372     }
 373     //Output image info
 374     if (m_Verbose)
 375     {
 376         std::cout<< "Setting output origin to "<<filter->GetOutputOrigin()<<"..."<<std::endl;
 377         std::cout<< "Setting output spacing to "<<filter->GetOutputSpacing()<<"..."<<std::endl;
 378         std::cout<< "Setting output size to "<<filter->GetOutputSize()<<"..."<<std::endl;
 379     }
 380
 381
 382     // -----------------------------------------------
 383     // Transform
 384     // -----------------------------------------------
 385     typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform< double, Dimension, Dimension > TransformType;
 386     TransformType::Pointer transform=TransformType::New();
 387     transform->SetTransformShape(m_ArgsInfo.shape_arg);
 388
 389     // Spline orders:  Default is cubic splines
 390     InputImageType::RegionType::SizeType splineOrders ;
 391     splineOrders.Fill(3);
 392     if (m_ArgsInfo.order_given)
 393         for (unsigned int i=0; i<Dimension;i++)
 394             splineOrders[i]=m_ArgsInfo.order_arg[i];
 395     if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 396
 397     // Mask
 398     typedef itk::ImageMaskSpatialObject<  Dimension >   MaskType;
 399     MaskType::Pointer  spatialObjectMask=NULL;
 400     if (m_ArgsInfo.mask_given)
 401     {
 402         typedef itk::Image< unsigned char, Dimension >   ImageMaskType;
 403         typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 404         MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 405         maskReader->SetFileName(m_ArgsInfo.mask_arg);
 406
 407         try
 408         {
 409             maskReader->Update();
 410         }
 411         catch ( itk::ExceptionObject & err )
 412         {
 413             std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 414             std::cerr << err << std::endl;
 415             return;
 416         }
 417         if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 418
 419         // Set the image to the spatialObject
 420         spatialObjectMask = MaskType::New();
 421         spatialObjectMask->SetImage( maskReader->GetOutput() );
 422     }
 423
 424
 425     // Samplingfactors
 426     InputImageType::SizeType samplingFactors;
 427     for (unsigned int i=0; i< Dimension; i++)
 428     {
 429         samplingFactors[i]= (int) ( input->GetSpacing()[i]/ filter->GetOutputSpacing()[i]);
 430         if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 431     }
 432     if ( !(m_ArgsInfo.shape_arg%2) )samplingFactors[Dimension-1]=5;
 433
 434     // Set
 435     transform->SetSplineOrders(splineOrders);
 436     transform->SetMask(spatialObjectMask);
 437     transform->SetLUTSamplingFactors(samplingFactors);
 438     transform->SetCoefficientImage(input);
 439     filter->SetTransform(transform);
 440
 441
 442     // -----------------------------------------------
 443     // Update
 444     // -----------------------------------------------
 445     if (m_Verbose)std::cout<< "Converting the BSpline transform..."<<std::endl;
 446     try
 447     {
 448         filter->Update();
 449     }
 450     catch (itk::ExceptionObject)
 451     {
 452         std::cerr<<"Error: Exception thrown during execution convertion filter!"<<std::endl;
 453     }
 454
 455     OutputImageType::Pointer output=filter->GetOutput();
 456
 457
 458     // -----------------------------------------------
 459     // Output
 460     // -----------------------------------------------
 461     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 462     WriterType::Pointer writer = WriterType::New();
 463     writer->SetFileName(m_ArgsInfo.output_arg);
 464     writer->SetInput(output);
 465     writer->Update();
 466
 467 }
 468 */
 469   //-----------------------------------------------------------
 470   // Update
 471   //-----------------------------------------------------------
 472   void ConvertBSplineDeformableTransformToVFGenericFilter::Update()
 473   {
 474
 475     // Read the Dimension and PixelType
 476     int Dimension, Components;
 477     std::string PixelType;
 478     ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType, Components);
 479
 480     // Call UpdateWithDim
 481     //if(Dimension==2) UpdateWithDim<2>(PixelType, Components);
 482     if(Dimension==3) UpdateWithDim<3>(PixelType, Components);
 483     //else if (Dimension==4) UpdateWithDim<4>(PixelType, Components);
 484     else
 485       {
 486         std::cout<<"Error, Only for 3 Dimensions!!!"<<std::endl ;
 487         return;
 488       }
 489   }
 490
 491 } //end clitk
 492
 493 #endif  //#define clitkConvertBSplineDeformableTransformToVFGenericFilter_cxx