raytrace

 [t,u,v,idx]=raytrace(p0,v0,node,face)

 perform a Havel-styled ray tracing for a triangular surface

 author: Qianqian Fang, <q.fang at neu.edu>

 input:
   p0: starting point coordinate of the ray
   v0: directional vector of the ray
   node: a list of node coordinates (nn x 3)
   face: a surface mesh triangle list (ne x 3)

 output:
   t: signed distance from p to the intersection point for each surface
      triangle, if ray is parallel to the triangle, t is set to Inf
   u: bary-centric coordinate 1 of all intersection points
   v: bary-centric coordinate 2 of all intersection points
      the final bary-centric triplet is [u,v,1-u-v]
   idx: optional output, if requested, idx lists the IDs of the face
      elements that intersects the ray; users can manually calc idx by

      idx=find(u>=0 & v>=0 & u+v<=1.0 & ~isinf(t));

 Reference: 
  [1] J. Havel and A. Herout, "Yet faster ray-triangle intersection (using 
          SSE4)," IEEE Trans. on Visualization and Computer Graphics,
          16(3):434-438 (2010)
  [2] Q. Fang, "Comment on 'A study on tetrahedron-based inhomogeneous 
          Monte-Carlo optical simulation'," Biomed. Opt. Express, (in
          press)

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function [t,u,v,idx]=raytrace(p0,v0,node,face)
0002 %
0003 % [t,u,v,idx]=raytrace(p0,v0,node,face)
0004 %
0005 % perform a Havel-styled ray tracing for a triangular surface
0006 %
0007 % author: Qianqian Fang, <q.fang at neu.edu>
0008 %
0009 % input:
0010 %   p0: starting point coordinate of the ray
0011 %   v0: directional vector of the ray
0012 %   node: a list of node coordinates (nn x 3)
0013 %   face: a surface mesh triangle list (ne x 3)
0014 %
0015 % output:
0016 %   t: signed distance from p to the intersection point for each surface
0017 %      triangle, if ray is parallel to the triangle, t is set to Inf
0018 %   u: bary-centric coordinate 1 of all intersection points
0019 %   v: bary-centric coordinate 2 of all intersection points
0020 %      the final bary-centric triplet is [u,v,1-u-v]
0021 %   idx: optional output, if requested, idx lists the IDs of the face
0022 %      elements that intersects the ray; users can manually calc idx by
0023 %
0024 %      idx=find(u>=0 & v>=0 & u+v<=1.0 & ~isinf(t));
0025 %
0026 % Reference:
0027 %  [1] J. Havel and A. Herout, "Yet faster ray-triangle intersection (using
0028 %          SSE4)," IEEE Trans. on Visualization and Computer Graphics,
0029 %          16(3):434-438 (2010)
0030 %  [2] Q. Fang, "Comment on 'A study on tetrahedron-based inhomogeneous
0031 %          Monte-Carlo optical simulation'," Biomed. Opt. Express, (in
0032 %          press)
0033 %
0034 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0035 %
0036 
0037 p0=p0(:)';
0038 v0=v0(:)';
0039 
0040 AB=node(face(:,2),1:3)-node(face(:,1),1:3);
0041 AC=node(face(:,3),1:3)-node(face(:,1),1:3);
0042 
0043 N=cross(AB',AC')';
0044 d=-dot(N',node(face(:,1),1:3)')';
0045 
0046 Rn2=1./sum((N.*N)')';
0047 
0048 N1=cross(AC',N')'.*repmat(Rn2,1,3);
0049 d1=-dot(N1',node(face(:,1),1:3)')';
0050 
0051 N2=cross(N',AB')'.*repmat(Rn2,1,3);
0052 d2=-dot(N2',node(face(:,1),1:3)')';
0053 
0054 den=(v0*N')';
0055 t=-(d+(p0*N')');
0056 P=(p0'*den'+v0'*t')';
0057 u=dot(P',N1')'+den.*d1;
0058 v=dot(P',N2')'+den.*d2;
0059 
0060 idx=find(den);
0061 den(idx)=1./den(idx);
0062 
0063 t=t.*den;
0064 u=u.*den;
0065 v=v.*den;
0066 
0067  % if den==0, ray is parallel to triangle, set t to infinity
0068 t(find(den==0))=Inf;
0069 
0070 if(nargout>=4)
0071     idx=find(u>=0 & v>=0 & u+v<=1.0 & ~isinf(t));
0072 end