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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 296 - (download) (as text) (annotate)
Sun Aug 15 23:06:00 2010 UTC (8 years, 11 months ago) by glk
File size: 2446 byte(s)
for gradient code
#include <stdio.h>
#include <math.h>

#define IMG_SZ 11

double V[IMG_SZ][IMG_SZ][IMG_SZ];		// image data
double 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

typedef double vec3[3];

void
probe (double grad[3], const double x[3])
{
  /* double x[3] = transform(p);	// image-space position */
  double nxd, nyd, nzd, fx, fy, fz;
  int nx, ny, nz;

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

  // compute kernel values for each axis
  double hx[4], dhx[4], hy[4], dhy[4], hz[4], dhz[4];
  for (int i = 1-s;  i <= s;  i++) {
    double 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]));
    dhx[i+s-1] = h[s-i][1] + t*(2*h[s-i][2] + t*3*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]));
    dhy[i+s-1] = h[s-i][1] + t*(2*h[s-i][2] + t*3*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]));
    dhz[i+s-1] = h[s-i][1] + t*(2*h[s-i][2] + t*3*h[s-i][3]);
  }

  vec3 vx = {0, 0, 0};
  for (int i = 1-s;  i <= s;  i++) {
    vec3 vy = {0, 0, 0};
    for (int j = 1-s;  j <= s;  j++) {
      vec3 vz = {0, 0, 0};
      for (int k = 1-s;  k <= s;  k++) {
        double v = V[nx+i][ny+j][nz+k];
        vz[0] += v * hz[k+s-1];
        vz[1] += v * hz[k+s-1];
        vz[2] += v * dhz[k+s-1];
      }
      vy[0] += vz[0] * hy[j+s-1];
      vy[1] += vz[1] * dhy[j+s-1];
      vy[2] += vz[2] * hy[j+s-1];
    }
    vx[0] += vy[0] * dhx[i+s-1];
    vx[1] += vy[1] * hx[i+s-1];
    vx[2] += vy[2] * hx[i+s-1];
  }

  grad[0] = vx[0];
  grad[1] = vx[1];
  grad[2] = vx[2];
  return;
}

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

int
main() {
  int I, J, K, N;
  double P[3], grad[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);
        probe(grad, P);
        printf("%g %g %g\n", grad[0], grad[1], grad[2]);
      }
    }
  }

  return 0;
}

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