| 28 |
|
COLOR acoef; /* division contribution coefficient */ |
| 29 |
|
struct s_ambsamp { |
| 30 |
|
COLOR v; /* hemisphere sample value */ |
| 31 |
< |
float p[3]; /* intersection point */ |
| 31 |
> |
FVECT p; /* intersection point */ |
| 32 |
|
} sa[1]; /* sample array (extends struct) */ |
| 33 |
|
} AMBHEMI; /* ambient sample hemisphere */ |
| 34 |
|
|
| 35 |
|
#define ambsamp(h,i,j) (h)->sa[(i)*(h)->ns + (j)] |
| 36 |
|
|
| 37 |
|
typedef struct { |
| 38 |
< |
FVECT r_i, r_i1, e_i; |
| 39 |
< |
double nf, I1, I2, J2; |
| 38 |
> |
FVECT r_i, r_i1, e_i, rI2_eJ2; |
| 39 |
> |
double nf, I1, I2; |
| 40 |
|
} FFTRI; /* vectors and coefficients for Hessian calculation */ |
| 41 |
|
|
| 42 |
|
|
| 103 |
|
setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL); |
| 104 |
|
else |
| 105 |
|
copycolor(ar.rcoef, hp->acoef); |
| 106 |
< |
if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) { |
| 107 |
< |
setcolor(ap->v, 0., 0., 0.); |
| 108 |
< |
VCOPY(ap->p, hp->rp->rop); |
| 109 |
< |
return(NULL); /* no sample taken */ |
| 110 |
< |
} |
| 106 |
> |
if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) |
| 107 |
> |
goto badsample; |
| 108 |
|
if (ambacc > FTINY) { |
| 109 |
|
multcolor(ar.rcoef, hp->acoef); |
| 110 |
|
scalecolor(ar.rcoef, 1./AVGREFL); |
| 121 |
|
dimlist[ndims++] = i*hp->ns + j + 90171; |
| 122 |
|
rayvalue(&ar); /* evaluate ray */ |
| 123 |
|
ndims--; |
| 124 |
+ |
if (ar.rt > 20.0*maxarad) /* limit vertex distance */ |
| 125 |
+ |
ar.rt = 20.0*maxarad; |
| 126 |
+ |
else if (ar.rt <= FTINY) /* should never happen! */ |
| 127 |
+ |
goto badsample; |
| 128 |
+ |
VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); |
| 129 |
|
multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ |
| 130 |
|
copycolor(ap->v, ar.rcol); |
| 129 |
– |
if (ar.rt > 20.0*maxarad) /* limit vertex distance */ |
| 130 |
– |
VSUM(ap->p, ar.rorg, ar.rdir, 20.0*maxarad); |
| 131 |
– |
else |
| 132 |
– |
VCOPY(ap->p, ar.rop); |
| 131 |
|
return(ap); |
| 132 |
+ |
badsample: |
| 133 |
+ |
setcolor(ap->v, 0., 0., 0.); |
| 134 |
+ |
VCOPY(ap->p, hp->rp->rop); |
| 135 |
+ |
return(NULL); |
| 136 |
|
} |
| 137 |
|
|
| 138 |
|
|
| 139 |
|
/* Compute vectors and coefficients for Hessian/gradient calcs */ |
| 140 |
|
static void |
| 141 |
< |
comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop) |
| 141 |
> |
comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop) |
| 142 |
|
{ |
| 143 |
< |
FVECT v1; |
| 144 |
< |
double dot_e, dot_er, dot_r, dot_r1; |
| 143 |
> |
FVECT vcp; |
| 144 |
> |
double dot_e, dot_er, dot_r, dot_r1, J2; |
| 145 |
> |
int i; |
| 146 |
|
|
| 147 |
|
VSUB(ftp->r_i, ap0, rop); |
| 148 |
|
VSUB(ftp->r_i1, ap1, rop); |
| 149 |
|
VSUB(ftp->e_i, ap1, ap0); |
| 150 |
< |
VCROSS(v1, ftp->e_i, ftp->r_i); |
| 151 |
< |
ftp->nf = 1.0/DOT(v1,v1); |
| 150 |
> |
VCROSS(vcp, ftp->e_i, ftp->r_i); |
| 151 |
> |
ftp->nf = 1.0/DOT(vcp,vcp); |
| 152 |
|
dot_e = DOT(ftp->e_i,ftp->e_i); |
| 153 |
|
dot_er = DOT(ftp->e_i, ftp->r_i); |
| 154 |
|
dot_r = DOT(ftp->r_i,ftp->r_i); |
| 157 |
|
sqrt( ftp->nf ); |
| 158 |
|
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r + |
| 159 |
|
dot_e*ftp->I1 )*0.5*ftp->nf; |
| 160 |
< |
ftp->J2 = 0.5/dot_e*( 1.0/dot_r - 1.0/dot_r1 ) - |
| 161 |
< |
dot_er/dot_e*ftp->I2; |
| 160 |
> |
J2 = 0.5/dot_e*( 1.0/dot_r - 1.0/dot_r1 ) - dot_er/dot_e*ftp->I2; |
| 161 |
> |
for (i = 3; i--; ) |
| 162 |
> |
ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; |
| 163 |
|
} |
| 164 |
|
|
| 165 |
|
|
| 180 |
|
static void |
| 181 |
|
comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) |
| 182 |
|
{ |
| 183 |
< |
FVECT v1, v2; |
| 183 |
> |
FVECT vcp; |
| 184 |
|
FVECT m1[3], m2[3], m3[3], m4[3]; |
| 185 |
|
double d1, d2, d3, d4; |
| 186 |
|
double I3, J3, K3; |
| 195 |
|
J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3; |
| 196 |
|
K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3); |
| 197 |
|
/* intermediate matrices */ |
| 198 |
< |
VCROSS(v1, nrm, ftp->e_i); |
| 199 |
< |
for (j = 3; j--; ) |
| 196 |
< |
v2[j] = ftp->I2*ftp->r_i[j] + ftp->J2*ftp->e_i[j]; |
| 197 |
< |
compose_matrix(m1, v1, v2); |
| 198 |
> |
VCROSS(vcp, nrm, ftp->e_i); |
| 199 |
> |
compose_matrix(m1, vcp, ftp->rI2_eJ2); |
| 200 |
|
compose_matrix(m2, ftp->r_i, ftp->r_i); |
| 201 |
|
compose_matrix(m3, ftp->e_i, ftp->e_i); |
| 202 |
|
compose_matrix(m4, ftp->r_i, ftp->e_i); |
| 203 |
< |
VCROSS(v1, ftp->r_i, ftp->e_i); |
| 204 |
< |
d1 = DOT(nrm, v1); |
| 203 |
> |
VCROSS(vcp, ftp->r_i, ftp->e_i); |
| 204 |
> |
d1 = DOT(nrm, vcp); |
| 205 |
|
d2 = -d1*ftp->I2; |
| 206 |
|
d1 *= 2.0; |
| 207 |
|
for (i = 3; i--; ) /* final matrix sum */ |
| 253 |
|
f1 = 2.0*DOT(nrm, vcp); |
| 254 |
|
VCROSS(vcp, nrm, ftp->e_i); |
| 255 |
|
for (i = 3; i--; ) |
| 256 |
< |
grad[i] = (0.5/PI)*( ftp->I1*vcp[i] + |
| 255 |
< |
f1*(ftp->I2*ftp->r_i[i] + ftp->J2*ftp->e_i[i]) ); |
| 256 |
> |
grad[i] = (0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] ); |
| 257 |
|
} |
| 258 |
|
|
| 259 |
|
|