# Calculate area under a polyline

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

• ComponentMeasurements[ColorNegate@Binarize@ListLinePlot[tmp, Filling -> Axis, Axes -> False], "Area"] Feb 7, 2015 at 20:23

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 GraphicsPolygonUtilsSimplePolygonPartition. 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 *)
GraphicsPolygonUtilsInPolygonQ[Polygon[poly], Mean[pts]],
point = Mean[pts];
Break[]]],
{t, Length[poly] - 2}];
If[VectorQ[point],
point,
Message[findPtInPoly::nopt, poly];
First[poly]]]

SeedRandom;
npts = 300;
tmp = RandomReal[100, {npts, 2}];
plot = ListLinePlot[tmp, Filling -> Axis]
With[{poly = First@Cases[Normal@plot, _Polygon, Infinity]},
Total[If[GraphicsPolygonUtilsInPolygonQ[poly, findPtInPoly[#]],
GraphicsPolygonUtilsPolygonArea[#], 0] & /@
GraphicsPolygonUtilsSimplePolygonPartition[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}.

• 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. Feb 7, 2015 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? Feb 7, 2015 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]]] Feb 7, 2015 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). Feb 7, 2015 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). Feb 8, 2015 at 1:29
tmp = {{0, 0}, {1, 1}, {2, 1}, {3, 2}, {1, 0.5}};
Graphics@ Polygon[tmp]


yields this: Its area is

Area[Polygon[tmp]]

(*  0.583333   *)


Have fun!

• 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 Feb 7, 2015 at 21:00