Compute area between two ListPlot joined list (area filled by Filling -> {1 -> {2}})

Question

Given two lists, how can I compute the area between their plots? (the blue region in the image below).

The horizontal distance between points in the list can be assumed to be 1, as done by ListPlot automatically in the example image.

For my use case, performance is really important so I would prefer a method that does not require doing any plotting.

Edit, solution chosen:

I finally went with this compilated version of the accepted version, which is very fast, at least on my tests:

BlocksAreaC = Compile[{{x, _Real, 2}},
   Block[
    {
     diff = x[[2]] - x[[1]],
     c = (-x[[1, 1]] + x[[2, 1]])
     },
    If[
     Times @@ Sign[diff] == 1,
     Total[Abs[diff]]/2,
     If[c == 0, 0,
      With[
       {h = c/(-x[[1, 1]] + x[[1, 2]] + x[[2, 1]] - x[[2, 2]])},
       Total[Abs[diff]*{h, 1 - h}/2]
       ]
      ]
     ]
    ],
   Parallelization -> True, RuntimeAttributes -> {Listable}, 
   CompilationTarget -> "C", RuntimeOptions -> "Speed", 
   CompilationOptions -> {"ExpressionOptimization" -> True, 
     "InlineCompiledFunctions" -> True, 
     "InlineExternalDefinitions" -> True}
   ];
areaBetweenLists7C = Compile[{{l1, _Real, 1}, {l2, _Real, 1}},
   Module[
    {
     s1 = Partition[l1, 2, 1],
     s2 = Partition[l2, 2, 1]
     },
    Total[Map[BlocksAreaC, Transpose[{s1, s2}]]]
    ],
   Parallelization -> True, RuntimeAttributes -> {Listable}, 
   CompilationTarget -> "C", RuntimeOptions -> "Speed", 
   CompilationOptions -> {"ExpressionOptimization" -> True, 
     "InlineCompiledFunctions" -> True, 
     "InlineExternalDefinitions" -> True}
   ];

Answer 1

One could also go with a fully analytical approach:

l1 = {1, 2, 4, 3, 1};
l2 = {3, 6, 5, 1, 0};

s1 = Subsequences[l1, {2}];
s2 = Subsequences[l2, {2}];
s = Transpose[{s1, s2}];

Edit: to avoid the code breaking when a polygon has area = 0, one can replace s with:

s = Select[Transpose[{s1, s2}], #[[1]] != #[[2]] &]

To speed up computation, one might want to add this check in a nested if within the BlocksArea function, but it would need some testing for checking which solution is the fastest.

BlocksArea[x_] := Block[{diff = x[[2]] - x[[1]], h},
  If[Times @@ Sign[diff] == 1,
   Total[Abs[diff]]/2,
   h = (-x[[1, 1]] + x[[2, 1]])/(-x[[1, 1]] + x[[1, 2]] + x[[2, 1]] - x[[2, 2]]);
   Total[Abs[diff]*{h, 1 - h}/2]
   ]
  ]

Total[BlocksArea /@ s]

47/6

what the code does is computing the area of the regions between 3 point piece by piece: if the points form trapezoid shape the area is computed as the sum of the bases, times the height, divided by 2. if the points form two triangles, one can find the height of the triangles and then the corresponding areas.

It looks like the code in my answer runs faster than the ones above (see a non-exhaustive benchmark below), and I'm pretty sure that one can speed up my code by rewriting it (keeping the same underline concept) for efficiently compiling it to C.

I think it's also suitable for parallelisation, as the lists of points can be split in the number of cores/kernel available for computing the area of each piece on a different thread.

Here a speed comparison:

kglr, code in point 3:

RepeatedTiming[
 {if1, if2} = 
  Interpolation[#, InterpolationOrder -> 1] & /@ 
   lists;
 NIntegrate[Abs[if1[t] - if2[t]], {t, 1, 5}]
 ]

{0.011, 7.83333}

J. M. is in limbo♦ code:

RepeatedTiming[
 Integrate[Abs[Apply[Subtract, makePW[#, t] & /@ lists]], {t, 1, 5}]
 ]

{0.020, 47/6}

my code:

RepeatedTiming[s1 = Subsequences[l1, {2}];
 s2 = Subsequences[l2, {2}];
 s = Transpose[{s1, s2}];
 Total[BlocksArea /@ s]]

{0.00033, 47/6}

EDIT: Difference when comparing two long lists (10000 samples)

generating the lists:

lists = {l1, l2} = RandomInteger[{0, 10}, {2, 10000}];

output of the timing (in the same order as above):

kglr: 0.66

J. M. is in limbo♦: 4.4

Fraccalo: 0.19

Answer 2

1. You can extract the polygons in llp using Cases:

llp = ListLinePlot[{{1, 2, 4, 3, 1}, {3, 6, 5, 1, 0}}, Filling -> {1 -> {2}}];

polygons = Cases[Normal@llp, _Polygon, All]

{Polygon[{{1.,1.},{2.,2.},{3.,4.},{3.33333,3.66667},{3.33333,3.66667},{3.,5.},{2.,6.},{1.,3.}}],
Polygon[{{3.33333,3.66667},{4.,1.},{5.,0.},{5.,1.},{4.,3.},{3.33333,3.66667}}]}

Area /@ polygons

{5.66667, 2.16667}

Total @ %

7.83334

2. You can use BoundaryDiscretizeGraphics to get a MeshRegion and get its Area or RegionMeasure:

Area @ BoundaryDiscretizeGraphics @ llp

7.833333333333334

RegionMeasure @ BoundaryDiscretizeGraphics[llp]

7.833333333333334

3. Use Interpolation on the two lists to get two functions and NIntegrate to get area between the two:

lists = {{1, 2, 4, 3, 1}, {3, 6, 5, 1, 0}};
{if1, if2} = Interpolation[#, InterpolationOrder -> 1] & /@ lists; (*thanks: Mr.Wizard*)

NIntegrate[Abs[if1[t] - if2[t]], {t, 1, 5}]

7.833333171170218

Answer 3

In fact, you can compute the exact answer for this case if you explicitly assemble the piecewise linear function representing the two "connect-the-dots" plots in the OP, and then feed the integrand to Integrate[]. Here's one way to derive the required piecewise linear interpolant:

makePW[ya_?VectorQ, t_] := 
    Piecewise[MapIndexed[{InterpolatingPolynomial[{{#2, #1[[1]]}, {#2 + 1, #1[[2]]}}, t],
                          First[#2] <= t <= First[#2] + 1} &, Partition[ya, 2, 1]]]

Then,

lists = {{1, 2, 4, 3, 1}, {3, 6, 5, 1, 0}};

Integrate[Abs[Apply[Subtract, makePW[#, t] & /@ lists]], {t, 1, 5}]
   47/6

---

One can also reformulate an equivalent procedure in terms of UnitStep[]:

makePW2[ya_?VectorQ, t_] := 
    With[{lp = Most[ya] + Differences[ya] (t - Range[Length[ya] - 1])}, 
         First[lp] + Differences[lp].UnitStep[t - Range[2, Length[lp]]]]

and Integrate[Abs[Apply[Subtract, makePW2[#, t] & /@ lists]], {t, 1, 5}] should yield the same answer.

Answer 4

EDIT

This approach doesn't work, as the method of choosing points to include in the polygon (counting crossings of rays extending from point) results in points sometimes in, sometimes out. With enough points, you are guaranteed for it to fail.

Keeping this answer as it has useful info, and to head off any future ventures down this path.

---

This is easy to state

upper = Transpose@{Range[5], {1, 2, 4, 3, 1}}
lower = Transpose@{Range[5], {3, 6, 5, 1, 0}}

polyPts = Join@@{upper,Reverse@lower};
ptOrder = {1,2,3,4,5,6,7,8,9,10,1};

poly = Polygon[polyPts, ptOrder];
Graphics[poly]

Integrate[1, {x, y} ∈ poly]
(* 7.83333 *)

Note: works in Mma 12 but not in 11, it appears.