Skip to content
Snippets Groups Projects
Commit f7ae6e7d authored by Julien Malik's avatar Julien Malik
Browse files

ENH: split Hyperspectral application in two application : VCA and Unmixing

parent 2cceab98
Branches
Tags
No related merge requests found
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)
......@@ -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;
......
/*=========================================================================
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)
......@@ -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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment