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

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: warp3d.c,v 3.7 2004/03/28 20:33:14 schorsch Exp $";
3 #endif
4 /*
5 * 3D warping routines.
6 */
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <math.h>
11
12 #include "rterror.h"
13 #include "rtio.h"
14 #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 W3VEC v;
22 } KEYDP; /* key/data allocation pair */
23
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 24 /* number of points to allocate at a time */
28
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( /* compute square of distance between points */
44 register W3VEC p1,
45 register W3VEC p2
46 )
47 {
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 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, ov;
67
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(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 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;
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 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;
120 KEYDP *kd;
121 int rval = W3OK;
122 register int i;
123
124 le = lu_find(&wp->grid.gtab, ndx);
125 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 rval = warp3dex(kd->v, gpt, wp);
139 le->key = (char *)kd->n;
140 le->data = (char *)kd->v;
141 }
142 W3VCPY(ov, (float *)le->data);
143 return(rval);
144 }
145
146
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];
156 GNDX gi;
157 int rval = W3OK;
158 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 l3interp(ov, cv, rem, 3); /* perform trilinear interpolation */
169 return(rval);
170 }
171
172
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;
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 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 vt;
205 register int i;
206
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 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 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( /* allocate and initialize WARP3D struct */
229 int flgs
230 )
231 {
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 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;
246 return(wp);
247 }
248
249
250 WARP3D *
251 load3dw( /* load 3D warp from file */
252 char *fn,
253 WARP3D *wp
254 )
255 {
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 if (wp == NULL && (wp = new3dw(0)) == NULL)
265 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 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( /* 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;
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 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((void *)wp->ip,
323 (wp->npts+AHUNK)*sizeof(W3VEC));
324 if (na == NULL) return(0);
325 wp->ip = na;
326 na = (W3VEC *)realloc((void *)wp->ov,
327 (wp->npts+AHUNK)*sizeof(W3VEC));
328 if (na == NULL) return(0);
329 wp->ov = na;
330 }
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 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 extern void
356 free3dw( /* free WARP3D data */
357 register WARP3D *wp
358 )
359 {
360 done3dgrid(&wp->grid);
361 free((void *)wp->ip);
362 free((void *)wp->ov);
363 free((void *)wp);
364 }
365
366
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]);
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 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);
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 W3VCPY(gmax, wp->ulim);
397 gridstep = sqrt(wp->d2min);
398 for (i = 0; i < 3; i++) {
399 wp->grid.gmin[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 = 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;
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 static void
428 done3dgrid( /* free interpolating grid for warp */
429 register struct grid3d *gp
430 )
431 {
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 }