 % Template for ICASSP-2004 paper; to be used with:
%          spconf.sty  - ICASSP/ICIP LaTeX style file, and
%          IEEEbib.bst - IEEE bibliography style file.
% --------------------------------------------------------------------------
\documentclass{article}
\usepackage{spconf,amsmath,epsfig,amssymb,mathrsfs}
%\usepackage{psfrag}
%\psfrag{+x}{+ \sf j}
%\psfrag{-x}{- \sf j}

% Example definitions.
% --------------------
\def\x{{\mathbf x}}
\def\L{{\cal L}}

% Title.
% ------
\title{\Large DIRECTIONAL HYPERCOMPLEX WAVELETS \\[0mm]
FOR MULTIDIMENSIONAL SIGNAL ANALYSIS AND PROCESSING}
%
% Single address.
% ---------------
\name{Wai Lam Chan, Hyeokho Choi, Richard G.\ Baraniuk
\thanks{Email: \{wailam, choi, richb\}@rice.edu. Web:
dsp.rice.edu. This work was supported by ONR grant
N00014--02--1--0353, AFOSR grant F49620-01-1-0378, DARPA, and the
Texas Instruments Leadership University Program.}}
\address{%\small
Department of Electrical and Computer Engineering, Rice
University, Houston, Texas, USA}

%\vspace{-0.5in}
%
% 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}
\ninept
\maketitle
%
\begin{abstract}
\sloppy We extend the wavelet transform to handle multidimensional
signals that are smooth save for singularities along
lower-dimensional manifolds. We first generalize the complex
wavelet transform to higher dimensions using a multidimensional
Hilbert transform. Then, using the resulting {\em hypercomplex
wavelet transform} (HWT) as a building block, we construct new
classes of nearly shift-invariant wavelet frames that are oriented
along lower-dimensional subspaces. The HWT can be computed
conveniently using a 1-D dual-tree complex wavelet transform along
each signal axis. A series of denoising experiments demonstrates
that the HWT significantly outperforms both standard tensor
product wavelets and complex wavelets when the signals of interest
sport a local manifold structure.
\end{abstract}
%
\vspace{-0.15in}

\section{INTRODUCTION}
\label{sec:intro}

\vspace{-0.05in}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%key concept:
%     1-d               2-d                   3-d
%dwt   SV          SV, non-directional      SV, non-directional
%cwt   SI          SI, directional, but     SI, but 2-dim manifold dim
%                  1 manifold dim
%
%goal: higher dim, shift invariant, directional wavelets with different
%      manifold dimension
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\sloppy The multiscale analysis of the discrete wavelet transform
(DWT) has proved natural and powerful for dealing with
singularity-rich one-dimensional (1-D) signals. In contrast to the
sinusoidal basis functions underlying Fourier analysis, wavelet
basis functions are localized in both time and frequency and thus
yield very sparse and structured representations of {\em piecewise
smooth signals} (signals that are smooth except for a finite
number of jump discontinuities). This sparsity and structure boost
the performance of algorithms such as denoising by shrinkage
\cite{donoho95denoising} and compression
\cite{Vetterli95wavelets}.


In this paper, we study how to best extend wavelet representations
to handle piecewise smooth {\em multidimensional signals} (2-D,
3-D, and beyond).  Such signals are smooth almost everywhere over
their $n$-D domain but feature singularities on {\em manifold}
structures of lower dimension $m\leq n$.\footnote{A
manifold is a topological space that is locally Euclidean
\cite{manifold}.} For example, an edge discontinuity in a 2-D
image lies along a 1-D curve; a bee or aircraft flying through 3-D
space tracks out a smooth 1-D curve in 4-D space-time; a video of
a smooth object moving smoothly in time lies on a 2-D manifold in
3-D space-time; and points in space can be viewed as lying on 0-D
manifolds in 3-D space.


Two fundamental problems have hindered the application of wavelets
to multidimensional piecewise smooth signals. First, manifold
singularities have a preferred {\em orientation}: they are smooth
in the ``direction'' of the manifold (along the $m$-D tangent
space) and contain a singularity (rapid change) in the
``direction'' of the normal (along the $(n-m)$-D normal space).
Thus, to obtain a sparse representation, wavelets intersecting the
manifold should be smooth (lowpass) in the direction of the
manifold and oscillatory (bandpass) in the direction of the
normal. Unfortunately, standard separable wavelets constructed
from tensor products of 1-D wavelets have poor orientation
selectivity. In 2-D, for example, one of the three tensor
wavelets is simultaneously oscillatory at {\em both} $+45^\circ$
and $-45^\circ$ in the image plane. Figure \ref{fig:basis}(a)
illustrates that the situation is even more grim in 3-D; this tensor wavelet is
simultaneously oscillatory along the 4 diagonal directions within the cube. In
general, tensor wavelets are only appropriate for representing
$0$-D structures (points) in $n$-D space.


Second, unlike the complex Fourier representation, the standard
real-valued wavelet transform in any dimension is highly {\em
shift varying}, meaning that a small shift in the signal
drastically changes the values of the wavelet coefficients.  This
complicates wavelet coefficient modeling and processing
considerably.


\sloppy Recently, a promising approach for generating directional
wavelets has emerged via the {\em complex wavelet transform} (CWT)
\cite{lina,kingsbury,ivan}.  Kingsbury's 1-D {\em dual-tree} CWT
expansion employs two independent wavelet filterbanks
corresponding to a ``real'' wavelet and an ``imaginary'' wavelet
that form an approximate Hilbert transform pair \cite{kingsbury}.
The resulting $2\times$ redundant tight frame enables a
magnitude/phase coherent representation akin to that of Fourier
analysis.  Moreover, in 2-D a clever combination of the
independent 1-D CWT filterbanks creates a directional set of six
complex wavelets oriented at multiples of $30^\circ$; these
wavelets are very natural for representing 1-D edge structure in
images \cite{wakin}. As an added bonus, the smooth magnitude
envelope of the complex wavelet yields a nearly shift-invariant
transform.


In \cite{ivan}, Selesnick and Li extended this 2-D construction to a
3-D CWT that is natural for studying video sequences and other data with 
($n=3$, $m=2$). As we see in Fig.\ \ref{fig:basis}(b), each 3-D CWT
wavelet oscillates along a 1-D line (its orientation) and is
smooth in the orthogonal 2-D plane. Extended further to $n$-D,
Selesnick and Li's CWT is natural for representing signals always
lying on $m=n-1$ dimensional manifolds.


\begin{figure*}[htb]
\begin{center}
\begin{tabular}{ccccc}
(a) \epsfig{figure=DWT3Dbasis.eps,width=2.5cm} &\hspace*{5mm}&
%
(b) \epsfig{figure=CWT3Dbasis.eps,width=2.5cm} &\hspace*{5mm}&
%
(c) \epsfig{figure=QWT3Dbasis.eps,width=2.5cm}
\end{tabular}
\end{center}

\vspace{-0.2in}

\caption{\sl Isosurfaces (points taking the same function value)
of wavelet bases in 3-d space ($n=3)$.  Red is large positive;
blue is large negative. (a) Standard DWT tensor-product wavelet
lacks any strong directionality; (b) $m=2$ DHWT wavelet (real
part) is bandpass along a 1-D line and lowpass along the
orthogonal 2-D plane \protect\cite{ivan}; (c) $m=1$ DHWT wavelet
(real part) is bandpass on a 2-D plane and lowpass on the
orthogonal 1-D line.}\label{fig:basis}
%
\end{figure*}



\begin{figure*}[htb]
\begin{center}
\begin{tabular}{ccccc}
(a) \epsfig{figure=DWT3D.eps,width=2.5cm} &\hspace*{5mm}&
%
(b) \epsfig{figure=CWT3D.eps,width=2.5cm} &\hspace*{5mm}&
%
(c)~~ \epsfig{figure=QWT3D.eps,width=2.5cm}
\end{tabular}

\end{center}

\vspace{-0.2in}

\caption{\sl Fourier-domain support of the 3-D wavelets from Fig.\
\protect\ref{fig:basis} ($n=3$). Part (a) also represents the
Fourier support of the corresponding HWT wavelet.}
\label{fig:cube}
%
\vspace{-0.2in}
\end{figure*}


But what about signals with singularities on $m<n-1$ dimensional manifolds
(the 1-D trajectory of a bee in 3-D, for example)?  In this paper,
we develop a general theory for directional wavelets that are
matched to any signal dimension $n$ and manifold dimension $m\leq
n$. We take a two-step approach.  First, we introduce a new {\em
hypercomplex wavelet transform} (HWT) that generalizes the CWT via
an analogy to the hypercomplex Fourier transform and its
associated generalized $n$-D Hilbert transform \cite{bulow_HFT}.
The HWT is conveniently obtained by a simple non-separable
combination of 1-D dual-tree CWTs. Second, we show that the HWT
serves as a basic building block for generating new classes of
{\em directional hypercomplex wavelet transforms} (DHWTs) for
handling different manifold dimensions $m$; we merely form a
simple linear combinations of the HWT components(orthogonal
projections) to obtain each new directional transform.  As an added
bonus, the magnitude envelopes of the HWT and DHWT wavelets are smooth
bumps, which makes the transforms nearly shift-invariant.


For different choices of $n$ and $m$ the HWT corresponds to some
interesting known and new wavelets, including dual-tree complex
wavelets ($n=1$) \cite{ivan}, new {\em quaternion wavelets}
($n=2$), new {\em octonion wavelets} ($n=3$), and so on. In 3-D,
octonion wavelets are the basic building blocks for constructing
both Selesnick and Li's CWT for $m=2$ (see Fig.\ref{fig:basis}(b))
and a new DHWT for $m=1$ (see Fig.\ref{fig:basis}(c)).


This paper is organized as follows.  In Sections \ref{sec:dwt} and
\ref{sec:cwt} we briefly review the DWT and CWT.  Section
\ref{sec:hwt} introduces the HWT and Section \ref{sec:dhwt}
develops the idea of orthogonal projections to obtain DHWTs.
Section \ref{sec:denoising} presents the results of  a number of
denoising experiments that demonstrate significant gains using
a DHWT with directionality matched to the manifold dimension at
hand. In the following, we will mainly focus on the cases $n=2$
and $n=3$ for clarity; however, all results generalize to
arbitrary $n$ and $m$.


\vspace{-0.1in}

\section{Discrete Wavelet Transform (DWT)}
\label{sec:dwt}

\vspace{-0.05in}

The discrete wavelet transform (DWT) represents a 1-D real signal
$f(t)$ in terms of shifted versions of a {\em scaling function}
$\phi(t)$ and shifted and scaled versions of a {\em wavelet
function} $\psi(t)$ \cite{Vetterli95wavelets}.  When
$\phi_{L,p}(t)=2^L\phi(2^L t -p)$ and
$\psi_{\ell,p}(t)=2^\ell\psi(2^\ell t -p)$ form an orthonormal
basis, we can represent any $f(t) \in L_2$ as
%
\vspace{-0.05in}
\begin{eqnarray}
f(t) & = & \sum_{p\in \mathbb{Z}} c_{L,p} \phi_{L,p}(t) +
\sum_{\ell\geq L,p\in \mathbb{Z}}
d_{\ell,p} \psi_{\ell,p}(t) \label{eqn:DWT1}
\end{eqnarray}
%\vspace{-0.1in}
%
where $c_{L,p} = \int f(t) \phi_{L,p}(t) dt$ and $d_{\ell,p} =
\int f(t) \psi_{\ell,p}(t) dt$ are the scaling and wavelet
coefficients.   $L$ sets the coarsest scale space that is spanned
by $\phi_{L,p}(t)$.  Behind each wavelet transform is a filterbank
based on a lowpass filter; we will use the notation $\phi_h(t)$,
$\psi_h(t)$ to denote the scaling and wavelet functions
corresponding to a particular filter $h$.

DWTs in higher dimensions are typically obtained using tensor
products of 1-D DWTs over each dimension.  In 2-D, for example, we
use the scaling function $\phi(x)\phi(y)$ and three wavelets
$\psi(x)\psi(y)$, $\phi(x)\psi(y)$, and $\psi(x)\phi(y)$ that are
oriented in the diagonal, horizontal, and vertical directions,
respectively \cite{ivan}.


\vspace{-0.1in}

\section{Complex Wavelet Transform (CWT)}
\label{sec:cwt}

\vspace{-0.05in}

The 1-D dual-tree CWT expands a real signal in terms of {\em two}
sets of wavelets and scaling functions obtained from two
independent filterbanks based on real filters $h$ and $g$
\cite{kingsbury,ivan}. $\psi_h(t)$ and $\psi_g(t)$ play the roles
of the real and imaginary parts of a complex analytic wavelet
$\psi^c(t) = \psi_h(t)+j\psi_g(t)$ that is supported only on
positive frequencies in the Fourier domain. Equivalently, the
imaginary wavelet is the {\em Hilbert transform} of the real
wavelet. The combined system is a $2\times$ redundant frame that,
by virtue of the fact that $|\psi^c(t)|$ is a smooth bump, is
shift-invariant.\footnote{In practice, in order to have
finite-length wavelets, the Hilbert transform is only
approximately satisfied, $\psi^c(t)$ is only approximately
analytic, and the CWT is only approximately shift-invariant
\cite{kingsbury,ivan}.}

