| 141 |
|
comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop) |
| 142 |
|
{ |
| 143 |
|
FVECT vcp; |
| 144 |
< |
double dot_e, dot_er, dot_r, dot_r1, J2; |
| 144 |
> |
double dot_e, dot_er, rdot_r, rdot_r1, J2; |
| 145 |
|
int i; |
| 146 |
|
|
| 147 |
|
VSUB(ftp->r_i, ap0, rop); |
| 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); |
| 155 |
< |
dot_r1 = DOT(ftp->r_i1,ftp->r_i1); |
| 156 |
< |
ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) / sqrt(dot_r*dot_r1) ) * |
| 154 |
> |
rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i); |
| 155 |
> |
rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1); |
| 156 |
> |
ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) * |
| 157 |
|
sqrt( ftp->nf ); |
| 158 |
< |
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r + |
| 158 |
> |
ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r + |
| 159 |
|
dot_e*ftp->I1 )*0.5*ftp->nf; |
| 160 |
< |
J2 = 0.5/dot_e*( 1.0/dot_r - 1.0/dot_r1 ) - dot_er/dot_e*ftp->I2; |
| 160 |
> |
J2 = ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e; |
| 161 |
|
for (i = 3; i--; ) |
| 162 |
|
ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; |
| 163 |
|
} |
| 209 |
|
hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] + |
| 210 |
|
2.0*J3*m4[i][j] ); |
| 211 |
|
hess[i][j] += d2*(i==j); |
| 212 |
< |
hess[i][j] *= -1.0/PI; |
| 212 |
> |
hess[i][j] *= 1.0/PI; |
| 213 |
|
} |
| 214 |
|
} |
| 215 |
|
|
| 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] + f1*ftp->rI2_eJ2[i] ); |
| 256 |
> |
grad[i] = (-0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] ); |
| 257 |
|
} |
| 258 |
|
|
| 259 |
|
|
| 455 |
|
|
| 456 |
|
if (ra != NULL) /* extract eigenvectors & radii */ |
| 457 |
|
eigenvectors(uv, ra, hessian); |
| 458 |
< |
if (pg != NULL) { /* tangential position gradient/PI */ |
| 459 |
< |
pg[0] = DOT(gradient, uv[0]) / PI; |
| 460 |
< |
pg[1] = DOT(gradient, uv[1]) / PI; |
| 458 |
> |
if (pg != NULL) { /* tangential position gradient */ |
| 459 |
> |
pg[0] = DOT(gradient, uv[0]); |
| 460 |
> |
pg[1] = DOT(gradient, uv[1]); |
| 461 |
|
} |
| 462 |
|
} |
| 463 |
|
|
| 534 |
|
free(hp); |
| 535 |
|
return(-1); /* no radius or gradient calc. */ |
| 536 |
|
} |
| 537 |
< |
multcolor(acol, hp->acoef); /* normalize Y values */ |
| 538 |
< |
if ((d = bright(acol)) > FTINY) |
| 539 |
< |
d = 1.0/d; |
| 537 |
> |
if (bright(acol) > FTINY) /* normalize Y values */ |
| 538 |
> |
d = cnt/bright(acol); |
| 539 |
|
else |
| 540 |
|
d = 0.0; |
| 541 |
|
ap = hp->sa; /* relative Y channel from here on... */ |
| 542 |
|
for (i = hp->ns*hp->ns; i--; ap++) |
| 543 |
< |
colval(ap->v,CIEY) = bright(ap->v)*d + 0.0314; |
| 543 |
> |
colval(ap->v,CIEY) = bright(ap->v)*d + 0.01; |
| 544 |
|
|
| 545 |
|
if (uv == NULL) /* make sure we have axis pointers */ |
| 546 |
|
uv = my_uv; |
