三次元空間上の三角形
(geometry-3d/triangle.hpp)

• View this file on GitHub
• Last update: 2026-07-25 02:01:37+09:00
• Include: #include "geometry-3d/triangle.hpp"

三次元空間上の三角形に関する面積，包含判定，距離計算を扱う．

三角形を同一直線上にない三点 a, b, c を用いた Triangle3D で表す．

• triangle_normal(t)：三角形 t の法線ベクトル $(b-a)\times(c-a)$ を返す．
• area(t)：三角形 t の面積を返す．
• signed_volume(a, b, c, d)：四面体 $abcd$ の符号付き体積を返す．
• volume(a, b, c, d)：四面体 $abcd$ の体積を返す．
• barycentric_coordinates(t, p)：点 p の重心座標を返す．三角形が退化している場合は nullopt．p が三角形の平面上にあるかは判定しない．
• is_intersect_tp(t, p)：三角形 t が点 p を含むか判定する．辺上と頂点も含む．
• closest_point_tp(t, p)：点 p に最も近い三角形 t 上の点を返す．
• distance_tp(t, p)：三角形 t と点 p の距離を返す．

Depends on

• 三次元幾何の基本要素 (geometry-3d/geometry-base.hpp)
• 三次元空間上の直線 (geometry-3d/line.hpp)
• 三次元空間上の平面 (geometry-3d/plane.hpp)
• 三次元空間上の線分 (geometry-3d/segment.hpp)
• 二次元幾何の基本要素 (geometry/geometry-base.hpp)

Verified with

• verify/geometry-3d/UNIT_geometry_3d.test.cpp

Code

#pragma once

#include "geometry-3d/plane.hpp"
#include "geometry-3d/segment.hpp"

struct Triangle3D {
  Point3D a, b, c;

  Triangle3D() = default;
  Triangle3D(const Point3D& _a, const Point3D& _b, const Point3D& _c)
      : a(_a), b(_b), c(_c) {
    assert(abs(cross(b - a, c - a)) > EPS);
  }
};

Point3D triangle_normal(const Triangle3D& t) {
  return cross(t.b - t.a, t.c - t.a);
}

Real area(const Triangle3D& t) { return abs(triangle_normal(t)) / 2; }

Real signed_volume(const Point3D& a, const Point3D& b, const Point3D& c,
                   const Point3D& d) {
  return triple(b - a, c - a, d - a) / 6;
}

Real volume(const Point3D& a, const Point3D& b, const Point3D& c,
            const Point3D& d) {
  return abs(signed_volume(a, b, c, d));
}

optional<array<Real, 3>> barycentric_coordinates(const Triangle3D& t,
                                                  const Point3D& p) {
  Point3D u = t.b - t.a, v = t.c - t.a, w = p - t.a;
  Real uu = norm(u), uv = dot(u, v), vv = norm(v);
  Real wu = dot(w, u), wv = dot(w, v);
  Real d = uu * vv - uv * uv;
  if (abs(d) < EPS * EPS) return nullopt;
  Real y = (vv * wu - uv * wv) / d;
  Real z = (uu * wv - uv * wu) / d;
  return array<Real, 3>{1 - y - z, y, z};
}

bool is_intersect_tp(const Triangle3D& t, const Point3D& p) {
  Plane3D pl{t.a, t.b, t.c};
  if (distance_pp(pl, p) >= EPS) return false;
  auto b = barycentric_coordinates(t, p);
  if (!b) return false;
  return (*b)[0] >= -EPS && (*b)[1] >= -EPS && (*b)[2] >= -EPS;
}

Point3D closest_point_tp(const Triangle3D& t, const Point3D& p) {
  Plane3D pl{t.a, t.b, t.c};
  assert(norm(pl.n) > EPS * EPS);
  Point3D q = projection(pl, p);
  if (is_intersect_tp(t, q)) return q;
  array<Point3D, 3> ps = {closest_point_sp(Segment3D{t.a, t.b}, p),
                         closest_point_sp(Segment3D{t.b, t.c}, p),
                         closest_point_sp(Segment3D{t.c, t.a}, p)};
  return *min_element(begin(ps), end(ps), [&](const Point3D& x, const Point3D& y) {
    return norm(x - p) < norm(y - p);
  });
}

Real distance_tp(const Triangle3D& t, const Point3D& p) {
  return abs(p - closest_point_tp(t, p));
}

/**
 * @brief 三次元空間上の三角形
 * @docs docs/geometry-3d/triangle.md
 */

