ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkBSplineInterpolateImageFunction.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 /*=========================================================================
00019  *
00020  *  Portions of this file are subject to the VTK Toolkit Version 3 copyright.
00021  *
00022  *  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
00023  *
00024  *  For complete copyright, license and disclaimer of warranty information
00025  *  please refer to the NOTICE file at the top of the ITK source tree.
00026  *
00027  *=========================================================================*/
00028 #ifndef __itkBSplineInterpolateImageFunction_h
00029 #define __itkBSplineInterpolateImageFunction_h
00030 
00031 #include <vector>
00032 
00033 #include "itkInterpolateImageFunction.h"
00034 #include "vnl/vnl_matrix.h"
00035 
00036 #include "itkBSplineDecompositionImageFilter.h"
00037 #include "itkConceptChecking.h"
00038 #include "itkCovariantVector.h"
00039 
00040 namespace itk
00041 {
00080 template<
00081   class TImageType,
00082   class TCoordRep = double,
00083   class TCoefficientType = double >
00084 class ITK_EXPORT BSplineInterpolateImageFunction:
00085   public InterpolateImageFunction< TImageType, TCoordRep >
00086 {
00087 public:
00089   typedef BSplineInterpolateImageFunction                   Self;
00090   typedef InterpolateImageFunction< TImageType, TCoordRep > Superclass;
00091   typedef SmartPointer< Self >                              Pointer;
00092   typedef SmartPointer< const Self >                        ConstPointer;
00093 
00095   itkTypeMacro(BSplineInterpolateImageFunction, InterpolateImageFunction);
00096 
00098   itkNewMacro(Self);
00099 
00101   typedef typename Superclass::OutputType OutputType;
00102 
00104   typedef typename Superclass::InputImageType InputImageType;
00105 
00107   itkStaticConstMacro(ImageDimension, unsigned int, Superclass::ImageDimension);
00108 
00110   typedef typename Superclass::IndexType IndexType;
00111 
00113   typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
00114 
00116   typedef typename Superclass::PointType PointType;
00117 
00119   typedef ImageLinearIteratorWithIndex< TImageType > Iterator;
00120 
00122   typedef TCoefficientType CoefficientDataType;
00123   typedef Image< CoefficientDataType,
00124                  itkGetStaticConstMacro(ImageDimension) >
00125   CoefficientImageType;
00126 
00128   typedef BSplineDecompositionImageFilter< TImageType, CoefficientImageType > CoefficientFilter;
00129   typedef typename CoefficientFilter::Pointer                                 CoefficientFilterPointer;
00130 
00132   typedef CovariantVector< OutputType,
00133                            itkGetStaticConstMacro(ImageDimension) >
00134   CovariantVectorType;
00135 
00144   virtual OutputType Evaluate(const PointType & point) const
00145   {
00146     ContinuousIndexType index;
00147 
00148     this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
00149                                                                    index);
00150     // No thread info passed in, so call method that doesn't need thread ID.
00151     return ( this->EvaluateAtContinuousIndex(index) );
00152   }
00153 
00154   virtual OutputType Evaluate(const PointType & point,
00155                               ThreadIdType threadID) const
00156   {
00157     ContinuousIndexType index;
00158 
00159     this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
00160                                                                    index);
00161     return ( this->EvaluateAtContinuousIndex(index, threadID) );
00162   }
00163 
00164   virtual OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &
00165                                                index) const
00166   {
00167     // Don't know thread information, make evaluateIndex, weights on the stack.
00168     // Slower, but safer.
00169     vnl_matrix< long >   evaluateIndex( ImageDimension, ( m_SplineOrder + 1 ) );
00170     vnl_matrix< double > weights( ImageDimension, ( m_SplineOrder + 1 ) );
00171 
00172     // Pass evaluateIndex, weights by reference. They're only good as long
00173     // as this method is in scope.
00174     return this->EvaluateAtContinuousIndexInternal(index,
00175                                                    evaluateIndex,
00176                                                    weights);
00177   }
00178 
00179   virtual OutputType EvaluateAtContinuousIndex(const ContinuousIndexType &
00180                                                index,
00181                                                ThreadIdType threadID) const;
00182 
00183   CovariantVectorType EvaluateDerivative(const PointType & point) const
00184   {
00185     ContinuousIndexType index;
00186 
00187     this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
00188                                                                    index);
00189     // No thread info passed in, so call method that doesn't need thread ID.
00190     return ( this->EvaluateDerivativeAtContinuousIndex(index) );
00191   }
00192 
00193   CovariantVectorType EvaluateDerivative(const PointType & point,
00194                                          ThreadIdType threadID) const
00195   {
00196     ContinuousIndexType index;
00197 
00198     this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
00199                                                                    index);
00200     return ( this->EvaluateDerivativeAtContinuousIndex(index, threadID) );
00201   }
00202 
00203   CovariantVectorType EvaluateDerivativeAtContinuousIndex(
00204     const ContinuousIndexType & x) const
00205   {
00206     // Don't know thread information, make evaluateIndex, weights,
00207     // weightsDerivative
00208     // on the stack.
00209     // Slower, but safer.
00210     vnl_matrix< long >   evaluateIndex( ImageDimension, ( m_SplineOrder + 1 ) );
00211     vnl_matrix< double > weights( ImageDimension, ( m_SplineOrder + 1 ) );
00212     vnl_matrix< double > weightsDerivative( ImageDimension, ( m_SplineOrder + 1 ) );
00213 
00214     // Pass evaluateIndex, weights, weightsDerivative by reference. They're only
00215     // good
00216     // as long as this method is in scope.
00217     return this->EvaluateDerivativeAtContinuousIndexInternal(x,
00218                                                              evaluateIndex,
00219                                                              weights,
00220                                                              weightsDerivative);
00221   }
00222 
00223   CovariantVectorType EvaluateDerivativeAtContinuousIndex(
00224     const ContinuousIndexType & x,
00225     ThreadIdType threadID) const;
00226 
00227   void EvaluateValueAndDerivative(const PointType & point,
00228                                   OutputType & value,
00229                                   CovariantVectorType & deriv) const
00230   {
00231     ContinuousIndexType index;
00232 
00233     this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
00234                                                                    index);
00235 
00236     // No thread info passed in, so call method that doesn't need thread ID.
00237     this->EvaluateValueAndDerivativeAtContinuousIndex(index,
00238                                                       value,
00239                                                       deriv);
00240   }
00241 
00242   void EvaluateValueAndDerivative(const PointType & point,
00243                                   OutputType & value,
00244                                   CovariantVectorType & deriv,
00245                                   ThreadIdType threadID) const
00246   {
00247     ContinuousIndexType index;
00248 
00249     this->GetInputImage()->TransformPhysicalPointToContinuousIndex(point,
00250                                                                    index);
00251     this->EvaluateValueAndDerivativeAtContinuousIndex(index,
00252                                                       value,
00253                                                       deriv,
00254                                                       threadID);
00255   }
00256 
00257   void EvaluateValueAndDerivativeAtContinuousIndex(
00258     const ContinuousIndexType & x,
00259     OutputType & value,
00260     CovariantVectorType & deriv
00261     ) const
00262   {
00263     // Don't know thread information, make evaluateIndex, weights,
00264     // weightsDerivative
00265     // on the stack.
00266     // Slower, but safer.
00267     vnl_matrix< long >   evaluateIndex( ImageDimension, ( m_SplineOrder + 1 ) );
00268     vnl_matrix< double > weights( ImageDimension, ( m_SplineOrder + 1 ) );
00269     vnl_matrix< double > weightsDerivative( ImageDimension, ( m_SplineOrder + 1 ) );
00270 
00271     // Pass evaluateIndex, weights, weightsDerivative by reference. They're only
00272     // good
00273     // as long as this method is in scope.
00274     this->EvaluateValueAndDerivativeAtContinuousIndexInternal(x,
00275                                                               value,
00276                                                               deriv,
00277                                                               evaluateIndex,
00278                                                               weights,
00279                                                               weightsDerivative);
00280   }
00281 
00282   void EvaluateValueAndDerivativeAtContinuousIndex(
00283     const ContinuousIndexType & x,
00284     OutputType & value,
00285     CovariantVectorType & deriv,
00286     ThreadIdType threadID) const;
00287 
00290   void SetSplineOrder(unsigned int SplineOrder);
00291 
00292   itkGetConstMacro(SplineOrder, int);
00293 
00294   void SetNumberOfThreads(ThreadIdType numThreads);
00295 
00296   itkGetConstMacro(NumberOfThreads, ThreadIdType);
00297 
00299   virtual void SetInputImage(const TImageType *inputData);
00300 
00311   itkSetMacro(UseImageDirection, bool);
00312   itkGetConstMacro(UseImageDirection, bool);
00313   itkBooleanMacro(UseImageDirection);
00314 protected:
00316 
00335   virtual OutputType EvaluateAtContinuousIndexInternal(const ContinuousIndexType & index,
00336                                                        vnl_matrix< long > & evaluateIndex,
00337                                                        vnl_matrix< double > & weights) const;
00338 
00339   virtual void EvaluateValueAndDerivativeAtContinuousIndexInternal(const ContinuousIndexType & x,
00340                                                                    OutputType & value,
00341                                                                    CovariantVectorType & derivativeValue,
00342                                                                    vnl_matrix< long > & evaluateIndex,
00343                                                                    vnl_matrix< double > & weights,
00344                                                                    vnl_matrix< double > & weightsDerivative
00345                                                                    ) const;
00346 
00347   virtual CovariantVectorType EvaluateDerivativeAtContinuousIndexInternal(const ContinuousIndexType & x,
00348                                                                           vnl_matrix< long > & evaluateIndex,
00349                                                                           vnl_matrix< double > & weights,
00350                                                                           vnl_matrix< double > & weightsDerivative
00351                                                                           ) const;
00352 
00353   BSplineInterpolateImageFunction();
00354   ~BSplineInterpolateImageFunction();
00355   void PrintSelf(std::ostream & os, Indent indent) const;
00356 
00357   // These are needed by the smoothing spline routine.
00358   // temp storage for processing of Coefficients
00359   std::vector< CoefficientDataType >    m_Scratch;
00360   // Image size
00361   typename TImageType::SizeType m_DataLength;
00362   // User specified spline order (3rd or cubic is the default)
00363   unsigned int m_SplineOrder;
00364 
00365   // Spline coefficients
00366   typename CoefficientImageType::ConstPointer m_Coefficients;
00367 private:
00368   BSplineInterpolateImageFunction(const Self &); //purposely not implemented
00369   void operator=(const Self &);                  //purposely not implemented
00370 
00372   void SetInterpolationWeights(const ContinuousIndexType & x,
00373                                const vnl_matrix< long > & EvaluateIndex,
00374                                vnl_matrix< double > & weights,
00375                                unsigned int splineOrder) const;
00376 
00378   void SetDerivativeWeights(const ContinuousIndexType & x,
00379                             const vnl_matrix< long > & EvaluateIndex,
00380                             vnl_matrix< double > & weights,
00381                             unsigned int splineOrder) const;
00382 
00385   void GeneratePointsToIndex();
00386 
00388   void DetermineRegionOfSupport(vnl_matrix< long > & evaluateIndex,
00389                                 const ContinuousIndexType & x,
00390                                 unsigned int splineOrder) const;
00391 
00394   void ApplyMirrorBoundaryConditions(vnl_matrix< long > & evaluateIndex,
00395                                      unsigned int splineOrder) const;
00396 
00397   Iterator m_CIterator;                                    // Iterator for
00398                                                            // traversing spline
00399                                                            // coefficients.
00400   unsigned long m_MaxNumberInterpolationPoints;            // number of
00401                                                            // neighborhood
00402                                                            // points used for
00403                                                            // interpolation
00404   std::vector< IndexType > m_PointsToIndex;                // Preallocation of
00405                                                            // interpolation
00406                                                            // neighborhood
00407                                                            // indicies
00408 
00409   CoefficientFilterPointer m_CoefficientFilter;
00410 
00411   // flag to take or not the image direction into account when computing the
00412   // derivatives.
00413   bool m_UseImageDirection;
00414 
00415   ThreadIdType          m_NumberOfThreads;
00416   vnl_matrix< long > *  m_ThreadedEvaluateIndex;
00417   vnl_matrix< double > *m_ThreadedWeights;
00418   vnl_matrix< double > *m_ThreadedWeightsDerivative;
00419 };
00420 } // namespace itk
00421 
00422 #ifndef ITK_MANUAL_INSTANTIATION
00423 #include "itkBSplineInterpolateImageFunction.hxx"
00424 #endif
00425 
00426 #endif
00427