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

Diff Legend

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