Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /trunk/doc/probe/code/c3d-accumf.c
 [diderot] / trunk / doc / probe / code / c3d-accumf.c

# View of /trunk/doc/probe/code/c3d-accumf.c

Sun Aug 15 18:09:49 2010 UTC (11 years, 10 months ago) by jhr
File size: 1928 byte(s)
```  Working on debugging
```
```/* c3d-accumf.c
*
* COPYRIGHT (c) 2010 The Diderot Project (http://diderot.cs.uchicago.edu)
*
* Single-precision accumulator version.
*/

#include <stdio.h>
#include <math.h>

#define IMG_SZ 11

float V[IMG_SZ][IMG_SZ][IMG_SZ];		// image data
float h[4][4] = {		// bspln3
{ 1.33333,   2.0,  1.0,  0.166667 },   // -2 .. -1
{ 0.666667,  0.0, -1.0, -0.5 },        // -1 .. 0
{ 0.666667,  0.0, -1.0,  0.5 },        //  0 .. 1
{ 1.33333,  -2.0,  1.0, -0.166667 },   //  1 .. 2
};
const int s = 2;       // kernel support

float probe (float x[3]) {
float nxd, nyd, nzd, fx, fy, fz;
int nx, ny, nz;

fx = modff (x[0], &nxd); nx = (int)nxd;
fy = modff (x[1], &nyd); ny = (int)nyd;
fz = modff (x[2], &nzd); nz = (int)nzd;

// compute kernel values for each axis
float hx[4], hy[4], hz[4];
for (int i = 1-s;  i <= s;  i++) {
float t;
t = fx - i;
hx[i+s-1] = h[s-i][0] + t*(h[s-i][1] + t*(h[s-i][2] + t*h[s-i][3]));
t = fy - i;
hy[i+s-1] = h[s-i][0] + t*(h[s-i][1] + t*(h[s-i][2] + t*h[s-i][3]));
t = fz - i;
hz[i+s-1] = h[s-i][0] + t*(h[s-i][1] + t*(h[s-i][2] + t*h[s-i][3]));
}

float vx = 0.0;
for (int i = 1-s;  i <= s;  i++) {
float vy = 0.0;
for (int j = 1-s;  j <= s;  j++) {
float vz = 0.0;
for (int k = 1-s;  k <= s;  k++) {
vz += V[nx+i][ny+j][nz+k] * hz[k+s-1];
}
vy += vz * hy[j+s-1];
}
vx += vy * hx[i+s-1];
}

return vx;
}

#define AIR_AFFINE(i,x,I,o,O) ( \
((float)(O)-(o))*((float)(x)-(i)) / ((float)(I)-(i)) + (o))

int
main() {
int I, J, K, N;
float P[3];

V[5][5][5] = 1.0;
N=60;
for (I=0; I<N; I++) {
P[0] = AIR_AFFINE(0, I, N, 2.0001, 7.9999);
for (J=0; J<N; J++) {
P[1] = AIR_AFFINE(0, J, N, 2.0001, 7.9999);
for (K=0; K<N; K++) {
P[2] = AIR_AFFINE(0, K, N, 2.0001, 7.9999);
printf("%g\n", probe(P));
}
}
}

return 0;
}
```