1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
include <imhdr.h>
include "pvol.h"
$for (silrdx)
# VTRANSMIT -- Compute the intensities of each output image pixel in the
# current line as a function of its existing intensity plus the emission
# and absorption from each contributing voxel.
procedure vtransmit$t (inbuf,nx,ny,nz,nh, px1,px2, iline,iband,nvox, oline, vp)
PIXEL inbuf[nx,ny,nz,nh] # Input data buffer for current set of yz slices
int nx,ny,nz,nh # Dimensions of current input buffer
int px1,px2 # Range of columns in current yz slice set
int iline[nvox] # Input image lines for current projection ray
int iband[nvox] # Input image bands for current projection ray
int nvox # Number of voxels in current projection column
real oline[ARB] # output image line buffer
pointer vp # Volume projection descriptor
bool use_both
int i, vox, opelem, intelem, frontvox, backvox
real amin, amax, vox_op, vox_int, ival, attenuate, vimin, ifac, ofac, distwt
begin
# Dereference most frequently used structure elements.
amin = AMIN(vp)
amax = AMAX(vp)
vimin = VIMIN(vp)
intelem = 1
opelem = OPACELEM(vp)
if (nh > 1) {
use_both = true
if (opelem == 1)
intelem = 2
else if (IS_INDEFI(opelem))
opelem = 2
} else {
use_both = false
opelem = 1
}
# Set up for opacity, intensity, or both.
ifac = (IIMAX(vp) - IIMIN(vp)) / (VIMAX(vp) - vimin)
if (PTYPE(vp) == P_ATTENUATE || use_both)
ofac = (amax - amin) / (OMAX(vp) - OMIN(vp))
# Since we are in memory anyway, it is more convenient to traverse
# the columns in the outer loop and the voxels from different bands
# and lines in the inner loop. This is necessary when distance
# weighting and the distance cutoff option is on (we need to know
# the range of usable voxels in a given column before projecting).
if (PTYPE(vp) == P_INVDISPOW || PTYPE(vp) == P_MODN)
do i = px1, px2 {
if (DISCUTOFF(vp) == NO) {
frontvox = nvox
backvox = 1
} else {
frontvox = 1
backvox = nvox
do vox = 1, nvox {
vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],intelem]
if (vimin <= vox_int && vox_int < VIMAX(vp)) {
frontvox = max (frontvox, vox)
backvox = min (backvox, vox)
}
}
}
if (frontvox - backvox < 0)
next
do vox = backvox, frontvox {
distwt = (real(vox-backvox+1) /
real(frontvox-backvox+1)) ** DISPOWER(vp)
# Opacity transformation function.
if (use_both) {
vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
opelem] * OSCALE(vp)
if (vox_op < OMIN(vp))
attenuate = amax
else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
attenuate = amax - (vox_op - OMIN(vp)) * ofac
else
attenuate = amin
oline[i] = oline[i] * attenuate
}
# Intensity transformation function.
vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
intelem]
if (vox_int < vimin)
ival = IIMIN(vp)
else if (vimin <= vox_int && vox_int < VIMAX(vp))
ival = IIMIN(vp) + (vox_int - vimin) * ifac
else
ival = IIMAX(vp)
if (PTYPE(vp) == P_INVDISPOW)
oline[i] = oline[i] + ival * distwt
else if (mod (int (ival/(IIMAX(vp)-IIMIN(vp)) * 100.0),
MODN(vp)) == 0)
oline[i] = oline[i] + ival * distwt
}
}
else
do i = px1, px2
do vox = 1, nvox {
# Opacity transformation function.
if (PTYPE(vp) == P_ATTENUATE || use_both) {
vox_op = inbuf[(i-px1+1), iline[vox], iband[vox],
opelem] * OSCALE(vp)
if (vox_op < OMIN(vp))
attenuate = amax
else if (OMIN(vp) <= vox_op && vox_op < OMAX(vp))
attenuate = amax - (vox_op - OMIN(vp)) * ofac
else
attenuate = amin
oline[i] = oline[i] * attenuate
}
# Intensity transformation function.
if (PTYPE(vp) != P_ATTENUATE || use_both) {
vox_int = inbuf[(i-px1+1),iline[vox],iband[vox],
intelem]
if (vox_int < vimin)
ival = IIMIN(vp)
else if (vimin <= vox_int && vox_int < VIMAX(vp))
ival = IIMIN(vp) + (vox_int - vimin) * ifac
else
ival = IIMAX(vp)
if (PTYPE(vp) == P_AVERAGE)
oline[i] = oline[i] + ival * 1.0 / real(nvox)
else if (PTYPE(vp) == P_SUM)
oline[i] = oline[i] + ival
else if (PTYPE(vp) == P_LASTONLY)
if (ival > 0.0)
oline[i] = ival
}
}
end
$endfor
|
