otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.h

Go to the documentation of this file.

/*=========================================================================

 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
{
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());
    std::ostringstream oss;
    oss << std::fixed << std::setprecision(3);

    // Loop over each component
    const unsigned int size = outPixel.GetSize();
    for (unsigned int j = 0; j < size; ++j)
      {
      oss.str("");
      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();
        if (RowIdx != 0)
          {
          ColIdx = (i + 1) - RowIdx*TempChannelWeighting.Cols() - 1;
          }
        else
          {
          ColIdx = i;
          }
        // 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;

        }
      double temp = static_cast<double>(it.GetCenterPixel()[j]) * m_UpwardTransmittanceRatio[j] + contribution *
                    m_DiffuseRatio[j];
      oss << temp;
      outPixel[j] = static_cast<RealValueType>(atof(oss.str().c_str()));

      //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;
};

}

template <class TInputImage, class TOutputImage>
class ITK_EXPORT SurfaceAdjacencyEffect6SCorrectionSchemeFilter :
  public UnaryFunctorNeighborhoodImageFilter<TInputImage,
      TOutputImage,
      typename Functor::ComputeNeighborhoodContributionFunctor<itk::
          ConstNeighborhoodIterator
          <TInputImage>,
          typename
          TOutputImage
          ::PixelType> >
{
public:
  typedef Functor::ComputeNeighborhoodContributionFunctor<itk::ConstNeighborhoodIterator<TInputImage>,
      typename TOutputImage::PixelType> FunctorType;

  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;
  itkNewMacro(Self);

  itkTypeMacro(SurfaceAdjacencyEffect6SCorrectionSchemeFilter, UnaryFunctorNeighborhoodImageFilter);

  itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
  itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension);

  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;

  typedef itk::VariableSizeMatrix<double>           WeightingMatrixType;
  typedef typename std::vector<WeightingMatrixType> WeightingValuesContainerType;

  typedef typename itk::ConstNeighborhoodIterator<InputImageType> NeighborIterType;

  void SetWindowRadius(unsigned int rad)
  {
    this->SetRadius(rad);
    m_WindowRadius = rad;
    this->Modified();
  }
  itkGetConstReferenceMacro(WindowRadius, unsigned int);

  itkSetMacro(PixelSpacingInKilometers, double);
  itkGetMacro(PixelSpacingInKilometers, double);
  itkSetMacro(ZenithalViewingAngle, double);
  itkGetMacro(ZenithalViewingAngle, double);

  void SetAtmosphericRadiativeTerms(AtmosphericRadiativeTermsPointerType atmo)
  {
    m_AtmosphericRadiativeTerms = atmo;
    this->SetNthInput(1, m_AtmosphericRadiativeTerms);
    m_IsSetAtmosphericRadiativeTerms = true;
    this->Modified();
  }
  itkGetObjectMacro(AtmosphericRadiativeTerms, AtmosphericRadiativeTerms);

  itkSetObjectMacro(CorrectionParameters, AtmosphericCorrectionParameters);
  itkGetObjectMacro(CorrectionParameters, AtmosphericCorrectionParameters);

  itkSetMacro(AeronetFileName, std::string);
  itkGetMacro(AeronetFileName, std::string);

  itkSetMacro(FilterFunctionValuesFileName, std::string);
  itkGetMacro(FilterFunctionValuesFileName, std::string);

  void SetFilterFunctionCoef(FilterFunctionCoefVectorType vect)
  {
    m_FilterFunctionCoef = vect;
    this->Modified();
  }
  FilterFunctionCoefVectorType GetFilterFunctionCoef()
  {
    return m_FilterFunctionCoef;
  }

  void ComputeParameters();
  void GenerateParameters();
  void UpdateAtmosphericRadiativeTerms();

protected:
  SurfaceAdjacencyEffect6SCorrectionSchemeFilter();
  virtual ~SurfaceAdjacencyEffect6SCorrectionSchemeFilter() {}
  void PrintSelf(std::ostream& os, itk::Indent indent) const;

  virtual void GenerateOutputInformation();

  //virtual void BeforeThreadedGenerateData();

  virtual void Modified();

private:
  unsigned int m_WindowRadius;
  WeightingValuesContainerType m_WeightingValues;
  bool m_ParametersHaveBeenComputed;
  double m_PixelSpacingInKilometers;
  double m_ZenithalViewingAngle;
  AtmosphericRadiativeTermsPointerType m_AtmosphericRadiativeTerms;
  bool                                 m_IsSetAtmosphericRadiativeTerms;
  CorrectionParametersPointerType m_CorrectionParameters;
  std::string m_AeronetFileName;
  std::string m_FilterFunctionValuesFileName;
  FilterFunctionCoefVectorType m_FilterFunctionCoef;
  bool m_UseGenerateParameters;
};

} // end namespace otb

#ifndef OTB_MANUAL_INSTANTIATION
#include "otbSurfaceAdjacencyEffect6SCorrectionSchemeFilter.txx"
#endif

#endif