% Template for ICIP-2017 paper; to be used with:
%          spconf.sty  - ICASSP/ICIP LaTeX style file, and
%          IEEEbib.bst - IEEE bibliography style file.
% --------------------------------------------------------------------------
\documentclass{article}
\usepackage{spconf,amsmath,amssymb,eucal,graphicx}
\usepackage{cite}


% Local definitions.
% --------------------
\newcommand{\blf}{{\bf f}}
\newcommand{\bx}{{\mathbf x}}
\newcommand{\bk}{{\mathbf k}}
\newcommand{\be}{{\mathbf e}}
\newcommand{\bF}{{\mathbf F}}
\newcommand{\bH}{{\mathbf H}}
\newcommand{\bn}{{\mathbf n}}
\newcommand{\bZero}{{\mathbf 0}}
\newcommand{\bv}{{\mathbf v}}
\newcommand{\bw}{\mbox{\boldmath $\omega$}}
\newcommand{\bz}{{\mathbf z}}
\newcommand{\Hi}{{\cal H}}
\newcommand{\del}[1]{\nabla{#1}}
\newcommand{\Eqref}[1]{{(\ref{#1})}}
\renewcommand{\Im}[1]{\hbox{Im}\left[{#1}\right]}
\newcommand{\lra}{{\longrightarrow}}
\newcommand{\R}{{\mathbb R}}
\newcommand{\N}{{\mathbb N}}
\newcommand{\ra}{{\rightarrow}}
\newcommand{\Z}{{\mathbb Z}}
\newcommand{\C}{{\mathbb C}}
\renewcommand{\Re}[1]{\hbox{Re}\left[{#1}\right]}
\newcommand{\sgn}{\mathop{\rm sgn}\nolimits}
\newcommand{\vphi}{\varphi}
\newcommand{\w}{\omega}


% Title.
% ------
\title{AM-FM Image Filters}
%
% Single address.
% ---------------
\name{Chuong T. Nguyen and Joseph P. Havlicek
\thanks{This work was supported in part by the U.S.~Army Research Laboratory
and the U.S.~Army Research Office under grant W911NF-04-1-0221.}}
\address{School of Electrical and Computer Engineering,
   University of Oklahoma, Norman, OK, USA }

%
% For example:
% ------------
%\address{School\\
%	Department\\
%	Address}
%
% Two addresses (uncomment and modify for two-address case).
% ----------------------------------------------------------
%\twoauthors
%  {A. Author-one, B. Author-two\sthanks{Thanks to XYZ agency for funding.}}
%	{School A-B\\
%	Department A-B\\
%	Address A-B}
%  {C. Author-three, D. Author-four\sthanks{The fourth author performed the work
%	while at ...}}
%	{School C-D\\
%	Department C-D\\
%	Address C-D}
%
\begin{document}
%
% The following line forces 9-point fonts to gain some space.
% Comment it out for the default 10-point font size.
\ninept
%
\maketitle
%
\begin{abstract}
We introduce a multicomponent invertible AM-FM image transform and use
it to define new nonlinear AM-FM filters for performing modulation domain image
processing.  The key elements of the transform are analysis and synthesis
filterbanks based on the steerable image pyramid and perfect reconstruction
demodulation algorithms based on analytic differentiation of continuous
cubic tensor spline models fit to the unwrapped phase samples of a digital image.
We demonstrate spatially and spectrally localized orientation and frequency
selective filtering, simple image restoration, and image fusion in the
modulation domain.  These results are also among the first to demonstrate high
fidelity image reconstructions from computed multicomponent AM-FM models.
\end{abstract}
%
\begin{keywords}
AM-FM image models, AM-FM image filters, modulation domain signal processing,
multicomponent models
\end{keywords}
%
\section{Introduction}
\label{sec:intro}

Multicomponent AM-FM image models~\cite{MariosPAMI07,JOSA-Sep95,IP-Multi,IT-92}
represent an image $t\colon\R^2\ra\R$ as a sum
\begin{equation}
\label{eq:RealModel}
t(\bx) = \sum_{k=1}^K t_k(\bx) = \sum_{k=1}^K a_k(\bx)\cos[\vphi_k(\bx)]
\end{equation}
of nonstationary AM-FM functions
\begin{equation}
\label{eq:t_kReal}
t_k(\bx) = a_k(\bx)\cos[\vphi_k(\bx)],
\end{equation}
where $\bx\in\R^2$, $a_k(\bx) \geq 0$, and $\vphi_k(\bx)\colon\R^2\ra\R$.
Given the image $t(\bx)$, a computed AM-FM model consists of estimates
of the $K$ AM functions $a_k(\bx)$ which provide a dense local
characterization of the local texture contrast and the $K$ FM functions
$\del{\vphi_k(\bx)}$ which provide a dense characterization of the local
texture orientation and pattern spacing.

Such models have been used with great success
in a variety of image {\em analysis\/}
applications, including texture
segmentation,
3-D shape from texture,
texture-based stereopsis,
fingerprint classification,
content-based retrieval,
and regeneration of occluded and damaged textures~\cite{AP982ED},
as well as for infrared target tracking~\cite{ChuongOTCBVS07} and
in the analysis of (2-D) spectrograms of
human speech signals~\cite{EzzatPoggioICASSP07}.
To date, however, they have been considerably less successful in
applications requiring image {\em synthesis\/} in addition to 
analysis (to the best of our knowledge, reconstruction from
a computed AM-FM model has been attempted previously
only in~\cite{DoerschukSP96,IP-Jun96,IP-Multi}).
The reason is that some means must be devised for decomposing
the image into a sum of components~\Eqref{eq:t_kReal}
that are {\em locally coherent}~\cite{IT-92,IP-Multi} and
for isolating these components from one another on a jointly
localized basis in space and spatial frequency prior to
demodulation.  Because of their excellent joint localization
properties, 
banks of Gabor or Gabor-like bandpass filters have been used
almost universally for this purpose.
Unfortunately, the very properties that make these filters
attractive, {\em viz.}, localization and smoothness of the 
impulse and frequency responses simultaneously, also
imply that they cannot provide perfect 
reconstruction~\cite{AP982ED}.






\begin{figure}[t]
\centering
% CCAMD 
\includegraphics[width=3.0in]{lena}
\caption{Perfect reconstruction AM-FM model for performing
image processing in the modulation domain.}
\label{fig:CCAMD}
\end{figure}

Our goal is to formulate a new perfect
reconstruction (or approximately perfect reconstruction)
AM-FM image model and use it to perform signal processing
in the modulation domain, as depicted schematically in
Fig.~\ref{fig:CCAMD}.  Here, the AM-FM model plays a role
analogous to that of the Fourier transform in traditional
LTI filter theory; in particular, it is a clear requirement
that the input and output images be visually indistinguishable
when the signal processing block in Fig.~\ref{fig:CCAMD}
is set to the ``do nothing" operation.

Previously, we achieved perfect reconstruction
in~\cite{SivleyICIP06} by replacing the typical Gabor
filterbank with a nonseparable, orientation selective
perfect reconstruction wavelet filterbank.  But there were
two inherent difficulties with using this approach to
formulate practical AM-FM image filters.  First, a very
large number of filterbank channels was required to
generate locally coherent image components $t_k$, on the
order of $K=500$ for a $256 \times 256$ image.  Second,
due to the aliasing that invariably occurs in a critically
sampled wavelet filterbank~\cite{SimoncelliIT92}, it proved
difficult to define signal processing operations
on the modulating functions $a_k$ and $\del{\vphi_k}$ that would
lead to predictable, visually meaningful changes in the output
image relative to the input.

In this paper, we introduce a new invertible AM-FM image transform
where the analysis and synthesis filterbanks are based on an
appropriate adaptation of the steerable image pyramid
decomposition~\cite{SimoncelliIT92,SimoncelliICIP95}.  The
steerable pyramid admits several desirable properties that
make it attractive for computing AM-FM image models.  The bandpass
channels are Gabor-like and orientation selective, leading to
locally coherent image components that correspond well to visual
perception of the significant image structure.  Moreover,
because the steerable
pyramid is essentially alias free and self-inverting, it
naturally facilitates the intuitive design of AM-FM filters
capable of achieving visually meaningful, perceptually
motivated signal processing goals.  
We briefly describe computation of the invertible AM-FM transform
in Section~\ref{sec:AMFMTransform}, while practical examples of
AM-FM image filters are given in Section~\ref{sec:Examples}.





%\section{Partial Rotational Hilbert Transform}
%\label{sec:RotationalHilbert}
%Let $f(t)$ be a real signal, the 1-D analytic signal $z(t)$ constructed as
%\begin{equation}
%z(t) = f(t) + j\Hi\left[f(t)\right].
%\end{equation}
%
%$g(t)$ .The construction of the 1-D analytic signal $z(t)$ is well-defined. 
%The real part of $z(t)$ comes from the original real signal $f(t)$, while
%the imaginary part of $z(t)$ is computed by taking the Hilbert transform 
%of $f(t)$. When extended to n-D signal, the definition of the Hilbert
%transform is no longer unique. There exist different approaches proposed 
%in the literature such as hypercomplex signal by Bulow and 
%Sommer~\cite{BulowSP01}, signal orthant spectra by Haln~\cite{Haln92}, and
%partial Hilbert transform by Havlicek et. al.~\cite{ICIP97}.
%
%We restrict to the case of 2-D image signal. The imaginary part of the 
%complex signal $z(t)$ is constructed by the partial Hilbert transform
%approach. This construction, however, results in jumps between left half
%and right half of the 2-D frequency plane. Therefore, it is undesirable
%to have these jumps to occur within the passband of directional filters.
%To alliviate this problem, we rotate the direction of the Hilbert 
%transform such that it is perpendicular to the direction of directional
%filters. As a result, the spectrum of the complex image $z(t)$ is defined as
%\begin{equation}
%Z(\bw) = \begin{cases}
%2F(\bw) & \bw^T.\bn > 0 \\
%F(\bw) & \bw^T.\bn = 0 \\
%0   & otherwise
%\end{cases},
%\label{eq:RotationHilbert}
%\end{equation}
%where $\bn$ is the unit vector parallel to the direction of 
%directional filters.


\section{Invertible AM-FM Transform}
\label{sec:AMFMTransform}

A number of good approaches have been proposed for computing 
continuous and discrete domain multicomponent AM-FM image
models, including approximate techniques based on the multidimensional 
Teager-Kaiser energy operator~\cite{JOSA-Sep95} and
on a multidimensional extension of the 1-D analytic signal~\cite{IP-Multi}, as
well as emerging new techniques based on Max-Gabor analysis
over local patches~\cite{EzzatPoggioICASSP07} and on
the empirical mode decomposition~\cite{HuangEMD98}.  Here we focus on the
analytic signal-based approach because, for a complex extension
$z_k(\bx) = t_k(\bx) + j\Hi[t_k(\bx)] = a_k(\bx)\exp[j\vphi_k(\bx)]$ of a
continuous domain AM-FM
image component $t_k$ given by~\Eqref{eq:t_kReal}, where $\Hi[\cdot]$ is the
partial Hilbert transform, it leads to an
{\em exact\/} demodulation algorithm~\cite{IP-Multi}
\begin{eqnarray}
\del{\vphi_k(\bx)} &=&
    \Re{\frac{\del{z_k(\bx)}}{j z_k(\bx)}} \label{eq:FMalg}\\
a_k(\bx) &=& \left|z_k(\bx)\right|\label{eq:AMalg}
\end{eqnarray}
that is free of approximation errors.  Moreover, provided that boundary
conditions on the phase of each component are saved as part of the AM-FM
model, the component $t_k$ can be recovered exactly from the
modulations~\Eqref{eq:FMalg} and~\Eqref{eq:AMalg} by integrating to recover
the phase and then reconstructing according to~\Eqref{eq:t_kReal}.

To isolate the individual components $t_k(\bx)$, we implement the steerable
pyramid as described in~\cite{SimoncelliIT92,SimoncelliICIP95} to a depth of
four levels with a non-orientation selective high-pass filter $H_0$ and
non-orientation selective low-pass filters $L_n$ for $0 \leq n \leq 3$.  For each
of $L_0$, $L_1$, and $L_2$, the response is partitioned into eight orientation
selective channels using  oriented bandpass filters $B_p$, $0 \leq p \leq 7$.  
Although this
construction completes the usual implementation of the steerable pyramid, it is
not optimal for AM-FM modeling because the high-pass and low-pass residuals
delivered by the channels $H_0$ and $L_3$ are not orientation selective and
hence not locally coherent.  Therefore, we adapt the steerable pyramid by
further decomposing the responses of $H_0$ and $L_3$ into eight
orientation selective channels using oriented filters
$[-j\cos(\theta - p\pi/8)]^7$, $0 \leq p \leq 7$, obtained by discarding
the magnitude characteristic $B(\bw)$ from the right-hand side of
equation (5) in~\cite{SimoncelliICIP95}.  The residual remaining
in channels $H_0$ and $L_3$ after this decomposition tends to be small
in magnitude and generally contributes little to visual perception of the image;
we therefore discard it.  This results in an adapted steerable pyramid
decomposition with 40 channels delivering $K=40$ AM-FM image
components $t_k$.

We apply AM-FM demodulation between the analysis and synthesis filterbanks.  
For each component, the discrete versions of~\Eqref{eq:FMalg} 
and~\Eqref{eq:AMalg} given
in~\cite{IP-Multi} are first applied to obtain the AM functions
$a_k$ and estimates of the FM functions 
$\del{\vphi_k}$.   The estimated FM functions are then used to
guide a multidimensional phase unwrapping algorithm as described
in~\cite{SivleyICIP05}.  The unwrapped phase is fit with cubic
tensor product splines that we differentiate analytically
to obtain the frequency field $\del{\vphi_k}$~\cite{SivleyICIP06,SivleySSIAI06}.
















\section{Examples}
\label{sec:Examples}

To define AM-FM image filters, we apply signal processing to the
computed modulating functions expressed in polar coordinates
according to 
$r_k = |\del{\vphi_k}|$ and $\psi_k = \arg\del{\vphi_k}$.
The filtered image is then obtained by first reconstructing the
filtered image components from the processed modulations
according to~\Eqref{eq:t_kReal} and then subjecting the
filtered components to the synthesis filterbank.
By effectively realizing the modulation domain signal processing
paradigm depicted in Fig.~\ref{fig:CCAMD}, this approach opens
the door to an entirely new class of nonlinear image filters that
we have only begun to investigate.  We present a few examples
in this section, where we denote the processed modulations by
$\widehat a_k(\bx)$ and $\del{\widehat\vphi_k(\bx)}$.


%
\begin{figure*}[tp]
\centering
% filterbank
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(a)}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(b)}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(c)} 
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(d)} \\
\vspace{5pt}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(e)}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(f)}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(g)} 
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(h)} \\
\vspace{5pt}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(i)}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(j)} 
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(k)}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(l)} \\
\vspace{5pt}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(m)}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(n)}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(o)}
\includegraphics[width=1.53in,clip=true]
  {lena}\makebox[1.25em]{(p)}
