🌐 An Efficient Algorithm for Computing the Area of Composite Spherical Regions

Zhou Renfei, Hangzhou Xuejun High School

📄 Abstract

On a sphere, by repeatedly performing Boolean operations such as intersection, union, and difference on spherical cap regions, many complex regions can be represented; these are referred to in this paper as "composite regions". Computing the area of a spherical composite region is an important problem with wide applications in reality. This paper proposes a precise definition of the composite region area problem given in the form of a Boolean expression and presents an excellent algorithm. This algorithm calculates the exact value of the area of the spherical composite region with a worst-case running time of $O(nm \log m)$, where $n$ is the number of spherical cap regions and $m$ is the length of the Boolean expression.

1. 🚀 Introduction

Given $n$ spherical cap regions $D_1, D_2, \dots, D_n$ on a sphere and a Boolean expression of length $m$ using them as parameters, a composite region $D^*$ is defined. Here, the spherical caps are distinct; the Boolean expression is given explicitly, and each region can appear multiple times in the expression. The main issue discussed in this paper is how to calculate the area of this composite region—the Spherical Composite Region Area Problem.

This type of problem has wide applications in reality, such as:

Calculating the raw material usage for spherical crafts.
Calculating the area of specific geographical regions.
Calculating the coverage of communication base stations.
Calculating the effective service range of the "BeiDou" satellite navigation system on the Earth's surface.

These applications are particularly important in Computer-Aided Design (CAD) and Geographic Information Systems (GIS).

A special case of this problem on the plane has long appeared in algorithm competitions. As early as NOI 2004, the problem "Rainfall" appeared, requiring the calculation of the area of the union of several parallelograms on a plane. In subsequent years, "Area Union" problems on the plane appeared several times in NOI, such as the famous NOI 2005 "Moonlit Lemon Tree". As for the spherical composite region area problem, it appeared once in an ICPC contest.

Previously, there was no general method for solving such problems. There is a traditional numerical integration method that can calculate an approximate value of the composite region's area. Additionally, there is a scan-line method that calculates the exact value of the plane composite region area within a worst-case time complexity of $O(n^2 m \log m)$ by cutting the composite region. This paper proposes an excellent algorithm that calculates the exact area of planar or spherical composite regions within a worst-case time complexity of $O(nm \log m)$, featuring high efficiency, high precision, good numerical stability, and strong extensibility.

2. 🟢 Plane Area Union Problem

Problem 1 (Plane Area Union Problem): Given $n$ circular regions $D_1, \dots, D_n$ on a plane, define $D^* = D_1 \cup D_2 \cup \dots \cup D_n$, where each circular region $D_i$ refers to the interior of a circle. Calculate the area of $D^*$.

This chapter proposes an algorithm based on Green's formula to solve the above problem.

Green's Formula: Let $\Omega \subset \mathbb{R}^2$ be a closed region bounded by a finite number of piecewise smooth curves. If functions $P(x,y)$ and $Q(x,y)$ are continuous on $\Omega$ and have continuous partial derivatives $\partial Q / \partial x$ and $\partial P / \partial y$, then: $$ \oint_{\partial \Omega} P dx + Q dy = \iint_{\Omega} \left(\frac{\partial Q}{\partial x} - \frac{\partial P}{\partial y}\right) dx~dy, $$ where $\partial \Omega$ is the boundary of region $\Omega$. Its orientation is determined such that when a person walks along the positive direction of $\partial \Omega$, the region is always on their left.

Let $S_{D^*}$ denote the area of $D^*$, then $S_{D^*} = \iint_{D^*} dx~dy$. If we let $\Omega = D^*$, $P(x,y) = -y$, and $Q(x,y) = 0$, substituting into Green's formula yields: $$ \oint_{\partial D^{\ast}} (-y) dx = \iint_{D^{\ast}} dx~dy = S_{D^{\ast}} $$ Thus, calculating the area is transformed into calculating the line integral of the boundary $\partial D^*$.

As shown in Figure 1, $\partial D^*$ consists of several circular arcs $A_1, \dots, A_m$. According to the additivity of integrals: $$ \oint_{\partial D^*} (-y) dx = \sum_{i=1}^{m} \int_{A_i} (-y) dx, $$ where the direction of each $A_i$ is counter-clockwise. Thus, Problem 1 is transformed into two sub-problems:

