diff --git a/Code/ChangeDetection/otbCBAMI.h b/Code/ChangeDetection/otbCBAMI.h
new file mode 100644
index 0000000000000000000000000000000000000000..3980902a4251b8a2c91fce7f92447b2e11b8a327
--- /dev/null
+++ b/Code/ChangeDetection/otbCBAMI.h
@@ -0,0 +1,281 @@
+/*=========================================================================
+
+ 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 __otbCBAMI_h
+#define __otbCBAMI_h
+
+#include <vector>
+#include <stdlib.h>
+#include <math.h>
+
+
+namespace otb
+{
+
+// #define EPSILON_VALUE_CBAMI 0.01
+#define epsilon 0.01
+
+namespace Functor
+{
+
+template< class TInput1, class TInput2, class TOutput>
+class CBAMI
+{
+public:
+
+ typedef typename std::vector<TOutput> VectorType;
+ typedef typename VectorType::iterator IteratorType;
+ typedef typename std::vector<VectorType> VectorOfVectorType;
+ typedef typename VectorOfVectorType::iterator VecOfVecIteratorType;
+
+ CBAMI() {};
+ virtual ~CBAMI() {};
+ inline TOutput operator()( const TInput1 & itA,
+ const TInput2 & itB)
+ {
+ //const double epsilon = 0.01;
+ VectorType vecA;
+ VectorType vecB;
+
+ for (unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+
+ vecA.push_back(static_cast<double>(itA.GetPixel(pos)));
+ vecB.push_back(static_cast<double>(itB.GetPixel(pos)));
+
+ }
+
+ normalizeInPlace(vecA);
+ normalizeInPlace(vecB);
+
+ return static_cast<TOutput>( - vcl_log(PhiMI(vecA, vecB)+epsilon) );
+ }
+
+protected:
+
+ inline void normalizeInPlace(VectorType vx)
+ {
+
+ TOutput Ex = 0.0;
+
+ IteratorType itx;
+
+ for ( itx = vx.begin(); itx < vx.end(); ++itx)
+ {
+ Ex += static_cast<TOutput>(*itx);
+ }
+
+ Ex /= (vx.size());
+
+ TOutput Vx = 0.0;
+
+ for ( itx = vx.begin(); itx < vx.end(); ++itx)
+ {
+ Vx += static_cast<TOutput>(vcl_pow(static_cast<double>((*itx)-Ex),2));
+ }
+
+ Vx /= (vx.size());
+
+ for ( itx = vx.begin(); itx < vx.end(); ++itx)
+ {
+ (*itx) = ((*itx)-Ex)/vcl_sqrt(Vx);
+ }
+
+
+ }
+ inline TOutput Exyc(VectorType vx, VectorType vy)
+ {
+
+ TOutput Exy = 0.0;
+ TOutput Ex = 0.0;
+ TOutput Ey = 0.0;
+
+ IteratorType itx;
+ IteratorType ity;
+
+ for ( itx = vx.begin(), ity = vy.begin(); itx < vx.end(); ++itx, ++ity)
+ {
+ //Ex += (*itx);
+ //Ey += (*ity);
+ Exy += (*itx)*(*ity);
+
+ }
+
+ //Ex /= (vx.size());
+ //Ey /= (vy.size());
+ Exy /= (vx.size());
+
+ return Exy-Ex*Ey;
+ }
+
+ inline TOutput Exyztc(VectorType vx, VectorType vy, VectorType vz, VectorType vt)
+ {
+
+ TOutput Exyzt = 0.0;
+
+ TOutput Exyz = 0.0;
+ TOutput Exyt = 0.0;
+ TOutput Exzt = 0.0;
+ TOutput Eyzt = 0.0;
+
+ TOutput Exy = 0.0;
+ TOutput Exz = 0.0;
+ TOutput Ext = 0.0;
+ TOutput Eyz = 0.0;
+ TOutput Eyt = 0.0;
+ TOutput Ezt = 0.0;
+
+ TOutput Ex = 0.0;
+ TOutput Ey = 0.0;
+ TOutput Ez = 0.0;
+ TOutput Et = 0.0;
+
+
+ IteratorType itx;
+ IteratorType ity;
+ IteratorType itz;
+ IteratorType itt;
+
+ for ( itx = vx.begin(),
+ ity = vy.begin(),
+ itz = vz.begin(),
+ itt = vt.begin();
+ itx < vx.end();
+ ++itx,
+ ++ity,
+ itz++,
+ itt++)
+ {
+ //Ex += (*itx);
+ //Ey += (*ity);
+ //Ez += (*itz);
+ //Et += (*itt);
+
+ Exy += (*itx)*(*ity);
+ Exz += (*itx)*(*itz);
+ Ext += (*itx)*(*itt);
+ Eyz += (*ity)*(*itz);
+ Eyt += (*ity)*(*itt);
+ Ezt += (*itz)*(*itt);
+
+ Exyz += (*itx)*(*ity)*(*itz);
+ Exyt += (*itx)*(*ity)*(*itt);
+ Exzt += (*itx)*(*itz)*(*itt);
+ Eyzt += (*ity)*(*itz)*(*itt);
+
+ Exyzt += (*itx)*(*ity)*(*itz)*(*itt);
+
+ }
+
+ /*Ex /= (vx.size());
+ Ey /= (vx.size());
+ Ez /= (vx.size());
+ Et /= (vx.size()); */
+
+ Exy /= (vx.size());
+ Exz /= (vx.size());
+ Ext /= (vx.size());
+ Eyz /= (vx.size());
+ Eyt /= (vx.size());
+ Ezt /= (vx.size());
+
+ Exyz /= (vx.size());
+ Exyt /= (vx.size());
+ Exzt /= (vx.size());
+ Eyzt /= (vx.size());
+
+
+ TOutput result = Exyzt - Exyz*Et- Exyt*Ez- Exzt*Ey- Eyzt*Ex +
+ Exy*Ez*Et + Exz*Et*Ey + Ext*Ey*Ez + Eyz*Et*Ex + Eyt*Ex*Ez + Ezt*Ex*Ey -
+ 3*Ex*Ey*Ez*Et;
+
+ return result;
+ }
+
+ inline TOutput Rxy(VectorType va, VectorType vb)
+ {
+
+ return Exyc(va, vb);
+
+ }
+
+ inline TOutput Qxijkl(VectorType va, VectorType vb, VectorType vc, VectorType vd)
+ {
+// IteratorType ita;
+// IteratorType itb;
+// IteratorType itc;
+// IteratorType itd;
+
+
+ TOutput Eabcd_c = Exyztc(va, vb, vc, vd);
+
+
+ TOutput Eab_c = Exyc(va, vb);
+ TOutput Eac_c = Exyc(va, vc);
+ TOutput Ead_c = Exyc(va, vd);
+ TOutput Ecd_c = Exyc(vc, vd);
+ TOutput Ebd_c = Exyc(vb, vd);
+ TOutput Ebc_c = Exyc(vb, vc);
+
+ return Eabcd_c - Eab_c*Ecd_c - Eac_c*Ebd_c - Ead_c*Ebc_c;
+
+
+
+
+ }
+
+ inline TOutput PhiMI(VectorType v1, VectorType v2)
+ {
+
+
+ VectorOfVectorType donnees;
+ donnees.push_back(v1);
+ donnees.push_back(v2);
+
+ VecOfVecIteratorType iti;
+ VecOfVecIteratorType itj;
+ VecOfVecIteratorType itk;
+ VecOfVecIteratorType itl;
+
+ TOutput termeR = 0.0;
+ TOutput termeQ = 0.0;
+
+ for ( iti = donnees.begin(); iti < donnees.end(); ++iti )
+ for ( itj = donnees.begin(); itj < donnees.end(); ++itj )
+ {
+ if (iti != itj)
+ termeR += static_cast<TOutput>(vcl_pow(static_cast<double>(Rxy((*iti),(*itj))),2));
+
+ for ( itk = donnees.begin(); itk < donnees.end(); ++itk )
+ for ( itl = donnees.begin(); itl < donnees.end(); itl++ )
+ {
+ if ((iti != itj) || (iti != itk) || (iti != itl))
+ termeQ += static_cast<TOutput>(vcl_pow( static_cast<double>(Qxijkl((*iti),(*itj),(*itk),(*itl))),2));
+ }
+ }
+
+
+ return 1.0/4.0*termeR + 1.0/48.0*termeQ;
+
+ }
+
+};
+}
+
+}
+
+#endif
diff --git a/Code/ChangeDetection/otbCBAMIChangeDetector.h b/Code/ChangeDetection/otbCBAMIChangeDetector.h
index 339f1b79529ef3a7a76fd5878c2d23772e4c3bdf..e077237c092c8d901b25abc844dd37a4dc96894c 100644
--- a/Code/ChangeDetection/otbCBAMIChangeDetector.h
+++ b/Code/ChangeDetection/otbCBAMIChangeDetector.h
@@ -20,9 +20,7 @@
#include "otbBinaryFunctorNeighborhoodImageFilter.h"
-#include <vector>
-#include <stdlib.h>
-#include <math.h>
+#include "otbCBAMI.h"
namespace otb
@@ -53,255 +51,6 @@ namespace otb
* \ingroup IntensityImageFilters Multithreaded
*/
-// #define EPSILON_VALUE_CBAMI 0.01
-
-namespace Functor
-{
-
-template< class TInput1, class TInput2, class TOutput>
-class CBAMI
-{
-public:
-
- typedef typename std::vector<TOutput> VectorType;
- typedef typename VectorType::iterator IteratorType;
- typedef typename std::vector<VectorType> VectorOfVectorType;
- typedef typename VectorOfVectorType::iterator VecOfVecIteratorType;
-
- CBAMI() {};
- virtual ~CBAMI() {};
- inline TOutput operator()( const TInput1 & itA,
- const TInput2 & itB)
- {
- const double epsilon = 0.01;
- VectorType vecA;
- VectorType vecB;
-
- for (unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
-
- vecA.push_back(static_cast<TOutput>(itA.GetPixel(pos)));
- vecB.push_back(static_cast<TOutput>(itB.GetPixel(pos)));
-
- }
-
- normalizeInPlace(vecA);
- normalizeInPlace(vecB);
-
- return static_cast<TOutput>( - vcl_log(PhiMI(vecA, vecB)+epsilon) );
- }
-
-protected:
-
- inline void normalizeInPlace(VectorType vx)
- {
-
- TOutput Ex = 0.0;
-
- IteratorType itx;
-
- for ( itx = vx.begin(); itx < vx.end(); ++itx)
- {
- Ex += (*itx);
- }
-
- Ex /= (vx.size());
-
- TOutput Vx = 0.0;
-
- for ( itx = vx.begin(); itx < vx.end(); ++itx)
- {
- Vx += vcl_pow((*itx)-Ex,2);
- }
-
- Vx /= (vx.size());
-
- for ( itx = vx.begin(); itx < vx.end(); ++itx)
- {
- (*itx) = ((*itx)-Ex)/vcl_sqrt(Vx);
- }
-
-
- }
- inline TOutput Exyc(VectorType vx, VectorType vy)
- {
-
- TOutput Exy = 0.0;
- TOutput Ex = 0.0;
- TOutput Ey = 0.0;
-
- IteratorType itx;
- IteratorType ity;
-
- for ( itx = vx.begin(), ity = vy.begin(); itx < vx.end(); ++itx, ++ity)
- {
- //Ex += (*itx);
- //Ey += (*ity);
- Exy += (*itx)*(*ity);
-
- }
-
- //Ex /= (vx.size());
- //Ey /= (vy.size());
- Exy /= (vx.size());
-
- return Exy-Ex*Ey;
- }
-
- inline TOutput Exyztc(VectorType vx, VectorType vy, VectorType vz, VectorType vt)
- {
-
- TOutput Exyzt = 0.0;
-
- TOutput Exyz = 0.0;
- TOutput Exyt = 0.0;
- TOutput Exzt = 0.0;
- TOutput Eyzt = 0.0;
-
- TOutput Exy = 0.0;
- TOutput Exz = 0.0;
- TOutput Ext = 0.0;
- TOutput Eyz = 0.0;
- TOutput Eyt = 0.0;
- TOutput Ezt = 0.0;
-
- TOutput Ex = 0.0;
- TOutput Ey = 0.0;
- TOutput Ez = 0.0;
- TOutput Et = 0.0;
-
-
- IteratorType itx;
- IteratorType ity;
- IteratorType itz;
- IteratorType itt;
-
- for ( itx = vx.begin(),
- ity = vy.begin(),
- itz = vz.begin(),
- itt = vt.begin();
- itx < vx.end();
- ++itx,
- ++ity,
- itz++,
- itt++)
- {
- //Ex += (*itx);
- //Ey += (*ity);
- //Ez += (*itz);
- //Et += (*itt);
-
- Exy += (*itx)*(*ity);
- Exz += (*itx)*(*itz);
- Ext += (*itx)*(*itt);
- Eyz += (*ity)*(*itz);
- Eyt += (*ity)*(*itt);
- Ezt += (*itz)*(*itt);
-
- Exyz += (*itx)*(*ity)*(*itz);
- Exyt += (*itx)*(*ity)*(*itt);
- Exzt += (*itx)*(*itz)*(*itt);
- Eyzt += (*ity)*(*itz)*(*itt);
-
- Exyzt += (*itx)*(*ity)*(*itz)*(*itt);
-
- }
-
- /*Ex /= (vx.size());
- Ey /= (vx.size());
- Ez /= (vx.size());
- Et /= (vx.size()); */
-
- Exy /= (vx.size());
- Exz /= (vx.size());
- Ext /= (vx.size());
- Eyz /= (vx.size());
- Eyt /= (vx.size());
- Ezt /= (vx.size());
-
- Exyz /= (vx.size());
- Exyt /= (vx.size());
- Exzt /= (vx.size());
- Eyzt /= (vx.size());
-
-
- TOutput result = Exyzt - Exyz*Et- Exyt*Ez- Exzt*Ey- Eyzt*Ex +
- Exy*Ez*Et + Exz*Et*Ey + Ext*Ey*Ez + Eyz*Et*Ex + Eyt*Ex*Ez + Ezt*Ex*Ey -
- 3*Ex*Ey*Ez*Et;
-
- return result;
- }
-
- inline TOutput Rxy(VectorType va, VectorType vb)
- {
-
- return Exyc(va, vb);
-
- }
-
- inline TOutput Qxijkl(VectorType va, VectorType vb, VectorType vc, VectorType vd)
- {
-// IteratorType ita;
-// IteratorType itb;
-// IteratorType itc;
-// IteratorType itd;
-
-
- TOutput Eabcd_c = Exyztc(va, vb, vc, vd);
-
-
- TOutput Eab_c = Exyc(va, vb);
- TOutput Eac_c = Exyc(va, vc);
- TOutput Ead_c = Exyc(va, vd);
- TOutput Ecd_c = Exyc(vc, vd);
- TOutput Ebd_c = Exyc(vb, vd);
- TOutput Ebc_c = Exyc(vb, vc);
-
- return Eabcd_c - Eab_c*Ecd_c - Eac_c*Ebd_c - Ead_c*Ebc_c;
-
-
-
-
- }
-
- inline TOutput PhiMI(VectorType v1, VectorType v2)
- {
-
-
- VectorOfVectorType donnees;
- donnees.push_back(v1);
- donnees.push_back(v2);
-
- VecOfVecIteratorType iti;
- VecOfVecIteratorType itj;
- VecOfVecIteratorType itk;
- VecOfVecIteratorType itl;
-
- TOutput termeR = 0.0;
- TOutput termeQ = 0.0;
-
- for ( iti = donnees.begin(); iti < donnees.end(); ++iti )
- for ( itj = donnees.begin(); itj < donnees.end(); ++itj )
- {
- if (iti != itj)
- termeR += vcl_pow(Rxy((*iti),(*itj)),2);
-
- for ( itk = donnees.begin(); itk < donnees.end(); ++itk )
- for ( itl = donnees.begin(); itl < donnees.end(); itl++ )
- {
- if ((iti != itj) || (iti != itk) || (iti != itl))
- termeQ += vcl_pow( Qxijkl((*iti),(*itj),(*itk),(*itl)),2);
- }
- }
-
-
- return 1.0/4.0*termeR + 1.0/48.0*termeQ;
-
- }
-
-};
-}
-
template <class TInputImage1, class TInputImage2, class TOutputImage>
class ITK_EXPORT CBAMIChangeDetector :
public BinaryFunctorNeighborhoodImageFilter<
diff --git a/Code/ChangeDetection/otbCorrelationChangeDetector.h b/Code/ChangeDetection/otbCorrelationChangeDetector.h
index b2b252938ba637fbd6a99fde329035d24e5b518e..147e51849e931c22a67487bcdd0bc07b9c120c14 100644
--- a/Code/ChangeDetection/otbCorrelationChangeDetector.h
+++ b/Code/ChangeDetection/otbCorrelationChangeDetector.h
@@ -19,6 +19,7 @@
#define __otbCorrelationChangeDetector_h
#include "otbBinaryFunctorNeighborhoodImageFilter.h"
+#include "otbCrossCorrelation.h"
namespace otb
{
@@ -46,65 +47,6 @@ namespace otb
*
* \ingroup IntensityImageFilters Multithreaded
*/
-namespace Functor
-{
-
-template< class TInput1, class TInput2, class TOutput>
-class CrossCorrelation
-{
-public:
- CrossCorrelation() {};
- virtual ~CrossCorrelation() {};
- inline TOutput operator()( const TInput1 & itA,
- const TInput2 & itB)
- {
-
- TOutput meanA = itk::NumericTraits<TOutput>::Zero;
- TOutput meanB = itk::NumericTraits<TOutput>::Zero;
-
- for (unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
-
- meanA += static_cast<TOutput>(itA.GetPixel(pos));
- meanB += static_cast<TOutput>(itB.GetPixel(pos));
-
-
- }
-
- meanA /= itA.Size();
- meanB /= itB.Size();
-
- TOutput varA = itk::NumericTraits<TOutput>::Zero;
- TOutput varB = itk::NumericTraits<TOutput>::Zero;
-
- for (unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
-
- varA += static_cast<TOutput>( vcl_pow( static_cast<double>(itA.GetPixel(pos))-static_cast<double>(meanA),static_cast<double>(2.0)));
- varB += static_cast<TOutput>( vcl_pow( static_cast<double>(itB.GetPixel(pos))-static_cast<double>(meanB),static_cast<double>(2.0)));
-
- }
-
- varA /= itA.Size();
- varB /= itB.Size();
-
- TOutput crossCorrel = itk::NumericTraits<TOutput>::Zero;
-
- if (varA!= itk::NumericTraits<TOutput>::Zero && varB!= itk::NumericTraits<TOutput>::Zero)
- {
- for (unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
- crossCorrel += (static_cast<TOutput>(itA.GetPixel(pos))-meanA)*(static_cast<TOutput>(itB.GetPixel(pos))-meanB)/(itA.Size()*vcl_sqrt(varA*varB));
- }
- }
- else if (varA==itk::NumericTraits<TOutput>::Zero && varB==itk::NumericTraits<TOutput>::Zero)
- {
- crossCorrel = itk::NumericTraits<TOutput>::One;
- }
- return static_cast<TOutput>( itk::NumericTraits<TOutput>::One - crossCorrel );
- }
-};
-}
template <class TInputImage1, class TInputImage2, class TOutputImage>
class ITK_EXPORT CorrelationChangeDetector :
diff --git a/Code/ChangeDetection/otbCrossCorrelation.h b/Code/ChangeDetection/otbCrossCorrelation.h
new file mode 100644
index 0000000000000000000000000000000000000000..7d179ff353815c322f3d62f8a65587082e84be59
--- /dev/null
+++ b/Code/ChangeDetection/otbCrossCorrelation.h
@@ -0,0 +1,101 @@
+/*=========================================================================
+
+ 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 __otbCrossCorrelation_h
+#define __otbCrossCorrelationChangeDetector_h
+
+
+namespace otb
+{
+
+
+namespace Functor
+{
+
+/** \class Functor::CrossCorrelation
+ *
+ * - cast the input 1 pixel value to \c double
+ * - cast the input 2 pixel value to \c double
+ * - compute the difference of the two pixel values
+ * - compute the value of the cross-correlation
+ * - cast the \c double value resulting to the pixel type of the output image
+ * - store the casted value into the output image.
+ *
+ */
+
+
+template< class TInput1, class TInput2, class TOutput>
+class CrossCorrelation
+{
+public:
+ CrossCorrelation() {};
+ virtual ~CrossCorrelation() {};
+ inline TOutput operator()( const TInput1 & itA,
+ const TInput2 & itB)
+ {
+
+ TOutput meanA = itk::NumericTraits<TOutput>::Zero;
+ TOutput meanB = itk::NumericTraits<TOutput>::Zero;
+
+ for (unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+
+ meanA += static_cast<TOutput>(itA.GetPixel(pos));
+ meanB += static_cast<TOutput>(itB.GetPixel(pos));
+
+
+ }
+
+ meanA /= itA.Size();
+ meanB /= itB.Size();
+
+ TOutput varA = itk::NumericTraits<TOutput>::Zero;
+ TOutput varB = itk::NumericTraits<TOutput>::Zero;
+
+ for (unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+
+ varA += static_cast<TOutput>( vcl_pow( static_cast<double>(itA.GetPixel(pos))-static_cast<double>(meanA),static_cast<double>(2.0)));
+ varB += static_cast<TOutput>( vcl_pow( static_cast<double>(itB.GetPixel(pos))-static_cast<double>(meanB),static_cast<double>(2.0)));
+
+ }
+
+ varA /= itA.Size();
+ varB /= itB.Size();
+
+ TOutput crossCorrel = itk::NumericTraits<TOutput>::Zero;
+
+ if (varA!= itk::NumericTraits<TOutput>::Zero && varB!= itk::NumericTraits<TOutput>::Zero)
+ {
+ for (unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+ crossCorrel += (static_cast<TOutput>(itA.GetPixel(pos))-meanA)*(static_cast<TOutput>(itB.GetPixel(pos))-meanB)/(itA.Size()*vcl_sqrt(varA*varB));
+ }
+ }
+ else if (varA==itk::NumericTraits<TOutput>::Zero && varB==itk::NumericTraits<TOutput>::Zero)
+ {
+ crossCorrel = itk::NumericTraits<TOutput>::One;
+ }
+ return static_cast<TOutput>( itk::NumericTraits<TOutput>::One - crossCorrel );
+ }
+};
+}
+
+} // end namespace otb
+
+
+#endif
diff --git a/Code/ChangeDetection/otbJoinHistogramMI.h b/Code/ChangeDetection/otbJoinHistogramMI.h
new file mode 100644
index 0000000000000000000000000000000000000000..35eeef25f1ef98a4c0f39cadca2d3a45d97b8e1e
--- /dev/null
+++ b/Code/ChangeDetection/otbJoinHistogramMI.h
@@ -0,0 +1,136 @@
+/*=========================================================================
+
+ 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 __otbJoinHistogramMI_h
+#define __otbJoinHistogramMI_h
+
+#include "itkHistogram.h"
+
+namespace otb
+{
+
+namespace Functor
+{
+
+template< class TInput1, class TInput2, class TOutput>
+class JoinHistogramMI
+{
+public:
+ typedef double HistogramFrequencyType;
+ typedef typename itk::Statistics::Histogram<HistogramFrequencyType, 2> HistogramType;
+ JoinHistogramMI() {};
+ virtual ~JoinHistogramMI() {};
+ inline TOutput operator()( const TInput1 & itA,
+ const TInput2 & itB, const HistogramType* histogram)
+ {
+ TOutput jointEntropy = itk::NumericTraits<TOutput>::Zero;
+ HistogramFrequencyType totalFreq = histogram->GetTotalFrequency();
+
+ /* for(unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+ double value = static_cast<double>(itA.GetPixel(pos));
+
+ unsigned int bin =
+ HistogramFrequencyType freq = histogram.GetFrequency(, 0);
+ if (freq > 0)
+ {
+ entropyX += freq*vcl_log(freq);
+ }
+ }
+
+ entropyX = -entropyX/static_cast<TOutput>(totalFreq) + vcl_log(totalFreq);
+
+ for (unsigned int i = 0; i < this->GetHistogramSize()[1]; ++i)
+ {
+ HistogramFrequencyType freq = histogram.GetFrequency(i, 1);
+ if (freq > 0)
+ {
+ entropyY += freq*vcl_log(freq);
+ }
+ }
+
+ entropyY = -entropyY/static_cast<TOutput>(totalFreq) + vcl_log(totalFreq);
+
+ HistogramIteratorType it = histogram.Begin();
+ HistogramIteratorType end = histogram.End();
+ while (it != end)
+ {
+ HistogramFrequencyType freq = it.GetFrequency();
+ if (freq > 0)
+ {
+ jointEntropy += freq*vcl_log(freq);
+ }
+ ++it;
+ }
+
+ jointEntropy = -jointEntropy/static_cast<TOutput>(totalFreq) +
+ vcl_log(totalFreq);
+
+ return entropyX + entropyY - jointEntropy;*/
+
+
+ typename HistogramType::MeasurementVectorType sample;
+ for (unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+ double valueA = static_cast<double>(itA.GetPixel(pos));
+ double valueB = static_cast<double>(itB.GetPixel(pos));
+
+ sample[0] = valueA;
+ sample[1] = valueB;
+
+
+ HistogramFrequencyType freq = histogram->GetFrequency(
+ histogram->GetIndex(sample));
+ if (freq > 0)
+ {
+ jointEntropy += freq*vcl_log(freq);
+ }
+
+ }
+
+ jointEntropy = -jointEntropy/static_cast<TOutput>(totalFreq) +
+ vcl_log(totalFreq);
+
+ return jointEntropy;
+
+ /* TOutput meanA = 0.0;
+ TOutput meanB = 0.0;
+
+ for(unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+
+ meanA += static_cast<TOutput>(itA.GetPixel(pos));
+ meanB += static_cast<TOutput>(itB.GetPixel(pos));
+
+
+ }*/
+ return static_cast<TOutput>( 0 );
+ }
+
+ /* void SetHistogram(HistogramType* histo)
+ {
+ m_Histogram = histo;
+ }
+
+ protected:
+ HistogramType::Pointer m_Histogram;*/
+};
+}
+} // end namespace otb
+
+
+#endif
diff --git a/Code/ChangeDetection/otbJoinHistogramMIImageFilter.h b/Code/ChangeDetection/otbJoinHistogramMIImageFilter.h
index 65e59e853836f85d84adcc86d076055562feff99..3117a2bba829983312f52e7d0acd65a5fb252312 100644
--- a/Code/ChangeDetection/otbJoinHistogramMIImageFilter.h
+++ b/Code/ChangeDetection/otbJoinHistogramMIImageFilter.h
@@ -20,6 +20,7 @@
#include "otbBinaryFunctorNeighborhoodJoinHistogramImageFilter.h"
#include "itkHistogram.h"
+#include "otbJoinHistogramMI.h"
namespace otb
{
@@ -48,113 +49,7 @@ namespace otb
*
* \ingroup IntensityImageFilters Multithreaded
*/
-namespace Functor
-{
-template< class TInput1, class TInput2, class TOutput>
-class JoinHistogramMI
-{
-public:
- typedef double HistogramFrequencyType;
- typedef typename itk::Statistics::Histogram<HistogramFrequencyType, 2> HistogramType;
- JoinHistogramMI() {};
- virtual ~JoinHistogramMI() {};
- inline TOutput operator()( const TInput1 & itA,
- const TInput2 & itB, const HistogramType* histogram)
- {
- TOutput jointEntropy = itk::NumericTraits<TOutput>::Zero;
- HistogramFrequencyType totalFreq = histogram->GetTotalFrequency();
-
- /* for(unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
- double value = static_cast<double>(itA.GetPixel(pos));
-
- unsigned int bin =
- HistogramFrequencyType freq = histogram.GetFrequency(, 0);
- if (freq > 0)
- {
- entropyX += freq*vcl_log(freq);
- }
- }
-
- entropyX = -entropyX/static_cast<TOutput>(totalFreq) + vcl_log(totalFreq);
-
- for (unsigned int i = 0; i < this->GetHistogramSize()[1]; ++i)
- {
- HistogramFrequencyType freq = histogram.GetFrequency(i, 1);
- if (freq > 0)
- {
- entropyY += freq*vcl_log(freq);
- }
- }
-
- entropyY = -entropyY/static_cast<TOutput>(totalFreq) + vcl_log(totalFreq);
-
- HistogramIteratorType it = histogram.Begin();
- HistogramIteratorType end = histogram.End();
- while (it != end)
- {
- HistogramFrequencyType freq = it.GetFrequency();
- if (freq > 0)
- {
- jointEntropy += freq*vcl_log(freq);
- }
- ++it;
- }
-
- jointEntropy = -jointEntropy/static_cast<TOutput>(totalFreq) +
- vcl_log(totalFreq);
-
- return entropyX + entropyY - jointEntropy;*/
-
-
- typename HistogramType::MeasurementVectorType sample;
- for (unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
- double valueA = static_cast<double>(itA.GetPixel(pos));
- double valueB = static_cast<double>(itB.GetPixel(pos));
-
- sample[0] = valueA;
- sample[1] = valueB;
-
-
- HistogramFrequencyType freq = histogram->GetFrequency(
- histogram->GetIndex(sample));
- if (freq > 0)
- {
- jointEntropy += freq*vcl_log(freq);
- }
-
- }
-
- jointEntropy = -jointEntropy/static_cast<TOutput>(totalFreq) +
- vcl_log(totalFreq);
-
- return jointEntropy;
-
- /* TOutput meanA = 0.0;
- TOutput meanB = 0.0;
-
- for(unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
-
- meanA += static_cast<TOutput>(itA.GetPixel(pos));
- meanB += static_cast<TOutput>(itB.GetPixel(pos));
-
-
- }*/
- return static_cast<TOutput>( 0 );
- }
-
- /* void SetHistogram(HistogramType* histo)
- {
- m_Histogram = histo;
- }
-
- protected:
- HistogramType::Pointer m_Histogram;*/
-};
-}
template <class TInputImage1, class TInputImage2, class TOutputImage>
class ITK_EXPORT JoinHistogramMIImageFilter :
diff --git a/Code/ChangeDetection/otbLHMI.h b/Code/ChangeDetection/otbLHMI.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2fbae9b313d4fc238dbba2512a08f38587d08d1
--- /dev/null
+++ b/Code/ChangeDetection/otbLHMI.h
@@ -0,0 +1,179 @@
+/*=========================================================================
+
+ 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 __otbLHMI_h
+#define __otbLHMI_h
+
+#include <vector>
+#include <stdlib.h>
+#include <math.h>
+#include "itkHistogram.h"
+
+
+namespace otb
+{
+
+
+#define epsilon 0.01
+
+namespace Functor
+{
+
+template< class TInput1, class TInput2, class TOutput>
+class LHMI
+{
+
+public:
+
+ typedef typename std::vector<TOutput> VectorType;
+ typedef typename VectorType::iterator IteratorType;
+ typedef typename std::vector<VectorType> VectorOfVectorType;
+ typedef typename VectorOfVectorType::iterator VecOfVecIteratorType;
+ typedef double HistogramFrequencyType;
+ typedef typename itk::Statistics::Histogram<HistogramFrequencyType, 2> HistogramType;
+ typedef typename HistogramType::MeasurementVectorType
+ MeasurementVectorType;
+ typedef typename HistogramType::SizeType HistogramSizeType;
+ typedef typename HistogramType::Iterator HistogramIteratorType;
+
+
+ LHMI() {};
+ virtual ~LHMI() {};
+ inline TOutput operator()( const TInput1 & itA,
+ const TInput2 & itB)
+ {
+
+
+ HistogramType::Pointer histogram;
+
+ /** The histogram size. */
+ HistogramSizeType histogramSize;
+ /** The lower bound for samples in the histogram. */
+ MeasurementVectorType lowerBound;
+ /** The upper bound for samples in the histogram. */
+ MeasurementVectorType upperBound;
+
+ double upperBoundIncreaseFactor = 0.001;
+
+ histogramSize.Fill(256);
+
+ TOutput maxA = itA.GetPixel(0);
+ TOutput minA = itA.GetPixel(0);
+ TOutput maxB = itB.GetPixel(0);
+ TOutput minB = itB.GetPixel(0);
+
+ for (unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+
+ TOutput value = static_cast<TOutput>(itA.GetPixel(pos));
+
+ if (value > maxA)
+ maxA = value;
+ else if (value < minA)
+ minA = value;
+
+ value = static_cast<TOutput>(itB.GetPixel(pos));
+
+ if (value > maxB)
+ maxB = value;
+ else if (value < minB)
+ minB = value;
+
+
+ }
+
+
+ // Initialize the upper and lower bounds of the histogram.
+ lowerBound[0] = minA;
+ lowerBound[1] = minB;
+ upperBound[0] =
+ maxA + (maxA - minA ) * upperBoundIncreaseFactor;
+ upperBound[1] =
+ maxB + (maxB - minB ) * upperBoundIncreaseFactor;
+
+
+ histogram = HistogramType::New();
+
+ histogram->Initialize(histogramSize, lowerBound, upperBound);
+
+ for (unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+
+ typename HistogramType::MeasurementVectorType sample;
+ sample[0] = itA.GetPixel(pos);
+ sample[1] = itB.GetPixel(pos);
+ /*if(sample[0]!=NumericTraits<TOutput>::Zero &&
+ sample[1]!=NumericTraits<TOutput>::Zero)*/
+ histogram->IncreaseFrequency(sample, 1);
+
+ }
+
+
+
+ TOutput entropyX = itk::NumericTraits<TOutput>::Zero;
+ TOutput entropyY = itk::NumericTraits<TOutput>::Zero;
+ TOutput jointEntropy = itk::NumericTraits<TOutput>::Zero;
+ HistogramFrequencyType totalFreq = histogram->GetTotalFrequency();
+
+ for (unsigned int i = 0; i < histogram->GetSize()[0]; ++i)
+ {
+ HistogramFrequencyType freq = histogram->GetFrequency(i, 0);
+ if (freq > 0)
+ {
+ entropyX += freq*vcl_log(freq);
+ }
+ }
+
+ entropyX = -entropyX/static_cast<TOutput>(totalFreq) + vcl_log(totalFreq);
+
+ for (unsigned int i = 0; i < histogram->GetSize()[1]; ++i)
+ {
+ HistogramFrequencyType freq = histogram->GetFrequency(i, 1);
+ if (freq > 0)
+ {
+ entropyY += freq*vcl_log(freq);
+ }
+ }
+
+ entropyY = -entropyY/static_cast<TOutput>(totalFreq) + vcl_log(totalFreq);
+
+ HistogramIteratorType it = histogram->Begin();
+ HistogramIteratorType end = histogram->End();
+ while (it != end)
+ {
+ HistogramFrequencyType freq = it.GetFrequency();
+ if (freq > 0)
+ {
+ jointEntropy += freq*vcl_log(freq);
+ }
+ ++it;
+ }
+
+ jointEntropy = -jointEntropy/static_cast<TOutput>(totalFreq) +
+ vcl_log(totalFreq);
+
+ return static_cast<TOutput>( jointEntropy/(entropyX + entropyY) );
+ }
+
+
+
+};
+}
+
+}
+#endif
diff --git a/Code/ChangeDetection/otbLHMIChangeDetector.h b/Code/ChangeDetection/otbLHMIChangeDetector.h
index 9383c213de62536f80c274ef9c70663532dae198..621618edb9bea87cfc754f21af93d52240a24cd8 100644
--- a/Code/ChangeDetection/otbLHMIChangeDetector.h
+++ b/Code/ChangeDetection/otbLHMIChangeDetector.h
@@ -19,11 +19,7 @@
#define __otbLHMIChangeDetector_h
#include "otbBinaryFunctorNeighborhoodImageFilter.h"
-#include <vector>
-#include <stdlib.h>
-#include <math.h>
-#include "itkHistogram.h"
-
+#include "otbLHMI.h"
namespace otb
{
@@ -53,152 +49,6 @@ namespace otb
* \ingroup IntensityImageFilters Multithreaded
*/
-#define epsilon 0.01
-
-namespace Functor
-{
-
-template< class TInput1, class TInput2, class TOutput>
-class LHMI
-{
-
-public:
-
- typedef typename std::vector<TOutput> VectorType;
- typedef typename VectorType::iterator IteratorType;
- typedef typename std::vector<VectorType> VectorOfVectorType;
- typedef typename VectorOfVectorType::iterator VecOfVecIteratorType;
- typedef double HistogramFrequencyType;
- typedef typename itk::Statistics::Histogram<HistogramFrequencyType, 2> HistogramType;
- typedef typename HistogramType::MeasurementVectorType
- MeasurementVectorType;
- typedef typename HistogramType::SizeType HistogramSizeType;
- typedef typename HistogramType::Iterator HistogramIteratorType;
-
-
- LHMI() {};
- virtual ~LHMI() {};
- inline TOutput operator()( const TInput1 & itA,
- const TInput2 & itB)
- {
-
-
- HistogramType::Pointer histogram;
-
- /** The histogram size. */
- HistogramSizeType histogramSize;
- /** The lower bound for samples in the histogram. */
- MeasurementVectorType lowerBound;
- /** The upper bound for samples in the histogram. */
- MeasurementVectorType upperBound;
-
- double upperBoundIncreaseFactor = 0.001;
-
- histogramSize.Fill(256);
-
- TOutput maxA = itA.GetPixel(0);
- TOutput minA = itA.GetPixel(0);
- TOutput maxB = itB.GetPixel(0);
- TOutput minB = itB.GetPixel(0);
-
- for (unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
-
- TOutput value = static_cast<TOutput>(itA.GetPixel(pos));
-
- if (value > maxA)
- maxA = value;
- else if (value < minA)
- minA = value;
-
- value = static_cast<TOutput>(itB.GetPixel(pos));
-
- if (value > maxB)
- maxB = value;
- else if (value < minB)
- minB = value;
-
-
- }
-
-
- // Initialize the upper and lower bounds of the histogram.
- lowerBound[0] = minA;
- lowerBound[1] = minB;
- upperBound[0] =
- maxA + (maxA - minA ) * upperBoundIncreaseFactor;
- upperBound[1] =
- maxB + (maxB - minB ) * upperBoundIncreaseFactor;
-
-
- histogram = HistogramType::New();
-
- histogram->Initialize(histogramSize, lowerBound, upperBound);
-
- for (unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
-
- typename HistogramType::MeasurementVectorType sample;
- sample[0] = itA.GetPixel(pos);
- sample[1] = itB.GetPixel(pos);
- /*if(sample[0]!=NumericTraits<TOutput>::Zero &&
- sample[1]!=NumericTraits<TOutput>::Zero)*/
- histogram->IncreaseFrequency(sample, 1);
-
- }
-
-
-
- TOutput entropyX = itk::NumericTraits<TOutput>::Zero;
- TOutput entropyY = itk::NumericTraits<TOutput>::Zero;
- TOutput jointEntropy = itk::NumericTraits<TOutput>::Zero;
- HistogramFrequencyType totalFreq = histogram->GetTotalFrequency();
-
- for (unsigned int i = 0; i < histogram->GetSize()[0]; ++i)
- {
- HistogramFrequencyType freq = histogram->GetFrequency(i, 0);
- if (freq > 0)
- {
- entropyX += freq*vcl_log(freq);
- }
- }
-
- entropyX = -entropyX/static_cast<TOutput>(totalFreq) + vcl_log(totalFreq);
-
- for (unsigned int i = 0; i < histogram->GetSize()[1]; ++i)
- {
- HistogramFrequencyType freq = histogram->GetFrequency(i, 1);
- if (freq > 0)
- {
- entropyY += freq*vcl_log(freq);
- }
- }
-
- entropyY = -entropyY/static_cast<TOutput>(totalFreq) + vcl_log(totalFreq);
-
- HistogramIteratorType it = histogram->Begin();
- HistogramIteratorType end = histogram->End();
- while (it != end)
- {
- HistogramFrequencyType freq = it.GetFrequency();
- if (freq > 0)
- {
- jointEntropy += freq*vcl_log(freq);
- }
- ++it;
- }
-
- jointEntropy = -jointEntropy/static_cast<TOutput>(totalFreq) +
- vcl_log(totalFreq);
-
- return static_cast<TOutput>( jointEntropy/(entropyX + entropyY) );
- }
-
-
-
-};
-}
-
template <class TInputImage1, class TInputImage2, class TOutputImage>
class ITK_EXPORT LHMIChangeDetector :
public BinaryFunctorNeighborhoodImageFilter<
diff --git a/Code/ChangeDetection/otbMeanDifference.h b/Code/ChangeDetection/otbMeanDifference.h
new file mode 100644
index 0000000000000000000000000000000000000000..c908f315607f9120ac48cd7d98bcc295593d34ac
--- /dev/null
+++ b/Code/ChangeDetection/otbMeanDifference.h
@@ -0,0 +1,78 @@
+/*=========================================================================
+
+ 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 __otbMeanDifference_h
+#define __otbMeanDifference_h
+
+
+namespace otb
+{
+
+/** \class MeanDifferenceImageFilter
+ * \brief Implements neighborhood-wise the computation of mean difference.
+ *
+ * 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 difference of the two pixel values
+ * - compute the value of the difference of means
+ * - 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)
+ *
+ * \ingroup IntensityImageFilters Multithreaded
+ */
+namespace Functor
+{
+
+template< class TInput1, class TInput2, class TOutput>
+class MeanDifference
+{
+public:
+ MeanDifference() {};
+ virtual ~MeanDifference() {};
+ inline TOutput operator()( const TInput1 & itA,
+ const TInput2 & itB)
+ {
+
+ TOutput meanA = 0.0;
+ TOutput meanB = 0.0;
+
+ for (unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+
+ meanA += static_cast<TOutput>(itA.GetPixel(pos));
+ meanB += static_cast<TOutput>(itB.GetPixel(pos));
+
+
+ }
+ return static_cast<TOutput>( (meanA-meanB)/itA.Size() );
+ }
+};
+}
+}
+
+#endif
diff --git a/Code/ChangeDetection/otbMeanDifferenceImageFilter.h b/Code/ChangeDetection/otbMeanDifferenceImageFilter.h
index c009eeb5893e9b50586ed208b3d3fdfec999e27e..9e182c82ecb9e17a3067ff37cfaf6d9869e3e374 100644
--- a/Code/ChangeDetection/otbMeanDifferenceImageFilter.h
+++ b/Code/ChangeDetection/otbMeanDifferenceImageFilter.h
@@ -19,62 +19,11 @@
#define __otbMeanDifferenceImageFilter_h
#include "otbBinaryFunctorNeighborhoodImageFilter.h"
+#include "otbMeanDifference.h"
namespace otb
{
-/** \class MeanDifferenceImageFilter
- * \brief Implements neighborhood-wise the computation of mean difference.
- *
- * 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 difference of the two pixel values
- * - compute the value of the difference of means
- * - 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)
- *
- * \ingroup IntensityImageFilters Multithreaded
- */
-namespace Functor
-{
-
-template< class TInput1, class TInput2, class TOutput>
-class MeanDifference
-{
-public:
- MeanDifference() {};
- virtual ~MeanDifference() {};
- inline TOutput operator()( const TInput1 & itA,
- const TInput2 & itB)
- {
-
- TOutput meanA = 0.0;
- TOutput meanB = 0.0;
-
- for (unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
-
- meanA += static_cast<TOutput>(itA.GetPixel(pos));
- meanB += static_cast<TOutput>(itB.GetPixel(pos));
-
-
- }
- return static_cast<TOutput>( (meanA-meanB)/itA.Size() );
- }
-};
-}
-
template <class TInputImage1, class TInputImage2, class TOutputImage>
class ITK_EXPORT MeanDifferenceImageFilter :
public BinaryFunctorNeighborhoodImageFilter<
diff --git a/Code/ChangeDetection/otbMeanRatio.h b/Code/ChangeDetection/otbMeanRatio.h
new file mode 100644
index 0000000000000000000000000000000000000000..dbe0ba5a91be1e0892d57bd56d5f4871b996be6b
--- /dev/null
+++ b/Code/ChangeDetection/otbMeanRatio.h
@@ -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.
+
+=========================================================================*/
+#ifndef __otbMeanRatio_h
+#define __otbMeanRatio_h
+
+#include "otbBinaryFunctorNeighborhoodImageFilter.h"
+
+namespace otb
+{
+
+/** \class Functor::MeanRatio
+ *
+ * - compute the ratio of the two pixel values
+ * - compute the value of the ratio of means
+ * - cast the \c double value resulting to the pixel type of the output image
+ * - store the casted value into the output image.
+ *
+
+ * \ingroup Functor
+ */
+namespace Functor
+{
+
+template< class TInput1, class TInput2, class TOutput>
+class MeanRatio
+{
+public:
+ MeanRatio() {};
+ virtual ~MeanRatio() {};
+ inline TOutput operator()( const TInput1 & itA,
+ const TInput2 & itB)
+ {
+
+ TOutput meanA = 0.0;
+ TOutput meanB = 0.0;
+
+ for (unsigned long pos = 0; pos< itA.Size(); ++pos)
+ {
+
+ meanA += static_cast<TOutput>(itA.GetPixel(pos));
+ meanB += static_cast<TOutput>(itB.GetPixel(pos));
+
+
+ }
+
+ meanA /= itA.Size();
+ meanB /= itB.Size();
+
+ //std::cout<<"meanA= "<<meanA<<", meanB= "<<meanB<<std::endl;
+
+ TOutput ratio;
+
+ if(meanA == meanB)
+ ratio = 0.;
+ else if (meanA>meanB)
+ ratio = static_cast<TOutput>(1.0 - meanB/meanA);
+ else ratio = static_cast<TOutput>(1.0 - meanA/meanB);
+
+ return ratio;
+ }
+};
+}
+} // end namespace otb
+
+
+#endif
diff --git a/Code/ChangeDetection/otbMeanRatioImageFilter.h b/Code/ChangeDetection/otbMeanRatioImageFilter.h
index 209aa3c400bca249fb2ddba8111f48d32e66d32a..9e77191fcbb302620fd3cacee2f96696ec909b0a 100644
--- a/Code/ChangeDetection/otbMeanRatioImageFilter.h
+++ b/Code/ChangeDetection/otbMeanRatioImageFilter.h
@@ -19,6 +19,7 @@
#define __otbMeanRatioImageFilter_h
#include "otbBinaryFunctorNeighborhoodImageFilter.h"
+#include "otbMeanRatio.h"
namespace otb
{
@@ -44,48 +45,6 @@ namespace otb
*
* \ingroup IntensityImageFilters Multithreaded
*/
-namespace Functor
-{
-
-template< class TInput1, class TInput2, class TOutput>
-class MeanRatio
-{
-public:
- MeanRatio() {};
- virtual ~MeanRatio() {};
- inline TOutput operator()( const TInput1 & itA,
- const TInput2 & itB)
- {
-
- TOutput meanA = 0.0;
- TOutput meanB = 0.0;
-
- for (unsigned long pos = 0; pos< itA.Size(); ++pos)
- {
-
- meanA += static_cast<TOutput>(itA.GetPixel(pos));
- meanB += static_cast<TOutput>(itB.GetPixel(pos));
-
-
- }
-
- meanA /= itA.Size();
- meanB /= itB.Size();
-
- //std::cout<<"meanA= "<<meanA<<", meanB= "<<meanB<<std::endl;
-
- TOutput ratio;
-
- if(meanA == meanB)
- ratio = 0.;
- else if (meanA>meanB)
- ratio = static_cast<TOutput>(1.0 - meanB/meanA);
- else ratio = static_cast<TOutput>(1.0 - meanA/meanB);
-
- return ratio;
- }
-};
-}
template <class TInputImage1, class TInputImage2, class TOutputImage>
class ITK_EXPORT MeanRatioImageFilter :