1 ## Copyright (C) 2010 Alex Opie <lx_op@orcon.net.nz>
3 ## This program is free software; you can redistribute it and/or modify it
4 ## under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 3 of the License, or (at
6 ## your option) any later version.
8 ## This program is distributed in the hope that it will be useful, but
9 ## WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 ## General Public License for more details.
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; see the file COPYING. If not, see
15 ## <http://www.gnu.org/licenses/>.
19 ## @defun {@var{P} =} phantom ('Shepp-Logan', @var{n})
21 ## Produces the Shepp-Logan phantom, with size @var{n} x @var{n}.
22 ## If @var{n} is omitted, 256 is used.
24 ## @defunx {@var{P} =} phantom ('Modified Shepp-Logan', @var{n})
26 ## Produces a modified version of the Shepp-Logan phantom which has
27 ## higher contrast than the original, with size @var{n} x @var{n}.
28 ## If @var{n} is omitted, 256 is used.
30 ## @defunx {@var{P} =} phantom (@var{ellipses}, @var{n})
32 ## Produces a custom phantom using the ellipses described in @var{ellipses}.
33 ## Each row of @var{ellipses} describes one ellipse, and must have 6 columns:
34 ## @{I, a, b, x0, y0, phi@}:
37 ## is the additive intensity of the ellipse
40 ## is the length of the major axis
43 ## is the length of the minor axis
46 ## is the horizontal offset of the centre of the ellipse
49 ## is the vercal offset of the centre of the ellipse
52 ## is the counterclockwise rotation of the ellipse in degrees,
53 ## measured as the angle between the x axis and the ellipse major axis.
57 ## The image bounding box in the algorithm is @{[-1, -1], [1, 1]@}, so the
58 ## values of a, b, x0, y0 should all be specified with this in mind.
59 ## If @var{n} is omitted, 256 is used.
61 ## @defunx {@var{P} =} phantom (@var{n})
63 ## Creates a modified Shepp-Logan phantom with size @var{n} x @var{n}.
65 ## @defunx {@var{P} = } phantom ()
67 ## Creates a modified Shepp-Logan phantom with size 256 x 256.
70 ## Create a Shepp-Logan or modified Shepp-Logan phantom.
72 ## A phantom is a known object (either real or purely mathematical) that
73 ## is used for testing image reconstruction algorithms. The Shepp-Logan
74 ## phantom is a popular mathematical model of a cranial slice, made up
75 ## of a set of ellipses. This allows rigorous testing of computed
76 ## tomography (CT) algorithms as it can be analytically transformed with
77 ## the radon transform (see the function @command{radon}).
88 ## Shepp, L. A.; Logan, B. F.; Reconstructing Interior Head Tissue
89 ## from X-Ray Transmissions, IEEE Transactions on Nuclear Science,
92 ## Toft, P.; "The Radon Transform - Theory and Implementation", Ph.D. thesis,
93 ## Department of Mathematical Modelling, Technical University
94 ## of Denmark, June 1996.
96 function p = phantom (varargin)
98 [n, ellipses] = read_args (varargin {:});
103 # Create the pixel grid
104 xvals = (-1 : 2 / (n - 1) : 1);
105 xgrid = repmat (xvals, n, 1);
107 for i = 1:size (ellipses, 1)
109 a2 = ellipses (i, 2)^2;
110 b2 = ellipses (i, 3)^2;
111 x0 = ellipses (i, 4);
112 y0 = ellipses (i, 5);
113 phi = ellipses (i, 6) * pi / 180; # Rotation angle in radians
115 # Create the offset x and y values for the grid
117 y = rot90 (xgrid) - y0;
122 # Find the pixels within the ellipse
123 locs = find (((x .* cos_p + y .* sin_p).^2) ./ a2 ...
124 + ((y .* cos_p - x .* sin_p).^2) ./ b2 <= 1);
126 # Add the ellipse intensity to those pixels
127 p (locs) = p (locs) + I;
131 function [n, ellip] = read_args (varargin)
133 ellip = mod_shepp_logan ();
136 if (ischar (varargin {1}))
137 ellip = select_phantom (varargin {1});
138 elseif (numel (varargin {1}) == 1)
141 if (size (varargin {1}, 2) != 6)
142 error ("Wrong number of columns in user phantom");
144 ellip = varargin {1};
148 if (ischar (varargin {1}))
149 ellip = select_phantom (varargin {1});
151 if (size (varargin {1}, 2) != 6)
152 error ("Wrong number of columns in user phantom");
154 ellip = varargin {1};
157 warning ("Extra arguments passed to phantom were ignored");
161 function e = select_phantom (name)
162 if (strcmpi (name, 'shepp-logan'))
164 elseif (strcmpi (name, 'modified shepp-logan'))
165 e = mod_shepp_logan ();
167 error ("Unknown phantom type: %s", name);
171 function e = shepp_logan
172 # Standard head phantom, taken from Shepp & Logan
173 e = [ 2, .69, .92, 0, 0, 0;
174 -.98, .6624, .8740, 0, -.0184, 0;
175 -.02, .1100, .3100, .22, 0, -18;
176 -.02, .1600, .4100, -.22, 0, 18;
177 .01, .2100, .2500, 0, .35, 0;
178 .01, .0460, .0460, 0, .1, 0;
179 .02, .0460, .0460, 0, -.1, 0;
180 .01, .0460, .0230, -.08, -.605, 0;
181 .01, .0230, .0230, 0, -.606, 0;
182 .01, .0230, .0460, .06, -.605, 0];
185 function e = mod_shepp_logan
186 # Modified version of Shepp & Logan's head phantom,
187 # adjusted to improve contrast. Taken from Toft.
188 e = [ 1, .69, .92, 0, 0, 0;
189 -.80, .6624, .8740, 0, -.0184, 0;
190 -.20, .1100, .3100, .22, 0, -18;
191 -.20, .1600, .4100, -.22, 0, 18;
192 .10, .2100, .2500, 0, .35, 0;
193 .10, .0460, .0460, 0, .1, 0;
194 .10, .0460, .0460, 0, -.1, 0;
195 .10, .0460, .0230, -.08, -.605, 0;
196 .10, .0230, .0230, 0, -.606, 0;
197 .10, .0230, .0460, .06, -.605, 0];
201 # # Add any further phantoms of interest here
202 # e = [ 0, 0, 0, 0, 0, 0;
207 %! P = phantom (512);
212 % Fixed typo that prevented the use of custom ellipses
215 % First commit to repo