
\section{Introduction}

Histogram-based thresholding is a ubiquitous tool in image processing, medical imaging, and document analysis: The grayscale intensities of an input image are used to compute a histogram, and some algorithm is then applied to that histogram to identify an optimal threshold (corresponding to a bin location along the histogram's $x$-axis) with which the histogram is to be ``split'' into two parts. That threshold is then used to binarize the input image, under the assumption that this binarized image will then be used for some downstream task such as classification or segmentation.
Thresholding has been the focus of a half-century of research, and is well documented in several survey papers \cite{sahoo1988survey,zhang2008image}.
One might assume that a global thresholding approach has little value in a modern machine learning context: preprocessing an image via binarization discards potentially-valuable information in the input image, and training a neural network to perform binarization is a straightforward alternative.
However, training requires significant amounts of training data, and in many contexts (particularly medical imaging) such data may be prohibitively expensive or difficult to obtain. As such, there continues to be value in ``automatic'' thresholding algorithms that do not require training data, as is evidenced by the active use of these algorithms in the medical imaging literature.

\section{Preliminaries}

Given an input image $\mathrm{I}$ (Figure~\ref{fig:teaser1}) that is usually assumed to contain 8-bit or 16-bit quantized grayscale intensity values, a histogram-based automatic thresholding technique first constructs a histogram of that image consisting of a vector of histogram bin locations $\myvec{x}$ and a corresponding vector of histogram counts or weights $\myvec{n}$  (Figure~\ref{fig:teaser2}).
With this histogram, the algorithm attempts to find some threshold $t$ that separates the two halves of the histogram according to some criteria (\textit{e}.\textit{g}.  maximum likelihood, minimal distortion, \textit{etc}.), and this threshold is used to produce a binarized image $\mathrm{I} > t$ (Figures~\ref{fig:teaser3}-\ref{fig:teaser5}).
Many prior works assume that $\at{x}{i} = i$, which addresses the common case in which histogram bins exactly correspond to quantized pixel intensities, but we will make no such assumption.
This allows our algorithm (and our descriptions and implementations of baseline algorithms) to be applied to arbitrary sets of sorted points, and to histograms whose bins have arbitrary locations. For example, it is equivalent to binarize some quantized image $\mathrm{I}$ using a histogram:
\begin{equation}
(\myvec{n}, \myvec{x}) = \operatorname{hist}(\mathrm{I})\,,
\end{equation}
or using a vector of sorted values each with a ``count'' of 1:
\begin{equation}
 \myvec{n} = \vv{1}\,, \quad\quad \myvec{x} = \operatorname{sort}(\mathrm{I})\,.
\label{eq:sorted}
\end{equation}
This equivalence is not used in any result presented here, but is useful when binarizing a set of continuous values.

Most histogram-based thresholding techniques work by considering each possible ``split'' of the histogram: each value in $\myvec{x}$ is considered as a candidate for $t$, and two quantities are produced that reflect the surface statistics of $\at{\myvec{n}}{\myvec{x}{\,\leq\,}t}$ and $\at{\myvec{n}}{\myvec{x}{\,>\,}t}$.
A critical insight of many classic histogram-based thresholding techniques is that some of these quantities can be computed recursively. 
For example, the sum of all histogram counts up through some index $i$ ($\at{w^{(0)}}{i} = \sum_{j=0}^i \at{n}{j}$) need not be recomputed from scratch if that sum has already been previously computed for all the previous histogram element, and can instead just be updated with a single addition ($\at{w^{(0)}}{i} = \at{w^{(0)}}{i\,\mbox{--}1} + \at{n}{i}$).
Here we construct pairs of vectors of intermediate values that measure surface statistics of the histogram above and below each ``split'', to be used by GHT and our baselines:
\begin{align}
\myvec{\intweightscalar}^{(0)} &= \cumsum{\myvec{n}} &
\myvec{\intweightscalar}^{(1)} &= \decsum{\myvec{n}} \label{eq:prelim} \\
\myvec{\mixturescalar}^{(0)} &= \myvec{\intweightscalar}^{(0)} / \sumvec{\nvec} &
\myvec{\mixturescalar}^{(1)} &= \myvec{\intweightscalar}^{(1)} / \sumvec{\nvec} = 1-\myvec{\mixturescalar}^{(0)} \nonumber \\
\myvec{\meanscalar}^{(0)} &= \cumsum{\myvec{n} \myvec{x}} / \myvec{\intweightscalar}^{(0)} \nonumber & 
\myvec{\meanscalar}^{(1)} &= \decsum{\myvec{n} \myvec{x}} / \myvec{\intweightscalar}^{(1)} \nonumber\\
\myvec{d}^{(0)} &= \cumsum{\myvec{n} \myvec{x}^2} - \myvec{\intweightscalar}^{(0)} \left(\myvec{\meanscalar}^{(0)}\right)^2 \nonumber &
\myvec{d}^{(1)} &= \decsum{\myvec{n} \myvec{x}^2} - \myvec{\intweightscalar}^{(1)} \left(\myvec{\meanscalar}^{(1)}\right)^2 \nonumber
\end{align}
All multiplications, divisions, and exponentiations are element-wise, and $\cumsum{}$ and $\decsum{}$ are cumulative and ``reverse cumulative'' sums respectively.
$\sumvec{\nvec}$ is the sum of all histogram counts in $\myvec{n}$, while $\myvec{\intweightscalar}^{(0)}$ and $\myvec{\intweightscalar}^{(1)}$ are the sums of histogram counts in $\myvec{n}$ below and above each split of the histogram, respectively.
Similarly, each $\myvec{\mixturescalar}^{(k)}$ respectively represent normalized histogram counts (or mixture weights) above and below each split.
Each $\myvec{\meanscalar}^{(k)}$ and $\myvec{d}^{(k)}$ represent the mean and distortion of all elements of $\myvec{x}$ below and above each split, weighted by their histogram counts $\myvec{n}$.
Here ``distortion'' is used in the context of k-means, where it refers to the sum of squared distances of a set of points to their mean.
The computation of $\myvec{d}^{(k)}$ follows from three observations:
First, $\cumsum{\cdot}$ or $\decsum{\cdot}$ can be used (by weighting by $\myvec{n}$ and normalizing by $\myvec{\intweightscalar}$) to compute a vector of expected values. Second, the sample variance of a quantity is a function of the expectation of its values and its values squared ($\operatorname {Var} (X) = \operatorname {E}[X^{2}]-\operatorname {E} [X]^{2}$). Third, the sample variance of a set of points is proportional to the total distortion of those points.

Formally, $\at{\cumsum{\cdot}}{i}$ is the inclusive sum of every vector element at or before index $i$, and $\at{\decsum{\cdot}}{i}$ is the exclusive sum of everything after index $i$:
\begin{align}
    \at{\cumsum{\myvec{n}}}{i} &= \sum_{j=0}^{i} \at{n}{j} = \at{\cumsum{\myvec{n}}}{i-1} + \at{n}{i}\,, \\
    \at{\decsum{\myvec{n}}}{i} &= \!\!\!\!\!\sum_{j=i+1}^{\operatorname{len}(\nvec)\mbox{-}1}\!\!\!\! \at{n}{j} = \at{\decsum{\myvec{n}}}{i+1} + \at{n}{i+1}\,. \nonumber 
\end{align}
Note that the outputs of these cumulative and decremental sums have a length that is one less than that of $\myvec{x}$ and $\myvec{n}$, as these values measure ``splits'' of the data, and we only consider splits that contain at least one histogram element.
The subscripts in our vectors of statistics correspond to a threshold \emph{after} each index, which means that binarization should be performed by a greater-than comparison with the returned threshold ($\myvec{x} > t$).
In practice, $\decsum{\myvec{n}}$ can be computed using just the cumulative and total sums of $\myvec{n}$:
\begin{align}
    \decsum{\myvec{n}} &= \sumvec{\nvec} - \cumsum{\myvec{n}}\,.
\end{align}
As such, these quantities (and all other vector quantities that will be described later) can be computed efficiently with just one or two passes over the histogram.

\section{Algorithm}

Our Generalized Histogram Thresholding (GHT) algorithm is motivated by a straightforward Bayesian treatment of histogram thresholding.
We assume that each pixel in the input image is generated from a mixture of two probability distributions corresponding to intensities below and above some threshold, and we maximize the joint probability of all pixels (which are assumed to be i.i.d.):
\begin{gather}
    \prod_{x,y} \sum_{k=1}^2 p\Big(\mathrm{I}_{x,y} \, \Big| \, z_{x,y}=k \Big) p(z_{x,y}=k) \,.
    \label{eq:bayesianprob}
\end{gather}
The probability of each pixel intensity $\mathrm{I}_{x,y}$ conditioned on its (unknown) mixture component $z_{x,y}$ is:
\begin{equation}
\resizebox{\textwidth}{!}{$
    \displaystyle
    p\Big(\mathrm{I}_{x,y} \, \Big| \, z_{x,y}=k \Big) = \pi^{(k)} \mathcal{N}\!\left( \mathrm{I}_{x,y}  {\,\mid\,}  \mu^{(k)}, \sigma^{(k)} \right) 
    \sichisq\!\left(\sigma^{(k)}  {\,\mid\,} \pi^{(k)} \nu, \tau^2 \right)
     \operatorname{Beta}\!\left( \pi^{(k)}  {\,\mid\,}  \kappa, \omega \right)\,.
     \label{eq:probform}
     $}
\end{equation}
This is similar to a straightforward mixture of Gaussians: we have a mixture probability ($\pi^{(0)}$, $\pi^{(1)} = 1 - \pi^{(0)}$), two means ($\mu^{(0)}$, $\mu^{(1)}$), and two standard deviations ($\sigma^{(0)}$, $\sigma^{(1)}$).
But unlike a standard mixture of Gaussians, we place conjugate priors on model parameters $\pi^{(k)}$ and $\sigma^{(k)}$. We assume each standard deviation is the mode (\textit{i}.\textit{e}.,  the ``updated'' $\tau$ value) of a scaled inverse chi-squared distribution, which is the conjugate prior of a Gaussian with a known mean.
\begin{equation}
\sichisq(x {\,\mid\,} \nu, \tau^2) = {\frac {(\tau ^{2}\nu /2)^{\nu /2}}{\Gamma (\nu /2)}}~{\frac {\exp\!\left({\frac {-\nu \tau ^{2}}{2 x}}\right)}{x^{1+\nu /2}}}\,.
\end{equation}
And we assume the mixture probability is drawn from a beta distribution (the conjugate prior of a Bernoulli distribution).
\begin{equation}
\operatorname{Beta}(x {\,\mid\,} \kappa, \omega) = {\frac {x^{\alpha -1}(1-x)^{\beta -1}}{\operatorname{B} (\alpha ,\beta )}}\,, \quad\,\,
  \alpha = \kappa \omega + 1\,, \quad\,\,
  \beta = \kappa(1 - \omega) + 1\,.
\end{equation}
Where $\operatorname{B}$ is the beta function.
Note that our beta distribution uses a slightly unconventional parameterization in terms of its concentration $\kappa$ and mode $\omega$, which will be explained later. Also note that in Equation~\ref{eq:probform} the degrees-of-freedom parameter $\nu$ of the scaled inverse chi-squared distribution is rescaled according to the mixture probability, as will also be discussed later.
This Bayesian mixture of Gaussians has four hyperparameters that will determine the behavior of GHT: $\nu \geq 0, \tau \geq 0, \kappa \geq 0, \omega \in [0, 1]$.

From this probabilistic framework we can derive our histogram thresholding algorithm.
Let us rewrite Equation~\ref{eq:bayesianprob} in terms of histogram locations $\myvec{x}$ and histogram counts $\myvec{n}$ with the summation written explicitly, and then take its logarithm to produce a log-likelihood:
\begin{equation}
    \sum_i \at{n}{i} \log\!\bigg( p\Big( \at{x}{i} \, \Big| \, z_i=0 \Big) p(z_i=0)\,+\, p\Big( \at{x}{i} \, \Big| \, z_i=1 \Big) p(z_i=1) \bigg)\,.
    \label{eq:bayesianloghist}
\end{equation}
Now we make a simplifying model assumption that enables a single-pass histogram thresholding algorithm.
We will assume that, at each potential ``split'' of the data, each of the two splits of the histogram is generated entirely by one of the two Gaussians in our mixture.
More formally, we consider all possible sorted assignments of each histogram bin to either of the two Gaussians and maximize the expected complete log-likelihood (ECLL) of that assignment.
Jensen's inequality ensures that this ECLL is a lower bound on the marginal likelihood of the data that we would actually like to maximize.
This is the same basic approach (though not the same justification) as in ``Minimum Error Thresholding'' (MET)~\cite{KittlerIllingworth} and other similar approaches, but where our technique differs is in how the likelihood and the parameters of these Gaussians are estimated.
For each split of the data, we assume the posterior distribution of the latent variables determining the ownership of each histogram bin by each Gaussian (the unknown $z$ values) are wholly assigned to one of the two Gaussians according to the split.
This gives us the following ECLL as a function of the assumed split location's array index $i$:
\newcommand{\mathcal{L}_i}{\mathcal{L}_i}
\begin{equation}
    \mathcal{L}_i = \sum_{j=0}^i \at{n}{j} \log\!\bigg( p\Big( \at{x}{j} \, \Big| \, z_j=0 \Big) \bigg) + \sum_{j=i+1}^{\operatorname{len}(\nvec)\mbox{-}1} \at{n}{j} \log\!\bigg(  p\Big( \at{x}{j} \, \Big| \, z_j=1 \Big) \bigg)\,.
    \label{eq:ecll}
\end{equation}
Our proposed algorithm is to simply iterate over all possible values of $i$ (of which there are $255$ in our experiments) and return the value that maximizes $\mathcal{L}_i$.
As such, our algorithm can be viewed as a kind of inverted expectation-maximization in which the sweep over threshold resembles an M-step and the assignment of latent variables according to that threshold resembles an E-step \cite{SalakhutdinovRG03}.
Our GHT algorithm will be defined as:
\begin{equation}
    \operatorname{GHT}(\myvec{x},  \myvec{n}; \nu, \tau, \kappa, \omega) = \at{x}{\argmax{\mathcal{L}_i}} \,. \label{eq:threshprob}
\end{equation}
This definition of GHT is a bit unwieldy, but using the preliminary math of 
Equation~\ref{eq:prelim} and ignoring global shifts and scales of $\mathcal{L}_i$ that do not affect the argmax it can be simplified substantially:
\begin{align}
  \myvec{\varscalar}^{(0)} &= \frac{\myvec{\mixturescalar}^{(0)} \nu \tau^2 + \myvec{d}^{(0)}}{\myvec{\mixturescalar}^{(0)} \nu + \myvec{\intweightscalar}^{(0)}} \quad\quad\quad
  \myvec{\varscalar}^{(1)} = \frac{\myvec{\mixturescalar}^{(1)} \nu \tau^2 + \myvec{d}^{(1)}}{\myvec{\mixturescalar}^{(1)} \nu + \myvec{\intweightscalar}^{(1)}} \nonumber \\
  \myvec{\myscorescalar}^{(0)} &= -\frac{\myvec{d}^{(0)}}{\myvec{\varscalar}^{(0)}} - \myvec{\intweightscalar}^{(0)} \log\!\left(\!\myvec{\varscalar}^{(0)}\!\right)  + 2\!\left( \myvec{\intweightscalar}^{(0)}\!+\!\kappa\omega \right) \log\!\left(\!\myvec{\intweightscalar}^{(0)}\!\right)  \nonumber \\
  \myvec{\myscorescalar}^{(1)} &= -\frac{\myvec{d}^{(1)}}{\myvec{\varscalar}^{(1)}}  - \myvec{\intweightscalar}^{(1)} \log\!\left(\!\myvec{\varscalar}^{(1)}\!\right)  + 2\!\left( \myvec{\intweightscalar}^{(1)}\!+\!\kappa(1\!-\!\omega) \right) \log\!\left(\!\myvec{\intweightscalar}^{(1)}\!\right) \nonumber  \\
  & \operatorname{GHT}(\myvec{x},  \myvec{n}; \nu, \tau, \kappa, \omega) = \at{x}{\argmax{\myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)}}} \label{eq:thresh}
