otbImageToReflectanceImageFilter.h

Go to the documentation of this file.

/*
 * Copyright (C) 1999-2011 Insight Software Consortium
 * Copyright (C) 2005-2024 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 otbImageToReflectanceImageFilter_h
#define otbImageToReflectanceImageFilter_h
 
#include "otbImageToRadianceImageFilter.h"
#include "otbRadianceToReflectanceImageFilter.h"
 
namespace otb
{
namespace Functor
{
template <class TInput, class TOutput>
class ImageToReflectanceImageFunctor
{
public:
  ImageToReflectanceImageFunctor()
  {
  }
  virtual ~ImageToReflectanceImageFunctor()
  {
  }
 
  typedef Functor::ImageToRadianceImageFunctor<TInput, TOutput>       ImToLumFunctorType;
  typedef Functor::RadianceToReflectanceImageFunctor<TInput, TOutput> LumToReflecFunctorType;
 
  void SetAlpha(double alpha)
  {
    m_ImToLumFunctor.SetAlpha(alpha);
  }
  void SetBeta(double beta)
  {
    m_ImToLumFunctor.SetBeta(beta);
  }
  void SetSolarIllumination(double solarIllumination)
  {
    m_LumToReflecFunctor.SetSolarIllumination(solarIllumination);
  }
  void SetIlluminationCorrectionCoefficient(double coef)
  {
    m_LumToReflecFunctor.SetIlluminationCorrectionCoefficient(coef);
  }
  void SetUseClamp(bool useClamp)
  {
    m_LumToReflecFunctor.SetUseClamp(useClamp);
  }
 
  double GetAlpha()
  {
    return m_ImToLumFunctor.GetAlpha();
  }
  double GetBeta()
  {
    return m_ImToLumFunctor.GetBeta();
  }
  double GetSolarIllumination()
  {
    return m_LumToReflecFunctor.GetSolarIllumination();
  }
  double GetIlluminationCorrectionCoefficient()
  {
    return m_LumToReflecFunctor.GetIlluminationCorrectionCoefficient();
  }
  bool GetUseClamp()
  {
    return m_LumToReflecFunctor.GetUseClamp();
  }
 
  inline TOutput operator()(const TInput& inPixel) const
  {
    TOutput outPixel;
    TOutput tempPix;
    tempPix  = m_ImToLumFunctor(inPixel);
    outPixel = m_LumToReflecFunctor(tempPix);
 
    return outPixel;
  }
 
private:
  ImToLumFunctorType     m_ImToLumFunctor;
  LumToReflecFunctorType m_LumToReflecFunctor;
};
}
 
template <class TInputImage, class TOutputImage>
class ITK_EXPORT ImageToReflectanceImageFilter
    : public UnaryImageFunctorWithVectorImageFilter<
          TInputImage, TOutputImage,
          typename Functor::ImageToReflectanceImageFunctor<typename TInputImage::InternalPixelType, typename TOutputImage::InternalPixelType>>
{
public:
  itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
  itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension);
 
  typedef TInputImage  InputImageType;
  typedef TOutputImage OutputImageType;
  typedef typename Functor::ImageToReflectanceImageFunctor<typename InputImageType::InternalPixelType, typename OutputImageType::InternalPixelType> FunctorType;
 
  typedef ImageToReflectanceImageFilter Self;
  typedef UnaryImageFunctorWithVectorImageFilter<InputImageType, OutputImageType, FunctorType> Superclass;
  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;
 
  itkNewMacro(Self);
 
  itkTypeMacro(ImageToReflectanceImageFilter, UnaryImageFunctorWithVectorImageFiltermageFilter);
 
  typedef typename InputImageType::PixelType          InputPixelType;
  typedef typename InputImageType::InternalPixelType  InputInternalPixelType;
  typedef typename InputImageType::RegionType         InputImageRegionType;
  typedef typename OutputImageType::PixelType         OutputPixelType;
  typedef typename OutputImageType::InternalPixelType OutputInternalPixelType;
  typedef typename OutputImageType::RegionType        OutputImageRegionType;
 
  typedef typename itk::VariableLengthVector<double> VectorType;
 
  typedef typename InputImageType::SizeType SizeType;
 
  itkSetMacro(Alpha, VectorType);
 
  itkGetConstReferenceMacro(Alpha, VectorType);
 
  itkSetMacro(Beta, VectorType);
 
  itkGetConstReferenceMacro(Beta, VectorType);
 
  itkSetMacro(SolarIllumination, VectorType);
 
  itkGetConstReferenceMacro(SolarIllumination, VectorType);
 
  itkSetMacro(ZenithalSolarAngle, double);
 
  itkGetConstReferenceMacro(ZenithalSolarAngle, double);
 
  itkSetMacro(UseClamp, bool);
 
  itkGetConstReferenceMacro(UseClamp, bool);
 
  virtual void SetElevationSolarAngle(double elevationAngle)
  {
    double zenithalAngle = 90.0 - elevationAngle;
    if (this->m_ZenithalSolarAngle != zenithalAngle)
    {
      this->m_ZenithalSolarAngle = zenithalAngle;
      this->Modified();
    }
  }
 
  virtual double GetElevationSolarAngle() const
  {
    return 90.0 - this->m_ZenithalSolarAngle;
  }
 
  void SetFluxNormalizationCoefficient(double coef)
  {
    m_FluxNormalizationCoefficient      = coef;
    m_IsSetFluxNormalizationCoefficient = true;
    this->Modified();
  }
 
  void SetSolarDistance(double value)
  {
    m_SolarDistance      = value;
    m_IsSetSolarDistance = true;
    this->Modified();
  }
 
  itkGetConstReferenceMacro(SolarDistance, double);
 
  itkSetMacro(IsSetSolarDistance, bool);
 
  itkGetConstReferenceMacro(IsSetSolarDistance, bool);
 
  itkSetClampMacro(Day, int, 1, 31);
 
  itkGetConstReferenceMacro(Day, int);
 
  itkSetClampMacro(Month, int, 1, 12);
 
  itkGetConstReferenceMacro(Month, int);
 
protected:
  ImageToReflectanceImageFilter()
    : m_ZenithalSolarAngle(120.), // invalid value which will lead to negative radiometry
      m_FluxNormalizationCoefficient(1.),
      m_UseClamp(true),
      m_IsSetFluxNormalizationCoefficient(false),
      m_Day(0),
      m_Month(0),
      m_SolarDistance(1.0),
      m_IsSetSolarDistance(false)
  {
    m_Alpha.SetSize(0);
    m_Beta.SetSize(0);
    m_SolarIllumination.SetSize(0);
  };
 
  ~ImageToReflectanceImageFilter() override
  {
  }
 
  void BeforeThreadedGenerateData(void) override
  {
    const auto & metadata = this->GetInput()->GetImageMetadata();
 
    if (m_Alpha.GetSize() == 0 && metadata.HasBandMetadata(MDNum::PhysicalGain))
    {
      m_Alpha = metadata.GetAsVector(MDNum::PhysicalGain);
    }
 
    if (m_Beta.GetSize() == 0 && metadata.HasBandMetadata(MDNum::PhysicalBias))
    {
      m_Beta = metadata.GetAsVector(MDNum::PhysicalBias);
    }
 
    if (m_Day == 0 && (!m_IsSetFluxNormalizationCoefficient) && (!m_IsSetSolarDistance) 
        && metadata.Has(MDTime::AcquisitionDate))
    {
      m_Day = metadata[MDTime::AcquisitionDate].GetDay();
    }
 
    if (m_Month == 0 && (!m_IsSetFluxNormalizationCoefficient) && (!m_IsSetSolarDistance) 
        && metadata.Has(MDTime::AcquisitionDate))
    {
      m_Month = metadata[MDTime::AcquisitionDate].GetMonth();
    }
 
    if (m_SolarIllumination.GetSize() == 0 && metadata.HasBandMetadata(MDNum::SolarIrradiance))
    {
      m_SolarIllumination = metadata.GetAsVector(MDNum::SolarIrradiance);
    }
 
    if (m_ZenithalSolarAngle == 120.0 && metadata.Has(MDNum::SunElevation))
    {
      // the zenithal angle is the complementary of the elevation angle
      m_ZenithalSolarAngle = 90.0 - metadata[MDNum::SunElevation];
    }
 
    otbMsgDevMacro(<< "Using correction parameters: ");
    otbMsgDevMacro(<< "Alpha (gain): " << m_Alpha);
    otbMsgDevMacro(<< "Beta (bias):  " << m_Beta);
    otbMsgDevMacro(<< "Day:               " << m_Day);
    otbMsgDevMacro(<< "Month:             " << m_Month);
    otbMsgDevMacro(<< "Solar irradiance:  " << m_SolarIllumination);
    otbMsgDevMacro(<< "Zenithal angle:    " << m_ZenithalSolarAngle);
 
    if ((m_Alpha.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel()) || (m_Beta.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel()) ||
        (m_SolarIllumination.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel()))
    {
      itkExceptionMacro(<< "Alpha, Beta and SolarIllumination parameters should have the same size as the number of bands");
    }
 
    this->GetFunctorVector().clear();
    for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i)
    {
      FunctorType functor;
      double      coefTemp = 0.;
 
      if (m_IsSetFluxNormalizationCoefficient)
      {
        coefTemp = std::cos(m_ZenithalSolarAngle * CONST_PI_180) * m_FluxNormalizationCoefficient * m_FluxNormalizationCoefficient;
      }
      else if (m_IsSetSolarDistance)
      {
        coefTemp = std::cos(m_ZenithalSolarAngle * CONST_PI_180) / (m_SolarDistance * m_SolarDistance);
      }
      else if (m_Day * m_Month != 0 && m_Day < 32 && m_Month < 13)
      {
        coefTemp = std::cos(m_ZenithalSolarAngle * CONST_PI_180) * VarSol::GetVarSol(m_Day, m_Month);
      }
      else
      {
        itkExceptionMacro(<< "Day has to be included between 1 and 31, Month between 1 and 12.");
      }
 
      functor.SetIlluminationCorrectionCoefficient(1. / coefTemp);
      functor.SetAlpha(m_Alpha[i]);
      functor.SetBeta(m_Beta[i]);
      functor.SetSolarIllumination(m_SolarIllumination[i]);
      functor.SetUseClamp(m_UseClamp);
      this->GetFunctorVector().push_back(functor);
    }
  }
 
private:
  VectorType m_Alpha;
  VectorType m_Beta;
 
  double m_ZenithalSolarAngle;
 
  double m_FluxNormalizationCoefficient;
 
  VectorType m_SolarIllumination;
 
  bool m_UseClamp;
 
  bool m_IsSetFluxNormalizationCoefficient;
 
  int m_Day;
 
  int m_Month;
 
  double m_SolarDistance;
 
  bool m_IsSetSolarDistance;
};
 
} // end namespace otb
 
#endif