[libcamera-devel,v4,2/4] WIP: ipa: ipu3: Add an histogram class
diff mbox series

Message ID 20210330211210.194806-3-jeanmichel.hautbois@ideasonboard.com
State Superseded
Headers show
Series
  • Implement IPA algorithms and demo with IPU3
Related show

Commit Message

Jean-Michel Hautbois March 30, 2021, 9:12 p.m. UTC
This class will be used at least by AGC algorithm when quantiles are
needed for example.

Signed-off-by: Jean-Michel Hautbois <jeanmichel.hautbois@ideasonboard.com>
---
 src/ipa/libipa/histogram.cpp                 | 154 +++++++++++++++++++
 src/ipa/libipa/histogram.h                   |  40 +++++
 src/ipa/libipa/meson.build                   |   2 +
 src/ipa/raspberrypi/controller/histogram.hpp |   9 +-
 4 files changed, 197 insertions(+), 8 deletions(-)
 create mode 100644 src/ipa/libipa/histogram.cpp
 create mode 100644 src/ipa/libipa/histogram.h

Comments

Kieran Bingham April 12, 2021, 12:21 p.m. UTC | #1
Hi Jean-Michel,

in $SUBJECT

s/an/a/



On 30/03/2021 22:12, Jean-Michel Hautbois wrote:
> This class will be used at least by AGC algorithm when quantiles are
> needed for example.
> 

Is this really a Histogram class? Or is is a CumulativeFrequency class?

I may likely be mis-interpretting or just wrong - but I thought a
histogram stored values like:


value 0 1 2 3 4 5 6 7 8 9
count 5 0 3 5 6 9 2 0 2 5

Meaning that ... you can look up in the table to see that there are 5
occurrences of the value 9 in the dataset... whereas this stores:


value 0 1 2  3  4  5  6  7  8  9
count 5 8 8 13 19 28 30 30 32 37


Ok, so now I've drawn that out, I can see that the same information is
still there, as you can get the occurrences of any specific value by
subtracting from the previous 'bin'?


> Signed-off-by: Jean-Michel Hautbois <jeanmichel.hautbois@ideasonboard.com>
> ---
>  src/ipa/libipa/histogram.cpp                 | 154 +++++++++++++++++++
>  src/ipa/libipa/histogram.h                   |  40 +++++
>  src/ipa/libipa/meson.build                   |   2 +
>  src/ipa/raspberrypi/controller/histogram.hpp |   9 +-
>  4 files changed, 197 insertions(+), 8 deletions(-)
>  create mode 100644 src/ipa/libipa/histogram.cpp
>  create mode 100644 src/ipa/libipa/histogram.h
> 
> diff --git a/src/ipa/libipa/histogram.cpp b/src/ipa/libipa/histogram.cpp
> new file mode 100644
> index 00000000..46fbb940
> --- /dev/null
> +++ b/src/ipa/libipa/histogram.cpp
> @@ -0,0 +1,154 @@
> +/* SPDX-License-Identifier: BSD-2-Clause */
> +/*
> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
> + *
> + * histogram.cpp - histogram calculations
> + */
> +#include "histogram.h"
> +
> +#include <cmath>
> +
> +#include "libcamera/internal/log.h"
> +
> +/**
> + * \file histogram.h
> + * \brief Class to represent Histograms and manipulate them
> + */
> +
> +namespace libcamera {
> +
> +namespace ipa {
> +
> +/**
> + * \class Histogram
> + * \brief The base class for creating histograms
> + *
> + * The Histogram class defines a standard interface for IPA algorithms. By
> + * abstracting histograms, it makes it possible to implement generic code
> + * to manage histograms regardless of their specific type.

I don't think this defines a standard interface for IPA algorithms ;)

I think we can drop that paragraph and keep the one below.

> + *
> + * This class stores a cumulative frequency histogram, which is a mapping that
> + * counts the cumulative number of observations in all of the bins up to the
> + * specified bin. It can be used to find quantiles and averages between quantiles.
> + */
> +
> +/**
> + * \brief Create a cumulative histogram with a bin number of intervals

I suspect 'bin' was a previous value passed in here.

I think we should describe how data should be passed in.
I.e. I think perhaps it should be a pre-sorted histogram or something?


> + * \param[in] data a reference to the histogram

I don't think this is a reference anymore.

> + */
> +Histogram::Histogram(Span<uint32_t> data)
> +{
> +	cumulative_.reserve(data.size());
> +	cumulative_.push_back(0);
> +	for (const uint32_t &value : data)
> +		cumulative_.push_back(cumulative_.back() + value);
> +}
> +/**
> + * \fn Histogram::bins()
> + * \brief getter for number of bins

We don't normally reference them as 'getter's even if that's what they do.

To be consistent with other briefs, we should put something like:

  \brief Retrieve the number of bins currently stored in the Histogram


(stored by, might be 'used by'?)


> + * \return Number of bins
> + */
> +/**
> + * \fn Histogram::total()
> + * \brief getter for number of values

Same here,

 \brief Retrieve the total number of values in the data set

> + * \return Number of values
> + */
> +
> +/**
> + * \brief Cumulative frequency up to a (fractional) point in a bin.
> + * \param[in] bin the bin upon which to cumulate

s/the/The/
(Capital letter to start after the parameter itself.)

Reading below, does it also mean this should say "up to which" rather
than "upon which"?

i.e.
   upon which: to me, counts the values 'in' that bin.
   up to which: to me, counts the values up to that one...

> + *
> + * With F(p) the cumulative frequency of the histogram, the value is 0 at

What is F(p) here?



> + * the bottom of the histogram, and the maximum is the number of bins.
> + * The pixels are spread evenly throughout the “bin” in which they lie, so that
> + * F(p) is a continuous (monotonically increasing) function.
> + *
> + * \return The cumulated frequency from 0 up to the specified bin
> + */
> +uint64_t Histogram::cumulativeFreq(double bin) const
> +{
> +	if (bin <= 0)
> +		return 0;
> +	else if (bin >= bins())
> +		return total();
> +	int b = static_cast<int32_t>(bin);
> +	return cumulative_[b] +
> +	       (bin - b) * (cumulative_[b + 1] - cumulative_[b]);
> +}
> +
> +/**
> + * \brief Return the (fractional) bin of the point through the histogram
> + * \param[in] q the desired point (0 <= q <= 1)
> + * \param[in] first low limit (default is 0)
> + * \param[in] last high limit (default is UINT_MAX)
> + *
> + * A quantile gives us the point p = Q(q) in the range such that a proportion
> + * q of the pixels lie below p. A familiar quantile is Q(0.5) which is the median
> + * of a distribution.
> + *
> + * \return The fractional bin of the point
> + */
> +double Histogram::quantile(double q, uint32_t first, uint32_t last) const
> +{
> +	if (last == UINT_MAX)
> +		last = cumulative_.size() - 2;
> +	ASSERT(first <= last);
> +
> +	uint64_t item = q * total();
> +	/* Binary search to find the right bin */
> +	while (first < last) {
> +		int middle = (first + last) / 2;
> +		/* Is it between first and middle ? */
> +		if (cumulative_[middle + 1] > item)
> +			last = middle;
> +		else
> +			first = middle + 1;
> +	}
> +	ASSERT(item >= cumulative_[first] && item <= cumulative_[last + 1]);
> +
> +	double frac;
> +	if (cumulative_[first + 1] == cumulative_[first])
> +		frac = 0;
> +	else
> +		frac = (item - cumulative_[first]) / (cumulative_[first + 1] - cumulative_[first]);
> +	return first + frac;
> +}
> +
> +/**
> + * \brief Calculate the mean between two quantiles
> + * \param[in] lowQuantile low Quantile
> + * \param[in] highQuantile high Quantile
> + *
> + * Quantiles are not ideal for metering as they suffer several limitations.
> + * Instead, a concept is introduced here: inter-quantile mean.
> + * It returns the mean of all pixels between lowQuantile and highQuantile.
> + *
> + * \return The mean histogram bin value between the two quantiles
> + */
> +double Histogram::interQuantileMean(double lowQuantile, double highQuantile) const
> +{
> +	ASSERT(highQuantile > lowQuantile);
> +	/* Proportion of pixels which lies below lowQuantile */
> +	double lowPoint = quantile(lowQuantile);
> +	/* Proportion of pixels which lies below highQuantile */
> +	double highPoint = quantile(highQuantile, static_cast<uint32_t>(lowPoint));
> +	double sumBinFreq = 0, cumulFreq = 0;
> +
> +	for (double p_next = floor(lowPoint) + 1.0;
> +	     p_next <= ceil(highPoint);
> +	     lowPoint = p_next, p_next += 1.0) {
> +		int bin = floor(lowPoint);
> +		double freq = (cumulative_[bin + 1] - cumulative_[bin]) * (std::min(p_next, highPoint) - lowPoint);
> +
> +		/* Accumulate weigthed bin */
> +		sumBinFreq += bin * freq;
> +		/* Accumulate weights */
> +		cumulFreq += freq;
> +	}
> +	/* add 0.5 to give an average for bin mid-points */
> +	return sumBinFreq / cumulFreq + 0.5;
> +}
> +
> +} /* namespace ipa */
> +
> +} /* namespace libcamera */
> diff --git a/src/ipa/libipa/histogram.h b/src/ipa/libipa/histogram.h
> new file mode 100644
> index 00000000..dc7451aa
> --- /dev/null
> +++ b/src/ipa/libipa/histogram.h
> @@ -0,0 +1,40 @@
> +/* SPDX-License-Identifier: BSD-2-Clause */
> +/*
> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
> + *
> + * histogram.h - histogram calculation interface
> + */
> +#ifndef __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
> +#define __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
> +
> +#include <assert.h>
> +#include <limits.h>
> +#include <stdint.h>
> +
> +#include <vector>
> +
> +#include <libcamera/span.h>
> +
> +namespace libcamera {
> +
> +namespace ipa {
> +
> +class Histogram
> +{
> +public:
> +	Histogram(Span<uint32_t> data);
> +	size_t bins() const { return cumulative_.size() - 1; }
> +	uint64_t total() const { return cumulative_[cumulative_.size() - 1]; }
> +	uint64_t cumulativeFreq(double bin) const;

I'd make this cumulativeFrequency()

> +	double quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
> +	double interQuantileMean(double lowQuantile, double hiQuantile) const;
> +
> +private:
> +	std::vector<uint64_t> cumulative_;

I see this came from the RPi code.
Do we really store 64 bit values in here? or are they only 32bit? (or
smaller?)

In fact, I see we now construct with a Span<uint32_t> data, so
presumably this vector can be uint32_t.

We could template it to accept different sizes if that would help ...
but I think supporting 32 bits is probably fine if that's the span we
currently define....

> +};
> +
> +} /* namespace ipa */
> +
> +} /* namespace libcamera */
> +
> +#endif /* __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__ */
> diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
> index 1819711d..038fc490 100644
> --- a/src/ipa/libipa/meson.build
> +++ b/src/ipa/libipa/meson.build
> @@ -2,10 +2,12 @@
>  
>  libipa_headers = files([
>      'algorithm.h',
> +    'histogram.h'
>  ])
>  
>  libipa_sources = files([
>      'algorithm.cpp',
> +    'histogram.cpp'
>  ])
>  
>  libipa_includes = include_directories('..')
> diff --git a/src/ipa/raspberrypi/controller/histogram.hpp b/src/ipa/raspberrypi/controller/histogram.hpp
> index 90f5ac78..fc236416 100644
> --- a/src/ipa/raspberrypi/controller/histogram.hpp
> +++ b/src/ipa/raspberrypi/controller/histogram.hpp

I'm not sure if the modifications to the RPi histogram below should be
in this patch.

It looks like they're unrelated to the actual code move?

Perhaps leave this out, and we can 'convert' RPi controller to use the
'generic' histogram separately?


> @@ -10,9 +10,6 @@
>  #include <vector>
>  #include <cassert>
>  
> -// A simple histogram class, for use in particular to find "quantiles" and
> -// averages between "quantiles".
> -
>  namespace RPiController {
>  
>  class Histogram
> @@ -29,12 +26,8 @@ public:
>  	}
>  	uint32_t Bins() const { return cumulative_.size() - 1; }
>  	uint64_t Total() const { return cumulative_[cumulative_.size() - 1]; }
> -	// Cumulative frequency up to a (fractional) point in a bin.
>  	uint64_t CumulativeFreq(double bin) const;
> -	// Return the (fractional) bin of the point q (0 <= q <= 1) through the
> -	// histogram. Optionally provide limits to help.
> -	double Quantile(double q, int first = -1, int last = -1) const;
> -	// Return the average histogram bin value between the two quantiles.
> +	Histogram::quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
>  	double InterQuantileMean(double q_lo, double q_hi) const;
>  
>  private:
>
David Plowman April 12, 2021, 1:29 p.m. UTC | #2
Hi Jean-Michel, Kieran

Happy to see this code getting some use elsewhere! I thought I could
just add a couple of comments on some of our original thinking...

On Mon, 12 Apr 2021 at 13:21, Kieran Bingham
<kieran.bingham@ideasonboard.com> wrote:
>
> Hi Jean-Michel,
>
> in $SUBJECT
>
> s/an/a/
>
>
>
> On 30/03/2021 22:12, Jean-Michel Hautbois wrote:
> > This class will be used at least by AGC algorithm when quantiles are
> > needed for example.
> >
>
> Is this really a Histogram class? Or is is a CumulativeFrequency class?
>
> I may likely be mis-interpretting or just wrong - but I thought a
> histogram stored values like:
>
>
> value 0 1 2 3 4 5 6 7 8 9
> count 5 0 3 5 6 9 2 0 2 5
>
> Meaning that ... you can look up in the table to see that there are 5
> occurrences of the value 9 in the dataset... whereas this stores:
>
>
> value 0 1 2  3  4  5  6  7  8  9
> count 5 8 8 13 19 28 30 30 32 37
>
>
> Ok, so now I've drawn that out, I can see that the same information is
> still there, as you can get the occurrences of any specific value by
> subtracting from the previous 'bin'?

I think you've pretty much convinced yourself why it stores cumulative
frequency. Going from cumulative frequency back to per-bin values is a
single subtraction. Going the other way is a loop. Given that finding
"quantiles" is a fairly useful operation, cumulative frequencies make
sense.

>
>
> > Signed-off-by: Jean-Michel Hautbois <jeanmichel.hautbois@ideasonboard.com>
> > ---
> >  src/ipa/libipa/histogram.cpp                 | 154 +++++++++++++++++++
> >  src/ipa/libipa/histogram.h                   |  40 +++++
> >  src/ipa/libipa/meson.build                   |   2 +
> >  src/ipa/raspberrypi/controller/histogram.hpp |   9 +-
> >  4 files changed, 197 insertions(+), 8 deletions(-)
> >  create mode 100644 src/ipa/libipa/histogram.cpp
> >  create mode 100644 src/ipa/libipa/histogram.h

Will this header be includable by application code? I've sometimes
found myself getting an image and then wanting to do histogram
operations. (Of course, this class doesn't make the Histogram from
scratch, you have to start by totting up the bins yourself, but that's
not difficult and then you get means, quantiles and stuff for free...)

> >
> > diff --git a/src/ipa/libipa/histogram.cpp b/src/ipa/libipa/histogram.cpp
> > new file mode 100644
> > index 00000000..46fbb940
> > --- /dev/null
> > +++ b/src/ipa/libipa/histogram.cpp
> > @@ -0,0 +1,154 @@
> > +/* SPDX-License-Identifier: BSD-2-Clause */
> > +/*
> > + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
> > + *
> > + * histogram.cpp - histogram calculations
> > + */
> > +#include "histogram.h"
> > +
> > +#include <cmath>
> > +
> > +#include "libcamera/internal/log.h"
> > +
> > +/**
> > + * \file histogram.h
> > + * \brief Class to represent Histograms and manipulate them
> > + */
> > +
> > +namespace libcamera {
> > +
> > +namespace ipa {
> > +
> > +/**
> > + * \class Histogram
> > + * \brief The base class for creating histograms
> > + *
> > + * The Histogram class defines a standard interface for IPA algorithms. By
> > + * abstracting histograms, it makes it possible to implement generic code
> > + * to manage histograms regardless of their specific type.
>
> I don't think this defines a standard interface for IPA algorithms ;)
>
> I think we can drop that paragraph and keep the one below.
>
> > + *
> > + * This class stores a cumulative frequency histogram, which is a mapping that
> > + * counts the cumulative number of observations in all of the bins up to the
> > + * specified bin. It can be used to find quantiles and averages between quantiles.
> > + */
> > +
> > +/**
> > + * \brief Create a cumulative histogram with a bin number of intervals
>
> I suspect 'bin' was a previous value passed in here.
>
> I think we should describe how data should be passed in.
> I.e. I think perhaps it should be a pre-sorted histogram or something?
>
>
> > + * \param[in] data a reference to the histogram
>
> I don't think this is a reference anymore.
>
> > + */
> > +Histogram::Histogram(Span<uint32_t> data)
> > +{
> > +     cumulative_.reserve(data.size());
> > +     cumulative_.push_back(0);
> > +     for (const uint32_t &value : data)
> > +             cumulative_.push_back(cumulative_.back() + value);
> > +}
> > +/**
> > + * \fn Histogram::bins()
> > + * \brief getter for number of bins
>
> We don't normally reference them as 'getter's even if that's what they do.
>
> To be consistent with other briefs, we should put something like:
>
>   \brief Retrieve the number of bins currently stored in the Histogram
>
>
> (stored by, might be 'used by'?)
>
>
> > + * \return Number of bins
> > + */
> > +/**
> > + * \fn Histogram::total()
> > + * \brief getter for number of values
>
> Same here,
>
>  \brief Retrieve the total number of values in the data set
>
> > + * \return Number of values
> > + */
> > +
> > +/**
> > + * \brief Cumulative frequency up to a (fractional) point in a bin.
> > + * \param[in] bin the bin upon which to cumulate
>
> s/the/The/
> (Capital letter to start after the parameter itself.)
>
> Reading below, does it also mean this should say "up to which" rather
> than "upon which"?
>
> i.e.
>    upon which: to me, counts the values 'in' that bin.
>    up to which: to me, counts the values up to that one...

Also, I don't think there's a verb "to cumulate", is there?

>
> > + *
> > + * With F(p) the cumulative frequency of the histogram, the value is 0 at
>
> What is F(p) here?

I think this is defining it to be the cumulative frequency.

>
>
>
> > + * the bottom of the histogram, and the maximum is the number of bins.
> > + * The pixels are spread evenly throughout the “bin” in which they lie, so that
> > + * F(p) is a continuous (monotonically increasing) function.
> > + *
> > + * \return The cumulated frequency from 0 up to the specified bin

s/cumulated/cumulative/

> > + */
> > +uint64_t Histogram::cumulativeFreq(double bin) const
> > +{
> > +     if (bin <= 0)
> > +             return 0;
> > +     else if (bin >= bins())
> > +             return total();
> > +     int b = static_cast<int32_t>(bin);
> > +     return cumulative_[b] +
> > +            (bin - b) * (cumulative_[b + 1] - cumulative_[b]);
> > +}
> > +
> > +/**
> > + * \brief Return the (fractional) bin of the point through the histogram
> > + * \param[in] q the desired point (0 <= q <= 1)
> > + * \param[in] first low limit (default is 0)
> > + * \param[in] last high limit (default is UINT_MAX)
> > + *
> > + * A quantile gives us the point p = Q(q) in the range such that a proportion
> > + * q of the pixels lie below p. A familiar quantile is Q(0.5) which is the median
> > + * of a distribution.
> > + *
> > + * \return The fractional bin of the point
> > + */
> > +double Histogram::quantile(double q, uint32_t first, uint32_t last) const
> > +{
> > +     if (last == UINT_MAX)
> > +             last = cumulative_.size() - 2;
> > +     ASSERT(first <= last);
> > +
> > +     uint64_t item = q * total();
> > +     /* Binary search to find the right bin */
> > +     while (first < last) {
> > +             int middle = (first + last) / 2;
> > +             /* Is it between first and middle ? */
> > +             if (cumulative_[middle + 1] > item)
> > +                     last = middle;
> > +             else
> > +                     first = middle + 1;
> > +     }
> > +     ASSERT(item >= cumulative_[first] && item <= cumulative_[last + 1]);
> > +
> > +     double frac;
> > +     if (cumulative_[first + 1] == cumulative_[first])
> > +             frac = 0;
> > +     else
> > +             frac = (item - cumulative_[first]) / (cumulative_[first + 1] - cumulative_[first]);
> > +     return first + frac;
> > +}
> > +
> > +/**
> > + * \brief Calculate the mean between two quantiles
> > + * \param[in] lowQuantile low Quantile
> > + * \param[in] highQuantile high Quantile
> > + *
> > + * Quantiles are not ideal for metering as they suffer several limitations.
> > + * Instead, a concept is introduced here: inter-quantile mean.
> > + * It returns the mean of all pixels between lowQuantile and highQuantile.
> > + *
> > + * \return The mean histogram bin value between the two quantiles
> > + */
> > +double Histogram::interQuantileMean(double lowQuantile, double highQuantile) const
> > +{
> > +     ASSERT(highQuantile > lowQuantile);
> > +     /* Proportion of pixels which lies below lowQuantile */
> > +     double lowPoint = quantile(lowQuantile);
> > +     /* Proportion of pixels which lies below highQuantile */
> > +     double highPoint = quantile(highQuantile, static_cast<uint32_t>(lowPoint));
> > +     double sumBinFreq = 0, cumulFreq = 0;
> > +
> > +     for (double p_next = floor(lowPoint) + 1.0;
> > +          p_next <= ceil(highPoint);
> > +          lowPoint = p_next, p_next += 1.0) {
> > +             int bin = floor(lowPoint);
> > +             double freq = (cumulative_[bin + 1] - cumulative_[bin]) * (std::min(p_next, highPoint) - lowPoint);
> > +
> > +             /* Accumulate weigthed bin */
> > +             sumBinFreq += bin * freq;
> > +             /* Accumulate weights */
> > +             cumulFreq += freq;
> > +     }
> > +     /* add 0.5 to give an average for bin mid-points */
> > +     return sumBinFreq / cumulFreq + 0.5;
> > +}

I always meant to rewrite this so that the highPoint is found during
the loop, not by a separate call to quantile(). Job for another day...

> > +
> > +} /* namespace ipa */
> > +
> > +} /* namespace libcamera */
> > diff --git a/src/ipa/libipa/histogram.h b/src/ipa/libipa/histogram.h
> > new file mode 100644
> > index 00000000..dc7451aa
> > --- /dev/null
> > +++ b/src/ipa/libipa/histogram.h
> > @@ -0,0 +1,40 @@
> > +/* SPDX-License-Identifier: BSD-2-Clause */
> > +/*
> > + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
> > + *
> > + * histogram.h - histogram calculation interface
> > + */
> > +#ifndef __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
> > +#define __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
> > +
> > +#include <assert.h>
> > +#include <limits.h>
> > +#include <stdint.h>
> > +
> > +#include <vector>
> > +
> > +#include <libcamera/span.h>
> > +
> > +namespace libcamera {
> > +
> > +namespace ipa {
> > +
> > +class Histogram
> > +{
> > +public:
> > +     Histogram(Span<uint32_t> data);
> > +     size_t bins() const { return cumulative_.size() - 1; }
> > +     uint64_t total() const { return cumulative_[cumulative_.size() - 1]; }
> > +     uint64_t cumulativeFreq(double bin) const;
>
> I'd make this cumulativeFrequency()
>
> > +     double quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
> > +     double interQuantileMean(double lowQuantile, double hiQuantile) const;
> > +
> > +private:
> > +     std::vector<uint64_t> cumulative_;
>
> I see this came from the RPi code.
> Do we really store 64 bit values in here? or are they only 32bit? (or
> smaller?)
>
> In fact, I see we now construct with a Span<uint32_t> data, so
> presumably this vector can be uint32_t.
>
> We could template it to accept different sizes if that would help ...
> but I think supporting 32 bits is probably fine if that's the span we
> currently define....

I guess 32 bits (rather than 64-bit values) is OK, though it does
limit you to images no larger than 4GP (Gigapixels)! Whilst I don't
suppose that matters for any current hardware, you might want to
consider what future hardware might come along (there are Gigapixel
images out there...). The constructor was chosen to work with the
statistics that come out of the Pi ISP, which gives you 32-bit bins.

Best regards

David

>
> > +};
> > +
> > +} /* namespace ipa */
> > +
> > +} /* namespace libcamera */
> > +
> > +#endif /* __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__ */
> > diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
> > index 1819711d..038fc490 100644
> > --- a/src/ipa/libipa/meson.build
> > +++ b/src/ipa/libipa/meson.build
> > @@ -2,10 +2,12 @@
> >
> >  libipa_headers = files([
> >      'algorithm.h',
> > +    'histogram.h'
> >  ])
> >
> >  libipa_sources = files([
> >      'algorithm.cpp',
> > +    'histogram.cpp'
> >  ])
> >
> >  libipa_includes = include_directories('..')
> > diff --git a/src/ipa/raspberrypi/controller/histogram.hpp b/src/ipa/raspberrypi/controller/histogram.hpp
> > index 90f5ac78..fc236416 100644
> > --- a/src/ipa/raspberrypi/controller/histogram.hpp
> > +++ b/src/ipa/raspberrypi/controller/histogram.hpp
>
> I'm not sure if the modifications to the RPi histogram below should be
> in this patch.
>
> It looks like they're unrelated to the actual code move?
>
> Perhaps leave this out, and we can 'convert' RPi controller to use the
> 'generic' histogram separately?
>
>
> > @@ -10,9 +10,6 @@
> >  #include <vector>
> >  #include <cassert>
> >
> > -// A simple histogram class, for use in particular to find "quantiles" and
> > -// averages between "quantiles".
> > -
> >  namespace RPiController {
> >
> >  class Histogram
> > @@ -29,12 +26,8 @@ public:
> >       }
> >       uint32_t Bins() const { return cumulative_.size() - 1; }
> >       uint64_t Total() const { return cumulative_[cumulative_.size() - 1]; }
> > -     // Cumulative frequency up to a (fractional) point in a bin.
> >       uint64_t CumulativeFreq(double bin) const;
> > -     // Return the (fractional) bin of the point q (0 <= q <= 1) through the
> > -     // histogram. Optionally provide limits to help.
> > -     double Quantile(double q, int first = -1, int last = -1) const;
> > -     // Return the average histogram bin value between the two quantiles.
> > +     Histogram::quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
> >       double InterQuantileMean(double q_lo, double q_hi) const;
> >
> >  private:
> >
>
> --
> Regards
> --
> Kieran
> _______________________________________________
> libcamera-devel mailing list
> libcamera-devel@lists.libcamera.org
> https://lists.libcamera.org/listinfo/libcamera-devel
Kieran Bingham April 12, 2021, 6:53 p.m. UTC | #3
Hi David,

On 12/04/2021 14:29, David Plowman wrote:
> Hi Jean-Michel, Kieran
> 
> Happy to see this code getting some use elsewhere! I thought I could
> just add a couple of comments on some of our original thinking...
> 
> On Mon, 12 Apr 2021 at 13:21, Kieran Bingham
> <kieran.bingham@ideasonboard.com> wrote:
>>
>> Hi Jean-Michel,
>>
>> in $SUBJECT
>>
>> s/an/a/
>>
>>
>>
>> On 30/03/2021 22:12, Jean-Michel Hautbois wrote:
>>> This class will be used at least by AGC algorithm when quantiles are
>>> needed for example.
>>>
>>
>> Is this really a Histogram class? Or is is a CumulativeFrequency class?
>>
>> I may likely be mis-interpretting or just wrong - but I thought a
>> histogram stored values like:
>>
>>
>> value 0 1 2 3 4 5 6 7 8 9
>> count 5 0 3 5 6 9 2 0 2 5
>>
>> Meaning that ... you can look up in the table to see that there are 5
>> occurrences of the value 9 in the dataset... whereas this stores:
>>
>>
>> value 0 1 2  3  4  5  6  7  8  9
>> count 5 8 8 13 19 28 30 30 32 37
>>
>>
>> Ok, so now I've drawn that out, I can see that the same information is
>> still there, as you can get the occurrences of any specific value by
>> subtracting from the previous 'bin'?
> 
> I think you've pretty much convinced yourself why it stores cumulative
> frequency. Going from cumulative frequency back to per-bin values is a
> single subtraction. Going the other way is a loop. Given that finding
> "quantiles" is a fairly useful operation, cumulative frequencies make
> sense.

Picture and a thousand words ...

As soon as I drew out the tables I could see it... ;-)

> 
>>
>>
>>> Signed-off-by: Jean-Michel Hautbois <jeanmichel.hautbois@ideasonboard.com>
>>> ---
>>>  src/ipa/libipa/histogram.cpp                 | 154 +++++++++++++++++++
>>>  src/ipa/libipa/histogram.h                   |  40 +++++
>>>  src/ipa/libipa/meson.build                   |   2 +
>>>  src/ipa/raspberrypi/controller/histogram.hpp |   9 +-
>>>  4 files changed, 197 insertions(+), 8 deletions(-)
>>>  create mode 100644 src/ipa/libipa/histogram.cpp
>>>  create mode 100644 src/ipa/libipa/histogram.h
> 
> Will this header be includable by application code? I've sometimes
> found myself getting an image and then wanting to do histogram
> operations. (Of course, this class doesn't make the Histogram from
> scratch, you have to start by totting up the bins yourself, but that's
> not difficult and then you get means, quantiles and stuff for free...)

libipa is not currently exposed as as a public API.

But there's lots of code that I'd like to be more re-usable.

The question is how do we make all the nice helpers available without
committing to never changing them as a public interface ...

I know Laurent has previously discussed moving helpers to separate
library sections or such.

We need to work towards being able to expose a stable API/ABI for
libcamera - but perhaps we could expose 'helper' libraries which don't
need to offer that guarantee?


>>>
>>> diff --git a/src/ipa/libipa/histogram.cpp b/src/ipa/libipa/histogram.cpp
>>> new file mode 100644
>>> index 00000000..46fbb940
>>> --- /dev/null
>>> +++ b/src/ipa/libipa/histogram.cpp
>>> @@ -0,0 +1,154 @@
>>> +/* SPDX-License-Identifier: BSD-2-Clause */
>>> +/*
>>> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
>>> + *
>>> + * histogram.cpp - histogram calculations
>>> + */
>>> +#include "histogram.h"
>>> +
>>> +#include <cmath>
>>> +
>>> +#include "libcamera/internal/log.h"
>>> +
>>> +/**
>>> + * \file histogram.h
>>> + * \brief Class to represent Histograms and manipulate them
>>> + */
>>> +
>>> +namespace libcamera {
>>> +
>>> +namespace ipa {
>>> +
>>> +/**
>>> + * \class Histogram
>>> + * \brief The base class for creating histograms
>>> + *
>>> + * The Histogram class defines a standard interface for IPA algorithms. By
>>> + * abstracting histograms, it makes it possible to implement generic code
>>> + * to manage histograms regardless of their specific type.
>>
>> I don't think this defines a standard interface for IPA algorithms ;)
>>
>> I think we can drop that paragraph and keep the one below.
>>
>>> + *
>>> + * This class stores a cumulative frequency histogram, which is a mapping that
>>> + * counts the cumulative number of observations in all of the bins up to the
>>> + * specified bin. It can be used to find quantiles and averages between quantiles.
>>> + */
>>> +
>>> +/**
>>> + * \brief Create a cumulative histogram with a bin number of intervals
>>
>> I suspect 'bin' was a previous value passed in here.
>>
>> I think we should describe how data should be passed in.
>> I.e. I think perhaps it should be a pre-sorted histogram or something?
>>
>>
>>> + * \param[in] data a reference to the histogram
>>
>> I don't think this is a reference anymore.
>>
>>> + */
>>> +Histogram::Histogram(Span<uint32_t> data)
>>> +{
>>> +     cumulative_.reserve(data.size());
>>> +     cumulative_.push_back(0);
>>> +     for (const uint32_t &value : data)
>>> +             cumulative_.push_back(cumulative_.back() + value);
>>> +}
>>> +/**
>>> + * \fn Histogram::bins()
>>> + * \brief getter for number of bins
>>
>> We don't normally reference them as 'getter's even if that's what they do.
>>
>> To be consistent with other briefs, we should put something like:
>>
>>   \brief Retrieve the number of bins currently stored in the Histogram
>>
>>
>> (stored by, might be 'used by'?)
>>
>>
>>> + * \return Number of bins
>>> + */
>>> +/**
>>> + * \fn Histogram::total()
>>> + * \brief getter for number of values
>>
>> Same here,
>>
>>  \brief Retrieve the total number of values in the data set
>>
>>> + * \return Number of values
>>> + */
>>> +
>>> +/**
>>> + * \brief Cumulative frequency up to a (fractional) point in a bin.
>>> + * \param[in] bin the bin upon which to cumulate
>>
>> s/the/The/
>> (Capital letter to start after the parameter itself.)
>>
>> Reading below, does it also mean this should say "up to which" rather
>> than "upon which"?
>>
>> i.e.
>>    upon which: to me, counts the values 'in' that bin.
>>    up to which: to me, counts the values up to that one...
> 
> Also, I don't think there's a verb "to cumulate", is there?

I wondered if it should say to 'accumulate' ... but then I doubted
myself ...

>>> + *
>>> + * With F(p) the cumulative frequency of the histogram, the value is 0 at
>>
>> What is F(p) here?
> 
> I think this is defining it to be the cumulative frequency.
> 
>>
>>
>>
>>> + * the bottom of the histogram, and the maximum is the number of bins.
>>> + * The pixels are spread evenly throughout the “bin” in which they lie, so that
>>> + * F(p) is a continuous (monotonically increasing) function.
>>> + *
>>> + * \return The cumulated frequency from 0 up to the specified bin
> 
> s/cumulated/cumulative/
> 
>>> + */
>>> +uint64_t Histogram::cumulativeFreq(double bin) const
>>> +{
>>> +     if (bin <= 0)
>>> +             return 0;
>>> +     else if (bin >= bins())
>>> +             return total();
>>> +     int b = static_cast<int32_t>(bin);
>>> +     return cumulative_[b] +
>>> +            (bin - b) * (cumulative_[b + 1] - cumulative_[b]);
>>> +}
>>> +
>>> +/**
>>> + * \brief Return the (fractional) bin of the point through the histogram
>>> + * \param[in] q the desired point (0 <= q <= 1)
>>> + * \param[in] first low limit (default is 0)
>>> + * \param[in] last high limit (default is UINT_MAX)
>>> + *
>>> + * A quantile gives us the point p = Q(q) in the range such that a proportion
>>> + * q of the pixels lie below p. A familiar quantile is Q(0.5) which is the median
>>> + * of a distribution.
>>> + *
>>> + * \return The fractional bin of the point
>>> + */
>>> +double Histogram::quantile(double q, uint32_t first, uint32_t last) const
>>> +{
>>> +     if (last == UINT_MAX)
>>> +             last = cumulative_.size() - 2;
>>> +     ASSERT(first <= last);
>>> +
>>> +     uint64_t item = q * total();
>>> +     /* Binary search to find the right bin */
>>> +     while (first < last) {
>>> +             int middle = (first + last) / 2;
>>> +             /* Is it between first and middle ? */
>>> +             if (cumulative_[middle + 1] > item)
>>> +                     last = middle;
>>> +             else
>>> +                     first = middle + 1;
>>> +     }
>>> +     ASSERT(item >= cumulative_[first] && item <= cumulative_[last + 1]);
>>> +
>>> +     double frac;
>>> +     if (cumulative_[first + 1] == cumulative_[first])
>>> +             frac = 0;
>>> +     else
>>> +             frac = (item - cumulative_[first]) / (cumulative_[first + 1] - cumulative_[first]);
>>> +     return first + frac;
>>> +}
>>> +
>>> +/**
>>> + * \brief Calculate the mean between two quantiles
>>> + * \param[in] lowQuantile low Quantile
>>> + * \param[in] highQuantile high Quantile
>>> + *
>>> + * Quantiles are not ideal for metering as they suffer several limitations.
>>> + * Instead, a concept is introduced here: inter-quantile mean.
>>> + * It returns the mean of all pixels between lowQuantile and highQuantile.
>>> + *
>>> + * \return The mean histogram bin value between the two quantiles
>>> + */
>>> +double Histogram::interQuantileMean(double lowQuantile, double highQuantile) const
>>> +{
>>> +     ASSERT(highQuantile > lowQuantile);
>>> +     /* Proportion of pixels which lies below lowQuantile */
>>> +     double lowPoint = quantile(lowQuantile);
>>> +     /* Proportion of pixels which lies below highQuantile */
>>> +     double highPoint = quantile(highQuantile, static_cast<uint32_t>(lowPoint));
>>> +     double sumBinFreq = 0, cumulFreq = 0;
>>> +
>>> +     for (double p_next = floor(lowPoint) + 1.0;
>>> +          p_next <= ceil(highPoint);
>>> +          lowPoint = p_next, p_next += 1.0) {
>>> +             int bin = floor(lowPoint);
>>> +             double freq = (cumulative_[bin + 1] - cumulative_[bin]) * (std::min(p_next, highPoint) - lowPoint);
>>> +
>>> +             /* Accumulate weigthed bin */
>>> +             sumBinFreq += bin * freq;
>>> +             /* Accumulate weights */
>>> +             cumulFreq += freq;
>>> +     }
>>> +     /* add 0.5 to give an average for bin mid-points */
>>> +     return sumBinFreq / cumulFreq + 0.5;
>>> +}
> 
> I always meant to rewrite this so that the highPoint is found during
> the loop, not by a separate call to quantile(). Job for another day...
> 
>>> +
>>> +} /* namespace ipa */
>>> +
>>> +} /* namespace libcamera */
>>> diff --git a/src/ipa/libipa/histogram.h b/src/ipa/libipa/histogram.h
>>> new file mode 100644
>>> index 00000000..dc7451aa
>>> --- /dev/null
>>> +++ b/src/ipa/libipa/histogram.h
>>> @@ -0,0 +1,40 @@
>>> +/* SPDX-License-Identifier: BSD-2-Clause */
>>> +/*
>>> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
>>> + *
>>> + * histogram.h - histogram calculation interface
>>> + */
>>> +#ifndef __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
>>> +#define __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
>>> +
>>> +#include <assert.h>
>>> +#include <limits.h>
>>> +#include <stdint.h>
>>> +
>>> +#include <vector>
>>> +
>>> +#include <libcamera/span.h>
>>> +
>>> +namespace libcamera {
>>> +
>>> +namespace ipa {
>>> +
>>> +class Histogram
>>> +{
>>> +public:
>>> +     Histogram(Span<uint32_t> data);
>>> +     size_t bins() const { return cumulative_.size() - 1; }
>>> +     uint64_t total() const { return cumulative_[cumulative_.size() - 1]; }
>>> +     uint64_t cumulativeFreq(double bin) const;
>>
>> I'd make this cumulativeFrequency()
>>
>>> +     double quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
>>> +     double interQuantileMean(double lowQuantile, double hiQuantile) const;
>>> +
>>> +private:
>>> +     std::vector<uint64_t> cumulative_;
>>
>> I see this came from the RPi code.
>> Do we really store 64 bit values in here? or are they only 32bit? (or
>> smaller?)
>>
>> In fact, I see we now construct with a Span<uint32_t> data, so
>> presumably this vector can be uint32_t.
>>
>> We could template it to accept different sizes if that would help ...
>> but I think supporting 32 bits is probably fine if that's the span we
>> currently define....
> 
> I guess 32 bits (rather than 64-bit values) is OK, though it does
> limit you to images no larger than 4GP (Gigapixels)! Whilst I don't
> suppose that matters for any current hardware, you might want to
> consider what future hardware might come along (there are Gigapixel
> images out there...). The constructor was chosen to work with the
> statistics that come out of the Pi ISP, which gives you 32-bit bins.

Aha, I saw that the class also supported templates, but that didn't map
to the underlying storage... I had wondered if the template type should
determine that to be able to match the storage type to the type passed
in at construction...



> Best regards
> 
> David
> 
>>
>>> +};
>>> +
>>> +} /* namespace ipa */
>>> +
>>> +} /* namespace libcamera */
>>> +
>>> +#endif /* __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__ */
>>> diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
>>> index 1819711d..038fc490 100644
>>> --- a/src/ipa/libipa/meson.build
>>> +++ b/src/ipa/libipa/meson.build
>>> @@ -2,10 +2,12 @@
>>>
>>>  libipa_headers = files([
>>>      'algorithm.h',
>>> +    'histogram.h'
>>>  ])
>>>
>>>  libipa_sources = files([
>>>      'algorithm.cpp',
>>> +    'histogram.cpp'
>>>  ])
>>>
>>>  libipa_includes = include_directories('..')
>>> diff --git a/src/ipa/raspberrypi/controller/histogram.hpp b/src/ipa/raspberrypi/controller/histogram.hpp
>>> index 90f5ac78..fc236416 100644
>>> --- a/src/ipa/raspberrypi/controller/histogram.hpp
>>> +++ b/src/ipa/raspberrypi/controller/histogram.hpp
>>
>> I'm not sure if the modifications to the RPi histogram below should be
>> in this patch.
>>
>> It looks like they're unrelated to the actual code move?
>>
>> Perhaps leave this out, and we can 'convert' RPi controller to use the
>> 'generic' histogram separately?
>>
>>
>>> @@ -10,9 +10,6 @@
>>>  #include <vector>
>>>  #include <cassert>
>>>
>>> -// A simple histogram class, for use in particular to find "quantiles" and
>>> -// averages between "quantiles".
>>> -
>>>  namespace RPiController {
>>>
>>>  class Histogram
>>> @@ -29,12 +26,8 @@ public:
>>>       }
>>>       uint32_t Bins() const { return cumulative_.size() - 1; }
>>>       uint64_t Total() const { return cumulative_[cumulative_.size() - 1]; }
>>> -     // Cumulative frequency up to a (fractional) point in a bin.
>>>       uint64_t CumulativeFreq(double bin) const;
>>> -     // Return the (fractional) bin of the point q (0 <= q <= 1) through the
>>> -     // histogram. Optionally provide limits to help.
>>> -     double Quantile(double q, int first = -1, int last = -1) const;
>>> -     // Return the average histogram bin value between the two quantiles.
>>> +     Histogram::quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
>>>       double InterQuantileMean(double q_lo, double q_hi) const;
>>>
>>>  private:
>>>
>>
>> --
>> Regards
>> --
>> Kieran
>> _______________________________________________
>> libcamera-devel mailing list
>> libcamera-devel@lists.libcamera.org
>> https://lists.libcamera.org/listinfo/libcamera-devel
Jean-Michel Hautbois April 13, 2021, 5:42 a.m. UTC | #4
Hi David, Kieran,

