diff --git a/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.h b/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..486691e62315ba239ba4278a2bf4e9b3398c4fbd
--- /dev/null
+++ b/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.h
@@ -0,0 +1,148 @@
+/*=========================================================================
+
+ 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 __otbKullbackLeiblerDistanceImageFilter_h
+#define __otbKullbackLeiblerDistanceImageFilter_h
+
+#include "itkVariableLengthVector.h"
+#include "otbBinaryFunctorNeighborhoodImageFilter.h"
+
+namespace otb {
+
+template <class TInput>
+class CumulantsForEdgeworth
+{
+ public :
+ CumulantsForEdgeworth ( const TInput & input );
+ virtual ~CumulantsForEdgeworth () { };
+
+ // Divergence de KL
+ template <class TInput2>
+ double KL_divergence ( const CumulantsForEdgeworth<TInput2> & cumulants );
+
+ // access to attributes
+ inline double GetMean () const { return this->fMean; }
+ inline double GetVariance () const { return this->fVariance; }
+ inline double GetSkewness () const { return this->fSkewness; }
+ inline double GetKurtosis () const { return this->fKurtosis; }
+
+ protected :
+
+ // estimation des moments à partir du voisinage initial
+ int MakeSumAndMoments ( const TInput & input );
+ // transformation moment -> cumulants (pour Edgeworth)
+ int MakeCumulants();
+
+ // Attributs internes à la classe
+ double fSum0, fSum1, fSum2, fSum3, fSum4;
+ double fMu1, fMu2, fMu3, fMu4;
+ double fMean, fVariance, fSkewness, fKurtosis;
+};
+
+
+
+namespace Functor {
+
+ template < class TInput1, class TInput2, class TOutput >
+ class KullbackLeiblerDistance
+ {
+ public :
+ KullbackLeiblerDistance () { }
+ virtual ~KullbackLeiblerDistance () { }
+ TOutput operator () ( const TInput1 & it1, const TInput2 & it2 ) {
+ CumulantsForEdgeworth<TInput1> cum1 ( it1 );
+ CumulantsForEdgeworth<TInput2> cum2 ( it2 );
+ return static_cast<TOutput> ( cum1.KL_divergence( cum2 )
+ + cum2.KL_divergence( cum1 ) );
+ }
+ };
+
+} // Functor
+
+/** \class KullbackLeiblerDistanceImageFilter
+ * \brief Implements neighborhood-wise the computation of KullbackLeibler distance over Edgeworth approximation.
+ *
+ * This filter is parametrized over the types of the two
+ * input images and the type of the output image.
+ *
+ * Numeric conversions (castings) are done by the C++ defaults.
+ *
+ * The filter will walk over all the pixels in the two input images, and for
+ * each one of them it will do the following:
+ *
+ * - cast the input 1 pixel value to \c double
+ * - cast the input 2 pixel value to \c double
+ * - compute the first four cumulants of the two pixel values
+ * - compute the value of the Edgeorth approximation of the KL distance
+ * - cast the \c double value resulting to the pixel type of the output image
+ * - store the casted value into the output image.
+ *
+ * The filter expect all images to have the same dimension
+ * (e.g. all 2D, or all 3D, or all ND)
+ *
+ * See article of Lin Saito et Levine
+ * "Edgeworth Approximation of the Kullback-Leibler Distance Towards Problems in Image Analysis"
+ * and
+ * "Edgeworth Expansions of the Kullback-Leibler Information" (submitted to JASA, nov 25, 1999)
+ * http://www.math.ucdavis.edu/~saito/publications
+ *
+ * \ingroup IntensityImageFilters Multithreaded
+ */
+template <class TInputImage1, class TInputImage2, class TOutputImage>
+class ITK_EXPORT KullbackLeiblerDistanceImageFilter :
+ public otb::BinaryFunctorNeighborhoodImageFilter<
+ TInputImage1,TInputImage2,TOutputImage,
+ Functor::KullbackLeiblerDistance<
+ typename itk::ConstNeighborhoodIterator<TInputImage1>,
+ typename itk::ConstNeighborhoodIterator<TInputImage2>,
+ typename TOutputImage::PixelType> >
+{
+ public:
+ /** Standard class typedefs. */
+ typedef KullbackLeiblerDistanceImageFilter Self;
+ typedef typename otb::BinaryFunctorNeighborhoodImageFilter<
+ TInputImage1,TInputImage2,TOutputImage,
+ Functor::KullbackLeiblerDistance<
+ typename itk::ConstNeighborhoodIterator<TInputImage1>,
+ typename itk::ConstNeighborhoodIterator<TInputImage2>,
+ typename TOutputImage::PixelType>
+ > Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ protected:
+ KullbackLeiblerDistanceImageFilter() {}
+ virtual ~KullbackLeiblerDistanceImageFilter() {}
+
+ private:
+ KullbackLeiblerDistanceImageFilter(const Self&); //purposely not implemented
+ void operator=(const Self&); //purposely not implemented
+
+};
+
+#ifndef ITK_MANUAL_INSTANTIATION
+#include "otbKullbackLeiblerDistanceImageFilter.txx"
+#endif
+
+} // namespace otb
+
+#endif
+
+
diff --git a/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.txx b/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..7ec38194c7aec9c6cb8ad15ac482c3d5593635d0
--- /dev/null
+++ b/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.txx
@@ -0,0 +1,168 @@
+/*=========================================================================
+
+ 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 __otbKullbackLeiblerDistanceImageFilter_txx
+#define __otbKullbackLeiblerDistanceImageFilter_txx
+
+#include "otbKullbackLeiblerDistanceImageFilter.h"
+
+
+/* *******************************************************************
+ * Classe CumulantsForEdgeworth
+ *********************************************************************
+ */
+template <class TInput>
+CumulantsForEdgeworth<TInput>
+::CumulantsForEdgeworth ( const TInput & input )
+{
+ MakeSumAndMoments( input );
+ MakeCumulants();
+}
+
+/* ========================== Divergence de KL ======================= */
+
+template <class TInput>
+template <class TInput2>
+double
+CumulantsForEdgeworth<TInput>
+::KL_divergence ( const CumulantsForEdgeworth<TInput2> & cumulants )
+{
+ double cum1 = GetMean();
+ double cum2 = GetVariance();
+ double cum3 = GetSkewness();
+
+ double sqrt_cum2 = sqrt( cum2 );
+ double cum2_3 = cum2 * cum2 * cum2;
+ double cum3_2 = cum3 * cum3;
+
+ double tilde_cum1 = cumulants.GetMean();
+ double tilde_cum2 = cumulants.GetVariance();
+ double tilde_cum3 = cumulants.GetSkewness();
+ double tilde_cum4 = cumulants.GetKurtosis();
+
+ double tilde_cum2_2 = cum2 * cum2;
+ double tilde_cum2_3 = cum2 * tilde_cum2_2;
+ double tilde_cum2_6 = tilde_cum2_3 * tilde_cum2_3;
+ double tilde_cum3_2 = tilde_cum3 * tilde_cum3;
+
+ double beta = sqrt_cum2 / tilde_cum2;
+ double alpha = ( cum1 - tilde_cum1 ) / tilde_cum2;
+
+ double alpha_2 = alpha * alpha;
+ double alpha_4 = alpha_2 * alpha_2;
+ double alpha_6 = alpha_2 * alpha_4;
+
+ double beta_2 = beta * beta;
+ double beta_4 = beta_2 * beta_2;
+ double beta_6 = beta_2 * beta_4;
+
+ double c2 = alpha_2 + beta_2;
+ double c3 = alpha * ( alpha_2 + 3.0 * beta_2 );
+ double c4 = alpha_4 + 6.0 * alpha_2 * beta_2 + 3.0 * beta_4;
+ double c6 = alpha_6 + 15.0 * alpha_4 * beta_2 + 45.0 * alpha_2 * beta_4 + 15.0 * beta_6;
+
+ double a1 = c3 - 3.0 * alpha / tilde_cum2;
+ double a2 = c4 - 6.0 * c2 / tilde_cum2 + 3.0 / tilde_cum2_2; // Watch for tilde_cum2_2 mistake in the article!
+ double a3 = c6 - 15.0 * c4 / tilde_cum2 + 45.0 * c2 / tilde_cum2_2 - 15.0 / tilde_cum2_3;
+
+ double tmp = cum1 - tilde_cum1 + sqrt_cum2;
+ double resu = cum3_2 / ( 12.0 * cum2_3 )
+ + ( log( tilde_cum2 / cum2 )
+ - 1.0 + tmp * tmp / tilde_cum2 ) / 2.0
+ - ( tilde_cum3 * a1 / 6.0 + tilde_cum4 * a2 / 24.0 + tilde_cum3_2 * a3 / 72.0 )
+ - tilde_cum3_2 * ( c6 - 6.0 * c4 / cum2 + 9.0 * c2 / tilde_cum2_2 ) / 72.0
+ - 10.0 * cum3 * tilde_cum3 * ( cum1 - tilde_cum1 ) * ( cum2 - tilde_cum2 ) / tilde_cum2_6;
+
+ return resu < 0.0 ? 0.0 : resu;
+}
+
+/* ====== Estimation des moments à partir du voisinage initial ======= */
+
+template <class TInput>
+int
+CumulantsForEdgeworth<TInput>
+::MakeSumAndMoments ( const TInput & input )
+{
+
+ fSum0 = fSum1 = fSum2 = fSum3 = fSum4 = 0.0;
+ double pixel,pixel_2;
+
+ for ( unsigned long i = 0; i < input.Size(); i++ )
+ {
+ pixel = static_cast<double> ( input.GetPixel(i) );
+ pixel_2 = pixel * pixel;
+
+ fSum0 += 1.0;
+ fSum1 += pixel;
+ fSum2 += pixel_2;
+ fSum3 += pixel_2 * pixel;
+ fSum4 += pixel_2 * pixel_2;
+ }
+
+ fMu1 = fSum1 / fSum0;
+ fMu2 = fSum2 / fSum0 - fMu1 * fMu1;
+
+ if ( fMu2 <= 0.0 )
+ {
+ std::cerr << "Potential NAN detected in " << __PRETTY_FUNCTION__ << "\n";
+ fMu3 = 0.0;
+ fMu4 = 4.0;
+ return 1;
+ }
+
+ double sigma = sqrt( fMu2 );
+
+ itk::VariableLengthVector<double> tab ( input.Size() );
+ double * x = const_cast<double *> ( tab.GetDataPointer() );
+ for ( unsigned long i = 0; i < input.Size(); i++ )
+ *x++ = ( static_cast<double> ( input.GetPixel(i) ) - fMu1 ) / sigma;
+
+ fMu3 = fMu4 = 0.0;
+ x = const_cast<double *> ( tab.GetDataPointer() );
+ for ( unsigned long i = 0; i < input.Size(); i++ )
+ {
+ pixel = *x++;
+ pixel_2 = pixel * pixel;
+
+ fMu3 += pixel * pixel_2;
+ fMu4 += pixel_2 * pixel_2;
+ }
+
+ fMu3 /= fSum0;
+ fMu4 /= fSum0;
+
+ return 0;
+}
+
+/* =========== transformation moment -> cumulants ==================== */
+
+template <class TInput>
+int
+CumulantsForEdgeworth<TInput>
+::MakeCumulants ()
+{
+ fMean = fMu1;
+ fVariance = fMu2;
+ fSkewness = fMu3;
+ fKurtosis = fMu4 - 3.0;
+ return 0;
+}
+
+
+#endif
+
+
diff --git a/Code/Common/otbBinaryFunctorNeighborhoodImageFilter.txx b/Code/Common/otbBinaryFunctorNeighborhoodImageFilter.txx
index 59842a04b02b4492dadaaa4f757959f5def1e54b..176d3475fcc7558c2bbd959dfe0d7b8657dce85e 100755
--- a/Code/Common/otbBinaryFunctorNeighborhoodImageFilter.txx
+++ b/Code/Common/otbBinaryFunctorNeighborhoodImageFilter.txx
@@ -125,7 +125,7 @@ BinaryFunctorNeighborhoodImageFilter<TInputImage1, TInputImage2, TOutputImage, T
for (fit1=faceList1.begin(), fit2=faceList2.begin(); fit1 != faceList1.end(), fit2 != faceList2.end(); ++fit1, ++fit2)
{
neighInputIt1 = itk::ConstNeighborhoodIterator<TInputImage1>(r1, inputPtr1, *fit1);
- neighInputIt2 = itk::ConstNeighborhoodIterator<TInputImage2>(r1, inputPtr2, *fit2);
+ neighInputIt2 = itk::ConstNeighborhoodIterator<TInputImage2>(r2, inputPtr2, *fit2);
// outputIt = itk::ImageRegionIterator<TOutputImage>(outputPtr, outputRegionForThread);
outputIt = itk::ImageRegionIterator<TOutputImage>(outputPtr, *fit1);
diff --git a/Testing/Code/ChangeDetection/CMakeLists.txt b/Testing/Code/ChangeDetection/CMakeLists.txt
index cdb4b9e3c348d6961a25cb78ee90ee1a11561014..2ffdae80376af237dd34e53f059e2ada0e2ed0bb 100755
--- a/Testing/Code/ChangeDetection/CMakeLists.txt
+++ b/Testing/Code/ChangeDetection/CMakeLists.txt
@@ -69,8 +69,9 @@ ADD_TEST(cdTvLHMI ${CHANGEDETECTION_TESTS}
)
ADD_TEST(cdTvJHMI ${CHANGEDETECTION_TESTS}
- --compare-image ${TOL} ${BASELINE}/cdJHMIImage.png
- ${TEMP}/cdJHMIImage.png
+ --compare-image ${TOL}
+ ${BASELINE}/cdJHMIImage.png
+ ${TEMP}/cdJHMIImage.png
otbJHMIChangeDetectionTest
${INPUTDATA}/GomaAvant.png
${INPUTDATA}/GomaApres.png
@@ -78,6 +79,18 @@ ADD_TEST(cdTvJHMI ${CHANGEDETECTION_TESTS}
${TEMP}/cdJHMIImage.png
)
+ADD_TEST(cdTuKullbackLeiblerDistanceImageFilterNew ${CHANGEDETECTION_TESTS}
+ otbKullbackLeiblerDistanceImageFilterNew)
+
+ADD_TEST(cdTvKullbackLeiblerDistanceImageFilter ${CHANGEDETECTION_TESTS}
+ --compare-image ${EPSILON}
+ ${BASELINE}/cdTVKullbackLeiblerDistanceImageFilterOutput.hdr
+ ${TEMP}/cdTVKullbackLeiblerDistanceImageFilterOutput.hdr
+ otbKullbackLeiblerDistanceImageFilter
+ ${INPUTDATA}/GomaAvant.png
+ ${INPUTDATA}/GomaApres.png
+ ${TEMP}/cdTVKullbackLeiblerDistanceImageFilterOutput.hdr
+ 35)
# ------- Fichiers sources CXX -----------------------------------
@@ -88,6 +101,8 @@ otbMeanDiffChangeDetectionTest.cxx
otbMeanRatioChangeDetectionTest.cxx
otbLHMIChangeDetectionTest.cxx
otbJHMIChangeDetectionTest.cxx
+otbKullbackLeiblerDistanceImageFilterNew.cxx
+otbKullbackLeiblerDistanceImageFilter.cxx
)
diff --git a/Testing/Code/ChangeDetection/otbChangeDetectionTests.cxx b/Testing/Code/ChangeDetection/otbChangeDetectionTests.cxx
index 1fc6c9e7a1a47f57fc42c3023c2b05305f7c134f..ec33ff32046f2bd5f7c76fd148883fde6af1b9cc 100755
--- a/Testing/Code/ChangeDetection/otbChangeDetectionTests.cxx
+++ b/Testing/Code/ChangeDetection/otbChangeDetectionTests.cxx
@@ -32,4 +32,6 @@ REGISTER_TEST(otbMeanDiffChangeDetectionTest );
REGISTER_TEST(otbMeanRatioChangeDetectionTest );
REGISTER_TEST(otbLHMIChangeDetectionTest );
REGISTER_TEST(otbJHMIChangeDetectionTest );
+REGISTER_TEST(otbKullbackLeiblerDistanceImageFilterNew);
+REGISTER_TEST(otbKullbackLeiblerDistanceImageFilter);
}
diff --git a/Testing/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.cxx b/Testing/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..cb1b7e68c89b6b96f66f200a6c06c5add38a7104
--- /dev/null
+++ b/Testing/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.cxx
@@ -0,0 +1,81 @@
+/*=========================================================================
+
+ 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 "otbImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbKullbackLeiblerDistanceImageFilter.h"
+
+
+int otbKullbackLeiblerDistanceImageFilter(int argc, char * argv[])
+{
+ try
+ {
+ if(argc != 5)
+ {
+ std::cerr<<"Detection de changements par mesure de Kullback-Leibler, optimisee par un developpement de Edgeworth\n";
+ std::cerr << argv[0] << " imgAv imgAp imgResu winSize\n";
+ return 1;
+ }
+
+ char * fileName1 = argv[1];
+ char * fileName2 = argv[2];
+ char * fileNameOut = argv[3];
+ int winSize = atoi(argv[4]);
+
+
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+
+ typedef otb::Image<PixelType,Dimension> ImageType;
+ typedef otb::KullbackLeiblerDistanceImageFilter<ImageType,ImageType,ImageType> FilterType;
+ typedef otb::ImageFileReader<ImageType> ReaderType;
+ typedef otb::ImageFileWriter<ImageType> WriterType;
+
+ ReaderType::Pointer reader1 = ReaderType::New();
+ reader1->SetFileName( fileName1 );
+
+ ReaderType::Pointer reader2 = ReaderType::New();
+ reader2->SetFileName( fileName2 );
+
+ FilterType::Pointer filter = FilterType::New();
+ filter->SetRadius( (winSize-1)/2 );
+ filter->SetInput1( reader1->GetOutput() );
+ filter->SetInput2( reader2->GetOutput() );
+
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetFileName(fileNameOut);
+ writer->SetInput(filter->GetOutput());
+ writer->Update();
+
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject thrown !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ catch( ... )
+ {
+ std::cout << "Unknown exception thrown !" << std::endl;
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilterNew.cxx b/Testing/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilterNew.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7f40408688888d9723f1dc55335a1bc48942f4b9
--- /dev/null
+++ b/Testing/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilterNew.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 "otbImage.h"
+#include "otbKullbackLeiblerDistanceImageFilter.h"
+
+
+int otbKullbackLeiblerDistanceImageFilterNew(int argc, char * argv[])
+{
+ try
+ {
+ const unsigned int Dimension = 2;
+ typedef double PixelType;
+
+ typedef otb::Image<PixelType,Dimension> ImageType;
+ typedef otb::KullbackLeiblerDistanceImageFilter<ImageType,ImageType,ImageType> FilterType;
+
+ FilterType::Pointer filter = FilterType::New();
+
+ }
+
+ catch( itk::ExceptionObject & err )
+ {
+ std::cout << "Exception itk::ExceptionObject thrown !" << std::endl;
+ std::cout << err << std::endl;
+ return EXIT_FAILURE;
+ }
+
+ catch( ... )
+ {
+ std::cout << "Unknown exception thrown !" << std::endl;
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}