diff --git a/Code/BasicFilters/otbBCOInterpolateImageFunction.h b/Code/BasicFilters/otbBCOInterpolateImageFunction.h
new file mode 100644
index 0000000000000000000000000000000000000000..b717ea2cb297232777ff96dc2bf63383bd6ea0c8
--- /dev/null
+++ b/Code/BasicFilters/otbBCOInterpolateImageFunction.h
@@ -0,0 +1,157 @@
+/*=========================================================================
+
+ 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 __otbBCOInterpolateImageFunction_h
+#define __otbBCOInterpolateImageFunction_h
+
+#include "itkInterpolateImageFunction.h"
+#include "vnl/vnl_vector.h"
+
+namespace otb
+{
+
+/** \class
+ * \brief
+ *
+ *
+ *
+ *
+ *
+ *
+ *
+ *
+ * \warning
+ *
+ *
+ * \sa
+ *
+ * \ingroup ImageFunctions ImageInterpolators
+ */
+template <
+ class TInputImage,
+ class TCoordRep = double>
+class ITK_EXPORT BCOInterpolateImageFunction :
+ public itk::InterpolateImageFunction<TInputImage,TCoordRep>
+{
+public:
+ /** Standard class typedefs. */
+ typedef BCOInterpolateImageFunction Self;
+ typedef itk::InterpolateImageFunction<TInputImage,TCoordRep> Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(BCOInterpolateImageFunction, InterpolateImageFunction);
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** OutputType typedef support. */
+ typedef typename Superclass::OutputType OutputType;
+
+ /** InputImageType typedef support. */
+ typedef typename Superclass::InputImageType InputImageType;
+
+ /** InputPixelType typedef support. */
+ typedef typename Superclass::InputPixelType InputPixelType;
+
+ /** RealType typedef support. */
+ //typedef typename Superclass::RealType RealType;
+ typedef typename itk::NumericTraits<typename TInputImage::PixelType>::RealType RealType;
+
+ /** Dimension underlying input image. */
+ itkStaticConstMacro(ImageDimension, unsigned int,Superclass::ImageDimension);
+
+ /** Index typedef support. */
+ typedef typename Superclass::IndexType IndexType;
+ typedef typename Superclass::IndexValueType IndexValueType;
+
+ /** Point typedef support. */
+ typedef typename Superclass::PointType PointType;
+
+ /** ContinuousIndex typedef support. */
+ typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
+
+ /** Coeficients container type.*/
+ typedef vnl_vector<double> CoefContainerType;
+
+ /** Interpolate the image at a point position
+ *
+ * Returns the 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 Evaluate( const PointType& point ) const;//////////////////////////////////////////////////
+
+ /** Evaluate the function at a ContinuousIndex position
+ *
+ * Returns the linearly 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;/////////////////////
+
+ /** Interpolate the image at an index position.
+ *
+ * Simply returns the image value at the
+ * specified index 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 EvaluateAtIndex( const IndexType & index ) const;//////////////////////////////////////////
+
+ /** Set/Get the window radius */
+ virtual void SetRadius(unsigned int radius);
+ virtual unsigned int GetRadius() const;
+
+ /** Set/Get the optimisation coefficient (Common values are -0.5, -0.75 or -1.0) */
+ virtual void SetAlpha(double alpha);
+ virtual double GetAlpha() const;
+
+ /** Compute the BCO coefficients. */
+ void EvaluateCoef();
+
+protected:
+ BCOInterpolateImageFunction();
+ ~BCOInterpolateImageFunction();
+ void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+private:
+ BCOInterpolateImageFunction( const Self& ); //purposely not implemented
+ void operator=( const Self& ); //purposely not implemented
+
+ /** Used radius for the BCO */
+ double m_Radius;
+ /** Optimisation Coefficient */
+ double m_Alpha;
+ /** Used BCO coefficiet */
+ CoefContainerType m_BCOCoef;
+};
+
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbBCOInterpolateImageFunction.txx"
+#endif
+
+#endif
diff --git a/Code/BasicFilters/otbBCOInterpolateImageFunction.txx b/Code/BasicFilters/otbBCOInterpolateImageFunction.txx
new file mode 100644
index 0000000000000000000000000000000000000000..929a50d7503fd3b4812af7a4f971dbdcadf4aa94
--- /dev/null
+++ b/Code/BasicFilters/otbBCOInterpolateImageFunction.txx
@@ -0,0 +1,225 @@
+/*=========================================================================
+
+ 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 __otbProlateInterpolateImageFunction_txx
+#define __otbProlateInterpolateImageFunction_txx
+
+#include "otbBCOInterpolateImageFunction.h"
+
+#include "itkNumericTraits.h"
+
+namespace otb
+{
+
+/** Constructor */
+template <class TInputImage, class TCoordRep>
+BCOInterpolateImageFunction<TInputImage, TCoordRep>
+::BCOInterpolateImageFunction()
+{
+ m_Radius = 1;
+ m_Alpha = -0.5;
+}
+
+/** Destructor */
+template <class TInputImage, class TCoordRep>
+BCOInterpolateImageFunction<TInputImage, TCoordRep>
+::~BCOInterpolateImageFunction()
+{
+
+}
+
+template <class TInputImage, class TCoordRep>
+void BCOInterpolateImageFunction<TInputImage, TCoordRep>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf(os, indent);
+}
+
+template <class TInputImage, class TCoordRep>
+typename BCOInterpolateImageFunction< TInputImage, TCoordRep >
+::OutputType
+BCOInterpolateImageFunction<TInputImage, TCoordRep>
+::EvaluateAtContinuousIndex( const ContinuousIndexType & index ) const
+{
+ unsigned int dim; // index over dimension
+ unsigned int neighborNumber = (2*m_Radius+1) * (2*m_Radius+1);
+
+
+ /**
+ * Compute base index = closet index below point
+ * Compute distance from point to base index
+ */
+ IndexType baseIndex;
+ double distance[ImageDimension];
+
+ for( dim = 0; dim < ImageDimension; dim++ )
+ {
+ baseIndex[dim] = itk::Math::Floor< IndexValueType >( index[dim] );
+ }
+
+ /*
+ std::cout << "index : " << index << std::endl;
+ std::cout << "baseIndex : " << baseIndex << std::endl;
+ */
+
+ RealType value = 0;
+
+
+ IndexType neighIndex;
+ neighIndex.Fill(0);
+
+ vnl_vector<RealType> lineRes(2*m_Radius+1, 0.);
+
+ for( unsigned int counter = 0; counter < neighborNumber; counter++ )
+ {
+
+ unsigned int upper = counter; // each bit indicates upper/lower neighbour
+ double BCOCoefIndex; //link beetween counter and BCO Coeff.
+
+
+ // get neighbor index
+ for( dim = 0; dim < ImageDimension; dim++ )
+ {
+ if ( upper & 1 )
+ {
+ neighIndex[dim] = baseIndex[dim] + 1;
+#ifdef ITK_USE_CENTERED_PIXEL_COORDINATES_CONSISTENTLY
+ // Take care of the case where the pixel is just
+ // in the outer upper boundary of the image grid.
+ if( neighIndex[dim] > this->m_EndIndex[dim] )
+ {
+ neighIndex[dim] = this->m_EndIndex[dim];
+ }
+#endif
+ }
+ else
+ {
+ neighIndex[dim] = baseIndex[dim];
+#ifdef ITK_USE_CENTERED_PIXEL_COORDINATES_CONSISTENTLY
+ // Take care of the case where the pixel is just
+ // in the outer lower boundary of the image grid.
+ if( neighIndex[dim] < this->m_StartIndex[dim] )
+ {
+ neighIndex[dim] = this->m_StartIndex[dim];
+ }
+#endif
+ }
+
+ upper >>= 1;
+ }
+ // Proceed column
+ BCOCoefIndex = static_cast<int>(counter - floor(counter/(2*m_Radius+1))*(2*m_Radius+1));
+ lineRes[BCOCoefIndex] = lineRes[BCOCoefIndex]
+ + static_cast<RealType>( this->GetInputImage()->GetPixel( neighIndex ) ) * m_BCOCoef(BCOCoefIndex) ;
+ }
+
+ //Proceed line
+ for(unsigned int l=0; l<lineRes.size(); l++)
+ {
+ value = value + lineRes[l]*m_BCOCoef(l);
+ }
+
+ std::cout << "Fin " << value / (2*m_Radius * 2*m_Radius) << std::endl;
+ return ( static_cast<OutputType>( value / (2*m_Radius * 2*m_Radius) ) );
+}
+
+
+
+template <class TInputImage, class TCoordRep>
+void BCOInterpolateImageFunction<TInputImage, TCoordRep>
+::SetRadius(unsigned int radius)
+{
+ if (radius == 0)
+ {
+ itkExceptionMacro(<< "Radius Must Be Strictly Superior to 0");
+ }
+ else
+ {
+ m_Radius = radius;
+ }
+}
+
+template <class TInputImage, class TCoordRep>
+unsigned int BCOInterpolateImageFunction<TInputImage, TCoordRep>
+::GetRadius() const
+{
+ return m_Radius;
+}
+
+template <class TInputImage, class TCoordRep>
+void BCOInterpolateImageFunction<TInputImage, TCoordRep>
+::SetAlpha(double alpha)
+{
+ m_Alpha = alpha;
+}
+
+template <class TInputImage, class TCoordRep>
+double BCOInterpolateImageFunction<TInputImage, TCoordRep>
+::GetAlpha() const
+{
+ return m_Alpha;
+}
+
+template<class TInputImage, class TCoordRep>
+void
+BCOInterpolateImageFunction<TInputImage, TCoordRep>
+::EvaluateCoef()
+{
+ // Init BCO coefficient container
+ unsigned int size = 2*m_Radius+1;
+ m_BCOCoef = CoefContainerType(size, 0.);
+ double val = 0.;
+ //double step = 1./static_cast<double>(m_Radius);
+ double step = 1./static_cast<double>(size);
+
+
+ // Compute BCO coefficients
+ for (unsigned int i = 0; i <= m_Radius; i++)
+ {
+ // Compute the distance according to alpha.
+ double d = val;
+ if( d <= 2. )
+ {
+ if (d <= 1.)
+ {
+ m_BCOCoef[m_Radius+i] = (m_Alpha + 2.)*vcl_abs(vcl_pow(d, 3)) - (m_Alpha + 3.)*vcl_pow(d, 2) + 1;
+ }
+ else
+ {
+ m_BCOCoef[m_Radius+i] = m_Alpha*vcl_abs(vcl_pow(d, 3)) - 5*m_Alpha*vcl_pow(d, 2) + 8*m_Alpha*vcl_abs(d) - 4*m_Alpha;
+ }
+ }
+ else
+ {
+ m_BCOCoef[m_Radius+i] = 0;
+ }
+
+ if(i>0)
+ {
+ // the filter is symtric
+ m_BCOCoef[m_Radius-i] = m_BCOCoef[m_Radius+i];
+ }
+ val = val + step;
+ }
+ std::cout << m_BCOCoef << std::endl;
+
+ }
+
+
+} //namespace otb
+
+#endif
diff --git a/Testing/Code/BasicFilters/CMakeLists.txt b/Testing/Code/BasicFilters/CMakeLists.txt
index 568e3be66c5effc42df5d96109321ae2d16323f4..339410a01fc696436144b8f6867e0666e4315db0 100644
--- a/Testing/Code/BasicFilters/CMakeLists.txt
+++ b/Testing/Code/BasicFilters/CMakeLists.txt
@@ -673,7 +673,7 @@ ADD_TEST(bfTvProlateInterpolateImageFunction ${BASICFILTERS_TESTS6}
${TEMP}/bfProlateInterpolateImageFunctionCosOutput.tif
${TEMP}/bfProlateInterpolateImageFunctionITKCosOutput.tif
${TEMP}/bfProlateInterpolateImageFunctionProOutput.tif
- 3 # radius
+ 50 # radius
0.5 0.5
127.33 44.9
259.67 21.43
@@ -733,6 +733,28 @@ ADD_TEST(bfTvWindowedSincInterpolateImageHammingFunction ${BASICFILTERS_TESTS6}
-1 -1
)
+# ------- otb::BCOInterpolateImageFunction ------------------------------
+
+ADD_TEST(bfTuBCOInterpolateImageFunctionNew ${BASICFILTERS_TESTS6}
+ otbBCOInterpolateImageFunctionNew)
+
+ADD_TEST(bfTvBCOInterpolateImageFunction ${BASICFILTERS_TESTS6}
+ otbBCOInterpolateImageFunction
+ ${INPUTDATA}/poupees.tif
+ ${TEMP}/bfBCOInterpolateImageFunctionOutput.txt
+ 50 # radius
+ -0.5 # optimised bicubic
+ 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
+ )
+
+
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbBasicFiltersTests7 ~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1797,6 +1819,7 @@ otbWindowedSincInterpolateImageGaussianFunctionNew.cxx
otbWindowedSincInterpolateImageGaussianFunction.cxx
otbWindowedSincInterpolateImageHammingFunctionNew.cxx
otbWindowedSincInterpolateImageHammingFunction.cxx
+otbBCOInterpolateImageFunction.cxx
)
# A enrichir
diff --git a/Testing/Code/BasicFilters/otbBCOInterpolateImageFunction.cxx b/Testing/Code/BasicFilters/otbBCOInterpolateImageFunction.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a1ffa9cc88c202c4f476e1b55062827769caa010
--- /dev/null
+++ b/Testing/Code/BasicFilters/otbBCOInterpolateImageFunction.cxx
@@ -0,0 +1,100 @@
+/*=========================================================================
+
+ 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 "otbBCOInterpolateImageFunction.h"
+#include "otbVectorImage.h"
+#include "otbImage.h"
+#include "otbImageFileReader.h"
+#include <fstream>
+
+int otbBCOInterpolateImageFunctionNew(int argc, char * argv[])
+{
+ typedef otb::Image<double, 2> ImageType;
+ typedef otb::BCOInterpolateImageFunction<ImageType, double> InterpolatorType;
+
+ // Instantiating object
+ InterpolatorType::Pointer filter = InterpolatorType::New();
+
+ return EXIT_SUCCESS;
+}
+
+int otbBCOInterpolateImageFunction(int argc, char * argv[])
+{
+ const char * infname = argv[1];
+ const char * outfname = argv[2];
+ const unsigned int radius = atoi(argv[3]);
+ const double alpha = atof(argv[4]);
+
+ //typedef otb::VectorImage<double, 2> ImageType;
+
+ typedef otb::Image<double, 2> ImageType;
+ typedef otb::BCOInterpolateImageFunction<ImageType, double> InterpolatorType;
+ typedef InterpolatorType::ContinuousIndexType ContinuousIndexType;
+ typedef otb::ImageFileReader<ImageType> ReaderType;
+
+ int i = 5;
+
+ 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 filter = InterpolatorType::New();
+
+ ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(infname);
+ reader->Update();//TODO check if necessary
+
+ std::cout << "Alpha Checking : " << std::endl << filter->GetAlpha() << std::endl;
+ filter->SetAlpha(-1.0);
+ std::cout << filter->GetAlpha() << std::endl;
+ filter->SetAlpha(alpha);
+ std::cout << filter->GetAlpha() << std::endl;
+ std::cout << "Radius Checking : " << std::endl << filter->GetRadius() << std::endl;
+ filter->SetRadius(radius);
+ std::cout << filter->GetRadius() << std::endl;
+ std::cout << "Coef Checking : " << std::endl;
+ filter->EvaluateCoef();
+ std::cout << "Coef Generated ! " << std::endl;
+
+ filter->SetInputImage(reader->GetOutput());
+
+ std::ofstream file;
+ file.open(outfname);
+
+ for (std::vector<ContinuousIndexType>::iterator it = indicesList.begin(); it != indicesList.end(); ++it)
+ {
+ file << (*it) << " -> " << filter->EvaluateAtContinuousIndex((*it)) << std::endl;
+ }
+
+ file.close();
+
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/BasicFilters/otbBasicFiltersTests6.cxx b/Testing/Code/BasicFilters/otbBasicFiltersTests6.cxx
index bd0cd5acc36215ac2925d8b6b7afab2d9370b2c4..3540b7dfd630bee59b4373e7cd25b0787fcfd79a 100644
--- a/Testing/Code/BasicFilters/otbBasicFiltersTests6.cxx
+++ b/Testing/Code/BasicFilters/otbBasicFiltersTests6.cxx
@@ -37,4 +37,6 @@ void RegisterTests()
REGISTER_TEST(otbWindowedSincInterpolateImageGaussianFunction);
REGISTER_TEST(otbWindowedSincInterpolateImageHammingFunctionNew);
REGISTER_TEST(otbWindowedSincInterpolateImageHammingFunction);
+ REGISTER_TEST(otbBCOInterpolateImageFunctionNew);
+ REGISTER_TEST(otbBCOInterpolateImageFunction);
}