ADD: keep old radiometry classes, flagged as deprecated

Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.cxx

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#include "otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.h"
+#include "otbSIXSTraits.h"
+
+namespace otb
+{
+/**
+ * Constructor.
+ */
+
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms()
+{
+  this->ProcessObject::SetNumberOfRequiredInputs(1);
+  this->ProcessObject::SetNumberOfRequiredOutputs(1);
+  // Create the output. We use static_cast<> here because we know the default
+  // output must be of type TOutputPointSet
+  AtmosphericRadiativeTermsPointer output
+    = static_cast<AtmosphericRadiativeTermsType*>(this->MakeOutput(0).GetPointer());
+
+  this->ProcessObject::SetNthOutput(0, output.GetPointer());
+
+}
+/**
+ *
+ */
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms::DataObjectPointer
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::MakeOutput(unsigned int)
+{
+  return static_cast<itk::DataObject*>(AtmosphericRadiativeTermsType::New().GetPointer());
+}
+
+/**
+ *
+ */
+void
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::GraftOutput(itk::DataObject *graft)
+{
+  this->GraftNthOutput(0, graft);
+}
+
+/**
+ *
+ */
+void
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::GraftNthOutput(unsigned int idx, itk::DataObject *graft)
+{
+  if (idx >= this->GetNumberOfOutputs())
+    {
+    itkExceptionMacro(<< "Requested to graft output " << idx <<
+                      " but this filter only has " << this->GetNumberOfOutputs() << " Outputs.");
+    }
+
+  if (!graft)
+    {
+    itkExceptionMacro(<< "Requested to graft output that is a NULL pointer");
+    }
+
+  itk::DataObject * output = this->GetOutput(idx);
+
+  // Call Graft on the PointSet in order to copy meta-information, and containers.
+  output->Graft(graft);
+}
+
+/**
+ * Get the output data
+ * \return The data produced.
+ */
+
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms::AtmosphericRadiativeTermsType *
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::GetOutput(void)
+{
+  if (this->GetNumberOfOutputs() < 1)
+    {
+    return 0;
+    }
+  return static_cast<AtmosphericRadiativeTermsType *> (this->ProcessObject::GetOutput(0));
+}
+
+/**
+ *
+ */
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms::AtmosphericRadiativeTermsType *
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::GetOutput(unsigned int idx)
+{
+  return static_cast<AtmosphericRadiativeTermsType*>
+         (this->itk::ProcessObject::GetOutput(idx));
+}
+
+void
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::SetInput(const AtmosphericCorrectionParametersType *object)
+{
+  // A single input image
+  this->itk::ProcessObject::SetNthInput(0, const_cast<AtmosphericCorrectionParametersType*>(object));
+}
+
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms::AtmosphericCorrectionParametersType *
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::GetInput(void)
+{
+  // If there is no input
+  if (this->GetNumberOfInputs() != 1)
+    {
+    // exit
+    return 0;
+    }
+  // else return the first input
+  return static_cast<AtmosphericCorrectionParametersType *>
+         (this->itk::ProcessObject::GetInput(0));
+}
+
+void
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::GenerateData()
+{
+
+  AtmosphericCorrectionParametersPointer input = this->GetInput();
+  AtmosphericRadiativeTermsPointer       output = this->GetOutput();
+
+  output->GetValues().clear();
+  typedef AtmosphericCorrectionParameters::WavelengthSpectralBandVectorType WavelengthSpectralBandVectorType;
+  WavelengthSpectralBandVectorType WavelengthSpectralBandVector = input->GetWavelengthSpectralBand();
+  unsigned int                     NbBand = WavelengthSpectralBandVector->Size();
+
+  double atmosphericReflectance(0.);
+  double atmosphericSphericalAlbedo(0.);
+  double totalGaseousTransmission(0.);
+  double downwardTransmittance(0.);
+  double upwardTransmittance(0.);
+  double upwardDiffuseTransmittance(0.);
+  double upwardDirectTransmittance(0.);
+  double upwardDiffuseTransmittanceForRayleigh(0.);
+  double upwardDiffuseTransmittanceForAerosol(0.);
+
+  for (unsigned int i = 0; i < NbBand; ++i)
+    {
+    atmosphericReflectance = 0.;
+    atmosphericSphericalAlbedo = 0.;
+    totalGaseousTransmission = 0.;
+    downwardTransmittance = 0.;
+    upwardTransmittance = 0.;
+    upwardDiffuseTransmittance = 0.;
+    upwardDirectTransmittance = 0.;
+    upwardDiffuseTransmittanceForRayleigh = 0.;
+    upwardDiffuseTransmittanceForAerosol = 0.;
+    SIXSTraits::ComputeAtmosphericParameters(
+      input->GetSolarZenithalAngle(),                  /** The Solar zenithal angle */
+      input->GetSolarAzimutalAngle(),                  /** The Solar azimutal angle */
+      input->GetViewingZenithalAngle(),                /** The Viewing zenithal angle */
+      input->GetViewingAzimutalAngle(),                /** The Viewing azimutal angle */
+      input->GetMonth(),                               /** The Month */
+      input->GetDay(),                                 /** The Day (in the month) */
+      input->GetAtmosphericPressure(),                 /** The Atmospheric pressure */
+      input->GetWaterVaporAmount(),                    /** The Water vapor amount (Total water vapor content over vertical atmospheric column) */
+      input->GetOzoneAmount(),                         /** The Ozone amount (Stratospheric ozone layer content) */
+      input->GetAerosolModel(),                        /** The Aerosol model */
+      input->GetAerosolOptical(),                      /** The Aerosol optical (radiative impact of aerosol for the reference wavelength 550-nm) */
+      input->GetWavelengthSpectralBand()->GetNthElement(i), /** Wavelength for the spectral band definition */
+      /** Note : The Max wavelength spectral band value must be updated ! */
+      atmosphericReflectance,                          /** Atmospheric reflectance */
+      atmosphericSphericalAlbedo,                      /** atmospheric spherical albedo */
+      totalGaseousTransmission,                        /** Total gaseous transmission */
+      downwardTransmittance,                           /** downward transmittance */
+      upwardTransmittance,                             /** upward transmittance */
+      upwardDiffuseTransmittance,                      /** Upward diffuse transmittance */
+      upwardDirectTransmittance,                       /** Upward direct transmittance */
+      upwardDiffuseTransmittanceForRayleigh,           /** Upward diffuse transmittance for rayleigh */
+      upwardDiffuseTransmittanceForAerosol             /** Upward diffuse transmittance for aerosols */
+      );
+
+    output->SetIntrinsicAtmosphericReflectance(i, atmosphericReflectance);
+    output->SetSphericalAlbedo(i, atmosphericSphericalAlbedo);
+    output->SetTotalGaseousTransmission(i, totalGaseousTransmission);
+    output->SetDownwardTransmittance(i, downwardTransmittance);
+    output->SetUpwardTransmittance(i, upwardTransmittance);
+    output->SetUpwardDiffuseTransmittance(i, upwardDiffuseTransmittance);
+    output->SetUpwardDirectTransmittance(i, upwardDirectTransmittance);
+    output->SetUpwardDiffuseTransmittanceForRayleigh(i, upwardDiffuseTransmittanceForRayleigh);
+    output->SetUpwardDiffuseTransmittanceForAerosol(i, upwardDiffuseTransmittanceForAerosol);
+    output->SetWavelengthSpectralBand(i, input->GetWavelengthSpectralBand()->GetNthElement(i)->GetCenterSpectralValue());
+    }
+}
+
+/**
+ * PrintSelf Method
+ */
+void
+AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  Superclass::PrintSelf(os, indent);
+}
+
+} // end namespace otb

Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.h

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms_h
+#define __otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms_h
+
+#include "otbMacro.h"
+#include "itkProcessObject.h"
+#include "otbAtmosphericCorrectionParameters.h"
+#include "otbAtmosphericRadiativeTerms.h"
+
+namespace otb
+{
+/**
+ * \class AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+ * \brief This class computes the atmospheric radiative terms with 6S.
+ *
+ * It enables to compute a AtmosphericRadiativeTerms from a AtmosphericCorrectionParameters,
+ * which is used in the ReflectanceToSurfaceReflectanceImageFilter.
+ *
+ * \sa AtmosphericRadiativeTerms
+ * \sa AtmosphericCorrectionParameters
+ * \sa ReflectanceToSurfaceReflectanceImageFilter
+ * \ingroup DataSources
+ * \ingroup Radiometry
+ * \deprecated This class has been replaced by
+ *  otb::RadiometryCorrectionParametersToAtmosphericRadiativeTerms,
+ *  which only contains a static function that computes the radiative terms.
+ */
+class ITK_EXPORT AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms
+  : public itk::ProcessObject
+{
+public:
+  /** Standard typedefs */
+  typedef AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms Self;
+  typedef itk::ProcessObject                                           Superclass;
+  typedef itk::SmartPointer<Self>                                      Pointer;
+  typedef itk::SmartPointer<const Self>                                ConstPointer;
+
+  /** Runtime information */
+  itkTypeMacro(AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms, itk::ProcessObject);
+  /** Creation through the object factory */
+  itkNewMacro(Self);
+  /** Template parameters typedefs */
+  typedef AtmosphericCorrectionParameters              AtmosphericCorrectionParametersType;
+  typedef AtmosphericCorrectionParametersType::Pointer AtmosphericCorrectionParametersPointer;
+  typedef AtmosphericRadiativeTerms                    AtmosphericRadiativeTermsType;
+  typedef AtmosphericRadiativeTermsType::Pointer       AtmosphericRadiativeTermsPointer;
+
+  /** Set the Atmospheric Correction Parameters input of this process object */
+  void SetInput(const AtmosphericCorrectionParametersType *object);
+
+  /** Get the Atmospheric Correction Parameters input of this process object */
+  AtmosphericCorrectionParametersType * GetInput(void);
+
+  DataObjectPointer MakeOutput(unsigned int);
+  void GraftOutput(itk::DataObject *graft);
+  void GraftNthOutput(unsigned int idx, itk::DataObject *graft);
+
+  /** Get the Atmospheric Radiative Terms output of this process object.  */
+  virtual AtmosphericRadiativeTermsType * GetOutput(void);
+  virtual AtmosphericRadiativeTermsType * GetOutput(unsigned int idx);
+
+  /** Generate the output.*/
+  virtual void GenerateData();
+
+protected:
+  /** Constructor */
+  AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms();
+  /** Destructor */
+  virtual ~AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms() {}
+  /** PrintSelf method */
+  virtual void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+private:
+  AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms(const Self &); //purposely not implemented
+  void operator =(const Self&); //purposely not implemented
+
+};
+
+} // end namespace otb
+
+#endif