It is useful to recall that the Fourier transforms of the Hilbert
pair of wavelets are related by $\Psi_g(f)=-j{\rm
sgn}(f)\Psi_h(f)$, where ${\rm sgn}(f)$ equals $1$ when $f>0$ and
$-1$ when $f<0$.  Thus, summing $\Psi_h(f)+j\Psi_g(f)=\Psi^c(f)$
cancels the negative frequency part of the complex wavelet,
rendering it analytic.

This notion of using 1-D Hilbert transforms to cancel out negative
frequencies has been used to create directional wavelets in 2-D
\cite{kingsbury}, 3-D, and beyond \cite{ivan}.  However, as we
will see, these CWTs are matched only to manifolds of intrinsic
dimension $m=n-1$.  To best analyze manifolds of lower intrinsic
dimension, we must turn to the theory of multidimensional Hilbert
transforms.


\vspace{-0.1in}

\section{Hypercomplex Wavelet Transform (HWT)}
\label{sec:hwt}

\begin{figure}[htb]
%\centerline{\epsfig{figure=QWT2D.eps,width=6.5cm}}

\begin{minipage}[b]{.48\linewidth}
  \centering
  \centerline{\epsfig{figure=QWT2D_a.eps,width=1.8cm}}
%  \vspace{1.5cm}
  \centerline{(a) $\Psi_h(f_x)\Psi_h(f_y)$}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\epsfig{figure=QWT2D_b.eps,width=1.8cm}}
%  \vspace{1.5cm}
  \centerline{(b) $\Psi_g(f_x)\Psi_h(f_y)$ }\medskip
\end{minipage}
\begin{minipage}[b]{.48\linewidth}
  \centering
  \centerline{\epsfig{figure=QWT2D_c.eps,width=1.8cm}}
% \centerline{\epsfig{figure=image3.ps,width=4.0cm}}
%  \vspace{1.5cm}
  \centerline{(c) $\Psi_h(f_x)\Psi_g(f_y)$}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\epsfig{figure=QWT2D_d.eps,width=1.8cm}}
% \centerline{\epsfig{figure=image4.ps,width=4.0cm}}
%  \vspace{1.5cm}
  \centerline{(d) $\Psi_g(f_x)\Psi_g(f_y)$}\medskip
\end{minipage}


\vspace{-0.2in}

\caption{\sl Fourier-domain relationships among the four
components of a 2-D HWT wavelet $\Psi(f_x,f_y)$ in the diagonal
subband.} \label{fig:QWT2D} \vspace{-0.22in}
\end{figure}


\vspace{-0.05in}


There are no unique definitions of Hilbert transform and analytic
signal in 2-D and higher. To illustrate two possibilities,
consider a 2-D wavelet Fourier transform $\Psi(f_x,f_y)$.  The
straightforward extension of the 1-D definition sets
$\Psi(f_x,f_y)=0$ on all but a half-plane ($f_x>0$, for example)
and uses standard complex algebra for manipulation.  Another
natural definition sets $\Psi(f_x,f_y)=0$ on all but a quadrant
($f_x,f_y>0$, for example) and uses {\em quaternion} algebra for
manipulation \cite{bulow_book11}.\footnote{The set of quaternions
$\mathbb{H} = \{a+j_1b+j_2c+j_3d|a,b,c,d \in \mathbb{R}\}$ with
multiplication rules $j_1j_2=-j_2j_1=j_3$ and $j_1^2=j_2^2=-1$
\cite{kantor}.}

Adopting the latter rather than the former definition, we will now
introduce a new multiscale signal representation: the dual-tree
{\em hypercomplex wavelet transform} (HWT). For clarity, we will
stick to 2-D, but the HWT extends easily to any dimension $n$.


Each 2-D hypercomplex wavelet consists of a standard DWT tensor
wavelet plus three additional real wavelets obtained from the
three combinations of 1-D Hilbert transforms along the two
coordinates. More specifically, denote the 1-D Hilbert transform
operators along the $x$ and $y$ coordinates as ${\cal H}_x$ and
${\cal H}_y$, respectively. Then given the diagonal tensor product
wavelet $\psi_h(x)\psi_h(y)$ from Section \ref{sec:dwt}, we
complement it with
%
\begin{eqnarray}
{\cal H}_x \{\psi_h(x)\psi_h(y)\} &=& \psi_g(x)\psi_h(y), \label{eqn:Hx} \\
{\cal H}_y \{\psi_h(x)\psi_h(y)\} &=& \psi_h(x)\psi_g(y), \label{eqn:Hy} \\
{\cal H}_y {\cal H}_x \{\psi_h(x)\psi_h(y)\} &=&
\psi_g(x)\psi_g(y). \label{eqn:Hxy}  \vspace{-0.1in}
\end{eqnarray}
%
Conveniently, each component can be computed as a combination of
1-D dual-tree complex wavelets. These three components closely
resemble $\psi_h(x)\psi_h(y)$ but are phase-shifted by $90^\circ$
in the horizontal, vertical, and diagonal directions,
respectively; Fig.\ \ref{fig:QWT2D} interprets these relationships
in the Fourier domain. Using quaternion algebra, we can compactly
summarize the above process in a {\em quaternion wavelet}
$\psi^q(x,y) = \psi_h(x)\psi_h(y) - j_1\psi_g(x)\psi_h(y) -
j_2\psi_h(x)\psi_g(y) + j_3\psi_g(x)\psi_g(y)$.  The construction
of the other two quaternion wavelets based on $\phi_h(x)\psi_h(y)$
and $\psi_h(x)\phi_h(y)$ is identical. The resulting HWT, which we
term the dual-tree {\em quaternion wavelet transform} (QWT) in
this special 2-D case, is a $4\times$ redundant tight frame.


Generalization to 3-D and higher is straightforward. In $n$-D
there are $2^n-1$ different hypercomplex wavelets. Each wavelet
consists of a standard DWT tensor wavelet plus $2^n-1$ additional
real wavelets obtained from the $2^n-1$ combinations of 1-D
Hilbert transforms along the $n$ coordinates. The total redundancy
is $2^n\times$.

\begin{figure}[htb]
%\centerline{\epsfig{figure=CWT2D.eps,width=5.5cm}}
\begin{minipage}[b]{.48\linewidth}
  \centering
  \centerline{\epsfig{figure=CWT2D_a.eps,width=1.8cm}}
%  \vspace{1.5cm}
  \centerline{(a) $\mathscr{F}({\rm Re}(\psi^c(x,y)))$}\medskip
\end{minipage}
\hfill
\begin{minipage}[b]{0.48\linewidth}
  \centering
  \centerline{\epsfig{figure=CWT2D_b.eps,width=1.8cm}}
%  \vspace{1.5cm}
  \centerline{(b) $\mathscr{F}({\rm Im}(\psi^c(x,y)))$}\medskip
\end{minipage}
\vspace{-0.2in} \caption{\sl Fourier-domain relationships between
the real and imaginary parts of a 2-D diagonally oriented complex
wavelet. $\mathscr{F}(\cdot)$ is the Fourier Transform.}\label{fig:CWT2D} \vspace{-0.25in}
\end{figure}

Since we construct the HWT basis by Hilbert transforming a
standard DWT basis, the HWT basis has no better directional
selectivity than the DWT basis. However, just as the 1-D Hilbert
transform enables us to zero-out the negative frequency axis to
form an analytic signal, the $n$-D hypercomplex Hilbert transform
allows us zero-out all but one ``$n$-tant'' \footnote{The natural
$n$-D generalization of a quadrant in 2-D and octant in 3-D.} in
the Fourier domain For example, to obtain a single diagonally
oriented 2-D wavelet whose frequency support is focused in the
second quadrant of the Fourier plane, we combine the four
quaternion components as
$(\psi_h(x)\psi_h(y)+\psi_g(x)\psi_g(y))+j(\psi_h(x)\psi_g(y)-\psi_g(x)\psi_h(y))$.
(See Fig.\ \ref{fig:QWT2D}, and also note that this is an
orthogonal projection.) This capability will prove extremely
useful for constructing directional wavelets in the next section.

%\vspace{-0.1in}

\vspace{-.1in}

\section{Directional HWT (DHWT)}
\label{sec:dhwt}

\vspace{-.05in}

The HWT provides a convenient set of building blocks for
constructing new directional wavelet transforms.  The idea is
simple: we can generate a new wavelet with localized Fourier
support, and hence directionality \cite{ivan}, by
carefully recombining the $2^n$ components of an $n$-D HWT
wavelet. The idea is also powerful: with the HWT we can generate
$n$-D wavelets that are bandpass in subspaces of arbitrary
dimension $m$.


\vspace{-0.1in}

\subsection{DHWT in 2-D}
\label{subsec:2Ddhwt}

\vspace{-0.05in}

We illustrate the procedure by first re-deriving the 2-D CWT of
\cite{kingsbury,ivan} from the QWT of last section.  The
Fourier-domain support of each 2-D QWT wavelet consists of four
dyadic squares symmetrically placed with respect to the origin.
This limits its directionality. For example, the Fourier support
of the diagonally oriented QWT wavelet lies at the four corners of
a square in the Fourier plane (see Fig.~\ref{fig:QWT2D}(a)); hence
this wavelet is localized simultaneously along two lines of slope
$\pm 45^\circ$. To obtain a wavelet directed along a single line
(at either $+45^\circ$ or $-45^\circ$) we need to partition the
Fourier plane in half to isolate pairs of squares that lie on a
straight line through the origin. Since we can isolate individual
quadrants by recombining the QWT components (recall the end of the
Section \ref{sec:hwt}), we can generate a wavelet with
two-quadrant Fourier support as in Fig.~\ref{fig:CWT2D}(a).  In
the signal domain, this wavelet will be lowpass along the
$+45^\circ$ direction and bandpass along the orthogonal
$-45^\circ$ direction \cite{kingsbury,ivan}.

Repeating the same procedure for all three QWT wavelets with two
directions per subband, we obtain six directional basis functions,
each of which is a complex analytic signal. This is precisely the
2-D CWT \cite{kingsbury,ivan}.  Since the basis functions are
analytic, the CWT is shift-invariant.

\vspace{-0.1in}

\subsection{DHWT in 3-D}
\label{subsec:3Ddhwt}

\vspace{-0.05in}

The 3-D case is more interesting because there exist two entirely
different DHWTs corresponding to $m=1,2$, and one is new.  Let us
focus on the ``diagonal'' HWT wavelet whose Fourier support is
located in the eight corner cubes in 3-D Fourier space (see Fig.\
\ref{fig:cube}(a)). This dual-tree HWT basis function has eight
{\em octonion} components:\vspace{-0.1in}
\begin{eqnarray}
\psi_{1,1}=\psi_h(x)\psi_h(y)\psi_h(z), &
\psi_{2,1}=\psi_g(x)\psi_g(y)\psi_g(z), \nonumber \\
\psi_{1,2}=\psi_g(x)\psi_g(y)\psi_h(z), &
\psi_{2,2}=\psi_h(x)\psi_h(y)\psi_g(z), \nonumber \\
\psi_{1,3}=\psi_h(x)\psi_g(y)\psi_g(z), &
\psi_{2,3}=\psi_g(x)\psi_h(y)\psi_h(z), \nonumber \\
\psi_{1,4}=\psi_g(x)\psi_h(y)\psi_g(z), &
\psi_{2,4}=\psi_h(x)\psi_g(y)\psi_h(z). \label{eqn:OWTmatrix}
%\vspace{-0.5in}
\end{eqnarray}
%\vspace{-0.2in}
%
%\begin{equation}
%\begin{bmatrix}
%\psi_h(x)\psi_h(y)\psi_h(z) \\
%\psi_g(x)\psi_g(y)\psi_h(z) \\
%\psi_g(x)\psi_h(y)\psi_g(z) \\
%\psi_h(x)\psi_g(y)\psi_g(z) \\
%\psi_g(x)\psi_h(y)\psi_h(z) \\
%\psi_h(x)\psi_g(y)\psi_h(z) \\
%\psi_h(x)\psi_h(y)\psi_g(z) \\
%\psi_g(x)\psi_g(y)\psi_g(z)
%\end{bmatrix}
%\label{eqn:OWTmatrix}
%\end{equation}
%\vspace{-0.05in}
%
As in 2-D, we can linearly combine these eight components to
obtain a wavelet with spectral support in a single octant. We now
have two choices, depending on whether $m=1$ or 2.

\noindent {\bf Case $m=2$:}  To form each DHWT wavelet, we group a
pair of cubes located on a 1-D line running through the origin
(see Fig.\ \ref{fig:cube}(b)).  The resulting wavelet is complex
and bandpass along the 1-D line\footnote{Due to the rotation
property of the Fourier transform, the wavelet will actually be
bandpass along a 1-D line in the space domain that is rotated
$90^\circ$ with respect to the Fourier-domain line.} and lowpass
in the orthogonal direction, which is a 2-D plane through the
origin (see Fig.\ \ref{fig:basis}(b)).  This coincides with the
3-D CWT of Selesnick and Li \cite{ivan}.

\noindent {\bf Case $m=1$:}  Alternatively, we group sets of four
cubes located on a 2-D plane passing through the origin (see Fig.\
\ref{fig:cube}(c)). The resulting wavelet is bandpass along the
2-D plane and lowpass in the orthogonal direction, which is a 1-D
line through the origin (see Fig.\ \ref{fig:basis}(c)).

Due to the way the four cubes are combined -- for example as
$\psi_{1,1} \pm \psi_{1,3}$ and $\psi_{2,1} \pm \psi_{2,3}$ -- the
wavelets behave as a quaternion in the direction parallel to the
2-D plane.  This yields a smooth magnitude envelope and
shift-invariant transform.  A similar recipe works for the other
directions.

Starting from the HWT wavelet in Fig.\ \ref{fig:basis}(a), there
are six ways to group four corner cubes on 2-D planes through the
origin.  This result holds for each of the three HWT wavelets.
Therefore, the $m=1$ DHWT basis that covers all possible 3-D
directions has 18 such directions and thus three times the
redundancy of the $m=2$ DHWT ($24\times$ redundant in total).

The same construction applies for arbitrary $n$ and $m$.

\vspace{-.1in}

\section{DENOISING EXPERIMENTS}
\label{sec:denoising}

\vspace{-.05in}

To evaluate the performance of the DHWT, we performed a number of
denoising experiments in 3-D using data of different intrinsic
manifold dimension and wavelets matched to $m=0,1,2$.

The first three test signals have $m=1$:  {\em Stars} is a
synthetic video that pans across a field of stars and galaxies
(size: $64\times64\times 64$); {\em Rods} contains solid rods of
various radii and intensities (size: $64\times 64\times 64$);
{\em Helix} maps out a smooth 1-D helical curve (size: $64\times
64\times 64$).  The final test signal has $m=2$: {\em Sheet} is
a thin tilted plane (size: $64\times 64\times
64$).\footnote{Test signals are available from
dsp.rice.edu/Software/HWT.}


We added white Gaussian noise of $\sigma=15$ to each data set and
applied optimal hard thresholding following the program in
\cite{ivan} (see taco.poly.edu/WaveletSoftware).  We used the
dual-tree wavelet transform from \cite{ivan} with length-10
wavelets.  We measured the log-mean-square (LMS) error of the
denoised data; LMS errors before denoising were 65.58dB, 76.43dB,
69.67dB, and 77.71dB, respectively.


\begin{table}
\begin{center}
\caption{\sl 3-D denoising results in log-MSE.}
\label{tab:denoise} \vspace{.02in}
\begin{tabular}{|c|c|c|c|} \hline
%Test Data & SER (3-D QWT) & SER (3-D CWT) \\ \hline
test & CWT  &  DHWT & DWT  \\[-1mm]
%
signal &  ($m=2$) & ($m=1$) & ($m=0$) \\ \hline
Stars & 60.96 dB & \bf{53.82 dB} & 62.84 dB \\
Rods & 62.77 dB & \bf{60.24 dB} & 64.97 dB \\
Helix & 63.68 dB & \bf{61.02 dB} & 64.98 dB \\
Sheet & 59.35 dB & \bf{52.35 dB} & 61.67 dB \\ \hline
%1 & 12.92 dB & 5.78 dB & 3.90 dB\\
%2 & 6.50 dB & 3.97 dB & 1.77 dB \\
%3 & 14.39 dB & 7.39 dB & 5.07 dB \\
%4 & 5.72 dB & 3.06 dB & 1.76 dB \\ \hline
\end{tabular}
\end{center}
\vspace{-0.2in}
\end{table}

From Table \ref{tab:denoise} we observe the superiority of the
$m=1$ DHWT over the $m=2$ CWT and the usual tensor product DWT. We
are surprised that the $m=2$ CWT does not dominate on the sheet
data, since it should theoretically be better matched to the
intrinsic manifold dimension.  Perhaps the difference is due to
the extra averaging available in the $3\times$ more redundant
$m=1$ transform.

\vspace{-.1in}

\section{CONCLUSIONS}
\label{sec:conclusions}

\vspace{-.05in}

We have introduced a new family of multiscale signal
representations, the directional hypercomplex wavelet transforms,
which show considerable promise for multidimensional signal
processing.  DHWTs are especially designed for signals with a
local manifold structure, precisely where standard tensor product
wavelets have often disappointed. Their superior performance on a
number of 3-D denoising experiments suggests that they have
excellent sparsification capabilities.  Perhaps the main challenge
to both CWT and DHWT methods will be their redundancy, since
multidimensional data processing problems are already often
computationally challenging.  Certainly controlling redundancy
while preserving directionality is a fertile area for future
research.


\nocite{bulow_HFT}

% 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{mybiblio}


%\newpage
%\section{Final stuff to finish it off}

%\begin{itemize}

%\item careful spell check

%\item richb added ``real part'' of wavelets in figs 1 and 2.  is
%this correct?

%\item fix the citations

%\item make sure all [?] and (?), if any, are fixed

%\item make sure color figures (eg: 1, 2) are legible when printed
%5on a black and white printer.

%\item is there a word for what i call ``$n$-tant'' -- see bulow
%book.

%\item how BIG are the data test sets?

%\item get all the figures so they are at the top of the page

%\item columns of the table should be reordered ideally to go
%$m=0,1,2$ or $m=2,1,0$.  not a show stopper tho.

%\item check length of wavelets used in denoising expts.  i guessed
%10.

%\item comment end of section 5.1

%\item make sure that NOTHING (words, figure, ...) sticks out past
%the edge of the column.  if words do, then put \sloppy at the
%start of the para.

%\item references to fix:

%4% 5 10 put a space between the initials, ie:  N.G. becomes N. G.

%4 capitalize {D}aubechies

%6 isn't that conference wavelets X?  I would also put that volume
%number on [7] i think

%7 no capitals in the title of the paper

%2 3 9 book titles should have each word capitalized

%9 no need for "Springer series in ... sciences".  also chapter
%becomes ch.

%\end{itemize}

\end{document}






%The quaternion wavelet
%transform (QWT) of a 2-D signal $f(x,y)$ is given as
%\begin{equation}
%d_\ell(p) = \int_{\mathbb{R}^2} \psi_{\ell,p}^i(x) f(x,y)
%\psi_{\ell,p}^j(y) dx dy \label{eq:qwt}
%\end{equation}
%where $\psi_{\ell,p}^i(x) = \psi_{\ell,p}^h(x) - i\psi_{\ell,p}^g(x)$ and
%$\psi_{\ell,p}^j(y) = \psi_{\ell,p}^h(y) - j\psi_{\ell,p}^g(y)$, since
%$\psi^q(x,y)=\psi^i(x)\psi^j(y)$.

%Figure~\ref{fig:QWT2Dbasis} shows the 4 components of QWT in the
%diagonal subband and its magnitude. Note that the 4 functions are
%obtained by $90^o$ phase shift of the leftmost function in $x$
%direction, in $y$ direction, and both. The magnitude of the QWT
%basis is not wiggly as shown in the rightmost plot, indicating
%shift invariance.
%\begin{figure}[htb]
%\centerline{\epsfig{figure=QWT2Dbasis.eps,width=8.5cm}}
%\caption{\sl One set of 2-D quaternionic wavelet basis. Left to right:
%$\psi_h(x)\psi_h(y)$, $\psi_h(x)\psi_g(y)$, $\psi_g(x)\psi_h(y)$,
%$\psi_g(x)\psi_g(y)$, $|\psi^q(x,y)|$.} \label{fig:QWT2Dbasis}
%\end{figure}


%[*** refine, reorganize following paragraphs ***]





%The inversion of the 2-D hypercomplex wavelet transform is stable
%because the transform basis functions form a tight frame with
%frame constant 4 when $\psi_h(t)$ and $\psi_g(t)$ form orthonormal
%basis in 1-D. In addition, for each quaternion wavelet
%coefficient, we can obtain 3 phase angles for a localized image
%block. Two of those phase angles are directly related to the local
%signal shift in the two coordinates. The third phase angle is
%related to the local signal structure \cite{bulow_book11}.





Just as the two definitions of Hilbert transform and analytic
signal are at extremes (half-spaces vs.\ quadrants, octants, and
so on), there is a range of different hypercomplex wavelets
corresponding to these in bewteen.  However, the hypercomplex is
the most general, since you can construct any other by tiling in
frequency.  In 2-D for example, you obtain the right half-plane
(traditional definition of Hilbert transform) by tiling it with
two quadrants.  This corresponds to building orthogonal projection
of the HWT wavelets.
