otbRadianceToReflectanceImageFilter.h

Go to the documentation of this file.

/*
 * Copyright (C) 1999-2011 Insight Software Consortium
 * 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 otbRadianceToReflectanceImageFilter_h
#define otbRadianceToReflectanceImageFilter_h
 
#include "otbVarSol.h"
#include "otbUnaryImageFunctorWithVectorImageFilter.h"
#include "otbMacro.h"
#include <iomanip>
 
namespace otb
{
namespace Functor
{
template <class TInput, class TOutput>
class RadianceToReflectanceImageFunctor
{
public:
  RadianceToReflectanceImageFunctor() : m_SolarIllumination(1.0), m_IlluminationCorrectionCoefficient(1.0), m_UseClamp(true)
  {
  }
 
  virtual ~RadianceToReflectanceImageFunctor()
  {
  }
 
  void SetSolarIllumination(double solarIllumination)
  {
    m_SolarIllumination = solarIllumination;
  }
  void SetIlluminationCorrectionCoefficient(double coef)
  {
    m_IlluminationCorrectionCoefficient = coef;
  }
  void SetUseClamp(bool useClamp)
  {
    m_UseClamp = useClamp;
  }
 
  double GetSolarIllumination()
  {
    return m_SolarIllumination;
  }
  double GetIlluminationCorrectionCoefficient()
  {
    return m_IlluminationCorrectionCoefficient;
  }
  bool GetUseClamp()
  {
    return m_UseClamp;
  }
 
  inline TOutput operator()(const TInput& inPixel) const
  {
    TOutput outPixel;
    double  temp;
    temp = static_cast<double>(inPixel) * static_cast<double>(CONST_PI) * m_IlluminationCorrectionCoefficient / m_SolarIllumination;
 
    if (m_UseClamp)
    {
      temp = std::max(temp, 0.);
      temp = std::min(temp, 1.);
    }
    outPixel = static_cast<TOutput>(temp);
 
    return outPixel;
  }
 
private:
  double m_SolarIllumination;
  double m_IlluminationCorrectionCoefficient;
  bool   m_UseClamp;
};
}
 
template <class TInputImage, class TOutputImage>
class ITK_EXPORT RadianceToReflectanceImageFilter
    : public UnaryImageFunctorWithVectorImageFilter<
          TInputImage, TOutputImage,
          typename Functor::RadianceToReflectanceImageFunctor<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::RadianceToReflectanceImageFunctor<typename InputImageType::InternalPixelType, typename OutputImageType::InternalPixelType>
      FunctorType;
 
  typedef RadianceToReflectanceImageFilter Self;
  typedef UnaryImageFunctorWithVectorImageFilter<InputImageType, OutputImageType, FunctorType> Superclass;
  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;
 
  itkNewMacro(Self);
 
  itkTypeMacro(RadianceToReflectanceImageFilter, 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(SolarIllumination, VectorType);
 
  itkGetConstReferenceMacro(SolarIllumination, VectorType);
 
  itkSetMacro(ZenithalSolarAngle, double);
 
  itkGetConstReferenceMacro(ZenithalSolarAngle, double);
 
  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;
  }
 
  itkSetClampMacro(Day, int, 1, 31);
 
  itkGetConstReferenceMacro(Day, int);
 
  itkSetClampMacro(Month, int, 1, 12);
 
  itkGetConstReferenceMacro(Month, int);
 
  void SetFluxNormalizationCoefficient(double coef)
  {
    m_FluxNormalizationCoefficient      = coef;
    m_IsSetFluxNormalizationCoefficient = true;
    this->Modified();
  }
 
  itkGetConstReferenceMacro(FluxNormalizationCoefficient, double);
 
  void SetSolarDistance(double value)
  {
    m_SolarDistance      = value;
    m_IsSetSolarDistance = true;
    this->Modified();
  }
 
  itkGetConstReferenceMacro(SolarDistance, double);
 
  itkSetMacro(IsSetSolarDistance, bool);
 
  itkGetConstReferenceMacro(IsSetSolarDistance, bool);
 
  itkSetMacro(IsSetFluxNormalizationCoefficient, bool);
 
  itkGetConstReferenceMacro(IsSetFluxNormalizationCoefficient, bool);
 
  itkSetMacro(UseClamp, bool);
 
  itkGetConstReferenceMacro(UseClamp, bool);
 
protected:
  RadianceToReflectanceImageFilter()
    : m_ZenithalSolarAngle(120.0), // invalid value which will lead to negative radiometry
      m_FluxNormalizationCoefficient(1.),
      m_Day(0),
      m_Month(0),
      m_SolarDistance(1.0),
      m_IsSetFluxNormalizationCoefficient(false),
      m_IsSetSolarDistance(false),
      m_UseClamp(true)
  {
    m_SolarIllumination.SetSize(0);
  };
 
  ~RadianceToReflectanceImageFilter() override
  {
  }
 
  void BeforeThreadedGenerateData(void) override
  {
    const auto & metadata = this->GetInput()->GetImageMetadata();
 
    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(<< "Day:               " << m_Day);
    otbMsgDevMacro(<< "Month:             " << m_Month);
    otbMsgDevMacro(<< "Solar irradiance:  " << m_SolarIllumination);
    otbMsgDevMacro(<< "Zenithal angle:    " << m_ZenithalSolarAngle);
 
    if ((m_SolarIllumination.GetSize() != this->GetInput()->GetNumberOfComponentsPerPixel()))
    {
      itkExceptionMacro(<< "SolarIllumination parameter 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.SetSolarIllumination(static_cast<double>(m_SolarIllumination[i]));
      functor.SetUseClamp(m_UseClamp);
 
      this->GetFunctorVector().push_back(functor);
    }
  }
 
private:
  double m_ZenithalSolarAngle;
 
  double m_FluxNormalizationCoefficient;
 
  int m_Day;
 
  int m_Month;
 
  double m_SolarDistance;
 
  VectorType m_SolarIllumination;
 
  bool m_IsSetFluxNormalizationCoefficient;
 
  bool m_IsSetSolarDistance;
 
  bool m_UseClamp;
};
 
} // end namespace otb
#endif