Code/Radiometry/otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.h

+/*=========================================================================
+
+ Program:   ORFEO Toolbox
+ Language:  C++
+ Date:      $Date$
+ Version:   $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ Some parts of this code are derived from ITK. See ITKCopyright.txt
+ for details.
+
+
+    This software is distributed WITHOUT ANY WARRANTY; without even
+    the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+    PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter_h
+#define __otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter_h
+
+#include "itkNumericTraits.h"
+#include <vector>
+#include "itkConstNeighborhoodIterator.h"
+#include "otbUnaryFunctorNeighborhoodImageFilter.h"
+#include "itkVariableSizeMatrix.h"
+#include "otbAtmosphericRadiativeTerms.h"
+#include "otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.h"
+#include <iomanip>
+
+namespace otb
+{
+namespace Functor
+{
+/** \class ComputeNeighborhoodContributionFunctor
+*  \brief Unary neighborhood functor to compute the value of a pixel which is a sum
+*   of the surrounding pixels value ponderated by a coefficient.
+*
+*  \ingroup Functor
+* \ingroup Radiometry
+*/
+template <class TNeighIter, class TOutput>
+class ComputeNeighborhoodContributionFunctor
+{
+public:
+  ComputeNeighborhoodContributionFunctor() {}
+  virtual ~ComputeNeighborhoodContributionFunctor() {}
+
+  typedef itk::VariableSizeMatrix<double>           WeightingMatrixType;
+  typedef typename std::vector<WeightingMatrixType> WeightingValuesContainerType;
+  typedef typename TOutput::RealValueType           RealValueType;
+  typedef std::vector<double>                       DoubleContainerType;
+
+  void SetWeightingValues(const WeightingValuesContainerType& cont)
+  {
+    m_WeightingValues = cont;
+  }
+  void SetUpwardTransmittanceRatio(DoubleContainerType upwardTransmittanceRatio)
+  {
+    m_UpwardTransmittanceRatio = upwardTransmittanceRatio;
+  }
+  void SetDiffuseRatio(DoubleContainerType diffuseRatio)
+  {
+    m_DiffuseRatio = diffuseRatio;
+  }
+  WeightingValuesContainerType GetWeightingValues()
+  {
+    return m_WeightingValues;
+  }
+  DoubleContainerType GetUpwardTransmittanceRatio()
+  {
+    return m_UpwardTransmittanceRatio;
+  }
+  DoubleContainerType GetDiffuseRatio()
+  {
+    return m_DiffuseRatio;
+  }
+
+  inline TOutput operator ()(const TNeighIter& it)
+  {
+    unsigned int neighborhoodSize = it.Size();
+    double       contribution = 0.;
+    TOutput      outPixel;
+    outPixel.SetSize(it.GetCenterPixel().Size());
+
+    // Loop over each component
+    const unsigned int size = outPixel.GetSize();
+    for (unsigned int j = 0; j < size; ++j)
+    {
+      contribution = 0;
+      // Load the current channel ponderation value matrix
+      WeightingMatrixType TempChannelWeighting = m_WeightingValues[j];
+      // Loop over the neighborhood
+      for (unsigned int i = 0; i < neighborhoodSize; ++i)
+      {
+        // Current neighborhood pixel index calculation
+        unsigned int RowIdx = 0;
+        unsigned int ColIdx = 0;
+        RowIdx = i / TempChannelWeighting.Cols();
+        ColIdx = i - RowIdx*TempChannelWeighting.Cols();
+   
+        // Extract the current neighborhood pixel ponderation
+        double idVal = TempChannelWeighting(RowIdx, ColIdx);
+        // Extract the current neighborhood pixel value
+        TOutput tempPix = it.GetPixel(i);
+
+        contribution += static_cast<double>(tempPix[j]) * idVal;
+
+        }
+        
+        outPixel[j] = static_cast<RealValueType>(it.GetCenterPixel()[j]) *m_UpwardTransmittanceRatio[j]
+                        + contribution * m_DiffuseRatio[j];
+      }
+    return outPixel;
+  }
+
+private:
+  WeightingValuesContainerType m_WeightingValues;
+  DoubleContainerType          m_UpwardTransmittanceRatio;
+  DoubleContainerType          m_DiffuseRatio;
+};
+
+}
+
+/** \class SurfaceAdjacencyEffect6SCorrectionSchemeFilter
+ *  \brief Correct the scheme taking care of the surrounding pixels.
+ *
+ *   The SurfaceAdjacencyEffect6SCorrectionSchemeFilter class allows to introduce a neighbor correction to the
+ *   reflectance estimation. The satelite signal is considered as to be a combinaison of the signal coming from
+ *   the target pixel and a weighting of the siganls coming from the neighbor pixels.
+ *
+ * \ingroup Radiometry
+ * \deprecated This class has been replaced by otb::SurfaceAdjacencyEffectCorrectionSchemeFilter.
+ */
+template <class TInputImage, class TOutputImage>
+class ITK_EXPORT SurfaceAdjacencyEffect6SCorrectionSchemeFilter :
+  public UnaryFunctorNeighborhoodImageFilter<
+      TInputImage,
+      TOutputImage,
+      typename Functor::ComputeNeighborhoodContributionFunctor<itk::
+          ConstNeighborhoodIterator <TInputImage>,
+      typename TOutputImage::PixelType> >
+{
+public:
+  /** "typedef" to simplify the variables definition and the declaration. */
+  typedef Functor::ComputeNeighborhoodContributionFunctor<itk::ConstNeighborhoodIterator<TInputImage>,
+      typename TOutputImage::PixelType> FunctorType;
+
+  /** "typedef" for standard classes. */
+  typedef SurfaceAdjacencyEffect6SCorrectionSchemeFilter                              Self;
+  typedef UnaryFunctorNeighborhoodImageFilter<TInputImage, TOutputImage, FunctorType> Superclass;
+  typedef itk::SmartPointer<Self>                                                     Pointer;
+  typedef itk::SmartPointer<const Self>                                               ConstPointer;
+
+  typedef typename Superclass::InputImageType  InputImageType;
+  typedef typename Superclass::OutputImageType OutputImageType;
+
+  typedef std::vector<double> DoubleContainerType;
+  /** object factory method. */
+  itkNewMacro(Self);
+
+  /** return class name. */
+  itkTypeMacro(SurfaceAdjacencyEffect6SCorrectionSchemeFilter, UnaryFunctorNeighborhoodImageFilter);
+
+  /**   Extract input and output images dimensions.*/
+  itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
+  itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension);
+
+  /** Supported images definition. */
+  typedef typename InputImageType::PixelType          InputPixelType;
+  typedef typename InputImageType::InternalPixelType  InputInternalPixelType;
+  typedef typename InputImageType::RegionType         InputImageRegionType;
+  typedef typename InputImageType::SizeType           SizeType;
+  typedef typename OutputImageType::PixelType         OutputPixelType;
+  typedef typename OutputImageType::InternalPixelType OutputInternalPixelType;
+  typedef typename OutputImageType::RegionType        OutputImageRegionType;
+
+  typedef AtmosphericRadiativeTerms                                    AtmosphericRadiativeTermsType;
+  typedef typename AtmosphericRadiativeTermsType::Pointer              AtmosphericRadiativeTermsPointerType;
+  typedef AtmosphericCorrectionParameters::Pointer                     CorrectionParametersPointerType;
+  typedef AtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms Parameters2RadiativeTermsType;
+  typedef Parameters2RadiativeTermsType::Pointer                       Parameters2RadiativeTermsPointerType;
+
+  typedef FilterFunctionValues                       FilterFunctionValuesType;
+  typedef FilterFunctionValuesType::ValuesVectorType CoefVectorType;
+  typedef std::vector<CoefVectorType>                FilterFunctionCoefVectorType;
+
+  typedef itk::MetaDataDictionary MetaDataDictionaryType;
+
+  /** Storage ponderation values types*/
+  typedef itk::VariableSizeMatrix<double>           WeightingMatrixType;
+  typedef typename std::vector<WeightingMatrixType> WeightingValuesContainerType;
+
+  /** typedef for calculation*/
+  typedef typename itk::ConstNeighborhoodIterator<InputImageType> NeighborIterType;
+
+  /** Set/Get the Size of the neighbor window. */
+  void SetWindowRadius(unsigned int rad)
+  {
+    this->SetRadius(rad);
+    m_WindowRadius = rad;
+    this->Modified();
+  }
+  itkGetConstReferenceMacro(WindowRadius, unsigned int);
+
+  /** Set/Get the pixel spacing in kilometers */
+  itkSetMacro(PixelSpacingInKilometers, double);
+  itkGetMacro(PixelSpacingInKilometers, double);
+  /** Set/Get the viewing angle */
+  itkSetMacro(ZenithalViewingAngle, double);
+  itkGetMacro(ZenithalViewingAngle, double);
+
+  /** Get/Set Atmospheric Radiative Terms. */
+  void SetAtmosphericRadiativeTerms(AtmosphericRadiativeTermsPointerType atmo)
+  {
+    m_AtmosphericRadiativeTerms = atmo;
+    this->SetNthInput(1, m_AtmosphericRadiativeTerms);
+    m_IsSetAtmosphericRadiativeTerms = true;
+    this->Modified();
+  }
+  itkGetObjectMacro(AtmosphericRadiativeTerms, AtmosphericRadiativeTerms);
+
+  /** Get/Set Atmospheric Correction Parameters. */
+  itkSetObjectMacro(CorrectionParameters, AtmosphericCorrectionParameters);
+  itkGetObjectMacro(CorrectionParameters, AtmosphericCorrectionParameters);
+
+  /** Get/Set Aeronet file name. */
+  itkSetMacro(AeronetFileName, std::string);
+  itkGetMacro(AeronetFileName, std::string);
+
+  /** Get/Set Aeronet file name. */
+  itkSetMacro(FilterFunctionValuesFileName, std::string);
+  itkGetMacro(FilterFunctionValuesFileName, std::string);
+
+  /** Get/Set Filter function coef. */
+  void SetFilterFunctionCoef(FilterFunctionCoefVectorType vect)
+  {
+    m_FilterFunctionCoef = vect;
+    this->Modified();
+  }
+  FilterFunctionCoefVectorType GetFilterFunctionCoef()
+  {
+    return m_FilterFunctionCoef;
+  }
+
+  /** Compute the functor parameters */
+  void ComputeParameters();
+  /** Compute radiative terms if necessary and then update functors attibuts. */
+  void GenerateParameters();
+  /** Fill AtmosphericRadiativeTerms using image metadata*/
+  void UpdateAtmosphericRadiativeTerms();
+
+protected:
+  SurfaceAdjacencyEffect6SCorrectionSchemeFilter();
+  virtual ~SurfaceAdjacencyEffect6SCorrectionSchemeFilter() {}
+  void PrintSelf(std::ostream& os, itk::Indent indent) const;
+
+  /** Initialize the parameters of the functor before the threads run. */
+  virtual void BeforeThreadedGenerateData();
+
+  /** If modified, we need to compute the parameters again */
+  virtual void Modified();
+
+private:
+  /** Size of the window. */
+  unsigned int m_WindowRadius;
+  /** Weighting values for the neighbor pixels.*/
+  WeightingValuesContainerType m_WeightingValues;
+  /** True if the parameters have been generated */
+  bool m_ParametersHaveBeenComputed;
+  /** Pixel spacing in kilometers */
+  double m_PixelSpacingInKilometers;
+  /** Viewing angle in degree */
+  double m_ZenithalViewingAngle;
+  /** Radiative terms object */
+  AtmosphericRadiativeTermsPointerType m_AtmosphericRadiativeTerms;
+  bool                                 m_IsSetAtmosphericRadiativeTerms;
+  /** Atmospheric Correction Parameters. */
+  CorrectionParametersPointerType m_CorrectionParameters;
+  /** Path to an Aeronet data file, allows to compute aerosol optical and water vapor amounts. */
+  std::string m_AeronetFileName;
+  /** Path to an filter function values file. */
+  std::string m_FilterFunctionValuesFileName;
+  /** Contains the filter function values (each element is a vector and reprsents the values for each channel) */
+  FilterFunctionCoefVectorType m_FilterFunctionCoef;
+};
+
+} // end namespace otb
+
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.txx"
+#endif
+
+#endif

