ITK  4.4.0
Insight Segmentation and Registration Toolkit
itkBSplineInterpolateImageFunction.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 /*=========================================================================
19  *
20  * Portions of this file are subject to the VTK Toolkit Version 3 copyright.
21  *
22  * Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
23  *
24  * For complete copyright, license and disclaimer of warranty information
25  * please refer to the NOTICE file at the top of the ITK source tree.
26  *
27  *=========================================================================*/
28 #ifndef __itkBSplineInterpolateImageFunction_h
29 #define __itkBSplineInterpolateImageFunction_h
30 
31 #include <vector>
32 
34 #include "vnl/vnl_matrix.h"
35 
37 #include "itkConceptChecking.h"
38 #include "itkCovariantVector.h"
39 
40 namespace itk
41 {
80 template<
81  class TImageType,
82  class TCoordRep = double,
83  class TCoefficientType = double >
85  public InterpolateImageFunction< TImageType, TCoordRep >
86 {
87 public:
93 
96 
98  itkNewMacro(Self);
99 
101  typedef typename Superclass::OutputType OutputType;
102 
104  typedef typename Superclass::InputImageType InputImageType;
105 
107  itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
108 
110  typedef typename Superclass::IndexType IndexType;
111 
113  typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
114 
116  typedef typename Superclass::PointType PointType;
117 
120 
122  typedef TCoefficientType CoefficientDataType;
123  typedef Image< CoefficientDataType,
124  itkGetStaticConstMacro(ImageDimension) >
126 
130 
132  typedef CovariantVector< OutputType,
133  itkGetStaticConstMacro(ImageDimension) >
135 
144  virtual OutputType Evaluate(const PointType & point) const
145  {
146  ContinuousIndexType index;
147 
148  this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
149  index);
150  // No thread info passed in, so call method that doesn't need thread ID.
151  return ( this->EvaluateAtContinuousIndex(index) );
152  }
153 
154  virtual OutputType Evaluate(const PointType & point,
155  ThreadIdType threadID) const
156  {
157  ContinuousIndexType index;
158 
159  this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
160  index);
161  return ( this->EvaluateAtContinuousIndex(index, threadID) );
162  }
163 
164  virtual OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &
165  index) const
166  {
167  // Don't know thread information, make evaluateIndex, weights on the stack.
168  // Slower, but safer.
169  vnl_matrix< long > evaluateIndex( ImageDimension, ( m_SplineOrder + 1 ) );
170  vnl_matrix< double > weights( ImageDimension, ( m_SplineOrder + 1 ) );
171 
172  // Pass evaluateIndex, weights by reference. They're only good as long
173  // as this method is in scope.
174  return this->EvaluateAtContinuousIndexInternal(index,
175  evaluateIndex,
176  weights);
177  }
178 
179  virtual OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &
180  index,
181  ThreadIdType threadID) const;
182 
183  CovariantVectorType EvaluateDerivative(const PointType & point) const
184  {
185  ContinuousIndexType index;
186 
187  this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
188  index);
189  // No thread info passed in, so call method that doesn't need thread ID.
190  return ( this->EvaluateDerivativeAtContinuousIndex(index) );
191  }
192 
193  CovariantVectorType EvaluateDerivative(const PointType & point,
194  ThreadIdType threadID) const
195  {
196  ContinuousIndexType index;
197 
198  this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
199  index);
200  return ( this->EvaluateDerivativeAtContinuousIndex(index, threadID) );
201  }
202 
203  CovariantVectorType EvaluateDerivativeAtContinuousIndex(
204  const ContinuousIndexType & x) const
205  {
206  // Don't know thread information, make evaluateIndex, weights,
207  // weightsDerivative
208  // on the stack.
209  // Slower, but safer.
210  vnl_matrix< long > evaluateIndex( ImageDimension, ( m_SplineOrder + 1 ) );
211  vnl_matrix< double > weights( ImageDimension, ( m_SplineOrder + 1 ) );
212  vnl_matrix< double > weightsDerivative( ImageDimension, ( m_SplineOrder + 1 ) );
213 
214  // Pass evaluateIndex, weights, weightsDerivative by reference. They're only
215  // good
216  // as long as this method is in scope.
217  return this->EvaluateDerivativeAtContinuousIndexInternal(x,
218  evaluateIndex,
219  weights,
220  weightsDerivative);
221  }
222 
223  CovariantVectorType EvaluateDerivativeAtContinuousIndex(
224  const ContinuousIndexType & x,
225  ThreadIdType threadID) const;
226 
227  void EvaluateValueAndDerivative(const PointType & point,
228  OutputType & value,
229  CovariantVectorType & deriv) const
230  {
231  ContinuousIndexType index;
232 
233  this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
234  index);
235 
236  // No thread info passed in, so call method that doesn't need thread ID.
237  this->EvaluateValueAndDerivativeAtContinuousIndex(index,
238  value,
239  deriv);
240  }
241 
242  void EvaluateValueAndDerivative(const PointType & point,
243  OutputType & value,
244  CovariantVectorType & deriv,
245  ThreadIdType threadID) const
246  {
247  ContinuousIndexType index;
248 
249  this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
250  index);
251  this->EvaluateValueAndDerivativeAtContinuousIndex(index,
252  value,
253  deriv,
254  threadID);
255  }
256 
257  void EvaluateValueAndDerivativeAtContinuousIndex(
258  const ContinuousIndexType & x,
259  OutputType & value,
260  CovariantVectorType & deriv
261  ) const
262  {
263  // Don't know thread information, make evaluateIndex, weights,
264  // weightsDerivative
265  // on the stack.
266  // Slower, but safer.
267  vnl_matrix< long > evaluateIndex( ImageDimension, ( m_SplineOrder + 1 ) );
268  vnl_matrix< double > weights( ImageDimension, ( m_SplineOrder + 1 ) );
269  vnl_matrix< double > weightsDerivative( ImageDimension, ( m_SplineOrder + 1 ) );
270 
271  // Pass evaluateIndex, weights, weightsDerivative by reference. They're only
272  // good
273  // as long as this method is in scope.
274  this->EvaluateValueAndDerivativeAtContinuousIndexInternal(x,
275  value,
276  deriv,
277  evaluateIndex,
278  weights,
279  weightsDerivative);
280  }
281 
282  void EvaluateValueAndDerivativeAtContinuousIndex(
283  const ContinuousIndexType & x,
284  OutputType & value,
285  CovariantVectorType & deriv,
286  ThreadIdType threadID) const;
287 
290  void SetSplineOrder(unsigned int SplineOrder);
291 
292  itkGetConstMacro(SplineOrder, int);
293 
294  void SetNumberOfThreads(ThreadIdType numThreads);
295 
296  itkGetConstMacro(NumberOfThreads, ThreadIdType);
297 
299  virtual void SetInputImage(const TImageType *inputData);
300 
311  itkSetMacro(UseImageDirection, bool);
312  itkGetConstMacro(UseImageDirection, bool);
313  itkBooleanMacro(UseImageDirection);
315 
316 protected:
317 
336  virtual OutputType EvaluateAtContinuousIndexInternal(const ContinuousIndexType & index,
337  vnl_matrix< long > & evaluateIndex,
338  vnl_matrix< double > & weights) const;
339 
340  virtual void EvaluateValueAndDerivativeAtContinuousIndexInternal(const ContinuousIndexType & x,
341  OutputType & value,
342  CovariantVectorType & derivativeValue,
343  vnl_matrix< long > & evaluateIndex,
344  vnl_matrix< double > & weights,
345  vnl_matrix< double > & weightsDerivative
346  ) const;
347 
348  virtual CovariantVectorType EvaluateDerivativeAtContinuousIndexInternal(const ContinuousIndexType & x,
349  vnl_matrix< long > & evaluateIndex,
350  vnl_matrix< double > & weights,
351  vnl_matrix< double > & weightsDerivative
352  ) const;
353 
356  void PrintSelf(std::ostream & os, Indent indent) const;
357 
358  // These are needed by the smoothing spline routine.
359  // temp storage for processing of Coefficients
360  std::vector< CoefficientDataType > m_Scratch;
361  // Image size
362  typename TImageType::SizeType m_DataLength;
363  // User specified spline order (3rd or cubic is the default)
364  unsigned int m_SplineOrder;
365 
366  // Spline coefficients
368 
369 private:
370  BSplineInterpolateImageFunction(const Self &); //purposely not implemented
371  void operator=(const Self &); //purposely not implemented
372 
374  void SetInterpolationWeights(const ContinuousIndexType & x,
375  const vnl_matrix< long > & EvaluateIndex,
376  vnl_matrix< double > & weights,
377  unsigned int splineOrder) const;
378 
380  void SetDerivativeWeights(const ContinuousIndexType & x,
381  const vnl_matrix< long > & EvaluateIndex,
382  vnl_matrix< double > & weights,
383  unsigned int splineOrder) const;
384 
387  void GeneratePointsToIndex();
388 
390  void DetermineRegionOfSupport(vnl_matrix< long > & evaluateIndex,
391  const ContinuousIndexType & x,
392  unsigned int splineOrder) const;
393 
396  void ApplyMirrorBoundaryConditions(vnl_matrix< long > & evaluateIndex,
397  unsigned int splineOrder) const;
398 
399  Iterator m_CIterator; // Iterator for
400  // traversing spline
401  // coefficients.
402  unsigned long m_MaxNumberInterpolationPoints; // number of
403  // neighborhood
404  // points used for
405  // interpolation
406  std::vector< IndexType > m_PointsToIndex; // Preallocation of
407  // interpolation
408  // neighborhood
409  // indices
410 
412 
413  // flag to take or not the image direction into account when computing the
414  // derivatives.
416 
418  vnl_matrix< long > * m_ThreadedEvaluateIndex;
419  vnl_matrix< double > *m_ThreadedWeights;
420  vnl_matrix< double > *m_ThreadedWeightsDerivative;
421 };
422 } // namespace itk
423 
424 #ifndef ITK_MANUAL_INSTANTIATION
425 #include "itkBSplineInterpolateImageFunction.hxx"
426 #endif
427 
428 #endif
429