How do I use 'intersection' to get the overlapping line between two polyshapes?

Hi - how do I get the line intersection between two polyshapes that overlap along a line? A super-simplified example is below. For the real case, the shapes are more complicated, but I think the simplified case covers it. I'm maybe not understanding how the intersection function works! Thank you!

main()

function main()

%clear all;     % completely superfluous in function workspace
x1 = [1 4 4 1];
y1 = [1 1 4 4];
x2 = [4 6 6 4];
y2 = [2 2 3 3];
poly1 = polyshape(x1,y1);
% the poly2 'grows' from poly1
poly2 = polyshape(x2,y2);
poly2 = union(poly1,poly2);
growthPoly = subtract(poly2,poly1,"KeepCollinearPoints",true);
close all;
hold on;    
plot(poly1,'FaceColor','none','EdgeColor','k','LineWidth',12);
plot(poly2,'FaceColor','none','EdgeColor','g','LineWidth',8);
plot(growthPoly,'FaceColor','none','EdgeColor','r','LineWidth',4);
% this comes back empty and nothing gets plotted
growthLine = intersect(poly1,growthPoly,"KeepCollinearPoints",true);
plot(growthLine,'FaceColor','none','EdgeColor','b','LineWidth',2);
legend({'The initial poly1',...
    'poly2 grew from poly1',...
    'The growth area is the difference between poly2 and poly1'},'Box','off','Location','south','FontSize',18);
title({'How do I get the ''growthLine'' part that is red-on-black?';'I thought it would be the intersection of the black and red, but no dice';'I thought it would show up as a blue line'},'FontSize',22);
end

Answer 1

[Steven Lord]

These two polyshapes touch but do not intersect.

x1 = [1 4 4 1];
y1 = [1 1 4 4];
x2 = [4 6 6 4];
y2 = [2 2 3 3];
poly1 = polyshape(x1,y1);
% the poly2 'grows' from poly1
poly2 = polyshape(x2,y2);
poly2 = union(poly1,poly2);
growthPoly = subtract(poly2,poly1,"KeepCollinearPoints",true);
close all;
hold on;    
plot(poly1,'FaceColor','none','EdgeColor','k','LineWidth',12);
plot(poly2,'FaceColor','none','EdgeColor','g','LineWidth',8);
plot(growthPoly,'FaceColor','none','EdgeColor','r','LineWidth',4);

% this comes back empty and nothing gets plotted
growthLine = intersect(poly1,growthPoly,"KeepCollinearPoints",true)

growthLine = 
  polyshape with properties:

      Vertices: [0×2 double]
    NumRegions: 0
      NumHoles: 0

It looks like they intersect due to the LineWidths, but if we plot them with the default line widths:

figure
plot(poly1,'FaceColor','none','EdgeColor','k');
hold on
plot(growthPoly,'FaceColor','none','EdgeColor','g');

You could add a small polybuffer to the polyshapes as Loren Shure and Matt Tearle did to solve a "Car Talk" puzzle back in 2017.

d = 1/64;
poly1a = polybuffer(poly1, d);
growthPoly2 = polybuffer(growthPoly, d);
intersectRegion = intersect(poly1a, growthPoly2)

intersectRegion = 
  polyshape with properties:

      Vertices: [94×2 double]
    NumRegions: 1
      NumHoles: 0

figure
plot(poly1a, EdgeColor = "g", FaceColor = "none")
hold on
plot(growthPoly2, EdgeColor = "r", FaceColor = "none")
plot(intersectRegion, EdgeColor = "k", FaceColor = "none")

[Robert]

That is a cool idea, thank you. I'll give it a shot.

[dpb]

"...the two polyshapes touch but do not intersect.... but if we plot them with the default line widths"

x1 = [1 4 4 1];
y1 = [1 1 4 4];
x2 = [4 6 6 4];
y2 = [2 2 3 3];
poly1 = polyshape(x1,y1);
% the poly2 'grows' from poly1
poly2 = polyshape(x2,y2);
poly2 = union(poly1,poly2);
growthPoly = subtract(poly2,poly1,"KeepCollinearPoints",true);
plot(poly1,'FaceColor','none','EdgeColor','k');
hold on
plot(growthPoly,'FaceColor','none','EdgeColor','r');
xlim([3.99995 4.00005]), ylim([1.5 3.5])

No matter how much you expand the x axis, the two definitely do coincide along the x==4 edges as can be seen by

c=intersect(growthPoly.Vertices(:,1),poly1.Vertices(:,1))

c = 4

ia=(growthPoly.Vertices(:,1)==c);
ib=(poly1.Vertices(:,1)==c);
[growthPoly.Vertices(ia,:) poly1.Vertices(ib,:)]

ans = 2×4

4 3 4 4 4 2 4 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The definition of "intersect" appears to be too limited in my view or there should be another set of functions for the polyshape family that will find the overlaying boundaries without kludging the dimensions. It's possible and I suppose the reason for the chosen implementation is that the general floating point comparison may be an issue such that ismembertol was introduced such that rounding tolerances could be accomodated. That's essentially what they polybuffer does but it would be far cleaner and more to a user's expectations I think if the builtin functions accomodated the case. At least, the documentation should make a much more clear definition of what "intersect" really means and show an example of the workaround with the buffer.

However, the polyshape class doesn't know anything about just a line; would have to just retain the x,y pairs from growthPoly.Vertices above as the defining end points.

growthLine=polyshape(growthPoly.Vertices(ia,:))

Warning: Boundaries with fewer than 3 points were removed.

growthLine = 
  polyshape with properties:

      Vertices: [0×2 double]
    NumRegions: 0
      NumHoles: 0

growthLine=growthPoly.Vertices(ia,:);
hL=plot(growthLine(:,1),growthLine(:,2),'b-','linewidth',2);

[Steven Lord]

If you're sitting in a chair, what's the intersection between you and the chair?

Or to use the "determine which states are adjacent" scenario from the Car Talk puzzle I mentioned, what's the intersection between the states of Massachusetts and Connecticut? Can you point at a piece of soil and say that it is in both states simultaneously? They share an edge, but don't overlap.

The question "What boundaries do these two (or more) polyshapes share?" is a legitimate question. There's currently no object function / method included in MATLAB to answer that question directly. I think that's worth submitting an enhancement request to Support, ideally with your use case, so the developers understand what your requirements are. For example, if the two polyshapes share only a single point, would that count as a "boundary" for purposes of that question?

p1 = polyshape([0 1 1 0], [0 0 1 1]);
p2 = translate(p1, [1 1]);
plot(p1, FaceColor = 'r');
hold on
plot(p2, FaceColor = 'k');

Would this hypothetical sharedBoundary function return the point (1, 1)?

[Robert]

Would this hypothetical sharedBoundary function return the point (1, 1)?

I would say yes, for sure. Maybe as a shared vertex.

I have submitted an enhancement request pointing to this case!

[dpb]

@Steven Lord

A, The person and the chair are two separate objects and while in contact don't actually intersect, agreed...the two shapes might be modeled as initially in the same location, but only by retaining the two shapes could they be modeled to represent shifting or rising from the chair.

B. There are an innumberable number particles of sand/earth that are on the theoretical line dividing the two states that share ownership ... there's no black hole or empty canyon between the two states.

I'll agree that calling it a boundary rather than an intersection may be better terminology, but the toolset as implemented/documented leads to confusion and unexpected results the presented workaround seeming very klunky. I'd not found reason for use of the feature ere this Q? and, like @Robert found the result less than understandable from the documentation and only could observe that the implementation doesn't count overlapping boundaries as intersections and that isn't clearly documented.

@Robert

I would agree on the single point being the expected result as well. If they're coincident, they're still a joining point.

Answer 2

[dpb]

%function main()

x1 = [1 4 4 1];
y1 = [1 1 4 4];
x2 = [4 6 6 4];
y2 = [2 2 3 3];
poly1 = polyshape(x1,y1);
% the poly2 'grows' from poly1
poly2 = polyshape(x2,y2);
poly3 = union(poly1,poly2);
growthPoly = subtract(poly3,poly1,"KeepCollinearPoints",true);
hold on;    
plot(poly1,'FaceColor','none','EdgeColor','k','LineWidth',12);
plot(poly3,'FaceColor','none','EdgeColor','g','LineWidth',8);
plot(growthPoly,'FaceColor','none','EdgeColor','r','LineWidth',4);
% this comes back empty and nothing gets plotted
growthLine = intersect(poly1,growthPoly,"KeepCollinearPoints",true);
plot(growthLine,'FaceColor','none','EdgeColor','b','LineWidth',2);
legend({'The initial poly1',...
    'poly2 grew from poly1',...
    'The growth area is the difference between poly2 and poly1'}, ...
    'Box','off','Location','south','FontSize',10);

%end

How do I get the ''growthLine'' part that is red-on-black?

I thought it would be the intersection of the black and red, but no dice;

I thought it would show up as a blue line

%main
poly1.Vertices

ans = 4×2

1 1 1 4 4 4 4 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%poly2.Vertices
%poly3.Vertices
growthPoly.Vertices

ans = 4×2

4 3 6 3 6 2 4 2

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

overlaps(poly1,growthPoly)

ans = logical
   0

I'm not sure about the use of inbuilt functions for the polyshape objects, but it appears that intersect only considers one within the boundaries of the other and doesn't count the edges as actually intersecting.

I dunno if that would be considered a bug or is working as designed, but looks like a support ticket item to me. I agree it seems the line section should be found. The only example with a line is one drawn diagonally across the edge of a shape and it finds the in/out sections; doesn't address the crossing point itself.

[Robert]

That makes sense, thank you. I'll write a support ticket to see if they can add the option to consider edges as intersecting!

[Matt J]

No, it is not an implementable support request to have edges included in the intersection. Mathematically, it is true that two polygons can intersect at their edges in a line, but it is not a numerically stable problem to solve for such an intersection.

In your example, you have simple integer-valued vertices, but in general the vertex data will be floating point values with finite precision. The precision uncertainty makes it impossible to determine reproducibly whether 2 edges intersect.

[dpb]

I don't see it as being impossible, just have to explicitly set a fuzziness for what is considered a match. One would have to do something similar to what ismembertol does for the purpose in general arrays that aren't polyshapes. That's sorta' what the workaround @Steven Lord showed does with the fudging of the boundaries.

As noted, I'd allow as how the "intersection" maybe shouldn't include boundary-matching but I definitely think the existing toolset is lacking in not having the facility to find them.