From f7ae6e7ddd70ac6de5b31c7fab4d517d30eb02a6 Mon Sep 17 00:00:00 2001
From: Julien Malik <julien.malik@c-s.fr>
Date: Thu, 3 Nov 2011 19:22:37 +0100
Subject: [PATCH] ENH: split Hyperspectral application in two application : VCA
 and Unmixing

---
 Applications/Hyperspectral/CMakeLists.txt     |   1 +
 .../otbHyperspectralUnmixing.cxx              | 255 ++----------------
 .../otbVertexComponentAnalysis.cxx            | 123 +++++++++
 .../Applications/Hyperspectral/CMakeLists.txt |  35 +--
 4 files changed, 158 insertions(+), 256 deletions(-)
 create mode 100644 Applications/Hyperspectral/otbVertexComponentAnalysis.cxx

diff --git a/Applications/Hyperspectral/CMakeLists.txt b/Applications/Hyperspectral/CMakeLists.txt
index 5d19ce017b..dae81dfbd1 100644
--- a/Applications/Hyperspectral/CMakeLists.txt
+++ b/Applications/Hyperspectral/CMakeLists.txt
@@ -1,2 +1,3 @@
 
 OTB_CREATE_APPLICATION(NAME HyperspectralUnmixing SOURCES otbHyperspectralUnmixing.cxx LINK_LIBRARIES OTBIO;OTBCommon;OTBBasicFilters)
+OTB_CREATE_APPLICATION(NAME VertexComponentAnalysis SOURCES otbVertexComponentAnalysis.cxx LINK_LIBRARIES OTBIO;OTBCommon;OTBBasicFilters)
diff --git a/Applications/Hyperspectral/otbHyperspectralUnmixing.cxx b/Applications/Hyperspectral/otbHyperspectralUnmixing.cxx
index 0042ffe3d8..096281e00a 100644
--- a/Applications/Hyperspectral/otbHyperspectralUnmixing.cxx
+++ b/Applications/Hyperspectral/otbHyperspectralUnmixing.cxx
@@ -42,10 +42,10 @@ typedef otb::StreamingStatisticsVectorImageFilter<DoubleVectorImageType> Streami
 typedef otb::EigenvalueLikelihoodMaximisation<double> ELMType;
 typedef otb::VCAImageFilter<DoubleVectorImageType> VCAFilterType;
 
-typedef otb::UnConstrainedLeastSquareImageFilter<DoubleVectorImageType, DoubleVectorImageType, double> UCLSUnmixingFilterType;
-typedef otb::ISRAUnmixingImageFilter<DoubleVectorImageType, DoubleVectorImageType, double>             ISRAUnmixingFilterType;
-typedef otb::NCLSUnmixingImageFilter<DoubleVectorImageType, DoubleVectorImageType, double>             NCLSUnmixingFilterType;
-typedef otb::FCLSUnmixingImageFilter<DoubleVectorImageType, DoubleVectorImageType, double>             FCLSUnmixingFilterType;
+typedef otb::UnConstrainedLeastSquareImageFilter<DoubleVectorImageType,DoubleVectorImageType,double> UCLSUnmixingFilterType;
+typedef otb::ISRAUnmixingImageFilter<DoubleVectorImageType,DoubleVectorImageType,double>             ISRAUnmixingFilterType;
+typedef otb::NCLSUnmixingImageFilter<DoubleVectorImageType,DoubleVectorImageType,double>             NCLSUnmixingFilterType;
+typedef otb::FCLSUnmixingImageFilter<DoubleVectorImageType,DoubleVectorImageType,double>             FCLSUnmixingFilterType;
 
 typedef otb::VectorImageToMatrixImageFilter<DoubleVectorImageType> VectorImageToMatrixImageFilterType;
 
@@ -72,14 +72,13 @@ enum EndmembersEstimationMethod
 
 enum UnMixingMethod
 {
-  UnMixingMethod_NONE,
   UnMixingMethod_UCLS,
   UnMixingMethod_ISRA,
   UnMixingMethod_NCLS,
   UnMixingMethod_FCLS,
 };
 
-const char* UnMixingMethodNames [] = { "NONE", "UCLS", "ISRA", "NCLS", "FCLS", };
+const char* UnMixingMethodNames [] = { "UCLS", "ISRA", "NCLS", "FCLS", };
 
 
 class HyperspectralUnmixing : public Application
@@ -104,13 +103,12 @@ private:
     SetDescription("Unmix an hyperspectral image");
 
     // Documentation
-    SetDocName("Hyper spectral unmixing application");
-    SetDocLongDescription("This application \"unmix\" an hyperspectral image.");
+    SetDocName("Hyperspectral data unmixing");
+    SetDocLongDescription("Applies an unmixing algorithm to an hyperspectral data cube");
     SetDocLimitations("None");
     SetDocAuthors("OTB-Team");
     SetDocSeeAlso(" ");
-    SetDocCLExample(" ");
-    AddDocTag("Image manipulation");
+    SetDocCLExample("otbApplicationLauncherCommandLine VertexComponentAnalysis ${OTB-BIN}/bin --in ${INPUTDATA}/Hyperspectral/synthetic/hsi_cube.tif --ie ${INPUTDATA}/Hyperspectral/synthetic/endmembers.tif --out ${TEMP}/apTvHyHyperspectralUnmixing_UCLS.tif double --ua ucls");
     AddDocTag("Hyperspectral");
   }
 
@@ -121,45 +119,31 @@ private:
   void DoCreateParameters()
   {
     AddParameter(ParameterType_InputImage,  "in",   "Input Image");
-    AddParameter(ParameterType_OutputImage, "out",  "Output Image");
-    MandatoryOff("out");
-
-    AddParameter(ParameterType_Choice, "dr", "Dimension reduction");
-    AddChoice("dr.pca", "PCA");
-    AddChoice("dr.mnf", "MNF");
-    AddChoice("dr.napca", "NAPCA");
-    MandatoryOff("dr");
-
-    AddParameter(ParameterType_Int, "ne", "Number of endmembers");
-    SetDefaultParameterInt("ne", 1);
-    MandatoryOff("ne");
-
-    AddParameter(ParameterType_Choice, "de", "Dimensionnality estimation");
-    MandatoryOff("de");
-    AddChoice("de.elm", "ELM (Eigen Value Likelihood Estimation)");
+    SetParameterDescription("in","The hyperspectral data cube to unmix");
 
-    AddParameter(ParameterType_Choice, "ee", "Endmembers estimation method");
-    MandatoryOff("ee");
-    AddChoice("ee.vca", "VCA (Vertex Component Analysis)");
+    AddParameter(ParameterType_OutputImage, "out",  "Output Image");
+    SetParameterDescription("out","The output abundance map");
 
-    AddParameter(ParameterType_InputImage, "ie", "Input endmembers image");
-    MandatoryOff("ie");
+    AddParameter(ParameterType_InputImage,  "ie",   "Input endmembers");
+    SetParameterDescription("ie","The endmembers to use for unmixing. Must be stored as a multispectral image, where each pixel is interpreted as an endmember");
 
     AddParameter(ParameterType_Choice, "ua", "Unmixing algorithm");
+    SetParameterDescription("ua", "The algorithm to use for unmixing");
     MandatoryOff("ua");
-    AddChoice("ua.none", "None");
     AddChoice("ua.ucls", "UCLS");
+    SetParameterDescription("ua.ucls", "Unconstrained Least Square");
+
     AddChoice("ua.isra", "ISRA");
+    SetParameterDescription("ua.isra", "Image Space Reconstruction Algorithm");
+
     AddChoice("ua.ncls", "NCLS");
+    SetParameterDescription("ua.ncls", "Non-negative constrained Least Square");
+
     AddChoice("ua.fcls", "FCLS");
-    //SetParameterString("ua", "none");
+    SetParameterDescription("ua.fcls", "Fully constrained Least Square");
 
-    AddParameter(ParameterType_OutputImage, "oe", "Output Endmembers");
-    MandatoryOff("oe");
+    SetParameterString("ua", "ucls");
 
-    AddParameter(ParameterType_RAM, "ram", "Available RAM");
-    SetDefaultParameterInt("ram", 256);
-    MandatoryOff("ram");
   }
 
   void DoUpdateParameters()
@@ -171,177 +155,8 @@ private:
   {
     m_ProcessObjects.clear();
 
-    /*
-     *
-     * EXTRACT PARAMETERS
-     *
-     */
-    // Dimension Reduction
-    DimReductionMethod dimReduction = DimReductionMethod_NONE;
-    if ( IsParameterEnabled("dr") && HasValue("dr") )
-      {
-      std::string dimReductionStr = GetParameterString("dr");
-      if ( boost::to_upper_copy(dimReductionStr) == "PCA" )
-        {
-        dimReduction = DimReductionMethod_PCA;
-        }
-      else if ( boost::to_upper_copy(dimReductionStr) == "MNF" )
-        {
-        dimReduction = DimReductionMethod_MNF;
-        }
-      }
-    otbAppLogDEBUG(<< "Using "
-                   << (dimReduction == DimReductionMethod_NONE ? "NONE" : (dimReduction == DimReductionMethod_PCA ? "PCA" : "MNF") )
-                   << " dimensionality reduction method")
-
-    // Number of endmembers
-    unsigned int nbEndmembers = 0;
-    if ( IsParameterEnabled("ne") && HasValue("ne") )
-      {
-      nbEndmembers = GetParameterInt("ne");
-      otbAppLogDEBUG(<< "nbEndmembers: " << nbEndmembers)
-      }
-
-    // Dimensionnality estimation
-    if ( IsParameterEnabled("de") && HasValue("de") )
-      {
-      std::string dimEstMethodStr = GetParameterString("de");
-      if (boost::to_upper_copy(dimEstMethodStr) != "ELM")
-        {
-        std::cerr << "Only ELM is supported for parameter DimensionalityEstimationMethod" << std::endl;
-        }
-      otbAppLogDEBUG(<< "dimEstMethodStr: " << dimEstMethodStr)
-
-      }
-
-    if ( IsParameterEnabled("ee") && HasValue("ee") )
-      {
-      std::string eeEstMethodStr = GetParameterString("ee");
-      if (boost::to_upper_copy(eeEstMethodStr) != "VCA")
-        {
-        std::cerr << "Only VCA is supported for parameter EndmembersEstimationMethod" << std::endl;
-        }
-      otbAppLogDEBUG(<< "eeEstMethodStr: " << eeEstMethodStr)
-      }
-
-    std::string inputEndmembers;
-    if ( IsParameterEnabled("ie") && HasValue("ie") )
-      {
-      inputEndmembers = GetParameterString("ie");
-      otbAppLogDEBUG("inputEndmembers: " << inputEndmembers)
-      }
-
-    UnMixingMethod unmixingAlgo = UnMixingMethod_NONE;
-    if ( IsParameterEnabled("ua") && HasValue("ua") )
-      {
-      std::string unmixingAlgoStr = GetParameterString("ua");
-
-      if ( boost::to_upper_copy(unmixingAlgoStr) == "UCLS" )
-        {
-        unmixingAlgo = UnMixingMethod_UCLS;
-        }
-      else if ( boost::to_upper_copy(unmixingAlgoStr) == "NCLS" )
-        {
-        unmixingAlgo = UnMixingMethod_NCLS;
-        }
-      else if ( boost::to_upper_copy(unmixingAlgoStr) == "FCLS" )
-        {
-        unmixingAlgo = UnMixingMethod_FCLS;
-        }
-      else if ( boost::to_upper_copy(unmixingAlgoStr) == "ISRA" )
-        {
-        unmixingAlgo = UnMixingMethod_ISRA;
-        }
-      }
-    otbAppLogDEBUG( << "Using "
-                    << UnMixingMethodNames[unmixingAlgo]
-                    << " unmixing algorithm")
-
-    std::string outputEndmembers;
-    if ( IsParameterEnabled("oe") && HasValue("oe") && inputEndmembers.empty() )
-      {
-      outputEndmembers = GetParameterString("oe");
-      }
-    otbAppLogDEBUG( << "IsParameterEnabled(oe) " << IsParameterEnabled("oe") )
-    otbAppLogDEBUG( << "HasValue(oe) " << HasValue("oe") )
-    otbAppLogDEBUG( << "inputEndmembers.empty() " << inputEndmembers.empty() )
-    otbAppLogDEBUG( << "outputEndmembers: " << outputEndmembers )
-
-
-    /*
-     *
-     * PROCESSING
-     *
-     */
-
     DoubleVectorImageType::Pointer inputImage = GetParameterImage<DoubleVectorImageType>("in");
-
-    DoubleVectorImageType::Pointer endmembersImage;
-    if ( inputEndmembers.empty() )
-      {
-      if( nbEndmembers == 0 )
-        {
-        /*
-         * Compute stats of input image
-         */
-
-        std::cout << "Computing stats" << std::endl;
-
-        StreamingStatisticsVectorImageFilterType::Pointer stats = StreamingStatisticsVectorImageFilterType::New();
-
-        stats->SetInput(inputImage);
-        otb::StandardWriterWatcher watcher(stats->GetStreamer(), stats->GetFilter(), "Computing image stats");
-        stats->Update();
-
-        VectorType mean (stats->GetMean().GetDataPointer(), stats->GetMean().GetSize());
-        MatrixType covariance  = stats->GetCovariance().GetVnlMatrix();
-        MatrixType correlation = stats->GetCorrelation().GetVnlMatrix();
-
-        /*
-         * Estimate Endmembers Numbers
-         */
-        std::cout << "Estimate Endmembers Numbers by ELM" << std::endl;
-
-        ELMType::Pointer elm = ELMType::New();
-        elm->SetCovariance(covariance);
-        elm->SetCorrelation(correlation);
-        elm->SetNumberOfPixels(inputImage->GetLargestPossibleRegion().GetNumberOfPixels());
-        elm->Compute();
-
-        nbEndmembers = elm->GetNumberOfEndmembers();
-
-        std::cout << "ELM : " << nbEndmembers << " estimated endmembers" << std::endl;
-
-        std::cout << "ELM likelihood values: " << elm->GetLikelihood() << std::endl;
-
-        }
-      else
-        {
-        std::cout << "Using " << nbEndmembers << " endmembers" << std::endl;
-        }
-
-      /*
-       * Estimate Endmembers
-       */
-      std::cout << "Estimate Endmembers by VCA" << std::endl;
-
-      VCAFilterType::Pointer vca = VCAFilterType::New();
-      vca->SetNumberOfEndmembers(nbEndmembers);
-      vca->SetInput(inputImage);
-
-      endmembersImage = vca->GetOutput();
-      m_ProcessObjects.push_back(vca.GetPointer());
-      }
-    else
-      {
-      /*
-       * Read input endmembers
-       */
-      std::cout << "Read Endmembers " << inputEndmembers << std::endl;
-      endmembersImage = GetParameterImage<DoubleVectorImageType>("ie");
-
-      }
-  //  endmembersRef->Update();
+    DoubleVectorImageType::Pointer endmembersImage = GetParameterImage<DoubleVectorImageType>("ie");
 
     /*
      * Transform Endmembers image to matrix representation
@@ -359,10 +174,9 @@ private:
      * Unmix
      */
     DoubleVectorImageType::Pointer abundanceMap;
-    switch (unmixingAlgo)
+
+    switch ( static_cast<UnMixingMethod>(GetParameterInt("ua")) )
     {
-    case UnMixingMethod_NONE:
-      break;
     case UnMixingMethod_UCLS:
       {
       std::cout << "UCLS Unmixing" << std::endl;
@@ -425,25 +239,8 @@ private:
       break;
     }
 
-    if ( !outputEndmembers.empty() )
-      {
-      /*
-       * Write endmembers
-       */
-      std::cout << "Write endmembers " << outputEndmembers << std::endl;
-
-      SetParameterOutputImage<DoubleVectorImageType>("oe", endmembersImage);
-      }
+    SetParameterOutputImage<DoubleVectorImageType>("out", abundanceMap);
 
-    if ( unmixingAlgo != UnMixingMethod_NONE )
-      {
-      /*
-       * Write abundance map
-       */
-      //std::cout << "Write abundance map" << outputImageName << std::endl;
-
-      SetParameterOutputImage<DoubleVectorImageType>("out", abundanceMap);
-      }
   }
 
   std::vector<itk::ProcessObject::Pointer> m_ProcessObjects;
diff --git a/Applications/Hyperspectral/otbVertexComponentAnalysis.cxx b/Applications/Hyperspectral/otbVertexComponentAnalysis.cxx
new file mode 100644
index 0000000000..a8e24171b6
--- /dev/null
+++ b/Applications/Hyperspectral/otbVertexComponentAnalysis.cxx
@@ -0,0 +1,123 @@
+/*=========================================================================
+
+  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 "otbWrapperApplication.h"
+#include "otbWrapperApplicationFactory.h"
+
+#include <boost/algorithm/string.hpp>
+
+#include "otbStreamingStatisticsVectorImageFilter.h"
+#include "otbEigenvalueLikelihoodMaximisation.h"
+#include "otbVcaImageFilter.h"
+#include "otbUnConstrainedLeastSquareImageFilter.h"
+#include "otbISRAUnmixingImageFilter.h"
+#include "otbNCLSUnmixingImageFilter.h"
+#include "otbFCLSUnmixingImageFilter.h"
+
+#include "otbVectorImageToMatrixImageFilter.h"
+
+
+namespace otb
+{
+namespace Wrapper
+{
+
+const unsigned int Dimension = 2;
+
+typedef otb::VCAImageFilter<DoubleVectorImageType> VCAFilterType;
+typedef otb::VectorImageToMatrixImageFilter<DoubleVectorImageType> VectorImageToMatrixImageFilterType;
+
+
+class VertexComponentAnalysis : public Application
+{
+public:
+  /** Standard class typedefs. */
+  typedef VertexComponentAnalysis         Self;
+  typedef Application                   Superclass;
+  typedef itk::SmartPointer<Self>       Pointer;
+  typedef itk::SmartPointer<const Self> ConstPointer;
+
+  /** Standard macro */
+  itkNewMacro(Self);
+
+  itkTypeMacro(VertexComponentAnalysis, otb::Application);
+
+private:
+
+  VertexComponentAnalysis()
+  {
+    SetName("VertexComponentAnalysis");
+    SetDescription("Find endmembers in hyperspectral images with Vertex Component Analysis");
+
+    // Documentation
+    SetDocName("Vertex Component Analysis");
+    SetDocLongDescription("Applies the Vertex Component Analysis to an hyperspectral image to extract endmembers");
+    SetDocLimitations("None");
+    SetDocAuthors("OTB-Team");
+    SetDocSeeAlso(" ");
+    SetDocCLExample("otbApplicationLauncherCommandLine VertexComponentAnalysis ${OTB-BIN}/bin --in ${OTB_DATA_ROOT}/Input/Hyperspectral/synthetic/hsi_cube.tif --ne 5 --outendm ${TEMP}/apTvHyHyperspectralUnmixing_VCA.tif double");
+    AddDocTag("Hyperspectral");
+
+  }
+
+  virtual ~VertexComponentAnalysis()
+  {
+  }
+
+  void DoCreateParameters()
+  {
+    AddParameter(ParameterType_InputImage,  "in",   "Input Image");
+    SetParameterDescription("out","Input hyperspectral data cube");
+
+    AddParameter(ParameterType_Int, "ne", "Number of endmembers");
+    SetParameterDescription("ne","The number of endmembers to extract from the data cube");
+    SetParameterInt("ne", 1);
+    MandatoryOn("ne");
+
+    AddParameter(ParameterType_OutputImage, "outendm", "Output Endmembers");
+    SetParameterDescription("outendm","The endmebers, stored in a one-line multi-spectral image, each pixel representing an endmember");
+    MandatoryOn("outendm");
+  }
+
+  void DoUpdateParameters()
+  {
+    // Nothing to do here : all parameters are independent
+  }
+
+  void DoExecute()
+  {
+    DoubleVectorImageType::Pointer inputImage = GetParameterImage<DoubleVectorImageType>("in");
+    DoubleVectorImageType::Pointer endmembersImage;
+
+    const unsigned int nbEndmembers = GetParameterInt("ne");
+    VCAFilterType::Pointer vca = VCAFilterType::New();
+    vca->SetNumberOfEndmembers(nbEndmembers);
+    vca->SetInput(inputImage);
+
+    endmembersImage = vca->GetOutput();
+    m_Ref = vca.GetPointer();
+
+    SetParameterOutputImage<DoubleVectorImageType>("outendm", endmembersImage);
+  }
+
+  itk::LightObject::Pointer m_Ref;
+};
+
+}
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::VertexComponentAnalysis)
diff --git a/Testing/Applications/Hyperspectral/CMakeLists.txt b/Testing/Applications/Hyperspectral/CMakeLists.txt
index cdb0d3c0a2..5893421db1 100644
--- a/Testing/Applications/Hyperspectral/CMakeLists.txt
+++ b/Testing/Applications/Hyperspectral/CMakeLists.txt
@@ -52,37 +52,18 @@ add_test(NAME apTvHyHyperspectralUnmixing_FCLS
                  --ie ${INPUTDATA}/Hyperspectral/synthetic/endmembers.tif
                  --out ${TEMP}/apTvHyHyperspectralUnmixing_FCLS.tif double
                  --ua fcls )
