diff --git a/Code/BasicFilters/otbBSplineDecompositionImageFilter.h b/Code/BasicFilters/otbBSplineDecompositionImageFilter.h
new file mode 100755
index 0000000000000000000000000000000000000000..030d76edb021158b99699b0b707fe3dffca199a3
--- /dev/null
+++ b/Code/BasicFilters/otbBSplineDecompositionImageFilter.h
@@ -0,0 +1,129 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbBSplineDecompositionImageFilter_h
+#define __otbBSplineDecompositionImageFilter_h
+
+#include <vector>
+
+#include "itkImageLinearIteratorWithIndex.h"
+#include "vnl/vnl_matrix.h"
+
+#include "itkImageToImageFilter.h"
+
+namespace otb
+{
+/** \class otbBSplineDecompositionImageFilter
+ * \brief This class is an evolution of the itk::BSplineDecompositionImageFilter to handle
+ * huge images with this interpolator. For more documentation, please refer to the original
+ * class.
+ * \sa itk::BSplineDecompositionImageFilter
+ * \ingroup ImageFilters
+ */
+template <class TInputImage, class TOutputImage>
+class ITK_EXPORT BSplineDecompositionImageFilter :
+ public itk::ImageToImageFilter<TInputImage,TOutputImage>
+{
+public:
+ /** Standard class typedefs. */
+ typedef BSplineDecompositionImageFilter Self;
+ typedef itk::ImageToImageFilter<TInputImage,TOutputImage> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(BSplineDecompositionImageFilter, ImageToImageFilter);
+
+ /** New macro for creation of through a Smart Pointer */
+ itkNewMacro( Self );
+
+ /** Inherit input and output image types from Superclass. */
+ typedef typename Superclass::InputImageType InputImageType;
+ typedef typename Superclass::InputImagePointer InputImagePointer;
+ typedef typename Superclass::InputImageConstPointer InputImageConstPointer;
+ typedef typename Superclass::OutputImagePointer OutputImagePointer;
+
+ /** Dimension underlying input image. */
+ itkStaticConstMacro(ImageDimension, unsigned int,TInputImage::ImageDimension);
+ itkStaticConstMacro(OutputImageDimension, unsigned int,
+ TOutputImage::ImageDimension);
+
+ /** Iterator typedef support */
+ typedef itk::ImageLinearIteratorWithIndex<TOutputImage> OutputLinearIterator;
+
+ /** Get/Sets the Spline Order, supports 0th - 5th order splines. The default
+ * is a 3rd order spline. */
+ void SetSplineOrder(unsigned int SplineOrder);
+ itkGetMacro(SplineOrder, int);
+
+protected:
+ BSplineDecompositionImageFilter();
+ virtual ~BSplineDecompositionImageFilter() {};
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+ void GenerateData( );
+
+ /** These are needed by the smoothing spline routine. */
+ std::vector<double> m_Scratch; // temp storage for processing of Coefficients
+ typename TInputImage::SizeType m_DataLength; // Image size
+ unsigned int m_SplineOrder; // User specified spline order (3rd or cubic is the default)
+ double m_SplinePoles[3];// Poles calculated for a given spline order
+ int m_NumberOfPoles; // number of poles
+ double m_Tolerance; // Tolerance used for determining initial causal coefficient
+ unsigned int m_IteratorDirection; // Direction for iterator incrementing
+
+
+private:
+ BSplineDecompositionImageFilter( const Self& ); //purposely not implemented
+ void operator=( const Self& ); //purposely not implemented
+
+ /** Determines the poles given the Spline Order. */
+ virtual void SetPoles();
+
+ /** Converts a vector of data to a vector of Spline coefficients. */
+ virtual bool DataToCoefficients1D();
+
+ /** Converts an N-dimension image of data to an equivalent sized image
+ * of spline coefficients. */
+ void DataToCoefficientsND();
+
+ /** Determines the first coefficient for the causal filtering of the data. */
+ virtual void SetInitialCausalCoefficient(double z);
+
+ /** Determines the first coefficient for the anti-causal filtering of the data. */
+ virtual void SetInitialAntiCausalCoefficient(double z);
+
+ /** Used to initialize the Coefficients image before calculation. */
+ void CopyImageToImage();
+
+ /** Copies a vector of data from the Coefficients image to the m_Scratch vector. */
+ void CopyCoefficientsToScratch( OutputLinearIterator & );
+
+ /** Copies a vector of data from m_Scratch to the Coefficients image. */
+ void CopyScratchToCoefficients( OutputLinearIterator & );
+
+};
+
+
+} // namespace itk
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbBSplineDecompositionImageFilter.txx"
+#endif
+
+#endif
+
diff --git a/Code/BasicFilters/otbBSplineDecompositionImageFilter.txx b/Code/BasicFilters/otbBSplineDecompositionImageFilter.txx
new file mode 100755
index 0000000000000000000000000000000000000000..c94b2d293de23a62d2f6d7f452689cddc9baad88
--- /dev/null
+++ b/Code/BasicFilters/otbBSplineDecompositionImageFilter.txx
@@ -0,0 +1,369 @@
+#ifndef _itkBSplineDecompositionImageFilter_txx
+#define _itkBSplineDecompositionImageFilter_txx
+#include "itkBSplineDecompositionImageFilter.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIterator.h"
+#include "itkProgressReporter.h"
+#include "itkVector.h"
+
+namespace otb
+{
+
+/**
+ * Constructor
+ */
+template <class TInputImage, class TOutputImage>
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::BSplineDecompositionImageFilter()
+{
+ m_SplineOrder = 0;
+ int SplineOrder = 3;
+ m_Tolerance = 1e-10; // Need some guidance on this one...what is reasonable?
+ m_IteratorDirection = 0;
+ this->SetSplineOrder(SplineOrder);
+}
+
+
+/**
+ * Standard "PrintSelf" method
+ */
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::PrintSelf(
+ std::ostream& os,
+ itk::Indent indent) const
+{
+ Superclass::PrintSelf( os, indent );
+ os << indent << "Spline Order: " << m_SplineOrder << std::endl;
+
+}
+
+
+template <class TInputImage, class TOutputImage>
+bool
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::DataToCoefficients1D()
+{
+
+ // See Unser, 1993, Part II, Equation 2.5,
+ // or Unser, 1999, Box 2. for an explaination.
+
+ double c0 = 1.0;
+
+ if (m_DataLength[m_IteratorDirection] == 1) //Required by mirror boundaries
+ {
+ return false;
+ }
+
+ // Compute overall gain
+ for (int k = 0; k < m_NumberOfPoles; k++)
+ {
+ // Note for cubic splines lambda = 6
+ c0 = c0 * (1.0 - m_SplinePoles[k]) * (1.0 - 1.0 / m_SplinePoles[k]);
+ }
+
+ // apply the gain
+ for (unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; n++)
+ {
+ m_Scratch[n] *= c0;
+ }
+
+ // loop over all poles
+ for (int k = 0; k < m_NumberOfPoles; k++)
+ {
+ // causal initialization
+ this->SetInitialCausalCoefficient(m_SplinePoles[k]);
+ // causal recursion
+ for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; n++)
+ {
+ m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1];
+ }
+
+ // anticausal initialization
+ this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]);
+ // anticausal recursion
+ for ( int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--)
+ {
+ m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]);
+ }
+ }
+ return true;
+
+}
+
+
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::SetSplineOrder(unsigned int SplineOrder)
+{
+ if (SplineOrder == m_SplineOrder)
+ {
+ return;
+ }
+ m_SplineOrder = SplineOrder;
+ this->SetPoles();
+ this->Modified();
+
+}
+
+
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::SetPoles()
+{
+ /* See Unser, 1997. Part II, Table I for Pole values */
+ // See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman,
+ // 2000, pg. 416.
+ switch (m_SplineOrder)
+ {
+ case 3:
+ m_NumberOfPoles = 1;
+ m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0;
+ break;
+ case 0:
+ m_NumberOfPoles = 0;
+ break;
+ case 1:
+ m_NumberOfPoles = 0;
+ break;
+ case 2:
+ m_NumberOfPoles = 1;
+ m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0;
+ break;
+ case 4:
+ m_NumberOfPoles = 2;
+ m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
+ m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
+ break;
+ case 5:
+ m_NumberOfPoles = 2;
+ m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
+ - 13.0 / 2.0;
+ m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0)
+ - 13.0 / 2.0;
+ break;
+ default:
+ // SplineOrder not implemented yet.
+ itk::ExceptionObject err(__FILE__, __LINE__);
+ err.SetLocation( ITK_LOCATION);
+ err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
+ throw err;
+ break;
+ }
+}
+
+
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::SetInitialCausalCoefficient(double z)
+{
+ /* begining InitialCausalCoefficient */
+ /* See Unser, 1999, Box 2 for explaination */
+
+ double sum, zn, z2n, iz;
+ unsigned long horizon;
+
+ /* this initialization corresponds to mirror boundaries */
+ horizon = m_DataLength[m_IteratorDirection];
+ zn = z;
+ if (m_Tolerance > 0.0)
+ {
+ horizon = (long)vcl_ceil(log(m_Tolerance) / vcl_log(fabs(z)));
+ }
+ if (horizon < m_DataLength[m_IteratorDirection])
+ {
+ /* accelerated loop */
+ sum = m_Scratch[0]; // verify this
+ for (unsigned int n = 1; n < horizon; n++)
+ {
+ sum += zn * m_Scratch[n];
+ zn *= z;
+ }
+ m_Scratch[0] = sum;
+ }
+ else {
+ /* full loop */
+ iz = 1.0 / z;
+ z2n = vcl_pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
+ sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
+ z2n *= z2n * iz;
+ for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++)
+ {
+ sum += (zn + z2n) * m_Scratch[n];
+ zn *= z;
+ z2n *= iz;
+ }
+ m_Scratch[0] = sum / (1.0 - zn * zn);
+ }
+}
+
+
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::SetInitialAntiCausalCoefficient(double z)
+{
+ // this initialization corresponds to mirror boundaries
+ /* See Unser, 1999, Box 2 for explaination */
+ // Also see erratum at http://bigwww.epfl.ch/publications/unser9902.html
+ m_Scratch[m_DataLength[m_IteratorDirection] - 1] =
+ (z / (z * z - 1.0)) *
+ (z * m_Scratch[m_DataLength[m_IteratorDirection] - 2] + m_Scratch[m_DataLength[m_IteratorDirection] - 1]);
+}
+
+
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::DataToCoefficientsND()
+{
+ OutputImagePointer output = this->GetOutput();
+
+ itk::Size<ImageDimension> size = output->GetBufferedRegion().GetSize();
+
+ unsigned int count = output->GetBufferedRegion().GetNumberOfPixels() / size[0] * ImageDimension;
+
+ itk::ProgressReporter progress(this, 0, count, 10);
+
+ // Initialize coeffient array
+ this->CopyImageToImage(); // Coefficients are initialized to the input data
+
+ for (unsigned int n=0; n < ImageDimension; n++)
+ {
+ m_IteratorDirection = n;
+ // Loop through each dimension
+
+ // Initialize iterators
+ OutputLinearIterator CIterator( output, output->GetBufferedRegion() );
+ CIterator.SetDirection( m_IteratorDirection );
+ // For each data vector
+ while ( !CIterator.IsAtEnd() )
+ {
+ // Copy coefficients to scratch
+ this->CopyCoefficientsToScratch( CIterator );
+
+
+ // Perform 1D BSpline calculations
+ this->DataToCoefficients1D();
+
+ // Copy scratch back to coefficients.
+ // Brings us back to the end of the line we were working on.
+ CIterator.GoToBeginOfLine();
+ this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
+ CIterator.NextLine();
+ progress.CompletedPixel();
+ }
+ }
+}
+
+
+/**
+ * Copy the input image into the output image
+ */
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::CopyImageToImage()
+{
+
+ typedef itk::ImageRegionConstIteratorWithIndex< TInputImage > InputIterator;
+ typedef itk::ImageRegionIterator< TOutputImage > OutputIterator;
+ typedef typename TOutputImage::PixelType OutputPixelType;
+
+ InputIterator inIt( this->GetInput(), this->GetInput()->GetBufferedRegion() );
+ OutputIterator outIt( this->GetOutput(), this->GetOutput()->GetBufferedRegion() );
+
+ inIt = inIt.Begin();
+ outIt = outIt.Begin();
+
+ while ( !outIt.IsAtEnd() )
+ {
+ outIt.Set( static_cast<OutputPixelType>( inIt.Get() ) );
+ ++inIt;
+ ++outIt;
+ }
+
+}
+
+
+/**
+ * Copy the scratch to one line of the output image
+ */
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::CopyScratchToCoefficients(OutputLinearIterator & Iter)
+{
+ typedef typename TOutputImage::PixelType OutputPixelType;
+ unsigned long j = 0;
+ while ( !Iter.IsAtEndOfLine() )
+ {
+ Iter.Set( static_cast<OutputPixelType>( m_Scratch[j] ) );
+ ++Iter;
+ ++j;
+ }
+
+}
+
+
+/**
+ * Copy one line of the output image to the scratch
+ */
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::CopyCoefficientsToScratch(OutputLinearIterator & Iter)
+{
+ unsigned long j = 0;
+ while ( !Iter.IsAtEndOfLine() )
+ {
+ m_Scratch[j] = static_cast<double>( Iter.Get() ) ;
+ ++Iter;
+ ++j;
+ }
+}
+
+/**
+ * Generate data
+ */
+template <class TInputImage, class TOutputImage>
+void
+BSplineDecompositionImageFilter<TInputImage, TOutputImage>
+::GenerateData()
+{
+
+ // Allocate scratch memory
+ InputImageConstPointer inputPtr = this->GetInput();
+ m_DataLength = inputPtr->GetBufferedRegion().GetSize();
+
+ unsigned long maxLength = 0;
+ for ( unsigned int n = 0; n < ImageDimension; n++ )
+ {
+ if ( m_DataLength[n] > maxLength )
+ {
+ maxLength = m_DataLength[n];
+ }
+ }
+ m_Scratch.resize( maxLength );
+
+ // Allocate memory for output image
+ OutputImagePointer outputPtr = this->GetOutput();
+ outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
+ outputPtr->Allocate();
+
+ // Calculate actual output
+ this->DataToCoefficientsND();
+
+ // Clean up
+ m_Scratch.clear();
+
+}
+
+
+} // namespace otb
+
+#endif
diff --git a/Code/BasicFilters/otbBSplineInterpolateImageFunction.h b/Code/BasicFilters/otbBSplineInterpolateImageFunction.h
new file mode 100755
index 0000000000000000000000000000000000000000..cf8c5e4868e43369d1e17ed94ae64789a109676e
--- /dev/null
+++ b/Code/BasicFilters/otbBSplineInterpolateImageFunction.h
@@ -0,0 +1,199 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbBSplineInterpolateImageFunction_h
+#define __otbBSplineInterpolateImageFunction_h
+
+#include <vector>
+
+#include "itkImageLinearIteratorWithIndex.h"
+#include "itkInterpolateImageFunction.h"
+#include "vnl/vnl_matrix.h"
+
+#include "otbBSplineDecompositionImageFilter.h"
+#include "itkConceptChecking.h"
+#include "itkCovariantVector.h"
+
+namespace otb
+{
+/** \class BSplineInterpolateImageFunction
+ * \brief This class is an evolution of the itk::BSplineInterpolateImageFunction to handle
+ * huge images with this interpolator. For more documentation, please refer to the original
+ * class.
+ * \sa itk::BSplineInterpolateImageFunction
+ * \sa itk::BSplineDecompositionImageFilter
+ * \sa otb::BSplineDecompositionImageFilter
+ *
+ * \ingroup ImageFunctions
+ */
+template <
+ class TImageType,
+ class TCoordRep = double,
+ class TCoefficientType = double >
+class ITK_EXPORT BSplineInterpolateImageFunction :
+ public itk::InterpolateImageFunction<TImageType,TCoordRep>
+{
+public:
+ /** Standard class typedefs. */
+ typedef BSplineInterpolateImageFunction Self;
+ typedef itk::InterpolateImageFunction<TImageType,TCoordRep> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(BSplineInterpolateImageFunction, InterpolateImageFunction);
+
+
+ /** New macro for creation of through a Smart Pointer */
+ itkNewMacro( Self );
+
+ /** OutputType typedef support. */
+ typedef typename Superclass::OutputType OutputType;
+
+ /** InputImageType typedef support. */
+ typedef typename Superclass::InputImageType InputImageType;
+
+ /** Dimension underlying input image. */
+ itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension);
+
+ /** Index typedef support. */
+ typedef typename Superclass::IndexType IndexType;
+
+/** Region typedef support */
+typedef typename InputImageType::RegionType RegionType;
+
+ /** ContinuousIndex typedef support. */
+ typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
+
+ /** PointType typedef support */
+ typedef typename Superclass::PointType PointType;
+
+ /** Iterator typedef support */
+ typedef itk::ImageLinearIteratorWithIndex<TImageType> Iterator;
+
+ /** Internal Coefficient typedef support */
+ typedef TCoefficientType CoefficientDataType;
+ typedef itk::Image<CoefficientDataType,
+ itkGetStaticConstMacro(ImageDimension)
+ > CoefficientImageType;
+
+ /** Define filter for calculating the BSpline coefficients */
+ typedef otb::BSplineDecompositionImageFilter<TImageType, CoefficientImageType>
+ CoefficientFilter;
+ typedef typename CoefficientFilter::Pointer CoefficientFilterPointer;
+
+ /** Evaluate the function at a ContinuousIndex position.
+ *
+ * Returns the B-Spline interpolated image intensity at a
+ * specified point position. No bounds checking is done.
+ * The point is assume to lie within the image buffer.
+ *
+ * ImageFunction::IsInsideBuffer() can be used to check bounds before
+ * calling the method. */
+ virtual OutputType EvaluateAtContinuousIndex(
+ const ContinuousIndexType & index ) const;
+
+ /** Derivative typedef support */
+ typedef itk::CovariantVector<OutputType,
+ itkGetStaticConstMacro(ImageDimension)
+ > CovariantVectorType;
+
+ CovariantVectorType EvaluateDerivative( const PointType & point ) const
+ {
+ ContinuousIndexType index;
+ this->GetInputImage()->TransformPhysicalPointToContinuousIndex( point, index );
+ return ( this->EvaluateDerivativeAtContinuousIndex( index ) );
+ }
+
+ CovariantVectorType EvaluateDerivativeAtContinuousIndex(
+ const ContinuousIndexType & x ) const;
+
+
+ /** Get/Sets the Spline Order, supports 0th - 5th order splines. The default
+ * is a 3rd order spline. */
+ void SetSplineOrder(unsigned int SplineOrder);
+ itkGetMacro(SplineOrder, int);
+
+
+ /** Set the input image. This must be set by the user. */
+ virtual void SetInputImage(const TImageType * inputData);
+
+
+ /** Update coefficients filter. Coefficient filter are computed over the buffered
+ region of the input image. */
+virtual void UpdateCoefficientsFilter(void);
+
+protected:
+ BSplineInterpolateImageFunction();
+ virtual ~BSplineInterpolateImageFunction() {};
+ void operator=( const Self& ); //purposely not implemented
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+ // These are needed by the smoothing spline routine.
+ std::vector<CoefficientDataType> m_Scratch; // temp storage for processing of Coefficients
+ typename TImageType::SizeType m_DataLength; // Image size
+ unsigned int m_SplineOrder; // User specified spline order (3rd or cubic is the default)
+
+ typename CoefficientImageType::ConstPointer m_Coefficients; // Spline coefficients
+
+private:
+ BSplineInterpolateImageFunction( const Self& ); //purposely not implemented
+ /** Determines the weights for interpolation of the value x */
+ void SetInterpolationWeights( const ContinuousIndexType & x,
+ const vnl_matrix<long> & EvaluateIndex,
+ vnl_matrix<double> & weights,
+ unsigned int splineOrder ) const;
+
+ /** Determines the weights for the derivative portion of the value x */
+ void SetDerivativeWeights( const ContinuousIndexType & x,
+ const vnl_matrix<long> & EvaluateIndex,
+ vnl_matrix<double> & weights,
+ unsigned int splineOrder ) const;
+
+ /** Precomputation for converting the 1D index of the interpolation neighborhood
+ * to an N-dimensional index. */
+ void GeneratePointsToIndex( );
+
+ /** Determines the indicies to use give the splines region of support */
+ void DetermineRegionOfSupport( vnl_matrix<long> & evaluateIndex,
+ const ContinuousIndexType & x,
+ unsigned int splineOrder ) const;
+
+ /** Set the indicies in evaluateIndex at the boundaries based on mirror
+ * boundary conditions. */
+ void ApplyMirrorBoundaryConditions(vnl_matrix<long> & evaluateIndex,
+ unsigned int splineOrder) const;
+
+
+ Iterator m_CIterator; // Iterator for traversing spline coefficients.
+ unsigned long m_MaxNumberInterpolationPoints; // number of neighborhood points used for interpolation
+ std::vector<IndexType> m_PointsToIndex; // Preallocation of interpolation neighborhood indicies
+
+ CoefficientFilterPointer m_CoefficientFilter;
+
+ RegionType m_CurrentBufferedRegion;
+
+};
+
+} // namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbBSplineInterpolateImageFunction.txx"
+#endif
+
+#endif
+
diff --git a/Code/BasicFilters/otbBSplineInterpolateImageFunction.txx b/Code/BasicFilters/otbBSplineInterpolateImageFunction.txx
new file mode 100755
index 0000000000000000000000000000000000000000..38855492a3ae5a3f3a5efb6225b9a9281b5b6e63
--- /dev/null
+++ b/Code/BasicFilters/otbBSplineInterpolateImageFunction.txx
@@ -0,0 +1,579 @@
+/*=========================================================================
+
+ Program: Insight Segmentation & Registration Toolkit
+ Module: $RCSfile: itkBSplineInterpolateImageFunction.txx,v $
+ Language: C++
+ Date: $Date: 2006/03/19 04:36:55 $
+ Version: $Revision: 1.13 $
+
+ Copyright (c) Insight Software Consortium. All rights reserved.
+ See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
+
+ Portions of this code are covered under the VTK copyright.
+ See VTKCopyright.txt or http://www.kitware.com/VTKCopyright.htm for details.
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef _otbBSplineInterpolateImageFunction_txx
+#define _otbBSplineInterpolateImageFunction_txx
+#include "itkBSplineInterpolateImageFunction.h"
+#include "itkImageLinearIteratorWithIndex.h"
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIterator.h"
+
+#include "itkVector.h"
+
+#include "itkMatrix.h"
+
+namespace otb
+{
+
+/**
+ * Constructor
+ */
+template <class TImageType, class TCoordRep, class TCoefficientType>
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::BSplineInterpolateImageFunction()
+{
+ m_SplineOrder = 0;
+ unsigned int SplineOrder = 3;
+ m_CoefficientFilter = CoefficientFilter::New();
+ // ***TODO: Should we store coefficients in a variable or retrieve from filter?
+ m_Coefficients = CoefficientImageType::New();
+ this->SetSplineOrder(SplineOrder);
+}
+
+/**
+ * Standard "PrintSelf" method
+ */
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::PrintSelf(
+ std::ostream& os,
+ itk::Indent indent) const
+{
+ Superclass::PrintSelf( os, indent );
+ os << indent << "Spline Order: " << m_SplineOrder << std::endl;
+
+}
+
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::UpdateCoefficientsFilter(void)
+{
+ if(!m_CoefficientFilter->GetInput()->GetBufferedRegion().IsInside(m_CurrentBufferedRegion))
+ {
+ //std::cout<<"Updating coefs filter"<<std::endl;
+ m_CoefficientFilter->GetOutput()->UpdateOutputInformation();
+ m_CoefficientFilter->GetOutput()->SetRequestedRegion(m_CoefficientFilter->GetInput()->GetBufferedRegion());
+ m_CoefficientFilter->GetOutput()->PropagateRequestedRegion();
+ m_CoefficientFilter->GetOutput()->UpdateOutputData();
+ m_Coefficients = m_CoefficientFilter->GetOutput();
+ //std::cout<<"Updated region: "<<m_CoefficientFilter->GetInput()->GetBufferedRegion()<<std::endl;
+ }
+ m_CurrentBufferedRegion =m_CoefficientFilter->GetInput()->GetBufferedRegion();
+}
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::SetInputImage(const TImageType * inputData)
+{
+ if ( inputData )
+ {
+ m_CoefficientFilter->SetInput(inputData);
+
+ // the Coefficient Filter requires that the spline order and the input data be set.
+ // TODO: We need to ensure that this is only run once and only after both input and
+ // spline order have been set. Should we force an update after the
+ // splineOrder has been set also?
+
+
+
+ UpdateCoefficientsFilter();
+
+ // Call the Superclass implementation after, in case the filter
+ // pulls in more of the input image
+ Superclass::SetInputImage(inputData);
+
+ m_DataLength = inputData->GetBufferedRegion().GetSize();
+ }
+ else
+ {
+ m_Coefficients = NULL;
+ }
+}
+
+
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::SetSplineOrder(unsigned int SplineOrder)
+{
+ if (SplineOrder == m_SplineOrder)
+ {
+ return;
+ }
+ m_SplineOrder = SplineOrder;
+ m_CoefficientFilter->SetSplineOrder( SplineOrder );
+
+ //this->SetPoles();
+ m_MaxNumberInterpolationPoints = 1;
+ for (unsigned int n=0; n < ImageDimension; n++)
+ {
+ m_MaxNumberInterpolationPoints *= ( m_SplineOrder + 1);
+ }
+ this->GeneratePointsToIndex( );
+}
+
+
+template <class TImageType, class TCoordRep, class TCoefficientType>
+typename
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::OutputType
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::EvaluateAtContinuousIndex( const ContinuousIndexType & x ) const
+{
+ //UpdateCoefficientsFilter();
+ vnl_matrix<long> EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
+
+ // compute the interpolation indexes
+ this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
+
+ // Determine weights
+ vnl_matrix<double> weights(ImageDimension, ( m_SplineOrder + 1 ));
+ SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
+
+ //std::cout<<"EvaluateIndex: "<<std::endl;
+ //std::cout<<EvaluateIndex[0][0]<<"\t"<<EvaluateIndex[0][1]<<"\t"<<EvaluateIndex[0][2]<<"\t"<<EvaluateIndex[0][3]<<std::endl;
+ //std::cout<<EvaluateIndex[1][0]<<"\t"<<EvaluateIndex[1][1]<<"\t"<<EvaluateIndex[1][2]<<"\t"<<EvaluateIndex[1][3]<<std::endl;
+
+ // Modify EvaluateIndex at the boundaries using mirror boundary conditions
+ this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
+
+ // perform interpolation
+ double interpolated = 0.0;
+ IndexType coefficientIndex;
+ // Step through eachpoint in the N-dimensional interpolation cube.
+ for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
+ {
+ // translate each step into the N-dimensional index.
+ // IndexType pointIndex = PointToIndex( p );
+
+ double w = 1.0;
+ for (unsigned int n = 0; n < ImageDimension; n++ )
+ {
+ w *= weights[n][ m_PointsToIndex[p][n] ];
+ coefficientIndex[n] = EvaluateIndex[n][m_PointsToIndex[p][n]]; // Build up ND index for coefficients.
+ //std::cout<<"From inside: "<<n<<" "<<p<<" "<<m_PointsToIndex[p][n]<<" "<< EvaluateIndex[n][m_PointsToIndex[p][n]]<<std::endl;
+ }
+ //std::cout<<"CoefficientIndex: "<<coefficientIndex<<std::endl;
+ // Convert our step p to the appropriate point in ND space in the
+ // m_Coefficients cube.
+ interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
+ }
+
+/* double interpolated = 0.0;
+ IndexType coefficientIndex;
+ // Step through eachpoint in the N-dimensional interpolation cube.
+ for (unsigned int sp = 0; sp <= m_SplineOrder; sp++)
+ {
+ for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++)
+ {
+
+ double w = 1.0;
+ for (unsigned int n1 = 0; n1 < ImageDimension; n1++ )
+ {
+ w *= weights[n1][ sp1 ];
+ coefficientIndex[n1] = EvaluateIndex[n1][sp]; // Build up ND index for coefficients.
+ }
+
+ interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
+ }
+ }
+*/
+ return(interpolated);
+
+}
+
+
+template <class TImageType, class TCoordRep, class TCoefficientType>
+typename
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+:: CovariantVectorType
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::EvaluateDerivativeAtContinuousIndex( const ContinuousIndexType & x ) const
+{
+ UpdateCoefficientsFilter();
+ vnl_matrix<long> EvaluateIndex(ImageDimension, ( m_SplineOrder + 1 ));
+
+ // compute the interpolation indexes
+ // TODO: Do we need to revisit region of support for the derivatives?
+ this->DetermineRegionOfSupport(EvaluateIndex, x, m_SplineOrder);
+
+ // Determine weights
+ vnl_matrix<double> weights(ImageDimension, ( m_SplineOrder + 1 ));
+ SetInterpolationWeights( x, EvaluateIndex, weights, m_SplineOrder );
+
+ vnl_matrix<double> weightsDerivative(ImageDimension, ( m_SplineOrder + 1));
+ SetDerivativeWeights( x, EvaluateIndex, weightsDerivative, ( m_SplineOrder ) );
+
+ // Modify EvaluateIndex at the boundaries using mirror boundary conditions
+ this->ApplyMirrorBoundaryConditions(EvaluateIndex, m_SplineOrder);
+
+ // Calculate derivative
+ CovariantVectorType derivativeValue;
+ double tempValue;
+ IndexType coefficientIndex;
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ derivativeValue[n] = 0.0;
+ for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
+ {
+ tempValue = 1.0 ;
+ for (unsigned int n1 = 0; n1 < ImageDimension; n1++)
+ {
+ //coefficientIndex[n1] = EvaluateIndex[n1][sp];
+ coefficientIndex[n1] = EvaluateIndex[n1][m_PointsToIndex[p][n1]];
+
+ if (n1 == n)
+ {
+ //w *= weights[n][ m_PointsToIndex[p][n] ];
+ tempValue *= weightsDerivative[n1][ m_PointsToIndex[p][n1] ];
+ }
+ else
+ {
+ tempValue *= weights[n1][ m_PointsToIndex[p][n1] ];
+ }
+ }
+ derivativeValue[n] += m_Coefficients->GetPixel(coefficientIndex) * tempValue ;
+ }
+ derivativeValue[n] /= this->GetInputImage()->GetSpacing()[n]; // take spacing into account
+ }
+
+ return(derivativeValue);
+
+}
+
+
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::SetInterpolationWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
+ vnl_matrix<double> & weights, unsigned int splineOrder ) const
+{
+ // For speed improvements we could make each case a separate function and use
+ // function pointers to reference the correct weight order.
+ // Left as is for now for readability.
+ double w, w2, w4, t, t0, t1;
+
+ switch (splineOrder)
+ {
+ case 3:
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ w = x[n] - (double) EvaluateIndex[n][1];
+ weights[n][3] = (1.0 / 6.0) * w * w * w;
+ weights[n][0] = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - weights[n][3];
+ weights[n][2] = w + weights[n][0] - 2.0 * weights[n][3];
+ weights[n][1] = 1.0 - weights[n][0] - weights[n][2] - weights[n][3];
+ }
+ break;
+ case 0:
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ weights[n][0] = 1; // implements nearest neighbor
+ }
+ break;
+ case 1:
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ w = x[n] - (double) EvaluateIndex[n][0];
+ weights[n][1] = w;
+ weights[n][0] = 1.0 - w;
+ }
+ break;
+ case 2:
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ /* x */
+ w = x[n] - (double)EvaluateIndex[n][1];
+ weights[n][1] = 0.75 - w * w;
+ weights[n][2] = 0.5 * (w - weights[n][1] + 1.0);
+ weights[n][0] = 1.0 - weights[n][1] - weights[n][2];
+ }
+ break;
+ case 4:
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ /* x */
+ w = x[n] - (double)EvaluateIndex[n][2];
+ w2 = w * w;
+ t = (1.0 / 6.0) * w2;
+ weights[n][0] = 0.5 - w;
+ weights[n][0] *= weights[n][0];
+ weights[n][0] *= (1.0 / 24.0) * weights[n][0];
+ t0 = w * (t - 11.0 / 24.0);
+ t1 = 19.0 / 96.0 + w2 * (0.25 - t);
+ weights[n][1] = t1 + t0;
+ weights[n][3] = t1 - t0;
+ weights[n][4] = weights[n][0] + t0 + 0.5 * w;
+ weights[n][2] = 1.0 - weights[n][0] - weights[n][1] - weights[n][3] - weights[n][4];
+ }
+ break;
+ case 5:
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ /* x */
+ w = x[n] - (double)EvaluateIndex[n][2];
+ w2 = w * w;
+ weights[n][5] = (1.0 / 120.0) * w * w2 * w2;
+ w2 -= w;
+ w4 = w2 * w2;
+ w -= 0.5;
+ t = w2 * (w2 - 3.0);
+ weights[n][0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - weights[n][5];
+ t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
+ t1 = (-1.0 / 12.0) * w * (t + 4.0);
+ weights[n][2] = t0 + t1;
+ weights[n][3] = t0 - t1;
+ t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
+ t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
+ weights[n][1] = t0 + t1;
+ weights[n][4] = t0 - t1;
+ }
+ break;
+ default:
+ // SplineOrder not implemented yet.
+ itk::ExceptionObject err(__FILE__, __LINE__);
+ err.SetLocation( ITK_LOCATION );
+ err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
+ throw err;
+ break;
+ }
+
+}
+
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::SetDerivativeWeights( const ContinuousIndexType & x, const vnl_matrix<long> & EvaluateIndex,
+ vnl_matrix<double> & weights, unsigned int splineOrder ) const
+{
+ // For speed improvements we could make each case a separate function and use
+ // function pointers to reference the correct weight order.
+ // Another possiblity would be to loop inside the case statement (reducing the number
+ // of switch statement executions to one per routine call.
+ // Left as is for now for readability.
+ double w, w1, w2, w3, w4, w5, t, t0, t1, t2;
+ int derivativeSplineOrder = (int) splineOrder -1;
+
+ switch (derivativeSplineOrder)
+ {
+
+ // Calculates B(splineOrder) ( (x + 1/2) - xi) - B(splineOrder -1) ( (x - 1/2) - xi)
+ case -1:
+ // Why would we want to do this?
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ weights[n][0] = 0.0;
+ }
+ break;
+ case 0:
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ weights[n][0] = -1.0;
+ weights[n][1] = 1.0;
+ }
+ break;
+ case 1:
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ w = x[n] + 0.5 - (double)EvaluateIndex[n][1];
+ // w2 = w;
+ w1 = 1.0 - w;
+
+ weights[n][0] = 0.0 - w1;
+ weights[n][1] = w1 - w;
+ weights[n][2] = w;
+ }
+ break;
+ case 2:
+
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ w = x[n] + .5 - (double)EvaluateIndex[n][2];
+ w2 = 0.75 - w * w;
+ w3 = 0.5 * (w - w2 + 1.0);
+ w1 = 1.0 - w2 - w3;
+
+ weights[n][0] = 0.0 - w1;
+ weights[n][1] = w1 - w2;
+ weights[n][2] = w2 - w3;
+ weights[n][3] = w3;
+ }
+ break;
+ case 3:
+
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ w = x[n] + 0.5 - (double)EvaluateIndex[n][2];
+ w4 = (1.0 / 6.0) * w * w * w;
+ w1 = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - w4;
+ w3 = w + w1 - 2.0 * w4;
+ w2 = 1.0 - w1 - w3 - w4;
+
+ weights[n][0] = 0.0 - w1;
+ weights[n][1] = w1 - w2;
+ weights[n][2] = w2 - w3;
+ weights[n][3] = w3 - w4;
+ weights[n][4] = w4;
+ }
+ break;
+ case 4:
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ w = x[n] + .5 - (double)EvaluateIndex[n][3];
+ t2 = w * w;
+ t = (1.0 / 6.0) * t2;
+ w1 = 0.5 - w;
+ w1 *= w1;
+ w1 *= (1.0 / 24.0) * w1;
+ t0 = w * (t - 11.0 / 24.0);
+ t1 = 19.0 / 96.0 + t2 * (0.25 - t);
+ w2 = t1 + t0;
+ w4 = t1 - t0;
+ w5 = w1 + t0 + 0.5 * w;
+ w3 = 1.0 - w1 - w2 - w4 - w5;
+
+ weights[n][0] = 0.0 - w1;
+ weights[n][1] = w1 - w2;
+ weights[n][2] = w2 - w3;
+ weights[n][3] = w3 - w4;
+ weights[n][4] = w4 - w5;
+ weights[n][5] = w5;
+ }
+ break;
+
+ default:
+ // SplineOrder not implemented yet.
+ itk::ExceptionObject err(__FILE__, __LINE__);
+ err.SetLocation( ITK_LOCATION );
+ err.SetDescription( "SplineOrder (for derivatives) must be between 1 and 5. Requested spline order has not been implemented yet." );
+ throw err;
+ break;
+ }
+
+}
+
+
+// Generates m_PointsToIndex;
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::GeneratePointsToIndex( )
+{
+ // m_PointsToIndex is used to convert a sequential location to an N-dimension
+ // index vector. This is precomputed to save time during the interpolation routine.
+ m_PointsToIndex.resize(m_MaxNumberInterpolationPoints);
+ for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
+ {
+ int pp = p;
+ unsigned long indexFactor[ImageDimension];
+ indexFactor[0] = 1;
+ for (int j=1; j< static_cast<int>(ImageDimension); j++)
+ {
+ indexFactor[j] = indexFactor[j-1] * ( m_SplineOrder + 1 );
+ }
+ for (int j = (static_cast<int>(ImageDimension) - 1); j >= 0; j--)
+ {
+ m_PointsToIndex[p][j] = pp / indexFactor[j];
+ pp = pp % indexFactor[j];
+ }
+ }
+}
+
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::DetermineRegionOfSupport( vnl_matrix<long> & evaluateIndex,
+ const ContinuousIndexType & x,
+ unsigned int splineOrder ) const
+{
+ long indx;
+
+// compute the interpolation indexes
+ for (unsigned int n = 0; n< ImageDimension; n++)
+ {
+ if (splineOrder & 1) // Use this index calculation for odd splineOrder
+ {
+ indx = (long)vcl_floor((float)x[n]) - splineOrder / 2;
+ //std::cout<<"x: "<<x<<std::endl;
+ //std::cout<<"splineOrder: "<<splineOrder<<std::endl;
+ //std::cout<<"indx: "<<indx<<std::endl;
+ for (unsigned int k = 0; k <= splineOrder; k++)
+ {
+ evaluateIndex[n][k] = indx++;
+ }
+ }
+ else // Use this index calculation for even splineOrder
+ {
+
+ indx = (long)vcl_floor((float)(x[n] + 0.5)) - splineOrder / 2;
+ //std::cout<<"x: "<<x<<std::endl;
+ //std::cout<<"splineOrder: "<<splineOrder<<std::endl;
+ //std::cout<<"indx: "<<indx<<std::endl;
+ for (unsigned int k = 0; k <= splineOrder; k++)
+ {
+ evaluateIndex[n][k] = indx++;
+ }
+ }
+ }
+}
+
+template <class TImageType, class TCoordRep, class TCoefficientType>
+void
+BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
+::ApplyMirrorBoundaryConditions(vnl_matrix<long> & evaluateIndex,
+ unsigned int splineOrder) const
+{
+
+ for (unsigned int n = 0; n < ImageDimension; n++)
+ {
+ long dataLength = m_DataLength[n];
+ long dataOffset = m_CurrentBufferedRegion.GetIndex()[n];
+
+ // apply the mirror boundary conditions
+ // TODO: We could implement other boundary options beside mirror
+ if (m_DataLength[n] == 1)
+ {
+ for (unsigned int k = 0; k <= splineOrder; k++)
+ {
+ evaluateIndex[n][k] = 0;
+ }
+ }
+ else
+ {
+ for (unsigned int k = 0; k <= splineOrder; k++)
+ {
+ // btw - Think about this couldn't this be replaced with a more elagent modulus method?
+ evaluateIndex[n][k] = (evaluateIndex[n][k] < dataOffset) ? (dataOffset+(dataOffset-evaluateIndex[n][k])%dataLength)
+ : (evaluateIndex[n][k]);
+ if ((long) dataLength+dataOffset <= evaluateIndex[n][k])
+ {
+ evaluateIndex[n][k] = dataOffset + dataLength - (evaluateIndex[n][k]-dataOffset - dataLength)%dataLength;
+ }
+ }
+
+ }
+
+ }
+}
+
+} // namespace otb
+
+#endif
diff --git a/Code/Visu/otbImageAlternateViewer.h b/Code/Visu/otbImageAlternateViewer.h
index 480479aea8e31dad1c1af1f505c3d0a0149332b2..f6e518f23f08d8a8ea61ce622ef0bde00734d955 100644
--- a/Code/Visu/otbImageAlternateViewer.h
+++ b/Code/Visu/otbImageAlternateViewer.h
@@ -23,7 +23,7 @@ PURPOSE. See the above copyright notices for more information.
#include "otbImageList.h"
#include "itkInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"
-#include "itkBSplineInterpolateImageFunction.h"
+#include "otbBSplineInterpolateImageFunction.h"
#include "itkVectorImageToImageAdaptor.h"
#include "otbVectorImageToImageListFilter.h"
#include "otbImage.h"
@@ -76,7 +76,7 @@ class ITK_EXPORT ImageAlternateViewer
typedef itk::InterpolateImageFunction<SingleImageType,double> InterpolatorType;
typedef typename InterpolatorType::Pointer InterpolatorPointerType;
typedef itk::LinearInterpolateImageFunction<SingleImageType,double> DefaultInterpolatorType;
- typedef itk::BSplineInterpolateImageFunction<SingleImageType,double,double> BSplineInterpolatorType;
+ typedef BSplineInterpolateImageFunction<SingleImageType,double,double> BSplineInterpolatorType;
typedef itk::Function::HammingWindowFunction<4> WindowFunctionType;
typedef itk::ZeroFluxNeumannBoundaryCondition<SingleImageType> ConditionType;
diff --git a/Code/Visu/otbImageAlternateViewer.txx b/Code/Visu/otbImageAlternateViewer.txx
index 87f6585f5962ea5ed961d3dfd48de94885b60332..7df6bf217d0d3e755bdcce31a99bf5040a634bfc 100644
--- a/Code/Visu/otbImageAlternateViewer.txx
+++ b/Code/Visu/otbImageAlternateViewer.txx
@@ -56,10 +56,10 @@ namespace otb
// m_ZoomInInterpolator = WindowedSincInterpolatorType::New();
- m_ZoomInInterpolator=DefaultInterpolatorType::New();
+ m_ZoomInInterpolator=bsplineInterpolator;
// m_ZoomInInterpolator=bsplineInterpolator;
- m_ZoomOutInterpolator = DefaultInterpolatorType::New();
+ m_ZoomOutInterpolator = bsplineInterpolator;
IndexType index;
SizeType size;
diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt
index d854eee6e09dcbd8f728462948d9d39ff9a3d930..8ab2d36d8b15376cc6f5b30f41a4e50023226f5f 100755
--- a/Testing/Code/BasicFilters/CMakeLists.txt
+++ b/Testing/Code/BasicFilters/CMakeLists.txt
@@ -469,6 +469,43 @@ ADD_TEST(bfTvPrintableImageFilter ${BASICFILTERS_TESTS5}
${TEMP}/bfTvPrintableImageFilter.png
)
+# ------- otb::BSplineDecompositionImageFilter ------------------------------
+
+ADD_TEST(bfTuBSplineDecompositionImageFilterNew ${BASICFILTERS_TESTS6}
+ otbBSplineDecompositionImageFilterNew)
+
+ADD_TEST(bfTvBSplineDecompositionImageFilter ${BASICFILTERS_TESTS6}
+ --compare-image ${TOL}
+ ${BASELINE}/bfBSplineDecompositionImageFilterOutput.tif
+ ${TEMP}/bfBSplineDecompositionImageFilterOutput.tif
+ otbBSplineDecompositionImageFilter
+ ${INPUTDATA}/poupees.tif
+ ${TEMP}/bfBSplineDecompositionImageFilterOutput.tif
+ )
+
+# ------- otb::BSplineInterpolateImageFunction ------------------------------
+
+ADD_TEST(bfTuBSplineInterpolateImageFunctionNew ${BASICFILTERS_TESTS6}
+ otbBSplineInterpolateImageFunctionNew)
+
+ADD_TEST(bfTvBSplineInterpolateImageFunction ${BASICFILTERS_TESTS6}
+ --compare-ascii ${TOL}
+ ${BASELINE_FILES}/bfBSplineInterpolateImageFunctionOutput.txt
+ ${TEMP}/bfBSplineInterpolateImageFunctionOutput.txt
+ otbBSplineInterpolateImageFunction
+ ${INPUTDATA}/poupees.tif
+ ${TEMP}/bfBSplineInterpolateImageFunctionOutput.txt
+ 0.5 0.5
+ 127.33 44.9
+ 259.67 21.43
+ 12.13 61.79
+ 89.5 11
+ 128 128
+ 127.255 128.73
+ -1 -1
+ )
+
+
# A enrichir
SET(BasicFilters_SRCS1
otbLeeFilter.cxx
@@ -541,17 +578,10 @@ otbPrintableImageFilter.cxx
# A enrichir
SET(BasicFilters_SRCS6
-otbStreamingResampleImageFilterNew.cxx
-otbStreamingResampleImageFilter.cxx
-otbStreamingResampleImageFilterCompareWithITK.cxx
-otbStreamingStatisticsVectorImageFilterNew.cxx
-otbStreamingStatisticsVectorImageFilter.cxx
-otbMatrixTransposeMatrixImageFilterNew.cxx
-otbMatrixTransposeMatrixImageFilter.cxx
-otbUnaryImageFunctorWithVectorImageFilterNew.cxx
-otbUnaryImageFunctorWithVectorImageFilter.cxx
-otbPrintableImageFilterNew.cxx
-otbPrintableImageFilter.cxx
+otbBSplineDecompositionImageFilterNew.cxx
+otbBSplineDecompositionImageFilter.cxx
+otbBSplineInterpolateImageFunctionNew.cxx
+otbBSplineInterpolateImageFunction.cxx
)
INCLUDE_DIRECTORIES("${OTBTesting_BINARY_DIR}")
diff --git a/Testing/Code/BasicFilters/otbBSplineDecompositionImageFilter.cxx b/Testing/Code/BasicFilters/otbBSplineDecompositionImageFilter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..928dac2c5cb81c193075d966f6f2f30392d26821
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbBSplineDecompositionImageFilter.cxx
@@ -0,0 +1,50 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#include "itkExceptionObject.h"
+
+#include "otbBSplineDecompositionImageFilter.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+
+int otbBSplineDecompositionImageFilter(int argc, char * argv[])
+{
+ const char * infname = argv[1];
+ const char * outfname = argv[2];
+
+ typedef otb::Image<double,2> ImageType;
+ typedef otb::BSplineDecompositionImageFilter<ImageType,ImageType> BSplineDecompositionImageFilterType;
+ typedef otb::ImageFileReader<ImageType> ReaderType;
+ typedef otb::ImageFileWriter<ImageType> WriterType;
+
+ // Instantiating object
+ BSplineDecompositionImageFilterType::Pointer filter = BSplineDecompositionImageFilterType::New();
+
+ ReaderType::Pointer reader = ReaderType::New();
+ WriterType::Pointer writer = WriterType::New();
+
+ reader->SetFileName(infname);
+ writer->SetFileName(outfname);
+
+ filter->SetInput(reader->GetOutput());
+ writer->SetInput(filter->GetOutput());
+
+ writer->Update();
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/BasicFilters/otbBSplineDecompositionImageFilterNew.cxx b/Testing/Code/BasicFilters/otbBSplineDecompositionImageFilterNew.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..995a7110a5182ed004e5f01a89c8692ac9ebe440
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbBSplineDecompositionImageFilterNew.cxx
@@ -0,0 +1,32 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#include "itkExceptionObject.h"
+
+#include "otbBSplineDecompositionImageFilter.h"
+#include "otbImage.h"
+
+int otbBSplineDecompositionImageFilterNew(int argc, char * argv[])
+{
+ typedef otb::Image<double,2> ImageType;
+ typedef otb::BSplineDecompositionImageFilter<ImageType,ImageType> BSplineDecompositionImageFilterType;
+
+ // Instantiating object
+ BSplineDecompositionImageFilterType::Pointer filter = BSplineDecompositionImageFilterType::New();
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/BasicFilters/otbBSplineInterpolateImageFunction.cxx b/Testing/Code/BasicFilters/otbBSplineInterpolateImageFunction.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..248c432d61789882c7160c2cdef8d2e1d160c4a2
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbBSplineInterpolateImageFunction.cxx
@@ -0,0 +1,74 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#include "itkExceptionObject.h"
+
+#include "otbBSplineInterpolateImageFunction.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include <fstream>
+
+int otbBSplineInterpolateImageFunction(int argc, char * argv[])
+{
+ const char * infname = argv[1];
+ const char * outfname = argv[2];
+
+ typedef otb::Image<double,2> ImageType;
+ typedef otb::BSplineInterpolateImageFunction<ImageType,double,double> InterpolatorType;
+ typedef InterpolatorType::ContinuousIndexType ContinuousIndexType;
+ typedef otb::ImageFileReader<ImageType> ReaderType;
+
+ unsigned int i = 3;
+
+ std::vector<ContinuousIndexType> indicesList;
+
+ while(i<argc && (i+1)<argc)
+ {
+ ContinuousIndexType idx;
+
+ idx[0]=atof(argv[i]);
+ idx[1]=atof(argv[i+1]);
+
+ indicesList.push_back(idx);
+
+ i+=2;
+ }
+
+ // Instantiating object
+ InterpolatorType::Pointer interpolator = InterpolatorType::New();
+
+ ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(infname);
+ reader->Update();
+
+ interpolator->SetInputImage(reader->GetOutput());
+
+ std::ofstream file;
+ file.open(outfname);
+
+ file<<"Interpolated points in image "<<infname<<":"<<std::endl;
+ file<<std::endl;
+
+ for(std::vector<ContinuousIndexType>::iterator it = indicesList.begin();it!=indicesList.end();++it)
+ {
+ file<<(*it)<<" -> "<<interpolator->EvaluateAtContinuousIndex((*it))<<std::endl;
+ }
+
+ file.close();
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/BasicFilters/otbBSplineInterpolateImageFunctionNew.cxx b/Testing/Code/BasicFilters/otbBSplineInterpolateImageFunctionNew.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2ca067e9e8fc2c87007aa313b8c603180c554248
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbBSplineInterpolateImageFunctionNew.cxx
@@ -0,0 +1,32 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#include "itkExceptionObject.h"
+
+#include "otbBSplineInterpolateImageFunction.h"
+#include "otbImage.h"
+
+int otbBSplineInterpolateImageFunctionNew(int argc, char * argv[])
+{
+ typedef otb::Image<double,2> ImageType;
+ typedef otb::BSplineInterpolateImageFunction<ImageType,double,double> InterpolatorType;
+
+ // Instantiating object
+ InterpolatorType::Pointer filter = InterpolatorType::New();
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/BasicFilters/otbBasicFiltersTests6.cxx b/Testing/Code/BasicFilters/otbBasicFiltersTests6.cxx
index 1bbdb44d1f80cd90ed7bc1bcaedc81bff6b7ce00..b0c93c86fb6d12c64acb6fe265a566a97a0bc798 100755
--- a/Testing/Code/BasicFilters/otbBasicFiltersTests6.cxx
+++ b/Testing/Code/BasicFilters/otbBasicFiltersTests6.cxx
@@ -27,4 +27,8 @@
void RegisterTests()
{
+ REGISTER_TEST(otbBSplineDecompositionImageFilterNew);
+ REGISTER_TEST(otbBSplineDecompositionImageFilter);
+ REGISTER_TEST(otbBSplineInterpolateImageFunctionNew);
+ REGISTER_TEST(otbBSplineInterpolateImageFunction);
}