Commit f15c213f authored by Julien Malik's avatar Julien Malik

ENH: remove partially implemented classes and add missing tests

parent 4f0e9cd5
/*=========================================================================
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 __otbCLSPSTOUnmixingImageFilter_h
#define __otbCLSPSTOUnmixingImageFilter_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 CLSPSTOUnmixingFunctor
*
* \brief TODO
*
*/
template<class TInput, class TOutput, class TPrecision>
class CLSPSTOUnmixingFunctor
{
public:
typedef CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision> Self;
typedef TInput InputType;
typedef TOutput OutputType;
typedef TPrecision PrecisionType;
typedef vnl_vector<PrecisionType> VectorType;
typedef vnl_matrix<PrecisionType> MatrixType;
CLSPSTOUnmixingFunctor();
virtual ~CLSPSTOUnmixingFunctor();
unsigned int GetOutputSize() const;
bool operator !=(const CLSPSTOUnmixingFunctor& other) const;
bool operator ==(const CLSPSTOUnmixingFunctor& 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;
SVDPointerType m_Svd; // SVD of U
unsigned int m_OutputSize;
unsigned int m_MaxIteration;
};
}
/** \class CLSPSTOUnmixingImageFilter
*
* \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 CLSPSTOUnmixingImageFilter :
public otb::UnaryFunctorImageFilter<TInputImage, TOutputImage,
Functor::CLSPSTOUnmixingFunctor<typename TInputImage::PixelType,
typename TOutputImage::PixelType, TPrecision> >
{
public:
/** Standard class typedefs. */
typedef CLSPSTOUnmixingImageFilter Self;
typedef otb::UnaryFunctorImageFilter
<TInputImage,
TOutputImage,
Functor::CLSPSTOUnmixingFunctor<
typename TInputImage::PixelType,
typename TOutputImage::PixelType,
TPrecision>
> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef Functor::CLSPSTOUnmixingFunctor<
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(CLSPSTOUnmixingImageFilter, 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:
CLSPSTOUnmixingImageFilter();
virtual ~CLSPSTOUnmixingImageFilter();
void PrintSelf(std::ostream& os, itk::Indent indent) const;
private:
CLSPSTOUnmixingImageFilter(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbCLSPSTOUnmixingImageFilter.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 __otbCLSPSTOUnmixingImageFilter_txx
#define __otbCLSPSTOUnmixingImageFilter_txx
#include "otbCLSPSTOUnmixingImageFilter.h"
#include <algorithm>
namespace otb
{
namespace Functor
{
template <class TInput, class TOutput, class TPrecision>
CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>
::CLSPSTOUnmixingFunctor()
: m_OutputSize(0),
m_MaxIteration(100)
{
}
template <class TInput, class TOutput, class TPrecision>
CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>
::~CLSPSTOUnmixingFunctor()
{
}
template <class TInput, class TOutput, class TPrecision>
unsigned int
CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>
::GetOutputSize() const
{
return m_OutputSize;
}
template <class TInput, class TOutput, class TPrecision>
bool
CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>
::operator !=(const Self& other) const
{
return true;
}
template <class TInput, class TOutput, class TPrecision>
bool
CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>
::operator ==(const Self& other) const
{
return false;
}
template <class TInput, class TOutput, class TPrecision>
void
CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>
::SetEndmembersMatrix(const MatrixType& U)
{
m_U = U;
m_OutputSize = m_U.cols();
m_Svd.reset( new SVDType(m_U) );
}
template <class TInput, class TOutput, class TPrecision>
const typename CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>::MatrixType&
CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>
::GetEndmembersMatrix() const
{
return m_U;
}
template <class TInput, class TOutput, class TPrecision>
typename CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>::OutputType
CLSPSTOUnmixingFunctor<TInput, TOutput, TPrecision>
::operator ()(const InputType& in) const
{
// TODO : support different types between input and output ?
VectorType inVector(in.Size());
for (int i = 0; i < in.GetSize(); i++ )
{
inVector[i] = in[i];
}
// Initialize with Unconstrained Least Square solution
VectorType outVector = m_Svd->solve(inVector);
unsigned int nbEndmembers = m_OutputSize;
unsigned int nbBands = in.Size();
// Apply CLSPSTO iterations
for (int i = 0; i < m_MaxIteration; ++i)
{
// NOT IMPLEMENTED
}
OutputType out(outVector.size());
for (int i = 0; i < out.GetSize(); i++ )
{
out[i] = outVector[i];
}
return out;
}
}
template <class TInputImage, class TOutputImage, class TPrecision>
CLSPSTOUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::CLSPSTOUnmixingImageFilter()
{
}
template <class TInputImage, class TOutputImage, class TPrecision>
CLSPSTOUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::~CLSPSTOUnmixingImageFilter()
{
}
template <class TInputImage, class TOutputImage, class TPrecision>
void
CLSPSTOUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::SetEndmembersMatrix(const MatrixType& m)
{
this->GetFunctor().SetEndmembersMatrix(m);
this->Modified();
}
template <class TInputImage, class TOutputImage, class TPrecision>
const typename CLSPSTOUnmixingImageFilter<TInputImage,TOutputImage,TPrecision>::MatrixType&
CLSPSTOUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::GetEndmembersMatrix() const
{
return this->GetFunctor().GetEndmembersMatrix();
}
template <class TInputImage, class TOutputImage, class TPrecision>
void
CLSPSTOUnmixingImageFilter<TInputImage, TOutputImage, TPrecision>
::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf(os, indent);
}
} // end namespace
#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 __otbFullyConstrainedLeastSquareImageFilter_h
#define __otbFullyConstrainedLeastSquareImageFilter_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 FullyConstrainedLeastSquareFunctor
*
* \brief TODO
*
*/
template<class TInput, class TOutput, class TPrecision>
class FullyConstrainedLeastSquareFunctor
{
public:
typedef FullyConstrainedLeastSquareFunctor<TInput, TOutput, TPrecision> Self;
typedef TInput InputType;
typedef TOutput OutputType;
typedef TPrecision PrecisionType;
typedef vnl_vector<PrecisionType> VectorType;
typedef vnl_matrix<PrecisionType> MatrixType;
FullyConstrainedLeastSquareFunctor();
virtual ~FullyConstrainedLeastSquareFunctor();
unsigned int GetOutputSize() const;
bool operator !=(const FullyConstrainedLeastSquareFunctor& other) const;
bool operator ==(const FullyConstrainedLeastSquareFunctor& other) const;
void SetMatrix(const MatrixType& U);
OutputType operator ()(const InputType& in) const;
private:
static bool IsNonNegative(PrecisionType val)
{
return val >= 0;
}
typedef vnl_svd<PrecisionType> SVDType;
MatrixType m_U;
unsigned int m_OutputSize;
};
}
/** \class FullyConstrainedLeastSquareImageFilter
*
* \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 FullyConstrainedLeastSquareImageFilter :
public otb::UnaryFunctorImageFilter<TInputImage, TOutputImage,
Functor::FullyConstrainedLeastSquareFunctor<typename TInputImage::PixelType,
typename TOutputImage::PixelType, TPrecision> >
{
public:
/** Standard class typedefs. */
typedef FullyConstrainedLeastSquareImageFilter Self;
typedef otb::UnaryFunctorImageFilter
<TInputImage,
TOutputImage,
Functor::FullyConstrainedLeastSquareFunctor<
typename TInputImage::PixelType,
typename TOutputImage::PixelType,
TPrecision>
> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef Functor::FullyConstrainedLeastSquareFunctor<
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(FullyConstrainedLeastSquareImageFilter, otb::UnaryFunctorImageFilter);
/** Pixel types. */
typedef typename TInputImage::PixelType InputPixelType;
typedef typename TOutputImage::PixelType OutputPixelType;
void SetMatrix(const MatrixType& m);
protected:
FullyConstrainedLeastSquareImageFilter();
virtual ~FullyConstrainedLeastSquareImageFilter();
void PrintSelf(std::ostream& os, itk::Indent indent) const;
private:
FullyConstrainedLeastSquareImageFilter(const Self &); //purposely not implemented
void operator =(const Self&); //purposely not implemented
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbFullyConstrainedLeastSquareImageFilter.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 __otbFullyConstrainedLeastSquareImageFilter_txx
#define __otbFullyConstrainedLeastSquareImageFilter_txx
#include "otbFullyConstrainedLeastSquareImageFilter.h"
#include <algorithm>
namespace otb
{
namespace Functor
{
template <class TInput, class TOutput, class TPrecision>
FullyConstrainedLeastSquareFunctor<TInput, TOutput, TPrecision>
::FullyConstrainedLeastSquareFunctor()
: m_OutputSize(0)
{
}
template <class TInput, class TOutput, class TPrecision>
FullyConstrainedLeastSquareFunctor<TInput, TOutput, TPrecision>
::~FullyConstrainedLeastSquareFunctor()
{
}
template <class TInput, class TOutput, class TPrecision>
unsigned int
FullyConstrainedLeastSquareFunctor<TInput, TOutput, TPrecision>
::GetOutputSize() const
{
return m_OutputSize;
}
template <class TInput, class TOutput, class TPrecision>
bool
FullyConstrainedLeastSquareFunctor<TInput, TOutput, TPrecision>
::operator !=(const Self& other) const
{
return true;
}
template <class TInput, class TOutput, class TPrecision>
bool
FullyConstrainedLeastSquareFunctor<TInput, TOutput, TPrecision>
::operator ==(const Self& other) const
{
return false;
}
template <class TInput, class TOutput, class TPrecision>
void
FullyConstrainedLeastSquareFunctor<TInput, TOutput, TPrecision>
::SetMatrix(const MatrixType& U)
{
m_U = U;
m_OutputSize = m_U.cols();
}
template <class TInput, class TOutput, class TPrecision>
typename FullyConstrainedLeastSquareFunctor<TInput, TOutput, TPrecision>::OutputType
FullyConstrainedLeastSquareFunctor<TInput, TOutput, TPrecision>
::operator ()(const InputType& in) const
{
// TODO : support different types between input and output ?
VectorType inVector(in.GetDataPointer(), in.Size());
VectorType out(m_OutputSize);
MatrixType U = m_U;
unsigned int count = m_U.cols();
typedef std::vector<unsigned int> IndiciesType;
IndiciesType indicies;
for (unsigned int i = 0; i < m_OutputSize; ++i)
{
indicies.push_back(i);
}
for (count = m_U.cols(); count > 0; --count)
{
// First, solve the unconstrained least square system
VectorType als_hat = SVDType(U).solve(inVector);
MatrixType Z(1, count);
Z.fill(itk::NumericTraits<PrecisionType>::One);
MatrixType Zt = Z.transpose();
MatrixType UtU = U.transpose() * U;
MatrixType UtUinv = SVDType(UtU).inverse();
MatrixType S = UtUinv * Zt;
// Compute correction term needed for full additivity constraint
PrecisionType correctionTerm = (S * SVDType( Z * S ).inverse())(0,0) * ((Z * als_hat)[0] - 1);
// Remove the correction term from the unconstrained solution
VectorType afcls_hat = als_hat - correctionTerm;
// If everybody is positive, it is finished
if ( std::count_if(afcls_hat.begin(), afcls_hat.end(), Self::IsNonNegative) == count )
{
out.fill(0);
for (unsigned int j = 0; j < indicies.size(); ++j)
{
out[ indicies[j] ] = afcls_hat[ j ];
}
break;
}
// Multiply negative elements by their counterpart in the S vector
// and find max(abs(afcls_hat)) for indicies where elements are negative
unsigned int maxIdx = 0;
PrecisionType maxi = itk::NumericTraits<PrecisionType>::min();
for (unsigned int j = 0; j < afcls_hat.size(); ++j)
{
if (afcls_hat[j] < 0)