ray_intersect

(Toolbox/ray_intersect.m in BrainStorm 2.0 (Alpha))

Function Synopsis

[intersect,indx,t,u,v] = ray_intersect(fv,r,d,cull);

Help Text

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

Cross-Reference Information

colnorm                            C:\BrainStorm_2001\Toolbox\colnorm.m

bem_xfer                           C:\BrainStorm_2001\Toolbox\bem_xfer.m

Listing of function C:\BrainStorm_2001\Toolbox\ray_intersect.m

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

%<autobegin> ---------------------- 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.
%<autoend> ------------------------ 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>
% <copyright>
% 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