DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkSlowPolyLineParametricPath.hxx
Go to the documentation of this file.
1 #ifndef __itkSlowPolyLineParametricPath_hxx
2 #define __itkSlowPolyLineParametricPath_hxx
3 
5 #include <math.h>
6 
7 
8 
9 namespace itk
10 {
11 
12 //template<unsigned int VDimension>
13 //typename SlowPolyLineParametricPath<VDimension>::VectorType
14 //SlowPolyLineParametricPath<VDimension>
15 //::EvaluateDerivative(const InputType & input) const
16 //{
17 //}
18 
19 
20 
24 template <unsigned int VDimension>
27 {
28  this->SetDefaultInputStepSize( 0.3 );
29  m_UseCentralDifference=true;
30 }
31 
32 template< unsigned int VDimension >
33 typename PolyLineParametricPath< VDimension >::OutputType
35 ::Evaluate(const InputType & input) const
36 {
37  // Handle the endpoint carefully, since there is no following vertex
38  const VertexListType* vertexList = this->GetVertexList();
39  const InputType endPoint = static_cast< InputType >( vertexList->Size() - 1 );
40  if ( input > endPoint || itk::Math::FloatAlmostEqual( input, endPoint ) )
41  {
42  return vertexList->ElementAt(vertexList->Size() - 1); // the last vertex
43  }
44  if (input<0 || itk::Math::FloatAlmostEqual( input, 0.0 ))
45  {
46  return vertexList->ElementAt(0); // the first vertex
47  }
48 
49  const VertexType vertex0 = vertexList->ElementAt( int(input)<0.0 ? 0 : int(input) );
50  const VertexType vertex1 = vertexList->ElementAt( int(input)+1>=endPoint ? endPoint-1 : int(input)+1);
51 
52  const double fractionOfLineSegment = input - int(input);
53 
54  const PointType outputPoint = vertex0 + ( vertex1 - vertex0 ) * fractionOfLineSegment;
55 
56  // For some stupid reason, there is no easy way to cast from a point to a
57  // continuous index.
58  OutputType output;
59  for ( unsigned int i = 0; i < VDimension; i++ )
60  {
61  output[i] = outputPoint[i];
62  }
63 
64  return output;
65 }
66 
67 template<unsigned int VDimension>
68 typename PolyLineParametricPath<VDimension>::VectorType
70 ::EvaluateDerivative(const InputType & input, bool isDerivativeNormalizedByDistance) const
71 {
72  if (m_UseCentralDifference)
73  {
74  //Get next integral time-point
75  const InputType nextTimepoint = std::floor(input + 1.0);
76 
77  //Get previous integral time-point
78  const InputType previousTimepoint = std::floor(input - 1.0);
79 
80  //Calculate the continuous index for both points
81  const ContinuousIndexType nextIndex = this->Evaluate(nextTimepoint);
82  const ContinuousIndexType prevIndex = this->Evaluate(previousTimepoint);
83 
84  //For some reason, there's no way to convert ContinuousIndexType to VectorType
85  VectorType partialDerivatives;
86  for (unsigned int i = 0; i < VDimension; ++i)
87  {
88  partialDerivatives[i] = (nextIndex[i] - prevIndex[i]);
89  }
90 
91  if (isDerivativeNormalizedByDistance)
92  {
93  double dist = nextIndex.EuclideanDistanceTo(prevIndex);
94  if (dist>1e-10)
95  partialDerivatives /= dist;
96  }
97 
98  return partialDerivatives;
99  }
100  else
101  return Superclass::EvaluateDerivative(input);
102 }
103 
104 template<unsigned int VDimension>
108 {
109  int iterationCount;
110  bool tooSmall;
111  bool tooBig;
112  InputType inputStepSize;
113  InputType finalInputValue;
114  OffsetType offset;
115  IndexType currentImageIndex;
116  IndexType nextImageIndex;
117  IndexType finalImageIndex;
118 
119  iterationCount = 0;
120  inputStepSize = this->GetDefaultInputStepSize();
121 
122  // Are we already at (or past) the end of the input?
123  finalInputValue = this->EndOfInput();
124  currentImageIndex = this->EvaluateToIndex( input );
125  finalImageIndex = this->EvaluateToIndex( finalInputValue );
126  offset = finalImageIndex - currentImageIndex;
127  if( ( offset == this->GetZeroOffset() && input != this->StartOfInput() ) ||
128  ( input >=finalInputValue ) )
129  {
130  return this->GetZeroOffset();
131  }
132 
133  do
134  {
135  if( iterationCount++ > 10000 ) {itkExceptionMacro(<<"Too many iterations");}
136 
137  nextImageIndex = this->EvaluateToIndex( input + inputStepSize );
138  offset = nextImageIndex - currentImageIndex;
139 
140  tooBig = false;
141  tooSmall = ( offset == this->GetZeroOffset() );
142  if( tooSmall )
143  {
144  // increase the input step size, but don't go past the end of the input
145  inputStepSize *= 2;
146  if( (input + inputStepSize) >= finalInputValue ){
147  //inputStepSize = finalInputValue - input;
148  inputStepSize += this->GetDefaultInputStepSize();
149  }
150  }
151  else
152  {
153  // Search for an offset dimension that is too big
154  for( unsigned int i=0; i<VDimension && !tooBig; i++ )
155  {
156  tooBig = ( offset[i] >= 2 || offset[i] <= -2 );
157  }
158 
159  if( tooBig ){
160  //inputStepSize /= 1.5;
161  inputStepSize -= (this->GetDefaultInputStepSize()/0.5);
162  }
163  }
164  }
165  while( tooSmall || tooBig );
166 
167  input += inputStepSize;
168  return offset;
169 }
170 
171 
175 template <unsigned int VDimension>
176 void
178 ::PrintSelf( std::ostream& os, Indent indent) const
179 {
180  Superclass::PrintSelf( os, indent );
181 }
182 
183 
184 
185 } // end namespaceitk
186 
187 #endif
void PrintSelf(std::ostream &os, Indent indent) const ITK_OVERRIDE
virtual OffsetType IncrementInput(InputType &input) const ITK_OVERRIDE
virtual OutputType Evaluate(const InputType &input) const ITK_OVERRIDE
Superclass::ContinuousIndexType ContinuousIndexType
VectorType EvaluateDerivative(const InputType &input, bool isDerivativeNormalizedByDistance=false) const