Skip to content
Snippets Groups Projects
Commit 58e31bd3 authored by Antoine Regimbeau's avatar Antoine Regimbeau
Browse files

Remove false NCLS filter from OTB

parent bac38cc5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Modules/Hyperspectral/Unmixing/include/otbNCLSUnmixingImageFilter.h deleted 100644 → 0
+ 0
− 131
View file @ bac38cc5
/*
* Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
*
* This file is part of Orfeo Toolbox
*
* https://www.orfeo-toolbox.org/
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef otbNCLSUnmixingImageFilter_h
#define otbNCLSUnmixingImageFilter_h
#include "itkMacro.h"
#include "itkNumericTraits.h"
#include "otbFunctorImageFilter.h"
#include "vnl/algo/vnl_svd.h"
#include <boost/shared_ptr.hpp>
namespace otb
{
namespace Functor {
/** \class NCLSUnmixingFunctor
*
* \brief Performs fully constrained least squares on a pixel
*
* \sa NCLSUnmixingImageFilter
*
* \ingroup OTBUnmixing
*/
template<class TInput, class TOutput, class TPrecision>
class NCLSUnmixingFunctor
{
public:
typedef NCLSUnmixingFunctor<TInput, TOutput, TPrecision> Self;
typedef TInput InputType;
typedef TOutput OutputType;
typedef TPrecision PrecisionType;
typedef vnl_vector<PrecisionType> VectorType;
typedef vnl_matrix<PrecisionType> MatrixType;
NCLSUnmixingFunctor();
virtual ~NCLSUnmixingFunctor() = default;
size_t OutputSize(const std::array<size_t, 1>& nbBands) 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_Ut;
MatrixType m_UtUinv;
SVDPointerType m_Svd; // SVD of U
unsigned int m_OutputSize;
unsigned int m_MaxIteration;
};
}
/** \class NCLSUnmixingImageFilter
*
* \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 \f$A\f$, it solves, for each pixel \f$p\f$, the system
* \f$A \cdot x = p\f$ in the least square sense, with additional 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
*
* \ingroup OTBUnmixing
*/
template <typename TInputImage, typename TOutputImage, typename TPrecision>
using NCLSUnmixingImageFilter = FunctorImageFilter<Functor::NCLSUnmixingFunctor<typename TInputImage::PixelType, typename TOutputImage::PixelType, TPrecision>>;
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbNCLSUnmixingImageFilter.hxx"
#endif
#endif
Modules/Hyperspectral/Unmixing/include/otbNCLSUnmixingImageFilter.hxx deleted 100644 → 0
+ 0
− 112
View file @ bac38cc5
/*
* Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
*
* This file is part of Orfeo Toolbox
*
* https://www.orfeo-toolbox.org/
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef otbNCLSUnmixingImageFilter_hxx
#define otbNCLSUnmixingImageFilter_hxx
#include "otbNCLSUnmixingImageFilter.h"
#include <algorithm>
namespace otb
{
namespace Functor
{
template <class TInput, class TOutput, class TPrecision>
NCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::NCLSUnmixingFunctor()
: m_OutputSize(0),
m_MaxIteration(100)
{
}
template <class TInput, class TOutput, class TPrecision>
size_t
NCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::OutputSize(const std::array<size_t,1> &) const
{
return m_OutputSize;
}
template <class TInput, class TOutput, class TPrecision>
void
NCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::SetEndmembersMatrix(const MatrixType& U)
{
m_U = U;
m_Ut = m_U.transpose();
m_UtUinv = SVDType(m_Ut * m_U).inverse();
m_OutputSize = m_U.cols();
m_Svd.reset( new SVDType(m_U) );
}
template <class TInput, class TOutput, class TPrecision>
const typename NCLSUnmixingFunctor<TInput, TOutput, TPrecision>::MatrixType&
NCLSUnmixingFunctor<TInput, TOutput, TPrecision>
::GetEndmembersMatrix() const
{
return m_U;
}
template <class TInput, class TOutput, class TPrecision>
typename NCLSUnmixingFunctor<TInput, TOutput, TPrecision>::OutputType
NCLSUnmixingFunctor<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;
// Apply NCLS iterations
VectorType lambda(nbEndmembers);
VectorType nclsVector = uclsVector;
VectorType correction(uclsVector.size());
for (unsigned int i = 0; i < m_MaxIteration; ++i)
{
// Error in original paper : divergence
// lambda = m_Ut * (inVector - m_U * nclsVector);
lambda = m_Ut * (m_U * nclsVector - inVector);
correction = m_UtUinv * lambda;
nclsVector -= correction;
}
OutputType out(nclsVector.size());
for (unsigned int i = 0; i < out.GetSize(); ++i )
{
out[i] = nclsVector[i];
}
return out;
}
} // end namespace Functor
} // end namespace otb
#endif
Modules/Hyperspectral/Unmixing/test/CMakeLists.txt
+ 0
− 11
View file @ 58e31bd3
......@@ -23,7 +23,6 @@ otb_module_test()
set(OTBUnmixingTests
otbUnmixingTestDriver.cxx
otbMDMDNMFImageFilter.cxx
otbNCLSUnmixingImageFilter.cxx
otbISRAUnmixingImageFilter.cxx
otbUnConstrainedLeastSquareImageFilter.cxx
otbSparseUnmixingImageFilter.cxx
......@@ -51,16 +50,6 @@ otb_add_test(NAME hyTvMDMDNMFImageFilterTest COMMAND otbUnmixingTestDriver
${TEMP}/hyTvMDMDNMFImageFilterTest.tif
100)
otb_add_test(NAME hyTvNCLSUnmixingImageFilterTest COMMAND otbUnmixingTestDriver
--compare-image ${EPSILON_9}
${BASELINE}/hyTvNCLSImageFilterTest.tif
${TEMP}/hyTvNCLSImageFilterTest.tif
otbNCLSUnmixingImageFilterTest
${INPUTDATA}/Hyperspectral/synthetic/hsi_cube.tif
${INPUTDATA}/Hyperspectral/synthetic/endmembers.tif
${TEMP}/hyTvNCLSImageFilterTest.tif
10)
otb_add_test(NAME hyTvISRAUnmixingImageFilterTest COMMAND otbUnmixingTestDriver
--compare-image ${EPSILON_9}
${BASELINE}/hyTvISRAUnmixingImageFilterTest.tif
......
Modules/Hyperspectral/Unmixing/test/otbNCLSUnmixingImageFilter.cxx deleted 100644 → 0
+ 0
− 76
View file @ bac38cc5
/*
* Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
*
* This file is part of Orfeo Toolbox
*
* https://www.orfeo-toolbox.org/
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include "otbNCLSUnmixingImageFilter.h"
#include "otbVectorImage.h"
#include "otbImageFileReader.h"
#include "otbImageFileWriter.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::NCLSUnmixingImageFilter<ImageType, ImageType, PixelType> UnmixingImageFilterType;
typedef otb::VectorImageToMatrixImageFilter<ImageType> VectorImageToMatrixImageFilterType;
typedef otb::ImageFileWriter<ImageType> WriterType;
int otbNCLSUnmixingImageFilterTest(int itkNotUsed(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->GetModifiableFunctor().SetMaxIteration(maxIter);
//unmixer->SetNumberOfThreads(1);
unmixer->GetModifiableFunctor().SetEndmembersMatrix(endMember2Matrix->GetMatrix());
WriterType::Pointer writer = WriterType::New();
writer->SetFileName(outputImage);
writer->SetInput(unmixer->GetOutput());
writer->SetNumberOfDivisionsStrippedStreaming(10);
otb::StandardWriterWatcher w4(writer, unmixer,"NCLSUnmixingImageFilter");
writer->Update();
return EXIT_SUCCESS;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment