otbBayesianFusionFilter.h

Go to the documentation of this file.

/*
 * Copyright (C) 2007 Julien Radoux
 *
 * 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 otbBayesianFusionFilter_h
#define otbBayesianFusionFilter_h
 
#include "itkImageToImageFilter.h"
#include "otbVectorImage.h"
#include "otbStreamingStatisticsVectorImageFilter.h"
#include "otbFunctorImageFilter.h"
#include "otbStreamingMatrixTransposeMatrixImageFilter.h"
#include "otbImageToVectorImageCastFilter.h"
 
namespace otb
{
 
namespace Functor
{
template <class TInputMultiSpectral, class TInputMultiSpectralInterp, class TInputPanchro, class TOutput>
class BayesianFunctor
{
public:
  BayesianFunctor()
  {
  }
  virtual ~BayesianFunctor()
  {
  }
  typedef typename TInputMultiSpectral::RealValueType RealType;
  typedef typename itk::VariableSizeMatrix<RealType>  MatrixType;
 
  void SetLambda(float lambda)
  {
    m_Lambda = lambda;
  }
  void SetS(float S)
  {
    m_S = S;
  }
  void SetAlpha(float alpha)
  {
    m_Alpha = alpha;
  }
  void SetBeta(MatrixType matrix)
  {
    m_Beta = matrix;
  }
  void SetCovarianceInvMatrix(MatrixType matrix)
  {
    m_CovarianceInvMatrix = matrix;
  }
  void SetVcondopt(MatrixType matrix)
  {
    m_Vcondopt = matrix;
  }
  float GetLambda()
  {
    return m_Lambda;
  }
  float GetAlpha()
  {
    return m_Alpha;
  }
  float GetS()
  {
    return m_S;
  }
  MatrixType GetBeta()
  {
    return m_Beta;
  }
  MatrixType GetCovarianceInvMatrix()
  {
    return m_CovarianceInvMatrix;
  }
  MatrixType GetVcondopt()
  {
    return m_Vcondopt;
  }
 
  void operator()(TOutput& obs, const TInputMultiSpectral& itkNotUsed(ms), const TInputMultiSpectralInterp& msi, const TInputPanchro& p)
  {
    MatrixType obsMat, msiVect;
    obsMat.SetSize(1, obs.GetSize());
    msiVect.SetSize(1, msi.GetSize());
    for (unsigned int i = 0; i < msi.GetSize(); ++i)
    {
      msiVect(0, i) = msi[i];
    }
    obsMat = msiVect * m_CovarianceInvMatrix;
    obsMat *= 2 * (1 - m_Lambda);
    MatrixType PanVect;
    PanVect = m_Beta.GetTranspose();
    PanVect *= (p - m_Alpha);
    PanVect /= m_S;
    PanVect *= 2 * m_Lambda;
 
    if ((obsMat.Cols() != PanVect.Cols()) || (obsMat.Rows() != PanVect.Rows()))
    {
      itkGenericExceptionMacro(<< "Matrix with size (" << obsMat.Rows() << "," << obsMat.Cols() << ") cannot be subtracted from matrix with size ("
                               << PanVect.Rows() << "," << PanVect.Cols() << " )");
    }
 
    for (unsigned int r = 0; r < obsMat.Rows(); ++r)
    {
      for (unsigned int c = 0; c < obsMat.Cols(); ++c)
      {
        obsMat(r, c) += PanVect(r, c);
      }
    }
    obsMat *= m_Vcondopt;
    for (unsigned int i = 0; i < obs.GetSize(); ++i)
    {
      obs[i] = static_cast<typename TOutput::ValueType>(obsMat(0U, i));
    }
  }
 
  constexpr size_t OutputSize(const std::array<size_t, 3> inputsNbBands) const
  {
    return inputsNbBands[1];
  }
 
private:
  float      m_Lambda;
  float      m_S;
  float      m_Alpha;
  MatrixType m_CovarianceInvMatrix;
  MatrixType m_Beta;
  MatrixType m_Vcondopt;
};
}
 
/***** TODO ***
 * Complete the description with J. Radoux text
 */
 
/***** END TODO ***/
 
template <class TInputMultiSpectralImage, class TInputMultiSpectralInterpImage, class TInputPanchroImage, class TOutputImage>
class ITK_EXPORT BayesianFusionFilter
    : public FunctorImageFilter<Functor::BayesianFunctor<typename TInputMultiSpectralImage::PixelType, typename TInputMultiSpectralInterpImage::PixelType,
                                                         typename TInputPanchroImage::PixelType, typename TOutputImage::PixelType>>
{
public:
  itkStaticConstMacro(InputImageDimension, unsigned int, TInputMultiSpectralImage::ImageDimension);
  itkStaticConstMacro(OutputImageDimension, unsigned int, TOutputImage::ImageDimension);
 
  typedef TInputMultiSpectralImage       InputMultiSpectralImageType;
  typedef TInputMultiSpectralInterpImage InputMultiSpectralInterpImageType;
  typedef TInputPanchroImage             InputPanchroImageType;
  typedef TOutputImage                   OutputImageType;
 
  typedef BayesianFusionFilter Self;
  using BayesianFunctorType = Functor::BayesianFunctor<typename TInputMultiSpectralImage::PixelType, typename TInputMultiSpectralInterpImage::PixelType,
                                                       typename TInputPanchroImage::PixelType, typename TOutputImage::PixelType>;
  typedef FunctorImageFilter<BayesianFunctorType> Superclass;
  typedef itk::SmartPointer<Self>                 Pointer;
  typedef itk::SmartPointer<const Self>           ConstPointer;
 
  itkNewMacro(Self);
 
  itkTypeMacro(BayesianFusionFilter, FunctorImageFilter);
 
  typedef typename InputMultiSpectralImageType::PixelType               InputMultiSpectralPixelType;
  typedef typename InputMultiSpectralImageType::InternalPixelType       InputMultiSpectralInternalPixelType;
  typedef typename InputMultiSpectralInterpImageType::PixelType         InputMultiSpectralInterpPixelType;
  typedef typename InputMultiSpectralInterpImageType::InternalPixelType InputMultiSpectralInterpInternalPixelType;
  typedef typename InputPanchroImageType::PixelType                     InputPanchroPixelType;
  typedef typename OutputImageType::PixelType                           OutputPixelType;
  typedef typename OutputImageType::InternalPixelType                   OutputInternalPixelType;
 
  typedef typename itk::NumericTraits<InputPanchroPixelType>::RealType                     InputPanchroRealType;
  typedef typename itk::NumericTraits<InputMultiSpectralInternalPixelType>::RealType       InputMultiSpectralRealType;
  typedef typename InputMultiSpectralImageType::RegionType                                 InputMultiSpectralImageRegionType;
  typedef typename itk::NumericTraits<InputMultiSpectralInterpInternalPixelType>::RealType InputMultiSpectralInterpRealType;
  typedef typename InputMultiSpectralInterpImageType::RegionType                           InputMultiSpectralInterpImageRegionType;
  typedef typename InputPanchroImageType::RegionType                                       InputPanchroImageRegionType;
  typedef typename OutputImageType::RegionType                                             OutputImageRegionType;
 
  typedef typename InputMultiSpectralImageType::SizeType SizeType;
 
  typedef StreamingStatisticsVectorImageFilter<InputMultiSpectralInterpImageType> StreamingStatisticsVectorImageFilterType;
  typedef typename StreamingStatisticsVectorImageFilterType::MatrixType           MatrixType;
  typedef StreamingMatrixTransposeMatrixImageFilter<InputMultiSpectralImageType, InputMultiSpectralImageType> MSTransposeMSType;
  typedef ImageToVectorImageCastFilter<InputPanchroImageType, InputMultiSpectralImageType>                    CasterType;
 
  void SetMultiSpect(const InputMultiSpectralImageType* multiSpect)
  {
    this->template SetInput<0>(multiSpect);
  }
 
  void SetMultiSpectInterp(const InputMultiSpectralInterpImageType* multiSpectInterp)
  {
    this->template SetInput<1>(multiSpectInterp);
  }
 
  void SetPanchro(const InputPanchroImageType* panchro)
  {
    this->template SetInput<2>(panchro);
  }
 
  const InputMultiSpectralImageType* GetMultiSpect()
  {
    return this->template GetInput<0>();
  }
 
  const InputMultiSpectralInterpImageType* GetMultiSpectInterp()
  {
    return this->template GetInput<1>();
  }
 
  const InputPanchroImageType* GetPanchro()
  {
    return this->template GetInput<2>();
  }
 
  itkSetMacro(Lambda, float);
 
  itkGetConstReferenceMacro(Lambda, float);
 
  itkSetMacro(Beta, MatrixType);
 
  itkGetConstReferenceMacro(Beta, MatrixType);
 
  itkSetMacro(CovarianceMatrix, MatrixType);
 
  itkGetConstReferenceMacro(CovarianceMatrix, MatrixType);
 
  itkSetMacro(CovarianceInvMatrix, MatrixType);
 
  itkGetConstReferenceMacro(CovarianceInvMatrix, MatrixType);
 
  itkSetMacro(Vcondopt, MatrixType);
 
  itkGetConstReferenceMacro(Vcondopt, MatrixType);
 
  itkSetMacro(S, float);
 
  itkGetConstReferenceMacro(S, float);
 
protected:
  BayesianFusionFilter() : Superclass(BayesianFunctorType{}, {{0, 0}})
  {
    m_Lambda                      = 0.9999;
    m_S                           = 1;
    m_StatisticsHaveBeenGenerated = false;
  };
 
  ~BayesianFusionFilter() override{};
  void BeforeThreadedGenerateData() override;
 
  void ComputeInternalStatistics(void);
 
  void Modified(void) const override;
 
private:
  float m_Lambda;
  float m_S;
 
  MatrixType m_CovarianceMatrix;
 
  MatrixType m_CovarianceInvMatrix;
 
  MatrixType m_Beta;
 
  MatrixType m_Vcondopt;
 
  mutable bool m_StatisticsHaveBeenGenerated;
};
 
} // end namespace otb
 
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbBayesianFusionFilter.hxx"
#endif
 
#endif