%!TEX root = paper.tex
%

\section{A DSL for image analysis}
\label{sec:design}

One of the major contributions of this research will be the design of a domain-specific
language (DSL), which we call \lang{}, for image-analysis algorithms.
The goals of this design are as follows:
\begin{enumerate}
  \item
    Abstract away from the discrete nature of the image data and present the programmer
    with a continuous model of the image.
  \item
    Support the programming of image-analysis algorithms using the mathematical
    properties of the continuous data and its derived attributes.
  \item
    Expose the inherent parallelism of the algorithms, while abstracting away from the
    details of how computations are mapped onto parallel hardware.
\end{enumerate}%

In the remainder of this section, we describe our preliminary thoughts about what the programming
model of this language will be.
It is important to understand that this discussion represents very preliminary ideas
about the design space.
In particular, there is a tradeoff between the flexibility provided by lower-level
general-purpose language constructs and the expressiveness advantages of more abstract
higher-level mechanisms.
In the discussion below, we have biased the design toward more general-purpose
constructs, but with the expectation that the design will get higher level as
it evolves.

\subsection{Diderot by example}

\begin{figure}[t]
  \begin{quote}
    \input{vr-phong-a}
  \end{quote}%
\end{figure}%
\begin{figure}[p]
  \begin{quote}
    \input{vr-phong-b}
  \end{quote}%
\end{figure}%

\subsection{Design sketch}
In this section we sketch a preliminary design for \lang{} using the examples from
the previous section as motivation.

A \lang{} program has three logical parts.
The first are global declarations of the datasets and fields that the
analysis algorithm computes over, as well as other global variables that
control the analysis.
The datasets and global variables (but not the derived fields) are the
inputs to the program.
The second part defines \emph{actors}, which are independent computational
agents that define the algorithm.
The third part specifies the coordination of the actors --- how they interact and
what the global termination conditions are.

\subsubsection{Types}
\lang{} is a statically-typed language.
It has the usual scalar types, including booleans, integer types of various sizes,
and single and double-precision floating-point values.
Similar to graphics-shader languages~\cite{opengl-shading:3ed} and
OpenCL~\cite{opencl-spec-v1}, it also includes small-dimension vector and matrix types.
In addition to these basic computational types, \lang{} also includes \emph{datasets},
\emph{fields}, and \emph{actor} types, which are discussed below.
Our initial design does not include user-defined data structures, but we expect to
add them at some point in the language evolution.

\subsubsection{Datasets}
An image-analysis program begins with the image data that is to be analyzed.
In \lang{}, the image data is represented by a dataset, which is a
multi-dimensional array of scalars, vectors, or tensors.
These data sets come from scientific or medical imaging instruments.
We provide a simple form for declaring a data set.
For example, the following declares \texttt{data} to be a three-dimensional
array of floating-point values.
\begin{quote}\begin{lstlisting}
dataset float[][][] data;
\end{lstlisting}\end{quote}%
Datasets also come with their own coordinate system, which defines how
integer indices map to world coordinates.

Datasets are not limited to arrays of scalars; they may also have vector or
tensor elements.
On disk, datasets are stored using the NRRD (Nearly Raw Raster Data) format~\cite{nrrd},
which supports multidimensional arrays of multidimensional data, such as those produced
by modern imaging instruments.

\subsubsection{Fields and kernels}
As discussed in \secref{sec:image}, we are interested in image analysis algorithms that
work with the continuous field that is reconstructed from discrete data.
Our DSL supports these algorithms by including continuous fields as a primitive notion.
Fields can be defined by applying a convolution kernel to a dataset.
For example, here we define a three-dimensional field generated from the dataset \text{data}
using a convolution kernel named ``\texttt{cubic}''
\begin{quote}\begin{lstlisting}
field F = cubic (data);
\end{lstlisting}\end{quote}%
The dimension of the field is inferred from the dimension of the dataset.
In our initial design, we expect to provide a predefined set of convolution kernels
(essentially those already supported by Teem~\cite{teem}).

The main operation of fields is the \emph{probe}, which extracts the value at some
location.
This operation is similar to indexing an array, except that the indices are
floating-point values.
For example,
\begin{quote}\begin{lstlisting}
F@(x, y, z)
\end{lstlisting}\end{quote}%
evaluates to the scalar value of the field \texttt{F} at the coordinate $(\mathtt{x},\mathtt{y},\mathtt{z})$.
We can also define new fields by applying operators to them.
For example,
\begin{quote}\begin{lstlisting}
(D (D F)) @ v
\end{lstlisting}\end{quote}%
probes the Hessian field derived by taking the second derivative of \texttt{F} at the
location specified by the vector \texttt{v}.
An important feature of the design is that operations, such as $\mathrm{FA}$
(fractional anisotropy of a tensor), are treated as
higher-order operators that map fields to fields (\ie{}, a field of tensors to a field of
scalars).
% implies other fields:
%	D F  --- gradient of F (has type float[3])
%	D (D F) --- Hessian of F (has type float[3][3])
%	
% QUESTION: is the domain of a field always in the interval [0..1] (for each dimension)
% or should we allow the programmer to specify the domain?
% ANSWER: domain is part of nrrd file format.

\subsection{Other global definitions}
In addition to datasets and fields, a \lang{} program may have other global
definitions.
These include auxiliary functions, such as the color composition and transfer
functions for direct-volume rendering (\secref{sec:dvr}),
and program inputs, which are denoted using the \kw{in} keyword:
\begin{quote}\begin{lstlisting}
in float timeStep;  // simulation timestep
\end{lstlisting}\end{quote}%
Other examples of globals include constants, such as the matrix that maps
pixel positions to image-plane positions in direct-volume rendering.
Note that global variables are immutable; \lang{} does not have an explicit notion
of global state.

\subsubsection{Actors}
Actors are an abstraction of the computations that the image analysis performs.
Examples include the rays in direct-volume rendering (\secref{sec:dvr}),
the fibers in fiber tractography (\secref{sec:fibers}), and the
particles in a particle system (\secref{sec:particles}).
A given program may have multiple types of actors, but we expect that most algorithms
will involve only a single actor definition.

An actor definition is similar to a class declaration in an object-oriented
language, and includes declarations of its local state, initialization,
and an \kw{update} method.
The \kw{update} methods of actors are written using a simple C-style language that is
similar to the graphics shader language GLSL~\cite{opengl-shading:3ed}.
For example, the direct volume rendering algorithm (Algorithm~\ref{alg:volrend})
is implemented as the following actor definition:
\begin{quote}\begin{lstlisting}
actor Ray (vec2i p)
{
  // map pixel location to 3D position in image plane
    vec3f pos = ToImagePlane * ((float)p.x, (float)p.y, 1.0);
  // ray direction
    vec3f dir = normalize(pos - eyePosition);
    float t = 0.0;
    vec4f color = initialColor;

    update () {
        vec3f newPos = pos + t*dir;
        if (newPos.z < farClip) {
            color = compose(color, transferFn(F@newPos));
            if (color.a == 1.0)
                stable;  // color is opaque, so we are done
            else
                t += timeStep;
        }
        else
            stable;
    }
}
\end{lstlisting}\end{quote}%
In this example, the \emph{instance variables} \texttt{pos}, \texttt{dir},
\texttt{t}, and \texttt{color} comprise the local state of an actor.
The \kw{update} method evolves this state until the actor reaches an
\emph{stable} state (denoted by executing the \kw{stable} statement).
An actor may also declare itself dead using the \kw{die} statement and actors may
create new actors using the \kw{new} statement.

For some algorithms, such as those described in \secref{sec:particles}, it is
necessary for actors to examine the state of nearby actors.
This behavior has two aspects: the definition of an actor's \emph{neighborhood},
which is defined declaratively in the global-coordination part of the program (see below),
and the computation on neighbors.
An actor's neighborhood can be represented as array of actors that are passed
as a parameter to the \kw{update} method.
The \kw{update} method can query, but not modify, the instance variables of these actors.
% actor states: active, stable, dead
% actors can create new actors

While the actor \kw{update} methods have some similarity to shader-language programs
and OpenCL kernels, there are some important differences.
First, actors compute over the continuous fields defined in the program\footnote{
  Shader programs treat texture data as continuous, but the interpolation mechanisms
  used in graphics applications are not suitable for image-analysis problems.
} instead of discrete data.
Second, actors can interact with other actors.
Third, actors can die or become stable.
Lastly, actors can create new actors.

% note: system could have multiple kinds of actors?

% QUESTION: how are actor interactions specified?

\subsubsection{Global coordination}
The last part of a \lang{} specification is the definition of the global properties of the
algorithm.
These include the initial set of actors, which for the direct volume rendering
application might be defined using a set-comprehension notation
\begin{quote}\begin{lstlisting}
initially { { Ray (vec2i(i, j)) | i in 0..1023 } | j in 0..1023 };
\end{lstlisting}\end{quote}%
The global properties also include global termination conditions,
such a limit on the number of iterations or a predicate defined on a
global property, such as global energy, and a mechanism for reading out
the result of the computation from the local actor states.
Lastly, global coordination includes properties, such as the neighbors
of a particle, that define inter-agent interactions.

\subsubsection{Semantics}

The semantics of a \lang{} program is based on an iterative model that is similar
to bulk-synchronous parallelism.
Informally, the state of an actor is represented by a triple: $\langle{}a, \sigma, S\rangle{}$,
where $a$ is a unique actor name, $\sigma$ is a mapping from instance variables to values,
and $S$ is the actor's status, which is either $\mathbf{ACTIVE}$ or $\mathbf{STABLE}$.
The state of a program $\mathcal{P}$ is then defined to be a set of actor states.
The \kw{update} method is then modeled as a curried function that takes the program state
and an actor name and produces a set of actor states.
The \kw{update} method returns a set, since the actor may die, in which case the result is the
empty set, or may create new actors, which will also be included in the resulting set.
A single iteration of the program can then be defined by
\begin{displaymath}
  \mathcal{P}_{i+1}
    = \mathrm{Stable}(\mathcal{P}_i) \cup \left( \bigcup_{\langle{}a,\sigma,S\rangle{} \in \mathrm{Active}(\mathcal{P}_i)}{U\,\mathcal{P}_i\,a} \right)
\end{displaymath}%
where $U$ is the function denoted by the \kw{update} method and
\begin{eqnarray*}
  \mathrm{Active}(\mathcal{P}) & = &
    \{ \langle{}a,\sigma, \mathbf{ACTIVE}\rangle{} | \langle{}a, \sigma, \mathbf{ACTIVE}\rangle \in \mathcal{P} \} \\
  \mathrm{Stable}(\mathcal{P}) & = &
    \{ \langle{}a,\sigma, \mathbf{STABLE}\rangle{} | \langle{}a, \sigma, \mathbf{STABLE}\rangle \in \mathcal{P} \}
\end{eqnarray*}%
Note that this rule has the effect of discarding the dead actors.
A program has reached a terminal state when there are no active actors left.

Global coordination is modeled by a program state to program state function
that is applied between iterations.
At this time, we also check for global termination conditions.

Because the \kw{update} method is given the program state $\mathcal{P}_i$ to
compute the actor's $i+1$ state, each actor's update is independent of the others (even if
they depend on a neighboring actor's state).
This property makes the program amenable for efficient parallel execution.