REFAC: KullbackLeibler based change detection

Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.h

@@ -58,19 +58,24 @@ public :
   {
     return this->fKurtosis;
   }
+  inline bool IsDataAvailable () const {
+    return this->fDataAvailable;
+  }
 
 protected :
 
   /** Moment estimation from intial neighborhood */
-  int  MakeSumAndMoments  ( const TInput & input );
+  void  MakeSumAndMoments  ( const TInput & input );
   /** Moment estimation from raw data */
-  int  MakeSumAndMoments  ( const itk::Image< typename TInput::ImageType::PixelType, 1 > * input );
+  void  MakeSumAndMoments  ( const itk::Image< typename TInput::ImageType::PixelType, 1 > * input );
   /** transformation moment -> cumulants (for Edgeworth) */
-  int MakeCumulants();
+  void MakeCumulants();
 
   double  fSum0, fSum1, fSum2, fSum3, fSum4;
   double  fMu1, fMu2, fMu3, fMu4;
   double  fMean, fVariance, fSkewness, fKurtosis;
+
+  bool  fDataAvailable;
 };
 
 
@@ -90,7 +95,13 @@ public :
   TOutput operator () ( const TInput1 & it1, const TInput2 & it2 )
   {
     CumulantsForEdgeworth<TInput1> cum1 ( it1 );
+    if ( !cum1.IsDataAvailable() )
+      return static_cast<TOutput>( 0. );
+
     CumulantsForEdgeworth<TInput2> cum2 ( it2 );
+    if ( !cum2.IsDataAvailable() )
+      return static_cast<TOutput>( 0. );
+
     return static_cast<TOutput> ( cum1.Divergence( cum2 )
                                   + cum2.Divergence( cum1 ) );
   }

Code/ChangeDetection/otbKullbackLeiblerDistanceImageFilter.txx

@@ -109,7 +109,7 @@ CumulantsForEdgeworth< TInput >
 /* ========== Moment estimation from initial neighborhood ============ */
 
 template < class TInput >
-int
+void
 CumulantsForEdgeworth< TInput >
 ::MakeSumAndMoments ( const TInput & input )
 {
@@ -134,10 +134,11 @@ CumulantsForEdgeworth< TInput >
 
   if ( fMu2 <= 0.0 )
   {
-    otbGenericMsgDebugMacro( << "Potential NAN detected in function MakeSumAndMoments.");
+    //otbGenericMsgDebugMacro( << "Potential NAN detected in function MakeSumAndMoments.");
     fMu3 = 0.0;
     fMu4 = 4.0;
-    return 1;
+    fDataAvailable = false;
+    //return 1;
   }
 
   double sigma = sqrt( fMu2 );
@@ -161,13 +162,13 @@ CumulantsForEdgeworth< TInput >
   fMu3 /= fSum0;
   fMu4 /= fSum0;
 
-  return 0;
+  fDataAvailable = true;
 }
 
 /* ================= Moment estimation from raw data ================= */
 
 template < class TInput >
-int
+void
 CumulantsForEdgeworth< TInput >
 ::MakeSumAndMoments ( const itk::Image< typename TInput::ImageType::PixelType, 1 > * input )
 {
@@ -202,10 +203,11 @@ CumulantsForEdgeworth< TInput >
 
   if ( fMu2 <= 0.0 )
   {
-    otbGenericMsgDebugMacro( << "Potential NAN detected in function MakeSumAndMoments.");
+    //otbGenericMsgDebugMacro( << "Potential NAN detected in function MakeSumAndMoments.");
     fMu3 = 0.0;
     fMu4 = 4.0;
-    return 1;
+    //return 1;
+    fDataAvailable = false;
   }
 
   double sigma = sqrt( fMu2 );
@@ -238,22 +240,26 @@ CumulantsForEdgeworth< TInput >
 
   otbGenericMsgDebugMacro( << "Moments: " << fMu1 << ",\t" << fMu2 << ",\t" << fMu3 << ",\t" << fMu4 );
 
-  return 0;
+  fDataAvailable = true;
+  //return 0;
 }
 
 
 /* ================= moments -> cumulants transformation ============= */
 
 template < class TInput >
-int
+void
 CumulantsForEdgeworth< TInput >
 ::MakeCumulants ()
 {
-  fMean = fMu1;
-  fVariance = fMu2;
-  fSkewness = fMu3;
-  fKurtosis = fMu4 - 3.0;
-  return 0;
+  if ( fDataAvailable )
+  {
+    fMean = fMu1;
+    fVariance = fMu2;
+    fSkewness = fMu3;
+    fKurtosis = fMu4 - 3.0;
+  }
+  //return 0;
 }
 
 } // end of namespace otb

Code/ChangeDetection/otbKullbackLeiblerProfileImageFilter.h

@@ -57,26 +57,32 @@ public :
   {
     return this->fCum[i];
   }
+  // Check data availability
+  inline bool IsDataAvailable () const {
+    return this->fDataAvailable;
+  }
 
   // debug
   int m_debug;
 
 protected :
 
-  // Estimation des moments � partir de voisinnages emboit�s
+  // Momentum Estimation from encapsulated neighborhood
   int  MakeSumAndMoments  ( const TInput & input, std::vector< itk::Array2D<int> > & mask );
-  // Estimation des moments � partir de la petite taille de fenetre
+  // momentum estimation from the smaller window
   int InitSumAndMoments ( const TInput & input, itk::Array2D<int> & mask );
   //
   int ReInitSumAndMoments ( const TInput & input, itk::Array2D<int> & mask, int level );
