diff --git a/Data/Baseline/OTB/Images/apTvHySpectralAngleClassification_sam.tif b/Data/Baseline/OTB/Images/apTvHySpectralAngleClassification_sam.tif
new file mode 100644
index 0000000000000000000000000000000000000000..5b40c8cb58ca17168eb0b64bd9b258172603014f
--- /dev/null
+++ b/Data/Baseline/OTB/Images/apTvHySpectralAngleClassification_sam.tif
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:41ba89a71a9bb2af17ba1bc66f0e1197fb843b59a207a7eae6bffb44f2578dc1
+size 400
diff --git a/Data/Baseline/OTB/Images/apTvHySpectralAngleClassification_sid.tif b/Data/Baseline/OTB/Images/apTvHySpectralAngleClassification_sid.tif
new file mode 100644
index 0000000000000000000000000000000000000000..6e4a1386201eff63effa649d782c66ce280c8f18
--- /dev/null
+++ b/Data/Baseline/OTB/Images/apTvHySpectralAngleClassification_sid.tif
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e86975e7d76ce475074fef8fd6d276fd8bb3c5bc091b296402cbf146154dfb7f
+size 684
diff --git a/Data/Baseline/OTB/Images/apTvHySpectralAngleClassification_sid_measure.tif b/Data/Baseline/OTB/Images/apTvHySpectralAngleClassification_sid_measure.tif
new file mode 100644
index 0000000000000000000000000000000000000000..68c604fc5baf6e60fbd81d6bc57ca8b24d68e4ca
--- /dev/null
+++ b/Data/Baseline/OTB/Images/apTvHySpectralAngleClassification_sid_measure.tif
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2c6285dfdf032e26076dfbde20b23d99931a49ca392a72b7e0671792d3dfa691
+size 6332
diff --git a/Documentation/Cookbook/Art/HyperspectralImages/classification.png b/Documentation/Cookbook/Art/HyperspectralImages/classification.png
new file mode 100644
index 0000000000000000000000000000000000000000..c84e24bea344beca46782e3dee9b9af5a29948e9
Binary files /dev/null and b/Documentation/Cookbook/Art/HyperspectralImages/classification.png differ
diff --git a/Documentation/Cookbook/rst/recipes/hyperspectral.rst b/Documentation/Cookbook/rst/recipes/hyperspectral.rst
index 9e89886fdf2fe302a99b22bee85e1650a60d712f..ef6fd1b451c757f5bface174383cd542542b7844 100644
--- a/Documentation/Cookbook/rst/recipes/hyperspectral.rst
+++ b/Documentation/Cookbook/rst/recipes/hyperspectral.rst
@@ -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
+
diff --git a/Documentation/Cookbook/rst/recipes/pbclassif.rst b/Documentation/Cookbook/rst/recipes/pbclassif.rst
index dbeeadd8abb8b9b775db2ffb3f65b834fbdea7b0..8a13d970b9ec03ff6d301b350060e1db8f166ab6 100644
--- a/Documentation/Cookbook/rst/recipes/pbclassif.rst
+++ b/Documentation/Cookbook/rst/recipes/pbclassif.rst
@@ -1,3 +1,5 @@
+.. _classif:
+
Classification
==============
diff --git a/Modules/Applications/AppHyperspectral/app/CMakeLists.txt b/Modules/Applications/AppHyperspectral/app/CMakeLists.txt
index 062aecb339119c6aa2a0154f23b45900a75c0503..9ab5697e1ec3c38314903f24b0debdafb1bf9574 100644
--- a/Modules/Applications/AppHyperspectral/app/CMakeLists.txt
+++ b/Modules/Applications/AppHyperspectral/app/CMakeLists.txt
@@ -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})
+
diff --git a/Modules/Applications/AppHyperspectral/app/otbSpectralAngleClassification.cxx b/Modules/Applications/AppHyperspectral/app/otbSpectralAngleClassification.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0755543feb9be40718321500af4b61bc2df8f5f0
--- /dev/null
+++ b/Modules/Applications/AppHyperspectral/app/otbSpectralAngleClassification.cxx
@@ -0,0 +1,217 @@
+/*
+ * 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)
diff --git a/Modules/Applications/AppHyperspectral/test/CMakeLists.txt b/Modules/Applications/AppHyperspectral/test/CMakeLists.txt
index 1a9b6ae1e8eaaff28027a0b93ca92c8b7057dbd7..7165de5422174e7f9d819cedafa11be0280c6658 100644
--- a/Modules/Applications/AppHyperspectral/test/CMakeLists.txt
+++ b/Modules/Applications/AppHyperspectral/test/CMakeLists.txt
@@ -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
+)
diff --git a/Modules/Filtering/ImageManipulation/include/otbBinarySpectralAngleFunctor.h b/Modules/Filtering/ImageManipulation/include/otbBinarySpectralAngleFunctor.h
index cf87c9b597712ceae6cd600920812a03dbd9a6c6..cc483ab61aad7b6cad8b711e828072afe1ebd660 100644
--- a/Modules/Filtering/ImageManipulation/include/otbBinarySpectralAngleFunctor.h
+++ b/Modules/Filtering/ImageManipulation/include/otbBinarySpectralAngleFunctor.h
@@ -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);
}
};
diff --git a/Modules/Filtering/ImageManipulation/include/otbSpectralAngleDistanceImageFilter.h b/Modules/Filtering/ImageManipulation/include/otbSpectralAngleDistanceImageFilter.h
index 5c3d09538ccc9eaf5b28293f171800e078113409..e22e7192d9cf36e2638d97e3def9810ff9045773 100644
--- a/Modules/Filtering/ImageManipulation/include/otbSpectralAngleDistanceImageFilter.h
+++ b/Modules/Filtering/ImageManipulation/include/otbSpectralAngleDistanceImageFilter.h
@@ -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.
*
diff --git a/Modules/Filtering/ImageManipulation/include/otbSpectralAngleFunctor.h b/Modules/Filtering/ImageManipulation/include/otbSpectralAngleFunctor.h
index 1ca9a87501aaa7f646b638b5b5384c92ecdf19af..6ad8f3a6733a522bc49eabfc0976edfafc0afd1b 100644
--- a/Modules/Filtering/ImageManipulation/include/otbSpectralAngleFunctor.h
+++ b/Modules/Filtering/ImageManipulation/include/otbSpectralAngleFunctor.h
@@ -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
diff --git a/Modules/Filtering/ImageManipulation/include/otbSpectralInformationDivergenceFunctor.h b/Modules/Filtering/ImageManipulation/include/otbSpectralInformationDivergenceFunctor.h
new file mode 100644
index 0000000000000000000000000000000000000000..5e52e0b4a23d9b0e3d61bc633aaf78606c0901ab
--- /dev/null
+++ b/Modules/Filtering/ImageManipulation/include/otbSpectralInformationDivergenceFunctor.h
@@ -0,0 +1,118 @@
+/*
+ * 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.
+ */
+
+#ifndef otbSpectralInformationDivergenceFunctor_h
+#define otbSpectralInformationDivergenceFunctor_h
+
+#include "otbMath.h"
+#include <algorithm>
+#include <vector>
+#include <numeric>
+
+namespace otb
+{
+namespace Functor
+{
+
+/** \class SpectralAngleFunctor
+ * \brief This functor computes the spectral information divergence according to a reference pixel.
+ *
+ * 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.
+ *
+ * \ingroup OTBImageManipulation
+ */
+template <class TInput, class TReference, class TOutput>
+class SpectralInformationDivergenceFunctor
+{
+public:
+ SpectralInformationDivergenceFunctor() = default;
+ virtual ~SpectralInformationDivergenceFunctor() = default;
+
+ using OutputValueType = typename TOutput::ValueType;
+
+ // Binary operator
+ inline TOutput operator()(const TInput& input) const
+ {
+ TOutput res;
+ res.SetSize(m_ReferenceProbabilities.size());
+
+ auto inputProbability = ComputeProbabilityMassFunction(input);
+
+ for (unsigned int i = 0; i< m_ReferenceProbabilities.size(); i++)
+ {
+ res[i] = ComputeSpectralInformationDivergence(inputProbability, m_ReferenceProbabilities[i]);
+ }
+
+ return res;
+ }
+
+ size_t OutputSize(...) const
+ {
+ return m_ReferenceProbabilities.size();
+ }
+
+ void SetReferencePixels(std::vector<TReference> const & ref)
+ {
+ // We only store the probability mass function associated with the reference pixels, are the latter are not needed
+ // in the sid computation.
+ m_ReferenceProbabilities.clear();
+ for (auto const & pixel : ref)
+ {
+ m_ReferenceProbabilities.push_back(ComputeProbabilityMassFunction(pixel));
+ }
+ }
+
+private:
+ inline TInput ComputeProbabilityMassFunction(TInput const & input) const
+ {
+ for (unsigned int i = 0; i < input.Size(); i++)
+ {
+ // Input pixel should be non negative (e.g. reflectance, radiance)
+ if (input[i] <= 0)
+ {
+ throw std::runtime_error("Input pixels of the spectral information divergence algorithm should be strictly positive.");
+ }
+ }
+
+ return input / std::accumulate(&input[0], &input[input.Size()], 0.0);
+ }
+
+ inline OutputValueType ComputeSpectralInformationDivergence(TInput const & p, TInput const & q) const
+ {
+ assert(p.Size() == q.Size());
+ OutputValueType sid = 0.0;
+ for (unsigned int i = 0; i < p.Size(); i++)
+ {
+ // Compute SID : p[i] * std::log(p[i]/q[i]) + q[i] * std::log(q[i]/p[i]);
+ sid += (p[i] - q[i]) * std::log(p[i]/q[i]);
+ }
+ return sid;
+ }
+
+ /** Probability mass function associated with the reference pixel */
+ std::vector<TReference> m_ReferenceProbabilities;
+};
+
+} // end namespace functor
+} // end namespace otb
+
+#endif //otbSpectralInformationDivergenceFunctor_h
diff --git a/Modules/Filtering/ImageManipulation/include/otbSqrtSpectralAngleFunctor.h b/Modules/Filtering/ImageManipulation/include/otbSqrtSpectralAngleFunctor.h
index e5b2350c335f4a39be909e96b964469faac87d71..b708ba2ee32d7af937de34785ba8b779f720e047 100644
--- a/Modules/Filtering/ImageManipulation/include/otbSqrtSpectralAngleFunctor.h
+++ b/Modules/Filtering/ImageManipulation/include/otbSqrtSpectralAngleFunctor.h
@@ -35,25 +35,44 @@ namespace Functor
*
* \ingroup OTBImageManipulation
*/
-template <class TInputVectorPixel, class TOutputPixel>
-class SqrtSpectralAngleFunctor : public SpectralAngleFunctor<TInputVectorPixel, TOutputPixel>
+
+/** \class SpectralAngleFunctor
+ * \brief This functor computes the spectral angle according to a reference pixel.
+ *
+ * \ingroup OTBImageManipulation
+ */
+template <class TInput, class TOutputValue>
+class SqrtSpectralAngleFunctor
{
public:
- typedef SqrtSpectralAngleFunctor Self;
- typedef SpectralAngleFunctor<TInputVectorPixel, TOutputPixel> Superclass;
-
SqrtSpectralAngleFunctor()
{
+ m_ReferencePixel.SetSize(4);
+ m_ReferencePixel.Fill(1);
}
- ~SqrtSpectralAngleFunctor() override
+
+ virtual ~SqrtSpectralAngleFunctor() = default;
+
+ // Binary operator
+ inline TOutputValue operator()(TInput const & inPix) const
{
+ return std::sqrt(SpectralAngleDetails::ComputeSpectralAngle<TInput, TInput, TOutputValue>(inPix, inPix.GetNorm(), m_ReferencePixel, m_RefNorm));
}
-protected:
- TOutputPixel Evaluate(const TInputVectorPixel& inPix) const override
+ void SetReferencePixel(TInput const & ref)
{
- return static_cast<TOutputPixel>(std::sqrt(Superclass::Evaluate(inPix)));
+ m_ReferencePixel = ref;
+ m_RefNorm = ref.GetNorm();
}
+
+ TInput GetReferencePixel() const
+ {
+ return m_ReferencePixel;
+ }
+
+private :
+ TInput m_ReferencePixel;
+ double m_RefNorm;
};
} // end namespace Functor
diff --git a/Modules/Radiometry/Indices/include/otbWaterSqrtSpectralAngleImageFilter.h b/Modules/Radiometry/Indices/include/otbWaterSqrtSpectralAngleImageFilter.h
index ed9e8e56f4f767962a6ea4fff8358702467a5a00..57b9c94febaf66dc08d9726be7e91fdce8161970 100644
--- a/Modules/Radiometry/Indices/include/otbWaterSqrtSpectralAngleImageFilter.h
+++ b/Modules/Radiometry/Indices/include/otbWaterSqrtSpectralAngleImageFilter.h
@@ -23,6 +23,7 @@
#include "otbSqrtSpectralAngleFunctor.h"
#include "itkUnaryFunctorImageFilter.h"
+#include "otbRadiometricIndex.h"
namespace otb
{
@@ -38,100 +39,102 @@ namespace Functor
*
* \ingroup OTBIndices
*/
-template <class TInputVectorPixel, class TOutputPixel>
-class WaterSqrtSpectralAngleFunctor : public SqrtSpectralAngleFunctor<TInputVectorPixel, TOutputPixel>
+template <class TInput, class TOutput>
+class WaterSqrtSpectralAngleFunctor : public RadiometricIndex<TInput, TOutput>
{
public:
- typedef WaterSqrtSpectralAngleFunctor Self;
- typedef SqrtSpectralAngleFunctor<TInputVectorPixel, TOutputPixel> Superclass;
- typedef TInputVectorPixel InputVectorPixelType;
- WaterSqrtSpectralAngleFunctor()
- {
+ using typename RadiometricIndex<TInput, TOutput>::PixelType;
- // Set the channels indices
- m_BlueIndex = 0;
- m_GreenIndex = 1;
- m_RedIndex = 2;
- m_NIRIndex = 3;
-
- // Set reference water value
- InputVectorPixelType reference;
- reference.SetSize(4);
- reference[0] = 136.0;
- reference[1] = 132.0;
- reference[2] = 47.0;
- reference[3] = 24.0;
- this->SetReferenceWaterPixel(reference);
+ WaterSqrtSpectralAngleFunctor()
+ : RadiometricIndex<TInput, TOutput>({CommonBandNames::BLUE, CommonBandNames::GREEN, CommonBandNames::RED, CommonBandNames::NIR} )
+ {
+ // Set default reference water value
+ m_ReferencePixel.SetSize(4);
+ m_ReferencePixel[0] = 136.0;
+ m_ReferencePixel[1] = 132.0;
+ m_ReferencePixel[2] = 47.0;
+ m_ReferencePixel[3] = 24.0;
+ m_RefNorm = m_ReferencePixel.GetNorm();
}
- ~WaterSqrtSpectralAngleFunctor() override
+
+ virtual ~WaterSqrtSpectralAngleFunctor() = default;
+
+ /** Set Reference Pixel
+ * \param ref : The new reference pixel, the band indices will be used.
+ * */
+ void SetReferenceWaterPixel(PixelType ref)
{
+ assert(m_ReferencePixel.Size() == 4);
+ m_ReferencePixel[0] = this->Value(CommonBandNames::BLUE, ref);
+ m_ReferencePixel[1] = this->Value(CommonBandNames::GREEN, ref);
+ m_ReferencePixel[2] = this->Value(CommonBandNames::RED, ref);
+ m_ReferencePixel[3] = this->Value(CommonBandNames::NIR, ref);
+ m_RefNorm = m_ReferencePixel.GetNorm();
}
- /** Set Reference Pixel */
- void SetReferenceWaterPixel(InputVectorPixelType ref)
+ // Binary operator
+ inline TOutput operator()(PixelType const & inPix) const override
{
- if (ref.GetSize() != 4)
- {
- }
- InputVectorPixelType reference;
- reference.SetSize(4);
- reference[m_BlueIndex] = ref[0];
- reference[m_GreenIndex] = ref[1];
- reference[m_RedIndex] = ref[2];
- reference[m_NIRIndex] = ref[3];
- this->SetReferencePixel(reference);
+ PixelType pix(4);
+ pix[0] = this->Value(CommonBandNames::BLUE, inPix);
+ pix[1] = this->Value(CommonBandNames::GREEN, inPix);
+ pix[2] = this->Value(CommonBandNames::RED, inPix);
+ pix[3] = this->Value(CommonBandNames::NIR, inPix);
+
+ return std::sqrt(SpectralAngleDetails::ComputeSpectralAngle<PixelType, PixelType, TOutput>
+ (pix, pix.GetNorm(),
+ m_ReferencePixel, m_RefNorm));
}
- /** Getters and setters */
+ /**@{*/
+ /** Legacy getters and setters (use SetBandIndex/GetBandIndex instead)
+ * \deprecated
+ * */
void SetBlueChannel(unsigned int channel)
{
- m_BlueIndex = channel;
+ this->SetBandIndex(CommonBandNames::BLUE, channel + 1);
}
unsigned int GetBlueChannel() const
{
- return m_BlueIndex;
+ return this->GetBandIndex(CommonBandNames::BLUE) - 1;
}
void SetGreenChannel(unsigned int channel)
{
- m_GreenIndex = channel;
+ this->SetBandIndex(CommonBandNames::GREEN, channel + 1);
}
unsigned int GetGreenChannel() const
{
- return m_GreenIndex;
+ return this->GetBandIndex(CommonBandNames::GREEN) - 1;
}
void SetRedChannel(unsigned int channel)
{
- m_RedIndex = channel;
+ this->SetBandIndex(CommonBandNames::RED, channel + 1);
}
unsigned int GetRedChannel() const
{
- return m_RedIndex;
+ return this->GetBandIndex(CommonBandNames::RED) - 1;
}
void SetNIRChannel(unsigned int channel)
{
- m_NIRIndex = channel;
+ this->SetBandIndex(CommonBandNames::NIR, channel + 1);
}
unsigned int GetNIRChannel() const
{
- return m_NIRIndex;
+ return this->GetBandIndex(CommonBandNames::NIR) - 1;
}
+ /**@}*/
+
-protected:
- inline TOutputPixel Evaluate(const TInputVectorPixel& inPix) const override
- {
- return static_cast<TOutputPixel>(Superclass::Evaluate(inPix));
- }
+private:
- /** Channels */
- int m_BlueIndex;
- int m_GreenIndex;
- int m_RedIndex;
- int m_NIRIndex;
+ PixelType m_ReferencePixel;
+ double m_RefNorm;
};
} // End namespace Functor
/** \class WaterSqrtSpectralAngleImageFilter
+ * \deprecated
* \brief Compute a radiometric water indice
*
* This filter calculates a pixel wise water indice by calculating
@@ -147,11 +150,11 @@ protected:
* probability of water.
*
* \sa SpectralAngleDistanceImageFilter
- *
+ *
* \ingroup OTBIndices
*/
template <class TInputVectorImage, class TOutputImage,
- class TFunction = Functor::WaterSqrtSpectralAngleFunctor<typename TInputVectorImage::PixelType, typename TOutputImage::PixelType>>
+ class TFunction = Functor::WaterSqrtSpectralAngleFunctor<double, typename TOutputImage::PixelType>>
class ITK_EXPORT WaterSqrtSpectralAngleImageFilter : public itk::UnaryFunctorImageFilter<TInputVectorImage, TOutputImage, TFunction>
{
public:
diff --git a/Modules/Radiometry/Indices/test/otbWaterSqrtSpectralAngleImageFilter.cxx b/Modules/Radiometry/Indices/test/otbWaterSqrtSpectralAngleImageFilter.cxx
index 589a168fd9c12e9b6ac7826e670384efc9f87111..83f81eebade31d3016f5e4bbb85f34e3be27c1b1 100644
--- a/Modules/Radiometry/Indices/test/otbWaterSqrtSpectralAngleImageFilter.cxx
+++ b/Modules/Radiometry/Indices/test/otbWaterSqrtSpectralAngleImageFilter.cxx
@@ -39,6 +39,7 @@ int otbWaterSqrtSpectralAngleImageFilter(int itkNotUsed(argc), char* argv[])
// Instantiating objects
WaterSqrtSpectralAngleImageFilterType::Pointer filter = WaterSqrtSpectralAngleImageFilterType::New();
+
ReaderType::Pointer reader = ReaderType::New();
WriterType::Pointer writer = WriterType::New();
@@ -60,7 +61,7 @@ int otbWaterSqrtSpectralAngleImageFilter(int itkNotUsed(argc), char* argv[])
filter->SetInput(reader->GetOutput());
writer->SetInput(filter->GetOutput());
+
writer->Update();
-
return EXIT_SUCCESS;
}
diff --git a/Modules/Segmentation/CCOBIA/include/otbConnectedComponentMuParserFunctor.h b/Modules/Segmentation/CCOBIA/include/otbConnectedComponentMuParserFunctor.h
index 26e6b2034bfbda61e6e1235f0d05a6d489142d1b..f9e3ed36cb2e022746c436403bc136e4e216ff72 100644
--- a/Modules/Segmentation/CCOBIA/include/otbConnectedComponentMuParserFunctor.h
+++ b/Modules/Segmentation/CCOBIA/include/otbConnectedComponentMuParserFunctor.h
@@ -25,7 +25,7 @@
#include "otbParser.h"
#include "otbMacro.h"
-
+#include "otbBinarySpectralAngleFunctor.h"
#include <vnl/algo/vnl_lsqr.h>
#include <vnl/vnl_sparse_matrix_linear_system.h>
@@ -115,31 +115,9 @@ public:
m_IntensityP2 = m_IntensityP2 / (static_cast<double>(m_NbOfBands));
m_Distance = std::sqrt(m_Distance);
-
- // compute spectralAngle
- double scalarProd = 0.0;
- double normProd = 0.0;
- double normProd1 = 0.0;
- double normProd2 = 0.0;
-
- for (unsigned int i = 0; i < p1.Size(); ++i)
- {
- scalarProd += p1[i] * p2[i];
- normProd1 += p1[i] * p1[i];
- normProd2 += p2[i] * p2[i];
- }
- normProd = normProd1 * normProd2;
-
- if (normProd == 0.0)
- {
- m_SpectralAngle = 0.0;
- }
- else
- {
- m_SpectralAngle = std::acos(scalarProd / std::sqrt(normProd));
- }
-
- //
+
+ BinarySpectralAngleFunctor<TInput, TInput, double> spectralAngleFunctor;
+ m_SpectralAngle = spectralAngleFunctor(p1, p2);
value = m_Parser->Eval();