diff --git a/Application/CMakeLists.txt b/Application/CMakeLists.txt
index 7d3117c16683e7bc28e99b408e022834e33746f6..0ae28d26fe4b492b3fd3fb097630a50c1befb728 100644
--- a/Application/CMakeLists.txt
+++ b/Application/CMakeLists.txt
@@ -1,3 +1,32 @@
ADD_EXECUTABLE(otbAvirisBandSelection otbAvirisBandSelection.cxx)
TARGET_LINK_LIBRARIES(otbAvirisBandSelection OTBCommon OTBIO OTBBasicFilters )
+
+ADD_EXECUTABLE(otbHyperspectralUnmixing otbHyperspectralUnmixing.cxx)
+TARGET_LINK_LIBRARIES(otbHyperspectralUnmixing OTBCommon OTBIO OTBBasicFilters )
+
+
+
+
+SET(TESTEXE_DIR ${Hyperspectral_BINARY_DIR}/bin)
+SET(TEMP ${Hyperspectral_BINARY_DIR}/Testing/Temporary)
+SET(DATA ${Hyperspectral_SOURCE_DIR}/Data)
+
+# parser->AddInputImage();
+# parser->AddOption( "--DimReduceMethod", "Dimensionnality reduction method (NONE,PCA,MCF)", "-dr", 1, false );
+# parser->AddOption( "--NumEndmembers", "Number of endmembers (if not given or 0, will be estimated with DimensionalityEstimationMethod)", "-ne", 1, false );
+# parser->AddOption( "--DimensionalityEstimationMethod", "Dimensionality estimation method (ELM)", "-de", 1, false );
+# parser->AddOption( "--EndmembersEstimationMethod", "Endmembers estimation method (VCA)", "-ee", 1, false );
+# parser->AddOption( "--InputEndmembers", "Input endmembers image (must have NumComponents components or fit the input image nb of components)", "-ie", 1, false );
+# parser->AddOption( "--UnmixingAlgorithm", "Unmixing algorithm (NONE, UCLS, ISRA)", "-ua", 1, false );
+# parser->AddOption( "--OutputEndmembers", "Output estimated endmembers image", "-oe", 1, false );
+# parser->AddOutputImage();
+
+ADD_TEST(hyTvUnmixApplication
+ ${TESTEXE_DIR}/otbHyperspectralUnmixing
+ otbHyperspectralUnmixing
+ -in ${DATA}/AVIRIS/extract128/extract128_goodbands.hdr
+
+ 8)
+
+
\ No newline at end of file
diff --git a/Application/otbHyperspectralUnmixing.cxx b/Application/otbHyperspectralUnmixing.cxx
index 9168b17500af2b7751c5c873ed2e7fe5e40d0c31..7e43add3a27bd3d6dd872b4d83f1525e696ac38f 100644
--- a/Application/otbHyperspectralUnmixing.cxx
+++ b/Application/otbHyperspectralUnmixing.cxx
@@ -27,6 +27,8 @@
#include "otbVcaImageFilter.h"
#include "otbUnConstrainedLeastSquareImageFilter.h"
#include "otbISRAUnmixingImageFilter.h"
+#include "otbUnConstrainedLeastSquareImageFilter.h"
+#include "otbVectorImageToMatrixImageFilter.h"
const unsigned int Dimension = 2;
typedef float PixelType;
@@ -36,25 +38,51 @@ typedef otb::VectorImage<PixelType, Dimension> VectorImageType;
typedef otb::ImageFileReader<VectorImageType> ReaderType;
typedef otb::StreamingImageFileWriter<VectorImageType> WriterType;
-typedef otb::StreamingStatisticsVectorImageFilter2<VectorImageType> StreamingStatisticsVectorImageFilterType;
+typedef otb::StreamingStatisticsVectorImageFilter2<VectorImageType, PixelType> StreamingStatisticsVectorImageFilterType;
typedef otb::EigenvalueLikelihoodMaximisation<PixelType> ELMType;
typedef otb::VCAImageFilter<VectorImageType> VCAFilterType;
-typedef enum DimReduction
+typedef otb::UnConstrainedLeastSquareImageFilter<VectorImageType,VectorImageType,PixelType> UCLSUnmixingFilterType;
+typedef otb::ISRAUnmixingImageFilter<VectorImageType,VectorImageType,PixelType> ISRAUnmixingFilterType;
+
+typedef otb::VectorImageToMatrixImageFilter<VectorImageType> VectorImageToMatrixImageFilterType;
+
+typedef vnl_vector<PixelType> VectorType;
+typedef vnl_matrix<PixelType> MatrixType;
+
+
+
+enum DimReductionMethod
+{
+ DimReductionMethod_NONE,
+ DimReductionMethod_PCA,
+ DimReductionMethod_MNF
+};
+
+enum DimensionalityEstimationMethod
+{
+ DimensionalityEstimationMethod_ELM
+};
+
+enum EndmembersEstimationMethod
{
- DimReduction_NONE,
- DimReduction_PCA,
- DimReduction_MNF
+ EndmembersEstimationMethod_VCA
};
-typedef enum UnMixingAlgo
+enum UnMixingMethod
{
- UnMixingAlgo_UCLS,
- UnMixingAlgo_ISRA,
+ UnMixingMethod_NONE,
+ UnMixingMethod_UCLS,
+ UnMixingMethod_ISRA
};
int main(int argc, char * argv[])
{
+ /*
+ *
+ * DEFINE APPLICATION
+ *
+ */
typedef otb::CommandLineArgumentParser ParserType;
ParserType::Pointer parser = ParserType::New();
@@ -65,17 +93,15 @@ int main(int argc, char * argv[])
parser->AddOption( "--DimensionalityEstimationMethod", "Dimensionality estimation method (ELM)", "-de", 1, false );
parser->AddOption( "--EndmembersEstimationMethod", "Endmembers estimation method (VCA)", "-ee", 1, false );
parser->AddOption( "--InputEndmembers", "Input endmembers image (must have NumComponents components or fit the input image nb of components)", "-ie", 1, false );
- parser->AddOption( "--UnmixingAlgorithm", "Unmixing algorithm (UCLS, ISRA)", "-ua", 1, false );
+ parser->AddOption( "--UnmixingAlgorithm", "Unmixing algorithm (NONE, UCLS, ISRA)", "-ua", 1, false );
parser->AddOption( "--OutputEndmembers", "Output estimated endmembers image", "-oe", 1, false );
parser->AddOutputImage();
-
/*
*
* EXTRACT PARAMETERS
*
*/
-
typedef otb::CommandLineArgumentParseResult ParserResultType;
ParserResultType::Pointer parseResult = ParserResultType::New();
@@ -85,7 +111,6 @@ int main(int argc, char * argv[])
}
catch( itk::ExceptionObject & err )
{
- //std::cerr << argv[0] << " applies MNF transformations\n";
std::string descriptionException = err.GetDescription();
if ( descriptionException.find("ParseCommandLine(): Help Parser")
!= std::string::npos )
@@ -99,23 +124,22 @@ int main(int argc, char * argv[])
const char * inputImageName = parseResult->GetInputImage().c_str();
const char * outputImageName = parseResult->GetOutputImage().c_str();
- DimReduction dimReduction = DimReduction_NONE;
+ DimReductionMethod dimReduction = DimReductionMethod_NONE;
if ( parseResult->IsOptionPresent("--DimReduceMethod") )
{
std::string dimReductionStr = parseResult->GetParameterString("--DimReduceMethod");
- if ( boost::to_upper(dimReductionStr) == "PCA" )
+ if ( boost::to_upper_copy(dimReductionStr) == "PCA" )
{
- dimReduction = DimReduction_PCA;
+ dimReduction = DimReductionMethod_PCA;
}
- else if ( boost::to_upper(dimReductionStr) == "MNF" )
+ else if ( boost::to_upper_copy(dimReductionStr) == "MNF" )
{
- dimReduction = DimReduction_MNF;
+ dimReduction = DimReductionMethod_MNF;
}
}
- otbMsgDevMacro(<< "Using "
- << (dimReduction == DimReduction_NONE ? "NONE" : (dimReduction == DimReduction_PCA ? : "PCA" : DimReduction_MNF) )
- << " dimensionality reduction method"
- );
+ otbMsgDevMacro( << "Using "
+ << (dimReduction == DimReductionMethod_NONE ? "NONE" : (dimReduction == DimReductionMethod_PCA ? : "PCA" : "MNF") )
+ << " dimensionality reduction method" );
unsigned int nbEndmembers = 0;
if ( parseResult->IsOptionPresent("--NumEndmembers") )
@@ -126,7 +150,7 @@ int main(int argc, char * argv[])
if ( parseResult->IsOptionPresent("--DimensionalityEstimationMethod") )
{
std::string dimEstMethodStr = parseResult->GetParameterString("--DimensionalityEstimationMethod");
- if (boost::to_upper(dimEstMethodStr) != "ELM")
+ if (boost::to_upper_copy(dimEstMethodStr) != "ELM")
{
std::cerr << "Only ELM is supported for parameter DimensionalityEstimationMethod" << std::endl;
return EXIT_FAILURE;
@@ -136,7 +160,7 @@ int main(int argc, char * argv[])
if ( parseResult->IsOptionPresent("--EndmembersEstimationMethod") )
{
std::string eeEstMethodStr = parseResult->GetParameterString("--EndmembersEstimationMethod");
- if (boost::to_upper(eeEstMethodStr) != "VCA")
+ if (boost::to_upper_copy(eeEstMethodStr) != "VCA")
{
std::cerr << "Only VCA is supported for parameter EndmembersEstimationMethod" << std::endl;
return EXIT_FAILURE;
@@ -146,23 +170,22 @@ int main(int argc, char * argv[])
const char * inputEndmembers = parseResult->IsOptionPresent("--InputEndmembers") ?
parseResult->GetParameterString("--InputEndmembers").c_str() : NULL;
- UnMixingAlgo unmixingAlgo = UnMixingAlgo_UCLS;
- if ( parseResult->IsOptionPresent("--DimReduceMethod") )
+ UnMixingMethod unmixingAlgo = UnMixingMethod_NONE;
+ if ( parseResult->IsOptionPresent("--UnmixingAlgorithm") )
{
- std::string unmixingAlgoStr = parseResult->GetParameterString("--DimReduceMethod");
- if ( boost::to_upper(unmixingAlgoStr) == "UCLS" )
+ std::string unmixingAlgoStr = parseResult->GetParameterString("--UnmixingAlgorithm");
+ if ( boost::to_upper_copy(unmixingAlgoStr) == "UCLS" )
{
- unmixingAlgo = UnMixingAlgo_UCLS;
+ unmixingAlgo = UnMixingMethod_UCLS;
}
- else if ( boost::to_upper(unmixingAlgoStr) == "ISRA" )
+ else if ( boost::to_upper_copy(unmixingAlgoStr) == "ISRA" )
{
- unmixingAlgo = UnMixingAlgo_ISRA;
+ unmixingAlgo = UnMixingMethod_ISRA;
}
}
- otbMsgDevMacro(<< "Using "
- << (unmixingAlgo == UnMixingAlgo_UCLS ? "UCLS" : "ISRA" )
- << " unmixing algorithm"
- );
+ otbMsgDevMacro( << "Using "
+ << (unmixingAlgo == UnMixingMethod_NONE ? "NONE" : (unmixingAlgo == UnMixingMethod_UCLS ? : "UCLS" : "ISRA") )
+ << " unmixing algorithm" );
const char * outputEndmembers = parseResult->IsOptionPresent("--OutputEndmembers") ?
parseResult->GetParameterString("--OutputEndmembers").c_str() : NULL;
@@ -178,10 +201,27 @@ int main(int argc, char * argv[])
*
*/
+ ReaderType::Pointer readerImage = ReaderType::New();
+ readerImage->SetFileName(inputImageName);
+
+ VectorImageType::Pointer inputImage = readerImage->GetOutput();
/*
* Compute stats of input image, we will need this all along the road
*/
+ std::cout << "Computing stats" << std::endl;
+
+ StreamingStatisticsVectorImageFilterType::Pointer stats = StreamingStatisticsVectorImageFilterType::New();
+
+ stats->SetInput(inputImage);
+ stats->Update();
+
+ VectorType mean (stats->GetMean().GetDataPointer(), stats->GetMean().GetSize());
+ MatrixType covariance = stats->GetCovariance().GetVnlMatrix();
+ MatrixType correlation = stats->GetCorrelation().GetVnlMatrix();
+
+ VectorImageType::Pointer endmembersImage;
+ itk::ProcessObject::Pointer endmembersRef;
if ( inputEndmembers == NULL )
{
@@ -190,49 +230,119 @@ int main(int argc, char * argv[])
/*
* 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;
}
/*
* Estimate Endmembers
*/
+ std::cout << "Estimate Endmembers by VCA" << std::endl;
+ VCAFilterType::Pointer vca = VCAFilterType::New();
+ vca->SetNumberOfEndmembers(nbEndmembers);
+ vca->SetInput(inputImage);
+
+ endmembersImage = vca->GetOutput();
+ endmembersRef = vca;
}
else
{
/*
* Read input endmembers
*/
+ std::cout << "Read Endmembers " << inputEndmembers << std::endl;
+ ReaderType::Pointer readerEndmembers = ReaderType::New();
+ readerEndmembers->SetFileName(inputEndmembers);
+ endmembersImage = readerEndmembers->GetOutput();
+ endmembersRef = readerEndmembers;
}
-
/*
- * Unmix
+ * Transform Endmembers image to matrix representation
*/
+ VectorImageToMatrixImageFilterType::Pointer endMember2Matrix = VectorImageToMatrixImageFilterType::New();
+ endMember2Matrix->SetInput(endmembersImage);
+ endMember2Matrix->Update();
+ MatrixType endMembersMatrix = endMember2Matrix->GetMatrix();
/*
- * Write abundance map
+ * Unmix
*/
- WriterType::Pointer writer = WriterType::New();
- writer->SetFileName(outputImageName);
- //writer->SetInput(scale->GetOutput());
- writer->Update();
+ VectorImageType::Pointer abundanceMap;
+ itk::ProcessObject::Pointer unmixerRef;
+
+ switch (unmixingAlgo)
+ {
+ case UnMixingMethod_NONE:
+ break;
+ case UnMixingMethod_UCLS:
+ {
+ std::cout << "UCLS Unmixing" << std::endl;
+
+ UCLSUnmixingFilterType::Pointer unmixer =
+ UCLSUnmixingFilterType::New();
+ unmixer->SetInput(readerImage->GetOutput());
+ unmixer->SetMatrix(endMembersMatrix);
+ unmixer->SetNumberOfThreads(1); // FIXME : currently buggy
+
+ unmixerRef = unmixer;
+ abundanceMap = unmixer->GetOutput();
+ }
+ break;
+ case UnMixingMethod_ISRA:
+ {
+ std::cout << "ISRA Unmixing" << std::endl;
+
+ ISRAUnmixingFilterType::Pointer unmixer =
+ ISRAUnmixingFilterType::New();
+
+ unmixer->SetInput(readerImage->GetOutput());
+ unmixer->SetEndmembersMatrix(endMembersMatrix);
+ unmixerRef = unmixer;
+ abundanceMap = unmixer->GetOutput();
+ }
+ break;
+ default:
+ break;
+ }
if ( outputEndmembers != NULL )
{
/*
* Write endmembers
*/
+ std::cout << "Write endmembers " << outputEndmembers << std::endl;
WriterType::Pointer endmembersWriter = WriterType::New();
endmembersWriter->SetFileName(outputEndmembers);
- //endmembersWriter->SetInput(scale->GetOutput());
+ endmembersWriter->SetInput(endmembersImage);
endmembersWriter->Update();
}
+ if ( unmixingAlgo != UnMixingMethod_NONE )
+ {
+ /*
+ * Write abundance map
+ */
+ std::cout << "Write abundance map" << outputImageName << std::endl;
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetFileName(outputImageName);
+ writer->SetInput(abundanceMap);
+ writer->Update();
+ }
return EXIT_SUCCESS;
}