X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Fimage-1.0.15%2Fphantom.m;fp=octave_packages%2Fimage-1.0.15%2Fphantom.m;h=efd51a4675780059966e0d772aceb2adddce7139;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05 diff --git a/octave_packages/image-1.0.15/phantom.m b/octave_packages/image-1.0.15/phantom.m new file mode 100644 index 0000000..efd51a4 --- /dev/null +++ b/octave_packages/image-1.0.15/phantom.m @@ -0,0 +1,215 @@ +## Copyright (C) 2010 Alex Opie +## +## This program is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## This program is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## +## @defun {@var{P} =} phantom ('Shepp-Logan', @var{n}) +## +## Produces the Shepp-Logan phantom, with size @var{n} x @var{n}. +## If @var{n} is omitted, 256 is used. +## +## @defunx {@var{P} =} phantom ('Modified Shepp-Logan', @var{n}) +## +## Produces a modified version of the Shepp-Logan phantom which has +## higher contrast than the original, with size @var{n} x @var{n}. +## If @var{n} is omitted, 256 is used. +## +## @defunx {@var{P} =} phantom (@var{ellipses}, @var{n}) +## +## Produces a custom phantom using the ellipses described in @var{ellipses}. +## Each row of @var{ellipses} describes one ellipse, and must have 6 columns: +## @{I, a, b, x0, y0, phi@}: +## @table @abbr +## @item I +## is the additive intensity of the ellipse +## +## @item a +## is the length of the major axis +## +## @item b +## is the length of the minor axis +## +## @item x0 +## is the horizontal offset of the centre of the ellipse +## +## @item y0 +## is the vercal offset of the centre of the ellipse +## +## @item phi +## is the counterclockwise rotation of the ellipse in degrees, +## measured as the angle between the x axis and the ellipse major axis. +## +## @end table +## +## The image bounding box in the algorithm is @{[-1, -1], [1, 1]@}, so the +## values of a, b, x0, y0 should all be specified with this in mind. +## If @var{n} is omitted, 256 is used. +## +## @defunx {@var{P} =} phantom (@var{n}) +## +## Creates a modified Shepp-Logan phantom with size @var{n} x @var{n}. +## +## @defunx {@var{P} = } phantom () +## +## Creates a modified Shepp-Logan phantom with size 256 x 256. +## @end defun +## +## Create a Shepp-Logan or modified Shepp-Logan phantom. +## +## A phantom is a known object (either real or purely mathematical) that +## is used for testing image reconstruction algorithms. The Shepp-Logan +## phantom is a popular mathematical model of a cranial slice, made up +## of a set of ellipses. This allows rigorous testing of computed +## tomography (CT) algorithms as it can be analytically transformed with +## the radon transform (see the function @command{radon}). +## +## Example: +## +## @example +## P = phantom (512); +## imshow (P, []); +## @end example +## +## References: +## +## Shepp, L. A.; Logan, B. F.; Reconstructing Interior Head Tissue +## from X-Ray Transmissions, IEEE Transactions on Nuclear Science, +## Feb. 1974, p. 232. +## +## Toft, P.; "The Radon Transform - Theory and Implementation", Ph.D. thesis, +## Department of Mathematical Modelling, Technical University +## of Denmark, June 1996. + +function p = phantom (varargin) + + [n, ellipses] = read_args (varargin {:}); + + # Blank image + p = zeros (n); + + # Create the pixel grid + xvals = (-1 : 2 / (n - 1) : 1); + xgrid = repmat (xvals, n, 1); + + for i = 1:size (ellipses, 1) + I = ellipses (i, 1); + a2 = ellipses (i, 2)^2; + b2 = ellipses (i, 3)^2; + x0 = ellipses (i, 4); + y0 = ellipses (i, 5); + phi = ellipses (i, 6) * pi / 180; # Rotation angle in radians + + # Create the offset x and y values for the grid + x = xgrid - x0; + y = rot90 (xgrid) - y0; + + cos_p = cos (phi); + sin_p = sin (phi); + + # Find the pixels within the ellipse + locs = find (((x .* cos_p + y .* sin_p).^2) ./ a2 ... + + ((y .* cos_p - x .* sin_p).^2) ./ b2 <= 1); + + # Add the ellipse intensity to those pixels + p (locs) = p (locs) + I; + endfor +endfunction + +function [n, ellip] = read_args (varargin) + n = 256; + ellip = mod_shepp_logan (); + + if (nargin == 1) + if (ischar (varargin {1})) + ellip = select_phantom (varargin {1}); + elseif (numel (varargin {1}) == 1) + n = varargin {1}; + else + if (size (varargin {1}, 2) != 6) + error ("Wrong number of columns in user phantom"); + endif + ellip = varargin {1}; + endif + elseif (nargin == 2) + n = varargin {2}; + if (ischar (varargin {1})) + ellip = select_phantom (varargin {1}); + else + if (size (varargin {1}, 2) != 6) + error ("Wrong number of columns in user phantom"); + endif + ellip = varargin {1}; + endif + elseif (nargin > 2) + warning ("Extra arguments passed to phantom were ignored"); + endif +endfunction + +function e = select_phantom (name) + if (strcmpi (name, 'shepp-logan')) + e = shepp_logan (); + elseif (strcmpi (name, 'modified shepp-logan')) + e = mod_shepp_logan (); + else + error ("Unknown phantom type: %s", name); + endif +endfunction + +function e = shepp_logan + # Standard head phantom, taken from Shepp & Logan + e = [ 2, .69, .92, 0, 0, 0; + -.98, .6624, .8740, 0, -.0184, 0; + -.02, .1100, .3100, .22, 0, -18; + -.02, .1600, .4100, -.22, 0, 18; + .01, .2100, .2500, 0, .35, 0; + .01, .0460, .0460, 0, .1, 0; + .02, .0460, .0460, 0, -.1, 0; + .01, .0460, .0230, -.08, -.605, 0; + .01, .0230, .0230, 0, -.606, 0; + .01, .0230, .0460, .06, -.605, 0]; +endfunction + +function e = mod_shepp_logan + # Modified version of Shepp & Logan's head phantom, + # adjusted to improve contrast. Taken from Toft. + e = [ 1, .69, .92, 0, 0, 0; + -.80, .6624, .8740, 0, -.0184, 0; + -.20, .1100, .3100, .22, 0, -18; + -.20, .1600, .4100, -.22, 0, 18; + .10, .2100, .2500, 0, .35, 0; + .10, .0460, .0460, 0, .1, 0; + .10, .0460, .0460, 0, -.1, 0; + .10, .0460, .0230, -.08, -.605, 0; + .10, .0230, .0230, 0, -.606, 0; + .10, .0230, .0460, .06, -.605, 0]; +endfunction + +#function e = ?? +# # Add any further phantoms of interest here +# e = [ 0, 0, 0, 0, 0, 0; +# 0, 0, 0, 0, 0, 0]; +#endfunction + +%!demo +%! P = phantom (512); +%! imshow (P, []); + +% $Log$ +% 2010/03/31 lxop +% Fixed typo that prevented the use of custom ellipses +% +% 2010/03/09 lxop +% First commit to repo