histogram_analyzer.hpp

#pragma once
 
#include "../_cexports/c_light_meter.h"
 
#include "../global.hpp"
#include "../image.hpp"
#include "../exception.hpp"
#include "../rect.hpp"
#include "../area_2d.hpp"
 
#include <vector>
#include <cstdint>
#include <cmath>
#include <iterator>
#include <algorithm>
 
 
namespace Cvb
{
CVB_BEGIN_INLINE_NS
 
 
 
namespace Foundation
{
 
 
namespace HistogramAnalyzer
{
 
// forward declarations
class Histogram;
namespace Private
{
  std::vector<Histogram> CreateHistograms (const Image & image, Area2D aoi, double density, bool ignoreImageCs);
}
 
 
class Histogram
{
public:
  Histogram () noexcept
    : numPixels_(0), median_(0), mean_(0.0), mode_ (0.0), min_(0.0), max_(0.0), variance_(0.0), standardDeviation_(0.0), values_()
  {}
 
private:
  Histogram (CExports::LMH lmHandle, int plane)
    : numPixels_(0), median_(0), mean_(0.0), mode_ (0.0), min_(0.0), max_(0.0), variance_(0.0), standardDeviation_(0.0), values_()
  {
    // retrieve the histogram data
    size_t numEntries;
    CVB_CALL_CAPI_CHECKED (LMGetDataBufferEntries (lmHandle, plane, &numEntries));
    values_.resize (numEntries);
    for (size_t i = 0; i < numEntries; ++i)
    {
      double val;
      CVB_CALL_CAPI_CHECKED (LMGetSingleHistogramEntry(lmHandle, plane, i, &val));
      values_[i] = static_cast<uint64_t> (std::llround (val));
    }
 
    // retrieve the statistic results
    double out;
    CVB_CALL_CAPI_CHECKED (LMGetStatisticTotal (lmHandle, plane, &out));
    numPixels_ = static_cast<uint64_t> (std::llround (out));
    CVB_CALL_CAPI_CHECKED (LMGetStatisticMean (lmHandle, plane, &mean_));
    CVB_CALL_CAPI_CHECKED (LMGetStatisticMode (lmHandle, plane, &mode_));
    CVB_CALL_CAPI_CHECKED (LMGetStatisticMin (lmHandle, plane, &min_));
    CVB_CALL_CAPI_CHECKED (LMGetStatisticMax (lmHandle, plane, &max_));
    CVB_CALL_CAPI_CHECKED (LMGetStatisticVariance (lmHandle, plane, &variance_));
    CVB_CALL_CAPI_CHECKED (LMGetStatisticStdDev (lmHandle, plane, &standardDeviation_));
    median_ = CalculateMedianIndex();
  }
 
  friend std::vector<Histogram> Private::CreateHistograms (const Image & image, Area2D aoi, double density, bool ignoreImageCs);
 
  int CalculateMedianIndex() const noexcept
  {
    std::size_t idx = 0;
    for (uint64_t countBelow = 0; (idx < values_.size ()) && (countBelow < (numPixels_ / 2)); ++idx)
    {
      countBelow += values_[idx];
    }
    return static_cast<int>(idx);
  }
 
public:
 
