-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathphantom3d.m
224 lines (187 loc) · 8.1 KB
/
phantom3d.m
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
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
function [p,ellipse]=phantom3d(varargin)
%PHANTOM3D Three-dimensional analogue of MATLAB Shepp-Logan phantom
% P = PHANTOM3D(DEF,N) generates a 3D head phantom that can
% be used to test 3-D reconstruction algorithms.
%
% DEF is a string that specifies the type of head phantom to generate.
% Valid values are:
%
% 'Shepp-Logan' A test image used widely by researchers in
% tomography
% 'Modified Shepp-Logan' (default) A variant of the Shepp-Logan phantom
% in which the contrast is improved for better
% visual perception.
%
% N is a scalar that specifies the grid size of P.
% If you omit the argument, N defaults to 64.
%
% P = PHANTOM3D(E,N) generates a user-defined phantom, where each row
% of the matrix E specifies an ellipsoid in the image. E has ten columns,
% with each column containing a different parameter for the ellipsoids:
%
% Column 1: A the additive intensity value of the ellipsoid
% Column 2: a the length of the x semi-axis of the ellipsoid
% Column 3: b the length of the y semi-axis of the ellipsoid
% Column 4: c the length of the z semi-axis of the ellipsoid
% Column 5: x0 the x-coordinate of the center of the ellipsoid
% Column 6: y0 the y-coordinate of the center of the ellipsoid
% Column 7: z0 the z-coordinate of the center of the ellipsoid
% Column 8: phi phi Euler angle (in degrees) (rotation about z-axis)
% Column 9: theta theta Euler angle (in degrees) (rotation about x-axis)
% Column 10: psi psi Euler angle (in degrees) (rotation about z-axis)
%
% For purposes of generating the phantom, the domains for the x-, y-, and
% z-axes span [-1,1]. Columns 2 through 7 must be specified in terms
% of this range.
%
% [P,E] = PHANTOM3D(...) returns the matrix E used to generate the phantom.
%
% Class Support
% -------------
% All inputs must be of class double. All outputs are of class double.
%
% Remarks
% -------
% For any given voxel in the output image, the voxel's value is equal to the
% sum of the additive intensity values of all ellipsoids that the voxel is a
% part of. If a voxel is not part of any ellipsoid, its value is 0.
%
% The additive intensity value A for an ellipsoid can be positive or negative;
% if it is negative, the ellipsoid will be darker than the surrounding pixels.
% Note that, depending on the values of A, some voxels may have values outside
% the range [0,1].
%
% Example
% -------
% ph = phantom3d(128);
% figure, imshow(squeeze(ph(64,:,:)))
%
% Copyright 2005 Matthias Christian Schabel (matthias @ stanfordalumni . org)
% University of Utah Department of Radiology
% Utah Center for Advanced Imaging Research
% 729 Arapeen Drive
% Salt Lake City, UT 84108-1218
%
% This code is released under the Gnu Public License (GPL). For more information,
% see : http://www.gnu.org/copyleft/gpl.html
%
% Portions of this code are based on phantom.m, copyrighted by the Mathworks
%
[ellipse,n] = parse_inputs(varargin{:});
p = zeros([n n n]);
rng = ( (0:n-1)-(n-1)/2 ) / ((n-1)/2);
[x,y,z] = meshgrid(rng,rng,rng);
coord = [flatten(x); flatten(y); flatten(z)];
p = flatten(p);
for k = 1:size(ellipse,1)
A = ellipse(k,1); % Amplitude change for this ellipsoid
asq = ellipse(k,2)^2; % a^2
bsq = ellipse(k,3)^2; % b^2
csq = ellipse(k,4)^2; % c^2
x0 = ellipse(k,5); % x offset
y0 = ellipse(k,6); % y offset
z0 = ellipse(k,7); % z offset
phi = ellipse(k,8)*pi/180; % first Euler angle in radians
theta = ellipse(k,9)*pi/180; % second Euler angle in radians
psi = ellipse(k,10)*pi/180; % third Euler angle in radians
cphi = cos(phi);
sphi = sin(phi);
ctheta = cos(theta);
stheta = sin(theta);
cpsi = cos(psi);
spsi = sin(psi);
% Euler rotation matrix
alpha = [cpsi*cphi-ctheta*sphi*spsi cpsi*sphi+ctheta*cphi*spsi spsi*stheta;
-spsi*cphi-ctheta*sphi*cpsi -spsi*sphi+ctheta*cphi*cpsi cpsi*stheta;
stheta*sphi -stheta*cphi ctheta];
% rotated ellipsoid coordinates
coordp = alpha*coord;
idx = find((coordp(1,:)-x0).^2./asq + (coordp(2,:)-y0).^2./bsq + (coordp(3,:)-z0).^2./csq <= 1);
p(idx) = p(idx) + A;
end
p = reshape(p,[n n n]);
return;
function out = flatten(in)
out = reshape(in,[1 prod(size(in))]);
return;
function [e,n] = parse_inputs(varargin)
% e is the m-by-10 array which defines ellipsoids
% n is the size of the phantom brain image
n = 128; % The default size
e = [];
defaults = {'shepp-logan', 'modified shepp-logan', 'yu-ye-wang'};
for i=1:nargin
if ischar(varargin{i}) % Look for a default phantom
def = lower(varargin{i});
idx = strmatch(def, defaults);
if isempty(idx)
eid = sprintf('Images:%s:unknownPhantom',mfilename);
msg = 'Unknown default phantom selected.';
error(eid,'%s',msg);
end
switch defaults{idx}
case 'shepp-logan'
e = shepp_logan;
case 'modified shepp-logan'
e = modified_shepp_logan;
case 'yu-ye-wang'
e = yu_ye_wang;
end
elseif numel(varargin{i})==1
n = varargin{i}; % a scalar is the image size
elseif ndims(varargin{i})==2 && size(varargin{i},2)==10
e = varargin{i}; % user specified phantom
else
eid = sprintf('Images:%s:invalidInputArgs',mfilename);
msg = 'Invalid input arguments.';
error(eid,'%s',msg);
end
end
% ellipse is not yet defined
if isempty(e)
e = modified_shepp_logan;
end
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Default head phantoms: %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function e = shepp_logan
e = modified_shepp_logan;
e(:,1) = [1 -.98 -.02 -.02 .01 .01 .01 .01 .01 .01];
return;
function e = modified_shepp_logan
%
% This head phantom is the same as the Shepp-Logan except
% the intensities are changed to yield higher contrast in
% the image. Taken from Toft, 199-200.
%
% A a b c x0 y0 z0 phi theta psi
% -----------------------------------------------------------------
e = [ 1 .6900 .920 .810 0 0 0 0 0 0
-.8 .6624 .874 .780 0 -.0184 0 0 0 0
-.2 .1100 .310 .220 .22 0 0 -18 0 10
-.2 .1600 .410 .280 -.22 0 0 18 0 10
.1 .2100 .250 .410 0 .35 -.15 0 0 0
.1 .0460 .046 .050 0 .1 .25 0 0 0
.1 .0460 .046 .050 0 -.1 .25 0 0 0
.1 .0460 .023 .050 -.08 -.605 0 0 0 0
.1 .0230 .023 .020 0 -.606 0 0 0 0
.1 .0230 .046 .020 .06 -.605 0 0 0 0 ];
return;
function e = yu_ye_wang
%
% Yu H, Ye Y, Wang G, Katsevich-Type Algorithms for Variable Radius Spiral Cone-Beam CT
%
% A a b c x0 y0 z0 phi theta psi
% -----------------------------------------------------------------
e = [ 1 .6900 .920 .900 0 0 0 0 0 0
-.8 .6624 .874 .880 0 0 0 0 0 0
-.2 .4100 .160 .210 -.22 0 -.25 108 0 0
-.2 .3100 .110 .220 .22 0 -.25 72 0 0
.2 .2100 .250 .500 0 .35 -.25 0 0 0
.2 .0460 .046 .046 0 .1 -.25 0 0 0
.1 .0460 .023 .020 -.08 -.65 -.25 0 0 0
.1 .0460 .023 .020 .06 -.65 -.25 90 0 0
.2 .0560 .040 .100 .06 -.105 .625 90 0 0
-.2 .0560 .056 .100 0 .100 .625 0 0 0 ];
return;