ray2surf.m

function [p,e0]=ray2surf(node,elem,p0,v0,e0)
%
% [p,e0]=ray2surf(node,elem,p0,v0,e0)
%
% Determine the entry position and element for a ray to intersect a mesh
%
% author: Qianqian Fang (q.fang <at> neu.edu)
%
% input:
%     node: the mesh coordinate list
%     elem: the tetrahedral mesh element list, 4 columns
%     p0: origin of the ray
%     v0: direction vector of the ray
%     e0: search direction: '>' forward search, '<' backward search, '-' bidirectional
%
% output:
%     p: the intersection position
%     e0: if found, the index of the intersecting element ID
%
% this file is part of Mesh-based Monte Carlo (MMC)
%
% License: GPLv3, see http://mcx.sf.net/mmc/ for details
%

p=p0;
if(size(elem,2)==3)
    face=elem;
else
    face=volface(elem);
end
[t,u,v,idx]=raytrace(p0,v0,node,face);
if(isempty(idx))
    error('ray does not intersect with the mesh');
else
    t=t(idx);
    if(e0=='>')
%         idx1=find(t>=0);
        idx1=find(t>=1e-10);
    elseif(e0=='<')
        idx1=find(t<=0);
    elseif(isnan(e0) || e0=='-')
        idx1=1:length(t);
    else
        error('ray direction specifier is not recognized');
    end
    if(isempty(idx1))
        error('no intersection is found along the ray direction');
    end
    t0=abs(t(idx1));
    [tmin,loc]=min(t0);
    faceidx=idx(idx1(loc));

    % update source position
    p=p0+t(idx1(loc))*v0;

    if(nargout<2)
        return;
    end

    % find initial element id
    if(size(elem,2)==3)
        e0=faceidx;
    else
       felem=sort(face(faceidx,:));
       f=elem;
       f=[f(:,[1,2,3]);
          f(:,[2,1,4]);
          f(:,[1,3,4]);
          f(:,[2,4,3])];
       [tf,loc]=ismember(felem,sort(f,2),'rows');
       loc=mod(loc,size(elem,1));
       if(loc==0) loc=size(elem,1); end
       e0=loc;
    end
end