5
$\begingroup$

Consider the following code:

tmp = {{0, 0}, {1, 1}, {2, 1}, {3, 2}, {1, 0.5}};
ListLinePlot[tmp, Filling -> Axis]

Is there any easy way to compute filled area?

Edit: Actual data consists of thousands of points, so some analytical solution would be the best

$\endgroup$
1
  • $\begingroup$ ComponentMeasurements[ColorNegate@Binarize@ListLinePlot[tmp, Filling -> Axis, Axes -> False], "Area"] $\endgroup$
    – Sektor
    Feb 7, 2015 at 20:23

2 Answers 2

7
$\begingroup$
tmp = {{0, 0}, {1, 1}, {2, 1}, {3, 2}, {1, 0.5}};
plot = ListLinePlot[tmp, Filling -> Axis];

Area@DiscretizeGraphics[plot]
(*  0.833333  *)

I figured out how to use Graphics`PolygonUtils`SimplePolygonPartition. It subdivides the polygon from the plot into non self-intersecting, possibly nonconvex polygons, but some of the polygons it creates lie outside the original polygon. One difficulty is finding a point inside a polygon. We do that by searching for a point, the average of three consecutive vertices, that is inside a convex angle. The subdivision also creates vertices along edges. Numerical round-off error makes it difficult to detect when point is inside or outside the polygon, so we skip those.

ClearAll[findPtInPoly];
SetAttributes[findPtInPoly, Listable];
findPtInPoly::nopt = "Warning: Could not find point inside polygon ``; returning a vertex";
tolerance = 1*^-10;
findPtInPoly[Polygon[poly_]] := 
 Module[{point}, 
  Do[With[{pts = poly[[t ;; t + 2]]}, 
    If[VectorAngle @@ Differences[pts] > tolerance && (* == not collinear *)
      Graphics`PolygonUtils`InPolygonQ[Polygon[poly], Mean[pts]],
     point = Mean[pts];
     Break[]]],
   {t, Length[poly] - 2}];
  If[VectorQ[point],
   point,
   Message[findPtInPoly::nopt, poly];
   First[poly]]]

SeedRandom[0];
npts = 300;
tmp = RandomReal[100, {npts, 2}];
plot = ListLinePlot[tmp, Filling -> Axis]
With[{poly = First@Cases[Normal@plot, _Polygon, Infinity]},
  Total[If[Graphics`PolygonUtils`InPolygonQ[poly, findPtInPoly[#]],
      Graphics`PolygonUtils`PolygonArea[#], 0] & /@
    Graphics`PolygonUtils`SimplePolygonPartition[poly]]
  ] // AbsoluteTiming
Area@DiscretizeGraphics[plot] // AbsoluteTiming
(*
  {2.931620, 4368.75}
  {32.840698, 4368.75}
*)

It seems to have better time complexity than Area @* DiscretizeGraphics. On inputs of size npts equal to 100, 200, 300, the timing ratios of the two methods are {1.8112, 4.16505, 11.2022}.

$\endgroup$
9
  • $\begingroup$ Interesting trick. Unfortunately, it is really slow, there is no way to run this code on big polylines (thousands of points)... I think that there is some obvious analytical solution with computing areas of trapeze. $\endgroup$
    – user25235
    Feb 7, 2015 at 20:52
  • $\begingroup$ Trying to debug your new code. It gives a wagon of errors, first being "First::normal: Nonatomic expression expected at position 1 in First[poly]. >>". BTW, in what version of mathematica did you run it? $\endgroup$
    – user25235
    Feb 7, 2015 at 22:07
  • $\begingroup$ Figured what you did. A really cunning approach! Works fast and nice, thank you. In case somebody will try to reproduce this, step by step version ("With" expression did not work for me): p = First[Cases[Normal@plot, _Polygon, Infinity]] g = GraphicsPolygonUtilsSimplePolygonPartition[p] Total[GraphicsPolygonUtilsPolygonArea /@ Pick[g, GraphicsPolygonUtilsInPolygonQ[p, Mean@*First /@ g]]] $\endgroup$
    – user25235
    Feb 7, 2015 at 22:21
  • $\begingroup$ Ummm... your npts is used for the range of the reals, not for the number of them. You always use 100 reals. When I try 300 instead of 100 (points) it already takes 3.5 seconds for me. Also, our answers seem to disagree quite a bit for larger numbers of points. E.g. for 200 points, seed 1, yours gets 4164.96 while mine gets 4553.61 (yours is a fair bit faster though and also seems to scale better than mine). $\endgroup$ Feb 7, 2015 at 22:43
  • $\begingroup$ @MichaelE2 I guess, you'd have to implement an exact solution manually to get a reference result. I also noticed that your code can return a negative result, while mine will always report a positive area. (Just try tmp = {{1, 1}, {-1, 1}}.) I don't know which is desired in that case. Also, mine won't really handle more than 6 (or so) integer coordinates (but doesn't have any trouble with large numbers of reals). $\endgroup$ Feb 8, 2015 at 1:29
0
$\begingroup$
tmp = {{0, 0}, {1, 1}, {2, 1}, {3, 2}, {1, 0.5}}; 
Graphics@ Polygon[tmp]

yields this:

enter image description here

Its area is

Area[Polygon[tmp]]

(*  0.583333   *)

Have fun!

$\endgroup$
1
  • $\begingroup$ Thank you for the answer but indeed, i need something different. Right now i am trying code like this: N@Total[Partition[tmp,2,1] /.{{x1_, y1_}, {x2_, y2_}} -> (x2 - x1)*(y1 + (y2 - y1)/2)] The key is to somehow deal with overlaps $\endgroup$
    – user25235
    Feb 7, 2015 at 21:00

Your Answer

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

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