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

SCM Repository

[diderot] Annotation of /trunk/doc/probe/code/c3dgrad.c
ViewVC logotype

Annotation of /trunk/doc/probe/code/c3dgrad.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : glk 296
2 :     #include <stdio.h>
3 :     #include <math.h>
4 :    
5 :     #define IMG_SZ 11
6 :    
7 :     double V[IMG_SZ][IMG_SZ][IMG_SZ]; // image data
8 :     double h[4][4] = { // bspln3
9 :     { 1.33333, 2.0, 1.0, 0.166667 }, // -2 .. -1
10 :     { 0.666667, 0.0, -1.0, -0.5 }, // -1 .. 0
11 :     { 0.666667, 0.0, -1.0, 0.5 }, // 0 .. 1
12 :     { 1.33333, -2.0, 1.0, -0.166667 }, // 1 .. 2
13 :     };
14 :     const int s = 2; // kernel support
15 :    
16 :     typedef double vec3[3];
17 :    
18 :     void
19 :     probe (double grad[3], const double x[3])
20 :     {
21 :     /* double x[3] = transform(p); // image-space position */
22 :     double nxd, nyd, nzd, fx, fy, fz;
23 :     int nx, ny, nz;
24 :    
25 :     fx = modf (x[0], &nxd); nx = (int)nxd;
26 :     fy = modf (x[1], &nyd); ny = (int)nyd;
27 :     fz = modf (x[2], &nzd); nz = (int)nzd;
28 :    
29 :     // compute kernel values for each axis
30 :     double hx[4], dhx[4], hy[4], dhy[4], hz[4], dhz[4];
31 :     for (int i = 1-s; i <= s; i++) {
32 :     double t;
33 :     t = fx - i;
34 :     hx[i+s-1] = h[s-i][0] + t*(h[s-i][1] + t*(h[s-i][2] + t*h[s-i][3]));
35 :     dhx[i+s-1] = h[s-i][1] + t*(2*h[s-i][2] + t*3*h[s-i][3]);
36 :     t = fy - i;
37 :     hy[i+s-1] = h[s-i][0] + t*(h[s-i][1] + t*(h[s-i][2] + t*h[s-i][3]));
38 :     dhy[i+s-1] = h[s-i][1] + t*(2*h[s-i][2] + t*3*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 :     dhz[i+s-1] = h[s-i][1] + t*(2*h[s-i][2] + t*3*h[s-i][3]);
42 :     }
43 :    
44 :     vec3 vx = {0, 0, 0};
45 :     for (int i = 1-s; i <= s; i++) {
46 :     vec3 vy = {0, 0, 0};
47 :     for (int j = 1-s; j <= s; j++) {
48 :     vec3 vz = {0, 0, 0};
49 :     for (int k = 1-s; k <= s; k++) {
50 :     double v = V[nx+i][ny+j][nz+k];
51 :     vz[0] += v * hz[k+s-1];
52 :     vz[1] += v * hz[k+s-1];
53 :     vz[2] += v * dhz[k+s-1];
54 :     }
55 :     vy[0] += vz[0] * hy[j+s-1];
56 :     vy[1] += vz[1] * dhy[j+s-1];
57 :     vy[2] += vz[2] * hy[j+s-1];
58 :     }
59 :     vx[0] += vy[0] * dhx[i+s-1];
60 :     vx[1] += vy[1] * hx[i+s-1];
61 :     vx[2] += vy[2] * hx[i+s-1];
62 :     }
63 :    
64 :     grad[0] = vx[0];
65 :     grad[1] = vx[1];
66 :     grad[2] = vx[2];
67 :     return;
68 :     }
69 :    
70 :     #define AIR_AFFINE(i,x,I,o,O) ( \
71 :     ((double)(O)-(o))*((double)(x)-(i)) / ((double)(I)-(i)) + (o))
72 :    
73 :     int
74 :     main() {
75 :     int I, J, K, N;
76 :     double P[3], grad[3];
77 :    
78 :     V[5][5][5] = 1.0;
79 :     N=60;
80 :     for (I=0; I<N; I++) {
81 :     P[0] = AIR_AFFINE(0, I, N, 2.0001, 7.9999);
82 :     for (J=0; J<N; J++) {
83 :     P[1] = AIR_AFFINE(0, J, N, 2.0001, 7.9999);
84 :     for (K=0; K<N; K++) {
85 :     P[2] = AIR_AFFINE(0, K, N, 2.0001, 7.9999);
86 :     probe(grad, P);
87 :     printf("%g %g %g\n", grad[0], grad[1], grad[2]);
88 :     }
89 :     }
90 :     }
91 :    
92 :     return 0;
93 :     }

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