%  Make data for a circular phantom for CT reconstruction
%
%  Calling: ct_data = pht_circle (K);
%
%  Input:  K - Number of projection to be calculated
%
%  Output: ct_data - CT data where first index is xm and second
%                    index is the projection angle
%
%  Version 1.1, 30/6-1997, Joergen Arendt Jensen
%  Version 1.2, 4/12-2006, Joergen Arendt Jensen
%      Error in boundary calculation fixed (floor replaced by round)

function ct_data = pht_circle (K)

%  Setup parameters

N=128;       %  Total number of samples
Ndiv2=N/2;   %  Half the total number
Nr=30;       %  Radius of circle

%  Calculate the projection of the circle

A=Nr;   %  Major axis
B=Nr;     %  Minor axis
rho=1;    %  Values inside ellipse
x1=0;     %  Offset in x-direction
y1=0;     %  Offset in y-direction
alpha=0;  %  Rotation with x-axis

%  Preallocate the data structure

ct_data=zeros(N,K);

%  Generate the data

index=1;
s=sqrt(x1^2+y1^2);
gamma=atan2(y1,x1);
for theta1=0:pi/K:(pi-pi/(2*K))
  theta = theta1 - alpha;
  a2theta = (A*cos(theta))^2 + (B*sin(theta))^2;

  toffset = s*cos(gamma-theta1);
  tboundary = round(sqrt(a2theta));
%   tboundary = floor(sqrt(a2theta));
  tvalues = (1:N) - Ndiv2 + toffset;
  tinside = abs(tvalues) < tboundary;

  P = 2*rho*A*B/a2theta*sqrt(a2theta-tvalues.^2);
  data = P .* tinside;

  ct_data(:,index) = data';
  index = index + 1;
  end





