ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/warp3d.c
Revision: 3.4
Committed: Wed Sep 2 18:41:44 1998 UTC (25 years, 7 months ago) by gwlarson
Content type: text/plain
Branch: MAIN
Changes since 3.3: +3 -3 lines
Log Message:
changed lookup routines to use unsigned long hash values

File Contents

# User Rev Content
1 gwlarson 3.4 /* Copyright (c) 1998 Silicon Graphics, Inc. */
2 greg 3.1
3     #ifndef lint
4 gwlarson 3.4 static char SCCSid[] = "$SunId$ SGI";
5 greg 3.1 #endif
6    
7     /*
8     * 3D warping routines.
9     */
10    
11     #include <stdio.h>
12     #include <math.h>
13     #include "fvect.h"
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 greg 3.2 W3VEC v;
21 greg 3.1 } KEYDP; /* key/data allocation pair */
22    
23 greg 3.2 #define fgetvec(p,v) (fgetval(p,'f',v) > 0 && fgetval(p,'f',v+1) > 0 \
24     && fgetval(p,'f',v+2) > 0)
25 greg 3.1
26 greg 3.2 #define AHUNK 24 /* number of points to allocate at a time */
27 greg 3.1
28     #ifndef malloc
29     extern char *malloc(), *realloc();
30     #endif
31     extern void free();
32    
33    
34     double
35     wpdist2(p1, p2) /* compute square of distance between points */
36     register W3VEC p1, p2;
37     {
38     double d, d2;
39    
40     d = p1[0] - p2[0]; d2 = d*d;
41     d = p1[1] - p2[1]; d2 += d*d;
42     d = p1[2] - p2[2]; d2 += d*d;
43     return(d2);
44     }
45    
46    
47     int
48     warp3d(po, pi, wp) /* warp 3D point according to the given map */
49     W3VEC po, pi;
50     register WARP3D *wp;
51     {
52     int rval = W3OK;
53     GNDX gi;
54 greg 3.2 W3VEC gd, ov;
55 greg 3.1
56 greg 3.2 if (wp->grid.gn[0] == 0 && (rval = new3dgrid(wp)) != W3OK)
57 greg 3.1 return(rval);
58     rval |= gridpoint(gi, gd, pi, &wp->grid);
59 greg 3.2 if (wp->grid.flags & W3EXACT)
60     rval |= warp3dex(ov, pi, wp);
61     else if (wp->grid.flags & W3FAST)
62     rval |= get3dgpt(ov, gi, wp);
63     else
64     rval |= get3dgin(ov, gi, gd, wp);
65     po[0] = pi[0] + ov[0];
66     po[1] = pi[1] + ov[1];
67     po[2] = pi[2] + ov[2];
68 greg 3.1 return(rval);
69     }
70    
71    
72     int
73     gridpoint(ndx, rem, ipt, gp) /* compute grid position for ipt */
74     GNDX ndx;
75     W3VEC rem, ipt;
76     register struct grid3d *gp;
77     {
78     int rval = W3OK;
79     register int i;
80    
81     for (i = 0; i < 3; i++) {
82     rem[i] = (ipt[i] - gp->gmin[i])/gp->gstep[i];
83     if (rem[i] < 0.) {
84     ndx[i] = 0;
85     rval = W3GAMUT;
86     } else if ((int)rem[i] >= gp->gn[i]) {
87     ndx[i] = gp->gn[i] - 1;
88     rval = W3GAMUT;
89     } else
90     ndx[i] = (int)rem[i];
91     rem[i] -= (double)ndx[i];
92     }
93     return(rval);
94     }
95    
96    
97     int
98 greg 3.2 get3dgpt(ov, ndx, wp) /* get value for voxel */
99     W3VEC ov;
100 greg 3.1 GNDX ndx;
101     register WARP3D *wp;
102     {
103     W3VEC gpt;
104     register LUENT *le;
105     KEYDP *kd;
106     int rval = W3OK;
107     register int i;
108    
109 greg 3.2 le = lu_find(&wp->grid.gtab, ndx);
110 greg 3.1 if (le == NULL)
111     return(W3ERROR);
112     if (le->data == NULL) {
113     if (le->key != NULL)
114     kd = (KEYDP *)le->key;
115     else if ((kd = (KEYDP *)malloc(sizeof(KEYDP))) == NULL)
116     return(W3ERROR);
117     for (i = 0; i < 3; i++) {
118     kd->n[i] = ndx[i];
119     gpt[i] = wp->grid.gmin[i] + ndx[i]*wp->grid.gstep[i];
120     if (wp->grid.flags & W3FAST) /* on centers */
121     gpt[i] += .5*wp->grid.gstep[i];
122     }
123 greg 3.2 rval = warp3dex(kd->v, gpt, wp);
124 greg 3.1 le->key = (char *)kd->n;
125 greg 3.2 le->data = (char *)kd->v;
126 greg 3.1 }
127 greg 3.2 W3VCPY(ov, (float *)le->data);
128 greg 3.1 return(rval);
129     }
130    
131    
132     int
133 greg 3.2 get3dgin(ov, ndx, rem, wp) /* interpolate from warp grid */
134     W3VEC ov, rem;
135 greg 3.1 GNDX ndx;
136     WARP3D *wp;
137     {
138     W3VEC cv[8];
139     GNDX gi;
140 greg 3.3 int rval = W3OK;
141 greg 3.1 register int i;
142     /* get corner values */
143     for (i = 0; i < 8; i++) {
144     gi[0] = ndx[0] + (i & 1);
145     gi[1] = ndx[1] + (i>>1 & 1);
146     gi[2] = ndx[2] + (i>>2);
147     rval |= get3dgpt(cv[i], gi, wp);
148     }
149     if (rval & W3ERROR)
150     return(rval);
151 greg 3.2 l3interp(ov, cv, rem, 3); /* perform trilinear interpolation */
152 greg 3.1 return(rval);
153     }
154    
155    
156     l3interp(vo, cl, dv, n) /* trilinear interpolation (recursive) */
157     W3VEC vo, *cl, dv;
158     int n;
159     {
160     W3VEC v0, v1;
161     register int i;
162    
163     if (--n) {
164     l3interp(v0, cl, dv, n);
165     l3interp(v1, cl+(1<<n), dv, n);
166     } else {
167     W3VCPY(v0, cl[0]);
168     W3VCPY(v1, cl[1]);
169     }
170     for (i = 0; i < 3; i++)
171     vo[i] = (1.-dv[n])*v0[i] + dv[n]*v1[i];
172     }
173    
174    
175     int
176 greg 3.2 warp3dex(ov, pi, wp) /* compute warp using 1/r^2 weighted avg. */
177     W3VEC ov, pi;
178 greg 3.1 register WARP3D *wp;
179     {
180     double d2, w, wsum;
181 greg 3.2 W3VEC vt;
182 greg 3.1 register int i;
183    
184 greg 3.2 vt[0] = vt[1] = vt[2] = 0.;
185 greg 3.1 wsum = 0.;
186     for (i = wp->npts; i--; ) {
187     d2 = wpdist2(pi, wp->ip[i]);
188     if (d2 <= MIND*MIND) w = 1./(MIND*MIND);
189     else w = 1./d2;
190 greg 3.2 vt[0] += w*wp->ov[i][0];
191     vt[1] += w*wp->ov[i][1];
192     vt[2] += w*wp->ov[i][2];
193 greg 3.1 wsum += w;
194     }
195     if (wsum > 0.) {
196 greg 3.2 ov[0] = vt[0]/wsum;
197     ov[1] = vt[1]/wsum;
198     ov[2] = vt[2]/wsum;
199 greg 3.1 }
200     return(W3OK);
201     }
202    
203    
204     WARP3D *
205     new3dw(flgs) /* allocate and initialize WARP3D struct */
206     int flgs;
207     {
208     register WARP3D *wp;
209    
210     if ((flgs & (W3EXACT|W3FAST)) == (W3EXACT|W3FAST)) {
211     eputs("new3dw: only one of W3EXACT or W3FAST\n");
212     return(NULL);
213     }
214     if ((wp = (WARP3D *)malloc(sizeof(WARP3D))) == NULL) {
215     eputs("new3dw: no memory\n");
216     return(NULL);
217     }
218 greg 3.2 wp->ip = wp->ov = NULL;
219 greg 3.1 wp->npts = 0;
220     wp->grid.flags = flgs;
221     wp->grid.gn[0] = wp->grid.gn[1] = wp->grid.gn[2] = 0;
222     return(wp);
223     }
224    
225    
226     WARP3D *
227     load3dw(fn, wp) /* load 3D warp from file */
228     char *fn;
229     WARP3D *wp;
230     {
231     FILE *fp;
232     W3VEC inp, outp;
233    
234     if ((fp = fopen(fn, "r")) == NULL) {
235     eputs(fn);
236     eputs(": cannot open\n");
237     return(NULL);
238     }
239 greg 3.2 if (wp == NULL && (wp = new3dw(0)) == NULL)
240 greg 3.1 goto memerr;
241     while (fgetvec(fp, inp) && fgetvec(fp, outp))
242     if (!add3dpt(wp, inp, outp))
243     goto memerr;
244     if (!feof(fp)) {
245     wputs(fn);
246     wputs(": non-number in 3D warp file\n");
247     }
248     goto cleanup;
249     memerr:
250     eputs("load3dw: out of memory\n");
251     cleanup:
252     fclose(fp);
253     return(wp);
254     }
255    
256    
257     int
258 greg 3.2 set3dwfl(wp, flgs) /* reset warp flags */
259     register WARP3D *wp;
260     int flgs;
261     {
262     if (flgs == wp->grid.flags)
263     return(0);
264     if ((flgs & (W3EXACT|W3FAST)) == (W3EXACT|W3FAST)) {
265     eputs("set3dwfl: only one of W3EXACT or W3FAST\n");
266     return(-1);
267     }
268     wp->grid.flags = flgs;
269     done3dgrid(&wp->grid); /* old grid is invalid */
270     return(0);
271     }
272    
273    
274     int
275 greg 3.1 add3dpt(wp, pti, pto) /* add 3D point pair to warp structure */
276     register WARP3D *wp;
277     W3VEC pti, pto;
278     {
279     double d2;
280     register W3VEC *na;
281     register int i;
282    
283     if (wp->npts == 0) { /* initialize */
284     wp->ip = (W3VEC *)malloc(AHUNK*sizeof(W3VEC));
285     if (wp->ip == NULL) return(0);
286 greg 3.2 wp->ov = (W3VEC *)malloc(AHUNK*sizeof(W3VEC));
287     if (wp->ov == NULL) return(0);
288 greg 3.1 wp->d2min = 1e10;
289     wp->d2max = 0.;
290     W3VCPY(wp->llim, pti);
291     W3VCPY(wp->ulim, pti);
292     } else {
293     if (wp->npts % AHUNK == 0) { /* allocate another hunk */
294     na = (W3VEC *)realloc((char *)wp->ip,
295     (wp->npts+AHUNK)*sizeof(W3VEC));
296     if (na == NULL) return(0);
297     wp->ip = na;
298 greg 3.2 na = (W3VEC *)realloc((char *)wp->ov,
299 greg 3.1 (wp->npts+AHUNK)*sizeof(W3VEC));
300     if (na == NULL) return(0);
301 greg 3.2 wp->ov = na;
302 greg 3.1 }
303     for (i = 0; i < 3; i++) /* check boundaries */
304     if (pti[i] < wp->llim[i])
305     wp->llim[i] = pti[i];
306     else if (pti[i] > wp->ulim[i])
307     wp->ulim[i] = pti[i];
308     for (i = wp->npts; i--; ) { /* check distances */
309     d2 = wpdist2(pti, wp->ip[i]);
310     if (d2 < MIND*MIND) /* value is too close */
311     return(wp->npts);
312     if (d2 < wp->d2min)
313     wp->d2min = d2;
314     if (d2 > wp->d2max)
315     wp->d2max = d2;
316     }
317     }
318     W3VCPY(wp->ip[wp->npts], pti); /* add new point */
319 greg 3.2 wp->ov[wp->npts][0] = pto[0] - pti[0];
320     wp->ov[wp->npts][1] = pto[1] - pti[1];
321     wp->ov[wp->npts][2] = pto[2] - pti[2];
322 greg 3.1 done3dgrid(&wp->grid); /* old grid is invalid */
323     return(++wp->npts);
324     }
325    
326    
327     free3dw(wp) /* free WARP3D data */
328     register WARP3D *wp;
329     {
330     done3dgrid(&wp->grid);
331     free((char *)wp->ip);
332 greg 3.2 free((char *)wp->ov);
333 greg 3.1 free((char *)wp);
334     }
335    
336    
337 gwlarson 3.4 unsigned long
338 greg 3.1 gridhash(gp) /* hash a grid point index */
339     GNDX gp;
340     {
341     return(((unsigned long)gp[0]<<GNBITS | gp[1])<<GNBITS | gp[2]);
342     }
343    
344    
345     int
346 greg 3.2 new3dgrid(wp) /* initialize interpolating grid for warp */
347 greg 3.1 register WARP3D *wp;
348     {
349 greg 3.2 W3VEC gmax;
350 greg 3.1 double gridstep, d;
351 greg 3.2 int n;
352 greg 3.1 register int i;
353     /* free old grid (if any) */
354     done3dgrid(&wp->grid);
355     /* check parameters */
356     if (wp->npts < 2)
357     return(W3BADMAP); /* undefined for less than 2 points */
358     if (wp->d2max < MIND)
359     return(W3BADMAP); /* not enough disparity */
360     /* compute gamut */
361     W3VCPY(wp->grid.gmin, wp->llim);
362 greg 3.2 W3VCPY(gmax, wp->ulim);
363 greg 3.1 gridstep = sqrt(wp->d2min);
364     for (i = 0; i < 3; i++) {
365     wp->grid.gmin[i] -= .501*gridstep;
366 greg 3.2 gmax[i] += .501*gridstep;
367 greg 3.1 }
368     if (wp->grid.flags & W3EXACT) {
369     wp->grid.gn[0] = wp->grid.gn[1] = wp->grid.gn[2] = 1;
370 greg 3.3 wp->grid.gstep[0] = gmax[0] - wp->grid.gmin[0];
371     wp->grid.gstep[1] = gmax[1] - wp->grid.gmin[1];
372     wp->grid.gstep[2] = gmax[2] - wp->grid.gmin[2];
373 greg 3.1 return(W3OK); /* no interpolation, so no grid */
374     }
375     /* create grid */
376     for (i = 0; i < 3; i++) {
377 greg 3.2 d = gmax[i] - wp->grid.gmin[i];
378     n = d/gridstep + .5;
379     if (n >= MAXGN-1)
380     n = MAXGN-2;
381     wp->grid.gstep[i] = d / n;
382     wp->grid.gn[i] = n;
383 greg 3.1 }
384     /* initialize lookup table */
385     wp->grid.gtab.hashf = gridhash;
386     wp->grid.gtab.keycmp = NULL;
387     wp->grid.gtab.freek = free;
388     wp->grid.gtab.freed = NULL;
389     return(lu_init(&wp->grid.gtab, 1024) ? W3OK : W3ERROR);
390     }
391    
392    
393     done3dgrid(gp) /* free interpolating grid for warp */
394     register struct grid3d *gp;
395     {
396     if (gp->gn[0] == 0)
397     return;
398     lu_done(&gp->gtab);
399     gp->gn[0] = gp->gn[1] = gp->gn[2] = 0;
400     }