New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

geometry : Support for implicit intersecting polygons #12931

Closed
wants to merge 3 commits into
base: master
from

Conversation

Projects
None yet
4 participants
@ArifAhmed1995
Contributor

ArifAhmed1995 commented Jul 10, 2017

Intersecting polygons with both implicit and explicit intersections are allowed. Remaining XFAILs for intpoly module pass as a result.

@certik

This comment has been minimized.

Show comment
Hide comment
@certik

certik Jul 10, 2017

Member

@ArifAhmed1995 nice! I'll review it soon. It looks like you are able to reproduce results from @ndotsu's paper.

Member

certik commented Jul 10, 2017

@ArifAhmed1995 nice! I'll review it soon. It looks like you are able to reproduce results from @ndotsu's paper.

@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 10, 2017

Member

Does this do the right thing when there are multiple intersections, e.g. a W with a strikethrough made by dropping from the top right down half height and then going left the width of the W at the top.

>>> assert len(Polygon((0,2),(1,0),(2,2),(3,0),(4,2),(4,1),(0,1)).sides) == 15

Regarding directed_sides...I have to think about that a bit. I would prefer another method, either Polygon.directed_sides (a method, not a function that allows you to pass True/False) or Segment.canonical.

📝 why are directed sides needed?

Member

smichr commented Jul 10, 2017

Does this do the right thing when there are multiple intersections, e.g. a W with a strikethrough made by dropping from the top right down half height and then going left the width of the W at the top.

>>> assert len(Polygon((0,2),(1,0),(2,2),(3,0),(4,2),(4,1),(0,1)).sides) == 15

Regarding directed_sides...I have to think about that a bit. I would prefer another method, either Polygon.directed_sides (a method, not a function that allows you to pass True/False) or Segment.canonical.

📝 why are directed sides needed?

ArifAhmed1995 added some commits Jul 10, 2017

geometry : Support for implicit intersecting polygons
There is an optional argument to ignore the implicit intersections
which is required for the algorithm in `intpoly` module to work correctly.
@ArifAhmed1995

This comment has been minimized.

Show comment
Hide comment
@ArifAhmed1995

ArifAhmed1995 Jul 11, 2017

Contributor

@smichr I've added the test you mentioned and changed directed_sides to a getter method.
Direction is important when a line segment is decomposed into it's constituents. Because the traversal of points would be incorrect otherwise.

Contributor

ArifAhmed1995 commented Jul 11, 2017

@smichr I've added the test you mentioned and changed directed_sides to a getter method.
Direction is important when a line segment is decomposed into it's constituents. Because the traversal of points would be incorrect otherwise.

Show outdated Hide outdated sympy/geometry/polygon.py
sides = {}
for i, side in enumerate(rv.directed_sides):
sides[i] = [side]
for i in sides:

This comment has been minimized.

@smichr
@smichr

smichr Jul 11, 2017

Member
Show outdated Hide outdated sympy/geometry/polygon.py
consider a polygon shaped like a Z with the top left connecting to the
bottom right of the Z. The point in the middle of the Z may or may not be
explicitly given. However it should be noted that not all characteristics
of the implicit case be accurate(such as area):

This comment has been minimized.

@smichr

smichr Jul 12, 2017

Member

Suggested rewrite to the effect of:

When a segment connecting two of the provided vertices crosses
one or more other segments of the polygon, the points of intersection
will be automatically inserted into the arguments. Note that two
polygons may have the same shape but have different area since area
depends on the *path* taken to create the polygon. For example, 
creating a Z with a line connecting the bottom right to the top left
gives an overall shape of two triangles meeting at common vertex. 
Defining/drawing the shape this way, however, is consistent with
moving in a counter-clockwise direction around the lower triangle and 
a clockwise direction around the top triangle so their areas cancel:

>>> mid, tl, tr, br, bl = (0,0), (-1,1), (1,1), (1,-1), (-1, -1)
>>> cwccw = Polygon(mid, tl, tr, mid, bl, br, mid)
>>> cwccw.area
0

If the same shape is made by defining a clockwise path around
both triangles, the areas do not cancel:

>>> cw = Polygon(mid, tl, tr, mid, br, bl, mid)
>>> cw.area
-2

When allowing the intersections to be computed automatically,
the path is consistent with the path specified by the vertices:

>>> Polygon(tl, tr, bl, br, tl).area
0

To disallow intersections...[blah blah]...then an error will be raised.  <-- is that right?

The caution sounds a little scary in saying that the area may not be accurate.

@smichr

smichr Jul 12, 2017

Member

Suggested rewrite to the effect of:

When a segment connecting two of the provided vertices crosses
one or more other segments of the polygon, the points of intersection
will be automatically inserted into the arguments. Note that two
polygons may have the same shape but have different area since area
depends on the *path* taken to create the polygon. For example, 
creating a Z with a line connecting the bottom right to the top left
gives an overall shape of two triangles meeting at common vertex. 
Defining/drawing the shape this way, however, is consistent with
moving in a counter-clockwise direction around the lower triangle and 
a clockwise direction around the top triangle so their areas cancel:

>>> mid, tl, tr, br, bl = (0,0), (-1,1), (1,1), (1,-1), (-1, -1)
>>> cwccw = Polygon(mid, tl, tr, mid, bl, br, mid)
>>> cwccw.area
0

