shepp-logan模型

function [pellipse]=phantom3d（varargin）

%PHANTOM3D Three-dimensional analogue of MATLAB Shepp-Logan phantom
%   P = PHANTOM3D（DEFN） 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（EN） 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 [-11].  Columns 2 through 7 must be specified in terms
%   of this range.
%
%   [PE] = 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 [01].
%    
%   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