Detection de changement par distance de Kullback-leibler.

Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.h

+/*=========================================================================
+
+  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
+
+

Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.txx

+/*=========================================================================
+
+  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
+
+

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);

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
 )
 
 

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);
 }

Testing/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.cxx

+/*=========================================================================
+
+  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;
+}

Testing/Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilterNew.cxx

+/*=========================================================================
+
+  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;
+}