octave_packages/image-1.0.15/phantom.m

   1 ## Copyright (C) 2010  Alex Opie  <lx_op@orcon.net.nz>
   2 ##
   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.
   7 ##
   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.
  12 ##
  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/>.
  16
  17 ## -*- texinfo -*-
  18 ##
  19 ## @defun {@var{P} =} phantom ('Shepp-Logan', @var{n})
  20 ##
  21 ## Produces the Shepp-Logan phantom, with size @var{n} x @var{n}.
  22 ## If @var{n} is omitted, 256 is used.
  23 ##
  24 ## @defunx {@var{P} =} phantom ('Modified Shepp-Logan', @var{n})
  25 ##
  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.
  29 ##
  30 ## @defunx {@var{P} =} phantom (@var{ellipses}, @var{n})
  31 ##
  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@}:
  35 ##  @table @abbr
  36 ##  @item I
  37 ##    is the additive intensity of the ellipse
  38 ##
  39 ##  @item a
  40 ##    is the length of the major axis
  41 ##
  42 ##  @item b
  43 ##    is the length of the minor axis
  44 ##
  45 ##  @item x0
  46 ##    is the horizontal offset of the centre of the ellipse
  47 ##
  48 ##  @item y0
  49 ##    is the vercal offset of the centre of the ellipse
  50 ##
  51 ##  @item phi
  52 ##    is the counterclockwise rotation of the ellipse in degrees,
  53 ##    measured as the angle between the x axis and the ellipse major axis.
  54 ##
  55 ## @end table
  56 ##
  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.
  60 ##
  61 ## @defunx {@var{P} =} phantom (@var{n})
  62 ##
  63 ## Creates a modified Shepp-Logan phantom with size @var{n} x @var{n}.
  64 ##
  65 ## @defunx {@var{P} = } phantom ()
  66 ##
  67 ## Creates a modified Shepp-Logan phantom with size 256 x 256.
  68 ## @end defun
  69 ##
  70 ## Create a Shepp-Logan or modified Shepp-Logan phantom.
  71 ##
  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}).
  78 ##
  79 ## Example:
  80 ##
  81 ## @example
  82 ##   P = phantom (512);
  83 ##   imshow (P, []);
  84 ## @end example
  85 ##
  86 ## References:
  87 ##
  88 ##  Shepp, L. A.; Logan, B. F.; Reconstructing Interior Head Tissue
  89 ##  from X-Ray Transmissions, IEEE Transactions on Nuclear Science,
  90 ##  Feb. 1974, p. 232.
  91 ##
  92 ##  Toft, P.; "The Radon Transform - Theory and Implementation", Ph.D. thesis,
  93 ##  Department of Mathematical Modelling, Technical University
  94 ##  of Denmark, June 1996.
  95
  96 function p = phantom (varargin)
  97
  98   [n, ellipses] = read_args (varargin {:});
  99
 100   # Blank image
 101   p = zeros (n);
 102
 103   # Create the pixel grid
 104   xvals = (-1 : 2 / (n - 1) : 1);
 105   xgrid = repmat (xvals, n, 1);
 106
 107   for i = 1:size (ellipses, 1)
 108     I   = ellipses (i, 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
 114
 115     # Create the offset x and y values for the grid
 116     x = xgrid - x0;
 117     y = rot90 (xgrid) - y0;
 118
 119     cos_p = cos (phi);
 120     sin_p = sin (phi);
 121
 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);
 125
 126     # Add the ellipse intensity to those pixels
 127     p (locs) = p (locs) + I;
 128   endfor
 129 endfunction
 130
 131 function [n, ellip] = read_args (varargin)
 132   n = 256;
 133   ellip = mod_shepp_logan ();
 134
 135   if (nargin == 1)
 136     if (ischar (varargin {1}))
 137       ellip = select_phantom (varargin {1});
 138     elseif (numel (varargin {1}) == 1)
 139       n = varargin {1};
 140     else
 141       if (size (varargin {1}, 2) != 6)
 142         error ("Wrong number of columns in user phantom");
 143       endif
 144       ellip = varargin {1};
 145     endif
 146   elseif (nargin == 2)
 147     n = varargin {2};
 148     if (ischar (varargin {1}))
 149       ellip = select_phantom (varargin {1});
 150     else
 151       if (size (varargin {1}, 2) != 6)
 152         error ("Wrong number of columns in user phantom");
 153       endif
 154       ellip = varargin {1};
 155     endif
 156   elseif (nargin > 2)
 157     warning ("Extra arguments passed to phantom were ignored");
 158   endif
 159 endfunction
 160
 161 function e = select_phantom (name)
 162   if (strcmpi (name, 'shepp-logan'))
 163     e = shepp_logan ();
 164   elseif (strcmpi (name, 'modified shepp-logan'))
 165     e = mod_shepp_logan ();
 166   else
 167     error ("Unknown phantom type: %s", name);
 168   endif
 169 endfunction
 170
 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];
 183 endfunction
 184
 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];
 198 endfunction
 199
 200 #function e = ??
 201 #  # Add any further phantoms of interest here
 202 #  e = [ 0, 0, 0, 0, 0, 0;
 203 #        0, 0, 0, 0, 0, 0];
 204 #endfunction
 205
 206 %!demo
 207 %! P = phantom (512);
 208 %! imshow (P, []);
 209
 210 % $Log$
 211 % 2010/03/31 lxop
 212 % Fixed typo that prevented the use of custom ellipses
 213 %
 214 % 2010/03/09 lxop
 215 % First commit to repo