HyperspectralUnmixingExample.cxxΒΆ
Example usage:
./HyperspectralUnmixingExample Input/AVIRIS/Indian_pines_corrected.tif \
Output/Indian_pines_Unmixing_Output.tif \
Output/UnmixingbandOutput1.png \
Output/UnmixingbandOutput2.png \
Output/UnmixingbandOutput3.png \
16
Example source code (HyperspectralUnmixingExample.cxx):
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "otbVectorRescaleIntensityImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "otbVectorRescaleIntensityImageFilter.h"
#include "otbMultiToMonoChannelExtractROI.h"
// This example illustrates the use of the
// \doxygen{otb}{VcaImageFilter} and \doxygen{otb}{UnConstrainedLeastSquareImageFilter}.
// The VCA filter computes endmembers using the Vertex Component Analysis
// and UCLS performs unmixing on the input image using these endmembers.
//
// The first step required to use these filters is to include its header files.
#include "otbVcaImageFilter.h"
#include "otbUnConstrainedLeastSquareImageFilter.h"
#include "itkMersenneTwisterRandomVariateGenerator.h"
int main(int argc, char* argv[])
{
if (argc != 7)
{
std::cerr << "Usage: " << argv[0];
std::cerr << "inputFileName outputFileName outputPrettyFilename1 outputPrettyFilename2 outputPrettyFilename3 EstimatenumberOfEndmembers" << std::endl;
return EXIT_FAILURE;
}
const char* infname = argv[1];
const char* outfname = argv[2];
const unsigned int estimateNumberOfEndmembers = atoi(argv[6]);
const unsigned int Dimension = 2;
using PixelType = double;
// We start by defining the types for the images and the reader and
// the writer. We choose to work with a \doxygen{otb}{VectorImage},
// since we will produce a multi-channel image (the principal
// components) from a multi-channel input image.
using ImageType = otb::VectorImage<PixelType, Dimension>;
using ReaderType = otb::ImageFileReader<ImageType>;
using WriterType = otb::ImageFileWriter<ImageType>;
using RescalerType = otb::VectorRescaleIntensityImageFilter<ImageType, ImageType>;
using VectorImageToMatrixImageFilterType = otb::VectorImageToMatrixImageFilter<ImageType>;
// We set the seed of the random number generator.
itk::Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed(121212);
// We instantiate now the image reader and we set the image file name.
ReaderType::Pointer reader = ReaderType::New();
reader->SetFileName(infname);
reader->UpdateOutputInformation();
// For now we need to rescale the input image between 0 and 1 to
// perform the unmixing algorithm. We use the
// \doxygen{otb}{VectorRescaleIntensityImageFilter}.
RescalerType::Pointer rescaler = RescalerType::New();
rescaler->SetInput(reader->GetOutput());
ImageType::PixelType minPixel, maxPixel;
minPixel.SetSize(reader->GetOutput()->GetNumberOfComponentsPerPixel());
maxPixel.SetSize(reader->GetOutput()->GetNumberOfComponentsPerPixel());
minPixel.Fill(0.);
maxPixel.Fill(1.);
rescaler->SetOutputMinimum(minPixel);
rescaler->SetOutputMaximum(maxPixel);
// We define the type for the VCA filter. It is templated over the input
// image type. The only parameter is the number of endmembers which
// needs to be estimated.
// We can now the instantiate the filter.
using VCAFilterType = otb::VCAImageFilter<ImageType>;
VCAFilterType::Pointer vca = VCAFilterType::New();
vca->SetNumberOfEndmembers(estimateNumberOfEndmembers);
vca->SetInput(rescaler->GetOutput());
// We transform the output estimate endmembers to a matrix representation
VectorImageToMatrixImageFilterType::Pointer endMember2Matrix = VectorImageToMatrixImageFilterType::New();
endMember2Matrix->SetInput(vca->GetOutput());
endMember2Matrix->Update();
// We can now procedd to the unmixing algorithm.Parameters
// needed are the input image and the endmembers matrix.
using UCLSUnmixingFilterType = otb::UnConstrainedLeastSquareImageFilter<ImageType, ImageType, double>;
UCLSUnmixingFilterType::Pointer unmixer = UCLSUnmixingFilterType::New();
unmixer->SetInput(rescaler->GetOutput());
unmixer->GetModifiableFunctor().SetMatrix(endMember2Matrix->GetMatrix());
unmixer->SetNumberOfThreads(1); // FIXME : currently buggy
// We now instantiate the writer and set the file name for the
// output image and trigger the unmixing computation with the \code{Update()} of the writer.
WriterType::Pointer writer = WriterType::New();
writer->SetInput(unmixer->GetOutput());
writer->SetFileName(outfname);
writer->Update();
// Figure~\ref{fig:UNMIXING_FILTER} shows the result of applying unmixing
// to an AVIRIS image (220 bands). This dataset is freely available
// at \url{http://www.ehu.es/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes#Indian_Pines}
// \begin{figure}
// \center
// \includegraphics[width=0.25\textwidth]{UnmixingbandOutput1.eps}
// \includegraphics[width=0.25\textwidth]{UnmixingbandOutput2.eps}
// \includegraphics[width=0.25\textwidth]{UnmixingbandOutput3.eps}
// \itkcaption[Unmixing Filter]{Result of applying the
// \doxygen{otb}{VcaImageFilter} and \doxygen{otb}{UnConstrainedLeastSquareImageFilter} to an image. From left
// to right and top to bottom:
// first abundance map band, second abundance map band, third abundance map band.}
// \label{fig:UNMIXING_FILTER}
// \end{figure}
using MonoImageType = otb::Image<PixelType, Dimension>;
using ExtractROIFilterType = otb::MultiToMonoChannelExtractROI<PixelType, PixelType>;
using OutputImageType = otb::Image<unsigned char, 2>;
using PrettyRescalerType = itk::RescaleIntensityImageFilter<MonoImageType, OutputImageType>;
using WriterType2 = otb::ImageFileWriter<OutputImageType>;
for (unsigned int cpt = 0; cpt < 3; ++cpt)
{
ExtractROIFilterType::Pointer extractROIFilter = ExtractROIFilterType::New();
PrettyRescalerType::Pointer prettyRescaler = PrettyRescalerType::New();
WriterType2::Pointer writer2 = WriterType2::New();
extractROIFilter->SetInput(unmixer->GetOutput());
extractROIFilter->SetChannel(cpt + 1);
prettyRescaler->SetInput(extractROIFilter->GetOutput());
prettyRescaler->SetOutputMinimum(0);
prettyRescaler->SetOutputMaximum(255);
writer2->SetInput(prettyRescaler->GetOutput());
writer2->SetFileName(argv[cpt + 3]);
writer2->Update();
}
return EXIT_SUCCESS;
}