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