| 158 |
|
} |
| 159 |
|
|
| 160 |
|
|
| 161 |
< |
int |
| 162 |
< |
weightf(hp, x0, x1, x2) /* voxel weighting function */ |
| 163 |
< |
register HOLO *hp; |
| 164 |
< |
register int x0, x1, x2; |
| 165 |
< |
{ |
| 166 |
< |
switch (vlet(OCCUPANCY)) { |
| 167 |
< |
case 'U': /* uniform weighting */ |
| 168 |
< |
return(1); |
| 169 |
< |
case 'C': /* center weighting (crude) */ |
| 170 |
< |
x0 += x0 - hp->grid[0] + 1; |
| 171 |
< |
x0 = abs(x0)*hp->grid[1]*hp->grid[2]; |
| 172 |
< |
x1 += x1 - hp->grid[1] + 1; |
| 173 |
< |
x1 = abs(x1)*hp->grid[0]*hp->grid[2]; |
| 174 |
< |
x2 += x2 - hp->grid[2] + 1; |
| 175 |
< |
x2 = abs(x2)*hp->grid[0]*hp->grid[1]; |
| 176 |
< |
return(hp->grid[0]*hp->grid[1]*hp->grid[2] - |
| 177 |
< |
(x0+x1+x2)/3); |
| 178 |
< |
default: |
| 179 |
< |
badvalue(OCCUPANCY); |
| 180 |
< |
} |
| 181 |
< |
} |
| 182 |
< |
|
| 183 |
< |
|
| 184 |
< |
/* The following is by Daniel Cohen, taken from Graphics Gems IV, p. 368 */ |
| 185 |
< |
long |
| 186 |
< |
lineweight(hp, x, y, z, dx, dy, dz) /* compute weights along a line */ |
| 161 |
> |
double |
| 162 |
> |
beamvolume(hp, bi) /* compute approximate volume of a beam */ |
| 163 |
|
HOLO *hp; |
| 164 |
< |
int x, y, z, dx, dy, dz; |
| 164 |
> |
int bi; |
| 165 |
|
{ |
| 166 |
< |
long wres = 0; |
| 167 |
< |
int n, sx, sy, sz, exy, exz, ezy, ax, ay, az, bx, by, bz; |
| 168 |
< |
|
| 169 |
< |
sx = sgn(dx); sy = sgn(dy); sz = sgn(dz); |
| 170 |
< |
ax = abs(dx); ay = abs(dy); az = abs(dz); |
| 171 |
< |
bx = 2*ax; by = 2*ay; bz = 2*az; |
| 172 |
< |
exy = ay-ax; exz = az-ax; ezy = ay-az; |
| 173 |
< |
n = ax+ay+az + 1; /* added increment to visit last */ |
| 174 |
< |
while (n--) { |
| 175 |
< |
wres += weightf(hp, x, y, z); |
| 176 |
< |
if (exy < 0) { |
| 177 |
< |
if (exz < 0) { |
| 178 |
< |
x += sx; |
| 179 |
< |
exy += by; exz += bz; |
| 180 |
< |
} else { |
| 181 |
< |
z += sz; |
| 182 |
< |
exz -= bx; ezy += by; |
| 183 |
< |
} |
| 184 |
< |
} else { |
| 185 |
< |
if (ezy < 0) { |
| 186 |
< |
z += sz; |
| 211 |
< |
exz -= bx; ezy += by; |
| 212 |
< |
} else { |
| 213 |
< |
y += sy; |
| 214 |
< |
exy -= bx; ezy -= bz; |
| 215 |
< |
} |
| 216 |
< |
} |
| 166 |
> |
GCOORD gc[2]; |
| 167 |
> |
FVECT cp[4], edgeA, edgeB, cent[2]; |
| 168 |
> |
FVECT v, crossp[2], diffv; |
| 169 |
> |
double vol[2]; |
| 170 |
> |
register int i; |
| 171 |
> |
/* get grid coordinates */ |
| 172 |
> |
if (!hdbcoord(gc, hp, bi)) |
| 173 |
> |
error(CONSISTENCY, "bad beam index in beamvolume"); |
| 174 |
> |
for (i = 0; i < 2; i++) { /* compute cell area vectors */ |
| 175 |
> |
hdcell(cp, hp, gc+i); |
| 176 |
> |
VSUM(edgeA, cp[1], cp[0], -1.0); |
| 177 |
> |
VSUM(edgeB, cp[3], cp[1], -1.0); |
| 178 |
> |
fcross(crossp[i], edgeA, edgeB); |
| 179 |
> |
VSUM(edgeA, cp[2], cp[3], -1.0); |
| 180 |
> |
VSUM(edgeB, cp[0], cp[2], -1.0); |
| 181 |
> |
fcross(v, edgeA, edgeB); |
| 182 |
> |
VSUM(crossp[i], crossp[i], v, 1.0); |
| 183 |
> |
/* compute center */ |
| 184 |
> |
cent[i][0] = 0.5*(cp[0][0] + cp[2][0]); |
| 185 |
> |
cent[i][1] = 0.5*(cp[0][1] + cp[2][1]); |
| 186 |
> |
cent[i][2] = 0.5*(cp[0][2] + cp[2][2]); |
| 187 |
|
} |
| 188 |
< |
return(wres); |
| 188 |
> |
/* compute difference vector */ |
| 189 |
> |
VSUM(diffv, cent[1], cent[0], -1.0); |
| 190 |
> |
for (i = 0; i < 2; i++) { /* compute volume contributions */ |
| 191 |
> |
vol[i] = 0.25*DOT(crossp[i], diffv); |
| 192 |
> |
if (vol[i] < 0.) vol[i] = -vol[i]; |
| 193 |
> |
} |
| 194 |
> |
return(vol[0] + vol[1]); /* return total volume */ |
| 195 |
|
} |
| 196 |
|
|
| 197 |
|
|
| 198 |
|
init_global() /* initialize global ray computation */ |
| 199 |
|
{ |
| 200 |
|
long wtotal = 0; |
| 225 |
– |
int i, j; |
| 226 |
– |
int lseg[2][3]; |
| 201 |
|
double frac; |
| 202 |
< |
register int k; |
| 202 |
> |
int i; |
| 203 |
> |
register int j, k; |
| 204 |
|
/* free old list */ |
| 205 |
|
if (complen > 0) |
| 206 |
|
free((char *)complist); |
| 213 |
|
error(SYSTEM, "out of memory in init_global"); |
| 214 |
|
/* compute beam weights */ |
| 215 |
|
k = 0; |
| 216 |
< |
for (j = 0; hdlist[j] != NULL; j++) |
| 216 |
> |
for (j = 0; hdlist[j] != NULL; j++) { |
| 217 |
> |
frac = 512. * hdlist[j]->wg[0] * |
| 218 |
> |
hdlist[j]->wg[1] * hdlist[j]->wg[2]; |
| 219 |
|
for (i = nbeams(hdlist[j]); i > 0; i--) { |
| 243 |
– |
hdlseg(lseg, hdlist[j], i); |
| 220 |
|
complist[k].hd = j; |
| 221 |
|
complist[k].bi = i; |
| 222 |
< |
complist[k].nr = lineweight( hdlist[j], |
| 247 |
< |
lseg[0][0], lseg[0][1], lseg[0][2], |
| 248 |
< |
lseg[1][0] - lseg[0][0], |
| 249 |
< |
lseg[1][1] - lseg[0][1], |
| 250 |
< |
lseg[1][2] - lseg[0][2] ); |
| 222 |
> |
complist[k].nr = frac*beamvolume(hdlist[j], i) + 0.5; |
| 223 |
|
wtotal += complist[k++].nr; |
| 224 |
|
} |
| 225 |
+ |
} |
| 226 |
|
/* adjust weights */ |
| 227 |
|
if (vdef(DISKSPACE)) |
| 228 |
|
frac = 1024.*1024.*vflt(DISKSPACE) / (wtotal*sizeof(RAYVAL)); |