\end{align}
Each $\myvec{\myscorescalar}^{(k)}$ can be thought of as a log-likelihood of the data for each Gaussian at each split, and each $\myvec{\varscalar}^{(k)}$ can be thought of as an estimate of the variance of each Gaussian at each split.
GHT simply scores each split of the histogram and returns the value of $\myvec{x}$ that maximizes the $\myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)}$ score (in the case of ties we return the mean of all $\myvec{x}$ values that maximize that score).
Because this only requires element-wise computation using the previously-described quantities in Equation~\ref{eq:prelim}, GHT requires just a single linear scan over the histogram.

The definition of each $\myvec{\varscalar}^{(k)}$ follows from our choice to model each Gaussian's variance with a scaled inverse chi-squared distribution: the $\nu$ and $\tau^2$ hyperparameters of the scaled inverse chi-squared distribution are updated according to the scaled sample variance (shown here as distortion $\myvec{d}^{(k)}$) to produce an updated posterior hyperparameter that we use as the Gaussian's estimated variance  $\myvec{\varscalar}^{(k)}$.
Using a conjugate prior update instead of the sample variance has little effect when the subset of $\myvec{n}$ being processed has significant mass (\textit{i}.\textit{e}.,  when $\myvec{\intweightscalar}^{(k)}$ is large) but has a notable effect when a Gaussian is being fit to a sparsely populated subset of the histogram.
Note that $\myvec{\varscalar}^{(k)}$ slightly deviates from the traditional conjugate prior update, which would omit the $\myvec{\mixturescalar}^{(k)}$ scaling on $\nu$ in the numerator and denominator (according to our decision to rescale each degrees-of-freedom parameter by each mixture probability).
This additional scaling counteracts the fact that, at different splits of our histogram, the total histogram count on either side of the split will vary. A ``correct'' update would assume a constant number of degrees-of-freedom regardless of where the split is being performed, which would result in the conjugate prior here having a substantially different effect at each split, thereby making the resulting $\myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)}$ scores not comparable to each other and making the argmax over these scores not meaningful.

