function phantom_2(ProjectFile, rho1, rho2, nproj, npos, PNG, EPS, MAT)
%---------------------------------------
%
%  To generate images and projections 
%  using ellipses for testing
%  tomography reconstruction codes
%
%  E.g. input:
%
%  rho1 = -1;        % minimum rho
%  rho2 = 1;         % maximum rho
%  nproj = 110;      % number of angles
%  npos  = 127;      % number of rays
%  PNG = 0;          % 1 save png, 0 don't save
%  EPS = 0;          % 1 save eps, 0 don't save
%  MAT = 0;          % 1 save mat, 0 don't save
%
%  phantom_2('sino', -1, 1, 110, 127, 0, 0, 0);
%
%  Kai Hock
%  Cockcroft Institute
%  10 June 2010
%---------------------------------------

%--describe  image -----
%
% Each row describes an ellipse.
% 'Refractive index' refers to intensity.
% At points where ellipses,overlap, the intensities are cumulated.
%
%-----------------------------------------------------------
%    center          major   minor   rotation   refractive
%    coordinate      axis    axis    angle      index
%     x      y                       (deg)

%  Shepp-Logan, Kak and Slaney (1988), p. 55
AI1 = [0      0       0.92    0.69     90        2.0
      0     -0.0184  0.874   0.6624   90       -0.98
      0.22   0       0.31    0.11     72       -0.02
     -0.22   0       0.41    0.16    108       -0.02
      0      0.35    0.25    0.21     90        0.01
      0      0.1     0.046   0.046    0         0.01
      0     -0.1     0.046   0.046    0         0.01
     -0.08  -0.605   0.046   0.023    0         0.01
      0     -0.605   0.023   0.023    0         0.01
      0.06  -0.605   0.046   0.023    90        0.01];

% Shepp-Logan, 
% Murrell, "Computer-Aided Tomography," 
% The Mathematical Journal (1996), p. 60
AI2 = [0      0       0.92    0.69     90        2.0
      0     -0.0184  0.874   0.6624   90       -0.9
      0.22   0       0.31    0.11     72       -0.1
     -0.22   0       0.41    0.16    108       -0.1
      0      0.35    0.25    0.21     90        0.3
      0      0.1     0.046   0.046    0         0.3
      0     -0.1     0.046   0.046    0         0.3
     -0.08  -0.605   0.046   0.023    0         0.3
      0     -0.605   0.023   0.023    0         0.3
      0.06  -0.605   0.046   0.023    90        0.3];

%AI = AI1;
AI = AI2;

%---image size -----------

xlow  = -1; % image left
xhigh =  1; % image right
ylow  = -1; % image bottom
yhigh =  1; % image top
xstep = 0.01; % x resolution
ystep = 0.01; % y resolution

%------------------------------------------------------------

x1 = AI(:, 1);   % centre coordinate x
y1 = AI(:, 2);  % centre coordinate y
A  = AI(:, 3);   % semi major axis
B  = AI(:, 4);   % semi minor axis
a1 = AI(:, 5);   % (deg) rotation angle
ri = AI(:, 6);   % refractive index
n1 = length(x1); % number of ellipses

%----intensity matrix -------

x0 = xlow:xstep:xhigh;
y0 = ylow:ystep:yhigh;
nx = length(x0);
ny = length(y0);
[yy, xx] = ndgrid(-y0, x0);
f0 = zeros(ny, nx);

for i1 = 1:n1,
  alpha1 = a1(i1)/180*pi;
  x =  (xx-x1(i1))*cos(alpha1) + (yy-y1(i1))*sin(alpha1);
  y = -(xx-x1(i1))*sin(alpha1) + (yy-y1(i1))*cos(alpha1);
  f0 = f0 + ri(i1)*heaviside(1 - (x/A(i1)) .^2 - (y/B(i1)) .^2);
end
fmax = max(max(f0));

close all
figure
imagesc(x0, -y0, f0/fmax, [0 1])
set(gca,'YDir','Normal')
colormap(gray)
axis square
axis([xlow xhigh ylow yhigh])
set(gca, 'FontSize', 24, 'FontWeight', 'b')
xlabel('x','FontSize', 24, 'FontWeight', 'b')
ylabel('y','FontSize', 24, 'FontWeight', 'b')

if PNG == 1,
    fprintf('Saving image to png\n');
    print('-dpng', 'phantom','-r360');
end
if EPS == 1,
    fprintf('Saving image to eps\n');
    print('-depsc', 'phantom','-r360');
end

%--------------------------------------------------------------------
%
% ProjectFile data format:
%
% nproj     - number of projections;
% angles    - column vector [nproj x 1] of angles
% weights   - column vector [nproj x 1] of weights for each projection.
%             For linear interpolation, use the sum of the angular intervals
%              on the two sides of each projection as the weight.
%             (Set to ones by default.) 
% npos      - number of positions, must be the same for each projection
% positions - matrix [nproj x npos] for positions of each and every projection
%             (for each projection, the positions must have uniform intervals, 
%              and the same number of points)
% centre    - column vector [nproj x 1], 
%             the position of the centre of rotation for each projection
%             (set to zeros by default) 
% sinogram  - matrix of projection data [nproj x npos], 
%             where npos = number of positions across the projections
%
%--------------------------------------------------------------------

%-- set projection parameters ------------

centre  = zeros(nproj, 1);
weights = ones(nproj, 1) * pi/nproj;
angles = (1:nproj)' / nproj * pi;

drho = (rho2 - rho1) / (npos-1);
rho  = rho1 + (0:(npos-1)) * drho;
for i = 1:nproj,
  positions(i, :) = rho;
end

%----- generate projections ------------------------
fprintf('Generating projections\n');

sinogram = zeros(nproj, npos);

for iproj = 1:nproj,
    p1 = zeros(1, npos);
    
    for i1 = 1:n1,
        theta = angles(iproj) - a1(i1)/180*pi;
        a2 = sqrt((A(i1)*cos(theta))^2 + (B(i1)*sin(theta))^2);
        
        s1 = sqrt(x1(i1)^2 + y1(i1)^2);
        gamma1 = atan2(y1(i1), x1(i1));
        
        for ipos = 1:npos,
            t1 = rho(ipos) - s1 * cos(gamma1 - angles(iproj));
            if (abs(t1) < a2),
                p1(ipos) = p1(ipos) + 2*ri(i1)*A(i1)*B(i1)/a2^2 * sqrt(a2^2 - t1^2);
            end
        end
    end
    sinogram(iproj, :) = p1;
end

%--------------------------------

figure
imagesc(rho, angles/pi*180, sinogram/max(max(sinogram)), [0 1]);
set(gca,'YDir','Normal')
colormap(gray)
axis square
set(gca, 'FontSize', 24, 'FontWeight', 'b')
xlabel('\rho','FontSize', 24, 'FontWeight', 'b')
ylabel('\theta (deg)','FontSize', 24, 'FontWeight', 'b')

if PNG == 1,
    fprintf('Saving image to png\n');
    print('-dpng', 'sino','-r360');
end
if EPS == 1,
    fprintf('Saving image to eps\n');
    print('-depsc', 'sino','-r360');
end
%------------------------------------------

if MAT == 1,
fprintf('Saving sinogram data\n');
save(ProjectFile, 'nproj', 'angles', 'weights', ...
                  'npos', 'positions', 'centre', ...
                  'sinogram');
end