polyglint2d: Integrating a 2D linear function over convex polygonal bounds

The package provides functions which compute integrals of the form $$\int_{Ax \leq b} ( c^{\rm T} x + d ) \, {\rm d}x,$$ where $c^{\rm T} x + d$ is an arbitrary two-dimensional linear function, i.e., $c \in {\mathbb R}^2$, $d \in {\mathbb R}$, and the region of integration is an arbitrary closed, convex, two-dimensional polygon, i.e., $A x \leq b$ for $A \in {\mathbb R}^{n \times 2}$, $b \in {\mathbb R}^n$.

Equivalently, given a set $P$ of points in ${\mathbb R}^2$, we can compute the integral $$\int_{{\rm ConvHull}(P)} ( c^{\rm T} x + d ) \, {\rm d}x,$$ where ${\rm ConvHull}(P)$ is the convex hull of $P$.

Computing the area of a shape

The simplest example is to compute the area of a polytope, which we get when $c = {\bf 0}$ and $d = 1$.

For example, suppose we wish to compute the area of the following shape.

# Define the corners.
x1, x2 = .3, .7
points = ((0, 0), (1, 0), (x2, 1), (x1, 1))

import scipy as sp
import scipy.spatial

# Get the convex hull.
hull = sp.spatial.ConvexHull(points)

# Display it.
%matplotlib inline

import matplotlib.pyplot as plt
import numpy as np

points_ = np.array(points)

plt.plot(points_[:, 0], points_[:, 1], 'o')  # show original points
plt.fill(points_[hull.vertices, 0], points_[hull.vertices, 1],
         color='lightgreen', edgecolor='black', alpha=0.5)
plt.xlabel("X"); plt.ylabel("Y")
plt.gca().set_aspect('equal'); plt.grid(True)
plt.show()

A common approach is to split the shape into two triangles and a rectangle, and compute the area of each separately.

shapes = (
    ((0, 0), (x1, 0), (x1, 1)),
    ((x1, 0), (x2, 0), (x2, 1), (x1, 1)),
    ((x2, 0), (1, 0), (x2, 1)),
)

for ps, color in zip(shapes, ("red", "green", "blue")):
    points_ = np.array(ps)
    plt.plot(points_[:, 0], points_[:, 1], 'o', c="k")  # show original points
    plt.fill(points_[:, 0], points_[:, 1],
         color=color, edgecolor='black', alpha=0.4)
plt.xlabel("X"); plt.ylabel("Y")
plt.gca().set_aspect('equal'); plt.grid(True)

plt.annotate("$\\frac{.3}{2}$", (.2, .4), ha="center", va="center", fontsize=20)
plt.annotate("$(.7-.3)$", (.5, .4), ha="center", va="center", fontsize=13)
plt.annotate("$\\frac{.3}{2}$", (.8, .4), ha="center", va="center", fontsize=20)
plt.annotate("$\\Sigma A = .7$", (.9, .9), ha="center", va="center", fontsize=16)
plt.show()

Indeed this is more or less what happens under the hood when we make the following call:

from setiptah.polyglint2d import *

c, d = np.zeros(2), 1.  # f(x) = 1
area = integrate_over_convexhull(c, d, points)

# print(points)
print(f"The area of our trapezoid is: {area:.4f}")

assert abs(area - .7) < 1e-12

The area of our trapezoid is: 0.7000

We can obtain an identical result using an inequality form of the same region:

def below_line(x1, x2, y1, y2):
    m = (y2 - y1) / (x2 - x1)
    b = y1 - m * x1
    return [-m, 1, b]

Ab = np.array([
    [0, -1, 0],  # y >= 0
    [0, 1, 1],  # y <= 1
    below_line(0., .3, 0., 1.),  # to the right of the left edge
    below_line(.7, 1., 1., 0.),  # to the left of the right edge
])
A, b = Ab[:, :-1], Ab[:, -1]

def matrix_to_latex(A, env="bmatrix", fmt="{}"):
    rows = [" & ".join(fmt.format(val) for val in row) for row in A]
    body = " \\\\\n".join(rows)
    return f"\\begin{{{env}}}\n{body}\n\\end{{{env}}}"

def array_to_latex(vec, env="bmatrix", fmt="{}"):
    rows = [fmt.format(val) for val in vec]
    body = " \\\\\n".join(rows)
    return f"\\begin{{{env}}}\n{body}\n\\end{{{env}}}"

from IPython.display import display, Latex
array_to_latex(b)
display(Latex(f'${{\\rm Region}} = {matrix_to_latex(A)} {array_to_latex(["x", "y"])} \\leq {array_to_latex(b)}$'))

${\rm Region} = \begin{bmatrix} 0.0 & -1.0 \\ 0.0 & 1.0 \\ -3.3333333333333335 & 1.0 \\ 3.333333333333333 & 1.0 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} \leq \begin{bmatrix} 0.0 \\ 1.0 \\ 0.0 \\ 3.333333333333333 \end{bmatrix}$

area = integrate(c, d, A, b)
display(Latex(f"Area of the described region: {area}"))

assert abs(area - .7) < 1e-12

Area of the described region: 0.7000000000000001

