function [isinside,pt,coord]=linextriangle(p0,p1,plane)
%  [isinside,pt,coord]=linextriangle(p0,p1,plane)
%
%  calculate the intersection of a 3d line (passing two points)
%  with a plane (determined by 3 points)
%
%  author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
%  date: 12/12/2008
%
% parameters: 
%      p0: a 3d point in form of (x,y,z)
%      p1: another 3d point in form of (x,y,z), p0 and p1 determins the line
%      plane: a 3x3 matrix, each row is a 3d point in form of (x,y,z)
%             this is used to define a plane
% outputs:
%      isinside: a boolean variable, 1 for the intersection is within the 
%               3d triangle determined by the 3 points in plane; 0 is outside
%      pt: the coordinates of the intersection pint
%      coord: 1x3 vector, if isinside=1, coord will record the barycentric 
%          coordinate for the intersection point within the triangle; 
%          otherwise it will be all zeros.
%
% for degenerated lines or triangles, this will stop
%
% Please find more information at http://iso2mesh.sf.net/cgi-bin/index.cgi?metch
%
% this function is part of "metch" toobox, see COPYING for license

[a,b,c,d]=getplanefrom3pt(plane);

if(a*a+b*b+c*c==0.0)
        error('degenerated plane');
end

dl_n=sum([a b c].*(p1-p0));

if(dl_n==0.0)
    error('degenerated line');
end

% solve for the intersection point
t=-(a*p0(1)+b*p0(2)+c*p0(3)+d)/dl_n;
pt=p0+(p1-p0)*t;


dist=sum(abs(diff(plane)));
[md,imax]=sort(dist);
if(md(2)==0.0)
    error('degenerated triangle');
end
goodidx=imax(2:end);

ptproj=pt(goodidx);
mat0=[plane(:,goodidx)',ptproj';1 1 1 1];

isinside=0;
coord=[0 0 0];

det1=det(mat0(:,[4 2 3],:));
det2=det(mat0(:,[1 4 3],:));
if(det1*det2<0) 
    return; 
end
det3=det(mat0(:,[1 2 4],:));
if(det2*det3<0) 
    return;
end
if(det1*det3<0) 
    return;
end
isinside=1;
det0=det(mat0(:,1:3));

coord=[det1 det2 det3]/det0;
