\documentclass[11pt,twocolumn]{article}
\usepackage{times,amsmath,amsfonts}
\usepackage[dvips]{graphicx,graphics}



\setlength{\textheight}{10in} \setlength{\textwidth}{7.1in}
\setlength{\columnsep}{0.3in} \setlength{\topmargin}{0in}
\setlength{\headheight}{-0.3in} \setlength{\headsep}{0in}
\setlength{\oddsidemargin}{-0.6in}  % Centers text.


% Add the period after section numbers.  Adjust spacing.
%\newcommand{\Section}[1]{\vspace{-5pt}\section{\hskip -1em.~~#1}\vspace{-3pt}}
%\newcommand{\SubSection}[1]{\vspace{-2pt}\subsection{\hskip -1em.~~#1}
%        \vspace{-3pt}}
\newcommand{\bqn}{\begin{eqnarray*}}
\newcommand{\eqn}{\end{eqnarray*}}
\newcommand{\bq}{\begin{eqnarray}}
\newcommand{\eq}{\end{eqnarray}}

\begin{document}

\title{Stat 992: Lecture 20\\
Bivariate smoothing on a sphere}
\author{Moo K. Chung \tt{mchung@stat.wisc.edu}}
\maketitle \thispagestyle{empty}

%\begin{figure}
%\centering
%\renewcommand{\baselinestretch}{1}
%\includegraphics[scale=0.8]{homework05-1.eps}
%\end{figure}

\begin{enumerate}

\item {\bf Problem 15.} Find spherical Laplacian.\\
{\em Solution.} (Shubing Wang).
Parametrize unit sphere $S^2$ by spherical coordinates. Let $p =(x,y,z) \in S^2$.
$$x=\sin \theta \cos \psi, \; y=\sin \theta \sin \psi, \; z = \cos \theta.$$There are slightly different parameterizations but this is the usual standard. The tangent vectors which form basis in the tangent space $T_p$ is $$\frac{\partial p}{\partial \theta} = (\cos \theta \cos \psi, \cos \theta \sin \psi, -\sin \theta)$$  and $$\frac{\partial p}{\partial \psi} = (-\sin \theta \sin \psi, \sin \theta \cos \psi, 0)$$
Note that $$\langle\frac{\partial p}{\partial \theta},
\frac{\partial p}{\partial \psi} \rangle = 0, \langle
\frac{\partial p}{\partial \theta},  \frac{\partial p}{\partial
\theta} \rangle = 1, \langle \frac{\partial p}{\partial \psi},
\frac{\partial p}{\partial \psi} \rangle = \sin ^2\theta.$$ So we
have orthogonal coordinates and this is why you can construct
orghonal grid system on a sphere. Then the Riemmanian metric
tensors are $$g_{11}=1, g_{12}=g_{21}=0, g_{22}=\sin^2 \theta.$$
Note that $\det g = \sin^2 \theta$ and $g^{11}=1, g^{12}=g^{21}=0,
g^{22}=1/\sin^2 \theta$. From previous lectures, the
Laplace-Beltrami operator on manifolds is given by
$$ \Delta = \frac{1}{\det G^{1/2}}\sum_{i,j=1}^2
\frac{\partial }{\partial u^i}\Big(\det
G^{1/2}g^{ij}\frac{\partial}{\partial
 u^j}\Big)$$
where $(u^1,u^2)=(\theta,\psi)$. Because of the orthogonality,
terms cancel out and we have \bqn \Delta_{S^2} &=& \frac{1}{\sin
\theta}\Big[ \frac{\partial }{\partial \theta}\Big(\sin \theta
\frac{\partial}{\partial
 \theta}\Big) + \frac{\partial }{\partial \psi}\Big(\frac{1}{\sin \theta} \frac{\partial}{\partial
 \psi}\Big)\Big]\\
 &=& \frac{1}{\sin \theta}\frac{\partial
}{\partial \theta}\Big(\sin \theta \frac{\partial}{\partial
 \theta}\Big) +\frac{1}{\sin^2 \theta} \frac{\partial^2 }{\partial^2 \psi}
 \eqn



\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima20-1.eps}
\includegraphics[scale=0.25]{ima20-2.eps}
\caption{Left: Coordinates $(\theta_i,\psi_i)$ where data are
observed. Sampling density is lower in polar regions. Right: North
pole regions are magnified to show the nonconformity of sampling
density. If one tries to formulate nonparametric regression in
these nonregular data points, one encounters numerical
implementation problem.}
\end{figure}


\item {\em Diffusion smoothing on a sphere.} Based on the spherical Laplacian given above, we will
construct bivariate diffusion smoothing in $S^2$. Suppose we have
observations $Y(x), x \in S^2$. Then we need to solve
$$\frac{\partial f}{\partial t} = \Delta_{S^2} f$$
with the initial condition $f(\theta,\psi,0)=Y(\theta,\psi)$.

Suppose we have polygonal mesh data available. Assume
$p_1,\cdots,p_m$ to be vertices are around $p_0$. The
corresponding spherical coordinates of vertex $p_i=(x_i,y_i)$ is
denoted by $(\theta_i,\psi_i)$. Then following the idea of
estimating Laplacian in polygonal surfaces, we fit a quardratic
polynomials locally via the least-squares estimation to estimate
derivatives $\frac{\partial}{\partial \theta},
\frac{\partial^2}{\partial \theta^2}, \frac{\partial^2}{\partial
\psi^2}$. However, the problem with the spherical coordinates is
that near the both north and south poles $\theta=0,\pi$, the
number of vertices of polygonal surfaces decreases. So when we
flatten the data defined on a sphere to a rectangle $$R=\{ 0 \leq
\theta \leq \pi, 0 \leq \psi < 2\pi\}$$ there is the problem of
lower sampling density on both top and bottom of the rectangular
domain. Thus the local polynomial fitting may perform badly in
these regions. To avoid this parameterization problem, estimation
of the Laplacian near the polar regions are done by rotating the
sphere by 90 degree, i.e. by adding $\theta$ by $\pi/2$. Note that
the Laplacian should be invariant under such coordinate
transformation. Flattening polygonal surface of $S^2$ to rectangle
$R$ destroys the topology so we need to construct a new vertex
connectivity map.

Download data from\\
http://tezpur.keck.waisman.wisc.edu/~chung/
teaching/data/MNImesh/
\begin{verbatim}
[tri,coord,nbr,n]=MNIgetmesh('white.obj');
tem=Scoord.^2; for
i=1:40962
  coord(:,i)=coord(:,i)
              /sqrt(sum(tem(:,i)));
end;
[theta, psi,r] =cart2sph(coord(1,:),
             coord(2,:),coord(3,:));
\end{verbatim}


\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima20-3.eps}
\includegraphics[scale=0.25]{ima20-5.eps}
\caption{Left: Polar regions are magnified to show the
nonuniformity of sampling density. It is numerically nontrivial to
implement most of nonparametric regressions if data are spaced
nonregularly in higher dimensions. Right: Voronoi diagram. The
dual of Voronoi is the Delaunay triangulation.}
\end{figure}


\item {\em Delaunay triangulation.} Given $n$ points
$P=\{p_1,\cdots, p_n \} \subset \mathbb{R}^2$, triangulation
$\mathcal{T}(P)$ of $P$ is a planar subdivision of the convex hull
of $P$ into non overlapping triangles with vertices from $P$.
$\mathcal{T}(P)$ is maximal in a sense that adding additional age
connecting two vertices in $P$ will make $\mathcal{T}(P)$
nonplanar. There are only finitely many different triangulations.

One unique triangulation is given by {\em Delaunay diagram}
$\mathcal{DT}(P)$ which is a dual graph of {\em Voronoi diagram}
$\mathcal{V}(P)$. A Delaunay diagram or Delaunay triangulation is
defined as a triangulation of $V$ such that the circumcircle of
any triangle in $\mathcal{DT}(P)$ does not contain points of $P$
in its interior. Assuming no four or more points of $P$ are
cocircular, it can be shown to be unique. The probability of
finite random points forming a circle is zero so most of Delaunay
triangulation encountered in practice should be unique unless the
coordinates of data are biased somehow. Algorithms for generating
Delaunay diagrams are well known in computional geometry and also
available in {\em MATLAB}. See M.De Berg's Computational Geometry:
Algorithms and Applications, 2nd rev. ed. Berlin: Springer-Verlag,
pp. 147-163, 2000 and H. Edelsbrunner's  Geometry and Topology for
Mesh Generation,

\begin{verbatim}
DT=delaunay(theta,psi);
figure;trimesh(DT,theta,psi,thickness);

voronoi(theta,psi);
\end{verbatim}

{\bf Problem 30.} If $n$ is the number of points in $P$ and $k$
points lie on the boundary of $\mathcal{DT}(P)$, there are $A$
number of triangles and $B$ number of edges in $\mathcal{T}(P)$.
Find $A$ and $B$. Prove your answer and demonstrate for simple
examples using built-in functions in {\tt MATLAB}.

{\bf Problem 31.} Based on Delaunay triangulation, estimate the
Laplacian $\Delta_{S^2}$ via the least squares estimation and do
diffusion smoothing on a sphere. Use {\tt thickness.vv} in the
same directory.
\end{enumerate}
\end{document}