The beta distribution used in our probabilistic formulation (parameterized by a concentration $\kappa$ and a mode $\omega$) has a similar effect as the ``anisotropy coefficients'' used by other models \cite{kapur1985new,pun1981entropic}: setting $\omega$ near $0$ biases the algorithm towards thresholding near the start of the histogram, and setting it near $1$ biases the algorithm towards thresholding near the end of the histogram. The concentration $\kappa$ parameter determines the strength of this effect, and setting $\kappa=0$ disables this regularizer entirely.
Note that our parameterization of concentration $\kappa$ differs slightly from the traditional formulation, where setting $\kappa=1$ has no effect and setting $\kappa < 1$ moves the predicted threshold \emph{away} from the mode. We instead assume the practitioner will only want to consider biasing the algorithm's threshold \emph{towards} the mode, and parameterize the distribution accordingly.
These hyperparameters allow for GHT to be biased towards outputs where a particular fraction of the image lies above or below the output threshold. For example, in an OCR / digit recognition context, it may be useful to inform the algorithm that the majority of the image is expected to not be ink, so as to prevent output thresholds that erroneously separate the non-ink background of the page into two halves.

Note that, unlike many histogram based thresholding algorithms, we do not require or assume that our histogram counts $\myvec{n}$ are normalized. In all experiments, we use raw counts as our histogram values. This is important to our approach, as it means that the behavior of the magnitude of $\myvec{n}$ and the behavior induced by varying our $\nu$ and $\kappa$ hyperparameters is consistent:
\newcommand{a}{a}
\begin{equation}
\forall a\in\mathbb{R}_{>0}\,\, \operatorname{GHT}(\myvec{x}, a \myvec{n}; a \nu, \tau, a \kappa, \omega) = \operatorname{GHT}(\myvec{n}, \myvec{x}; \nu, \tau, \kappa, \omega) \,.
\end{equation}
This means that, for example, doubling the number of pixels in an image is equivalent to halving the values of the two hyperparameters that control the ``strength'' of our two Bayesian model components.

