From 4199553b51884ff3607b3074781adfabef46d4ff Mon Sep 17 00:00:00 2001
From: Ludovic Hussonnois <ludovic.hussonnois@c-s.fr>
Date: Thu, 2 Mar 2017 17:33:19 +0100
Subject: [PATCH] ENH: Add ContingencyTable and Calculator for image and
 MeasurementVector

---
 .../include/otbContingencyTable.h             | 140 ++++++++++++++++++
 .../include/otbContingencyTableCalculator.h   |  96 ++++++++++++
 .../include/otbContingencyTableCalculator.txx | 120 +++++++++++++++
 3 files changed, 356 insertions(+)
 create mode 100644 Modules/Learning/Unsupervised/include/otbContingencyTable.h
 create mode 100644 Modules/Learning/Unsupervised/include/otbContingencyTableCalculator.h
 create mode 100644 Modules/Learning/Unsupervised/include/otbContingencyTableCalculator.txx

diff --git a/Modules/Learning/Unsupervised/include/otbContingencyTable.h b/Modules/Learning/Unsupervised/include/otbContingencyTable.h
new file mode 100644
index 0000000000..70178e45cb
--- /dev/null
+++ b/Modules/Learning/Unsupervised/include/otbContingencyTable.h
@@ -0,0 +1,140 @@
+/*=========================================================================
+
+  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 otbContingencyTable_h
+#define otbContingencyTable_h
+
+#include <vector>
+#include <iostream>
+#include <iomanip>
+
+namespace otb
+{
+template<class TClassLabel>
+class ContingencyTable
+{
+public:
+  typedef itk::VariableSizeMatrix<unsigned long> MatrixType;
+  typedef std::vector<TClassLabel> LabelList;
+
+  ContingencyTable(LabelList referenceLabels, LabelList producedLabels) : refLabels( referenceLabels ),
+                                                                          prodLabels( producedLabels )
+  {
+    unsigned int rows = static_cast<unsigned int>(refLabels.size());
+    unsigned int cols = static_cast<unsigned int>(prodLabels.size());
+    matrix.SetSize( rows, cols );
+    matrix.Fill( 0 );
+  }
+
+  MatrixType matrix;
+
+  friend std::ostream &operator<<(std::ostream &o, const ContingencyTable<TClassLabel> &contingencyTable)
+  {
+
+    // Retrieve the maximal width from the matrix and the labels
+    size_t maxWidth = 6;
+    for( size_t i = 0; i << contingencyTable.prodLabels.size(); ++i )
+      {
+      std::ostringstream oss;
+      oss << contingencyTable.prodLabels[i];
+      size_t length = oss.str().length();
+      if( length > maxWidth )
+        maxWidth = length;
+      }
+
+    for( size_t i = 0; i << contingencyTable.refLabels.size(); ++i )
+      {
+      std::ostringstream oss;
+      oss << contingencyTable.refLabels[i];
+      size_t length = oss.str().length();
+      if( length > maxWidth )
+        maxWidth = length;
+      }
+
+    for( unsigned int i = 0; i < contingencyTable.matrix.Rows(); ++i )
+      {
+      for( unsigned int j = 0; j < contingencyTable.matrix.Cols(); ++j )
+        {
+        std::ostringstream oss;
+        oss << contingencyTable.matrix( i, j );
+        size_t length = oss.str().length();
+        if( length > maxWidth )
+          maxWidth = length;
+        }
+      }
+
+    int width = static_cast<int>(maxWidth);
+
+    // Write the first line of the matrix (produced labels)
+    o << std::setfill(' ') << std::setw( width ) << "labels";
+    for( size_t i = 0; i < contingencyTable.prodLabels.size(); ++i )
+      {
+      o << std::setfill(' ') << std::setw( width ) << contingencyTable.prodLabels[i];
+      }
+    o << std::endl;
+
+    // For each line write the reference label, then the count value
+    for( unsigned int i = 0; i < contingencyTable.matrix.Rows(); ++i )
+      {
+      o << std::setfill(' ') << std::setw( width ) << contingencyTable.refLabels[i];
+      for( unsigned int j = 0; j < contingencyTable.matrix.Cols(); ++j )
+        {
+        o << std::setfill(' ') << std::setw( width ) << contingencyTable.matrix( i, j );
+        }
+      o << std::endl;
+      }
+
+    return o;
+  }
+
+  std::string to_csv() const
+  {
+    const char separator = ',';
+
+    std::ostringstream oss;
+    oss << "labels";
+    for( size_t i = 0; i < prodLabels.size(); ++i )
+      {
+      oss << separator << prodLabels[i];
+      }
+    oss << std::endl;
+
+    // For each line write the reference label, then the count value
+    for( unsigned int i = 0; i < matrix.Rows(); ++i )
+      {
+      oss << refLabels[i];
+      for( unsigned int j = 0; j < matrix.Cols(); ++j )
+        {
+        oss << separator << matrix( i, j );
+        }
+      oss << std::endl;
+      }
+    oss << std::endl;
+
+    return oss.str();
+  }
+
+private:
+  LabelList refLabels;
+  LabelList prodLabels;
+
+
+};
+}
+
+#endif //otbContingencyTable_h
diff --git a/Modules/Learning/Unsupervised/include/otbContingencyTableCalculator.h b/Modules/Learning/Unsupervised/include/otbContingencyTableCalculator.h
new file mode 100644
index 0000000000..f6f6df9731
--- /dev/null
+++ b/Modules/Learning/Unsupervised/include/otbContingencyTableCalculator.h
@@ -0,0 +1,96 @@
+/*=========================================================================
+
+ 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 otbContingencyTableCalculator_h
+#define otbContingencyTableCalculator_h
+
+#include "itkObject.h"
+#include "itkObjectFactory.h"
+#include "otbContingencyTable.h"
+#include <map>
+
+namespace otb
+{
+/**
+ * \brief ContingencyTableCalculator provide facilities to compute ContingencyTable with different structure type
+ * The size of the label list should be the same for both list.
+ * \tparam TRefListLabel data structure type which contain the reference labels
+ * \tparam TProdListLabel data structure type which contain the produced labels
+ * \tparam TClassLabel the label data type
+ *
+ * \ingroup OTBUnsupervised
+ */
+template<class TClassLabel>
+class ITK_EXPORT ContingencyTableCalculator : public itk::Object
+{
+public:
+  /** Standard class typedefs */
+  typedef ContingencyTableCalculator    Self;
+  typedef itk::Object                   Superclass;
+  typedef itk::SmartPointer <Self>      Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  /** Run-time type information (and related methods). */
+  itkTypeMacro(ContingencyTableCalculator, itk::Object);
+
+  /** Method for creation through the object factory. */
+  itkNewMacro(Self);
+
+  typedef ContingencyTable<TClassLabel> ContingencyTableType;
+
+  typedef typename std::map<TClassLabel, unsigned long> CountMapType;
+  typedef typename std::map<TClassLabel, CountMapType > MapOfClassesType;
+
+  /** Populate the confusion Matrix for a image iteration */
+  template<class TRefIterator, class TProdIterator>
+  void Compute(TRefIterator itRef, TProdIterator itProd);
+
+  /** Populate the confusion Matrix with input which provide GetMeasurementVector()[0] access */
+  template<class TRefIterator, class TProdIterator>
+  void Compute(TRefIterator refBegin, TRefIterator refEnd, TProdIterator prodBegin, TProdIterator prodEnd);
+
+  itkGetConstMacro(NumberOfRefClasses, unsigned long);
+  itkGetConstMacro(NumberOfProdClasses, unsigned long);
+  itkGetConstMacro(NumberOfSamples, unsigned long);
+
+  void Clear();
+  ContingencyTableType GetContingencyTable();
+
+protected:
+  ContingencyTableCalculator();
+  ~ContingencyTableCalculator() ITK_OVERRIDE {}
+  //void PrintSelf(std::ostream& os, itk::Indent indent) const ITK_OVERRIDE;
+
+private:
+  ContingencyTableCalculator(const Self &); //purposely not implemented
+  void operator=(const Self &); //purposely not implemented
+
+  MapOfClassesType m_LabelCount;
+  unsigned long m_NumberOfRefClasses;
+  unsigned long m_NumberOfProdClasses;
+  unsigned long m_NumberOfSamples;
+
+
+};
+}
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbContingencyTableCalculator.txx"
+#endif
+
+#endif //otbContingencyTableCalculator_h
diff --git a/Modules/Learning/Unsupervised/include/otbContingencyTableCalculator.txx b/Modules/Learning/Unsupervised/include/otbContingencyTableCalculator.txx
new file mode 100644
index 0000000000..22f0e04d5d
--- /dev/null
+++ b/Modules/Learning/Unsupervised/include/otbContingencyTableCalculator.txx
@@ -0,0 +1,120 @@
+/*=========================================================================
+
+ 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 otbContingencyTableCalculator_txx
+#define otbContingencyTableCalculator_txx
+
+#include "otbContingencyTableCalculator.h"
+#include "itkVariableLengthVector.h"
+#include "itkListSample.h"
+
+#include <set>
+
+namespace otb
+{
+template<class TClassLabel>
+ContingencyTableCalculator<TClassLabel>::ContingencyTableCalculator()
+: m_NumberOfRefClasses(0), m_NumberOfProdClasses(0), m_NumberOfSamples(0)
+{}
+
+template<class TClassLabel>
+void
+ContingencyTableCalculator<TClassLabel>
+::Clear()
+{
+  m_LabelCount.clear();
+  m_NumberOfRefClasses = 0;
+  m_NumberOfProdClasses = 0;
+  m_NumberOfSamples = 0;
+}
+
+template<class TClassLabel>
+template<class TRefIterator, class TProdIterator>
+void
+ContingencyTableCalculator<TClassLabel>
+::Compute(TRefIterator refBegin, TRefIterator refEnd, TProdIterator prodBegin, TProdIterator prodEnd)
+{
+    while( refBegin != refEnd && prodBegin != prodEnd )
+      {
+      ++m_LabelCount[refBegin.GetMeasurementVector()[0]][prodBegin.GetMeasurementVector()[0]];
+      ++refBegin;
+      ++prodBegin;
+      ++m_NumberOfSamples;
+      }
+}
+
+template<class TClassLabel>
+template<class TRefIterator, class TProdIterator>
+void
+ContingencyTableCalculator<TClassLabel>
+::Compute(TRefIterator itRef, TProdIterator itProd)
+{
+  while( !itRef.IsAtEnd() && !itProd.IsAtEnd() )
+    {
+    ++m_LabelCount[itRef.Get()][itProd.Get()];
+    ++itRef;
+    ++itRef;
+    ++m_NumberOfSamples;
+    }
+}
+
+
+template<class TClassLabel>
+typename ContingencyTableCalculator<TClassLabel>::ContingencyTableType
+ContingencyTableCalculator<TClassLabel>
+::GetContingencyTable()
+{
+  std::set<TClassLabel> refLabels;
+  std::set<TClassLabel> prodLabels;
+
+  // Retrieve all labels needed to iterate over all labelCount
+  for(typename MapOfClassesType::const_iterator refIt = m_LabelCount.begin(); refIt != m_LabelCount.end(); ++refIt)
+    {
+      refLabels.insert(refIt->first);
+      CountMapType cmt = refIt->second;
+      for(typename CountMapType::const_iterator prodIt = cmt.begin(); prodIt != cmt.end(); ++prodIt)
+        {
+        prodLabels.insert(prodIt->first);
+        }
+    }
+
+  m_NumberOfRefClasses = refLabels.size();
+  m_NumberOfProdClasses = prodLabels.size();
+
+  unsigned int rows = static_cast<unsigned int>(m_NumberOfRefClasses);
+  unsigned int cols = static_cast<unsigned int>(m_NumberOfProdClasses);
+
+  std::vector<TClassLabel> referenceLabels;
+  std::vector<TClassLabel> producedLabels;
+
+  std::copy(refLabels.begin(), refLabels.end(), std::back_inserter(referenceLabels));
+  std::copy(prodLabels.begin(), prodLabels.end(), std::back_inserter(producedLabels));
+
+  ContingencyTableType contingencyTable(referenceLabels, producedLabels);
+
+  for( unsigned int i = 0; i < rows; ++i )
+    for( unsigned int j = 0; j < cols; ++j )
+      contingencyTable.matrix(i,j) = m_LabelCount[referenceLabels[i]][producedLabels[j]];
+
+
+  return contingencyTable;
+}
+
+}
+
+#endif
-- 
GitLab