\caption{
Examples.
(a) radial chirp image. 
(b) LTI notch filter frequency response.
(c) Output of LTI notch filter.
(d) Difference image between (a) and (c) showing lack of spatial localization
    in the LTI filtering operation.
(e) Orientation selective attenuation performed by AM-FM notch filter.
(f) Frequency selective enhancement performed by AM-FM bandpass filter.
(g) Original Lena image.
(h) Corrupted by linear blur and additive noise.
(i) Result of na\"{i}ve LTI high-pass restoration filter.
(j) Result of elementary AM-FM enhancement/restoration filter.
(k),(l) Clocks input image pair.
(m) AM-FM image fusion result.
(n),(o) Registered CT and MR images ({\tt imagefusion.org}).
(p) AM-FM image fusion result.
}
\label{fig:Examples}
\end{figure*}

A synthetic radial chirp image is shown in Fig.~\ref{fig:Examples}(a).  By
construction, the amplitude is constant and the phase is quadratic along
radials emanating from the center of the image.  Our signal processing goal
is to attenuate nonstationary structure that is oriented at odd multiples
of $\pi/4$.  As a baseline for comparison, we implemented an LTI notch
filter with frequency response given in Fig.~\ref{fig:Examples}(b).  The
output of this filter is shown in Fig.~\ref{fig:Examples}(c) and
exhibits undesirable artifacts as expected:  the LTI filter attenuates
Fourier components on a spatially global scale, which achieves the
desired result but also degrades the subtle constructive and destructive
interference between Fourier components that creates the image structure
at orientations other than odd multiples of $\pi/4$.  This is demonstrated
by Fig.~\ref{fig:Examples}(d), which gives the residual between the
original image in Fig.~\ref{fig:Examples}(a) and the LTI processed image
in Fig.~\ref{fig:Examples}(c).  We define an AM-FM notch filter by 
\begin{equation}
\label{eq:chirpFilterAM}
\widehat a_k(\bx) = \left\{\begin{array}{ll}
16\delta_k(\bx)a_k(\bx) / \pi, &
   \mbox{$\delta_k(\bx) < \pi/16$,}\\
a_k(\bx), & \mbox{otherwise}
\end{array}\right.
\end{equation}
and $\del{\widehat\vphi_k(\bx)} = \del{\vphi_k(\bx)}$, where
$\delta_k(\bx) = |(|\psi_k(\bx)| - \pi/4)|$ is an amplitude
scaling factor equal to the radian angular distance between
$\psi_k(\bx)$ and $\pm\pi/4$ (because the partial Hilbert transform
has action in the horizontal direction, $\psi_k(\bx)$ is restricted
to quadrants I and IV of the frequency plane).
The resulting processed image is given in Fig.~\ref{fig:Examples}(e),
where it may be seen that the signal processing goal has been achieved.
Moreover, because the AM-FM filter is capable of attenuating oriented
structure on a spatially local basis, it delivers a result that is
free of the undesirable artifacts seen in the LTI filter output.


In addition to orientation selective processing,
AM-FM filters can also be used to perform spatially local
amplification and attenuation based on magnitude frequency.  
To amplify a band of radial frequencies, we
implemented a frequency selective filter according to
\begin{equation}
\label{eq:chirpFilterRing}
\widehat a_k(\bx) = \left\{\begin{array}{ll}
2a_k(\bx), & 0.2 < \mbox{$r_k(\bx) < 0.35$,}\\
a_k(\bx), & \mbox{otherwise}
\end{array}\right.
\end{equation}
and $\del{\widehat\vphi_k(\bx)} = \del{\vphi_k(\bx)}$, where
$r_k(\bx)$ is given in units of cycles per pixel.  The AM-FM
filtering result is given in Fig.~\ref{fig:Examples}(f).

Simplistic AM-FM image restoration is illustrated in 
Fig.~\ref{fig:Examples}(g)-(j).  
The original image is
given in Fig.~\ref{fig:Examples}(g), while the image in
Fig.~\ref{fig:Examples}(h) was obtained by applying a
low-pass linear blur and adding Gaussian white noise.
As a baseline comparison, the result of a na\"{i}ve high-pass 
LTI filter approximating the pseudo-inverse is shown in
Fig.~\ref{fig:Examples}(i).  The design concept for the AM-FM 
restoration filter combines elements similar to both
wavelet shrinkage and unsharp masking.  The noise power is
distributed widely throughout the steerable pyramid channels resulting
in a relatively small contribution to the individual 
AM functions $a_k$, whereas the coherent image structure tends
to be jointly localized resulting in strong contributions
to the AM functions, particularly in the vicinity of edges.
We applied a simple threshold to the amplitude modulations
computed from Fig.~\ref{fig:Examples}(h) to attenuate the noise,
reconstructed, and then applied the same high-pass filter that
was used in Fig.~\ref{fig:Examples}(i) to generate a high-pass
mask.  This mask was added back to the degraded image in
Fig.~\ref{fig:Examples}(h) to obtain the enhanced/restored
result shown in Fig.~\ref{fig:Examples}(j).

Finally, elementary examples of AM-FM image fusion based on 
local contrast are given Fig.~\ref{fig:Examples}(k)-(m) and
Fig.~\ref{fig:Examples}(n)-(p).  The main idea is that
sharp intensity edges are associated with local
AM values $a_k(\bx)$ that are relatively large, whereas
defocused surfaces are associated with AM values that are
relatively smaller than those of their in-focus counterparts.  
Thus,
for a problem like the well-known pair of clock images shown in
Fig.~\ref{fig:Examples}(k) and~(l), a fused image showing
both clocks in focus 
can be obtained as follows.  We first compute a multicomponent
AM-FM model for each of the two input images shown in
Fig.~\ref{fig:Examples}(k) and~(l).  For each pixel 
$\widehat t_k(\bx)$ of component $k$ in the fused image, we 
take $\widehat a_k(\bx)$ and
$\del{\widehat\vphi_k(\bx)}$ directly from the input image for
which $a_k(\bx)$ is larger on a pixel-by-pixel basis.
The fused image result is given in Fig.~\ref{fig:Examples}(m).
An identical AM-FM algorithm was used to obtain the result
shown in Fig.~\ref{fig:Examples}(p) by fusing the CT image
of Fig.~\ref{fig:Examples}(n) and the MR image
of Fig.~\ref{fig:Examples}(o) (CT and MR images
courtesy of {\tt imagefusion.org}).




\section{Conclusion}
\label{sec:Conclusion}
In this paper, we introduced an invertible AM-FM image transform based on
an adapted steerable image pyramid in combination with powerful
perfect reconstruction joint demodulation algorithms.  We used this transform
to develop new nonlinear AM-FM filters for performing image processing
in the modulation domain.  Compared to LTI filtering, the advantages of
AM-FM filters are that they are capable of performing spatially and
spectrally localized processing directly in terms of the visually
important nonstationary image structure.  The results given here
represent an entirely new class of nonlinear image filters that
we have only begun to investigate.  We demonstrated several simple
AM-FM filters, including notch and bandpass filters for amplifying
and attenuating spatially local structure based on orientation and
granularity ({\em i.e.}, magnitude frequency) and an elementary AM-FM image 
restoration filter.
We also applied AM-FM filtering to perform image fusion based on
local contrast in two typical scenarios.  Compared to previous AM-FM
image analysis techniques, these results are significant in that they
are among the first to demonstrate high fidelity AM-FM image reconstructions
and intuitive, practical modulation domain filter design.

Most of the filters considered here relied heavily on AM-based processing.
FM-based processing is generally more difficult due to the line singularity
associated with the partial Hilbert transform, which tends to perturb the
computed frequencies moreso than the amplitudes.
We are currently developing
strategies for rotating the Hilbert transform on a component-by-component
basis to steer the singularity away from the instantaneous frequency
vectors $\del{\vphi_k}$.  It will also be interesting to investigate a
wide variety of AM-FM filtering techniques, including the application of
LTI and nonlinear spatial filters directly to the computed modulations.








% References should be produced using the bibtex program from suitable
% BiBTeX files (here: strings, refs, manuals). The IEEEbib.bst bibliography
% style file from IEEE produces unsorted bibliography list.
% -------------------------------------------------------------------------
\bibliographystyle{IEEEbib}
\bibliography{refs}

\end{document}
