Commit ff2b2c5d authored by Cédric Traizet's avatar Cédric Traizet

Merge branch 'spectral_angle_classif' into 'develop'

Spectral angle classification

Closes #2012

See merge request !708
parents e0e030d5 f223a1df
Pipeline #4881 passed with stages
in 39 minutes and 12 seconds
......@@ -106,6 +106,68 @@ image is shown below:
Resulting unmixed image, here the first three bands are displayed.
Classification
--------------
A common task in hyperspectral image processing is to classify a data cube. For example
one might want to match pixels to reference endmembers. The Spectral Angle Mapper (SAM) algorithm
[1] does that by computing a spectral measure between each pixel and the references.
The spectral angle of a pixel `x` with a reference pixel `r` is defined by :
.. math::
sam[x, r] = cos^{-1}(\frac{<x,r>}{\|x\| * \|r\| } )
where `<x,r>` denotes the scalar product between x and r.
This is also called the spectral angle between `x` and `r`.
In SAM classification the spectral angle is computed for each pixel with a set of reference pixels,
the class associated with the pixel is the index of the reference pixel that has the lowest spectral angle.
::
otbcli_SpectralAngleClassification -in inputImage.tif
-ie endmembers.tif
-out classification.tif
-mode sam
The result of the classification is shown below :
.. figure:: ../Art/HyperspectralImages/classification.png
:width: 70%
:align: center
Resulting unmixed classified image.
Another algorithm is available in this application based on spectral information divergence [1],
if we define the probability mass function p of a pixel x by :
.. math::
x = [x_1, x_2 ,..., x_L ]^T \\
p = [p_1, p_2 ,..., p_L ]^T \\
p_i = \frac{x_i}{\sum_{j=1}^L x_j}
The spectral information divergence between a pixel `x` and a reference `r` is defined by :
.. math::
sid[x, r] = \sum_{j=1}^{L} p_j *log(\frac{p_j}{q_j}) + \sum_{j=1}^{L} q_j * log(\frac{q_j}{p_j})
where p and q are respectively the probability mass function of x and r. As with the SAM algorithm,
the class associated with the pixel is the index of the reference pixel that has the lowest spectral information
divergence.
Note that the framework described in the
:ref:`classification recipe<classif>` can also be applied to hyperspectral data.
Anomaly detection
-----------------
......@@ -183,6 +245,10 @@ The value of the threshold depends on how sensitive the anomaly detector should
Left: Computed Rx score, right: detected anomalies (in red)
*[1] Du, Yingzi & Chang, Chein-I & Ren, Hsuan & Chang, Chein-Chi & Jensen, James & D'Amico, Francis. (2004).
New Hyperspectral Discrimination Measure for Spectral Characterization. Optical Engineering - OPT ENG. 43.*
.. _here: http://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes#Cuprite
.. _AVIRIS: https://aviris.jpl.nasa.gov/
.. _Pavia: http://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes#Pavia_University_scene
.. _classif:
Classification
==============
......
......@@ -38,3 +38,8 @@ otb_create_application(
SOURCES otbEndmemberNumberEstimation.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
otb_create_application(
NAME SpectralAngleClassification
SOURCES otbSpectralAngleClassification.cxx
LINK_LIBRARIES ${${otb-module}_LIBRARIES})
/*
* Copyright (C) 2005-2020 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 "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbSpectralAngleFunctor.h"
#include "otbFunctorImageFilter.h"
#include "otbSpectralInformationDivergenceFunctor.h"
namespace otb
{
namespace Wrapper
{
class SpectralAngleClassification : public Application
{
public:
/** @name Standard class typedefs
* @{
*/
using Self = SpectralAngleClassification;
using Superclass = Application;
using Pointer = itk::SmartPointer<Self>;
using ConstPointer = itk::SmartPointer<const Self>;
/** @} */
/** @name Standard macro
* @{
*/
itkNewMacro(Self);
itkTypeMacro(SpectralAngleClassification, otb::Application);
/** @} */
using ValueType = float;
using ImageType = otb::VectorImage<ValueType>;
using PixelType = ImageType::PixelType;
using SAMFilterType = otb::FunctorImageFilter<otb::Functor::SpectralAngleMapperFunctor<PixelType, PixelType, PixelType>>;
using SIDFilterType = otb::FunctorImageFilter<otb::Functor::SpectralInformationDivergenceFunctor<PixelType, PixelType, PixelType>>;
private:
void DoInit() override
{
SetName("SpectralAngleClassification");
SetDescription("Classifies an image using a spectral measure.");
// Documentation
SetDocLongDescription("This application computes a spectral measure on a hyperspectral data cube using a set of "
"reference pixels. The image is then classified by finding for each pixel the index of its closest endmember, "
"i.e. the endmember corresponding to the lowest measurement value. The output classification is "
"labeled from 1 to L, L being the number of endmembers. A threshold can be set for the "
"classification step, only measurement values lower than this threshold will be considered, other "
"will be classified as `background values` (default value of 0). \n\n"
"Two measures are available : the spectral angle mapper and the spectral information divergence. See [1] for details");
SetDocLimitations("In sid mode, the pixels of the input image and the input endmembers should be strictly positive. \n"
"The endmember image is fully loaded in memory.");
SetDocAuthors("OTB-Team");
SetDocSeeAlso("VertexComponentAnalysis \n"
"[1] Du, Yingzi & Chang, Chein-I & Ren, Hsuan & Chang, Chein-Chi & Jensen, James & D'Amico, Francis. (2004). "
"New Hyperspectral Discrimination Measure for Spectral Characterization. Optical Engineering - OPT ENG. 43."
" 1777-1786. 10.1117/1.1766301. ");
AddDocTag(Tags::Hyperspectral);
AddParameter(ParameterType_InputImage, "in", "Input Image Filename");
SetParameterDescription("in", "The hyperspectral data cube input");
AddParameter(ParameterType_InputImage, "ie", "Input endmembers");
SetParameterDescription("ie",
"The endmembers (estimated pure pixels) to "
"use for unmixing. Must be stored as a multispectral image, where "
"each pixel is interpreted as an endmember.");
AddParameter(ParameterType_OutputImage, "measure", "Output spectral angle values");
SetParameterDescription("measure",
"Output image containing for each pixel from the input image the computed measure relative to each endmember");
MandatoryOff("measure");
AddParameter(ParameterType_OutputImage, "out", "Output classified image");
SetParameterDescription("out",
"Output classified image, classified pixels are labeled from 1 to L, L being the number of endmember in the image.");
MandatoryOff("out");
AddParameter(ParameterType_Choice, "mode", "Measure used for classification");
SetParameterDescription("mode", "Measure used for classification");
MandatoryOff("mode");
AddChoice("mode.sam", "Spectral angle mapper");
SetParameterDescription("mode.sam", "Spectral angle mapper (SAM) measure.");
AddChoice("mode.sid", "Spectral information divergence");
SetParameterDescription("mode.sid", "Spectral information divergence (SID) measure. "
"Input pixel values should be strictly positive.");
AddParameter(ParameterType_Float, "threshold", "Classification threshold");
SetParameterDescription("threshold",
"Pixel with a measurement greater than this threshold relatively to "
"a reference pixel are not classified. The same threshold is used for all classes.");
MandatoryOff("threshold");
AddParameter(ParameterType_Int, "bv", "Background value");
SetParameterDescription("bv",
"Value of unclassified pixels in the classification image "
"(this parameter is only used if a threshold is set).");
MandatoryOff("bv");
SetDefaultParameterInt("bv", 0);
AddRAMParameter();
SetMultiWriting(true);
// Doc example parameter settings
SetDocExampleParameterValue("in", "cupriteSubHsi.tif");
SetDocExampleParameterValue("ie", "cupriteEndmembers.tif");
SetDocExampleParameterValue("out", "classification.tif");
SetDocExampleParameterValue("measure", "measure.tif");
SetDocExampleParameterValue("mode", "sam");
SetDocExampleParameterValue("threshold", "0.1");
SetOfficialDocLink();
}
void DoUpdateParameters() override
{
// Nothing to do here : all parameters are independent
}
void DoExecute() override
{
auto inputEndmemberImage = GetParameterImage("ie");
inputEndmemberImage->Update();
itk::ImageRegionConstIterator<ImageType> it(inputEndmemberImage, inputEndmemberImage->GetLargestPossibleRegion());
std::vector<PixelType> endmembers;
for (it.GoToBegin(); !it.IsAtEnd(); ++it)
{
endmembers.push_back(it.Get());
}
// Process objects to keep references on the actual filter instanciated below and its output
itk::LightObject::Pointer filter;
ImageType::Pointer filterOutput;
auto mode = GetParameterString("mode");
if (mode == "sam")
{
auto SAMFilter = SAMFilterType::New();
SAMFilter->GetModifiableFunctor().SetReferencePixels(endmembers);
SAMFilter->SetInput(GetParameterImage("in"));
filter = SAMFilter;
filterOutput = SAMFilter->GetOutput();
}
else if (mode == "sid")
{
auto SIDFilter = SIDFilterType::New();
SIDFilter->GetModifiableFunctor().SetReferencePixels(endmembers);
SIDFilter->SetInput(GetParameterImage("in"));
filter = SIDFilter;
filterOutput = SIDFilter->GetOutput();
}
if (HasValue("measure"))
{
SetParameterOutputImage("measure", filterOutput);
}
if (HasValue("out"))
{
auto threshold = HasValue("threshold") ? GetParameterFloat("threshold")
: std::numeric_limits<ValueType>::max();
auto bv = GetParameterInt("bv");
// This lambda return the index of the minimum value in a pixel, values above threshold are classified as background values.
auto minIndexLambda = [threshold, bv](PixelType const & pixel)
{
auto minElem = std::min_element(&pixel[0], &pixel[pixel.Size()]);
return static_cast<int>(*minElem < threshold ? std::distance(&pixel[0], minElem) + 1 : bv);
};
auto classificationFilter = NewFunctorFilter(minIndexLambda);
classificationFilter->SetInput(filterOutput);
SetParameterOutputImage("out", classificationFilter->GetOutput());
RegisterPipeline();
}
else
{
RegisterPipeline();
}
}
};
}
}
OTB_APPLICATION_EXPORT(otb::Wrapper::SpectralAngleClassification)
......@@ -75,3 +75,34 @@ otb_test_application(NAME apTvHyLocalRxDetection
${BASELINE}/apTvHyLocalRxDetection.tif
${TEMP}/apTvHyLocalRxDetection.tif
)
#----------- SpectralAngleClassification TESTS ------------------------
otb_test_application(NAME apTvHySpectralAngleClassification_sam
APP SpectralAngleClassification
OPTIONS -in ${OTB_DATA_ROOT}/Input/Hyperspectral/synthetic/hsi_cube.tif
-ie ${INPUTDATA}/Hyperspectral/synthetic/endmembers.tif
-out ${TEMP}/apTvHySpectralAngleClassification_sam.tif
-mode sam
VALID --compare-image ${EPSILON_9}
${BASELINE}/apTvHySpectralAngleClassification_sam.tif
${TEMP}/apTvHySpectralAngleClassification_sam.tif
)
#----------- LocalRxDetection TESTS ------------------------
otb_test_application(NAME apTvHySpectralAngleClassification_sid
APP SpectralAngleClassification
OPTIONS -in ${OTB_DATA_ROOT}/Input/Hyperspectral/synthetic/hsi_cube.tif
-ie ${INPUTDATA}/Hyperspectral/synthetic/endmembers.tif
-out ${TEMP}/apTvHySpectralAngleClassification_sid.tif
-mode sid
-measure ${TEMP}/apTvHySpectralAngleClassification_sid_measure.tif
-threshold 0.1
-bv -1
VALID --compare-n-images ${EPSILON_6} 2
${BASELINE}/apTvHySpectralAngleClassification_sid.tif
${TEMP}/apTvHySpectralAngleClassification_sid.tif
${BASELINE}/apTvHySpectralAngleClassification_sid_measure.tif
${TEMP}/apTvHySpectralAngleClassification_sid_measure.tif
)
......@@ -23,12 +23,16 @@
#include <algorithm>
#include "otbMath.h"
#include "otbSpectralAngleFunctor.h"
namespace otb
{
/** \class BinarySpectralAngleFunctor
* \brief This functor computes the spectral angle between two pixels.
*
* If the pixels have different sizes, only the first components of the
* largest pixel will be considered.
*
* It can be used as a functor in a BinaryFunctorImageFilter to
* compute the pixel-by-pixel spectral angles.
*
......@@ -40,40 +44,24 @@ template <class TInput1, class TInput2, class TOutputValue>
class BinarySpectralAngleFunctor
{
public:
BinarySpectralAngleFunctor()
{
}
BinarySpectralAngleFunctor() = default;
virtual ~BinarySpectralAngleFunctor()
{
}
virtual ~BinarySpectralAngleFunctor() = default;
// Binary operator
inline TOutputValue operator()(const TInput1& a, const TInput2& b) const
TOutputValue operator()(const TInput1& in1, const TInput2& in2) const
{
const double Epsilon = 1E-10;
double dist = 0.0;
double scalarProd = 0.0;
double norma = 0.0;
double normb = 0.0;
double sqrtNormProd = 0.0;
for (unsigned int i = 0; i < std::min(a.Size(), b.Size()); ++i)
// Compute norms.
auto in1Norm = 0;
auto in2Norm = 0;
auto nbIter = std::min(in1.Size(), in2.Size());
for (unsigned int i = 0; i < nbIter; ++i)
{
scalarProd += a[i] * b[i];
norma += a[i] * a[i];
normb += b[i] * b[i];
in1Norm += in1[i] * in1[i];
in2Norm += in2[i] * in2[i];
}
sqrtNormProd = std::sqrt(norma * normb);
if (std::abs(sqrtNormProd) < Epsilon || scalarProd / sqrtNormProd > 1)
{
dist = 0.0;
}
else
{
dist = std::acos(scalarProd / sqrtNormProd);
}
return static_cast<TOutputValue>(dist);
return SpectralAngleDetails::ComputeSpectralAngle<TInput1, TInput2, TOutputValue>(in1, in1Norm, in2, in2Norm);
}
};
......
......@@ -26,6 +26,7 @@
namespace otb
{
/** \class SpectralAngleDistanceImageFilter
* \deprecated
* \brief This filter implements the computation of the spectral angle
* distance with respect to a reference pixel.
*
......
......@@ -23,16 +23,45 @@
#include "otbMath.h"
#include <algorithm>
#include <vector>
#include <numeric>
namespace otb
{
namespace Functor
{
namespace SpectralAngleDetails
{
/** \fn
* \brief This function computes spectral angle between a pixel and a reference, the norm of these inputs.
* should also be given as parameter of this function.
* */
template <class TInput, class TReference, class TOutput>
TOutput ComputeSpectralAngle(TInput const & input, typename TInput ::ValueType const & inputNorm,
TReference const & reference, typename TReference::ValueType refNorm)
{
auto minSize = std::min(input.Size(), reference.Size());
double scalarProduct = std::inner_product(&input[0], &input[minSize], &reference[0],0. );
auto normProd = inputNorm * refNorm;
if ((normProd < 1.e-12) || (scalarProduct / normProd > 1))
{
return static_cast<TOutput>(0.0);
}
else
{
return static_cast<TOutput>(std::acos(scalarProduct / normProd));
}
}
} // end namespace SpectralAngleDetails
/** \class SpectralAngleFunctor
* \brief This functor computes the spectral angle according to a reference pixel.
*
* \ingroup OTBImageManipulation
*/
namespace Functor
{
template <class TInput, class TOutputValue>
class SpectralAngleFunctor
{
......@@ -43,64 +72,82 @@ public:
m_ReferencePixel.Fill(1);
}
virtual ~SpectralAngleFunctor()
{
}
~SpectralAngleFunctor() = default;
// Binary operator
inline TOutputValue operator()(const TInput& inPix) const
inline TOutputValue operator()(TInput const & inPix) const
{
return this->Evaluate(inPix);
return SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, TOutputValue>(inPix, inPix.GetNorm(), m_ReferencePixel, m_RefNorm);
}
void SetReferencePixel(TInput ref)
void SetReferencePixel(TInput const & ref)
{
m_ReferencePixel = ref;
m_RefNorm = 0.0;
for (unsigned int i = 0; i < ref.Size(); ++i)
{
m_RefNorm += ref[i] * ref[i];
}
m_RefNorm = std::sqrt(static_cast<double>(m_RefNorm));
m_RefNorm = ref.GetNorm();
}
TInput GetReferencePixel() const
{
return m_ReferencePixel;
}
protected:
// This method can be reimplemented in subclasses to actually
// compute the index value
virtual TOutputValue Evaluate(const TInput& inPix) const
private :
TInput m_ReferencePixel;
double m_RefNorm;
};
/** \class SpectralAngleMapperFunctor
* \brief This functor computes the spectral angle according to a vector of reference pixel.
*
* \ingroup OTBImageManipulation
*/
template <class TInput, class TReference, class TOutput>
class SpectralAngleMapperFunctor
{
public:
SpectralAngleMapperFunctor() = default;
virtual ~SpectralAngleMapperFunctor() = default;
// Binary operator
inline TOutput operator()(const TInput& inPix) const
{
TOutputValue out;
double dist = 0.0;
double scalarProd = 0.0;
double normProd = 0.0;
double normProd1 = 0.0;
double sqrtNormProd = 0.0;
for (unsigned int i = 0; i < std::min(inPix.Size(), m_ReferencePixel.Size()); ++i)
{
scalarProd += inPix[i] * m_ReferencePixel[i];
normProd1 += inPix[i] * inPix[i];
}
normProd = normProd1 * m_RefNorm * m_RefNorm;
sqrtNormProd = std::sqrt(normProd);
if ((sqrtNormProd == 0.0) || (scalarProd / sqrtNormProd > 1))
TOutput res(m_ReferencePixels.size());
auto inputNorm = inPix.GetNorm();
for (unsigned int i = 0; i< m_ReferencePixels.size(); i++)
{
dist = 0.0;
res[i] = SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, typename TOutput::ValueType>
(inPix, inputNorm, m_ReferencePixels[i], m_ReferenceNorm[i]);
}
else
return res;
}
size_t OutputSize(...) const
{
return m_ReferencePixels.size();
}
void SetReferencePixels(std::vector<TReference> ref)
{
m_ReferencePixels = std::move(ref);
m_ReferenceNorm.clear();
// Precompute the norm of reference pixels
for (auto const & pixel : m_ReferencePixels)
{
dist = std::acos(scalarProd / sqrtNormProd);
m_ReferenceNorm.push_back(pixel.GetNorm());
}
out = static_cast<TOutputValue>(dist);
return out;
}
std::vector<TReference> const & GetReferencePixels() const
{
return m_ReferencePixels;
}
TInput m_ReferencePixel;
double m_RefNorm;
private:
std::vector<TReference> m_ReferencePixels;
std::vector<double> m_ReferenceNorm;
};
} // end namespace functor
......
/*
* Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
*
* This file is part of Orfeo Toolbox
*