利用http://blog.sina.com.cn/s/blog_6163bdeb0102dvcs.html提供的思想，用matlab编写了程序如下

% 判断点与三角形的位置关系
function flag = PointInTriangle(P, Triangle, method)
if nargin == 2
    method = 's';
end

flag = zeros(size(P, 1), 1);
switch method
    case 'a'
        for i = 1:size(P, 1)
            flag(i) = AreaMth(P(i, :), Triangle);
        end
    case 's'
        for i = 1:size(P, 1)
            flag(i) = SyntropyMth(P(i, :), Triangle);
        end
    otherwise
        error('没有这种方法！');
end
end

% 同向法
function flag = SyntropyMth(P, Triangle)
P(3) = 0;
Triangle(:,3) = 0;
A = Triangle(1, :);
B = Triangle(2, :);
C = Triangle(3, :);

Pa = SameSide(P,A,B,C);
Pb = SameSide(P,B,A,C);
Pc = SameSide(P,C,A,B);

if Pa>0 && Pb>0 && Pc>0
    flag = 1;
elseif Pa*Pb*Pc == 0
    flag = 0;
else
    flag = -1;
end

% 两点在线的同一侧
    function flag = SameSide(p1, p2, A, B)
        cp1 = cross(B-A, p1-A);
        cp2 = cross(B-A, p2-A);
        flag = sign(dot(cp1, cp2)); % 1-同侧，0-线上，-1-异侧
    end
end

% 重心法或面积法
function flag = AreaMth(P, Triangle)
A = Triangle(1, :);
B = Triangle(2, :);
C = Triangle(3, :);

% Compute vectors
v0 = C - A;
v1 = B - A;
v2 = P - A;

% Compute dot products
dot00 = dot(v0, v0);
dot01 = dot(v0, v1);
dot02 = dot(v0, v2);
dot11 = dot(v1, v1);
dot12 = dot(v1, v2);

% Compute barycentric coordinates
invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
u = (dot11 * dot02 - dot01 * dot12) * invDenom;
v = (dot00 * dot12 - dot01 * dot02) * invDenom;

% Check if point is in triangle
if u > 0 && v > 0 && u + v < 1
    flag = 1;
elseif u == 0 || v == 0 || u + v == 1
    flag = 0;
else
    flag = -1;
end
end

测试如下

P = [-1 2;
    0 1;
    0.5 0.5];
Triangle = [0 0;
    0 2;
    2 0];

flag = PointInTriangle(P, Triangle, 'a')

结果为

flag =
    -1
     0
     1