If the same shape is made by defining a clockwise path around
both triangles, the areas do not cancel:

>>> cw = Polygon(mid, tl, tr, mid, br, bl, mid)
>>> cw.area
-2

When allowing the intersections to be computed automatically,
the path is consistent with the path specified by the vertices:

>>> Polygon(tl, tr, bl, br, tl).area
0

To disallow intersections...[blah blah]...then an error will be raised.  <-- is that right?

The caution sounds a little scary in saying that the area may not be accurate.

@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 12, 2017

Member

OK, I've looked at this more closely. I largely think the modifications are not necessary. To make the polytope integration work we just have to not raise an error and then not compute the intersections so that's easy. When the interersections are computed, they give the same result (for area) as if they were not computed. So I don't see the value added for computing the intersections. The only thing I could see the intersections being useful for is if someone wanted to draw a non-overlapping version of the polygon by using overlapping segments which are then "deconstructed" into the desired polygon. e.g. the vertical bowtie (Z with connected endpoints).

Though the name would be longer, I would prefer 'ignore_intersections' to be used instead of 'ignore_int' since I always think of python integers when I see 'int'.

Questions

Is there a name for the polygon when it is intended to be interpreted sharing vertices between portions, not crossing sides?
What properties besides area depend on this behavior?
If only area depends on it, we might just include a property "footprint_area" to compute the area covered by the polygon rather than the more abstract area that depends on the direction in which the sides are drawn. And maybe there is an easier way to compute that than computing all the intersections. My suspicion is that there is not. I would love to be wrong, though!

Ideas

Maybe your code could be added to a Polygon.hull property which would return the non-crossing outline of a polygon. A convex polygon would return the convex hull; an intersecting one would return a ccw traversed polygon (whose area would then be positive) whose segments never cross each other but only meet at vertices (which your code currently detects).

Member

smichr commented Jul 12, 2017

OK, I've looked at this more closely. I largely think the modifications are not necessary. To make the polytope integration work we just have to not raise an error and then not compute the intersections so that's easy. When the interersections are computed, they give the same result (for area) as if they were not computed. So I don't see the value added for computing the intersections. The only thing I could see the intersections being useful for is if someone wanted to draw a non-overlapping version of the polygon by using overlapping segments which are then "deconstructed" into the desired polygon. e.g. the vertical bowtie (Z with connected endpoints).

Though the name would be longer, I would prefer 'ignore_intersections' to be used instead of 'ignore_int' since I always think of python integers when I see 'int'.

Questions

Is there a name for the polygon when it is intended to be interpreted sharing vertices between portions, not crossing sides?
What properties besides area depend on this behavior?
If only area depends on it, we might just include a property "footprint_area" to compute the area covered by the polygon rather than the more abstract area that depends on the direction in which the sides are drawn. And maybe there is an easier way to compute that than computing all the intersections. My suspicion is that there is not. I would love to be wrong, though!

Ideas

Maybe your code could be added to a Polygon.hull property which would return the non-crossing outline of a polygon. A convex polygon would return the convex hull; an intersecting one would return a ccw traversed polygon (whose area would then be positive) whose segments never cross each other but only meet at vertices (which your code currently detects).

@ArifAhmed1995

This comment has been minimized.

Show comment
Hide comment
@ArifAhmed1995

ArifAhmed1995 Jul 12, 2017

Contributor

@smichr

When the interersections are computed, they give the same result (for area) as if they were not computed.

Unfortunately they do not :

>>> polytope_integrate(Polygon((0, 2), (2, 2), (0, 0), (2, 0)), 1)
0
>>> polytope_integrate(Polygon((0, 2), (2, 2), (1, 1), (2, 0), (0, 0), (1, 1)), 1)
2

This was on master with only the error throwing block commented out.
This discrepancy in calculating area or other integrals on such polytopes requires a traversal which avoids decomposing segments. This can be achieved by first computing the intersections by the BO algorithm and the correct traversal by Fleury's algorithm with an extra constraint.
The extra constraint is that during traversal travelling between three collinear points is disallowed.
Due to lack of better name as of now, I'll refer to this as : outer traversal, since one does not cross any other segment.

Is there a name for the polygon when it is intended to be interpreted sharing vertices between portions, not crossing sides?

Would the term complex polygon suffice ?

What properties besides area depend on this behavior?

Any integral computed on such a polytope without outer traversal would be incorrect.

So, I think computing intersection using BO algorithm and then traversal using Fleury's algorithm would be the correct solution. Also, the default value of ignore_intersections should be changed to True.
Waiting for your thoughts.

Contributor

ArifAhmed1995 commented Jul 12, 2017

@smichr

When the interersections are computed, they give the same result (for area) as if they were not computed.

Unfortunately they do not :

>>> polytope_integrate(Polygon((0, 2), (2, 2), (0, 0), (2, 0)), 1)
0
>>> polytope_integrate(Polygon((0, 2), (2, 2), (1, 1), (2, 0), (0, 0), (1, 1)), 1)
2

This was on master with only the error throwing block commented out.
This discrepancy in calculating area or other integrals on such polytopes requires a traversal which avoids decomposing segments. This can be achieved by first computing the intersections by the BO algorithm and the correct traversal by Fleury's algorithm with an extra constraint.
The extra constraint is that during traversal travelling between three collinear points is disallowed.
Due to lack of better name as of now, I'll refer to this as : outer traversal, since one does not cross any other segment.

Is there a name for the polygon when it is intended to be interpreted sharing vertices between portions, not crossing sides?

Would the term complex polygon suffice ?

What properties besides area depend on this behavior?

Any integral computed on such a polytope without outer traversal would be incorrect.

So, I think computing intersection using BO algorithm and then traversal using Fleury's algorithm would be the correct solution. Also, the default value of ignore_intersections should be changed to True.
Waiting for your thoughts.

@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 12, 2017

Member

Leave a comment

Unfortunately they do not :

sorry -- can only leave a brief comment. I thought that in the passing polytope tests you shut off intersections. And I was actually referring to the regular area (not the polytope integration) so pimplicitintersections.area == pexplicitintersections.area`

Would the term complex polygon suffice ?

Or maybe simply hull (instead of convex hull).

default value...should be changed to True

I think so, given that the area (not polytope integration) is the same -- or am I wrong about that?

Member

smichr commented Jul 12, 2017

Leave a comment

Unfortunately they do not :

sorry -- can only leave a brief comment. I thought that in the passing polytope tests you shut off intersections. And I was actually referring to the regular area (not the polytope integration) so pimplicitintersections.area == pexplicitintersections.area`

Would the term complex polygon suffice ?

Or maybe simply hull (instead of convex hull).

default value...should be changed to True

I think so, given that the area (not polytope integration) is the same -- or am I wrong about that?

@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 12, 2017

Member

📝 note to self, see what happens with p=Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0)) and is this a way to put a hole in a polygon Polygon((0,0),(1,1),(2,0),(1,0),(1,1/8),(3/2,1/4),(1/2,1/4),(1,1/8),(1,0))

Member

smichr commented Jul 12, 2017

📝 note to self, see what happens with p=Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0)) and is this a way to put a hole in a polygon Polygon((0,0),(1,1),(2,0),(1,0),(1,1/8),(3/2,1/4),(1/2,1/4),(1,1/8),(1,0))

@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 12, 2017

Member

I see that the hull would not be unique for a collection of points. See this for example.

Member

smichr commented Jul 12, 2017

I see that the hull would not be unique for a collection of points. See this for example.

@ArifAhmed1995

This comment has been minimized.

Show comment
Hide comment
@ArifAhmed1995

ArifAhmed1995 Jul 12, 2017

Contributor

@smichr

I thought that in the passing polytope tests you shut off intersections. And I was actually referring to the regular area (not the polytope integration) so pimplicitintersections.area == pexplicitintersections.area
If only area depends on it, we might just include a property "footprint_area" to compute the area covered by the polygon rather than the more abstract area that depends on the direction in which the sides are drawn.

Oops ! I was mistaken. The results by polytope_integrate are correct.
Now, I don't think the abstract area which depends on direction of sides is required at all. Users will expect area to mean footprint_area. So, when regardless of the value of ignore_intersections the area method should look at only initial vertices supplied and ignore intersections.

see what happens with p=Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0)) and is this a way to put a hole in a polygon Polygon((0,0),(1,1),(2,0),(1,0),(1,1/8),(3/2,1/4),(1/2,1/4),(1,1/8),(1,0))

It seems that neither of those cases work as expected. I'll look at it again when I wake up :

>>> Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0))
Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0), Point2D(1, 0), Point2D(1, 1), Point2D(2, 1), Point2D(1, 0))
>>> Polygon((0,0),(1,1),(2,0),(1,0),(1,1/8),(3/2,1/4),(1/2,1/4),(1,1/8),(1,0))
Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0), Point2D(1, 0), Point2D(0, 0), Point2D(1, 0))
Contributor

ArifAhmed1995 commented Jul 12, 2017

@smichr

I thought that in the passing polytope tests you shut off intersections. And I was actually referring to the regular area (not the polytope integration) so pimplicitintersections.area == pexplicitintersections.area
If only area depends on it, we might just include a property "footprint_area" to compute the area covered by the polygon rather than the more abstract area that depends on the direction in which the sides are drawn.

Oops ! I was mistaken. The results by polytope_integrate are correct.
Now, I don't think the abstract area which depends on direction of sides is required at all. Users will expect area to mean footprint_area. So, when regardless of the value of ignore_intersections the area method should look at only initial vertices supplied and ignore intersections.

see what happens with p=Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0)) and is this a way to put a hole in a polygon Polygon((0,0),(1,1),(2,0),(1,0),(1,1/8),(3/2,1/4),(1/2,1/4),(1,1/8),(1,0))

It seems that neither of those cases work as expected. I'll look at it again when I wake up :

>>> Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0))
Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0), Point2D(1, 0), Point2D(1, 1), Point2D(2, 1), Point2D(1, 0))
>>> Polygon((0,0),(1,1),(2,0),(1,0),(1,1/8),(3/2,1/4),(1/2,1/4),(1,1/8),(1,0))
Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0), Point2D(1, 0), Point2D(0, 0), Point2D(1, 0))
@ArifAhmed1995

This comment has been minimized.

Show comment
Hide comment
@ArifAhmed1995

ArifAhmed1995 Jul 13, 2017

Contributor

