\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 24\\
Maxima of random fields 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}

\begin{enumerate}
\item {\em Intrinsic volumes.}. Read Worsley, K.J. (2003).
Detecting activation in fMRI data. Statistical Methods in Medical
Research, 12:401-418 for a fast overview. Cao, J. and Worsley,
K.J. (2001). Applications of random fields in human brain mapping.
In M. Moore (Ed.) Spatial Statistics: Methodological Aspects and
Applications, Springer Lecture Notes in Statistics, 159:169-182.

For any smooth random fields $Y$, define excursion set $A_h = \{ x
: Y(x) > h\}$. Then it was shown by Adler that
$$ P(\sup_{x \in \Omega} Y(x) > h) \approx \mathbb{E} \chi(\Omega \cup A_h) = \sum_{d=0}^N
\mu_d(\Omega)\rho_d(h)$$ where $\mu_d$ is $d$-dimensional
Minkowski functional or {\em intrinsic volume} of $\Omega$ and
$\rho_d$ is $d$-dimentional Euler characteristic (EC) density of
$Y$. Most inference in imaging is done in 3D so $N=3$. For details
on intrinsic volume, read V. Schmidt and E. Spodarev's Joint
Estimators for the Specific Intrinsic Volumes of Stationary Random
Sets (2003).


For search region, we can always contain $\Omega$ in a sphere with
radius $r$. Then the intrinsic volumes are $\mu_0 = \chi(\Omega)
=1$, $\mu_1 = 4r$, $\mu_2=\|\partial \Omega\| = 4\pi r^2$, $\mu_3
= \|\Omega\| = \frac{4}{3}\pi r^3$. If we want to use a 3D
rectangle of size $a \times b \times c$, the intrinsic volumes are
$\mu_0=1, \mu_1= a+ b+c, \mu_2 = ab + bc + ac, \mu_3=abc$. For
regular image grid system with equal voxel size $\delta$ in all
direction, the intrinsic volume can be computed in the following
fashion. Let $V$ be the total number of vertices that forms the
corners of voxels, $E$ be the total number of edges connecting
each vertices, $F$ be the number of faces formed by 4 connected
edges. Then $\mu_0=\chi(\Omega) = V -E+F, \mu_1=(E-2F)\delta,
\mu_2= F\delta^2$.

\item {\em Example: corpus callosum of the brain.} The corpus
callosum (CC) $\Omega$ is the white matter region that connects
the left hemisphere to the right hemisphere. The average shape of
midsagittal cross-section is given in Figure 1. To find the number
of edges and pixels contained in $\Omega$, we start from a seed
pixel $p_0=(40,40)$ that is the splenium of CC and implement
random search algorithm and found $V=328$, $E=599$ and $F=272$.
The resolution of image is $\delta=2 \mbox{ mm }$. So
$\mu_0=1,\mu_1=55\delta,\mu_2=272\delta^2$.

\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.25]{ima24-0.eps}
\includegraphics[scale=0.25]{ima24-1.eps}
\includegraphics[scale=0.25]{ima24-2.eps}
\includegraphics[scale=0.25]{ima24-3.eps}
\caption{Top left: original white matter density of the corpus
callosum. Search region $\Omega$ is defined as the white matter
region with density higher than 0.1. It will be a little bit
larger than the actual corpus callosum shape.}
\end{figure}

\begin{verbatim}
omega=zeros(dim1,dim2);
x=40;y=40;
omega(x,y)=1;

Emap=[3 2 1 0];
E=4;F=1;
for i=1:10000
  cur_x = x+ (-1)^binornd(1,0.5);
  cur_y = y+ (-1)^binornd(1,0.5);
  if (image(cur_x,cur_y) >0.1)
    if (omega(cur_x,cur_y) ~=1)
      n_nbr=omega(cur_x-1,cur_y-1)
            +omega(cur_x+1,cur_y+1)
            +omega(cur_x-1,cur_y+1)
            +omega(cur_x+1,cur_y-1);
      E=E+Emap(n_nbr);
      F=F+1;
    end;
    omega(cur_x,cur_y)=1;
    x=cur_x;
    y=cur_y;
  end;
end;
\end{verbatim}


\item {\em Gaussian fields.} The EC density is given by
$$\rho_d(h) = \mathbb{E} \big[(Y> h)\det( -\ddot{Y_d})|\dot{Y_d}\big]
P(\dot{Y_d}=0).$$ where dot notation indicates differentiation
with respect to the first $d$ components. Computation of this
conditional expectation is nontrivial other than Gaussian fields.
For zero mean Gaussian field $Y$,
$$\rho_0 = 1 - \Phi(h),$$
$$\rho_1 = \lambda^{1/2}\frac{e^{-h^2/2}}{\sqrt{2\pi}},$$
$$\rho_2 = \lambda h\frac{e^{-h^2/2}}{\sqrt{2\pi}},$$
$$\rho_3 = \lambda^{3/2}(h^2-1) \frac{e^{-h^2/2}}{\sqrt{2\pi}}$$
where $\lambda$ will be defined later.

\item {\em $F$-fields.} It is defined in Worsley, K.J. (1994).
Local maxima and the expected Euler characteristic of excursion
sets of $\chi2$, $F$ and $t$ fields. Advances in Applied
Probability, 26:13-42. Let $X_1,\cdots,X_{\alpha}, Y_1,\cdots,
Y_{\beta}$ be i.i.d. stationary zero mean Gaussian fields. Then
$F$-field with $\alpha$ and $\beta$ degrees of freedom is given by
$$F(x) = \frac{\sum_{j=1}^{\alpha} X_j^2(x)/\alpha}{\sum_{j=1}^{\beta} Y_j^2(x)/\beta}.$$
To avoid singularity, we need to assume the total degrees of
freedom $\alpha + \beta \gg N$ to be fairly larger than the
dimension of space. See Worsley (1994) for detail. The EC-density
for $F$-field is given by \bqn \rho_0 &=&\int_h^\infty {{\Gamma
({{\alpha + \beta}\over 2})}\over {\Gamma ({\alpha \over 2})\Gamma
({\beta \over 2})}} {\alpha\over \beta}\left({{\alpha x}\over
\beta }\right)^{(\alpha-2) \over 2} \left(1+{{\alpha x}\over \beta
}\right)^{-{(\alpha+\beta)\over 2}} d x,\\
\rho_1&=&\lambda^{1/2} {{\Gamma ({{\alpha+\beta-1}\over
2})2^{1\over 2}}\over{ \Gamma ({\alpha\over 2})\Gamma ({\beta
\over 2})}} \left({{\alpha h}\over \beta}\right)^{(\alpha-1) \over
2} \left(1+{{\alpha h}\over \beta}\right)^{-{(\alpha+\beta-2)\over
2}},\\
\rho_2&=&\lambda {{\Gamma ({{\alpha +\beta-2}\over 2})}\over
{\Gamma ({\alpha \over 2})\Gamma ({\beta \over 2})}}
\left({{\alpha h}\over \beta }\right)^{(\alpha-2) \over 2}
\left(1+{{\alpha h}\over \beta }\right)^{-{(\alpha +\beta-2)\over 2}}\\
&\times&\left[(\beta -1){{\alpha h}\over \beta}-(\alpha-1)\right].
\eqn

%\rho_3(t)&=&{{(4 \log_e 2)^{3\over 2}}\over{(2\pi)^{3\over 2}}}
%{{\Gamma\left({{\beta+\alpha-3}\over 2}\right)2^{-{1\over
%2}}}\over{ \Gamma\left({\beta\over
%2}\right)\Gamma\left({\alpha\over 2}\right)}}
%\left({{\alphat}\over \beta}\right)^{{1\over 2}(\alpha-3)}
%\left(1+{{\alphat}\over \beta}\right)^{-{1\over 2}(\beta+\alpha-2)}\\
%&\times&\left[(\beta-1)(\beta-2)\left({{\alphat}\over
%\beta}\right)^2-(2\beta \alpha-\beta -\alpha-1)
%\left({{\alphat}\over \beta }\right)+(\alpha-1)(\alpha-2)\right]\\
%\end{eqnarray*}

\item {\em Amount of smoothing.} $\lambda$ measures the smoothness
of fields, defined as the variance of the derivative of component
of $Y$. Consider zero mean unit variance field $Y$, assume that it
can be constructed from white noise $W$ by convolving it with
Gaussian kernel $W$, i.e. $Y=K_{\sigma}*W$. Obviously this
generates isotropic fields with covariance function $R$. Then we
can show that
$${\bf Var} \dot Y =\lambda I$$ where
$\lambda = \frac{4\ln 2}{2\pi \sigma^2}.$ Note that $R_Y(x,y)=
\int K_{\sigma}(x-z)K_{\sigma}(y-z) \;dz$. So \bqn \lambda_{ij}
&=& \mathbb{E} \Big(\frac{\partial Y}{\partial x_i}\frac{\partial
Y}{\partial x_j}\Big)\\ &=& \frac{\partial^2}{\partial x_i
\partial x_j}R_Y(x,x) = \int \frac{\partial
K_{\sigma}(z)}{\partial x_i}\frac{\partial K_{\sigma}(z)}{\partial
x_j} \;dz\eqn


\begin{figure}[t]
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.28]{grid.eps}
\caption{There are only 4 possible configuration of adding
additional pixel to a search region $\Omega$. Each case
corresponds to adding 3, 2 1 and 0 edges and it has been
implemented as a table called {\tt Emap} in the algorithm and
figuring out how many edges to add is implemented as a look up
table.}
\end{figure}


{\bf Problem 37.} Simplify the above expression and show $\lambda
\propto \sigma^{-2}$. If possible generalize the result for the
general anisotropic kernel smoothed image $K_H*W$.


\end{enumerate}
\end{document}