Additionally, GHT's output and hyperparameters (excluding $\tau$, which exists solely to encode absolute scale) are invariant to positive affine transformations of the histogram bin centers:
\newcommand{b}{b}
\begin{equation}
\forall a\in\mathbb{R}_{>0}\,\, \forall b\in\mathbb{R}\,\, \operatorname{GHT}(a \myvec{x} + b, \myvec{n}; \nu, a \tau, \kappa, \omega) = \operatorname{GHT}(\myvec{n}, \myvec{x}; \nu, \tau, \kappa, \omega) \,.
\end{equation}
One could extend GHT to be sensitive to absolute intensity values by including conjugate priors over the means of the Gaussians, but we found little experimental value in that level of control. This is likely because the dataset we evaluate on does not have a standardized notion of brightness or exposure (as is often the case for binarization tasks).

In the following subsections we demonstrate how GHT generalizes three standard approaches to histogram thresholding by including them as special cases of our hyperparameter settings: Minimum Error Thresholding~\cite{KittlerIllingworth}, Otsu's method~\cite{Otsu}, and weighted percentile thresholding~\cite{percentile}.


\subsection{Special Case: Minimum Error Thresholding}
\label{subsec:met}

Because Minimum Error Thresholding (MET)~\cite{KittlerIllingworth}, like our approach, works by maximizing the expected complete log-likelihood of the input histogram under a mixture of two Gaussians, it is straightforward to express it using the already-defined quantities in Equation~\ref{eq:prelim}:
\begin{align}
    \myvec{\metscorescalar} = 1 + \myvec{\intweightscalar}^{(0)}\!\log\!\left(\!\frac{\myvec{d}^{(0)}}{\myvec{\intweightscalar}^{(0)}}\!\right) + \myvec{\intweightscalar}^{(1)}\!\log\!\left(\!\frac{\myvec{d}^{(1)}}{\myvec{\intweightscalar}^{(1)}}\!\right) \nonumber - 2\!\left(\!\myvec{\intweightscalar}^{(0)}\!\log\!\left(\!\myvec{\intweightscalar}^{(0)}\!\right) + \myvec{\intweightscalar}^{(1)}\!\log\!\left(\!\myvec{\intweightscalar}^{(1)}\!\right)\!\right)\,, \nonumber \\
    \operatorname{MET}(\myvec{x},  \myvec{n}) = \at{x}{\argmin{ \myvec{\metscorescalar}}}\,. \span 
\end{align}
Because GHT is simply a Bayesian treatment of this same process, it includes MET as a special case. If we set $\nu=0$ and $\kappa=0$ (under these conditions the values of $\tau$ and $\omega$ are irrelevant) in Equation~\ref{eq:thresh}, we see that $\myvec{\varscalar}^{(k)}$ reduces to $\myvec{d}^{(k)} / \myvec{\intweightscalar}^{(k)}$, which causes each score to simplify dramatically:
\begin{equation}
    \nu=\kappa=0 \implies \myvec{\myscorescalar}^{(k)} = -\myvec{\intweightscalar}^{(k)} -\myvec{\intweightscalar}^{(k)} \log\!\left(\!\frac{\myvec{d}^{(k)}}{\myvec{\intweightscalar}^{(k)}}\!\right) + 2 \myvec{\intweightscalar}^{(k)} \log\!\left(\!\myvec{\intweightscalar}^{(k)}\!\right)\,.
\end{equation}
Because $\myvec{\intweightscalar}^{(0)} + \myvec{\intweightscalar} = \sumvec{\nvec}$, the score maximized by GHT can be simplified:
\begin{equation}
    \nu=\kappa=0 \implies \myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)} = 1 - \sumvec{\nvec} - \myvec{\metscorescalar}\,.
\end{equation}
We see that, for these particular hyperparameter settings, the total score $\myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)}$ that GHT maximizes is simply an affine transformation (with a negative scale) of the score  $\myvec{\metscorescalar}$ that MET maximizes. This means that the index that optimizes either quantity is identical, and so the algorithms (under these hyperparameter settings) are equivalent:
\begin{equation}
\operatorname{MET}(\myvec{n}, \myvec{x}) = \operatorname{GHT}(\myvec{n}, \myvec{x}; 0, \cdot, 0, \cdot)\,.
\end{equation}