  uint64_t NumPixels () const noexcept
  {
    return numPixels_;
  }
 
 
  uint64_t Median () const noexcept
  {
    // we decided not to call it median index
    return median_;
  }
 
 
  double Mean () const noexcept
  {
    return mean_;
  }
 
 
  double Mode () const noexcept
  {
    return mode_;
  }
 
 
  double Min () const noexcept
  {
    return min_;
  }
 
 
  double Max () const noexcept
  {
    return max_;
  }
 
 
  double Variance () const noexcept
  {
    return variance_;
  }
 
 
  double StandardDeviation () const noexcept
  {
    return standardDeviation_;
  }
 
 
  int Count () const noexcept
  {
    return static_cast<int>(values_.size());
  }
 
 
  std::vector<uint64_t> Values () const
  {
    return values_;
  }
 
 
  uint64_t operator[](int index) const
  {
    if (static_cast<std::size_t>(index) >= values_.size())
    {
      throw std::out_of_range ("histogram index out of range");
    }
 
    return values_[index];
  }
 
private:
  uint64_t numPixels_;
  int median_;
  double mean_;
  double mode_;
  double min_;
  double max_;
  double variance_;
  double standardDeviation_;
  std::vector<uint64_t> values_;
};
 
namespace Private
{
typedef HandleGuard<void, CVB_CALL_CAPI(ReleaseLMVoid)> ReleaseLMGuard;
 
/* Worker of the histogram creation functions */
inline std::vector<Histogram> CreateHistograms (const Image & image, Area2D aoi, double density, bool ignoreImageCs)
{
  std::vector<Histogram> histograms;
 
  CExports::LMH hLM;
  CVB_CALL_CAPI_CHECKED (LMCreate (&hLM));
  ReleaseLMGuard hLMHolder (reinterpret_cast<void*>(hLM));
 
  CVB_CALL_CAPI_CHECKED (LMSetImage(hLM, image.Handle()));
  for (CExports::cvbdim_t i = 0; i < image.PlanesCount (); ++i)
  {
    CVB_CALL_CAPI_CHECKED (LMSetEntireImageFlag (hLM, i, false));
    CVB_CALL_CAPI_CHECKED (LMSetProcessFlag (hLM, i, true));
  }
 
  for (CExports::cvbdim_t i = 0; i < image.PlanesCount(); ++i)
  {
    CVB_CALL_CAPI_CHECKED (LMSetDensity(hLM, i, std::lround(density * 1000.0)));
    CVB_CALL_CAPI_CHECKED (LMSetUseCSFlag(hLM, i, !ignoreImageCs));
    CVB_CALL_CAPI_CHECKED (LMSetArea(hLM, i, reinterpret_cast<CExports::TArea&>(aoi)));
  }
  CVB_CALL_CAPI_CHECKED (LMExecute(hLM));
 
  for (CExports::cvbdim_t i = 0; i < image.PlanesCount (); ++i)
  {
    histograms.push_back (Histogram (hLM, i));
  }
 
  return histograms;
}
 
} /* namespace Private */
 
 
inline std::vector<Histogram> CreateImageHistograms (const Image & image, Rect<int> aoi, double density = 1.0)
{
  return Private::CreateHistograms (image, Area2D(static_cast<Rect<double>>(aoi)), density, true);
}
 
 
inline std::vector<Histogram> CreateImageHistograms (const Image & image, Area2D aoi, double density = 1.0)
{
  return Private::CreateHistograms (image, aoi, density, false);
}
 
 
inline std::vector<Histogram> CreateImageHistograms (const Image & image, double density = 1.0)
{
  return CreateImageHistograms (image, image.Bounds(), density);
}
 
 
inline Histogram CreatePlaneHistogram (const ImagePlane & imagePlane, Rect<int> aoi, double density = 1.0)
{
  return CreateImageHistograms (*imagePlane.Map(), aoi, density) [0];
}
 
 
inline Histogram CreatePlaneHistogram (const ImagePlane & imagePlane, Area2D aoi, double density = 1.0)
{
  return CreateImageHistograms (*imagePlane.Map(), aoi, density) [0];
}
 
 
inline Histogram CreatePlaneHistogram (const ImagePlane & imagePlane, double density = 1.0)
{
  return CreateImageHistograms (*imagePlane.Map(), density) [0];
}
 
namespace Private
{
 
struct PeakData
{
  int GrayValue;
  double Quality;
};
 
inline std::vector<double> MakeDerivationFilter ()
{
  return std::vector<double> { -1.0, 0.0, 1.0 };
}
 
inline std::vector<double> MakeBlurFilter (int blurSize)
{
  return std::vector<double> (blurSize, 1.0 / blurSize);
}
 
template <class HISTVAL>
inline std::vector<double> MakeBorderedHistogram (const std::vector<HISTVAL> &histogram, size_t kernelSize)
{
  std::vector<double> result (histogram.size() + kernelSize - 1);
 
  size_t j = 0;
  auto firstValue = histogram[0];
  for (size_t i = 0, bckOffset = kernelSize / 2; i < bckOffset; ++i)
  {
    result[j++] = static_cast<double> (firstValue);
  }
 
  for (size_t i = 0; i < histogram.size(); ++i)
  {
    result[j++] = static_cast<double> (histogram[i]);
  }
 
  auto lastValue = histogram[histogram.size() - 1];
  for (; j < result.size (); ++j)
  {
    result[j] = static_cast<double> (lastValue);
  }
 
  return result;
}
 
template <class RANGE>
inline typename TypedRange<std::vector<double>, double, RANGE>::type ApplyKernel (const std::vector<double> &borderedHistogram, const RANGE &kernel)
{
  auto kernelRange = MakeRangeAdapter<double> (kernel, 1);
  std::vector<double> result (borderedHistogram.size() - kernelRange.Size() + 1);
 
  for (size_t i = 0; i < result.size(); ++i)
  {
    double val = 0.0;
    for (size_t k = 0; k < kernelRange.Size (); ++k)
    {
      val += borderedHistogram[i + k] * kernel[k];
    }
    result[i] = val;
  }
  return result;
}
 
inline std::vector<PeakData> FindZeroCrossings (const std::vector<double> &blurredHistogram, const std::vector<double> &derivedHistogram)
{
  // for initialization we need to make sure that a value != 0 is loaded
  auto itNonZero = std::find_if (derivedHistogram.begin (), derivedHistogram.end (), [](double val) {return val != 0.0; });
  if (itNonZero == derivedHistogram.end())
  {
    return std::vector<PeakData> ();
  }
  // turn the iterator to index
  size_t i = std::distance (derivedHistogram.begin(), itNonZero);
 
  auto sign = [](double val) {return (val == 0.0) ? 0 : ((val < 0.0) ? -1 : 1); };
 
  std::vector<PeakData> positions;
  auto lastStart = i++;
  auto lastSign = sign (derivedHistogram[lastStart]);
  for (; i < derivedHistogram.size(); i++)
  {
    // three cases sign before == sign after; sign before != sign after;
    // sign after == 0
    auto thisSign = sign (derivedHistogram[i]);
    if (thisSign != 0)
    {
      if (thisSign != lastSign)
      {
        auto peakPos = (i + lastStart) / 2;
        positions.push_back(PeakData{static_cast<int>(peakPos), blurredHistogram[peakPos]});
      }
      lastSign = thisSign;
      lastStart = i;
    }
  }
  return positions;
}
 
inline std::vector<int> FilterByMinDistance (const std::vector<PeakData> &zeroCrossings, int minDiff)
{
  std::vector<PeakData> descendingCrossings (zeroCrossings);
  std::sort (descendingCrossings.begin (), descendingCrossings.end (), [](PeakData first, PeakData second) { return first.Quality > second.Quality; });
  std::vector<int> values (descendingCrossings.size());
  std::transform (descendingCrossings.begin (), descendingCrossings.end (), values.begin (), [](PeakData peak) { return peak.GrayValue; });
 
  for (size_t i = 0; i < values.size() - 1; ++i)
  {
    for (size_t k = i + 1; k < values.size(); ++k)
    {
      auto diff = std::abs (values[i] - values[k]);
      if (diff < minDiff)
      {
        values.erase (values.begin() + k--);
      }
    }
  }
 
  return values;
}
 
} /* namespace Private */
 
 
 
template <class HISTVAL, class RANGE>
inline typename TypedRange<std::vector<double>, double, RANGE>::type FilterHistogram (const std::vector<HISTVAL> &histogram, const RANGE &kernel)
{
  auto kernelSize = std::distance(std::begin(kernel), std::end(kernel));
  if ((kernelSize < 1) || (static_cast<size_t>(kernelSize) > histogram.size()))
  {
    throw std::out_of_range ("histogram filter kernel size out of range");
  }
 
  auto borderedHistogram = Private::MakeBorderedHistogram (histogram, kernelSize);
  return Private::ApplyKernel (borderedHistogram, kernel);
}
 
 
inline std::vector<int> FindHistogramPeaks (const std::vector<uint64_t> &histogram, int blurSize, int minDiff)
{
  if ((blurSize < 1) || (static_cast<size_t>(blurSize) > histogram.size()) || (minDiff < 1))
  {
    throw std::out_of_range ("histogram peak search parameters out of range");
  }
 
  auto blurredHistogram = FilterHistogram (histogram, Private::MakeBlurFilter(static_cast<size_t>(blurSize)));
  auto derivedHistogram = FilterHistogram (blurredHistogram, Private::MakeDerivationFilter());
  auto zeroCrossings = Private::FindZeroCrossings (blurredHistogram, derivedHistogram);
  return Private::FilterByMinDistance (zeroCrossings, minDiff);
}
 
 
inline uint64_t SumHistogramBetween (const std::vector<uint64_t> &histogram, int lowerLimit, int upperLimit)
{
  if ((lowerLimit < 0) || (static_cast<size_t>(upperLimit) >= histogram.size()))
  {
    throw std::out_of_range ("histogram peak search parameters out of range");
  }
 
  uint64_t nres = 0;
  for (int i = lowerLimit; i <= upperLimit; ++i)
  {
    nres += histogram[i];
  }
  return nres;
}
 
} /* namespace HistogramAnalyzer */
 
using HistogramAnalyzer::Histogram;
 
using HistogramAnalyzer::CreateImageHistograms;
using HistogramAnalyzer::CreatePlaneHistogram;
using HistogramAnalyzer::FilterHistogram;
using HistogramAnalyzer::FindHistogramPeaks;
using HistogramAnalyzer::SumHistogramBetween;
 
 
} /* namespace Foundation */
CVB_END_INLINE_NS
} /* namespace Cvb */