otbBayesianFusionFilter.hxx

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 otbBayesianFusionFilter_hxx
#define otbBayesianFusionFilter_hxx
 
#include "otbBayesianFusionFilter.h"
 
namespace otb
{
 
template <class TInputMultiSpectralImage, class TInputMultiSpectralInterpImage, class TInputPanchroImage, class TOutputImage>
void BayesianFusionFilter<TInputMultiSpectralImage, TInputMultiSpectralInterpImage, TInputPanchroImage, TOutputImage>::Modified() const
{
  Superclass::Modified();
  m_StatisticsHaveBeenGenerated = false;
}
 
template <class TInputMultiSpectralImage, class TInputMultiSpectralInterpImage, class TInputPanchroImage, class TOutputImage>
void BayesianFusionFilter<TInputMultiSpectralImage, TInputMultiSpectralInterpImage, TInputPanchroImage, TOutputImage>::BeforeThreadedGenerateData()
{
  if (!m_StatisticsHaveBeenGenerated)
  {
    this->ComputeInternalStatistics();
    m_StatisticsHaveBeenGenerated = true;
  }
}
 
template <class TInputMultiSpectralImage, class TInputMultiSpectralInterpImage, class TInputPanchroImage, class TOutputImage>
void BayesianFusionFilter<TInputMultiSpectralImage, TInputMultiSpectralInterpImage, TInputPanchroImage, TOutputImage>::ComputeInternalStatistics()
{
  OutputImageRegionType msiRequestedRegion     = this->GetMultiSpectInterp()->GetRequestedRegion();
  OutputImageRegionType msRequestedRegion      = this->GetMultiSpect()->GetRequestedRegion();
  OutputImageRegionType panchroRequestedRegion = this->GetPanchro()->GetRequestedRegion();
 
  // Allocate output
  typename OutputImageType::Pointer                   output          = this->GetOutput();
  typename InputMultiSpectralImageType::Pointer       multiSpec       = const_cast<InputMultiSpectralImageType*>(this->GetMultiSpect());
  typename InputMultiSpectralInterpImageType::Pointer multiSpecInterp = const_cast<InputMultiSpectralInterpImageType*>(this->GetMultiSpectInterp());
  typename InputPanchroImageType::Pointer             panchro         = const_cast<InputPanchroImageType*>(this->GetPanchro());
 
  m_Beta.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel() + 1, 1);
  m_Beta.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
 
  m_CovarianceMatrix.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel(), multiSpecInterp->GetNumberOfComponentsPerPixel());
  m_CovarianceMatrix.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
 
  m_CovarianceInvMatrix.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel(), multiSpecInterp->GetNumberOfComponentsPerPixel());
  m_CovarianceInvMatrix.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
 
  m_Vcondopt.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel(), multiSpecInterp->GetNumberOfComponentsPerPixel());
  m_Vcondopt.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
  typename StreamingStatisticsVectorImageFilterType::Pointer covComputefilter = StreamingStatisticsVectorImageFilterType::New();
 
  covComputefilter->SetInput(multiSpecInterp);
  covComputefilter->Update();
 
  m_CovarianceMatrix = covComputefilter->GetCovariance();
  otbMsgDebugMacro(<< "Covariance: " << m_CovarianceMatrix);
 
  m_CovarianceInvMatrix = m_CovarianceMatrix.GetInverse();
  // MatrixTransform only support vectorimage input
  typename CasterType::Pointer caster = CasterType::New();
  caster->SetInput(panchro);
  // caster->Update();
  // Compute the transpose multispectral image multiplied by itself
  typename MSTransposeMSType::Pointer msTransposeMs = MSTransposeMSType::New();
  // Compute the transpose multispectral image multiplied by the panchromatic one
  typename MSTransposeMSType::Pointer msTransposePan = MSTransposeMSType::New();
  // Add a dimension filled with ones to the images
  msTransposeMs->SetUsePadFirstInput(true);
  msTransposeMs->SetUsePadSecondInput(true);
  msTransposePan->SetUsePadFirstInput(true);
 
  msTransposeMs->SetFirstInput(multiSpec);
  msTransposeMs->SetSecondInput(multiSpec);
 
  msTransposePan->SetFirstInput(multiSpec);
  msTransposePan->SetSecondInput(caster->GetOutput());
 
  msTransposeMs->Update();
  otbMsgDebugMacro(<< "MsTMs: " << msTransposeMs->GetResultOutput()->Get());
  msTransposePan->Update();
  otbMsgDebugMacro(<< "MsTPan: " << msTransposePan->GetResultOutput()->Get());
 
  MatrixType temp;
  temp   = msTransposeMs->GetResultOutput()->Get().GetInverse();
  m_Beta = temp * msTransposePan->GetResultOutput()->Get();
 
  // S computation : quadratique mean of the regression residue
  // Compute the transpose panchromatic image multiplied by itself
  typename MSTransposeMSType::Pointer panTransposePan = MSTransposeMSType::New();
  panTransposePan->SetFirstInput(caster->GetOutput());
  panTransposePan->SetSecondInput(caster->GetOutput());
  panTransposePan->Update();
  otbMsgDebugMacro(<< "PanTPan: " << msTransposePan->GetResultOutput()->Get());
  MatrixType S, tempS, tempS2;
  S     = panTransposePan->GetResultOutput()->Get();
  tempS = msTransposePan->GetResultOutput()->Get().GetTranspose();
  tempS = tempS * m_Beta;
 
  if ((S.Rows() != tempS.Rows()) || (S.Cols() != tempS.Cols()))
  {
    itkExceptionMacro(<< "Matrix with size (" << S.Rows() << "," << S.Cols() << ") cannot be subtracted from matrix with size (" << tempS.Rows() << ","
                      << tempS.Cols() << " )");
  }
  for (unsigned int r = 0; r < S.Rows(); ++r)
  {
    for (unsigned int c = 0; c < S.Cols(); ++c)
    {
      S(r, c) -= tempS(r, c);
    }
  }
 
/**** END TODO ****/
 
  tempS  = m_Beta.GetTranspose();
  tempS2 = msTransposePan->GetResultOutput()->Get();
  tempS  = tempS * tempS2;
  if ((S.Rows() != tempS.Rows()) || (S.Cols() != tempS.Cols()))
  {
    itkExceptionMacro(<< "Matrix with size (" << S.Rows() << "," << S.Cols() << ") cannot be subtracted from matrix with size (" << tempS.Rows() << ","
                      << tempS.Cols() << " )");
  }
  for (unsigned int r = 0; r < S.Rows(); ++r)
  {
    for (unsigned int c = 0; c < S.Cols(); ++c)
    {
      S(r, c) -= tempS(r, c);
    }
  }
 
/**** END TODO ****/
 
  MatrixType xxT, xxTb, xxTbT, xxTbTb;
  xxT    = msTransposeMs->GetResultOutput()->Get().GetTranspose();
  xxTb   = xxT * m_Beta;
  xxTbT  = xxTb.GetTranspose();
  xxTbTb = xxTbT * m_Beta;
  if ((S.Rows() != xxTbTb.Rows()) || (S.Cols() != xxTbTb.Cols()))
  {
    itkExceptionMacro(<< "Matrix with size (" << S.Rows() << "," << S.Cols() << ") cannot be subtracted from matrix with size (" << xxTbTb.Rows() << ","
                      << xxTbTb.Cols() << " )");
  }
 
  for (unsigned int r = 0; r < S.Rows(); ++r)
  {
    for (unsigned int c = 0; c < S.Cols(); ++c)
    {
      S(r, c) += xxTbTb(r, c);
    }
  }
/**** END TODO ****/
 
  unsigned int size1 = multiSpec->GetLargestPossibleRegion().GetSize()[0] * multiSpec->GetLargestPossibleRegion().GetSize()[1];
  unsigned int size2 = multiSpec->GetNumberOfComponentsPerPixel() + 1;
  m_S                = S(0, 0);
  m_S /= static_cast<float>(size1 - size2);
 
  // cutBeta is the N-1 last m_Beta element matrix.
  // varPan contains transpose(cutBeta)*cutBeta/S
  MatrixType varPan, cutBeta;
  varPan.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel(), 1);
  varPan.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
  cutBeta.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel(), 1);
  cutBeta.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
  // Take the N-1 m_Beta last elements
  for (unsigned int r = 1; r < m_Beta.Rows(); ++r)
  {
    cutBeta(r - 1, 0) = m_Beta(r, 0);
  }
  varPan = cutBeta;
 
  MatrixType tempvarPan;
  tempvarPan = varPan.GetTranspose();
  varPan *= tempvarPan;
  varPan /= m_S;
  // Compute the optimization matrix : m_Vcondopt
  // eye is the identical matrix which size is the number of components of the multispectral image
  MatrixType eye;
  eye.SetSize(multiSpecInterp->GetNumberOfComponentsPerPixel(), multiSpecInterp->GetNumberOfComponentsPerPixel());
  eye.Fill(itk::NumericTraits<InputMultiSpectralInterpRealType>::Zero);
  for (unsigned int r = 1; r < eye.Rows(); ++r)
  {
    eye(r, r) = std::pow(10., -12.);
  }
 
  if ((m_Vcondopt.Rows() != varPan.Rows()) || (m_Vcondopt.Cols() != varPan.Cols()) || (m_Vcondopt.Rows() != m_CovarianceInvMatrix.Rows()) ||
      (m_Vcondopt.Cols() != m_CovarianceInvMatrix.Cols()))
  {
    itkExceptionMacro(<< "Matrix with size (" << m_Vcondopt.Rows() << "," << m_Vcondopt.Cols() << ") cannot be subtracted from matrix with size ("
                      << varPan.Rows() << "," << varPan.Cols() << " ) or ( " << m_CovarianceInvMatrix.Rows() << "," << m_CovarianceInvMatrix.Cols() << ")");
  }
  for (unsigned int r = 0; r < m_Vcondopt.Rows(); ++r)
  {
    for (unsigned int c = 0; c < m_Vcondopt.Cols(); ++c)
    {
      m_Vcondopt(r, c) = 2 * m_Lambda * varPan(r, c) + 2 * m_CovarianceInvMatrix(r, c) * (1 - m_Lambda) + eye(r, c);
    }
  }
 
  m_Vcondopt = m_Vcondopt.GetInverse();
  // Functor initialization
  this->GetModifiableFunctor().SetVcondopt(m_Vcondopt);
  this->GetModifiableFunctor().SetBeta(cutBeta);
  this->GetModifiableFunctor().SetAlpha(m_Beta(0, 0));
  this->GetModifiableFunctor().SetCovarianceInvMatrix(m_CovarianceInvMatrix);
  this->GetModifiableFunctor().SetLambda(m_Lambda);
  this->GetModifiableFunctor().SetS(m_S);
 
  // Restore the previous buffered data
  multiSpecInterp->SetRequestedRegion(msiRequestedRegion);
  multiSpecInterp->PropagateRequestedRegion();
  multiSpecInterp->UpdateOutputData();
 
  multiSpec->SetRequestedRegion(msRequestedRegion);
  multiSpec->PropagateRequestedRegion();
  multiSpec->UpdateOutputData();
 
  panchro->SetRequestedRegion(panchroRequestedRegion);
  panchro->PropagateRequestedRegion();
  panchro->UpdateOutputData();
}
} // end namespace otb
 
#endif