\subsection{Special Case: Otsu's Method}
\label{subsec:otsu}

Otsu's method for histogram thresholding~\cite{Otsu} works by directly maximizing inter-class variance for the two sides of the split histogram, which is equivalent to indirectly minimizing the total intra-class variance of those two sides~\cite{Otsu}.
Otsu's method can also be expressed using the quantities in Equation~\ref{eq:prelim}:
\begin{equation}
  \myvec{\otsuscorescalar} = \myvec{\intweightscalar}^{(0)} \myvec{\intweightscalar}^{(1)} \left(\myvec{\meanscalar}^{(0)} - \myvec{\meanscalar}^{(1)} \right)^2, \quad \operatorname{Otsu}(\myvec{n}, \myvec{x}) = \at{x}{\argmax{ \myvec{\otsuscorescalar} }}\,.
\end{equation}
To clarify the connection between Otsu's method and GHT, we rewrite the score $\myvec{\otsuscorescalar}$ of Otsu's method as an explicit sum of intra-class variances:
\begin{align}
\myvec{\otsuscorescalar} &= \sumvec{\nvec} \dotprod{\myvec{n}}{\myvec{x}^2} - \dotprod{\myvec{n}}{\myvec{x}}^2 - \sumvec{\nvec}\!\left(\myvec{d}^{(0)} + \myvec{d}^{(1)}\right)\,. \label{eq:distortion_equiv}
\end{align}
Now let us take the limit of the $\myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)}$ score that is maximized by GHT as $\nu$ approaches infinity:
\begin{equation}
\lim_{\nu \rightarrow \infty} \myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)} = -\frac{\myvec{d}^{(0)} + \myvec{d}^{(1)}}{\tau^2} + 2 \myvec{\intweightscalar}^{(0)} \log\!\left(\!\frac{\myvec{\intweightscalar}^{(0)}}{\tau}\!\right) + 2 \myvec{\intweightscalar}^{(1)} \log\!\left(\!\frac{\myvec{\intweightscalar}^{(1)}}{\tau}\!\right)\,.
\label{eq:otsulim}
\end{equation}
With this we can set the $\tau$ hyperparameter to be infinitesimally close to zero, in which case the score in Equation~\ref{eq:otsulim} becomes dominated by its first term.
Therefore, as $\nu$ approaches infinity and $\tau$ approaches zero, the score maximized by GHT is proportional to the (negative) score that is indirectly minimized by Otsu's method:
\begin{equation}
\nu \gg 0, \, \tau = \epsilon \implies \myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)} \approx -\frac{\myvec{d}^{(0)} + \myvec{d}^{(1)}}{\tau^2}\,,
\end{equation}
where $\epsilon$ is a small positive number.
This observation fits naturally with the well-understood equivalence of k-means (which works by minimizing distortion) and the asymptotic behavior of maximum likelihood under a mixture of Gaussians model as the variance of each Gaussian approaches zero~\cite{KulisJ12}.
Given this equivalence between the scores that are optimized by MET and GHT (subject to these specific hyperparameters) we can conclude that Otsu's method is a special case of GHT:
\begin{equation}
\operatorname{Otsu}(\myvec{n}, \myvec{x}) = \lim_{(\nu, \tau) \rightarrow (\infty, 0)} \operatorname{GHT}(\myvec{n}, \myvec{x}; \nu, \tau, 0, \cdot)\,.
\end{equation}



\begin{figure}[t]
    \centering
    \includegraphics[width=0.85\linewidth]{images/nu_sweep.png}
    \caption{
    On the right we have a toy input histogram (shown rotated), and on the left we have a plot showing the predicted threshold ($y$-axis) of GHT ($\tau=0.01$, $\kappa=0$, shown in black) relative to the hyperparameter value $\nu$ ($x$-axis).
    Note that the $y$-axis is shared across the two plots.
    Both MET and Otsu's method predict incorrect thresholds (shown in red and blue, respectively) despite the evident gap in the input histogram.
    GHT, which includes MET and Otsu's method as special cases ($\nu=0$ and $\nu=\infty$, respectively) allows us to interpolate between these two extremes by varying $\nu$, and produces the correct threshold for a wide range of values $\nu \in [\sim\!2^2, \sim\!2^8]$. 
    }
    \label{fig:nu_sweep}
\end{figure}

We have demonstrated that GHT is a generalization of Otsu's method (when $\nu = \infty$) and MET (when $\nu = 0$), but for this observation to be of any practical importance, there must be value in GHT when $\nu$ is set to some value in between those two extremes.
To demonstrate this, in Figure~\ref{fig:nu_sweep} we visualize the effect of $\nu$ on GHT's behavior. Similar to the example shown in Figure~\ref{fig:teaser}, we see that both Otsu's method and MET both perform poorly when faced with a simple histogram that contains two separate and uneven uniform distributions: MET prefers to express all histogram elements using a single Gaussian by selecting a threshold at one end of the histogram, and Otsu's method selects a threshold that splits the larger of the two modes in half instead of splitting the two modes apart.
But by varying $\nu$ from $0$ to $\infty$ (while also setting $\tau$ to a small value) we see that a wide range of values of $\nu$ results in GHT correctly separating the two modes of this histogram. Additionally, we see that the range of hyperparameters that reproduces this behavior is wide, demonstrating the robustness of performance with respect to the specific value of this parameter.


\subsection{Relationship to Histogram Bin Width}

