From dabc50b1d24e60870c1098d299cc9720f58fd23d Mon Sep 17 00:00:00 2001
From: Julien Malik <julien.malik@c-s.fr>
Date: Fri, 18 Feb 2011 08:36:59 +0100
Subject: [PATCH] ADD: add FullyConstrainedLeastSquare image filter
---
...tbFullyConstrainedLeastSquareImageFilter.h | 159 ++++++++++++++
...FullyConstrainedLeastSquareImageFilter.txx | 201 ++++++++++++++++++
Testing/CMakeLists.txt | 20 ++
...FullyConstrainedLeastSquareImageFilter.cxx | 75 +++++++
Testing/otbHyperTests1.cxx | 2 +
5 files changed, 457 insertions(+)
create mode 100644 Code/Hyperspectral/otbFullyConstrainedLeastSquareImageFilter.h
create mode 100644 Code/Hyperspectral/otbFullyConstrainedLeastSquareImageFilter.txx
create mode 100644 Testing/otbFullyConstrainedLeastSquareImageFilter.cxx
diff --git a/Code/Hyperspectral/otbFullyConstrainedLeastSquareImageFilter.h b/Code/Hyperspectral/otbFullyConstrainedLeastSquareImageFilter.h
new file mode 100644
index 0000000000..e8be0dcc3b
--- /dev/null
+++ b/Code/Hyperspectral/otbFullyConstrainedLeastSquareImageFilter.h
@@ -0,0 +1,159 @@
+/*=========================================================================
+
+ 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 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
+
diff --git a/Code/Hyperspectral/otbFullyConstrainedLeastSquareImageFilter.txx b/Code/Hyperspectral/otbFullyConstrainedLeastSquareImageFilter.txx
new file mode 100644
index 0000000000..7c457a2f12
--- /dev/null
+++ b/Code/Hyperspectral/otbFullyConstrainedLeastSquareImageFilter.txx
@@ -0,0 +1,201 @@
+/*=========================================================================
+
+ 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)
+ {
+ afcls_hat[j] = afcls_hat[j] / S(j,0);
+
+ if (vcl_abs(afcls_hat[j]) > maxi)
+ {
+ maxi = vcl_abs(afcls_hat[j]);
+ maxIdx = j;
+ }
+ }
+ }
+
+ // Remove column maxIdx from U
+ MatrixType Unew(U.rows(), U.cols() - 1);
+ for (unsigned int Ucol = 0, UcolNew = 0; Ucol < U.cols(); Ucol++)
+ {
+ if (Ucol != maxIdx)
+ {
+ Unew.set_column(UcolNew, U.get_column(Ucol));
+ UcolNew++;
+ }
+ }
+
+ // Remove maxIdx-th element from indicies
+ IndiciesType::iterator itToRemove = indicies.begin() + maxIdx;
+ indicies.erase(itToRemove);
+
+ // Update U
+ U = Unew;
+ }
+
+ // Convert to OutputType (vnl_vector -> VariableLengthVector)
+ return OutputType(out.data_block(), out.size());
+}
+
+}
+
+template <class TInputImage, class TOutputImage, class TPrecision>
+FullyConstrainedLeastSquareImageFilter<TInputImage, TOutputImage, TPrecision>
+::FullyConstrainedLeastSquareImageFilter()
+{
+}
+
+template <class TInputImage, class TOutputImage, class TPrecision>
+FullyConstrainedLeastSquareImageFilter<TInputImage, TOutputImage, TPrecision>
+::~FullyConstrainedLeastSquareImageFilter()
+{
+}
+
+template <class TInputImage, class TOutputImage, class TPrecision>
+void
+FullyConstrainedLeastSquareImageFilter<TInputImage, TOutputImage, TPrecision>
+::SetMatrix(const MatrixType& m)
+{
+ this->GetFunctor().SetMatrix(m);
+ this->Modified();
+}
+
+template <class TInputImage, class TOutputImage, class TPrecision>
+void
+FullyConstrainedLeastSquareImageFilter<TInputImage, TOutputImage, TPrecision>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+ Superclass::PrintSelf(os, indent);
+}
+
+} // end namespace
+
+#endif
diff --git a/Testing/CMakeLists.txt b/Testing/CMakeLists.txt
index a8cff678b4..b34eaa27bf 100644
--- a/Testing/CMakeLists.txt
+++ b/Testing/CMakeLists.txt
@@ -10,6 +10,7 @@ SET( otbHyperTests1_SRC
otbStreamingStatisticsVectorImageFilter2.cxx
otbMatrixMultiplyImageFilter.cxx
otbUnConstrainedLeastSquareImageFilter.cxx
+ otbFullyConstrainedLeastSquareImageFilter.cxx
otbVectorImageToMatrixImageFilter.cxx
otbProjectiveProjection.cxx
otbVCA.cxx
@@ -99,6 +100,25 @@ ADD_TEST(bfTvUnConstrainedLeastSquareImageFilter
${TEMP}/bfUnConstrainedLeastSquareImageFilter.tif)
+ADD_TEST(bfTuFullyConstrainedLeastSquareImageFilterNew
+ ${TESTEXE_DIR}/otbHyperTests1
+ otbFullyConstrainedLeastSquareImageFilterNewTest)
+
+ADD_TEST(bfTvFullyConstrainedLeastSquareImageFilter_SYNTHETIC
+ ${TESTEXE_DIR}/otbHyperTests1
+ otbFullyConstrainedLeastSquareImageFilterTest
+ ${DATA}/SYNTHETIC100/hsi_cube.tif
+ ${DATA}/SYNTHETIC100/endmembers.tif
+ ${TEMP}/bfFullyConstrainedLeastSquareImageFilter.tif)
+
+
+ADD_TEST(bfTvFullyConstrainedLeastSquareImageFilter_AVIRIS512
+ ${TESTEXE_DIR}/otbHyperTests1
+ otbFullyConstrainedLeastSquareImageFilterTest
+ ${DATA}/AVIRIS/extract512/f970619t01p02_r02_sc01_goodbands.hdr
+ ${TEMP}/aviris_vcaimagefilter_13.hdr
+ ${TEMP}/bfFullyConstrainedLeastSquareImageFilter_AVIRIS512.tif)
+
ADD_TEST(bfTuVectorImageToMatrixNew
diff --git a/Testing/otbFullyConstrainedLeastSquareImageFilter.cxx b/Testing/otbFullyConstrainedLeastSquareImageFilter.cxx
new file mode 100644
index 0000000000..2010f3c9ce
--- /dev/null
+++ b/Testing/otbFullyConstrainedLeastSquareImageFilter.cxx
@@ -0,0 +1,75 @@
+/*=========================================================================
+
+ 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 "otbFullyConstrainedLeastSquareImageFilter.h"
+
+#include "otbVectorImage.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+#include "otbVectorImageToMatrixImageFilter.h"
+
+const unsigned int Dimension = 2;
+typedef double PixelType;
+
+typedef otb::VectorImage<PixelType, Dimension> ImageType;
+typedef otb::ImageFileReader<ImageType> ReaderType;
+typedef otb::FullyConstrainedLeastSquareImageFilter<ImageType,ImageType,double> FullyConstrainedLeastSquareSolverType;
+typedef otb::VectorImageToMatrixImageFilter<ImageType> VectorImageToMatrixImageFilterType;
+typedef otb::ImageFileWriter<ImageType> WriterType;
+
+int otbFullyConstrainedLeastSquareImageFilterNewTest(int argc, char * argv[])
+{
+ FullyConstrainedLeastSquareSolverType::Pointer filter = FullyConstrainedLeastSquareSolverType::New();
+ std::cout << filter << std::endl;
+ return EXIT_SUCCESS;
+}
+
+int otbFullyConstrainedLeastSquareImageFilterTest(int argc, char * argv[])
+{
+ const char * inputImage = argv[1];
+ const char * inputEndmembers = argv[2];
+ const char * outputImage = argv[3];
+
+ 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<double>(endMembers);
+
+ FullyConstrainedLeastSquareSolverType::Pointer unmixer = \
+ FullyConstrainedLeastSquareSolverType::New();
+
+ unmixer->SetInput(readerImage->GetOutput());
+ unmixer->SetMatrix(endMember2Matrix->GetMatrix());
+
+ unmixer->Update();
+
+ WriterType::Pointer writer = WriterType::New();
+ writer->SetFileName(outputImage);
+ writer->SetInput(unmixer->GetOutput());
+ writer->Update();
+
+ return EXIT_SUCCESS;
+}
diff --git a/Testing/otbHyperTests1.cxx b/Testing/otbHyperTests1.cxx
index 4cf869af24..d06c4c65ff 100644
--- a/Testing/otbHyperTests1.cxx
+++ b/Testing/otbHyperTests1.cxx
@@ -32,6 +32,8 @@ void RegisterTests()
REGISTER_TEST(otbMatrixMultiplyImageFilterTest);
REGISTER_TEST(otbUnConstrainedLeastSquareImageFilterNewTest);
REGISTER_TEST(otbUnConstrainedLeastSquareImageFilterTest);
+ REGISTER_TEST(otbFullyConstrainedLeastSquareImageFilterNewTest);
+ REGISTER_TEST(otbFullyConstrainedLeastSquareImageFilterTest);
REGISTER_TEST(otbVectorImageToMatrixNewTest);
REGISTER_TEST(otbVectorImageToMatrixTest);
REGISTER_TEST(vahineVCA);
--
GitLab