Commit fe9aa991 authored by Julien Malik's avatar Julien Malik

BUG: remove CLSPSTO and FCLS unmixing algorithms, not ready yet for release

parent 1798508e
......@@ -22,7 +22,7 @@
#include "otbUnConstrainedLeastSquareImageFilter.h"
#include "otbISRAUnmixingImageFilter.h"
#include "otbNCLSUnmixingImageFilter.h"
#include "otbFCLSUnmixingImageFilter.h"
//#include "otbFCLSUnmixingImageFilter.h"
#include "otbMDMDNMFImageFilter.h"
#include "otbVectorImageToMatrixImageFilter.h"
......@@ -40,7 +40,7 @@ typedef otb::StreamingStatisticsVectorImageFilter<DoubleVectorImageType> Streami
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::FCLSUnmixingImageFilter<DoubleVectorImageType, DoubleVectorImageType, double> FCLSUnmixingFilterType;
typedef otb::MDMDNMFImageFilter<DoubleVectorImageType, DoubleVectorImageType> MDMDNMFUnmixingFilterType;
typedef otb::VectorImageToMatrixImageFilter<DoubleVectorImageType> VectorImageToMatrixImageFilterType;
......@@ -69,7 +69,7 @@ enum EndmembersEstimationMethod
enum UnMixingMethod
{
UnMixingMethod_UCLS,
UnMixingMethod_FCLS,
//UnMixingMethod_FCLS,
UnMixingMethod_NCLS,
UnMixingMethod_ISRA,
UnMixingMethod_MDMDNMF,
......@@ -122,8 +122,8 @@ private:
AddChoice("ua.ucls", "UCLS");
SetParameterDescription("ua.ucls", "Unconstrained Least Square");
AddChoice("ua.fcls", "FCLS");
SetParameterDescription("ua.fcls", "Fully constrained Least Square");
// AddChoice("ua.fcls", "FCLS");
// SetParameterDescription("ua.fcls", "Fully constrained Least Square");
AddChoice("ua.ncls", "NCLS");
SetParameterDescription("ua.ncls", "Non-negative constrained Least Square");
......@@ -216,6 +216,7 @@ private:
}
break;
/*
case UnMixingMethod_FCLS:
{
std::cout << "FCLS Unmixing" << std::endl;
......@@ -230,6 +231,7 @@ private:
}
break;
*/
case UnMixingMethod_MDMDNMF:
{
std::cout << "MDMD-NMF Unmixing" << std::endl;
......
/*=========================================================================
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.
=========================================================================*/
#ifndef __otbFCLSUnmixingImageFilter_h
#define __otbFCLSUnmixingImageFilter_h
#include "itkMacro.h"
#include "itkNumericTraits.h"
#include "otbUnaryFunctorImageFilter.h"
#include "vnl/algo/vnl_svd.h"
#include <boost/shared_ptr.hpp>
namespace otb
{
namespace Functor {
/** \class FCLSUnmixingFunctor
*
* \brief TODO
*
*/
template<class TInput, class TOutput, class TPrecision>
class FCLSUnmixingFunctor
{
public:
typedef FCLSUnmixingFunctor<TInput, TOutput, TPrecision> Self;
typedef TInput InputType;
typedef TOutput OutputType;
typedef TPrecision PrecisionType;
typedef vnl_vector<PrecisionType> VectorType;
typedef vnl_matrix<PrecisionType> MatrixType;
FCLSUnmixingFunctor();
virtual ~FCLSUnmixingFunctor();
unsigned int GetOutputSize() const;
bool operator !=(const FCLSUnmixingFunctor& other) const;
bool operator ==(const FCLSUnmixingFunctor& other) const;
void SetEndmembersMatrix(const MatrixType& U);
const MatrixType& GetEndmembersMatrix(void) const;
void SetMaxIteration(unsigned int val)
{
m_MaxIteration = val;
}
unsigned int GetMaxIteration() const
{
return m_MaxIteration;
}
OutputType operator ()(const InputType& in) const;
private:
static bool IsNonNegative(PrecisionType val)
{
return val >= 0;
}
typedef vnl_svd<PrecisionType> SVDType;
typedef boost::shared_ptr<SVDType> SVDPointerType;
MatrixType m_U;
MatrixType m_N;
MatrixType m_Nt;
MatrixType m_NtNinv;
SVDPointerType m_Svd; // SVD of U
unsigned int m_OutputSize;
unsigned int m_MaxIteration;
};
}
/** \class FCLSUnmixingImageFilter
*
* \brief Performs fully constrained least squares on each pixel of a VectorImage
*
* This filter takes as input a multiband image and a matrix.
* If the matrix is called $A$, it solves, for each pixel $p$, the system
* \f$A \cdot x = p\f$ in the least square sense, with additionnal constraints on the solution
* \f$\hat{x}\f$ ensuring positivity (each component is positive) and additivity (the sum of
* all components is 1).
*
* The main use of this filter is to unmix an hyperspectral dataset,
* where \f$A\f$ is the mixing matrix, in which each row corresponds to an endmember signature.
*
* The number of rows in \f$A\f$ must match the input image number of bands.
* The number of bands in the output image will be the number of columns of $A$
*
* References
* "Fully Constrained Least-Squares Based Linear Unmixing." Daniel Heinz,
* Chein-I Chang, and Mark L.G. Althouse. IEEE. 1999.
*
* \ingroup Hyperspectral
* \ingroup Streamed
* \ingroup Threaded
*/
template <class TInputImage, class TOutputImage, class TPrecision>
class ITK_EXPORT FCLSUnmixingImageFilter :
public otb::UnaryFunctorImageFilter<TInputImage, TOutputImage,
Functor::FCLSUnmixingFunctor<typename TInputImage::PixelType,
typename TOutputImage::PixelType, TPrecision> >
{
public:
/** Standard class typedefs. */
typedef FCLSUnmixingImageFilter Self;
typedef otb::UnaryFunctorImageFilter
<TInputImage,
TOutputImage,
Functor::FCLSUnmixingFunctor<
typename TInputImage::PixelType,
typename TOutputImage::PixelType,
TPrecision>
> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef Functor::FCLSUnmixingFunctor<
typename TInputImage::PixelType,
typename TOutputImage::PixelType,
TPrecision> FunctorType;
typedef typename FunctorType::MatrixType MatrixType;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Run-time type information (and related methods). */
itkTypeMacro(FCLSUnmixingImageFilter, otb::UnaryFunctorImageFilter);
/** Pixel types. */
typedef typename TInputImage::PixelType InputPixelType;
typedef typename TOutputImage::PixelType OutputPixelType;
void SetEndmembersMatrix(const MatrixType& m);
const MatrixType& GetEndmembersMatrix() const;
void SetMaxIteration( unsigned int val )
{
this->GetFunctor().SetMaxIteration(val);
this->Modified();
}
unsigned int GetMaxIteration() const
{
return this->GetFunctor().GetMaxIteration();
}
protected:
FCLSUnmixingImageFilter();
virtual ~FCLSUnmixingImageFilter();
void PrintSelf(std::ostream& os, itk::Indent indent) const;
private:
FCLSUnmixingImageFilter(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbFCLSUnmixingImageFilter.txx"
#endif
#endif
/*=========================================================================
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.
=========================================================================*/
#ifndef __otbFCLSUnmixingImageFilter_txx
#define __otbFCLSUnmixingImageFilter_txx
#include "otbFCLSUnmixingImageFilter.h"
#include <algorithm>
namespace otb
{
namespace Functor
{
template <class TInput, class TOutput, class TPrecision>
FCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::FCLSUnmixingFunctor()
: m_OutputSize(0),
m_MaxIteration(10)
{
}
template <class TInput, class TOutput, class TPrecision>
FCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::~FCLSUnmixingFunctor()
{
}
template <class TInput, class TOutput, class TPrecision>
unsigned int
FCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::GetOutputSize() const
{
return m_OutputSize;
}
template <class TInput, class TOutput, class TPrecision>
bool
FCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::operator != (const Self& other) const
{
return true;
}
template <class TInput, class TOutput, class TPrecision>
bool
FCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::operator == (const Self& other) const
{
return false;
}
template <class TInput, class TOutput, class TPrecision>
void
FCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::SetEndmembersMatrix(const MatrixType& U)
{
//const double Delta = 1.0E-6;
m_U = U;
m_N.set_size(U.rows() + 1, U.cols());
for (unsigned int r = 0; r < U.rows(); ++r)
{
m_N.set_row(r, U.get_row(r));
}
m_N.set_row(U.rows(), 1.0E6);
m_Nt = m_N.transpose();
MatrixType NtN = m_Nt * m_N;
m_NtNinv = SVDType(NtN).inverse();
m_OutputSize = m_N.cols();
m_Svd.reset( new SVDType(m_U) );
}
template <class TInput, class TOutput, class TPrecision>
const typename FCLSUnmixingFunctor<TInput, TOutput, TPrecision>::MatrixType&
FCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::GetEndmembersMatrix() const
{
return m_U;
}
template <class TInput, class TOutput, class TPrecision>
typename FCLSUnmixingFunctor<TInput, TOutput, TPrecision>::OutputType
FCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::operator ()(const InputType& in) const
{
// TODO : support different types between input and output ?
VectorType inVector(in.Size());
for (unsigned int i = 0; i < in.GetSize(); ++i )
{
inVector[i] = in[i];
}
// Initialize with Unconstrained Least Square solution
VectorType uclsVector = m_Svd->solve(inVector);
unsigned int nbEndmembers = m_OutputSize;
unsigned int nbBands = in.Size();
// Apply FCLS iterations
//const double Delta = 1.0E-6;
VectorType s(nbBands + 1);
for (unsigned int r = 0; r < nbBands; ++r)
{
s[r] = inVector[r];
}
s[nbBands] = 1.0E6;
VectorType lambda(nbEndmembers);
VectorType fclsVector = uclsVector;
VectorType correction(uclsVector.size());
for (unsigned int i = 0; i < m_MaxIteration; ++i)
{
// Error in original paper : divergence
// lambda = m_Nt * (s - m_N * fclsVector);
lambda = m_Nt * (m_N * fclsVector - s);
correction = m_NtNinv * lambda;
fclsVector -= correction;
}
OutputType out(fclsVector.size());
for (unsigned int i = 0; i < out.GetSize(); ++i )
{
out[i] = fclsVector[i];
}
return out;
}
}
template <class TInputImage, class TOutputImage, class TPrecision>
FCLSUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::FCLSUnmixingImageFilter()
{
}
template <class TInputImage, class TOutputImage, class TPrecision>
FCLSUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::~FCLSUnmixingImageFilter()
{
}
template <class TInputImage, class TOutputImage, class TPrecision>
void
FCLSUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::SetEndmembersMatrix(const MatrixType& m)
{
this->GetFunctor().SetEndmembersMatrix(m);
this->Modified();
}
template <class TInputImage, class TOutputImage, class TPrecision>
const typename FCLSUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>::MatrixType&
FCLSUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::GetEndmembersMatrix() const
{
return this->GetFunctor().GetEndmembersMatrix();
}
template <class TInputImage, class TOutputImage, class TPrecision>
void
FCLSUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
} // end namespace
#endif
......@@ -121,22 +121,6 @@ ADD_TEST(hyTvISRAUnmixingImageFilterTest
${TEMP}/hyTvISRAUnmixingImageFilterTest.tif
10)
ADD_TEST(hyTuFCLSUnmixingImageFilterNew
${Hyperspectral_TESTS1}
otbFCLSUnmixingImageFilterNewTest)
ADD_TEST(hyTvFCLSUnmixingImageFilterTest
${Hyperspectral_TESTS1}
--compare-image ${EPSILON_9}
${BASELINE_FILES}/hyTvFCLSUnmixingImageFilterTest.tif
${TEMP}/hyTvFCLSUnmixingImageFilterTest.tif
otbFCLSUnmixingImageFilterTest
${INPUTDATA}/Hyperspectral/synthetic/hsi_cube.tif
${INPUTDATA}/Hyperspectral/synthetic/endmembers.tif
${TEMP}/hyTvFCLSUnmixingImageFilterTest.tif
10)
# HyperspectralTests3
ADD_TEST(hyTuMDMDNMFImageFilterNew
${Hyperspectral_TESTS3}
......@@ -165,7 +149,6 @@ SET(Hyperspectral_SRCS1
otbHyperspectralTests1.cxx
otbEigenvalueLikelihoodMaximization.cxx
otbUnConstrainedLeastSquareImageFilter.cxx
otbFCLSUnmixingImageFilter.cxx
otbISRAUnmixingImageFilter.cxx
otbNCLSUnmixingImageFilter.cxx
otbSparseUnmixingImageFilterNew.cxx
......
/*=========================================================================
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 "otbCLSPSTOUnmixingImageFilter.h"
#include "otbVectorImage.h"
#include "otbImageFileReader.h"
#include "otbStreamingImageFileWriter.h"
#include "otbVectorImageToMatrixImageFilter.h"
#include "otbStandardWriterWatcher.h"
const unsigned int Dimension = 2;
typedef float PixelType;
typedef otb::VectorImage<PixelType, Dimension> ImageType;
typedef otb::ImageFileReader<ImageType> ReaderType;
typedef otb::CLSPSTOUnmixingImageFilter<ImageType, ImageType, float> UnmixingImageFilterType;
typedef otb::VectorImageToMatrixImageFilter<ImageType> VectorImageToMatrixImageFilterType;
typedef otb::StreamingImageFileWriter<ImageType> WriterType;
int otbCLSPSTOUnmixingImageFilterNewTest(int argc, char * argv[])
{
UnmixingImageFilterType::Pointer filter = UnmixingImageFilterType::New();
std::cout << filter << std::endl;
return EXIT_SUCCESS;
}
int otbCLSPSTOUnmixingImageFilterTest(int argc, char * argv[])
{
const char * inputImage = argv[1];
const char * inputEndmembers = argv[2];
const char * outputImage = argv[3];
int maxIter = atoi(argv[4]);
ReaderType::Pointer readerImage = ReaderType::New();
readerImage->SetFileName(inputImage);
ReaderType::Pointer readerEndMembers = ReaderType::New();
readerEndMembers->SetFileName(inputEndmembers);
VectorImageToMatrixImageFilterType::Pointer endMember2Matrix = VectorImageToMatrixImageFilterType::New();
endMember2Matrix->SetInput(readerEndMembers->GetOutput());
endMember2Matrix->Update();
typedef VectorImageToMatrixImageFilterType::MatrixType MatrixType;
MatrixType endMembers = endMember2Matrix->GetMatrix();
MatrixType pinv = vnl_matrix_inverse<PixelType>(endMembers);
UnmixingImageFilterType::Pointer unmixer = UnmixingImageFilterType::New();
unmixer->SetInput(readerImage->GetOutput());
unmixer->SetMaxIteration(maxIter);
//unmixer->SetNumberOfThreads(1);
unmixer->SetEndmembersMatrix(endMember2Matrix->GetMatrix());
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outputImage);
writer->SetInput(unmixer->GetOutput());
writer->SetNumberOfDivisionsStrippedStreaming(10);
otb::StandardWriterWatcher w4(writer, unmixer,"CLSPSTOUnmixingImageFilter");
writer->Update();
return EXIT_SUCCESS;
}
/*=========================================================================
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 "otbFCLSUnmixingImageFilter.h"
#include "otbVectorImage.h"
#include "otbImageFileReader.h"
#include "otbStreamingImageFileWriter.h"
#include "otbVectorImageToMatrixImageFilter.h"
#include "otbStandardWriterWatcher.h"
const unsigned int Dimension = 2;
typedef double PixelType;
typedef otb::VectorImage<PixelType, Dimension> ImageType;
typedef otb::ImageFileReader<ImageType> ReaderType;
typedef otb::FCLSUnmixingImageFilter<ImageType, ImageType, PixelType> UnmixingImageFilterType;
typedef otb::VectorImageToMatrixImageFilter<ImageType> VectorImageToMatrixImageFilterType;
typedef otb::StreamingImageFileWriter<ImageType> WriterType;
int otbFCLSUnmixingImageFilterNewTest(int argc, char * argv[])
{
UnmixingImageFilterType::Pointer filter = UnmixingImageFilterType::New();
std::cout << filter << std::endl;
return EXIT_SUCCESS;
}
int otbFCLSUnmixingImageFilterTest(int argc, char * argv[])
{
const char * inputImage = argv[1];
const char * inputEndmembers = argv[2];
const char * outputImage = argv[3];
int maxIter = atoi(argv[4]);
ReaderType::Pointer readerImage = ReaderType::New();
readerImage->SetFileName(inputImage);
ReaderType::Pointer readerEndMembers = ReaderType::New();
readerEndMembers->SetFileName(inputEndmembers);
VectorImageToMatrixImageFilterType::Pointer endMember2Matrix = VectorImageToMatrixImageFilterType::New();
endMember2Matrix->SetInput(readerEndMembers->GetOutput());
endMember2Matrix->Update();
typedef VectorImageToMatrixImageFilterType::MatrixType MatrixType;
MatrixType endMembers = endMember2Matrix->GetMatrix();
MatrixType pinv = vnl_matrix_inverse<PixelType>(endMembers);
UnmixingImageFilterType::Pointer unmixer = UnmixingImageFilterType::New();
unmixer->SetInput(readerImage->GetOutput());
unmixer->SetMaxIteration(maxIter);
//unmixer->SetNumberOfThreads(1);
unmixer->SetEndmembersMatrix(endMember2Matrix->GetMatrix());
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outputImage);
writer->SetInput(unmixer->GetOutput());
writer->SetNumberOfDivisionsStrippedStreaming(10);
otb::StandardWriterWatcher w4(writer, unmixer,"FC