Most histogram thresholding techniques (including GHT as it is used in this paper) construct a histogram with as many bins as there are pixel intensities --- an image of 8-bit integers results in a histogram with 256 bins, each with a bin width of 1. However, the performance of histogram thresholding techniques often depends on the selection of the input histogram's bin width, with a coarse binning resulting in more ``stable'' performance and a fine binning resulting in a more precisely localized threshold. Many histogram thresholding techniques therefore vary the bin width of their histograms, either by treating it as a user-defined hyperparameter~\cite{coudray2010robust} or by determining it automatically~\cite{Kadhim,onumanyi2017modified} using classical statistical techniques~\cite{doane1976aesthetic,sturges1926choice}. Similarly, one could instead construct a fine histogram that is then convolved with a Gaussian before thresholding~\cite{lindblad2000histogram}, which is equivalent to constructing the histogram using a Parzen window. Blurring or coarsening histograms both serve the same purpose of filtering out high frequency variation in the histogram, as coarsening a histogram is equivalent blurring (with a box filter) and then decimating a fine histogram.

This practice of varying histogram bins or blurring a histogram can be viewed through the lens of GHT.
Consider a histogram $(\myvec{n}, \myvec{x})$, and let us assume (contrary to the equivalence described in Equation~\ref{eq:sorted}) that the spacing between the bins centers in $\myvec{x}$ is constant.
Consider a discrete Gaussian blur filter $\boldsymbol{f}(\blursigma)$ with a standard deviation of $\sigma$, where that filter is normalized ($\norm{\boldsymbol{f}(\blursigma)}_1=1$).
Let us consider $\myvec{n} * \boldsymbol{f}(\blursigma)$, which is the convolution of histogram counts with our Gaussian blur (assuming a ``full'' convolution of $\myvec{n}$ and $\myvec{x}$ with zero boundary conditions).
This significantly affects the sample variance of the histogram $v$, which (because convolution is linear) is:
\begin{equation}
v = \frac{\sum_i \at{(\myvec{n} * \boldsymbol{f}(\blursigma))}{i} ( \at{x}{i} - \mu)^2}{\sumvec{\nvec}} = \frac{\sumvec{\nvec} \sigma^2 + d}{\sumvec{\nvec}}\,, 
\quad\quad d = \sum_i \at{n}{i} ( \at{x}{i} - \mu)^2 \,. 
\end{equation}
We use $v$ and $d$ to describe sample variance and distortion as before, though here these values are scalars as we compute a single estimate of both for the entire histogram.
This sample variance $v$ resembles the definition of $\myvec{\varscalar}^{(k)}$ in Equation~\ref{eq:thresh}, and we can make the two identical in the limit by setting GHT's hyperparameters to $\nu = \epsilon \sumvec{\nvec}$ and $\tau = \sigma / \sqrt{\epsilon}$, where $\epsilon$ is a small positive number.
With this equivalence we see that GHT's $\tau$ hyperparameter serves a similar purpose as coarsening/blurring the input histogram --- blurring the input histogram with a Gaussian filter with standard deviation of $\sigma$ is roughly equivalent to setting GHT's $\tau$ parameter proportionally to $\sigma$. Or more formally:
\begin{equation}
\operatorname{MET}\!\left( \myvec{n}, \myvec{x} * \boldsymbol{f}(\blursigma) \right) \approx \operatorname{GHT}\!\left(\myvec{n}, \myvec{x}; \sumvec{\nvec} \epsilon, \sigma / \sqrt{\epsilon}, 0, \cdot \right)\,.
\end{equation}
Unlike the other algorithmic equivalences stated in this work, this equivalence is approximate: changing $\sigma$ is not exactly equivalent to blurring the input histogram. Instead, it is equivalent to, at each split of $\myvec{n}$, blurring the left and right halves of $\myvec{n}$ independently --- the histogram cannot be blurred ``across'' the split. Still, we observed that these two approaches produce identical threshold results in the vast majority of instances. 

This relationship between GHT's $\tau$ hyperparameter and the common practice of varying a histogram's bin width provides a theoretical grounding of both GHT and other work: GHT may perform well because varying $\tau$ implicitly blurs the input histogram, or prior techniques based on variable bin widths may work well because they are implicitly imposing a conjugate prior on sample variance. This relationship may also shed some light on other prior work exploring the limitations of MET thresholding due to sample variance estimates not being properly ``blurred'' across the splitting threshold~\cite{cho1989improvement}.


\subsection{Special Case: Weighted Percentile}
\label{sec:wprctile}

A simple approach for binarizing some continuous signal is to use the nth percentile of the input data as the threshold, such as its median.
This simple baseline is also expressible as a special case of GHT.
If we set $\kappa$ to a large value we see that the score $\myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)}$ being maximized by GHT becomes dominated by the terms of the score related to $\log(\myvec{\mixturescalar}^{(k)})$:
\begin{equation}
  \resizebox{\textwidth}{!}{$
  \displaystyle
  \kappa \gg 0\!\implies\!\frac{\myvec{\myscorescalar}^{(0)} + \myvec{\myscorescalar}^{(1)}}{2} \approx \left( \sumvec{\nvec} + \kappa \right)\log\left(\sumvec{\nvec}\right) + \kappa \left( \omega \log\!\left(\!\myvec{\mixturescalar}^{(0)}\!\right) + (1 - \omega) \log\!\left(\!1 - \myvec{\mixturescalar}^{(0)}\!\right) \right)\,.
  $}
\end{equation}
By setting its derivative to zero we see that this score is maximized at the split location $i$ where $\at{\pi^{(0)}}{i}$ is as close as possible to $\omega$. This condition is also satisfied by the $(100\omega)$th percentile of the histogrammed data in $\mathrm{I}$, or equivalently, by taking the $(100\omega)$th weighted percentile of $\myvec{x}$ (where each bin is weighted by $\myvec{n}$).
From this we can conclude that the weighted percentile of the histogram (or equivalently, the percentile of the image) is a special case of GHT:
\begin{equation}
\operatorname{wprctile}(\myvec{x},  \myvec{n}, 100\omega) = \lim_{\kappa \rightarrow \infty} \operatorname{GHT}(\myvec{n}, \myvec{x}, 0, \cdot, \kappa, \omega)\,.
\end{equation}
This also follows straightforwardly from how GHT uses a beta distribution: if this beta distribution's concentration $\kappa$ is set to a large value, the mode of the posterior distribution will approach the mode of the prior $\omega$.

