ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/px/warp3d.c
Revision: 3.10
Committed: Fri Feb 18 00:40:25 2011 UTC (13 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2, rad4R2P2, rad5R0, rad5R1, rad4R2, rad4R1, rad4R2P1
Changes since 3.9: +2 -2 lines
Log Message:
Major code reorg, moving mgflib to common and introducing BSDF material

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: warp3d.c,v 3.9 2010/09/03 21:18:15 greg 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 lut_hashf_t gridhash;
37 static int new3dgrid(WARP3D *wp);
38 static void done3dgrid(struct grid3d *gp);
39
40
41 double
42 wpdist2( /* compute square of distance between points */
43 register W3VEC p1,
44 register W3VEC p2
45 )
46 {
47 double d, d2;
48
49 d = p1[0] - p2[0]; d2 = d*d;
50 d = p1[1] - p2[1]; d2 += d*d;
51 d = p1[2] - p2[2]; d2 += d*d;
52 return(d2);
53 }
54
55
56 int
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, ov;
66
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(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 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;
93
94 for (i = 0; i < 3; i++) {
95 rem[i] = (ipt[i] - gp->gmin[i])/gp->gstep[i];
96 if (rem[i] < 0.) {
97 ndx[i] = 0;
98 rval = W3GAMUT;
99 } else if ((int)rem[i] >= gp->gn[i]) {
100 ndx[i] = gp->gn[i] - 1;
101 rval = W3GAMUT;
102 } else
103 ndx[i] = (int)rem[i];
104 rem[i] -= (double)ndx[i];
105 }
106 return(rval);
107 }
108
109
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;
119 KEYDP *kd;
120 int rval = W3OK;
121 register int i;
122
123 le = lu_find(&wp->grid.gtab, ndx);
124 if (le == NULL)
125 return(W3ERROR);
126 if (le->data == NULL) {
127 if (le->key != NULL)
128 kd = (KEYDP *)le->key;
129 else if ((kd = (KEYDP *)malloc(sizeof(KEYDP))) == NULL)
130 return(W3ERROR);
131 for (i = 0; i < 3; i++) {
132 kd->n[i] = ndx[i];
133 gpt[i] = wp->grid.gmin[i] + ndx[i]*wp->grid.gstep[i];
134 if (wp->grid.flags & W3FAST) /* on centers */
135 gpt[i] += .5*wp->grid.gstep[i];
136 }
137 rval = warp3dex(kd->v, gpt, wp);
138 le->key = (char *)kd->n;
139 le->data = (char *)kd->v;
140 }
141 W3VCPY(ov, (float *)le->data);
142 return(rval);
143 }
144
145
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];
155 GNDX gi;
156 int rval = W3OK;
157 register int i;
158 /* get corner values */
159 for (i = 0; i < 8; i++) {
160 gi[0] = ndx[0] + (i & 1);
161 gi[1] = ndx[1] + (i>>1 & 1);
162 gi[2] = ndx[2] + (i>>2);
163 rval |= get3dgpt(cv[i], gi, wp);
164 }
165 if (rval & W3ERROR)
166 return(rval);
167 l3interp(ov, cv, rem, 3); /* perform trilinear interpolation */
168 return(rval);
169 }
170
171
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;
182
183 if (--n) {
184 l3interp(v0, cl, dv, n);
185 l3interp(v1, cl+(1<<n), dv, n);
186 } else {
187 W3VCPY(v0, cl[0]);
188 W3VCPY(v1, cl[1]);
189 }
190 for (i = 0; i < 3; i++)
191 vo[i] = (1.-dv[n])*v0[i] + dv[n]*v1[i];
192 }
193
194
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 vt;
204 register int i;
205
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 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 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( /* allocate and initialize WARP3D struct */
228 int flgs
229 )
230 {
231 register WARP3D *wp;
232
233 if ((flgs & (W3EXACT|W3FAST)) == (W3EXACT|W3FAST)) {
234 eputs("new3dw: only one of W3EXACT or W3FAST\n");
235 return(NULL);
236 }
237 if ((wp = (WARP3D *)malloc(sizeof(WARP3D))) == NULL) {
238 eputs("new3dw: no memory\n");
239 return(NULL);
240 }
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;
245 return(wp);
246 }
247
248
249 WARP3D *
250 load3dw( /* load 3D warp from file */
251 char *fn,
252 WARP3D *wp
253 )
254 {
255 FILE *fp;
256 W3VEC inp, outp;
257
258 if ((fp = fopen(fn, "r")) == NULL) {
259 eputs(fn);
260 eputs(": cannot open\n");
261 return(NULL);
262 }
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))
267 goto memerr;
268 if (!feof(fp)) {
269 wputs(fn);
270 wputs(": non-number in 3D warp file\n");
271 }
272 goto cleanup;
273 memerr:
274 eputs("load3dw: out of memory\n");
275 cleanup:
276 fclose(fp);
277 return(wp);
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( /* 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;
308 register int i;
309
310 if (wp->npts == 0) { /* initialize */
311 wp->ip = (W3VEC *)malloc(AHUNK*sizeof(W3VEC));
312 if (wp->ip == 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((void *)wp->ip,
322 (wp->npts+AHUNK)*sizeof(W3VEC));
323 if (na == NULL) return(0);
324 wp->ip = na;
325 na = (W3VEC *)realloc((void *)wp->ov,
326 (wp->npts+AHUNK)*sizeof(W3VEC));
327 if (na == NULL) return(0);
328 wp->ov = na;
329 }
330 for (i = 0; i < 3; i++) /* check boundaries */
331 if (pti[i] < wp->llim[i])
332 wp->llim[i] = pti[i];
333 else if (pti[i] > wp->ulim[i])
334 wp->ulim[i] = pti[i];
335 for (i = wp->npts; i--; ) { /* check distances */
336 d2 = wpdist2(pti, wp->ip[i]);
337 if (d2 < MIND*MIND) /* value is too close */
338 return(wp->npts);
339 if (d2 < wp->d2min)
340 wp->d2min = d2;
341 if (d2 > wp->d2max)
342 wp->d2max = d2;
343 }
344 }
345 W3VCPY(wp->ip[wp->npts], pti); /* add new point */
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 extern void
355 free3dw( /* free WARP3D data */
356 register WARP3D *wp
357 )
358 {
359 done3dgrid(&wp->grid);
360 free((void *)wp->ip);
361 free((void *)wp->ov);
362 free((void *)wp);
363 }
364
365
366 static unsigned long
367 gridhash( /* hash a grid point index */
368 const char *gp
369 )
370 {
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 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);
386 /* check parameters */
387 if (wp->npts < 2)
388 return(W3BADMAP); /* undefined for less than 2 points */
389 if (wp->d2max < MIND)
390 return(W3BADMAP); /* not enough disparity */
391 /* compute gamut */
392 W3VCPY(wp->grid.gmin, wp->llim);
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 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 = 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;
417 wp->grid.gtab.keycmp = NULL;
418 wp->grid.gtab.freek = free;
419 wp->grid.gtab.freed = NULL;
420 return(lu_init(&wp->grid.gtab, 1024) ? W3OK : W3ERROR);
421 }
422
423
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;
431 lu_done(&gp->gtab);
432 gp->gn[0] = gp->gn[1] = gp->gn[2] = 0;
433 }