Add 3d contents

content/geometry/3dFaces.h

@@ -0,0 +1,70 @@
+/**
+ * Author: Victor Lecomte
+ * Date: 2009-04-14
+ * License: CC0
+ * Source: https://victorlecomte.com/cp-geo.pdf
+ * Description: Get polyhedron volume
+ * Usage: reorient the faces first
+ * Status: untested
+ */
+#pragma once
+#include "Point3D.h"
+struct edge {
+  int v; bool same; // = is the common edge in the same order?
+};
+// Given a series of faces (lists of points), reverse some of them
+// so that their orientations are consistent
+void reorient(vector<vector<p3>> &fs) {
+  int n = fs.size();
+// Find the common edges and create the resulting graph
+  vector<vector<edge>> g(n);
+  map<pair<p3, p3>, int> es;
+  for (int u = 0; u < n; u++) {
+    for (int i = 0, m = fs[u].size(); i < m; i++) {
+      p3 a = fs[u][i], b = fs[u][(i + 1) % m];
+// Let’s look at edge [AB]
+      if (es.count({a, b})) { // seen in same order
+        int v = es[ {a, b}];
+        g[u].push_back({v, true});
+        g[v].push_back({u, true});
+      } else if (es.count({b, a})) { // seen in different order
+        int v = es[ {b, a}];
+        g[u].push_back({v, false});
+        g[v].push_back({u, false});
+      } else { // not seen yet
+        es[ {a, b}] = u;
+      }
+    }
+  }
+  vector<bool> vis(n, false), flip(n);
+  flip[0] = false;
+  queue<int> q;
+  q.push(0);
+  while (!q.empty()) {
+    int u = q.front();
+    q.pop();
+    for (edge e : g[u]) {
+      if (!vis[e.v]) {
+        vis[e.v] = true;
+        flip[e.v] = (flip[u] ^ e.same);
+        q.push(e.v);
+      }
+    }
+  }
+  for (int u = 0; u < n; u++)
+    if (flip[u])
+      reverse(fs[u].begin(), fs[u].end());
+}
+
+p3 vectorArea2(vector<p3> p) {
+  p3 S = p3(0, 0);
+  for (int i = 0, n = p.size(); i < n; i++)
+    S = S + p[i] * p[(i + 1) % n];
+  return S;
+}
+double area(vector<p3> p) { return vectorArea2(p).dist() / 2.0; }
+double volume(vector<vector<p3>> fs) {
+  double vol6 = 0.0;
+  for (vector<p3> f : fs) vol6 += (vectorArea2(f) | f[0]);
+  return abs(vol6) / 6.0;
+}

content/geometry/Point3D.h

@@ -1,55 +1,114 @@
 /**
- * Author: Ulf Lundstrom with inspiration from tinyKACTL
+ * Author: Ulf Lundstrom, tinyKACTL, hocky, Victor Lecomte
  * Date: 2009-04-14
  * License: CC0
  * Source:
  * Description: Class to handle points in 3D space.
  * 	T can be e.g. double or long long.
  * Usage:
- * Status: tested, except for phi and theta
+ * Status: untested
  */
 #pragma once
 
 template<class T> struct Point3D {
-	typedef Point3D P;
-	typedef const P& R;
-	T x, y, z;
-	explicit Point3D(T x=0, T y=0, T z=0) : x(x), y(y), z(z) {}
-	bool operator<(R p) const {
-		return tie(x, y, z) < tie(p.x, p.y, p.z); }
-	bool operator==(R p) const {
-		return tie(x, y, z) == tie(p.x, p.y, p.z); }
-	P operator+(R p) const { return P(x+p.x, y+p.y, z+p.z); }
-	P operator-(R p) const { return P(x-p.x, y-p.y, z-p.z); }
-	P operator*(T d) const { return P(x*d, y*d, z*d); }
-	P operator/(T d) const { return P(x/d, y/d, z/d); }
-	T dot(R p) const { return x*p.x + y*p.y + z*p.z; }
-	P cross(R p) const {
-		return P(y*p.z - z*p.y, z*p.x - x*p.z, x*p.y - y*p.x);
-	}
-	T dist2() const { return x*x + y*y + z*z; }
-	double dist() const { return sqrt((double)dist2()); }
-	//Azimuthal angle (longitude) to x-axis in interval [-pi, pi]
-	double phi() const { return atan2(y, x); } 
-	//Zenith angle (latitude) to the z-axis in interval [0, pi]
-	double theta() const { return atan2(sqrt(x*x+y*y),z); }
-	P unit() const { return *this/(T)dist(); } //makes dist()=1
-	//returns unit vector normal to *this and p
-	P normal(P p) const { return cross(p).unit(); }
-	//returns point rotated 'angle' radians ccw around axis
-	P rotate(double angle, P axis) const {
-		double s = sin(angle), c = cos(angle); P u = axis.unit();
-		return u*dot(u)*(1-c) + (*this)*c - cross(u)*s;
-	}
+  typedef Point3D P;
+  typedef const P& R;
+  T x, y, z;
+  bool operator<(R p) const {
+    return tie(x, y, z) < tie(p.x, p.y, p.z);
+  }
+  bool operator==(R p) const {
+    return tie(x, y, z) == tie(p.x, p.y, p.z);
+  }
+  P operator+(R p) const { return P(x + p.x, y + p.y, z + p.z); }
+  P operator-(R p) const { return P(x - p.x, y - p.y, z - p.z); }
+  P operator*(T d) const { return P(x * d, y * d, z * d); }
+  P operator/(T d) const { return P(x / d, y / d, z / d); }
+  T dot(R p) const { return x * p.x + y * p.y + z * p.z; }
+  T operator|(R p) const { return (*this).dot(p); }
+  // v * w = (v.dist() * w.dist() * sin(theta)) n;
+  // n is the vector perp to the plane made by v and w
+  P cross(R p) const {return P(y * p.z - z * p.y, z * p.x - x * p.z, x * p.y - y * p.x);}
+  P operator*(R p) const { return (*this).cross(p); }
+  // given plane, check whether this point is below or above it
+  T orient(R p, R q, R r) {return (q - p) * (r - p) | ((*this) - p);}
+  // (*this) is a normal vector of plane p,q,r
+  // determine the 2d orientation of p.ccw(q, r). normal got by cross
+  T orientByNormal(R p, R q, R r) {return (q - p) * (r - p) | (*this);}
+  T dist2() const { return x * x + y * y + z * z; }
+  double dist() const { return sqrt((double)dist2()); }
+  //Azimuthal angle (longitude) to x-axis in interval [-pi, pi]
+  double phi() const { return atan2(y, x); }
+  //Zenith angle (latitude) to the z-axis in interval [0, pi]
+  double theta() const { return atan2(sqrt(x * x + y * y), z); }
+  P unit() const { return *this / (T)dist(); } //makes dist()=1
+  //returns unit vector normal to *this and p
+  P normal(P p) const { return cross(p).unit(); }
+  //returns point rotated 'angle' radians ccw around axis
+  P rotate(double angle, P axis) const {
+    double s = sin(angle), c = cos(angle); P u = axis.unit();
+    return u * dot(u) * (1 - c) + (*this) * c - cross(u) * s;
+  }
+  double angle(const P &other) const {
+    return acos(min(fabs((*this) | other) / dist() / other.dist(), 1.0));
+  }
 };
 
-pair<bool, PD> intersectPoint(PD rayPoint, PD rayVector, Plane plane) {
-	// plane.fi is point passing plane
-	// plane.se is normal vector from O
-  PD diff = rayPoint - plane.fi;
-  LD prod1 = diff.dot(plane.se);
-  LD prod2 = rayVector.dot(plane.se);
-  if (prod2 == 0) return {0, PD()};
-  LD prod3 = prod1 / prod2;
-  return {1, rayPoint - rayVector * prod3};
-}
+
+typedef double T;
+typedef Point3D<T> p3;
+
+
+template <class T> int sgn(T x) { return (x > 0) - (x < 0); }
+struct Plane {
+  // normal vector is better if its a unit vector.
+  p3 n; T d;
+  // From normal n and offset d
+  Plane(p3 _n, T _d) : n(_n.unit()), d(_d) {}
+  // From normal n and point P
+  Plane(p3 _n, p3 p) : n(_n.unit()), d(_n | p) {}
+  // From three non-collinear points P,Q,R
+  Plane(p3 p, p3 q, p3 r) : Plane((q - p) * (r - p), p) {}
+  // weighted distance of point p, if n is unit vector, its signed distance
+  T side(p3 p) {return (n | p) - d;}
+  double dist(p3 p) {return fabs(side(p));}
+  // translate a plane by vector t
+  Plane translate(p3 t) {return {n, d + (n | t)};}
+  Plane shiftUp(double shift) {return {n, d + shift};}
+  p3 proj(p3 p) {return p - n * side(p);}
+  p3 refl(p3 p) {return p - n * 2 * side(p);}
+  double angle(const Plane &p2) const { return n.angle(p2.n); }
+};
+
+struct Line3d {
+  // d is a unit shift vector
+  p3 d, o;
+  Line3d(p3 p, p3 q) : d((q - p).unit()), o(p) {}
+  // From two planes p1, p2 (requires T = double)
+  // if d == 0, plane is parallel
+  Line3d(Plane p1, Plane p2) {
+    d = (p1.n * p2.n).unit();
+    o = (p2.n * p1.d - p1.n * p2.d) * d;
+  }
+  double dist2(p3 p) const {return (d * (p - o)).dist2();}
+  double dist(p3 p) const {return sqrt(dist2(p));}
+  bool cmpProj(p3 p, p3 q) {return (d | p) < (d | q);}
+  p3 proj(p3 p) {return o + d * (d | (p - o));}
+  p3 refl(p3 p) {return proj(p) * 2 - p;}
+  // when d dot p.n is 0, line is parallel with plane
+  p3 inter(Plane p) {return o - d * p.side(o) / (d | p.n);}
+  double dist(const Line3d &other) const {
+    p3 n = d * other.d;
+    if (n == p3(0, 0)) // parallel
+      return dist(other.o);
+    return fabs((other.o - o) | n);
+  }
+  p3 closestFromLine(const Line3d &other) const {
+    p3 n2 = other.d * (d * other.d);
+    return o + d * ((other.o - o) | n2) / (d | n2);
+  }
+  double angle(const Line3d &other) {
+    return d.angle(other.d);
+  }
+  double angle(const Plane &p) { return PI / 2 - p.n.angle(d); }
+};

content/geometry/Sphere.h

@@ -0,0 +1,87 @@
+/**
+ * Author: Victor Lecomte, chilli
+ * Date: 2019-04-26
+ * License: CC0
+ * Source: https://vlecomte.github.io/cp-geo.pdf
+ * Description: Spheres utility function
+ * Status: untested
+ */
+
+#pragma once
+#include "Point3D.h"
+constexpr p3 zero(0, 0, 0);
+
+p3 sph(double r, double lat, double lon) {
+  lat *= M_PI / 180, lon *= M_PI / 180;
+  return {r * cos(lat)*cos(lon), r * cos(lat)*sin(lon), r * sin(lat)};
+}
+int sphereLine(p3 o, double r, Line3d l, pair<p3, p3> &out) {
+  double h2 = r * r - l.dist2(o);
+  if (h2 < 0) return 0;
+  p3 p = l.proj(o);
+  p3 h = l.d * sqrt(h2) / (l.d).dist();
+  out = {p - h, p + h};
+  return 1 + (h2 > 0);
+}
+double greatCircleDist(p3 o, double r, p3 a, p3 b) {
+  return r * (a - o).angle(b - o);
+}
+bool validSegment(p3 a, p3 b) {
+  return !(a * b == zero) || (a | b) > 0;
+}
+
+bool properInter(p3 a, p3 b, p3 c, p3 d, p3 &out) {
+  p3 ab = a * b, cd = c * d; // normals of planes OAB and OCD
+  int oa = sgn(cd | a), ob = sgn(cd | b), oc = sgn(ab | c), od = sgn(ab | d);
+  out = ab * cd * od;
+  return (oa != ob && oc != od && oa != oc);
+}
+
+bool onSphSegment(p3 a, p3 b, p3 p) {
+  p3 n = a * b;
+  if (n == zero)
+    return a * p == zero && (a | p) > 0;
+  return (n | p) == 0 && (n | a * p) >= 0 && (n | b * p) <= 0;
+}
+
+struct directionSet : vector<p3> {
+  using vector::vector; // import constructors
+  void insert(p3 p) {
+    for (p3 q : *this) if (p * q == zero) return;
+    push_back(p);
+  }
+};
+directionSet intersSph(p3 a, p3 b, p3 c, p3 d) {
+  assert(validSegment(a, b) && validSegment(c, d));
+  p3 out;
+  if (properInter(a, b, c, d, out)) return {out};
+  directionSet s;
+  if (onSphSegment(c, d, a)) s.insert(a);
+  if (onSphSegment(c, d, b)) s.insert(b);
+  if (onSphSegment(a, b, c)) s.insert(c);
+  if (onSphSegment(a, b, d)) s.insert(d);
+  return s;
+}
+double angleSph(p3 a, p3 b, p3 c) {
+  return (a * b).angle(a * c);
+}
+double orientedAngleSph(p3 a, p3 b, p3 c) {
+  if ((a * b | c) >= 0)
+    return angleSph(a, b, c);
+  else
+    return 2 * M_PI - angleSph(a, b, c);
+}
+double areaOnSphere(double r, vector<p3> p) {
+  int n = p.size();
+  double sum = -(n - 2) * M_PI;
+  for (int i = 0; i < n; i++)
+    sum += orientedAngleSph(p[(i + 1) % n], p[(i + 2) % n], p[i]);
+  return r * r * sum;
+}
+
+int windingNumber3D(vector<vector<p3>> fs) {
+  double sum = 0;
+  for (vector<p3> f : fs)
+    sum += remainder(areaOnSphere(1, f), 4 * M_PI);
+  return round(sum / (4 * M_PI));
+}

content/geometry/chapter.tex

@@ -49,11 +49,13 @@ \section{Misc. Point Set Problems}
 	\kactlimport{ClosestPair.h}
 	\kactlimport{ManhattanMST.h}
 	\kactlimport{kdTree.h}
-	% \kactlimport{DelaunayTriangulation.h}
 	\kactlimport{FastDelaunay.h}
 
 \section{3D}
-	\kactlimport{PolyhedronVolume.h}
+	% \kactlimport{PolyhedronVolume.h}
 	\kactlimport{Point3D.h}
+	\kactlimport{DelaunayTriangulation.h}
+	\kactlimport{3dFaces.h}
+	\kactlimport{Sphere.h}
 	\kactlimport{3dHull.h}
-	\kactlimport{sphericalDistance.h}
+	% \kactlimport{sphericalDistance.h}