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

Diff Legend

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