8 |
|
* for Irradiance Caching" by Schwarzhaupt, Wann Jensen, & Jarosz |
9 |
|
* from ACM SIGGRAPH Asia 2012 conference proceedings. |
10 |
|
* |
11 |
+ |
* Added book-keeping optimization to avoid calculations that would |
12 |
+ |
* cancel due to traversal both directions on edges that are adjacent |
13 |
+ |
* to same-valued triangles. This cuts about half of Hessian math. |
14 |
+ |
* |
15 |
|
* Declarations of external symbols in ambient.h |
16 |
|
*/ |
17 |
|
|
23 |
|
|
24 |
|
#ifdef NEWAMB |
25 |
|
|
26 |
+ |
/* #define HEM_MULT 4.0 /* hem multiplier (bigger => sparser cache) */ |
27 |
+ |
|
28 |
|
extern void SDsquare2disk(double ds[2], double seedx, double seedy); |
29 |
|
|
30 |
+ |
/* vertex direction bit positions */ |
31 |
+ |
#define VDB_xy 0 |
32 |
+ |
#define VDB_y 01 |
33 |
+ |
#define VDB_x 02 |
34 |
+ |
#define VDB_Xy 03 |
35 |
+ |
#define VDB_xY 04 |
36 |
+ |
#define VDB_X 05 |
37 |
+ |
#define VDB_Y 06 |
38 |
+ |
#define VDB_XY 07 |
39 |
+ |
/* get opposite vertex direction bit */ |
40 |
+ |
#define VDB_OPP(f) (~(f) & 07) |
41 |
+ |
/* adjacent triangle vertex flags */ |
42 |
+ |
static const int adjacent_trifl[8] = { |
43 |
+ |
0, /* forbidden diagonal */ |
44 |
+ |
1<<VDB_x|1<<VDB_y|1<<VDB_Xy, |
45 |
+ |
1<<VDB_y|1<<VDB_x|1<<VDB_xY, |
46 |
+ |
1<<VDB_y|1<<VDB_Xy|1<<VDB_X, |
47 |
+ |
1<<VDB_x|1<<VDB_xY|1<<VDB_Y, |
48 |
+ |
1<<VDB_Xy|1<<VDB_X|1<<VDB_Y, |
49 |
+ |
1<<VDB_xY|1<<VDB_Y|1<<VDB_X, |
50 |
+ |
0, /* forbidden diagonal */ |
51 |
+ |
}; |
52 |
+ |
|
53 |
|
typedef struct { |
54 |
+ |
COLOR v; /* hemisphere sample value */ |
55 |
+ |
float d; /* reciprocal distance (1/rt) */ |
56 |
+ |
FVECT p; /* intersection point */ |
57 |
+ |
} AMBSAMP; /* sample value */ |
58 |
+ |
|
59 |
+ |
typedef struct { |
60 |
|
RAY *rp; /* originating ray sample */ |
61 |
|
FVECT ux, uy; /* tangent axis unit vectors */ |
62 |
|
int ns; /* number of samples per axis */ |
63 |
|
COLOR acoef; /* division contribution coefficient */ |
64 |
< |
struct s_ambsamp { |
30 |
< |
COLOR v; /* hemisphere sample value */ |
31 |
< |
FVECT p; /* intersection point */ |
32 |
< |
} sa[1]; /* sample array (extends struct) */ |
64 |
> |
AMBSAMP sa[1]; /* sample array (extends struct) */ |
65 |
|
} AMBHEMI; /* ambient sample hemisphere */ |
66 |
|
|
67 |
< |
typedef struct s_ambsamp AMBSAMP; |
67 |
> |
#define ambndx(h,i,j) ((i)*(h)->ns + (j)) |
68 |
> |
#define ambsam(h,i,j) (h)->sa[ambndx(h,i,j)] |
69 |
|
|
37 |
– |
#define ambsamp(h,i,j) (h)->sa[(i)*(h)->ns + (j)] |
38 |
– |
|
70 |
|
typedef struct { |
71 |
|
FVECT r_i, r_i1, e_i, rcp, rI2_eJ2; |
72 |
|
double I1, I2; |
73 |
+ |
int valid; |
74 |
|
} FFTRI; /* vectors and coefficients for Hessian calculation */ |
75 |
|
|
76 |
|
|
77 |
+ |
/* Get index for adjacent vertex */ |
78 |
+ |
static int |
79 |
+ |
adjacent_verti(AMBHEMI *hp, int i, int j, int dbit) |
80 |
+ |
{ |
81 |
+ |
int i0 = i*hp->ns + j; |
82 |
+ |
|
83 |
+ |
switch (dbit) { |
84 |
+ |
case VDB_y: return(i0 - hp->ns); |
85 |
+ |
case VDB_x: return(i0 - 1); |
86 |
+ |
case VDB_Xy: return(i0 - hp->ns + 1); |
87 |
+ |
case VDB_xY: return(i0 + hp->ns - 1); |
88 |
+ |
case VDB_X: return(i0 + 1); |
89 |
+ |
case VDB_Y: return(i0 + hp->ns); |
90 |
+ |
/* the following should never occur */ |
91 |
+ |
case VDB_xy: return(i0 - hp->ns - 1); |
92 |
+ |
case VDB_XY: return(i0 + hp->ns + 1); |
93 |
+ |
} |
94 |
+ |
return(-1); |
95 |
+ |
} |
96 |
+ |
|
97 |
+ |
|
98 |
+ |
/* Get vertex direction bit for the opposite edge to complete triangle */ |
99 |
+ |
static int |
100 |
+ |
vdb_edge(int db1, int db2) |
101 |
+ |
{ |
102 |
+ |
switch (db1) { |
103 |
+ |
case VDB_x: return(db2==VDB_y ? VDB_Xy : VDB_Y); |
104 |
+ |
case VDB_y: return(db2==VDB_x ? VDB_xY : VDB_X); |
105 |
+ |
case VDB_X: return(db2==VDB_Xy ? VDB_y : VDB_xY); |
106 |
+ |
case VDB_Y: return(db2==VDB_xY ? VDB_x : VDB_Xy); |
107 |
+ |
case VDB_xY: return(db2==VDB_x ? VDB_y : VDB_X); |
108 |
+ |
case VDB_Xy: return(db2==VDB_y ? VDB_x : VDB_Y); |
109 |
+ |
} |
110 |
+ |
error(INTERNAL, "forbidden diagonal in vdb_edge()"); |
111 |
+ |
return(-1); |
112 |
+ |
} |
113 |
+ |
|
114 |
+ |
|
115 |
|
static AMBHEMI * |
116 |
|
inithemi( /* initialize sampling hemisphere */ |
117 |
|
COLOR ac, |
158 |
|
} |
159 |
|
|
160 |
|
|
161 |
< |
/* Prepare ambient division sample */ |
161 |
> |
/* Sample ambient division and apply weighting coefficient */ |
162 |
|
static int |
163 |
< |
prepambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n) |
163 |
> |
getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n) |
164 |
|
{ |
165 |
|
int hlist[3], ii; |
166 |
|
double spt[2], zd; |
176 |
|
scalecolor(arp->rcoef, 1./AVGREFL); |
177 |
|
} |
178 |
|
hlist[0] = hp->rp->rno; |
179 |
< |
hlist[1] = i; |
180 |
< |
hlist[2] = j; |
179 |
> |
hlist[1] = j; |
180 |
> |
hlist[2] = i; |
181 |
|
multisamp(spt, 2, urand(ilhash(hlist,3)+n)); |
182 |
|
if (!n) { /* avoid border samples for n==0 */ |
183 |
< |
if ((spt[0] < 0.1) | (spt[0] > 0.9)) |
183 |
> |
if ((spt[0] < 0.1) | (spt[0] >= 0.9)) |
184 |
|
spt[0] = 0.1 + 0.8*frandom(); |
185 |
< |
if ((spt[1] < 0.1) | (spt[1] > 0.9)) |
185 |
> |
if ((spt[1] < 0.1) | (spt[1] >= 0.9)) |
186 |
|
spt[1] = 0.1 + 0.8*frandom(); |
187 |
|
} |
188 |
< |
SDsquare2disk(spt, (i+spt[0])/hp->ns, (j+spt[1])/hp->ns); |
188 |
> |
SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns); |
189 |
|
zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]); |
190 |
|
for (ii = 3; ii--; ) |
191 |
|
arp->rdir[ii] = spt[0]*hp->ux[ii] + |
192 |
|
spt[1]*hp->uy[ii] + |
193 |
|
zd*hp->rp->ron[ii]; |
194 |
|
checknorm(arp->rdir); |
195 |
+ |
dimlist[ndims++] = ambndx(hp,i,j) + 90171; |
196 |
+ |
rayvalue(arp); /* evaluate ray */ |
197 |
+ |
ndims--; /* apply coefficient */ |
198 |
+ |
multcolor(arp->rcol, arp->rcoef); |
199 |
|
return(1); |
200 |
|
} |
201 |
|
|
202 |
|
|
203 |
|
static AMBSAMP * |
204 |
< |
ambsample( /* sample an ambient direction */ |
204 |
> |
ambsample( /* initial ambient division sample */ |
205 |
|
AMBHEMI *hp, |
206 |
|
int i, |
207 |
|
int j |
208 |
|
) |
209 |
|
{ |
210 |
< |
AMBSAMP *ap = &ambsamp(hp,i,j); |
210 |
> |
AMBSAMP *ap = &ambsam(hp,i,j); |
211 |
|
RAY ar; |
212 |
|
/* generate hemispherical sample */ |
213 |
< |
if (!prepambsamp(&ar, hp, i, j, 0)) |
214 |
< |
goto badsample; |
215 |
< |
dimlist[ndims++] = i*hp->ns + j + 90171; |
216 |
< |
rayvalue(&ar); /* evaluate ray */ |
217 |
< |
ndims--; |
144 |
< |
/* limit vertex distance */ |
213 |
> |
if (!getambsamp(&ar, hp, i, j, 0) || ar.rt <= FTINY) { |
214 |
> |
memset(ap, 0, sizeof(AMBSAMP)); |
215 |
> |
return(NULL); |
216 |
> |
} |
217 |
> |
ap->d = 1.0/ar.rt; /* limit vertex distance */ |
218 |
|
if (ar.rt > 10.0*thescene.cusize) |
219 |
|
ar.rt = 10.0*thescene.cusize; |
147 |
– |
else if (ar.rt <= FTINY) /* should never happen! */ |
148 |
– |
goto badsample; |
220 |
|
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
150 |
– |
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
221 |
|
copycolor(ap->v, ar.rcol); |
222 |
|
return(ap); |
153 |
– |
badsample: |
154 |
– |
setcolor(ap->v, 0., 0., 0.); |
155 |
– |
VCOPY(ap->p, hp->rp->rop); |
156 |
– |
return(NULL); |
223 |
|
} |
224 |
|
|
225 |
|
|
227 |
|
static float * |
228 |
|
getambdiffs(AMBHEMI *hp) |
229 |
|
{ |
230 |
< |
float *earr = calloc(hp->ns*hp->ns, sizeof(float)); |
230 |
> |
float *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float)); |
231 |
|
float *ep; |
232 |
|
AMBSAMP *ap; |
233 |
|
double b, d2; |
269 |
|
} |
270 |
|
|
271 |
|
|
272 |
< |
/* Perform super-sampling on hemisphere */ |
272 |
> |
/* Perform super-sampling on hemisphere (introduces bias) */ |
273 |
|
static void |
274 |
|
ambsupersamp(double acol[3], AMBHEMI *hp, int cnt) |
275 |
|
{ |
292 |
|
int nss = *ep/e2sum*cnt + frandom(); |
293 |
|
setcolor(asum, 0., 0., 0.); |
294 |
|
for (n = 1; n <= nss; n++) { |
295 |
< |
if (!prepambsamp(&ar, hp, i, j, n)) { |
295 |
> |
if (!getambsamp(&ar, hp, i, j, n)) { |
296 |
|
nss = n-1; |
297 |
|
break; |
298 |
|
} |
233 |
– |
dimlist[ndims++] = i*hp->ns + j + 90171; |
234 |
– |
rayvalue(&ar); /* evaluate super-sample */ |
235 |
– |
ndims--; |
236 |
– |
multcolor(ar.rcol, ar.rcoef); |
299 |
|
addcolor(asum, ar.rcol); |
300 |
|
} |
301 |
|
if (nss) { /* update returned ambient value */ |
311 |
|
} |
312 |
|
|
313 |
|
|
314 |
+ |
/* Compute vertex flags, indicating farthest in each direction */ |
315 |
+ |
static uby8 * |
316 |
+ |
vertex_flags(AMBHEMI *hp) |
317 |
+ |
{ |
318 |
+ |
uby8 *vflags = (uby8 *)calloc(hp->ns*hp->ns, sizeof(uby8)); |
319 |
+ |
uby8 *vf; |
320 |
+ |
AMBSAMP *ap; |
321 |
+ |
int i, j; |
322 |
+ |
|
323 |
+ |
if (vflags == NULL) |
324 |
+ |
error(SYSTEM, "out of memory in vertex_flags()"); |
325 |
+ |
vf = vflags; |
326 |
+ |
ap = hp->sa; /* compute farthest along first row */ |
327 |
+ |
for (j = 0; j < hp->ns-1; j++, vf++, ap++) |
328 |
+ |
if (ap[0].d <= ap[1].d) |
329 |
+ |
vf[0] |= 1<<VDB_X; |
330 |
+ |
else |
331 |
+ |
vf[1] |= 1<<VDB_x; |
332 |
+ |
++vf; ++ap; |
333 |
+ |
/* flag subsequent rows */ |
334 |
+ |
for (i = 1; i < hp->ns; i++) { |
335 |
+ |
for (j = 0; j < hp->ns-1; j++, vf++, ap++) { |
336 |
+ |
if (ap[0].d <= ap[-hp->ns].d) /* row before */ |
337 |
+ |
vf[0] |= 1<<VDB_y; |
338 |
+ |
else |
339 |
+ |
vf[-hp->ns] |= 1<<VDB_Y; |
340 |
+ |
if (ap[0].d <= ap[1-hp->ns].d) /* diagonal we care about */ |
341 |
+ |
vf[0] |= 1<<VDB_Xy; |
342 |
+ |
else |
343 |
+ |
vf[1-hp->ns] |= 1<<VDB_xY; |
344 |
+ |
if (ap[0].d <= ap[1].d) /* column after */ |
345 |
+ |
vf[0] |= 1<<VDB_X; |
346 |
+ |
else |
347 |
+ |
vf[1] |= 1<<VDB_x; |
348 |
+ |
} |
349 |
+ |
if (ap[0].d <= ap[-hp->ns].d) /* final column edge */ |
350 |
+ |
vf[0] |= 1<<VDB_y; |
351 |
+ |
else |
352 |
+ |
vf[-hp->ns] |= 1<<VDB_Y; |
353 |
+ |
++vf; ++ap; |
354 |
+ |
} |
355 |
+ |
return(vflags); |
356 |
+ |
} |
357 |
+ |
|
358 |
+ |
|
359 |
+ |
/* Return brightness of farthest ambient sample */ |
360 |
+ |
static double |
361 |
+ |
back_ambval(AMBHEMI *hp, int i, int j, int dbit1, int dbit2, const uby8 *vflags) |
362 |
+ |
{ |
363 |
+ |
const int v0 = ambndx(hp,i,j); |
364 |
+ |
const int tflags = (1<<dbit1 | 1<<dbit2); |
365 |
+ |
int v1, v2; |
366 |
+ |
|
367 |
+ |
if ((vflags[v0] & tflags) == tflags) /* is v0 the farthest? */ |
368 |
+ |
return(colval(hp->sa[v0].v,CIEY)); |
369 |
+ |
v1 = adjacent_verti(hp, i, j, dbit1); |
370 |
+ |
if (vflags[v0] & 1<<dbit2) /* v1 farthest if v0>v2 */ |
371 |
+ |
return(colval(hp->sa[v1].v,CIEY)); |
372 |
+ |
v2 = adjacent_verti(hp, i, j, dbit2); |
373 |
+ |
if (vflags[v0] & 1<<dbit1) /* v2 farthest if v0>v1 */ |
374 |
+ |
return(colval(hp->sa[v2].v,CIEY)); |
375 |
+ |
/* else check if v1>v2 */ |
376 |
+ |
if (vflags[v1] & 1<<vdb_edge(dbit1,dbit2)) |
377 |
+ |
return(colval(hp->sa[v1].v,CIEY)); |
378 |
+ |
return(colval(hp->sa[v2].v,CIEY)); |
379 |
+ |
} |
380 |
+ |
|
381 |
+ |
|
382 |
|
/* Compute vectors and coefficients for Hessian/gradient calcs */ |
383 |
|
static void |
384 |
< |
comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop) |
384 |
> |
comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int dbit, const uby8 *vflags) |
385 |
|
{ |
386 |
< |
double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2; |
387 |
< |
int i; |
386 |
> |
const int i0 = ambndx(hp,i,j); |
387 |
> |
double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2; |
388 |
> |
int i1, ii; |
389 |
|
|
390 |
< |
VSUB(ftp->r_i, ap0, rop); |
391 |
< |
VSUB(ftp->r_i1, ap1, rop); |
392 |
< |
VSUB(ftp->e_i, ap1, ap0); |
390 |
> |
ftp->valid = 0; /* check if we can skip this edge */ |
391 |
> |
ii = adjacent_trifl[dbit]; |
392 |
> |
if ((vflags[i0] & ii) == ii) /* cancels if vertex used as value */ |
393 |
> |
return; |
394 |
> |
i1 = adjacent_verti(hp, i, j, dbit); |
395 |
> |
ii = adjacent_trifl[VDB_OPP(dbit)]; |
396 |
> |
if ((vflags[i1] & ii) == ii) /* on either end (for both triangles) */ |
397 |
> |
return; |
398 |
> |
/* else go ahead with calculation */ |
399 |
> |
VSUB(ftp->r_i, hp->sa[i0].p, hp->rp->rop); |
400 |
> |
VSUB(ftp->r_i1, hp->sa[i1].p, hp->rp->rop); |
401 |
> |
VSUB(ftp->e_i, hp->sa[i1].p, hp->sa[i0].p); |
402 |
|
VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1); |
403 |
|
rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp); |
404 |
|
dot_e = DOT(ftp->e_i,ftp->e_i); |
410 |
|
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r + |
411 |
|
dot_e*ftp->I1 )*0.5*rdot_cp; |
412 |
|
J2 = ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e; |
413 |
< |
for (i = 3; i--; ) |
414 |
< |
ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; |
413 |
> |
for (ii = 3; ii--; ) |
414 |
> |
ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii]; |
415 |
> |
ftp->valid++; |
416 |
|
} |
417 |
|
|
418 |
|
|
438 |
|
double d1, d2, d3, d4; |
439 |
|
double I3, J3, K3; |
440 |
|
int i, j; |
441 |
+ |
|
442 |
+ |
if (!ftp->valid) { /* preemptive test */ |
443 |
+ |
memset(hess, 0, sizeof(FVECT)*3); |
444 |
+ |
return; |
445 |
+ |
} |
446 |
|
/* compute intermediate coefficients */ |
447 |
|
d1 = 1.0/DOT(ftp->r_i,ftp->r_i); |
448 |
|
d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1); |
466 |
|
hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] + |
467 |
|
2.0*J3*m4[i][j] ); |
468 |
|
hess[i][j] += d2*(i==j); |
469 |
< |
hess[i][j] *= 1.0/PI; |
469 |
> |
hess[i][j] *= -1.0/PI; |
470 |
|
} |
471 |
|
} |
472 |
|
|
488 |
|
/* Add to radiometric Hessian from the given triangle */ |
489 |
|
static void |
490 |
|
add2hessian(FVECT hess[3], FVECT ehess1[3], |
491 |
< |
FVECT ehess2[3], FVECT ehess3[3], COLORV v) |
491 |
> |
FVECT ehess2[3], FVECT ehess3[3], double v) |
492 |
|
{ |
493 |
|
int i, j; |
494 |
|
|
506 |
|
double f1; |
507 |
|
int i; |
508 |
|
|
509 |
+ |
if (!ftp->valid) { /* preemptive test */ |
510 |
+ |
memset(grad, 0, sizeof(FVECT)); |
511 |
+ |
return; |
512 |
+ |
} |
513 |
|
f1 = 2.0*DOT(nrm, ftp->rcp); |
514 |
|
VCROSS(ncp, nrm, ftp->e_i); |
515 |
|
for (i = 3; i--; ) |
516 |
< |
grad[i] = (-0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] ); |
516 |
> |
grad[i] = (0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] ); |
517 |
|
} |
518 |
|
|
519 |
|
|
529 |
|
|
530 |
|
/* Add to displacement gradient from the given triangle */ |
531 |
|
static void |
532 |
< |
add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v) |
532 |
> |
add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v) |
533 |
|
{ |
534 |
|
int i; |
535 |
|
|
538 |
|
} |
539 |
|
|
540 |
|
|
391 |
– |
/* Return brightness of furthest ambient sample */ |
392 |
– |
static COLORV |
393 |
– |
back_ambval(AMBSAMP *ap1, AMBSAMP *ap2, AMBSAMP *ap3, FVECT orig) |
394 |
– |
{ |
395 |
– |
COLORV vback; |
396 |
– |
FVECT vec; |
397 |
– |
double d2, d2best; |
398 |
– |
|
399 |
– |
VSUB(vec, ap1->p, orig); |
400 |
– |
d2best = DOT(vec,vec); |
401 |
– |
vback = colval(ap1->v,CIEY); |
402 |
– |
VSUB(vec, ap2->p, orig); |
403 |
– |
d2 = DOT(vec,vec); |
404 |
– |
if (d2 > d2best) { |
405 |
– |
d2best = d2; |
406 |
– |
vback = colval(ap2->v,CIEY); |
407 |
– |
} |
408 |
– |
VSUB(vec, ap3->p, orig); |
409 |
– |
d2 = DOT(vec,vec); |
410 |
– |
if (d2 > d2best) |
411 |
– |
return(colval(ap3->v,CIEY)); |
412 |
– |
return(vback); |
413 |
– |
} |
414 |
– |
|
415 |
– |
|
541 |
|
/* Compute anisotropic radii and eigenvector directions */ |
542 |
|
static int |
543 |
|
eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3]) |
598 |
|
static char memerrmsg[] = "out of memory in ambHessian()"; |
599 |
|
FVECT (*hessrow)[3] = NULL; |
600 |
|
FVECT *gradrow = NULL; |
601 |
+ |
uby8 *vflags; |
602 |
|
FVECT hessian[3]; |
603 |
|
FVECT gradient; |
604 |
|
FFTRI fftr; |
620 |
|
error(SYSTEM, memerrmsg); |
621 |
|
memset(gradient, 0, sizeof(gradient)); |
622 |
|
} |
623 |
+ |
/* get vertex position flags */ |
624 |
+ |
vflags = vertex_flags(hp); |
625 |
|
/* compute first row of edges */ |
626 |
|
for (j = 0; j < hp->ns-1; j++) { |
627 |
< |
comp_fftri(&fftr, ambsamp(hp,0,j).p, |
500 |
< |
ambsamp(hp,0,j+1).p, hp->rp->rop); |
627 |
> |
comp_fftri(&fftr, hp, 0, j, VDB_X, vflags); |
628 |
|
if (hessrow != NULL) |
629 |
|
comp_hessian(hessrow[j], &fftr, hp->rp->ron); |
630 |
|
if (gradrow != NULL) |
634 |
|
for (i = 0; i < hp->ns-1; i++) { |
635 |
|
FVECT hesscol[3]; /* compute first vertical edge */ |
636 |
|
FVECT gradcol; |
637 |
< |
comp_fftri(&fftr, ambsamp(hp,i,0).p, |
511 |
< |
ambsamp(hp,i+1,0).p, hp->rp->rop); |
637 |
> |
comp_fftri(&fftr, hp, i, 0, VDB_Y, vflags); |
638 |
|
if (hessrow != NULL) |
639 |
|
comp_hessian(hesscol, &fftr, hp->rp->ron); |
640 |
|
if (gradrow != NULL) |
642 |
|
for (j = 0; j < hp->ns-1; j++) { |
643 |
|
FVECT hessdia[3]; /* compute triangle contributions */ |
644 |
|
FVECT graddia; |
645 |
< |
COLORV backg; |
646 |
< |
backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1), |
521 |
< |
&ambsamp(hp,i+1,j), hp->rp->rop); |
645 |
> |
double backg; |
646 |
> |
backg = back_ambval(hp, i, j, VDB_X, VDB_Y, vflags); |
647 |
|
/* diagonal (inner) edge */ |
648 |
< |
comp_fftri(&fftr, ambsamp(hp,i,j+1).p, |
524 |
< |
ambsamp(hp,i+1,j).p, hp->rp->rop); |
648 |
> |
comp_fftri(&fftr, hp, i, j+1, VDB_xY, vflags); |
649 |
|
if (hessrow != NULL) { |
650 |
|
comp_hessian(hessdia, &fftr, hp->rp->ron); |
651 |
|
rev_hessian(hesscol); |
657 |
|
add2gradient(gradient, gradrow[j], graddia, gradcol, backg); |
658 |
|
} |
659 |
|
/* initialize edge in next row */ |
660 |
< |
comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p, |
537 |
< |
ambsamp(hp,i+1,j).p, hp->rp->rop); |
660 |
> |
comp_fftri(&fftr, hp, i+1, j+1, VDB_x, vflags); |
661 |
|
if (hessrow != NULL) |
662 |
|
comp_hessian(hessrow[j], &fftr, hp->rp->ron); |
663 |
|
if (gradrow != NULL) |
664 |
|
comp_gradient(gradrow[j], &fftr, hp->rp->ron); |
665 |
|
/* new column edge & paired triangle */ |
666 |
< |
backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1), |
667 |
< |
&ambsamp(hp,i+1,j), hp->rp->rop); |
545 |
< |
comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p, |
546 |
< |
hp->rp->rop); |
666 |
> |
backg = back_ambval(hp, i+1, j+1, VDB_x, VDB_y, vflags); |
667 |
> |
comp_fftri(&fftr, hp, i, j+1, VDB_Y, vflags); |
668 |
|
if (hessrow != NULL) { |
669 |
|
comp_hessian(hesscol, &fftr, hp->rp->ron); |
670 |
|
rev_hessian(hessdia); |
684 |
|
/* release row buffers */ |
685 |
|
if (hessrow != NULL) free(hessrow); |
686 |
|
if (gradrow != NULL) free(gradrow); |
687 |
+ |
free(vflags); |
688 |
|
|
689 |
|
if (ra != NULL) /* extract eigenvectors & radii */ |
690 |
|
eigenvectors(uv, ra, hessian); |
720 |
|
} |
721 |
|
|
722 |
|
|
723 |
+ |
/* Make sure radii don't extend beyond what we see in our periphery */ |
724 |
+ |
static void |
725 |
+ |
hem_radii(AMBHEMI *hp, FVECT uv[2], float ra[2]) |
726 |
+ |
{ |
727 |
+ |
#ifdef HEM_MULT |
728 |
+ |
double udsum = 0, vdsum = 0; |
729 |
+ |
double uwsum = 0, vwsum = 0; |
730 |
+ |
int i, j; |
731 |
+ |
/* circle around perimeter */ |
732 |
+ |
for (i = 0; i < hp->ns; i++) |
733 |
+ |
for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) { |
734 |
+ |
AMBSAMP *ap = &ambsam(hp,i,j); |
735 |
+ |
FVECT vec; |
736 |
+ |
double us2, vs2; |
737 |
+ |
VSUB(vec, ap->p, hp->rp->rop); |
738 |
+ |
us2 = DOT(vec, uv[0]) * ap->d; |
739 |
+ |
us2 *= us2; |
740 |
+ |
vs2 = DOT(vec, uv[1]) * ap->d; |
741 |
+ |
vs2 *= vs2; |
742 |
+ |
udsum += us2 * ap->d; |
743 |
+ |
uwsum += us2; |
744 |
+ |
vdsum += vs2 * ap->d; |
745 |
+ |
vwsum += vs2; |
746 |
+ |
} |
747 |
+ |
uwsum *= HEM_MULT; /* adjust effective hem size */ |
748 |
+ |
vwsum *= HEM_MULT; |
749 |
+ |
/* cap radii (recall d=1/rt) */ |
750 |
+ |
if (ra[0]*udsum > uwsum) |
751 |
+ |
ra[0] = uwsum/udsum; |
752 |
+ |
if (ra[1]*vdsum > vwsum) |
753 |
+ |
ra[1] = vwsum/vdsum; |
754 |
+ |
#endif |
755 |
+ |
} |
756 |
+ |
|
757 |
+ |
|
758 |
|
int |
759 |
|
doambient( /* compute ambient component */ |
760 |
|
COLOR rcol, /* input/output color */ |
767 |
|
) |
768 |
|
{ |
769 |
|
AMBHEMI *hp = inithemi(rcol, r, wt); |
770 |
< |
int cnt = 0; |
770 |
> |
int cnt; |
771 |
|
FVECT my_uv[2]; |
772 |
|
double d, K, acol[3]; |
773 |
|
AMBSAMP *ap; |
785 |
|
dg[0] = dg[1] = 0.0; |
786 |
|
/* sample the hemisphere */ |
787 |
|
acol[0] = acol[1] = acol[2] = 0.0; |
788 |
+ |
cnt = 0; |
789 |
|
for (i = hp->ns; i--; ) |
790 |
|
for (j = hp->ns; j--; ) |
791 |
|
if ((ap = ambsample(hp, i, j)) != NULL) { |
810 |
|
free(hp); |
811 |
|
return(-1); /* no radius or gradient calc. */ |
812 |
|
} |
813 |
< |
if (bright(acol) > FTINY) { /* normalize Y values */ |
814 |
< |
d = 0.99*cnt/bright(acol); |
813 |
> |
if ((d = bright(acol)) > FTINY) { /* normalize Y values */ |
814 |
> |
d = 0.99*(hp->ns*hp->ns)/d; |
815 |
|
K = 0.01; |
816 |
< |
} else { /* geometric Hessian fall-back */ |
659 |
< |
d = 0.0; |
816 |
> |
} else { /* or fall back on geometric Hessian */ |
817 |
|
K = 1.0; |
818 |
|
pg = NULL; |
819 |
|
dg = NULL; |
836 |
|
ra[0] = 1.0/d; |
837 |
|
if (ra[1]*(d = fabs(pg[1])) > 1.0) |
838 |
|
ra[1] = 1.0/d; |
682 |
– |
if (ra[0] > ra[1]) |
683 |
– |
ra[0] = ra[1]; |
839 |
|
} |
840 |
+ |
hem_radii(hp, uv, ra); |
841 |
+ |
if (ra[0] > ra[1]) |
842 |
+ |
ra[0] = ra[1]; |
843 |
|
if (ra[0] < minarad) { |
844 |
|
ra[0] = minarad; |
845 |
|
if (ra[1] < minarad) |