Code/Radiometry/otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.txx

+/*=========================================================================
+
+  Program:   ORFEO Toolbox
+  Language:  C++
+  Date:      $Date$
+  Version:   $Revision$
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+  See OTBCopyright.txt for details.
+
+
+  Some parts of this code are derived from ITK. See ITKCopyright.txt
+  for details.
+
+
+     This software is distributed WITHOUT ANY WARRANTY; without even
+     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+     PURPOSE.  See the above copyright notices for more information.
+
+=========================================================================*/
+
+#ifndef __otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter_txx
+#define __otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter_txx
+
+#include "otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.h"
+
+#include "itkConstNeighborhoodIterator.h"
+#include "itkNeighborhoodInnerProduct.h"
+#include "itkImageRegionIterator.h"
+#include "itkNeighborhoodAlgorithm.h"
+#include "itkZeroFluxNeumannBoundaryCondition.h"
+#include "itkOffset.h"
+#include "itkProgressReporter.h"
+#include "otbImage.h"
+#include "otbSIXSTraits.h"
+#include "otbMath.h"
+#include "otbOpticalImageMetadataInterfaceFactory.h"
+#include "otbOpticalImageMetadataInterface.h"
+
+namespace otb
+{
+
+template <class TInputImage, class TOutputImage>
+SurfaceAdjacencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
+::SurfaceAdjacencyEffect6SCorrectionSchemeFilter()
+{
+  m_WindowRadius = 1;
+  m_ParametersHaveBeenComputed=false;
+  m_PixelSpacingInKilometers = 1.;
+  m_ZenithalViewingAngle = 361.;
+  m_AtmosphericRadiativeTerms = AtmosphericRadiativeTermsType::New();
+  m_CorrectionParameters      = AtmosphericCorrectionParameters::New();
+  m_IsSetAtmosphericRadiativeTerms = false;
+  m_AeronetFileName = "";
+  m_FilterFunctionValuesFileName = "";
+  m_FilterFunctionCoef.clear();
+}
+
+template <class TInputImage, class TOutputImage>
+void
+SurfaceAdjacencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
+::BeforeThreadedGenerateData()
+{
+  Superclass::BeforeThreadedGenerateData();
+
+  this->GenerateParameters();
+
+  if (!m_ParametersHaveBeenComputed)
+    {
+    this->ComputeParameters();
+    m_ParametersHaveBeenComputed = true;
+    }
+}
+
+template <class TInputImage, class TOutputImage>
+void
+SurfaceAdjacencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
+::Modified()
+{
+  Superclass::Modified();
+  m_ParametersHaveBeenComputed  = false;
+}
+
+template <class TInputImage, class TOutputImage>
+void
+SurfaceAdjacencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
+::UpdateAtmosphericRadiativeTerms()
+{
+  MetaDataDictionaryType dict = this->GetInput()->GetMetaDataDictionary();
+
+  OpticalImageMetadataInterface::Pointer imageMetadataInterface = OpticalImageMetadataInterfaceFactory::CreateIMI(dict);
+
+  if ((m_CorrectionParameters->GetDay() == 0))
+    {
+    m_CorrectionParameters->SetDay(imageMetadataInterface->GetDay());
+    }
+
+  if ((m_CorrectionParameters->GetMonth() == 0))
+    {
+    m_CorrectionParameters->SetMonth(imageMetadataInterface->GetMonth());
+    }
+
+  if ((m_CorrectionParameters->GetSolarZenithalAngle() == 361.))
+    {
+    m_CorrectionParameters->SetSolarZenithalAngle(90. - imageMetadataInterface->GetSunElevation());
+    }
+
+  if ((m_CorrectionParameters->GetSolarAzimutalAngle() == 361.))
+    {
+    m_CorrectionParameters->SetSolarAzimutalAngle(imageMetadataInterface->GetSunAzimuth());
+    }
+
+  if ((m_CorrectionParameters->GetViewingZenithalAngle() == 361.))
+    {
+    m_CorrectionParameters->SetViewingZenithalAngle(90. - imageMetadataInterface->GetSatElevation());
+    }
+
+  if (m_ZenithalViewingAngle == 361.)
+    {
+    this->SetZenithalViewingAngle(90. - imageMetadataInterface->GetSatElevation());
+    }
+
+  if ((m_CorrectionParameters->GetViewingAzimutalAngle() == 361.))
+    {
+    m_CorrectionParameters->SetViewingAzimutalAngle(imageMetadataInterface->GetSatAzimuth());
+    }
+
+  if (m_AeronetFileName != "")
+    {
+    m_CorrectionParameters->UpdateAeronetData(m_AeronetFileName,
+                                              imageMetadataInterface->GetYear(),
+                                              imageMetadataInterface->GetHour(),
+                                              imageMetadataInterface->GetMinute());
+    }
+
+  // load filter function values
+  if (m_FilterFunctionValuesFileName != "")
+    {
+    m_CorrectionParameters->LoadFilterFunctionValue(m_FilterFunctionValuesFileName);
+    }
+  // the user has set the filter function values
+  else
+    {
+    if (m_FilterFunctionCoef.size() != this->GetInput()->GetNumberOfComponentsPerPixel())
+      {
+      itkExceptionMacro(<< "Filter Function and image channels mismatch.");
+      }
+    for (unsigned int i = 0; i < this->GetInput()->GetNumberOfComponentsPerPixel(); ++i)
+      {
+      FilterFunctionValuesType::Pointer functionValues = FilterFunctionValuesType::New();
+      functionValues->SetFilterFunctionValues(m_FilterFunctionCoef[i]);
+      functionValues->SetMinSpectralValue(imageMetadataInterface->GetFirstWavelengths()[i]);
+      functionValues->SetMaxSpectralValue(imageMetadataInterface->GetLastWavelengths()[i]);
+
+      m_CorrectionParameters->SetWavelengthSpectralBandWithIndex(i, functionValues);
+      }
+    }
+
+  Parameters2RadiativeTermsPointerType param2Terms = Parameters2RadiativeTermsType::New();
+
+  param2Terms->SetInput(m_CorrectionParameters);
+  param2Terms->Update();
+
+  m_AtmosphericRadiativeTerms = param2Terms->GetOutput();
+}
+
+template <class TInputImage, class TOutputImage>
+void
+SurfaceAdjacencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
+::GenerateParameters()
+{
+  if (m_IsSetAtmosphericRadiativeTerms == false)
+    {
+    this->UpdateAtmosphericRadiativeTerms();
+    }
+}
+
+template <class TInputImage, class TOutputImage>
+void
+SurfaceAdjacencyEffect6SCorrectionSchemeFilter<TInputImage, TOutputImage>
+::ComputeParameters()
+{
+  // get pointers to the input and output
+  typename InputImageType::Pointer  inputPtr = const_cast<TInputImage *>(this->GetInput());
+  typename OutputImageType::Pointer outputPtr = const_cast<TOutputImage *>(this->GetOutput());
+
+  WeightingMatrixType radiusMatrix(2*m_WindowRadius + 1, 2*m_WindowRadius + 1);
+  radiusMatrix.Fill(0.);
+
+  double center = static_cast<double>(m_WindowRadius);
+
+  for (unsigned int i = 0; i < m_WindowRadius + 1; ++i)
+    {
+    for (unsigned int j = 0; j < m_WindowRadius + 1; ++j)
+      {
+      double id = static_cast<double>(i);
+      double jd = static_cast<double>(j);
+      double currentRadius = m_PixelSpacingInKilometers * vcl_sqrt(vcl_pow(id - center, 2) + vcl_pow(jd - center, 2));
+      radiusMatrix(i, j) = currentRadius;
+      radiusMatrix(2 * m_WindowRadius - i, j) = currentRadius;
+      radiusMatrix(2 * m_WindowRadius - i, 2 * m_WindowRadius - j) = currentRadius;
+      radiusMatrix(i, 2 * m_WindowRadius - j) = currentRadius;
+      }
+    }
+
+  for (unsigned int band = 0; band < inputPtr->GetNumberOfComponentsPerPixel(); ++band)
+    {
+    double rayleigh = m_AtmosphericRadiativeTerms->GetUpwardDiffuseTransmittanceForRayleigh(band);
+    double aerosol =  m_AtmosphericRadiativeTerms->GetUpwardDiffuseTransmittanceForAerosol(band);
+    
+    WeightingMatrixType currentWeightingMatrix(2*m_WindowRadius + 1, 2*m_WindowRadius + 1);
+    currentWeightingMatrix.Fill(0.);
+
+    for (unsigned int i = 0; i < 2 * m_WindowRadius + 1; ++i)
+      {
+      for (unsigned int j = 0; j < 2 * m_WindowRadius + 1; ++j)
+        {
+        double notUsed1, notUsed2;
+        double factor = 1;
+        double palt = 1000.;
+        SIXSTraits::ComputeEnvironmentalContribution(rayleigh,
+                                                     aerosol,
+                                                     radiusMatrix(i,j),
+                                                     palt,
+                                                     vcl_cos(m_ZenithalViewingAngle * CONST_PI_180),
+                                                     notUsed1,
+                                                     notUsed2,
+                                                     factor);                                                                                                        //Call to 6S
+        currentWeightingMatrix(i, j) = factor;
+        }
+      }
+    m_WeightingValues.push_back(currentWeightingMatrix);
+    }
+
+  DoubleContainerType upwardTransmittanceRatio, diffuseRatio;
+
+  for (unsigned int band = 0; band < inputPtr->GetNumberOfComponentsPerPixel(); ++band)
+    {
+    upwardTransmittanceRatio.push_back(m_AtmosphericRadiativeTerms->GetUpwardTransmittance(
+                                         band) / m_AtmosphericRadiativeTerms->GetUpwardDirectTransmittance(band));
+    diffuseRatio.push_back(m_AtmosphericRadiativeTerms->GetUpwardDiffuseTransmittance(
+                             band) / m_AtmosphericRadiativeTerms->GetUpwardDirectTransmittance(band));
+    }
+  this->GetFunctor().SetUpwardTransmittanceRatio(upwardTransmittanceRatio);
+  this->GetFunctor().SetDiffuseRatio(diffuseRatio);
+  this->GetFunctor().SetWeightingValues(m_WeightingValues);
+}
+/**
+ * Standard "PrintSelf" method
+ */
+template <class TInputImage, class TOutput>
+void
+SurfaceAdjacencyEffect6SCorrectionSchemeFilter<TInputImage, TOutput>
+::PrintSelf(std::ostream& os, itk::Indent indent) const
+{
+  os << indent << "Radius : " << m_WindowRadius << std::endl;
+  os << indent << "Pixel spacing in kilometers: " << m_PixelSpacingInKilometers << std::endl;
+  os << indent << "Zenithal viewing angle in degree: " << m_ZenithalViewingAngle << std::endl;
+}
+
+} // end namespace otb
+
+#endif