Thanks for your comments !

On 12/04/2021 20:53, Kieran Bingham wrote:
> Hi David,
> 
> On 12/04/2021 14:29, David Plowman wrote:
>> Hi Jean-Michel, Kieran
>>
>> Happy to see this code getting some use elsewhere! I thought I could
>> just add a couple of comments on some of our original thinking...

Indeed that is your work as a start !

>> On Mon, 12 Apr 2021 at 13:21, Kieran Bingham
>> <kieran.bingham@ideasonboard.com> wrote:
>>>
>>> Hi Jean-Michel,
>>>
>>> in $SUBJECT
>>>
>>> s/an/a/
>>>
>>>
>>>
>>> On 30/03/2021 22:12, Jean-Michel Hautbois wrote:
>>>> This class will be used at least by AGC algorithm when quantiles are
>>>> needed for example.
>>>>
>>>
>>> Is this really a Histogram class? Or is is a CumulativeFrequency class?
>>>
>>> I may likely be mis-interpretting or just wrong - but I thought a
>>> histogram stored values like:
>>>
>>>
>>> value 0 1 2 3 4 5 6 7 8 9
>>> count 5 0 3 5 6 9 2 0 2 5
>>>
>>> Meaning that ... you can look up in the table to see that there are 5
>>> occurrences of the value 9 in the dataset... whereas this stores:
>>>
>>>
>>> value 0 1 2  3  4  5  6  7  8  9
>>> count 5 8 8 13 19 28 30 30 32 37
>>>
>>>
>>> Ok, so now I've drawn that out, I can see that the same information is
>>> still there, as you can get the occurrences of any specific value by
>>> subtracting from the previous 'bin'?
>>
>> I think you've pretty much convinced yourself why it stores cumulative
>> frequency. Going from cumulative frequency back to per-bin values is a
>> single subtraction. Going the other way is a loop. Given that finding
>> "quantiles" is a fairly useful operation, cumulative frequencies make
>> sense.
> 
> Picture and a thousand words ...
> 
> As soon as I drew out the tables I could see it... ;-)

OK, as Kieran wondered and had to draw the tables, I may use your
comment David in the commit message ?

"This class stores a cumulative frequency histogram. Going from
cumulative frequency back to per-bin values is a single subtraction,
while going the other way is a loop."

>>
>>>
>>>
>>>> Signed-off-by: Jean-Michel Hautbois <jeanmichel.hautbois@ideasonboard.com>
>>>> ---
>>>>  src/ipa/libipa/histogram.cpp                 | 154 +++++++++++++++++++
>>>>  src/ipa/libipa/histogram.h                   |  40 +++++
>>>>  src/ipa/libipa/meson.build                   |   2 +
>>>>  src/ipa/raspberrypi/controller/histogram.hpp |   9 +-
>>>>  4 files changed, 197 insertions(+), 8 deletions(-)
>>>>  create mode 100644 src/ipa/libipa/histogram.cpp
>>>>  create mode 100644 src/ipa/libipa/histogram.h
>>
>> Will this header be includable by application code? I've sometimes
>> found myself getting an image and then wanting to do histogram
>> operations. (Of course, this class doesn't make the Histogram from
>> scratch, you have to start by totting up the bins yourself, but that's
>> not difficult and then you get means, quantiles and stuff for free...)
> 
> libipa is not currently exposed as as a public API.
> 
> But there's lots of code that I'd like to be more re-usable.
> 
> The question is how do we make all the nice helpers available without
> committing to never changing them as a public interface ...
> 
> I know Laurent has previously discussed moving helpers to separate
> library sections or such.
> 
> We need to work towards being able to expose a stable API/ABI for
> libcamera - but perhaps we could expose 'helper' libraries which don't
> need to offer that guarantee?
> 

I wish I could display the histogram in qcam too, for example.
But I suppose this class would only be a helper for a higher level one
in Qt ?

