Point inside triangle(s)?

Hi all, here's a little puzzle that I've confused myself over for the last hour... it's time to clear my head.

I have 4 triangles and one query point. I want to determine which of these triangles contain my query point. A secondary goal will be to identify which triangles have an edge that crosses the query point. This is basically inpolygon(), but for multiple-geometries-per-point, rather than multiple-points-per-geometry (and I'm sure I can make something more efficient than inpolygon since I'm only dealing with triangles).

I thought I'd go about it this way:

• Each triangle has 3 edges.
• Work around the triangles (starting at edge A), get the equation of the line (ie, y=mx+b) that makes up that edge.
• Plug my query point into that equation (ie, result = mx+b-y)
• If result is positive, my query point is below that edge.
• If result is negative, my query point is above that edge.
• If result is zero, my query point is coincident with that edge
• If my query point is above one line but below the other two, it should be inside the triangle

All this is neatly encapsulated in the (commented... running) code below:

    qPtXY = [-0.9414   -0.2134];
    candVerts = [
        -0.9360   -0.2325
        -0.9511   -0.1625
        -0.9511   -0.1681
        -0.9612   -0.1918
        -0.9511   -0.1625
        -0.9239   -0.2012
        -0.9444   -0.2616
        -0.9450   -0.2000
        -0.9239   -0.2012
        -0.9444   -0.2616
        -0.9612   -0.1918
        -0.9511   -0.1625];

    vertxsA = candVerts(1:3:end,:); vertxsB = candVerts(2:3:end,:); vertxsC = candVerts(3:3:end,:);
    % Work out the equation y = mx + b, for the lines from vert1 to vert2
    % First, get m as ydiff/xdiff
    edgeVecs = vertxsB - vertxsA;
    m = edgeVecs(:,2)./edgeVecs(:,1);
    % Next, get b as b=y-mx
    b = vertxsA(:,2) - vertxsA(:,1).*m;
    % Now get the "other side" of 0 = mx + b - y, using the query pt's xy
    testA = m*qPtXY(1) + b -qPtXY(2);

    % Do the same for the second edgeB
    edgeVecs = vertxsC - vertxsB;
    m = edgeVecs(:,2)./edgeVecs(:,1);
    b = vertxsB(:,2) - vertxsB(:,1).*m;
    testB = m*qPtXY(1) + b -qPtXY(2);
    % Do the same for edgeC
    edgeVecs = vertxsA - vertxsC;
    m = edgeVecs(:,2)./edgeVecs(:,1);
    b = vertxsC(:,2) - vertxsC(:,1).*m;
    testC = m*qPtXY(1) + b -qPtXY(2);

    % If the qPt is "above" one line and "below" two lines, it should be
    % inside the triangle
    testSet = [testA testB testC];
    inMask = sum(testSet>0,2)==1 & sum(testSet<0,2)==2;

And the result can be visualised with the following code. Bold lines if the triangle encloses the query point, thin lines if it does not:

    figure, plot(qPtXY(1),qPtXY(2),'r*'), hold on
    triaCols = lines(length(testC));
    lineWidths = 4*inMask + 0.5;
    for ptNo = 1:length(inMask)
        pts3 = [vertxsA(ptNo,:); vertxsB(ptNo,:); vertxsC(ptNo,:)];
        plot(pts3([1:end 1],1), pts3([1:end 1],2), '.b-','Color',triaCols(ptNo,:),'LineWidth',lineWidths(ptNo))
        waitforbuttonpress
    end

It seemed to work, up until the last triangle which obviously doesn't enclose the query point but still meets my criteria for doing so. Obviously my criteria is wrong (or there's a typo somewhere), but I've just gotten myself in a muddle about it. Can someone pick out the errors in my assumptions?

I would also be open to alternative approaches (such as using dot products rather than line equations... that may overcome the issue my method has with vertical edges).

For context, I want to do things in an efficient way so bonus votes for keeping things vectorised. This code will run for between 2 to 20 triangles each time, but looped many many times.

Thanks, Sven.

[per isakson]

Did you google "point inside triangle"?

Answer 1

[the cyclist]

Would the inpolygon() function help?

[Sven]

Ha, sorry cyclist... inpolygon() is set up well for one geometry and many query points... however I am pretty sure it will be hugely inefficient for my "one query point tested on 20 triangles, looped many thousands of times".

Thanks though :)

Answer 2

[Image Analyst]

So you mean like it might take 2 seconds instead of 2 milliseconds? Have you actually tried inpolygon to see how long it would take? I did it (used inpolygon) to check if a point was inside of 10 thousand triangles and it took only 0.55 seconds. What kind of speed must you have to declare success? Copy, paste, and run the following demo.

    clc;    % Clear the command window.
    close all;  % Close all figures (except those of imtool.)
    clear;  % Erase all existing variables.
    workspace;  % Make sure the workspace panel is showing.
    format longg;
    format compact;
    fontSize = 20;

    qPtXY = [0.7 0.8];
    numberOfTriangles = 10000;
    trianglesX = rand(numberOfTriangles, 3);
    trianglesY = rand(numberOfTriangles, 3);

    % START OF ALGORITHM
    tic;
    itsInside = false(numberOfTriangles, 1);
    for t = 1 : numberOfTriangles
      oneTriangleX = [trianglesX(t, 1), trianglesX(t, 2), trianglesX(t, 3), trianglesX(t, 1)];
      oneTriangleY = [trianglesY(t, 1), trianglesY(t, 2), trianglesY(t, 3), trianglesY(t, 1)];
      if inpolygon(qPtXY(1), qPtXY(2), oneTriangleX, oneTriangleY)
        itsInside(t) = true;
      end
    end
    toc
    % ALL DONE!

    % Now plot them all
    % First set up the plots.
    figure;
    subplot(1, 2, 1);
    plot(qPtXY(1), qPtXY(2), 'r*');
    hold on;
    grid on;
    axis on;
    xlim([0 1]);
    ylim([0 1]);
    title('Enclosing Triangles', 'FontSize', fontSize);
    subplot(1, 2, 2);
    plot(qPtXY(1), qPtXY(2), 'r*');
    hold on;
    grid on;
    axis on;
    xlim([0 1]);
    ylim([0 1]);
    title('Non-Enclosing Triangles', 'FontSize', fontSize);
    % Enlarge figure to full screen.
    set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
    % Give a name to the title bar.
    set(gcf,'name','Demo by ImageAnalyst','numbertitle','off') 
    drawnow;

    % Now draw individual triangles, in groups of 100.
    inCount = 0;
    outCount = 0;
    for t = 1 : numberOfTriangles
      oneTriangleX = [trianglesX(t, 1), trianglesX(t, 2), trianglesX(t, 3), trianglesX(t, 1)];
      oneTriangleY = [trianglesY(t, 1), trianglesY(t, 2), trianglesY(t, 3), trianglesY(t, 1)];
      if itsInside(t)
        subplot(1, 2, 1);
        plot(oneTriangleX, oneTriangleY, 'g-');
        if mod(t, 100)
          inCount = inCount + 1;
          inPlotTitle = sprintf('Showing %d enclosing triangles', inCount);
        end
      else
        subplot(1, 2, 2);
        plot(oneTriangleX, oneTriangleY, 'b-');
        if mod(t, 100)
          outCount = outCount + 1;
          outPlotTitle = sprintf('Showing %d non-enclosing triangles', outCount);
        end
        outCount = outCount + 1;
      end
      if mod(t, 100) == 0
        subplot(1, 2, 1);
        title(inPlotTitle, 'FontSize', fontSize);
        subplot(1, 2, 2);
        title(outPlotTitle, 'FontSize', fontSize);
        promptMessage = sprintf('Do you want to Continue displaying triangles,\nor Quit?');
        titleBarCaption = 'Continue?';
        button = questdlg(promptMessage, titleBarCaption, 'Continue', 'Quit', 'Continue');
        if strcmp(button, 'Quit')
          return;
        end
      end
    end

Answer 3

[Anton Semechko]

Hi Sven,

I have encountered this problem many times in the past. I am aware of at least two approaches to solve it. The first one is outlined in this paper: Fast, Minimum Storage Ray-Triangle Intersection (Moller & Thumbore, 1997). If the point lies in the same plane as the triangle, the problem can be reduced to a 2D case. Let {x_ji}, i=1,2,3 be the vertices of the triangle T_j, such that x_ji is an element of R2 space. Using barycentric coordinates an arbitrary point inside the triangle can be defined as: x=x_j1+u*(x_j2-x_j1)+v*(x_j3-x_j1), where 0>=u,v>=1 and u+v<=1. So to check if some point p is inside T_j, solve x-x_j1=[(x_j2-x_j1)+(x_j3-x_j1)][u;v], for u and v and check that 0>=u,v>=1 and u+v<=1.

The second approach is a little different. Suppose that d_12, d_23, and d_31 are the direction vectors perpendicular to the triangle edges, lying in the same plane as the triangle. For example, let d_12 be perpendicular to the edge passing trough vertices x_1 and x_2, and let the remaining two direction vectors be defined similarly. Then to test if an arbitrary point p is inside the triangle, compute three dot products:

dot(p-x_1,d_12)

dot(p-x_2,d_23)

dot(p-x_3,d_31)

If all three are negative then the point is inside the triangle.

I haven't read any of the replies above. So I apologize if the is any redundancy is my response.

Hope this helps.

[Sven]

Thanks Anton. I went with the second (dot product) method (see my answer above). The first paper/implementation looks interesting, but I think it's particular strength lies in the fact it works equally well for aribitrarily oriented 3D rays. In what I'm doing (here), I'm considering only vertical rays so my problem reduces nicely to 2D where (I think) the dot product method is pretty efficient. I know at least that it's easier to grasp :)

You get a vote for this though ;)

Answer 4

[Azzi Abdelmalek]

instead "waitforbuttonpress" use "pause"; i executed your code, it works with "pause"

[Sven]

Azzi, the waitforbuttonpress() is purely so that you can see each triangle one at a time (before they're thrown onto the screen together). This displaying of the result is not the code itself, just a visualisation. Remove waitforbuttonpress() entirely if you'd like.

The code "runs" just fine, but "works" is an entirely different concept... and it doesn't do that. This is demonstrated by the last of the 4 triangles being detected as encompassing the query point when it clearly does not.

Answer 5

[Sven]

So here's an answer (that works!) using dot products.

Imagine one vertex of a triangle. The two edges stemming from that vertex make two vectors (a "V" shape). Normalise these two vectors, and take their dot product. If the dot product is near to 1, the "V" shape is very acute, if the dot product is near -1, the "V" shape is very obtuse.

Now make another unit vector from the shared vertex (the apex of the "V" shape) to the query point, and calculate the two dot products between this unit vector and the two edge unit vectors.

If the two dot products using the query point are larger than the dot product from the edges of the "V" shape, the query point must lie "inside" the V.

If this is the case for all three vertexes on the triangle, the point is inside the triangle. Code below:

    % Get my N triangles' vertices into a 3-by-2-by-N shape (easier to work with)
    myVerts = permute(reshape(candVerts',2,3,[]),[2 1 3]);
    % Get unit vectors pointing from each triangle vertex to my query point
    vert2ptVecs = bsxfun(@minus, qPtXY, myVerts);
    vert2ptUVecs = bsxfun(@rdivide, vert2ptVecs, sqrt(sum(vert2ptVecs.^2,2)));
    % Get unit vectors pointing around each triangle (along edge A, edge B, edge C)
    edgeVecs = myVerts([2 3 1],:,:) - myVerts;
    edgeUVecs = bsxfun(@rdivide, edgeVecs, sqrt(sum(edgeVecs.^2,2)));
    % Get the inner product between edgeA.edgeC, edgeB.edgeA, edgeC.edgeB
    edgeEdgeDotPs = sum(edgeUVecs .* -edgeUVecs([3 1 2],:,:),2);
    % Get the inner product between each edge unit vec and the unit vect from qPt to vertex
    edgeQPntDotPs = sum(edgeUVecs .* vert2ptUVecs,2);
    qPntEdgeDotPs = sum(vert2ptUVecs .* -edgeUVecs([3 1 2],:,:),2);
    % If both inner products 2 edges to the query point are greater than the inner product between
    % the two edges themselves, the query point is between the V shape made by the two edges. If
    % this is true for all 3 edge pair, the query point is inside the triangle.
    result = squeeze(all(edgeQPntDotPs>edgeEdgeDotPs & qPntEdgeDotPs>edgeEdgeDotPs,1));

Elapsed time is 0.002080 seconds.

Now for comparison, using inpolygon():

    nTrias = size(candVerts,1)/3;
    result = false(nTrias,1);
    for tNo = 1:nTrias
        thisTria = candVerts((1:3)+(tNo-1)*3,:);
        result(tNo) = inpolygon(qPtXY(1),qPtXY(2), thisTria(:,1), thisTria(:,2));
    end

For speed comparison, my example above (looped 10000 times) takes 1.267937 seconds. The inpolygon() version takes 4.494493 seconds. I'll throw in any other solution times here.

[Image Analyst]

4.5 seconds for 10,000 triangles? My code is virtually identical and it took one tenth that time. But whatever, as long as you have something you're happy with (though it seems a lot harder to follow/understand than the simple inpolygon approach).

[Sven]

Thanks IA.

Here's a comparison using the same points setup:

numberOfTriangles = 10000;
qPtXY = [0.7 0.8];
myVerts = rand(3,2,numberOfTriangles);

% START OF ALGORITHM 1
tic
vert2ptVecs = bsxfun(@minus, qPtXY, myVerts);
vert2ptUVecs = bsxfun(@rdivide, vert2ptVecs, sqrt(sum(vert2ptVecs.^2,2)));
edgeVecs = myVerts([2 3 1],:,:) - myVerts;
edgeUVecs = bsxfun(@rdivide, edgeVecs, sqrt(sum(edgeVecs.^2,2)));
edgeEdgeDotPs = sum(edgeUVecs .* -edgeUVecs([3 1 2],:,:),2);
edgeQPntDotPs = sum(edgeUVecs .* vert2ptUVecs,2);
qPntEdgeDotPs = sum(vert2ptUVecs .* -edgeUVecs([3 1 2],:,:),2);
result = squeeze(all(edgeQPntDotPs>edgeEdgeDotPs & qPntEdgeDotPs>edgeEdgeDotPs,1));
toc

% START OF ALGORITHM 2
tic;
itsInside = false(numberOfTriangles, 1);
for t = 1 : numberOfTriangles
    itsInside(t) = inpolygon(qPtXY(1), qPtXY(2), myVerts(:,1,t), myVerts(:,2,t));
end
toc

Elapsed time is 0.010236 seconds. Elapsed time is 0.821573 seconds.

So I get an 80x difference. This is of course the "best case scenario" for the vectorised version over the looped inpolygon() version - here it's one run of 10000 triangles rather than 10000 runs of 4-to-20 triangles as my original description stated (which showed a 3.5x speedup). But with some pre-computation of triangle edge unit vectors, I think I can get closer to the 80x difference than the 3.5x difference, which starts to get noticeable. You get a vote though for the very, very, very thorough graphical display ;)

My answer doesn't yet handles ON triangles instead of IN triangles (so there's a minor difference to inpolygon()), but that should be a reasonably simple addition.

For reference, I'm trying to write a 3d version of inpolygon()... inpolyhedral()... which I'll put on the FEX when done. There have been others put up there (John D'Errico's inhull is great but only for convex hulls), but I'm trying to optimise speed as we often work with big meshes and previous versions have us sitting for ~1minute when I think it can be seconds.

[Sean de Wolski]

@Sven, vert2con on the FEX also works well for (2:n)dimensional polygons, but once again they must be convex.

[Image Analyst]

You won't be able to handle ON an edge of a triangle since the numbers are quantized and it's virtually impossible for two floating point numbers to match (I assume you've read the FAQ, right: http://matlab.wikia.com/wiki/FAQ#Why_is_0.3_-_0.2_-_0.1_.28or_similar.29_not_equal_to_zero.3F). But you can determine if it's within some tolerance distance of the edge.

[Sven]

Yup, got that. Tolerance has to be used to get ON... that's why it's not in my first pass solution :)

[Matt Kindig]

Hi Sven,

I have used an entry on the File Exchange named 'inpoly' as a replacement for inpolygon. It is MUCH faster than inpolygon in my experience; also it is able to detect whether a point is on the edge. I would image that looping through all triangles with this function would be fairly fast.

Matt

[Sven]

Hi all, thanks for the tips.

The FEX entry I was writing is now up here: http://www.mathworks.com/matlabcentral/fileexchange/37856-inpolyhedron-are-points-inside-a-volume.

It uses the ray-casting algorithm along the Z-dimension, which is where the point-in-polygon part of the algorithm is used to test if a point inside the projection of a surface onto the XY plane.