]> Creatis software - clitk.git/blob - registration/clitkPointListTransform.txx
Moved from repository clitk to clitk.private/tests_dav
[clitk.git] / 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