diff --git a/Examples/Hyperspectral/HyperspectralUnmixing.cxx b/Examples/Hyperspectral/HyperspectralUnmixing.cxx index e44701971f64389c3336690fd0f62a06c217aa5f..26476ec34b8a20b7165874a9da79d9cb9de6864f 100644 --- a/Examples/Hyperspectral/HyperspectralUnmixing.cxx +++ b/Examples/Hyperspectral/HyperspectralUnmixing.cxx @@ -19,118 +19,186 @@ #include "otbImageFileReader.h" #include "otbImageFileWriter.h" #include "otbStreamingStatisticsVectorImageFilter.h" -#include "otbVcaImageFilter.h" -#include "otbUnConstrainedLeastSquareImageFilter.h" #include "otbVectorRescaleIntensityImageFilter.h" #include "otbVectorImageToMatrixImageFilter.h" #include "itkRescaleIntensityImageFilter.h" #include "otbVectorRescaleIntensityImageFilter.h" #include "otbMultiToMonoChannelExtractROI.h" // Software Guide : BeginCommandLineArgs -// INPUTS: {ROI_QB_MUL_1.png} -// OUTPUTS: {PCAOutput.tif}, {InversePCAOutput.tif}, {InversePCAOutput1.png}, {PCAOutput1.png}, {PCAOutput2.png}, {PCAOutput3.png} -// 3 +// INPUTS: {Indian_pines_corrected.tif} +// OUTPUTS: {Indian_pines_Unmixing_Output.tif}, {UnmixingbandOutput1.png}, {UnmixingbandOutput2.png}, {UnmixingbandOutput3.png} +// 16 // Software Guide : EndCommandLineArgs + +// Software Guide : BeginLatex +// +// 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. +// +// Software Guide : EndLatex + +// Software Guide : BeginCodeSnippet +#include "otbVcaImageFilter.h" +#include "otbUnConstrainedLeastSquareImageFilter.h" +// Software Guide : EndCodeSnippet int main(int argc, char * argv[]) { - if (argc != 5) + if (argc != 7) { std::cerr << "Usage: " << argv[0]; - std::cerr << "inputFileName outputFileName EstimatenumberOfEndmembers outputPrettyFilename1 outputPrettyFilename2 outputPrettyFilename3 " << std::endl; + 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[3]); + const unsigned int estimateNumberOfEndmembers = atoi(argv[6]); const unsigned int Dimension = 2; - typedef double PixelType; - typedef otb::VectorImage<PixelType, Dimension> ImageType; - typedef otb::VectorImage<unsigned char, 2> OutputPrettyImageType; + typedef double PixelType; + typedef otb::VectorImage<unsigned char, 2> OutputPrettyImageType; - typedef otb::ImageFileReader<ImageType> ReaderType; - typedef otb::ImageFileWriter<ImageType> WriterType; - - typedef otb::VectorRescaleIntensityImageFilter<ImageType,ImageType> RescalerType; - typedef otb::VectorRescaleIntensityImageFilter<ImageType,OutputPrettyImageType> RescalerType2; - - typedef otb::StreamingStatisticsVectorImageFilter<ImageType> StreamingStatisticsVectorImageFilterType; - typedef otb::VCAImageFilter<ImageType> VCAFilterType; + // Software Guide : BeginLatex + // + // 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. + // + // Software Guide : EndLatex - typedef otb::UnConstrainedLeastSquareImageFilter<ImageType,ImageType,double> UCLSUnmixingFilterType; - - typedef otb::VectorImageToMatrixImageFilter<ImageType> VectorImageToMatrixImageFilterType; + // Software Guide : BeginCodeSnippet + typedef otb::VectorImage<PixelType, Dimension> ImageType; + typedef otb::ImageFileReader<ImageType> ReaderType; + typedef otb::ImageFileWriter<ImageType> WriterType; + // Software Guide : EndCodeSnippet + + typedef otb::VectorRescaleIntensityImageFilter<ImageType,ImageType> RescalerType; + typedef otb::StreamingStatisticsVectorImageFilter<ImageType> StreamingStatisticsVectorImageFilterType; + typedef otb::VectorImageToMatrixImageFilter<ImageType> VectorImageToMatrixImageFilterType; typedef vnl_vector<double> VectorType; typedef vnl_matrix<double> MatrixType; - /// / Noise filtering - // typedef otb::LocalActivityVectorImageFilter< ImageType, ImageType > NoiseFilterType; + // Software Guide : EndCodeSnippet + // Software Guide : BeginLatex + // + // We instantiate now the image reader and we set the image file name. + // + // Software Guide : EndLatex - // // Image filtering - // typedef otb::MNFImageFilter< ImageType, ImageType, - // NoiseFilterType, otb::Transform::FORWARD > FilterType; - - //read image + // Software Guide : BeginCodeSnippet ReaderType::Pointer reader = ReaderType::New(); reader->SetFileName(infname); + // Software Guide : EndCodeSnippet + reader->UpdateOutputInformation(); + // Software Guide : BeginLatex + // + // For now we need to rescale the input image between 0 and 1 to + // perform the unmixing algorithm. We use the + // \doxygen{otb}{VectorRescaleIntensityImageFilter}. + // + // Software Guide : EndLatex //rescale input image between 0 and 1 + + // Software Guide : BeginCodeSnippet RescalerType::Pointer rescaler = RescalerType::New(); rescaler->SetInput(reader->GetOutput()); - rescaler->SetOutputMinimum(0); - rescaler->SetOutputMaximum(1.); - /* - * Estimate Endmembers - */ - std::cout << "Estimate Endmembers by VCA" << std::endl; - + + 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); + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // 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. + // + // Software Guide : EndLatex + + // Software Guide : BeginCodeSnippet + typedef otb::VCAImageFilter<ImageType> VCAFilterType; VCAFilterType::Pointer vca = VCAFilterType::New(); + vca->SetNumberOfEndmembers(estimateNumberOfEndmembers); vca->SetInput(rescaler->GetOutput()); + // Software Guide : EndCodeSnippet - ImageType::Pointer endmembersImage; - endmembersImage = vca->GetOutput(); - - /* - * Transform Endmembers image to matrix representation - */ - std::cout << "Endmembers extracted" << std::endl; - std::cout << "Converting endmembers to matrix" << std::endl; + + // Software Guide : BeginLatex + // + // We transform the output estimate endmembers to a matrix representation + // + // Software Guide : EndLatex + + // Software Guide : BeginCodeSnippet VectorImageToMatrixImageFilterType::Pointer endMember2Matrix = VectorImageToMatrixImageFilterType::New(); - endMember2Matrix->SetInput(endmembersImage); + endMember2Matrix->SetInput(vca->GetOutput()); endMember2Matrix->Update(); - - MatrixType endMembersMatrix = endMember2Matrix->GetMatrix(); - std::cout << "Endmembers matrix : " << endMembersMatrix << std::endl; - - /* - * Unmix - */ - UCLSUnmixingFilterType::Pointer unmixer = - UCLSUnmixingFilterType::New(); - + + // Software Guide : EndCodeSnippet + + // Software Guide : BeginLatex + // + // We can now procedd to the unmixing algorithm.Parameters + // needed are the input image and the endmembers matrix. + // + // Software Guide : EndLatex + + // Software Guide : BeginCodeSnippet + typedef otb::UnConstrainedLeastSquareImageFilter<ImageType,ImageType,double> UCLSUnmixingFilterType; + UCLSUnmixingFilterType::Pointer unmixer = UCLSUnmixingFilterType::New(); unmixer->SetInput(rescaler->GetOutput()); - unmixer->SetMatrix(endMembersMatrix); - unmixer->SetNumberOfThreads(1); // FIXME : currently buggy - - ImageType::Pointer abundanceMap; - abundanceMap = unmixer->GetOutput(); - - /* - * Write abundance map - */ - std::cout << "Write abundance map " << outfname << std::endl; + unmixer->SetMatrix(endMember2Matrix->GetMatrix()); + // Software Guide : EndCodeSnippet + unmixer->SetNumberOfThreads(1); // FIXME : currently buggy + + // Software Guide : BeginLatex + // + // 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. + // + // Software Guide : EndLatex + + // Software Guide : BeginCodeSnippet WriterType::Pointer writer = WriterType::New(); - writer->SetInput(abundanceMap); + writer->SetInput(unmixer->GetOutput()); writer->SetFileName(outfname); - writer->Update(); + // Software Guide : EndCodeSnippet + // Software Guide : BeginLatex + // Figure~\ref{fig:UNMIXING_FILTER} shows the result of applying unmixing + // to a 220 band AVIRIS image. + // \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} + // + // Software Guide : EndLatex typedef otb::Image<PixelType, Dimension> MonoImageType; typedef otb::MultiToMonoChannelExtractROI<PixelType, PixelType> ExtractROIFilterType; typedef otb::Image<unsigned char, 2> OutputImageType; @@ -144,7 +212,7 @@ int main(int argc, char * argv[]) PrettyRescalerType::Pointer prettyRescaler = PrettyRescalerType::New(); WriterType2::Pointer writer2 = WriterType2::New(); - extractROIFilter->SetInput(abundanceMap); + extractROIFilter->SetInput(unmixer->GetOutput()); extractROIFilter->SetChannel(cpt + 1); prettyRescaler->SetInput(extractROIFilter->GetOutput()); @@ -152,7 +220,7 @@ int main(int argc, char * argv[]) prettyRescaler->SetOutputMaximum(255); writer2->SetInput(prettyRescaler->GetOutput()); - writer2->SetFileName(argv[cpt + 4]); + writer2->SetFileName(argv[cpt + 3]); writer2->Update(); }