MRG: Merge remote-tracking branch 'remotes/origin/contingency_table' into develop

Modules/Applications/AppClassification/app/otbComputeConfusionMatrix.cxx

@@ -28,6 +28,8 @@
 #include "otbRAMDrivenAdaptativeStreamingManager.h"
 
 #include "otbConfusionMatrixMeasurements.h"
+#include "otbContingencyTableCalculator.h"
+#include "otbContingencyTable.h"
 
 namespace otb
 {
@@ -77,8 +79,27 @@ public:
   typedef ConfusionMatrixMeasurementsType::MapOfClassesType                     MapOfClassesType;
   typedef ConfusionMatrixMeasurementsType::MeasurementType                      MeasurementType;
 
+  typedef ContingencyTable<ClassLabelType>  ContingencyTableType;
+  typedef ContingencyTableType::Pointer     ContingencyTablePointerType;
+
 
 private:
+  Int32ImageType* m_Input;
+  Int32ImageType::Pointer m_Reference;
+  RAMDrivenAdaptativeStreamingManagerType::Pointer m_StreamingManager;
+  otb::ogr::DataSource::Pointer m_OgrRef;
+  RasterizeFilterType::Pointer m_RasterizeReference;
+
+  struct StreamingInitializationData
+  {
+    bool refhasnodata;
+    bool prodhasnodata;
+    int  prodnodata;
+    int  refnodata;
+    unsigned long numberOfStreamDivisions;
+  };
+
+
   void DoInit() ITK_OVERRIDE
   {
   SetName("ComputeConfusionMatrix");
@@ -103,6 +124,12 @@ private:
   AddParameter(ParameterType_OutputFilename, "out", "Matrix output");
   SetParameterDescription("out", "Filename to store the output matrix (csv format)");
 
+  AddParameter(ParameterType_Choice,"format","set the output format to contingency table or confusion matrix");
+  SetParameterDescription("format","Choice of the output format as a contingency table for unsupervised algorithms"
+          "or confusion matrix for supervised ones.");
+  AddChoice("format.confusionmatrix","Choice of a confusion matrix as output.");
+  AddChoice("format.contingencytable","Choice of a contingency table as output.");
+
   AddParameter(ParameterType_Choice,"ref","Ground truth");
   SetParameterDescription("ref","Choice of ground truth format");
 
@@ -120,10 +147,23 @@ private:
   SetParameterDescription("ref.vector.field","Field name containing the label values");
   SetListViewSingleSelectionMode("ref.vector.field",true);
   
-  AddParameter(ParameterType_Int,"nodatalabel","Value for nodata pixels");
-  SetParameterDescription("nodatalabel", "Label for the NoData class. Such input pixels will be discarded from the "
-      "ground truth and from the input classification map. By default, 'nodatalabel = 0'.");
+  AddParameter(ParameterType_Int,"ref.raster.nodata","Value for nodata pixels in ref raster");
+  SetDefaultParameterInt("ref.raster.nodata",0);
+  SetParameterDescription("ref.raster.nodata","Label to be treated as no data in ref raster.");
+  MandatoryOff("ref.raster.nodata");
+  DisableParameter("ref.raster.nodata");
+
+  AddParameter(ParameterType_Int,"ref.vector.nodata","Value for nodata pixels in ref vector");
+  SetDefaultParameterInt("ref.vector.nodata",0);
+  SetParameterDescription("ref.vector.nodata","Label to be treated as no data in ref vector. Please note that this value is always used in vector mode, to generate default values. Please set it to a value that does not correspond to a class label.");
+  MandatoryOff("ref.vector.nodata");
+  DisableParameter("ref.vector.nodata");
+
+
+  AddParameter(ParameterType_Int,"nodatalabel","Value for nodata pixels in input image");
+  SetParameterDescription("nodatalabel","Label to be treated as no data in input image");
   SetDefaultParameterInt("nodatalabel",0);
+  
   MandatoryOff("nodatalabel");
   DisableParameter("nodatalabel");
 
@@ -135,7 +175,7 @@ private:
   SetDocExampleParameterValue("ref", "vector");
   SetDocExampleParameterValue("ref.vector.in","VectorData_QB1_bis.shp");
   SetDocExampleParameterValue("ref.vector.field","Class");
-  SetDocExampleParameterValue("nodatalabel","255");
+  SetDocExampleParameterValue("ref.vector.nodata","255");
   }
 
   void DoUpdateParameters() ITK_OVERRIDE
@@ -168,6 +208,19 @@ private:
       }    
   }
 
+  void LogContingencyTable(const ContingencyTablePointerType& contingencyTable)
+  {
+    otbAppLogINFO("Contingency table: reference labels (rows) vs. produced labels (cols)\n" << (*contingencyTable.GetPointer()));
+  }
+
+  void m_WriteContingencyTable(const ContingencyTablePointerType& contingencyTable)
+  {
+    std::ofstream outFile;
+    outFile.open( this->GetParameterString( "out" ).c_str() );
+    outFile << contingencyTable->ToCSV();
+    outFile.close();
+  }
+
   std::string LogConfusionMatrix(MapOfClassesType* mapOfClasses, ConfusionMatrixType* matrix)
     {
       // Compute minimal width
@@ -240,61 +293,118 @@ private:
     }
 
 
-  void DoExecute() ITK_OVERRIDE
+  StreamingInitializationData InitStreamingData()
   {
-    Int32ImageType* input = this->GetParameterInt32Image("in");
 
-    std::string field;
-    int nodata = this->GetParameterInt("nodatalabel");
+    StreamingInitializationData sid;
 
+    m_Input = this->GetParameterInt32Image("in");
 
-    Int32ImageType::Pointer reference;
-    otb::ogr::DataSource::Pointer ogrRef;
-    RasterizeFilterType::Pointer rasterizeReference = RasterizeFilterType::New();
+    std::string field;
+
+    sid.prodnodata = this->GetParameterInt("nodatalabel");
+    sid.prodhasnodata = this->IsParameterEnabled("nodatalabel");
 
     if (GetParameterString("ref") == "raster")
       {
-      reference = this->GetParameterInt32Image("ref.raster.in");
+      sid.refnodata = this->GetParameterInt("ref.raster.nodata");
+      sid.refhasnodata = this->IsParameterEnabled("ref.raster.nodata");
+      m_Reference = this->GetParameterInt32Image("ref.raster.in");
       }
     else
       {
-      ogrRef = otb::ogr::DataSource::New(GetParameterString("ref.vector.in"), otb::ogr::DataSource::Modes::Read);
-
-    // Get field name
-    std::vector<int> selectedCFieldIdx = GetSelectedItems("ref.vector.field");
-    
-    if(selectedCFieldIdx.empty())
-      {
-      otbAppLogFATAL(<<"No field has been selected for data labelling!");
-      }
+      // Force nodata to true since it will be generated during rasterization
+      sid.refhasnodata = true;
+      sid.refnodata = this->GetParameterInt("ref.vector.nodata");
       
-      std::vector<std::string> cFieldNames = GetChoiceNames("ref.vector.field");  
+      m_OgrRef = otb::ogr::DataSource::New(GetParameterString("ref.vector.in"), otb::ogr::DataSource::Modes::Read);
+
+      // Get field name
+      std::vector<int> selectedCFieldIdx = GetSelectedItems("ref.vector.field");
+
+      if(selectedCFieldIdx.empty())
+        {
+        otbAppLogFATAL(<<"No field has been selected for data labelling!");
+        }
+
+      std::vector<std::string> cFieldNames = GetChoiceNames("ref.vector.field");
       field = cFieldNames[selectedCFieldIdx.front()];
-      
-      rasterizeReference->AddOGRDataSource(ogrRef);
-      rasterizeReference->SetOutputParametersFromImage(input);
-      rasterizeReference->SetBackgroundValue(nodata);
-      rasterizeReference->SetBurnAttribute(field.c_str());
 
-      reference = rasterizeReference->GetOutput();
-      reference->UpdateOutputInformation();
+      m_RasterizeReference = RasterizeFilterType::New();
+      m_RasterizeReference->AddOGRDataSource(m_OgrRef);
+      m_RasterizeReference->SetOutputParametersFromImage(m_Input);
+      m_RasterizeReference->SetBackgroundValue(sid.refnodata);
+      m_RasterizeReference->SetBurnAttribute(field.c_str());
+
+      m_Reference = m_RasterizeReference->GetOutput();
+      m_Reference->UpdateOutputInformation();
       }
 
     // Prepare local streaming
-    
-    RAMDrivenAdaptativeStreamingManagerType::Pointer
-      streamingManager = RAMDrivenAdaptativeStreamingManagerType::New();
+
+
+    m_StreamingManager = RAMDrivenAdaptativeStreamingManagerType::New();
     int availableRAM = GetParameterInt("ram");
-    streamingManager->SetAvailableRAMInMB(availableRAM);
+    m_StreamingManager->SetAvailableRAMInMB( static_cast<unsigned int>( availableRAM ) );
     float bias = 2.0; // empiric value;
-    streamingManager->SetBias(bias);
-    
-    streamingManager->PrepareStreaming(input, input->GetLargestPossibleRegion());
+    m_StreamingManager->SetBias(bias);
+
+    m_StreamingManager->PrepareStreaming(m_Input, m_Input->GetLargestPossibleRegion());
+
+    sid.numberOfStreamDivisions = m_StreamingManager->GetNumberOfSplits();
+
+    otbAppLogINFO("Number of stream divisions : "<<sid.numberOfStreamDivisions);
+
+    return sid;
+  }
 
-    unsigned long numberOfStreamDivisions = streamingManager->GetNumberOfSplits();
+  void DoExecute() ITK_OVERRIDE
+  {
+    StreamingInitializationData sid = InitStreamingData();
 
-    otbAppLogINFO("Number of stream divisions : "<<numberOfStreamDivisions);
+    if(GetParameterString("format") == "contingencytable")
+      {
+      DoExecuteContingencyTable( sid );
+      }
+    else
+      {
+      DoExecuteConfusionMatrix( sid );
+      }
+  }
 
+  void DoExecuteContingencyTable(const StreamingInitializationData& sid)
+  {
+    typedef ContingencyTableCalculator<ClassLabelType> ContingencyTableCalculatorType;
+    ContingencyTableCalculatorType::Pointer calculator = ContingencyTableCalculatorType::New();
+
+    for (unsigned int index = 0; index < sid.numberOfStreamDivisions; index++)
+      {
+      RegionType streamRegion = m_StreamingManager->GetSplit( index );
+
+      m_Input->SetRequestedRegion( streamRegion );
+      m_Input->PropagateRequestedRegion();
+      m_Input->UpdateOutputData();
+
+      m_Reference->SetRequestedRegion( streamRegion );
+      m_Reference->PropagateRequestedRegion();
+      m_Reference->UpdateOutputData();
+
+      ImageIteratorType itInput( m_Input, streamRegion );
+      itInput.GoToBegin();
+
+      ImageIteratorType itRef( m_Reference, streamRegion );
+      itRef.GoToBegin();
+
+      calculator->Compute( itRef, itInput,sid.refhasnodata,sid.refnodata,sid.prodhasnodata,sid.prodnodata);
+      }
+
+    ContingencyTablePointerType contingencyTable = calculator->BuildContingencyTable();
+    LogContingencyTable(contingencyTable);
+    m_WriteContingencyTable(contingencyTable);
+  }
+
+  void DoExecuteConfusionMatrix(const StreamingInitializationData& sid)
+  {
 
     // Extraction of the Class Labels from the Reference image/rasterized vector data + filling of m_Matrix
     MapOfClassesType  mapOfClassesRef, mapOfClassesProd;
@@ -302,22 +412,22 @@ private:
     ClassLabelType labelRef = 0, labelProd = 0;
     int itLabelRef = 0, itLabelProd = 0;
 
-    for (unsigned int index = 0; index < numberOfStreamDivisions; index++)
+    for (unsigned int index = 0; index < sid.numberOfStreamDivisions; index++)
       {
-      RegionType streamRegion = streamingManager->GetSplit(index);
+      RegionType streamRegion = m_StreamingManager->GetSplit(index);
 
-      input->SetRequestedRegion(streamRegion);
-      input->PropagateRequestedRegion();
-      input->UpdateOutputData();
+      m_Input->SetRequestedRegion(streamRegion);
+      m_Input->PropagateRequestedRegion();
+      m_Input->UpdateOutputData();
 
-      reference->SetRequestedRegion(streamRegion);
-      reference->PropagateRequestedRegion();
-      reference->UpdateOutputData();
+      m_Reference->SetRequestedRegion(streamRegion);
+      m_Reference->PropagateRequestedRegion();
+      m_Reference->UpdateOutputData();
 
-      ImageIteratorType itInput(input, streamRegion);
+      ImageIteratorType itInput(m_Input, streamRegion);
       itInput.GoToBegin();
 
-      ImageIteratorType itRef(reference, streamRegion);
+      ImageIteratorType itRef(m_Reference, streamRegion);
       itRef.GoToBegin();
 
       while (!itRef.IsAtEnd())
@@ -326,7 +436,7 @@ private:
         labelProd = static_cast<ClassLabelType> (itInput.Get());
 
         // Extraction of the reference/produced class labels
-        if ((labelRef != nodata) && (labelProd != nodata))
+        if ((!sid.refhasnodata || labelRef != sid.refnodata) && (!sid.prodhasnodata || labelProd != sid.prodnodata))
           {
           // If the current labels have not been added to their respective mapOfClasses yet
           if (mapOfClassesRef.insert(MapOfClassesType::value_type(labelRef, itLabelRef)).second)
@@ -425,8 +535,8 @@ private:
 
 
     // Initialization of the Confusion Matrix for the application LOG and for measurements
-    int nbClassesRef = mapOfClassesRef.size();
-    int nbClassesProd = mapOfClassesProd.size();
+    int nbClassesRef = static_cast<int>(mapOfClassesRef.size());
+    int nbClassesProd = static_cast<int>(mapOfClassesProd.size());
 
     // Formatting m_MatrixLOG from m_Matrix in order to make m_MatrixLOG a square matrix
     // from the reference labels in mapOfClassesRef

Modules/Applications/AppClassification/app/otbTrainVectorClassifier.cxx

@@ -22,6 +22,7 @@
 
 // Validation
 #include "otbConfusionMatrixCalculator.h"
+#include "otbContingencyTableCalculator.h"
 
 namespace otb
 {
@@ -48,6 +49,9 @@ public:
   typedef ConfusionMatrixCalculatorType::ConfusionMatrixType ConfusionMatrixType;
   typedef ConfusionMatrixCalculatorType::MapOfIndicesType MapOfIndicesType;
   typedef ConfusionMatrixCalculatorType::ClassLabelType ClassLabelType;
+  
+  typedef ContingencyTable<ClassLabelType> ContingencyTableType;
+  typedef ContingencyTableType::Pointer    ContingencyTablePointerType;
 
 protected:
   void DoInit()
@@ -72,32 +76,57 @@ protected:
 
   void DoExecute()
   {
-    // Enforce the need of class field name in supervised mode
-    if (GetClassifierCategory() == Supervised)
+    m_FeaturesInfo.SetClassFieldNames( GetChoiceNames( "cfield" ), GetSelectedItems( "cfield" ) );
+
+    if( m_FeaturesInfo.m_SelectedCFieldIdx.empty() && GetClassifierCategory() == Supervised )
       {
-      m_FeaturesInfo.SetClassFieldNames( GetChoiceNames( "cfield" ), GetSelectedItems( "cfield" ) );
+      otbAppLogFATAL( << "No field has been selected for data labelling!" );
+      }
 
-      if( m_FeaturesInfo.m_SelectedCFieldIdx.empty() )
-        {
-        otbAppLogFATAL( << "No field has been selected for data labelling!" );
-        }
+    Superclass::DoExecute();
+
+    if (GetClassifierCategory() == Supervised)
+      {
+      ConfusionMatrixCalculatorType::Pointer confMatCalc = ComputeConfusionMatrix( m_PredictedList,
+                                                                                   m_ClassificationSamplesWithLabel.labeledListSample );
+      WriteConfusionMatrix( confMatCalc );
+      }
+    else
+      {
+      ContingencyTablePointerType table = ComputeContingencyTable( m_PredictedList,
+                                                                   m_ClassificationSamplesWithLabel.labeledListSample );
+      WriteContingencyTable( table );
       }
+  }
 
-      Superclass::DoExecute();
+  ContingencyTablePointerType ComputeContingencyTable(const TargetListSampleType::Pointer &predictedListSample,
+                                                      const TargetListSampleType::Pointer &performanceLabeledListSample)
+  {
+    typedef ContingencyTableCalculator<ClassLabelType> ContigencyTableCalcutaltorType;
 
-      if (GetClassifierCategory() == Supervised)
-        {
-        ConfusionMatrixCalculatorType::Pointer confMatCalc = ComputeConfusionMatrix( m_PredictedList,
-                                                                                     m_ClassificationSamplesWithLabel.labeledListSample );
-        WriteConfusionMatrix( confMatCalc );
-        }
-      else
-        {
-        // TODO Compute Contingency Table
-        }
+    ContigencyTableCalcutaltorType::Pointer contingencyTableCalculator = ContigencyTableCalcutaltorType::New();
+    contingencyTableCalculator->Compute(performanceLabeledListSample->Begin(),
+                                        performanceLabeledListSample->End(),predictedListSample->Begin(), predictedListSample->End());
+
+    otbAppLogINFO( "Training performances:" );
+
+    otbAppLogINFO(<<"Contingency table: reference labels (rows) vs. produced labels (cols)\n"<<contingencyTableCalculator->BuildContingencyTable());
+
+    return contingencyTableCalculator->BuildContingencyTable();
   }
 
 
+  void WriteContingencyTable(const ContingencyTablePointerType& table)
+  {
+    if(IsParameterEnabled("io.confmatout"))
+      {
+      // Write contingency table
+      std::ofstream outFile;
+      outFile.open( this->GetParameterString( "io.confmatout" ).c_str() );
+      outFile << table->ToCSV();
+      }
+  }
+
 
   ConfusionMatrixCalculatorType::Pointer
   ComputeConfusionMatrix(const TargetListSampleType::Pointer &predictedListSample,
@@ -111,7 +140,7 @@ protected:
     confMatCalc->SetProducedLabels( predictedListSample );
     confMatCalc->Compute();
 
-    otbAppLogINFO( "training performances" );
+    otbAppLogINFO( "Training performances:" );
     LogConfusionMatrix( confMatCalc );
 
     for( unsigned int itClasses = 0; itClasses < confMatCalc->GetNumberOfClasses(); itClasses++ )

Modules/Applications/AppClassification/include/otbLearningApplicationBase.txx

@@ -33,6 +33,7 @@ namespace Wrapper
 template <class TInputValue, class TOutputValue>
 LearningApplicationBase<TInputValue,TOutputValue>
 ::LearningApplicationBase() : m_RegressionFlag(false)
+
 {
 }
 
@@ -59,7 +60,9 @@ LearningApplicationBase<TInputValue,TOutputValue>
 
   InitUnsupervisedClassifierParams();
   std::vector<std::string> allClassifier = GetChoiceKeys("classifier");
-  m_UnsupervisedClassifier.assign(allClassifier.begin() + m_SupervisedClassifier.size(), allClassifier.end());
+  // Check for empty unsupervised classifier
+  if( allClassifier.size() > m_UnsupervisedClassifier.size() )
+    m_UnsupervisedClassifier.assign( allClassifier.begin() + m_SupervisedClassifier.size(), allClassifier.end() );
 }
 
 template <class TInputValue, class TOutputValue>
@@ -67,10 +70,16 @@ typename LearningApplicationBase<TInputValue,TOutputValue>::ClassifierCategory
 LearningApplicationBase<TInputValue,TOutputValue>
 ::GetClassifierCategory()
 {
-  bool foundUnsupervised =
-          std::find(m_UnsupervisedClassifier.begin(), m_UnsupervisedClassifier.end(),
-                    GetParameterString("classifier")) != m_UnsupervisedClassifier.end();
-  return foundUnsupervised ? Unsupervised : Supervised;
+  if( m_UnsupervisedClassifier.empty() )
+    {
+    return Supervised;
+    }
+  else
+    {
+    bool foundUnsupervised = std::find( m_UnsupervisedClassifier.begin(), m_UnsupervisedClassifier.end(),
+                                        GetParameterString( "classifier" ) ) != m_UnsupervisedClassifier.end();
+    return foundUnsupervised ? Unsupervised : Supervised;
+    }
 }
 
 template <class TInputValue, class TOutputValue>

Modules/Applications/AppClassification/include/otbTrainVectorBase.h

@@ -125,7 +125,7 @@ protected:
     {
       m_SelectedCFieldIdx = selectedCFieldIdx;
       // Handle only one class field name, if several are provided only the first one is used.
-      m_SelectedCFieldName = cFieldNames[selectedCFieldIdx.front()];
+      m_SelectedCFieldName = selectedCFieldIdx.empty() ? cFieldNames.front() : cFieldNames[selectedCFieldIdx.front()];
     }
   };
 

Modules/Applications/AppClassification/test/CMakeLists.txt

@@ -392,6 +392,7 @@ otb_test_application(NAME apTvComputeConfusionMatrixExtraReferenceLabelsR
   OPTIONS -in ${INPUTDATA}/Classification/QB_1_ortho_C7.tif
   -ref raster
   -ref.raster.in ${INPUTDATA}/Classification/clLabeledImageQB456_1_NoData_255.tif
+  -ref.raster.nodata 255
   -nodatalabel 255
   -out ${TEMP}/apTvComputeConfusionMatrixExtraRefLabelsROut.csv
   VALID   --compare-ascii ${NOTOL}
@@ -403,6 +404,7 @@ otb_test_application(NAME apTvComputeConfusionMatrixR
   OPTIONS -in ${OTBAPP_BASELINE}/clLabeledImageQB123_1.tif
   -ref raster
   -ref.raster.in ${INPUTDATA}/Classification/clLabeledImageQB456_1_NoData_255.tif
+  -ref.raster.nodata 255
   -nodatalabel 255
   -out ${TEMP}/apTvComputeConfusionMatrixTconfusionROut.csv
   VALID   --compare-ascii ${NOTOL}
@@ -414,12 +416,27 @@ otb_test_application(NAME apTvComputeConfusionMatrixExtraProducedLabelsR
   OPTIONS -in ${INPUTDATA}/Classification/clLabeledImageQB456_1_NoData_255.tif
   -ref raster
   -ref.raster.in ${INPUTDATA}/Classification/QB_1_ortho_C8.tif
+  -ref.raster.nodata 255
   -nodatalabel 255
   -out ${TEMP}/apTvComputeConfusionMatrixExtraProdLabelsROut.csv
   VALID   --compare-ascii ${NOTOL}
   ${OTBAPP_BASELINE_FILES}/apTvComputeConfusionMatrixExtraProdLabelsROut.csv
   ${TEMP}/apTvComputeConfusionMatrixExtraProdLabelsROut.csv)
 
+#----------- ComputeContingencyTable TESTS ----------------
+
+otb_test_application(NAME apTvComputeContingencyTableExtraProducedLabelsR
+  APP  ComputeConfusionMatrix
+  OPTIONS -in ${INPUTDATA}/Classification/clLabeledImageQB456_1_NoData_255.tif
+  -ref raster
+  -ref.raster.in ${INPUTDATA}/Classification/QB_1_ortho_C8.tif
+  -ref.raster.nodata 255
+  -nodatalabel 255
+  -format contingencytable
+  -out ${TEMP}/apTvComputeContingencyTableExtraProdLabelsROut.csv
+  VALID   --compare-ascii ${NOTOL}
+  ${OTBAPP_BASELINE_FILES}/apTvComputeContingencyTableExtraProdLabelsROut.csv
+  ${TEMP}/apTvComputeContingencyTableExtraProdLabelsROut.csv)
 
 #----------- FusionOfClassifications TESTS ----------------
 otb_test_application(NAME apTvFusionOfClassificationsDSPrecision6Inputs

Modules/Learning/Unsupervised/include/otbContingencyTable.h

+/*
+ * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part of Orfeo Toolbox
+ *
+ *     https://www.orfeo-toolbox.org/
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef otbContingencyTable_h
+#define otbContingencyTable_h
+
+#include <vector>
+#include <iostream>
+#include <iomanip>
+#include <itkObject.h>
+#include <itkObjectFactory.h>
+#include <itkVariableSizeMatrix.h>
+
+namespace otb
+{
+template<class TClassLabel>
+class ContingencyTable : public itk::Object
+{
+public:
+  /** Standard class typedefs */
+  typedef ContingencyTable              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 itk::VariableSizeMatrix<unsigned long> MatrixType;
+  typedef std::vector<TClassLabel>               LabelList;
+
+  MatrixType matrix;
+
+  void SetLabels(const LabelList &referenceLabels, const LabelList &producedLabels)
+  {
+    m_RefLabels = referenceLabels;
+    m_ProdLabels = producedLabels;
+    unsigned int rows = static_cast<unsigned int>(m_RefLabels.size());
+    unsigned int cols = static_cast<unsigned int>(m_ProdLabels.size());
+    matrix.SetSize( rows, cols );
+    matrix.Fill( 0 );
+  }
+
+  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;
+    maxWidth = GetLabelsMaximumLength(contingencyTable.m_ProdLabels, maxWidth);
+    maxWidth = GetLabelsMaximumLength(contingencyTable.m_RefLabels, maxWidth);
+
+    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)+1;
+
+    // Write the first line of the matrix (produced labels)
+    o << std::setfill(' ') << std::setw( width ) << "labels";
+    for( size_t i = 0; i < contingencyTable.m_ProdLabels.size(); ++i )
+      {
+      o << std::setfill(' ') << std::setw( width ) << contingencyTable.m_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.m_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;