Quickhull
Quickhull is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
- Problem Description
Quickhull is a computational geometry algorithm used to compute the convex hull of a set of points. The convex hull is the smallest convex set that contains all the points in the set. In three dimensions, this is the smallest convex polyhedron that encloses all the points, essentially forming the "outer boundary" of the points.
The Quickhull algorithm, also known as the "divide-and-conquer" convex hull algorithm, efficiently constructs the convex hull by recursively partitioning the point set. It is named for its similarity to the Quicksort algorithm due to its use of pivoting and partitioning.
- Input
- A set of points in 3D space, where each point is represented by its coordinates .
- Output
- The surface area of the convex hull of the given set of points.
- Algorithm Overview (for 3D)
1. Find Initial Tetrahedron: Identify a set of four points that form a tetrahedron and are likely to be part of the convex hull. A common approach is to find the two points with minimum and maximum x-coordinates, and then find a point farthest from the line segment connecting them. Then, find a fourth point farthest from the plane defined by these three points. These four points form an initial tetrahedron that is guaranteed to be inside the convex hull.
2. Initialize Faces and Point Sets: The initial tetrahedron has four triangular faces. For each face, maintain a list of points that are *outside* the plane of that face.
3. Recursively Build the Hull: While there are faces with points outside:
- Select a face that has points outside.
- Find the point farthest from the plane of this face. This point is guaranteed to be on the convex hull.
- Identify all faces that are "visible" from this farthest point. A face is visible if the point is on the "outside" of the face's plane.
- Remove the visible faces and the points associated with them (those previously outside those faces).
- Create new triangular faces by connecting the farthest point to the edges of the horizon. The horizon is the boundary of the visible faces on the existing hull.
- For each new face, determine which of the previously outside points are now outside the plane of this new face.
4. Termination: The process terminates when no faces have points outside their respective planes.
5. Calculate Surface Area: Once the convex hull (represented by a set of triangular faces) is constructed, calculate the area of each triangular face and sum them up to get the total surface area of the convex hull. The area of a triangle in 3D can be calculated using the magnitude of the cross product of two of its edge vectors.
- Applications
- Collision Detection: Used in computer graphics and game development to simplify bounding shapes for efficient collision checks.
- Shape Analysis: Helps in analyzing the shape and extent of 3D datasets, useful in fields like medical imaging and scientific simulations.
- Path Planning: Assists in robotics and 3D mapping for navigating complex environments by defining the outer boundaries of obstacles.
- Data Visualization: Simplifies the representation of complex 3D data by outlining their extremities.
- Task
Given n points in 3D space, calculate the surface area of the convex hull formed by these points.
C#
using System;
using System.Collections.Generic;
namespace QuickHull3D
{
class Program
{
const int MAXN = 2500;
const double EPS = 1e-8;
// Globals
static List<Facet> FAC = new List<Facet>();
static List<List<Vect>> pts = new List<List<Vect>>();
static int TIME = 0;
static Edge[,] e = new Edge[2, MAXN];
static int[] vistime = new int[MAXN];
static List<int> queue = new List<int>();
static List<int> resfnew = new List<int>();
static List<int> resfdel = new List<int>();
static List<Vect> respt = new List<Vect>();
static void Main(string[] args)
{
PreConvexHulls();
// Example: unit tetrahedron
int n = 4;
var point = new Vect[n + 1];
point[1] = new Vect(0, 0, 0, 1);
point[2] = new Vect(1, 0, 0, 2);
point[3] = new Vect(0, 1, 0, 3);
point[4] = new Vect(0, 0, 1, 4);
var hull = QuickHull3D(point, n);
Console.WriteLine("{0:F3}", hull.GetSurfaceArea());
}
static void PreConvexHulls()
{
// Reserve index 0
pts.Add(new List<Vect>());
FAC.Add(new Facet()); // dummy facet[0]
// Initialize edge array
for (int i = 0; i < 2; i++)
for (int j = 0; j < MAXN; j++)
e[i, j] = new Edge();
}
class Vect
{
public double X, Y, Z;
public int Id;
public Vect(double x = 0, double y = 0, double z = 0, int id = 0)
{
X = x; Y = y; Z = z; Id = id;
}
public Vect Sub(Vect b) => new Vect(X - b.X, Y - b.Y, Z - b.Z);
public Vect Cross(Vect b) => new Vect(
Y * b.Z - Z * b.Y,
Z * b.X - X * b.Z,
X * b.Y - Y * b.X);
public double Dot(Vect b) => X * b.X + Y * b.Y + Z * b.Z;
public double Mag() => Math.Sqrt(X * X + Y * Y + Z * Z);
public bool Eq(Vect b) => eq(X, b.X) && eq(Y, b.Y) && eq(Z, b.Z);
}
class Line
{
public Vect U, V;
public Line(Vect u, Vect v) { U = u; V = v; }
}
class Plane
{
public Vect U, V, W;
public Plane() { }
public Plane(Vect u, Vect v, Vect w) { U = u; V = v; W = w; }
public Vect Normal() => V.Sub(U).Cross(W.Sub(U));
public Vect VecAt(int i)
{
if (i == 0) return U;
if (i == 1) return V;
return W;
}
public int VecId(int i) => VecAt(i).Id;
}
class Facet
{
public int[] N = new int[3];
public int Id;
public int VisitTime;
public bool IsDeleted;
public Plane P = new Plane();
public Facet() { }
public Facet(int id, Plane p)
{
Id = id; P = p;
}
}
class Edge
{
public int NetId, FacetId;
public Edge() { }
}
class ConvexHulls3d
{
public int Index;
double surfaceArea = 0.0;
public ConvexHulls3d(int idx) { Index = idx; }
void DfsArea(int f)
{
if (FAC[f].VisitTime == TIME) return;
FAC[f].VisitTime = TIME;
var nrm = FAC[f].P.Normal();
surfaceArea += nrm.Mag() / 2.0;
for (int i = 0; i < 3; i++)
DfsArea(FAC[f].N[i]);
}
public double GetSurfaceArea()
{
if (gtr(surfaceArea, 0.0)) return surfaceArea;
TIME++;
DfsArea(Index);
return surfaceArea;
}
public int GetHorizon(int f, Vect p, List<int> resDel)
{
if (!isAbove(p, FAC[f].P)) return 0;
if (FAC[f].VisitTime == TIME) return -1;
FAC[f].VisitTime = TIME;
FAC[f].IsDeleted = true;
resDel.Add(FAC[f].Id);
int ret = -2;
for (int i = 0; i < 3; i++)
{
int r = GetHorizon(FAC[f].N[i], p, resDel);
if (r == 0)
{
int a = FAC[f].P.VecId(i);
int b = FAC[f].P.VecId((i + 1) % 3);
foreach (var pair in new[]{(pt:a, idx:0),(pt:b, idx:1)})
{
int pt = pair.pt; int idx = pair.idx;
if (vistime[pt] != TIME)
{
vistime[pt] = TIME;
e[0, pt].NetId = idx == 0 ? b : a;
e[0, pt].FacetId = FAC[f].N[i];
}
else
{
e[1, pt].NetId = idx == 0 ? b : a;
e[1, pt].FacetId = FAC[f].N[i];
}
}
ret = a;
}
else if (r != -1 && r != -2)
{
ret = r;
}
}
return ret;
}
}
// Utilities
static bool eq(double a, double b) => Math.Abs(a - b) < EPS;
static bool gtr(double a, double b) => a - b > EPS;
static double Abs(double x) => x < 0 ? -x : x;
static double DistPointPlane(Vect v, Plane p)
{
double num = v.Sub(p.U).Dot(p.Normal());
double den = p.Normal().Mag();
return num / den;
}
static double DistPointLine(Vect v, Line l)
{
double d = v.Sub(l.U).Mag();
if (d == 0) return 0;
return l.V.Sub(l.U).Cross(v.Sub(l.U)).Mag() / l.V.Sub(l.U).Mag();
}
static double DistPointPoint(Vect a, Vect b) => a.Sub(b).Mag();
static bool isAbove(Vect v, Plane p) =>
gtr(v.Sub(p.U).Dot(p.Normal()), 0);
static ConvexHulls3d GetStart(Vect[] point, int totp)
{
Vect[] extremes = new Vect[6];
for (int i = 0; i < 6; i++) extremes[i] = point[1];
for (int i = 1; i <= totp; i++)
{
var v = point[i];
if (gtr(v.X, extremes[0].X)) extremes[0] = v;
if (gtr(extremes[1].X, v.X)) extremes[1] = v;
if (gtr(v.Y, extremes[2].Y)) extremes[2] = v;
if (gtr(extremes[3].Y, v.Y)) extremes[3] = v;
if (gtr(v.Z, extremes[4].Z)) extremes[4] = v;
if (gtr(extremes[5].Z, v.Z)) extremes[5] = v;
}
// Furthest pair
Vect s0 = extremes[0], s1 = extremes[1];
for (int i = 0; i < 6; i++)
for (int j = i + 1; j < 6; j++)
{
var d = DistPointPoint(extremes[i], extremes[j]);
if (gtr(d, DistPointPoint(s0, s1)))
{
s0 = extremes[i];
s1 = extremes[j];
}
}
// Furthest from line
var line = new Line(s0, s1);
Vect s2 = extremes[0];
for (int i = 0; i < 6; i++)
{
if (gtr(DistPointLine(extremes[i], line),
DistPointLine(s2, line)))
s2 = extremes[i];
}
// Furthest from plane
Vect s3 = point[1];
var basePlane = new Plane(s0, s1, s2);
for (int i = 1; i <= totp; i++)
{
var d1 = Abs(DistPointPlane(point[i], basePlane));
var d2 = Abs(DistPointPlane(s3, basePlane));
if (gtr(d1, d2)) s3 = point[i];
}
if (gtr(0, DistPointPlane(s3, basePlane)))
(s1, s2) = (s2, s1);
// Create 4 initial facets
int[] f = new int[4];
for (int i = 0; i < 4; i++)
{
FAC.Add(new Facet { Id = FAC.Count });
f[i] = FAC.Count - 1;
}
FAC[f[0]].P = new Plane(s0, s2, s1);
FAC[f[1]].P = new Plane(s0, s1, s3);
FAC[f[2]].P = new Plane(s1, s2, s3);
FAC[f[3]].P = new Plane(s2, s0, s3);
FAC[f[0]].N = new[] { f[3], f[2], f[1] };
FAC[f[1]].N = new[] { f[0], f[2], f[3] };
FAC[f[2]].N = new[] { f[0], f[3], f[1] };
FAC[f[3]].N = new[] { f[0], f[1], f[2] };
// Prepare pts lists
for (int i = 0; i < 4; i++)
pts.Add(new List<Vect>());
// Assign points
for (int i = 1; i <= totp; i++)
{
var v = point[i];
if (v.Eq(s0) || v.Eq(s1) || v.Eq(s2) || v.Eq(s3))
continue;
for (int j = 0; j < 4; j++)
if (isAbove(v, FAC[f[j]].P))
{
pts[f[j]].Add(v);
break;
}
}
return new ConvexHulls3d(f[0]);
}
static ConvexHulls3d QuickHull3D(Vect[] point, int totp)
{
var hull = GetStart(point, totp);
// Initialize queue
queue.Clear();
queue.Add(hull.Index);
for (int i = 0; i < 3; i++)
queue.Add(FAC[hull.Index].N[i]);
int snew = 0;
while (queue.Count > 0)
{
int nf = queue[0];
queue.RemoveAt(0);
if (FAC[nf].IsDeleted || pts[nf].Count == 0)
{
if (!FAC[nf].IsDeleted)
snew = nf;
continue;
}
// Farthest point
var p = pts[nf][0];
foreach (var v in pts[nf])
{
if (gtr(DistPointPlane(v, FAC[nf].P),
DistPointPlane(p, FAC[nf].P)))
p = v;
}
// Find horizon
TIME++;
resfdel.Clear();
int s = hull.GetHorizon(nf, p, resfdel);
// Build new faces
resfnew.Clear();
TIME++;
int from = -1, lastf = -1, fstf = -1;
while (vistime[s] != TIME)
{
vistime[s] = TIME;
int net, f, fnew;
if (e[0, s].NetId == from)
{
net = e[1, s].NetId;
f = e[1, s].FacetId;
}
else
{
net = e[0, s].NetId;
f = e[0, s].FacetId;
}
// Find indices on facet f
int pt1 = 0, pt2 = 0;
for (int i = 0; i < 3; i++)
{
if (point[s].Eq(FAC[f].P.VecAt(i)))
pt1 = i;
if (point[net].Eq(FAC[f].P.VecAt(i)))
pt2 = i;
}
if ((pt1 + 1) % 3 != pt2)
(pt1, pt2) = (pt2, pt1);
// Create new facet
FAC.Add(new Facet
{
Id = FAC.Count,
P = new Plane(FAC[f].P.VecAt(pt2),
FAC[f].P.VecAt(pt1),
p)
});
fnew = FAC.Count - 1;
pts.Add(new List<Vect>());
resfnew.Add(fnew);
FAC[fnew].N[0] = f;
FAC[f].N[pt1] = fnew;
if (lastf >= 0)
{
// Link with previous new facet
if (FAC[fnew].P.VecAt(1).Eq(FAC[lastf].P.U))
{
FAC[fnew].N[1] = lastf;
FAC[lastf].N[2] = fnew;
}
else
{
FAC[fnew].N[2] = lastf;
FAC[lastf].N[1] = fnew;
}
}
else
{
fstf = fnew;
}
lastf = fnew;
from = s;
s = net;
}
// Close the loop
if (FAC[fstf].P.VecAt(1).Eq(FAC[lastf].P.U))
{
FAC[fstf].N[1] = lastf;
FAC[lastf].N[2] = fstf;
}
else
{
FAC[fstf].N[2] = lastf;
FAC[lastf].N[1] = fstf;
}
// Collect deleted points
respt.Clear();
foreach (var fid in resfdel)
{
respt.AddRange(pts[fid]);
pts[fid].Clear();
}
// Reassign
foreach (var v in respt)
{
if (v == p) continue;
foreach (var fid in resfnew)
{
if (isAbove(v, FAC[fid].P))
{
pts[fid].Add(v);
break;
}
}
}
// Enqueue new faces
foreach (var fid in resfnew)
queue.Add(fid);
}
hull.Index = snew;
return hull;
}
}
}
- Output:
2.366
C++
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <vector>
constexpr uint32_t MAX_SIZE = 2'500;
constexpr double EPSILON = 0.000'000'01;
// Numerical utilities
bool is_greater_than(const double& a, const double& b) {
return ( a - b ) > EPSILON;
}
bool is_equal(const double& a, const double& b) {
return std::abs(a - b) < EPSILON;
}
class Vector {
public:
Vector(const double& x = 0, const double& y = 0, const double& z = 0, const uint32_t& id = 0)
: x(x), y(y), z(z), id(id) {}
Vector subtract(const Vector& other) const {
return Vector(x - other.x, y - other.y, z - other.z);
}
Vector vector_product(const Vector& other) const {
return Vector(y * other.z - z * other.y, z * other.x - x * other.z, x * other.y - y * other.x);
}
double scalar_product(const Vector& other) const {
return x * other.x + y * other.y + z * other.z;
}
double magnitude() const {
return std::sqrt(x * x + y * y + z * z);
}
bool equals(const Vector& b) const {
return is_equal(x, b.x) && is_equal(y, b.y) && is_equal(z, b.z);
}
double x, y, z;
uint32_t id;
};
class Line {
public:
Line(const Vector& u, const Vector& v) : u(u), v(v) {}
Vector u, v;
};
class Plane {
public:
Plane(const Vector& u, const Vector& v, const Vector& w) : u(u), v(v), w(w) {}
Plane() {}
Vector normal() const {
return v.subtract(u).vector_product(w.subtract(u));
}
Vector vector_at_index(const uint32_t& i) const {
if ( i == 0 ) { return u; }
if ( i == 1 ) { return v; }
if ( i == 2 ) { return w; }
return 0;
}
uint32_t vector_id(const uint32_t& i) const {
return vector_at_index(i).id;
}
Vector u, v, w;
};
class Facet {
public:
Facet(const uint32_t& id, const Plane& p) : id(id), plane(p) {}
Facet(const uint32_t& id) : id(id) {}
Facet() {}
std::vector<uint32_t> N;
uint32_t id = 0;
uint32_t visited_time = 0;
bool is_deleted = false;
Plane plane;
};
class Edge {
public:
Edge() {}
uint32_t netID = 0;
uint32_t faceID = 0;
};
// Global variables
std::vector<Facet> facets;
std::vector<std::vector<Vector>> hull_points;
uint32_t time_step = 0;
std::vector<std::vector<Edge>> edges = { 2, std::vector<Edge>(MAX_SIZE) };
std::vector<uint32_t> visit_time(MAX_SIZE, 0);
std::vector<uint32_t> queue;
std::vector<uint32_t> resfnew;
std::vector<uint32_t> resfdel;
std::vector<Vector> respt;
// Geometric utilities
double distance_point_plane(const Vector& vec, const Plane& plane) {
const double num = vec.subtract(plane.u).scalar_product(plane.normal());
const double den = plane.normal().magnitude();
return num / den;
}
double distance_point_line(const Vector& vec, const Line& line) {
const double length = vec.subtract(line.u).magnitude();
return ( length == 0.0 ) ? 0.0 :
line.v.subtract(line.u).vector_product(vec.subtract(line.u)).magnitude() / line.v.subtract(line.u).magnitude();
}
double distance_point_point(const Vector& a, const Vector& b) {
return a.subtract(b).magnitude();
}
bool is_above(const Vector& point, const Plane& plane) {
return is_greater_than(point.subtract(plane.u).scalar_product(plane.normal()), 0.0);
}
class ConvexHulls3d {
public:
ConvexHulls3d(const uint32_t& index) : index(index) {}
double get_surface_area() {
if ( is_greater_than(surface_area, 0.0) ) {
return surface_area;
}
time_step++;
dfs_area(index);
return surface_area;
}
uint32_t get_horizon(const uint32_t& f, const Vector& point, std::vector<uint32_t> resDel) const {
Facet Ff = facets[f];
if ( ! is_above(point, Ff.plane) ) {
return 0;
}
if ( Ff.visited_time == time_step ) {
return -1;
}
Ff.visited_time = time_step;
Ff.is_deleted = true;
resDel.emplace_back(Ff.id);
uint32_t result = -2;
for ( uint32_t i = 0; i < 3; ++i ) {
const uint32_t ni = Ff.N[i];
const int32_t horizon = get_horizon(ni, point, resDel);
if ( horizon == 0 ) {
const uint32_t a = facets[f].plane.vector_id(i);
const uint32_t b = facets[f].plane.vector_id((i + 1) % 3);
for ( uint32_t idx = 0; idx < 2; ++idx ) {
const uint32_t pt = ( idx == 0 ) ? a : b;
const uint32_t facet = ni;
if ( visit_time[pt] != time_step ) {
visit_time[pt] = time_step;
edges[0][pt].netID = ( idx == 0 ) ? b : a;
edges[0][pt].faceID = facet;
} else {
edges[1][pt].netID = ( idx == 0 ) ? b : a;
edges[1][pt].faceID = facet;
}
}
result = a;
} else if ( horizon != -1 && horizon != -2 ) {
result = horizon;
}
}
return result;
}
uint32_t index;
double surface_area = 0.0;
private:
void dfs_area(const uint32_t& f) {
if ( facets[f].visited_time == time_step ) {
return;
}
facets[f].visited_time = time_step;
const Vector normal = facets[f].plane.normal();
surface_area += normal.magnitude() / 2.0;
for ( uint32_t i = 0; i < 3; ++i ) {
dfs_area(facets[f].N[i]);
}
}
};
void prepareConvexHulls() {
// Reserve index 0
hull_points.emplace_back(std::vector<Vector>());
facets.emplace_back(Facet());
// Initialize edge vector
for ( uint32_t i = 0; i < 2; ++i ) {
for ( uint32_t j = 0; j < MAX_SIZE; ++j ) {
edges[i][j] = Edge();
}
}
}
ConvexHulls3d get_initial_hull(const std::vector<Vector>& points, const uint32_t& total_points) {
std::vector<Vector> extremes(6);
for ( uint32_t i = 0; i < 6; ++i ) { extremes[i] = points[1]; }
for ( uint32_t i = 1; i <= total_points; ++i ) {
Vector point = points[i];
if ( is_greater_than(point.x, extremes[0].x) ) { extremes[0] = point; }
if ( is_greater_than(extremes[1].x, point.x) ) { extremes[1] = point; }
if ( is_greater_than(point.y, extremes[2].y) ) { extremes[2] = point; }
if ( is_greater_than(extremes[3].y, point.y) ) { extremes[3] = point; }
if ( is_greater_than(point.z, extremes[4].z) ) { extremes[4] = point; }
if ( is_greater_than(extremes[5].z, point.z) ) { extremes[5] = point; }
}
// Furthest pair
Vector extreme0 = extremes[0];
Vector extreme1 = extremes[1];
for ( uint32_t i = 0; i < 6; ++i ) {
for ( uint32_t j = i + 1; j < 6; ++j ) {
const double distance = distance_point_point(extremes[i], extremes[j]);
if ( is_greater_than(distance, distance_point_point(extreme0, extreme1)) ) {
extreme0 = extremes[i];
extreme1 = extremes[j];
}
}
}
// Furthest from line
const Line line(extreme0, extreme1);
Vector extreme2 = extremes[0];
for ( uint32_t i = 0; i < 6; ++i ) {
if ( is_greater_than(distance_point_line(extremes[i], line), distance_point_line(extreme2, line)) ) {
extreme2 = extremes[i];
}
}
// Furthest from plane
Vector extreme3 = points[1];
const Plane basePlane(extreme0, extreme1, extreme2);
for ( uint32_t i = 1; i <= total_points; ++i ) {
const double distance1 = std::fabs(distance_point_plane(points[i], basePlane));
const double distance2 = std::fabs(distance_point_plane(extreme3, basePlane));
if ( is_greater_than(distance1, distance2) ) {
extreme3 = points[i];
}
}
if ( is_greater_than(0, distance_point_plane(extreme3, basePlane)) ) {
std::swap(extreme1, extreme2);
}
// Create 4 initial facets
std::vector<uint32_t> f(4);
for ( uint32_t i = 0; i < 4; ++i ) {
facets.emplace_back(Facet(facets.size()));
f[i] = facets.size() - 1;
}
facets[f[0]].plane = Plane(extreme0, extreme2, extreme1);
facets[f[1]].plane = Plane(extreme0, extreme1, extreme3);
facets[f[2]].plane = Plane(extreme1, extreme2, extreme3);
facets[f[3]].plane = Plane(extreme2, extreme0, extreme3);
facets[f[0]].N = { f[3], f[2], f[1] };
facets[f[1]].N = { f[0], f[2], f[3] };
facets[f[2]].N = { f[0], f[3], f[1] };
facets[f[3]].N = { f[0], f[1], f[2] };
// Prepare hull_points vector
for ( uint32_t i = 0; i < 4; ++i ) {
hull_points.emplace_back(std::vector<Vector>());
}
// Assign points
for ( uint32_t i = 1; i <= total_points; ++i ) {
const Vector point = points[i];
if ( point.equals(extreme0) || point.equals(extreme1) || point.equals(extreme2) || point.equals(extreme3) ) {
continue;
}
for ( uint32_t j = 0; j < 4; ++j ) {
if ( is_above(point, facets[f[j]].plane) ) {
hull_points[f[j]].emplace_back(point);
break;
}
}
}
return ConvexHulls3d(f[0]);
}
ConvexHulls3d QuickHull3D(const std::vector<Vector>& points, const uint32_t& total_points) {
ConvexHulls3d hull = get_initial_hull(points, total_points);
// Initialize queue
queue.clear();
queue.emplace_back(hull.index);
for ( uint32_t i = 0; i < 3; ++i ) {
queue.emplace_back(facets[hull.index].N[i]);
}
uint32_t new_horizon = 0;
while ( ! queue.empty() ) {
uint32_t nf = queue.front();
queue.erase(queue.begin());
if ( facets[nf].is_deleted || hull_points[nf].empty()) {
if ( ! facets[nf].is_deleted ) {
new_horizon = nf;
}
continue;
}
// Farthest point
Vector point = hull_points[nf][0];
for ( const Vector& vec : hull_points[nf] ) {
if ( is_greater_than(distance_point_plane(vec, facets[nf].plane),
distance_point_plane(point, facets[nf].plane)) ) {
point = vec;
}
}
// Find horizon
time_step++;
resfdel.clear();
uint32_t horizon = hull.get_horizon(nf, point, resfdel);
// Build new faces
resfnew.clear();
time_step++;
uint32_t from = -1;
uint32_t last_f = -1;
uint32_t first_f = -1;
while ( visit_time[horizon] != time_step ) {
visit_time[horizon] = time_step;
uint32_t net;
uint32_t f;
uint32_t new_f;
if ( edges[0][horizon].netID == from ) {
net = edges[1][horizon].netID;
f = edges[1][horizon].faceID;
} else {
net = edges[0][horizon].netID;
f = edges[0][horizon].faceID;
}
// Find indices on facet f
uint32_t pt1 = 0;
uint32_t pt2 = 0;
for ( uint32_t i = 0; i < 3; ++i ) {
if ( points[horizon].equals(facets[f].plane.vector_at_index(i)) ) {
pt1 = i;
}
if ( points[net].equals(facets[f].plane.vector_at_index(i)) ) {
pt2 = i;
}
}
if ( ( pt1 + 1 ) % 3 != pt2 ) {
std::swap(pt1, pt2);
}
// Create new facet
facets.emplace_back(Facet(
facets.size(),
Plane(facets[f].plane.vector_at_index(pt2),
facets[f].plane.vector_at_index(pt1), point)));
new_f = facets.size() - 1;
hull_points.emplace_back(std::vector<Vector>());
resfnew.emplace_back(new_f);
facets[new_f].N[0] = f;
facets[f].N[pt1] = new_f;
if ( last_f >= 0 ) {
// Link with previous new facet
if ( facets[new_f].plane.vector_at_index(1).equals(facets[last_f].plane.u) ) {
facets[new_f].N[1] = last_f;
facets[last_f].N[2] = new_f;
} else {
facets[new_f].N[2] = last_f;
facets[last_f].N[1] = new_f;
}
} else {
first_f = new_f;
}
last_f = new_f;
from = horizon;
horizon = net;
}
// Close the loop
if ( facets[first_f].plane.vector_at_index(1).equals(facets[last_f].plane.u) ) {
facets[first_f].N[1] = last_f;
facets[last_f].N[2] = first_f;
} else {
facets[first_f].N[2] = last_f;
facets[last_f].N[1] = first_f;
}
// Collect deleted points
respt.clear();
for ( const uint32_t& f_id : resfdel ) {
respt.insert(respt.end(), hull_points[f_id].begin(), hull_points[f_id].end());
hull_points[f_id].clear();
}
// Reassign
for ( const Vector& vec : respt ) {
if ( vec.equals(point) ) {
continue;
}
for ( const uint32_t& f_id : resfnew ) {
if ( is_above(vec, facets[f_id].plane) ) {
hull_points[f_id].emplace_back(vec);
break;
}
}
}
// Enqueue new faces
for ( const uint32_t& f_id : resfnew ) {
queue.emplace_back(f_id);
}
}
hull.index = new_horizon;
return hull;
}
int main() {
prepareConvexHulls();
// Example: a tetrahedron
constexpr uint32_t n = 4;
std::vector<Vector> points(n + 1);
// An empty or placeholder value for points[0] is required for all examples
points[1] = Vector(0, 0, 0, 1); // The last argument of each vector is its index
points[2] = Vector(1, 0, 0, 2);
points[3] = Vector(0, 1, 0, 3);
points[4] = Vector(0, 0, 1, 4);
ConvexHulls3d hull = QuickHull3D(points, n);
std::cout << std::fixed << std::setprecision(3) << hull.get_surface_area() << std::endl;
}
- Output:
2.366
FreeBASIC
Const NULL As Any Ptr = 0
Const MAXN = 2500
Const EPS = 1e-8
' Vect represents a 3D point/vector
Type Vect
x As Single
y As Single
z As Single
id As Integer
End Type
' Linea is the segment u->v
Type Linea
u As Vect Ptr
v As Vect Ptr
End Type
' Plano is defined by three points u, v, w
Type Plano
u As Vect Ptr
v As Vect Ptr
w As Vect Ptr
End Type
' Faceta is a face in the hull
Type Faceta
n(2) As Integer
id As Integer
vistime As Integer
isdel As Integer ' 0 = false, 1 = true
p As Plano
End Type
' Edge on the horizon
Type Edge
netid As Integer
facetid As Integer
End Type
' ConvexHulls3d manages hull construction
Type ConvexHulls3d
index As Integer
surfacearea As Single
End Type
' Global variables
Dim Shared FAC() As Faceta Ptr
Dim Shared pts(MAXN, MAXN) As Vect Ptr
Dim Shared G_TIME As Integer
Dim Shared e(1, MAXN-1) As Edge
Dim Shared vistimeArr(MAXN-1) As Integer
Dim Shared que() As Integer
Dim Shared resfnew() As Integer
Dim Shared resfdel() As Integer
Dim Shared respt() As Vect Ptr
'' Auxiliary vector functions
Function Vect_Sub(a As Vect Ptr, b As Vect Ptr) As Vect
Dim res As Vect
res.x = a->x - b->x
res.y = a->y - b->y
res.z = a->z - b->z
res.id = 0
Return res
End Function
Function Vect_Cross(a As Vect Ptr, b As Vect Ptr) As Vect
Dim res As Vect
res.x = a->y * b->z - a->z * b->y
res.y = a->z * b->x - a->x * b->z
res.z = a->x * b->y - a->y * b->x
res.id = 0
Return res
End Function
Function Vect_Dot(a As Vect Ptr, b As Vect Ptr) As Single
Return a->x * b->x + a->y * b->y + a->z * b->z
End Function
Function Vect_Mag(a As Vect Ptr) As Single
Return Sqr(a->x * a->x + a->y * a->y + a->z * a->z)
End Function
Function Vect_Eq(a As Vect Ptr, b As Vect Ptr) As Integer
Return Abs(a->x - b->x) < EPS And Abs(a->y - b->y) < EPS And Abs(a->z - b->z) < EPS
End Function
'' Operations with planes
Function Plano_Normal(p As Plano Ptr) As Vect
Dim ab As Vect = Vect_Sub(p->v, p->u)
Dim ac As Vect = Vect_Sub(p->w, p->u)
Return Vect_Cross(@ab, @ac)
End Function
Function distPuntoPlano(v As Vect Ptr, p As Plano) As Single
Dim tmp_sub As Vect = Vect_Sub(v, p.u)
Dim n As Vect = Plano_Normal(@p)
Dim num As Single = Vect_Dot(@tmp_sub, @n)
Dim den As Single = Vect_Mag(@n)
Return num / den
End Function
Function distPuntoLinea(v As Vect Ptr, f As Linea) As Single
Dim sub_fu As Vect = Vect_Sub(f.v, f.u)
Dim sub_vu As Vect = Vect_Sub(v, f.u)
Dim cross As Vect = Vect_Cross(@sub_fu, @sub_vu)
Return Vect_Mag(@cross) / Vect_Mag(@sub_fu)
End Function
Function distPuntoPunto(u As Vect Ptr, v As Vect Ptr) As Single
Dim tmp_sub As Vect = Vect_Sub(u, v)
Return Vect_Mag(@tmp_sub)
End Function
Function isAbove(v As Vect Ptr, p As Plano) As Integer
Dim tmp_sub As Vect = Vect_Sub(v, p.u)
Dim n As Vect = Plano_Normal(@p)
Return (Vect_Dot(@tmp_sub, @n) > EPS)
End Function
' Initialize globals
Sub preConvexHulls()
Redim FAC(0)
FAC(0) = Callocate(Sizeof(Faceta)) ' dummy facet
G_TIME = 0
End Sub
' DFS to accumulate area
Sub dfsArea(h As ConvexHulls3d Ptr, nf As Integer)
If FAC(nf)->vistime = G_TIME Then Exit Sub
FAC(nf)->vistime = G_TIME
Dim nrm As Vect = Plano_Normal(@FAC(nf)->p)
h->surfacearea += Vect_Mag(@nrm) / 2.0
For i As Integer = 0 To 2
dfsArea(h, FAC(nf)->n(i))
Next
End Sub
' GetSurfaceArea returns (and caches) the hull surface area
Function GetSurfaceArea(h As ConvexHulls3d Ptr) As Single
If h->surfacearea > EPS Then Return h->surfacearea
G_TIME += 1
dfsArea(h, h->index)
Return h->surfacearea
End Function
' Recursively find the horizon around point p
Function getHorizon(nf As Integer, p As Vect Ptr, resfdel() As Integer) As Integer
FAC(nf)->vistime = G_TIME
FAC(nf)->isdel = 1
Redim Preserve resfdel(Ubound(resfdel) + 1)
resfdel(Ubound(resfdel)) = nf
Dim s As Integer = 0
For i As Integer = 0 To 2
Dim nei As Integer = FAC(nf)->n(i)
If FAC(nei)->vistime = G_TIME Then Continue For
If isAbove(p, FAC(nei)->p) Then
Dim result As Integer = getHorizon(nei, p, resfdel())
s += result
Else
s += 1
Redim Preserve FAC(Ubound(FAC) + 1)
FAC(Ubound(FAC)) = Callocate(Sizeof(Faceta))
FAC(Ubound(FAC))->id = Ubound(FAC)
Dim newf As Integer = Ubound(FAC)
Redim Preserve resfnew(Ubound(resfnew) + 1)
resfnew(Ubound(resfnew)) = newf
End If
Next
Return s
End Function
Sub AsignarPuntoFacetaa(nf As Integer, v As Vect Ptr)
Dim i As Integer = 0
While pts(nf, i) <> 0 And i < MAXN
i += 1
Wend
If i < MAXN Then pts(nf, i) = v
End Sub
Sub ActualizarCola()
If Ubound(que) >= 1 Then
Dim tmp() As Integer
Redim tmp(Ubound(que) - 1)
For i As Integer = 0 To Ubound(tmp)
tmp(i) = que(i + 1)
Next
Redim que(Ubound(tmp))
For i As Integer = 0 To Ubound(tmp)
que(i) = tmp(i)
Next
Else
Redim que(-1)
End If
End Sub
' Build the initial tetrahedron
Function getStart(punto() As Vect Ptr, totp As Integer) As ConvexHulls3d Ptr
Dim As Vect Ptr pt(5), s(3)
Dim As Integer i, j
For i = 0 To 5
pt(i) = punto(1)
Next
' pick extreme coords
For i = 1 To totp
Dim v As Vect Ptr = punto(i)
If v->x > pt(0)->x Then pt(0) = v
If v->x < pt(1)->x Then pt(1) = v
If v->y > pt(2)->y Then pt(2) = v
If v->y < pt(3)->y Then pt(3) = v
If v->z > pt(4)->z Then pt(4) = v
If v->z < pt(5)->z Then pt(5) = v
Next
' Take the two points with the largest distance
s(0) = pt(0): s(1) = pt(1)
For i = 0 To 5
For j = i + 1 To 5
If distPuntoPunto(pt(i), pt(j)) > distPuntoPunto(s(0), s(1)) Then
s(0) = pt(i): s(1) = pt(j)
End If
Next
Next
' Take the point farthest from the line connecting the two points
Dim L As Linea
L.u = s(0): L.v = s(1)
s(2) = pt(0)
For i = 0 To 5
If distPuntoLinea(pt(i), L) > distPuntoLinea(s(2), L) Then s(2) = pt(i)
Next
' Take the point farthest from the face
Dim basePlano As Plano
basePlano.u = s(0): basePlano.v = s(1): basePlano.w = s(2)
s(3) = punto(1)
For i = 1 To totp
If Abs(distPuntoPlano(punto(i), basePlano)) > Abs(distPuntoPlano(s(3), basePlano)) Then
s(3) = punto(i)
End If
Next
' Ensure that the constructed face faces outwards
If distPuntoPlano(s(3), basePlano) < 0 Then Swap s(1), s(2)
' Make 4 facets
Dim f(3) As Integer
For i = 0 To 3
Redim Preserve FAC(Ubound(FAC) + 1)
FAC(Ubound(FAC)) = Callocate(Sizeof(Faceta))
FAC(Ubound(FAC))->id = Ubound(FAC)
f(i) = Ubound(FAC)
Next
FAC(f(0))->p.u = s(0): FAC(f(0))->p.v = s(2): FAC(f(0))->p.w = s(1)
FAC(f(1))->p.u = s(0): FAC(f(1))->p.v = s(1): FAC(f(1))->p.w = s(3)
FAC(f(2))->p.u = s(1): FAC(f(2))->p.v = s(2): FAC(f(2))->p.w = s(3)
FAC(f(3))->p.u = s(2): FAC(f(3))->p.v = s(0): FAC(f(3))->p.w = s(3)
FAC(f(0))->n(0) = f(3): FAC(f(0))->n(1) = f(2): FAC(f(0))->n(2) = f(1)
FAC(f(1))->n(0) = f(0): FAC(f(1))->n(1) = f(2): FAC(f(1))->n(2) = f(3)
FAC(f(2))->n(0) = f(0): FAC(f(2))->n(1) = f(3): FAC(f(2))->n(2) = f(1)
FAC(f(3))->n(0) = f(0): FAC(f(3))->n(1) = f(1): FAC(f(3))->n(2) = f(2)
' Assign point set space to four faces
For i = 0 To MAXN
For j = 0 To MAXN
pts(i, j) = NULL
Next
Next
' Assign points to four faces
For i = 1 To totp
Dim v As Vect Ptr = punto(i)
If Vect_Eq(v, s(0)) Or Vect_Eq(v, s(1)) Or Vect_Eq(v, s(2)) Or Vect_Eq(v, s(3)) Then Continue For
For j = 0 To 3
If isAbove(v, FAC(f(j))->p) Then
AsignarPuntoFacetaa(f(j), v)
Exit For
End If
Next
Next
' Return the initial simplex, using a face as index
Dim hull As ConvexHulls3d Ptr = Callocate(Sizeof(ConvexHulls3d))
hull->index = f(0)
Return hull
End Function
'' The main QuickHull3D loop
Function quickHull3d(punto() As Vect Ptr, totp As Integer) As ConvexHulls3d Ptr
Dim As Integer i, j, k
Dim hull As ConvexHulls3d Ptr = getStart(punto(), totp)
' Add the face of initial simplex to queue
Redim que(0)
que(0) = hull->index
For i = 0 To 2
Redim Preserve que(Ubound(que) + 1)
que(Ubound(que)) = FAC(hull->index)->n(i)
Next
' snew saves index face of the final convex hull
Dim snew As Integer = 0
While Ubound(que) >= 0
Dim nf As Integer = que(0)
ActualizarCola()
' Skip if the current face has been deleted
If FAC(nf)->isdel Then Continue While
' Find the farthest point from the face
Dim p As Vect Ptr = NULL
For i = 0 To MAXN
If pts(nf, i) <> 0 Then
p = pts(nf, i)
Exit For
End If
Next
If p = NULL Then
If FAC(nf)->isdel = 0 Then snew = nf
Continue While
End If
' Find the horizon
G_TIME += 1
Redim resfdel(0)
resfdel(0) = 0
Redim resfdel(-1)
Redim resfnew(0)
resfnew(0) = 0
Redim resfnew(-1)
' The current face must be deleted, so start dfs from current face
Dim s As Integer = getHorizon(nf, p, resfdel())
For i = 0 To Ubound(resfdel)
Dim fid As Integer = resfdel(i)
For j = 0 To MAXN
Dim pt As Vect Ptr = pts(fid, j)
If pt = NULL Then Continue For
If Vect_Eq(pt, p) Then Continue For
' reassign them
For k = 0 To Ubound(resfnew)
Dim nfid As Integer = resfnew(k)
If isAbove(pt, FAC(nfid)->p) Then
AsignarPuntoFacetaa(nfid, pt)
Exit For
End If
Next
Next
Next
' enqueue new faces
For i = 0 To Ubound(resfnew)
Redim Preserve que(Ubound(que) + 1)
que(Ubound(que)) = resfnew(i)
Next
Wend
If snew = 0 And Ubound(FAC) > 0 Then
For i = 1 To Ubound(FAC)
If FAC(i) <> NULL And FAC(i)->isdel = 0 Then
snew = i
Exit For
End If
Next
End If
hull->index = snew
Return hull
End Function
Sub main()
Dim As Integer i
preConvexHulls()
' Example: a unit tetrahedron
Dim As Vect Ptr punto(1 To 4)
For i = 1 To 4
punto(i) = Callocate(Sizeof(Vect))
Next
punto(1)->x = 0: punto(1)->y = 0: punto(1)->z = 0
punto(2)->x = 1: punto(2)->y = 0: punto(2)->z = 0
punto(3)->x = 0: punto(3)->y = 1: punto(3)->z = 0
punto(4)->x = 0: punto(4)->y = 0: punto(4)->z = 1
Dim As ConvexHulls3d Ptr hull = quickHull3d(punto(), 4)
Print Using "##.###"; GetSurfaceArea(hull)
' Free memory
For i = 1 To 4
Deallocate(punto(i))
Next
For i = 0 To Ubound(FAC)
If FAC(i) <> NULL Then Deallocate(FAC(i))
Next
Deallocate(hull)
End Sub
main()
Sleep
- Output:
2.366
Go
package main
import (
"fmt"
"math"
)
const (
MAXN = 2500
EPS = 1e-8
)
// Globals
var (
FAC []*Facet
pts [][]*Vect
TIME int
e [2][MAXN]Edge
vistimeArr [MAXN]int
que []int
resfnew []int
resfdel []int
respt []*Vect
)
// Vect represents a 3D point/vector
type Vect struct {
x, y, z float64
id int
}
func (a *Vect) Sub(b *Vect) *Vect {
return &Vect{a.x - b.x, a.y - b.y, a.z - b.z, 0}
}
func (a *Vect) Cross(b *Vect) *Vect {
return &Vect{
a.y*b.z - a.z*b.y,
a.z*b.x - a.x*b.z,
a.x*b.y - a.y*b.x,
0,
}
}
func (a *Vect) Dot(b *Vect) float64 {
return a.x*b.x + a.y*b.y + a.z*b.z
}
func (a *Vect) Mag() float64 {
return math.Sqrt(a.x*a.x + a.y*a.y + a.z*a.z)
}
func (a *Vect) Eq(b *Vect) bool {
return eq(a.x, b.x) && eq(a.y, b.y) && eq(a.z, b.z)
}
// Line is the segment u->v
type Line struct {
u, v *Vect
}
// Plane is defined by three points u, v, w
type Plane struct {
u, v, w *Vect
}
// Normal vector of the plane
func (p *Plane) Normal() *Vect {
return p.v.Sub(p.u).Cross(p.w.Sub(p.u))
}
// vecAt returns the i-th defining point of the plane (0=u,1=v,2=w)
func (p *Plane) vecAt(i int) *Vect {
switch i {
case 0:
return p.u
case 1:
return p.v
default:
return p.w
}
}
// vecID returns the id of the i-th plane point
func (p *Plane) vecID(i int) int {
return p.vecAt(i).id
}
// Facet is a face in the hull
type Facet struct {
n [3]int // neighbor facet indices
id int
vistime int
isdel bool
p Plane
}
// Edge on the horizon
type Edge struct {
netid, facetid int
}
// ConvexHulls3d manages hull construction
type ConvexHulls3d struct {
index int
surfacearea float64
}
// Basic floating comparisons
func eq(a, b float64) bool { return math.Abs(a-b) < EPS }
func gtr(a, b float64) bool { return a-b > EPS }
func Abs(x float64) float64 { if x < 0 { return -x }; return x }
// Distances
func distPointPlane(v *Vect, p Plane) float64 {
num := v.Sub(p.u).Dot(p.Normal())
den := p.Normal().Mag()
return num / den
}
func distPointLine(v *Vect, f Line) float64 {
d := v.Sub(f.u).Mag()
if d == 0 {
return 0
}
return f.v.Sub(f.u).Cross(v.Sub(f.u)).Mag() / f.v.Sub(f.u).Mag()
}
func distPointPoint(u, v *Vect) float64 {
return u.Sub(v).Mag()
}
func isAbove(v *Vect, p Plane) bool {
return gtr(v.Sub(p.u).Dot(p.Normal()), 0)
}
// Initialize globals
func preConvexHulls() {
pts = append(pts, nil) // reserve index 0
FAC = append(FAC, &Facet{}) // dummy facet at 0
}
// DFS to accumulate area
func (h *ConvexHulls3d) dfsArea(nf int) {
if FAC[nf].vistime == TIME {
return
}
FAC[nf].vistime = TIME
nrm := FAC[nf].p.Normal()
h.surfacearea += nrm.Mag() / 2
for i := 0; i < 3; i++ {
h.dfsArea(FAC[nf].n[i])
}
}
// GetSurfaceArea returns (and caches) the hull surface area
func (h *ConvexHulls3d) GetSurfaceArea() float64 {
if gtr(h.surfacearea, 0) {
return h.surfacearea
}
TIME++
h.dfsArea(h.index)
return h.surfacearea
}
// Recursively find the horizon around point p
func (h *ConvexHulls3d) getHorizon(f int, p *Vect, resfdel *[]int) int {
if !isAbove(p, FAC[f].p) {
return 0
}
if FAC[f].vistime == TIME {
return -1
}
FAC[f].vistime = TIME
FAC[f].isdel = true
*resfdel = append(*resfdel, FAC[f].id)
ret := -2
for i := 0; i < 3; i++ {
r := h.getHorizon(FAC[f].n[i], p, resfdel)
if r == 0 {
a := FAC[f].p.vecID(i)
b := FAC[f].p.vecID((i + 1) % 3)
for j, pt := range []int{a, b} {
if vistimeArr[pt] != TIME {
vistimeArr[pt] = TIME
e[0][pt] = Edge{netid: []int{b, a}[j], facetid: FAC[f].n[i]}
} else {
e[1][pt] = Edge{netid: []int{b, a}[j], facetid: FAC[f].n[i]}
}
}
ret = a
} else if r != -1 && r != -2 {
ret = r
}
}
return ret
}
// Build the initial tetrahedron
func getStart(point []*Vect, totp int) *ConvexHulls3d {
var pt [6]*Vect
var s [4]*Vect
for i := range pt {
pt[i] = point[1]
}
// pick extreme coords
for i := 1; i <= totp; i++ {
v := point[i]
if gtr(v.x, pt[0].x) {
pt[0] = v
}
if gtr(pt[1].x, v.x) {
pt[1] = v
}
if gtr(v.y, pt[2].y) {
pt[2] = v
}
if gtr(pt[3].y, v.y) {
pt[3] = v
}
if gtr(v.z, pt[4].z) {
pt[4] = v
}
if gtr(pt[5].z, v.z) {
pt[5] = v
}
}
// take furthest pair
s[0], s[1] = pt[0], pt[1]
for i := 0; i < 6; i++ {
for j := i + 1; j < 6; j++ {
d := distPointPoint(pt[i], pt[j])
if gtr(d, distPointPoint(s[0], s[1])) {
s[0], s[1] = pt[i], pt[j]
}
}
}
// furthest from that line
L := Line{s[0], s[1]}
s[2] = pt[0]
for i := 0; i < 6; i++ {
if gtr(distPointLine(pt[i], L), distPointLine(s[2], L)) {
s[2] = pt[i]
}
}
// furthest from the plane
s[3] = point[1]
base := Plane{s[0], s[1], s[2]}
for i := 1; i <= totp; i++ {
d1 := Abs(distPointPlane(point[i], base))
d2 := Abs(distPointPlane(s[3], base))
if gtr(d1, d2) {
s[3] = point[i]
}
}
if gtr(0, distPointPlane(s[3], base)) {
s[1], s[2] = s[2], s[1]
}
// make 4 facets
f := make([]int, 4)
for i := 0; i < 4; i++ {
FAC = append(FAC, &Facet{id: len(FAC)})
f[i] = len(FAC) - 1
}
FAC[f[0]].p = Plane{s[0], s[2], s[1]}
FAC[f[1]].p = Plane{s[0], s[1], s[3]}
FAC[f[2]].p = Plane{s[1], s[2], s[3]}
FAC[f[3]].p = Plane{s[2], s[0], s[3]}
FAC[f[0]].n = [3]int{f[3], f[2], f[1]}
FAC[f[1]].n = [3]int{f[0], f[2], f[3]}
FAC[f[2]].n = [3]int{f[0], f[3], f[1]}
FAC[f[3]].n = [3]int{f[0], f[1], f[2]}
// prepare point‐lists
for i := 0; i < 4; i++ {
pts = append(pts, nil)
}
// assign remaining points
for i := 1; i <= totp; i++ {
v := point[i]
if v.Eq(s[0]) || v.Eq(s[1]) || v.Eq(s[2]) || v.Eq(s[3]) {
continue
}
for j := 0; j < 4; j++ {
if isAbove(v, FAC[f[j]].p) {
pts[f[j]] = append(pts[f[j]], v)
break
}
}
}
return &ConvexHulls3d{index: f[0]}
}
// The main QuickHull3D loop
func quickHull3d(point []*Vect, totp int) *ConvexHulls3d {
hull := getStart(point, totp)
que = []int{hull.index}
for i := 0; i < 3; i++ {
que = append(que, FAC[hull.index].n[i])
}
snew := 0
for len(que) > 0 {
nf := que[0]
que = que[1:]
if FAC[nf].isdel || len(pts[nf]) == 0 {
if !FAC[nf].isdel {
snew = nf
}
continue
}
// find farthest point
p := pts[nf][0]
for _, v := range pts[nf][1:] {
if gtr(distPointPlane(v, FAC[nf].p), distPointPlane(p, FAC[nf].p)) {
p = v
}
}
// find horizon
TIME++
resfdel = resfdel[:0]
s := hull.getHorizon(nf, p, &resfdel)
// build new faces around the horizon
resfnew = resfnew[:0]
TIME++
from := -1
lastf := -1
fstf := -1
for vistimeArr[s] != TIME {
vistimeArr[s] = TIME
var net, f, fnew int
if e[0][s].netid == from {
net, f = e[1][s].netid, e[1][s].facetid
} else {
net, f = e[0][s].netid, e[0][s].facetid
}
// find indices of s and net on facet f
pt1, pt2 := 0, 0
for i := 0; i < 3; i++ {
if point[s].Eq(FAC[f].p.vecAt(i)) {
pt1 = i
}
if point[net].Eq(FAC[f].p.vecAt(i)) {
pt2 = i
}
}
if (pt1+1)%3 != pt2 {
pt1, pt2 = pt2, pt1
}
// create new facet facing outwards
FAC = append(FAC, &Facet{
id: len(FAC),
p: Plane{FAC[f].p.vecAt(pt2), FAC[f].p.vecAt(pt1), p},
})
fnew = len(FAC) - 1
pts = append(pts, nil)
resfnew = append(resfnew, fnew)
FAC[fnew].n[0] = f
FAC[f].n[pt1] = fnew
if lastf >= 0 {
// link with previous new facet
if FAC[fnew].p.vecAt(1).Eq(FAC[lastf].p.u) {
FAC[fnew].n[1], FAC[lastf].n[2] = lastf, fnew
} else {
FAC[fnew].n[2], FAC[lastf].n[1] = lastf, fnew
}
} else {
fstf = fnew
}
lastf = fnew
from = s
s = net
}
// close the loop
if FAC[fstf].p.v.Eq(FAC[lastf].p.u) {
FAC[fstf].n[1], FAC[lastf].n[2] = lastf, fstf
} else {
FAC[fstf].n[2], FAC[lastf].n[1] = lastf, fstf
}
// collect points from deleted faces
respt = respt[:0]
for _, fid := range resfdel {
respt = append(respt, pts[fid]...)
pts[fid] = nil
}
// reassign them
for _, v := range respt {
if v == p {
continue
}
for _, fid := range resfnew {
if isAbove(v, FAC[fid].p) {
pts[fid] = append(pts[fid], v)
break
}
}
}
// enqueue new faces
for _, fid := range resfnew {
que = append(que, fid)
}
}
hull.index = snew
return hull
}
func main() {
preConvexHulls()
// Example: a unit tetrahedron
point := []*Vect{
nil,
{0, 0, 0, 1},
{1, 0, 0, 2},
{0, 1, 0, 3},
{0, 0, 1, 4},
}
hull := quickHull3d(point, 4)
fmt.Printf("%.3f\n", hull.GetSurfaceArea())
}
- Output:
2.366
Java
import java.util.*;
public class QuickHull3D {
static final int MAXN = 2500;
static final double EPS = 1e-8;
// Globals
static List<Facet> FAC = new ArrayList<>();
static List<List<Vect>> pts = new ArrayList<>();
static int TIME = 0;
static Edge[][] e = new Edge[2][MAXN];
static int[] vistime = new int[MAXN];
static List<Integer> queue = new ArrayList<>();
static List<Integer> resfnew = new ArrayList<>();
static List<Integer> resfdel = new ArrayList<>();
static List<Vect> respt = new ArrayList<>();
public static void main(String[] args) {
preConvexHulls();
// Example: unit tetrahedron
int n = 4;
Vect[] point = new Vect[n + 1];
point[1] = new Vect(0, 0, 0, 1);
point[2] = new Vect(1, 0, 0, 2);
point[3] = new Vect(0, 1, 0, 3);
point[4] = new Vect(0, 0, 1, 4);
ConvexHulls3d hull = quickHull3D(point, n);
System.out.printf("%.3f%n", hull.getSurfaceArea());
}
static void preConvexHulls() {
// Reserve index 0
pts.add(new ArrayList<>()); // dummy point list[0]
FAC.add(new Facet()); // dummy facet[0]
// Initialize edge array
for (int i = 0; i < 2; i++) {
for (int j = 0; j < MAXN; j++) {
e[i][j] = new Edge();
}
}
}
static class Vect {
double X, Y, Z;
int Id;
Vect(double x, double y, double z, int id) { X = x; Y = y; Z = z; Id = id; }
Vect sub(Vect b) { return new Vect(X - b.X, Y - b.Y, Z - b.Z, 0); }
Vect cross(Vect b) {
return new Vect(
Y * b.Z - Z * b.Y,
Z * b.X - X * b.Z,
X * b.Y - Y * b.X,
0
);
}
double dot(Vect b) { return X * b.X + Y * b.Y + Z * b.Z; }
double mag() { return Math.sqrt(X * X + Y * Y + Z * Z); }
boolean eq(Vect b) {
// Qualify eq to the outer class's static method:
return QuickHull3D.eq(X, b.X)
&& QuickHull3D.eq(Y, b.Y)
&& QuickHull3D.eq(Z, b.Z);
}
}
static class Line {
Vect U, V;
Line(Vect u, Vect v) { U = u; V = v; }
}
static class Plane {
Vect U, V, W;
Plane() {}
Plane(Vect u, Vect v, Vect w) { U = u; V = v; W = w; }
Vect normal() { return V.sub(U).cross(W.sub(U)); }
Vect vecAt(int i) {
if (i == 0) return U;
if (i == 1) return V;
return W;
}
int vecId(int i) { return vecAt(i).Id; }
}
static class Facet {
int[] N = new int[3];
int Id;
int visitTime;
boolean isDeleted;
Plane P = new Plane();
}
static class Edge {
int netId, facetId;
}
static class ConvexHulls3d {
int Index;
double surfaceArea = 0.0;
ConvexHulls3d(int idx) { Index = idx; }
void dfsArea(int f) {
Facet F = FAC.get(f);
if (F.visitTime == TIME) return;
F.visitTime = TIME;
Vect nrm = F.P.normal();
surfaceArea += nrm.mag() / 2.0;
for (int i = 0; i < 3; i++) {
dfsArea(F.N[i]);
}
}
public double getSurfaceArea() {
if (gtr(surfaceArea, 0.0)) return surfaceArea;
TIME++;
dfsArea(Index);
return surfaceArea;
}
public int getHorizon(int f, Vect p, List<Integer> resDel) {
Facet Ff = FAC.get(f);
if (!isAbove(p, Ff.P)) return 0;
if (Ff.visitTime == TIME) return -1;
Ff.visitTime = TIME;
Ff.isDeleted = true;
resDel.add(Ff.Id);
int ret = -2;
for (int i = 0; i < 3; i++) {
int ni = Ff.N[i];
int r = getHorizon(ni, p, resDel);
if (r == 0) {
int a = Ff.P.vecId(i);
int b = Ff.P.vecId((i + 1) % 3);
for (int idx = 0; idx < 2; idx++) {
int pt = (idx == 0 ? a : b);
int facet = ni;
if (vistime[pt] != TIME) {
vistime[pt] = TIME;
e[0][pt].netId = (idx == 0 ? b : a);
e[0][pt].facetId = facet;
} else {
e[1][pt].netId = (idx == 0 ? b : a);
e[1][pt].facetId = facet;
}
}
ret = a;
} else if (r != -1 && r != -2) {
ret = r;
}
}
return ret;
}
}
// Utilities
static boolean eq(double a, double b) { return Math.abs(a - b) < EPS; }
static boolean gtr(double a, double b) { return a - b > EPS; }
static double abs(double x) { return x < 0 ? -x : x; }
static double distPointPlane(Vect v, Plane p) {
double num = v.sub(p.U).dot(p.normal());
double den = p.normal().mag();
return num / den;
}
static double distPointLine(Vect v, Line l) {
double d = v.sub(l.U).mag();
if (d == 0) return 0;
return l.V.sub(l.U).cross(v.sub(l.U)).mag() / l.V.sub(l.U).mag();
}
static double distPointPoint(Vect a, Vect b) { return a.sub(b).mag(); }
static boolean isAbove(Vect v, Plane p) {
return gtr(v.sub(p.U).dot(p.normal()), 0);
}
static ConvexHulls3d getStart(Vect[] point, int totp) {
Vect[] extremes = new Vect[6];
for (int i = 0; i < 6; i++) extremes[i] = point[1];
for (int i = 1; i <= totp; i++) {
Vect v = point[i];
if (gtr(v.X, extremes[0].X)) extremes[0] = v;
if (gtr(extremes[1].X, v.X)) extremes[1] = v;
if (gtr(v.Y, extremes[2].Y)) extremes[2] = v;
if (gtr(extremes[3].Y, v.Y)) extremes[3] = v;
if (gtr(v.Z, extremes[4].Z)) extremes[4] = v;
if (gtr(extremes[5].Z, v.Z)) extremes[5] = v;
}
// Furthest pair
Vect s0 = extremes[0], s1 = extremes[1];
for (int i = 0; i < 6; i++) {
for (int j = i + 1; j < 6; j++) {
double d = distPointPoint(extremes[i], extremes[j]);
if (gtr(d, distPointPoint(s0, s1))) {
s0 = extremes[i];
s1 = extremes[j];
}
}
}
// Furthest from line
Line line = new Line(s0, s1);
Vect s2 = extremes[0];
for (int i = 0; i < 6; i++) {
if (gtr(distPointLine(extremes[i], line),
distPointLine(s2, line)))
s2 = extremes[i];
}
// Furthest from plane
Vect s3 = point[1];
Plane basePlane = new Plane(s0, s1, s2);
for (int i = 1; i <= totp; i++) {
double d1 = abs(distPointPlane(point[i], basePlane));
double d2 = abs(distPointPlane(s3, basePlane));
if (gtr(d1, d2)) s3 = point[i];
}
if (gtr(0, distPointPlane(s3, basePlane))) {
Vect tmp = s1; s1 = s2; s2 = tmp;
}
// Create 4 initial facets
int[] f = new int[4];
for (int i = 0; i < 4; i++) {
Facet F = new Facet();
F.Id = FAC.size();
FAC.add(F);
f[i] = FAC.size() - 1;
}
FAC.get(f[0]).P = new Plane(s0, s2, s1);
FAC.get(f[1]).P = new Plane(s0, s1, s3);
FAC.get(f[2]).P = new Plane(s1, s2, s3);
FAC.get(f[3]).P = new Plane(s2, s0, s3);
FAC.get(f[0]).N = new int[]{f[3], f[2], f[1]};
FAC.get(f[1]).N = new int[]{f[0], f[2], f[3]};
FAC.get(f[2]).N = new int[]{f[0], f[3], f[1]};
FAC.get(f[3]).N = new int[]{f[0], f[1], f[2]};
// Prepare pts lists for the 4 facets
for (int i = 0; i < 4; i++) {
pts.add(new ArrayList<>());
}
// Assign points to facets
for (int i = 1; i <= totp; i++) {
Vect v = point[i];
if (v.eq(s0) || v.eq(s1) || v.eq(s2) || v.eq(s3)) continue;
for (int j = 0; j < 4; j++) {
if (isAbove(v, FAC.get(f[j]).P)) {
pts.get(f[j]).add(v);
break;
}
}
}
return new ConvexHulls3d(f[0]);
}
static ConvexHulls3d quickHull3D(Vect[] point, int totp) {
ConvexHulls3d hull = getStart(point, totp);
queue.clear();
queue.add(hull.Index);
for (int i = 0; i < 3; i++) {
queue.add(FAC.get(hull.Index).N[i]);
}
int snew = 0;
while (!queue.isEmpty()) {
int nf = queue.remove(0);
Facet Fnf = FAC.get(nf);
if (Fnf.isDeleted || pts.get(nf).isEmpty()) {
if (!Fnf.isDeleted) snew = nf;
continue;
}
// Farthest point from facet plane
Vect p = pts.get(nf).get(0);
for (Vect v : pts.get(nf)) {
if (gtr(distPointPlane(v, Fnf.P), distPointPlane(p, Fnf.P))) {
p = v;
}
}
// Find horizon
TIME++;
resfdel.clear();
int s = hull.getHorizon(nf, p, resfdel);
// Build new faces
resfnew.clear();
TIME++;
int from = -1, lastf = -1, fstf = -1;
while (vistime[s] != TIME) {
vistime[s] = TIME;
int net, fidx;
if (e[0][s].netId == from) {
net = e[1][s].netId;
fidx = e[1][s].facetId;
} else {
net = e[0][s].netId;
fidx = e[0][s].facetId;
}
// Find indices on facet fidx
Facet Ff = FAC.get(fidx);
int pt1 = 0, pt2 = 0;
for (int i = 0; i < 3; i++) {
if (point[s].eq(Ff.P.vecAt(i))) pt1 = i;
if (point[net].eq(Ff.P.vecAt(i))) pt2 = i;
}
if ((pt1 + 1) % 3 != pt2) {
int tmp = pt1; pt1 = pt2; pt2 = tmp;
}
// Create new facet
Facet newF = new Facet();
newF.Id = FAC.size();
newF.P = new Plane(
Ff.P.vecAt(pt2),
Ff.P.vecAt(pt1),
p
);
FAC.add(newF);
pts.add(new ArrayList<>());
int fnew = FAC.size() - 1;
resfnew.add(fnew);
// Link neighborhoods
newF.N[0] = fidx;
Ff.N[pt1] = fnew;
if (lastf >= 0) {
Plane Pf1 = newF.P, Plast = FAC.get(lastf).P;
if (Pf1.vecAt(1).eq(Plast.U)) {
newF.N[1] = lastf;
FAC.get(lastf).N[2] = fnew;
} else {
newF.N[2] = lastf;
FAC.get(lastf).N[1] = fnew;
}
} else {
fstf = fnew;
}
lastf = fnew;
from = s;
s = net;
}
// Close the loop
Facet Ffst = FAC.get(fstf), Flast = FAC.get(lastf);
if (Ffst.P.vecAt(1).eq(Flast.P.U)) {
Ffst.N[1] = lastf;
Flast.N[2] = fstf;
} else {
Ffst.N[2] = lastf;
Flast.N[1] = fstf;
}
// Collect deleted points
respt.clear();
for (int fid : resfdel) {
respt.addAll(pts.get(fid));
pts.get(fid).clear();
}
// Reassign points
for (Vect v : respt) {
if (v == p) continue;
for (int fid : resfnew) {
if (isAbove(v, FAC.get(fid).P)) {
pts.get(fid).add(v);
break;
}
}
}
// Enqueue new faces
for (int fid : resfnew) {
queue.add(fid);
}
}
hull.Index = snew;
return hull;
}
}
- Output:
2.366
JavaScript
// Constants
const MAXN = 2500;
const EPS = 1e-8;
function gtr(a, b) {
return a - b > EPS;
}
function eq(a, b) {
return Math.abs(a - b) < EPS;
}
function Abs(x) {
return x < 0 ? -x : x;
}
// 3D Vector
class Vect {
constructor(x = 0.0, y = 0.0, z = 0.0, id = 0) {
this.x = x;
this.y = y;
this.z = z;
this.id = id;
}
subtract(o) {
return new Vect(this.x - o.x, this.y - o.y, this.z - o.z);
}
// cross product
cross(o) {
return new Vect(
this.y * o.z - this.z * o.y,
this.z * o.x - this.x * o.z,
this.x * o.y - this.y * o.x
);
}
// dot product
dot(o) {
return this.x * o.x + this.y * o.y + this.z * o.z;
}
magnitude() {
return Math.sqrt(this.x * this.x + this.y * this.y + this.z * this.z);
}
equals(o) {
return eq(this.x, o.x) && eq(this.y, o.y) && eq(this.z, o.z);
}
toString() {
return `Vect(${this.x}, ${this.y}, ${this.z}, ${this.id})`;
}
}
// Line from u to v
class Line {
constructor(u = null, v = null) {
this.u = u;
this.v = v;
}
}
// Plane through three points u,v,w
class Plane {
constructor(u = null, v = null, w = null) {
this.vec = (u && v && w) ? [u, v, w] : [null, null, null];
}
// normal vector = (v - u) × (w - u)
normal() {
const [u, v, w] = this.vec;
return v.subtract(u).cross(w.subtract(u));
}
// one point on plane
u() {
return this.vec[0];
}
}
// Signed distance point-plane
function distPointPlane(pt, pl) {
const n = pl.normal();
return pt.subtract(pl.u()).dot(n) / n.magnitude();
}
// Unsigned distance point-line
function distPointLine(pt, ln) {
const u_to_pt = pt.subtract(ln.u);
const dir = ln.v.subtract(ln.u);
if (u_to_pt.magnitude() === 0) return 0;
return u_to_pt.cross(dir).magnitude() / dir.magnitude();
}
// Distance point-point
function distPointPoint(a, b) {
return a.subtract(b).magnitude();
}
// is pt strictly above plane?
function isAbove(pt, pl) {
const n = pl.normal();
return gtr(pt.subtract(pl.u()).dot(n), 0);
}
// Facet (triangle face of hull)
class Facet {
constructor(id = 0, p = null) {
this.n = [0, 0, 0]; // neighbors
this.id = id;
this.vistime = 0;
this.isdel = false;
this.p = p || new Plane();
}
setNeighbors(n1, n2, n3) {
this.n = [n1, n2, n3];
}
toString() {
return `Facet(id=${this.id}, isdel=${this.isdel}, vistime=${this.vistime})`;
}
}
// Edge info for horizon detection
class Edge {
constructor() {
this.netid = 0;
this.facetid = 0;
}
}
// Global variables
let FAC = [];
let pts = []; // points assigned to each face
let TIME = 0;
// Convex hull wrapper
class ConvexHulls3d {
constructor(index) {
this.index = index;
this.surfaceArea = 0.0;
}
dfsArea(fidx) {
if (FAC[fidx].vistime === TIME) return;
FAC[fidx].vistime = TIME;
const n = FAC[fidx].p.normal();
this.surfaceArea += n.magnitude() / 2;
for (let i = 0; i < 3; i++) {
this.dfsArea(FAC[fidx].n[i]);
}
}
getSurfaceArea() {
if (gtr(this.surfaceArea, 0)) {
return this.surfaceArea;
}
TIME++;
this.dfsArea(this.index);
return this.surfaceArea;
}
getHorizon(f, p, vistime, e1, e2, resfdel) {
if (!isAbove(p, FAC[f].p)) return 0;
if (FAC[f].vistime === TIME) return -1;
FAC[f].vistime = TIME;
FAC[f].isdel = true;
resfdel.push(FAC[f].id);
let ret = -2;
for (let i = 0; i < 3; i++) {
const res = this.getHorizon(FAC[f].n[i], p, vistime, e1, e2, resfdel);
if (res === 0) {
const ptIds = [
FAC[f].p.vec[i].id,
FAC[f].p.vec[(i + 1) % 3].id
];
for (let j = 0; j < 2; j++) {
const pid = ptIds[j];
if (vistime[pid] !== TIME) {
vistime[pid] = TIME;
e1[pid].netid = ptIds[(j + 1) % 2];
e1[pid].facetid = FAC[f].n[i];
} else {
e2[pid].netid = ptIds[(j + 1) % 2];
e2[pid].facetid = FAC[f].n[i];
}
}
ret = ptIds[0];
} else if (res !== -1 && res !== -2) {
ret = res;
}
}
return ret;
}
}
// Initialize convex hull system
function preConvexHulls() {
pts = [[]];
FAC = [ new Facet() ];
}
// Build initial tetrahedron (simplex)
function getStart(point, totp) {
// pick extreme points along axes
let pt = Array(6).fill(point[1]);
let s = Array(4).fill(point[1]);
for (let i = 2; i <= totp; i++) {
if (gtr(point[i].x, pt[0].x)) pt[0] = point[i];
if (gtr(pt[1].x, point[i].x)) pt[1] = point[i];
if (gtr(point[i].y, pt[2].y)) pt[2] = point[i];
if (gtr(pt[3].y, point[i].y)) pt[3] = point[i];
if (gtr(point[i].z, pt[4].z)) pt[4] = point[i];
if (gtr(pt[5].z, point[i].z)) pt[5] = point[i];
}
// farthest pair among these 6
for (let i = 0; i < 6; i++) {
for (let j = i + 1; j < 6; j++) {
if (gtr(distPointPoint(pt[i], pt[j]), distPointPoint(s[0], s[1]))) {
s[0] = pt[i];
s[1] = pt[j];
}
}
}
// farthest point from line s[0]-s[1]
for (let i = 0; i < 6; i++) {
if (gtr(
distPointLine(pt[i], new Line(s[0], s[1])),
distPointLine(s[2], new Line(s[0], s[1]))
)) {
s[2] = pt[i];
}
}
// farthest point from plane through s[0],s[1],s[2]
for (let i = 1; i <= totp; i++) {
if (gtr(
Abs(distPointPlane(point[i], new Plane(s[0], s[1], s[2]))),
Abs(distPointPlane(s[3], new Plane(s[0], s[1], s[2])))
)) {
s[3] = point[i];
}
}
// ensure correct orientation
if (gtr(0, distPointPlane(s[3], new Plane(s[0], s[1], s[2])))) {
[s[1], s[2]] = [s[2], s[1]];
}
// create 4 faces
let f = new Array(4);
for (let i = 0; i < 4; i++) {
FAC.push(new Facet(FAC.length));
f[i] = FAC.length - 1;
}
FAC[f[0]].p = new Plane(s[0], s[2], s[1]);
FAC[f[1]].p = new Plane(s[0], s[1], s[3]);
FAC[f[2]].p = new Plane(s[1], s[2], s[3]);
FAC[f[3]].p = new Plane(s[2], s[0], s[3]);
FAC[f[0]].setNeighbors(f[3], f[2], f[1]);
FAC[f[1]].setNeighbors(f[0], f[2], f[3]);
FAC[f[2]].setNeighbors(f[0], f[3], f[1]);
FAC[f[3]].setNeighbors(f[0], f[1], f[2]);
// allocate point lists
for (let i = 0; i < 4; i++) pts.push([]);
// assign remaining points to faces
for (let i = 1; i <= totp; i++) {
const pi = point[i];
if (pi.equals(s[0]) || pi.equals(s[1]) || pi.equals(s[2]) || pi.equals(s[3])) continue;
for (let j = 0; j < 4; j++) {
if (isAbove(pi, FAC[f[j]].p)) {
pts[f[j]].push(pi);
break;
}
}
}
return new ConvexHulls3d(f[0]);
}
// QuickHull main routine
function quickHull3d(point, totp) {
let hull = getStart(point, totp);
let queue = [hull.index,
...FAC[hull.index].n];
let snew = 0;
// horizon data
const e = [ Array(MAXN).fill().map(() => new Edge()),
Array(MAXN).fill().map(() => new Edge()) ];
const vistime = Array(MAXN).fill(0);
while (queue.length > 0) {
const nf = queue.shift();
if (FAC[nf].isdel) continue;
if (pts[nf].length === 0) {
snew = nf;
continue;
}
// farthest point from face
let p = pts[nf][0];
for (let pt of pts[nf]) {
if (gtr(distPointPlane(pt, FAC[nf].p), distPointPlane(p, FAC[nf].p))) {
p = pt;
}
}
// find horizon
TIME++;
let resfdel = [];
const startEdge = hull.getHorizon(nf, p, vistime, e[0], e[1], resfdel);
// build new faces around horizon
let resfnew = [];
TIME++;
let from = 0, lastf = 0, fstf = 0, s = startEdge;
while (vistime[s] !== TIME) {
vistime[s] = TIME;
let net, adjF;
if (e[0][s].netid === from) {
net = e[1][s].netid; adjF = e[1][s].facetid;
} else {
net = e[0][s].netid; adjF = e[0][s].facetid;
}
// find indices of s and net in adj face
let pt1 = -1, pt2 = -1;
for (let i = 0; i < 3; i++) {
if (point[s] === FAC[adjF].p.vec[i]) pt1 = i;
if (point[net] === FAC[adjF].p.vec[i]) pt2 = i;
}
if ((pt1 + 1) % 3 !== pt2) [pt1, pt2] = [pt2, pt1];
// new face [pt2, pt1, p]
FAC.push(new Facet(FAC.length, new Plane(
FAC[adjF].p.vec[pt2],
FAC[adjF].p.vec[pt1],
p
)));
const fnew = FAC.length - 1;
pts.push([]);
resfnew.push(fnew);
// link adjacency
FAC[fnew].n[0] = adjF;
FAC[adjF].n[pt1] = fnew;
if (lastf) {
if (FAC[fnew].p.vec[1] === FAC[lastf].p.vec[0]) {
FAC[fnew].n[1] = lastf;
FAC[lastf].n[2] = fnew;
} else {
FAC[fnew].n[2] = lastf;
FAC[lastf].n[1] = fnew;
}
} else {
fstf = fnew;
}
lastf = fnew;
from = s;
s = net;
}
// close the loop
if (FAC[fstf].p.vec[1] === FAC[lastf].p.vec[0]) {
FAC[fstf].n[1] = lastf;
FAC[lastf].n[2] = fstf;
} else {
FAC[fstf].n[2] = lastf;
FAC[lastf].n[1] = fstf;
}
// collect and reassign points
let respt = [];
for (let fid of resfdel) {
respt.push(...pts[fid]);
pts[fid] = [];
}
for (let qpt of respt) {
if (qpt === p) continue;
for (let nfnew of resfnew) {
if (isAbove(qpt, FAC[nfnew].p)) {
pts[nfnew].push(qpt);
break;
}
}
}
// enqueue new faces
queue.push(...resfnew);
}
hull.index = snew;
return hull;
}
// Example usage
(function main() {
preConvexHulls();
const n = 4;
const input = [
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0]
];
const point = Array(n + 1);
for (let i = 1; i <= n; i++) {
const [x, y, z] = input[i - 1];
point[i] = new Vect(x, y, z, i);
}
const hull = quickHull3d(point, n);
console.log(hull.getSurfaceArea().toFixed(3));
})();
- Output:
2.366
Julia
# QuickHull3D.jl
using Printf
const MAXN = 2500
const EPS = 1e-8
# Tell Julia we intend to extend these Base methods
import Base: -, ==
# 3D vector with an ID
mutable struct Vect
x::Float64; y::Float64; z::Float64; id::Int
end
# Subtract two Vect’s
function -(a::Vect, b::Vect)
Vect(a.x - b.x, a.y - b.y, a.z - b.z, 0)
end
# Equality of two Vect’s (approximate)
function ==(a::Vect, b::Vect)
abs(a.x-b.x)<EPS && abs(a.y-b.y)<EPS && abs(a.z-b.z)<EPS
end
# Cross, dot, magnitude
function cross(a::Vect, b::Vect)
Vect(
a.y*b.z - a.z*b.y,
a.z*b.x - a.x*b.z,
a.x*b.y - a.y*b.x,
0
)
end
dot(a::Vect, b::Vect) = a.x*b.x + a.y*b.y + a.z*b.z
mag(a::Vect) = sqrt(a.x^2 + a.y^2 + a.z^2)
# Line and Plane types
struct Line
u::Vect; v::Vect
end
struct Plane
u::Vect; v::Vect; w::Vect
end
normal(p::Plane) = cross(p.v - p.u, p.w - p.u)
vecAt(p::Plane, i::Int) = (i==1 ? p.u : i==2 ? p.v : p.w)
vecId(p::Plane, i::Int) = vecAt(p, i).id
# Floating‐point comparison
gtr(a,b) = a - b > EPS
# Distances
dist_point_plane(v::Vect, p::Plane) = dot(v - p.u, normal(p)) / mag(normal(p))
dist_point_line(v::Vect, l::Line) = begin
d = mag(v - l.u)
d>0 ? mag(cross(l.v - l.u, v - l.u)) / mag(l.v - l.u) : 0.0
end
dist_point_point(a::Vect, b::Vect) = mag(a - b)
isabove(v::Vect, p::Plane) = gtr(dot(v - p.u, normal(p)), 0)
# Facet and Edge
mutable struct Facet
n::NTuple{3,Int} # neighbor facet indices
id::Int # facet ID
vistime::Int # timestamp for DFS
isdel::Bool # deleted?
p::Plane # supporting plane
end
Facet() = Facet((0,0,0),0,0,false,Plane(Vect(0,0,0,0),Vect(0,0,0,0),Vect(0,0,0,0)))
mutable struct Edge
netid::Int
facetid::Int
end
Edge() = Edge(0,0)
mutable struct ConvexHulls3d
index::Int
surfacearea::Float64
end
ConvexHulls3d(idx::Int) = ConvexHulls3d(idx, 0.0)
# Globals
FAC = Facet[]
pts = Vector{Vect}[]
TIME = 0
e = [ Edge() for i in 1:2, j in 1:MAXN ]
vistime = zeros(Int, MAXN)
queue = Int[]
resfnew = Int[]
resfdel = Int[]
respt = Vect[]
# Initialize
function preConvexHulls()
push!(pts, Vect[]) # reserve index 1
push!(FAC, Facet()) # dummy facet[1]
end
# DFS for surface area
function dfsArea(h::ConvexHulls3d, fidx::Int)
f = FAC[fidx]
if f.vistime == TIME return end
f.vistime = TIME
nrm = normal(f.p)
h.surfacearea += mag(nrm)/2
for i in 1:3
dfsArea(h, f.n[i])
end
end
function getSurfaceArea(h::ConvexHulls3d)
if h.surfacearea > 0
return h.surfacearea
end
global TIME
TIME += 1
dfsArea(h, h.index)
return h.surfacearea
end
# Horizon search
function getHorizon(h::ConvexHulls3d, fidx::Int, p::Vect, resdel::Vector{Int})
f = FAC[fidx]
if !isabove(p, f.p) return 0 end
if f.vistime == TIME return -1 end
f.vistime = TIME
f.isdel = true
push!(resdel, fidx)
ret = -2
for i in 1:3
r = getHorizon(h, f.n[i], p, resdel)
if r == 0
a = vecId(f.p, i)
b = vecId(f.p, (i % 3) + 1)
for (j, pt) in enumerate((a,b))
slot = vistime[pt] != TIME ? 1 : 2
vistime[pt] = TIME
e[slot, pt] = Edge(j==1 ? b : a, f.n[i])
end
ret = a
elseif r != -1 && r != -2
ret = r
end
end
return ret
end
# Build initial tetrahedron
function getStart(point::Vector{Vect}, totp::Int)
# pick extremes
pt = [ point[1] for _ in 1:6 ]
for i in 1:totp
v = point[i]
if gtr(v.x, pt[1].x) pt[1] = v end
if gtr(pt[2].x, v.x) pt[2] = v end
if gtr(v.y, pt[3].y) pt[3] = v end
if gtr(pt[4].y, v.y) pt[4] = v end
if gtr(v.z, pt[5].z) pt[5] = v end
if gtr(pt[6].z, v.z) pt[6] = v end
end
# furthest pair
s = [pt[1], pt[2], point[1], point[1]]
for i in 1:6, j in i+1:6
d = dist_point_point(pt[i], pt[j])
if gtr(d, dist_point_point(s[1], s[2]))
s[1], s[2] = pt[i], pt[j]
end
end
# furthest from line
L = Line(s[1], s[2])
s[3] = pt[1]
for i in 1:6
if gtr(dist_point_line(pt[i], L), dist_point_line(s[3], L))
s[3] = pt[i]
end
end
# furthest from plane
s[4] = point[1]
base = Plane(s[1], s[2], s[3])
for i in 1:totp
if gtr(abs(dist_point_plane(point[i], base)),
abs(dist_point_plane(s[4], base)))
s[4] = point[i]
end
end
if gtr(0, dist_point_plane(s[4], base))
s[2], s[3] = s[3], s[2]
end
# create 4 facets
fidx = Int[]
for i in 1:4
push!(FAC, Facet())
FAC[end].id = length(FAC)
push!(fidx, length(FAC))
end
FAC[fidx[1]].p = Plane(s[1], s[3], s[2])
FAC[fidx[2]].p = Plane(s[1], s[2], s[4])
FAC[fidx[3]].p = Plane(s[2], s[3], s[4])
FAC[fidx[4]].p = Plane(s[3], s[1], s[4])
FAC[fidx[1]].n = (fidx[4], fidx[3], fidx[2])
FAC[fidx[2]].n = (fidx[1], fidx[3], fidx[4])
FAC[fidx[3]].n = (fidx[1], fidx[4], fidx[2])
FAC[fidx[4]].n = (fidx[1], fidx[2], fidx[3])
# prepare pts lists
for _ in 1:4
push!(pts, Vect[])
end
# assign points
for i in 1:totp
v = point[i]
if v == s[1] || v == s[2] || v == s[3] || v == s[4]
continue
end
for j in 1:4
if isabove(v, FAC[fidx[j]].p)
push!(pts[fidx[j]], v)
break
end
end
end
return ConvexHulls3d(fidx[1])
end
# Main QuickHull3D
function quickHull3d(point::Vector{Vect}, totp::Int)
hull = getStart(point, totp)
empty!(queue)
push!(queue, hull.index)
for i in 1:3
push!(queue, FAC[hull.index].n[i])
end
snew = hull.index
while !isempty(queue)
nf = popfirst!(queue)
if FAC[nf].isdel || isempty(pts[nf])
if !FAC[nf].isdel
snew = nf
end
continue
end
# farthest point
p = pts[nf][1]
for v in pts[nf]
if gtr(dist_point_plane(v, FAC[nf].p),
dist_point_plane(p, FAC[nf].p))
p = v
end
end
# find horizon
TIME += 1
empty!(resfdel)
s = getHorizon(hull, nf, p, resfdel)
# build around horizon
empty!(resfnew)
TIME += 1
from = -1; lastf = -1; fstf = -1
while vistime[s] != TIME
vistime[s] = TIME
e1, e2 = e[1,s], e[2,s]
net, adj = (e1.netid==from ? (e2.netid, e2.facetid)
: (e1.netid, e1.facetid))
# find edge orientation
pt1 = findfirst(i->vecAt(FAC[adj].p,i)==point[s], 1:3)
pt2 = findfirst(i->vecAt(FAC[adj].p,i)==point[net], 1:3)
if (pt1 % 3) + 1 != pt2
pt1, pt2 = pt2, pt1
end
newp = Plane(vecAt(FAC[adj].p,pt2),
vecAt(FAC[adj].p,pt1),
p)
push!(FAC, Facet(length(FAC)+1, newp))
fnew = length(FAC)
push!(pts, Vect[])
push!(resfnew, fnew)
FAC[fnew].n = (adj,0,0)
FAC[adj].n = setindex(FAC[adj].n, fnew, pt1)
if lastf != -1
if vecAt(FAC[fnew].p,2) == vecAt(FAC[lastf].p,1)
FAC[fnew].n = (FAC[fnew].n[1], lastf, FAC[fnew].n[3])
FAC[lastf].n = (FAC[lastf].n[1], FAC[lastf].n[2], fnew)
else
FAC[fnew].n = (FAC[fnew].n[1], FAC[fnew].n[2], lastf)
FAC[lastf].n = (FAC[lastf].n[1], fnew, FAC[lastf].n[3])
end
else
fstf = fnew
end
lastf, from, s = fnew, s, net
end
# close loop
if vecAt(FAC[fstf].p,2) == vecAt(FAC[lastf].p,1)
FAC[fstf].n = (FAC[fstf].n[1], lastf, FAC[fstf].n[3])
FAC[lastf].n = (FAC[lastf].n[1], FAC[lastf].n[2], fstf)
else
FAC[fstf].n = (FAC[fstf].n[1], FAC[fstf].n[2], lastf)
FAC[lastf].n = (FAC[lastf].n[1], fstf, FAC[lastf].n[3])
end
# collect and reassign points
empty!(respt)
for fid in resfdel
append!(respt, pts[fid])
empty!(pts[fid])
end
for v in respt
if v != p
for fid in resfnew
if isabove(v, FAC[fid].p)
push!(pts[fid], v)
break
end
end
end
end
for fid in resfnew
push!(queue, fid)
end
end
hull.index = snew
return hull
end
# Example: unit tetrahedron
preConvexHulls()
points = [
Vect(0.0,0.0,0.0,1),
Vect(1.0,0.0,0.0,2),
Vect(0.0,1.0,0.0,3),
Vect(0.0,0.0,1.0,4),
]
hull = quickHull3d(points, length(points))
@printf("%.3f\n", getSurfaceArea(hull))
- Output:
2.366
Lua
-- QuickHull3D in Lua
local MAXN = 2500
local EPS = 1e-8
-- Globals
local FAC = {} -- list of facets
local pts = {} -- points assigned to each facet
local TIME = 0
local e = { {}, {} } -- horizon edges: e[1][ptid], e[2][ptid]
local vistime = {} -- per‐point visit timestamp
local queue = {} -- facet queue
local resfnew, resfdel, respt = {}, {}, {} -- temp lists
-- Initialize edge table and vistime
for k = 1,2 do
e[k] = {}
for i = 1, MAXN do
e[k][i] = { netid = 0, facetid = 0 }
end
end
for i = 1, MAXN do vistime[i] = 0 end
-- Vector class
local Vect = {}
Vect.__index = Vect
function Vect.new(x,y,z,id)
return setmetatable({ x=x, y=y, z=z, id=id }, Vect)
end
function Vect:sub(o)
return Vect.new(self.x-o.x, self.y-o.y, self.z-o.z, 0)
end
function Vect:cross(o)
return Vect.new(
self.y*o.z - self.z*o.y,
self.z*o.x - self.x*o.z,
self.x*o.y - self.y*o.x,
0
)
end
function Vect:dot(o)
return self.x*o.x + self.y*o.y + self.z*o.z
end
function Vect:mag()
return math.sqrt(self.x*self.x + self.y*self.y + self.z*self.z)
end
function Vect:eq(o)
local function eq(a,b) return math.abs(a-b)<EPS end
return eq(self.x,o.x) and eq(self.y,o.y) and eq(self.z,o.z)
end
-- Line and Plane
local Line = {}
function Line.new(u,v) return { u=u, v=v } end
local Plane = {}
Plane.__index = Plane
function Plane.new(u,v,w)
return setmetatable({ u=u, v=v, w=w }, Plane)
end
function Plane:normal()
return self.v:sub(self.u):cross(self.w:sub(self.u))
end
function Plane:vecAt(i)
-- i = 1,2,3
if i==1 then return self.u
elseif i==2 then return self.v
else return self.w
end
end
function Plane:vecId(i)
return self:vecAt(i).id
end
-- Facet class
local Facet = {}
Facet.__index = Facet
function Facet.new()
return setmetatable({ n={0,0,0}, id=0, vistime=0, isdel=false, p=Plane.new(nil,nil,nil) }, Facet)
end
-- Convex hull manager
local ConvexHulls3d = {}
ConvexHulls3d.__index = ConvexHulls3d
function ConvexHulls3d.new(idx)
return setmetatable({ index=idx, surfacearea=0.0 }, ConvexHulls3d)
end
-- Utilities
local function gtr(a,b) return a-b>EPS end
local function ABS(x) return x<0 and -x or x end
-- Distances
local function distPointPlane(v,p)
local n = p:normal()
return (v:sub(p.u):dot(n)) / n:mag()
end
local function distPointLine(v,l)
local d = v:sub(l.u):mag()
if d==0 then return 0 end
return l.v:sub(l.u):cross(v:sub(l.u)):mag() / l.v:sub(l.u):mag()
end
local function distPointPoint(a,b)
return a:sub(b):mag()
end
local function isAbove(v,p)
return gtr(v:sub(p.u):dot(p:normal()), 0)
end
-- Pre‐initialize
local function preConvexHulls()
table.insert(pts,{})
table.insert(FAC,Facet.new()) -- dummy at index 1
end
-- DFS for surface area
function ConvexHulls3d:dfsArea(fidx)
local f = FAC[fidx]
if f.vistime == TIME then return end
f.vistime = TIME
local n = f.p:normal()
self.surfacearea = self.surfacearea + n:mag()/2
for i=1,3 do self:dfsArea(f.n[i]) end
end
function ConvexHulls3d:getSurfaceArea()
if self.surfacearea>0 then return self.surfacearea end
TIME = TIME + 1
self:dfsArea(self.index)
return self.surfacearea
end
-- Find the horizon from facet fidx looking toward point p
function ConvexHulls3d:getHorizon(fidx,p,resfdel)
local f = FAC[fidx]
if not isAbove(p,f.p) then return 0 end
if f.vistime==TIME then return -1 end
f.vistime = TIME
f.isdel = true
table.insert(resfdel, fidx)
local ret = -2
for i=1,3 do
local r = self:getHorizon(f.n[i], p, resfdel)
if r==0 then
local a = f.p:vecId(i)
local b = f.p:vecId((i%3)+1)
for j,pt in ipairs{a,b} do
if vistime[pt]~=TIME then
vistime[pt]=TIME
e[1][pt] = { netid = (j==1 and b or a), facetid = f.n[i] }
else
e[2][pt] = { netid = (j==1 and b or a), facetid = f.n[i] }
end
end
ret = a
elseif r~=-1 and r~=-2 then
ret = r
end
end
return ret
end
-- Build initial tetrahedron
local function getStart(point, totp)
-- pick 6 extremes
local pt = {}
for i=1,6 do pt[i]=point[1] end
for i=1,totp do
local v = point[i]
if gtr(v.x,pt[1].x) then pt[1]=v end
if gtr(pt[2].x,v.x) then pt[2]=v end
if gtr(v.y,pt[3].y) then pt[3]=v end
if gtr(pt[4].y,v.y) then pt[4]=v end
if gtr(v.z,pt[5].z) then pt[5]=v end
if gtr(pt[6].z,v.z) then pt[6]=v end
end
-- furthest pair
local s = { pt[1], pt[2], point[1], point[1] }
for i=1,6 do
for j=i+1,6 do
if gtr(distPointPoint(pt[i],pt[j]), distPointPoint(s[1],s[2])) then
s[1],s[2]=pt[i],pt[j]
end
end
end
-- furthest from line
local L = Line.new(s[1],s[2])
s[3]=pt[1]
for i=1,6 do
if gtr(distPointLine(pt[i],L), distPointLine(s[3],L)) then
s[3]=pt[i]
end
end
-- furthest from plane
s[4]=point[1]
local base = Plane.new(s[1],s[2],s[3])
for i=1,totp do
if gtr(ABS(distPointPlane(point[i],base)), ABS(distPointPlane(s[4],base))) then
s[4]=point[i]
end
end
if gtr(0, distPointPlane(s[4],base)) then
s[2],s[3]=s[3],s[2]
end
-- make 4 facets
local fidx = {}
for i=1,4 do
local f = Facet.new()
table.insert(FAC,f)
f.id = #FAC
fidx[i] = #FAC
end
FAC[fidx[1]].p = Plane.new(s[1],s[3],s[2])
FAC[fidx[2]].p = Plane.new(s[1],s[2],s[4])
FAC[fidx[3]].p = Plane.new(s[2],s[3],s[4])
FAC[fidx[4]].p = Plane.new(s[3],s[1],s[4])
-- neighbors
FAC[fidx[1]].n = { fidx[4], fidx[3], fidx[2] }
FAC[fidx[2]].n = { fidx[1], fidx[3], fidx[4] }
FAC[fidx[3]].n = { fidx[1], fidx[4], fidx[2] }
FAC[fidx[4]].n = { fidx[1], fidx[2], fidx[3] }
-- prepare pts lists
for i=1,4 do pts[fidx[i]] = {} end
-- assign points above each face
for i=1,totp do
local v = point[i]
if not (v:eq(s[1]) or v:eq(s[2]) or v:eq(s[3]) or v:eq(s[4])) then
for j=1,4 do
if isAbove(v,FAC[fidx[j]].p) then
table.insert(pts[fidx[j]], v)
break
end
end
end
end
return ConvexHulls3d.new(fidx[1])
end
-- Main QuickHull3D
local function quickHull3d(point, totp)
local hull = getStart(point, totp)
-- init queue
queue = { hull.index }
for _,nbr in ipairs(FAC[hull.index].n) do
table.insert(queue, nbr)
end
local snew = hull.index
while #queue>0 do
local nf = table.remove(queue,1)
local face = FAC[nf]
if face.isdel or #pts[nf]==0 then
if not face.isdel then snew = nf end
else
-- farthest point from face
local p = pts[nf][1]
for i=2,#pts[nf] do
local v = pts[nf][i]
if gtr(distPointPlane(v,face.p), distPointPlane(p,face.p)) then
p = v
end
end
-- find horizon
TIME = TIME + 1
resfdel = {}
local s = hull:getHorizon(nf, p, resfdel)
-- build new faces around horizon
resfnew = {}
TIME = TIME + 1
local from = -1
local lastf, fstf = -1, -1
while vistime[s] ~= TIME do
vistime[s] = TIME
local net, adj = 0,0
if e[1][s].netid == from then
net = e[2][s].netid; adj = e[2][s].facetid
else
net = e[1][s].netid; adj = e[1][s].facetid
end
-- find indices of s and net on face adj
local pt1,pt2 = 1,1
for i=1,3 do
if point[s]:eq(FAC[adj].p:vecAt(i)) then pt1=i end
if point[net]:eq(FAC[adj].p:vecAt(i)) then pt2=i end
end
if (pt1%3)+1 ~= pt2 then pt1,pt2 = pt2,pt1 end
-- create new facet
local nfnew = Facet.new()
nfnew.p = Plane.new(FAC[adj].p:vecAt(pt2),
FAC[adj].p:vecAt(pt1), p)
table.insert(FAC,nfnew)
nfnew.id = #FAC
table.insert(pts,{})
table.insert(resfnew, nfnew.id)
-- link adjacency
nfnew.n[1] = adj
FAC[adj].n[pt1] = nfnew.id
if lastf~=-1 then
local A1,A2 = nfnew.p:vecAt(2), FAC[lastf].p.u
if A1:eq(A2) then
nfnew.n[2], FAC[lastf].n[3] = lastf, nfnew.id
else
nfnew.n[3], FAC[lastf].n[2] = lastf, nfnew.id
end
else
fstf = nfnew.id
end
lastf = nfnew.id
from, s = s, net
end
-- close loop
local A1,A2 = FAC[fstf].p:vecAt(2), FAC[lastf].p.u
if A1:eq(A2) then
FAC[fstf].n[2], FAC[lastf].n[3] = lastf, fstf
else
FAC[fstf].n[3], FAC[lastf].n[2] = lastf, fstf
end
-- collect deleted points
respt = {}
for _,fid in ipairs(resfdel) do
for _,v in ipairs(pts[fid]) do
table.insert(respt, v)
end
pts[fid] = {}
end
-- reassign
for _,v in ipairs(respt) do
if v~=p then
for _,fid in ipairs(resfnew) do
if isAbove(v,FAC[fid].p) then
table.insert(pts[fid], v)
break
end
end
end
end
-- enqueue new faces
for _,fid in ipairs(resfnew) do
table.insert(queue, fid)
end
end
end
hull.index = snew
return hull
end
-- Example usage
preConvexHulls()
local n = 4
local point = {}
local coords = {
{0,0,0}, {1,0,0}, {0,1,0}, {0,0,1}
}
for i=1,n do
local c = coords[i]
point[i] = Vect.new(c[1],c[2],c[3], i)
end
local hull = quickHull3d(point, n)
print(string.format("%.3f", hull:getSurfaceArea()))
- Output:
2.366
Phix
constant MAXN = 2500,
EPS = 1e-8
// Globals
integer TIME = 0
sequence Facets = {},
pts = {}, // vect
e = repeat(repeat(0,MAXN),2), // Edge
vistimeArr = repeat(0,MAXN),
resfnew,
resfdel
// Basic floating comparisons
function eq(atom a, b) return abs(a-b) < EPS end function
function gtr(atom a, b) return a-b > EPS end function
// Vect represents a 3D point/vector
enum X, Y, Z, VECID
type vect(sequence s)
return length(s)=VECID
and atom(s[X])
and atom(s[Y])
and atom(s[Z])
and integer(s[VECID])
end type
function Sub(vect a, b)
vect res = {a[X]-b[X], a[Y]-b[Y], a[Z]-b[Z], 0}
return res
end function
function Cross(vect a, b)
vect res = {a[Y]*b[Z] - a[Z]*b[Y],
a[Z]*b[X] - a[X]*b[Z],
a[X]*b[Y] - a[Y]*b[X],
0}
return res
end function
function Dot(vect a, b)
return a[X]*b[X] + a[Y]*b[Y] + a[Z]*b[Z]
end function
function Mag(vect a)
return sqrt(a[X]*a[X] + a[Y]*a[Y] + a[Z]*a[Z])
end function
function Eq(vect a, b)
return eq(a[X], b[X])
and eq(a[Y], b[Y])
and eq(a[Z], b[Z])
end function
// line is the segment u->v
enum U, V, W
type line(sequence l)
return length(l)=V
and vect(l[U])
and vect(l[V])
end type
// Plane is defined by three points u, v, w
type plane(sequence p)
return length(p)=W
and vect(p[U])
and vect(p[V])
and vect(p[W])
end type
// Normal vector of the plane
function Normal(plane p)
vect res = Cross(Sub(p[V],p[U]),Sub(p[W],p[U]))
return res
end function
// Facet is a face in the hull
enum NEIGH, /*FACET_ID,*/ VISTIME, ISDEL, P
type Facet(sequence f)
return length(f)=P
and sequence(f[NEIGH]) -- 3 ints
-- and integer(f[FACET_ID])
and integer(f[VISTIME])
and bool(f[ISDEL])
and plane(f[P])
end type
// Edge on the horizon
enum NETID, FACET_ID
type Edge(sequence e)
return length(e)=FACET_ID
and integer(e[NETID])
and integer(e[FACET_ID])
end type
// ConvexHulls3d manages hull construction
enum HULLID, SURFACEAREA
type ConvexHulls3d(sequence h)
return length(h)=SURFACEAREA
and integer(h[HULLID])
and atom(h[SURFACEAREA])
end type
// Distances
function distPointPlane(vect v, plane p)
vect np = Normal(p)
atom num := Dot(Sub(v,p[U]),np),
den := Mag(np)
return num / den
end function
function distPointLine(vect v, line f)
atom d := Mag(Sub(v,f[U]))
if d == 0 then return 0 end if
return Mag(Cross(Sub(f[V],f[U]),Sub(v,f[U]))) / Mag(Sub(f[V],f[U]))
end function
function distPointPoint(vect u, v)
return Mag(Sub(u,v))
end function
function isAbove(vect v, plane p)
return gtr(Dot(Sub(v,p[U]),Normal(p)), 0)
end function
// DFS to accumulate area
function dfsArea(ConvexHulls3d h, integer nf)
if Facets[nf][VISTIME] != TIME then
Facets[nf][VISTIME] = TIME
vect nrm := Normal(Facets[nf][P])
h[SURFACEAREA] += Mag(nrm) / 2
for i=1 to 3 do
h = dfsArea(h,Facets[nf][NEIGH][i])
end for
end if
return h
end function
// GetSurfaceArea returns (and caches) the hull surface area
function GetSurfaceArea(ConvexHulls3d h)
if not gtr(h[SURFACEAREA], 0) then
TIME += 1
h = dfsArea(h,h[HULLID])
end if
return h[SURFACEAREA]
end function
// Recursively find the horizon around point p
function getHorizon(ConvexHulls3d h, integer f, vect p)
if not isAbove(p, Facets[f][P]) then
return 0
end if
if Facets[f][VISTIME] == TIME then
return -1
end if
Facets[f][VISTIME] = TIME
Facets[f][ISDEL] = true
resfdel = append(resfdel, f)
integer ret := -2
for i=1 to 3 do
integer r := getHorizon(h,Facets[f][NEIGH][i], p)
if r == 0 then
integer j = {2,3,1}[i],
a := Facets[f][P][i][VECID],
b := Facets[f][P][j][VECID]
for j, pt in {a, b} do
integer edgedx = 2
if vistimeArr[pt] != TIME then
vistimeArr[pt] = TIME
edgedx = 1
end if
e[edgedx][pt] = {{b,a}[j], Facets[f][NEIGH][i]}
end for
ret = a
elsif r != -1 and r != -2 then
ret = r
end if
end for
return ret
end function
// Build the initial tetrahedron
function getStart(sequence points)
sequence pt = repeat(points[1],6)
vect s1,s2,s3,s4
// pick extreme coords
for v in points do
if gtr(v[X], pt[1][X]) then pt[1] = v end if
if gtr(pt[2][X], v[X]) then pt[2] = v end if
if gtr(v[Y], pt[3][Y]) then pt[3] = v end if
if gtr(pt[4][Y], v[Y]) then pt[4] = v end if
if gtr(v[Z], pt[5][Z]) then pt[5] = v end if
if gtr(pt[6][Z], v[Z]) then pt[6] = v end if
end for
// take furthest pair
{s1,s2} = pt[1..2]
for i,pi in pt do
for pj in pt from i+1 do
if gtr(distPointPoint(pi, pj),
distPointPoint(s1, s2)) then
{s1, s2} = {pi, pj}
end if
end for
end for
// furthest from that line
line L := {s1, s2}
s3 = pt[1]
for pi in pt from 2 do
if gtr(distPointLine(pi, L),
distPointLine(s3, L)) then
s3 = pi
end if
end for
// furthest from the plane
s4 = points[1]
plane base := {s1, s2, s3}
for p in points do
if gtr(abs(distPointPlane(p, base)),
abs(distPointPlane(s4, base))) then
s4 = p
end if
end for
if gtr(0, distPointPlane(s4, base)) then
{s2, s3} = {s3, s2}
end if
// make 4 facets
Facets = {{{4, 3, 2},0,0,{s1, s3, s2}},
{{1, 3, 4},0,0,{s1, s2, s4}},
{{1, 4, 2},0,0,{s2, s3, s4}},
{{1, 2, 3},0,0,{s3, s1, s4}}}
// prepare point-lists
for i=1 to 4 do
pts = append(pts, {})
end for
// assign remaining points
for v in points do
if not Eq(v,s1)
and not Eq(v,s2)
and not Eq(v,s3)
and not Eq(v,s4) then
for j=1 to 4 do
if isAbove(v, Facets[j][P]) then
pts[j] = append(pts[j], v)
exit
end if
end for
end if
end for
ConvexHulls3d res = {1,0}
return res
end function
// The main QuickHull3D loop
function quickHull3d(sequence points)
ConvexHulls3d hull := getStart(points)
sequence queue = hull[HULLID]&Facets[hull[HULLID]][NEIGH]
integer snew := 0
while length(queue) do
integer nf := queue[1]
queue = queue[2..$]
if Facets[nf][ISDEL] or length(pts[nf]) == 0 then
if not Facets[nf][ISDEL] then
snew = nf
end if
else
// find farthest point
vect p := pts[nf][1]
for v in pts[nf] from 2 do
if gtr(distPointPlane(v, Facets[nf][P]),
distPointPlane(p, Facets[nf][P])) then
p = v
end if
end for
// find horizon
TIME += 1
resfdel = {}
integer s := getHorizon(hull,nf, p)
// build new faces around the horizon
resfnew = {}
TIME += 1
integer frm := -1,
lastf := -1,
fstf := -1
while vistimeArr[s] != TIME do
vistimeArr[s] = TIME
integer net, f, fnew
if e[0][s][NETID] == frm then
{net, f} = {e[1][s][NETID], e[1][s][FACET_ID]}
else
{net, f} = {e[0][s][NETID], e[0][s][FACET_ID]}
end if
// find indices of s and net on facet f
integer pt1, pt2
for i=1 to 3 do
if Eq(points[s],Facets[f][P][i]) then
pt1 = i
end if
if Eq(points[net],Facets[f][P][i]) then
pt2 = i
end if
end for
if {2,3,1}[pt1] != pt2 then
{pt1, pt2} = {pt2, pt1}
end if
// create new facet facing outwards
Facets = append(Facets, {{0,0,0},0,0,{Facets[f][P][pt2], Facets[f][P][pt1], p}})
fnew = length(Facets)
pts = append(pts, {})
resfnew = append(resfnew, fnew)
Facets[fnew][NEIGH][1] = f
Facets[f][NEIGH][pt1] = fnew
if lastf >= 0 then
// link with previous new facet
if Eq(Facets[fnew][P][2],Facets[lastf][P][U]) then
{Facets[fnew][NEIGH][1], Facets[lastf][NEIGH][2]} = {lastf, fnew}
else
{Facets[fnew][NEIGH][2], Facets[lastf][NEIGH][1]} = {lastf, fnew}
end if
else
fstf = fnew
end if
lastf = fnew
frm = s
s = net
end while
// close the loop
if Eq(Facets[fstf][P][V],Facets[lastf][P][U]) then
{Facets[fstf][NEIGH][1], Facets[lastf][NEIGH][2]} = {lastf, fstf}
else
{Facets[fstf][NEIGH][2], Facets[lastf][NEIGH][1]} = {lastf, fstf}
end if
// collect points from deleted faces
sequence respt = {}
for fid in resfdel do
respt &= pts[fid]
pts[fid] = {}
end for
// reassign them
for v in respt do
if v != p then
for fid in resfnew do
if isAbove(v, Facets[fid][P]) then
pts[fid] = append(pts[fid], v)
exit
end if
end for
end if
end for
// enqueue new faces
queue &= resfnew
end if
end while
hull[HULLID] = snew
return hull
end function
// Example: a unit tetrahedron
constant unit_tetrahedron = {{0, 0, 0, 1},
{1, 0, 0, 2},
{0, 1, 0, 3},
{0, 0, 1, 4}}
printf(1,"%.3f\n", GetSurfaceArea(quickHull3d(unit_tetrahedron)))
- Output:
2.366
Python
import math
# Constants
MAXN = 2500
EPS = 1e-8
class Vect:
def __init__(self, x=0.0, y=0.0, z=0.0, id=0):
self.x = x
self.y = y
self.z = z
self.id = id
def __sub__(self, other):
return Vect(self.x - other.x, self.y - other.y, self.z - other.z)
def __truediv__(self, other):
return Vect(self.y * other.z - self.z * other.y,
self.z * other.x - self.x * other.z,
self.x * other.y - self.y * other.x)
def __mul__(self, other):
return self.x * other.x + self.y * other.y + self.z * other.z
def m(self):
return math.sqrt(self.x * self.x + self.y * self.y + self.z * self.z)
def __eq__(self, other):
return eq(self.x, other.x) and eq(self.y, other.y) and eq(self.z, other.z)
def __repr__(self):
return f"Vect({self.x}, {self.y}, {self.z}, {self.id})"
class Line:
def __init__(self, u=None, v=None):
self.u = u
self.v = v
class Plane:
def __init__(self, u=None, v=None, w=None):
self.vec = [u, v, w] if u is not None and v is not None and w is not None else [None, None, None]
def normal(self):
return (self.vec[1] - self.vec[0]) / (self.vec[2] - self.vec[0])
def u(self):
return self.vec[0]
def gtr(a, b):
return a - b > EPS
def eq(a, b):
return -EPS < a - b < EPS
def Abs(x):
return -x if gtr(0, x) else x
# Signed distance
def dist_point_plane(v, p):
return (v - p.u()) * p.normal() / p.normal().m()
# Unsigned distance
def dist_point_line(v, f):
return (((f.v - f.u) / (v - f.u)).m() / (f.v - f.u).m()) if (v - f.u).m() > 0 else 0
def dist_point_point(u, v):
return (u - v).m()
def isabove(v, p):
return gtr((v - p.u()) * p.normal(), 0)
# Convex Hull Structures
TIME = 0
class Facet:
def __init__(self, id=0, p=None):
# neighbor，correspond to point (u->v, v->w, w->u)
self.n = [0, 0, 0]
self.id = id
# access timestamp
self.vistime = 0
self.isdel = False
self.p = p if p is not None else Plane()
def in_(self, n1, n2, n3):
self.n = [n1, n2, n3]
def __repr__(self):
return f"Facet(id={self.id}, isdel={self.isdel}, vistime={self.vistime})"
# Edge of the horizon
class Edge:
def __init__(self):
self.netid = 0
self.facetid = 0
# Store all faces
FAC = []
class ConvexHulls3d:
def __init__(self, index):
# Index face
self.index = index
self.surfacearea = 0.0
def dfsArea(self, nf):
global TIME, FAC
# Already visited in current timestamp
if FAC[nf].vistime == TIME:
return
FAC[nf].vistime = TIME
if FAC[nf].p.normal() is not None: # check if plane is initialized
self.surfacearea += FAC[nf].p.normal().m() / 2
for i in range(3):
self.dfsArea(FAC[nf].n[i])
def getSurfaceArea(self):
global TIME, FAC
if gtr(self.surfacearea, 0):
return self.surfacearea
TIME += 1
self.dfsArea(self.index)
return self.surfacearea
def getHorizon(self, f, p, vistime, e1, e2, resfdel):
global TIME, FAC
if not isabove(p, FAC[f].p):
return 0
if FAC[f].vistime == TIME:
return -1
FAC[f].vistime = TIME
# Mark the deleted face
FAC[f].isdel = True
resfdel.append(FAC[f].id)
ret = -2
for i in range(3):
res = self.getHorizon(FAC[f].n[i], p, vistime, e1, e2, resfdel)
if res == 0:
pt = [FAC[f].p.vec[i].id, FAC[f].p.vec[(i + 1) % 3].id]
for j in range(2):
if vistime[pt[j]] != TIME:
vistime[pt[j]] = TIME
e1[pt[j]].netid = pt[(j + 1) % 2]
e1[pt[j]].facetid = FAC[f].n[i]
else:
e2[pt[j]].netid = pt[(j + 1) % 2]
e2[pt[j]].facetid = FAC[f].n[i]
ret = pt[0]
elif res != -1 and res != -2:
# The face is enclosed in the middle
ret = res
return ret
# Global Variables
# Global points
pts = []
# Construct initial simplex
def getStart(point, totp):
global FAC, pts
pt = [point[1]] * 6
s = [point[1]] * 4
# Find the maximum point of the coordinate axis
for i in range(2, totp + 1):
if gtr(point[i].x, pt[0].x):
pt[0] = point[i]
if gtr(pt[1].x, point[i].x):
pt[1] = point[i]
if gtr(point[i].y, pt[2].y):
pt[2] = point[i]
if gtr(pt[3].y, point[i].y):
pt[3] = point[i]
if gtr(point[i].z, pt[4].z):
pt[4] = point[i]
if gtr(pt[5].z, point[i].z):
pt[5] = point[i]
# Take the two points with the largest distance
for i in range(6):
for j in range(i + 1, 6):
if gtr(dist_point_point(pt[i], pt[j]), dist_point_point(s[0], s[1])):
s[0] = pt[i]
s[1] = pt[j]
# Take the point farthest from the line connecting the two points
for i in range(6):
if gtr(dist_point_line(pt[i], Line(s[0], s[1])), dist_point_line(s[2], Line(s[0], s[1]))):
s[2] = pt[i]
# Take the point farthest from the face
for i in range(1, totp + 1): # !!
if gtr(Abs(dist_point_plane(point[i], Plane(s[0], s[1], s[2]))),
Abs(dist_point_plane(s[3], Plane(s[0], s[1], s[2])))):
s[3] = point[i]
# Ensure that the constructed face faces outwards
if gtr(0, dist_point_plane(s[3], Plane(s[0], s[1], s[2]))):
s[1], s[2] = s[2], s[1]
# Construct simplex
f = [0] * 4
for i in range(4):
FAC.append(Facet(id=len(FAC)))
f[i] = len(FAC) - 1
FAC[f[0]].p = Plane(s[0], s[2], s[1]) # Bottom face
FAC[f[1]].p = Plane(s[0], s[1], s[3])
FAC[f[2]].p = Plane(s[1], s[2], s[3])
FAC[f[3]].p = Plane(s[2], s[0], s[3])
FAC[f[0]].in_(f[3], f[2], f[1])
FAC[f[1]].in_(f[0], f[2], f[3])
FAC[f[2]].in_(f[0], f[3], f[1])
FAC[f[3]].in_(f[0], f[1], f[2])
# Assign point set space to four faces
for i in range(4):
pts.append([])
# Assign points to four faces
for i in range(1, totp + 1):
if point[i] == s[0] or point[i] == s[1] or point[i] == s[2] or point[i] == s[3]:
continue
for j in range(4):
if isabove(point[i], FAC[f[j]].p):
pts[f[j]].append(point[i])
break
# Return the initial simplex, using a face as index
return ConvexHulls3d(f[0])
# Border line graph information
e = [[Edge() for _ in range(MAXN)] for _ in range(2)]
# Timestamp of each point access
vistime = [0] * MAXN
que = []
# Save the newly constructed face
resfnew = []
# Save the deleted face
resfdel = []
# Save the point to be allocated
respt = []
def quickHull3d(point, totp):
global FAC, pts, TIME, e, vistime, que, resfnew, resfdel, respt
hull = getStart(point, totp)
# Add the face of initial simplex to queue
que = [hull.index]
for i in range(3):
que.append(FAC[hull.index].n[i])
# snew saves index face of the final convex hull
snew = 0
while que:
nf = que.pop(0)
# Skip if the current face has been deleted
if FAC[nf].isdel:
continue
# Skip if no vertices are allocated to the current face
if not pts[nf]:
snew = nf
continue
# Find the farthest point from the face
p = pts[nf][0]
for i in range(1, len(pts[nf])):
if gtr(dist_point_plane(pts[nf][i], FAC[nf].p), dist_point_plane(p, FAC[nf].p)):
p = pts[nf][i]
# Find the horizon
TIME += 1
resfdel = []
# The current face must be deleted, so start dfs from current face
s = hull.getHorizon(nf, p, vistime, e[0], e[1], resfdel)
# Iterate over horizon(go around a circle), construct new face
# When finding horizon, we can't know whether an edge is clockwise or counterclockwise, so it needs to be judged here
resfnew = []
TIME += 1
from_ = 0 # The previous visited point
lastf = 0 # The last created face
fstf = 0 # The first created face
while vistime[s] != TIME:
# Record whether the current point has been visited with timestamp
vistime[s] = TIME
net = 0 # Next point
f = 0 # The unseen face connected to the current edge on horizon
fnew = 0 # New face
# Make sure the traversal direction is correct
if e[0][s].netid == from_:
net = e[1][s].netid
f = e[1][s].facetid
else:
net = e[0][s].netid
f = e[0][s].facetid
# Find the counterclockwise information of these two points on the adjacent face
pt1 = -1
pt2 = -1
for i in range(3):
if point[s] == FAC[f].p.vec[i]:
pt1 = i
if point[net] == FAC[f].p.vec[i]:
pt2 = i
# Make sure pt1->pt2 is arranged counterclockwise by adjacent face points
if (pt1 + 1) % 3 != pt2:
pt1, pt2 = pt2, pt1 # swap
# The face constructed in this way faces outward
FAC.append(Facet(id=len(FAC), p=Plane(FAC[f].p.vec[pt2], FAC[f].p.vec[pt1], p)))
fnew = len(FAC) - 1
pts.append([])
resfnew.append(fnew)
# Maintain adjacency information
FAC[fnew].n[0] = f
FAC[f].n[pt1] = fnew
if lastf:
# Can't determine whether to traverse clockwise or counterclockwise in advance
# Maintain adjacency information between new faces
if FAC[fnew].p.vec[1] == FAC[lastf].p.vec[0]:
FAC[fnew].n[1] = lastf
FAC[lastf].n[2] = fnew
else:
FAC[fnew].n[2] = lastf
FAC[lastf].n[1] = fnew
else:
fstf = fnew # No new face yet
lastf = fnew
from_ = s
s = net
# Give the new face head and tail maintenance critical information
if FAC[fstf].p.vec[1] == FAC[lastf].p.vec[0]:
FAC[fstf].n[1] = lastf
FAC[lastf].n[2] = fstf
else:
FAC[fstf].n[2] = lastf
FAC[lastf].n[1] = fstf
# Get all the points to be assigned
respt = []
for i in range(len(resfdel)):
for j in range(len(pts[resfdel[i]])):
respt.append(pts[resfdel[i]][j])
pts[resfdel[i]] = []
# Assign points
for i in range(len(respt)):
if respt[i] == p:
continue # Skip the points used to create the new face
for j in range(len(resfnew)):
if isabove(respt[i], FAC[resfnew[j]].p):
pts[resfnew[j]].append(respt[i])
break # Make sure the points are not reassigned
# Add the new face to queue
for i in range(len(resfnew)):
que.append(resfnew[i])
hull.index = snew
return hull
def preConvexHulls():
global pts, FAC
# 0 for reservation
pts.append([])
FAC.append(Facet())
if __name__ == "__main__":
preConvexHulls()
n = 4 # number of points
point = [None] * (n + 1) # 1-indexed to match original code
my_input=[[0.0,0.0,0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0] ]
for i in range(1, n + 1):
x, y, z = my_input[i-1]
point[i] = Vect(x, y, z, id=i)
hull = quickHull3d(point, n)
print(f"{hull.getSurfaceArea():.3f}")
- Output:
2.366
Rust
use std::f64;
const MAXN: usize = 2500;
const EPS: f64 = 1e-8;
#[derive(Clone, Copy, Debug)]
struct Vect {
x: f64,
y: f64,
z: f64,
id: usize,
}
impl Vect {
fn new(x: f64, y: f64, z: f64, id: usize) -> Self {
Vect { x, y, z, id }
}
fn sub(&self, other: &Vect) -> Vect {
Vect::new(self.x - other.x, self.y - other.y, self.z - other.z, 0)
}
fn cross(&self, other: &Vect) -> Vect {
Vect::new(
self.y * other.z - self.z * other.y,
self.z * other.x - self.x * other.z,
self.x * other.y - self.y * other.x,
0,
)
}
fn dot(&self, other: &Vect) -> f64 {
self.x * other.x + self.y * other.y + self.z * other.z
}
fn m(&self) -> f64 {
f64::sqrt(self.x * self.x + self.y * self.y + self.z * self.z)
}
fn eq(&self, other: &Vect) -> bool {
eq(self.x, other.x) && eq(self.y, other.y) && eq(self.z, other.z)
}
}
struct Line {
u: Vect,
v: Vect,
}
impl Line {
fn new(u: Vect, v: Vect) -> Self {
Line { u, v }
}
}
#[derive(Clone)]
struct Plane {
vec: [Vect; 3],
}
impl Plane {
fn new(u: Vect, v: Vect, w: Vect) -> Self {
Plane { vec: [u, v, w] }
}
fn normal(&self) -> Vect {
self.vec[1].sub(&self.vec[0]).cross(&self.vec[2].sub(&self.vec[0]))
}
fn u(&self) -> Vect {
self.vec[0]
}
}
fn gtr(a: f64, b: f64) -> bool {
a - b > EPS
}
fn eq(a: f64, b: f64) -> bool {
-EPS < a - b && a - b < EPS
}
fn abs(x: f64) -> f64 {
if gtr(0.0, x) { -x } else { x }
}
// Signed distance
fn dist_point_plane(v: &Vect, p: &Plane) -> f64 {
let normal = p.normal();
v.sub(&p.u()).dot(&normal) / normal.m()
}
// Unsigned distance
fn dist_point_line(v: &Vect, f: &Line) -> f64 {
if v.sub(&f.u).m() > 0.0 {
f.v.sub(&f.u).cross(&v.sub(&f.u)).m() / f.v.sub(&f.u).m()
} else {
0.0
}
}
fn dist_point_point(u: &Vect, v: &Vect) -> f64 {
u.sub(v).m()
}
fn isabove(v: &Vect, p: &Plane) -> bool {
gtr(v.sub(&p.u()).dot(&p.normal()), 0.0)
}
// Convex Hull Structures
#[derive(Clone)]
struct Facet {
n: [usize; 3], // neighbors, correspond to point (u->v, v->w, w->u)
id: usize,
vistime: usize, // access timestamp
isdel: bool,
p: Plane,
}
impl Facet {
fn new(id: usize, p: Plane) -> Self {
Facet {
n: [0, 0, 0],
id,
vistime: 0,
isdel: false,
p,
}
}
fn set_neighbors(&mut self, n1: usize, n2: usize, n3: usize) {
self.n = [n1, n2, n3];
}
}
// Edge of the horizon
#[derive(Clone)]
struct Edge {
netid: usize,
facetid: usize,
}
impl Edge {
fn new() -> Self {
Edge { netid: 0, facetid: 0 }
}
}
struct ConvexHulls3d {
index: usize, // Index face
surfacearea: f64,
}
impl ConvexHulls3d {
fn new(index: usize) -> Self {
ConvexHulls3d {
index,
surfacearea: 0.0,
}
}
fn dfs_area(&mut self, nf: usize, fac: &mut Vec<Facet>, time: &mut usize) {
// Already visited in current timestamp
if fac[nf].vistime == *time {
return;
}
fac[nf].vistime = *time;
self.surfacearea += fac[nf].p.normal().m() / 2.0;
// Make a copy of neighbors to avoid borrow issues
let neighbors = fac[nf].n;
for i in 0..3 {
self.dfs_area(neighbors[i], fac, time);
}
}
fn get_surface_area(&mut self, fac: &mut Vec<Facet>, time: &mut usize) -> f64 {
if gtr(self.surfacearea, 0.0) {
return self.surfacearea;
}
*time += 1;
self.dfs_area(self.index, fac, time);
self.surfacearea
}
fn get_horizon(
&self,
f: usize,
p: Vect,
vistime: &mut [usize],
e1: &mut [Edge],
e2: &mut [Edge],
resfdel: &mut Vec<usize>,
fac: &mut Vec<Facet>,
time: usize,
) -> isize {
if !isabove(&p, &fac[f].p) {
return 0;
}
if fac[f].vistime == time {
return -1;
}
fac[f].vistime = time;
// Mark the deleted face
fac[f].isdel = true;
resfdel.push(fac[f].id);
let mut ret: isize = -2;
// Make a copy of neighbors to avoid borrow issues
let neighbors = fac[f].n;
for i in 0..3 {
let res = self.get_horizon(
neighbors[i],
p,
vistime,
e1,
e2,
resfdel,
fac,
time,
);
if res == 0 {
let pt = [
fac[f].p.vec[i].id,
fac[f].p.vec[(i + 1) % 3].id
];
for j in 0..2 {
if vistime[pt[j]] != time {
vistime[pt[j]] = time;
e1[pt[j]].netid = pt[(j + 1) % 2];
e1[pt[j]].facetid = neighbors[i];
} else {
e2[pt[j]].netid = pt[(j + 1) % 2];
e2[pt[j]].facetid = neighbors[i];
}
}
ret = pt[0] as isize;
} else if res != -1 && res != -2 {
// The face is enclosed in the middle
ret = res;
}
}
ret
}
}
// Construct initial simplex
fn get_start(point: &[Vect], totp: usize) -> (ConvexHulls3d, Vec<Facet>, Vec<Vec<Vect>>) {
let mut fac = Vec::new();
let mut pts = Vec::new();
// Add a dummy facet at index 0
fac.push(Facet::new(0, Plane::new(
Vect::new(0.0, 0.0, 0.0, 0),
Vect::new(0.0, 0.0, 0.0, 0),
Vect::new(0.0, 0.0, 0.0, 0)
)));
let mut pt = [point[1]; 6];
let mut s = [point[1]; 4];
// Find the maximum point of the coordinate axis
for i in 2..=totp {
if gtr(point[i].x, pt[0].x) {
pt[0] = point[i];
}
if gtr(pt[1].x, point[i].x) {
pt[1] = point[i];
}
if gtr(point[i].y, pt[2].y) {
pt[2] = point[i];
}
if gtr(pt[3].y, point[i].y) {
pt[3] = point[i];
}
if gtr(point[i].z, pt[4].z) {
pt[4] = point[i];
}
if gtr(pt[5].z, point[i].z) {
pt[5] = point[i];
}
}
// Take the two points with the largest distance
for i in 0..6 {
for j in (i+1)..6 {
if gtr(dist_point_point(&pt[i], &pt[j]), dist_point_point(&s[0], &s[1])) {
s[0] = pt[i];
s[1] = pt[j];
}
}
}
// Take the point farthest from the line connecting the two points
for i in 0..6 {
if gtr(
dist_point_line(&pt[i], &Line::new(s[0], s[1])),
dist_point_line(&s[2], &Line::new(s[0], s[1]))
) {
s[2] = pt[i];
}
}
// Take the point farthest from the face
for i in 1..=totp {
if gtr(
abs(dist_point_plane(&point[i], &Plane::new(s[0], s[1], s[2]))),
abs(dist_point_plane(&s[3], &Plane::new(s[0], s[1], s[2])))
) {
s[3] = point[i];
}
}
// Ensure that the constructed face faces outwards
if gtr(0.0, dist_point_plane(&s[3], &Plane::new(s[0], s[1], s[2]))) {
// Use the array swap method instead of std::mem::swap
s.swap(1, 2);
}
// Construct simplex
let mut f = [0; 4];
for i in 0..4 {
fac.push(Facet::new(
fac.len(),
Plane::new(
Vect::new(0.0, 0.0, 0.0, 0),
Vect::new(0.0, 0.0, 0.0, 0),
Vect::new(0.0, 0.0, 0.0, 0)
)
));
f[i] = fac.len() - 1;
}
fac[f[0]].p = Plane::new(s[0], s[2], s[1]); // Bottom face
fac[f[1]].p = Plane::new(s[0], s[1], s[3]);
fac[f[2]].p = Plane::new(s[1], s[2], s[3]);
fac[f[3]].p = Plane::new(s[2], s[0], s[3]);
fac[f[0]].set_neighbors(f[3], f[2], f[1]);
fac[f[1]].set_neighbors(f[0], f[2], f[3]);
fac[f[2]].set_neighbors(f[0], f[3], f[1]);
fac[f[3]].set_neighbors(f[0], f[1], f[2]);
// Assign point set space
for _ in 0..5 {
pts.push(Vec::new());
}
// Assign points to four faces
for i in 1..=totp {
if point[i].eq(&s[0]) || point[i].eq(&s[1]) || point[i].eq(&s[2]) || point[i].eq(&s[3]) {
continue;
}
for j in 0..4 {
if isabove(&point[i], &fac[f[j]].p) {
pts[f[j]].push(point[i]);
break;
}
}
}
// Return the initial simplex, using a face as index
(ConvexHulls3d::new(f[0]), fac, pts)
}
fn quick_hull_3d(point: &[Vect], totp: usize) -> (ConvexHulls3d, Vec<Facet>) {
let (mut hull, mut fac, mut pts) = get_start(point, totp);
// Add the face of initial simplex to queue
let mut que = Vec::new();
que.push(hull.index);
for i in 0..3 {
que.push(fac[hull.index].n[i]);
}
// snew saves index face of the final convex hull
let mut snew = 0;
let mut time = 0;
// Border line graph information
let mut e1 = vec![Edge::new(); MAXN];
let mut e2 = vec![Edge::new(); MAXN];
// Timestamp of each point access
let mut vistime = vec![0; MAXN];
while !que.is_empty() {
let nf = que.remove(0);
// Skip if the current face has been deleted
if fac[nf].isdel {
continue;
}
// Skip if no vertices are allocated to the current face
if pts[nf].is_empty() {
snew = nf;
continue;
}
// Find the farthest point from the face
let mut p = pts[nf][0];
for i in 1..pts[nf].len() {
if gtr(
dist_point_plane(&pts[nf][i], &fac[nf].p),
dist_point_plane(&p, &fac[nf].p)
) {
p = pts[nf][i];
}
}
// Find the horizon
time += 1;
let mut resfdel = Vec::new();
// The current face must be deleted, so start dfs from current face
let s = hull.get_horizon(
nf,
p,
&mut vistime,
&mut e1,
&mut e2,
&mut resfdel,
&mut fac,
time,
);
// Iterate over horizon(go around a circle), construct new face
let mut resfnew = Vec::new();
time += 1;
let mut from = 0; // The previous visited point
let mut lastf = 0; // The last created face
let mut fstf = 0; // The first created face
let mut s = s as usize;
while vistime[s] != time {
// Record whether the current point has been visited with timestamp
vistime[s] = time;
let net ; // Next point
let f ; // The unseen face connected to the current edge on horizon
let fnew; // New face
// Make sure the traversal direction is correct
if e1[s].netid == from {
net = e2[s].netid;
f = e2[s].facetid;
} else {
net = e1[s].netid;
f = e1[s].facetid;
}
// Find the counterclockwise information of these two points on the adjacent face
let mut pt1 = -1;
let mut pt2 = -1;
for i in 0..3 {
if point[s].eq(&fac[f].p.vec[i]) {
pt1 = i as isize;
}
if point[net].eq(&fac[f].p.vec[i]) {
pt2 = i as isize;
}
}
// Make sure pt1->pt2 is arranged counterclockwise by adjacent face points
if (pt1 + 1) % 3 != pt2 {
let temp = pt1;
pt1 = pt2;
pt2 = temp;
}
// The face constructed in this way faces outward
fac.push(Facet::new(
fac.len(),
Plane::new(
fac[f].p.vec[pt2 as usize],
fac[f].p.vec[pt1 as usize],
p
),
));
fnew = fac.len() - 1;
pts.push(Vec::new());
resfnew.push(fnew);
// Maintain adjacency information
fac[fnew].n[0] = f;
fac[f].n[pt1 as usize] = fnew;
if lastf != 0 {
// Can't determine whether to traverse clockwise or counterclockwise in advance
// Maintain adjacency information between new faces
if fac[fnew].p.vec[1].eq(&fac[lastf].p.vec[0]) {
fac[fnew].n[1] = lastf;
fac[lastf].n[2] = fnew;
} else {
fac[fnew].n[2] = lastf;
fac[lastf].n[1] = fnew;
}
} else {
fstf = fnew; // No new face yet
}
lastf = fnew;
from = s;
s = net;
}
// Give the new face head and tail maintenance critical information
if fac[fstf].p.vec[1].eq(&fac[lastf].p.vec[0]) {
fac[fstf].n[1] = lastf;
fac[lastf].n[2] = fstf;
} else {
fac[fstf].n[2] = lastf;
fac[lastf].n[1] = fstf;
}
// Get all the points to be assigned
let mut respt = Vec::new();
for i in 0..resfdel.len() {
for j in 0..pts[resfdel[i]].len() {
respt.push(pts[resfdel[i]][j]);
}
pts[resfdel[i]].clear();
}
// Assign points
for i in 0..respt.len() {
if respt[i].eq(&p) {
continue; // Skip the points used to create the new face
}
for j in 0..resfnew.len() {
if isabove(&respt[i], &fac[resfnew[j]].p) {
pts[resfnew[j]].push(respt[i]);
break; // Make sure the points are not reassigned
}
}
}
// Add the new face to queue
for i in 0..resfnew.len() {
que.push(resfnew[i]);
}
}
hull.index = snew;
(hull, fac)
}
fn main() {
let n = 4; // number of points
let mut point = vec![Vect::new(0.0, 0.0, 0.0, 0)]; // 0th element is a placeholder
let my_input = vec![
[0.0, 0.0, 0.0],
[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0]
];
for i in 1..=n {
let (x, y, z) = (my_input[i-1][0], my_input[i-1][1], my_input[i-1][2]);
point.push(Vect::new(x, y, z, i));
}
let (mut hull, mut fac) = quick_hull_3d(&point, n);
let mut time = 0;
println!("{:.3}", hull.get_surface_area(&mut fac, &mut time));
}
- Output:
2.366
Wren
var MAXN = 2500
var EPS = 1e-8
class Vect {
construct new(x, y, z, id) {
_x = x
_y = y
_z = z
_id = id
}
static new(x, y, z) { new(x, y, z, 0) }
x { _x }
x=(v) { _x = v }
y { _y }
y=(v) { _y = v }
z { _z }
z=(v) { _z = v }
id { _id }
id=(v) { _id = v }
-(other) { Vect.new(_x - other.x, _y - other.y, _z - other.z) }
/(other) {
return Vect.new(_y * other.z - _z * other.y,
_z * other.x - _x * other.z,
_x * other.y - _y * other.x)
}
*(other) { _x * other.x + _y * other.y + _z * other.z }
m { (this * this).sqrt }
==(other) { _x == other.x && _y == other.y && _z == other.z }
toString { "Vect(%(_x), %(_y), %(_z), %(_id))" }
}
class Line {
construct new(u, v) {
_u = u
_v = v
}
u { _u }
u=(x) { _u = x }
v { _v }
v=(x) { _v = x }
}
class Plane {
construct new(u, v, w) {
_vec = (u && v && w) ? [u, v, w] : [null, null, null]
}
static new() { new(null, null, null) }
normal { (_vec[1] - _vec[0]) / (_vec[2] - _vec[0]) }
u { _vec[0] }
v { _vec[1] }
w { _vec[2] }
}
var Gtr = Fn.new { |a, b| a - b > EPS }
var eq = Fn.new { |a, b| -EPS < a - b && a - b < EPS }
var abs = Fn.new { |x| Gtr.call(0, x) ? -x : x }
// Signed distance.
var distPointPlane = Fn.new { |v, p| (v - p.u) * p.normal / p.normal.m }
// Unsigned distance
var distPointLine = Fn.new { |v, f|
if ((v - f.u).m > 0) {
return ((f.v - f.u) / (v - f.u)).m / (f.v - f.u).m
}
return 0
}
var distPointPoint = Fn.new { |u, v| (u - v).m }
var isAbove = Fn.new { |v, p| Gtr.call((v - p.u) * p.normal, 0) }
// Convex Hull Structures
var TIME = 0
class Facet {
construct new(id, p) {
// neighbor，corresponds to point (u->v, v->w, w->u)
_n = [0, 0, 0]
_id = id
// access timestamp
_vistime = 0
_isdel = false
_p = p ? p : Plane.new()
}
static new() { new(0, null) }
n { _n }
id { _id }
id =(i) { _id = i }
p { _p }
p =(v) { _p = v }
vistime { _vistime }
vistime=(v) { _vistime = v }
isdel { _isdel }
isdel=(v) { _isdel = v }
setN(n1, n2, n3) { _n = [n1, n2, n3] }
toString { "Facet(id=%(_id), isdel=%(_isdel), vistime=%(_vistime))" }
}
// Edge of the horizon.
class Edge {
construct new() {
_netid = 0
_facetid = 0
}
netid { _netid }
netid=(v) { _netid = v }
facetid { _facetid }
facetid=(v) { _facetid = v }
}
// Store all faces.
var FAC = []
class ConvexHulls3d {
construct new(index) {
// Index face.
_index = index
_surfaceArea = 0
}
index { _index }
index=(i) { _index = i }
dfsArea(nf) {
// Already visited in current timestamp.
if (FAC[nf].vistime == TIME) return
FAC[nf].vistime = TIME
// Check if plane is intialized.
if (FAC[nf].p.normal) _surfaceArea = _surfaceArea + FAC[nf].p.normal.m / 2
for (i in 0..2) dfsArea(FAC[nf].n[i])
}
getSurfaceArea() {
if (Gtr.call(_surfaceArea, 0)) return _surfaceArea
TIME = TIME + 1
dfsArea(_index)
return _surfaceArea
}
getHorizon(f, p, vistime, e1, e2, resfdel) {
if (!isAbove.call(p, FAC[f].p)) return 0
if (FAC[f].vistime == TIME) return -1
FAC[f].vistime = TIME
// Mark the deleted face.
FAC[f].isdel = true
resfdel.add(FAC[f].id)
var ret = -2
for (i in 0..2) {
var res = getHorizon(FAC[f].n[i], p, vistime, e1, e2, resfdel)
if (res == 0) {
var pt = [FAC[f].p.vec[i].id, FAC[f].p.vec[(i + 1) % 3].id]
for (j in 0..1) {
if (vistime[pt[j]] != TIME) {
vistime[pt[j]] = TIME
e1[pt[j]].netid = pt[(j + 1) % 2]
e1[pt[j]].facetid = FAC[f].n[i]
} else {
e2[pt[j]].netid = pt[(j + 1) % 2]
e2[pt[j]].facetid = FAC[f].n[i]
}
}
ret = pt[0]
} else if (res != -1 && res != -2) {
// The face is enclosed in the middle.
ret = res
}
}
return ret
}
}
// Global points.
var pts = []
// Construct initial simplex.
var getStart = Fn.new { |point, totp|
var pt = [point[1]] * 6
var s = [point[1]] * 4
// Find the maximum point of the coordinate axis.
var i = 2
while (i < totp + 1) {
if (Gtr.call(point[i].x, pt[0].x)) pt[0] = point[i]
if (Gtr.call(pt[1].x, point[i].x)) pt[1] = point[i]
if (Gtr.call(point[i].y, pt[2].y)) pt[2] = point[i]
if (Gtr.call(pt[3].y, point[i].y)) pt[3] = point[i]
if (Gtr.call(point[i].z, pt[4].z)) pt[4] = point[i]
if (Gtr.call(pt[5].z, point[i].z)) pt[5] = point[i]
i = i + 1
}
// Take the two points with the largest distance.
for (i in 0..5) {
var j = i + 1
while (j < 6) {
if (Gtr.call(distPointPoint.call(pt[i], pt[j]), distPointPoint.call(s[0], s[1]))) {
s[0] = pt[i]
s[1] = pt[j]
}
j = j + 1
}
}
// Take the point farthest from the line connecting the two points.
for (i in 0..5) {
if (Gtr.call(distPointLine.call(pt[i], Line.new(s[0], s[1])),
distPointLine.call(s[2], Line.new(s[0], s[1])))) {
s[2] = pt[i]
}
}
// Take the point farthest from the face.
for (i in 1...totp + 1){
if (Gtr.call(abs.call(distPointPlane.call(point[i], Plane.new(s[0], s[1], s[2]))),
abs.call(distPointPlane.call(s[3], Plane.new(s[0], s[1], s[2]))))) {
s[3] = point[i]
}
}
// Ensure that the constructed face faces outwards.
if (Gtr.call(0, distPointPlane.call(s[3], Plane.new(s[0], s[1], s[2])))) {
s.swap(1, 2)
}
// Construct simplex
var f = List.filled(4, 0)
for (i in 0..3) {
FAC.add(Facet.new(FAC.count, null))
f[i] = FAC.count - 1
}
FAC[f[0]].p = Plane.new(s[0], s[2], s[1]) // Bottom face
FAC[f[1]].p = Plane.new(s[0], s[1], s[3])
FAC[f[2]].p = Plane.new(s[1], s[2], s[3])
FAC[f[3]].p = Plane.new(s[2], s[0], s[3])
FAC[f[0]].setN(f[3], f[2], f[1])
FAC[f[1]].setN(f[0], f[2], f[3])
FAC[f[2]].setN(f[0], f[3], f[1])
FAC[f[3]].setN(f[0], f[1], f[2])
// Assign point set space to four faces.
for (i in 0..3) pts.add([])
// Assign points to four faces.
for (i in 1...totp + 1) {
if (point[i] == s[0] || point[i] == s[1] || point[i] == s[2] || point[i] == s[3]) {
continue
}
for (j in 0..3) {
if (isAbove.call(point[i], FAC[f[j]].p)) {
pts[f[j]].add(point[i])
break
}
}
}
// Return the initial simplex, using a face as index.
return ConvexHulls3d.new(f[0])
}
// Border line graph information.
var e = List.filled(2, null)
for (i in 0..1) {
e[i] = List.filled(MAXN, null)
for (j in 0...MAXN) e[i][j] = Edge.new()
}
// Timestamp of each point access.
var vistime = List.filled(MAXN, 0)
var que = []
// Save the newly constructed face.
var resfnew = []
// Save the deleted face
var resfdel = []
// Save the point to be allocated
var respt = []
var quickHull3d = Fn.new { |point, totp|
var hull = getStart.call(point, totp)
// Add the face of initial simplex to queue.
var que = [hull.index]
for (i in 0..2) que.add(FAC[hull.index].n[i])
// snew saves index face of the final convex hull.
var snew = 0
while (!que.isEmpty) {
var nf = que.removeAt(0)
// Skip if the current face has been deleted.
if (FAC[nf].isdel) continue
// Skip if no vertices are allocated to the current face.
if (pts[nf].isEmpty) {
snew = nf
continue
}
// Find the farthest point from the face.
var p = pts[nf][0]
for (i in 1...pts[nf].count) {
if (Gtr.call(distPointPlane.call(pts[nf][i], FAC[nf].p), distPointPlane.call(p, FAC[nf].p))) {
p = pts[nf][i]
}
}
// Find the horizon
TIME = TIME + 1
resfdel = []
// The current face must be deleted, so start dfs from current face
var s = hull.getHorizon(nf, p, vistime, e[0], e[1], resfdel)
// Iterate over horizon(go around a circle), construct new face.
// When finding horizon, we can't know whether an edge is clockwise
// or counterclockwise, so it needs to be judged here.
resfnew = []
TIME = TIME + 1
var from = 0 // the previous visited point
var lastf = 0 // the last created face
var fstf = 0 // the first created face
while (vistime[s] != TIME) {
// Record whether the current point has been visited with timestamp.
vistime[s] = TIME
var net = 0 // next point
var f = 0 // the unseen face connected to the current edge on horizon
var fnew = 0 // new face
// Make sure the traversal direction is correct.
if (e[0][s].netid == from) {
net = e[1][s].netid
f = e[1][s].facetid
} else {
net = e[0][s].netid
f = e[0][s].facetid
}
// Find the counterclockwise information of these two points on the adjacent face.
var pt1 = -1
var pt2 = -1
for (i in 0..2) {
if (point[s] == FAC[f].p.vec[i]) pt1 = i
if (point[net] == FAC[f].p.vec[i]) pt2 = i
}
// Make sure pt1->pt2 is arranged counterclockwise by adjacent face points.
if ((pt1 + 1) % 3 != pt2) {
var t = pt1
pt1 = pt2
pt2 = t
}
// The face constructed in this way faces outward
FAC.add(Facet.new(FAC.count, Plane.new(FAC[f].p.vec[pt2], FAC[f].p.vec[pt1], p)))
fnew = FAC.count - 1
pts.add([])
resfnew.add(fnew)
// Maintain adjacency information.
FAC[fnew].n[0] = f
FAC[f].n[pt1] = fnew
if (lastf != 0) {
// Can't determine whether to traverse clockwise or counterclockwise in advance.
// Maintain adjacency information between new faces.
if (FAC[fnew].p.vec[1] == FAC[lastf].p.vec[0]) {
FAC[fnew].n[1] = lastf
FAC[lastf].n[2] = fnew
} else {
FAC[fnew].n[2] = lastf
FAC[lastf].n[1] = fnew
}
} else {
fstf = fnew // no new face yet
}
lastf = fnew
from = s
s = net
}
// Give the new face head and tail maintenance critical information.
if (FAC[fstf].p.vec[1] == FAC[lastf].p.vec[0]) {
FAC[fstf].n[1] = lastf
FAC[lastf].n[2] = fstf
} else {
FAC[fstf].n[2] = lastf
FAC[lastf].n[1] = fstf
}
// Get all the points to be assigned.
respt = []
for (i in 0...resfdel.count) {
for (j in 0...pts[resfdel[i]].count) respt.add(pts[resfdel[i]][j])
pts[resfdel[i]] = []
}
// Assign points.
for (i in 0...respt.count) {
if (respt[i] == p) continue // skip the points used to create the new face
for (j in 0...resfnew.count) {
if (isAbove.call(respt[i], FAC[resfnew[j]].p)) {
pts[resfnew[j]].add(respt[i])
break // make sure the points are not reassigned
}
}
}
// Add the new face to queue.
for (i in 0...resfnew.count) que.add(resfnew[i])
}
hull.index = snew
return hull
}
var preConvexHulls = Fn.new {
// 0 for reservation
pts.add([])
FAC.add(Facet.new())
}
preConvexHulls.call()
var n = 4 // number of points
var point = List.filled(n + 1, null) // 1-indexed to match original code
var myInput = [ [0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1] ]
for (i in 1..n) {
var x = myInput[i-1][0]
var y = myInput[i-1][1]
var z = myInput[i-1][2]
point[i] = Vect.new(x, y, z, i)
}
var hull = quickHull3d.call(point, n)
var sa = hull.getSurfaceArea()
System.print((sa * 1000).round / 1000)
- Output:
2.366