Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /trunk/doc/probe/paper.tex
ViewVC logotype

View of /trunk/doc/probe/paper.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 156 - (download) (as text) (annotate)
Mon Jul 12 02:25:43 2010 UTC (9 years, 7 months ago) by jhr
File size: 5676 byte(s)
  Adding notes on code for probing a field.
\documentclass[11pt]{article}

\input{defs}
\usepackage{graphicx}
\usepackage{listings}
\lstset{
  basicstyle=\ttfamily,
  keywordstyle=\bfseries,
  showstringspaces=false,
}

\lstset{
  language=C,
}

\setlength{\textwidth}{6in}
\setlength{\oddsidemargin}{0.25in}
\setlength{\evensidemargin}{0.25in}
\setlength{\parskip}{5pt}

\newcommand{\matM}{\mathbf{M}}
\newcommand{\vecp}{\mathbf{p}}
\newcommand{\vecn}{\mathbf{n}}
\newcommand{\vecf}{\mathbf{f}}
\newcommand{\VEC}[1]{\left\langle{#1}\right\rangle}
\newcommand{\FLOOR}[1]{\left\lfloor{#1}\right\rfloor}

\title{Compiling probe operations for Diderot}
\author{
 John Reppy \\
  University of Chicago \\
  {\small\tt{}jhr@cs.uchicago.edu} \\
}
\date{\today}

\begin{document}

\maketitle
\thispagestyle{empty}

\bibliographystyle{../common/alpha}
\bibliography{../common/strings-short,../common/manticore}

\section{Introduction}

This note describes the code needed to implement a probe of a field operation.

In the discussion below, we use $\matM^{-1}$ for the homogeneous matrix that maps from world
coordinates to image coordinates and use $\vecp$ for the position vector.

\section{Probing a 1D scalar field}
The simplest case is probing a 1D scalar field $F = V\circledast{}h$, where $s$ is the support
of $h$.
The probe $F\mkw{@}x$ is computed as follows:
\begin{eqnarray*}
  \left[\begin{array}{c} x' \\ 1 \end{array}\right] & = & \matM^{-1} \left[\begin{array}{c} x \\ 1 \end{array}\right] \qquad \text{\textit{transform to image space}} \\
  n  & = & \FLOOR{x'} \qquad \text{\textit{integer part of position}} \\
  f  & = & x' - f \qquad \text{\textit{fractional part of position}} \\
  F\mkw{@}x & = & \sum_{i=1-s}^s {V(n+i) h(f - i)}
\end{eqnarray*}%
The convolution $h$ is represented as a symmetric piecewise polynomial formed
from $s$ polynomials $h_1,\,\ldots,\,h_s$.
\begin{displaymath}
  h(x) = \left\{\begin{array}{ll}
    0 & \text{$x < s$} \\
    h_i(x) & \text{$-s \leq -i \leq x < 1-i \leq 0$} \\
    h_i(x) & \text{$0 < i-1 \leq x < i \leq s$} \\
    0 & \text{$x > s$} \\
  \end{array}\right.
\end{displaymath}%
Thus, we can rewrite the probe computation as
\begin{displaymath}
  F\mkw{@}x = \sum_{i=1-s}^{0} {V(n+i) h_{1-i}(f - i)} + \sum_{i=1}^s {V(n+i) h_i(f - i)}
\end{displaymath}%
\figref{fig:1d-probe} illustrates the situation for a kernel $h$ with support $s = 2$,
and \figref{fig:1d-probe-code} gives the C code for the probe operation, assuming that
$h$ is represented by third-degree polynomials.
\begin{figure}[t]
  \begin{center}
    \includegraphics[scale=0.8]{pictures/convo}
  \end{center}%
  \caption{1D scalar probe}
  \label{fig:1d-probe}
\end{figure}%

\begin{figure}[t]
\begin{quote}
\begin{lstlisting}
double img[];			// image data
double h1[4], h2[4];		// kernel

double probe (double x)
{
  double imgX = transform(x);	// image-space position
  double n, f;	

  f = modf (imgPos, &n);

  double value = 0.0, x;

  w = f + 1.0;
  value += img[n-1] * (h2[0] + w*(h2[1] + w*(h2[2] + w*h2[3])));
  w = f;
  value += img[n] * (h1[0] + w*(h1[1] + w*(h1[2] + w*h1[3])));
  w = f - 1.0;
  value += img[n+1] * (h1[0] + w*(h1[1] + w*(h1[2] + w*h1[3])));
  w = f - 2.0;
  value += img[n+2] * (h2[0] + w*(h2[1] + w*(h2[2] + w*h2[3])));

  return value;
}
\end{lstlisting}%
\end{quote}%
\caption{Computing $F\mkw{@}x$ for 1D scalar field in C}
\label{fig:1d-probe-code}
\end{figure}

\begin{figure}[t]
\begin{quote}
\begin{lstlisting}
double4 d = (double4)(h2[0], h1[0], h1[0], h2[0]); // x^0 coeffs
double4 c = (double4)(h2[1], h1[1], h1[1], h2[1]); // x^1 coeffs
double4 b = (double4)(h2[2], h1[2], h1[2], h2[2]); // x^2 coeffs
double4 a = (double4)(h2[3], h1[3], h1[3], h2[3]); // x^3 coeffs
double4 offsets = (double4)(1.0, 0.0, -1.0, -2.0);

double probe (double x)
{
  double imgX = transform(x);	// image-space position
  double n, f;	

  f = modf (imgPos, &n);

  double4 v = (double4)(img[n-1], img[n], img[n+1], img[n+2]);
  double4 w = (double4)(f) + offsets;
  return dot(v, d + w*(c + w*(b + w*a)));
}
\end{lstlisting}%
\end{quote}%
\caption{Computing $F\mkw{@}x$ for 1D scalar field in OpenCL}
\label{fig:1d-probe-code-opencl}
\end{figure}

\section{Probing a 3D scalar field}

The more common case is when the field is a convolution of a scalar 3-dimensional
field ($F = V\circledast{}h$).
Let  $s$ be the support of $h$.
Then the probe $F\mkw{@}\vecp$ is computed as follows:
\begin{eqnarray*}
  \vecp' & = & \matM^{-1} \vecp \qquad \text{\textit{transform to image space}} \\
  \vecn  & = & \FLOOR{\vecp'} \qquad \text{\textit{integer part of position}} \\
  \vecf  & = & \vecp' - \vecn \qquad \text{\textit{fractional part of position}} \\
  F\mkw{@}\vecp & = & \sum_{i=1-s}^s {\sum_{j=1-s}^s {\sum_{k=1-s}^s {V(\vecn+\VEC{i,j,k}) h(\vecf_x - i) h(\vecf_y - j) h(\vecf_z - k)}}}
\end{eqnarray*}%

\begin{figure}[t]
\begin{quote}
\begin{lstlisting}
typedef double vec3[3];

typedef struct {
    int		degree;
    double	coeff[];
} polynomial;

typedef struct {
    int		support;
    polynomial	*segments[];
} kernel;

double probe (vec3 ***img, kernel *h, vec3 pos)
{
    vec3	imgPos, u, v;
    
    MultMv4 (Minv, pos, imgPos);
    
    for (int i = 0;  i < 3;  i++) v[i] = modf(imgPos[i], &(u[i]));
    
    double value = 0.0;
    int lo = 1 - h->support;
    int hi = h->support;
    for (int i = lo;  i < hi;  i++) {
    	for (int j = lo;  j < hi;  j++) {
    	    for (int k = lo;  k < hi;  k++) {
	        value += img[(int)u[0]+i][(int)u[1]+j][(int)u[2]+k]
	        	* eval(h->segments[abs(i)], v[0]-i)
			* eval(h->segments[abs(j)], v[1]-j)
			* eval(h->segments[abs(k)], v[2]-k);
	    }
	}
    }

}
\end{lstlisting}%
\end{quote}%
\end{figure}

\end{document}

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0