#line 2 "geometry-3d/triangle.hpp"

#line 2 "geometry-3d/plane.hpp"

#line 2 "geometry-3d/line.hpp"

#line 2 "geometry-3d/geometry-base.hpp"

#line 2 "geometry/geometry-base.hpp"

#include <bits/stdc++.h>

using Real = long double;
constexpr Real EPS = 1e-10;
constexpr Real PI = 3.141592653589793238462643383279L;
bool equals(Real x, Real y) { return fabs(x - y) < EPS; }
int sign(Real a) { return equals(a, 0) ? 0 : (a > 0 ? 1 : -1); }

template <class R>
struct PointBase {
  using P = PointBase;
  R x, y;
  PointBase() : x(0), y(0) {}
  PointBase(R _x, R _y) : x(_x), y(_y) {}
  template <typename T, typename U>
  PointBase(const pair<T, U>& p) : x(p.first), y(p.second) {}

  P operator+(const P& r) const { return P{x + r.x, y + r.y}; }
  P operator-(const P& r) const { return P{x - r.x, y - r.y}; }
  P operator-() const { return P{-x, -y}; }
  P operator*(R r) const { return P{x * r, y * r}; }
  P operator/(R r) const { return P{x / r, y / r}; }

  P& operator+=(const P& r) { return (*this) = (*this) + r; }
  P& operator-=(const P& r) { return (*this) = (*this) - r; }
  P& operator*=(R r) { return (*this) = (*this) * r; }
  P& operator/=(R r) { return (*this) = (*this) / r; }

  bool operator<(const P& r) const { return x != r.x ? x < r.x : y < r.y; }
  bool operator==(const P& r) const { return x == r.x and y == r.y; }
  bool operator!=(const P& r) const { return !((*this) == r); }

  P rotate(R rad) const {
    return {x * cos(rad) - y * sin(rad), x * sin(rad) + y * cos(rad)};
  }
  P rotate90() const { return {-y, x}; }

  R real() const { return x; }
  R imag() const { return y; }

  friend P operator*(R r, const P& p) { return p * r; }
  friend R real(const P& p) { return p.x; }
  friend R imag(const P& p) { return p.y; }
  friend R dot(const P& l, const P& r) { return l.x * r.x + l.y * r.y; }
  friend R cross(const P& l, const P& r) { return l.x * r.y - l.y * r.x; }
  friend R abs(const P& p) { return sqrt(p.x * p.x + p.y * p.y); }
  friend R norm(const P& p) { return p.x * p.x + p.y * p.y; }
  friend R arg(const P& p) { return atan2(p.y, p.x); }

  friend istream& operator>>(istream& is, P& p) {
    R a, b;
    is >> a >> b;
    p = P{a, b};
    return is;
  }
  friend ostream& operator<<(ostream& os, const P& p) {
    return os << p.x << " " << p.y;
  }
};
using Point = PointBase<Real>;
using Points = vector<Point>;

// relative position of c from a->b
int ccw(const Point& a, const Point& b, const Point& c) {
  Point x = b - a, y = c - a;
  if (cross(x, y) > EPS) return +1;        // counter-clockwise
  if (cross(x, y) < -EPS) return -1;       // clockwise
  if (dot(x, y) < -EPS) return +2;         // collinear in the order c-a-b
  if (norm(x) + EPS < norm(y)) return -2;  // collinear in the order a-b-c
  return 0;                                // collinear in the order a-c-b
}

/**
 * @brief 二次元幾何の基本要素
 * @docs docs/geometry/geometry-base.md
 */
#line 4 "geometry-3d/geometry-base.hpp"

template <class R>
struct Point3DBase {
  using P = Point3DBase;
  R x, y, z;

  Point3DBase() : x(0), y(0), z(0) {}
  Point3DBase(R _x, R _y, R _z) : x(_x), y(_y), z(_z) {}

  P operator+(const P& r) const { return {x + r.x, y + r.y, z + r.z}; }
  P operator-(const P& r) const { return {x - r.x, y - r.y, z - r.z}; }
  P operator-() const { return {-x, -y, -z}; }
  P operator*(R r) const { return {x * r, y * r, z * r}; }
  P operator/(R r) const { return {x / r, y / r, z / r}; }

  P& operator+=(const P& r) { return (*this) = (*this) + r; }
  P& operator-=(const P& r) { return (*this) = (*this) - r; }
  P& operator*=(R r) { return (*this) = (*this) * r; }
  P& operator/=(R r) { return (*this) = (*this) / r; }

  bool operator<(const P& r) const {
    if (x != r.x) return x < r.x;
    return y != r.y ? y < r.y : z < r.z;
  }
  bool operator==(const P& r) const {
    return x == r.x && y == r.y && z == r.z;
  }
  bool operator!=(const P& r) const { return !((*this) == r); }

  friend P operator*(R r, const P& p) { return p * r; }
  friend R dot(const P& l, const P& r) {
    return l.x * r.x + l.y * r.y + l.z * r.z;
  }
  friend P cross(const P& l, const P& r) {
    return {l.y * r.z - l.z * r.y, l.z * r.x - l.x * r.z,
            l.x * r.y - l.y * r.x};
  }
  friend R norm(const P& p) { return dot(p, p); }
  friend R abs(const P& p) { return sqrt(norm(p)); }

  friend istream& operator>>(istream& is, P& p) {
    return is >> p.x >> p.y >> p.z;
  }
  friend ostream& operator<<(ostream& os, const P& p) {
    return os << p.x << " " << p.y << " " << p.z;
  }
};

using Point3D = Point3DBase<Real>;
using Points3D = vector<Point3D>;

Real triple(const Point3D& a, const Point3D& b, const Point3D& c) {
  return dot(a, cross(b, c));
}

Point3D normalize(const Point3D& p) {
  Real len = abs(p);
  assert(len > EPS);
  return p / len;
}

Real angle(const Point3D& a, const Point3D& b) {
  Real d = abs(a) * abs(b);
  assert(d > EPS);
  return acos(clamp(dot(a, b) / d, (Real)-1, (Real)1));
}

bool equals(const Point3D& a, const Point3D& b) { return abs(a - b) < EPS; }

/**
 * @brief 三次元幾何の基本要素
 * @docs docs/geometry-3d/geometry-base.md
 */
#line 4 "geometry-3d/line.hpp"

struct Line3D {
  Point3D a, b;

  Line3D() = default;
  Line3D(const Point3D& _a, const Point3D& _b) : a(_a), b(_b) {
    assert(abs(b - a) > EPS);
  }

  friend istream& operator>>(istream& is, Line3D& l) { return is >> l.a >> l.b; }
  friend ostream& operator<<(ostream& os, const Line3D& l) {
    return os << l.a << " to " << l.b;
  }
};

using Lines3D = vector<Line3D>;

bool is_intersect_lp(const Line3D& l, const Point3D& p) {
  Point3D d = l.b - l.a;
  assert(abs(d) > EPS);
  return abs(cross(d, p - l.a)) < EPS * abs(d);
}

bool is_parallel(const Line3D& l, const Line3D& m) {
  Point3D u = l.b - l.a, v = m.b - m.a;
  assert(abs(u) > EPS && abs(v) > EPS);
  return abs(cross(u, v)) < EPS * abs(u) * abs(v);
}

bool is_orthogonal(const Line3D& l, const Line3D& m) {
  Point3D u = l.b - l.a, v = m.b - m.a;
  assert(abs(u) > EPS && abs(v) > EPS);
  return abs(dot(u, v)) < EPS * abs(u) * abs(v);
}

Point3D projection(const Line3D& l, const Point3D& p) {
  Point3D d = l.b - l.a;
  assert(norm(d) > EPS * EPS);
  return l.a + d * (dot(p - l.a, d) / norm(d));
}

Point3D reflection(const Line3D& l, const Point3D& p) {
  return projection(l, p) * 2 - p;
}

Real distance_lp(const Line3D& l, const Point3D& p) {
  return abs(p - projection(l, p));
}

pair<Point3D, Point3D> closest_points_ll(const Line3D& l, const Line3D& m) {
  Point3D u = l.b - l.a, v = m.b - m.a, w = l.a - m.a;
  Real a = norm(u), b = dot(u, v), c = norm(v);
  assert(a > EPS * EPS && c > EPS * EPS);
  Real d = dot(u, w), e = dot(v, w);
  Real det = a * c - b * b;
  if (abs(det) < EPS * EPS * a * c) {
    Point3D q = projection(m, l.a);
    return {l.a, q};
  }
  Real s = (b * e - c * d) / det;
  Real t = (a * e - b * d) / det;
  return {l.a + u * s, m.a + v * t};
}

Real distance_ll(const Line3D& l, const Line3D& m) {
  auto [p, q] = closest_points_ll(l, m);
  return abs(p - q);
}

bool is_intersect_ll(const Line3D& l, const Line3D& m) {
  return distance_ll(l, m) < EPS;
}

optional<Point3D> cross_point_ll(const Line3D& l, const Line3D& m) {
  if (is_parallel(l, m)) return nullopt;
  auto [p, q] = closest_points_ll(l, m);
  if (abs(p - q) >= EPS) return nullopt;
  return (p + q) / 2;
}

/**
 * @brief 三次元空間上の直線
 * @docs docs/geometry-3d/line.md
 */
#line 4 "geometry-3d/plane.hpp"

struct Plane3D {
  Point3D p, n;

  Plane3D() = default;
  Plane3D(const Point3D& _p, const Point3D& _n) : p(_p), n(_n) {
    assert(abs(n) > EPS);
  }
  Plane3D(const Point3D& a, const Point3D& b, const Point3D& c)
      : p(a), n(cross(b - a, c - a)) {
    assert(abs(n) > EPS);
  }
};

using Planes3D = vector<Plane3D>;

Real plane_value(const Plane3D& pl, const Point3D& p) {
  return dot(pl.n, p - pl.p);
}

bool is_intersect_pp(const Plane3D& pl, const Point3D& p) {
  return abs(plane_value(pl, p)) < EPS * abs(pl.n);
}

bool is_parallel(const Plane3D& a, const Plane3D& b) {
  return abs(cross(a.n, b.n)) < EPS * abs(a.n) * abs(b.n);
}

bool is_orthogonal(const Plane3D& a, const Plane3D& b) {
  return abs(dot(a.n, b.n)) < EPS * abs(a.n) * abs(b.n);
}

Point3D projection(const Plane3D& pl, const Point3D& p) {
  assert(norm(pl.n) > EPS * EPS);
  return p - pl.n * (plane_value(pl, p) / norm(pl.n));
}

Point3D reflection(const Plane3D& pl, const Point3D& p) {
  return projection(pl, p) * 2 - p;
}

Real signed_distance_pp(const Plane3D& pl, const Point3D& p) {
  return plane_value(pl, p) / abs(pl.n);
}

Real distance_pp(const Plane3D& pl, const Point3D& p) {
  return abs(signed_distance_pp(pl, p));
}

bool is_intersect_lp(const Line3D& l, const Plane3D& pl) {
  Point3D d = l.b - l.a;
  bool parallel = abs(dot(pl.n, d)) < EPS * abs(pl.n) * abs(d);
  return !parallel || is_intersect_pp(pl, l.a);
}

optional<Point3D> cross_point_lp(const Line3D& l, const Plane3D& pl) {
  Point3D v = l.b - l.a;
  Real d = dot(pl.n, v);
  if (abs(d) < EPS * abs(pl.n) * abs(v)) return nullopt;
  Real t = dot(pl.n, pl.p - l.a) / d;
  return l.a + v * t;
}

Real distance_lp(const Line3D& l, const Plane3D& pl) {
  return is_intersect_lp(l, pl) ? 0 : distance_pp(pl, l.a);
}

optional<Line3D> cross_line_pp(const Plane3D& a, const Plane3D& b) {
  Point3D d = cross(a.n, b.n);
  Real d2 = norm(d);
  if (d2 < EPS * EPS * norm(a.n) * norm(b.n)) return nullopt;
  Real da = dot(a.n, a.p), db = dot(b.n, b.p);
  Point3D p = cross(da * b.n - db * a.n, d) / d2;
  return Line3D{p, p + d};
}

Real distance_pp(const Plane3D& a, const Plane3D& b) {
  return is_parallel(a, b) ? distance_pp(a, b.p) : 0;
}

/**
 * @brief 三次元空間上の平面
 * @docs docs/geometry-3d/plane.md
 */
#line 2 "geometry-3d/segment.hpp"

#line 4 "geometry-3d/segment.hpp"

struct Segment3D : Line3D {
  Segment3D() = default;
  Segment3D(const Point3D& _a, const Point3D& _b) {
    a = _a;
    b = _b;
  }
};

using Segments3D = vector<Segment3D>;

Point3D closest_point_sp(const Segment3D& s, const Point3D& p) {
  Point3D d = s.b - s.a;
  if (norm(d) <= EPS * EPS) return s.a;
  Real t = clamp(dot(p - s.a, d) / norm(d), (Real)0, (Real)1);
  return s.a + d * t;
}

