\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 26\\
Curve modeling II.}
\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}
We lecture is based on {\bf Problem 25}.
\begin{enumerate}
\item {\em Direct curvature estimation.} (Shubing Wang)  The
curvature is given by
$$\kappa = \frac{4A({\bf p}_{i-1},{\bf p}_i, {\bf p}_{i+1})}
{|{\bf p}_{i-1} - {\bf p}_i||{\bf p}_i - {\bf p_{i+1}}||{\bf
p}_{i+1} - {\bf p}_{i-1}|}$$ where $A({\bf p}_{i-1},{\bf p}_i,
{\bf p}_{i+1})$ is the area of triangle with vertices ${\bf
p}_{i-1}, {\bf p}_i, {\bf p}_{i+1}$.
\begin{verbatim}
%Modified from Shubing's code
load CCcurve.data
CC=[CCcurve(171,:);CCcurve;CCcurve(1,:)]
for i=2:172
 kappa2(i-1)=4*polyarea(CC((i-1):(i+1)...
       ,1),CC((i-1):(i+1),2))...
       /norm(CC(i-1,:)-CC(i,:))...
       /norm(CC(i+1,:)-CC(i,:))...
       /norm(CC(i+1,:)-CC(i-1,:))
end;
\end{verbatim}
\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.24]{ima25-8.eps}
\includegraphics[scale=0.24]{ima26-1.eps}
\caption{Left: curvature function $\kappa(\lambda)$ via Shubing's
method. Right: finite difference method. Right: scatter plot of
curvature estimations. Two methods gives almost identical
estimation within $0.013\%$ difference.}
\end{figure}
\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
%\includegraphics[scale=0.24]{ima25-1.eps}
\includegraphics[scale=0.24]{ima25-2.eps}
\includegraphics[scale=0.24]{ima25-3.eps}
\includegraphics[scale=0.24]{ima25-4.eps}
\includegraphics[scale=0.24]{ima25-8.eps}
\caption{Top left: arc length function $\lambda$ given as a
function of the vertex index. Top right: arc length
parameterization of $x$-component. Bottom left: arc length
parameterization of $y$-component. Bottom right: curvature
estimation via arclength parameterization.}
\end{figure}

\item {\em Curvature estimation via acrlength parameterization.}
This is trivial once we figure out how to implement arclength
parameterization in {\tt MATLAB}. The first derivatives ${\bf
X}'(\lambda)$ are estimated via the finite difference and the
second derivative follows similarly. Finite difference for non
regular grid is

$$f'(\lambda) = \frac{f(\lambda+\Delta\lambda_1) - f(\lambda)}{\Delta
\lambda_1}$$ and $$f''(\lambda) =
\frac{\frac{f(\lambda+\Delta\lambda_1) - f(\lambda)}{\Delta
\lambda_1} - \frac{f(\lambda) - f(\lambda -\Delta
\lambda_2)}{\Delta \lambda_2}}{\frac{\Delta \lambda_1 + \Delta
\lambda_2}{2}}.$$

\begin{verbatim}
x=CCcurve(:,1);
y=CCcurve(:,2);
dx=diff(x); dy=diff(y);
dl=sqrt(dx.^2 + dy.^2);
%arclength function
for i=1:172
  l(i)=sum(dl(1:i))
end;
%fist difference
Dx=dx./dl(1:172);
figure;plot(l,Dx,'linewidth',2);
Dy=dy./dl(1:172);

%second difference
DDl= (dl(1:171)+dl(2:172))/2;
DDx = diff(Dx)./DDl;
DDy = diff(Dy)./DDl;
L= (l(1:171)+l(2:172))/2;
kappa=sqrt(DDx.^2 + DDy.^2)
figure;plot(L,kappa,'linewidth',2);
\end{verbatim}

\item {\em Diffusion smoothing on a curve}. Note that $\kappa(0) =
\kappa(L)$ where $L$ is the total arclength of $\partial \Omega$.
 Let us display the curvature function in 2D
using the following {\tt MATLAB} codes.
\begin{verbatim}
x=CCcurve(:,1) y=CCcurve(:,2)
map=hot(34)
colormap(map)
for i=1:171
  fc=min(ceil(kappa(i)+0.0001),34)
    plot(x(i),y(i),'s',...
     'MarkerEdgeColor','k',...
     'MarkerFaceColor',map(fc,:),...
     'MarkerSize',7)
    hold on;
end; colorbar
\end{verbatim}
If we smooth $\kappa(\lambda)$ in interval $[0,\lambda(2\pi))$, we
can not grantee the condition $\kappa(0) = \kappa(L)$. How can we
smooth periodic data? We can use periodic spline smoothing or
diffusion smoothing or iterated kernel smoothing on $\partial
\Omega$ directly. Since the diffusion smoothing should gives
identical result to the iterated kernel smoothing, we implement
iterated kernel smoothing ($\sigma = \sqrt 5 \cdot 0.2 =0.45$) as
follows:

\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.5]{ima26-2.eps}
\includegraphics[scale=0.5]{ima26-3.eps}
\includegraphics[scale=0.5]{ima26-4.eps}
\caption{Top: Curvature function $\kappa$. Middle: Smoothed
curvature function.Bottom: Comparison}
\end{figure}
\begin{verbatim}
kappa=[kappa(171) kappa' kappa(1)]
kernel=inline('exp(-x.^2/sigma^2)...
       /sum(exp(-x.^2/sigma^2))')

for j=1:5
  for i=2:172
    nbr_dist=[dl(i-1) 0 dl(i)];
    weight= kernel(0.2,nbr_dist);
    tem(i-1) = dot(kappa((i-1):(i+1)),...
                   weight);
  end;
  kappa=[tem(171) tem tem(1)]
end;
\end{verbatim}
\end{enumerate}
{\bf Problem 39.} Redo {\bf Problem 25} answering following
questions. Given parameterization $X(\lambda)$ of curve $\partial
\Omega$, what is the 1D Laplacian on $\partial \Omega$? Estimate
the Laplacian in our data numerically. Now solve the diffusion
equation on $\partial \Omega$ with initial condition $\kappa$.
Make sure you are not solving the diffusion equation in the
parameter space, which is a wrong way to do diffusion smoothing!


\end{document}

Suppose that ${\bf U}_j$ is the displacement field that maps point
${\bf x}$ in the corpus callosum template $\partial
\Omega_{atlas}$ to the corresponding homologous point ${\bf x} +
{\bf U}_j$ in the corpus callosum of subject $j$, $\partial
\Omega_j$ (Figure 2). Let $\lambda_j$ and $\lambda_{atlas}$ to be
the arc-length element of subject $j$ and template respectively.
Then the arc-length dilatation with respect to the template at
each point ${\bf x} \in \partial \Omega_{atlas}$ is defined as
$$\Delta \lambda_j({\bf x}) = \frac{\lambda_j({\bf x} + {\bf U}({\bf x}))
- \lambda_{atlas}({\bf x})}{\lambda_{atlas}({\bf x})}.$$ The
curvature dilatation $\Delta \kappa_j$ is defined similarly. As in
Chung {\em et al.} (2001) and Chung {\em et al.} (2003), we will
model metric dilatations as Gaussian random fields.