5

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

user25235
  • 53
  • 3
  • ComponentMeasurements[ColorNegate@Binarize@ListLinePlot[tmp, Filling -> Axis, Axes -> False], "Area"] – Sektor Feb 07 '15 at 20:23

2 Answers2

7
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}.

Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • 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. – user25235 Feb 07 '15 at 20:52
  • 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? – user25235 Feb 07 '15 at 22:07
  • 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]]] – user25235 Feb 07 '15 at 22:21
  • 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). – Martin Ender Feb 07 '15 at 22:43
  • @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). – Martin Ender Feb 08 '15 at 01:29
  • @user25235 I think I figured out how to use the PolygonUtils functions. The time complexity is still worse than n^2, but I supposed that's expected when the computation depends on the intersections of the lines. – Michael E2 Feb 09 '15 at 01:50
  • @MichaelE2 I just tested Area@Polygon on a very simple self-intersecting case, and it doesn't actually compute the area that's rendered by Graphics@Polygon. Instead it seems to double-count areas which are covered multiple times. I'm deleting my answer. – Martin Ender Feb 09 '15 at 11:42
  • @MichaelE2 If I change the code into tmp = {{0, 0}, {1, 1}, {2, 1}, {3, 2}, {1, 6}}; plot = ListLinePlot[tmp, Filling -> Axis]; Area@DiscretizeGraphics[plot] I got the answer is 8. My maple has the same result. But your code is 6. I don't understand why? – minthao_2011 Feb 09 '15 at 15:08
  • @minthao_2011 I get 6, if I add Gridlines and just count. – Michael E2 Feb 09 '15 at 15:17
0
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!

Alexei Boulbitch
  • 39,397
  • 2
  • 47
  • 96
  • 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 – user25235 Feb 07 '15 at 21:00