% each row is an ellipse.
ells = [1, 5, 3, 9, 0, 0, 0, 0, 0, 0
	1, 7, 17, 3, 10, 3, 0, 0, 0, 0
	3, 10, 1, 1, -10, -5, 0, pi/4, 0, 0
	2, 10, 4, 20, 0, 10, 0, pi/8, 2*pi/3, 0];

p = EllipsoidPhantom(ells);

xStart = [-25 -25 -25];
xStep = [.25 .25 .25];
xSize = [200 200 200];
f = p.eval(xStart, xStep, xSize);


for slice = [65, 101, 121]
	figure
	imagesc((0:xSize(1)-1) * xStep(1) + xStart(1), ...
		(0:xSize(2)-1) * xStep(2) + xStart(2), ...
		f(:,:,slice)', [0 5]);
	axis xy
	colorbar
	axis equal
	xlabel('x_0')
	ylabel('x_1')
	title(sprintf('phantom slice, x_2 = %g', xStep(3)*(slice-1) + xStart(3)));
end

%% take x-ray measurements
% setup the y grid, which is 2D
yStart = [-30 -30];
yStep = [.25 .25];
ySize = [300 300];

% set angles determining the projection directions
phi =   [0  0     0     pi/4  pi/2  pi/4];
theta = [0  -pi/2  pi/2  0     0     pi/4];
psi =   [0  0     0     0     0     0   ];
numProj = length(phi); % number of projection angles

% generate a basis orthogonal to the projection direction for
% each projection.
[r0, r1, r2] = Euler3D(phi, theta, psi);
P = reshape([r0(:)'; r1(:)'], 2, 3, numProj);

% take the measurements
g = p.xray(yStart, yStep, ySize, P);

% plot


for gInd = [1 2]
	figure
	imagesc( (0:ySize(1)-1)*yStep(1) + yStart(1), ...
		(0:ySize(2)-1)*yStep(2) + yStart(2), ...
		g(:,:,gInd)');
	axis xy
	axis equal
	colorbar
	xlabel(sprintf('y_0, (direction [%.0f %.0f %.0f])', r0(:,gInd)));
	ylabel(sprintf('y_1 = (direction [%.0f %.0f %.0f])', r1(:,gInd)));
	title(sprintf('projection along [%.0f %.0f %.0f]', r2(:,gInd)));
end

figure
volume_plot(xStart, xStep, f, .5);