diff --git a/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.cxx b/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0cef1e2a1a349f78f18ad1b50e287aae128ecd80
--- /dev/null
+++ b/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.cxx
@@ -0,0 +1,216 @@
+/*=========================================================================
+
+ 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
diff --git a/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.h b/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.h
new file mode 100644
index 0000000000000000000000000000000000000000..f37f35ac8a77095d24121a0d59b77746b9ef86be
--- /dev/null
+++ b/Code/Radiometry/otbAtmosphericCorrectionParametersTo6SAtmosphericRadiativeTerms.h
@@ -0,0 +1,97 @@
+/*=========================================================================
+
+ 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
diff --git a/Code/Radiometry/otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.h b/Code/Radiometry/otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..d7f1b75c71b9f2143dfe08cce52bc4023846efcb
--- /dev/null
+++ b/Code/Radiometry/otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.h
@@ -0,0 +1,297 @@
+/*=========================================================================
+
+ 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
diff --git a/Code/Radiometry/otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.txx b/Code/Radiometry/otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..b050c6d8a484df6da78010b610c219a5e8ef84a6
--- /dev/null
+++ b/Code/Radiometry/otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.txx
@@ -0,0 +1,264 @@
+/*=========================================================================
+
+ 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