\chapter{Midway image equalization}
\label{chapter_midway}

\begin{resume}
On désigne par égalisation mi-chemin toute méthode donnant à deux images un histogramme commun, tout en respectant au maximum leurs radiométries de départ. Ce chapitre présente une analyse axiomatique de l'égalisation mi-chemin, qui conduit à deux méthodes distinctes. Celles-ci sont comparées des points de vue théorique et pratique, en analysant leurs effets respectifs sur le bruit de quantification et le spectre des images. Les deux méthodes sont ensuite testées sur plusieurs exemples réels.
\end{resume}

\medskip

\noindent \begin{abstract}
  {\rm Midway image equalization} means any method giving
   to a pair of images the same histogram, while
maintaining  as much as possible their previous grey level
dynamics. In this chapter, we present an axiomatic analysis of image
equalization which leads us to derive two possible methods. Both 
methods are then compared in theory and in practice for two
reliability criteria, namely their effect on quantization noise
and on the support of the Fourier spectrum. A mathematical
analysis of the properties of these methods is
performed. Their algorithms are described and they are tested on such
typical pairs as satellite image stereo pairs and different
photographs of a same painting.
\end{abstract}

\medskip

This chapter corresponds to the article ``Midway image equalization''~\cite{delon04midway}.

\medskip

\section{Introduction}

The comparison of two images, in order to extract a mutual
information, is one of the main themes in computer vision. The
pair of  images can be obtained  in many ways: they can be a
stereo pair, two images of the same object (a painting for
example), multi-channel images of the same region, images of a
movie, etc. This comparison  is perceptually greatly improved if
both images have the same grey level dynamics. In addition, a lot
of image comparison algorithms, based on grey level, take as basic
assumption that intensities of corresponding points in both images
are equal. It has, however, been often emphasized that the
constant brightness hypothesis is nothing but a visual illusion 
\cite{Metz}. As it is well known by experts in stereo vision, this
assumption is generally false for stereo pairs and deviations from
this assumption cannot even be modelled by affine transforms
\cite{Cox}. Consequently, if we want to compare visually and
numerically two images, it is useful to give them first the same
dynamic range and luminance. For now, we will limit the study to grey level
images (the case of color images will be studied in chapter~\ref{chapter_transport}).

We are led to the following problem.

\begin{itemize} 

\item \indent From two images $u_0$ and $u_1$,
construct by contrast changes two images $\tilde{u}_0$ and
$\tilde{u}_1$, which have the same histogram.

\item This common histogram $h$ should stand midway between the previous histograms of
$u_0$ and $u_1$, and be as close as possible to each of them. This
treatment must avoid to favor one histogram rather than the other.
\end{itemize} With this intention, two additional aspects are
taken into consideration:
\begin{itemize}

\item Digital images being quantized,
the quantization noise must remain as uniform as possible.

\item  $u_0$ and $u_1$ must remain ``well sampled'' if they were,
namely, the support of their Fourier spectrum should not expand too
much by the applied contrast change.
\end{itemize}
Obviously, the involved contrast changes  have to be increasing
and continuous, for  reasons of reversibility and coherence. We will call any
method addressing the problem with the above mentioned  requirements a 
\textit{midway image equalization}.

As we shall see, the requirements yield naturally two methods. The
first one is to choose the midway histogram as an "inverse mean''
of the previous ones. We will see how this method catches up with
the work of Cox et al. \cite{Cox}, who addressed the same problem
and proposed a discrete algorithm. The second one we shall derive
is inspired by \textit{histogram specification}. It consists, if
$H_0$ and $H_1$ are the cumulative histograms of $u_0$ and $u_1$,
to find a contrast change $\psi$ such that $\psi \circ \psi
={H_0}^{-1} \circ H_1$ and to replace $u_1$ by $\psi(u_1)$ and
$u_0$ by ${\psi}^{-1}(u_0)$. This second method is more difficult
to approach, because there exist several solutions to the
``\textit{square root}'' problem, and their graphs can be, as we
shall see, oscillatory. We will introduce and describe the two
announced methods in the next section. The third part of the chapter
deals with the theoretical construction and mathematical
properties of the \textit{ square iterative root}. We will see how
to choose this square root optimally towards the quantization
noise. The final section will be devoted to
experiments, to the comparison of both methods with specification, 
and their effects on the images spectra and on the quantization noise.




\section{How to bring two histograms closer?}

\subsection{Some definitions and notations}

Let us begin with some notations that we will use all along the first chapters.
Let $\Omega \subseteq {\mathbb{R}}^2$ be the image domain, and $u:\Omega \rightarrow [0,1]$ \footnote{This definition is not
restrictive, the range of an image can be put on $[0,1]$ by an
affine transformation.}  an
image. If $\mu_u$ is the measure representing the distribution of its values on $[0,1]$, we note $H_u:[0,1]\rightarrow [0,1]$ the corresponding repartition function. In image analysis terms, $H_u$ is called the {\em cumulative histogram} of $u$. The histogram $h_u$ of an image $u$ is the distributional derivative of $H_u$.
%If the measure $\mu_u$ is absolutely continuous with respect to the Lebesgue measure,  its density $h_u$ is called the histogram of $u$. 
%The histogram and cumulative histogram of an image are linked by the relation $H_u(\lambda)=\int_0^{\lambda} h_u(t)dt$.

In practice, we have access to discrete images, defined on a finite set $\Omega$ of size $N$. Thus, the cumulative histogram of $u$ can be computed as $H_u(\lambda)= |\{x \in \Omega;\; u(x) \leq \lambda\}|/N$. An image, its histogram and its cumulative histogram can be observed on the first line of Figure~\ref{marais}.

\begin{figure}[h]
\centerline{
\subfigure[{\em Image $u$, the corresponding grey level histogram $h_u$, and the cumulative histogram $H_u$.}]{
\includegraphics[width=6cm]{marais_carre}
\includegraphics[width=4.5cm, height=4.5cm]{marais_histo}
\includegraphics[width=4.5cm, height=4.5cm]{marais_histo_cumul}}}
\centerline{
\subfigure[{\em  Equalized image $H_u(u)$, its histogram and its cumulative histogram. In the discrete case, histogram equalization flattens the histogram as much as possible. After equalization, if this histogram is smoothed with a small kernel, it quickly becomes constant.}]{
\includegraphics[width=6cm]{marais_eq_carre}
\includegraphics[width=4.5cm, height=4.5cm]{marais_eq_histo}
\includegraphics[width=4.5cm, height=4.5cm]{marais_eq_histo_cumul}}}
\caption{\em  We see on this example that image equalization can be visually very destructive. In this marshland image, after equalization, the water is almost not distinguishable from the vegetation, especially in the left part of the image (see in particular the part in the rectangular zone). }
\label{marais}
\end{figure}



\begin{proposition}
Let $\ph:[0,1]\rightarrow [0,1]$ be a continuous and strictly
increasing function (such a function will be called an
\textit{increasing homeomorphism} on $[0,1]$). If  $u$ is an image, then $\ph(u)$ is another image, whose cumulative histogram is $H_u \circ {\ph}^{-1}$.
\end{proposition}
We point out that such a transformation, called {\bf contrast change}, preserves the order of the grey levels in $v$ and does not create or cancel any image
information. However, if $\ph$ is  just nondecreasing, information can be lost. If the cumulative distribution $H_u$ is continuous and strictly increasing, $H_u(u)$ is an image whose distribution
function is uniform on $[0,1]$.  This transformation is called
\textit{histogram equalization}. General information on histogram
modification can be found, for example, in \cite{woods}.

\subsection{A description of classical methods.}

We consider two images $u_0$ and $u_1$, and their cumulative histograms
$H_0$ and $H_1$. In order to give them the same radiometry, 
we have to choose an intermediate cumulative histogram $H$, and to make $u_0$ and $u_1$ correspond to it by  contrast changes. In the light of
the above proposition, we see that the only way to perform such a
transformation is to replace each image $u_i$ by $\ph_i(u_i)$
where $\ph_i={H}^{-1} \circ H_i$ ($i \in \{0,1\}$). All the
methods described below are  based on this principle.


\medskip
\begin{description}

\item[\bf 1. ``Classical equalization'' method]

  In this case, $H$ is the identity function on $[0,1]$. Each image $u_i$ becomes
  $H_i(u_i)$, which means that each histogram is equalized. Now, if 
  both histograms are similar and far from
  being constant, it is not judicious to flatten them to make them
  correspond. Indeed, any contrast change affects 
the aspect of an image, especially when this contrast change is far from the identity. It results in two main artifacts. First, the images being discrete in practice, contrast changes can produce a loss of information. This loss is particularly strong when the contrast change is far from the identity. Secondly, equalization can also be perceptually bad because it enhances the quantization noise in some places, especially in low contrasted regions, and produces this way artificial edges or textures (see Figure~\ref{piano}). 

  \medskip

\begin{figure}
\centerline{
\includegraphics[width=7cm]{piano}
\includegraphics[width=7cm]{piano_eq}
}
\caption{\em Effect of histogram equalization on the quantization noise. On the left, the original image. On the right, the same image after histogram equalization. The effect of this equalization on the dark areas (the piano, the left part of the wall), which are low contrasted, is perceptually dramatic.}
\label{piano}
\end{figure}
 
\item[\bf 2. Histogram Specification]
  
  One of the images is chosen as a reference, and  $H$ is taken as
  the cumulative histogram of this reference. For example, if $u_0$
  is the reference, $u_0$ remains unchanged and $u_1$ becomes
  ${H_0}^{-1} \circ H_1 (u_1)$, which has the same histogram as
  $u_0$ (we say that $h_1$ is \textit{specified} on $h_0$). This
  method is called \textit{histogram specification} \cite{woods}. It
  can be judicious if both images have the same kind of dynamic
  range. For the same reason as in
  equalization,  this method can also product contouring artifacts.
  We can expect a little alteration of each image to be less affecting than a
  large modification of one of them.
 % \psfrag{n1}{$n_0$}
%  \psfrag{n2}{$n_1$} \psfrag{n12}{$\frac{n_0+n_1}{2}$}
%  \psfrag{h}{$h$} \psfrag{h1}{$h_1$} \psfrag{h2}{$h_2$}
%  \psfrag{h3}{$h_3$} \psfrag{h4}{$h_4$} \psfrag{H1}{$H_0$}
%  \psfrag{H2}{$H_1$} \psfrag{H3}{$H_3$} \psfrag{H}{$H$}

\begin{figure}
\begin{center}
\begin{tabular}{>{\centering}m{4cm} >{\centering}m{4cm}  >{\centering}m{4cm}}
\includegraphics[width=3.5cm]{degrade2} &
\includegraphics[width=3cm, angle=270]{degrade2_histo} &
\includegraphics[width=3cm, angle=270]{degrade2_histocum}
\end{tabular}
\vspace{0.2cm}

\begin{tabular}{>{\centering}m{4cm} >{\centering}m{4cm}  >{\centering}m{4cm}}
\includegraphics[width=3.5cm]{degrade1} &
\includegraphics[width=3cm, angle=270]{degrade1_histo} &
\includegraphics[width=3cm, angle=270]{degrade1_histocum}
\end{tabular}
\vspace{0.2cm}

\begin{tabular}{>{\centering}m{4cm} >{\centering}m{4cm}  >{\centering}m{4cm}}
\includegraphics[width=3.5cm]{degrade_moyenne} &
\includegraphics[width=3cm, angle=270]{moyenne_histo} &
\includegraphics[width=3cm, angle=270]{moyenne_histocum}
\end{tabular}
\vspace{0.2cm}

\begin{tabular}{>{\centering}m{4cm} >{\centering}m{4cm}  >{\centering}m{4cm}}
\includegraphics[width=3.5cm]{degrade_midway} &
\includegraphics[width=3cm, angle=270]{midway_histo} &
\includegraphics[width=3cm, angle=270]{midway_histocum}
\end{tabular}
\end{center}
%\begin{figure}
%\centerline{\subfigure{
%\includegraphics[width=3cm]{degrade1}
%\hspace{0.2cm}
%\includegraphics[width=9cm]{moyenne3}
%\hspace{0.2cm}
%\includegraphics[width=3cm]{degrade2}}}
%\centerline{
%\includegraphics[width=3cm]{degrade_moyenne}
%\hspace{0.2cm}
%\includegraphics[width=9cm]{moyenne4}
%\hspace{0.2cm}
%\includegraphics[width=3cm]{degrade_midway}}
\caption{\em First and second lines: two gradations images, their grey level histograms $h_0$ and $h_1$, and their cumulative histograms $H_0$ and $H_1$. The grey levels of the first gradation are around $n_0$, and the grey levels of the second one are around $n_1$. 
Third line:  intermediate image obtained when using the diret mean $\tilde{h}$ between $h_0$ and $h_1$ as a target histogram. Note the artificial edge that appears in the image when we use the mean histogram. 
Fourth line: intermediate image obtained when using as a target histogram the ``natural'' intermediate between $h_0$ and $h_1$, called $h$.}
\label{moyenne.fig}
\end{figure}

\medskip

\item[\bf 3. Mean]
  
  If we want $H$ as a compromise between $H_0$ and $H_1$, the
  most elementary function that we can imagine is the average of those
  two cumulative histograms. However, the following example proves
  that this idea is not always judicious: let $u_0$ and $u_1$ be two
  images, whose histograms $h_0$ and $h_1$ are translated versions of each other.
Such images are represented in the two first lines of Figure~\ref{moyenne.fig}. 
  In this example, the images have no common grey
  level.  Now,  their histograms $h_0$ and $h_1$ are
  two similar crenel functions, respectively centered at $n_0$
  and $n_1$. It follows that the most natural intermediary (in the
  sense stated in the introduction) for them should be the same type of crenel function, centered exactly at the middle of $n_0$ and $n_1$ (we note $h$
  this ``natural'' intermediary). 
  Now, let $\tilde{h}$ be the direct average of $h_0$ and $h_1$. 
We see that $\tilde{h}$ is made up of two
  successive modes instead of one. Consequently, $\tilde{h}$ is not a good
  compromise between $h_0$ and $h_1$. Indeed, the composition of
  $u_0$ by ${\tilde{H}}^{-1} \circ H_0$  spreads out the grey levels of
  the image on two regions centered at $n_0$ and $n_1$, in order
  to make its histogram correspond as much as possible to $\tilde{h}$.
  This ``stretching'' has at least one drawback: the grey levels of
  $u_0$, which were uniformly distributed around $n_0$, are
  se\-pa\-rated in two distant regions, and this separation creates
  some structure in $u_0$ which did not appear before (third line of Figure~\ref{moyenne.fig}). If we use $H$
  instead of $\tilde{H}$, this phenomenon  does not appear: the grey
  levels of $u_0$ and $u_1$ are just translated in order to be
  distributed around $\frac{n_0+n_1}{2}$ (fourth line of Figure~\ref{moyenne.fig}). Here, this ``natural'' intermediary is clearly related to $H_0$ and $H_1$ by the relation :
$$H^{-1}=\frac{H_0^{-1}+H_1^{-1}}{2}.$$
  
\end{description}



\subsection{Axioms of midway equalization}

We keep the notations of the previous section. We search for two
contrast changes $\ph_0$ and $\ph_1$ such that $\ph_0(u_0)$ and
$\ph_1(u_1)$ have the same histogram. Thus, $${\ph_0}^{-1} \circ
\ph_1={H_0}^{-1} \circ H_1.$$ In order to select the right method,
two axiomatic approaches are possible. A first one would enhance
invariance and conservation principles for all images. Now,
another more intuitive way to specify axiomatically a method is to
state what the expected result should be for a short list of key
examples. Here, we choose this second axiomatic specification
method by considering a very simple particular case: if $u_0$
and $u_1$ represent the same image up to a contrast
change, namely $u_0=f_0(u)$ and $u_1=f_1(u)$, where $f_0$ and $f_1$ are two increasing functions.
Under the assumption that $u_0$ and $u_1$ share the same geometry, we shall consider two alternative axioms, both natural in this context.



\medskip
\noindent \textit{Axiom 1:} if the intermediate image between $u_0$ and $u_1$ is an image $u$, there exists $f$  such that $u_0=f(u)$ and $u_1=f^{-1}(u)$.
\medskip

\noindent or,

\medskip
\noindent \textit{Axiom 2:} the  intermediate image between $u_0$
and $u_1$ is a direct average function of $f_0(u)$ and $f_1(u)$. The average function is  independent of $f_0$ and $f_1$.

\medskip
This axiom means that there exists a function $F:[0,1]^2\rightarrow [0,1]$, such that for every pair of images $(u_0,u_1)$, the inverse of their intermediate cumulative histogram will be $F(H_0^{-1},H_1^{-1})$.  It can also be interpreted by saying that the operator computing the intermediate histogram act independently on the level sets of rank $\lambda$ (see chapter~\ref{chapter_ste} for an intuitive explanation of this property).

A third general axiom can possibly be added, concerning the simple case where $u_0$ and $u_1$ differ by a constant.

\medskip
\noindent \textit{Axiom 3:}   If $u_1=u_0+\lambda$, then
$\ph_0(u_0)=\ph_1(u_1)=u_0+\frac{\lambda}{2}$.
\medskip

\noindent In the example  presented above (Figure~\ref{moyenne.fig}), this axiom  yields
exactly the only sound midway histogram. 



We shall now examine the axiom 1 and the pairs of axioms 2 \& 3, in turn.
Each choice leads to a different midway equalization method.

\medskip
{\bf Analysis of Axiom 1 : the Root method.}
%\noindent Notice that axiom $2$ implies the axiom $1$. Indeed, if
%$u_1=u_0+\lambda$, we can write
%$u_1=(u_0+\frac{\lambda}{2})+\frac{\lambda}{2}$ and
%$u_0=(u_0+\frac{\lambda}{2})-\frac{\lambda}{2}$. Then the
%intermediate image  is $u_0+\frac{\lambda}{2}$.
 Assume that $u_0$ and $u_1$ come from the same image (this intermediate image is unknown), up to two contrast changes $f$ and $f^{-1}$. 
 The function $f$ is unknown, but satisfies the equality $ H_0 \circ f =H_1 \circ f^{-1}$, thus
$${H_0}^{-1} \circ H_1=f \circ f.$$

It follows that the general method following this axiom consists precisely in
building an increasing function $\psi$ such that $\psi \circ \psi
= {H_0}^{-1} \circ H_1$ and replacing $u_1$ by  $\psi(u_1)$ and
$u_0$ by ${\psi}^{-1}(u_0)$.
In the next section, we will see that such a function always exists, and how to construct it. We can argue  that if such a \textit{square iterative
root} exists, it will be closer to the identity than the function
${H_0}^{-1} \circ H_1$. Consequently, the changes performed on
both images by using $\psi$ or $\psi^{-1}$ will be smaller than
the one which is performed on only one image by specification.

\medskip
{\bf Analysis of Axioms 2 \& 3: the Inverse Average Equalization.}
 \noindent If the axioms $2$ and $3$ are satisfied, the
axiom $3$ implies that:
$$\forall x,\; \forall \lambda,\;F(x+\lambda,x)=x+\frac{\lambda}{2},$$
which yields  $F(x,y)=\frac{x+y}{2},\; \forall (x,y) \in [0,1]^2$.
 Thus, we have:
\begin{proposition}
The only intermediate image between $f_0(u)$
and $f_1(u)$ satisfying the axioms $2$ and $3$ is
$\frac{f_0(u)+f_1(u)}{2}$. 
\end{proposition}

\noindent The intermediate cumulative histogram constructed in this way is
\begin{eqnarray*}
H_{\textit{new}}&=&H \circ (\frac{f_0+f_1}{2})^{-1}=(\frac{f_0 \circ H^{-1}+f_1 \circ H^{-1}}{2})^{-1}\\
&=&(\frac{{H_0}^{-1}+{H_1}^{-1}}{2})^{-1}.\\
\end{eqnarray*}
The general method deduced from the previous proposition
will be called \textit{Inverse Average Equalization} (i.a.e.). The
contrast changes involved are:
$$\ph_0=\frac{id+{H_0}^{-1}\circ H_1}{2} \;\text{  and  }\; \ph_1=\frac{id+{H_1}^{-1}\circ H_0}{2}.$$

\medskip 
{\bf The Cox et al. algorithm}.
\noindent We shall now prove that the Inverse Average Equalization is
nothing but a continuous version of the discrete algorithm
proposed by Cox et al. in \cite{Cox}. The dynamic histogram
warping (DHW) presented in \cite{Cox} uses dynamic programming for
matching directly histogram values. Only one-to-one, one-to-many
or many-to-one mappings are allowed.

Suppose the range $[n_0,\tilde{n}_0]$ in $u_0$ has been
mapped with the intensity $n_1$ in $u_1$. Cox et al. propose to
simultaneously warp both histograms by replacing the previous grey
levels in both images by \footnote{this is not specified in their
paper, but it can be deduced from their code.}
$n_0=\frac{n_1+\frac{n_0+\tilde{n}_0}{2}}{2}$. Now, we have
$H_0(\frac{n_0+\tilde{n}_0}{2}) \simeq H_1(n_1)$, then:
$$n_0 \simeq \frac{n_1+{H_0}^{-1} \circ H_1 (n_1)}{2}.$$
In the same way, if $n \in [n_0,\tilde{n}_0]$, we have
$$n_0 \simeq \frac{n+{H_1}^{-1} \circ H_0 (n)}{2}.$$
It follows that $u_0$ is actually composed by
$\psi_0=\frac{{H_0}^{-1}+{H_1}^{-1}}{2} \circ H_0$ and $u_1$ by
$\psi_1=\frac{{H_0}^{-1}+{H_1}^{-1}}{2} \circ H_1$. Thus, the
intermediate histogram constructed this way is exactly
$$H=H_0 \circ {\psi_0}^{-1} = H_1 \circ {\psi_1}^{-1} =(\frac{{H_0}^{-1}+{H_1}^{-1}}{2})^{-1}.$$
Their algorithm is  a very efficient way to perform this
transformation on discrete data.










\section{Root. Theory. }

Let us see how the   method deriving from Axioms 2 \& 3 can be
performed in practical terms. Any solution $\psi$ of the
functional equation
$$\psi^{(\circ N)}=\ph$$
is called an \textit{iterative root} of the given function $\ph$.
Iterative roots are studied in \textit{Iterative functional
equations} by M.Kuczma, B.Choczewski and R.Ger \cite{Kuc}. Some of
the results of this section (propositions~\ref{prop_construction}, \ref{prop_c1_condition} and \ref{prop_c1_existence}) can be found in this book. We reformulate them here in the image processing numerical context, which is new.

In our case, we are just interested in square iterative roots
($N=2$) of functions from an interval into itself. We shall even
limit ourselves to the case where $\ph: [0,1] \rightarrow [0,1]$
is an increasing homeomorphism \footnote{The case  ``$\Phi: [0,a]
\rightarrow [0,a]$'' with $a>0$ can be deduced from the previous
one via the function $\ph(x)=\frac{1}{a}\Phi(ax)$. If $\psi$ is a
s.i.r. of $\ph$, then $x \rightarrow a\psi(\frac{x}{a})$ is a
s.i.r. of $\Phi$.}. Consequently, we call \textit{square iterative
root} (s.i.r.) of  $\ph$ every increasing homeomorphism $\psi:
[0,1] \rightarrow [0,1]$ such that $\psi \circ \psi = \ph$. For
example, $x \rightarrow x^{\sqrt{2}}$ is a square iterative root
of $x \rightarrow x^{2}$.




This definition leads to several questions, as the existence or
uniqueness of such an iterative root. Let us start with some
remarks.

\begin{itemize}

\item We already know that $\varphi(0)=\psi(0)=0$ and $\varphi(1)=\psi(1)=1$.
If $\varphi$ and $\psi$ are $C^1$ functions, then
$\ph'(0)=\psi'(0)^2$ and  $\ph'(1)=\psi'(1)^2$.

\item If $x_o \in ]0,1[$ is such that $\ph(x_o)=x_o$, then $\psi_{[0,x_o]}$
(respectively $\psi_{]x_o,1]}$) just relies on $\ph_{[0,x_o]}$
(resp. $\ph_{]x_o,1]}$). Then, we can limit the search for a
s.i.r. to the case $\ph(x)>x \text{ for all } x \text{ in }
]0,1[$. The case  $\ph(x)<x$ follows immediately: $\psi$ is a
s.i.r of $\ph$ iff $1-\psi(1-x)$ is a s.i.r. of $1-\ph(1-x)$.

\item If $\ph(x)>x$ on $]0,1[$, and if $\psi$ is a s.i.r. of $\ph$, then $x<\psi(x)<\ph(x) \; \forall x \in ]0,1[$.
\end{itemize}


\subsection{Construction of a Square Iterative Root}
\begin{proposition}
Let $(x_0,y_0)$ be a point strictly located between the graph of
$\ph$ and the line $y=x$. Then $\ph$ has an infinite amount of
square iterative roots $\psi$ such that $\psi(x_0)=y_0$.
\label{prop_construction}
\end{proposition}

\begin{proof}
%\psfrag{x-3}{$x_{-3}$} \psfrag{x-2}{$x_{-2}$}
%\psfrag{x-1}{$x_{-1}$} \psfrag{x0}{$x_0$} \psfrag{x1}{$x_{1}$}
%\psfrag{x2}{$x_{2}$} \psfrag{x3}{$x_{3}$} \psfrag{y-3}{$y_{-3}$}
%\psfrag{y-2}{$y_{-2}$} \psfrag{y-1}{$y_{-1}$} \psfrag{y0}{$y_0$}
%\psfrag{y1}{$y_{1}$} \psfrag{y2}{$y_{2}$} \psfrag{y3}{$y_{3}$}
%\psfrag{y0'}{${y_0}'$} \psfrag{y2'}{${y_2}'$}
%\psfrag{x1'}{${x_1}'$} \psfrag{x3'}{${x_3}'$}
%\psfrag{psi(x0)}{$\psi(x_0)$} \psfrag{ps-(x0)}{${\psi}^{-1}(x_0)$}
%\psfrag{ps-3(x0)}{${\psi}^{-3}(x_0)$}
%\psfrag{psi3(x0)}{${\psi}^{\circ (3)}(x_0)$}
%\psfrag{ph(x0)}{$\ph(x_0)$} \psfrag{ph-(x0)}{${\ph}^{-1}(x_0)$}
%\psfrag{y=x}{$y=x$}

\begin{figure}[h]
\centerline{\includegraphics[width=8.0cm]{iterations}}
\caption{\em Construction of a Square Iterative Root.}
\label{construction.fig}
\end{figure}



Let $(x_0,y_0)$ be a point in the area between the graph of $\ph$ and the line $y=x$. We want to
build a square iterative root $\psi$ of $\ph$, such that
$\psi(x_0)=y_0$. We are already able to know the sequences of
positive and negative iterates of $x_0$ by $\psi$. Indeed, we know
$\ph$ and
$$\forall n \in {\mathbb{Z}}, \; {\psi}^{\circ (2n)}(x_0)={\ph}^{\circ n}(x_0), \; {\psi}^{\circ (2n+1)}(x_0)={\ph}^{\circ n}(y_0).$$
Next, since we want $\psi$ to be an increasing function, its graph
must remain in the sequence of rectangles defined by those
iterations (Figure~\ref{construction.fig}). Now, we remark that any continuous and increasing function $\psi$ defined on $[x_0,\psi(x_0)]$ and satisfying $\psi(\psi(x_0))=\ph(x_0)$, can be prolonged on $[0,1]$ as a square iterative root of $\ph$. Indeed, ${\psi}^{-1}$ is defined and continuous on
$[\psi(x_0),\psi^{\circ (2)}(x_0)]$ (because $\psi$ is an
homeomorphism from  $[x_0,\psi(x_0)]$ onto $[\psi(x_0),\psi^{\circ
(2)}(x_0)]$). Therefore $\psi=\ph \circ {\psi}^{-1}$ can also be
defined continuously on $[\psi(x_0),\psi^{\circ(2)}(x_0)]$. The
continuity of $\psi$ at $\psi(x_0)$ is ensured by the
construction. This process implies that the function $\psi$
constructed is continuous on  $]0,1[$, and by extension, on $[0,1]$.

\end{proof}


In order to illustrate the process described above, we have used
it to construct some square iterative roots of the function $x
\rightarrow \sqrt{x}$ in Figure~\ref{ex_racines.ps}. For each choice
of the point $(x_0,y_0)$, $\psi$ has been defined as an affine
function on the set $[x_0,y_0]$.

\begin{figure}[h]
\centerline{\includegraphics[width=8.0cm]{ex_racines2}}
\caption{\em Some square iterative roots of the function $x \rightarrow
x^{\frac{1}{2}}$. The darker curve is the explicit s.i.r. $x
\rightarrow x^{\sqrt(\frac{1}{2})}$. } 
\label{ex_racines.ps}
\end{figure}




Our problem is now to single out a root among all possible square
roots, which would be more judicious than the other ones. In what
follows, we analyze the properties of all square roots and discuss
criteria for a best choice.

\subsection{Properties of Square Iterative Roots}



\subsubsection{Oscillations}

\begin{proposition}
Let $\ph: [0,1] \rightarrow [0,1]$ be an increasing
homeomorphism, and $\psi_0$ and $\psi_1$ be two square iterative
roots of $\ph$. Then, for every positive number $\eps$, there is
an infinite amount of $x \in [0,\eps]$ (likewise in $[1-\eps,1]$)
such that $\psi_0$ and $\psi_1$ cross each other at $x$.
\end{proposition}
\begin{proof}
Assume that there exists a positive number $\eps$ such that
$$\forall x \in ]0,\eps[,\; \psi_0(x)\geq \psi_1(x).$$
 Since $\psi_0$ is continuous and $\psi_0(0)=0$, there exists a
positive number $\delta$ such that  $\forall x \in ]0,\delta[,\;
\psi_0(x) \in ]0,\eps[$. Thus,
$$\forall x \in ]0,\inf(\delta,\eps)[,\; \psi_0 \circ \psi_0(x) \geq \psi_1 \circ \psi_0(x) \geq \psi_1 \circ \psi_1(x) $$
(the second inequality comes from the fact that $\psi_1$ is
increasing). Now, the extremal terms must be equal to $\ph(x)$.
Then, all the terms are equal, $\psi_0$ and $\psi_1$ take the same
values on $[0,\inf(\delta,\eps)]$ and consequently $\psi_0$ and
$\psi_1$ are equal everywhere (see  the explicit construction).

\end{proof}


\subsubsection{Self-similarity}

\begin{proposition}
Let $\ph$ be linear on $[0,\eps]$ ($\eps > 0$), with a slope
$a>1$. Then, if $\psi$, a square iterative root of $\ph$, is $C^1$ on a
neighborhood of $0$, $\psi(x)=\sqrt{a}x$ on  $[0,\eps]$.
\end{proposition}

\begin{proof}
Let $\psi$ be a square iterative root of $\ph$ and $x_0$ be a
point in $[0,\eps]$. We define the sequences $(x_n)$, $(y_n)$ and
$(r_n)$ as $(x_n)=\psi^{-n}(x_0)$, $(y_n)=\psi(x_n)$, and
$$r_n=\frac{y_n-y_{n+1}}{x_n-x_{n+1}},$$
where $\psi^{-n}$ denotes the $n^{th}$ iterate of $\psi^{-1}$.
Then, the sequences $r_{2n}$ and $r_{2n+1}$ remain constant: odd
rectangles (resp. even) are similar (Figure~\ref{autosimilaire.fig}).
Moreover, the graph of $\psi$ in two odd (resp. two even)
rectangles is the same, up to a similarity. The points
$(x_{2n},y_{2n})$ (resp. the points $(x_{2n+1},y_{2n+1})$) are all
on a line going through $(0,0)$. Consequently, $\psi$ is a $C^1$
function in a neighborhood of $0$ if and only if both lines are
the same for all $x_0$ in $[0,\eps]$, {\em i.e.} iff $\psi(x)=\sqrt{a}x$ on $[0,\eps]$. We can also
add, following the same argument, that this function is the only
square root which is concave near $0$.


%\psfrag{x0}{$x_0$} \psfrag{x1}{$x_{1}$} \psfrag{x2}{$x_{2}$}
%\psfrag{x3}{$x_{3}$} \psfrag{y0}{$y_0$} \psfrag{y1}{$y_{1}$}
%\psfrag{y2}{$y_{2}$} \psfrag{y3}{$y_{3}$} \psfrag{y=ax}{$y=ax$}
%\psfrag{y=x}{$y=x$}
\begin{figure}[h]
\centerline{\includegraphics[width=8.0cm]{autosimilaire}}
\caption{\em Behaviour of a square iterative root in the linear case.} \label{autosimilaire.fig}
\end{figure}

\end{proof}


\subsection{$C^1$ Square Iterative Roots}

We focus our attention on the case when $\ph$ is a $C^1$ function
on the open set $]0,1[$, and when $\ph'$ does not reach $0$.

\begin{proposition}
Let $\ph$ be a  $C^1$ increasing homeomorphism on $]0,1[$, such
that $\ph'>0$. There exists an infinite amount of
square iterative roots $\psi$, which are $C^1$ functions on the
open set $]0,1[$. Moreover, in the previous construction, one only
has to build $\psi$ as a $C^1$ function on $]x_0,y_0[$ such that
$\psi'({y_0}^-)=\frac{\ph'(x_0)}{\psi'({x_0}^+)}$ to be sure that
$\psi$ is a $C^1$ function on $]0,1[$.
\label{prop_c1_condition}
\end{proposition}
\begin{proof}
The derivability and continuity of the derivative on the open sets are easy to prove. Let see what happens at the edges of the intervals.
Let $\psi$ be a square iterative root of $\ph$ going through
$(x_0,y_0)$. Then,
$$\forall y \in ]y_0,\ph(x_0)],\;\exists x \in ]x_0,y_0],\; y=\psi(x).$$

