%% define the 'x' grid
% These are the locations in the reconstruction (image) where you want to
% evaluate the phantom. Useful for generating ground truth and debugging.
xStart = [-25 -25];
xStep = [.25 .25];
xSize = [200 200];

%% define ellipses
% each row is an ellipse, specified by
% [intensity, radis a, radius b, centerX, centerY, rotation]
% just like 'help phantom'
ells = [1, 10, 3, 11, 5, pi/8
	.4, 5, 5, 0, 0, 0
	2, 10, 7, -1, -10, 2*pi/3
	2 .5 .5 -5 7 0];


%% create the phantom object
p = EllipsoidPhantom(ells);
f = p.eval(xStart, xStep, xSize);

figure
% rows of f correspond to x_0, cols to x_1 , pages to x_2, so we need to
% transpose to plot it in the conventional way, with "x" going
% left-right and "y" going up-down
imagesc( (0:xSize(1)-1) * xStep(1) + xStart(1), ...
	(0:xSize(2)-1) * xStep(2) + xStart(2), ...
	f');
axis xy
colorbar
axis equal
xlabel('x_0')
ylabel('x_1')
title('phantom');

%% calculate a sinogram
numProj = 100; % number of projection angles

% P is a 1 x 2 x numProj matrix specifying the projection directions.
% Specifically, it is a collection of numProj 1 x 2 unit vectors giving the
% direction perpendicular to the direction of integration. This vector also
% specifies the coordinate system for the projection domain.
P = zeros(1, 2, numProj);
t = (0:numProj-1) * pi/numProj;
for i = 1:numProj
	P(:,:,i) = [cos(t(i)) sin(t(i))];
end

% define the 'y' grid, which are coordinates in the projection domain
yStart = -25;
yStep = .25;
ySize = 200;
g = p.xray(yStart, yStep, ySize, P);

figure;
imagesc(g, 'XData', t, 'YData', (0:ySize-1)*yStep + yStart)
xlabel('angle perpendicular to projection direction')
ylabel('y')
axis xy
colorbar
title('sinogram');