In Figure~\ref{fig:omega_sweep} we demonstrate the value of this beta distribution regularization.
We see that Otsu's method and MET behave unpredictably when given multi-modal inputs, as the metrics they optimize have no reason to prefer a split that groups more of the histogram modes above the threshold versus one that groups more of the modes below the threshold.
This can be circumvented by using the weighted percentile of the histogram as a threshold, but this requires that the target percentile $\omega$ be precisely specified beforehand: For any particular input histogram, it is possible to identify a value for $\omega$ that produces the desired threshold, but this percentile target will not generalize to another histogram with the same overall appearance but with a slightly different distribution of mass across that desired threshold level.
GHT, in contrast, is able to produce sensible thresholds for all valid values of $\omega$ --- when the value is less than $\sfrac{1}{2}$ the recovered threshold precisely separates the first mode from the latter two, and when the value is greater than $\sfrac{1}{2}$ it precisely separates the first two modes from the latter one. As such, we see that this model component provides an effective means to bias GHT towards specific kinds of thresholds, while still causing it to respect the relative spread of histogram bin counts.



\begin{figure}[t]
    \centering
    \includegraphics[trim=80 10 70 40, clip, width=0.85\linewidth]{images/omega_sweep.png}
    \caption{
    On the right we have a toy input histogram (shown rotated), and on the left we have a plot showing the predicted threshold ($y$-axis) of GHT ($\nu = 200$, $\tau=0.01$, $\kappa = 0.1$, shown in black) relative to the hyperparameter value $\omega$ ($x$-axis).
    Note that the $y$-axis is shared across the two plots.
    Both MET and Otsu's method (shown in red and blue, respectively) predict thresholds that separate two of the three modes from each other arbitrarily. Computing the weighted percentile of the histogram (shown in yellow) can produce \emph{any} threshold, but reproducing the desired split requires exactly specifying the correct value of $\omega$, which likely differs across inputs. GHT (which includes these three other algorithms as special cases)
    produces thresholds that accurately separate the three modes, and the location of that threshold is robust to the precise value of $\omega$ and depends only on it being below or above $\sfrac{1}{2}$.
    }
    \label{fig:omega_sweep}
\end{figure}

\section{Experiments}