>>>>
>>>> diff --git a/src/ipa/libipa/histogram.cpp b/src/ipa/libipa/histogram.cpp
>>>> new file mode 100644
>>>> index 00000000..46fbb940
>>>> --- /dev/null
>>>> +++ b/src/ipa/libipa/histogram.cpp
>>>> @@ -0,0 +1,154 @@
>>>> +/* SPDX-License-Identifier: BSD-2-Clause */
>>>> +/*
>>>> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
>>>> + *
>>>> + * histogram.cpp - histogram calculations
>>>> + */
>>>> +#include "histogram.h"
>>>> +
>>>> +#include <cmath>
>>>> +
>>>> +#include "libcamera/internal/log.h"
>>>> +
>>>> +/**
>>>> + * \file histogram.h
>>>> + * \brief Class to represent Histograms and manipulate them
>>>> + */
>>>> +
>>>> +namespace libcamera {
>>>> +
>>>> +namespace ipa {
>>>> +
>>>> +/**
>>>> + * \class Histogram
>>>> + * \brief The base class for creating histograms
>>>> + *
>>>> + * The Histogram class defines a standard interface for IPA algorithms. By
>>>> + * abstracting histograms, it makes it possible to implement generic code
>>>> + * to manage histograms regardless of their specific type.
>>>
>>> I don't think this defines a standard interface for IPA algorithms ;)
>>>
>>> I think we can drop that paragraph and keep the one below.
>>>
>>>> + *
>>>> + * This class stores a cumulative frequency histogram, which is a mapping that
>>>> + * counts the cumulative number of observations in all of the bins up to the
>>>> + * specified bin. It can be used to find quantiles and averages between quantiles.
>>>> + */
>>>> +
>>>> +/**
>>>> + * \brief Create a cumulative histogram with a bin number of intervals
>>>
>>> I suspect 'bin' was a previous value passed in here.
>>>
>>> I think we should describe how data should be passed in.
>>> I.e. I think perhaps it should be a pre-sorted histogram or something?

Mmh, you are right, it should be pre-sorted...

>>>
>>>> + * \param[in] data a reference to the histogram
>>>
>>> I don't think this is a reference anymore.
>>>
>>>> + */
>>>> +Histogram::Histogram(Span<uint32_t> data)
>>>> +{
>>>> +     cumulative_.reserve(data.size());
>>>> +     cumulative_.push_back(0);
>>>> +     for (const uint32_t &value : data)
>>>> +             cumulative_.push_back(cumulative_.back() + value);
>>>> +}
>>>> +/**
>>>> + * \fn Histogram::bins()
>>>> + * \brief getter for number of bins
>>>
>>> We don't normally reference them as 'getter's even if that's what they do.
>>>
>>> To be consistent with other briefs, we should put something like:
>>>
>>>   \brief Retrieve the number of bins currently stored in the Histogram
>>>
>>>
>>> (stored by, might be 'used by'?)
>>>
>>>
>>>> + * \return Number of bins
>>>> + */
>>>> +/**
>>>> + * \fn Histogram::total()
>>>> + * \brief getter for number of values
>>>
>>> Same here,
>>>
>>>  \brief Retrieve the total number of values in the data set
>>>
>>>> + * \return Number of values
>>>> + */
>>>> +
>>>> +/**
>>>> + * \brief Cumulative frequency up to a (fractional) point in a bin.
>>>> + * \param[in] bin the bin upon which to cumulate
>>>
>>> s/the/The/
>>> (Capital letter to start after the parameter itself.)
>>>
>>> Reading below, does it also mean this should say "up to which" rather
>>> than "upon which"?
>>>
>>> i.e.
>>>    upon which: to me, counts the values 'in' that bin.
>>>    up to which: to me, counts the values up to that one...
>>
>> Also, I don't think there's a verb "to cumulate", is there?
> 
> I wondered if it should say to 'accumulate' ... but then I doubted
> myself ...

:-/ Sorry for making you doubting on your vocabulary. I have to say I
like to translate directly from French :-p.

>>>> + *
>>>> + * With F(p) the cumulative frequency of the histogram, the value is 0 at
>>>
>>> What is F(p) here?
>>
>> I think this is defining it to be the cumulative frequency.

Well, I think this is the way to define a function, "With ... a ..." ?
Should I rephrase it ?

>>>
>>>
>>>
>>>> + * the bottom of the histogram, and the maximum is the number of bins.
>>>> + * The pixels are spread evenly throughout the “bin” in which they lie, so that
>>>> + * F(p) is a continuous (monotonically increasing) function.
>>>> + *
>>>> + * \return The cumulated frequency from 0 up to the specified bin
>>
>> s/cumulated/cumulative/
>>
>>>> + */
>>>> +uint64_t Histogram::cumulativeFreq(double bin) const
>>>> +{
>>>> +     if (bin <= 0)
>>>> +             return 0;
>>>> +     else if (bin >= bins())
>>>> +             return total();
>>>> +     int b = static_cast<int32_t>(bin);
>>>> +     return cumulative_[b] +
>>>> +            (bin - b) * (cumulative_[b + 1] - cumulative_[b]);
>>>> +}
>>>> +
>>>> +/**
>>>> + * \brief Return the (fractional) bin of the point through the histogram
>>>> + * \param[in] q the desired point (0 <= q <= 1)
>>>> + * \param[in] first low limit (default is 0)
>>>> + * \param[in] last high limit (default is UINT_MAX)
>>>> + *
>>>> + * A quantile gives us the point p = Q(q) in the range such that a proportion
>>>> + * q of the pixels lie below p. A familiar quantile is Q(0.5) which is the median
>>>> + * of a distribution.
>>>> + *
>>>> + * \return The fractional bin of the point
>>>> + */
>>>> +double Histogram::quantile(double q, uint32_t first, uint32_t last) const
>>>> +{
>>>> +     if (last == UINT_MAX)
>>>> +             last = cumulative_.size() - 2;
>>>> +     ASSERT(first <= last);
>>>> +
>>>> +     uint64_t item = q * total();
>>>> +     /* Binary search to find the right bin */
>>>> +     while (first < last) {
>>>> +             int middle = (first + last) / 2;
>>>> +             /* Is it between first and middle ? */
>>>> +             if (cumulative_[middle + 1] > item)
>>>> +                     last = middle;
>>>> +             else
>>>> +                     first = middle + 1;
>>>> +     }
>>>> +     ASSERT(item >= cumulative_[first] && item <= cumulative_[last + 1]);
>>>> +
>>>> +     double frac;
>>>> +     if (cumulative_[first + 1] == cumulative_[first])
>>>> +             frac = 0;
>>>> +     else
>>>> +             frac = (item - cumulative_[first]) / (cumulative_[first + 1] - cumulative_[first]);
>>>> +     return first + frac;
>>>> +}
>>>> +
>>>> +/**
>>>> + * \brief Calculate the mean between two quantiles
>>>> + * \param[in] lowQuantile low Quantile
>>>> + * \param[in] highQuantile high Quantile
>>>> + *
>>>> + * Quantiles are not ideal for metering as they suffer several limitations.
>>>> + * Instead, a concept is introduced here: inter-quantile mean.
>>>> + * It returns the mean of all pixels between lowQuantile and highQuantile.
>>>> + *
>>>> + * \return The mean histogram bin value between the two quantiles
>>>> + */
>>>> +double Histogram::interQuantileMean(double lowQuantile, double highQuantile) const
>>>> +{
>>>> +     ASSERT(highQuantile > lowQuantile);
>>>> +     /* Proportion of pixels which lies below lowQuantile */
>>>> +     double lowPoint = quantile(lowQuantile);
>>>> +     /* Proportion of pixels which lies below highQuantile */
>>>> +     double highPoint = quantile(highQuantile, static_cast<uint32_t>(lowPoint));
>>>> +     double sumBinFreq = 0, cumulFreq = 0;
>>>> +
>>>> +     for (double p_next = floor(lowPoint) + 1.0;
>>>> +          p_next <= ceil(highPoint);
>>>> +          lowPoint = p_next, p_next += 1.0) {
>>>> +             int bin = floor(lowPoint);
>>>> +             double freq = (cumulative_[bin + 1] - cumulative_[bin]) * (std::min(p_next, highPoint) - lowPoint);
>>>> +
>>>> +             /* Accumulate weigthed bin */
>>>> +             sumBinFreq += bin * freq;
>>>> +             /* Accumulate weights */
>>>> +             cumulFreq += freq;
>>>> +     }
>>>> +     /* add 0.5 to give an average for bin mid-points */
>>>> +     return sumBinFreq / cumulFreq + 0.5;
>>>> +}
>>
>> I always meant to rewrite this so that the highPoint is found during
>> the loop, not by a separate call to quantile(). Job for another day...
>>
>>>> +
>>>> +} /* namespace ipa */
>>>> +
>>>> +} /* namespace libcamera */
>>>> diff --git a/src/ipa/libipa/histogram.h b/src/ipa/libipa/histogram.h
>>>> new file mode 100644
>>>> index 00000000..dc7451aa
>>>> --- /dev/null
>>>> +++ b/src/ipa/libipa/histogram.h
>>>> @@ -0,0 +1,40 @@
>>>> +/* SPDX-License-Identifier: BSD-2-Clause */
>>>> +/*
>>>> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
>>>> + *
>>>> + * histogram.h - histogram calculation interface
>>>> + */
>>>> +#ifndef __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
>>>> +#define __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
>>>> +
>>>> +#include <assert.h>
>>>> +#include <limits.h>
>>>> +#include <stdint.h>
>>>> +
>>>> +#include <vector>
>>>> +
>>>> +#include <libcamera/span.h>
>>>> +
>>>> +namespace libcamera {
>>>> +
>>>> +namespace ipa {
>>>> +
>>>> +class Histogram
>>>> +{
>>>> +public:
>>>> +     Histogram(Span<uint32_t> data);
>>>> +     size_t bins() const { return cumulative_.size() - 1; }
>>>> +     uint64_t total() const { return cumulative_[cumulative_.size() - 1]; }
>>>> +     uint64_t cumulativeFreq(double bin) const;
>>>
>>> I'd make this cumulativeFrequency()
>>>
>>>> +     double quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
>>>> +     double interQuantileMean(double lowQuantile, double hiQuantile) const;
>>>> +
>>>> +private:
>>>> +     std::vector<uint64_t> cumulative_;
>>>
>>> I see this came from the RPi code.
>>> Do we really store 64 bit values in here? or are they only 32bit? (or
>>> smaller?)
>>>
>>> In fact, I see we now construct with a Span<uint32_t> data, so
>>> presumably this vector can be uint32_t.
>>>
>>> We could template it to accept different sizes if that would help ...
>>> but I think supporting 32 bits is probably fine if that's the span we
>>> currently define....
>>
>> I guess 32 bits (rather than 64-bit values) is OK, though it does
>> limit you to images no larger than 4GP (Gigapixels)! Whilst I don't
>> suppose that matters for any current hardware, you might want to
>> consider what future hardware might come along (there are Gigapixel
>> images out there...). The constructor was chosen to work with the
>> statistics that come out of the Pi ISP, which gives you 32-bit bins.
> 
> Aha, I saw that the class also supported templates, but that didn't map
> to the underlying storage... I had wondered if the template type should
> determine that to be able to match the storage type to the type passed
> in at construction...

I started with the template type, but could not find a good enough
reason not to specify a type.
Rockchip ISP also stores u32 for its histograms. I don't think lots of
ISPs are using u64 for the bins ?

> 
> 
>> Best regards
>>
>> David
>>
>>>
>>>> +};
>>>> +
>>>> +} /* namespace ipa */
>>>> +
>>>> +} /* namespace libcamera */
>>>> +
>>>> +#endif /* __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__ */
>>>> diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
>>>> index 1819711d..038fc490 100644
>>>> --- a/src/ipa/libipa/meson.build
>>>> +++ b/src/ipa/libipa/meson.build
>>>> @@ -2,10 +2,12 @@
>>>>
>>>>  libipa_headers = files([
>>>>      'algorithm.h',
>>>> +    'histogram.h'
>>>>  ])
>>>>
>>>>  libipa_sources = files([
>>>>      'algorithm.cpp',
>>>> +    'histogram.cpp'
>>>>  ])
>>>>
>>>>  libipa_includes = include_directories('..')
>>>> diff --git a/src/ipa/raspberrypi/controller/histogram.hpp b/src/ipa/raspberrypi/controller/histogram.hpp
>>>> index 90f5ac78..fc236416 100644
>>>> --- a/src/ipa/raspberrypi/controller/histogram.hpp
>>>> +++ b/src/ipa/raspberrypi/controller/histogram.hpp
>>>
>>> I'm not sure if the modifications to the RPi histogram below should be
>>> in this patch.
>>>
>>> It looks like they're unrelated to the actual code move?
>>>
>>> Perhaps leave this out, and we can 'convert' RPi controller to use the
>>> 'generic' histogram separately?
>>>
>>>
>>>> @@ -10,9 +10,6 @@
>>>>  #include <vector>
>>>>  #include <cassert>
>>>>
>>>> -// A simple histogram class, for use in particular to find "quantiles" and
>>>> -// averages between "quantiles".
>>>> -
>>>>  namespace RPiController {
>>>>
>>>>  class Histogram
>>>> @@ -29,12 +26,8 @@ public:
>>>>       }
>>>>       uint32_t Bins() const { return cumulative_.size() - 1; }
>>>>       uint64_t Total() const { return cumulative_[cumulative_.size() - 1]; }
>>>> -     // Cumulative frequency up to a (fractional) point in a bin.
>>>>       uint64_t CumulativeFreq(double bin) const;
>>>> -     // Return the (fractional) bin of the point q (0 <= q <= 1) through the
>>>> -     // histogram. Optionally provide limits to help.
>>>> -     double Quantile(double q, int first = -1, int last = -1) const;
>>>> -     // Return the average histogram bin value between the two quantiles.
>>>> +     Histogram::quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
>>>>       double InterQuantileMean(double q_lo, double q_hi) const;
>>>>
>>>>  private:
>>>>
>>>
>>> --
>>> Regards
>>> --
>>> Kieran
>>> _______________________________________________
>>> libcamera-devel mailing list
>>> libcamera-devel@lists.libcamera.org
>>> https://lists.libcamera.org/listinfo/libcamera-devel
>
Jean-Michel Hautbois April 13, 2021, 5:45 a.m. UTC | #5
On 13/04/2021 07:42, Jean-Michel Hautbois wrote:
> Hi David, Kieran,
> 
> Thanks for your comments !
> 
> On 12/04/2021 20:53, Kieran Bingham wrote:
>> Hi David,
>>
>> On 12/04/2021 14:29, David Plowman wrote:
>>> Hi Jean-Michel, Kieran
>>>
>>> Happy to see this code getting some use elsewhere! I thought I could
>>> just add a couple of comments on some of our original thinking...
> 
> Indeed that is your work as a start !
> 
>>> On Mon, 12 Apr 2021 at 13:21, Kieran Bingham
>>> <kieran.bingham@ideasonboard.com> wrote:
>>>>
>>>> Hi Jean-Michel,
>>>>
>>>> in $SUBJECT
>>>>
>>>> s/an/a/
>>>>
>>>>
>>>>
>>>> On 30/03/2021 22:12, Jean-Michel Hautbois wrote:
>>>>> This class will be used at least by AGC algorithm when quantiles are
>>>>> needed for example.
>>>>>
>>>>
>>>> Is this really a Histogram class? Or is is a CumulativeFrequency class?
>>>>
>>>> I may likely be mis-interpretting or just wrong - but I thought a
>>>> histogram stored values like:
>>>>
>>>>
>>>> value 0 1 2 3 4 5 6 7 8 9
>>>> count 5 0 3 5 6 9 2 0 2 5
>>>>
>>>> Meaning that ... you can look up in the table to see that there are 5
>>>> occurrences of the value 9 in the dataset... whereas this stores:
>>>>
>>>>
>>>> value 0 1 2  3  4  5  6  7  8  9
>>>> count 5 8 8 13 19 28 30 30 32 37
>>>>
>>>>
>>>> Ok, so now I've drawn that out, I can see that the same information is
>>>> still there, as you can get the occurrences of any specific value by
>>>> subtracting from the previous 'bin'?
>>>
>>> I think you've pretty much convinced yourself why it stores cumulative
>>> frequency. Going from cumulative frequency back to per-bin values is a
>>> single subtraction. Going the other way is a loop. Given that finding
>>> "quantiles" is a fairly useful operation, cumulative frequencies make
>>> sense.
>>
>> Picture and a thousand words ...
>>
>> As soon as I drew out the tables I could see it... ;-)
> 
> OK, as Kieran wondered and had to draw the tables, I may use your
> comment David in the commit message ?
> 
> "This class stores a cumulative frequency histogram. Going from
> cumulative frequency back to per-bin values is a single subtraction,
> while going the other way is a loop."
> 
>>>
>>>>
>>>>
>>>>> Signed-off-by: Jean-Michel Hautbois <jeanmichel.hautbois@ideasonboard.com>
>>>>> ---
>>>>>  src/ipa/libipa/histogram.cpp                 | 154 +++++++++++++++++++
>>>>>  src/ipa/libipa/histogram.h                   |  40 +++++
>>>>>  src/ipa/libipa/meson.build                   |   2 +
>>>>>  src/ipa/raspberrypi/controller/histogram.hpp |   9 +-
>>>>>  4 files changed, 197 insertions(+), 8 deletions(-)
>>>>>  create mode 100644 src/ipa/libipa/histogram.cpp
>>>>>  create mode 100644 src/ipa/libipa/histogram.h
>>>
>>> Will this header be includable by application code? I've sometimes
>>> found myself getting an image and then wanting to do histogram
>>> operations. (Of course, this class doesn't make the Histogram from
>>> scratch, you have to start by totting up the bins yourself, but that's
>>> not difficult and then you get means, quantiles and stuff for free...)
>>
>> libipa is not currently exposed as as a public API.
>>
>> But there's lots of code that I'd like to be more re-usable.
>>
>> The question is how do we make all the nice helpers available without
>> committing to never changing them as a public interface ...
>>
>> I know Laurent has previously discussed moving helpers to separate
>> library sections or such.
>>
>> We need to work towards being able to expose a stable API/ABI for
>> libcamera - but perhaps we could expose 'helper' libraries which don't
>> need to offer that guarantee?
>>
> 
> I wish I could display the histogram in qcam too, for example.
> But I suppose this class would only be a helper for a higher level one
> in Qt ?
> 
>>>>>
>>>>> diff --git a/src/ipa/libipa/histogram.cpp b/src/ipa/libipa/histogram.cpp
>>>>> new file mode 100644
>>>>> index 00000000..46fbb940
>>>>> --- /dev/null
>>>>> +++ b/src/ipa/libipa/histogram.cpp
>>>>> @@ -0,0 +1,154 @@
>>>>> +/* SPDX-License-Identifier: BSD-2-Clause */
>>>>> +/*
>>>>> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
>>>>> + *
>>>>> + * histogram.cpp - histogram calculations
>>>>> + */
>>>>> +#include "histogram.h"
>>>>> +
>>>>> +#include <cmath>
>>>>> +
>>>>> +#include "libcamera/internal/log.h"
>>>>> +
>>>>> +/**
>>>>> + * \file histogram.h
>>>>> + * \brief Class to represent Histograms and manipulate them
>>>>> + */
>>>>> +
>>>>> +namespace libcamera {
>>>>> +
>>>>> +namespace ipa {
>>>>> +
>>>>> +/**
>>>>> + * \class Histogram
>>>>> + * \brief The base class for creating histograms
>>>>> + *
>>>>> + * The Histogram class defines a standard interface for IPA algorithms. By
>>>>> + * abstracting histograms, it makes it possible to implement generic code
>>>>> + * to manage histograms regardless of their specific type.
>>>>
>>>> I don't think this defines a standard interface for IPA algorithms ;)
>>>>
>>>> I think we can drop that paragraph and keep the one below.
>>>>
>>>>> + *
>>>>> + * This class stores a cumulative frequency histogram, which is a mapping that
>>>>> + * counts the cumulative number of observations in all of the bins up to the
>>>>> + * specified bin. It can be used to find quantiles and averages between quantiles.
>>>>> + */
>>>>> +
>>>>> +/**
>>>>> + * \brief Create a cumulative histogram with a bin number of intervals
>>>>
>>>> I suspect 'bin' was a previous value passed in here.
>>>>
>>>> I think we should describe how data should be passed in.
>>>> I.e. I think perhaps it should be a pre-sorted histogram or something?
> 
> Mmh, you are right, it should be pre-sorted...
> 
>>>>
>>>>> + * \param[in] data a reference to the histogram
>>>>
>>>> I don't think this is a reference anymore.
>>>>
>>>>> + */
>>>>> +Histogram::Histogram(Span<uint32_t> data)
>>>>> +{
>>>>> +     cumulative_.reserve(data.size());
>>>>> +     cumulative_.push_back(0);
>>>>> +     for (const uint32_t &value : data)
>>>>> +             cumulative_.push_back(cumulative_.back() + value);
>>>>> +}
>>>>> +/**
>>>>> + * \fn Histogram::bins()
>>>>> + * \brief getter for number of bins
>>>>
>>>> We don't normally reference them as 'getter's even if that's what they do.
>>>>
>>>> To be consistent with other briefs, we should put something like:
>>>>
>>>>   \brief Retrieve the number of bins currently stored in the Histogram
>>>>
>>>>
>>>> (stored by, might be 'used by'?)
>>>>
>>>>
>>>>> + * \return Number of bins
>>>>> + */
>>>>> +/**
>>>>> + * \fn Histogram::total()
>>>>> + * \brief getter for number of values
>>>>
>>>> Same here,
>>>>
>>>>  \brief Retrieve the total number of values in the data set
>>>>
>>>>> + * \return Number of values
>>>>> + */
>>>>> +
>>>>> +/**
>>>>> + * \brief Cumulative frequency up to a (fractional) point in a bin.
>>>>> + * \param[in] bin the bin upon which to cumulate
>>>>
>>>> s/the/The/
>>>> (Capital letter to start after the parameter itself.)
>>>>
>>>> Reading below, does it also mean this should say "up to which" rather
>>>> than "upon which"?
>>>>
>>>> i.e.
>>>>    upon which: to me, counts the values 'in' that bin.
>>>>    up to which: to me, counts the values up to that one...
>>>
>>> Also, I don't think there's a verb "to cumulate", is there?
>>
>> I wondered if it should say to 'accumulate' ... but then I doubted
>> myself ...
> 
> :-/ Sorry for making you doubting on your vocabulary. I have to say I
> like to translate directly from French :-p.
> 
>>>>> + *
>>>>> + * With F(p) the cumulative frequency of the histogram, the value is 0 at
>>>>
>>>> What is F(p) here?
>>>
>>> I think this is defining it to be the cumulative frequency.
> 
> Well, I think this is the way to define a function, "With ... a ..." ?
> Should I rephrase it ?
> 
>>>>
>>>>
>>>>
>>>>> + * the bottom of the histogram, and the maximum is the number of bins.
>>>>> + * The pixels are spread evenly throughout the “bin” in which they lie, so that
>>>>> + * F(p) is a continuous (monotonically increasing) function.
>>>>> + *
>>>>> + * \return The cumulated frequency from 0 up to the specified bin
>>>
>>> s/cumulated/cumulative/
>>>
>>>>> + */
>>>>> +uint64_t Histogram::cumulativeFreq(double bin) const
>>>>> +{
>>>>> +     if (bin <= 0)
>>>>> +             return 0;
>>>>> +     else if (bin >= bins())
>>>>> +             return total();
>>>>> +     int b = static_cast<int32_t>(bin);
>>>>> +     return cumulative_[b] +
>>>>> +            (bin - b) * (cumulative_[b + 1] - cumulative_[b]);
>>>>> +}
>>>>> +
>>>>> +/**
>>>>> + * \brief Return the (fractional) bin of the point through the histogram
>>>>> + * \param[in] q the desired point (0 <= q <= 1)
>>>>> + * \param[in] first low limit (default is 0)
>>>>> + * \param[in] last high limit (default is UINT_MAX)
>>>>> + *
>>>>> + * A quantile gives us the point p = Q(q) in the range such that a proportion
>>>>> + * q of the pixels lie below p. A familiar quantile is Q(0.5) which is the median
>>>>> + * of a distribution.
>>>>> + *
>>>>> + * \return The fractional bin of the point
>>>>> + */
>>>>> +double Histogram::quantile(double q, uint32_t first, uint32_t last) const
>>>>> +{
>>>>> +     if (last == UINT_MAX)
>>>>> +             last = cumulative_.size() - 2;
>>>>> +     ASSERT(first <= last);
>>>>> +
>>>>> +     uint64_t item = q * total();
>>>>> +     /* Binary search to find the right bin */
>>>>> +     while (first < last) {
>>>>> +             int middle = (first + last) / 2;
>>>>> +             /* Is it between first and middle ? */
>>>>> +             if (cumulative_[middle + 1] > item)
>>>>> +                     last = middle;
>>>>> +             else
>>>>> +                     first = middle + 1;
>>>>> +     }
>>>>> +     ASSERT(item >= cumulative_[first] && item <= cumulative_[last + 1]);
>>>>> +
>>>>> +     double frac;
>>>>> +     if (cumulative_[first + 1] == cumulative_[first])
>>>>> +             frac = 0;
>>>>> +     else
>>>>> +             frac = (item - cumulative_[first]) / (cumulative_[first + 1] - cumulative_[first]);
>>>>> +     return first + frac;
>>>>> +}
>>>>> +
>>>>> +/**
>>>>> + * \brief Calculate the mean between two quantiles
>>>>> + * \param[in] lowQuantile low Quantile
>>>>> + * \param[in] highQuantile high Quantile
>>>>> + *
>>>>> + * Quantiles are not ideal for metering as they suffer several limitations.
>>>>> + * Instead, a concept is introduced here: inter-quantile mean.
>>>>> + * It returns the mean of all pixels between lowQuantile and highQuantile.
>>>>> + *
>>>>> + * \return The mean histogram bin value between the two quantiles
>>>>> + */
>>>>> +double Histogram::interQuantileMean(double lowQuantile, double highQuantile) const
>>>>> +{
>>>>> +     ASSERT(highQuantile > lowQuantile);
>>>>> +     /* Proportion of pixels which lies below lowQuantile */
>>>>> +     double lowPoint = quantile(lowQuantile);
>>>>> +     /* Proportion of pixels which lies below highQuantile */
>>>>> +     double highPoint = quantile(highQuantile, static_cast<uint32_t>(lowPoint));
>>>>> +     double sumBinFreq = 0, cumulFreq = 0;
>>>>> +
>>>>> +     for (double p_next = floor(lowPoint) + 1.0;
>>>>> +          p_next <= ceil(highPoint);
>>>>> +          lowPoint = p_next, p_next += 1.0) {
>>>>> +             int bin = floor(lowPoint);
>>>>> +             double freq = (cumulative_[bin + 1] - cumulative_[bin]) * (std::min(p_next, highPoint) - lowPoint);
>>>>> +
>>>>> +             /* Accumulate weigthed bin */
>>>>> +             sumBinFreq += bin * freq;
>>>>> +             /* Accumulate weights */
>>>>> +             cumulFreq += freq;
>>>>> +     }
>>>>> +     /* add 0.5 to give an average for bin mid-points */
>>>>> +     return sumBinFreq / cumulFreq + 0.5;
>>>>> +}
>>>
>>> I always meant to rewrite this so that the highPoint is found during
>>> the loop, not by a separate call to quantile(). Job for another day...
>>>
>>>>> +
>>>>> +} /* namespace ipa */
>>>>> +
>>>>> +} /* namespace libcamera */
>>>>> diff --git a/src/ipa/libipa/histogram.h b/src/ipa/libipa/histogram.h
>>>>> new file mode 100644
>>>>> index 00000000..dc7451aa
>>>>> --- /dev/null
>>>>> +++ b/src/ipa/libipa/histogram.h
>>>>> @@ -0,0 +1,40 @@
>>>>> +/* SPDX-License-Identifier: BSD-2-Clause */
>>>>> +/*
>>>>> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
>>>>> + *
>>>>> + * histogram.h - histogram calculation interface
>>>>> + */
>>>>> +#ifndef __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
>>>>> +#define __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
>>>>> +
>>>>> +#include <assert.h>
>>>>> +#include <limits.h>
>>>>> +#include <stdint.h>
>>>>> +
>>>>> +#include <vector>
>>>>> +
>>>>> +#include <libcamera/span.h>
>>>>> +
>>>>> +namespace libcamera {
>>>>> +
>>>>> +namespace ipa {
>>>>> +
>>>>> +class Histogram
>>>>> +{
>>>>> +public:
>>>>> +     Histogram(Span<uint32_t> data);
>>>>> +     size_t bins() const { return cumulative_.size() - 1; }
>>>>> +     uint64_t total() const { return cumulative_[cumulative_.size() - 1]; }
>>>>> +     uint64_t cumulativeFreq(double bin) const;
>>>>
>>>> I'd make this cumulativeFrequency()
>>>>
>>>>> +     double quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
>>>>> +     double interQuantileMean(double lowQuantile, double hiQuantile) const;
>>>>> +
>>>>> +private:
>>>>> +     std::vector<uint64_t> cumulative_;
>>>>
>>>> I see this came from the RPi code.
>>>> Do we really store 64 bit values in here? or are they only 32bit? (or
>>>> smaller?)
>>>>
>>>> In fact, I see we now construct with a Span<uint32_t> data, so
>>>> presumably this vector can be uint32_t.
>>>>
>>>> We could template it to accept different sizes if that would help ...
>>>> but I think supporting 32 bits is probably fine if that's the span we
>>>> currently define....
>>>
>>> I guess 32 bits (rather than 64-bit values) is OK, though it does
>>> limit you to images no larger than 4GP (Gigapixels)! Whilst I don't
>>> suppose that matters for any current hardware, you might want to
>>> consider what future hardware might come along (there are Gigapixel
>>> images out there...). The constructor was chosen to work with the
>>> statistics that come out of the Pi ISP, which gives you 32-bit bins.
>>
>> Aha, I saw that the class also supported templates, but that didn't map
>> to the underlying storage... I had wondered if the template type should
>> determine that to be able to match the storage type to the type passed
>> in at construction...
> 
> I started with the template type, but could not find a good enough
> reason not to specify a type.
> Rockchip ISP also stores u32 for its histograms. I don't think lots of
> ISPs are using u64 for the bins ?
> 
>>
>>
>>> Best regards
>>>
>>> David
>>>
>>>>
>>>>> +};
>>>>> +
>>>>> +} /* namespace ipa */
>>>>> +
>>>>> +} /* namespace libcamera */
>>>>> +
>>>>> +#endif /* __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__ */
>>>>> diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
>>>>> index 1819711d..038fc490 100644
>>>>> --- a/src/ipa/libipa/meson.build
>>>>> +++ b/src/ipa/libipa/meson.build
>>>>> @@ -2,10 +2,12 @@
>>>>>
>>>>>  libipa_headers = files([
>>>>>      'algorithm.h',
>>>>> +    'histogram.h'
>>>>>  ])
>>>>>
>>>>>  libipa_sources = files([
>>>>>      'algorithm.cpp',
>>>>> +    'histogram.cpp'
>>>>>  ])
>>>>>
>>>>>  libipa_includes = include_directories('..')
>>>>> diff --git a/src/ipa/raspberrypi/controller/histogram.hpp b/src/ipa/raspberrypi/controller/histogram.hpp
>>>>> index 90f5ac78..fc236416 100644
>>>>> --- a/src/ipa/raspberrypi/controller/histogram.hpp
>>>>> +++ b/src/ipa/raspberrypi/controller/histogram.hpp
>>>>
>>>> I'm not sure if the modifications to the RPi histogram below should be
>>>> in this patch.
>>>>
>>>> It looks like they're unrelated to the actual code move?
>>>>
>>>> Perhaps leave this out, and we can 'convert' RPi controller to use the
>>>> 'generic' histogram separately?
>>>>
>>>>
>>>>> @@ -10,9 +10,6 @@
>>>>>  #include <vector>
>>>>>  #include <cassert>
>>>>>
>>>>> -// A simple histogram class, for use in particular to find "quantiles" and
>>>>> -// averages between "quantiles".
>>>>> -
>>>>>  namespace RPiController {
>>>>>
>>>>>  class Histogram
>>>>> @@ -29,12 +26,8 @@ public:
>>>>>       }
>>>>>       uint32_t Bins() const { return cumulative_.size() - 1; }
>>>>>       uint64_t Total() const { return cumulative_[cumulative_.size() - 1]; }
>>>>> -     // Cumulative frequency up to a (fractional) point in a bin.
>>>>>       uint64_t CumulativeFreq(double bin) const;
>>>>> -     // Return the (fractional) bin of the point q (0 <= q <= 1) through the
>>>>> -     // histogram. Optionally provide limits to help.
>>>>> -     double Quantile(double q, int first = -1, int last = -1) const;
>>>>> -     // Return the average histogram bin value between the two quantiles.
>>>>> +     Histogram::quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
>>>>>       double InterQuantileMean(double q_lo, double q_hi) const;