-  // transformation moment -> cumulants (pour Edgeworth)
+  // transformation moment -> cumulants (for Edgeworth)
   int MakeCumulants();
 
-  // Attributs internes � la classe
+  // Internal variables
   double  fSum0, fSum1, fSum2, fSum3, fSum4;
   CumulantSet fMu;
   CumulantSet fCum;
 
+  bool fDataAvailable;
+
 private :
   CumulantsForEdgeworthProfile (); // Not implemented
   CumulantsForEdgeworthProfile ( const TInput & input ); // Not implemented

Code/ChangeDetection/otbKullbackLeiblerProfileImageFilter.txx

@@ -181,7 +181,10 @@ CumulantsForEdgeworthProfile<TInput>
     }
   }
   if ( fSum0 == 0.0 )
+  {
+    fDataAvailable = false;
     return 1;
+  }
 
   double mu1, mu2;
 
@@ -189,7 +192,10 @@ CumulantsForEdgeworthProfile<TInput>
   mu2 = fSum2 / fSum0 - mu1 * mu1;
 
   if ( mu2 == 0.0 )
+  {
+    fDataAvailable = false;
     return 1;
+  }
 
   double sigma = sqrt( mu2 );
 
@@ -224,13 +230,18 @@ CumulantsForEdgeworthProfile<TInput>
   mu4 /= fSum0;
 
   if ( vnl_math_isnan( mu3 ) || vnl_math_isnan( mu4 ) )
+  {
+    fDataAvailable = false;
     return 1;
+  }
 
   fMu[0][0] = mu1;
   fMu[0][1] = mu2;
   fMu[0][2] = mu3;
   fMu[0][3] = mu4;
 
+  fDataAvailable = true;
+
   return 0;
 }
 
@@ -245,10 +256,10 @@ CumulantsForEdgeworthProfile<TInput>
   fMu[level].Fill(0.0);
   // mise a jour du comptage...
   double sum0 = 0.0,
-                sum1 = 0.0,
-                       sum2 = 0.0,
-                              sum3 = 0.0,
-                                     sum4 = 0.0;
+          sum1 = 0.0,
+          sum2 = 0.0,
+          sum3 = 0.0,
+          sum4 = 0.0;
 
   double pixel,pixel_2;
 
@@ -310,6 +321,9 @@ int
 CumulantsForEdgeworthProfile<TInput>
 ::MakeCumulants()
 {
+  if ( !IsDataAvailable() )
+    return 1;
+
   fCum.resize( fMu.size() );
   fCum = fMu;
 
@@ -416,7 +430,6 @@ KullbackLeiblerProfile<TInput1,TInput2,TOutput>
 ( const TInput1 & it1, const TInput2 & it2 )
 {
   CumulantsForEdgeworthProfile<TInput1> cum1 ( it1, m_mask );
-  CumulantsForEdgeworthProfile<TInput2> cum2 ( it2, m_mask );
 
   if ( cum1.m_debug )
   {
@@ -425,6 +438,8 @@ KullbackLeiblerProfile<TInput1,TInput2,TOutput>
     return static_cast<TOutput> ( resu );
   }
 
+  CumulantsForEdgeworthProfile<TInput2> cum2 ( it2, m_mask );
+
   if ( cum2.m_debug )
   {
     itk::VariableLengthVector<double> resu ( m_RadiusMax - m_RadiusMin + 1 );

Code/ChangeDetection/otbKullbackLeiblerSupervizedDistanceImageFilter.txx

@@ -72,6 +72,11 @@ KullbackLeiblerSupervizedDistance< TInput1, TInput2, TInputROIImage, TOutput >
 
   m_CumROI1 = new CumulantsForEdgeworth< ROIInputType1 > ( convertion1->GetOutput() );
 
+  if ( !m_CumROI1->IsDataAvailable() )
+  {
+    throw itk::ExceptionObject( __FILE__, __LINE__, "Cumulants estimated from ROI in image 1 are not usable", ITK_LOCATION );
+  }
+
   typedef ROIdataConversion< typename TInput2::ImageType, TInputROIImage >
   ROIConversionType2;
 
@@ -87,6 +92,11 @@ KullbackLeiblerSupervizedDistance< TInput1, TInput2, TInputROIImage, TOutput >
     delete m_CumROI2;
 
   m_CumROI2 = new CumulantsForEdgeworth< ROIInputType2 > ( convertion2->GetOutput() );
+
+  if ( !m_CumROI2->IsDataAvailable() )
+  {
+    throw itk::ExceptionObject( __FILE__, __LINE__, "Cumulants estimated from ROI in image 2 are not usable", ITK_LOCATION );
+  }
 }
 
 template < class TInput1, class TInput2, class TInputROIImage, class TOutput >
@@ -95,7 +105,19 @@ KullbackLeiblerSupervizedDistance< TInput1, TInput2, TInputROIImage, TOutput >
 ::operator () ( const TInput1 & it1, const TInput2 & it2 )
 {
   CumulantsForEdgeworth<TInput1> cum1 ( it1 );
+
+  if ( !cum1.IsDataAvailable() )
+  {
+    return static_cast<TOutput> ( 0. );
+  }
+
   CumulantsForEdgeworth<TInput2> cum2 ( it2 );
+
+  if ( !cum2.IsDataAvailable() )
+  {
+    return static_cast<TOutput> ( 0. );
+  }
+
   return static_cast<TOutput> ( m_CumROI1->Divergence( cum1 )
                                 + m_CumROI2->Divergence( cum2 ) );
 }