raysurf

 [t,u,v,idx,xnode]=raysurf(p,v,node,face)

 perform a Havel-styled ray tracing for a triangular surface

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

 input:
   p0: list of starting points of the rays
   v0: directional vector of the rays, 
   node: a list of node coordinates (nn x 3)
   face: a surface mesh triangle list (ne x 3)

 output:
   t: distance from p0 to the intersection point for each surface
      triangle, if t(i)=NaN, no intersection was found for that ray
   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: idx lists the IDs of the face elements that intersects 
      each ray
   xnode: optional output, if requested, xnode gives the intersection
      point coordinates; to compute manually, xnode=p0+repmat(t,1,3).*v0

 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,xnode]=raysurf(p0,v0,node,face)
0002 %
0003 % [t,u,v,idx,xnode]=raysurf(p,v,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: list of starting points of the rays
0011 %   v0: directional vector of the rays,
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: distance from p0 to the intersection point for each surface
0017 %      triangle, if t(i)=NaN, no intersection was found for that ray
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: idx lists the IDs of the face elements that intersects
0022 %      each ray
0023 %   xnode: optional output, if requested, xnode gives the intersection
0024 %      point coordinates; to compute manually, xnode=p0+repmat(t,1,3).*v0
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 len=size(p0,1);
0038 if(len==0)
0039    error('p0 can not be empty');
0040 end
0041 if(size(node,2)<3)
0042    error('node must contain at least 3 columns');
0043 end
0044 if(size(face,2)<3)
0045    error('face must contain at least 3 columns');
0046 end
0047 
0048 if(size(v0,1)==1 | size(v0,2)==1 & len>1)
0049    v0=repmat(v0(:)',len,1);
0050 end
0051 
0052 t=zeros(len,1)*nan;
0053 u=t;
0054 v=t;
0055 idx=t;
0056 
0057 for i=1:len
0058    [ti,ui,vi,id]=raytrace(p0(i,:),v0(i,:),node,face);
0059    if(isempty(id)) continue; end
0060    ti=ti(id);
0061    tpid=find(ti>=0);
0062    if(isempty(tpid)) continue; end
0063    [tmin,tloc]=min(ti(find(ti>=0)));
0064    t(i)=tmin;
0065    u(i)=ui(id(tpid(tloc)));
0066    v(i)=vi(id(tpid(tloc)));
0067    idx(i)=id(tpid(tloc));
0068 end
0069 
0070 if(nargout>=5)
0071    xnode=p0+repmat(t,1,3).*v0;
0072 end
0073