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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : jhr 291 /* c3d-vec.c
2 :     *
3 :     * COPYRIGHT (c) 2010 The Diderot Project (http://diderot.cs.uchicago.edu)
4 :     * All rights reserved.
5 :     *
6 :     * Single-precision vector version.
7 :     */
8 :    
9 :     #include <stdio.h>
10 :     #include <math.h>
11 :     #include <assert.h>
12 :    
13 :     typedef float float4 __attribute__ ((__vector_size__ (16)));
14 :    
15 :     typedef union {
16 :     float f[4];
17 :     float4 v;
18 :     } f4u_t;
19 :    
20 :     #define IMG_SZ 11
21 :    
22 :     float V[IMG_SZ][IMG_SZ][IMG_SZ]; // image data
23 :     float h[4][4] = { // bspln3
24 :     { 1.33333, 2.0, 1.0, 0.166667 }, // -2 .. -1
25 :     { 0.666667, 0.0, -1.0, -0.5 }, // -1 .. 0
26 :     { 0.666667, 0.0, -1.0, 0.5 }, // 0 .. 1
27 :     { 1.33333, -2.0, 1.0, -0.166667 }, // 1 .. 2
28 :     };
29 :     const int s = 2; // kernel support
30 :    
31 :     float4 a, b, c, d;
32 :    
33 :     // fake vector dot product. There are SSE instructions for computing this.
34 :     //
35 :     static inline float dot (float4 a, float4 b)
36 :     {
37 :     f4u_t av, bv;
38 :     av.v = a;
39 :     bv.v = b;
40 :     float s = 0.0;
41 :     for (int i = 0; i < 4; i++)
42 :     s += av.f[i] * bv.f[i];
43 :    
44 :     return s;
45 :     }
46 :    
47 :     static inline void printv (float4 v)
48 :     {
49 :     f4u_t f;
50 :     f.v = v;
51 :     fprintf(stderr, "<%f, %f, %f, %f>", f.f[0], f.f[1], f.f[2], f.f[3]);
52 :     }
53 :    
54 :     float probe (float x[3]) {
55 :     float nxd, nyd, nzd, fx, fy, fz;
56 :     int nx, ny, nz;
57 :    
58 :     fx = modff (x[0], &nxd); nx = (int)nxd;
59 :     fy = modff (x[1], &nyd); ny = (int)nyd;
60 :     fz = modff (x[2], &nzd); nz = (int)nzd;
61 :    
62 :     // compute kernel values for each axis
63 :     float4 t, hx, hy, hz;
64 :     t = (float4){fx + 1, fx, fx - 1, fx - 2};
65 :     hx = d + t*(c + t*(b + t*a));
66 :     t = (float4){fy + 1, fy, fy - 1, fy - 2};
67 :     hy = d + t*(c + t*(b + t*a));
68 :     t = (float4){fz + 1, fz, fz - 1, fz - 2};
69 :     hz = d + t*(c + t*(b + t*a));
70 :    
71 :     f4u_t vy, vz;
72 :     for (int i = 1-s; i <= s; i++) {
73 :     for (int j = 1-s; j <= s; j++) {
74 :     float4 v = (float4){
75 :     V[nx+i][ny+j][nz-1],
76 :     V[nx+i][ny+j][nz],
77 :     V[nx+i][ny+j][nz+1],
78 :     V[nx+i][ny+j][nz+2]};
79 :     vz.f[j+s-1] = dot(v, hz);
80 :     }
81 :     vy.f[i+s-1] = dot(vz.v, hy);
82 :     }
83 :    
84 :     return dot(vy.v, hx);
85 :     }
86 :    
87 :     #define AIR_AFFINE(i,x,I,o,O) ( \
88 :     ((float)(O)-(o))*((float)(x)-(i)) / ((float)(I)-(i)) + (o))
89 :    
90 :     int
91 :     main() {
92 :     int I, J, K, N;
93 :     float P[3];
94 :    
95 :     // initialize coefficients
96 :     d = (float4){h[3][0], h[2][0], h[1][0], h[0][0]}; // x^0 coeffs
97 :     c = (float4){h[3][1], h[2][1], h[1][1], h[0][1]}; // x^1 coeffs
98 :     b = (float4){h[3][2], h[2][2], h[1][2], h[0][2]}; // x^2 coeffs
99 :     a = (float4){h[3][3], h[2][3], h[1][3], h[0][3]}; // x^3 coeffs
100 :    
101 :     V[5][5][5] = 1.0;
102 :     N=60;
103 :     for (I=0; I<N; I++) {
104 :     P[0] = AIR_AFFINE(0, I, N, 2.0001, 7.9999);
105 :     for (J=0; J<N; J++) {
106 :     P[1] = AIR_AFFINE(0, J, N, 2.0001, 7.9999);
107 :     for (K=0; K<N; K++) {
108 :     P[2] = AIR_AFFINE(0, K, N, 2.0001, 7.9999);
109 :     printf("%g\n", probe(P));
110 :     }
111 :     }
112 :     }
113 :    
114 :     return 0;
115 :     }

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