General linear integration

The area of a two-dimensional region $Ax \leq b$ is also equal to the volume of the three-dimensional region contained by the cylinder $Ax \leq b$ and between the $z=0$ and $z=1$ planes.

(Feel free to skip the plotting code.)

from mpl_toolkits.mplot3d import Axes3D  # sometimes needed for 3D
from mpl_toolkits.mplot3d import art3d
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
from IPython.display import Image

from dataclasses import dataclass
from typing import Any

from setiptah.polyglint2d._geom import *

@dataclass(frozen=True)
class SurfaceArt:
    """This is how we will draw surfaces on a plot."""
    points: Any
    faces: Any
    face_color: str
    face_alpha: float = .3
    wire_color: str = "k"
    line_style: str = "-"

    def draw(self, ax, coords = None):
        if coords is not None:
            coords.extend(self.points)
        
        # Compute face polygons
        try:
            face_polys = [[self.points[i] for i in face] for face in self.faces]
        except:
            print(len(self.points), self.faces)
            raise
        
        # Fill layer
        face_collection = Poly3DCollection(
            face_polys, facecolors=self.face_color, edgecolors='none', alpha=self.face_alpha
        )
        ax.add_collection3d(face_collection)
        
        # Wireframe layer
        edges = set()
        for face in self.faces:
            for i in range(len(face)):
                edge = tuple(sorted((face[i], face[(i + 1) % len(face)])))
                edges.add(edge)
        edge_lines = [[self.points[i], self.points[j]] for i, j in edges]
    
        if self.wire_color:
            edge_collection = Line3DCollection(edge_lines, colors=self.wire_color, linewidths=1)
            ax.add_collection3d(edge_collection)

def trapezoid_surfaces(points, c, d):
    trapezoid_dict = convex2d_to_traps(points)

    # Choose a colormap (e.g., 'viridis', 'plasma', 'tab10', etc.)
    cmap = plt.cm.viridis  # or plt.get_cmap('viridis')
    face_colors = [cmap(t) for t in np.linspace(0, 1, len(trapezoid_dict))]

    for trapezoid, face_color in zip(trapezoid_dict.items(), face_colors):
        (x1, x2), ((mtop, btop), (mbot, bbot)) = trapezoid
        traversal = traverse_trapezoid(x1, x2, mtop, btop, mbot, bbot)
        points_ = compute_points(traversal, c, d)
        faces = compute_faces(len(traversal))
        yield SurfaceArt(points_, faces, face_color=face_color, wire_color="black")

def bounding_surfaces(points, c, d):
    M = np.array(points)

    def scale_from_center(arr, scale):
        mu = np.mean(arr)
        return mu + scale * (arr - mu)

    xmin, xmax = scale_from_center(extents(M[:, 0]), 1.2)
    ymin, ymax = scale_from_center(extents(M[:, 1]), 1.2)
    
    window_points2 = [
        (xmin, ymin),
        (xmax, ymin),
        (xmax, ymax),
        (xmin, ymax),
    ]

    # planar surfaces
    cx, cy = c
    top = [
        (x, y, cx * x + cy * y + d)
        for x, y in window_points2
    ]
    ground = [
        (x, y, 0)
        for x, y in window_points2
    ]
    plane_faces = [(0, 1, 2, 3)]

    plane_opts = dict(wire_color=None, face_alpha=.2)
    yield SurfaceArt(ground, plane_faces, "g", **plane_opts)
    yield SurfaceArt(top, plane_faces, "r", **plane_opts)

    # cyclinder surfaces
    zs = np.array([
        z
        for points in (top, ground)
        for x, y, z in points
    ])
    zmin, zmax = scale_from_center(extents(zs), 1.2)
    
    cylinder = [
        (x, y, z)
        for x, y in points
        for z in (zmin, zmax)
    ]    
    cylinder_faces = compute_faces(len(points))
    yield SurfaceArt(cylinder, cylinder_faces, "b", wire_color="gray", face_alpha=.1)

def all_surfaces(points, c, d):
    yield from bounding_surfaces(points, c, d)
    yield from trapezoid_surfaces(points, c, d)

def traverse_trapezoid(x1, x2, mtop, btop, mbot, bbot):
    top = np.array([mtop, btop])
    bot = np.array([mbot, bbot])

    return [
        (x1, np.polyval(bot, x1)),
        (x2, np.polyval(bot, x2)),
        (x2, np.polyval(top, x2)),
        (x1, np.polyval(top, x1)),
    ]

def compute_points(traversal2, c, d):
    cx, cy = c
    return [
        (x, y, z)
        for x, y in traversal2
        for z in (0, cx * x + cy * y + d)
    ]

def compute_faces(n_traversal):
    ring_bottom = 2 * np.arange(n_traversal)
    ring_top = ring_bottom + 1
    cylinder_faces = [
        ring_bottom,
        ring_top,
    ]
    for k in ring_bottom:
        path = (k, k+2, k+3, k+1)
        path = tuple(k % (2 * n_traversal) for k in path)
        cylinder_faces.append(path)
    return cylinder_faces

def extents(arr):
    return arr.min(), arr.max()

