Commit d49c65be authored by Antoine Regimbeau's avatar Antoine Regimbeau

Merge branch 'develop' into cmake_undefined_var

parents 79856999 93111e34
......@@ -23,8 +23,6 @@
#include "otbUnConstrainedLeastSquareImageFilter.h"
#include "otbISRAUnmixingImageFilter.h"
#include "otbNCLSUnmixingImageFilter.h"
//#include "otbFCLSUnmixingImageFilter.h"
#include "otbMDMDNMFImageFilter.h"
......@@ -34,8 +32,6 @@ namespace Wrapper
{
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::MDMDNMFImageFilter<DoubleVectorImageType, DoubleVectorImageType> MDMDNMFUnmixingFilterType;
typedef otb::VectorImageToMatrixImageFilter<DoubleVectorImageType> VectorImageToMatrixImageFilterType;
......@@ -65,12 +61,12 @@ enum UnMixingMethod
{
UnMixingMethod_UCLS,
//UnMixingMethod_FCLS,
UnMixingMethod_NCLS,
// UnMixingMethod_NCLS,
UnMixingMethod_ISRA,
UnMixingMethod_MDMDNMF,
};
const char* UnMixingMethodNames [] = { "UCLS", "FCLS", "NCLS", "ISRA", "MDMDNMF", };
const char* UnMixingMethodNames [] = { "UCLS", "ISRA", "MDMDNMF", };
class HyperspectralUnmixing : public Application
......@@ -104,9 +100,7 @@ private:
"be estimated using the VertexComponentAnalysis application.\n\n"
"The application allows estimating the abundance maps with several algorithms:\n\n"
"* Unconstrained Least Square (ucls)\n"
// "* Fully Constrained Least Square (fcls)\n"
"* Image Space Reconstruction Algorithm (isra)\n"
"* Non-negative constrained\n"
"* Least Square (ncls)\n"
"* Minimum Dispersion Constrained Non Negative Matrix Factorization (MDMDNMF)."
);
......@@ -135,12 +129,6 @@ private:
AddChoice("ua.ucls", "UCLS");
SetParameterDescription("ua.ucls", "Unconstrained 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");
AddChoice("ua.isra", "ISRA");
SetParameterDescription("ua.isra", "Image Space Reconstruction Algorithm");
......@@ -217,36 +205,6 @@ private:
}
break;
case UnMixingMethod_NCLS:
{
otbAppLogINFO("NCLS Unmixing");
NCLSUnmixingFilterType::Pointer unmixer =
NCLSUnmixingFilterType::New();
unmixer->SetInput(inputImage);
unmixer->GetModifiableFunctor().SetEndmembersMatrix(endMembersMatrix);
abundanceMap = unmixer->GetOutput();
m_ProcessObjects.push_back(unmixer.GetPointer());
}
break;
/*
case UnMixingMethod_FCLS:
{
otbAppLogINFO("FCLS Unmixing");
FCLSUnmixingFilterType::Pointer unmixer =
FCLSUnmixingFilterType::New();
unmixer->SetInput(inputImage);
unmixer->SetEndmembersMatrix(endMembersMatrix);
abundanceMap = unmixer->GetOutput();
m_ProcessObjects.push_back(unmixer.GetPointer());
}
break;
*/
case UnMixingMethod_MDMDNMF:
{
otbAppLogINFO("MDMD-NMF Unmixing");
......
......@@ -30,17 +30,6 @@ otb_test_application(NAME apTvHyHyperspectralUnmixing_UCLS
${BASELINE}/apTvHyHyperspectralUnmixing_UCLS.tif
${TEMP}/apTvHyHyperspectralUnmixing_UCLS.tif)
otb_test_application(NAME apTvHyHyperspectralUnmixing_NCLS
APP HyperspectralUnmixing
OPTIONS -in ${INPUTDATA}/Hyperspectral/synthetic/hsi_cube.tif
-ie ${INPUTDATA}/Hyperspectral/synthetic/endmembers.tif
-out ${TEMP}/apTvHyHyperspectralUnmixing_NCLS.tif double
-ua ncls
VALID --compare-image ${EPSILON_9}
${BASELINE}/apTvHyHyperspectralUnmixing_NCLS.tif
${TEMP}/apTvHyHyperspectralUnmixing_NCLS.tif)
#----------- VertexComponentAnalysis TESTS ----------------
otb_test_application(NAME apTvHyVertexComponentAnalysis
APP VertexComponentAnalysis
......
......@@ -230,24 +230,9 @@ StreamingWarpImageFilter<TInputImage, TOutputImage, TDisplacementField>
{
inputPtr->SetRequestedRegion(inputRequestedRegion);
}
else
{
inputFinalSize.Fill(0);
inputRequestedRegion.SetSize(inputFinalSize);
inputFinalIndex.Fill(0);
inputRequestedRegion.SetIndex(inputFinalIndex);
// store what we tried to request (prior to trying to crop)
inputPtr->SetRequestedRegion(inputRequestedRegion);
// // build an exception
// itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
// e.SetLocation(ITK_LOCATION);
// e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
// e.SetDataObject(inputPtr);
// throw e;
}
// Here we do not throw an exception, we just do nothing since
// resampling filter can legitimately ask for a region outside of
// input's largest possible region.
}
......
/*
* 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
/*
* 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
......@@ -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
......
/*
* 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;
}
......@@ -24,7 +24,6 @@ void RegisterTests()
{
REGISTER_TEST(otbMDMDNMFImageFilterTest);
REGISTER_TEST(otbMDMDNMFImageFilterTest2);
REGISTER_TEST(otbNCLSUnmixingImageFilterTest);
REGISTER_TEST(otbISRAUnmixingImageFilterTest);
REGISTER_TEST(otbUnConstrainedLeastSquareImageFilterTest);
REGISTER_TEST(otbSparseUnmixingImageFilterTest);
......
Markdown is supported
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