geometry_3d.hpp¶
提供三维点向量、混合积、点线面距离、平面投影、线面交及四面体有向体积；退化直线和平面返回空结果。
Complexity: Time: O(1) per primitive. Space: O(1).
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(1) per primitive. Space: O(1).
#include <algorithm>
#include <cassert>
#include <cmath>
#include <optional>
namespace noya {
/// @brief Three-dimensional point/vector with lexicographic order.
/// Vector arithmetic and products use T; choose T wide enough for intermediate
/// products. Metric operations convert coordinates before subtracting them.
template <class T> struct point3 {
T x{};
T y{};
T z{};
point3() = default;
point3(T x_, T y_, T z_) : x(x_), y(y_), z(z_) {}
point3 &operator+=(const point3 &rhs) {
x += rhs.x;
y += rhs.y;
z += rhs.z;
return *this;
}
point3 &operator-=(const point3 &rhs) {
x -= rhs.x;
y -= rhs.y;
z -= rhs.z;
return *this;
}
point3 &operator*=(const T &scl) {
x *= scl;
y *= scl;
z *= scl;
return *this;
}
point3 &operator/=(const T &scl) {
x /= scl;
y /= scl;
z /= scl;
return *this;
}
friend point3 operator+(point3 lhs, const point3 &rhs) { return lhs += rhs; }
friend point3 operator-(point3 lhs, const point3 &rhs) { return lhs -= rhs; }
friend point3 operator*(point3 val, const T &scl) { return val *= scl; }
friend point3 operator*(const T &scl, point3 val) { return val *= scl; }
friend point3 operator/(point3 val, const T &scl) { return val /= scl; }
friend bool operator==(const point3 &, const point3 &) = default;
friend bool operator<(const point3 &lhs, const point3 &rhs) {
if (lhs.x != rhs.x)
return lhs.x < rhs.x;
if (lhs.y != rhs.y)
return lhs.y < rhs.y;
return lhs.z < rhs.z;
}
};
/// @brief Plane through p with normal n; n must be nonzero for a valid plane.
template <class T> struct plane3 {
point3<T> p;
point3<T> n;
};
template <class T> T dot(const point3<T> &a, const point3<T> &b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
template <class T> point3<T> cross(const point3<T> &a, const point3<T> &b) {
return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
}
template <class T> T norm2(const point3<T> &val) { return dot(val, val); }
/// @brief Signed volume of the parallelepiped spanned by a, b, and c.
template <class T>
T scalar_triple_product(const point3<T> &a, const point3<T> &b,
const point3<T> &c) {
return dot(a, cross(b, c));
}
namespace geometry_3d_internal {
template <class T> point3<long double> cast_point(const point3<T> &val) {
return {static_cast<long double>(val.x), static_cast<long double>(val.y),
static_cast<long double>(val.z)};
}
inline long double length(const point3<long double> &val) {
return std::hypot(val.x, val.y, val.z);
}
} // namespace geometry_3d_internal
template <class T>
long double distance(const point3<T> &a, const point3<T> &b) {
using geometry_3d_internal::cast_point;
return geometry_3d_internal::length(cast_point(a) - cast_point(b));
}
/// @brief Distance to the infinite line through a and b; nullopt if a == b.
template <class T>
std::optional<long double> point_line_distance(const point3<T> &p,
const point3<T> &a,
const point3<T> &b) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto dir = cast_point(b) - cast_point(a);
long double len = length(dir);
if (len == 0)
return std::nullopt;
return length(cross(cast_point(p) - cast_point(a), dir / len));
}
/// @brief Distance to the closed segment [a, b]; equal endpoints form a point.
template <class T>
long double point_segment_distance(const point3<T> &p, const point3<T> &a,
const point3<T> &b) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto rel = cast_point(p) - cast_point(a);
auto dir = cast_point(b) - cast_point(a);
long double len = length(dir);
if (len == 0)
return length(rel);
dir /= len;
long double pro = std::clamp(dot(rel, dir), 0.0L, len);
return length(rel - dir * pro);
}
/// @brief Plane through three points with a unit normal following (b-a)x(c-a).
/// Return nullopt for repeated points or when the sine of their angle is at
/// most eps. eps is dimensionless and must be nonnegative.
template <class T>
std::optional<plane3<long double>>
plane_through_points(const point3<T> &a, const point3<T> &b, const point3<T> &c,
long double eps = 1e-12L) {
assert(eps >= 0);
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto lhs = cast_point(b) - cast_point(a);
auto rhs = cast_point(c) - cast_point(a);
long double al = length(lhs), bl = length(rhs);
if (al == 0 || bl == 0)
return std::nullopt;
auto nm = cross(lhs / al, rhs / bl);
long double nl = length(nm);
if (nl <= eps)
return std::nullopt;
return plane3<long double>{cast_point(a), nm / nl};
}
/// @brief Perpendicular distance to a plane; nullopt for a zero normal.
template <class T, class U>
std::optional<long double> point_plane_distance(const point3<T> &p,
const plane3<U> &pl) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto nm = cast_point(pl.n);
long double len = length(nm);
if (len == 0)
return std::nullopt;
return std::abs(dot(cast_point(p) - cast_point(pl.p), nm / len));
}
/// @brief Orthogonal projection onto a plane; nullopt for a zero normal.
template <class T, class U>
std::optional<point3<long double>> plane_projection(const point3<T> &p,
const plane3<U> &pl) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto nm = cast_point(pl.n);
long double len = length(nm);
if (len == 0)
return std::nullopt;
nm /= len;
auto pt = cast_point(p);
return pt - nm * dot(pt - cast_point(pl.p), nm);
}
/// @brief Unique intersection of the infinite line (a, b) with a plane.
/// Return nullopt for a degenerate line/plane, a parallel line, or a line
/// contained in the plane. eps bounds the absolute dot product of the unit
/// line direction and unit plane normal, so scaling the normal has no effect.
template <class T, class U>
std::optional<point3<long double>>
line_plane_intersection(const point3<T> &a, const point3<T> &b,
const plane3<U> &pl, long double eps = 1e-12L) {
assert(eps >= 0);
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto dir = cast_point(b) - cast_point(a);
auto nm = cast_point(pl.n);
long double dl = length(dir), nl = length(nm);
if (dl == 0 || nl == 0)
return std::nullopt;
dir /= dl;
nm /= nl;
long double den = dot(dir, nm);
if (std::abs(den) <= eps)
return std::nullopt;
auto pt = cast_point(a);
return pt + dir * (dot(cast_point(pl.p) - pt, nm) / den);
}
/// @brief Signed tetrahedron volume det(b-a, c-a, d-a)/6, in long double.
template <class T>
long double signed_tetrahedron_volume(const point3<T> &a, const point3<T> &b,
const point3<T> &c, const point3<T> &d) {
using geometry_3d_internal::cast_point;
auto o = cast_point(a);
return scalar_triple_product(cast_point(b) - o, cast_point(c) - o,
cast_point(d) - o) /
6;
}
} // namespace noya
// NOYA_GEOMETRY_3D_HPP
#ifndef NOYA_GEOMETRY_3D_HPP
#define NOYA_GEOMETRY_3D_HPP 1
/// @complexity Time: O(1) per primitive. Space: O(1).
#include <algorithm>
#include <cassert>
#include <cmath>
#include <optional>
namespace noya {
/// @brief Three-dimensional point/vector with lexicographic order.
/// Vector arithmetic and products use T; choose T wide enough for intermediate
/// products. Metric operations convert coordinates before subtracting them.
template <class T> struct point3 {
T x{};
T y{};
T z{};
point3() = default;
point3(T x_, T y_, T z_) : x(x_), y(y_), z(z_) {}
point3 &operator+=(const point3 &rhs) {
x += rhs.x;
y += rhs.y;
z += rhs.z;
return *this;
}
point3 &operator-=(const point3 &rhs) {
x -= rhs.x;
y -= rhs.y;
z -= rhs.z;
return *this;
}
point3 &operator*=(const T &scl) {
x *= scl;
y *= scl;
z *= scl;
return *this;
}
point3 &operator/=(const T &scl) {
x /= scl;
y /= scl;
z /= scl;
return *this;
}
friend point3 operator+(point3 lhs, const point3 &rhs) { return lhs += rhs; }
friend point3 operator-(point3 lhs, const point3 &rhs) { return lhs -= rhs; }
friend point3 operator*(point3 val, const T &scl) { return val *= scl; }
friend point3 operator*(const T &scl, point3 val) { return val *= scl; }
friend point3 operator/(point3 val, const T &scl) { return val /= scl; }
friend bool operator==(const point3 &, const point3 &) = default;
friend bool operator<(const point3 &lhs, const point3 &rhs) {
if (lhs.x != rhs.x)
return lhs.x < rhs.x;
if (lhs.y != rhs.y)
return lhs.y < rhs.y;
return lhs.z < rhs.z;
}
};
/// @brief Plane through p with normal n; n must be nonzero for a valid plane.
template <class T> struct plane3 {
point3<T> p;
point3<T> n;
};
template <class T> T dot(const point3<T> &a, const point3<T> &b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
template <class T> point3<T> cross(const point3<T> &a, const point3<T> &b) {
return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
}
template <class T> T norm2(const point3<T> &val) { return dot(val, val); }
/// @brief Signed volume of the parallelepiped spanned by a, b, and c.
template <class T>
T scalar_triple_product(const point3<T> &a, const point3<T> &b,
const point3<T> &c) {
return dot(a, cross(b, c));
}
namespace geometry_3d_internal {
template <class T> point3<long double> cast_point(const point3<T> &val) {
return {static_cast<long double>(val.x), static_cast<long double>(val.y),
static_cast<long double>(val.z)};
}
inline long double length(const point3<long double> &val) {
return std::hypot(val.x, val.y, val.z);
}
} // namespace geometry_3d_internal
template <class T>
long double distance(const point3<T> &a, const point3<T> &b) {
using geometry_3d_internal::cast_point;
return geometry_3d_internal::length(cast_point(a) - cast_point(b));
}
/// @brief Distance to the infinite line through a and b; nullopt if a == b.
template <class T>
std::optional<long double> point_line_distance(const point3<T> &p,
const point3<T> &a,
const point3<T> &b) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto dir = cast_point(b) - cast_point(a);
long double len = length(dir);
if (len == 0)
return std::nullopt;
return length(cross(cast_point(p) - cast_point(a), dir / len));
}
/// @brief Distance to the closed segment [a, b]; equal endpoints form a point.
template <class T>
long double point_segment_distance(const point3<T> &p, const point3<T> &a,
const point3<T> &b) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto rel = cast_point(p) - cast_point(a);
auto dir = cast_point(b) - cast_point(a);
long double len = length(dir);
if (len == 0)
return length(rel);
dir /= len;
long double pro = std::clamp(dot(rel, dir), 0.0L, len);
return length(rel - dir * pro);
}
/// @brief Plane through three points with a unit normal following (b-a)x(c-a).
/// Return nullopt for repeated points or when the sine of their angle is at
/// most eps. eps is dimensionless and must be nonnegative.
template <class T>
std::optional<plane3<long double>>
plane_through_points(const point3<T> &a, const point3<T> &b, const point3<T> &c,
long double eps = 1e-12L) {
assert(eps >= 0);
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto lhs = cast_point(b) - cast_point(a);
auto rhs = cast_point(c) - cast_point(a);
long double al = length(lhs), bl = length(rhs);
if (al == 0 || bl == 0)
return std::nullopt;
auto nm = cross(lhs / al, rhs / bl);
long double nl = length(nm);
if (nl <= eps)
return std::nullopt;
return plane3<long double>{cast_point(a), nm / nl};
}
/// @brief Perpendicular distance to a plane; nullopt for a zero normal.
template <class T, class U>
std::optional<long double> point_plane_distance(const point3<T> &p,
const plane3<U> &pl) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto nm = cast_point(pl.n);
long double len = length(nm);
if (len == 0)
return std::nullopt;
return std::abs(dot(cast_point(p) - cast_point(pl.p), nm / len));
}
/// @brief Orthogonal projection onto a plane; nullopt for a zero normal.
template <class T, class U>
std::optional<point3<long double>> plane_projection(const point3<T> &p,
const plane3<U> &pl) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto nm = cast_point(pl.n);
long double len = length(nm);
if (len == 0)
return std::nullopt;
nm /= len;
auto pt = cast_point(p);
return pt - nm * dot(pt - cast_point(pl.p), nm);
}
/// @brief Unique intersection of the infinite line (a, b) with a plane.
/// Return nullopt for a degenerate line/plane, a parallel line, or a line
/// contained in the plane. eps bounds the absolute dot product of the unit
/// line direction and unit plane normal, so scaling the normal has no effect.
template <class T, class U>
std::optional<point3<long double>>
line_plane_intersection(const point3<T> &a, const point3<T> &b,
const plane3<U> &pl, long double eps = 1e-12L) {
assert(eps >= 0);
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto dir = cast_point(b) - cast_point(a);
auto nm = cast_point(pl.n);
long double dl = length(dir), nl = length(nm);
if (dl == 0 || nl == 0)
return std::nullopt;
dir /= dl;
nm /= nl;
long double den = dot(dir, nm);
if (std::abs(den) <= eps)
return std::nullopt;
auto pt = cast_point(a);
return pt + dir * (dot(cast_point(pl.p) - pt, nm) / den);
}
/// @brief Signed tetrahedron volume det(b-a, c-a, d-a)/6, in long double.
template <class T>
long double signed_tetrahedron_volume(const point3<T> &a, const point3<T> &b,
const point3<T> &c, const point3<T> &d) {
using geometry_3d_internal::cast_point;
auto o = cast_point(a);
return scalar_triple_product(cast_point(b) - o, cast_point(c) - o,
cast_point(d) - o) /
6;
}
} // namespace noya
#endif // NOYA_GEOMETRY_3D_HPP
/// @complexity Time: O(1) per primitive. Space: O(1).
#include <algorithm>
#include <cassert>
#include <cmath>
#include <optional>
namespace noya {
/// @brief Three-dimensional point/vector with lexicographic order.
/// Vector arithmetic and products use T; choose T wide enough for intermediate
/// products. Metric operations convert coordinates before subtracting them.
template <class T> struct point3 {
T x{};
T y{};
T z{};
point3() = default;
point3(T x_, T y_, T z_) : x(x_), y(y_), z(z_) {}
point3 &operator+=(const point3 &rhs) {
x += rhs.x;
y += rhs.y;
z += rhs.z;
return *this;
}
point3 &operator-=(const point3 &rhs) {
x -= rhs.x;
y -= rhs.y;
z -= rhs.z;
return *this;
}
point3 &operator*=(const T &scl) {
x *= scl;
y *= scl;
z *= scl;
return *this;
}
point3 &operator/=(const T &scl) {
x /= scl;
y /= scl;
z /= scl;
return *this;
}
friend point3 operator+(point3 lhs, const point3 &rhs) { return lhs += rhs; }
friend point3 operator-(point3 lhs, const point3 &rhs) { return lhs -= rhs; }
friend point3 operator*(point3 val, const T &scl) { return val *= scl; }
friend point3 operator*(const T &scl, point3 val) { return val *= scl; }
friend point3 operator/(point3 val, const T &scl) { return val /= scl; }
friend bool operator==(const point3 &, const point3 &) = default;
friend bool operator<(const point3 &lhs, const point3 &rhs) {
if (lhs.x != rhs.x)
return lhs.x < rhs.x;
if (lhs.y != rhs.y)
return lhs.y < rhs.y;
return lhs.z < rhs.z;
}
};
/// @brief Plane through p with normal n; n must be nonzero for a valid plane.
template <class T> struct plane3 {
point3<T> p;
point3<T> n;
};
template <class T> T dot(const point3<T> &a, const point3<T> &b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}
template <class T> point3<T> cross(const point3<T> &a, const point3<T> &b) {
return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
}
template <class T> T norm2(const point3<T> &val) { return dot(val, val); }
/// @brief Signed volume of the parallelepiped spanned by a, b, and c.
template <class T>
T scalar_triple_product(const point3<T> &a, const point3<T> &b,
const point3<T> &c) {
return dot(a, cross(b, c));
}
namespace geometry_3d_internal {
template <class T> point3<long double> cast_point(const point3<T> &val) {
return {static_cast<long double>(val.x), static_cast<long double>(val.y),
static_cast<long double>(val.z)};
}
inline long double length(const point3<long double> &val) {
return std::hypot(val.x, val.y, val.z);
}
} // namespace geometry_3d_internal
template <class T>
long double distance(const point3<T> &a, const point3<T> &b) {
using geometry_3d_internal::cast_point;
return geometry_3d_internal::length(cast_point(a) - cast_point(b));
}
/// @brief Distance to the infinite line through a and b; nullopt if a == b.
template <class T>
std::optional<long double> point_line_distance(const point3<T> &p,
const point3<T> &a,
const point3<T> &b) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto dir = cast_point(b) - cast_point(a);
long double len = length(dir);
if (len == 0)
return std::nullopt;
return length(cross(cast_point(p) - cast_point(a), dir / len));
}
/// @brief Distance to the closed segment [a, b]; equal endpoints form a point.
template <class T>
long double point_segment_distance(const point3<T> &p, const point3<T> &a,
const point3<T> &b) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto rel = cast_point(p) - cast_point(a);
auto dir = cast_point(b) - cast_point(a);
long double len = length(dir);
if (len == 0)
return length(rel);
dir /= len;
long double pro = std::clamp(dot(rel, dir), 0.0L, len);
return length(rel - dir * pro);
}
/// @brief Plane through three points with a unit normal following (b-a)x(c-a).
/// Return nullopt for repeated points or when the sine of their angle is at
/// most eps. eps is dimensionless and must be nonnegative.
template <class T>
std::optional<plane3<long double>>
plane_through_points(const point3<T> &a, const point3<T> &b, const point3<T> &c,
long double eps = 1e-12L) {
assert(eps >= 0);
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto lhs = cast_point(b) - cast_point(a);
auto rhs = cast_point(c) - cast_point(a);
long double al = length(lhs), bl = length(rhs);
if (al == 0 || bl == 0)
return std::nullopt;
auto nm = cross(lhs / al, rhs / bl);
long double nl = length(nm);
if (nl <= eps)
return std::nullopt;
return plane3<long double>{cast_point(a), nm / nl};
}
/// @brief Perpendicular distance to a plane; nullopt for a zero normal.
template <class T, class U>
std::optional<long double> point_plane_distance(const point3<T> &p,
const plane3<U> &pl) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto nm = cast_point(pl.n);
long double len = length(nm);
if (len == 0)
return std::nullopt;
return std::abs(dot(cast_point(p) - cast_point(pl.p), nm / len));
}
/// @brief Orthogonal projection onto a plane; nullopt for a zero normal.
template <class T, class U>
std::optional<point3<long double>> plane_projection(const point3<T> &p,
const plane3<U> &pl) {
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto nm = cast_point(pl.n);
long double len = length(nm);
if (len == 0)
return std::nullopt;
nm /= len;
auto pt = cast_point(p);
return pt - nm * dot(pt - cast_point(pl.p), nm);
}
/// @brief Unique intersection of the infinite line (a, b) with a plane.
/// Return nullopt for a degenerate line/plane, a parallel line, or a line
/// contained in the plane. eps bounds the absolute dot product of the unit
/// line direction and unit plane normal, so scaling the normal has no effect.
template <class T, class U>
std::optional<point3<long double>>
line_plane_intersection(const point3<T> &a, const point3<T> &b,
const plane3<U> &pl, long double eps = 1e-12L) {
assert(eps >= 0);
using geometry_3d_internal::cast_point;
using geometry_3d_internal::length;
auto dir = cast_point(b) - cast_point(a);
auto nm = cast_point(pl.n);
long double dl = length(dir), nl = length(nm);
if (dl == 0 || nl == 0)
return std::nullopt;
dir /= dl;
nm /= nl;
long double den = dot(dir, nm);
if (std::abs(den) <= eps)
return std::nullopt;
auto pt = cast_point(a);
return pt + dir * (dot(cast_point(pl.p) - pt, nm) / den);
}
/// @brief Signed tetrahedron volume det(b-a, c-a, d-a)/6, in long double.
template <class T>
long double signed_tetrahedron_volume(const point3<T> &a, const point3<T> &b,
const point3<T> &c, const point3<T> &d) {
using geometry_3d_internal::cast_point;
auto o = cast_point(a);
return scalar_triple_product(cast_point(b) - o, cast_point(c) - o,
cast_point(d) - o) /
6;
}
} // namespace noya
// NOYA_GEOMETRY_3D_HPP