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

ENH: fill up HyperspectralUnmixing application

parent 90d129b8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Application/CMakeLists.txt
+ 29
− 0
View file @ 1f9ed549
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
Application/otbHyperspectralUnmixing.cxx
+ 152
− 42
View file @ 1f9ed549
......@@ -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;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment