registration/clitkPointListTransform.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://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 __clitkPointListTransform_txx
  19 #define __clitkPointListTransform_txx
  20 #include "clitkPointListTransform.h"
  21
  22
  23 namespace clitk
  24 {
  25
  26   // Constructor
  27   template<class TScalarType, unsigned int NDimensions, unsigned int NOutputDimensions>
  28   PointListTransform<TScalarType, NDimensions, NOutputDimensions>
  29   ::PointListTransform():Superclass(NOutputDimensions,1)
  30   {
  31     m_PointLists.resize(0);
  32     m_PointList.resize(1);
  33   }
  34
  35   // Find the point list in the lists
  36   template<class TScalarType, unsigned int NDimensions,unsigned int NOutputDimensions>
  37   typename PointListTransform<TScalarType,  NDimensions, NOutputDimensions>::PointListType
  38   PointListTransform<TScalarType, NDimensions,NOutputDimensions>::
  39   GetCorrespondingPointList(const InputPointType &inputPoint) const
  40   {
  41     SpacePointType point;
  42     for(unsigned int j=0; j< SpaceDimension;j++)
  43       point[j]=inputPoint[j];
  44
  45     if(m_PointList[0]==point) return m_PointList;
  46     else
  47       {
  48         for (unsigned int i=0; i< m_PointLists.size();i++)
  49           if(m_PointLists[i][0]==point) return m_PointLists[i];
  50       }
  51
  52     itkExceptionMacro(<<"Point List not found");
  53   }
  54
  55
  56   // Transform a point
  57   template<class TScalarType, unsigned int NDimensions,unsigned int NOutputDimensions>
  58   typename PointListTransform<TScalarType,  NDimensions,NOutputDimensions>::OutputPointType
  59   PointListTransform<TScalarType, NDimensions,NOutputDimensions>::
  60   TransformPoint(const InputPointType &inputPoint) const
  61   {
  62
  63     // -------------------------------
  64     // Get the corresponding point list
  65     m_PointList = this->GetCorrespondingPointList(inputPoint);
  66
  67     // -------------------------------
  68     // Create 1D vector image
  69     typename PointListImageType::Pointer pointListImage=PointListImageType::New();
  70     typename PointListImageType::RegionType region;
  71     region.SetSize(0,m_PointList.size()+6);
  72     pointListImage->SetRegions(region);
  73     pointListImage->Allocate();
  74     typename PointListImageType::SpacingType spacing;
  75     spacing[0]=1;
  76     pointListImage->SetSpacing(spacing);
  77     typename PointListImageType::PointType origin;
  78     origin[0]=-2.;
  79     pointListImage->SetOrigin(origin);
  80
  81
  82     // -------------------------------
  83     // Copy Point list to image
  84     typedef itk::ImageRegionIterator<PointListImageType> IteratorType;
  85     IteratorType it(pointListImage, pointListImage->GetLargestPossibleRegion());
  86
  87     // First points are the last
  88     PointListImagePixelType pixel;
  89     for (unsigned int j=0; j<2;j++)
  90       {
  91         for (unsigned int i=0; i <SpaceDimension; i++)
  92           pixel[i]=m_PointList[m_PointList.size()-2+j][i];
  93         it.Set(pixel);
  94         ++it;
  95       }
  96
  97     // Copy the rest
  98     unsigned int position=0;
  99     while(position< m_PointList.size())
 100       {
 101         for (unsigned int i=0; i <SpaceDimension; i++)
 102           pixel[i]=m_PointList[position][i];
 103         it.Set(pixel);
 104         ++it;
 105         ++position;
 106       }
 107
 108     // last points are the first
 109     for (unsigned int j=0; j<4;j++)
 110       {
 111         for (unsigned int i=0; i <SpaceDimension; i++)
 112           pixel[i]=m_PointList[j][i];
 113         it.Set(pixel);
 114         ++it;
 115       }
 116
 117     // -------------------------------
 118     // Set 1D image to vectorInterpolator
 119     m_Interpolator->SetInputImage(pointListImage);
 120
 121
 122     // -------------------------------
 123     // Evaluate at phase value
 124     typename PointListImageType::PointType t;
 125     t[0]=inputPoint[ImageDimension-1];
 126
 127     // Inside valid region?
 128     if ( (t[0] >= 0) &&
 129          (t[0] < m_PointList.size()) )
 130       {
 131         pixel= m_Interpolator->Evaluate(t);
 132         OutputPointType outputPoint;
 133         for (unsigned int i=0; i < SpaceDimension; i++)
 134           outputPoint[i]=pixel[i];
 135         outputPoint[ImageDimension-1]=t[0];
 136         return outputPoint;
 137       }
 138     // No displacement
 139     else return inputPoint;
 140
 141   }
 142
 143
 144 } // namespace clitk
 145
 146 #endif