22 |
|
extern void SDsquare2disk(double ds[2], double seedx, double seedy); |
23 |
|
|
24 |
|
typedef struct { |
25 |
+ |
COLOR v; /* hemisphere sample value */ |
26 |
+ |
FVECT p; /* intersection point */ |
27 |
+ |
} AMBSAMP; /* sample value */ |
28 |
+ |
|
29 |
+ |
typedef struct { |
30 |
|
RAY *rp; /* originating ray sample */ |
31 |
|
FVECT ux, uy; /* tangent axis unit vectors */ |
32 |
|
int ns; /* number of samples per axis */ |
33 |
|
COLOR acoef; /* division contribution coefficient */ |
34 |
< |
struct s_ambsamp { |
30 |
< |
COLOR v; /* hemisphere sample value */ |
31 |
< |
FVECT p; /* intersection point */ |
32 |
< |
} sa[1]; /* sample array (extends struct) */ |
34 |
> |
AMBSAMP sa[1]; /* sample array (extends struct) */ |
35 |
|
} AMBHEMI; /* ambient sample hemisphere */ |
36 |
|
|
37 |
< |
typedef struct s_ambsamp AMBSAMP; |
37 |
> |
#define ambsam(h,i,j) (h)->sa[(i)*(h)->ns + (j)] |
38 |
|
|
37 |
– |
#define ambsamp(h,i,j) (h)->sa[(i)*(h)->ns + (j)] |
38 |
– |
|
39 |
|
typedef struct { |
40 |
|
FVECT r_i, r_i1, e_i, rcp, rI2_eJ2; |
41 |
|
double I1, I2; |
88 |
|
} |
89 |
|
|
90 |
|
|
91 |
< |
/* Prepare ambient division sample */ |
91 |
> |
/* Sample ambient division and apply weighting coefficient */ |
92 |
|
static int |
93 |
< |
prepambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n) |
93 |
> |
getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n) |
94 |
|
{ |
95 |
|
int hlist[3], ii; |
96 |
|
double spt[2], zd; |
122 |
|
spt[1]*hp->uy[ii] + |
123 |
|
zd*hp->rp->ron[ii]; |
124 |
|
checknorm(arp->rdir); |
125 |
+ |
dimlist[ndims++] = i*hp->ns + j + 90171; |
126 |
+ |
rayvalue(arp); /* evaluate ray */ |
127 |
+ |
ndims--; /* apply coefficient */ |
128 |
+ |
multcolor(arp->rcol, arp->rcoef); |
129 |
|
return(1); |
130 |
|
} |
131 |
|
|
132 |
|
|
133 |
|
static AMBSAMP * |
134 |
< |
ambsample( /* sample an ambient direction */ |
134 |
> |
ambsample( /* initial ambient division sample */ |
135 |
|
AMBHEMI *hp, |
136 |
|
int i, |
137 |
|
int j |
138 |
|
) |
139 |
|
{ |
140 |
< |
AMBSAMP *ap = &ambsamp(hp,i,j); |
140 |
> |
AMBSAMP *ap = &ambsam(hp,i,j); |
141 |
|
RAY ar; |
142 |
|
/* generate hemispherical sample */ |
143 |
< |
if (!prepambsamp(&ar, hp, i, j, 0)) |
143 |
> |
if (!getambsamp(&ar, hp, i, j, 0)) |
144 |
|
goto badsample; |
141 |
– |
dimlist[ndims++] = i*hp->ns + j + 90171; |
142 |
– |
rayvalue(&ar); /* evaluate ray */ |
143 |
– |
ndims--; |
145 |
|
/* limit vertex distance */ |
146 |
|
if (ar.rt > 10.0*thescene.cusize) |
147 |
|
ar.rt = 10.0*thescene.cusize; |
148 |
|
else if (ar.rt <= FTINY) /* should never happen! */ |
149 |
|
goto badsample; |
150 |
|
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
150 |
– |
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
151 |
|
copycolor(ap->v, ar.rcol); |
152 |
|
return(ap); |
153 |
|
badsample: |
161 |
|
static float * |
162 |
|
getambdiffs(AMBHEMI *hp) |
163 |
|
{ |
164 |
< |
float *earr = calloc(hp->ns*hp->ns, sizeof(float)); |
164 |
> |
float *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float)); |
165 |
|
float *ep; |
166 |
+ |
AMBSAMP *ap; |
167 |
|
double b, d2; |
168 |
|
int i, j; |
169 |
|
|
170 |
|
if (earr == NULL) /* out of memory? */ |
171 |
|
return(NULL); |
172 |
|
/* compute squared neighbor diffs */ |
173 |
< |
for (ep = earr, i = 0; i < hp->ns; i++) |
174 |
< |
for (j = 0; j < hp->ns; j++, ep++) { |
175 |
< |
b = bright(ambsamp(hp,i,j).v); |
173 |
> |
for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++) |
174 |
> |
for (j = 0; j < hp->ns; j++, ap++, ep++) { |
175 |
> |
b = bright(ap[0].v); |
176 |
|
if (i) { /* from above */ |
177 |
< |
d2 = b - bright(ambsamp(hp,i-1,j).v); |
177 |
> |
d2 = b - bright(ap[-hp->ns].v); |
178 |
|
d2 *= d2; |
179 |
|
ep[0] += d2; |
180 |
|
ep[-hp->ns] += d2; |
181 |
|
} |
182 |
|
if (j) { /* from behind */ |
183 |
< |
d2 = b - bright(ambsamp(hp,i,j-1).v); |
183 |
> |
d2 = b - bright(ap[-1].v); |
184 |
|
d2 *= d2; |
185 |
|
ep[0] += d2; |
186 |
|
ep[-1] += d2; |
203 |
|
} |
204 |
|
|
205 |
|
|
206 |
< |
/* Perform super-sampling on hemisphere */ |
206 |
> |
/* Perform super-sampling on hemisphere (introduces bias) */ |
207 |
|
static void |
208 |
|
ambsupersamp(double acol[3], AMBHEMI *hp, int cnt) |
209 |
|
{ |
226 |
|
int nss = *ep/e2sum*cnt + frandom(); |
227 |
|
setcolor(asum, 0., 0., 0.); |
228 |
|
for (n = 1; n <= nss; n++) { |
229 |
< |
if (!prepambsamp(&ar, hp, i, j, n)) { |
229 |
> |
if (!getambsamp(&ar, hp, i, j, n)) { |
230 |
|
nss = n-1; |
231 |
|
break; |
232 |
|
} |
232 |
– |
dimlist[ndims++] = i*hp->ns + j + 90171; |
233 |
– |
rayvalue(&ar); /* evaluate super-sample */ |
234 |
– |
ndims--; |
235 |
– |
multcolor(ar.rcol, ar.rcoef); |
233 |
|
addcolor(asum, ar.rcol); |
234 |
|
} |
235 |
|
if (nss) { /* update returned ambient value */ |
492 |
|
} |
493 |
|
/* compute first row of edges */ |
494 |
|
for (j = 0; j < hp->ns-1; j++) { |
495 |
< |
comp_fftri(&fftr, ambsamp(hp,0,j).p, |
496 |
< |
ambsamp(hp,0,j+1).p, hp->rp->rop); |
495 |
> |
comp_fftri(&fftr, ambsam(hp,0,j).p, |
496 |
> |
ambsam(hp,0,j+1).p, hp->rp->rop); |
497 |
|
if (hessrow != NULL) |
498 |
|
comp_hessian(hessrow[j], &fftr, hp->rp->ron); |
499 |
|
if (gradrow != NULL) |
503 |
|
for (i = 0; i < hp->ns-1; i++) { |
504 |
|
FVECT hesscol[3]; /* compute first vertical edge */ |
505 |
|
FVECT gradcol; |
506 |
< |
comp_fftri(&fftr, ambsamp(hp,i,0).p, |
507 |
< |
ambsamp(hp,i+1,0).p, hp->rp->rop); |
506 |
> |
comp_fftri(&fftr, ambsam(hp,i,0).p, |
507 |
> |
ambsam(hp,i+1,0).p, hp->rp->rop); |
508 |
|
if (hessrow != NULL) |
509 |
|
comp_hessian(hesscol, &fftr, hp->rp->ron); |
510 |
|
if (gradrow != NULL) |
513 |
|
FVECT hessdia[3]; /* compute triangle contributions */ |
514 |
|
FVECT graddia; |
515 |
|
COLORV backg; |
516 |
< |
backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1), |
517 |
< |
&ambsamp(hp,i+1,j), hp->rp->rop); |
516 |
> |
backg = back_ambval(&ambsam(hp,i,j), &ambsam(hp,i,j+1), |
517 |
> |
&ambsam(hp,i+1,j), hp->rp->rop); |
518 |
|
/* diagonal (inner) edge */ |
519 |
< |
comp_fftri(&fftr, ambsamp(hp,i,j+1).p, |
520 |
< |
ambsamp(hp,i+1,j).p, hp->rp->rop); |
519 |
> |
comp_fftri(&fftr, ambsam(hp,i,j+1).p, |
520 |
> |
ambsam(hp,i+1,j).p, hp->rp->rop); |
521 |
|
if (hessrow != NULL) { |
522 |
|
comp_hessian(hessdia, &fftr, hp->rp->ron); |
523 |
|
rev_hessian(hesscol); |
529 |
|
add2gradient(gradient, gradrow[j], graddia, gradcol, backg); |
530 |
|
} |
531 |
|
/* initialize edge in next row */ |
532 |
< |
comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p, |
533 |
< |
ambsamp(hp,i+1,j).p, hp->rp->rop); |
532 |
> |
comp_fftri(&fftr, ambsam(hp,i+1,j+1).p, |
533 |
> |
ambsam(hp,i+1,j).p, hp->rp->rop); |
534 |
|
if (hessrow != NULL) |
535 |
|
comp_hessian(hessrow[j], &fftr, hp->rp->ron); |
536 |
|
if (gradrow != NULL) |
537 |
|
comp_gradient(gradrow[j], &fftr, hp->rp->ron); |
538 |
|
/* new column edge & paired triangle */ |
539 |
< |
backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1), |
540 |
< |
&ambsamp(hp,i+1,j), hp->rp->rop); |
541 |
< |
comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p, |
539 |
> |
backg = back_ambval(&ambsam(hp,i,j+1), &ambsam(hp,i+1,j+1), |
540 |
> |
&ambsam(hp,i+1,j), hp->rp->rop); |
541 |
> |
comp_fftri(&fftr, ambsam(hp,i,j+1).p, ambsam(hp,i+1,j+1).p, |
542 |
|
hp->rp->rop); |
543 |
|
if (hessrow != NULL) { |
544 |
|
comp_hessian(hesscol, &fftr, hp->rp->ron); |
606 |
|
) |
607 |
|
{ |
608 |
|
AMBHEMI *hp = inithemi(rcol, r, wt); |
609 |
< |
int cnt = 0; |
609 |
> |
int cnt; |
610 |
|
FVECT my_uv[2]; |
611 |
|
double d, K, acol[3]; |
612 |
|
AMBSAMP *ap; |
624 |
|
dg[0] = dg[1] = 0.0; |
625 |
|
/* sample the hemisphere */ |
626 |
|
acol[0] = acol[1] = acol[2] = 0.0; |
627 |
+ |
cnt = 0; |
628 |
|
for (i = hp->ns; i--; ) |
629 |
|
for (j = hp->ns; j--; ) |
630 |
|
if ((ap = ambsample(hp, i, j)) != NULL) { |
649 |
|
free(hp); |
650 |
|
return(-1); /* no radius or gradient calc. */ |
651 |
|
} |
652 |
< |
if (bright(acol) > FTINY) { /* normalize Y values */ |
653 |
< |
d = 0.99*cnt/bright(acol); |
652 |
> |
if ((d = bright(acol)) > FTINY) { /* normalize Y values */ |
653 |
> |
d = 0.99*(hp->ns*hp->ns)/d; |
654 |
|
K = 0.01; |
655 |
< |
} else { /* geometric Hessian fall-back */ |
658 |
< |
d = 0.0; |
655 |
> |
} else { /* or fall back on geometric Hessian */ |
656 |
|
K = 1.0; |
657 |
|
pg = NULL; |
658 |
|
dg = NULL; |