Problem 1.1: Find $A_1, A_2, \dots, A_m$.
Problem 1.2: Calculate $\int_{A_i} (-y) dx$.

🔍 Discussing Problem 1.1

Without loss of generality, assume circles $D_1, D_2, \dots, D_n$ are distinct. To find the arcs constituting $\partial D^*$, we must calculate the parts of the circumference not covered by other circles for each circle respectively.

Define the "Polar Angle" operation: for a point $P$ on the circumference with center $O$, the polar angle is the angle rotated from the positive x-axis to the vector $OP$, with a range of $[-\pi, \pi)$. The problem is transformed into calculating the uncovered parts on $[-\pi, \pi)$.

Problem 1.1a: Given $n'$ intervals on $[-\pi, \pi)$, find the parts not covered by any interval, expressed as several disjoint intervals. Using the scan-line idea, let a moving point $P$ move from $-\pi$ to $\pi$, maintaining how many intervals contain $P$. This solves the problem with a running time upper bound of $O(n' \log n')$. Performing this for all circles solves Problem 1.1 in a worst-case time complexity of $O(n^2 \log n)$.

[Image of angular sweep line algorithm on a circle]

📐 Discussing Problem 1.2

Suppose we want to calculate $\int_{A_i} (-y) dx$ for an arc $A_i = \text{arc } PMQ$. Let $\Omega_1$ be the segment (bow shape) bounded by $A_i$ and vector $\vec{QP}$, and $\Omega_2$ be the region bounded by $\vec{QP}$, the x-axis, and vertical lines. Applying Green's formula to $\Omega_1$: $$ \int_{A_i} (-y) dx = S_{\Omega_1} + S_{\Omega_2} $$ where $S_{\Omega_2}$ is the signed area of region $\Omega_2$. Using this formula, the integral of an arc $A_i$ can be calculated in constant time.

Summary: We obtain an excellent algorithm solving Problem 1 with a running time upper bound of $O(n^2 \log n)$.

3. 🌍 Spherical Area Union Problem

Problem 2 (Spherical Area Union Problem): Given $n$ spherical cap regions $D_1, \dots, D_n$ on a sphere (caps not larger than a hemisphere), define $D^* = D_1 \cup \dots \cup D_n$. Calculate the area of $D^*$.

Basic Knowledge of Spherical Geometry:

Assume the sphere is a unit sphere, center $O=(0,0,0)$, radius 1.
A plane intersecting the sphere forms a circle; the part of the sphere on one side is a spherical cap.
Two points are antipodal if their connecting line passes through the sphere's center.
A "Great Circle" passes through the center; other circles are "Small Circles". The shortest path between non-antipodal points $A, B$ is the spherical line segment $AB$.
Spherical coordinates: $x = \cos \varphi \cdot \cos \lambda$, $y = \cos \varphi \cdot \sin \lambda$, $z = \sin \varphi$, where $\lambda \in [0, 2\pi)$ and $\varphi \in [-\pi/2, \pi/2]$.

Theorem 1: Assume $\Omega$ is a region on the sphere bounded by finite piecewise smooth curves, and the North Pole $P_N$ is not on the boundary of $\Omega$. Then: $$ S_{\Omega} = \oint_{\partial \Omega} (-\sin \varphi - 1) d\lambda + [P_N \in \Omega]4\pi. $$

Proof: Using the area element $da = \cos \varphi \cdot d\lambda \cdot d\varphi$ and mapping to the "Lat-Lon Plane" $\hat{\Omega}$ , we apply Green's formula with $P = -\sin \varphi - 1$ and $Q = 0$. The boundary integral accounts for the domain boundary ($Part-1$) and the region boundary $\partial \Omega$ ($Part-2$).

Corollary: $S_{\Omega}$ and $\oint_{\partial \Omega} (-\sin \varphi - 1) d\lambda$ can be derived from each other in constant time.

Solving Problem 2: First, randomly rotate the sphere so no boundary passes through the North Pole. The problem converts to calculating $\oint_{\partial D^*} (-\sin \varphi - 1) d\lambda$.

Problem 2.1: Find arcs of each circle not covered by others using the same angular sweep-line method as Chapter 2.
Problem 2.2: Calculate $\int_{C} (-\sin \varphi - 1) d\lambda$ for a directed arc $C$.
- Case 1: $C$ is part of a small circle. Form a closed region $A$ with chord $C_2$. Using the Corollary, relate the integral to the area of the spherical segment $A$ and the integral of the chord. The area of $A$ is derived from the area of a spherical sector ($S = 2\pi h \theta$) and a spherical triangle ($S_{\triangle ABC} = \angle A + \angle B + \angle C - \pi$).
- Case 2: $C$ is part of a great circle. Reduce to spherical line segments.

The structure is similar to Chapter 2, with time complexity $O(n^2 \log n)$.

4. 🔗 Boolean Operations

We define a composite region $D^*$ using an $n$-ary Boolean function $F(\lambda_1, \dots, \lambda_n)$ where $\lambda_i = [P \in D_i]$. $$ D^* = { P \mid F(\lambda_1, \dots, \lambda_n) = 1 } $$

Problem 3 (Spherical Composite Region Area Problem): Given spherical caps $D_1, \dots, D_n$ and a Boolean expression $\alpha_F$ of length $m$, calculate the area of $D^*$.

We need to calculate the boundary $\partial D^*$. An arc on $\partial D_i$ is part of $\partial D^*$ if its left side is in $D^*$ and right side is out of $D^*$ (or vice-versa for the reverse direction). This leads to Problem 3.1.

Problem 3.1: For each circle, map boundary to $[-\pi, \pi)$. Other circles cover intervals where $\lambda_j=1$. Find angular intervals where $F$ evaluates to 1. Using a scan-line, as point $P$ moves, there are $O(n)$ "events" where a $\lambda_j$ changes value.

Expression Tree Data Structure: Build a tree $T(\alpha_F)$ representing the expression.

Problem 3.2: Design a data structure for INITIALIZE-TREE, MODIFY-TREE, QUERY-TREE.
Algorithm 1 (Naive): Recompute values up to the root. Time per update is bounded by tree depth.
Algorithm 2 (Optimized): Use Heavy-Light Decomposition. Compose unary Boolean functions $g_u(x)$ along heavy chains. Use a Segment Tree to maintain these compositions. The update time is $O(\log^2 m)$, or $O(\log m)$ using Global Balanced Binary Trees.

Total Complexity: Initialization takes $O(m)$. The leaves of the expression tree are modified constant times per circle scan. Problem 3.1 is solved in $O(m \log m)$. The total time for Problem 3 is $O(nm \log m)$.

5. 💡 Applications and Discussion

Spherical Polygons: Can be represented by the intersection of hemispheres (triangles) and their unions.
Plane Composite Regions: By introducing Half-Planes and using dot products (similar to polar angles), the algorithm extends to planes with complexity $O(nm \log m)$. This solves NOI problems like "Rainfall" and "Moonlit Lemon Tree".
Easy Boundary Cases: For problems with special properties (e.g., "Moonlit Lemon Tree"), specialized logic can calculate boundaries in $O(n)$.
DAG Boolean Functions: If variables are reused (DAG instead of Tree), recomputing the DAG takes $O(m)$, leading to total complexity $O(n^2 m)$. For volume calculations (constructive solid geometry), this can be $O(n^3 I)$.
Surface Area of Sphere Union: The surface of a union of spheres is composed of spherical caps. This converts to the Spherical Composite Region Area problem, yielding an $O(n^3 \log n)$ algorithm.

6. 🏁 Conclusion

This paper introduced the Spherical Composite Region Area problem and a Green's formula-based solution. The algorithm is efficient ($O(nm \log m)$) and adaptable to plane geometry. Future work could involve Ellipsoidal composite regions or improving efficiency for DAG-based Boolean functions.

Acknowledgements Thanks to Zhou Hangrui and Xu Zhean for their assistance and discussions.

References

[1] Yang Zhe. Research on SPOJ375 QTREE Solution, 2007.

[2] Hu Weidong. Moonlit Lemon Tree, NOI 2005.

[3] Xu Mingkuan. 2D Computational Geometry Algorithms, NOI Winter Camp 2018.

[4] Spherical Geometry, People's Education Press, 2008.

[5] Chang Gengzhe, Shi Jihuai. Mathematical Analysis Tutorial, 2012.

[6] Area. ACM/ICPC Asia Regional Nanjing Online, 2013.

[7] Todhunter I. Spherical Trigonometry, 1863.