function [t,u,v,idx,xnode]=raysurf(p0,v0,node,face)
%
% [t,u,v,idx,xnode]=raysurf(p,v,node,face)
%
% perform a Havel-styled ray tracing for a triangular surface
%
% author: Qianqian Fang, <fangq at nmr.mgh.harvard.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)
%

len=size(p0,1);
if(len==0)
   error('p0 can not be empty');
end
if(size(node,2)<3)
   error('node must contain at least 3 columns');
end
if(size(face,2)<3)
   error('face must contain at least 3 columns');
end

if(size(v0,1)==1 | size(v0,2)==1 & len>1)
   v0=repmat(v0(:)',len,1);
end

t=zeros(len,1)*nan;
u=t;
v=t;
idx=t;

for i=1:len
   [ti,ui,vi,id]=raytrace(p0(i,:),v0(i,:),node,face);
   if(isempty(id)) continue; end
   ti=ti(id);
   tpid=find(ti>=0);
   if(isempty(tpid)) continue; end
   [tmin,tloc]=min(ti(find(ti>=0)));
   t(i)=tmin;
   u(i)=ui(id(tpid(tloc)));
   v(i)=vi(id(tpid(tloc)));
   idx(i)=id(tpid(tloc));
end

if(nargout>=5)
   xnode=p0+repmat(t,1,3).*v0;
end

