otbLandsatTMIndices.h

Go to the documentation of this file.

/*
 * Copyright (C) 2005-2022 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 otbLandsatTMIndices_h
#define otbLandsatTMIndices_h
 
#include "otbMath.h"
#include "otbBandName.h"
#include "itkVariableLengthVector.h"
#include "otbFuzzyVariable.h"
#include <vector>
#include <algorithm>
#include <string>
 
namespace otb
{
 
namespace Functor
{
 
namespace LandsatTM
{
 
enum SATType
{
  L5,
  L7
};
enum DegreeType
{
  Kelvin,
  HundredsKelvin,
  Celsius
};
enum ReflectanceType
{
  Thousands,
  Normalized
};
template <class TInput, class TOutput>
class LandsatTMIndexBase
{
public:
  // operators !=
  bool operator!=(const LandsatTMIndexBase&) const
  {
    return false;
  }
  // operator ==
  bool operator==(const LandsatTMIndexBase& other) const
  {
    return !(*this != other);
  }
 
  LandsatTMIndexBase()
    : m_EpsilonToBeConsideredAsZero(0.0000001),
      m_TM1(0),
      m_TM2(1),
      m_TM3(2),
      m_TM4(3),
      m_TM5(4),
      m_TM60(5),
      m_TM61(5),
      m_TM62(6),
      m_TM7(7),
      m_SAT(L7),
      m_Degree(Celsius),
      m_Reflectance(Normalized)
  {
  }
  virtual ~LandsatTMIndexBase()
  {
  }
 
  void SetIndex(BandName::LandsatTMBandNames band, unsigned int channel)
  {
    switch (band)
    {
    case BandName::TM1:
      m_TM1 = channel;
      break;
    case BandName::TM2:
      m_TM2 = channel;
      break;
    case BandName::TM3:
      m_TM3 = channel;
      break;
    case BandName::TM4:
      m_TM4 = channel;
      break;
    case BandName::TM5:
      m_TM5 = channel;
      break;
    case BandName::TM60:
      m_TM60 = channel;
      break;
    case BandName::TM61:
      m_TM61 = channel;
      break;
    case BandName::TM62:
      m_TM62 = channel;
      break;
    case BandName::TM7:
      m_TM7 = channel;
      break;
    }
  }
 
  unsigned int GetIndex(BandName::LandsatTMBandNames band) const
  {
    switch (band)
    {
    case BandName::TM1:
      return m_TM1;
      break;
    case BandName::TM2:
      return m_TM2;
      break;
    case BandName::TM3:
      return m_TM3;
      break;
    case BandName::TM4:
      return m_TM4;
      break;
    case BandName::TM5:
      return m_TM5;
      break;
    case BandName::TM60:
      return m_TM60;
      break;
    case BandName::TM61:
      return m_TM61;
      break;
    case BandName::TM62:
      return m_TM62;
      break;
    case BandName::TM7:
      return m_TM7;
      break;
    }
    return m_TM1;
  }
 
  unsigned int GetTM1() const
  {
    return this->m_TM1;
  }
 
  unsigned int GetTM2() const
  {
    return this->m_TM2;
  }
 
  unsigned int GetTM3() const
  {
    return this->m_TM3;
  }
 
  unsigned int GetTM4() const
  {
    return this->m_TM4;
  }
 
  unsigned int GetTM5() const
  {
    return this->m_TM5;
  }
 
  unsigned int GetTM60() const
  {
    return this->m_TM60;
  }
 
  unsigned int GetTM61() const
  {
    return this->m_TM61;
  }
 
  unsigned int GetTM62() const
  {
    return this->m_TM62;
  }
 
  unsigned int GetTM7() const
  {
    return this->m_TM7;
  }
 
  void SetSAT(SATType sat)
  {
 
    this->m_SAT = sat;
 
    if (sat == L5)
    {
      m_TM1  = 0;
      m_TM2  = 1;
      m_TM3  = 2;
      m_TM4  = 3;
      m_TM5  = 4;
      m_TM60 = 5;
      m_TM7  = 6;
      m_SAT  = L5;
    }
    else
    {
      m_TM1  = 0;
      m_TM2  = 1;
      m_TM3  = 2;
      m_TM4  = 3;
      m_TM5  = 4;
      m_TM61 = 5;
      m_TM62 = 6;
      m_TM7  = 7;
      m_SAT  = L7;
    }
  }
 
  SATType GetSAT() const
  {
    return this->m_SAT;
  }
 
  void SetDegree(DegreeType deg)
  {
    this->m_Degree = deg;
  }
 
  DegreeType GetDegre() const
  {
    return this->m_Degree;
  }
 
  void SetReflectance(ReflectanceType ref)
  {
    this->m_Reflectance = ref;
  }
 
  ReflectanceType GetReflectance() const
  {
    return this->m_Reflectance;
  }
 
  double GetEpsilonToBeConsideredAsZero() const
  {
    return this->m_EpsilonToBeConsideredAsZero;
  }
 
protected:
 
  TInput PrepareValues(const TInput& inputPixel)
  {
 
    TInput newPix(inputPixel);
    if (this->m_Degree == Kelvin)
    {
      if (this->m_SAT == L5)
      {
        newPix[this->m_TM60] = newPix[this->m_TM60] - 273.15;
      }
      else
      {
        newPix[this->m_TM61] = newPix[this->m_TM61] - 273.15;
        newPix[this->m_TM62] = newPix[this->m_TM62] - 273.15;
      }
    }
    else if (this->m_Degree == HundredsKelvin)
    {
      if (this->m_SAT == L5)
      {
        newPix[this->m_TM60] = newPix[this->m_TM60] / 100. - 273.15;
      }
      else
      {
        newPix[this->m_TM61] = newPix[this->m_TM61] / 100. - 273.15;
        newPix[this->m_TM62] = newPix[this->m_TM62] / 100. - 273.15;
      }
    }
    if (this->m_Reflectance == Thousands)
    {
      newPix[this->m_TM1] = newPix[this->m_TM1] / 1000.;
      newPix[this->m_TM2] = newPix[this->m_TM2] / 1000.;
      newPix[this->m_TM3] = newPix[this->m_TM3] / 1000.;
      newPix[this->m_TM4] = newPix[this->m_TM4] / 1000.;
      newPix[this->m_TM5] = newPix[this->m_TM5] / 1000.;
      newPix[this->m_TM7] = newPix[this->m_TM7] / 1000.;
    }
 
    return newPix;
  }
 
 
  double m_EpsilonToBeConsideredAsZero;
 
  unsigned int m_TM1;
  unsigned int m_TM2;
  unsigned int m_TM3;
  unsigned int m_TM4;
  unsigned int m_TM5;
  unsigned int m_TM60;
  unsigned int m_TM61;
  unsigned int m_TM62;
  unsigned int m_TM7;
 
  SATType         m_SAT;
  DegreeType      m_Degree;
  ReflectanceType m_Reflectance;
};
 
 
template <class TInput, class TOutput>
class LandsatTMIndex : public LandsatTMIndexBase<TInput, TOutput>
{
public:
  // Operator performing the work
  inline TOutput operator()(const TInput& inputPixel) const;
 
  virtual std::string GetName() const = 0;
 
  LandsatTMIndex()
  {
  }
  ~LandsatTMIndex() override
  {
  }
};
 
template <class TInput, class TOutput>
class Bright : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "Bright";
  }
 
  Bright()
  {
  }
  ~Bright() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
 
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = (newPixel[this->m_TM1] + newPixel[this->m_TM2] + 2 * newPixel[this->m_TM3] + 2 * newPixel[this->m_TM4] + newPixel[this->m_TM5] +
                     newPixel[this->m_TM7]) /
                    8.0;
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class Vis : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "Vis";
  }
 
  Vis()
  {
  }
  ~Vis() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = (newPixel[this->m_TM1] + newPixel[this->m_TM2] + newPixel[this->m_TM3]) / 3.0;
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class NIR : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "NIR";
  }
 
  NIR()
  {
  }
  ~NIR() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = newPixel[this->m_TM4];
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class MIR1 : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "MIR1";
  }
 
  MIR1()
  {
  }
  ~MIR1() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = newPixel[this->m_TM5];
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class MIR2 : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "MIR2";
  }
 
  MIR2()
  {
  }
  ~MIR2() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = newPixel[this->m_TM7];
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class TIR : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "TIR";
  }
 
  TIR()
  {
  }
  ~TIR() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = newPixel[this->m_TM62];
 
    if (this->m_SAT == L5)
      result = newPixel[this->m_TM60];
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class MIRTIR : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "MIRTIR";
  }
 
  MIRTIR()
  {
  }
  ~MIRTIR() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double tir  = newPixel[this->m_TM62];
    double mir1 = newPixel[this->m_TM5];
 
    if (this->m_SAT == L5)
      tir = newPixel[this->m_TM60];
 
    double result = 255 * (1 - mir1) * (tir + 100) / 100.;
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class NDVI : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "NDVI";
  }
 
  NDVI()
  {
  }
  ~NDVI() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = (newPixel[this->m_TM4] - newPixel[this->m_TM3]) / (newPixel[this->m_TM4] + newPixel[this->m_TM3] + this->m_EpsilonToBeConsideredAsZero);
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class NDBSI : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "NDBSI";
  }
 
  NDBSI()
  {
  }
  ~NDBSI() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = (newPixel[this->m_TM5] - newPixel[this->m_TM4]) / (newPixel[this->m_TM5] + newPixel[this->m_TM4] + this->m_EpsilonToBeConsideredAsZero);
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class BIO : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "BIO";
  }
 
  BIO()
  {
  }
  ~BIO() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = ((newPixel[this->m_TM5] + newPixel[this->m_TM3]) - (newPixel[this->m_TM4] + newPixel[this->m_TM1])) /
                    ((newPixel[this->m_TM5] + newPixel[this->m_TM3]) + (newPixel[this->m_TM4] + newPixel[this->m_TM1]));
 
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class NDSI : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "NDSI";
  }
 
  NDSI()
  {
  }
  ~NDSI() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = (newPixel[this->m_TM2] - newPixel[this->m_TM5]) / (newPixel[this->m_TM2] + newPixel[this->m_TM5] + this->m_EpsilonToBeConsideredAsZero);
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class NDSIVis : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "NDSIVis";
  }
 
  NDSIVis()
  {
  }
  ~NDSIVis() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double vis    = (newPixel[this->m_TM1] + newPixel[this->m_TM2] + newPixel[this->m_TM3]) / 3.0;
    double result = (vis - newPixel[this->m_TM5]) / (vis + newPixel[this->m_TM5] + this->m_EpsilonToBeConsideredAsZero);
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class NDBBBI : public LandsatTMIndex<TInput, TOutput>
{
public:
  std::string GetName() const override
  {
    return "NDBBBI";
  }
 
  NDBBBI()
  {
  }
  ~NDBBBI() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    double result = (newPixel[this->m_TM1] - newPixel[this->m_TM5]) / (newPixel[this->m_TM1] + newPixel[this->m_TM5] + this->m_EpsilonToBeConsideredAsZero);
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput>
class LinguisticVariables : public LandsatTMIndexBase<TInput, itk::FixedArray<unsigned int, 11>>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef typename itk::FixedArray<unsigned int, 11>        OutputPixelType;
  typedef otb::FuzzyVariable<unsigned short, PrecisionType> FuzzyVarType;
 
  enum LinguisticValues
  {
    MINlv = 0,
    Low   = MINlv,
    Medium,
    MAXlv = 2,
    High  = MAXlv
  };
  enum Indices
  {
    MINid  = 0,
    bright = MINid,
    vis,
    nir,
    mir1,
    mir2,
    tir,
    mirtir,
    ndsivis,
    ndbbbi,
    ndvi,
    MAXid = 10,
    ndbsi = MAXid
  };
 
  virtual std::string GetName() const
  {
    return "LandsatTM Linguistic Variables";
  }
 
  LinguisticVariables()
  {
    m_FvBright  = FuzzyVarType::New();
    m_FvVis     = FuzzyVarType::New();
    m_FvNIR     = FuzzyVarType::New();
    m_FvMIR1    = FuzzyVarType::New();
    m_FvMIR2    = FuzzyVarType::New();
    m_FvTIR     = FuzzyVarType::New();
    m_FvMIRTIR  = FuzzyVarType::New();
    m_FvNDSIVis = FuzzyVarType::New();
    m_FvNDBBBI  = FuzzyVarType::New();
    m_FvNDVI    = FuzzyVarType::New();
    m_FvNDBSI   = FuzzyVarType::New();
 
    PrecisionType maximumValue = itk::NumericTraits<PrecisionType>::max();
    PrecisionType minimumValue = itk::NumericTraits<PrecisionType>::NonpositiveMin();
 
    // the thresholds are computed wrt Baraldi's paper (normalized 0-255 values)
    m_FvBright->SetMembership(Low, minimumValue, minimumValue, 40 / 255., 40 / 255.);
    m_FvBright->SetMembership(Medium, 40 / 255., 40 / 255., 60 / 255., 60 / 255.);
    m_FvBright->SetMembership(High, 60 / 255., 60 / 255., maximumValue, maximumValue);
 
    m_FvVis->SetMembership(Low, minimumValue, minimumValue, 30 / 255., 30 / 255.);
    m_FvVis->SetMembership(Medium, 30 / 255., 30 / 255., 50 / 255., 50 / 255.);
    m_FvVis->SetMembership(High, 50 / 255., 50 / 255., maximumValue, maximumValue);
 
    m_FvNIR->SetMembership(Low, minimumValue, minimumValue, 40 / 255., 40 / 255.);
    m_FvNIR->SetMembership(Medium, 40 / 255., 40 / 255., 60 / 255., 60 / 255.);
    m_FvNIR->SetMembership(High, 60 / 255., 60 / 255., maximumValue, maximumValue);
 
    m_FvMIR1->SetMembership(Low, minimumValue, minimumValue, 40 / 255., 40 / 255.);
    m_FvMIR1->SetMembership(Medium, 40 / 255., 40 / 255., 60 / 255., 60 / 255.);
    m_FvMIR1->SetMembership(High, 60 / 255., 60 / 255., maximumValue, maximumValue);
 
    m_FvMIR2->SetMembership(Low, minimumValue, minimumValue, 30 / 255., 30 / 255.);
    m_FvMIR2->SetMembership(Medium, 30 / 255., 30 / 255., 50 / 255., 50 / 255.);
    m_FvMIR2->SetMembership(High, 50 / 255., 50 / 255., maximumValue, maximumValue);
 
    m_FvTIR->SetMembership(Low, minimumValue, minimumValue, 0, 0);
    m_FvTIR->SetMembership(Medium, 0, 0, 28, 28);
    m_FvTIR->SetMembership(High, 28, 28, maximumValue, maximumValue);
 
    m_FvMIRTIR->SetMembership(Low, minimumValue, minimumValue, 180, 180);
    m_FvMIRTIR->SetMembership(Medium, 180, 180, 220, 220);
    m_FvMIRTIR->SetMembership(High, 220, 220, maximumValue, maximumValue);
 
    m_FvNDSIVis->SetMembership(Low, minimumValue, minimumValue, 0, 0);
    m_FvNDSIVis->SetMembership(Medium, 0, 0, 0.5, 0.5);
    m_FvNDSIVis->SetMembership(High, 0.5, 0.5, maximumValue, maximumValue);
 
    m_FvNDBBBI->SetMembership(Low, minimumValue, minimumValue, -0.20, -0.20);
    m_FvNDBBBI->SetMembership(Medium, -0.20, -0.20, 0.10, 0.10);
    m_FvNDBBBI->SetMembership(High, 0.10, 0.10, maximumValue, maximumValue);
 
    m_FvNDVI->SetMembership(Low, minimumValue, minimumValue, 0.35, 0.35);
    m_FvNDVI->SetMembership(Medium, 0.35, 0.35, 0.6, 0.6);
    m_FvNDVI->SetMembership(High, 0.6, 0.6, maximumValue, maximumValue);
 
    m_FvNDBSI->SetMembership(Low, minimumValue, minimumValue, -0.20, -0.20);
    m_FvNDBSI->SetMembership(Medium, -0.20, -0.20, 0.0, 0.0);
    m_FvNDBSI->SetMembership(High, 0.0, 0.0, maximumValue, maximumValue);
  }
  ~LinguisticVariables() override
  {
  }
 
  inline itk::FixedArray<unsigned int, 11> operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    itk::FixedArray<unsigned int, 11> result;
 
    result[bright] = m_FvBright->GetMaxVar(Bright<TInput, PrecisionType>()(newPixel));
    result[vis]    = m_FvVis->GetMaxVar(Vis<TInput, PrecisionType>()(newPixel));
    result[nir]    = m_FvNIR->GetMaxVar(NIR<TInput, PrecisionType>()(newPixel));
    result[mir1]   = m_FvMIR1->GetMaxVar(MIR1<TInput, PrecisionType>()(newPixel));
 
    MIR2<TInput, PrecisionType> mir2F;
    mir2F.SetSAT(this->m_SAT);
    result[mir2] = m_FvMIR2->GetMaxVar(mir2F(newPixel));
 
    TIR<TInput, PrecisionType> tirF;
    tirF.SetSAT(this->m_SAT);
    result[tir] = m_FvTIR->GetMaxVar(tirF(newPixel));
 
    MIRTIR<TInput, PrecisionType> mirtirF;
    mirtirF.SetSAT(this->m_SAT);
    result[mirtir] = m_FvMIRTIR->GetMaxVar(mirtirF(newPixel));
 
    result[ndsivis] = m_FvNDSIVis->GetMaxVar(NDSIVis<TInput, PrecisionType>()(newPixel));
 
    result[ndbbbi] = m_FvNDBBBI->GetMaxVar(NDBBBI<TInput, PrecisionType>()(newPixel));
 
    result[ndvi] = m_FvNDVI->GetMaxVar(NDVI<TInput, PrecisionType>()(newPixel));
 
    result[ndbsi] = m_FvNDBSI->GetMaxVar(NDBSI<TInput, PrecisionType>()(newPixel));
 
    return result;
  }
 
protected:
  typename FuzzyVarType::Pointer m_FvBright;
  typename FuzzyVarType::Pointer m_FvVis;
  typename FuzzyVarType::Pointer m_FvNIR;
  typename FuzzyVarType::Pointer m_FvMIR1;
  typename FuzzyVarType::Pointer m_FvMIR2;
  typename FuzzyVarType::Pointer m_FvTIR;
  typename FuzzyVarType::Pointer m_FvMIRTIR;
  typename FuzzyVarType::Pointer m_FvNDSIVis;
  typename FuzzyVarType::Pointer m_FvNDBBBI;
  typename FuzzyVarType::Pointer m_FvNDVI;
  typename FuzzyVarType::Pointer m_FvNDBSI;
};
 
 
template <class TInput, class TOutput>
class KernelSpectralRule : public LandsatTMIndexBase<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  virtual std::string GetName() const
  {
    return "LandsatTM KernelSpectralRule";
  }
 
  KernelSpectralRule() : m_TV1(0.7), m_TV2(0.5)
  {
  }
  ~KernelSpectralRule() override
  {
  }
 
  void SetTV1(PrecisionType tv1)
  {
    this->m_TV1 = tv1;
  }
 
  void SetTV2(PrecisionType tv2)
  {
    this->m_TV2 = tv2;
  }
 
  PrecisionType GetTV1() const
  {
    return this->m_TV1;
  }
 
  PrecisionType GetTV2() const
  {
    return this->m_TV2;
  }
 
protected:
  PrecisionType m_TV1;
 
  PrecisionType m_TV2;
 
  void SetMinMax(const TInput& inputPixel, PrecisionType* max13, PrecisionType* min123, PrecisionType* max123, PrecisionType* min12347, PrecisionType* max12347,
                 PrecisionType* max234, PrecisionType* max45)
  {
    std::vector<PrecisionType> v13;
    v13.push_back(inputPixel[this->m_TM1]);
    v13.push_back(inputPixel[this->m_TM3]);
 
    *max13 = *(std::max_element(v13.begin(), v13.end()));
 
 
    std::vector<PrecisionType> v123;
    v123.push_back(inputPixel[this->m_TM1]);
    v123.push_back(inputPixel[this->m_TM2]);
    v123.push_back(inputPixel[this->m_TM3]);
 
    *max123 = *(std::max_element(v123.begin(), v123.end()));
    *min123 = *(std::min_element(v123.begin(), v123.end()));
 
 
    std::vector<PrecisionType> v12347;
    v12347.push_back(inputPixel[this->m_TM1]);
    v12347.push_back(inputPixel[this->m_TM2]);
    v12347.push_back(inputPixel[this->m_TM3]);
    v12347.push_back(inputPixel[this->m_TM4]);
    v12347.push_back(inputPixel[this->m_TM7]);
 
    *max12347 = *(std::max_element(v12347.begin(), v12347.end()));
    *min12347 = *(std::min_element(v12347.begin(), v12347.end()));
 
    std::vector<PrecisionType> v234;
    v234.push_back(inputPixel[this->m_TM2]);
    v234.push_back(inputPixel[this->m_TM3]);
    v234.push_back(inputPixel[this->m_TM4]);
 
    *max234 = *(std::max_element(v234.begin(), v234.end()));
 
    std::vector<PrecisionType> v45;
    v45.push_back(inputPixel[this->m_TM4]);
    v45.push_back(inputPixel[this->m_TM5]);
 
    *max45 = *(std::max_element(v45.begin(), v45.end()));
  }
};
 
template <class TInput, class TOutput>
class ThickCloudsSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM ThickCloudsSpectralRule";
  }
 
  ThickCloudsSpectralRule()
  {
  }
  ~ThickCloudsSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput        newPixel(this->PrepareValues(inputPixel));
    PrecisionType max13;
    PrecisionType max123;
    PrecisionType min123;
    PrecisionType max12347;
    PrecisionType min12347;
    PrecisionType max234;
    PrecisionType max45;
    this->SetMinMax(newPixel, &max13, &min123, &max123, &min12347, &max12347, &max234, &max45);
 
    bool result = (((min123 >= (this->m_TV1 * max123) && (max123 <= this->m_TV1 * newPixel[this->m_TM4])) ||
                    ((newPixel[this->m_TM2] >= this->m_TV1 * max13) && (max123 <= newPixel[this->m_TM4]))) &&
                   (newPixel[this->m_TM5] <= this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM5] >= this->m_TV1 * max123) &&
                   (newPixel[this->m_TM7] <= this->m_TV1 * newPixel[this->m_TM4]));
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class ThinCloudsSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM ThinCloudsSpectralRule";
  }
 
  ThinCloudsSpectralRule()
  {
  }
  ~ThinCloudsSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput        newPixel(this->PrepareValues(inputPixel));
    PrecisionType max13;
    PrecisionType max123;
    PrecisionType min123;
    PrecisionType max12347;
    PrecisionType min12347;
    PrecisionType max234;
    PrecisionType max45;
    this->SetMinMax(newPixel, &max13, &min123, &max123, &min12347, &max12347, &max234, &max45);
 
 
    bool result = (min123 >= (this->m_TV1 * max123)) && (newPixel[this->m_TM4] >= max123) &&
                  !((newPixel[this->m_TM1] <= newPixel[this->m_TM2] && newPixel[this->m_TM2] <= newPixel[this->m_TM3] &&
                     newPixel[this->m_TM3] <= newPixel[this->m_TM4]) &&
                    (newPixel[this->m_TM3] >= this->m_TV1 * newPixel[this->m_TM4])) &&
                  (newPixel[this->m_TM4] >= this->m_TV1 * newPixel[this->m_TM5]) && (newPixel[this->m_TM5] >= this->m_TV1 * newPixel[this->m_TM4]) &&
                  (newPixel[this->m_TM5] >= this->m_TV1 * max123) && (newPixel[this->m_TM5] >= this->m_TV1 * newPixel[this->m_TM7]);
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class SnowOrIceSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM SnowOrIceSpectralRule";
  }
 
  SnowOrIceSpectralRule()
  {
  }
  ~SnowOrIceSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput        newPixel(this->PrepareValues(inputPixel));
    PrecisionType max13;
    PrecisionType max123;
    PrecisionType min123;
    PrecisionType max12347;
    PrecisionType min12347;
    PrecisionType max234;
    PrecisionType max45;
    this->SetMinMax(newPixel, &max13, &min123, &max123, &min12347, &max12347, &max234, &max45);
 
 
    bool result = (min123 >= (this->m_TV1 * max123)) && (newPixel[this->m_TM4] >= (this->m_TV1 * max123)) &&
                  (newPixel[this->m_TM5] <= this->m_TV2 * newPixel[this->m_TM4]) && (newPixel[this->m_TM5] <= this->m_TV1 * min123) &&
                  (newPixel[this->m_TM5] <= this->m_TV1 * max123) && (newPixel[this->m_TM7] <= this->m_TV2 * newPixel[this->m_TM4]) &&
                  (newPixel[this->m_TM7] <= this->m_TV1 * min123);
 
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class WaterOrShadowSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM WaterOrShadowSpectralRule";
  }
 
  WaterOrShadowSpectralRule()
  {
  }
  ~WaterOrShadowSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    bool   result = (newPixel[this->m_TM1] >= newPixel[this->m_TM2]) && (newPixel[this->m_TM2] >= newPixel[this->m_TM3]) &&
                  (newPixel[this->m_TM3] >= newPixel[this->m_TM4]) && (newPixel[this->m_TM4] >= newPixel[this->m_TM5]) &&
                  (newPixel[this->m_TM4] >= newPixel[this->m_TM7]);
 
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class PitbogOrGreenhouseSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM PitbogOrGreenhouseSpectralRule";
  }
 
  PitbogOrGreenhouseSpectralRule()
  {
  }
  ~PitbogOrGreenhouseSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput        newPixel(this->PrepareValues(inputPixel));
    PrecisionType max13;
    PrecisionType max123;
    PrecisionType min123;
    PrecisionType max12347;
    PrecisionType min12347;
    PrecisionType max234;
    PrecisionType max45;
    this->SetMinMax(newPixel, &max13, &min123, &max123, &min12347, &max12347, &max234, &max45);
 
 
    bool result = (newPixel[this->m_TM3] >= this->m_TV1 * newPixel[this->m_TM1]) && (newPixel[this->m_TM1] >= this->m_TV1 * newPixel[this->m_TM3]) &&
                  (max123 <= this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM5] <= this->m_TV1 * newPixel[this->m_TM4]) &&
                  (newPixel[this->m_TM3] >= this->m_TV2 * newPixel[this->m_TM5]) && (min123 >= this->m_TV1 * newPixel[this->m_TM7]);
 
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class DominantBlueSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM DominantBlueSpectralRule";
  }
 
  DominantBlueSpectralRule()
  {
  }
  ~DominantBlueSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    bool   result = (newPixel[this->m_TM1] >= this->m_TV1 * newPixel[this->m_TM2]) && (newPixel[this->m_TM1] >= this->m_TV1 * newPixel[this->m_TM3]) &&
                  (newPixel[this->m_TM1] >= this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM1] >= this->m_TV1 * newPixel[this->m_TM5]) &&
                  (newPixel[this->m_TM1] >= this->m_TV1 * newPixel[this->m_TM7]);
 
 
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class VegetationSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM VegetationSpectralRule";
  }
 
  VegetationSpectralRule()
  {
  }
  ~VegetationSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput        newPixel(this->PrepareValues(inputPixel));
    PrecisionType max13;
    PrecisionType max123;
    PrecisionType min123;
    PrecisionType max12347;
    PrecisionType min12347;
    PrecisionType max234;
    PrecisionType max45;
    this->SetMinMax(newPixel, &max13, &min123, &max123, &min12347, &max12347, &max234, &max45);
 
 
    bool result = (newPixel[this->m_TM2] >= this->m_TV2 * newPixel[this->m_TM1]) && (newPixel[this->m_TM2] >= this->m_TV1 * newPixel[this->m_TM3]) &&
                  (newPixel[this->m_TM3] < this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM4] > max123) &&
                  (newPixel[this->m_TM5] < this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM5] >= this->m_TV1 * newPixel[this->m_TM3]) &&
                  (newPixel[this->m_TM7] < this->m_TV1 * newPixel[this->m_TM5]);
 
 
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class RangelandSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM RangelandSpectralRule";
  }
 
  RangelandSpectralRule()
  {
  }
  ~RangelandSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput        newPixel(this->PrepareValues(inputPixel));
    PrecisionType max13;
    PrecisionType max123;
    PrecisionType min123;
    PrecisionType max12347;
    PrecisionType min12347;
    PrecisionType max234;
    PrecisionType max45;
    this->SetMinMax(newPixel, &max13, &min123, &max123, &min12347, &max12347, &max234, &max45);
 
 
    bool result = (newPixel[this->m_TM2] >= this->m_TV2 * newPixel[this->m_TM1]) && (newPixel[this->m_TM2] >= this->m_TV1 * newPixel[this->m_TM3]) &&
                  (newPixel[this->m_TM4] > max123) && (newPixel[this->m_TM3] < this->m_TV1 * newPixel[this->m_TM4]) &&
                  (newPixel[this->m_TM4] >= this->m_TV1 * newPixel[this->m_TM5]) && (newPixel[this->m_TM5] >= this->m_TV1 * newPixel[this->m_TM4]) &&
                  (newPixel[this->m_TM5] > max123) && (newPixel[this->m_TM7] < this->m_TV1 * max45) && (newPixel[this->m_TM5] >= newPixel[this->m_TM7]);
 
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class BarrenLandOrBuiltUpOrCloudsSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM BarrenLandOrBuiltUpOrCloudsSpectralRule";
  }
 
  BarrenLandOrBuiltUpOrCloudsSpectralRule()
  {
  }
  ~BarrenLandOrBuiltUpOrCloudsSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput        newPixel(this->PrepareValues(inputPixel));
    PrecisionType max13;
    PrecisionType max123;
    PrecisionType min123;
    PrecisionType max12347;
    PrecisionType min12347;
    PrecisionType max234;
    PrecisionType max45;
    this->SetMinMax(newPixel, &max13, &min123, &max123, &min12347, &max12347, &max234, &max45);
 
 
    bool result = (newPixel[this->m_TM3] >= this->m_TV2 * newPixel[this->m_TM1]) && (newPixel[this->m_TM3] >= this->m_TV1 * newPixel[this->m_TM2]) &&
                  (newPixel[this->m_TM4] >= this->m_TV1 * max123) && (newPixel[this->m_TM5] >= max123) &&
                  (newPixel[this->m_TM5] >= this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM5] >= this->m_TV1 * newPixel[this->m_TM7]) &&
                  (newPixel[this->m_TM7] >= this->m_TV2 * max45);
 
    return static_cast<TOutput>(result);
  }
};
 
template <class TInput, class TOutput>
class FlatResponseBarrenLandOrBuiltUpSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM FlatResponseBarrenLandOrBuiltUpSpectralRule";
  }
 
  FlatResponseBarrenLandOrBuiltUpSpectralRule()
  {
  }
  ~FlatResponseBarrenLandOrBuiltUpSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput        newPixel(this->PrepareValues(inputPixel));
    PrecisionType max13;
    PrecisionType max123;
    PrecisionType min123;
    PrecisionType max12347;
    PrecisionType min12347;
    PrecisionType max234;
    PrecisionType max45;
    this->SetMinMax(newPixel, &max13, &min123, &max123, &min12347, &max12347, &max234, &max45);
 
 
    bool result = (newPixel[this->m_TM5] >= this->m_TV1 * max12347) && (min12347 >= this->m_TV2 * newPixel[this->m_TM5]);
 
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class ShadowWithBarrenLandSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM ShadowWithBarrenLandSpectralRule";
  }
 
  ShadowWithBarrenLandSpectralRule()
  {
  }
  ~ShadowWithBarrenLandSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    bool   result = (newPixel[this->m_TM1] >= newPixel[this->m_TM2]) && (newPixel[this->m_TM2] >= newPixel[this->m_TM3]) &&
                  (newPixel[this->m_TM3] >= this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM1] >= newPixel[this->m_TM5]) &&
                  (newPixel[this->m_TM5] >= this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM5] >= this->m_TV1 * newPixel[this->m_TM7]);
 
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class ShadowWithVegetationSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM ShadowWithVegetationSpectralRule";
  }
 
  ShadowWithVegetationSpectralRule()
  {
  }
  ~ShadowWithVegetationSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    bool   result = (newPixel[this->m_TM1] >= newPixel[this->m_TM2]) && (newPixel[this->m_TM2] >= newPixel[this->m_TM3]) &&
                  (newPixel[this->m_TM1] >= this->m_TV2 * newPixel[this->m_TM4]) && (newPixel[this->m_TM3] < this->m_TV1 * newPixel[this->m_TM4]) &&
                  (newPixel[this->m_TM5] < this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM3] >= this->m_TV2 * newPixel[this->m_TM5]) &&
                  (newPixel[this->m_TM7] < this->m_TV1 * newPixel[this->m_TM4]);
 
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class ShadowCloudOrSnowSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM ShadowCloudOrSnowSpectralRule";
  }
 
  ShadowCloudOrSnowSpectralRule()
  {
  }
  ~ShadowCloudOrSnowSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput        newPixel(this->PrepareValues(inputPixel));
    PrecisionType max13;
    PrecisionType max123;
    PrecisionType min123;
    PrecisionType max12347;
    PrecisionType min12347;
    PrecisionType max234;
    PrecisionType max45;
    this->SetMinMax(newPixel, &max13, &min123, &max123, &min12347, &max12347, &max234, &max45);
 
 
    bool result = (newPixel[this->m_TM1] >= this->m_TV1 * max234) && (max234 >= this->m_TV1 * newPixel[this->m_TM1]) &&
                  (newPixel[this->m_TM5] < newPixel[this->m_TM1]) && (newPixel[this->m_TM7] < this->m_TV1 * newPixel[this->m_TM1]);
 
    return static_cast<TOutput>(result);
  }
};
 
 
template <class TInput, class TOutput>
class WetlandSpectralRule : public KernelSpectralRule<TInput, TOutput>
{
public:
  typedef typename TInput::ValueType PrecisionType;
  typedef bool                       OutputPixelType;
 
  std::string GetName() const override
  {
    return "LandsatTM WetlandSpectralRule";
  }
 
  WetlandSpectralRule()
  {
  }
  ~WetlandSpectralRule() override
  {
  }
 
  inline TOutput operator()(const TInput& inputPixel)
  {
    TInput newPixel(this->PrepareValues(inputPixel));
    bool   result = (newPixel[this->m_TM1] >= newPixel[this->m_TM2]) && (newPixel[this->m_TM2] >= newPixel[this->m_TM3]) &&
                  (newPixel[this->m_TM1] >= this->m_TV1 * newPixel[this->m_TM4]) && (newPixel[this->m_TM3] < newPixel[this->m_TM4]) &&
                  (newPixel[this->m_TM4] >= this->m_TV1 * newPixel[this->m_TM5]) && (newPixel[this->m_TM5] >= this->m_TV1 * newPixel[this->m_TM4]) &&
                  (newPixel[this->m_TM3] >= this->m_TV2 * newPixel[this->m_TM5]) && (newPixel[this->m_TM5] >= newPixel[this->m_TM7]);
 
    return static_cast<TOutput>(result);
  }
};
 
 
} // namespace LandsatTM
} // namespace Functor
} // namespace otb
 
#endif