I just saw that I also modified the histogram.hpp file in RPi !
That was not intended at first... and it obviously should be removed
from the patch !

>>>>>  private:
>>>>>
>>>>
>>>> --
>>>> Regards
>>>> --
>>>> Kieran
>>>> _______________________________________________
>>>> libcamera-devel mailing list
>>>> libcamera-devel@lists.libcamera.org
>>>> https://lists.libcamera.org/listinfo/libcamera-devel
>>
Kieran Bingham April 13, 2021, 9:48 a.m. UTC | #6
Hi JM,

On 13/04/2021 06:42, Jean-Michel Hautbois wrote:
> Hi David, Kieran,
> 
> Thanks for your comments !
> 
> On 12/04/2021 20:53, Kieran Bingham wrote:
>> Hi David,
>>
>> On 12/04/2021 14:29, David Plowman wrote:
>>> Hi Jean-Michel, Kieran
>>>
>>> Happy to see this code getting some use elsewhere! I thought I could
>>> just add a couple of comments on some of our original thinking...
> 
> Indeed that is your work as a start !
> 
>>> On Mon, 12 Apr 2021 at 13:21, Kieran Bingham
>>> <kieran.bingham@ideasonboard.com> wrote:
>>>>
>>>> Hi Jean-Michel,
>>>>
>>>> in $SUBJECT
>>>>
>>>> s/an/a/
>>>>
>>>>
>>>>
>>>> On 30/03/2021 22:12, Jean-Michel Hautbois wrote:
>>>>> This class will be used at least by AGC algorithm when quantiles are
>>>>> needed for example.
>>>>>
>>>>
>>>> Is this really a Histogram class? Or is is a CumulativeFrequency class?
>>>>
>>>> I may likely be mis-interpretting or just wrong - but I thought a
>>>> histogram stored values like:
>>>>
>>>>
>>>> value 0 1 2 3 4 5 6 7 8 9
>>>> count 5 0 3 5 6 9 2 0 2 5
>>>>
>>>> Meaning that ... you can look up in the table to see that there are 5
>>>> occurrences of the value 9 in the dataset... whereas this stores:
>>>>
>>>>
>>>> value 0 1 2  3  4  5  6  7  8  9
>>>> count 5 8 8 13 19 28 30 30 32 37
>>>>
>>>>
>>>> Ok, so now I've drawn that out, I can see that the same information is
>>>> still there, as you can get the occurrences of any specific value by
>>>> subtracting from the previous 'bin'?
>>>
>>> I think you've pretty much convinced yourself why it stores cumulative
>>> frequency. Going from cumulative frequency back to per-bin values is a
>>> single subtraction. Going the other way is a loop. Given that finding
>>> "quantiles" is a fairly useful operation, cumulative frequencies make
>>> sense.
>>
>> Picture and a thousand words ...
>>
>> As soon as I drew out the tables I could see it... ;-)
> 
> OK, as Kieran wondered and had to draw the tables, I may use your
> comment David in the commit message ?
> 
> "This class stores a cumulative frequency histogram. Going from
> cumulative frequency back to per-bin values is a single subtraction,
> while going the other way is a loop."

That sounds good but I'd leave out the last
  ", while going the other way is a loop."



> 
>>>
>>>>
>>>>
>>>>> Signed-off-by: Jean-Michel Hautbois <jeanmichel.hautbois@ideasonboard.com>
>>>>> ---
>>>>>  src/ipa/libipa/histogram.cpp                 | 154 +++++++++++++++++++
>>>>>  src/ipa/libipa/histogram.h                   |  40 +++++
>>>>>  src/ipa/libipa/meson.build                   |   2 +
>>>>>  src/ipa/raspberrypi/controller/histogram.hpp |   9 +-
>>>>>  4 files changed, 197 insertions(+), 8 deletions(-)
>>>>>  create mode 100644 src/ipa/libipa/histogram.cpp
>>>>>  create mode 100644 src/ipa/libipa/histogram.h
>>>
>>> Will this header be includable by application code? I've sometimes
>>> found myself getting an image and then wanting to do histogram
>>> operations. (Of course, this class doesn't make the Histogram from
>>> scratch, you have to start by totting up the bins yourself, but that's
>>> not difficult and then you get means, quantiles and stuff for free...)
>>
>> libipa is not currently exposed as as a public API.
>>
>> But there's lots of code that I'd like to be more re-usable.
>>
>> The question is how do we make all the nice helpers available without
>> committing to never changing them as a public interface ...
>>
>> I know Laurent has previously discussed moving helpers to separate
>> library sections or such.
>>
>> We need to work towards being able to expose a stable API/ABI for
>> libcamera - but perhaps we could expose 'helper' libraries which don't
>> need to offer that guarantee?
> 
> I wish I could display the histogram in qcam too, for example.
> But I suppose this class would only be a helper for a higher level one
> in Qt ?

I think I asked about storing histograms in image metadata to send to
applications when available as I think they're useful data (and if
provided by an ISP, that's better than post-processing/counting in
software) - but Laurent stated it would be hard to define the type and
what it would represent. Plus it may be a lot of data to copy if it's
used as a control.

Any thoughts here? Is it easy to define a histogram for an image as a
metadata control? (i.e. a span...)?

What image metadata would we report to applications 'if we could'?
- Brightness histogram?
- per-channel histogram?


>>>>> diff --git a/src/ipa/libipa/histogram.cpp b/src/ipa/libipa/histogram.cpp
>>>>> new file mode 100644
>>>>> index 00000000..46fbb940
>>>>> --- /dev/null
>>>>> +++ b/src/ipa/libipa/histogram.cpp
>>>>> @@ -0,0 +1,154 @@
>>>>> +/* SPDX-License-Identifier: BSD-2-Clause */
>>>>> +/*
>>>>> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
>>>>> + *
>>>>> + * histogram.cpp - histogram calculations
>>>>> + */
>>>>> +#include "histogram.h"
>>>>> +
>>>>> +#include <cmath>
>>>>> +
>>>>> +#include "libcamera/internal/log.h"
>>>>> +
>>>>> +/**
>>>>> + * \file histogram.h
>>>>> + * \brief Class to represent Histograms and manipulate them
>>>>> + */
>>>>> +
>>>>> +namespace libcamera {
>>>>> +
>>>>> +namespace ipa {
>>>>> +
>>>>> +/**
>>>>> + * \class Histogram
>>>>> + * \brief The base class for creating histograms
>>>>> + *
>>>>> + * The Histogram class defines a standard interface for IPA algorithms. By
>>>>> + * abstracting histograms, it makes it possible to implement generic code
>>>>> + * to manage histograms regardless of their specific type.
>>>>
>>>> I don't think this defines a standard interface for IPA algorithms ;)
>>>>
>>>> I think we can drop that paragraph and keep the one below.
>>>>
>>>>> + *
>>>>> + * This class stores a cumulative frequency histogram, which is a mapping that
>>>>> + * counts the cumulative number of observations in all of the bins up to the
>>>>> + * specified bin. It can be used to find quantiles and averages between quantiles.
>>>>> + */
>>>>> +
>>>>> +/**
>>>>> + * \brief Create a cumulative histogram with a bin number of intervals
>>>>
>>>> I suspect 'bin' was a previous value passed in here.
>>>>
>>>> I think we should describe how data should be passed in.
>>>> I.e. I think perhaps it should be a pre-sorted histogram or something?
> 
> Mmh, you are right, it should be pre-sorted...
> 
>>>>
>>>>> + * \param[in] data a reference to the histogram
>>>>
>>>> I don't think this is a reference anymore.
>>>>
>>>>> + */
>>>>> +Histogram::Histogram(Span<uint32_t> data)
>>>>> +{
>>>>> +     cumulative_.reserve(data.size());
>>>>> +     cumulative_.push_back(0);
>>>>> +     for (const uint32_t &value : data)
>>>>> +             cumulative_.push_back(cumulative_.back() + value);
>>>>> +}
>>>>> +/**
>>>>> + * \fn Histogram::bins()
>>>>> + * \brief getter for number of bins
>>>>
>>>> We don't normally reference them as 'getter's even if that's what they do.
>>>>
>>>> To be consistent with other briefs, we should put something like:
>>>>
>>>>   \brief Retrieve the number of bins currently stored in the Histogram
>>>>
>>>>
>>>> (stored by, might be 'used by'?)
>>>>
>>>>
>>>>> + * \return Number of bins
>>>>> + */
>>>>> +/**
>>>>> + * \fn Histogram::total()
>>>>> + * \brief getter for number of values
>>>>
>>>> Same here,
>>>>
>>>>  \brief Retrieve the total number of values in the data set
>>>>
>>>>> + * \return Number of values
>>>>> + */
>>>>> +
>>>>> +/**
>>>>> + * \brief Cumulative frequency up to a (fractional) point in a bin.
>>>>> + * \param[in] bin the bin upon which to cumulate
>>>>
>>>> s/the/The/
>>>> (Capital letter to start after the parameter itself.)
>>>>
>>>> Reading below, does it also mean this should say "up to which" rather
>>>> than "upon which"?
>>>>
>>>> i.e.
>>>>    upon which: to me, counts the values 'in' that bin.
>>>>    up to which: to me, counts the values up to that one...
>>>
>>> Also, I don't think there's a verb "to cumulate", is there?
>>
>> I wondered if it should say to 'accumulate' ... but then I doubted
>> myself ...
> 
> :-/ Sorry for making you doubting on your vocabulary. I have to say I
> like to translate directly from French :-p.
> 
>>>>> + *
>>>>> + * With F(p) the cumulative frequency of the histogram, the value is 0 at
>>>>
>>>> What is F(p) here?
>>>
>>> I think this is defining it to be the cumulative frequency.
> 
> Well, I think this is the way to define a function, "With ... a ..." ?
> Should I rephrase it ?

I guess we're mixing maths and programming terms ... given that it's a
mathematical function ... maybe that's appropriate ;-)



>>>>> + * the bottom of the histogram, and the maximum is the number of bins.
>>>>> + * The pixels are spread evenly throughout the “bin” in which they lie, so that
>>>>> + * F(p) is a continuous (monotonically increasing) function.
>>>>> + *
>>>>> + * \return The cumulated frequency from 0 up to the specified bin

If I try to rewrite this in a way that I would understand better I get
the following: (That doesn't mean what I've written is correct though,
so take it with a pinch of salt)

"""
Determine the cumulative frequency of the histogram up to the desired
fractional bin position.

Bins are assumed to contain an evenly spread set of values, giving a
continuous and monotonically increasing frequency throughout the bin.

The histogram data is stored internally as a cumulative set to optimise
this procedure.
"""


>>>
>>> s/cumulated/cumulative/
>>>
>>>>> + */
>>>>> +uint64_t Histogram::cumulativeFreq(double bin) const
>>>>> +{
>>>>> +     if (bin <= 0)
>>>>> +             return 0;
>>>>> +     else if (bin >= bins())
>>>>> +             return total();
>>>>> +     int b = static_cast<int32_t>(bin);
>>>>> +     return cumulative_[b] +
>>>>> +            (bin - b) * (cumulative_[b + 1] - cumulative_[b]);
>>>>> +}
>>>>> +
>>>>> +/**
>>>>> + * \brief Return the (fractional) bin of the point through the histogram
>>>>> + * \param[in] q the desired point (0 <= q <= 1)
>>>>> + * \param[in] first low limit (default is 0)
>>>>> + * \param[in] last high limit (default is UINT_MAX)
>>>>> + *
>>>>> + * A quantile gives us the point p = Q(q) in the range such that a proportion
>>>>> + * q of the pixels lie below p. A familiar quantile is Q(0.5) which is the median
>>>>> + * of a distribution.
>>>>> + *
>>>>> + * \return The fractional bin of the point
>>>>> + */
>>>>> +double Histogram::quantile(double q, uint32_t first, uint32_t last) const
>>>>> +{
>>>>> +     if (last == UINT_MAX)
>>>>> +             last = cumulative_.size() - 2;
>>>>> +     ASSERT(first <= last);
>>>>> +
>>>>> +     uint64_t item = q * total();
>>>>> +     /* Binary search to find the right bin */
>>>>> +     while (first < last) {
>>>>> +             int middle = (first + last) / 2;
>>>>> +             /* Is it between first and middle ? */
>>>>> +             if (cumulative_[middle + 1] > item)
>>>>> +                     last = middle;
>>>>> +             else
>>>>> +                     first = middle + 1;
>>>>> +     }
>>>>> +     ASSERT(item >= cumulative_[first] && item <= cumulative_[last + 1]);
>>>>> +
>>>>> +     double frac;
>>>>> +     if (cumulative_[first + 1] == cumulative_[first])
>>>>> +             frac = 0;
>>>>> +     else
>>>>> +             frac = (item - cumulative_[first]) / (cumulative_[first + 1] - cumulative_[first]);
>>>>> +     return first + frac;
>>>>> +}
>>>>> +
>>>>> +/**
>>>>> + * \brief Calculate the mean between two quantiles
>>>>> + * \param[in] lowQuantile low Quantile
>>>>> + * \param[in] highQuantile high Quantile
>>>>> + *
>>>>> + * Quantiles are not ideal for metering as they suffer several limitations.
>>>>> + * Instead, a concept is introduced here: inter-quantile mean.
>>>>> + * It returns the mean of all pixels between lowQuantile and highQuantile.
>>>>> + *
>>>>> + * \return The mean histogram bin value between the two quantiles
>>>>> + */
>>>>> +double Histogram::interQuantileMean(double lowQuantile, double highQuantile) const
>>>>> +{
>>>>> +     ASSERT(highQuantile > lowQuantile);
>>>>> +     /* Proportion of pixels which lies below lowQuantile */
>>>>> +     double lowPoint = quantile(lowQuantile);
>>>>> +     /* Proportion of pixels which lies below highQuantile */
>>>>> +     double highPoint = quantile(highQuantile, static_cast<uint32_t>(lowPoint));
>>>>> +     double sumBinFreq = 0, cumulFreq = 0;
>>>>> +
>>>>> +     for (double p_next = floor(lowPoint) + 1.0;
>>>>> +          p_next <= ceil(highPoint);
>>>>> +          lowPoint = p_next, p_next += 1.0) {
>>>>> +             int bin = floor(lowPoint);
>>>>> +             double freq = (cumulative_[bin + 1] - cumulative_[bin]) * (std::min(p_next, highPoint) - lowPoint);
>>>>> +
>>>>> +             /* Accumulate weigthed bin */

s/weigthed/weighted/


>>>>> +             sumBinFreq += bin * freq;
>>>>> +             /* Accumulate weights */
>>>>> +             cumulFreq += freq;
>>>>> +     }
>>>>> +     /* add 0.5 to give an average for bin mid-points */
>>>>> +     return sumBinFreq / cumulFreq + 0.5;
>>>>> +}
>>>
>>> I always meant to rewrite this so that the highPoint is found during
>>> the loop, not by a separate call to quantile(). Job for another day...
>>>
>>>>> +
>>>>> +} /* namespace ipa */
>>>>> +
>>>>> +} /* namespace libcamera */
>>>>> diff --git a/src/ipa/libipa/histogram.h b/src/ipa/libipa/histogram.h
>>>>> new file mode 100644
>>>>> index 00000000..dc7451aa
>>>>> --- /dev/null
>>>>> +++ b/src/ipa/libipa/histogram.h
>>>>> @@ -0,0 +1,40 @@
>>>>> +/* SPDX-License-Identifier: BSD-2-Clause */
>>>>> +/*
>>>>> + * Copyright (C) 2019, Raspberry Pi (Trading) Limited
>>>>> + *
>>>>> + * histogram.h - histogram calculation interface
>>>>> + */
>>>>> +#ifndef __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
>>>>> +#define __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
>>>>> +
>>>>> +#include <assert.h>
>>>>> +#include <limits.h>
>>>>> +#include <stdint.h>
>>>>> +
>>>>> +#include <vector>
>>>>> +
>>>>> +#include <libcamera/span.h>
>>>>> +
>>>>> +namespace libcamera {
>>>>> +
>>>>> +namespace ipa {
>>>>> +
>>>>> +class Histogram
>>>>> +{
>>>>> +public:
>>>>> +     Histogram(Span<uint32_t> data);
>>>>> +     size_t bins() const { return cumulative_.size() - 1; }
>>>>> +     uint64_t total() const { return cumulative_[cumulative_.size() - 1]; }
>>>>> +     uint64_t cumulativeFreq(double bin) const;
>>>>
>>>> I'd make this cumulativeFrequency()
>>>>
>>>>> +     double quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
>>>>> +     double interQuantileMean(double lowQuantile, double hiQuantile) const;
>>>>> +
>>>>> +private:
>>>>> +     std::vector<uint64_t> cumulative_;
>>>>
>>>> I see this came from the RPi code.
>>>> Do we really store 64 bit values in here? or are they only 32bit? (or
>>>> smaller?)
>>>>
>>>> In fact, I see we now construct with a Span<uint32_t> data, so
>>>> presumably this vector can be uint32_t.
>>>>
>>>> We could template it to accept different sizes if that would help ...
>>>> but I think supporting 32 bits is probably fine if that's the span we
>>>> currently define....
>>>
>>> I guess 32 bits (rather than 64-bit values) is OK, though it does
>>> limit you to images no larger than 4GP (Gigapixels)! Whilst I don't
>>> suppose that matters for any current hardware, you might want to
>>> consider what future hardware might come along (there are Gigapixel
>>> images out there...). The constructor was chosen to work with the
>>> statistics that come out of the Pi ISP, which gives you 32-bit bins.
>>
>> Aha, I saw that the class also supported templates, but that didn't map
>> to the underlying storage... I had wondered if the template type should
>> determine that to be able to match the storage type to the type passed
>> in at construction...
> 
> I started with the template type, but could not find a good enough
> reason not to specify a type.
> Rockchip ISP also stores u32 for its histograms. I don't think lots of
> ISPs are using u64 for the bins ?

In fact, now I re-read this - the cumulative_ is not storing the same
type as the input span<T> is it.

It's changed, as it could be bigger than the input data (as it's additive).

So ... perhaps leaving it as 64 bit might be fine?

I guess we don't normally have huge amounts of bins, so the memory
consumption here is not going to be too exhaustive?


--
Kieran


>>> Best regards
>>>
>>> David
>>>
>>>>
>>>>> +};
>>>>> +
>>>>> +} /* namespace ipa */
>>>>> +
>>>>> +} /* namespace libcamera */
>>>>> +
>>>>> +#endif /* __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__ */
>>>>> diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
>>>>> index 1819711d..038fc490 100644
>>>>> --- a/src/ipa/libipa/meson.build
>>>>> +++ b/src/ipa/libipa/meson.build
>>>>> @@ -2,10 +2,12 @@
>>>>>
>>>>>  libipa_headers = files([
>>>>>      'algorithm.h',
>>>>> +    'histogram.h'
>>>>>  ])
>>>>>
>>>>>  libipa_sources = files([
>>>>>      'algorithm.cpp',
>>>>> +    'histogram.cpp'
>>>>>  ])
>>>>>
>>>>>  libipa_includes = include_directories('..')
>>>>> diff --git a/src/ipa/raspberrypi/controller/histogram.hpp b/src/ipa/raspberrypi/controller/histogram.hpp
>>>>> index 90f5ac78..fc236416 100644
>>>>> --- a/src/ipa/raspberrypi/controller/histogram.hpp
>>>>> +++ b/src/ipa/raspberrypi/controller/histogram.hpp
>>>>
>>>> I'm not sure if the modifications to the RPi histogram below should be
>>>> in this patch.
>>>>
>>>> It looks like they're unrelated to the actual code move?
>>>>
>>>> Perhaps leave this out, and we can 'convert' RPi controller to use the
>>>> 'generic' histogram separately?
>>>>
>>>>
>>>>> @@ -10,9 +10,6 @@
>>>>>  #include <vector>
>>>>>  #include <cassert>
>>>>>
>>>>> -// A simple histogram class, for use in particular to find "quantiles" and
>>>>> -// averages between "quantiles".
>>>>> -
>>>>>  namespace RPiController {
>>>>>
>>>>>  class Histogram
>>>>> @@ -29,12 +26,8 @@ public:
>>>>>       }
>>>>>       uint32_t Bins() const { return cumulative_.size() - 1; }
>>>>>       uint64_t Total() const { return cumulative_[cumulative_.size() - 1]; }
>>>>> -     // Cumulative frequency up to a (fractional) point in a bin.
>>>>>       uint64_t CumulativeFreq(double bin) const;
>>>>> -     // Return the (fractional) bin of the point q (0 <= q <= 1) through the
>>>>> -     // histogram. Optionally provide limits to help.
>>>>> -     double Quantile(double q, int first = -1, int last = -1) const;
>>>>> -     // Return the average histogram bin value between the two quantiles.
>>>>> +     Histogram::quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
>>>>>       double InterQuantileMean(double q_lo, double q_hi) const;
>>>>>
>>>>>  private:
>>>>>
>>>>
>>>> --
>>>> Regards
>>>> --
>>>> Kieran
>>>> _______________________________________________
>>>> libcamera-devel mailing list
>>>> libcamera-devel@lists.libcamera.org
>>>> https://lists.libcamera.org/listinfo/libcamera-devel
>>

