ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/warp3d.c
(Generate patch)

Comparing ray/src/px/warp3d.c (file contents):
Revision 3.1 by greg, Tue Feb 4 16:03:48 1997 UTC vs.
Revision 3.8 by greg, Tue May 25 22:04:14 2004 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1997 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   * 3D warping routines.
6   */
7  
8   #include <stdio.h>
9 <
9 > #include <stdlib.h>
10   #include <math.h>
11  
12 + #include "rterror.h"
13 + #include "rtio.h"
14   #include "fvect.h"
16 –
15   #include "warp3d.h"
16  
17   #define MIND    1e-5            /* minimum distance between input points */
18  
19   typedef struct {
20          GNDX    n;              /* index must be first */
21 <        W3VEC   p;
21 >        W3VEC   v;
22   } KEYDP;                /* key/data allocation pair */
23  
24 < #define fgetvec(f,v)    (fgetval(f,'f',v) > 0 && fgetval(f,'f',v+1) > 0 \
25 <                                && fgetval(f,'f',v+2) > 0)
24 > #define fgetvec(p,v)    (fgetval(p,'f',v) > 0 && fgetval(p,'f',v+1) > 0 \
25 >                                && fgetval(p,'f',v+2) > 0)
26  
27 < #define AHUNK   32              /* number of points to allocate at a time */
27 > #define AHUNK   24              /* number of points to allocate at a time */
28  
31 – #ifndef malloc
32 – extern char     *malloc(), *realloc();
33 – #endif
34 – extern void     free();
29  
30 + double wpdist2(W3VEC p1, W3VEC p2);
31 + static int gridpoint(GNDX ndx, W3VEC rem, W3VEC ipt, struct grid3d *gp);
32 + static int get3dgpt(W3VEC ov, GNDX ndx, WARP3D *wp);
33 + static int get3dgin(W3VEC ov, GNDX ndx, W3VEC rem, WARP3D *wp);
34 + static void l3interp(W3VEC vo, W3VEC *cl, W3VEC dv, int n);
35 + static int warp3dex(W3VEC ov, W3VEC pi, WARP3D *wp);
36 + //static unsigned long gridhash(const void *gp);
37 + static lut_hashf_t gridhash;
38 + static int new3dgrid(WARP3D *wp);
39 + static void done3dgrid(struct grid3d *gp);
40  
41 +
42   double
43 < wpdist2(p1, p2)                 /* compute square of distance between points */
44 < register W3VEC  p1, p2;
43 > wpdist2(                        /* compute square of distance between points */
44 >        register W3VEC  p1,
45 >        register W3VEC  p2
46 > )
47   {
48          double  d, d2;
49  
# Line 48 | Line 55 | register W3VEC p1, p2;
55  
56  
57   int
58 < warp3d(po, pi, wp)              /* warp 3D point according to the given map */
59 < W3VEC   po, pi;
60 < register WARP3D *wp;
58 > warp3d(         /* warp 3D point according to the given map */
59 >        W3VEC   po,
60 >        W3VEC   pi,
61 >        register WARP3D *wp
62 > )
63   {
64          int     rval = W3OK;
65          GNDX    gi;
66 <        W3VEC   gd;
66 >        W3VEC   gd, ov;
67  
68 <        if (wp->grid.gn[0] == 0 && (rval = new3dwgrid(wp)) != W3OK)
68 >        if (wp->grid.gn[0] == 0 && (rval = new3dgrid(wp)) != W3OK)
69                  return(rval);
70          rval |= gridpoint(gi, gd, pi, &wp->grid);
71 <        if (wp->grid.flags & W3EXACT) {
72 <                rval |= warp3dex(po, pi, wp);
73 <                return(rval);
74 <        }
75 <        if (wp->grid.flags & W3FAST) {
76 <                rval |= get3dgpt(po, gi, wp);
77 <                return(rval);
78 <        }
79 <        rval |= get3dgin(po, gi, gd, wp);
71 >        if (wp->grid.flags & W3EXACT)
72 >                rval |= warp3dex(ov, pi, wp);
73 >        else if (wp->grid.flags & W3FAST)
74 >                rval |= get3dgpt(ov, gi, wp);
75 >        else
76 >                rval |= get3dgin(ov, gi, gd, wp);
77 >        po[0] = pi[0] + ov[0];
78 >        po[1] = pi[1] + ov[1];
79 >        po[2] = pi[2] + ov[2];
80          return(rval);
81   }
82  
83  
84 < int
85 < gridpoint(ndx, rem, ipt, gp)    /* compute grid position for ipt */
86 < GNDX    ndx;
87 < W3VEC   rem, ipt;
88 < register struct grid3d  *gp;
84 > static int
85 > gridpoint(      /* compute grid position for ipt */
86 >        GNDX    ndx,
87 >        W3VEC   rem,
88 >        W3VEC   ipt,
89 >        register struct grid3d  *gp
90 > )
91   {
92          int     rval = W3OK;
93          register int    i;
# Line 97 | Line 108 | register struct grid3d *gp;
108   }
109  
110  
111 < int
112 < get3dgpt(po, ndx, wp)           /* get value for voxel */
113 < W3VEC   po;
114 < GNDX    ndx;
115 < register WARP3D *wp;
111 > static int
112 > get3dgpt(               /* get value for voxel */
113 >        W3VEC   ov,
114 >        GNDX    ndx,
115 >        register WARP3D *wp
116 > )
117   {
118          W3VEC   gpt;
119          register LUENT  *le;
# Line 109 | Line 121 | register WARP3D        *wp;
121          int     rval = W3OK;
122          register int    i;
123  
124 <        le = lu_find(&wp->grid.gtab, (char *)ndx);
124 >        le = lu_find(&wp->grid.gtab, ndx);
125          if (le == NULL)
126                  return(W3ERROR);
127          if (le->data == NULL) {
# Line 123 | Line 135 | register WARP3D        *wp;
135                          if (wp->grid.flags & W3FAST)    /* on centers */
136                                  gpt[i] += .5*wp->grid.gstep[i];
137                  }
138 <                rval |= warp3dex(kd->p, gpt, wp);
138 >                rval = warp3dex(kd->v, gpt, wp);
139                  le->key = (char *)kd->n;
140 <                le->data = (char *)kd->p;
140 >                le->data = (char *)kd->v;
141          }
142 <        W3VCPY(po, (float *)le->data);
142 >        W3VCPY(ov, (float *)le->data);
143          return(rval);
144   }
145  
146  
147 < int
148 < get3dgin(po, ndx, rem, wp)      /* interpolate from warp grid */
149 < W3VEC   po, rem;
150 < GNDX    ndx;
151 < WARP3D  *wp;
147 > static int
148 > get3dgin(       /* interpolate from warp grid */
149 >        W3VEC   ov,
150 >        GNDX    ndx,
151 >        W3VEC   rem,
152 >        WARP3D  *wp
153 > )
154   {
155          W3VEC   cv[8];
142 –        int     rval;
156          GNDX    gi;
157 +        int     rval = W3OK;
158          register int    i;
159                                          /* get corner values */
160          for (i = 0; i < 8; i++) {
# Line 151 | Line 165 | WARP3D *wp;
165          }
166          if (rval & W3ERROR)
167                  return(rval);
168 <        l3interp(po, cv, rem, 3);       /* perform trilinear interpolation */
168 >        l3interp(ov, cv, rem, 3);       /* perform trilinear interpolation */
169          return(rval);
170   }
171  
172  
173 < l3interp(vo, cl, dv, n)         /* trilinear interpolation (recursive) */
174 < W3VEC   vo, *cl, dv;
175 < int     n;
173 > static void
174 > l3interp(               /* trilinear interpolation (recursive) */
175 >        W3VEC   vo,
176 >        W3VEC   *cl,
177 >        W3VEC   dv,
178 >        int     n
179 > )
180   {
181          W3VEC   v0, v1;
182          register int    i;
# Line 175 | Line 193 | int    n;
193   }
194  
195  
196 < int
197 < warp3dex(po, pi, wp)            /* compute warp using 1/r^2 weighted avg. */
198 < W3VEC   po, pi;
199 < register WARP3D *wp;
196 > static int
197 > warp3dex(               /* compute warp using 1/r^2 weighted avg. */
198 >        W3VEC   ov,
199 >        W3VEC   pi,
200 >        register WARP3D *wp
201 > )
202   {
203          double  d2, w, wsum;
204 <        W3VEC   pt;
204 >        W3VEC   vt;
205          register int    i;
206  
207 <        pt[0] = pt[1] = pt[2] = 0.;
207 >        vt[0] = vt[1] = vt[2] = 0.;
208          wsum = 0.;
209          for (i = wp->npts; i--; ) {
210                  d2 = wpdist2(pi, wp->ip[i]);
211                  if (d2 <= MIND*MIND) w = 1./(MIND*MIND);
212                  else w = 1./d2;
213 <                pt[0] += w*wp->op[i][0];
214 <                pt[1] += w*wp->op[i][1];
215 <                pt[2] += w*wp->op[i][2];
213 >                vt[0] += w*wp->ov[i][0];
214 >                vt[1] += w*wp->ov[i][1];
215 >                vt[2] += w*wp->ov[i][2];
216                  wsum += w;
217          }
218          if (wsum > 0.) {
219 <                po[0] = pt[0]/wsum;
220 <                po[1] = pt[1]/wsum;
221 <                po[2] = pt[2]/wsum;
219 >                ov[0] = vt[0]/wsum;
220 >                ov[1] = vt[1]/wsum;
221 >                ov[2] = vt[2]/wsum;
222          }
223          return(W3OK);
224   }
225  
226  
227   WARP3D *
228 < new3dw(flgs)                    /* allocate and initialize WARP3D struct */
229 < int     flgs;
228 > new3dw(                 /* allocate and initialize WARP3D struct */
229 >        int     flgs
230 > )
231   {
232          register WARP3D  *wp;
233  
# Line 218 | Line 239 | int    flgs;
239                  eputs("new3dw: no memory\n");
240                  return(NULL);
241          }
242 <        wp->ip = wp->op = NULL;
242 >        wp->ip = wp->ov = NULL;
243          wp->npts = 0;
244          wp->grid.flags = flgs;
245          wp->grid.gn[0] = wp->grid.gn[1] = wp->grid.gn[2] = 0;
# Line 227 | Line 248 | int    flgs;
248  
249  
250   WARP3D *
251 < load3dw(fn, wp)                 /* load 3D warp from file */
252 < char    *fn;
253 < WARP3D  *wp;
251 > load3dw(                        /* load 3D warp from file */
252 >        char    *fn,
253 >        WARP3D  *wp
254 > )
255   {
256          FILE    *fp;
257          W3VEC   inp, outp;
# Line 239 | Line 261 | WARP3D *wp;
261                  eputs(": cannot open\n");
262                  return(NULL);
263          }
264 <        if (wp == NULL && (wp = new3dw()) == NULL)
264 >        if (wp == NULL && (wp = new3dw(0)) == NULL)
265                  goto memerr;
266          while (fgetvec(fp, inp) && fgetvec(fp, outp))
267                  if (!add3dpt(wp, inp, outp))
# Line 257 | Line 279 | cleanup:
279   }
280  
281  
282 + extern int
283 + set3dwfl(               /* reset warp flags */
284 +        register WARP3D *wp,
285 +        int     flgs
286 + )
287 + {
288 +        if (flgs == wp->grid.flags)
289 +                return(0);
290 +        if ((flgs & (W3EXACT|W3FAST)) == (W3EXACT|W3FAST)) {
291 +                eputs("set3dwfl: only one of W3EXACT or W3FAST\n");
292 +                return(-1);
293 +        }
294 +        wp->grid.flags = flgs;
295 +        done3dgrid(&wp->grid);          /* old grid is invalid */
296 +        return(0);
297 + }
298 +
299 +
300   int
301 < add3dpt(wp, pti, pto)           /* add 3D point pair to warp structure */
302 < register WARP3D *wp;
303 < W3VEC   pti, pto;
301 > add3dpt(                /* add 3D point pair to warp structure */
302 >        register WARP3D *wp,
303 >        W3VEC   pti,
304 >        W3VEC   pto
305 > )
306   {
307          double  d2;
308          register W3VEC  *na;
# Line 269 | Line 311 | W3VEC  pti, pto;
311          if (wp->npts == 0) {                    /* initialize */
312                  wp->ip = (W3VEC *)malloc(AHUNK*sizeof(W3VEC));
313                  if (wp->ip == NULL) return(0);
314 <                wp->op = (W3VEC *)malloc(AHUNK*sizeof(W3VEC));
315 <                if (wp->op == NULL) return(0);
314 >                wp->ov = (W3VEC *)malloc(AHUNK*sizeof(W3VEC));
315 >                if (wp->ov == NULL) return(0);
316                  wp->d2min = 1e10;
317                  wp->d2max = 0.;
318                  W3VCPY(wp->llim, pti);
319                  W3VCPY(wp->ulim, pti);
320          } else {
321                  if (wp->npts % AHUNK == 0) {    /* allocate another hunk */
322 <                        na = (W3VEC *)realloc((char *)wp->ip,
322 >                        na = (W3VEC *)realloc((void *)wp->ip,
323                                          (wp->npts+AHUNK)*sizeof(W3VEC));
324                          if (na == NULL) return(0);
325                          wp->ip = na;
326 <                        na = (W3VEC *)realloc((char *)wp->op,
326 >                        na = (W3VEC *)realloc((void *)wp->ov,
327                                          (wp->npts+AHUNK)*sizeof(W3VEC));
328                          if (na == NULL) return(0);
329 <                        wp->op = na;
329 >                        wp->ov = na;
330                  }
331                  for (i = 0; i < 3; i++)         /* check boundaries */
332                          if (pti[i] < wp->llim[i])
# Line 302 | Line 344 | W3VEC  pti, pto;
344                  }
345          }
346          W3VCPY(wp->ip[wp->npts], pti);  /* add new point */
347 <        W3VCPY(wp->op[wp->npts], pto);
347 >        wp->ov[wp->npts][0] = pto[0] - pti[0];
348 >        wp->ov[wp->npts][1] = pto[1] - pti[1];
349 >        wp->ov[wp->npts][2] = pto[2] - pti[2];
350          done3dgrid(&wp->grid);          /* old grid is invalid */
351          return(++wp->npts);
352   }
353  
354  
355 < free3dw(wp)                     /* free WARP3D data */
356 < register WARP3D *wp;
355 > extern void
356 > free3dw(                        /* free WARP3D data */
357 >        register WARP3D *wp
358 > )
359   {
360          done3dgrid(&wp->grid);
361 <        free((char *)wp->ip);
362 <        free((char *)wp->op);
363 <        free((char *)wp);
361 >        free((void *)wp->ip);
362 >        free((void *)wp->ov);
363 >        free((void *)wp);
364   }
365  
366  
367 < long
368 < gridhash(gp)                    /* hash a grid point index */
369 < GNDX    gp;
367 > static unsigned long
368 > gridhash(                       /* hash a grid point index */
369 >        //GNDX  gp
370 >        const void      *gp
371 > )
372   {
373 <        return(((unsigned long)gp[0]<<GNBITS | gp[1])<<GNBITS | gp[2]);
373 >        //return(((unsigned long)gp[0]<<GNBITS | gp[1])<<GNBITS | gp[2]);
374 >        return(((unsigned long)((const unsigned char*)gp)[0]<<GNBITS | ((const unsigned char*)gp)[1])<<GNBITS | ((const unsigned char*)gp)[2]);
375   }
376  
377  
378 < int
379 < new3dwgrid(wp)                  /* initialize interpolating grid for warp */
380 < register WARP3D *wp;
378 > static int
379 > new3dgrid(                      /* initialize interpolating grid for warp */
380 >        register WARP3D *wp
381 > )
382   {
383 +        W3VEC   gmax;
384          double  gridstep, d;
385 +        int     n;
386          register int    i;
387                                  /* free old grid (if any) */
388          done3dgrid(&wp->grid);
# Line 341 | Line 393 | register WARP3D        *wp;
393                  return(W3BADMAP);       /* not enough disparity */
394                                  /* compute gamut */
395          W3VCPY(wp->grid.gmin, wp->llim);
396 <        W3VCPY(wp->grid.gmax, wp->ulim);
396 >        W3VCPY(gmax, wp->ulim);
397          gridstep = sqrt(wp->d2min);
398          for (i = 0; i < 3; i++) {
399                  wp->grid.gmin[i] -= .501*gridstep;
400 <                wp->grid.gmax[i] += .501*gridstep;
400 >                gmax[i] += .501*gridstep;
401          }
402          if (wp->grid.flags & W3EXACT) {
403                  wp->grid.gn[0] = wp->grid.gn[1] = wp->grid.gn[2] = 1;
404 +                wp->grid.gstep[0] = gmax[0] - wp->grid.gmin[0];
405 +                wp->grid.gstep[1] = gmax[1] - wp->grid.gmin[1];
406 +                wp->grid.gstep[2] = gmax[2] - wp->grid.gmin[2];
407                  return(W3OK);           /* no interpolation, so no grid */
408          }
409                                  /* create grid */
410          for (i = 0; i < 3; i++) {
411 <                d = wp->grid.gmax[i] - wp->grid.gmin[i];
412 <                wp->grid.gn[i] = d/gridstep + .5;
413 <                if (wp->grid.gn[i] >= MAXGN)
414 <                        wp->grid.gn[i] = MAXGN-1;
415 <                wp->grid.gstep[i] = d / wp->grid.gn[i];
411 >                d = gmax[i] - wp->grid.gmin[i];
412 >                n = d/gridstep + .5;
413 >                if (n >= MAXGN-1)
414 >                        n = MAXGN-2;
415 >                wp->grid.gstep[i] = d / n;
416 >                wp->grid.gn[i] = n;
417          }
418                                  /* initialize lookup table */
419          wp->grid.gtab.hashf = gridhash;
# Line 368 | Line 424 | register WARP3D        *wp;
424   }
425  
426  
427 < done3dgrid(gp)                  /* free interpolating grid for warp */
428 < register struct grid3d  *gp;
427 > static void
428 > done3dgrid(                     /* free interpolating grid for warp */
429 >        register struct grid3d  *gp
430 > )
431   {
432          if (gp->gn[0] == 0)
433                  return;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines