tools/clitkBSplineCoefficientsToValuesGenericFilter.txx

   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://oncora1.lyon.fnclcc.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 clitkBSplineCoefficientsToValuesGenericFilter_txx
  19 #define clitkBSplineCoefficientsToValuesGenericFilter_txx
  20 #include <itkAffineTransform.h>
  21
  22 /* =================================================
  23  * @file   clitkBSplineCoefficientsToValuesGenericFilter.txx
  24  * @author
  25  * @date
  26  *
  27  * @brief
  28  *
  29  ===================================================*/
  30
  31
  32 namespace clitk
  33 {
  34
  35   //-------------------------------------------------------------------
  36   // Update with the number of dimensions
  37   //-------------------------------------------------------------------
  38   template<unsigned int Dimension>
  39   void
  40   BSplineCoefficientsToValuesGenericFilter::UpdateWithDim(std::string PixelType, unsigned int Components)
  41   {
  42     if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<<Components<<" component(s) of "<<  PixelType<<"..."<<std::endl;
  43
  44     if (Components==1)
  45       {
  46         //      if(PixelType == "short"){
  47         //        if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
  48         //        UpdateWithDimAndPixelType<Dimension, signed short>();
  49         //      }
  50         //    else if(PixelType == "unsigned_short"){
  51         //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
  52         //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
  53         //     }
  54
  55         //      else if (PixelType == "unsigned_char"){
  56         //        if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
  57         //        UpdateWithDimAndPixelType<Dimension, unsigned char>();
  58         //      }
  59
  60         //     else if (PixelType == "char"){
  61         //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
  62         //       UpdateWithDimAndPixelType<Dimension, signed char>();
  63         //     }
  64         //      else {
  65           if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
  66           UpdateWithDimAndPixelType<Dimension, float>();
  67           //    }
  68       }
  69
  70    else
  71     if (Components==3)
  72       {
  73         if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and 3D float (DVF)" << std::endl;
  74         UpdateWithDimAndVectorType<Dimension, itk::Vector<float, 3> >();
  75       }
  76
  77     else std::cerr<<"Number of components is "<<Components<<", not supported!"<<std::endl;
  78
  79   }
  80
  81
  82   //-------------------------------------------------------------------
  83   // Update with the number of dimensions and the pixeltype
  84   //-------------------------------------------------------------------
  85   template <unsigned int Dimension, class  PixelType>
  86   void
  87   BSplineCoefficientsToValuesGenericFilter::UpdateWithDimAndPixelType()
  88   {
  89
  90     // ImageTypes
  91     typedef itk::Image<PixelType, Dimension> InputImageType;
  92     typedef itk::Image<PixelType, Dimension> OutputImageType;
  93
  94     // Read the input
  95     typedef itk::ImageFileReader<InputImageType> InputReaderType;
  96     typename InputReaderType::Pointer reader = InputReaderType::New();
  97     reader->SetFileName( m_InputFileName);
  98     reader->Update();
  99     typename InputImageType::Pointer input= reader->GetOutput();
 100
 101     //Filter
 102     typedef  itk::ResampleImageFilter< InputImageType,OutputImageType >  ResampleFilterType;
 103     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 104
 105     // Properties
 106     if (m_ArgsInfo.like_given)
 107       {
 108         typename InputReaderType::Pointer likeReader=InputReaderType::New();
 109         likeReader->SetFileName(m_ArgsInfo.like_arg);
 110         likeReader->Update();
 111         resampler->SetOutputParametersFromImage(likeReader->GetOutput());
 112         if (m_Verbose)std::cout<<"Resampling output like "<<m_ArgsInfo.like_arg<<"..."<<std::endl;
 113       }
 114     else
 115       {
 116         // Size
 117         typename OutputImageType::SizeType outputSize;
 118         if (m_ArgsInfo.size_given)
 119           {
 120             for(unsigned int i=0; i< Dimension; i++)
 121               outputSize[i]=m_ArgsInfo.size_arg[i];
 122           }
 123         else outputSize=input->GetLargestPossibleRegion().GetSize();
 124         if (m_Verbose)std::cout<<"Setting the size to "<<outputSize<<"..."<<std::endl;
 125
 126         // Spacing
 127         typename OutputImageType::SpacingType outputSpacing;
 128         if (m_ArgsInfo.spacing_given)
 129           {
 130             for(unsigned int i=0; i< Dimension; i++)
 131               outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
 132           }
 133         else outputSpacing=input->GetSpacing();
 134         if (m_Verbose)std::cout<<"Setting the spacing to "<<outputSpacing<<"..."<<std::endl;
 135
 136         // Origin
 137         typename OutputImageType::PointType outputOrigin;
 138         if (m_ArgsInfo.origin_given)
 139           {
 140             for(unsigned int i=0; i< Dimension; i++)
 141               outputOrigin[i]=m_ArgsInfo.origin_arg[i];
 142           }
 143         else outputOrigin=input->GetOrigin();
 144         if (m_Verbose)std::cout<<"Setting the origin to "<<outputOrigin<<"..."<<std::endl;
 145
 146         // Set
 147         resampler->SetSize( outputSize );
 148         resampler->SetOutputSpacing( outputSpacing );
 149         resampler->SetOutputOrigin(  outputOrigin );
 150
 151       }
 152
 153     // Interp : coeff
 154     typedef itk::BSplineResampleImageFunction<InputImageType, double> InterpolatorType;
 155     typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
 156     interpolator->SetSplineOrder(m_ArgsInfo.order_arg);
 157
 158     // Set
 159     resampler->SetInterpolator( interpolator);
 160     resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 161     resampler->SetInput( input );
 162
 163     // Update
 164     try
 165       {
 166         resampler->Update();
 167       }
 168     catch(itk::ExceptionObject)
 169       {
 170         std::cerr<<"Error resampling the image"<<std::endl;
 171       }
 172
 173     typename OutputImageType::Pointer output = resampler->GetOutput();
 174
 175     // Output
 176     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 177     typename WriterType::Pointer writer = WriterType::New();
 178     writer->SetFileName(m_ArgsInfo.output_arg);
 179     writer->SetInput(output);
 180     writer->Update();
 181   }
 182
 183
 184   //-------------------------------------------------------------------
 185   // Update with the number of dimensions and the Vectortype
 186   //-------------------------------------------------------------------
 187   template <unsigned int Dimension, class  PixelType>
 188   void
 189   BSplineCoefficientsToValuesGenericFilter::UpdateWithDimAndVectorType()
 190   {
 191
 192     // ImageTypes
 193     typedef itk::Image<PixelType, Dimension> InputImageType;
 194     typedef itk::Image<PixelType, Dimension> OutputImageType;
 195
 196     // Read the input
 197     typedef itk::ImageFileReader<InputImageType> InputReaderType;
 198     typename InputReaderType::Pointer reader = InputReaderType::New();
 199     reader->SetFileName( m_InputFileName);
 200     reader->Update();
 201     typename InputImageType::Pointer input= reader->GetOutput();
 202
 203     //Filter
 204     typedef  itk::VectorResampleImageFilter< InputImageType,OutputImageType >  ResampleFilterType;
 205     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 206
 207     // Properties
 208     if (m_ArgsInfo.like_given)
 209       {
 210         typename InputReaderType::Pointer likeReader=InputReaderType::New();
 211         likeReader->SetFileName(m_ArgsInfo.like_arg);
 212         likeReader->Update();
 213
 214         // Set
 215         resampler->SetSize( likeReader->GetOutput()->GetLargestPossibleRegion().GetSize() );
 216         resampler->SetOutputSpacing(likeReader->GetOutput()->GetSpacing() );
 217         resampler->SetOutputOrigin(  likeReader->GetOutput()->GetOrigin() );
 218         if (m_Verbose)std::cout<<"Resampling output like "<<m_ArgsInfo.like_arg<<"..."<<std::endl;
 219       }
 220     else
 221       {
 222         // Size
 223         typename OutputImageType::SizeType outputSize;
 224         if (m_ArgsInfo.size_given)
 225           {
 226             for(unsigned int i=0; i< Dimension; i++)
 227               outputSize[i]=m_ArgsInfo.size_arg[i];
 228           }
 229         else outputSize=input->GetLargestPossibleRegion().GetSize();
 230         if (m_Verbose)std::cout<<"Setting the size to "<<outputSize<<"..."<<std::endl;
 231
 232         // Spacing
 233         typename OutputImageType::SpacingType outputSpacing;
 234         if (m_ArgsInfo.spacing_given)
 235           {
 236             for(unsigned int i=0; i< Dimension; i++)
 237               outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
 238           }
 239         else outputSpacing=input->GetSpacing();
 240         if (m_Verbose)std::cout<<"Setting the spacing to "<<outputSpacing<<"..."<<std::endl;
 241
 242         // Origin
 243         typename OutputImageType::PointType outputOrigin;
 244         if (m_ArgsInfo.origin_given)
 245           {
 246             for(unsigned int i=0; i< Dimension; i++)
 247               outputOrigin[i]=m_ArgsInfo.origin_arg[i];
 248           }
 249         else outputOrigin=input->GetOrigin();
 250         if (m_Verbose)std::cout<<"Setting the origin to "<<outputOrigin<<"..."<<std::endl;
 251
 252         // Set
 253         resampler->SetSize( outputSize );
 254         resampler->SetOutputSpacing( outputSpacing );
 255         resampler->SetOutputOrigin(  outputOrigin );
 256
 257       }
 258
 259     // Interp : coeff
 260     typedef clitk::VectorBSplineResampleImageFunction<InputImageType, double> InterpolatorType;
 261     typename InterpolatorType::Pointer interpolator=InterpolatorType::New();
 262     interpolator->SetSplineOrder(m_ArgsInfo.order_arg);
 263
 264
 265     // Set
 266     resampler->SetInterpolator( interpolator);
 267     resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
 268     resampler->SetInput( input );
 269
 270     // Update
 271     try
 272       {
 273         resampler->Update();
 274       }
 275     catch(itk::ExceptionObject)
 276       {
 277         std::cerr<<"Error resampling the image"<<std::endl;
 278       }
 279
 280     typename OutputImageType::Pointer output = resampler->GetOutput();
 281     //writeImage<OutputImageType>(output, "/home/jef/tmp/output1.mhd", true);
 282
 283
 284     // Matrix Transform
 285     if(m_ArgsInfo.matrix_given)
 286       {
 287         typedef itk::TransformToDeformationFieldSource<OutputImageType, double> ConvertorType;
 288         typename   ConvertorType::Pointer filter= ConvertorType::New();
 289         filter->SetOutputParametersFromImage(output);
 290
 291         typedef itk::AffineTransform<double, Dimension> TransformType;
 292         typename TransformType::Pointer transform =TransformType::New();
 293         itk::Matrix<double, Dimension+1, Dimension+1> homMatrix= ReadMatrix<Dimension>( m_ArgsInfo.matrix_arg);
 294         itk::Matrix<double, Dimension, Dimension> matrix =GetRotationalPartMatrix( homMatrix);
 295         itk::Vector<double, Dimension> offset= GetTranslationPartMatrix( homMatrix);
 296         transform->SetMatrix(matrix);
 297         transform->SetOffset(offset);
 298         filter->SetTransform(transform);
 299         filter->Update();
 300         typename OutputImageType::Pointer output2=filter->GetOutput();
 301         //writeImage<OutputImageType>(output2, "/home/jef/tmp/output2.mhd", true);
 302
 303         // Add
 304         typedef itk::AddImageFilter< OutputImageType, OutputImageType, OutputImageType > AddType;
 305         typename AddType::Pointer adder= AddType::New();
 306         adder->SetInput(output);
 307         adder->SetInput(1,output2);
 308         adder->Update();
 309         output=adder->GetOutput();
 310       }
 311
 312
 313     // Output
 314     typedef itk::ImageFileWriter<OutputImageType> WriterType;
 315     typename WriterType::Pointer writer = WriterType::New();
 316     writer->SetFileName(m_ArgsInfo.output_arg);
 317     writer->SetInput(output);
 318     writer->Update();
 319   }
 320
 321 }//end clitk
 322
 323 #endif //#define clitkBSplineCoefficientsToValuesGenericFilter_txx