| 156 |
|
} |
| 157 |
|
|
| 158 |
|
|
| 159 |
< |
/* Estimate errors based on ambient division differences */ |
| 159 |
> |
/* Estimate variance based on relative ambient division differences */ |
| 160 |
|
static float * |
| 161 |
|
getambdiffs(AMBHEMI *hp) |
| 162 |
|
{ |
| 164 |
|
float *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float)); |
| 165 |
|
float *ep; |
| 166 |
|
AMBSAMP *ap; |
| 167 |
< |
double b, d2; |
| 167 |
> |
double b, b1, d2; |
| 168 |
|
int i, j; |
| 169 |
|
|
| 170 |
|
if (earr == NULL) /* out of memory? */ |
| 171 |
|
return(NULL); |
| 172 |
< |
/* compute squared neighbor diffs */ |
| 172 |
> |
/* sum squared neighbor diffs */ |
| 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 = normf*(b - bright(ap[-hp->ns].v)); |
| 178 |
< |
d2 *= d2; |
| 177 |
> |
b1 = bright(ap[-hp->ns].v); |
| 178 |
> |
d2 = (b - b1)/(b + b1); |
| 179 |
> |
d2 *= d2*normf; |
| 180 |
|
ep[0] += d2; |
| 181 |
|
ep[-hp->ns] += d2; |
| 182 |
|
} |
| 183 |
|
if (!j) continue; |
| 184 |
|
/* from behind */ |
| 185 |
< |
d2 = normf*(b - bright(ap[-1].v)); |
| 186 |
< |
d2 *= d2; |
| 185 |
> |
b1 = bright(ap[-1].v); |
| 186 |
> |
d2 = (b - b1)/(b + b1); |
| 187 |
> |
d2 *= d2*normf; |
| 188 |
|
ep[0] += d2; |
| 189 |
|
ep[-1] += d2; |
| 190 |
|
if (!i) continue; |
| 191 |
|
/* diagonal */ |
| 192 |
< |
d2 = normf*(b - bright(ap[-hp->ns-1].v)); |
| 193 |
< |
d2 *= d2; |
| 192 |
> |
b1 = bright(ap[-hp->ns-1].v); |
| 193 |
> |
d2 = (b - b1)/(b + b1); |
| 194 |
> |
d2 *= d2*normf; |
| 195 |
|
ep[0] += d2; |
| 196 |
|
ep[-hp->ns-1] += d2; |
| 197 |
|
} |
| 218 |
|
{ |
| 219 |
|
float *earr = getambdiffs(hp); |
| 220 |
|
double e2rem = 0; |
| 218 |
– |
AMBSAMP *ap; |
| 221 |
|
float *ep; |
| 222 |
|
int i, j, n, nss; |
| 223 |
|
|
| 227 |
|
for (ep = earr + hp->ns*hp->ns; ep > earr; ) |
| 228 |
|
e2rem += *--ep; |
| 229 |
|
ep = earr; /* perform super-sampling */ |
| 230 |
< |
for (ap = hp->sa, i = 0; i < hp->ns; i++) |
| 231 |
< |
for (j = 0; j < hp->ns; j++, ap++) { |
| 230 |
> |
for (i = 0; i < hp->ns; i++) |
| 231 |
> |
for (j = 0; j < hp->ns; j++) { |
| 232 |
|
if (e2rem <= FTINY) |
| 233 |
|
goto done; /* nothing left to do */ |
| 234 |
|
nss = *ep/e2rem*cnt + frandom(); |
| 676 |
|
for (a1 = ang-ang_res; a1 <= ang+ang_res; a1 += ang_step) |
| 677 |
|
flgs |= 1L<<(int)(16/PI*(a1 + 2.*PI*(a1 < 0))); |
| 678 |
|
} |
| 677 |
– |
/* add low-angle incident (< 20deg) */ |
| 678 |
– |
if (fabs(hp->rp->rod) <= 0.342) { |
| 679 |
– |
u = -DOT(hp->rp->rdir, uv[0]); |
| 680 |
– |
v = -DOT(hp->rp->rdir, uv[1]); |
| 681 |
– |
if ((r0*r0*u*u + r1*r1*v*v) > hp->rp->rot*hp->rp->rot) { |
| 682 |
– |
ang = atan2a(v, u); |
| 683 |
– |
ang += 2.*PI*(ang < 0); |
| 684 |
– |
ang *= 16/PI; |
| 685 |
– |
if ((ang < .5) | (ang >= 31.5)) |
| 686 |
– |
flgs |= 0x80000001; |
| 687 |
– |
else |
| 688 |
– |
flgs |= 3L<<(int)(ang-.5); |
| 689 |
– |
} |
| 690 |
– |
} |
| 679 |
|
return(flgs); |
| 680 |
|
} |
| 681 |
|
|
