otbLHMI.h

Go to the documentation of this file.

/*
 * 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 otbLHMI_h
#define otbLHMI_h
 
#include <vector>
#include "itkHistogram.h"
 
namespace otb
{
 
namespace Functor
{
 
template <class TInput1, class TInput2, class TOutput>
class LHMI
{
 
public:
  typedef typename std::vector<TOutput>         VectorType;
  typedef typename VectorType::iterator         IteratorType;
  typedef typename std::vector<VectorType>      VectorOfVectorType;
  typedef typename VectorOfVectorType::iterator VecOfVecIteratorType;
  typedef double                                HistogramFrequencyType;
  typedef typename itk::Statistics::Histogram<HistogramFrequencyType, itk::Statistics::DenseFrequencyContainer2> HistogramType;
  typedef typename HistogramType::MeasurementVectorType MeasurementVectorType;
  typedef typename HistogramType::SizeType              HistogramSizeType;
  typedef typename HistogramType::Iterator              HistogramIteratorType;
 
  LHMI()
  {
  }
  virtual ~LHMI()
  {
  }
  inline TOutput operator()(const TInput1& itA, const TInput2& itB)
  {
    HistogramType::Pointer histogram;
 
    HistogramSizeType histogramSize(2);
 
    MeasurementVectorType lowerBound(2);
 
    MeasurementVectorType upperBound(2);
 
    double upperBoundIncreaseFactor = 0.001;
 
    histogramSize.Fill(256);
 
    TOutput maxA = itA.GetPixel(0);
    TOutput minA = itA.GetPixel(0);
    TOutput maxB = itB.GetPixel(0);
    TOutput minB = itB.GetPixel(0);
 
    for (unsigned long pos = 0; pos < itA.Size(); ++pos)
    {
 
      TOutput value = static_cast<TOutput>(itA.GetPixel(pos));
 
      if (value > maxA)
        maxA = value;
      else if (value < minA)
        minA = value;
 
      value = static_cast<TOutput>(itB.GetPixel(pos));
 
      if (value > maxB)
        maxB = value;
      else if (value < minB)
        minB = value;
    }
 
    // Initialize the upper and lower bounds of the histogram.
    lowerBound[0] = minA;
    lowerBound[1] = minB;
    upperBound[0] = maxA + (maxA - minA) * upperBoundIncreaseFactor;
    upperBound[1] = maxB + (maxB - minB) * upperBoundIncreaseFactor;
 
    histogram = HistogramType::New();
 
    histogram->SetMeasurementVectorSize(2);
    histogram->Initialize(histogramSize, lowerBound, upperBound);
 
    for (unsigned long pos = 0; pos < itA.Size(); ++pos)
    {
 
      typename HistogramType::IndexType sample(2);
      sample[0] = itA.GetPixel(pos);
      sample[1] = itB.GetPixel(pos);
      /*if(sample[0]!=NumericTraits<TOutput>::Zero &&
         sample[1]!=NumericTraits<TOutput>::Zero)*/
      histogram->IncreaseFrequencyOfIndex(sample, 1);
    }
 
    TOutput                entropyX     = itk::NumericTraits<TOutput>::Zero;
    TOutput                entropyY     = itk::NumericTraits<TOutput>::Zero;
    TOutput                jointEntropy = itk::NumericTraits<TOutput>::Zero;
    HistogramFrequencyType totalFreq    = histogram->GetTotalFrequency();
 
    for (unsigned int i = 0; i < histogram->GetSize()[0]; ++i)
    {
      HistogramFrequencyType freq = histogram->GetFrequency(i, 0);
      if (freq > 0)
      {
        entropyX += freq * std::log(freq);
      }
    }
 
    entropyX = -entropyX / static_cast<TOutput>(totalFreq) + std::log(totalFreq);
 
    for (unsigned int i = 0; i < histogram->GetSize()[1]; ++i)
    {
      HistogramFrequencyType freq = histogram->GetFrequency(i, 1);
      if (freq > 0)
      {
        entropyY += freq * std::log(freq);
      }
    }
 
    entropyY = -entropyY / static_cast<TOutput>(totalFreq) + std::log(totalFreq);
 
    HistogramIteratorType it  = histogram->Begin();
    HistogramIteratorType end = histogram->End();
    while (it != end)
    {
      HistogramFrequencyType freq = it.GetFrequency();
      if (freq > 0)
      {
        jointEntropy += freq * std::log(freq);
      }
      ++it;
    }
 
    jointEntropy = -jointEntropy / static_cast<TOutput>(totalFreq) + std::log(totalFreq);
 
    return static_cast<TOutput>(jointEntropy / (entropyX + entropyY));
  }
};
}
}
#endif