function [p,ellipse]=phantom(varargin)
%PHANTOM Generate a head phantom image.
%   P = PHANTOM(DEF,N) generates an image of a head phantom that can   
%   be used to test the numerical accuracy of RADON and IRADON or other  
%   2-D reconstruction algorithms.  P is a grayscale intensity image that
%   consists of one large ellipse (representing the brain) containing
%   several smaller ellipses (representing features in the brain).
%
%   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 number of rows and columns in P.
%   If you omit the argument, N defaults to 256.
% 
%   P = PHANTOM(E,N) generates a user-defined phantom, where each row
%   of the matrix E specifies an ellipse in the image.  E has six columns,
%   with each column containing a different parameter for the ellipses:
%   
%     Column 1:  A    the additive intensity value of the ellipse
%     Column 2:  a    the length of the horizontal semi-axis of the ellipse 
%     Column 3:  b    the length of the vertical semi-axis of the ellipse
%     Column 4:  x0   the x-coordinate of the center of the ellipse
%     Column 5:  y0   the y-coordinate of the center of the ellipse
%     Column 6:  phi  the angle (in degrees) between the horizontal semi-axis 
%                     of the ellipse and the x-axis of the image        
%
%   For purposes of generating the phantom, the domains for the x- and 
%   y-axes span [-1,1].  Columns 2 through 5 must be specified in terms
%   of this range.
%
%   [P,E] = PHANTOM(...) 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 pixel in the output image, the pixel's value is equal to the
%   sum of the additive intensity values of all ellipses that the pixel is a 
%   part of.  If a pixel is not part of any ellipse, its value is 0.  
%
%   The additive intensity value A for an ellipse can be positive or negative;
%   if it is negative, the ellipse will be darker than the surrounding pixels.
%   Note that, depending on the values of A, some pixels may have values outside
%   the range [0,1].
%    
%   See also RADON, IRADON.

%   Chris Griffin 4-22-97
%   Copyright 1993-1998 The MathWorks, Inc.  All Rights Reserved.
%   $Revision: 1.5 $  $Date: 1997/11/24 16:18:24 $

%   References: 
%      A. K. Jain, "Fundamentals of Digital Image Processing", p. 439.
%      P. A. Toft, "The Radon Transform, Theory and Implementation" (unpublished
%      dissertation), p. 199.

[ellipse,n] = parse_inputs(varargin{:});

p = zeros(n);

xax =  ( (0:n-1)-(n-1)/2 ) / ((n-1)/2); yax=xax; % x and y axes
xg = repmat(xax, n, 1);   % x coordinates, the y coordinates are rot90(xg)

for k = 1:size(ellipse,1)    
   asq = ellipse(k,2)^2;       % a^2
   bsq = ellipse(k,3)^2;       % b^2
   phi = ellipse(k,6)*pi/180;  % rotation angle in radians
   x0 = ellipse(k,4);          % x offset
   y0 = ellipse(k,5);          % y offset
   A = ellipse(k,1);           % Amplitude change for this ellipse
   x=xg-x0; y=rot90(xg)-y0;    % Center the ellipse
   cosp = cos(phi); sinp = sin(phi);
   idx=find(((x.*cosp + y.*sinp).^2)./asq + ((y.*cosp - x.*sinp).^2)./bsq <= 1); 
   p(idx) = p(idx) + A;
end
   
   
function [e,n] = parse_inputs(varargin)
%  e is the m-by-6 array which defines ellipses
%  n is the size of the phantom brain image

n=256;     % The default size
e = [];
defaults = {'shepp-logan', 'modified shepp-logan'};

for i=1:nargin
   if ischar(varargin{i})         % Look for a default phantom
      def = lower(varargin{i});
      idx = strmatch(def, defaults);
      if isempty(idx)
         error('Unknown default phantom selected.');
      end
      switch defaults{idx}
      case 'shepp-logan'
         e = shepp_logan;
      case 'modified shepp-logan'
         e = modified_shepp_logan;
      end
   elseif prod(size(varargin{i}))==1 
      n = varargin{i};            % a scalar is the image size
   elseif ndims(varargin{i})==2 & size(varargin{i},2)==6 
      e = varargin{i};            % user specified phantom
   else
      error('Invalid input arguments.');
   end
end

if isempty(e)                    % ellipse is not yet defined
   e = modified_shepp_logan;
end

   

      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Default head phantoms:   
%

function shep=shepp_logan
%
%  This is the default head phantom, taken from AK Jain, 439.
%
%         A    a     b    x0    y0    phi
%        ---------------------------------
shep = [  1   .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
         .01 .0460 .0460   0   -.1     0
         .01 .0460 .0230 -.08  -.605   0 
         .01 .0230 .0230   0   -.606   0
         .01 .0230 .0460  .06  -.605   0   ];
      
      
function toft=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    x0    y0    phi
%        ---------------------------------
toft = [  1   .69   .92    0     0     0   
        -.8  .6624 .8740   0  -.0184   0
        -.2  .1100 .3100  .22    0    -18
        -.2  .1600 .4100 -.22    0     18
         .1  .2100 .2500   0    .35    0
         .1  .0460 .0460   0    .1     0
         .1  .0460 .0460   0   -.1     0
         .1  .0460 .0230 -.08  -.605   0 
         .1  .0230 .0230   0   -.606   0
         .1  .0230 .0460  .06  -.605   0   ];