\documentclass[12pt,english]{article}
\usepackage{pslatex}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{geometry}
\geometry{verbose,letterpaper,tmargin=1in,bmargin=1in,lmargin=.75in,rmargin=1in}
\usepackage{setspace}
\onehalfspacing
\makeatletter
\usepackage{babel}
\usepackage{euscript}
\makeatother
\begin{document}
\title{Cratermatic Finder Routine Description}
\author{Michael Mendenhall}
\date{August 2006}
\maketitle
\newcommand{\rfr}[1]{\textbf{(#1)}}
\newcommand{\CM}{{\em Cratermatic}}
\newcommand{\sct}[1]{{\bf #1}}
\newcommand{\bbR}{\ensuremath{\mathbb{R}}}
\newcommand{\bbZ}{\ensuremath{\mathbb{Z}}}
\newcommand{\bbRR}{\ensuremath{\mathbb{R}^2}}
\newcommand{\vc}[1]{\ensuremath{\vec{\mathbf{#1}}}}
\newcommand{\vch}[1]{\ensuremath{\hat{\mathbf{#1}}}}
\newcommand{\grd}{\ensuremath{\vec{\nabla}}}
\newcommand{\mtx}[1]{\ensuremath{\mathrm{\mathbf{\left[#1\right]}}}}
\newcommand{\pd}[2]{\ensuremath{\frac{\partial #1}{\partial #2}}}
\newcommand{\sB}{\ensuremath{\EuScript{B}}}
\newcommand{\sW}{\ensuremath{\EuScript{W}}}
\section{Notation}
\subsection{The Topographic Landscape}
The dataset that \CM\ deals with is a rectangular $n \times m$ array of elevation values, which we will denote by $\mtx{Z}$, the $(i,j)$ element of which is a real number denoted by $Z_{i,j}$. It is often conceptually easier to think of the topography as a continuous, differentiable function $z:S\rightarrow\bbR$ over some domain $S\subset \bbRR$.
\subsection{Derivatives}
For actual computations on the discrete dataset, \CM\ replaces the continuous derivative $\pd{z}{x}$ by the discrete derivative $Dx(\mtx{Z})$,
\[ \left[ \begin{array}{ccc} -\frac{1}{2} & 0 & \frac{1}{2} \end{array} \right] \ast \mtx{Z} : Z_{i,j} \mapsto -\frac{1}{2}Z_{i-1,j} + \frac{1}{2}Z_{i+1,j} \]
(where $\ast$ denotes discrete convolution). Similarly, $\pd{z}{y}$ is replaced by $Dy(\mtx{Z})$,
\[ \left[ \begin{array}{c} -\frac{1}{2} \\ 0 \\ \frac{1}{2} \end{array} \right] \ast \mtx{Z} : Z_{i,j} \mapsto -\frac{1}{2}Z_{i,j-1} + \frac{1}{2}Z_{i,j+1} \]
For an analytic continuous function of two variables, these two partial derivatives are sufficient to specify the slope in any direction at any point. For the discrete dataset, however, it is useful to define additional derivatives in other directions, specifically the diagonal derivatives $Dxy(\mtx{Z})$
\[ \left[ \begin{array}{ccc}
-\frac{1}{2\sqrt{2}} & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & \frac{1}{2\sqrt{2}} \end{array} \right] \ast \mtx{Z} : Z_{i,j} \mapsto -\frac{1}{2\sqrt{2}}Z_{i-1,j-1} + \frac{1}{2\sqrt{2}}Z_{i+1,j+1} \]
and $Dxy'(\mtx{Z})$,
\[ \left[ \begin{array}{ccc}
0 & 0 & \frac{1}{2\sqrt{2}}\\
0 & 0 & 0 \\
-\frac{1}{2\sqrt{2}} & 0 & 0 \end{array} \right] \ast \mtx{Z} : Z_{i,j} \mapsto -\frac{1}{2\sqrt{2}}Z_{i-1,j+1} + \frac{1}{2\sqrt{2}}Z_{i+1,j-1} \]
\subsection{The C-Transform}
The craterfinding transform (C-Transform) is defined for continuous topography as
\[ C\circ z:\vc{x} \mapsto \iint_S e^{-\frac{\left|\vc{x}-\vc{x}'\right|^2}{2r^2}} \grd z(\vc{x}') \cdot \frac{\vc{x}'-\vc{x}}{\left|\vc{x}'-\vc{x}\right|} d\vc{x}' \]
To remove interference from unwanted areas of the terrain (for example, craters that have already been found), the region of integration may exclude those areas (this is implemented by setting the gradient to $0$ in the excluded areas).
\subsection{The Watershed Transform}
For an $m\times n$ topography dataset $\mtx{Z}$, let the watershed transform of $\mtx{Z}$ be denoted by $\sW$, a set of disjoint subsets of the $m\times n$ lattice whose union covers the $m\times n$ lattice.
\section{The "findcraters" routine}
Starting the procedure with target radius $r=5$,
\begin{enumerate}
\item Create a blurred copy of the topography $\mtx{G}$ by applying a Gaussian blur of radius $r/5$ to $\mtx{Z}$
\item Calculate the watershed transform $\sW$ of $\mtx{G}$
\item Calculate the gradient $\mtx{\nabla G}$ of $\mtx{G}$
\item Calculate the radius $r$ C-Transform $\mtx{C}$ of $\mtx{Z}$
\item Find crater candidate regions in $\mtx{C}$:
\begin{enumerate}
\item Calculate the discrete second derivatives of $\mtx{C}$ in 4 directions: $Dx(Dx(\mtx{C}))$, $Dy(Dy(\mtx{C}))$, $Dxy(Dxy(\mtx{C}))$, and $Dxy'(Dxy'(\mtx{C}))$
\item Create a binary image $\mtx{B}$ where $B_{i,j} = 1$ if all four second derivatives are positive at $(i,j)$, and $B_{i,j} = 0$ otherwise
\item Find the path-connected regions of 1's in $\mtx{B}$, each of which is an upward-convex region of $\mtx{C}$ which may be the interior of a crater. Denote this set of regions by $\sB$.
\end{enumerate}
\item Restrict each element $X \in \sB$ to those segments in $\sW$ whose minimum point is contained in $X$:
\begin{enumerate}
\item For each $W\in\sW$, find some minimum point $\vec{p}_W\in W$ such that $\mtx{G}(\vec{p}_W) \leq \mtx{G}(\vec{q})$ $\forall \vec{q}\in W$
\item For each $W\in\sW$, if $\vec{p}_W \notin X$, set $X = X\setminus W$
\end{enumerate}
As a result of this operation, no $W\in\sW$ can contain points from more than one $X\in\sB$. The future expansion of the regions $X\in\sB$ will be limited to the union of regions $W\in\sW$ which intersect $X$: for each $X\in\sB$, let $X'$ denote $\bigcup \{W\in\sW : W\cap X \neq \emptyset \}$.
\item Remove each element of $\sB$ for which 70\% or more of the points have been identified as a crater in a previous step (this prevents double-counting of craters)
\item Expand each region $X \in \sB$ to cover (we hope) the whole crater. This is done using a "cellular automata"-like approach, in which a set of "active points" expands to a new set of "active points" to cover more of the crater. Each "active point" $a$ is associated with a set of characteristics: a position $a_p$, a starting point $a_{p_0}$, an owner $a_X\in\sB$, a highest gradient value $a_g \in \bbR^+$, a distance $a_d\in\bbR^+$, and a "strength" $a_s \in \bbR^+$ which is used to determine when the point has travelled too far over too low slopes. For each region $X\in\sB$, we also track the highest gradient yet seen in the region, $X_g \in \bbR^+$, which starts out with the value $0$.
\begin{enumerate}
\item Create the set $A$ of initial active points
\begin{enumerate}
\item The initial active points are each $\vec{p}_W$ contained in some $X\in\sB$.
\item The owner of each active point is the $X\in\sB$ which contains it
\item The points start with $a_g = 0$, $a_d = 0$, $a_s = 1$
\end{enumerate}
\item Create a new set of active points $A'$ by expanding out from the current set of active points $A$. For each $a\in A$ with owner $X\in\sB$ for which $a_s \geq 10^{-r/5}$, we consider the possibility of creating a new point $a' \in A'$ for each of the 8 pixels surrounding $a$ which are also in $X'$. Let $\vc{dv} \in (\{-1,0,1\} \times \{-1,0,1\})\setminus \{(0,0)\}$ denote the direction of movement from $a$ to $a'$, and let $dz = (\mtx{G}(a'_p) - \mtx{G}(a_p))/|\vc{dv}|$ denote the slope of the move. Expansion to $a'$ is subject to the following restrictions:
\begin{enumerate}
\item If $a'$ has been expanded onto in a previous step, do not expand (no need to repeatedly visit the same area)
\item If $a_p \in a_X$ (i.e. we are still inside the owner region $X$), skip the rest of these restrictions
\item If $dz < 0$, do not expand (only travel uphill)
\item If $\frac{(a'_p - a_{p_0}) \cdot \mtx{\nabla G}(a'_p)}{|a'_p - a_{p_0}|\cdot |\mtx{\nabla G}(a'_p)|} > -0.7$, do not expand (only expand radially outward from starting point)
\end{enumerate}
If several $a\in A$ could potentially expand onto $a'$, the one with the minimum $a_d + |dv|$ is given priority. On expansion, the newly created $a'$ inherits the characteristics of $a$ with the following modifications:
\begin{itemize}
\item $a'_d = a_d + |\vc{dv}|$
\item $a'_g = \max(a_g,dz)$. We also update $X_g = \max(X_g,dz)$ for owner region $a_X$.
\item If $a'_p \notin a_X$ but $a_p \in a_X$ (i.e. we have just moved outside the owner region), set $a'_d = 1000$ and $a'_{p_0} = a_p$. Thus, for points outside the original region, distance and direction of travel is measured relative to the boundary of the region.
\item When $a'_p \notin a_X$ (outside of the original owner region), the strength $a'_s$ may be reduced by low slope $dz$ relative to previously recorded high values:
\begin{itemize}
\item if $dz < 0.6\cdot a'_g$ then $a'_s = a'_s \cdot dz / a'_g$
\item For owner region $a'_X$ of $a'$, if $dz < 0.2\cdot X_g$ then $a'_s = a'_s \cdot 3 dz / X_g$
\end{itemize}
\end{itemize}
\item So long as $A' \neq \emptyset$, repeat the preceding step with $A = A'$.
\end{enumerate}
\item Fill in any holes in each region $X$ so that each region is simply connected
\item Eliminate any regions which cover an area of less than $\frac{\pi r^2}{4}$ or more than $8\pi (r+5)^2$
\item Check the characteristics of each region $X$ to see whether it is indeed a crater. The crater identification process is based on 7 somewhat arbitrary tuning parameters $k_1,\cdots,k_7$ (described below), for which the optimal values are a matter for future research.
\begin{itemize}
\item find the center-of-mass of the region $\vc{c} = \frac{\iint_X \vc{x} d\vc{x}}{\iint_X d\vc{x}}$
\item Express the shape of $X$ in terms of the radius around $\vc{c}$ as a fourier series:
\[r(\theta) = r_0\cdot \left( 1 + \sum_{n=1}^\infty a_n \sin n\theta + b_n \cos n\theta \right) \]
\[ \pi r_0^2 \equiv \iint_X d\vc{x};\ \ a_n \equiv \frac{\iint_X \sin n\theta\,d\vc{x}}{\iint_X d\vc{x}};\ \ b_n \equiv \frac{\iint_X \cos n\theta\,d\vc{x}}{\iint_X d\vc{x}} \]
\item Let $r_i(\theta)$ denote the $i^{th}$ partial sum of the series, $r_i(\theta) = r_0 \cdot \left( 1 + \sum_{n=1}^i a_n \sin n\theta + b_n \cos n\theta \right)$.
\item Define $m_i \equiv \sqrt{a_i^2 + b_i^2}\ \forall i > 0$, the magnitude of the contribution from the $i^{th}$ component. For an ideal (perfectly round) crater, $m_i = 0 \ \forall i > 0$.
\item Eliminate as non-craters any regions for which $m_2 > k_1 = 0.25$ (region is too narrow; ideally 0)
\item Eliminate as non-craters any regions for which $m_3 > k_2 = 0.10$ (region is too lumpy; ideally 0)
\item To assess how well the first $n$ terms of the series reflect the shape of the region, we can measure how far inside or outside of the $n$-term sum boundary the actual region lies. With $X$ denoting, as usual, the actual crater region, let $X_n$ denote the region defined by the first $n$ terms of the radial fourier series. The "error region" is thus $E = (X\setminus X_n) \cup (X_n \setminus X)$. Define
\[d_n \equiv \frac{1}{\pi r_0^2} \sqrt{ \iint_{E} |\vc{x} - r_n(\theta)\vch{x}|^2 d\vc{x} } \]
which provides a (scale invariant) measure of deviation from the $n$-term sum (with perfect agreement indicated by $d_n = 0$).
\item Eliminate as non-craters any regions for which $d_2 > k_3 = 0.06$ (ideally 0)
\item Eliminate as non-craters any regions for which $d_3 > k_4 = 0.02$ (ideally 0)
\item Express the the average gradient $\vc{g}$ as a function of angle $\theta$ around the center point $\vc{c}$ in terms of a Fourier series
\[ \vc{g}(\theta) = \sum_{n=0}^\infty \left( a_n \sin n\theta + b_n \cos n\theta \right) \vch{x} +\sum_{n=0}^\infty \left( c_n \sin n\theta + d_n \cos n\theta \right) \vch{y} \]
\[ a_n \equiv \frac{\iint_X \frac{\partial}{\partial x}z(\vc{x}) \cdot \sin n\theta\,d\vc{x}}{\iint_X d\vc{x}},\ \mathrm{etc.}\]
For an ideal crater, $b_1 \approx c_1$ are the only nonzero components. $b_0$ and $d_0$ indicate the level of average background slope, and can mostly be ignored.
\item Eliminate as non-craters any regions for which $\min(\frac{|b_1|}{|c_1|},\frac{|c_1|}{|b_1|}) < k_5 = 0.5$ (balance between $x$ and $y$ direction slopes, ideally 1)
\item define $s_1 = \sqrt{b_1^2+c_1^2}$, $s_1 ' = \sqrt{a_1^2+d_1^2}$, $s_2 = \sqrt{a_2^2+b_2^2+c_2^2+d_2^2}$
\item Eliminate as non-craters any regions for which $\frac{s_1'}{s_1} > k_6 = 0.33$ (ideally 0)
\item Eliminate as non-craters any regions for which $\frac{s_2}{s_1} > k_7 = 0.33$ (ideally 0)
\end{itemize}
\item If $r<20$, restart the craterfinding process with new target radius $2r$. Otherwise, so long as both the image width and height of $\mtx{Z}$ are $>40$ pixels, scale the image $\mtx{Z}$ down by a factor of 2 in each direction to produce reduced image $\mtx{Z'}$, and apply the same craterfinding routine to $\mtx{Z'}$ with target radius $r=20$.
\end{enumerate}
\end{document}