pair<Point3D, Point3D> closest_points_ss(const Segment3D& s, const Segment3D& t) {
  Point3D d1 = s.b - s.a, d2 = t.b - t.a, r = s.a - t.a;
  Real a = norm(d1), e = norm(d2), x = 0, y = 0;
  if (a <= EPS * EPS && e <= EPS * EPS) return {s.a, t.a};
  if (a <= EPS * EPS) {
    y = clamp(dot(d2, r) / e, (Real)0, (Real)1);
  } else {
    Real c = dot(d1, r);
    if (e <= EPS * EPS) {
      x = clamp(-c / a, (Real)0, (Real)1);
    } else {
      Real b = dot(d1, d2), f = dot(d2, r);
      Real det = a * e - b * b;
      if (abs(det) > EPS * EPS * a * e) {
        x = clamp((b * f - c * e) / det, (Real)0, (Real)1);
      }
      y = (b * x + f) / e;
      if (y < 0) {
        y = 0;
        x = clamp(-c / a, (Real)0, (Real)1);
      } else if (y > 1) {
        y = 1;
        x = clamp((b - c) / a, (Real)0, (Real)1);
      }
    }
  }
  return {s.a + d1 * x, t.a + d2 * y};
}

Real distance_sp(const Segment3D& s, const Point3D& p) {
  return abs(p - closest_point_sp(s, p));
}

Real distance_ss(const Segment3D& s, const Segment3D& t) {
  auto [p, q] = closest_points_ss(s, t);
  return abs(p - q);
}

bool is_intersect_sp(const Segment3D& s, const Point3D& p) {
  return distance_sp(s, p) < EPS;
}

bool is_intersect_ss(const Segment3D& s, const Segment3D& t) {
  return distance_ss(s, t) < EPS;
}

/**
 * @brief 三次元空間上の線分
 * @docs docs/geometry-3d/segment.md
 */
#line 5 "geometry-3d/triangle.hpp"

struct Triangle3D {
  Point3D a, b, c;

  Triangle3D() = default;
  Triangle3D(const Point3D& _a, const Point3D& _b, const Point3D& _c)
      : a(_a), b(_b), c(_c) {
    assert(abs(cross(b - a, c - a)) > EPS);
  }
};

Point3D triangle_normal(const Triangle3D& t) {
  return cross(t.b - t.a, t.c - t.a);
}

Real area(const Triangle3D& t) { return abs(triangle_normal(t)) / 2; }

Real signed_volume(const Point3D& a, const Point3D& b, const Point3D& c,
                   const Point3D& d) {
  return triple(b - a, c - a, d - a) / 6;
}

Real volume(const Point3D& a, const Point3D& b, const Point3D& c,
            const Point3D& d) {
  return abs(signed_volume(a, b, c, d));
}

optional<array<Real, 3>> barycentric_coordinates(const Triangle3D& t,
                                                  const Point3D& p) {
  Point3D u = t.b - t.a, v = t.c - t.a, w = p - t.a;
  Real uu = norm(u), uv = dot(u, v), vv = norm(v);
  Real wu = dot(w, u), wv = dot(w, v);
  Real d = uu * vv - uv * uv;
  if (abs(d) < EPS * EPS) return nullopt;
  Real y = (vv * wu - uv * wv) / d;
  Real z = (uu * wv - uv * wu) / d;
  return array<Real, 3>{1 - y - z, y, z};
}

bool is_intersect_tp(const Triangle3D& t, const Point3D& p) {
  Plane3D pl{t.a, t.b, t.c};
  if (distance_pp(pl, p) >= EPS) return false;
  auto b = barycentric_coordinates(t, p);
  if (!b) return false;
  return (*b)[0] >= -EPS && (*b)[1] >= -EPS && (*b)[2] >= -EPS;
}

Point3D closest_point_tp(const Triangle3D& t, const Point3D& p) {
  Plane3D pl{t.a, t.b, t.c};
  assert(norm(pl.n) > EPS * EPS);
  Point3D q = projection(pl, p);
  if (is_intersect_tp(t, q)) return q;
  array<Point3D, 3> ps = {closest_point_sp(Segment3D{t.a, t.b}, p),
                         closest_point_sp(Segment3D{t.b, t.c}, p),
                         closest_point_sp(Segment3D{t.c, t.a}, p)};
  return *min_element(begin(ps), end(ps), [&](const Point3D& x, const Point3D& y) {
    return norm(x - p) < norm(y - p);
  });
}

Real distance_tp(const Triangle3D& t, const Point3D& p) {
  return abs(p - closest_point_tp(t, p));
}

/**
 * @brief 三次元空間上の三角形
 * @docs docs/geometry-3d/triangle.md
 */

Back to top page