linextriangle.m

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 <q.fang at neu.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;