So, if $y_0<y \leq \ph(x_0)$:
$$\frac{\psi(y)-\psi(y_0)}{y-y_0}=\frac{\ph(x)-\ph(x_0)}{x-x_0}\frac{x-x_0}{\psi(x)-\psi(x_0)}.$$
$\psi^{-1}$ is continuous, then $x \rightarrow x_0$  when $\psi(x)
\rightarrow \psi(x_0)$. Thus
$$lim_{y \stackrel{>}{\rightarrow} y_0}\frac{\psi(y)-\psi(y_0)}{y-y_0}=\frac{\ph'(x_0)}{\psi'({x_0}^+)}.$$

\end{proof}


\begin{proposition}
There exists a $C^1$ increasing homeomorphism $\ph: [0,1]
\rightarrow [0,1]$, with no $C^1$ square iterative root on $[0,1]$.
\label{prop_c1_existence}
\end{proposition}
\begin{proof}
It is clear that we can construct a $C^1$ increasing function $\ph$, satisfying:


$$\ph(x)=\left\{
\begin{array}{ll}
2x &  \text{ on } [0,\frac{1}{3}] \\
1+\frac{1}{3}(x-1)& \text{ on } [\frac{\sqrt{2}}{3},1].
\end{array}
\right.
$$

Let $\psi$ be a $C^1$ square iterative root of $\ph$. Then, $\psi$ must at least
be linear on $[0,\frac{1}{3}]$ and $[\frac{\sqrt{2}}{3},1]$, with
respective slopes$\sqrt{2}$ and $\sqrt{\frac{1}{3}}$ (see Figure~\ref{ctre_ex.fig}).
Hence, $\psi(\frac{\sqrt{2}}{3})=1+\frac{\sqrt{2}-3}{3 \sqrt{3}}$.
However, one iteration of the point
$(\frac{1}{3},\frac{\sqrt{2}}{3})$ give
$\psi(\frac{\sqrt{2}}{3})=\frac{2}{3}$. It follows that there is
no $C^1$ square iterative root of $\ph$.

%\psfrag{1/3}{$\frac{1}{3}$} \psfrag{2/3}{$\frac{2}{3}$}
%\psfrag{rac(2)/3}{$\frac{\sqrt{2}}{3}$} \psfrag{1}{$1$}
\begin{figure}[h]
\centerline{\includegraphics[width=8.0cm]{ctre_ex}}
\caption{\em Example of $C^1$ function $\ph$ without any $C^1$  square iterative root.}
\label{ctre_ex.fig}
\end{figure}
\end{proof}

The previous proposition implies that it is useless in general to search for
some $C^1$ square iterative root. In what follows, we will see
that the concavity of $\ph$ does not imply the existence of a
concave square root. However, it leads us to an interesting
property, which will help us to determine what kind of reasonable
criterion we can use, in order to choose among all these
solutions.




\subsection{Behavior towards the concavity of $\ph$}
The previous example is also suitable for showing that there is not
always a concave square root if $\ph$ is a concave function.
Indeed, in the previous construction, we can ensure the concavity
of $\ph$. Next, if $\psi$ is a concave s.i.r. of $\ph$, $\psi$ has
to be linear on $[0,\frac{1}{3}]$ and $[\frac{\sqrt{2}}{3},1]$
(see the part about self-similarity). Now, we have:
\begin{proposition}
If $\ph$ is concave near $0$, and if there exists a concave
square root $\psi$ for $\ph$, then $\psi$ is unique.
\label{prop_concave}
\end{proposition}
\begin{proof}
Assume that $\psi_0$ and $\psi_1$ are two distinct concave square
roots of $\ph$. Let $a>0$ be such that $\ph$ is concave on
$[0,a]$. It exists $x_0$ in $[0,a]$ such that
$y_0=\psi_0(x_0)>y'_0=\psi_1(x_0)$. Let  $x_k={\psi_0}^{-k}(x_0)$
and $y_k={\psi_0}^{-k}(y_0)$ be the iterates of $x_0$ and $y_0$ by
${\psi_0}^{-1}$, and $(x_k')$, $(y_k')$ those of $x_0$ and $y_0'$ by
${\psi_1}^{-1}$. All these values are of course in $[0,a]$. Then,
for each positive integer $k$:
$$\frac{\ph(y_{2k+2})-\ph({y'}_{2k+2})}{y_{2k+2}-{y'}_{2k+2}} \leq \frac{\ph(x_{2k+2})}{x_{2k+2}}.$$
Indeed, the first term is the slope of the chord between
$(y_{2k+2},\ph(y_{2k+2}))$ and $(y'_{2k+2},\ph(y'_{2k+2}))$, and
the second one is the slope of the chord between
$(x_{2k+2},\ph(x_{2k+2}))$ and $(0,0)$. Now,
$0<x_{2k+2}<y'_{2k+2}<y_{2k+2}$, thus the previous inequality
follows from the concavity of $\ph$. Consequently:
$$\forall k \geq 0,\;\frac{y_{2k}-{y'}_{2k}}{x_{2k}} \leq \frac{y_{2k+2}-{y'}_{2k+2}}{x_{2k+2}} .$$
The sequence $(\frac{y_{2k}-{y'}_{2k}}{x_{2k}})$ is increasing,
and then is always larger than its first term, which is positive.
Now, if $\psi_0$ and $\psi_1$ are concave, they have  right
derivatives at $0$, which must be equal to $\sqrt{\ph'(0)}$. It
follows that the sequences $ \frac{y_{2k}}{x_{2k}}$ and
$\frac{y'_{2k}}{x_{2k}}$ must have the same limit, and thus
$$\frac{y_{2k}-{y'}_{2k}}{x_{2k}} \stackrel{k \rightarrow
\infty}{\rightarrow} 0,$$ which contradicts the previous property.


%\psfrag{x0}{$x_0$} \psfrag{x1}{$x_{1}$} \psfrag{x2}{$x_{2}$}
%\psfrag{x3}{$x_{3}$} \psfrag{y0}{$y_0$} \psfrag{y1}{$y_{1}$}
%\psfrag{y2}{$y_{2}$} \psfrag{y3}{$y_{3}$} \psfrag{y=x}{$y=x$}
%\psfrag{y=ax}{$y=ax$} \psfrag{y0'}{${y_0}'$}
%\psfrag{y2'}{${y_2}'$} \psfrag{x1'}{${x_1}'$}
%\psfrag{x3'}{${x_3}'$}
\begin{figure}[h]
\centerline{\includegraphics[width=8.0cm]{concave_unicite}}
\caption{\em Iterations of ${\psi_0}^{-1}$ and ${\psi_1}^{-1}$.}
\label{concave_unicite.fig}
\end{figure}
\end{proof}


\section{The Choice of the Square Root}

\subsection{How to choose the ``best'' square iterative root ?}
Let $u_0$ and $u_1$ be the images, whose histograms have to be
brought closer. According to the above considerations, we have to
find a continuous function $\psi$, for which $\psi \circ \psi =
{H_0}^{-1} \circ H_1$. Then, $\psi(u_1)$ and ${\psi}^{-1}(u_0)$
will have the same histogram, which will be an intermediary
between $H_0$ and $H_1$. Now, we have seen that there is an
infinite amount of such functions. Even if we search $C^1$ square
roots on the open set $]0,1[$ (in case ${H_0}^{-1} \circ H_1$ is a
$C^1$ function itself), there exists an infinite amount of
solutions. On the other hand, there exist $C^1$ homeomorphisms on
$[0,1]$, none of whose square root is $C^1$ on the closed set
$[0,1]$, and one can find concave functions deprived of  any
concave s.i.r. Thus, no usual structural criterion seems at hand
to single out a right square root.

Let us come back to proposition~\ref{prop_concave}: if
a function $\ph$ has a concave square root $\psi$, then this
square root is the ``right'' one, the one we look for. Look at
Figure~\ref{ex_racines.ps}: it shows a variety of bad solutions for
$\ph(x)=\sqrt{x}$. We have seen that those solutions are
oscillating around each other, whatever $\ph$ is
(\textit{Prop.4}). In this case, each solution oscillates around
the concave one $\psi(x)=x^{\frac{1}{\sqrt{2}}}$, and changes its
convexity an infinite amount of times near  $0$ (and near $1$).
Thus, the concave solution is the one which ``oscillates''  as
little as possible. In the next section, we propose a variational
criterion based on a non oscillation requirement.




\subsection{Behavior towards the noise}

Assume that we have a uniform quantization noise $b$ added to $u$
($-0.5 \leq b \leq 0.5$). Then, after a composition by a function
$\psi: [0,255] \rightarrow [0,255]$:
$$\psi(u+b) \simeq \psi(u)+\psi'(u)b\;\;\text{(first order)}.$$
 The noise becomes $\psi'(u)b$, which depends on the grey
  level of the point. It follows that the noise is not uniform any more.
  The more the variations of $\psi'$ are important, the less the noise
  is uniform. Now, the  uniformity of  noise requirement is the most
  relevant one in applications. Therefore, if we want to perform
  a histogram modification, we must try to choose $\psi'$ with the smallest
   variations around its mean.
However, in our case, the contrast change used for the second image involves $\psi^{-1}$, we should consequently try to minimize:
$\int_0^1 (\psi'-1)^2+\int_0^1 (\frac{1}{\psi'(\psi^{-1})}-1)^2$.
Now, we can write 
$$\int_0^1 (\psi'-1)^2+\int_0^1 (\frac{1}{\psi'(\psi^{-1})}-1)^2=\int_0^1 (\psi'-1)^2+\int_0^1 (\frac{1}{\psi'}-1)^2\psi'=\int_0^1 (\psi'-1)^2(1+\frac{1}{\psi'}).$$
In practice, we observed that this minimization gives almost the same results as the minimization of $\int_0^1 (\psi'-1)^2$ (because $\psi'$ is close to $1$). Therefore, we shall make the approximation:
\\
\textit{Principle:} The square root which alters the least the uniformity of the quantization noise must minimize $\int_0^1 (\psi' -1)^2$.



\subsection{Results on examples}
We have seen and explained why the best square root for our
problem is the one which ``oscillates'' as little as possible. On
the other hand, a ``good'' solution, from the quantization noise
point of view, has to make $\int {\psi '}^2$ as small as possible.
In order to understand the link between both criteria, we
made the following experiment: let $\ph$ be a concave function,
whose concave square root $f$ we know. First, we compute the value
$y_0$ in $[\frac{1}{2},\ph(\frac{1}{2})]$, which minimizes the
discrete version of $\int {\psi '}^2$ (computed on the sequence of
iterates of the point $(\frac{1}{2},y_0)$). Next, we compare this
value to the value of the concave s.i.r. $f$ at $\frac{1}{2}$. The
results are in the following table (we use $30$ iterations of the
point, and the precision on $y_0$ is better than $10^{-3}$):

\begin{center}
\begin{tabular}[]{|l|c|c|r|}
\hline
$\ph$ & $f$ & $y_0$ & $f(\frac{1}{2})$ \\
\hline
$\sqrt{x}$ & $x^{\frac{1}{\sqrt{2}}}$ & $0.612553$ & $0.612547$\\
\hline
$f \circ f$ & $ln((e-1)x+1)$ & $0.6184$ & $0.62011$ \\
\hline
$x^{\frac{1}{4}}$ & $x^{\frac{1}{2}}$ & $0.7074$ & $0.7071$ \\
\hline
$f \circ f$ & $\frac{4}{\Pi} \arctan{x}$ & $0.5891$ & $0.5903$ \\
\hline
\end{tabular}
\end{center}

Thus, for those functions, the value $y_0$ computed is quite close
to the value at $\frac{1}{2}$ of the solution we would like to
construct.

\subsection{S.I.R. Algorithm}

Our algorithm is the following one (for the description, $\ph>Id$)
: first, let $x_o=\frac{1}{2}$. We choose $\psi(\frac{ 1}{2})$ in
the open set $]\frac{1}{2},\ph(\frac{1}{2})[$, such that the
sequence $({\psi}^{\circ k}(\frac{1}{2}),{\psi}^{\circ (k+1)}(\frac{1}{2}))_{k \in \Z}$ oscillates as little as possible. 
Now, the values of $\psi$ at all iterates (negative or positive) 
of $\frac{1}{2}$ are fixed (see Figure~\ref{construction.fig}). We 
already know that, once $\psi(\frac{1}{2})$ is chosen, the construction of $\psi$ as a continuous and increasing function on $[\frac{1}{2},\psi(\frac{1}{2})]$ is enough to 
completely define $\psi$ as continuous and increasing on $[0,1]$. Indeed, the
values of $\psi$ at the other points of the set $[0,1]$ will be
defined as iterates of those on $[\frac{1}{2},\psi(\frac{1}{2})]$.
So, we suggest to construct $\psi$ on a dichotomic subdivision of
$[\frac{1}{2},\psi(\frac{1}{2})]$.

The first step consists in the construction of $\psi(\frac{1}{2})$. This value must lie in $]\frac{1}{2},\ph(\frac{1}{2})[$, and is actually chosen in the open set such that the iterates of $(\frac{1}{2},\psi(\frac{1}{2}))$ minimize the discrete version of $\int \psi'^2$. We follow the same rule to construct by dichotomy the values of $\psi$ on $]\frac{1}{2},\psi(\frac{1}{2})[$.  


Then, for $k=0 \dots K$ ($k$ is the current depth of the dichotomy, and $K$ the final one), for each
$(\eps_0,...,\eps_k) \in {\{0,1\}}^{k+1}$, let
$$x_0=(\sum_{i=0}^{k} \eps_i {\left(\frac{1}{2}\right)}^i) (\psi(\frac{1}{2})-\frac{1}{2})+\frac{1}{2}.$$
We want to construct $\psi(x_0)$. Now, the value chosen must remain coherent with the points previously constructed in the dichotomy. Consequently, let $X_0$ and $X_1$ be the points on both sides of $x_0$ which have already been dealt with in previous steps of the dichotomy:
$$X_0=(\sum_{i=0}^{k-1} \eps_i {\left(\frac{1}{2}\right)}^i)(\psi(\frac{1}{2})-\frac{1}{2})+\frac{1}{2}.$$
$$X_1=X_0+(\frac{1}{2})^{k-1}(\psi(\frac{1}{2})-\frac{1}{2}).$$
$x_0$ is the middle of $[X_0,X_1]$ (see Figure~\ref{algo.fig}) and
$\psi(X_0)$ and $\psi(X_1)$ are already fixed.
Therefore, we choose $\psi(x_0)$ in the open set $]\psi(X_0),\psi(X_1)[$, such that the sequences of the iterates of $(x_0,\psi(x_0))$ minimizes the discrete version of $\int {\psi'}^2$.

%\psfrag{x0}{$x_0$} \psfrag{X0}{$X_0$} \psfrag{X1}{$X_1$}
%\psfrag{psi(x0)}{$\psi(x_0)$}
%\psfrag{(1/2,psi(1/2))}{$(\frac{1}{2},\psi(\frac{1}{2}))$}
%\psfrag{(psi(1/2),phi(1/2))}{$(\psi(\frac{1}{2}),\phi(\frac{1}{2}))$}
%\psfrag{phi(1/2)}{$\phi(\frac{1}{2})$}
%\psfrag{psi(X0)}{$\psi(X_0)$} \psfrag{psi(x0)}{$\psi(x_0)$}
%\psfrag{psi(X1)}{$\psi(X_1)$}
\begin{figure}[h]
\centerline{\includegraphics[width=8.0cm]{algo}}
\caption{\em Construction of $\psi$ on
$[\frac{1}{2},\psi(\frac{1}{2})]$.} \label{algo.fig}
\end{figure}


Thanks to this algorithm, we are able to construct a ``good''
square root for any function $\ph$ which does not cross the line
$y=x$ on $]0,1[$. If a function $\ph$ is an increasing
homeomorphism on $[0,1]$, we can use the same algorithm on each
interval $[x_0,x_1]$ such that $\ph(x_0)=x_0$, $\ph(x_1)=x_1$, and
$\ph(x) \neq x $ on $]x_0,x_1[$. This way, we can compute a
``good'' square root for any increasing homeomorphism of $[0,1]$,
and next for any increasing homeomorphism of a segment $[a,b]$.


\paragraph{Cost of the algorithm:} if $\ph$ is a discrete increasing function
on $[0,1]$, defined on a regular grid of $N$ bins, we know that the depth of the dichotomy should not exceed $\log_2(N)$. Moreover, the cost of construction of the image of a point of the dichotomy is constant. Hence, the algorithm is $O(N)$.








\section{Comparison of the methods. Experiments and applications}

\subsection{Theoretical Comparison of the methods}

First, let us remark that if $u_1=\lambda u_0$ then the s.i.r.
equalization builds as an intermediate image between $u_1$ and
$u_0$ the image $u=\sqrt{\lambda} u_0$, whereas the inverse
average equalization builds $u= \frac{\lambda+1}{2} u_0$. Now,
when $\lambda$ is close to $1$, $\sqrt{\lambda} \simeq
\frac{\lambda+1}{2}$, thus the results obtained with both methods
are quite similar. By the same arguments, we can show in a more
general way that when ${H_0}^{-1} \circ H_1$ is close to $Id$, the
results obtained with the two methods are almost identical.

However, the i.a.e. method is more conservative since it
maintains the mean of the two images constant. Indeed, if we note
$\bar{u_i}$ the mean of $u_i$, we have:
$$\bar{u_i}=\int_0^1 xh_i(x)dx = 1-\int_0^1 H_i(x)dx=\int_0^1 {H_i}^{-1}(x)dx.$$ Then the average
of the two means is:
$$\frac{\bar{u}_0+\bar{u}_1}{2}=\int \frac{{H_0}^{-1}+{H_1}^{-1}}{2}$$
and remains the same before and after the inverse average
equalization. This is not the case for the s.i.r method: if
$u_1=\lambda u_0$, then the average of their means before the
transformation is $\frac{\lambda+1}{2} \bar{u}_0$ and becomes
 $\sqrt{\lambda} \bar{u}_0$. This conservation principle leads us to 
assume that the use of the i.a.e. is maybe more relevant in iterative algorithms.


\subsection{Results on a stereo pair}


The top of Figure~\ref{toulouse} shows a pair of aerial images in the region
of Toulouse. Although the angle variation between both views is
small, and the photographs are taken at nearly the same time, we
see that the lightning conditions vary significantly (the radiometric differences can also come from a change in camera settings).
The second line shows the result of the specification of the
histogram of each image on the other one. The third line
shows both images after equalization with the s.i.r. method,
and the last line after inverse average equalization. 
Figure~\ref{4histo.ps} represents the cumulative histograms of the
images before and after the different transformations. We can see
that the histograms of the images after both equalizations are
close and really intermediate between the previous ones. All the
comparisons which follow are made between specification and the
s.i.r. equalization but qualitative results would be very similar
and conclusions identical with the inverse average equalization.

\begin{figure}
\centerline{
\subfigure[\em Stereo pair before equalization.]{
\includegraphics[width=8.0cm]{img11}
\includegraphics[width=8.0cm]{img12}}}
\centerline{
\subfigure[\em Same images after specification of their histograms on each other (left: the histogram of the first image has been specified on the second, and right: the histogram of the second image has been specified on the first).]{
\includegraphics[width=8.0cm]{img11sp12}
\includegraphics[width=8.0cm]{img1211}}}
\centerline{
\subfigure[\em Stereo pair after equalization via the s.i.r. algorithm.]{
\includegraphics[width=8.0cm]{img11rac}
\includegraphics[width=8.0cm]{img12rac}}}
\centerline{
\subfigure[\em Stereo pair after equalization via the i.a.e. algorithm.]{
\includegraphics[width=8.0cm]{img11ds}
\includegraphics[width=8.0cm]{img12ds}}}
\caption{Two pieces of aerial images of a region of Toulouse, before and after different kind of midway equalizations.}
\label{toulouse}
\end{figure}


If we scan  some image details, as illustrated on Figure~\ref{extraits}, 
the damages caused by a direct specification
become obvious. Let us specify the darker image on the brightest
one. Then the information loss, due to the reduction of dynamic
range, can be detected in the brightest areas. Look at the roof of
the bright building in the top left corner of the image
(first line of Figure~\ref{extraits}): the chimneys project horizontal shadows
on the roof. In the specified image, these shadows have almost
completely vanished, and we cannot even discern the presence of a
chimney anymore. In the same image after s.i.r. equalization, the
shadows are still entirely recognizable, and their size reduction
remains minimal. The second line of Figure~\ref{extraits} illustrates the same
phenomenon, observed in the bottom center of the image. The
structure present at the bottom of the image has completely
disappeared after specification and remains visible after midway
equalization. These examples show  how visual information can be
lost by specification and how  midway algorithms reduce
significantly this loss.



\begin{figure}[h]
\centerline{\includegraphics[width=8.0cm]{4histo}}
\caption{\em  Cumulative histograms of the stereo pair (Figure~\ref{toulouse}) before and after each transformation. 1 and 4: Cumulative histograms of the first and second image, 2: Common histogram after i.a.e., and 3: Common histogram after s.i.r. equalization. } \label{4histo.ps}
\end{figure}


\begin{figure}[h]
\centerline{\includegraphics[width=6.0cm]{histo_et_rac_1112}}
\caption{\em The first curve (1) represents ${H_0}^{-1} \circ H_1$ for the stereo pair and the second one (2) is the square iterative root computed in this case. }
\label{histo_et_rac_1112.ps}
\end{figure}


\begin{figure}[h]
\centerline{
\subfigure{\includegraphics[width=9.0cm]{extrait1}}}
\centerline{\subfigure{\includegraphics[width=9.0cm]{extrait2}}}
\caption{\em Two extracts of the stereo pair shown on Figure~\ref{toulouse}. For each line, from left to right: in the original image, in the specified one, in the original image after midway equalization.}
\label{extraits}
\end{figure}







\subsection{Multi-Channel images}
The top of Figure~\ref{multi} shows two pieces of multi-channel images of
Toulouse. The first one is extracted of the blue channel, and the
other from the infrared channel. The second and third line of the same figure show the same images after the s.i.r. and the i.a.e. algorithm.
We can see with this example one of the limitations of the method. The multichannel images have the particularity to present contrast inversions : for instance, the trees appear to be darker than the church in the blue channel, and are naturally brighter than the church in the infrared channel. The midway equalization being limited to increasing contrast changes, it obviously cannot handle these contrast inversions. Now, the main interest of matching the histograms of such images is to make their visual comparison easier and more relevant : we can see that in spite of these contrast inversions, the results remain visually good, which underlines the robustness of the method. Moreover, we observed that the results of correlation-based stereo algorihtms on multi-channel pairs were never worse and often improved after this kind of equalization.

\begin{figure}
\centerline{
\subfigure{
\includegraphics[width=4.0cm]{toul0}
\includegraphics[width=4.0cm]{toul3}}}
\centerline{
\subfigure{\includegraphics[width=4.0cm]{toul0v3_p512}
\includegraphics[width=4.0cm]{toul3v0_p512}}}
\centerline{
\subfigure{\includegraphics[width=4.0cm]{toulds0v3_p512}
\includegraphics[width=4.0cm]{toulds3v0_p512}}}
\caption{\em First line: two images of Toulouse (blue and infrared channel).
Second line:  same images after s.i.r. equalization.
Third line: same images after i.a.e. equalization.} 
\label{multi}
\end{figure}


%\begin{figure}
%\centerline{\includegraphics[width=4.0cm]{toulsp0v3_p512}
%\includegraphics[width=4.0cm]{toulsp3v0_p512}}
%\caption{\em  Blue and infrared channel after specification on each other.} \label{multi4}
%\end{figure}









\subsection{Different shots of the same picture}
The top of Figure~\ref{meduse} shows two different snapshots of the same painting, \textit{Le Radeau de la méduse}\footnote{Musée du Louvre, Paris.},
by Théodore Géricault, found on the internet. The second one is
brighter and seems to be damaged at the bottom left. The second line
 shows the same couple after s.i.r. equalization. The results of the inverse average equalization on this example are similar. Finally, the last line of Figure~\ref{meduse} shows the difference between both images after the midway equalization. We see clear differences around the edges, due to the fact that the original images are not completely similar from the geometric point of view.

\onecolumn
\begin{figure}
\centerline{
\subfigure[{\em Two shots of the \textit{Radeau de
la Méduse}, by Géricault}]{
\includegraphics[width=8.0cm]{meduse}
\includegraphics[width=8.0cm]{meduse5}}}
\centerline{
\subfigure[{\em The same images after s.i.r. equalization.}]{
\includegraphics[width=8.0cm]{meduse1_5}
\includegraphics[width=8.0cm]{meduse5_1}}}
\centerline{
\subfigure[{\em Image of the difference between both images after the midway equalization.}]{
\includegraphics[width=8.0cm]{diff_meduse}}}
\caption{\em Two different images of the same painting and their midway versions. The difference between both images after midway equalization is mostly geometric, due to the small distortions between the images.}
\label{meduse}
\end{figure}



\subsection{Effects on the quantization noise}

Let us see the effects of midway equalization on the quantization
noise. Let $u_0$ and $u_1$ be two images, and $\tilde{u}_0$ and
$\tilde{u}_1$ the same images after quantization on the integers
of $[0,255]$. We will assume as a good approximation that the quantization noise $b_i=u_i-\tilde{u}_i$ is uniformly distributed on $[-0.5,0.5]$. Thus, $E(b_i)=0$ and its
standard deviation is $\sqrt{\int_{-0.5}^{0.5}x^2
dx}=\sqrt{\frac{1}{12}} \simeq 0.289$. The following table shows
what becomes the standard deviation of $b_0$ and $b_1$ after the
different transformations \footnote{The non-quantized images used to compute these results have been obtained in the following way: for the multi-channel images and the stereo pair, we used the original images, which were coded in float values; for the painting, the images were already quantized, so we added a simulated noise, uniformly distributed on $[-0.5,0.5]$. }. For example, if $u_0$ after
specification on $u_1$ becomes $u_{0,1}$, and $\tilde{u}_0$ after
specification on $\tilde{u}_1$ becomes $\tilde{u}_{0,1}$,
we consider that the quantization noise $b_0$ is transformed in
$\tilde{b}_0=u_{0,1}-\tilde{u}_{0,1}$.



\begin{center}
\begin{tabular}[]{|l|c|c|c|c|r|}
\hline
$$ &  Specification & I.a.e. & S.i.r. \\
\hline
stereo left & $1.207121$ & $0.646560$ & $0.442230$ \\
\hline
stereo right &  $0.611246$ & $0.479748$ & $0.538039$\\
\hline
blue channel &  $1.848390$ & $1.069263$ & $1.040278$  \\
\hline
infrared channel  & $0.145176$ & $0.961488$ & $0.480899$\\
\hline
meduse left  & $0.598254$ & $0.458458$ & $0.540970$ \\
\hline
meduse right  & $0.393870$ & $0.354508$ &  $0.462411$\\
\hline
\end{tabular}
\end{center}

We observe that the effects of the I.a.e. and of the S.i.r. on the
noise are far more symmetrical than those of specification. With
these methods, the standard deviation of the noise is
multiplied by a number smaller than $2.5$, except in the case of
the multi-channel images. Indeed, those images are farther from
each other and consequently   subject to more drastic
transformations. In this case, the ratio is smaller than $3.7$.


\subsection{The Shannon point of view: effects on the spectrum}

In general, the composition of an image by a function does not maintain
the property ``to be well sampled, according to Shannon-Whittaker
sampling theorem''. Now, we argue that when we use correlation
technics on a stereo pair, in order to get the depth map of the
area, the images must be well sampled if we want a good precision
(see chapter~\ref{chapter_analytic}).

In theory, we understand that a modification of the histogram via
the s.i.r. algorithm must be smoother for the spectrum than direct
specification. Indeed, if $\psi:[0,1] \rightarrow [0,1]$ can be
approximated by a polynomial function:
$\psi(s)=s+a_2s^2+a_3s^3+...$, then:
$$\widehat{\psi(u)}=\hat{u}+a_2\hat{u}*\hat{u}+a_3\hat{u}*\hat{u}*\hat{u}+...$$
The values $a_i$ must be as small as possible in order for
$\hat{\psi(u)}$ to remain well sampled. This happens if $\psi$ remains as close to identity as possible and is smooth enough. Now, we have seen that the
functions used in the s.i.r. (see for example
Figure~\ref{histo_et_rac_1112.ps}) and the i.a.e. methods are indeed
closer to the identity than the one used in specification
(${H_0}^{-1} \circ H_1$). 
%However, it is not enough to conclude.

\begin{figure}[h]
\centerline{\includegraphics[width=3.0cm]{spectre1}
\includegraphics[width=3.0cm]{spectre2}
\includegraphics[width=3.0cm]{spectre3}
\includegraphics[width=3.0cm]{spectre4}
}
\caption{\em Spectra of $u_0$, $\bar{u}_0$, $\tilde{u}_0$ and $\breve{u}_0$. We can see that the spectra of $\tilde{u}_0$ and $\breve{u}_0$ remain more
condensed in the central square than the one of $\bar{u}_0$.}
\label{spectres.ps}
\end{figure}

In order to check experimentaly the effects of the different transformations on the spectrum, we construct the
following example: the left half of each image of the couple
at the top of Figure~\ref{toulouse} is zoomed by FFT, which gives $u_0$ and
$u_1$. Next, we specify $u_0$ on $u_1$ (result $\bar{u}_0$), and
we compute the mo\-dification of $u_0$ via the s.i.r. algorithm
(result $\tilde{u}_0$), or the i.a.e. algorithm (result $\breve{u}_0$).
Finally, we compute the spectra of $u_0$, $\bar{u}_0$,
$\tilde{u}_0$ and $\breve{u}_0$ (Figure~\ref{spectres.ps}). We see
that the spectra of $\tilde{u}_0$ and $\breve{u}_0$ remain more
in the central square than the spectrum of $\bar{u}_0$. Actually, if we
compute the rate of energy which goes out of the central square
after different modifications, we get $0.0107 $ for specification,
$0.0027$ for s.i.r., and $0.0013$ for i.a.e. The energy stays more
condensed in the spectral support of $u_0$ after an s.i.r. or an
i.a.e. modification. Consequently, $\tilde{u}_0$ and $\breve{u}_0$
remain better sampled than $\bar{u}_0$ \footnote{ Most natural
images have a spectrum really concentrated in the center of its
support. It follows that after a transform like a histogram
modification, the spectrum sprawls but does not go much out of its
support. The ratios of energy presented above illustrate this
fact. However, we win a factor larger than $5$ and this gain can
be really useful in satellite imaging applications to depth map
recovery.}.

\subsection{Comparison using the Kullback-Leibler distance}

Finally, we choose to compare experimentaly the results of both midway equalizations for a statistical distance. The following table contains the Kullback-Leibler distances \footnote{The distance between two histograms was computed symmetrically as $\frac{1}{2}(\int h_0\log \frac{h_0}{h_1}+\int h_1\log \frac{h_1}{h_0})$. The histograms were previously lightly smoothed in order to avoid them to vanish at some point.} between the histograms of the previous pairs, before, and after each kind of transformation.

\begin{center}
\begin{tabular}[]{|l|c|c|c|c|r|}
\hline
$$ &  before any transformation &  after I.a.e. & after S.i.r. \\
\hline
stereo pair & $1.622677$  & $0.005506$ & $0.001155$ \\
\hline
Different channels &  $1.586096$  & $0.008154$ &  $0.017710$  \\
\hline
\textit{Radeau de la Meduse}& $0.298163$  & $0.001402$ & $0.001435$ \\
\hline
\end{tabular}
\end{center}


We can observe that the algorithm giving the best performances depends on the input pair. So, these results do not allow us to conclude experimentally on the choice of the method.  


\section{Conclusion}
This chapter presented two ways of giving to a pair of images the
same histogram, while remaining coherent with their previous grey
level dynamics.  This intermediate histogram has been constructed
following some elementary axioms, which have
led to one method or the other. We have seen that they can be
applied to a wide range of image pairs, giving in each case a
better representation of the images for observation, comparison or
ulterior manipulations. It has been shown that those
transformations are smoother for the images than direct
specification: they reduce the loss of visual information
generated by this kind of alteration and are better towards the
quantization noise. They are also more respectful of the property
``to be well sampled''. The theoretical and experimental comparison of these midway equalizations was not enough to conclude on the choice of the best method. Both methods being acceptable, only practice allows to decide. From the structural point of view, the result given by i.a.e. is simpler and unique, given by an explicit formula. The s.i.r. equalization has to be considered as an alternative, and gives experimentally similar results, but its construction is certainly needlessly complicated, and more sensible to robustness problems. Consequently, we consider that the i.a.e. should be recommended as a standard method, which can be applied efficiently by the algorithm proposed by Cox et al.

We shall now study an extension of the midway equalization to image sequences. A generalization to color images will be proposed in chapter~\ref{chapter_transport}.






