ANIMA  4.0
animaVectorModelLinearInterpolateImageFunction.hxx
Go to the documentation of this file.
1 #pragma once
2 
4 #include <vnl/vnl_math.h>
5 
6 namespace anima
7 {
8 
12 template<class TInputImage, class TCoordRep>
13 const unsigned long
14 VectorModelLinearInterpolateImageFunction< TInputImage, TCoordRep >
15 ::m_Neighbors = 1 << TInputImage::ImageDimension;
16 
17 
21 template<class TInputImage, class TCoordRep>
24 {
25 
26 }
27 
31 template<class TInputImage, class TCoordRep>
36 {
37  IndexType baseIndex;
38  double distance[ImageDimension], oppDistance[ImageDimension];
39 
40  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
41  {
42  baseIndex[dim] = itk::Math::Floor< IndexValueType >( index[dim] );
43  distance[dim] = index[dim] - static_cast< double >( baseIndex[dim] );
44  oppDistance[dim] = 1.0 - distance[dim];
45  }
46 
47  unsigned int vectorDim = this->GetInputImage()->GetNumberOfComponentsPerPixel();
48  OutputType output(vectorDim);
49  this->InitializeZeroPixel(output);
50 
51  double totalOverlap = 0;
52 
53  for( unsigned int counter = 0; counter < m_Neighbors; ++counter)
54  {
55  double overlap = 1.0; // fraction overlap
56  unsigned int upper = counter; // each bit indicates upper/lower neighbour
57  IndexType neighIndex;
58 
59  // get neighbor index and overlap fraction
60  bool okValue = true;
61  for (unsigned int dim = 0; dim < ImageDimension; ++dim)
62  {
63  if (upper & 1)
64  {
65  neighIndex[dim] = baseIndex[dim] + 1;
66 
67  if (neighIndex[dim] > this->m_EndIndex[dim])
68  {
69  okValue = false;
70  break;
71  }
72 
73  overlap *= distance[dim];
74  }
75  else
76  {
77  neighIndex[dim] = baseIndex[dim];
78 
79  if (neighIndex[dim] < this->m_StartIndex[dim])
80  {
81  okValue = false;
82  break;
83  }
84 
85  overlap *= oppDistance[dim];
86  }
87 
88  upper >>= 1;
89  }
90 
91  // get neighbor value only if overlap is not zero
92  if ((overlap > 0) && okValue)
93  {
94  VectorPixelType input = static_cast <VectorPixelType> (this->GetInputImage()->GetPixel(neighIndex));
95 
96  if (!isZero(input))
97  {
98  this->AddValue(input,overlap,output);
99  totalOverlap += overlap;
100  }
101  }
102 
103  }
104 
105  if (totalOverlap >= 0.5)
106  output /= totalOverlap;
107  else
108  this->InitializeZeroPixel(output);
109 
110  return output;
111 }
112 
113 } // end namespace itk
virtual OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &index) const ITK_OVERRIDE
itk::VariableLengthVector< typename TInputImage::IOPixelType > VectorPixelType