11
$\begingroup$

So Mathematica uses an even-odd rule for self-intersecting polygons, both when rendering them and when calculating their area:

poly = Polygon@{{0., 0}, {3, 0}, {3, 2}, {1, 2}, {1, 1}, {2, 1}, {2, 3}, {0, 3}};
Graphics@poly
(* See image below *)
Area@poly
(* 7. *)
RegionMeasure@poly
(* 7. *)

enter image description here

But the area calculation doesn't seem to work when one of the vertices is repeated:

poly = Polygon@{{0., 0}, {3, 0}, {3, 3}, {1, 3}, {1, 1}, {3, 1}, {3, 3}, {0, 3}}
Graphics@poly
(* See image below *)
Area@poly
(* 13. *)
RegionMeasure@poly
(* 13. *)

enter image description here

I'd expect the area to be 5 as can be shown by slightly offsetting the points in the corner:

poly = Polygon@{{0., 0}, {3, 0}, {3, 2.999}, {1, 2.999}, {1, 1}, {2.999, 1}, {2.999, 3}, {0, 3}}
Area@poly
(* 5.004 *)
RegionMeasure@poly
(* 5.004 *)

It seems that the areas of the convex parts of the polygon are actually added instead of subtracted if I've got repeated vertices.

This looks likes a bug to me, but is it? Is there a fix other than applying random offsets on the order of the machine epsilon to the vertices?

Note that I am not interested in fixing the above example to a simpler non-intersecting polygon unless this can be done programmatically in the general case. I've only chosen such a simple example to illustrate the point. Assume that the polygon could be arbitrarily complex with repeated vertices.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ If you use exact coordinates, Polygon@{{0, 0}, {3, 0}, {3, 3}, {1, 3}, {1, 1}, {3, 1}, {3, 3}, {0, 3}}, the computed area is 4 instead of 5. $\endgroup$
    – Michael E2
    Commented Mar 11, 2015 at 13:55
  • 1
    $\begingroup$ I'm interpreting your problem like this: "Given an arbitrary closed polygonal path, how can one compute the area of the polygon according to the even-odd rule?" Is that correct? $\endgroup$
    – Michael E2
    Commented Mar 11, 2015 at 13:59
  • $\begingroup$ @MichaelE2 4 is even weirder than 5. If I call DiscretizeGraphics first I still get 13 by the way. As for your rephrasing, yes, that's fairly accurate, although I'd like to stress that I'm not looking for general polygon area algorithms, but rather for a way to make Area or RegionMeasure work for general polygons (or a good reason why they don't work in the presence of repeated vertices). But it looks like your answer provides that (at least the former). :) $\endgroup$
    – Martin Ender
    Commented Mar 11, 2015 at 16:40
  • $\begingroup$ I realize what you're looking for, but I don't know how to achieve it without some intervention (as in my answer). BTW, the 4 seems to be the area of the excluded square, as if it got the even-odd rule backwards; the 13 can be obtained by double-counting the excluded square. At this point, that seems like a bug to me. $\endgroup$
    – Michael E2
    Commented Mar 11, 2015 at 17:01

2 Answers 2

Reset to default
8
$\begingroup$

You could triangulate the polygon, pick out the triangles inside the polygon, and compute their area.

area[poly_Polygon] := 
 With[{mpoly = TriangulateMesh[DiscretizeGraphics[poly], 
                 "MaxCellMeasure" -> Area@Graphics`Mesh`BoundingAxisAlignedRectangle[poly]]},
  With[{pts = MeshCoordinates[mpoly],
        tri = MeshCells[mpoly, 2]},
   Area@MeshRegion[
     pts,
     Pick[
      tri, 
      tri /. Polygon[v_] :> 
        Graphics`PolygonUtils`InPolygonQ[poly, Mean[pts[[v]]]]
      ]
     ]
   ]]

poly = Polygon@{{0., 0}, {3, 0}, {3, 3}, {1, 3}, {1, 1}, {3, 1}, {3, 3}, {0, 3}};
area[poly]
(*  5.  *)

poly2 = Polygon@{{0., 0}, {3, 0}, {3, 2}, {1, 2}, {1, 1}, {2, 1}, {2, 3}, {0, 3}};
area[poly2]
(*  7.  *)

Response to comment: Bug indeed!!! (I'm using V10.0.2, OSX 10.10.2.; reported to Wolfram Support.)

DiscretizeRegion[
 Polygon@N@{{0, 0}, {2, 0}, {2, 2}, {1, 0}, {0, 2}}, 
 Method -> "RegionPlot"]

Mathematica graphics

Note that

RegionPlot[Polygon@N@{{0, 0}, {2, 0}, {2, 2}, {1, 0}, {0, 2}}]

works fine.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I'm accepting this, because it does indeed provide a general workaround. That being said, I'd still be interested in why Area and RegionMeasure can't handle this as they handle any other self-intersecting polygon, so if someone has an answer for that, please do post it. $\endgroup$
    – Martin Ender
    Commented Mar 11, 2015 at 16:49
  • 1
    $\begingroup$ It seems that this runs into problems if a vertex coincides with another edge of the polygon. Try Polygon@N@{{0, 0}, {2, 0}, {2, 2}, {1, 0}, {0, 2}}. $\endgroup$
    – Martin Ender
    Commented Mar 11, 2015 at 18:23
  • $\begingroup$ @MartinBüttner Please see my update. $\endgroup$
    – Michael E2
    Commented Mar 12, 2015 at 17:57
  • $\begingroup$ Re your bug report, 10.3 doesn't respond with an internal error any more but simply DiscretizeRegion::drf: DiscretizeRegion was unable to discretize the region [...] $\endgroup$
    – Martin Ender
    Commented Jan 22, 2016 at 14:31
4
$\begingroup$

I think it could probably issue a message. Perhaps you could report that to the technical service.

If you use the coordinates in an order 

mr = DiscretizeGraphics[
  Polygon[{{0., 0}, {3, 0}, {3, 1}, {1, 1}, {1, 3}, {0, 3}}]];

Area[mr]

things work as expected.

To check if coordinates have duplicates one can use 

{coords, map} = 
 Region`Mesh`DeleteDuplicateCoordinates[{{0., 0}, {3, 0}, {3, 3}, {1, 
    3}, {1, 1}, {3, 1}, {3, 3}, {0, 3}}]
{{{0.`, 0.`}, {3.`, 0.`}, {3.`, 3.`}, {1.`, 3.`}, {1.`, 1.`}, {3.`, 
   1.`}, {0.`, 3.`}}, {1, 2, 3, 4, 5, 6, 3, 7}}

This returns the coordinates without the duplicates and an incidence map. In the case there are no duplicates then the map is the Indentity. You could check for that. But again, the coordinate specification needs to be encircle the area in an systematic way.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Of course, the first example is no longer self-intersecting, whereas the second gives a completely different polygon than the one I'm interested. Ideally I'd be looking for a general way to deal with the self-intersecting ones without having to analyse them for an equivalent polygon without self-intersections (which may be a nontrivial task to begin with). $\endgroup$
    – Martin Ender
    Commented Mar 11, 2015 at 10:44
  • $\begingroup$ @MartinBüttner, I think it were good if you clarified your questions. Is it a bug: No, more of a messaging issue, what's a fix: give coordinates in a proper manner. If you want something else, it might help explaining that in the question. Concerning the second example: this was to demonstrate how you detect duplicate coordinates. The incidence map is not the new polygon. It's mapping incidences of, say, a GraphicsComplex to the new order. But you do not have that here. $\endgroup$
    – user21
    Commented Mar 11, 2015 at 11:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.