7.5.2 Polyhedron Volume

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

/*

Given an oriented triangular mesh of a closed polyhedron, compute its signed volume. Each triangle
contributes the signed volume of the tetrahedron formed by the triangle $abc$ and the origin
$((a \times b) \cdot c) / 6$. If faces are oriented outward, the result is positive for the
orientation used by `convex_hull_3d()` in 7.5.3; take `abs()` if only the magnitude is needed.

- `signed_polyhedron_volume(p, faces)` returns the signed volume. Each face triplet (`a`, `b`, `c`)
  represents the triangle with vertices `p[a]`, `p[b]`, and `p[c]`.

Overflow warning: for integral point types, each scalar triple product must fit in the point's
coordinate arithmetic type before it is converted to `double`.

Time Complexity:
- O(f) per call, where $f$ is the number of triangular faces.

Space Complexity:
- O(1) auxiliary.

*/

#include <cmath>
#include <tuple>
#include <vector>

struct Point3D {
  double x, y, z;
  Point3D(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
  double dot(const Point3D &p) const { return x * p.x + y * p.y + z * p.z; }

  Point3D cross(const Point3D &p) const {
    return {y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x};
  }
};

template<typename Pt, typename F>
double signed_polyhedron_volume(const std::vector<Pt> &p, const std::vector<F> &faces) {
  double volume = 0;
  for (const auto &[a, b, c] : faces) {
    volume += p[a].cross(p[b]).dot(p[c]);  // Overflow warning.
  }
  return volume / 6.0;
}

/*** Example Usage ***/

#include <cassert>
using namespace std;

bool EQ(double a, double b) {
  return fabs(a - b) < 1e-9;
}

int main() {
  using Face = tuple<int, int, int>;  // (a, b, c)
  vector<Point3D> p{{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
  vector<Face> faces{
      Face(0, 2, 1),
      Face(0, 1, 3),
      Face(0, 3, 2),
      Face(1, 2, 3),
  };
  assert(EQ(signed_polyhedron_volume(p, faces), 1.0 / 6.0));
  return 0;
}