registration/clitkPointTrajectoryGenericFilter.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 clitkPointTrajectoryGenericFilter_cxx
  19 #define clitkPointTrajectoryGenericFilter_cxx
  20
  21 /* =================================================
  22  * @file   clitkPointTrajectoryGenericFilter.cxx
  23  * @author
  24  * @date
  25  *
  26  * @brief
  27  *
  28  ===================================================*/
  29
  30 #include "clitkPointTrajectoryGenericFilter.h"
  31
  32
  33 namespace clitk
  34 {
  35
  36
  37   //-----------------------------------------------------------
  38   // Constructor
  39   //-----------------------------------------------------------
  40   PointTrajectoryGenericFilter::PointTrajectoryGenericFilter()
  41   {
  42     m_Verbose=false;
  43     //    m_InputFileName="";
  44   }
  45
  46
  47   //-----------------------------------------------------------
  48   // Update
  49   //-----------------------------------------------------------
  50   void PointTrajectoryGenericFilter::Update()
  51   {
  52     // ImageTypes
  53     const unsigned int ImageDimension=4;
  54     const unsigned int SpaceDimension=3;
  55     typedef itk::Vector<double, SpaceDimension> CoefficientPixelType;
  56     typedef itk::Vector<float, SpaceDimension> VectorPixelType;
  57     typedef itk::Image<CoefficientPixelType, ImageDimension> CoefficientImageType;
  58     typedef itk::Image<VectorPixelType, ImageDimension> VectorFieldType;
  59
  60
  61     // -----------------------------------------------
  62     // Reference Point List 3D
  63     // -----------------------------------------------
  64     typedef itk::Point<double, SpaceDimension> SpacePointType;
  65     typedef clitk::List<SpacePointType> PointListType;
  66     PointListType referencePointList;
  67     if (m_Verbose) std::cout<<"Reference point list:"<<std::endl;
  68     referencePointList=PointListType(m_ArgsInfo.ref_arg, m_Verbose);
  69
  70
  71     // -----------------------------------------------
  72     // Transform: based on points, 4DVF, spatio-Temporal transform
  73     // -----------------------------------------------
  74     typedef itk::Transform<double, 4, 4> TransformType;
  75     TransformType::Pointer transform;
  76     switch (m_ArgsInfo.transform_arg)
  77       {
  78         // ==========================
  79         // List of points
  80         // ==========================
  81       case 0:
  82         {
  83           //-----------------------------
  84           // Input point lists
  85           //-----------------------------
  86           typedef itk::Point<double, SpaceDimension> PointType;
  87           typedef clitk::List<PointType> PointListType;
  88           typedef clitk::Lists<PointType> PointListsType;
  89           PointListsType inputPointLists, sortedPointLists;
  90
  91           // Read the lists
  92           for (unsigned int i=0; i<m_ArgsInfo.points_given; i++)
  93             inputPointLists.push_back(PointListType(m_ArgsInfo.points_arg[i], m_Verbose) );
  94
  95           // Convert/sort the lists
  96           sortedPointLists.resize(inputPointLists[0].size());
  97           for (unsigned int i=0; i<inputPointLists[0].size(); i++)
  98             {
  99               sortedPointLists[i].push_back(referencePointList[i]);
 100               for (unsigned int j=0; j<inputPointLists.size(); j++)
 101                 {
 102                   sortedPointLists[i].push_back(inputPointLists[j][i]);
 103                 }
 104             }
 105
 106           // Point list Transform
 107           typedef clitk::PointListTransform<double, ImageDimension,ImageDimension> PointListTransformType;
 108           PointListTransformType::Pointer pointListTransform=PointListTransformType::New();
 109           pointListTransform->SetPointLists(sortedPointLists);
 110
 111           // Vector Interpolator
 112           typedef PointListTransformType::PointListImageType PointListImageType;
 113           typedef clitk::GenericVectorInterpolator<args_info_clitkPointTrajectory,PointListImageType, double> GenericVectorInterpolatorType;
 114           GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
 115           genericInterpolator->SetArgsInfo(m_ArgsInfo);
 116           typedef itk::VectorInterpolateImageFunction<PointListImageType, double> InterpolatorType;
 117           InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
 118           pointListTransform->SetInterpolator(interpolator);
 119           transform=pointListTransform;
 120
 121           break;
 122         }
 123
 124         // ==========================
 125         // 4D vector field
 126         // ==========================
 127       case 1:
 128         {
 129           // Deformation field transform
 130           typedef clitk::DeformationFieldTransform<double, ImageDimension,ImageDimension, SpaceDimension> DeformationFieldTransformType;
 131           DeformationFieldTransformType::Pointer deformationFieldTransform=DeformationFieldTransformType::New();
 132
 133           // The deformation field
 134           typedef DeformationFieldTransformType::DeformationFieldType DeformationFieldType;
 135           typedef itk::ImageFileReader<DeformationFieldType> InputReaderType;
 136           InputReaderType::Pointer reader = InputReaderType::New();
 137           reader->SetFileName( m_ArgsInfo.input_arg);
 138           reader->Update();
 139           DeformationFieldType::Pointer input= reader->GetOutput();
 140           deformationFieldTransform->SetDeformationField(input);
 141
 142           // Vector Interpolator
 143           typedef clitk::GenericVectorInterpolator<args_info_clitkPointTrajectory,DeformationFieldType, double> GenericVectorInterpolatorType;
 144           GenericVectorInterpolatorType::Pointer genericInterpolator=GenericVectorInterpolatorType::New();
 145           genericInterpolator->SetArgsInfo(m_ArgsInfo);
 146           typedef itk::VectorInterpolateImageFunction<DeformationFieldType, double> InterpolatorType;
 147           InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
 148           deformationFieldTransform->SetInterpolator(interpolator);
 149           transform=deformationFieldTransform;
 150
 151           break;
 152         }
 153
 154         // ==========================
 155         // Spatio-Temporal transform
 156         // ==========================
 157       case 2:
 158         {
 159           // S-T transform
 160           typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform< double, ImageDimension, ImageDimension > TransformType;
 161           TransformType::Pointer spatioTemporalTransform = TransformType::New();
 162
 163
 164           // Spline orders:  Default is cubic splines
 165           CoefficientImageType::RegionType::SizeType splineOrders ;
 166           splineOrders.Fill(3);
 167           if (m_ArgsInfo.order_given)
 168             for(unsigned int i=0; i<ImageDimension;i++)
 169               splineOrders[i]=m_ArgsInfo.order_arg[i];
 170           if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 171
 172           // Coefficient image
 173           typedef itk::ImageFileReader<CoefficientImageType> InputReaderType;
 174           InputReaderType::Pointer reader = InputReaderType::New();
 175           reader->SetFileName( m_ArgsInfo.input_arg);
 176           reader->Update();
 177           CoefficientImageType::Pointer input= reader->GetOutput();
 178           //      itk::Vector<double,3> vector;
 179           //      vector.Fill(0.);
 180           //      vector[2]=100;
 181           //      input->FillBuffer(vector);
 182
 183           // Mask
 184           typedef itk::ImageMaskSpatialObject<  ImageDimension >   MaskType;
 185           MaskType::Pointer  spatialObjectMask=NULL;
 186           if (m_ArgsInfo.mask_given)
 187             {
 188               typedef itk::Image< unsigned char, ImageDimension >   ImageMaskType;
 189               typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 190               MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 191               maskReader->SetFileName(m_ArgsInfo.mask_arg);
 192
 193               try
 194                 {
 195                   maskReader->Update();
 196                 }
 197               catch( itk::ExceptionObject & err )
 198                 {
 199                   std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 200                   std::cerr << err << std::endl;
 201                   return;
 202                 }
 203               if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 204
 205               // Set the image to the spatialObject
 206               spatialObjectMask = MaskType::New();
 207               spatialObjectMask->SetImage( maskReader->GetOutput() );
 208             }
 209
 210           // Samplingfactors
 211           CoefficientImageType::SizeType samplingFactors;
 212           for (unsigned int i=0; i< ImageDimension-1; i++)
 213             {
 214               samplingFactors[i]= (int) ( input->GetSpacing()[i]/ m_ArgsInfo.spacing_arg);
 215               if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 216             }
 217           samplingFactors[ImageDimension-1]= (int) ( input->GetSpacing()[ImageDimension-1]/ m_ArgsInfo.phaseIncrement_arg);
 218           if (m_Verbose) std::cout<<"Setting sampling factor "<<ImageDimension-1<<" to "<<samplingFactors[ImageDimension-1]<<"..."<<std::endl;
 219           if( (m_ArgsInfo.shape_arg==3) |
 220               (m_ArgsInfo.shape_arg==4) |
 221               (m_ArgsInfo.shape_arg==6) |
 222               (m_ArgsInfo.shape_arg==8)
 223               ) samplingFactors[ImageDimension-1]*=2.5;
 224
 225           // Set
 226           spatioTemporalTransform->SetTransformShape(m_ArgsInfo.shape_arg);
 227           spatioTemporalTransform->SetSplineOrders(splineOrders);
 228           spatioTemporalTransform->SetMask(spatialObjectMask);
 229           spatioTemporalTransform->SetLUTSamplingFactors(samplingFactors);
 230           spatioTemporalTransform->SetCoefficientImage(input);
 231           transform=spatioTemporalTransform;
 232
 233           break;
 234         }
 235
 236 //      // ==========================
 237 //      // Spatio-Temporal transform
 238 //      // ==========================
 239 //       case 3:
 240 //      {
 241 //        // S-T transform
 242 //        typedef clitk::BSplineSpatioTemporalDeformableTransform< double, ImageDimension, ImageDimension > TransformType;
 243 //        TransformType::Pointer spatioTemporalTransform = TransformType::New();
 244
 245
 246 //        // Spline orders:  Default is cubic splines
 247 //        CoefficientImageType::RegionType::SizeType splineOrders ;
 248 //        splineOrders.Fill(3);
 249 //        if (m_ArgsInfo.order_given)
 250 //          for(unsigned int i=0; i<ImageDimension;i++)
 251 //            splineOrders[i]=m_ArgsInfo.order_arg[i];
 252 //        if (m_Verbose) std::cout<<"Setting the spline orders  to "<<splineOrders<<"..."<<std::endl;
 253
 254 //        // Coefficient image
 255 //        typedef itk::ImageFileReader<CoefficientImageType> InputReaderType;
 256 //        InputReaderType::Pointer reader = InputReaderType::New();
 257 //        reader->SetFileName( m_ArgsInfo.input_arg);
 258 //        reader->Update();
 259 //        CoefficientImageType::Pointer input= reader->GetOutput();
 260 //        //      itk::Vector<double,3> vector;
 261 //        //      vector.Fill(0.);
 262 //        //      vector[2]=100;
 263 //        //      input->FillBuffer(vector);
 264
 265 //        // Mask
 266 //        typedef itk::ImageMaskSpatialObject<  ImageDimension >   MaskType;
 267 //        MaskType::Pointer  spatialObjectMask=NULL;
 268 //        if (m_ArgsInfo.mask_given)
 269 //          {
 270 //            typedef itk::Image< unsigned char, ImageDimension >   ImageMaskType;
 271 //            typedef itk::ImageFileReader< ImageMaskType >    MaskReaderType;
 272 //            MaskReaderType::Pointer  maskReader = MaskReaderType::New();
 273 //            maskReader->SetFileName(m_ArgsInfo.mask_arg);
 274
 275 //            try
 276 //              {
 277 //                maskReader->Update();
 278 //              }
 279 //            catch( itk::ExceptionObject & err )
 280 //              {
 281 //                std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
 282 //                std::cerr << err << std::endl;
 283 //                return;
 284 //              }
 285 //            if (m_Verbose)std::cout <<"Mask was read..." <<std::endl;
 286
 287 //            // Set the image to the spatialObject
 288 //            spatialObjectMask = MaskType::New();
 289 //            spatialObjectMask->SetImage( maskReader->GetOutput() );
 290 //          }
 291
 292 //        // Samplingfactors
 293 //        CoefficientImageType::SizeType samplingFactors;
 294 //        for (unsigned int i=0; i< ImageDimension-1; i++)
 295 //          {
 296 //            samplingFactors[i]= (int) ( input->GetSpacing()[i]/ m_ArgsInfo.spacing_arg);
 297 //            if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
 298 //          }
 299 //        samplingFactors[ImageDimension-1]= (int) ( input->GetSpacing()[ImageDimension-1]/ m_ArgsInfo.phaseIncrement_arg);
 300 //        if (m_Verbose) std::cout<<"Setting sampling factor "<<ImageDimension-1<<" to "<<samplingFactors[ImageDimension-1]<<"..."<<std::endl;
 301
 302 //        // Set
 303 //        //spatioTemporalTransform->SetTransformShape(m_ArgsInfo.shape_arg);
 304 //        spatioTemporalTransform->SetSplineOrders(splineOrders);
 305 //        spatioTemporalTransform->SetMask(spatialObjectMask);
 306 //        spatioTemporalTransform->SetLUTSamplingFactors(samplingFactors);
 307 //        spatioTemporalTransform->SetCoefficientImage(input);
 308 //        transform=spatioTemporalTransform;
 309
 310 //        break;
 311 //      }
 312       }
 313
 314
 315     // -----------------------------------------------
 316     // Construct Spatio-Temporal Point lists 4D
 317     // -----------------------------------------------
 318     typedef itk::Point<double, ImageDimension> SpaceTimePointType;
 319     typedef clitk::Lists<SpaceTimePointType> PointListsType;
 320     PointListsType pointLists(referencePointList.size());
 321     SpaceTimePointType spaceTimePoint;
 322     double phase;
 323     for (unsigned int i=0; i<referencePointList.size(); i++)
 324       {
 325         for (unsigned int d=0; d<SpaceDimension; d++)
 326           {
 327             spaceTimePoint[d]=referencePointList[i][d];
 328           }
 329         phase=0;
 330         while (phase<10.)
 331           {
 332             spaceTimePoint[ImageDimension-1]=phase;
 333             pointLists[i].push_back(spaceTimePoint);
 334             phase+=m_ArgsInfo.phaseIncrement_arg;
 335           }
 336       }
 337
 338
 339     // -----------------------------------------------
 340     // Transform Points
 341     // -----------------------------------------------
 342     typedef itk::Vector<double, ImageDimension> VectorType;
 343     typedef clitk::List<VectorType> VectorListType;
 344     typedef clitk::Lists<VectorType> VectorListsType;
 345     VectorListsType displacementLists(pointLists.size());
 346     VectorType displacement;
 347     for (unsigned int i=0; i<pointLists.size(); i++)
 348       {
 349         if (m_Verbose) std::cout<<"Transforming point "<<pointLists[i][0]<<"..."<<std::endl;
 350         for (unsigned int j=0; j<pointLists[i].size(); j++)
 351           {
 352             spaceTimePoint= transform->TransformPoint(pointLists[i][j]);
 353             if (m_Verbose) std::cout<<"Transformed point "<<spaceTimePoint<<"..."<<std::endl;
 354             displacement=spaceTimePoint-pointLists[i][j];
 355             displacementLists[i].push_back(displacement);
 356           }
 357       }
 358
 359
 360     // -----------------------------------------------
 361     // Write displacements
 362     // -----------------------------------------------
 363     std::vector<std::string> filenames;
 364     for (unsigned int i=0;i<displacementLists.size();i++)
 365       {
 366         std::ostringstream number_ostr;
 367         number_ostr << i;
 368         std::string number_str;
 369         if (i<10)  number_str=  "0"+number_ostr.str();
 370         else number_str =  number_ostr.str();
 371         filenames.push_back(m_ArgsInfo.trajectory_arg+number_str);
 372       }
 373     displacementLists.Write(filenames, m_Verbose );
 374
 375   }
 376
 377
 378 } //end clitk
 379
 380 #endif  //#define clitkPointTrajectoryGenericFilter_cxx