Actually, I ran the above in Python 2.7 so forgot about the integer division. The second one is fine.

>>> Polygon((0,0),(1,1),(2,0),(1,0),(1,S(1)/8),(S(3)/2,S(1)/4),(S(1)/2,S(1)/4),(1,S(1)/8),(1,0))
Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0), Point2D(1, 0), Point2D(1, 1/8), Point2D(3/2, 1/4), Point2D(1/2, 1/4), Point2D(1, 1/8), Point2D(1, 0))

The first one is incorrectly classified as being convex.
If I change that check then it gives correct result(incorrect traversal but correctly finds intersections) :
Change :

if not convex and not ignore_int:
    ...

to

if not ignore_int:
    ...

then result :

>>> Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0))
Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(3/2, 1/2), Point2D(2, 0), Point2D(1, 0), Point2D(1, 1), Point2D(2, 1), Point2D(3/2, 1/2), Point2D(1, 0))
Contributor

ArifAhmed1995 commented Jul 13, 2017

Actually, I ran the above in Python 2.7 so forgot about the integer division. The second one is fine.

>>> Polygon((0,0),(1,1),(2,0),(1,0),(1,S(1)/8),(S(3)/2,S(1)/4),(S(1)/2,S(1)/4),(1,S(1)/8),(1,0))
Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0), Point2D(1, 0), Point2D(1, 1/8), Point2D(3/2, 1/4), Point2D(1/2, 1/4), Point2D(1, 1/8), Point2D(1, 0))

The first one is incorrectly classified as being convex.
If I change that check then it gives correct result(incorrect traversal but correctly finds intersections) :
Change :

if not convex and not ignore_int:
    ...

to

if not ignore_int:
    ...

then result :

>>> Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0))
Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(3/2, 1/2), Point2D(2, 0), Point2D(1, 0), Point2D(1, 1), Point2D(2, 1), Point2D(3/2, 1/2), Point2D(1, 0))
@ArifAhmed1995

This comment has been minimized.

Show comment
Hide comment
@ArifAhmed1995

ArifAhmed1995 Jul 13, 2017

Contributor

@smichr
I think that the traversal should be done with BO algorithm and the area should be considering only vertices supplied. For example, if one calculates the area of the polygon

Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0))

using result as given by Fleury's traversal then the result will be incorrect cause then the Triangle((1, 0), (S(3)/2, S(1)/2), (1, 1)) is ignored.

Maybe your code could be added to a Polygon.hull property which would return the non-crossing outline of a polygon. A convex polygon would return the convex hull; an intersecting one would return a ccw traversed polygon (whose area would then be positive) whose segments never cross each other but only meet at vertices (which your code currently detects).

But on the other hand while a hull algorithm gives correct area it loses information about the polygon since then the Segment2D((1, 0), (1, 1)) is ignored.

So, in total would the following be a good solution ? :

  • Calculate area ignoring intersections(that is footprint_area), cause there is no use for the area that depends on segment direction.

  • Compute normal traversal by figuring out intersection points by BO algorithm and then place that intersection point between two endpoints of the segment. That is Segment2D((1, 1), (2, 0)) --> Segment2D((1, 1), (S(3)/2, S(1)/2)) + Segment2D((S(3)/2, S(1)/2), (2, 0)).

  • Add a property Polygon.hull as explained by you above.

Contributor

ArifAhmed1995 commented Jul 13, 2017

@smichr
I think that the traversal should be done with BO algorithm and the area should be considering only vertices supplied. For example, if one calculates the area of the polygon

Polygon((0,0),(1,1),(2,0),(1,0),(1,1),(2,1),(1,0))

using result as given by Fleury's traversal then the result will be incorrect cause then the Triangle((1, 0), (S(3)/2, S(1)/2), (1, 1)) is ignored.

Maybe your code could be added to a Polygon.hull property which would return the non-crossing outline of a polygon. A convex polygon would return the convex hull; an intersecting one would return a ccw traversed polygon (whose area would then be positive) whose segments never cross each other but only meet at vertices (which your code currently detects).

But on the other hand while a hull algorithm gives correct area it loses information about the polygon since then the Segment2D((1, 0), (1, 1)) is ignored.

So, in total would the following be a good solution ? :

  • Calculate area ignoring intersections(that is footprint_area), cause there is no use for the area that depends on segment direction.

  • Compute normal traversal by figuring out intersection points by BO algorithm and then place that intersection point between two endpoints of the segment. That is Segment2D((1, 1), (2, 0)) --> Segment2D((1, 1), (S(3)/2, S(1)/2)) + Segment2D((S(3)/2, S(1)/2), (2, 0)).

  • Add a property Polygon.hull as explained by you above.

@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 13, 2017

Member

I can look at this over the next few days, but will be a little slow in responding perhaps. I am interested but don't have a lot of free time at the moment. I'll try ring in this evening.

Member

smichr commented Jul 13, 2017

I can look at this over the next few days, but will be a little slow in responding perhaps. I am interested but don't have a lot of free time at the moment. I'll try ring in this evening.

@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 15, 2017

Member

cause there is no use for the area that depends on segment direction.

Ahh, but there is :-) That is the purpose of the polygon with a hole that I gave above and present again here:

This is the polygon with a "hole" defined as a triangle drawn in the ccw direction while the large triangle is drawn in the cw direction.