\begin{table*}[h!]
\centering
\resizebox{\textwidth}{!}{  
\begin{tabular}{@{}l|cccc||ccc@{}}
\quad Algorithm & $\nu$ & $\tau$ & $\kappa$ & $\omega$ & $F_1 \times 100 \uparrow$ & PSNR $\uparrow$ & DRD~\cite{lu2004distance} $\downarrow$ \\ \hline\hline
{\bf GHT (MET Case)} & -  & -  & -  & -                                                    &                    $60.40 \pm 20.65$ &                    $11.21 \pm 3.50$ &                    $45.32 \pm 41.35$ \\
Kefali \textit{et al}.  \cite{dibco2016,Sari2014TextEF} &&&&                                          &                    $76.10 \pm 13.81$ &                    $15.35 \pm 3.19$ &                    $9.16 \pm 4.87$ \\
Raza \cite{dibco2016} &&&&                                                                 &                    $76.28 \pm 9.71$ &                    $14.21 \pm 2.21$ &                    $15.14 \pm 9.42$ \\
{\bf GHT (wprctile Case)} & -  & -  & $10^{60}$  & $2^{-3.75}$                             &                    $76.77 \pm 14.50$ &                    $15.44 \pm 3.40$ &                    $12.91 \pm 17.19$ \\
Sauvola \cite{dibco2016,Sauvola2000AdaptiveDI} &&&&                                        &                    $82.52 \pm 9.65$ &                    $16.42 \pm 2.87$ &                    $7.49 \pm 3.97$ \\
Khan \& Mollah \cite{dibco2016} &&&&                                                       &                    $84.32 \pm 6.81$ &                    $16.59 \pm 2.99$ &                    $6.94 \pm 3.33$ \\
Tensmeyer \& Martinez \cite{FCN,tensmeyer2017,Wolf2002TextLE} &&&&                         &                    $85.57 \pm 6.75$ &                    $17.50 \pm 3.43$ &                    $5.00 \pm 2.60$ \\
de Almeida \& de Mello \cite{dibco2016} &&&&                                               &                    $86.24 \pm 5.79$ &                    $17.52 \pm 3.42$ &                    $5.25 \pm 2.88$ \\
Otsu's Method \cite{Otsu,dibco2016} &&&&                                                            &                    $86.61 \pm 7.26$ &                    $17.80 \pm 4.51$ &                    $5.56 \pm 4.44$ \\
{\bf GHT (No wprctile)} & $2^{50.5}$  & $2^{0.125}$  & -  & -                              &                    $87.16 \pm 6.32$ &                    $17.97 \pm 4.00$ &                    $5.04 \pm 3.17$ \\
{\bf GHT (Otsu Case)} & $10^{60}$  & $10^{-15}$  & -  & -                                  &                    $87.19 \pm 6.28$ &                    $17.97 \pm 4.01$ &                    $5.04 \pm 3.16$ \\
Otsu's Method (Our Impl.) \cite{Otsu} &&&&                                                 &                    $87.19 \pm 6.28$ &                    $17.97 \pm 4.01$ &                    $5.04 \pm 3.16$ \\
Nafchi \textit{et al}.  - 1 \cite{dibco2016,Nafchi} &&&&                                              &                    $87.60 \pm 4.85$ &                    $17.86 \pm 3.51$ &                    $4.51 \pm 1.62$ \\
Kligler \cite{Howe2013,Katz2007,Kligler2017} &&&&                                          &                    $87.61 \pm 6.99$ &                    $18.11 \pm 4.27$ &                    $5.21 \pm 5.28$ \\
Roe \& de Mello \cite{dibco2016} &&&&                                                      &                    $87.97 \pm 5.17$ &                    $18.00 \pm 3.68$ &                    $4.49 \pm 2.65$ \\
Nafchi \textit{et al}.  - 2 \cite{dibco2016,Nafchi} &&&&                                              &                    $88.11 \pm 4.63$ &                    $18.00 \pm 3.41$ &                    $4.38 \pm 1.65$ \\
Hassa{\"i}ne \textit{et al}. - 1 \cite{Hassane2011,dibco2016} &&&&                                    &                    $88.22 \pm 4.80$ &                    $18.22 \pm 3.41$ &                    $4.01 \pm 1.49$ \\
Hassa{\"i}ne \textit{et al}.  - 2 \cite{Hassane2012,dibco2016} &&&&                                   & \cellcolor{yellow} $88.47 \pm 4.45$ & \cellcolor{yellow} $18.29 \pm 3.35$ & \cellcolor{orange} $3.93 \pm 1.37$ \\
Hassa{\"i}ne \textit{et al}.  - 3 \cite{Hassane2012,Hassane2011,dibco2016} &&&&                       & \cellcolor{orange} $88.72 \pm 4.68$ & \cellcolor{orange} $18.45 \pm 3.41$ & \cellcolor{red}    $3.86 \pm 1.57$ \\
{\bf GHT} & $2^{29.5}$  & $2^{3.125}$  & $2^{22.25}$  & $2^{-3.25}$                        & \cellcolor{red}    $88.77 \pm 4.99$ & \cellcolor{red}    $18.55 \pm 3.46$ & \cellcolor{yellow} $3.99 \pm 1.77$ \\
\hline
\rowcolor{lightgray} Oracle Global Threshold &&&&                                          & $90.69 \pm 3.92$ & $19.17 \pm 3.29$ & $3.57 \pm 1.84$
\end{tabular}
}
\caption{
Results on the 2016 Handwritten Document Image Binarization Contest (H-DIBCO) challenge \cite{dibco2016}, in terms of the arithmetic mean and standard deviation of F1 score (multiplied by $100$ to match the conventions used by \cite{dibco2016}), PSNR, and Distance-Reciprocal Distortion (DRD)~\cite{lu2004distance}. The scores for all baseline algorithms are taken from \cite{dibco2016}.
Our GHT algorithm and it's ablations and special cases (some of which correspond to other algorithms) are indicated in bold. The settings of the four hyperparameters governing GHT's behavior are also provided.
}
\label{table:results}
\end{table*}

To evaluate GHT, we use the 2016 Handwritten Document Image Binarization Contest (H-DIBCO) challenge \cite{dibco2016}. This challenge consists of images of handwritten documents alongside ground-truth segmentations of those documents, and algorithms are evaluated by how well their binarizations match the ground truth.
We use the 2016 challenge because this is the most recent instantiation of this challenge where a global thresholding algorithm (one that selects a single threshold for use on all pixels) can be competitive:
Later versions of this challenge use input images that contain content outside of the page being binarized which, due to the background of the images often being black, renders any algorithm that produces a single global threshold for binarization ineffective.
GHT's four hyperparameters were tuned using coordinate descent to maximize the arithmetic mean of $F_1$ scores (a metric used by the challenge) over a small training dataset (or perhaps more accurately, a ``tuning'' dataset).
For training data we use the $8$ hand-written images from the 2013 version of the same challenge, which is consistent with the tuning procedure advocated by the challenge (any training data from past challenges may be used, though we found it sufficient to use only relevant data from the most recent challenge).
Input images are processed by taking the per-pixel $\operatorname{max}()$ across color channels to produce a grayscale image, computing a single 256-bin histogram of the resulting values, applying GHT to that histogram to produce a single global threshold, and binarizing the grayscale input image according to that threshold.

In Table~\ref{table:results} we report GHT's performance against all algorithms evaluated in the challenge~\cite{dibco2016}, using the $F_1$, PSNR, and Distance-Reciprocal Distortion~\cite{lu2004distance} metrics used by the challenge~\cite{dibco2016}.
We see that GHT produces the lowest-error of all entrants to the H-DIBCO 2016 challenge for two of the three metrics used. GHT outperforms Otsu's method and MET by a significant margin, and also outperforms or matches the performance of significantly more complicated techniques that rely on large neural networks with thousands or millions of learned parameters, and that produce an arbitrary per-pixel binary mask instead of GHT's single global threshold. This is despite GHT's simplicity: it requires roughly a dozen lines of code to implement (far fewer if an implementation of MET or Otsu's method is available), requires no training, and has only four parameters that were tuned on a small dataset of only 8 training images. 
See Figures~\ref{fig:results1}-\ref{fig:results4} for visualizations of GHT alongside MET and Otsu's method.

We augment Table~\ref{table:results} with some additional results not present in \cite{dibco2016}. We present an ``Oracle Global Threshold'' algorithm, which shows the performance of an oracle that selects the best-performing (according to $F_1$) global threshold individually for each test image.
We present the special cases of GHT that correspond to MET and Otsu's method, to verify that the special case corresponding to Otsu's method performs identically to our own implementation of Otsu's method and nearly identically to the implementation presented in \cite{dibco2016}.
We present additional ablations of GHT to demonstrate the contribution of each algorithm component: The ``wprctile Case'' model sets $\kappa$ to an extremely large value and tunes $\omega$ on our training set, and performs poorly.
The ``No wprctile'' model sets $\kappa=0$ and exhibits worse performance than complete GHT.

See Appendix A for reference implementations of GHT, as well as reference implementations of the other algorithms that were demonstrated to be special cases of GHT in Sections~\ref{subsec:met}-\ref{sec:wprctile}.

\section{Conclusion}

We have presented Generalized Histogram Thresholding, a simple, fast, and effective technique for histogram-based image thresholding. GHT includes several classic techniques as special cases (Otsu's method, Minimum Error Thresholding, and weighted percentile thresholding) and thereby serves as a unifying framework for those discrete algorithms, in addition to providing a theoretical grounding for the common practice of varying a histogram's bin width when thresholding.
GHT is exceedingly simple: it can be implemented in just a dozen lines of python (far fewer if an implementation of MET or Otsu's method is available) and has just four tunable hyperparameters.
Because it requires just a single sweep over a histogram of image intensities, GHT is fast to evaluate: its computational complexity is comparable to Otsu's method.
Despite its simplicity and speed (and its inherent limitations as a \emph{global} thresholding algorithm) GHT outperforms or matches the performance of all submitted techniques on the 2016 H-DIBCO image binarization challenge --- including deep neural networks that have been trained to produce arbitrary per-pixel binarizations.

\bibliographystyle{splncs04}