plt.close("all")
fig = plt.figure(figsize=(8, 6))
for k, (elev, azimuth) in enumerate([
    (70, -90), (20, -90), (20, -45), (20, -135)
]):
    plt.subplot(2,2, k+1, projection='3d')
    ax = plt.gca()

    ax.view_init(elev=elev, azim=azimuth)  # camera
    
    all_coords = []
    
    # Draw.
    for surf in all_surfaces(points, c, d):
        surf.draw(ax, all_coords)
    
    # ax.auto_scale_xyz(all_coords, all_coords, all_coords)
    all_coords = np.array(all_coords)
    ax.set_xlim(*extents(all_coords[:, 0]))
    ax.set_ylim(*extents(all_coords[:, 1]))
    ax.set_zlim(*extents(all_coords[:, 2]))
    ax.set_aspect("equal")
    # plt.tight_layout()
    
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_zlabel('$c^{\\rm T} x + d$')
    
fig.suptitle("A trapezoidal decomposition of the bounded region (with bounding planes);\n different viewpoints");

The real magic happens though when $c'x + d$ is not a constant function like area, but a general linear function in 2D.

Suppose now we are integrating over a region shown below.

ts = np.linspace(0, 2*np.pi, 10)[:-1]
points = np.vstack([np.cos(ts), np.sin(ts)]).T
points = [list(row) for row in points]

points_ = np.array(points)

plt.plot(points_[:, 0], points_[:, 1], 'o')  # show original points
plt.fill(points_[:, 0], points_[:, 1],
         color='lightgreen', edgecolor='black', alpha=0.5)
plt.xlabel("$x_1$"); plt.ylabel("$x_2$")
plt.gca().set_aspect('equal'); plt.grid(True)
plt.show()

Let's demonstrate the integral $\int_{{\rm ConvHull}(P)} ( \frac{1}{2} x_1 + \frac{1}{2} x_2 + \frac{1}{10} ) \, {\rm d}x$ over the shown region:

def trapezoid_surfaces(points, c, d):
    """

    We need to refine our method for trapezoid decomposition. The function plane could (does!) cross the z=0.

    """
    from setiptah.polyglint2d import vertices_to_hull_inequality
    
    A, b = vertices_to_hull_inequality(points)
    
    A_below = np.vstack((A, c))
    b_below = np.concatenate((b, [-d]))
    hull_below = enumerate_vertices_2d(A_below, b_below)

    A_above = np.vstack((A, -c))
    b_above = np.concatenate((b, [d]))
    hull_above = enumerate_vertices_2d(A_above, b_above)

    parts = [
        (trapezoid, wire_color)
        for hull, wire_color in (
            (hull_above, "black"),
            (hull_below, "red"),
        )
        if len(hull) > 0
        for trapezoid in convex2d_to_traps(hull).items()
    ]

    # Choose a colormap (e.g., 'viridis', 'plasma', 'tab10', etc.)
    cmap = plt.cm.viridis  # or plt.get_cmap('viridis')
    face_colors = [cmap(t) for t in np.linspace(0, 1, len(parts))]

    for (trapezoid, wire_color), face_color in zip(parts, face_colors):
        (x1, x2), ((mtop, btop), (mbot, bbot)) = trapezoid
        traversal = traverse_trapezoid(x1, x2, mtop, btop, mbot, bbot)
        points_ = compute_points(traversal, c, d)
        faces = compute_faces(len(traversal))
        yield SurfaceArt(points_, faces, face_color=face_color, wire_color=wire_color)

c, d = np.array([.5, .5]), 0.1


plt.close("all")
fig = plt.figure(figsize=(8, 6))

for k, (elev, azimuth) in enumerate([
    (0, -60), (20, -80), (20, -45), (20, -135)
]):
    plt.subplot(2,2, k+1, projection='3d')
    ax = plt.gca()

    ax.view_init(elev=elev, azim=azimuth)  # camera

    all_coords = []

    # Draw.
    for surf in all_surfaces(points, c, d):
        surf.draw(ax, all_coords)
    
    # ax.auto_scale_xyz(all_coords, all_coords, all_coords)
    all_coords = np.array(all_coords)
    ax.set_xlim(*extents(all_coords[:, 0]))
    ax.set_ylim(*extents(all_coords[:, 1]))
    ax.set_zlim(*extents(all_coords[:, 2]))
    ax.set_aspect("equal")
    # plt.tight_layout()
    
    ax.set_xlabel('$x_1$')
    ax.set_ylabel('$x_2$')
    ax.set_zlabel('$c^{\\rm T} x + d$')
    
fig.suptitle("A trapezoidal decomposition of the bounded region (with bounding surfaces);\n different viewpoints");

Because the function intersects the $z=0$ plane within the cylinder, then the integral is the difference of two volumes: It is the volume of the region with the black wire frame---lying above the $z=0$ plane---minus the volume of the region with the red wire frame.

Since we added a small positive constant $d=\frac{1}{10}$, we should expect a slightly positive net result.

ans = integrate_over_convexhull(c, d, points)

print(f"The value of the integral is: {ans:.4f}")

The value of the integral is: 0.2893