function [fan_proj,dat] = CTFanProjSim(zsource, ddetector_angle_degrees, varargin)
% function [fan_proj,dat] = CTFanProjSim(zsource, ddetector_angle_degrees, varargin)
%  Computes the fan beam 2D projections for an 'arc' detector
%  The source and detector are assumed to be a rigid body and rotate around the origin
%  The detector angle and the source motion increments are assumed to be equal
%  so we can use the fast fan-to-parallel algorithm of Kak and Slaney Sec. 3.4.3
%  The detector arc has a center at zsource and the number of angles is set to span the object
%  The source scans 360 degrees. 
%  inputs:
%     zsource: the initial source position (x+iy)
%     ddetector_angle: (degrees) the detector angle and the source motion increment--assumed to be equal
%     detector_angle_multiplier: (optional = 1.2) increase detector angle range so can span object
%     do_parallel: (optional = false) if present computes parallel projections and outputs in dat struct
%     the following optional parameters specify the object. At least one must be present
%        'ellipse': a matrix of 6 or more columns (density can have more than one value) : 
%            one row per matrix
%            [center_x, center_y, semi_major_axis,semi_minor_axis,angle(degrees), density(1 ...) ]
%        'gaussian': a matrix with same columns as for ellipse. Projections of exp(-pi*(x^2/a^2 + y^2/b^2))
%            range of object is set to nsigma*a nsigma*b
%        'nsigma': (optional=2.5) multiplier of a,b in gaussian to give detector length
%        'polygon': a convex polygon specified by its vertexes. Specified as a struct array with each struct having fields
%            vertex: a complex vector giving the polygon vertex coordinates. 
%               Assumed closed (i.e. last vertex = first) but will be closed if not
%            density: a vector giving the density values
%    outputs:
%      fan_proj: a matrix of the projections. (ndetectors, nangles, nbasis)
%      dat: intermediate calculation results
%
%  by Robert E. Alvarez http://aprendtech.com
%  9/13/10 - 2011-Sep-19 16:13

nreqargs = 2;
assert(nargin>=nreqargs);
assert(numel(zsource)==1);
ddetector_angle = ddetector_angle_degrees*pi/180; % convert to radians
nsigma = 2.5; 
do_parallel = false; 
detector_angle_multiplier = 1.2; % make sure can span object with projections
zcorners4 = []; % 4 corners of bounding rectangle

do_ellipse = false;
do_polygon = false;
do_gaussian = false;
nbasis = [];
  % process the optional arguments
if(nargin>nreqargs)
  i=1;
  while(i<=size(varargin,2))
     switch lower(varargin{i})
     case 'do_parallel';  do_parallel = true;  
     case 'nsigma';  nsigma = varargin{i+1};  i = i+1;
     case 'detector_angle_multiplier';  detector_angle_multiplier = varargin{i+1};  i = i+1;
     case 'zdetectors';  zdetectors = varargin{i+1};  i = i+1;
     case 'ndetectors';  ndetectors = varargin{i+1};  i = i+1;
     case 'ellipse'; % only ellipse implemented for now
%    ellipse.spec  =      [center_x, center_y, major_axis,minor_axis,angle(degrees), density(1 ...) ]
         if do_ellipse, error('CTprojSim: only one set of ellipses can be defined'); end
             % store data in struct to keep it together         
         ellipse.spec = varargin{i+1}; 
         ellipse.nshapes = size(ellipse.spec,1);
            % extract the vector 'density'. all shapes must have same number of dimensions
         ellipse.density = ellipse.spec(:,6:end);
         nb = size(ellipse.density,2);
         if isempty(nbasis)
            nbasis = nb;
         elseif nbasis ~= nb
            error('CTprojSim: nbasis changes from %d to %d in ellipse ellipse.spec',nbasis,nb);
         end
         ellipse.zcenters = zcomp(ellipse.spec(:,1:2));
         ellipse.as = ellipse.spec(:,3);
         ellipse.bs = ellipse.spec(:,4);
         ellipse.angles =ellipse.spec(:,5)*pi/180; % convert ellipse angles to radians

            % find bounding rectangle of all ellipses
         [ellipse.zrange,ellipse.sm2] = EllipseCornersLoc(ellipse,1);
         zcorners4 = zrange([ ellipse.zrange; zcorners4]); 
         dat.ellipse = ellipse;
         do_ellipse = true;
         i=i+1;
         
           
    otherwise
      error('CTprojSim: Unknown argument %s given',varargin{i});
     end
     i=i+1;
  end
end

if ~ (do_ellipse || do_polygon || do_gaussian)
   error('CTprojSim: Need to specify a set of shapes');
end

R = max(abs(zcorners4));  % radius of circle, centered on origin that includes all the objects
S = abs(zsource); % distance from origin to source
assert( R < S ); % source cannot be in object
detector_spacing = S*ddetector_angle; % use this to scale density

  % use detector_angle for detector angles and source_angle for source angles 