>>> p = Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0), Point2D(1, 0),
Point2D(1, 1/8), Point2D(3/2, 1/4), Point2D(1/2, 1/4), Point2D(1,
1/8), Point2D(1, 0))

Here is the hole alone:

>>> h = Triangle(Point2D(1, 1/8), Point2D(3/2, 1/4), Point2D(1/2, 1/4))

Here is the large triangle

b = Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0))

Now the areas:

>>> b.area
-1
>>> h.area
1/16
>>> p.area
-15/16
Member

smichr commented Jul 15, 2017

cause there is no use for the area that depends on segment direction.

Ahh, but there is :-) That is the purpose of the polygon with a hole that I gave above and present again here:

This is the polygon with a "hole" defined as a triangle drawn in the ccw direction while the large triangle is drawn in the cw direction.

>>> p = Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0), Point2D(1, 0),
Point2D(1, 1/8), Point2D(3/2, 1/4), Point2D(1/2, 1/4), Point2D(1,
1/8), Point2D(1, 0))

Here is the hole alone:

>>> h = Triangle(Point2D(1, 1/8), Point2D(3/2, 1/4), Point2D(1/2, 1/4))

Here is the large triangle

b = Polygon(Point2D(0, 0), Point2D(1, 1), Point2D(2, 0))

Now the areas:

>>> b.area
-1
>>> h.area
1/16
>>> p.area
-15/16
@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 15, 2017

Member

Something here may be of interest though it refers to even-odd rule area.

Member

smichr commented Jul 15, 2017

Something here may be of interest though it refers to even-odd rule area.

Show outdated Hide outdated sympy/geometry/line.py
p1, p2 = p2, p1
elif (p1.x == p2.x) == True and (p1.y > p2.y) == True:
p1, p2 = p2, p1
dir = kwargs.get('dir', False)

This comment has been minimized.

@smichr

smichr Jul 16, 2017

Member

The docstring should be edited to tell that a Segment can now be directed. Should the same thing be done for the 3D version of Segment?

Did you consider using the keyword 'evaluate' instead of 'dir'? That's a common way to get an object returned with the arguments in the order given, e.g. Mul(x, 2, evaluate=False).args == (x, 2)

@smichr

smichr Jul 16, 2017

Member

The docstring should be edited to tell that a Segment can now be directed. Should the same thing be done for the 3D version of Segment?

Did you consider using the keyword 'evaluate' instead of 'dir'? That's a common way to get an object returned with the arguments in the order given, e.g. Mul(x, 2, evaluate=False).args == (x, 2)

This comment has been minimized.

@ArifAhmed1995

ArifAhmed1995 Jul 17, 2017

Contributor

I'll do it once the issue with area is cleared up. evaluate seems a little strange at first but I guess it will be better since it's known to be used to keep args in order.

@ArifAhmed1995

ArifAhmed1995 Jul 17, 2017

Contributor

I'll do it once the issue with area is cleared up. evaluate seems a little strange at first but I guess it will be better since it's known to be used to keep args in order.

@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 17, 2017

Member

This fails, I think: (A: The first one is incorrectly classified as being convex.)

assert Point(1, 5/S(2)) in Polygon(
(0, 0), (5, 0),
(1, 5), (1, 1),
(2, 2), (0, 3), ignore_int=False).args

is_convex also needs to check whether there are any collinear points -- if a polygon is self intersecting, it isn't convex.

Here is a first attempt at footprint_area:

def footprint_area(p):
    pargs = list(p.args)
    if pargs[-1] != pargs[0]:
        pargs.append(pargs[0])
    for i in range(len(pargs)):
        if pargs[0].is_collinear(*pargs[1:3]):
            pargs.append(pargs.pop(0))
        else:
            p1, p2 = pargs[:2]
            area = (p2.x - p1.x)*(p1.y + p2.y)/2
            dir = 1
            narea = 0
            break
    else:
        return abs(p.area)
    for pt in pargs[2:]:
        if p1.is_collinear(p2, pt):
            dir = -dir
        da = (pt.x - p2.x)*(pt.y + p2.y)/2
        if dir == 1:
            area += da
        else:
            narea -= da
        p1, p2 = p2, pt
    return max(abs(area), abs(narea))

pts = [(0, 0), (5, 0), (1, 5), (1, 1), (2, 2), (0, 3)]
for i in range(len(pts)):
    pts.append(pts.pop(0))
    assert footprint_area(Polygon(*pts, ignore_int=False)) == 51/S(4)
Member

smichr commented Jul 17, 2017

This fails, I think: (A: The first one is incorrectly classified as being convex.)

assert Point(1, 5/S(2)) in Polygon(
(0, 0), (5, 0),
(1, 5), (1, 1),
(2, 2), (0, 3), ignore_int=False).args

is_convex also needs to check whether there are any collinear points -- if a polygon is self intersecting, it isn't convex.

Here is a first attempt at footprint_area:

def footprint_area(p):
    pargs = list(p.args)
    if pargs[-1] != pargs[0]:
        pargs.append(pargs[0])
    for i in range(len(pargs)):
        if pargs[0].is_collinear(*pargs[1:3]):
            pargs.append(pargs.pop(0))
        else:
            p1, p2 = pargs[:2]
            area = (p2.x - p1.x)*(p1.y + p2.y)/2
            dir = 1
            narea = 0
            break
    else:
        return abs(p.area)
    for pt in pargs[2:]:
        if p1.is_collinear(p2, pt):
            dir = -dir
        da = (pt.x - p2.x)*(pt.y + p2.y)/2
        if dir == 1:
            area += da
        else:
            narea -= da
        p1, p2 = p2, pt
    return max(abs(area), abs(narea))

pts = [(0, 0), (5, 0), (1, 5), (1, 1), (2, 2), (0, 3)]
for i in range(len(pts)):
    pts.append(pts.pop(0))
    assert footprint_area(Polygon(*pts, ignore_int=False)) == 51/S(4)
@ArifAhmed1995

This comment has been minimized.

Show comment
Hide comment
@ArifAhmed1995

ArifAhmed1995 Jul 17, 2017

Contributor

This fails, I think:
assert Point(1, 5/S(2)) in Polygon(
(0, 0), (5, 0),
(1, 5), (1, 1),
(2, 2), (0, 3), ignore_int=False).args

The test passes :

>>> Point(1, 5/S(2)) in Polygon((0, 0), (5, 0),(1, 5), (1, 1),(2, 2), (0, 3), ignore_int=False).args
True
>>> Polygon((0, 0), (5, 0),(1, 5), (1, 1),(2, 2), (0, 3), ignore_int=False)
Polygon(Point2D(0, 0), Point2D(5, 0), Point2D(1, 5), Point2D(1, 5/2), Point2D(1, 1), Point2D(2, 2), Point2D(1, 5/2), Point2D(0, 3))

I implemented one of the techniques mentioned on codegolf and it ran fine for the following two polygons :

>>> p1 = Polygon((0,-3), (-1,-2), (2,0), (-1,2), (0,3))
>>> p1.area
-13/3
>>> p2 = Polygon((0,0), (1,0), (0,1), (1,1))
>>> p2.area
1/2
Contributor

ArifAhmed1995 commented Jul 17, 2017

This fails, I think:
assert Point(1, 5/S(2)) in Polygon(
(0, 0), (5, 0),
(1, 5), (1, 1),
(2, 2), (0, 3), ignore_int=False).args

The test passes :

>>> Point(1, 5/S(2)) in Polygon((0, 0), (5, 0),(1, 5), (1, 1),(2, 2), (0, 3), ignore_int=False).args
True
>>> Polygon((0, 0), (5, 0),(1, 5), (1, 1),(2, 2), (0, 3), ignore_int=False)
Polygon(Point2D(0, 0), Point2D(5, 0), Point2D(1, 5), Point2D(1, 5/2), Point2D(1, 1), Point2D(2, 2), Point2D(1, 5/2), Point2D(0, 3))

I implemented one of the techniques mentioned on codegolf and it ran fine for the following two polygons :

>>> p1 = Polygon((0,-3), (-1,-2), (2,0), (-1,2), (0,3))
>>> p1.area
-13/3
>>> p2 = Polygon((0,0), (1,0), (0,1), (1,1))
>>> p2.area
1/2
@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 17, 2017

Member

I think there is a problem creating a polygon with these points:

pts = ((0,0),(3,0),(3,3),(0,3),(0,-1),(2,-1),(-1,2),(1,2),(-1,0))
Polygon(*pts, ignore_int=False)

A segment with equal points is created which collapses to that point. When accessed later as though it is a segment is when the error is raised:

Traceback (most recent call last):
  File "C:\Python34\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python34\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python34\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
    exec(cmd, globals, locals)
  File "jnk.py", line 18, in <module>
    Polygon(*pts, ignore_int=False)
  File "sympy\sympy\geometry\polygon.py", line 194, in __new__
    [Segment(pts[0], hit, dir=True),
  File "sympy\sympy\geometry\line.py", line 1417, in __new__
    p1, p2 = Point._normalize_dimension(Point(p1), Point(p2))
  File "sympy\sympy\geometry\point.py", line 119, in __new__
    if len(coords) == 0 and kwargs.get('dim', None):
TypeError: object of type 'Zero' has no len()
Member

smichr commented Jul 17, 2017

I think there is a problem creating a polygon with these points:

pts = ((0,0),(3,0),(3,3),(0,3),(0,-1),(2,-1),(-1,2),(1,2),(-1,0))
Polygon(*pts, ignore_int=False)

A segment with equal points is created which collapses to that point. When accessed later as though it is a segment is when the error is raised:

Traceback (most recent call last):
  File "C:\Python34\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 323, in RunScript
    debugger.run(codeObject, __main__.__dict__, start_stepping=0)
  File "C:\Python34\Lib\site-packages\pythonwin\pywin\debugger\__init__.py", line 60, in run
    _GetCurrentDebugger().run(cmd, globals,locals, start_stepping)
  File "C:\Python34\Lib\site-packages\pythonwin\pywin\debugger\debugger.py", line 654, in run
    exec(cmd, globals, locals)
  File "jnk.py", line 18, in <module>
    Polygon(*pts, ignore_int=False)
  File "sympy\sympy\geometry\polygon.py", line 194, in __new__
    [Segment(pts[0], hit, dir=True),
  File "sympy\sympy\geometry\line.py", line 1417, in __new__
    p1, p2 = Point._normalize_dimension(Point(p1), Point(p2))
  File "sympy\sympy\geometry\point.py", line 119, in __new__
    if len(coords) == 0 and kwargs.get('dim', None):
TypeError: object of type 'Zero' has no len()
@ArifAhmed1995

This comment has been minimized.

Show comment
Hide comment
@ArifAhmed1995

ArifAhmed1995 Jul 22, 2017

Contributor

@smichr

pts = ((0,0),(3,0),(3,3),(0,3),(0,-1),(2,-1),(-1,2),(1,2),(-1,0))
Polygon(*pts, ignore_int=False)

I've put a fix in, now this is the output :

>>> pts = ((0,0),(3,0),(3,3),(0,3),(0,-1),(2,-1),(-1,2),(1,2),(-1,0))
>>> p = Polygon(*pts, ignore_intersection=False)
>>> p
Polygon(Point2D(3, 0), Point2D(3, 3), Point2D(0, 3), Point2D(0, 2), Point2D(0, 1), Point2D(0, 0), Point2D(0, -1), Point2D(2, -1), Point2D(1, 0), Point2D(0, 1), Point2D(-1, 2), Point2D(0, 2), Point2D(1, 2), Point2D(0, 1), Point2D(-1, 0), Point2D(0, 0), Point2D(1, 0))
>>> p.sides
[Segment2D(Point2D(3, 0), Point2D(3, 3)), Segment2D(Point2D(0, 3), Point2D(3, 3)), Segment2D(Point2D(0, 2), Point2D(0, 3)), Segment2D(Point2D(0, 1), Point2D(0, 2)), Segment2D(Point2D(0, 0), Point2D(0, 1)), Segment2D(Point2D(0, -1), Point2D(0, 0)), Segment2D(Point2D(0, -1), Point2D(2, -1)), Segment2D(Point2D(1, 0), Point2D(2, -1)), Segment2D(Point2D(0, 1), Point2D(1, 0)), Segment2D(Point2D(-1, 2), Point2D(0, 1)), Segment2D(Point2D(-1, 2), Point2D(0, 2)), Segment2D(Point2D(0, 2), Point2D(1, 2)), Segment2D(Point2D(0, 1), Point2D(1, 2)), Segment2D(Point2D(-1, 0), Point2D(0, 1)), Segment2D(Point2D(-1, 0), Point2D(0, 0)), Segment2D(Point2D(0, 0), Point2D(1, 0)), Segment2D(Point2D(1, 0), Point2D(3, 0))]

But this technique for finding intersections points is quite inefficient and will have to be replaced by the BO algorithm. This is the reason why I didn't put all the tests mentioned on the codegolf site in the test_polygon file. The answers returned are correct but calculating intersection points takes too long. Once the points are found out the area calculation is fast.
However, the area algorithm doesn't work when there are more than two lines passing through same intersection point. For the above polygon the program seems to spiral into an infinite loop.

Here is a first attempt at footprint_area:

The code doesn't seem to work for the following polygon :

>>> p = Polygon((0,0), (1,0), (0,1), (1,1), ignore_intersection=False)
>>> p.footprint_area
1/4

Am I right in thinking that footprint_area should return the area given by the hull of the polygon ?
Can you point me to a source that explains the logic ? I could not understand it by looking at the code.

So, there seems to be two issues left :

  • Very slow algorithm for calculating the intersections.
  • Area calculation doesn't work for polygons having more than two sides passing through the same point.

Can I request that I try to fix these two issues in separate PRs and this current one be reviewed and merged ? The reason is to bring closure to the 2D use case of the intpoly module.

Contributor

ArifAhmed1995 commented Jul 22, 2017

@smichr

pts = ((0,0),(3,0),(3,3),(0,3),(0,-1),(2,-1),(-1,2),(1,2),(-1,0))
Polygon(*pts, ignore_int=False)

I've put a fix in, now this is the output :

>>> pts = ((0,0),(3,0),(3,3),(0,3),(0,-1),(2,-1),(-1,2),(1,2),(-1,0))
>>> p = Polygon(*pts, ignore_intersection=False)
>>> p
Polygon(Point2D(3, 0), Point2D(3, 3), Point2D(0, 3), Point2D(0, 2), Point2D(0, 1), Point2D(0, 0), Point2D(0, -1), Point2D(2, -1), Point2D(1, 0), Point2D(0, 1), Point2D(-1, 2), Point2D(0, 2), Point2D(1, 2), Point2D(0, 1), Point2D(-1, 0), Point2D(0, 0), Point2D(1, 0))
>>> p.sides
[Segment2D(Point2D(3, 0), Point2D(3, 3)), Segment2D(Point2D(0, 3), Point2D(3, 3)), Segment2D(Point2D(0, 2), Point2D(0, 3)), Segment2D(Point2D(0, 1), Point2D(0, 2)), Segment2D(Point2D(0, 0), Point2D(0, 1)), Segment2D(Point2D(0, -1), Point2D(0, 0)), Segment2D(Point2D(0, -1), Point2D(2, -1)), Segment2D(Point2D(1, 0), Point2D(2, -1)), Segment2D(Point2D(0, 1), Point2D(1, 0)), Segment2D(Point2D(-1, 2), Point2D(0, 1)), Segment2D(Point2D(-1, 2), Point2D(0, 2)), Segment2D(Point2D(0, 2), Point2D(1, 2)), Segment2D(Point2D(0, 1), Point2D(1, 2)), Segment2D(Point2D(-1, 0), Point2D(0, 1)), Segment2D(Point2D(-1, 0), Point2D(0, 0)), Segment2D(Point2D(0, 0), Point2D(1, 0)), Segment2D(Point2D(1, 0), Point2D(3, 0))]

But this technique for finding intersections points is quite inefficient and will have to be replaced by the BO algorithm. This is the reason why I didn't put all the tests mentioned on the codegolf site in the test_polygon file. The answers returned are correct but calculating intersection points takes too long. Once the points are found out the area calculation is fast.
However, the area algorithm doesn't work when there are more than two lines passing through same intersection point. For the above polygon the program seems to spiral into an infinite loop.

Here is a first attempt at footprint_area:

The code doesn't seem to work for the following polygon :

>>> p = Polygon((0,0), (1,0), (0,1), (1,1), ignore_intersection=False)
>>> p.footprint_area
1/4

Am I right in thinking that footprint_area should return the area given by the hull of the polygon ?
Can you point me to a source that explains the logic ? I could not understand it by looking at the code.

So, there seems to be two issues left :

  • Very slow algorithm for calculating the intersections.
  • Area calculation doesn't work for polygons having more than two sides passing through the same point.

Can I request that I try to fix these two issues in separate PRs and this current one be reviewed and merged ? The reason is to bring closure to the 2D use case of the intpoly module.

p1, p2 = p2, p1
elif (p1.x == p2.x) == True and (p1.y > p2.y) == True:
p1, p2 = p2, p1
evaluate = kwargs.get('evaluate', False)

This comment has been minimized.

@smichr

smichr Jul 22, 2017

Member

You got the sense wrong, I think: if evaluate=False you shouldn't try to order the points. Would anything break if Segments didn't try to reorder their points? Line doesn't do so. All that might need to be done is the equals method be modified to return isinstance(o, self.func) and list(ordered(self.args)) == list(ordered(o.args)) instead of self == o.

@smichr

smichr Jul 22, 2017

Member

You got the sense wrong, I think: if evaluate=False you shouldn't try to order the points. Would anything break if Segments didn't try to reorder their points? Line doesn't do so. All that might need to be done is the equals method be modified to return isinstance(o, self.func) and list(ordered(self.args)) == list(ordered(o.args)) instead of self == o.

# part of the convex hull can possibly intersect with other
# sides of the polygon...but for now we use the n**2 algorithm
# and check if any side intersects with any preceding side,
# excluding the ones it is connected to
try:

This comment has been minimized.

@smichr

smichr Jul 22, 2017

Member

convex is no longer used. But if intersecting polygons are allowed, that method needs to be modified to check that there are no collinear segments (which should only happen if a point was needed for a self intersecting polygon). It might be better to store this attribute (is_self_intersecting) at instantiation rather than have to compute it again for testing whether the polygon is convex or not. For that matter, should the convex attribute be stored at instantiation?

@smichr

smichr Jul 22, 2017

Member

convex is no longer used. But if intersecting polygons are allowed, that method needs to be modified to check that there are no collinear segments (which should only happen if a point was needed for a self intersecting polygon). It might be better to store this attribute (is_self_intersecting) at instantiation rather than have to compute it again for testing whether the polygon is convex or not. For that matter, should the convex attribute be stored at instantiation?

@@ -389,6 +455,49 @@ def centroid(self):
return Point(simplify(A*cx), simplify(A*cy))
@property
def directed_sides(self):

This comment has been minimized.

@smichr

smichr Jul 22, 2017

Member

if segments no longer modify the point order then this wouldn't be needed.

@smichr

smichr Jul 22, 2017

Member

if segments no longer modify the point order then this wouldn't be needed.

This comment has been minimized.

@ArifAhmed1995

ArifAhmed1995 Jul 26, 2017

Contributor

Allright, so evaluate is not needed ? Will it be better if all sides are designated to be directed ?

@ArifAhmed1995

ArifAhmed1995 Jul 26, 2017

Contributor

Allright, so evaluate is not needed ? Will it be better if all sides are designated to be directed ?

@smichr

This comment has been minimized.

Show comment
Hide comment
@smichr

smichr Jul 22, 2017

Member

current one be reviewed and merged

So as a first step (correct me if I am wrong):

  • change Segment to not order args
  • allow Polygon to have intersecting sides
  • track whether the polygon is self intersecting and raise a NotImplementedError if the area of such is requested
  • modify is_convex to return False if the polygon is self intersecting
Member

smichr commented Jul 22, 2017

current one be reviewed and merged

So as a first step (correct me if I am wrong):

  • change Segment to not order args
  • allow Polygon to have intersecting sides
  • track whether the polygon is self intersecting and raise a NotImplementedError if the area of such is requested
  • modify is_convex to return False if the polygon is self intersecting

@gxyd gxyd added the geometry label Oct 1, 2017

@smichr smichr closed this in #13846 Jan 6, 2018

@ArifAhmed1995 ArifAhmed1995 deleted the ArifAhmed1995:polygon_intersect branch Jan 6, 2018

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment