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

# SCM Repository

[diderot] Diff of /benchmarks/programs/vr-lite-cam/bmark-teem.c
 [diderot] / benchmarks / programs / vr-lite-cam / bmark-teem.c

# Diff of /benchmarks/programs/vr-lite-cam/bmark-teem.c

revision 1575, Sat Oct 29 16:50:18 2011 UTC revision 1606, Thu Nov 3 22:51:58 2011 UTC
# Line 1  Line 1
1  /*! \file vr-lite-cam.c  /*! \file bmark-teem.c
2   *   *
3   * \author Nick Seltzer   * \author Nick Seltzer & Gordon Kindlmann
4   */   */
5
6  /*  /*
# Line 14  Line 14
14  #include "teem/nrrd.h"  #include "teem/nrrd.h"
15  #include "teem/gage.h"  #include "teem/gage.h"
16
17  #include <sys/time.h>  #include "teem-defs.h"

static double GetTime ()
{
struct timeval tv;

gettimeofday (&tv, 0);

return (double)tv.tv_sec + 0.000001 * (double)tv.tv_usec;
}

#define STATIC_INLINE static inline
18
19  STATIC_INLINE double dot3(double vec1[3], double vec2[3])  STATIC_INLINE double dot3(double vec1[3], double vec2[3])
20  {  {
# Line 76  Line 65
65
66  STATIC_INLINE bool inside(double pos[3])  STATIC_INLINE bool inside(double pos[3])
67  {  {
68    // XXX - Hack    // XXX - Hack specific to ../../data/vfrhand-nohip.nhdr"
69    if(pos[0] < -0.5 || pos[0] > 175.5) return false;    // but for some reason still not matching what Diderot says
70    if(pos[1] < -0.5 || pos[1] > 186.5) return false;    if(pos[0] < 2*0.9375 || pos[0] > (176-3)*0.9375) return false;
71    if(pos[2] < -0.5 || pos[2] > 189.5) return false;    if(pos[1] < 2*0.9375 || pos[1] > (187-3)*0.9375) return false;
72      if(pos[2] < 2.0      || pos[2] > (190-3)) return false;
73    return true;    return true;
74  }  }
75
# Line 93  Line 83
83    return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));    return (unsigned char)((0.0f * (in_max - in) + 255.0f * (in - in_min)) / (in_max - in_min));
84  }  }
85
86  int main ()  int main (int argc, const char *argv[])
87  {  {
88      const char *me = argv[0];
89    char *infile = "../../data/vfrhand-nohip.nhdr";    char *infile = "../../data/vfrhand-nohip.nhdr";
90    char *outfile = "vr-lite-cam.pgm";    char *outfile = "bmark-teem.nrrd";
double Zero[3] = {0, 0, 0};
91    double camEye[3] = {127.331, -1322.05, 272.53};    double camEye[3] = {127.331, -1322.05, 272.53};
92    double camAt[3] = {63.0, 82.6536, 98.0};    double camAt[3] = {63.0, 82.6536, 98.0};
93    double camUp[3] = {0.9987, 0.0459166, -0.0221267};    double camUp[3] = {0.9987, 0.0459166, -0.0221267};
# Line 106  Line 96
96    double camFOV = 5.0;    double camFOV = 5.0;
97    int imgResU = 480;    int imgResU = 480;
98    int imgResV = 345;    int imgResV = 345;
99    double rayStep = 0.5;    double rayStep = 0.1;
100    double valOpacMin = 400.0;        // 400.0 for skin, 1150.0 for bone    double valOpacMin = 400.0;        // 400.0 for skin, 1150.0 for bone
101    double valOpacMax = 700.0;        // 700.0 for skin, 1450.0 for bone    double valOpacMax = 700.0;        // 700.0 for skin, 1450.0 for bone
102
103
104    double temp[3];    double temp[3];
105    double temp2[3];    double temp2[3];
106    double temp3[3];    double temp3[3];
# Line 123  Line 115
115    norm3(temp, camU);    norm3(temp, camU);
116    double camV[3];    double camV[3];
117    cross(camN, camU, camV);    cross(camN, camU, camV);
118    double camVmax = tan(camFOV*3.1415926536/360.0)*camDist;    double camVmax = tan(camFOV*AIR_PI/360.0)*camDist;
119    double camUmax = camVmax*(double)(imgResU)/(double)(imgResV);    double camUmax = camVmax*(double)(imgResU)/(double)(imgResV);
120    double lightVspDir[3] = {0.9, -1.0, -2.5};    double lightVspDir[3] = {0.9, -1.0, -2.5};
121    scale_vec3(lightVspDir[0], camU, temp);    scale_vec3(lightVspDir[0], camU, temp);
# Line 136  Line 128
128    double phongKa = 0.05;    double phongKa = 0.05;
129    double phongKd = 0.80;    double phongKd = 0.80;
130    double phongKs = 0.20;    double phongKs = 0.20;
131    double phongSp = 30.0;    double phongSp = 50.0;
132
133    // Read in the data    // Read in the data
134    char* err;    char* err;
135    Nrrd *nin;    Nrrd *nin;
136    int status;    int status;
137    nin = nrrdNew();    airArray *mop;
138    if(nin == NULL) {    mop = airMopNew();
err = biffGetDone(NRRD);
fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
free(err);
return -1;
}
139
140      nin = nrrdNew();
141      airMopAdd(mop, nin, (airMopper)nrrdNuke, airMopAlways);
142
143    status = nrrdLoad(nin, infile, NULL);    // Get data for mip
144    if (status) {    if (nrrdLoad(nin, infile, NULL)) {
145      err = biffGetDone(NRRD);      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
146      fprintf(stderr, "Trouble reading \"%s\":\n%s", infile, err);      fprintf(stderr, "%s: Trouble reading \"%s\":\n%s", me, infile, err);
147      free(err);      airMopError(mop);
nrrdNix(nin);
148      return -1;      return -1;
149    }    }
150
// Get data for mip
151    gageContext *ctx;    gageContext *ctx;
152    gagePerVolume *pvl;    gagePerVolume *pvl;
153    const double *val;    const double *val;
154    const double *grad;    const double *grad;
155    double kparm[NRRD_KERNEL_PARMS_NUM];    double kparm[NRRD_KERNEL_PARMS_NUM] = {0.0};
156
157    ctx = gageContextNew();    ctx = gageContextNew();
158    if(ctx == NULL) {    airMopAdd(mop, ctx, (airMopper)gageContextNix, airMopAlways);
159      err = biffGetDone(GAGE);    if (!( pvl = gagePerVolumeNew(ctx, nin, gageKindScl) )) {
160      fprintf(stderr, "Trouble allocating gage context:\n%s", err);      airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
161      free(err);      fprintf(stderr, "%s: couldn't create gage context:\n%s", me, err);
162      nrrdNuke(nin);      airMopError(mop);
163      return -1;      return -1;
164    }    }
165      if (gagePerVolumeAttach(ctx, pvl)
166    pvl = gagePerVolumeNew(ctx, nin, gageKindScl);        || gageKernelSet(ctx, gageKernel00, nrrdKernelBSpline3, kparm)
167    if(pvl == NULL) {        || gageKernelSet(ctx, gageKernel11, nrrdKernelBSpline3D, kparm)
168      err = biffGetDone(GAGE);        || gageQueryItemOn(ctx, pvl, gageSclValue)
169      fprintf(stderr, "Trouble allocating gage PerVolume:\n%s", err);        || gageQueryItemOn(ctx, pvl, gageSclGradVec)
170      free(err);        || gageUpdate(ctx)) {
171      gageContextNix(ctx);      airMopAdd(mop, err = biffGetDone(GAGE), airFree, airMopAlways);
172      nrrdNuke(nin);      fprintf(stderr, "%s: Trouble setting up gage context:\n%s", me, err);
173      return -1;      airMopError(mop);
}

status = gagePerVolumeAttach(ctx, pvl);
if(status) {
err = biffGetDone(GAGE);
fprintf(stderr, "Trouble attaching gage PerVolume:\n%s", err);
free(err);
gageContextNix(ctx);
nrrdNuke(nin);
return -1;
}

// It's okay to use kparm unitialized here: gageKernelSet ignores it.
status = gageKernelSet(ctx, gageKernel00, nrrdKernelBSpline3, kparm);
if(status) {
err = biffGetDone(GAGE);
fprintf(stderr, "Trouble setting kernel:\n%s", err);
free(err);
gageContextNix(ctx);
nrrdNuke(nin);
return -1;
}

// It's okay to use kparm unitialized here: gageKernelSet ignores it.
status = gageKernelSet(ctx, gageKernel11, nrrdKernelBSpline3D, kparm);
if(status) {
err = biffGetDone(GAGE);
fprintf(stderr, "Trouble setting kernel:\n%s", err);
free(err);
gageContextNix(ctx);
nrrdNuke(nin);
return -1;
}

status = gageQueryItemOn(ctx, pvl, gageSclValue);
if(status) {
err = biffGetDone(GAGE);
fprintf(stderr, "Trouble with gageQueryItemOn:\n%s", err);
free(err);
gageContextNix(ctx);
nrrdNuke(nin);
174      return -1;      return -1;
175    }    }
176
177    status = gageQueryItemOn(ctx, pvl, gageSclGradVec);    if (!( (val = gageAnswerPointer(ctx, pvl, gageSclValue))
178    if(status) {           && (grad = gageAnswerPointer(ctx, pvl, gageSclGradVec)) )) {
179      err = biffGetDone(GAGE);      fprintf(stderr, "%s: Trouble getting answer pointer\n", me);
180      fprintf(stderr, "Trouble with gageQueryItemOn:\n%s", err);      airMopError(mop);
free(err);
gageContextNix(ctx);
nrrdNuke(nin);
return -1;
}

status = gageUpdate(ctx);
if(status) {
err = biffGetDone(GAGE);
fprintf(stderr, "Trouble with gageUpdate:\n%s", err);
free(err);
gageContextNix(ctx);
nrrdNuke(nin);
return -1;
}

val = gageAnswerPointer(ctx, pvl, gageSclValue);
if(val == NULL) {
err = biffGetDone(GAGE);
fprintf(stderr, "Trouble getting answer pointer:\n%s", err);
free(err);
gageContextNix(ctx);
nrrdNuke(nin);
return -1;
}

grad = gageAnswerPointer(ctx, pvl, gageSclGradVec);
if(grad == NULL) {
err = biffGetDone(GAGE);
fprintf(stderr, "Trouble getting answer pointer:\n%s", err);
free(err);
gageContextNix(ctx);
nrrdNuke(nin);
181      return -1;      return -1;
182    }    }
183
184    int ui, vi;    int ui, vi;
185    double rayU, rayV, rayN;    double rayU, rayV, rayN;
186    double rayVec[3];    double rayVec[3];
187    double rayDir[3];    double toEye[3];
188    double rayPos[3];    double rayPos[3];
189    double transp, gray;    double transp, gray;
double output[4];
190    double norm[3];    double norm[3];
191    double alpha;    double alpha;
192    double ld;    double ld;
193    double hd;    double hd;
194    double mat;    double mat;
195    bool set = false;
196    double max_gray = 0;    // allocate output image
197    double min_gray = 0;    Nrrd *nimg = nrrdNew();
198    double *out_data = malloc(sizeof(double) * 4 * imgResU * imgResV);    airMopAdd(mop, nimg, (airMopper)nrrdNuke, airMopAlways);
199    if(out_data == NULL) {    if (nrrdAlloc_va(nimg, nrrdTypeDouble, 3,
200      fprintf(stderr, "Trouble with malloc.\n");                     AIR_CAST(size_t, 4),
201      gageContextNix(ctx);                     AIR_CAST(size_t, imgResU),
202      nrrdNuke(nin);                     AIR_CAST(size_t, imgResV))) {
203        airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
204        fprintf(stderr, "%s: Trouble allocating image output:\n%s", me, err);
205        airMopError(mop);
206      return -1;      return -1;
207    }    }
208      double *out_data = AIR_CAST(double *, nimg->data);
209    double t0 = GetTime(); // start timing    /*
210      fprintf(stderr, "hit return\n");
211      fgetc(stdin);
212      fprintf(stderr, "here we go\n");
213      */
214      double t0 = airTime(); // start timing
215    for(vi = 0; vi < imgResV; vi++) {    for(vi = 0; vi < imgResV; vi++) {
216      rayV = lerp(-camVmax, camVmax, -0.5, (double)(vi), (double)(imgResV) - 0.5);      rayV = lerp(-camVmax, camVmax, -0.5, (double)(vi), (double)(imgResV) - 0.5);
scale_vec3(rayV, camV, temp);
217      for(ui = 0; ui < imgResU; ui++) {      for(ui = 0; ui < imgResU; ui++) {
218          double len;
219        rayU = lerp(-camUmax, camUmax, -0.5, (double)(ui), (double)(imgResU) - 0.5);        rayU = lerp(-camUmax, camUmax, -0.5, (double)(ui), (double)(imgResU) - 0.5);
220        scale_vec3(rayU, camU, temp2);        ELL_3V_SCALE_ADD3(rayVec,   1.0, camN,   rayU/camDist, camU,   rayV/camDist, camV);
221        add3(temp, temp2, temp2);        ELL_3V_SCALE(toEye, -1, rayVec);
222        scale_vec3(camDist, camN, temp3);        ELL_3V_NORM(toEye, toEye, len);
add3(temp2, temp3, temp3);
scale_vec3(1 / camDist, temp3, rayVec);
norm3(rayVec, rayDir);
223        transp = 1;        transp = 1;
224        gray = 0;        gray = 0;
225        output[0] = 0;
output[1] = 0;
output[2] = 0;
output[3] = 0;
226        for(rayN = camVspNear; rayN <= camVspFar; rayN += rayStep) {        for(rayN = camVspNear; rayN <= camVspFar; rayN += rayStep) {
227          scale_vec3(rayN, rayDir, temp2);          ELL_3V_SCALE_ADD2(rayPos,   1.0, camEye,  rayN, rayVec);
228          add3(temp2, camEye, rayPos);
229          if(inside(rayPos)) {          if(inside(rayPos)) {
230            status = gageProbe(ctx, rayPos[0], rayPos[1], rayPos[2]);            if (gageProbeSpace(ctx, rayPos[0], rayPos[1], rayPos[2],
231            if(status) {                               AIR_FALSE /* index space */, AIR_TRUE /* clamp */)) {
232              err = ctx->errStr;              fprintf(stderr, "%s: Trouble with gageProbe:\n%s", me, ctx->errStr);
233              fprintf(stderr, "Trouble with gageProbe:\n%s", err);              airMopError(mop); return -1;
free(err);
free(out_data);
gageContextNix(ctx);
nrrdNuke(nin);
return -1;
234            }            }
235            if(*val > valOpacMin) {            if(*val > valOpacMin) {
236              temp2[0] = *(grad + 0);              ELL_3V_SCALE(temp2, -1, grad);
237              temp2[1] = *(grad + 1);              ELL_3V_NORM(norm, temp2, len);
temp2[2] = *(grad + 2);
sub3(Zero, temp2, temp2);
norm3(temp2, norm);
238              alpha = lerp(0.0, 1.0,  valOpacMin, *val, valOpacMax);              alpha = lerp(0.0, 1.0,  valOpacMin, *val, valOpacMax);
239              alpha = (alpha < 1.0) ? alpha : 1;              alpha = (alpha < 1.0) ? alpha : 1;
240              ld = dot3(norm, lightDir);              ld = dot3(norm, lightDir);
241              ld = (ld < 0) ? 0 : ld;              ld = (ld < 0) ? 0 : ld;
242              sub3(camEye, rayPos, temp2);              ELL_3V_ADD2(temp2, toEye, lightDir);
243              norm3(temp2, temp2);              ELL_3V_NORM(temp2, temp2, len);
add3(temp2, lightDir, temp2);
norm3(temp2, temp2);
244              hd = dot3(norm, temp2);              hd = dot3(norm, temp2);
245              hd = (hd < 0) ? 0 : hd;              hd = (hd < 0) ? 0 : hd;
246              mat = phongKa +              mat = (phongKa +
247                    phongKd * ((ld > 0) ? ld : 0) +                     phongKd * ld +
248                    phongKs * ((hd > 0) ? pow(hd, phongSp) : 0);                     phongKs * pow(hd, phongSp));
249              gray = gray + transp * alpha * mat;              gray = gray + transp * alpha * mat;
250              transp *= 1 - alpha;              transp *= 1 - alpha;
251            }            }
# Line 351  Line 255
255            break;            break;
256          }          }
257        }        }
*(out_data + ui * imgResV * 4 + vi * 4 + 0) = gray;
*(out_data + ui * imgResV * 4 + vi * 4 + 1) = gray;
*(out_data + ui * imgResV * 4 + vi * 4 + 2) = gray;
*(out_data + ui * imgResV * 4 + vi * 4 + 3) = 1.0f - transp;
if(!set) {
max_gray = gray;
min_gray = gray;
set = true;
}
max_gray = (gray > max_gray) ? gray : max_gray;
min_gray = (gray < min_gray) ? gray : min_gray;
}
}
258
259    unsigned char *uc_out_data = malloc(sizeof(unsigned char) * imgResU * imgResV);        ELL_4V_SET(out_data + 4*(ui + imgResU*vi), gray, gray, gray, 1.0-transp);
260    for(vi = 0; vi < imgResV; vi++) {
for(ui = 0; ui < imgResU; ui++) {
*(uc_out_data + vi * imgResU + ui) = ulerp(min_gray, *(out_data + ui * imgResV * 4 + vi * 4 + 0), max_gray);
261      }      }
262    }    }
263
264    double totalTime = GetTime() - t0;    double totalTime = airTime() - t0;
265    printf("usr=%f\n", totalTime);    printf("usr=%f\n", totalTime);
266
267    // Write out data    if (nrrdSave(outfile, nimg, NULL)) {
268    Nrrd *nout;      airMopAdd(mop, err = biffGetDone(NRRD), airFree, airMopAlways);
269    nout = nrrdNew();      fprintf(stderr, "%s: Trouble writing output:\n%s", me, err);
270    if(nout == NULL) {      airMopError(mop);
err = biffGetDone(NRRD);
fprintf(stderr, "Trouble allocating nrrd struct:\n%s", err);
free(err);
free(out_data);
free(uc_out_data);
gageContextNix(ctx);
nrrdNuke(nin);
271      return -1;      return -1;
272    }    }
273
274    status = nrrdWrap_va(nout, uc_out_data, nrrdTypeUChar, 2, imgResU, imgResV);    airMopOkay(mop);
if(status) {
err = biffGetDone(NRRD);
fprintf(stderr, "Trouble wrapping nrrd struct:\n%s", err);
free(err);
nrrdNix(nout);
free(out_data);
gageContextNix(ctx);
nrrdNuke(nin);
return -1;
}

status = nrrdSave(outfile, nout, NULL);
if(status) {
err = biffGetDone(NRRD);
fprintf(stderr, "Trouble saving nrrd struct:\n%s", err);
free(err);
nrrdNuke(nout);
free(out_data);
gageContextNix(ctx);
nrrdNuke(nin);
return -1;
}

free(out_data);
nrrdNuke(nout);
gageContextNix(ctx);
nrrdNuke(nin);
275    return 0;    return 0;
276
277  }  }

Legend:
 Removed from v.1575 changed lines Added in v.1606

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