detector_angle_max = detector_angle_multiplier*asin(R/S); % asin == arcsine
   % detector angles--center is at source
detector_angles = angle(-zsource) + (-detector_angle_max:ddetector_angle:detector_angle_max)'; 
    % note transpose for column matrix
    
ndetector_angles = numel(detector_angles);

source_angles = 0:ddetector_angle:(2*pi); % source positions + 1
source_angles = source_angles(1:end-1); % do not wrap around
nsource_angles = numel(source_angles);

fan_proj = zeros(ndetector_angles, nsource_angles,nbasis); % fan projections

if do_ellipse
%           compute the length of lines through ellipses using formula from my blog post:
%            http://aprendtech.com/wordpress/?p=89
%            T = 2*ab*sqrt(sm^2 - s^2)/sm^2  s<sm
%             ddetector_angle = detector_angle - ellipse_detector_angle --detector_angle = angle of perpendicular to line
%            sm^2 = a^2*(cos(ddetector_angle))^2 + b^2*(sin(ddetector_angle))^2, a,b, are major,minor semi-axes
%          compute projections for all ellipses and all projection angles
     abs_prod = 2*prod([ellipse.as ellipse.bs],2);
     ellipse.density = ellipse.density/detector_spacing; % scale to get accurate reconstruct values
     for kangle = 1:nsource_angles
       src_angle = source_angles(kangle);
       pj_angles_rot = detector_angles + src_angle;
       nhats = exp(1i*pj_angles_rot);
       nhats_perp = 1i*nhats;
       zsource_rot = zsource*exp(1i*src_angle);
       for k = 1:ellipse.nshapes
         ddetector_angles = pj_angles_rot - ellipse.angles(k); % ddetector_angles for this ellipse for all projection lines
         sm2 = EllipseSm2(ellipse.as(k), ellipse.bs(k), ddetector_angles); % sm^2 for this ellipse for all projection lines
              % distance of projection lines from ellipse center
         Rs = zdot(nhats_perp,zsource_rot - ellipse.zcenters(k)) ;
         ds2 = sm2 - Rs.^2;
         idxs = find(ds2>=0);
         if isempty(idxs) % put up debugging information
             warning('CTFanProjSim: object out of field');
             fh = figure;
             axis([-3 3 -3 3])
             hold on
             plot(0,0,'+g') % origin
             plot(zsource_rot + 1i*eps,'xk');
             plot(ellipse.zcenters(k) + 1i*eps,'or');
             for kline = 1:2:numel(nhats)
                 plot([zsource_rot; zsource_rot + nhats(kline)*2.5],'-b');
                 plot([ellipse.zcenters(k) + 1i*eps, ellipse.zcenters(k) + 1i*eps + nhats_perp(kline)*Rs(kline)],'-r');
             end
             hold off
         end
         Ts = abs_prod(k)*sqrt( ds2(idxs) )./sm2(idxs); % length of intersection of line with ellipse
         Ls = repmat(Ts(:),1,nbasis).*repmat(ellipse.density(k,:),numel(idxs),1); % length*density
         fan_proj(idxs,kangle,:) = fan_proj(idxs,kangle,:) + Ls;
       end    
    end
end % do_ellipse

if do_parallel
    sz = size(fan_proj);
    pp = NaN*ones(sz);
    ns = (1:nsource_angles) -1; 
    for kr = (1:ndetector_angles) -1 
        idx2para = sub2ind(sz, (kr+1)*ones(size(ns)),mod(kr + ns,nsource_angles) + 1);
        pp(idx2para) = fan_proj(kr+1,:);
    end % for kr = 1:ndetector_angles
    dat.para_proj = pp;
    para_angles = linspace(-detector_angle_max, -detector_angle_max + 2*pi, nsource_angles+1);
    dat.para_angles = para_angles(1:end-1)*180/pi; % in degrees
end % if do_parallel



%--------------  LOCAL FUNCTIONS ---------------
    function [zrng,sm2] = EllipseCornersLoc(s,nsigma) % find overall rectangle enclosing a set of ellipses. 
      % ellipse parameters are fields of s
             sm2 = EllipseSm2(nsigma*s.as,nsigma*s.bs,s.angles); % in x-dir
             sm2_perp =  EllipseSm2(nsigma*s.as,nsigma*s.bs,s.angles + pi/2); % in y-dir
             zcorners_rect = zcomp(sqrt([sm2 sm2_perp]));
             zcorners_rect = [zcorners_rect, -zcorners_rect];
                % find range after offset by zcenters
             zrng = zrange(repmat(s.zcenters,1,2) + zcorners_rect);

    function sm2 = EllipseSm2(a,b,ddetector_angle) % square of max offset of tangent line for ellipse
      sm2 = a.^2.*(cos(ddetector_angle)).^2 + b.^2.*(sin(ddetector_angle)).^2;  
  
     function d = zdot(z1,z2) % dot product for complex vectors
    d = real( z1(:) .* conj(z2(:)));
