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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 156 - (view) (download) (as text)

1 : jhr 156 \documentclass[11pt]{article}
2 :    
3 :     \input{defs}
4 :     \usepackage{graphicx}
5 :     \usepackage{listings}
6 :     \lstset{
7 :     basicstyle=\ttfamily,
8 :     keywordstyle=\bfseries,
9 :     showstringspaces=false,
10 :     }
11 :    
12 :     \lstset{
13 :     language=C,
14 :     }
15 :    
16 :     \setlength{\textwidth}{6in}
17 :     \setlength{\oddsidemargin}{0.25in}
18 :     \setlength{\evensidemargin}{0.25in}
19 :     \setlength{\parskip}{5pt}
20 :    
21 :     \newcommand{\matM}{\mathbf{M}}
22 :     \newcommand{\vecp}{\mathbf{p}}
23 :     \newcommand{\vecn}{\mathbf{n}}
24 :     \newcommand{\vecf}{\mathbf{f}}
25 :     \newcommand{\VEC}[1]{\left\langle{#1}\right\rangle}
26 :     \newcommand{\FLOOR}[1]{\left\lfloor{#1}\right\rfloor}
27 :    
28 :     \title{Compiling probe operations for Diderot}
29 :     \author{
30 :     John Reppy \\
31 :     University of Chicago \\
32 :     {\small\tt{}jhr@cs.uchicago.edu} \\
33 :     }
34 :     \date{\today}
35 :    
36 :     \begin{document}
37 :    
38 :     \maketitle
39 :     \thispagestyle{empty}
40 :    
41 :     \bibliographystyle{../common/alpha}
42 :     \bibliography{../common/strings-short,../common/manticore}
43 :    
44 :     \section{Introduction}
45 :    
46 :     This note describes the code needed to implement a probe of a field operation.
47 :    
48 :     In the discussion below, we use $\matM^{-1}$ for the homogeneous matrix that maps from world
49 :     coordinates to image coordinates and use $\vecp$ for the position vector.
50 :    
51 :     \section{Probing a 1D scalar field}
52 :     The simplest case is probing a 1D scalar field $F = V\circledast{}h$, where $s$ is the support
53 :     of $h$.
54 :     The probe $F\mkw{@}x$ is computed as follows:
55 :     \begin{eqnarray*}
56 :     \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}} \\
57 :     n & = & \FLOOR{x'} \qquad \text{\textit{integer part of position}} \\
58 :     f & = & x' - f \qquad \text{\textit{fractional part of position}} \\
59 :     F\mkw{@}x & = & \sum_{i=1-s}^s {V(n+i) h(f - i)}
60 :     \end{eqnarray*}%
61 :     The convolution $h$ is represented as a symmetric piecewise polynomial formed
62 :     from $s$ polynomials $h_1,\,\ldots,\,h_s$.
63 :     \begin{displaymath}
64 :     h(x) = \left\{\begin{array}{ll}
65 :     0 & \text{$x < s$} \\
66 :     h_i(x) & \text{$-s \leq -i \leq x < 1-i \leq 0$} \\
67 :     h_i(x) & \text{$0 < i-1 \leq x < i \leq s$} \\
68 :     0 & \text{$x > s$} \\
69 :     \end{array}\right.
70 :     \end{displaymath}%
71 :     Thus, we can rewrite the probe computation as
72 :     \begin{displaymath}
73 :     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)}
74 :     \end{displaymath}%
75 :     \figref{fig:1d-probe} illustrates the situation for a kernel $h$ with support $s = 2$,
76 :     and \figref{fig:1d-probe-code} gives the C code for the probe operation, assuming that
77 :     $h$ is represented by third-degree polynomials.
78 :     \begin{figure}[t]
79 :     \begin{center}
80 :     \includegraphics[scale=0.8]{pictures/convo}
81 :     \end{center}%
82 :     \caption{1D scalar probe}
83 :     \label{fig:1d-probe}
84 :     \end{figure}%
85 :    
86 :     \begin{figure}[t]
87 :     \begin{quote}
88 :     \begin{lstlisting}
89 :     double img[]; // image data
90 :     double h1[4], h2[4]; // kernel
91 :    
92 :     double probe (double x)
93 :     {
94 :     double imgX = transform(x); // image-space position
95 :     double n, f;
96 :    
97 :     f = modf (imgPos, &n);
98 :    
99 :     double value = 0.0, x;
100 :    
101 :     w = f + 1.0;
102 :     value += img[n-1] * (h2[0] + w*(h2[1] + w*(h2[2] + w*h2[3])));
103 :     w = f;
104 :     value += img[n] * (h1[0] + w*(h1[1] + w*(h1[2] + w*h1[3])));
105 :     w = f - 1.0;
106 :     value += img[n+1] * (h1[0] + w*(h1[1] + w*(h1[2] + w*h1[3])));
107 :     w = f - 2.0;
108 :     value += img[n+2] * (h2[0] + w*(h2[1] + w*(h2[2] + w*h2[3])));
109 :    
110 :     return value;
111 :     }
112 :     \end{lstlisting}%
113 :     \end{quote}%
114 :     \caption{Computing $F\mkw{@}x$ for 1D scalar field in C}
115 :     \label{fig:1d-probe-code}
116 :     \end{figure}
117 :    
118 :     \begin{figure}[t]
119 :     \begin{quote}
120 :     \begin{lstlisting}
121 :     double4 d = (double4)(h2[0], h1[0], h1[0], h2[0]); // x^0 coeffs
122 :     double4 c = (double4)(h2[1], h1[1], h1[1], h2[1]); // x^1 coeffs
123 :     double4 b = (double4)(h2[2], h1[2], h1[2], h2[2]); // x^2 coeffs
124 :     double4 a = (double4)(h2[3], h1[3], h1[3], h2[3]); // x^3 coeffs
125 :     double4 offsets = (double4)(1.0, 0.0, -1.0, -2.0);
126 :    
127 :     double probe (double x)
128 :     {
129 :     double imgX = transform(x); // image-space position
130 :     double n, f;
131 :    
132 :     f = modf (imgPos, &n);
133 :    
134 :     double4 v = (double4)(img[n-1], img[n], img[n+1], img[n+2]);
135 :     double4 w = (double4)(f) + offsets;
136 :     return dot(v, d + w*(c + w*(b + w*a)));
137 :     }
138 :     \end{lstlisting}%
139 :     \end{quote}%
140 :     \caption{Computing $F\mkw{@}x$ for 1D scalar field in OpenCL}
141 :     \label{fig:1d-probe-code-opencl}
142 :     \end{figure}
143 :    
144 :     \section{Probing a 3D scalar field}
145 :    
146 :     The more common case is when the field is a convolution of a scalar 3-dimensional
147 :     field ($F = V\circledast{}h$).
148 :     Let $s$ be the support of $h$.
149 :     Then the probe $F\mkw{@}\vecp$ is computed as follows:
150 :     \begin{eqnarray*}
151 :     \vecp' & = & \matM^{-1} \vecp \qquad \text{\textit{transform to image space}} \\
152 :     \vecn & = & \FLOOR{\vecp'} \qquad \text{\textit{integer part of position}} \\
153 :     \vecf & = & \vecp' - \vecn \qquad \text{\textit{fractional part of position}} \\
154 :     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)}}}
155 :     \end{eqnarray*}%
156 :    
157 :     \begin{figure}[t]
158 :     \begin{quote}
159 :     \begin{lstlisting}
160 :     typedef double vec3[3];
161 :    
162 :     typedef struct {
163 :     int degree;
164 :     double coeff[];
165 :     } polynomial;
166 :    
167 :     typedef struct {
168 :     int support;
169 :     polynomial *segments[];
170 :     } kernel;
171 :    
172 :     double probe (vec3 ***img, kernel *h, vec3 pos)
173 :     {
174 :     vec3 imgPos, u, v;
175 :    
176 :     MultMv4 (Minv, pos, imgPos);
177 :    
178 :     for (int i = 0; i < 3; i++) v[i] = modf(imgPos[i], &(u[i]));
179 :    
180 :     double value = 0.0;
181 :     int lo = 1 - h->support;
182 :     int hi = h->support;
183 :     for (int i = lo; i < hi; i++) {
184 :     for (int j = lo; j < hi; j++) {
185 :     for (int k = lo; k < hi; k++) {
186 :     value += img[(int)u[0]+i][(int)u[1]+j][(int)u[2]+k]
187 :     * eval(h->segments[abs(i)], v[0]-i)
188 :     * eval(h->segments[abs(j)], v[1]-j)
189 :     * eval(h->segments[abs(k)], v[2]-k);
190 :     }
191 :     }
192 :     }
193 :    
194 :     }
195 :     \end{lstlisting}%
196 :     \end{quote}%
197 :     \end{figure}
198 :    
199 :     \end{document}

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