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