cylinder_intersection.m

function [ in, rin ] = cylinder_intersection( r_in, e, sD, sh, surf )
%   c2 = e( :, 2 ).^2 + e( :, 3 ).^2; % d^2 coefficients
%   c1 = 2 * ( e( :, 2 ) .* r_in( :, 2 ) + e( :, 3 ) .* r_in( :, 3 ) ); % d^1 coefficients
%   c0 = r_in( :, 2 ).^2 + r_in( :, 3 ).^2 - ( sD / 2 ).^2; % d^0 coefficients
% 
%   % solve c2 * d^2 + c1 * d + c0 = 0 for d
%   D =  c1.^2 - 4 * c2 .* c0;
%   D( D < -1e-12 ) = Inf; % suppress negative values (no intersection with the half-cone) to avoid imaginary values
%   D( D < 0 ) = 0; % remove round off negative Ds
%   D = sqrt( D );
%   dm = 0.5 * ( -c1 - D ) ./ c2;
%   dp = 0.5 * ( -c1 + D ) ./ c2;

  a2 = ( sD/2 ).^2;
  [ ~, ~, dm, dp ] = quadric_intersection( 0, a2, a2, 1, r_in, e ); % orient the reference frame along the z direction instead of the x direction
  % discount the previous intersection, where the ray originates
  if isa( surf, 'Screen' ) || isa( surf, 'Retina' ) % allow negative ray directions to use with virtual images
      dm( abs( dm ) < 1e-12 ) = realmax;
      dp( abs( dp ) < 1e-12 ) = realmax;
  else
      dm( dm < 1e-12 ) = realmax;
      dp( dp < 1e-12 ) = realmax;
  end

  % form the intersection vectors, determine which is the first intersection, which is the second
  fst = dm < dp; % don't count the ray's origin on a surface;
  snd = ~fst;
  % first try the first intersection
  rin1 = r_in + repmat( dm .* fst + dp .* snd, 1, 3 ) .* e; % use the closest intersection with the surface
  in1 = rin1( :, 1 ) >= 0 & rin1( :, 1 ) <= sh;
  % then try the second intersection
  rin2 = r_in + repmat( dm .* snd + dp .* fst, 1, 3 ) .* e; % use the closest intersection with the surface
  in2 = rin2( :, 1 ) >= 0 & rin2( :, 1 ) <= sh;
  in = in1 | in2; % use either intersection
  rin = rin1;
  rin( ~in1 & in2, : ) = rin2( ~in1 & in2, : ); % only use second intersection if the first is invalid
end