Patch
diff mbox series

diff --git a/src/ipa/libipa/histogram.cpp b/src/ipa/libipa/histogram.cpp
new file mode 100644
index 00000000..46fbb940
--- /dev/null
+++ b/src/ipa/libipa/histogram.cpp
@@ -0,0 +1,154 @@ 
+/* SPDX-License-Identifier: BSD-2-Clause */
+/*
+ * Copyright (C) 2019, Raspberry Pi (Trading) Limited
+ *
+ * histogram.cpp - histogram calculations
+ */
+#include "histogram.h"
+
+#include <cmath>
+
+#include "libcamera/internal/log.h"
+
+/**
+ * \file histogram.h
+ * \brief Class to represent Histograms and manipulate them
+ */
+
+namespace libcamera {
+
+namespace ipa {
+
+/**
+ * \class Histogram
+ * \brief The base class for creating histograms
+ *
+ * The Histogram class defines a standard interface for IPA algorithms. By
+ * abstracting histograms, it makes it possible to implement generic code
+ * to manage histograms regardless of their specific type.
+ *
+ * This class stores a cumulative frequency histogram, which is a mapping that
+ * counts the cumulative number of observations in all of the bins up to the
+ * specified bin. It can be used to find quantiles and averages between quantiles.
+ */
+
+/**
+ * \brief Create a cumulative histogram with a bin number of intervals
+ * \param[in] data a reference to the histogram
+ */
+Histogram::Histogram(Span<uint32_t> data)
+{
+	cumulative_.reserve(data.size());
+	cumulative_.push_back(0);
+	for (const uint32_t &value : data)
+		cumulative_.push_back(cumulative_.back() + value);
+}
+/**
+ * \fn Histogram::bins()
+ * \brief getter for number of bins
+ * \return Number of bins
+ */
+/**
+ * \fn Histogram::total()
+ * \brief getter for number of values
+ * \return Number of values
+ */
+
+/**
+ * \brief Cumulative frequency up to a (fractional) point in a bin.
+ * \param[in] bin the bin upon which to cumulate
+ *
+ * With F(p) the cumulative frequency of the histogram, the value is 0 at
+ * the bottom of the histogram, and the maximum is the number of bins.
+ * The pixels are spread evenly throughout the “bin” in which they lie, so that
+ * F(p) is a continuous (monotonically increasing) function.
+ *
+ * \return The cumulated frequency from 0 up to the specified bin
+ */
+uint64_t Histogram::cumulativeFreq(double bin) const
+{
+	if (bin <= 0)
+		return 0;
+	else if (bin >= bins())
+		return total();
+	int b = static_cast<int32_t>(bin);
+	return cumulative_[b] +
+	       (bin - b) * (cumulative_[b + 1] - cumulative_[b]);
+}
+
+/**
+ * \brief Return the (fractional) bin of the point through the histogram
+ * \param[in] q the desired point (0 <= q <= 1)
+ * \param[in] first low limit (default is 0)
+ * \param[in] last high limit (default is UINT_MAX)
+ *
+ * A quantile gives us the point p = Q(q) in the range such that a proportion
+ * q of the pixels lie below p. A familiar quantile is Q(0.5) which is the median
+ * of a distribution.
+ *
+ * \return The fractional bin of the point
+ */
+double Histogram::quantile(double q, uint32_t first, uint32_t last) const
+{
+	if (last == UINT_MAX)
+		last = cumulative_.size() - 2;
+	ASSERT(first <= last);
+
+	uint64_t item = q * total();
+	/* Binary search to find the right bin */
+	while (first < last) {
+		int middle = (first + last) / 2;
+		/* Is it between first and middle ? */
+		if (cumulative_[middle + 1] > item)
+			last = middle;
+		else
+			first = middle + 1;
+	}
+	ASSERT(item >= cumulative_[first] && item <= cumulative_[last + 1]);
+
+	double frac;
+	if (cumulative_[first + 1] == cumulative_[first])
+		frac = 0;
+	else
+		frac = (item - cumulative_[first]) / (cumulative_[first + 1] - cumulative_[first]);
+	return first + frac;
+}
+
+/**
+ * \brief Calculate the mean between two quantiles
+ * \param[in] lowQuantile low Quantile
+ * \param[in] highQuantile high Quantile
+ *
+ * Quantiles are not ideal for metering as they suffer several limitations.
+ * Instead, a concept is introduced here: inter-quantile mean.
+ * It returns the mean of all pixels between lowQuantile and highQuantile.
+ *
+ * \return The mean histogram bin value between the two quantiles
+ */
+double Histogram::interQuantileMean(double lowQuantile, double highQuantile) const
+{
+	ASSERT(highQuantile > lowQuantile);
+	/* Proportion of pixels which lies below lowQuantile */
+	double lowPoint = quantile(lowQuantile);
+	/* Proportion of pixels which lies below highQuantile */
+	double highPoint = quantile(highQuantile, static_cast<uint32_t>(lowPoint));
+	double sumBinFreq = 0, cumulFreq = 0;
+
+	for (double p_next = floor(lowPoint) + 1.0;
+	     p_next <= ceil(highPoint);
+	     lowPoint = p_next, p_next += 1.0) {
+		int bin = floor(lowPoint);
+		double freq = (cumulative_[bin + 1] - cumulative_[bin]) * (std::min(p_next, highPoint) - lowPoint);
+
+		/* Accumulate weigthed bin */
+		sumBinFreq += bin * freq;
+		/* Accumulate weights */
+		cumulFreq += freq;
+	}
+	/* add 0.5 to give an average for bin mid-points */
+	return sumBinFreq / cumulFreq + 0.5;
+}
+
+} /* namespace ipa */
+
+} /* namespace libcamera */
diff --git a/src/ipa/libipa/histogram.h b/src/ipa/libipa/histogram.h
new file mode 100644
index 00000000..dc7451aa
--- /dev/null
+++ b/src/ipa/libipa/histogram.h
@@ -0,0 +1,40 @@ 
+/* SPDX-License-Identifier: BSD-2-Clause */
+/*
+ * Copyright (C) 2019, Raspberry Pi (Trading) Limited
+ *
+ * histogram.h - histogram calculation interface
+ */
+#ifndef __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
+#define __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__
+
+#include <assert.h>
+#include <limits.h>
+#include <stdint.h>
+
+#include <vector>
+
+#include <libcamera/span.h>
+
+namespace libcamera {
+
+namespace ipa {
+
+class Histogram
+{
+public:
+	Histogram(Span<uint32_t> data);
+	size_t bins() const { return cumulative_.size() - 1; }
+	uint64_t total() const { return cumulative_[cumulative_.size() - 1]; }
+	uint64_t cumulativeFreq(double bin) const;
+	double quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
+	double interQuantileMean(double lowQuantile, double hiQuantile) const;
+
+private:
+	std::vector<uint64_t> cumulative_;
+};
+
+} /* namespace ipa */
+
+} /* namespace libcamera */
+
+#endif /* __LIBCAMERA_IPA_LIBIPA_HISTOGRAM_H__ */
diff --git a/src/ipa/libipa/meson.build b/src/ipa/libipa/meson.build
index 1819711d..038fc490 100644
--- a/src/ipa/libipa/meson.build
+++ b/src/ipa/libipa/meson.build
@@ -2,10 +2,12 @@ 
 
 libipa_headers = files([
     'algorithm.h',
+    'histogram.h'
 ])
 
 libipa_sources = files([
     'algorithm.cpp',
+    'histogram.cpp'
 ])
 
 libipa_includes = include_directories('..')
diff --git a/src/ipa/raspberrypi/controller/histogram.hpp b/src/ipa/raspberrypi/controller/histogram.hpp
index 90f5ac78..fc236416 100644
--- a/src/ipa/raspberrypi/controller/histogram.hpp
+++ b/src/ipa/raspberrypi/controller/histogram.hpp
@@ -10,9 +10,6 @@ 
 #include <vector>
 #include <cassert>
 
-// A simple histogram class, for use in particular to find "quantiles" and
-// averages between "quantiles".
-
 namespace RPiController {
 
 class Histogram
@@ -29,12 +26,8 @@  public:
 	}
 	uint32_t Bins() const { return cumulative_.size() - 1; }
 	uint64_t Total() const { return cumulative_[cumulative_.size() - 1]; }
-	// Cumulative frequency up to a (fractional) point in a bin.
 	uint64_t CumulativeFreq(double bin) const;
-	// Return the (fractional) bin of the point q (0 <= q <= 1) through the
-	// histogram. Optionally provide limits to help.
-	double Quantile(double q, int first = -1, int last = -1) const;
-	// Return the average histogram bin value between the two quantiles.
+	Histogram::quantile(double q, uint32_t first = 0, uint32_t last = UINT_MAX) const;
 	double InterQuantileMean(double q_lo, double q_hi) const;
 
 private: