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

SCM Repository

[diderot] Annotation of /trunk/doc/probe/code/c3d-accumf.c
ViewVC logotype

Annotation of /trunk/doc/probe/code/c3d-accumf.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : jhr 291 /* c3d-accumf.c
2 :     *
3 :     * COPYRIGHT (c) 2010 The Diderot Project (http://diderot.cs.uchicago.edu)
4 :     * All rights reserved.
5 :     *
6 :     * Single-precision accumulator version.
7 :     */
8 :    
9 :     #include <stdio.h>
10 :     #include <math.h>
11 :    
12 :     #define IMG_SZ 11
13 :    
14 :     float V[IMG_SZ][IMG_SZ][IMG_SZ]; // image data
15 :     float h[4][4] = { // bspln3
16 :     { 1.33333, 2.0, 1.0, 0.166667 }, // -2 .. -1
17 :     { 0.666667, 0.0, -1.0, -0.5 }, // -1 .. 0
18 :     { 0.666667, 0.0, -1.0, 0.5 }, // 0 .. 1
19 :     { 1.33333, -2.0, 1.0, -0.166667 }, // 1 .. 2
20 :     };
21 :     const int s = 2; // kernel support
22 :    
23 :     float probe (float x[3]) {
24 :     float nxd, nyd, nzd, fx, fy, fz;
25 :     int nx, ny, nz;
26 :    
27 :     fx = modff (x[0], &nxd); nx = (int)nxd;
28 :     fy = modff (x[1], &nyd); ny = (int)nyd;
29 :     fz = modff (x[2], &nzd); nz = (int)nzd;
30 :    
31 :     // compute kernel values for each axis
32 :     float hx[4], hy[4], hz[4];
33 :     for (int i = 1-s; i <= s; i++) {
34 :     float t;
35 :     t = fx - i;
36 :     hx[i+s-1] = h[s-i][0] + t*(h[s-i][1] + t*(h[s-i][2] + t*h[s-i][3]));
37 :     t = fy - i;
38 :     hy[i+s-1] = h[s-i][0] + t*(h[s-i][1] + t*(h[s-i][2] + t*h[s-i][3]));
39 :     t = fz - i;
40 :     hz[i+s-1] = h[s-i][0] + t*(h[s-i][1] + t*(h[s-i][2] + t*h[s-i][3]));
41 :     }
42 :    
43 :     float vx = 0.0;
44 :     for (int i = 1-s; i <= s; i++) {
45 :     float vy = 0.0;
46 :     for (int j = 1-s; j <= s; j++) {
47 :     float vz = 0.0;
48 :     for (int k = 1-s; k <= s; k++) {
49 :     vz += V[nx+i][ny+j][nz+k] * hz[k+s-1];
50 :     }
51 :     vy += vz * hy[j+s-1];
52 :     }
53 :     vx += vy * hx[i+s-1];
54 :     }
55 :    
56 :     return vx;
57 :     }
58 :    
59 :     #define AIR_AFFINE(i,x,I,o,O) ( \
60 :     ((float)(O)-(o))*((float)(x)-(i)) / ((float)(I)-(i)) + (o))
61 :    
62 :     int
63 :     main() {
64 :     int I, J, K, N;
65 :     float P[3];
66 :    
67 :     V[5][5][5] = 1.0;
68 :     N=60;
69 :     for (I=0; I<N; I++) {
70 :     P[0] = AIR_AFFINE(0, I, N, 2.0001, 7.9999);
71 :     for (J=0; J<N; J++) {
72 :     P[1] = AIR_AFFINE(0, J, N, 2.0001, 7.9999);
73 :     for (K=0; K<N; K++) {
74 :     P[2] = AIR_AFFINE(0, K, N, 2.0001, 7.9999);
75 :     printf("%g\n", probe(P));
76 :     }
77 :     }
78 :     }
79 :    
80 :     return 0;
81 :     }

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