function [intersect,indx,t,u,v] = ray_intersect(fv,r,d,cull);
%RAY_INTERSECT - find intersection of a ray with a set of faces
% function [intersect,indx,t,u,v] = ray_intersect(fv,r,d,cull);
% function indx = ray_intersect(fv,r,d,cull);
% FV is structure of faces and vertices, in standard Matlab convention
% of m x 3 triangular faces and n x 3 corresponding vertices
% R is the observation point, with unit directional vector D.
% (D will be automatically scaled to unity in the program.)
% CULL is one of 'i', 'o', or 'b', for inside (backside),
% outside (frontside), or both sides
% to use in the calculation. Only first letter of CULL is
% considered, so 'inside', 'outside', or 'both' are acceptable.
% For example, 'i' means the ray passes from the observation r in the
% direction d through the backside and out the frontside of the triangle.
% In other words, the observation is inside a closed surface.
% Note that the program assumes right-handedness in the faces description,
% such that V0 to V1 to V2 points "out" or "front."
% OUTPUT:
% INDX is the index such that fv.faces(INDX,:) gives the
% triangles intersected.
% T,U,V gives distance t to the point(s) of
% intersection, with point(s) of intersection given by
% INTERSECT = (1-U-V) .* fv.vertices(fv.faces(indx,1),:) + ...
% U .* fv.vertices(fv.faces(indx,2),:) + V .* fv.vertices(fv.faces(indx,3),:);
% where U is repmat(u,1,3) and similarly V.
%
% See also delaunay_surface, surf_distance, ray_intersect_slow
% ---------------------- 14-Jun-2004 17:12:33 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Utility - Numeric
%
% Alphabetical list of external functions (non-Matlab):
% toolbox\colnorm.m
%
% At Check-in: $Author: Mosher $ $Revision: 4 $ $Date: 6/14/04 3:38p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-Jun-2004
%
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
% University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
% Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
% CNRS, Hopital de la Salpetriere, Paris, France
%
% See BrainStorm website at http://neuroimage.usc.edu for further information.
%
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
% ------------------------ 14-Jun-2004 17:12:33 -----------------------
% Uses Tomas Moeller and Ben Trumbore's "Fast, Minimum Storage
% Ray/Triangle Intersection" 1997approach. Given a line defined by
% L(t) = R + t * D and a triangle vertex V0 and edges E1 and E2
% emanating from VO. Then the point of intersection between the
% ray and the plane defined by the triangle is
% [-D, E1, E2] [t;u;v] = R - V0;
% If u and v are both 0 <= u,v <= 1, then the intersection
% is in the bounds of the triangle. The point on the triangle
% is given by T(u,v) = (1 - u - v)*V0 + u*V1 + v*V2, where
% E1 = V1 - V0, and E2 = V2 - V0. {u,v} are the barycentric
% coordinates of the triangle.
% The solution is simply Cramer's rule,
% [t] = {T,E1,E2} /
% [u] = {-D,T,E2} / {-D,E1,E2}.
% [v] = {-D,E1,T} /
% where T is the translation R - V0, and {a,b,c} is the
% triple scalar product.
% Moeller and Trumbore's trick is to first cull the set of triangles
% looking for determinants in the correct direction. In this reduced
% subset, they look for triangles whose u is properly bounded. In the
% greatly reduced subset, the look for triangles whose v is properly
% bounded. The finally apply the division to get the units correct.
%
%
% Copyright Dr. John C. Mosher, Los Alamos National Laboratory
% 28 December 2001
tol = 1e-6; % tolerance for parallel triangles,
% will be scaled below to maximum triangle length
EPS = 1e-6; % tolerance in the barycentric coordinates
if(~exist('cull','var')),
cull = [];
end
if(isempty(cull)),
cull = 'o'; % outside
end
d = d(:)/norm(d(:)); % force to column unity
r = r(:); % force to column
% culling can be 'o' outside (front),'i', inside (back), or 'b' both
cull = lower(cull(1));
mfaces = size(fv.faces,1); %number of triangles
% fv is structure of faces and vertices, in matlab convention of m x 3 and n x 3
% define the edges first
% edge from vertext 0 to vertex 1 first. then vertex 2
e1 = (fv.vertices(fv.faces(:,2),:) - fv.vertices(fv.faces(:,1),:))'; % 3 x m
e2 = (fv.vertices(fv.faces(:,3),:) - fv.vertices(fv.faces(:,1),:))'; % 3 x m
% want to economically set the tolerance to something relative to the scale of the
% triangles. Calculating the area would require expensive cross product
% area = e1 x e2 / 2. Just use side one as good enough
tol = max(colnorm(e1))*tol;
% partial calculation using Moeller's "p" and "q"
% p is the d x e2
p = cross(d(:,ones(1,mfaces)),e2);
% e1 dot d cross e2 is the determinant of each triangle
% (note determinant is twice the area of the triangle).
Det = sum(e1 .* p); % determinant of each triangle
% if the determinant is (almost) zero, then the
% ray is passing (almost) parallel to the surface of the
% triangle. So there is (almost) no intersection.
% The ray has a direction. If the direction of the ray hits
% the backside of the triangle, then the ray and the triangle
% are both pointing in the same direction, and we say that
% the triangle is passing from the inside to the outside or
% backside to frontside.
% If the ray is pointing in the opposite direction, then
% the ray is passing from the outside to the inside (frontside
% to backside).
% If the determinant is negative, then we are passing through
% the backside of the triangle. If positive, then the
% frontside.
switch cull
case 'b' % both sides are considered
ndx = find(Det < -tol | Det > tol);
case 'i' % only those for which we are passing from inside to outside
ndx = find(Det < -tol);
case 'o' % passing from outside to inside
ndx = find(Det > tol);
end
% so ndx gives us a subset of triangles for further processing
% calculate distance from V0 to the observation for each triangle
T = r(:,ones(1,length(ndx))) - fv.vertices(fv.faces(ndx,1),:)';
% now solve for the first barycentric coordinate, unscaled
u = sum(T .* p(:,ndx)); % dot product
u = u ./ Det(ndx); % handles negative cases
% find properly bounded u
ndx2 = find(u >= -EPS & u <= 1+EPS);
if(isempty(ndx2)),
indx = [];
t = [];
u = [];
v = [];
return
end
% ndx2 is double referenced, such that ndx(ndx2) refers
% to the original triangle numbering. Should be a greatly
% reduced set of indices
% the other Moeller variable is "q", q = T x E1
% form only for valid u
q = cross(T(:,ndx2),e1(:,ndx(ndx2)));
% form the other barycentric coordinate
v = sum(d(:,ones(1,size(q,2))) .* q); % dot product
v = v ./ Det(ndx(ndx2));
% lookfor valid v
ndx3 = find(v >= -EPS & (u(ndx2) + v) <= 1+EPS);
% ndx3 may be null or singleton
if(isempty(ndx3)),
indx = [];
t = [];
u = [];
v = [];
return
else
indx = ndx(ndx2(ndx3)); % relative indexing
end
% indx gives us the original absolution indexing
% scale each solution by the determinant
t = sum(e2(:,indx) .* q(:,ndx3))' ./ Det(indx)'; % dot product
u = u(ndx2(ndx3))';
v = v(ndx3)';
U = repmat(u,1,3);
V = repmat(v,1,3);
intersect = (1-U-V).*fv.vertices(fv.faces(indx,1),:) + ...
U .* fv.vertices(fv.faces(indx,2),:) + V .* fv.vertices(fv.faces(indx,3),:);
intersect = intersect'; % one intersection per column
return