-                 
-add_test(NAME apTvHyHyperspectralUnmixing_VCA
-         COMMAND otbTestDriver
-                 --compare-image ${EPSILON_9}
-                   ${BASELINE}/apTvHyHyperspectralUnmixing_VCA.tif
-                   ${TEMP}/apTvHyHyperspectralUnmixing_VCA.tif
-                 Execute $<TARGET_FILE:otbApplicationLauncherCommandLine>
-                 HyperspectralUnmixing
-                 $<TARGET_FILE_DIR:otbapp_HyperspectralUnmixing>
-                 --in ${INPUTDATA}/Hyperspectral/synthetic/hsi_cube.tif
-                 --ne 5 
-                 --oe ${TEMP}/apTvHyHyperspectralUnmixing_VCA.tif double )
 
-add_test(NAME apTvHyHyperspectralUnmixing_VCA_UCLS
-         COMMAND otbTestDriver
-                 Execute $<TARGET_FILE:otbApplicationLauncherCommandLine>
-                 HyperspectralUnmixing
-                 $<TARGET_FILE_DIR:otbapp_HyperspectralUnmixing>
-                 --in ${OTB_DATA_ROOT}/Input/Hyperspectral/synthetic/hsi_cube.tif
-                 --out ${TEMP}/apTvHyHyperspectralUnmixing_ucls.tif double
-                 --ne 5
-                 --ua ucls )
+#--- VertexComponentAnalysis ---#
 
-add_test(NAME apTvHyHyperspectralUnmixing_VCA_FCLS
+add_test(NAME apTvHyVertexComponentAnalysis
          COMMAND otbTestDriver
                  Execute $<TARGET_FILE:otbApplicationLauncherCommandLine>
-                 HyperspectralUnmixing
-                 $<TARGET_FILE_DIR:otbapp_HyperspectralUnmixing>
+                 VertexComponentAnalysis
+                 $<TARGET_FILE_DIR:otbapp_VertexComponentAnalysis>
                  --in ${OTB_DATA_ROOT}/Input/Hyperspectral/synthetic/hsi_cube.tif
-                 --out ${TEMP}/apTvHyHyperspectralUnmixing_fcls.tif double
                  --ne 5
-                 --ua fcls )
+                 --outendm ${TEMP}/apTvHyHyperspectralUnmixing_VCA.tif double )
+
+
 
-    
\ No newline at end of file
+    
-- 
GitLab