-
Notifications
You must be signed in to change notification settings - Fork 21
/
cone_intersection.m
60 lines (57 loc) · 2.33 KB
/
cone_intersection.m
1
function [ in, rin ] = cone_intersection( r_in, e, rad1, rad2, sag, the, surf )%% rad1 should be < rad2. Sag is at rad1. The cone is opening up here% % find the cone's vertex coordinates and cone slice limits minsag = sag; if abs( the - pi ) < 1e-10 vx = sag; else vx = sag - rad1 / tan( the ); end maxsag = sag + ( rad2 - rad1 ) / tan( the ); if minsag == maxsag maxsag = minsag + realmin; % make the two tiny different to avoid rays disappearing at frontoparallel surfaces elseif minsag > maxsag tmp = minsag; minsag = maxsag; maxsag = tmp; end v = [ vx 0 0 ]; % cone vertex coordinates dt = r_in - repmat( v, size( r_in, 1 ), 1 ); % rays in the reference frame of the cone vertex at [ 0 0 0 ] a2 = tan( the ).^2;% if a2 == 0% a2;% end [ ~, ~, dm, dp ] = quadric_intersection( -1, a2, a2, 0, dt, 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; 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 r2 = sum( rin1( :, 2 : 3 ).^2, 2 ); if rad1 < rad2 in1 = r2 <= rad2^2 & r2 >= rad1^2 & rin1( :, 1 ) >= minsag & rin1( :, 1 ) <= maxsag; else in1 = r2 <= rad1^2 & r2 >= rad2^2 & rin1( :, 1 ) >= minsag & rin1( :, 1 ) <= maxsag; end % then try the second intersection rin2 = r_in + repmat( dm .* snd + dp .* fst, 1, 3 ) .* e; % use the closest intersection with the surface r2 = sum( rin2( :, 2 : 3 ).^2, 2 ); if rad1 < rad2 in2 = r2 <= rad2^2 & r2 >= rad1^2 & rin2( :, 1 ) >= minsag & rin2( :, 1 ) <= maxsag; else in2 = r2 <= rad1^2 & r2 >= rad2^2 & rin2( :, 1 ) >= minsag & rin2( :, 1 ) <= maxsag; end in = in1 | in2; % use either intersection rin = rin1; rin( ~in1 & in2, : ) = rin2( ~in1 & in2, : ); % only use second intersection if the first is invalidend