Fast directional continuous spherical wavelet transform algorithms
J.~D.~McEwen
~M.~P.~Hobson
~D.~J.~Mortlock
~A.~N.~Lasenby% % Michael~Shell
~\IEEEmembership{Member,~IEEE,} % John~Doe
~\IEEEmembership{Fellow,~OSA,} % and~Jane~Doe
We describe the construction of a spherical wavelet analysis through the inverse stereographic projection of the Euclidean planar wavelet framework, introduced originally by Antoine and Vandergheynst and developed further by Wiaux \etal. Fast algorithms for performing the directional continuous wavelet analysis on the unit sphere are presented. The fast directional algorithm, based on the fast spherical convolution algorithm developed by Wandelt and G\'{o}rski, provides a saving of $\order(\sqrt{\npix})$ over a direct quadrature implementation for \npix\ pixels on the sphere, and allows one to perform a directional spherical wavelet analysis of a $10^6$ pixel map on a personal computer. %reduces typical execution time on a 3GHz laptop with 512MB %of memory from almost one hour to a fraction of a second for a %pixelisation containing approximately 50,000 pixels.
fcswt_abstract.tex
%------------------------------------------------------------------------------
% Abstract
%------------------------------------------------------------------------------

\begin{abstract}
We describe the construction of a spherical wavelet analysis through
the inverse stereographic projection of the Euclidean planar wavelet
framework, introduced originally by Antoine and Vandergheynst and developed further by
Wiaux \etal.  
Fast algorithms for performing the directional continuous wavelet
analysis on the unit sphere are presented.  The fast directional algorithm,
based on the fast spherical convolution algorithm developed by Wandelt
and G\'{o}rski,
provides a saving of $\order(\sqrt{\npix})$ over a direct quadrature
implementation for \npix\ pixels on the sphere, and 
allows one to perform a directional spherical wavelet analysis of a 
$10^6$ pixel map on a personal computer.
%reduces typical execution time on a 3GHz laptop with 512MB
%of memory from almost one hour to a fraction of a second for a
%pixelisation containing approximately 50,000 pixels.
\end{abstract}

fcswt_body.tex
%------------------------------------------------------------------------------
% Introduction
%------------------------------------------------------------------------------
\section{Introduction}
\label{sec:intro}

\PARstart{W}{avelet} analysis has proven useful in many
applications due to 
the ability of wavelets to resolve localised signal content in both
scale and space.  Many of these applications, however, are restricted
to data defined in Euclidean space: the 1-dimensional line, 
the 2-dimensional plane and,
occasionally, higher dimensions.  Nevertheless, data are often measured
or defined on other manifolds, such as the 2-sphere.  
Examples where data are measured on the sphere are found in 
astrophysics (\eg\ \cite{bennett:1996,bennett:2003}), 
planetary science (\eg\ \cite{wieczorek:2006,wieczorek:1998,turcotte:1981}), 
geophysics (\eg\ \cite{whaler:1994,swenson:2002,simons:2006}), 
computer vision (\eg\ \cite{ramamoorthi:2004}) and 
quantum chemistry (\eg\ \cite{choi:1999,ritchie:1999}).
% Some examples
% include astrophysical observations, geophysical data and reflectance
% functions for computer vision, all of which are defined on the sphere 
% (see \cite{healy:2003} for additional examples and references). 
% \fbox{give references and more examples}
% 
To realise the potential benefits that may be provided by wavelets
in such settings, ordinary Euclidean wavelet analysis must be extended
to spherical geometry.

A number of attempts have been made to extend wavelets to the unit sphere.  
Discrete second generation wavelets on the sphere that are based on a multiresolution analysis have been developed \cite{schroder:1995,sweldens:1996}.  Other authors have focused on the continuous wavelet transform on the sphere.  A number of works construct a solution using a harmonic approach \cite{narcowich:1996,potts:1996,freeden:1997a,freeden:1997b}, however these solutions suffer from the poor localisation of the spherical harmonic functions.  Others adopt a tangent bundle viewpoint \cite{torresani:1995,dahlke:1996}, thereby avoiding the necessity to define a dilation operator on the sphere.  A satisfactory extension of the continuous wavelet transform to the sphere is defined by \cite{holschneider:1996}, however this construction requires an abstract dilation parameter that must satisfy a number of {\it ad hoc} assumptions.
% 
% (a brief review of which is given in the introductions of
% \cite{antoine:1999,antoine:1998}).  However, many of these constructions
% are not fully satisfactory, missing important elements of the
% generic wavelet framework.  
% 
More recently, a consistent and
satisfactory framework for wavelets defined on the unit sphere has been
constructed and developed by 
\cite{antoine:1999,antoine:1998,antoine:2002,antoine:2004,%
bogdanova:2004,demanet:2003,wiaux:2005}.  Moreover, this construction is derived entirely from group theoretic principles and inherently satisfies a number of natural requirements.
We consider the continuous spherical wavelet transform (\cswt)
developed in these last works.  
For a more detailed review of the attempts made at constructing a wavelet transform on the unit sphere see \cite{antoine:1999,antoine:1998,mhaskar:2001}.

Current and future data-sets defined on the sphere are of
considerable size.  The current Wilkinson Microwave Anisotropy Probe (WMAP) data of the cosmic microwave
background (CMB) contain approximately $3\times10^6$ pixels on the
sphere, whereas the forthcoming Planck CMB mission will generate maps with  
approximately $50\times10^6$ pixels.  
Fast algorithms are therefore required to perform the \cswt\ on
practical data-sets.
A semi-fast algorithm to implement the \cswt\ is presented in
\cite{antoine:2002} and implemented in the 
\yawtb\footnote{\url{http://www.fyma.ucl.ac.be/projects/yawtb/}}
(Yet-Another-Wavelet-Toolbox) Matlab wavelet
toolbox (which also makes use of the SpharmonicKit\footnote{\url{http://www.cs.dartmouth.edu/~geelong/sphere/}}).
However, this implementation is restricted to an equi-angular
tessellation of the sphere.  The beauty of this tessellation is its
simplicity and ability to be easily represented in matrix form.
However, the pixels of an equi-angular tessellation are densely spaced
about the poles and do not have equal areas. 
Other tessellations of the sphere also exist, such as those constructed to minimise some energy measure \cite{saff:1997,hardin:2004,du:2003} or those constructed for more practical or numerical purposes (for example
the  IGLOO\footnote{\url{http://www.mrao.cam.ac.uk/projects/cpac/igloo/}}
\cite{crittenden:1998},
\healpix\footnote{\url{http://healpix.jpl.nasa.gov/}} 
\cite{gorski:2005}
and GLESP\footnote{\url{http://www.glesp.nbi.dk/}} \cite{doroshkevich:2005}
schemes).
%Ideally, we desire a fast implementation of the \cswt\ that is
%tessellation independent. 
There is thus a need for a fast implementation of the
\cswt\ that is 
not tied to any particular tessellation of the sphere.

We fill this void by presenting a fast algorithm for performing
the directional \cswt. 
The \cswt\ at a particular scale is essentially a
spherical convolution, hence we may apply the fast spherical convolution
algorithm developed
by \cite{wandelt:2001} to evaluate the wavelet transform.
The algorithm is posed in harmonic space and thus is
independent of the underlying tessellation of the sphere, (although an
iso-latitude tessellation does enable faster spherical harmonic
transforms, thereby increasing the speed of the algorithm). 
The framework supports both non-azimuthally symmetric spherical wavelets\footnote{Azimuthally symmetric spherical wavelets are also often referred to as axisymmetric wavelets.} and a decomposition that employs anisotropic dilations, however no synthesis is possible when anisotropic dilations are incorporated.
%We have implemented our fast algorithm (and others) in Fortran 90.
For an illustration of a spherical wavelet analysis on a practical
problem of considerable size we refer the reader to our recent works to
test the CMB for deviations from Gaussianity \cite{mcewen:2004,mcewen:2006b} and to detect the \iswtext\ (\isw) effect \cite{mcewen:2006}.
Both of these works involve performing 1000 Monte Carlo simulations and  
would not have been feasible without a fast directional \cswt\ algorithm.

The remainder of this paper is structured as follows.  
In \sectn{\ref{sec:cswt}} we describe the \cswt\ in the framework
presented by \cite{wiaux:2005}.  
For a more complete treatment of the spherical wavelet transform in
this framework and the correspondence between spherical and Euclidean
wavelets we recommend that the reader refer to \cite{wiaux:2005}.
We also present an extension to anisotropic dilations, however in this case the basis functions are not strictly wavelets hence perfect reconstruction is not possible.
% extend the spherical wavelet framework described in \cite{wiaux:2005} 
% to the more general case of anisotropic dilations in \sectn{\ref{sec:cswt}}.
%
% We also describe in \sectn{\ref{sec:cswt}} the extension of the spherical dilation 
% operator described in \cite{wiaux:2005} to the more general anisotropic case.
% although we do extend the spherical dilation operator to the more general
% anisotropic case.  
Various algorithms to perform the \cswt\ are
described and then compared in \sectn{\ref{sec:algorithms}}.
An application of our implementation is demonstrated in \sectn{\ref{sec:application}}.
Concluding remarks are made in \sectn{\ref{sec:conclusions}}.



%------------------------------------------------------------------------------
% The continuous spherical wavelet transform
%------------------------------------------------------------------------------
\section{The continuous spherical wavelet transform}
\label{sec:cswt}

To extend Euclidean wavelet analysis to spherical geometry a number of
requirements must be satisfied: (i) the signals and wavelets
must live fully on the unit sphere; (ii) the transform must involve local
dilations of some kind on the unit sphere; and (iii) the spherical wavelet
transform  should reduce locally to the Euclidean transform on the
tangent plane (\ie\ the Euclidean limit must be satisfied)
\cite{holschneider:1996}.
The final requirement is intuitively obvious; the sphere is asymptotically
flat, hence any spherical wavelet transform should match the planar
Euclidean transform on small scales, or equivalently, for a large
radius of curvature.
%  (by, for example, an analogous approximation to the spherical harmonics by Bessel functions, \eg\ \cite{simons:2006}). \fbox{Ok?}

The spherical wavelet transform developed in
\cite{antoine:1999,antoine:1998,antoine:2002,antoine:2004,bogdanova:2004,demanet:2003}
satisfies all of these requirements, moreover each requirement
naturally follows from the construction.  The construction of this
transform is derived entirely from group theoretical principles.
However, in a recent work by \cite{wiaux:2005} this formalism is
reintroduced independently of the original group theoretic
formalism, in an equivalent, practical and self-consistent approach. 
We adopt this approach herein.
%and extend it to anisotropic dilations on the
% sphere in two orthogonal directions.

The correspondence principle between spherical and Euclidean wavelets
is developed by \cite{wiaux:2005}, relating the concepts of planar
Euclidean wavelets to spherical wavelets through a stereographic
projection.  We use the stereographic projection to define affine
transformations on the unit sphere that facilitate the construction of a
wavelet basis on the unit sphere.  The spherical wavelet transform may then
be defined as the projection on to this basis, where the spherical
wavelets must satisfy the appropriate admissibility criterion to
ensure perfect reconstruction.


%------------------------------------------------------------------------------
\subsection{Stereographic projection}

In order to construct a correspondence between wavelets on the plane
($\realno^2$) and sphere (\sphere) a projection operator between the
two spaces must be chosen.  
It is shown in \cite{wiaux:2005} that the stereographic projection is the unique unitary, radial and conformal diffeomorphism between the sphere and the plane.
% The stereographic projection is not
% necessarily a unique choice, nevertheless it provides a natural
% definition of dilations on the unit sphere (as shown in 
% \sectn{\ref{sec:affine_transformations}}) and also preserves
% directionality.  

The stereographic projection is defined by projecting a point on the
unit sphere to a point on the tangent plane at the north pole, by casting a
ray though the point and the south pole.  The point on the unit sphere is
mapped on to the intersection of this ray and the tangent plane (see
\fig{\ref{fig:stereographic_projection}}).  
% Notice that under this
% projection 
% %the orientational \eulc\ angle on the unit sphere is identical
% %to the polar angle on the plane, hence 
% directionality is preserved.
Formally, we may define the stereographic projection operator as 
$\spo : \sa \rightarrow \xvect = \spo \sa = (\rad(\saa),\sab)$ where 
$\rad=2 \tan(\saa/2)$, $\sa \equiv (\sas) \in \sphere$ denotes
spherical coordinates with colatitude \saa\ and longitude \sab\
and $\xvect \in \realno^2$ is a point in the plane, denoted here by
the polar coordinates $(\rad,\sab)$.
The inverse operator is 
$\spo^{-1} : \xvect \rightarrow \sa = \spo^{-1} \xvect =
(\saa(\rad),\sab)$,
where $\saa(\rad)=2\tan^{-1}(\rad/2)$.

Following the formulation of \cite{wiaux:2005} again, we define the
action of the stereographic projection operator on functions on the
plane and sphere.  We consider the space of square integrable
functions in $L^2(\realno^2,\dx^2\xvect)$ on the plane and 
$L^2(\sphere,\dmu)$ on the unit sphere, where 
$\dmu = \sin\saa \dx\saa \dx\sab$ is the usual rotation invariant measure
on the sphere. 
The action of the stereographic projection operator 
$\spo : s \in L^2(\sphere,\dmu) \rightarrow p=\spo s \in L^2(\realno^2,\dx^2\xvect)$
on functions is defined as
\begin{equation}
\label{eqn:sp}
p(\rad,\sab) = (\spo s)(\rad,\sab) =
(1 + r^2/4)^{-1}
s(\saa(\rad),\sab)
\spcend .
\end{equation}
The inverse stereographic projection operator 
$\spo^{-1} : p \in L^2(\realno^2,\dx^2\xvect) \rightarrow s=\spo^{-1}
p \in L^2(\sphere,\dmu)$
on functions is then
\begin{equation}
\label{eqn:isp}
s(\sas) = (\spo^{-1} p)(\sas) =
[1 + \tan^2(\saa/2)]
p(\rad(\saa),\sab)
\spcend .
\end{equation}
The pre-factors introduced ensure that the $L^2$-norm of functions
through the forward and inverse projections are conserved. 
In the Euclidean limit, the stereographic projection and inverse
naturally reduce to the identity operator \cite{antoine:1999}. 

\begin{figure}
%\scriptsize
%\psfrag{x}{$x$}
%\psfrag{y}{$y$}
%\psfrag{z}{$z$}
%\psfrag{t}{$2\tan(\saa/2)$}
%\psfrag{a}{$\saa$}
\centerline{
  \includegraphics[width=86mm]{figures/stereographic_projection_3d_psfrag}}
\caption[Stereographic projection]{Stereographic projection of the
  sphere onto the plane.}
\label{fig:stereographic_projection}
\end{figure}

%------------------------------------------------------------------------------
\subsection{Affine transformations on the sphere}
%\subsection{Wavelet basis}
\label{sec:affine_transformations}

A wavelet basis is constructed on the unit sphere in
\sectn{\ref{sec:wavelet_transform}} by applying the 
spherical extension of Euclidean translations and dilations 
to mother wavelets defined on the unit sphere.  The extension of these
affine transformations to the sphere, facilitated by the stereographic
projection operator, are defined here.

The natural extension of Euclidean translations on the unit sphere are rotations.
These are characterised by the elements of the rotation group \sothree,
which we parameterise in terms of the three Euler angles
$\rho=(\euls)$.\footnote{We adopt the $zyz$ Euler convention corresponding to the rotation of a physical body in a \emph{fixed} co-ordinate system about the $z$, $y$ and $z$ axes by \eulc, \eulb\ and \eula\ respectively.}  The rotation of a square-integrable function 
$s\in L^2(\sphere,\dmu)$ is defined by
\begin{equation}
[\rot(\rho) s](\sa) = s(\rho^{-1} \sa), \; \; \rho \in \sothree 
\spcend .
\end{equation}

Dilations on the unit sphere are constructed by first
lifting the sphere to the plane by the stereographic projection,
\mbox{followed} by the usual Euclidean dilation in the plane, before
re-projecting back onto the sphere.  We generalise to anisotropic dilations on the sphere (a similar anisotropic dilation operator on the sphere has been independently proposed by \cite{tosic:2005}), however in this setting we do not achieve a wavelet basis and hence cannot synthesise our original signal.
We define the anisotropic
Euclidean dilation operator in $L^2(\realno^2,\dx^2\xvect)$ as
\begin{equation}
[\dilsmall(\scalea,\scaleb) p] (x,y) = \scalea^{-1/2}\scaleb^{-1/2} \:
p(\scalea^{-1}x,\scaleb^{-1}y)
\spcend ,
\end{equation}
for the non-zero positive scales 
$\scalea, \scaleb \in \realno^{+}_{\ast}$.
The $\scalea^{-1/2}\scaleb^{-1/2}$ normalisation factor ensures the
$L^2$-norm is preserved.  The spherical dilation operator 
$\dil(\scaleab) : s(\sas) \rightarrow [\dil(\scaleab)s](\sas)$
in $L^2(\sphere,\dmu)$ is defined as the conjugation by \spo\ of
the Euclidean dilation $d(\scaleab)$ in $L^2(\realno^2,\dx^2\xvect)$
on the tangent plane at the north pole:
\begin{equation}
\dil(\scaleab) = \spo^{-1} \, \dilsmall(\scaleab) \, \spo
\spcend .
\end{equation}
The norm of functions in $L^2(\sphere,\dmu)$ is preserved by the
spherical dilation as both the stereographic projection operator and
Euclidean dilations preserve the norm of functions.
Extending the isotropic spherical dilation operator defined by 
\cite{wiaux:2005} to anisotropic dilations, we obtain
\begin{equation}
[\dil(\scaleab) s](\sa)
= [\cocycle(\scaleab,\sas)]^{1/2} \: s(\sa_{1/\scalea,1/\scaleb})
\spcend ,
\end{equation}
where $\sa_\scaleab=(\saa_\scaleab,\sab_\scaleab)$,
\begin{displaymath}
\tan(\saa_\scaleab/2) = 
\tan(\saa/2)
\sqrt{\scalea^2 \cos^2{\sab} + \scaleb^2 \sin^2{\sab}}
\end{displaymath}
and 
\begin{displaymath}
\tan(\sab_\scaleab) = 
\frac{\scaleb}{\scalea}
\tan(\sab)
\spcend .
\end{displaymath}
For the case where $\scalea=\scaleb$ the anisotropic
dilation reduces to the usual isotropic case defined by 
$\tan(\saa_\scalea/2) = \scalea \tan(\saa/2)$ and
$\sab_\scalea=\sab$.  
The $\cocycle(\scaleab,\sas)$ cocycle term follows from the factors
introduced in the stereographic projection of functions to preserve
the $L^2$-norm.  Alternatively, the cocycle may be derived
explicitly to preserve the $L^2$-norm when the stereographic projection
of functions do not have these pre-factor terms.  The cocycle of an
anisotropic spherical dilation is defined by
\begin{equation}
\cocycle(\scaleab, \sas) =
\frac{4 \scalea^3 \scaleb^3}
{ 
\left(
A_{-}\cos\saa + A_{+}
\right)^2
}
\spcend ,
\end{equation}
where 
\begin{displaymath}
A_\pm = \scalea^2\scaleb^2 \pm \scalea^2 \sin^2\sab \pm \scaleb^2 \cos^2\sab
\spcend .
\end{displaymath}
For the case where $\scalea=\scaleb$ the anisotropic cocycle reduces
to the usual isotropic cocycle
\begin{displaymath}
\cocycle(\scalea, \scalea, \sas) =
\frac{4 \scalea^2}
{ [(\scalea^2-1)\cos{\saa} + \scalea^2+1]^2 }
\spcend .
\end{displaymath}
Although the ability to perform anisotropic dilations is of practical use, we do not achieve a wavelet basis in this setting. 
In the anisotropic setting a bounded admissibility integral cannot be determined (even in the plane), thus the synthesis of a signal from its coefficients cannot be performed. 
This results from there being no direct means of evaluating the proper measure in the absence of a group structure.
The projection of a signal onto basis functions undergoing anisotropic dilations may be performed in an analogous manner to the following discussion of the wavelet transform.  However, since these basis functions are not wavelets we restrict the following discussion to isotropic dilations.

%------------------------------------------------------------------------------
\subsection{Wavelet transform}
\label{sec:wavelet_transform}

A wavelet basis on the unit sphere may now be constructed from rotations and
isotropic dilations (where $\scalea=\scaleb$) of a mother spherical wavelet $\wav \in L^2(\sphere,\dmu)$. 
The corresponding wavelet family 
\mbox{$\{ \wav_{\scalea,\rho} \equiv \rot(\eul) \dil(\scalea,\scalea) \wav,
\; {\rho \in \sothree,} \; {\scalea \in \realno_{\ast}^{+}} \}$}
provides an over-complete set of functions in $L^2(\sphere,\dmu)$. 
%
The \cswt\ of $\sky \in L^2(\sphere,\dmu)$ is given
by the projection onto each wavelet basis function in the usual
manner,
\begin{equation}
\skywav(\scalea, \eul) 
\equiv
\int_{\sphere}
\dmu \:
\wav_{\scalea,\eul}^\cconj(\sa) \:
\sky(\sa)
=
\left< \wav_{\scalea,\eul} \mid \sky \right>
\spcend ,
\label{eqn:cswt}
\end{equation}
where the \cconj\ denotes complex conjugation.

The transform is general in the sense that all orientations in the
rotation group \sothree\ are considered, thus directional structure is
naturally incorporated.  It is important to note, however, that only
\emph{local} directions make any sense on \sphere.  There is no global
way of defining directions on the sphere\footnote{There is no
differentiable vector field of constant norm on the sphere and hence
no global way of defining directions.} -- there will always be some
singular point where the definition fails.  

The synthesis of a signal on the unit sphere from its wavelet coefficients
is given by
% \ifcol
%   \begin{eqnarray}
%   \sky(\sa) = &
%   \displaystyle
%   \!\!\!\!\!\!  \!\!\!\!\!\! \!\!\!
%   \int_{\sothree} \dx\eul
%   \int_0^\infty \frac{\dx\scalea}{\scalea^3} \nonumber \\
%   &\times \:
%   \skywav(\scalea, \eul) \:
% 	 [\rot(\eul) L_\wav \wav_\scalea](\sa)
% 	 \spcend ,
%   \end{eqnarray}
% \else
  \begin{equation}
  \sky(\sa) = 
  \int_{\sothree} \dx\eul
  \int_0^\infty \frac{\dx\scalea}{\scalea^3} \:\:
  \skywav(\scalea, \eul) \:
  [\rot(\eul) L_\wav \wav_\scalea](\sa)
  \spcend ,
  \end{equation}
% \fi
where $\dx\eul=\sin\eulb \dx\eula\dx\eulb\dx\eulc$.
The $L_\wav$ operator in $L^2(\sphere,\dmu)$ is defined by the
action
\begin{equation}
(\shcoeff{L_\wav f})_\elm = \shcoeff{f}_\elm / \admissc_\wav^\el 
\end{equation}
on the spherical harmonic coefficients of functions \mbox{$f \in L^2(\sphere,\dmu)$}, 
where  $\admissc_\wav^\el$ is defined below.
The hat denotes the spherical harmonic coefficients 
\begin{displaymath}
\shcoeff{f}_\elm 
= 
\int_{\sphere}
\dmu \:
\sh_\elm^\cconj(\sa) \:
f(\sa)
= \left< \sh_\elm \mid f \right> 
\end{displaymath}
of the decomposition
\begin{displaymath}
f(\sa) = \sum_{\el=0}^\infty \sum_{\m=-\el}^\el \shcoeff{f}_\elm \sh_\elm(\sa)
\spcend .
\end{displaymath}
We adopt the Condon-Shortley phase convention where the normalised 
spherical harmonics are defined by
\begin{displaymath}
\sh_{\elm}(\sa) = (-1)^\m \sqrt{\frac{2\el+1}{4\pi} 
\frac{(\el-\m)!}{(\el+\m)!}} \: 
\aleg{\el}{\m}{\cos\saa} \:
{ e}^{\img \m \sab}
\spcend ,
\end{displaymath}
where $\aleg{\el}{\m}{x}$ are the associated Legendre functions.  Using this normalisation the orthogonality of the spherical harmonic functions is given by
\begin{equation}
\label{eqn:shortho}
\int_\sphere
\dmu \:
\sh_{\elm}(\sa)
\sh_{\el\p \m\p}^\cconj(\sa)
=
\delta_{\el\el\p} \delta_{\m\m\p}
\spcend ,
\end{equation}
where $\delta_{ij}$ is Kronecker delta function.
In order to ensure the perfect reconstruction of a signal synthesised from its wavelet coefficients, one requires
the admissibility condition 
  \begin{equation}
  \label{eqn:admiss_full}
  0 <
  \admissc_\wav^\el \equiv
  \frac{8\pi^2}{2\el+1}
  \sum_{\m=-\el}^\el 
  \int_0^\infty
  \frac{\dx\scalea}{\scalea^3}
  \mid (\shcoeff{\wav_\scalea})_\elm \mid^2
  < \infty
  \end{equation}
to hold for all $\el \in \natno$. 
A proof of the admissibility criterion is given by \cite{wiaux:2005}.
%  The extension of the admissibility condition to
% anisotropic dilations follows directly from the proof given for the
% isotropic case by \cite{wiaux:2005} (appendix A).  
Practically, it is
difficult to apply \eqn{\ref{eqn:admiss_full}} directly, thus a
necessary (and almost sufficient) condition for admissibility is the
zero-mean \mbox{condition \cite{antoine:1999}}
\begin{equation}
\admissc_\wav \equiv \int_\sphere  \dmu \: \frac{\wav(\sa)}{1+\cos\saa}=0
\spcend .
\end{equation}
In the Euclidean limit this condition naturally reduces to the
necessary zero-mean condition for Euclidean wavelets
\cite{wiaux:2005}.


%------------------------------------------------------------------------------
\subsection{Correspondence principle and spherical wavelets}
\label{sec:wavelets}

The correspondence principle between spherical and Euclidean wavelets
states that the inverse stereographic projection of an
\emph{admissible} wavelet on the plane yields an 
\emph{admissible} wavelet on the unit sphere.  This result is proved by
\cite{wiaux:2005}.
% (appendix B); the proof remains essentially
% unchanged for our extension to anisotropic dilations.  
%
Hence, mother spherical wavelets may be constructed from the projection
of mother Euclidean wavelets on the plane: 
\begin{equation}
\label{eqn:wav_proj}
\wav(\sa) = [\spo^{-1}\wav_{\realno^2}](\sa)
\spcend ,
\end{equation}
where $\wav_{\realno^2} \in L^2(\realno^2,\dx^2\xvect)$ is an admissible wavelet in the plane.
Directional spherical wavelets may be naturally constructed in this 
setting -- they are simply the projection of directional Euclidean planar
wavelets on to the sphere.

We give examples of three spherical wavelets: the spherical Mexican
hat wavelet (\smhw); the spherical butterfly wavelet (\sbw); and
the spherical real Morlet wavelet (\smw).  These spherical wavelets are
illustrated in \fig{\ref{fig:mother_wavelets}}.
%
Each spherical wavelet is constructed by the stereographic projection
of the corresponding Euclidean wavelet onto the sphere, where the
Euclidean planar wavelets are defined by 
\begin{displaymath}
\wav_{\realno^2}^{\rm \smhw}(r,\sab) = 
\frac{1}{2}
(2 - r^2) \, e^{-r^2/2}
\spcend ,
\end{displaymath}
\begin{displaymath}
\wav_{\realno^2}^{\rm \sbw}(x,y) = 
x \, e^{-(x^2+y^2)/2}
\end{displaymath}
and
\begin{displaymath}
\wav_{\realno^2}^{\rm \smw}(\xvect; \kvect) = 
{\rm Re} \left(
e^{ \img \kvect \cdot \xvect / \sqrt{2}} \:
e^{ - \| \xvect \|^2 /2}
\right)
\end{displaymath}
respectively, where $\kvect$ is the wave vector of the \smw. 
%($\kvect$ specifies the direction orthogonal to the ripples).
The \smhw\ is proportional to the Laplacian of a Gaussian,
whereas the \sbw\ is proportional to the first partial derivative of a
Gaussian in the $x$-direction.  The \smw\ is a Gaussian modulated
sinusoid, or Gabor wavelet.

\begin{figure*}
\begin{minipage}{\textwidth}
\centering
\ifbw
  \mbox{
  \subfigure[Spherical Mexican hat wavelet (\smhw)]{\includegraphics[width=\mapwidth]{figures_bw/mexh_d0-2_moll_bw}} \quad
  \subfigure[Spherical butterfly wavelet (\sbw)]{\includegraphics[width=\mapwidth]{figures_bw/bfly_d0-2_moll_bw}} \quad
  \subfigure[Spherical real Morlet wavelet (\smw)]{\includegraphics[width=\mapwidth]{figures_bw/mrlt_d0-2_moll_bw}}
  }
\else
  \mbox{
  \subfigure[Spherical Mexican hat wavelet (\smhw)]{\includegraphics[width=\mapwidth]{figures/mexh_d0-2_moll}} \quad
  \subfigure[Spherical butterfly wavelet (\sbw)]{\includegraphics[width=\mapwidth]{figures/bfly_d0-2_moll}} \quad
  \subfigure[Spherical real Morlet wavelet (\smw)]{\includegraphics[width=\mapwidth]{figures/mrlt_d0-2_moll}}
  }
\fi 
\caption{Spherical wavelets at scale $\scalea=\scaleb=0.2$.  Wavelet
  maps are displayed in the Mollweide projection, where the wavelets
  have been rotated down from the north pole for ease of observation.
  The \smw\ is plotted for wave vector $\kvect=(10,0)^T$.}
\label{fig:mother_wavelets}
\end{minipage}
\end{figure*}

A full directional wavelet analysis on the unit sphere for large data sets
has previously been prohibited by the computational infeasibility of any
implementation. 
The computational burden of computing many orientations may be
reduced by using steerable wavelets, for which any continuous
orientation can be computed from a small number of basis
orientations \cite{wiaux:2005}.  This is achieved since steerable wavelets have a limited azimuthal band limit and may thus be represented as a finite sum of trigonometric exponentials \cite{wiaux:2005}.
However, in this case one must still
compute the initial transform for more than one orientation,
so although the computational burden is reduced, it is still
significant.  Moreover, we also require a fast approach for general 
non-steerable, directional wavelets.
We address this problem in the following section by
presenting a fast algorithm to perform the directional \cswt.


%------------------------------------------------------------------------------
% Algorithms
%------------------------------------------------------------------------------
\section{Algorithms}
\label{sec:algorithms}

A range of algorithms of varying computational efficiency and numerical
accuracy are presented to perform the \cswt\ described in
\sectn{\ref{sec:cswt}}.  
We implement these algorithms in Fortran 90 and subsequently compare
computational complexity and typical execution time. 
The synthesis of a signal from its wavelet coefficients is not
considered any further.
%Computational complexity and typical execution times for each
%algorithm are subsequently compared.  
Without loss of generality we consider only a single dilation
(\ie\ fixed \scalea\ and \scaleb).
 
%------------------------------------------------------------------------------
\subsection{Tessellation schemes}
\label{sec:pixelisations}

It is necessary to discretise both the spherical coordinates of a
function defined on the unit sphere and also the Euler angle representation
of the \sothree\ rotation group.
% Alternative first sentence:
%Both the Euler angle representation of the SO(3) rotation group and
%also the spherical coordinates of a function defined on the unit sphere
%must be discretised. 
%
The fast algorithms we present are performed in harmonic space and
hence are tessellation independent, provided an appropriate spherical
harmonic transform is defined for the tessellation.  However, the
semi-fast algorithm is restricted to an equi-angular
tessellation of the sphere.  
%
The various tessellation schemes adopted are defined below.
%We first define the various tessellation schemes adopted.

The equi-angular tessellation (also known as the equidistant
cylindrical projection (ECP)) of the spherical coordinates is defined
by 
%\begin{equation}
$
\gridecp = \{ \saa_\ith=\frac{\pi\ith}{\nth}, \,
\sab_\ith=\frac{2\pi\iph}{\nph} :
0 \leq \ith \leq \nth - 1, \,
0 \leq \iph \leq \nph - 1
 \}
$. %\end{equation}
Let $\npix=\nth\nph$ denote the number of pixels in the tessellation.

We also consider the \healpix\ tessellation scheme since it is
commonly used for astrophysical data-sets of the CMB.  Pixels in the
\healpix\ scheme are of equal area and are located on rings of
constant latitude (the latter feature enables fast spherical harmonic 
transforms on the pixelised grid).  We refer the reader to
\cite{gorski:2005} for details of the \healpix\ scheme and here just
define the \healpix\ grid in terms of pixel indices:
%\begin{equation}
$
\gridhpix = \{ 
(\sas)_\ipix = (\saa_\ipix, \sab_\ipix) :
0 \leq \ipix \leq \npix-1
\}
$.
%\spcend ,\end{equation}
The \healpix\ resolution is parameterised by \nside, where
$\npix=12\nside^2$.
It should be noted that an exact quadrature formula does not exist for the \healpix\ tessellation, thus spherical harmonic transforms are necessarily approximate.  This is not the case for the ECP or other practical tessellations (\eg\ IGLOO and GLESP) where exact quadrature may be performed.

The Euler angle domain of the spherical wavelet coefficients is in
general arbitrary, however we use the equi-angular discretisation
defined by
%\begin{eqnarray}
$
\grideula = \{ 
\eula_\ia=\frac{2\pi\ia}{\na}, \,
\eulb_\ib=\frac{\pi\ib}{\nb}, \,
\eulc_\ig=\frac{2\pi\ig}{\ngm} :
0 \leq \ia \leq \na - 1, \,
0 \leq \ib \leq \nb - 1, \,
0 \leq \ig \leq \ngm - 1 
\}
$.
%\spcend .\end{eqnarray}
Our fast algorithm, however, requires (for convenience) the
tessellation 
%\begin{eqnarray}
$
\grideulb = \{ 
\eula_\ia=\frac{2\pi\ia}{\na}, \,
\eulb_\ib=\frac{2\pi\ib}{\nb}, \,
\eulc_\ig=\frac{2\pi\ig}{\ngm} :
0 \leq \ia \leq \na - 1,  \,
0 \leq \ib \leq 2\nb - 1, \,
0 \leq \ig \leq \ngm - 1 
\}
$,
%\spcend ,\end{eqnarray}
where the \eulb\ sampling is repeated.
Evaluating \eulb\ over the range $0$ to $2\pi$ is redundant, covering
the \sothree\ manifold exactly twice.  Nonetheless, the use of our fast
algorithm requires this range.  
Such an approach is not uncommon, as \cite{driscoll:1994,healy:2003}
also oversample a function of the sphere in the \saa\ direction in
order to develop fast spherical harmonic transforms on equi-angular grids. 
The overhead associated with our inefficient discretisation is more than
offset by the fast algorithm it affords, as described in
\sectn{\ref{sec:alg_fast}}.   


%------------------------------------------------------------------------------
\subsection{Direct case}

The \cswt\ defined by \eqn{\ref{eqn:cswt}} may be implemented directly
by applying an appropriate quadrature rule.
Using index subscripts to denote sampled signals, the direct \cswt\
implementation is given by 
\ifcol
  \begin{eqnarray}
  \label{eqn:alg_direct}
  \lefteqn{(\skywav)_{\ia,\ib,\ig} =}  \nonumber \\
  &&\displaystyle \sum_{\ipix=0}^{\npix-1}
  \left[\rot\left(\frac{2\pi\ia}{\na},\frac{\pi\ib}{\nb},\frac{2\pi\ig}{\ngm}\right)%
  \wav\right]_{\ipix}^\cconj
  \sky_{\ipix} \:
  \weight_\ipix
  \spcend ,
  \end{eqnarray}
\else
  \begin{equation}
  \label{eqn:alg_direct}
  (\skywav)_{\ia,\ib,\ig} = \sum_{\ipix=0}^{\npix-1}
  \left[\rot\left(\frac{2\pi\ia}{\na},\frac{\pi\ib}{\nb},\frac{2\pi\ig}{\ngm}\right)%
  \wav\right]_{\ipix}^\cconj 
  \sky_{\ipix} \:
  \weight_\ipix
  \spcend ,
  \end{equation}
\fi
where the pixel sum is over all the pixels of any chosen tessellation.
%may be over either the ECP \gridecp\ grid or the \healpix\ \gridhpix\ grid.  
The weights for the ECP \gridecp\ grid are given by 
$\weight_\ipix=\weight_{\ith}=\frac{2\pi^2 \sin{\saa_\ith}}{\nth\nph}$,
whereas the equal pixel areas of the \healpix\ \gridhpix\ grid ensure
the pixel weights, given by $\weight_\ipix=\frac{4\pi}{\npix}$, are
independent of position.  

Discretisation techniques other than the plain Riemann sum used in
\eqn{\ref{eqn:alg_direct}} would be beneficial only if additional
regularity conditions are imposed on the signal \sky\
\cite{antoine:2002}.  It is also possible to choose other weights to
achieve a better approximation of the integral.  An example of a
different equi-angular discretisation and a different choice for the
weights is given by the sampling theorem for band-limited functions on
the sphere developed by \cite{healy:2003}. 

Evaluation of \eqn{\ref{eqn:alg_direct}} requires the computation of a
\mbox{2-dimensional} summation, evaluated over a \mbox{3-dimensional} grid.  We
assume the number of samples for each discretised angle, except
$\gamma$, is of the same order \n.  Typically the number of samples
in the $\gamma$ direction is of a much lower order, so we treat this
term separately.  The complexity of the direct algorithm is 
$\order(\ngm\n^4)$.

% All rotations are performed directly in real space in this algorithm.
% Arbitrary rotations cannot be performed on a discretised sphere in
% real space without introducing numerical errors,%
% \footnote{A finite point set on the sphere that is invariant under
%   rotations does not exist, thus arbitrary rotations cannot be
%   performed exactly in real space on a tessellated sphere.}
% hence the accuracy of the algorithm is compromised for low
% resolutions. 


%------------------------------------------------------------------------------
\subsection{Semi-fast case}

We rederive here the semi-fast implementation of the \cswt\ described
by \cite{antoine:2002} and implemented using the \yawtb\ Matlab wavelet toolbox and SpharmonicKit.
This algorithm involves performing a separation of variables so that
one rotation may be performed in Fourier space.  The algorithm is
restricted to the equi-angular grid \gridecp\ (in essence pixels must
only be defined on equal latitude rings, however some form of
interpolation and down-sampling is then required to extract samples
for equal longitudes). 

Firstly, the \eula\ rotation is represented by shifting the corresponding
wavelet samples to give 
\ifcol
  \begin{eqnarray*}
  \lefteqn{(\skywav)_{\ia,\ib,\ig} =} \nonumber \\
  &&
  \displaystyle
  \sum_{\ith=0}^{\nth-1}
  \sum_{\iph=0}^{\nph-1}
  \left[\rot\left(0,\frac{\pi\ib}{\nb},\frac{2\pi\ig}{\ngm}\right)
  \wav\right]_{\ith,\iph-\ia}^\cconj \!\!
  \sky_{\ith,\iph} 
  \weight_\ith
  \spcend ,
  \end{eqnarray*}
\else
  \begin{displaymath}
  (\skywav)_{\ia,\ib,\ig} = 
  \sum_{\ith=0}^{\nth-1}
  \sum_{\iph=0}^{\nph-1}
  \left[\rot\left(0,\frac{\pi\ib}{\nb},\frac{2\pi\ig}{\ngm}\right)
  \wav\right]_{\ith,\iph-\ia}^\cconj \:
  \sky_{\ith,\iph} \:
  \weight_\ith
  \spcend ,
  \end{displaymath}
\fi
where the index \iph\ is extended periodically with period \nph.  The
discrete space convolution theorem may then be applied to represent
the inner summation as the inverse discrete Fourier transform (\dft) of
the product of the wavelet and signal \dft\ samples (note that only a
\mbox{1-dimensional} \dft\ is performed in the azimuthal direction):
\ifcol
  \begin{eqnarray}
  \label{eqn:alg_semi}
  \lefteqn{(\skywav)_{\ia,\ib,\ig} =} \nonumber \\ 
  &&
  \sum_{\ith=0}^{\nth-1}
  \left\{
  \frac{1}{\nph} \right.
  \sum_{\ik=0}^{\nph-1}
  \mathcal{F}^\cconj
  \left[\rot\left(0,\frac{\pi\ib}{\nb},\frac{2\pi\ig}{\ngm}\right)
  \wav\right]_{\ith,\ik} \nonumber \\
  &&\left. \times \:
  {\mathcal{F}(\sky)}_{\ith,\ik} \:
  e^\frac{\img2\pi\ik\iph}{\nph}
  \right\}
  \weight_\ith
  \spcend ,
  \end{eqnarray}
\else
  \begin{equation}
  \label{eqn:alg_semi}
  (\skywav)_{\ia,\ib,\ig} = 
  \sum_{\ith=0}^{\nth-1}
  \left\{
  \frac{1}{\nph}
  \sum_{\ik=0}^{\nph-1}
  \mathcal{F}^\cconj
  \left[\rot\left(0,\frac{\pi\ib}{\nb},\frac{2\pi\ig}{\ngm}\right)
  \wav\right]_{\ith,\ik} \:
  {\mathcal{F}(\sky)}_{\ith,\ik} \:
  e^\frac{\img2\pi\ik\iph}{\nph}
  \right\}
  \weight_\ith
  \spcend ,
  \end{equation}
\fi
where $\mathcal{F}(\cdot)_{n,k}$ denotes \mbox{1-dimensional} \dft\ coefficients.
A fast Fourier transform (\fft) may then be applied to evaluate
simultaneously all of the $\ia$ terms of the expression enclosed in the 
curly braces in \eqn{\ref{eqn:alg_semi}}.  A final
summation (integral) over $\ith$ produces the spherical wavelet
coefficients for a given $\ib$ and $\ig$, for \emph{all} $\ia$.
%
Applying an \fft\ to evaluate simultaneously one summation rapidly,
reduces the complexity of the \cswt\ implementation to 
$\order(\ngm \n^3 \log_2 \n  )$.

% Now only two of the three rotations must be performed in real space.
% This reduces the numerical errors associated with rotating a function
% on the sphere in real space for low resolution tessellations.  To
% eliminate this problem altogether, we subsequently consider fast
% \cswt\ algorithms posed purely in harmonic space. 

% These algorithms
%are also independent of the tessellation of the sphere.
% provided a
%spherical harmonic transform algorithm is defined for the tessellation.


%------------------------------------------------------------------------------
\subsection{Fast azimuthally symmetric case}

The fast azimuthally symmetric \cswt\ algorithm is posed in harmonic space, where 
$\shcoeff{s}_\elm = \left< \sh_\elm \mid s \right>$ are the spherical harmonic
coefficients of a function $s \in L^2(\sphere,\dmu)$, as described in
\sectn{\ref{sec:wavelet_transform}}.  For the special case where the
wavelet is azimuthally symmetric (\ie\ invariant under azimuthal rotations), it is essentially only a function of \saa\ and may be represented in terms of its Legendre expansion.
% the \eulc\ angle drops out and the wavelet coefficients are defined over
% the sphere, \ie\ over only \eula\ and \eulb\ angles. 
In this case the 
harmonic representation of the wavelet coefficients is given by
the product of the signal and wavelet spherical harmonic coefficients:
\begin{equation}
\label{eqn:fast_isotropic}
(\shcoeff{\skywav})_\elm = 
\sqrt{\frac{4 \pi}{2\el+1}} \:
\shcoeff{\wav}_{\el 0}^\cconj \:
\shcoeff{\sky}_\elm
\spcend ,
\end{equation}
noting that the harmonic coefficients of an azimuthally symmetric wavelet are zero
for $\m\neq0$.  In practice, one requires that at least one of the
signals, usually the wavelet, has a finite band limit so that
negligible power is present in those coefficients above a certain
\elmax.  We then only need to consider $\el \leq \elmax$ (a detailed discussion of the determination of \elmax\ is presented in \cite{bogdanova:2004}).
Once the spherical harmonic representation of the wavelet coefficients
is calculated by \eqn{\ref{eqn:fast_isotropic}}, the inverse spherical
harmonic transform is applied to compute the wavelet coefficients in
the Euler domain. 
%The algorithm simply proceeds as follows: (i)
%compute the spherical harmonic coefficients of the wavelet and signal
%using a spherical harmonic transform;
%(ii) compute the spherical harmonic coefficients of the wavelet
%coefficients using \eqn{\ref{eqn:fast_isotropic}}; (iii) compute the
%wavelet coefficients from their spherical harmonic representation
%using an inverse spherical harmonic transform.
%
The complexity of the fast isotropic \cswt\ algorithm is dominated by
the spherical harmonic transforms.  
For a tessellation containing pixels on rings of constant latitude, a fast spherical harmonic transform may be performed (see \eg \cite{mohlenkamp:1997,mohlenkamp:1999,gorski:2005}).  This reduces the complexity of the spherical harmonic transform from $\order(\n^4)$ to $\order(\n^3)=\order(\npix^{3/2})$. 
For certain tessellation schemes fast spherical harmonic transforms of lower complexity are also available, however these are related directly to the tessellation (\eg\ \cite{mohlenkamp:1997,mohlenkamp:1999,driscoll:1994,healy:2003}).  In particular, the algorithm developed by \cite{driscoll:1994,healy:2003} for the ECP tessellation scales as $\order(\n^2\log\n)$.

The fast azimuthally symmetric \cswt\ algorithm is posed purely in harmonic space and consequently 
% any errors introduced by performing rotations in real space are eliminated, and moreover, 
the algorithm is tessellation independent.  However, we are
restricted to azimuthally symmetric wavelets and lose the ability to
perform the directional analysis inherent in the wavelet transform
construction.  %This is a severe limitation.


%------------------------------------------------------------------------------
\subsection{Fast directional case}
\label{sec:alg_fast}

We present the most general fast directional \cswt\ algorithm for
non-azimuthally symmetric wavelets, \ie\ steerable and directional wavelets, in this section.
Again, the algorithm is posed purely in harmonic space and so 
is tessellation independent.  We do, however, use the
equi-angular \grideulb\ discretisation of the wavelet coefficient
domain, although other discretisations may be used if \fft s are also 
defined on these grids.
%
The \cswt\ at a particular scale (\ie\ a particular \scalea\ and
\scaleb) is essentially a spherical 
{con\-volu\-tion}, hence we apply the fast spherical
convolution algorithm proposed by \cite{wandelt:2001} to evaluate
the wavelet transform.   
The algorithm proceeds by factoring the rotation into two separate
rotations, each of which involves only a constant polar rotation component. 
%that is independent of all Euler angles. 
Azimuthal rotations may then be performed in harmonic space at far
less computation expense than polar rotations. 
We subsequently rederive the fast spherical convolution algorithm
developed by \cite{wandelt:2001}, as applied to our application of
evaluating the \cswt.
The harmonic representation of the \cswt\ is first presented, followed
by the discretisation and fast implementation.

\subsubsection{Harmonic formulation}

Substituting the spherical harmonic expansions of the wavelet and signal
into the wavelet transform defined by \eqn{\ref{eqn:cswt}}, and noting the
orthogonality of the spherical harmonics described by \eqn{\ref{eqn:shortho}},
yields the harmonic representation
\ifcol
  \begin{eqnarray}
  \label{eqn:harmonic:1}
  \lefteqn{\skywav(\eula, \eulb, \eulc) =} \nonumber \\
  &&
  \sum_{\el =0}^{\elmax} \:
  \sum_{\m=-\el}^{\el} \:
  \sum_{\m\p=-\el}^{\el}
  \left [
  {\dmatbig_{\m\m\p}^{\el }(\euls) \,
  \shcoeff{\wav}_{\el \m\p}} \right ]^\cconj 
  \shcoeff{\sky}_{\el \m}
  \spcend .
  \end{eqnarray}
\else
  \begin{equation}
  \label{eqn:harmonic:1}
  \skywav(\eula, \eulb, \eulc) =
  \sum_{\el =0}^{\elmax} \:
  \sum_{\m=-\el}^{\el} \:
  \sum_{\m\p=-\el}^{\el}
  \left [
  {\dmatbig_{\m\m\p}^{\el }(\euls) \,
  \shcoeff{\wav}_{\el \m\p}} \right ]^\cconj \:
  \shcoeff{\sky}_{\el \m}
  \spcend .
  \end{equation}
\fi
Again, we assume negligible power above \elmax\ in at least one of the
signals, usually the wavelet, so that the outer summation is truncated
to \elmax. 
%
The additional summation over $\m\p$ and the $\dmatbig_{mm\p}^\el$
Wigner rotation matrices that are introduced
characterise the rotation of a spherical harmonic, noting that a rotated
spherical harmonic may be represented by a sum of 
harmonics of the same $\el$ \cite{brink:1993,ritchie:1999}:
\begin{displaymath}
\left[\rot(\euls)\sh_{\elm}\right](\sa) = 
\sum_{\m\p=-\el}^{\el} 
\dmatbig_{\m\p\m}^{\el}(\euls) \: \sh_{\el \m\p}(\sa)
\spcend .
\end{displaymath}
The Wigner rotation matrices may be decomposed as \cite{brink:1993,ritchie:1999}
\begin{equation}
\label{eqn:d_decomp}
\dmatbig_{\m\m\p}^{\el}(\euls)
= e^{-\img \m\eula} \:
\dmatsmall_{\m\m\p}^\el(\eulb) \:
e^{-\img \m\p\eulc}
\spcend ,
\end{equation}
where the real polar \dmatsmall-matrix is defined by \cite{brink:1993}
\ifcol
  \begin{eqnarray*}
  \lefteqn{\dmatsmall^\el_{\m \m\p}(\eulb)  = 
  \sum_{t = \max(0, \m - \m\p)}^{\min(\el + \m, \el - \m\p)} (-1)^t} \nonumber\\
  && \times \frac{
  \left[(\el+\m)! \, (\el-\m)! \, (\el+\m\p)! \, (\el-\m\p)! \, \right]^{1/2}}
  {(\el+\m-t)! \, (\el-\m\p-t)! \, (t+\m\p-\m)! \, t!} \nonumber \\
  && \times \left[\cos\!\left(\frac{\eulb}{2} \right)\right]^{2\el+\m-\m\p-2t}
  \left[\sin \! \left(\frac{\eulb}{2} \right)\right]^{\m\p-\m+2t}
  \spcend ,
  \end{eqnarray*}
\else
  \begin{eqnarray*}
  \dmatsmall^\el_{\m \m\p}(\eulb) & = & 
  \sum_{t = \max(0, \m - \m\p)}^{\min(\el + \m, \el - \m\p)} (-1)^t
  \frac{
  \left[(\el+\m)! \, (\el-\m)! \, (\el+\m\p)! \, (\el-\m\p)! \, \right]^{1/2}}
  {(\el+\m-t)! \, (\el-\m\p-t)! \, (t+\m\p-\m)! \, t!} \nonumber \\
  && \times \left[\cos\!\left(\frac{\eulb}{2} \right)\right]^{2\el+\m-\m\p-2t}
  \left[\sin \! \left(\frac{\eulb}{2} \right)\right]^{\m\p-\m+2t}
  \spcend ,
  \end{eqnarray*}
\fi
and the sum over $t$ is defined so that the arguments of the factorials
are non-negative.  Recursion formulae are available to compute rapidly the Wigner rotation matrices in the basis of either complex \cite{risbo:1996,choi:1999} or real \cite{ivanic:1996,blanco:1997} spherical harmonics.  We employ the recursion formulae described in \cite{risbo:1996} in our implementation.
%	
The decomposition shown in \eqn{\ref{eqn:d_decomp}} is exploited by
factoring the rotation $\rot(\euls)$ into two separate rotations,
both of which only contain 
a constant $\pm \pi/2$ polar rotation:
\ifcol
  \begin{eqnarray}
  \label{eqn:rot_factor}
  \rot(\euls)
  &=& \rot(\eula-\pi/2, \; -\pi/2, \; \eulb) \nonumber \\
  &&\times \:
  \rot(0, \; \pi/2, \; \eulc+\pi/2)
  \spcend .
  \end{eqnarray}
\else
  \begin{equation}
  \label{eqn:rot_factor}
  \rot(\euls)
  = \rot(\eula-\pi/2, \; -\pi/2, \; \eulb) \:\:
  \rot(0, \; \pi/2, \; \eulc+\pi/2)
  \spcend .
  \end{equation}
\fi
By factoring the rotation in this manner and applying the decomposition
described by \eqn{\ref{eqn:d_decomp}}, \eqn{\ref{eqn:harmonic:1}} 
can be rewritten as
\ifcol
  \begin{eqnarray}
  \label{eqn:harmonic:2}
  \lefteqn{\skywav(\euls) =} \nonumber \\
  &&
  \sum_{\el=0}^{\elmax} \:
  \sum_{\m=-\el}^{\el} \:
  \sum_{\m\p=-\el}^{\el}
  \sum_{\m\pp=-\min(\mmax,\el)}^{\min(\mmax,\el)}
  \dmatsmall_{\m\p \m}^\el(\pi/2) \:
  \nonumber \\
  && \times \:
  \dmatsmall_{\m\p \m\pp}^\el(\pi/2) \:
  \shcoeff{\wav}_{\el \m\pp}^\cconj \:
  \shcoeff{\sky}_{\el \m} \nonumber \\
  && \times \:
  e^{\img [\m(\eula-\pi/2) + \m\p\eulb + 
  \m\pp(\eulc+\pi/2) ]}
  \spcend ,
  \end{eqnarray}
\else
  \begin{eqnarray}
  \label{eqn:harmonic:2}
  \skywav(\euls) &=&
  \sum_{\el=0}^{\elmax} \:
  \sum_{\m=-\el}^{\el} \:
  \sum_{\m\p=-\el}^{\el} \:
  \sum_{\m\pp=-\min(\mmax,\el)}^{\min(\mmax,\el)}
  \dmatsmall_{\m\p \m}^\el(\pi/2) \:
  \dmatsmall_{\m\p \m\pp}^\el(\pi/2)
  \nonumber \\
  & & 
  \times \:\:
  \shcoeff{\wav}_{\el \m\pp}^\cconj \:
  \shcoeff{\sky}_{\el \m} \:
  e^{\img [\m(\eula-\pi/2) + \m\p\eulb + 
  \m\pp(\eulc+\pi/2) ]}
  \spcend ,
  \end{eqnarray}
\fi
where the symmetry relationship
\mbox{$\dmatsmall_{\m \m\p}^{\el}(-\eulb)=\dmatsmall_{\m\p \m}^{\el}(\eulb)$}
has been applied.  
(A similar factoring of the rotation and harmonic space representation has been independently performed in \cite{bulow:2001}.)
%
Steerable wavelets may have a low effective band limit 
$\mmax \ll \elmax$, in which case the the inner summation in \eqn{\ref{eqn:harmonic:2}} may be truncated to \mmax.
For general directional wavelets this is not the case and one must use $\mmax = \elmax$.
% In many cases it is likely that the wavelet will have minimal 
% azimuthal structure compared to the signal under analysis, in
% which case it may also have a lower effective azimuthal band limit
% $\mmax \ll \elmax$.

Evaluating the harmonic formulation given by
\eqn{\ref{eqn:harmonic:2}} directly would be no more efficient that 
approximating the initial integral \eqn{\ref{eqn:cswt}} using
simple quadrature.  
However, 
% rotations are elegantly represented in harmonic space and the  
% approximation and interpolation required in any real space
% discretisation is avoided.  Moreover, 
\eqn{\ref{eqn:harmonic:2}}
is represented in such a way that the presence of complex exponentials
may be exploited such that \fft s may be applied to evaluate rapidly
the three summations simultaneously. 

%The harmonic formulation presented replaces the continuous integral
%of \eqn{\ref{eqn:cswt}} by finite summations, although evaluating these
%summations directly would be no more efficient that approximating the
%initial integral using simple quadrature.  
%Rotations are
%elegantly represented in harmonic space, however, and the
%approximation and interpolation required in any real space
%discretisation is avoided.  Moreover, \eqn{\ref{eqn:harmonic:2}}
%is represented in such a way that the presence of complex exponentials
%may be exploited such that \fft s may be applied to evaluate rapidly
%the three summations simultaneously. 


\subsubsection{Fast implementation}

Azimuthal rotations may be applied with far less computational expense
than polar rotations since they appear within complex exponentials
in \eqn{\ref{eqn:harmonic:2}}.  Although the \dmatsmall-matrices can be
evaluated reasonably quickly and reliably using recursion formulae, the
basis for the fast implementation is to avoid these polar rotations as
much as possible and use \fft s to evaluate rapidly all of the azimuthal
rotations simultaneously. 
This is the motivation for factoring the rotation by
\eqn{\ref{eqn:rot_factor}} so that all Euler angles occur as azimuthal
rotations.

The discretisation of each Euler angle may in general be {arbi\-trary}.
However, to exploit standard \fft\ routines uniform sampling is
adopted, \ie\ grid \grideulb\ is used
(see \sectn{\ref{sec:pixelisations}}).  As mentioned, this
discretisation is inefficient since it covers the \sothree\ manifold
exactly twice, nevertheless it enables the use of standard \fft\
routines which significantly increases the speed of the algorithm. 
%
Discretising \eqn{\ref{eqn:harmonic:2}} in this manner and
interchanging the order of summation we obtain
\ifcol
  \begin{eqnarray*}
%   \label{eqn:harmonic:3}
  \lefteqn{(\skywav)_{\ia, \ib, \ig}  =} \nonumber \\
  && \displaystyle 
  \sum_{\m=-\elmax}^{\elmax} \:
  \sum_{\m\p=-\elmax}^{\elmax} \:
  \sum_{\m\pp=-\mmax}^{\mmax} 
  \nonumber \\
  &&
  \displaystyle
  \sum_{\el=\max(\mid \m \mid, \mid \m\p \mid, \mid \m\pp \mid)}^{\elmax}
  \dmatsmall_{\m\p \m}^\el(\pi/2) \:
  \dmatsmall_{\m\p \m\pp}^\el(\pi/2) \:
  \shcoeff{\wav}_{\elm\pp}^\cconj \:
  \shcoeff{\sky}_{\elm} \nonumber \\
  \nonumber \\
  &&\times \:
  e^{\img [\m(2\pi\ia/\na - \pi/2) 
  + \m\p 2\pi\ib/\nb
  + \m\pp (2\pi\ig/\ngm + \pi/2) ]}
  \spcend .
  \end{eqnarray*}
\else
  \begin{eqnarray*}
%   \label{eqn:harmonic:3}
  (\skywav)_{\ia, \ib, \ig} & = &
  \sum_{\m=-\elmax}^{\elmax} \:
  \sum_{\m\p=-\elmax}^{\elmax} \:
  \sum_{\m\pp=-\mmax}^{\mmax} 
  \sum_{\el=\max(\mid \m \mid, \mid \m\p \mid, \mid \m\pp \mid)}^{\elmax}
  \dmatsmall_{\m\p \m}^\el(\pi/2) \:
  \dmatsmall_{\m\p \m\pp}^\el(\pi/2)
  \nonumber \\
  & & 
  \times
  \shcoeff{\wav}_{\elm\pp}^\cconj \:
  \shcoeff{\sky}_{\elm} \:
  e^{\img [\m(2\pi\ia/\na - \pi/2) 
  + \m\p 2\pi\ib/\nb
  + \m\pp (2\pi\ig/\ngm + \pi/2) ]}
  \spcend .
  \end{eqnarray*}
\fi
%
Shifting the indices yields
\ifcol
  \begin{eqnarray}
  \label{eqn:cswt_fast}
  \lefteqn{(\skywav)_{\ia, \ib, \ig} =} \nonumber \\
  &&\displaystyle \sum_{\m=0}^{2\elmax} \:
  \sum_{\m\p=0}^{2\elmax} \:
  \sum_{\m\pp=0}^{2\mmax}
  e^{\img 2 \pi
  ( \ia\m/\na + \ib\m\p/\nb + \ig\m\pp/\ngm )} \nonumber \\
  &&\times \:
  \cswtfftterm_{\m,\m\p,\m\pp}
  \spcend ,
  \end{eqnarray}
\else
  \begin{equation}
  \label{eqn:cswt_fast}
  (\skywav)_{\ia, \ib, \ig} =
  \sum_{\m=0}^{2\elmax} \:
  \sum_{\m\p=0}^{2\elmax} \:
  \sum_{\m\pp=0}^{2\mmax}
  e^{\img 2 \pi
  ( \ia\m/\na + \ib\m\p/\nb + \ig\m\pp/\ngm )}
  \cswtfftterm_{\m,\m\p,\m\pp}
  \spcend ,
  \end{equation}
\fi
where
\ifcol
  \begin{eqnarray}
  \lefteqn{\cswtfftterm_{\m,\m\p,\m\pp} =
  e^{\img(\m\pp-\m)\pi/2}} \nonumber \\
  &&\times
  \displaystyle
  \sum_{\el=\max(\mid \m \mid, \mid \m\p \mid, \mid \m\pp \mid)}^{\elmax}
  \dmatsmall_{\m\p \m}^\el(\pi/2) \:
  \dmatsmall_{\m\p \m\pp}^\el(\pi/2) \nonumber \\
  &&\times \:
  \shcoeff{\wav}_{\el \m\pp}^\cconj \:
  \shcoeff{\sky}_{\el \m}
  \label{eqn:cswt_fast_term}
  \end{eqnarray}
\else
  \begin{equation}
  \cswtfftterm_{\m,\m\p,\m\pp} =
  e^{\img(\m\pp-\m)\pi/2}
  \sum_{\el=\max(\mid \m \mid, \mid \m\p \mid, \mid \m\pp \mid)}^{\elmax}
  \dmatsmall_{\m\p \m}^\el(\pi/2) \:
  \dmatsmall_{\m\p \m\pp}^\el(\pi/2) \:
  \shcoeff{\wav}_{\el \m\pp}^\cconj \:
  \shcoeff{\sky}_{\el \m}
  \label{eqn:cswt_fast_term}
  \end{equation}
\fi
is extended periodically.  
Note that the phase shift introduced by shifting the indices of the
summation in \eqn{\ref{eqn:cswt_fast}}, shifts the
$\cswtfftterm_{\m,\m\p,\m\pp}$ indices back.  
Making the associations \mbox{$\na=2\elmax+1$}, $\nb=2\elmax+1$ and
$\ngm=2\mmax+1$, one notices that \eqn{\ref{eqn:cswt_fast}} is 
the unnormalised \mbox{3-dimensional} inverse \dft\ of
\eqn{\ref{eqn:cswt_fast_term}}. 
Nyquist sampling in $(\euls)$ is inherently ensured by the
associations made with \elmax\ and \mmax.

The \cswt\ may be performed rapidly in spherical harmonic
space by constructing the \cswtfftterm-matrix of
\eqn{\ref{eqn:cswt_fast_term}} from spherical harmonic coefficients and
precomputed \mbox{\dmatsmall-matrices},
followed by the application of an \fft\ to 
evaluate rapidly all three Euler angles of the discretised \cswt\
simultaneously.  The complexity of the algorithm is dominated by the 
computation of the \cswtfftterm-matrix.  This involves performing a
\mbox{1-dimensional} summation over a \mbox{3-dimensional} grid, hence the
algorithm is of order $\order(\ngm\n^3)$.

%\subsubsection{Exploiting conjugate symmetry}
Additional benefits may be realised for real signals by exploiting the
resulting conjugate symmetry relationship.  
Memory and computational requirements may be reduced by a
further factor of two for real signals by noting the conjugate
symmetry relationship
$\cswtfftterm_{-m,-m\p,-m\pp}=\cswtfftterm_{m,m\p,m\pp}^\cconj$.
In our implementation we use the 
complex-to-real \fft\ routines of the FFTW\footnote{\url{http://www.fftw.org/}}
package, which are approximately twice as fast as the
equivalent complex-to-complex routines. 


%------------------------------------------------------------------------------
\subsection{Comparison}

We summarise the computational complexities of the various \cswt\
algorithms for a 
single pair of scales and single orientation in
\tbl{\ref{tbl:complexity}}.   
The complexity of each algorithm scales with the
number of dilations considered and, for those algorithms that
facilitate a directional analysis (\ie\ all but the fast azimuthally symmetric
algorithm), with the number of orientations considered.  
The most general fast directional algorithm provides a
saving of $\order(\n)$ over the direct case, where the number of
pixels on the sphere and the harmonic band limit are related to \n\ by 
$\order(\npix)=\order(\elmax^2)=\order(\n^2)$.

\begin{table}
%\ifcol
  \caption{Algorithm complexity for one scale and one orientation.}
%\else
%  \caption{Complexity of each algorithm for one scale and one
%    orientation.  The complexity of each algorithm simply scales with the
%    number of dilations considered and, for those algorithms that
%    facilitate a directional analysis, with the number of orientations
%    considered.  The number of pixels on the sphere and the harmonic
%    band limit is related to \n\ by
%    $\order(\npix)=\order(\elmax^2)=\order(\n^2)$.}
%\fi
\label{tbl:complexity}
\begin{center}
\begin{tabular}{ll} \cline{1-2}
Algorithm & Complexity \\ \cline{1-2}
Direct           & $\order(\n^4)$ \\
Semi-fast        & $\order(\n^3\log_2\n)$ \\
Fast azimuthally symmetric   & $\order(\n^3)$ \\
Fast directional & $\order(\n^3)$ \\ \cline{1-2}
\end{tabular}
\end{center}
\end{table}

We implement all algorithms in Fortran 90, adopting the \healpix\ tessellation of the sphere
% and being free to choose any tessellation scheme we adopt the
% \healpix\ tessellation of the sphere #
(which, incidentally, 
is the tessellation scheme of the WMAP CMB data \cite{bennett:2003}).
%using the \healpix\ tessellation of the sphere.  
Typical execution times for the 
algorithms are tabulated in \tbl{\ref{tbl:exec_time}} for a range of 
resolutions from 110 down to 3.4 arcminutes.
The improvements provided by the fast
algorithms are apparent.  Indeed, it is not feasible to run the direct
algorithm on data-sets with a resolution much greater than
$\npix \simeq 5\times 10^5$.
For data-sets of practical size, such as the WMAP 
($\npix \simeq 3\times10^6$) and
forthcoming Planck ($\npix \simeq 50\times10^6$) CMB data, the fast
implementations of the \cswt\ are essential.

The semi-fast algorithm is also implemented using the \healpix\
tessellation.  However, to perform the outer summation (integration)
continual interpolation followed by down-sampling is required on each
iso-latitude ring to essentially resample the data on an ECP tessellation. 
This increases the execution 
time of the implementation of the semi-fast algorithm on the \healpix\
grid to an extent that the semi-fast algorithm provides little
advantage over the direct algorithm.  To appreciate the advantages of
the semi-fast approach it must be implemented on an ECP tessellation,
hence we do not tabulate the execution times for our implementation of
this algorithm on the \healpix\ grid as it provides an unfair
comparison. 

It is also important to note that although complexity scales with the
number of dilations and orientations considered, 
execution time does not for the fast algorithms.
Execution time does scale in this manner for
the direct algorithm.
There are a number of additional overheads
associated with the fast algorithms, such as computing spherical
harmonic coefficients and \dmatsmall-matrices.  Consequently, the fast
algorithms provide additional execution time savings that are not
directly apparent in \tbl{\ref{tbl:exec_time}}.  For example, the
execution time of the fast azimuthally symmetric and directional algorithms for 10
dilations at a resolution of $\npix \simeq 8 \times 10^5$
($\nside=256$) are 3:08.20 
and 7:06.83 (minutes:seconds) respectively, which is considerably
faster than ten times the execution time of one dilation.

\begin{table}
%\ifcol
  \caption{Typical execution time (minutes:seconds) for each algorithm
    run on an Intel P4-M 3GHz laptop with 512MB of memory.}
%\else
%  \caption{Typical execution time for each algorithm run on an Intel P4-M
%    3GHz laptop with 512MB of memory.  It is not computationally
%    feasible to run the direct algorithm on resolutions with $\nside>64$.
%    Times are give in minutes:seconds.}
%\fi
\label{tbl:exec_time}
\begin{center}
\begin{tabular}{rrrrr} \cline{1-5}
\multicolumn{2}{c}{Resolution} & \multicolumn{3}{c}{Algorithm
  execution time} \\
\multicolumn{1}{c}{\nside} & \multicolumn{1}{c}{\npix} & Direct & Fast isotropic & Fast directional \\ \cline{1-5}
  32 &     12,288 & 3:25.37         & 0:00.06  &  0:00.10 \\
  64 &     49,152 & 54:31.75        & 0:00.38  &  0:00.74 \\
 256 &    786,432 & -- \hspace{3mm} & 0:28.00  &  0:52.55 \\
 512 &  3,257,292 & -- \hspace{3mm} & 3:43.69  &  7:57.75 \\
1024 & 12,582,912 & -- \hspace{3mm} & 28:23.85 & 71:31.68 \\\cline{1-5}
\end{tabular}
\end{center}
\end{table}

%------------------------------------------------------------------------------
%\subsection{Numerical admissibility}


%------------------------------------------------------------------------------
% Application
%------------------------------------------------------------------------------
\section{Application}
\label{sec:application}

We demonstrate in this section the application of our \cswt\ implementation to binary Earth data.  In \fig{\ref{fig:world}} the Earth data and the corresponding spherical wavelet coefficients are shown.  We use the SBW defined in \sectn{\ref{sec:wavelets}} to perform the analysis.  This is a steerable wavelet \cite{wiaux:2005}, however our implementation is in general valid for any directional wavelet.  Notice how the wavelet coefficient maps corresponding to different oriented wavelets pick out corresponding oriented structure in the data.  As the dilation scale is increased, the scale of the features extracted also increases accordingly.

The ability to probe oriented structure in data defined on a sphere is of important practical use.  Certain physical processes may be localised on the sphere in scale, position and orientation (\eg\ signatures of cosmic strings in the CMB \cite{kaiser:1984} or correlations induced in the CMB by the nearby galaxy distribution \cite{crittenden:1996}).  Thus, analysing the statistical properties of spherical wavelet coefficients individually for a range of scales, positions and orientations may allow one to detect such effects with greater significance.
Indeed, using a directional spherical wavelet analysis we have made very strong detections of non-Gaussianity in the CMB \cite{mcewen:2004,mcewen:2006b} and the strongest detection made to date of the ISW effect \cite{mcewen:2006}.


\begin{figure}
\centering
\ifbw
  \subfigure[Binary Earth map]{\includegraphics[width=\coeffwidth]{figures/worldmap_nside256_flip}}
  \subfigure[Spherical butterfly wavelet coefficients for scale $\scalea=\scaleb=0.03$ and orientation $\eulc=0^\circ$]{\includegraphics[width=\coeffwidth]{figures/wc_worldmap_n256_butterfly_sky_id01_ig01_flip}}
  \subfigure[Spherical butterfly wavelet coefficients for scale $\scalea=\scaleb=0.03$ and orientation $\eulc=144^\circ$]{\includegraphics[width=\coeffwidth]{figures/wc_worldmap_n256_butterfly_sky_id01_ig03_flip}}
  \subfigure[Spherical butterfly wavelet coefficients for scale $\scalea=\scaleb=0.12$ and orientation $\eulc=0^\circ$]{\includegraphics[width=\coeffwidth]{figures/wc_worldmap_n256_butterfly_sky_id02_ig01_flip}}
  \subfigure[Spherical butterfly wavelet coefficients for scale $\scalea=\scaleb=0.12$ and orientation $\eulc=144^\circ$]{\includegraphics[width=\coeffwidth]{figures/wc_worldmap_n256_butterfly_sky_id02_ig03_flip}}
\else
  \subfigure[Binary Earth map]{\includegraphics[width=\coeffwidth]{figures/worldmap_nside256_flip}}
  \subfigure[Spherical butterfly wavelet coefficients for scale $\scalea=\scaleb=0.03$ and orientation $\eulc=0^\circ$]{\includegraphics[width=\coeffwidth]{figures/wc_worldmap_n256_butterfly_sky_id01_ig01_flip}}
  \subfigure[Spherical butterfly wavelet coefficients for scale $\scalea=\scaleb=0.03$ and orientation $\eulc=144^\circ$]{\includegraphics[width=\coeffwidth]{figures/wc_worldmap_n256_butterfly_sky_id01_ig03_flip}}
  \subfigure[Spherical butterfly wavelet coefficients for scale $\scalea=\scaleb=0.12$ and orientation $\eulc=0^\circ$]{\includegraphics[width=\coeffwidth]{figures/wc_worldmap_n256_butterfly_sky_id02_ig01_flip}}
  \subfigure[Spherical butterfly wavelet coefficients for scale $\scalea=\scaleb=0.12$ and orientation $\eulc=144^\circ$]{\includegraphics[width=\coeffwidth]{figures/wc_worldmap_n256_butterfly_sky_id02_ig03_flip}}
\fi 
\caption{Binary Earth data and corresponding spherical butterfly wavelet coefficients.  The wavelet coefficients for combinations of only two scales and two orientations are shown.  Notice how the wavelet coefficient maps corresponding to different oriented wavelets pick out corresponding oriented structure in the data.  As the dilation scale is increased, the scale of the features extracted also increases accordingly.  Directional wavelets therefore allow one not only to probe signal characteristics localised in scale and space, but also in orientation.  However, the concept of orientation on the sphere is necessarily a local one (see text \sectn{\ref{sec:wavelet_transform}} for further details).
}
\label{fig:world}
\end{figure}

%------------------------------------------------------------------------------
% Conclusions
%------------------------------------------------------------------------------
\section{Concluding remarks}
\label{sec:conclusions}

The extension of Euclidean wavelet analysis to the sphere has been
described in the framework presented by \cite{wiaux:2005}, where the
stereographic projection is used to relate the spherical and Euclidean
constructions.  
% Like the plane, the sphere is a two-dimensional
% manifold, and thus a natural candidate for anisotropic dilations
% defined in two orthogonal directions.  
We extend the concept of the spherical dilation presented
by \cite{wiaux:2005} to anisotropic dilations.  Although anisotropic dilations are of practical use, the resulting basis one projects onto does not classify as a wavelet basis since perfect reconstruction is not possible.

% noting the extended
% admissibility condition that follows directly from the proof by
% \cite{wiaux:2005} (appendix A).  We also note that the proof of the
% correspondence principle given by \cite{wiaux:2005} (appendix B),
% relating planar Euclidean and spherical wavelets, also remains
% essentially unchanged for the anisotropic dilation.

Current and forthcoming data-sets on the sphere, of the CMB for
example, are of considerable size as higher resolutions are
necessary to test new physics.  Consequently, we present fast
algorithms to implement the \cswt\ as an analysis without such
algorithms is not computationally feasible.  A range of
algorithms are described, from the direct quadrature approximation, to
the semi-fast equi-angular algorithm where one rotation is performed
in Fourier space, to the fast azimuthally symmetric and directional algorithms
posed purely in spherical harmonic space.  Posing the \cswt\ purely in
harmonic space naturally ensures the resulting algorithms are tessellation
independent.
%  and (ii) enables an analytic representation of rotations,
% thereby ensuring any errors  introduced by performing rotations in real
% space are eliminated.
The most general fast directional algorithm provides a
saving of $\order(\sqrt{\npix})=\order(\elmax)$ over the direct
implementation and 
may be performed down to a few
arcminutes even with limited computational resources.
 % (the resolution of forthcoming Planck CMB data). 
%takes just over three and a half minutes to compute
%on a 3GHz laptop with 512MB of memory for a tessellation of
%approximately three million pixels (the resolution of the current WMAP
%CMB data). 

Data is observed on a sphere in a range of applications.
% astrophysical and geophysical applications, not to mention others. 
In many of these cases the ability to perform a wavelet analysis would
be useful.
For example,
spherical wavelets may be used to probe the CMB for deviations from
the standard assumption of Gaussianity, or to search for compact
objects embedded in the CMB, such as cosmic strings, a predicted but
as yet unobserved phenomenon.  
The extension of wavelet analysis to the sphere enables one to probe
new physics in many areas, and is facilitated on large practical
data-sets by our fast directional \cswt\ algorithm.


%------------------------------------------------------------------------------
% Acknowledgements
%------------------------------------------------------------------------------
\section*{Acknowledgements}

JDM is supported by a Commonwealth (Cambridge) Scholarship from
the Association of Commonwealth Universities and the Cambridge
Commonwealth Trust.
DJM is supported by the UK Particle Physics and Astronomy Research
Council (PPARC).
The implementations described in this paper use the \healpix\
\cite{gorski:2005} and FFTW packages.  
We also acknowledge use of the \yawtb\ Matlab toolbox for the binary Earth data defined on the sphere.

fcswt_commands.tex

\newcommand{\cswtfftterm}{\ensuremath{T}}

\newcommand{\el}{\ensuremath{\ell}}
\newcommand{\m}{\ensuremath{m}}
\newcommand{\elm}{\ensuremath{{\el \m}}}
\newcommand{\mmax}{\ensuremath{\m_{\rm max}}}
\newcommand{\elmax}{\ensuremath{{\el_{\rm max}}}}
\newcommand{\dmatbig}{\ensuremath{D}}
\newcommand{\dmatsmall}{\ensuremath{d}}
\newcommand{\admissc}{\ensuremath{C}}
\newcommand{\sh}{\ensuremath{Y}}
\newcommand{\shcoeff}[1]{\ensuremath{\widehat{#1}}}

% \newcommand{\exp}[1]{\ensuremath{{{\rm e}^{#1}}}
% \newcommand{\aleg}[2]{\ensuremath{{P_{#1}^{#1}}}
\newcommand{\aleg}[3]{\ensuremath{P_{#1}^{#2}({#3})}}


% --------------------------------------
% Section abbreviations
% --------------------------------------

\newcommand{\eqn}[1]
	{(#1)}
\newcommand{\Eqn}[1]
	{(#1)}
\newcommand{\tbl}[1]
	{Table~#1}
\newcommand{\Tbl}[1]
	{Table~#1}
\newcommand{\fig}[1]
	{Fig.~#1}
\newcommand{\Fig}[1]
	{Fig.~#1}
\newcommand{\sectn}[1]
	{section~#1}
\newcommand{\Sectn}[1]
	{Section~#1}
\newcommand{\ch}[1]
	{Chapter~#1}
\newcommand{\Ch}[1]
	{Chapter~#1}
\newcommand{\appn}[1]
	{Appendix~#1}
\newcommand{\Appn}[1]
	{Appendix~#1}


\newcommand{\eg}{\mbox{\it{e.g.}}}
\newcommand{\ie}{\mbox{\it{i.e.}}}
\newcommand{\cf}{\mbox{\it{cf.}}}
\newcommand{\urlsize}{\footnotesize}


\newcommand{\fft}{{FFT}}
\newcommand{\dft}{{DFT}}

\newcommand{\isw}
	{{ISW}}
\newcommand{\iswtext}
	{integrated Sachs-Wolfe}

\newcommand{\cswt}
	{{CSWT}}
\newcommand{\cswttext}
	{continuous spherical wavelet transform}
\newcommand{\Cswttext}
	{Continuous Spherical Wavelet Transform}
\newcommand{\smhw}{{SMHW}}
\newcommand{\smw}{{SMW}}
\newcommand{\sbw}{{SBW}}

\newcommand{\yawtb}{{YAWTb}}
\newcommand{\healpix}{{HEALPix}}
\newcommand{\healpixtext}{{Hierarchical Equal Area isoLatitude Pixelisation of the sphere}}
\newcommand{\xvect}{\ensuremath{\mathbf{x}}}
\newcommand{\kvect}{\ensuremath{\mathbf{k}}}

\newcommand{\sky}{\ensuremath{s}}
%\newcommand{\skywav}{\ensuremath{(\mathcal{W}_{\wav}\sky)}}
\newcommand{\skywav}{\ensuremath{{W_\wav^\sky}}}
%	{\ensuremath{W}}

\newcommand{\spcend}{\ensuremath{\:}}

\newcommand{\img}{\ensuremath{\mathit{i}}}
\newcommand{\dx}{\ensuremath{\mathrm{\,d}}}
%\newcommand{\dmu}{\ensuremath{\dx \mu(\sa)}}
\newcommand{\dmu}{\ensuremath{\dx \Omega}}

\newcommand{\cconj}{\ensuremath{\ast}}  % note command here is renew...
\newcommand{\sa}{\ensuremath{\omega}}
\newcommand{\saa}{\ensuremath{\theta}}
\newcommand{\sab}{\ensuremath{\phi}}
\newcommand{\sas}{\ensuremath{\saa, \sab}}
\newcommand{\eul}{\ensuremath{\mathbf{\rho}}}
\newcommand{\euls}{\ensuremath{\eula, \eulb, \eulc}}
\newcommand{\eula}{\ensuremath{\alpha}}
\newcommand{\eulb}{\ensuremath{\beta}}
\newcommand{\eulc}{\ensuremath{\gamma}}

\newcommand{\spo}{\ensuremath{\Pi}}

\newcommand{\cocycle}{\ensuremath{\lambda}}
\newcommand{\pnorm}{\ensuremath{p}}
\newcommand{\pnormtext}{\ensuremath{\pnorm}-norm}
\newcommand{\realno}{\ensuremath{\mathbb{R}}}
\newcommand{\natno}{\ensuremath{\mathbb{N}}}
\newcommand{\sothree}{\ensuremath{\mathrm{SO}(3)}}
\newcommand{\sphere}{\ensuremath{{S^2}}}
\newcommand{\dil}{\ensuremath{\mathcal{D}}}
\newcommand{\dilsmall}{\ensuremath{d}}
%\newcommand{\dil}{\ensuremath{D}}
\newcommand{\rot}{\ensuremath{R}}
\newcommand{\scalea}{\ensuremath{a}}
\newcommand{\scaleb}{\ensuremath{b}}
\newcommand{\scaleab}{\ensuremath{{\scalea,\scaleb}}}

\newcommand{\rad}{\ensuremath{r}}
\newcommand{\p}{\ensuremath{^\prime}}
\newcommand{\pp}{\ensuremath{^{\prime\prime}}}
\newcommand{\wav}{\ensuremath{\psi}}

\newcommand{\grideula}{\ensuremath{\mathcal{E}_1}}
\newcommand{\grideulb}{\ensuremath{\mathcal{E}_2}}
\newcommand{\gridhpix}{\ensuremath{\mathcal{H}}}
\newcommand{\gridecp}{\ensuremath{\mathcal{C}}}
\newcommand{\ia}{\ensuremath{{n_\eula}}}
\newcommand{\ib}{\ensuremath{{n_\eulb}}}
\newcommand{\ig}{\ensuremath{{n_\eulc}}}
\newcommand{\ik}{\ensuremath{{k}}}
\newcommand{\ith}{\ensuremath{{n_\saa}}}
\newcommand{\iph}{\ensuremath{{n_\sab}}}
\newcommand{\ipix}{\ensuremath{{p}}}
\newcommand{\na}{\ensuremath{{N_\eula}}}
\newcommand{\nb}{\ensuremath{{N_\eulb}}}
\newcommand{\ngm}{\ensuremath{{N_\eulc}}}
\newcommand{\nth}{\ensuremath{{N_\saa}}}
\newcommand{\nph}{\ensuremath{{N_\sab}}}
\newcommand{\n}{\ensuremath{{N}}}
\newcommand{\nside}{\ensuremath{{N_{\rm{side}}}}}
\newcommand{\npix}{\ensuremath{{N_{\rm{pix}}}}}
\newcommand{\weight}{\ensuremath{w}}
\newcommand{\order}{\ensuremath{\mathcal{O}}}

\newcommand{\astroph}
	{astro-ph}
\newcommand{\mnras}
	{\mbox{MNRAS}}

fcswt_ieee.tex
%------------------------------------------------------------------------------
% Fast directional continuous spherical wavelet transform algorithms
% Version: IEEE Transactions on Signal Processing
%------------------------------------------------------------------------------

%% Adapted from: bare_jrnl.tex
%%               V1.2
%%               2002/11/18
%%               by Michael Shell
%%               mshell@ece.gatech.edu

% Note that the a4paper option is mainly intended so that authors in
% countries using A4 can easily print to A4 and see how their papers will
% look in print. Authors are encouraged to use U.S. letter paper when 
% submitting to IEEE. Use the testflow package mentioned above to verify
% correct handling of both paper sizes by the author's LaTeX system.
%
% Also note that the "draftcls" or "draftclsnofoot", not "draft", option
% should be used if it is desired that the figures are to be displayed in
% draft mode.

%------------------------------------------------------------------------------

% Define constant options.

\newif\ifcol
\coltrue

\newif\ifbw
\bwfalse

\newif\ifdraft
\draftfalse
% \drafttrue

%\newcommand{\pagesize}{a4paper}
\newcommand{\pagesize}{letter}

% If draft set col status to false.
\ifdraft
  \colfalse
\fi


\newlength{\mapwidth}
\setlength{\mapwidth}{55mm}
\ifdraft
  \setlength{\mapwidth}{48mm}
\fi


\newlength{\coeffwidth}
\setlength{\coeffwidth}{62mm}


%------------------------------------------------------------------------------

\ifdraft
  \documentclass[11pt,draftcls,\pagesize,onecolumn,twoside]{IEEEtran}
\else
  \documentclass[10pt,journal,\pagesize,twoside]{IEEEtran}
\fi

\usepackage{url}

\usepackage{ifpdf}
\ifpdf
  \usepackage[pdftex,
    pdfauthor     = {Jason McEwen},
    pdftitle      = {Fast directional CSWT algorithms},
    pdfsubject    = {Continuous spherical wavelets},
    pdfkeywords   = {Wavelet transforms -- spheres -- convolution},
    pdfproducer   = {pdflatex with hyperref},
    pdfcreator    = {pdflatex}]{hyperref}
\else
  \usepackage[ps2pdf,colorlinks=true]{hyperref}
  \hypersetup{%
    pdfauthor     = {Jason McEwen},
    pdftitle      = {Fast directional CSWT algorithms},
    pdfsubject    = {Continuous spherical wavelets},
    pdfkeywords   = {Wavelet transforms -- spheres -- convolution},
    pdfauthor     = {ps2pdf with hyperref},
    pdfcreator    = {LateX},
  }
\fi

\usepackage{cite}      % Written by Donald Arseneau
                        % V1.6 and later of IEEEtran pre-defines the format
                        % of the cite.sty package \cite{} output to follow
                        % that of IEEE. Loading the cite package will
                        % result in citation numbers being automatically
                        % sorted and properly "ranged". i.e.,
                        % [1], [9], [2], [7], [5], [6]
                        % (without using cite.sty)
                        % will become:
                        % [1], [2], [5]--[7], [9] (using cite.sty)
                        % cite.sty's \cite will automatically add leading
                        % space, if needed. Use cite.sty's noadjust option
                        % (cite.sty V3.8 and later) if you want to turn this
                        % off. cite.sty is already installed on most LaTeX
                        % systems. The latest version can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/contrib/supported/cite/

%\usepackage{graphicx}  % Written by David Carlisle and Sebastian Rahtz
                        % Required if you want graphics, photos, etc.
                        % graphicx.sty is already installed on most LaTeX
                        % systems. The latest version and documentation can
                        % be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/required/graphics/
                        % Another good source of documentation is "Using
                        % Imported Graphics in LaTeX2e" by Keith Reckdahl
                        % which can be found as esplatex.ps and epslatex.pdf
                        % at: http://www.ctan.org/tex-archive/info/
% NOTE: for dual use with latex and pdflatex, instead load graphicx like:
%\ifx\pdfoutput\undefined
%  \usepackage[dvips]{graphicx}
%\else
%  \usepackage[pdftex]{graphicx}
%\fi

\ifpdf
  \usepackage[pdftex]{graphicx}
\else
  \usepackage[dvips]{graphicx}
\fi

% However, be warned that pdflatex will require graphics to be in PDF
% (not EPS) format and will preclude the use of PostScript based LaTeX
% packages such as psfrag.sty and pstricks.sty. IEEE conferences typically
% allow PDF graphics (and hence pdfLaTeX). However, IEEE journals do not
% (yet) allow image formats other than EPS or TIFF. Therefore, authors of
% journal papers should use traditional LaTeX with EPS graphics.

%\usepackage{psfrag}    % Written by Craig Barratt, Michael C. Grant,
                        % and David Carlisle
                        % This package allows you to substitute LaTeX
                        % commands for text in imported EPS graphic files.
                        % In this way, LaTeX symbols can be placed into
                        % graphics that have been generated by other
                        % applications. You must use latex->dvips->ps2pdf
                        % workflow (not direct pdf output from pdflatex) if
                        % you wish to use this capability because it works
                        % via some PostScript tricks. Alternatively, the
                        % graphics could be processed as separate files via
                        % psfrag and dvips, then converted to PDF for
                        % inclusion in the main file which uses pdflatex.
                        % Docs are in "The PSfrag System" by Michael C. Grant
                        % and David Carlisle. There is also some information 
                        % about using psfrag in "Using Imported Graphics in
                        % LaTeX2e" by Keith Reckdahl which documents the
                        % graphicx package (see above). The psfrag package
                        % and documentation can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/contrib/supported/psfrag/

\usepackage{subfigure} % Written by Steven Douglas Cochran
                        % This package makes it easy to put subfigures
                        % in your figures. i.e., "figure 1a and 1b"
                        % Docs are in "Using Imported Graphics in LaTeX2e"
                        % by Keith Reckdahl which also documents the graphicx
                        % package (see above). subfigure.sty is already
                        % installed on most LaTeX systems. The latest version
                        % and documentation can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/contrib/supported/subfigure/

%\usepackage{url}       % Written by Donald Arseneau
                        % Provides better support for handling and breaking
                        % URLs. url.sty is already installed on most LaTeX
                        % systems. The latest version can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/contrib/other/misc/
                        % Read the url.sty source comments for usage information.

%\usepackage{stfloats}  % Written by Sigitas Tolusis
                        % Gives LaTeX2e the ability to do double column
                        % floats at the bottom of the page as well as the top.
                        % (e.g., "\begin{figure*}[!b]" is not normally
                        % possible in LaTeX2e). This is an invasive package
                        % which rewrites many portions of the LaTeX2e output
                        % routines. It may not work with other packages that
                        % modify the LaTeX2e output routine and/or with other
                        % versions of LaTeX. The latest version and
                        % documentation can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/contrib/supported/sttools/
                        % Documentation is contained in the stfloats.sty
                        % comments as well as in the presfull.pdf file.
                        % Do not use the stfloats baselinefloat ability as
                        % IEEE does not allow \baselineskip to stretch.
                        % Authors submitting work to the IEEE should note
                        % that IEEE rarely uses double column equations and
                        % that authors should try to avoid such use.
                        % Do not be tempted to use the cuted.sty or
                        % midfloat.sty package (by the same author) as IEEE
                        % does not format its papers in such ways.

\usepackage{amssymb}    % Required for some math symbols uses, e.g. R for real numbers.

%\usepackage{amsmath}   % From the American Mathematical Society
                        % A popular package that provides many helpful commands
                        % for dealing with mathematics. Note that the AMSmath
                        % package sets \interdisplaylinepenalty to 10000 thus
                        % preventing page breaks from occurring within multiline
                        % equations. Use:
                        %\interdisplaylinepenalty=2500
                        % after loading amsmath to restore such page breaks
                        % as IEEEtran.cls normally does. amsmath.sty is already
                        % installed on most LaTeX systems. The latest version
                        % and documentation can be obtained at:
                        % http://www.ctan.org/tex-archive/macros/latex/required/amslatex/math/


% Other popular packages for formatting tables and equations include:

%\usepackage{array}
% Frank Mittelbach's and David Carlisle's array.sty which improves the
% LaTeX2e array and tabular environments to provide better appearances and
% additional user controls. array.sty is already installed on most systems.
% The latest version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/tools/

% Mark Wooding's extremely powerful MDW tools, especially mdwmath.sty and
% mdwtab.sty which are used to format equations and tables, respectively.
% The MDWtools set is already installed on most LaTeX systems. The lastest
% version and documentation is available at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/supported/mdwtools/


% V1.6 of IEEEtran contains the IEEEeqnarray family of commands that can
% be used to generate multiline equations as well as matrices, tables, etc.

% Also of notable interest:

% Scott Pakin's eqparbox package for creating (automatically sized) equal
% width boxes. Available:
% http://www.ctan.org/tex-archive/macros/latex/contrib/supported/eqparbox/

% Notes on hyperref:
% IEEEtran.cls attempts to be compliant with the hyperref package, written
% by Heiko Oberdiek and Sebastian Rahtz, which provides hyperlinks within
% a document as well as an index for PDF files (produced via pdflatex).
% However, it is a tad difficult to properly interface LaTeX classes and
% packages with this (necessarily) complex and invasive package. It is
% recommended that hyperref not be used for work that is to be submitted
% to the IEEE. Users who wish to use hyperref *must* ensure that their
% hyperref version is 6.72u or later *and* IEEEtran.cls is version 1.6b
% or later. The latest version of hyperref can be obtained at:
%
% http://www.ctan.org/tex-archive/macros/latex/contrib/supported/hyperref/
%
% Also, be aware that cite.sty (as of version 3.9, 11/2001) and hyperref.sty
% (as of version 6.72t, 2002/07/25) do not work optimally together.
% To mediate the differences between these two packages, IEEEtran.cls, as
% of v1.6b, predefines a command that fools hyperref into thinking that
% the natbib package is being used - causing it not to modify the existing
% citation commands, and allowing cite.sty to operate as normal. However,
% as a result, citation numbers will not be hyperlinked. Another side effect
% of this approach is that the natbib.sty package will not properly load
% under IEEEtran.cls. However, current versions of natbib are not capable
% of compressing and sorting citation numbers in IEEE's style - so this
% should not be an issue. If, for some strange reason, the user wants to
% load natbib.sty under IEEEtran.cls, the following code must be placed
% before natbib.sty can be loaded:
%
% \makeatletter
% \let\NAT@parse\undefined
% \makeatother
%
% Hyperref should be loaded differently depending on whether pdflatex
% or traditional latex is being used:
%
%\ifx\pdfoutput\undefined
%  \usepackage[dvips,hypertex]{hyperref}
%\else
%  \usepackage[pdftex,hypertexnames=false]{hyperref}
%\fi
%
% Pdflatex produces superior hyperref results and is the recommended
% compiler for such use.

%------------------------------------------------------------------------------

% Correct any bad hyphenation here
%\hyphenation{op-tical net-works semi-conduc-tor}

% Include commands.
\include{fcswt_commands}
\newcommand{\etal}{\mbox{\it{et al.}}}

%\ifx\pdfoutput\undefined
  \renewcommand{\urlsize}{\footnotesize}
%\else
%  \renewcommand{\urlsize}{\scriptsize}
%\fi

\begin{document}
%
%
\title{Fast directional continuous spherical wavelet transform algorithms}
%
%
% author names and IEEE memberships
% note positions of commas and nonbreaking spaces ( ~ ) LaTeX will not break
% a structure at a ~ so this keeps an author's name from being broken across
% two lines.
% use \thanks{} to gain access to the first footnote area
% a separate \thanks must be used for each paragraph as LaTeX2e's \thanks
% was not built to handle multiple paragraphs
\author{J.~D.~McEwen,~M.~P.~Hobson,~D.~J.~Mortlock,~A.~N.~Lasenby%
%        Michael~Shell,~\IEEEmembership{Member,~IEEE,}
%        John~Doe,~\IEEEmembership{Fellow,~OSA,}
%        and~Jane~Doe,~\IEEEmembership{Life~Fellow,~IEEE}% <-this % stops a space
\thanks{Manuscript received 14 June, 2005}% <-this % stops a space
\thanks{J.~D.~McEwen,~M.~P.~Hobson and A.~N.~Lasenby are with the 
        Astrophysics Group, Cavendish Laboratory, Cambridge, UK.}
\thanks{D.~J.~Mortlock is with the Institute of Astronomy, Cambridge, UK.}
\thanks{E-mail: mcewen@mrao.cam.ac.uk (J.~D.~McEwen)}}
% note the % following the last \IEEEmembership and also the first \thanks - 
% these prevent an unwanted space from occurring between the last author name
% and the end of the author line. i.e., if you had this:
% 
% \author{....lastname \thanks{...} \thanks{...} }
%                     ^------------^------------^----Do not want these spaces!
%
% a space would be appended to the last name and could cause every name on that
% line to be shifted left slightly. This is one of those "LaTeX things". For
% instance, "A\textbf{} \textbf{}B" will typeset as "A B" not "AB". If you want
% "AB" then you have to do: "A\textbf{}\textbf{}B"
% \thanks is no different in this regard, so shield the last } of each \thanks
% that ends a line with a % and do not let a space in before the next \thanks.
% Spaces after \IEEEmembership other than the last one are OK (and needed) as
% you are supposed to have spaces between the names. For what it is worth,
% this is a minor point as most people would not even notice if the said evil
% space somehow managed to creep in.
%
% The paper headers
\markboth{IEEE Transactions on Signal Processing,~Vol.~--, No.~--,~June~2005}%
{McEwen \MakeLowercase{\textit{et al.}}: Fast directional \cswt\ algorithms}
% The only time the second header will appear is for the odd numbered pages
% after the title page when using the twoside option.
% 
% *** Note that you probably will NOT want to include the author's name in ***
% *** the headers of peer review papers.                                   ***

% If you want to put a publisher's ID mark on the page
% (can leave text blank if you just want to see how the
% text height on the first page will be reduced by IEEE)
%\pubid{0000--0000/00\$00.00~\copyright~2002 IEEE}

% use only for invited papers
%\specialpapernotice{(Invited Paper)}

% make the title area
\maketitle

%------------------------------------------------------------------------------

% Abstract.
\input{fcswt_abstract}

% Keywords.
\begin{keywords}
Wavelet transforms, spheres, convolution.
\end{keywords}

% For peerreview papers, inserts a page break and creates the second title.
% Will be ignored for other modes.
\IEEEpeerreviewmaketitle

% Main body.
\input{fcswt_body}


% Reminder: the "draftcls" or "draftclsnofoot", not "draft", class option
% should be used if it is desired that the figures are to be displayed while
% in draft mode.


% An example of a double column floating figure using two subfigures.
% (The subfigure.sty package must be loaded for this to work.)
% The subfigure \label commands are set within each subfigure command, the
% \label for the overall fgure must come after \caption.
% \hfil must be used as a separator to get equal spacing
%
%\begin{figure*}
%\centerline{\subfigure[Case I]{\includegraphics[width=2.5in]{subfigcase1}
% where an .eps filename suffix will be assumed under latex, 
% and a .pdf suffix will be assumed for pdflatex
%\label{fig_first_case}}
%\hfil
%\subfigure[Case II]{\includegraphics[width=2.5in]{subfigcase2}
% where an .eps filename suffix will be assumed under latex, 
% and a .pdf suffix will be assumed for pdflatex
%\label{fig_second_case}}}
%\caption{Simulation results}
%\label{fig_sim}
%\end{figure*}

% An example of a floating table. Note that, for IEEE style tables, the 
% \caption command should come BEFORE the table. Table text will default to
% \footnotesize as IEEE normally uses this smaller font for tables.
% The \label must come after \caption as always.
%
%\begin{table}
%% increase table row spacing, adjust to taste
%\renewcommand{\arraystretch}{1.3}
%\caption{An Example of a Table}
%\label{table_example}
%\centering
%% Some packages, such as MDW tools, offer better commands for making tables
%% than the plain LaTeX2e tabular which is used here.
%\begin{tabular}{|c||c|}
%\hline
%One & Two\\
%\hline
%Three & Four\\
%\hline
%\end{tabular}
%\end{table}

% use section* for acknowledgement
%\section*{Acknowledgment}
% optional entry into table of contents (if used)
%\addcontentsline{toc}{section}{Acknowledgment}
%The authors would like to thank...

% trigger a \newpage just before the given reference
% number - used to balance the columns on the last page
% adjust value as needed - may need to be readjusted if
% the document is modified later
%\IEEEtriggeratref{8}
% The "triggered" command can be changed if desired:
%\IEEEtriggercmd{\enlargethispage{-5in}}

% references section
% NOTE: BibTeX documentation can be easily obtained at:
% http://www.ctan.org/tex-archive/biblio/bibtex/contrib/doc/

% can use a bibliography generated by BibTeX as a .bbl file
% standard IEEE bibliography style from:
% http://www.ctan.org/tex-archive/macros/latex/contrib/supported/IEEEtran/bibtex
\bibliographystyle{IEEEtran}
% argument is your BibTeX string definitions and bibliography database(s)
\bibliography{fcswt_bibabbrv,fcswt_bib}
%
% <OR> manually copy in the resultant .bbl file
% set second argument of \begin to the number of references
% (used to reserve space for the reference number labels box)
%\begin{thebibliography}{1}

%\bibitem{IEEEhowto:kopka}
%H.~Kopka and P.~W. Daly, \emph{A Guide to {\LaTeX}}, 3rd~ed.\hskip 1em plus
%  0.5em minus 0.4em\relax Harlow, England: Addison-Wesley, 1999.

%\end{thebibliography}

% biography section
% 
% If you have an EPS/PDF photo (graphicx package needed) extra braces are
% needed around the contents of the optional argument to biography to prevent
% the LaTeX parser from getting confused when it sees the complicated
% \includegraphics command within an optional argument. (You could create
% your own custom macro containing the \includegraphics command to make things
% simpler here.)
%\begin{biography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]{mshell}}]{Michael Shell}
% where an .eps filename suffix will be assumed under latex, and a .pdf suffix
% will be assumed for pdflatex; or if you just want to reserve a space for
% a photo:

\begin{biography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]%
{photos/jmcewen}}]{Jason McEwen}
was born in Wellington, New Zealand, in August 1979.
He received a B.E.\ (Hons) degree in Electrical and Computer Engineering
from the University of Canterbury, New Zealand, in 2002.  

Currently, he is working towards a Ph.D.\ degree at the Astrophysics Group,
Cavendish Laboratory, Cambridge.  His area of interests include
spherical wavelets and the cosmic microwave background.
\end{biography}

% if you will not have a photo at all:
\begin{biography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]%
{photos/mhobson}}]{Michael Hobson} 
was born in Birmingham, England, in September 1967. He
received the B.A.\ degree in natural sciences with honours and
the Ph.D.\ degree in astrophysics from the University of
Cambridge, England, in 1989 and 1993 respectively.

Since 1993, he has been a member of the Astrophysics Group of the
Cavendish Laboratory at the University of Cambridge, where he has been
a Reader in Astrophysics and Cosmology since 2003. His research
interests include theoretical and observational cosmology,
particularly anisotropies in the cosmic microwave background,
gravitation, Bayesian analysis techniques and theoretical optics. 
\end{biography}

% insert where needed to balance the two columns on the last page
%\newpage

\begin{biography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]%
{photos/dmortlock}}]{Daniel Mortlock}
was born in Melbourne, Australia, in June 1973.
He received the B.Sc.\ and Ph.D.\ degrees in science from The University
of Melbourne, Melbourne, Australia, in 1994 and 2000, respectively.

From 1999 to 2002 he was research associate at The Cavendish Laboratory,
Cambridge, UK.  In 2002 he moved to The Institute of Astronomy,
Cambridge, UK as a postdoctoral fellow.  His research interests include
cosmology, gravitational lensing and statistical analysis of large data-sets.
\end{biography}

\begin{biography}[{\includegraphics[width=1in,height=1.25in,clip,keepaspectratio]%
{photos/alasenby}}]{Anthony Lasenby}
was born in Malvern, England, in June 1954. He
received a B.A.\ then M.A.\ from the University of Cambridge in
1975 and 1979, an M.Sc.\ from Queen Mary College, London in 1978
and a Ph.D.\ from the University of Manchester in 1981. 

His Ph.D.\
work was carried out at the Jodrell Bank Radio Observatory
specializing in the Cosmic Microwave Background, which has been a
major subject of his research ever since. After a brief period at
the National Radio Astronomy Observatory in America, he moved
from Manchester to Cambridge in 1984, and has been at the
Cavendish Laboratory Cambridge since then. He is currently Head
of the Astrophysics Group and the Mullard Radio Astronomy
Observatory in the Cavendish Laboratory, and a Deputy Head of the
Laboratory. His other main interests include theoretical physics
and cosmology, the application of new geometric techniques in
computer graphics and electromagnetic modelling, and statistical
techniques in data analysis.
\end{biography}

% You can push biographies down or up by placing
% a \vfill before or after them. The appropriate
% use of \vfill depends on what kind of text is
% on the last page and whether or not the columns
% are being equalized.

%\vfill

% Can be used to pull up biographies so that the bottom of the last one
% is flush with the other column.
%\enlargethispage{-5in}


\end{document}



fcswt_ieee.bbl
\begin{thebibliography}{10}
\providecommand{\url}[1]{#1}
\csname url@rmstyle\endcsname
\providecommand{\newblock}{\relax}
\providecommand{\bibinfo}[2]{#2}
\providecommand\BIBentrySTDinterwordspacing{\spaceskip=0pt\relax}
\providecommand\BIBentryALTinterwordstretchfactor{4}
\providecommand\BIBentryALTinterwordspacing{\spaceskip=\fontdimen2\font plus
\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
  \fontdimen4\font\relax}
\providecommand\BIBforeignlanguage[2]{{%
\expandafter\ifx\csname l@#1\endcsname\relax
\typeout{** WARNING: IEEEtran.bst: No hyphenation pattern has been}%
\typeout{** loaded for the language `#1'. Using the pattern for}%
\typeout{** the default language instead.}%
\else
\language=\csname l@#1\endcsname
\fi
#2}}

\bibitem{bennett:1996}
C.~L. Bennett, A.~Banday, K.~M. Gorski, G.~Hinshaw, P.~Jackson, P.~Keegstra,
  A.~Kogut, G.~F. Smoot, D.~T. Wilkinson, and E.~L. Wright, ``4-year {COBE DMR}
  cosmic microwave background observations: {M}aps and basic results,''
  \emph{Astrophys.\ J.\ Lett.}, vol. 464, no.~1, pp. 1--4, 1996.

\bibitem{bennett:2003}
C.~L. Bennett, M.~Halpern, G.~Hinshaw, N.~Jarosik, A.~Kogut, M.~Limon, S.~S.
  Meyer, L.~Page, D.~N. Spergel, G.~S. Tucker, E.~Wollack, E.~L. Wright,
  C.~Barnes, M.~R. Greason, R.~S. Hill, E.~Komatsu, M.~R. Nolta, N.~Odegard,
  H.~V. Peiris, L.~Verde, and J.~L. Weiland, ``First-year {W}ilkinson
  {M}icrowave {A}nisotropy {P}robe ({WMAP}) observations: {P}reliminary maps
  and basic results,'' \emph{Astrophys.\ J.\ Supp.}, vol. 148, no.~1, pp.
  1--28, 2003.

\bibitem{wieczorek:2006}
M.~A. Wieczorek, ``Gravity and topology of terrestrial planets,''
  \emph{submitted to Treatise on Geophysics}, 2006.

\bibitem{wieczorek:1998}
M.~A. Wieczorek and R.~J. Phillips, ``Potential anomalies on a sphere:
  Applications to the thickness of the lunar crust,'' \emph{J.~Geophys.~Res.},
  vol. 103, pp. 383--395, 1998.

\bibitem{turcotte:1981}
D.~L. Turcotte, R.~J. Willemann, W.~F. Haxby, and J.~Norberry, ``Role of
  membrane stresses in the support of planetary topology,''
  \emph{J.~Geophys.~Res.}, vol.~86, pp. 3951--3959, 1981.

\bibitem{whaler:1994}
K.~A. Whaler, ``Downward continuation of {M}agsat lithosphere anomalies to the
  {E}arth's surface,'' \emph{Geophys.~J.~R. astr.~Soc.}, vol. 116, pp.
  267--278, 1994.

\bibitem{swenson:2002}
S.~Swenson and J.~Wahr, ``Methods for inferring regional surface-mass anomalies
  from {GRACE} measurements of time-variable gravity,''
  \emph{J.~Geophys.~Res.}, vol. 107, pp. 2193--278, 2002.

\bibitem{simons:2006}
F.~J. Simons, F.~A. Dahlen, and M.~A. Wieczorek, ``Spatiospectral concentration
  on a sphere,'' \emph{{SIAM} {R}ev.}, p. in press, 2006.

\bibitem{ramamoorthi:2004}
R.~Ramamoorthi and P.~Hanrahan, ``A signal processing framework for
  reflection,'' \emph{ACM Transactions on Graphics}, vol.~23, no.~4, pp.
  1004--1042, 2004.

\bibitem{choi:1999}
C.~H. Choi, J.~Ivanic, M.~S. Gordon, and K.~Ruedenberg, ``Rapid and stable
  determination of rotation matrices between spherical harmonics by direct
  recursion,'' \emph{J.~Chem.~Phys.}, vol. 111, no.~19, pp. 8825--8831, 1999.

\bibitem{ritchie:1999}
D.~W. Ritchie and G.~J.~L. Kemp, ``Fast computation, rotation and comparison of
  low resolution spherical harmonic molecular surfaces,''
  \emph{J.~Comput.~Chem.}, vol.~20, no.~4, pp. 383--395, 1999.

\bibitem{schroder:1995}
P.~Schr{\"o}der and W.~Sweldens, ``Spherical wavelets: {E}fficiently
  representing functions on the sphere,'' in \emph{Computer Graphics
  Proceedings ({SIGGRAPH} `95)}, 1995, pp. 161--172.

\bibitem{sweldens:1996}
W.~Sweldens, ``The lifting scheme: a custom-design constriction of biorthogonal
  wavelets,'' \emph{Applied Comput.\ Harm.\ Anal.}, vol.~3, no.~2, pp.
  186--200, 1996.

\bibitem{narcowich:1996}
F.~J. Narcowich and J.~D. Ward, ``Non-stationary wavelets on the {\it m}-sphere
  for scattered data,'' \emph{Applied Comput.\ Harm.\ Anal.}, vol.~3, pp.
  324--336, 1996.

\bibitem{potts:1996}
D.~Potts, G.~Steidl, and M.~Tasche, ``Kernels of spherical harmonics and
  spherical frames,'' in \emph{Advanced Topic in Multivariate Approximation},
  F.~Fontanella, K.~Jetter, and P.-J. Laurent, Eds.\hskip 1em plus 0.5em minus
  0.4em\relax Singapore: World Scientific, 1996, pp. 1--154.

\bibitem{freeden:1997a}
W.~Freeden and U.~Windheuser, ``Combined spherical harmonic and wavelet
  expansion -- a future concept in the {E}arth's gravitational determination,''
  \emph{Applied Comput.\ Harm.\ Anal.}, vol.~4, pp. 1--37, 1997.

\bibitem{freeden:1997b}
W.~Freeden, T.~Gervens, and M.~Schreiner, \emph{Constructive approximation on
  the sphere -- with application to geomathematics}.\hskip 1em plus 0.5em minus
  0.4em\relax Clarendon Press, Oxford, 1997.

\bibitem{torresani:1995}
B.~Torr\'{e}sani, ``Position-frequency analysis for signals defined on
  spheres,'' \emph{Signal Proc.}, vol.~43, pp. 341--346, 1995.

\bibitem{dahlke:1996}
S.~Dahlke and P.~Maass, ``Continuous wavelet transforms with applications to
  analyzing functions on sphere,'' \emph{J.\ Fourier Anal.\ and Appl.}, vol.~2,
  pp. 379--396, 1996.

\bibitem{holschneider:1996}
M.~Holschneider, ``Continuous wavelet transforms on the sphere,'' \emph{J.\
  Math.\ Phys.}, vol.~37, pp. 4156--4165, 1996.

\bibitem{antoine:1999}
J.-P. Antoine and P.~Vandergheynst, ``Wavelets on the 2-sphere : A group
  theoretical approach,'' \emph{Applied Comput.\ Harm.\ Anal.}, vol.~7, pp.
  1--30, 1999.

\bibitem{antoine:1998}
------, ``Wavelets on the n-sphere and related manifolds,'' \emph{J.\ Math.\
  Phys.}, vol.~39, no.~8, pp. 3987--4008, 1998.

\bibitem{antoine:2002}
J.-P. Antoine, L.~Demanet, L.~Jacques, and P.~Vandergheynst, ``Wavelets on the
  sphere: Implementation and approximations,'' \emph{Applied Comput.\ Harm.\
  Anal.}, vol.~13, no.~3, pp. 177--200, 2002.

\bibitem{antoine:2004}
J.-P. Antoine, R.~Murenzi, P.~Vandergheynst, and S.~T. Ali,
  \emph{Two-dimensional wavelets and their relatives}.\hskip 1em plus 0.5em
  minus 0.4em\relax Cambridge: Cambridge University Press, 2004.

\bibitem{bogdanova:2004}
I.~Bogdanova, P.~Vandergheynst, J.-P. Antoine, L.~Jacques, and M.~Morvidone,
  ``Stereographic wavelet frames on the sphere,'' \emph{Applied Comput.\ Harm.\
  Anal.}, vol.~19, no.~2, pp. 223--252, 2005.

\bibitem{demanet:2003}
L.~Demanet and P.~Vandergheynst, ``Gabor wavelets on the sphere,'' in
  \emph{SPIE}, 2003.

\bibitem{wiaux:2005}
Y.~Wiaux, L.~Jacques, , and P.~Vandergheynst, ``Correspondence principle
  between spherical and {E}uclidean wavelets,'' \emph{Astrophys.\ J.}, vol.
  632, no.~1, pp. 15--28, 2005.

\bibitem{mhaskar:2001}
H.~N. Mhasker, F.~J. Narcowich, and J.~D. Ward, ``Representing and analysing
  scattered data on spheres,'' in \emph{Multivariate Approximation and
  Applications}, N.~Dyn, D.~Leviaton, D.~Levin, and A.~Pinkus, Eds.\hskip 1em
  plus 0.5em minus 0.4em\relax Cambridge University Press, 2001.

\bibitem{saff:1997}
E.~B. Saff and A.~B.~J. Kuijlaars, ``Distributing many points on a sphere,''
  \emph{Math. Intelligencer}, vol.~19, pp. 5--11, 1997.

\bibitem{hardin:2004}
D.~P. Hardin and E.~B. Saff, ``Discretizing manifolds via minimum energy
  points,'' \emph{Notices of Amer. Math. Soc.}, vol.~51, no.~10, pp.
  1186--1194, 2004.

\bibitem{du:2003}
Q.~Du, M.~D. Gunzburger, and I.~Ju, ``Constrained centroidal voronoi
  tessellations for surfaces,'' \emph{SIAM J.\ Sci. Comput.}, vol.~24, pp.
  1488--1506, 2003.

\bibitem{crittenden:1998}
R.~G. Crittenden and N.~G. Turok, ``Exactly azimuthal pixelizations of the
  sky,'' \emph{preprint (astro-ph/9806374)}, 1998.

\bibitem{gorski:2005}
K.~M. G\'{o}rski, E.~Hivon, A.~J. Banday, B.~D. Wandelt, F.~K. Hansen,
  M.~Reinecke, and M.~Bartelmann, ``Healpix -- a framework for high resolution
  discretization and fast analysis of data distributed on the sphere,''
  \emph{Astrophys.\ J.}, vol. 622, pp. 759--771, 2005.

\bibitem{doroshkevich:2005}
A.~G. Doroshkevich, P.~D. Naselsky, O.~V. Verkhodanov, D.~I. Novikov, V.~I.
  Turchaninov, I.~D. Novikov, P.~R. Christensen, and L.~Y. Chiang,
  ``Gauss--{L}egendre {S}ky {P}ixelization ({GLESP}) for {CMB} maps,''
  \emph{Int. J. Mod. Phys. D.}, vol.~14, no.~2, pp. 275--290, 2005.

\bibitem{wandelt:2001}
B.~D. Wandelt and K.~M. G\'{o}rski, ``Fast convolution on the sphere,''
  \emph{Phys.\ Rev.\ D.}, vol.~63, no. 123002, pp. 1--6, 2001.

\bibitem{mcewen:2004}
J.~D. McEwen, M.~P. Hobson, A.~N. Lasenby, and D.~J. Mortlock, ``A
  high-significance detection of non-{G}aussianity in the {WMAP} 1-yr data
  using directional spherical wavelets,'' \emph{Mon.\ Not.\ Roy.\ Astr.\ Soc.},
  vol. 359, pp. 1583--1596, 2005.

\bibitem{mcewen:2006b}
------, ``A high-significance detection of non-{G}aussianity in the {WMAP} 3-yr
  data using directional spherical wavelets,'' \emph{submitted to Mon.\ Not.\
  Roy.\ Astr.\ Soc}, 2006.

\bibitem{mcewen:2006}
J.~D. McEwen, P.~Vielva, M.~P. Hobson, E.~Mart\'{\i}nez-Gonz\'{a}lez, and A.~N.
  Lasenby, ``Detection of the {ISW} effect and corresponding dark energy
  constraints made with directional spherical wavelets,'' \emph{submitted to
  Mon.\ Not.\ Roy.\ Astr.\ Soc}, 2006.

\bibitem{tosic:2005}
I.~Tosic, P.~Frossard, and P.~Vandergheynst, ``Progressive coding of 3{D}
  objects based on overcomplete decompositions,'' \emph{submitted to IEEE
  Trans. Circuits Syst. Video Technol.}, 2005.

\bibitem{driscoll:1994}
J.~R. Driscoll and D.~M.~J. Healy, ``Computing {F}ourier transforms and
  convolutions on the sphere,'' \emph{Advances in Applied Mathematics},
  vol.~15, pp. 202--250, 1994.

\bibitem{healy:2003}
D.~M.~J. Healy, D.~Rockmore, P.~J. Kostelec, and S.~S.~B. Moore, ``{FFT}s for
  the 2-sphere -- improvements and variations,'' \emph{J.\ Fourier Anal.\ and
  Appl.}, vol.~9, no.~4, pp. 341--385, 2003.

\bibitem{mohlenkamp:1997}
M.~J. Mohlenkamp, ``A fast transform for spherical harmonics,'' Ph.D.
  dissertation, Yale University, Department of Mathematics, 1997.

\bibitem{mohlenkamp:1999}
------, ``A fast transform for spherical harmonics,'' \emph{J.\ Fourier Anal.\
  and Appl.}, vol.~5, no. 2/3, pp. 159--184, 1999.

\bibitem{brink:1993}
D.~M. Brink and G.~R. Satchler, \emph{Angular Momentum}, 3rd~ed.\hskip 1em plus
  0.5em minus 0.4em\relax Oxford: Clarendon Press, 1993.

\bibitem{risbo:1996}
T.~Risbo, ``Fourier transform summation of {L}egendre series and d-functions,''
  \emph{J.\ Geodesy}, vol.~70, no.~7, pp. 383--396, 1996.

\bibitem{ivanic:1996}
J.~Ivanic and K.~Ruedenberg, ``Rotation matrices for real spherical harmonics.
  {D}irect determination by recursion,'' \emph{J.~Phys.~Chem.~A}, vol. 100,
  no.~15, pp. 6342--6347, 1996.

\bibitem{blanco:1997}
M.~A. Blanco, M.~Fl{\'o}rez, and M.~Bermejo, ``Evaluation of the rotation
  matrices in the basis of real spherical harmonics,'' \emph{J.~Mol.~Struct.
  (Theochem)}, vol. 419, pp. 19--27, 1997.

\bibitem{bulow:2001}
T.~Bulow and K.~Daniilidis, ``Surface representation using spherical harmoniscs
  and {G}abor wavelets on the sphere,'' University of Pennsylvania, Tech. Rep.
  MS-CIS-01-37, 2001.

\bibitem{kaiser:1984}
N.~Kaiser and A.~Stebbins, ``Microwave anisotropy due to cosmic strings,''
  \emph{Nature}, vol. 310, pp. 391--393, 1984.

\bibitem{crittenden:1996}
R.~G. Crittenden and N.~Turok, ``Looking for {L}ambda with the {R}ees-{S}ciama
  effect,'' \emph{Phys.\ Rev.\ Lett.}, vol.~76, pp. 575--578, 1996.

\end{thebibliography}
