ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/warp3d.c
Revision: 3.8
Committed: Tue May 25 22:04:14 2004 UTC (19 years, 10 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad3R7P2, rad3R7P1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9
Changes since 3.7: +4 -4 lines
Log Message:
Added const modifier to key and other parameters in lookup.h

File Contents

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