Quickhull
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.
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
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.
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
#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
import 'dart:math';
import 'dart:io';
// ---------------------------------------------------------------------
// Constants
// ---------------------------------------------------------------------
const int MAX_SIZE = 2500;
const double EPSILON = 1e-8;
// ---------------------------------------------------------------------
// Numerical utilities (exact copy of the C++ helpers)
// ---------------------------------------------------------------------
bool isGreaterThan(double a, double b) => (a - b) > EPSILON;
bool isEqual(double a, double b) => (a - b).abs() < EPSILON;
// ---------------------------------------------------------------------
// Geometry primitives
// ---------------------------------------------------------------------
class Vector {
double x, y, z;
int id; // (in C++ it is uint32_t, here we use int)
Vector([this.x = 0, this.y = 0, this.z = 0, this.id = 0]);
Vector subtract(Vector other) =>
Vector(x - other.x, y - other.y, z - other.z);
Vector vectorProduct(Vector other) => Vector(
y * other.z - z * other.y,
z * other.x - x * other.z,
x * other.y - y * other.x,
);
double scalarProduct(Vector other) =>
x * other.x + y * other.y + z * other.z;
double magnitude() => sqrt(x * x + y * y + z * z);
bool equals(Vector b) =>
isEqual(x, b.x) && isEqual(y, b.y) && isEqual(z, b.z);
}
// ---------------------------------------------------------------------
// Line and Plane
// ---------------------------------------------------------------------
class Line {
Vector u, v;
Line(this.u, this.v);
}
class Plane {
Vector u, v, w;
Plane([Vector? a, Vector? b, Vector? c])
: u = a ?? Vector(),
v = b ?? Vector(),
w = c ?? Vector();
Vector normal() => v.subtract(u).vectorProduct(w.subtract(u));
Vector vectorAtIndex(int i) {
switch (i) {
case 0:
return u;
case 1:
return v;
case 2:
return w;
default:
throw RangeError('Plane only has three vectors');
}
}
int vectorId(int i) => vectorAtIndex(i).id;
}
// ---------------------------------------------------------------------
// Facet and Edge
// ---------------------------------------------------------------------
class Facet {
List<int> N = List.filled(3, 0);
int id = 0;
int visitedTime = 0;
bool isDeleted = false;
Plane plane = Plane();
Facet([this.id = 0, Plane? p]) {
if (p != null) plane = p;
}
}
class Edge {
int netID = 0;
int faceID = 0;
}
// ---------------------------------------------------------------------
// Global containers (exactly the same layout as the C++ version)
// ---------------------------------------------------------------------
class _Globals {
static List<Facet> facets = <Facet>[];
static List<List<Vector>> hullPoints = <List<Vector>>[];
static int timeStep = 0;
// edges[2][MAX_SIZE]
static List<List<Edge>> edges = List.generate(
2, (_) => List.generate(MAX_SIZE, (_) => Edge()));
static List<int> visitTime = List.filled(MAX_SIZE, 0);
static List<int> queue = <int>[];
static List<int> resFNew = <int>[];
static List<int> resFDel = <int>[];
static List<Vector> resPt = <Vector>[];
}
// ---------------------------------------------------------------------
// Geometry helpers (distance, above‑test …)
// ---------------------------------------------------------------------
double distancePointPlane(Vector vec, Plane plane) {
final num = vec.subtract(plane.u).scalarProduct(plane.normal());
final den = plane.normal().magnitude();
return num / den;
}
double distancePointLine(Vector vec, Line line) {
final length = vec.subtract(line.u).magnitude();
if (length == 0.0) return 0.0;
final cross = line.v
.subtract(line.u)
.vectorProduct(vec.subtract(line.u))
.magnitude();
final denom = line.v.subtract(line.u).magnitude();
return cross / denom;
}
double distancePointPoint(Vector a, Vector b) => a.subtract(b).magnitude();
bool isAbove(Vector point, Plane plane) =>
isGreaterThan(point.subtract(plane.u).scalarProduct(plane.normal()), 0.0);
// ---------------------------------------------------------------------
// Convex‑Hull wrapper (the class that contains the surface‑area routine)
// ---------------------------------------------------------------------
class ConvexHulls3d {
int index; // the facet that currently holds the hull
double surfaceArea = 0.0;
ConvexHulls3d(this.index);
double getSurfaceArea() {
if (isGreaterThan(surfaceArea, 0.0)) return surfaceArea;
_Globals.timeStep++;
_dfsArea(index);
return surfaceArea;
}
// -----------------------------------------------------------------
// Recursive horizon search – returns the first point of the horizon
// (or –1 / –2 as sentinel values, exactly like the C++ code)
// -----------------------------------------------------------------
int getHorizon(int f, Vector point, List<int> resDel) {
final Facet Ff = _Globals.facets[f];
if (!isAbove(point, Ff.plane)) return 0;
if (Ff.visitedTime == _Globals.timeStep) return -1;
Ff.visitedTime = _Globals.timeStep;
Ff.isDeleted = true;
resDel.add(Ff.id);
int result = -2; // “not decided yet”
for (int i = 0; i < 3; ++i) {
final int ni = Ff.N[i];
final int horizon = getHorizon(ni, point, resDel);
if (horizon == 0) {
// Horizon edge – store it in the global “edges” table
final int a = _Globals.facets[f].plane.vectorId(i);
final int b = _Globals.facets[f].plane.vectorId((i + 1) % 3);
for (int idx = 0; idx < 2; ++idx) {
final int pt = (idx == 0) ? a : b;
final int facet = ni;
if (_Globals.visitTime[pt] != _Globals.timeStep) {
_Globals.visitTime[pt] = _Globals.timeStep;
_Globals.edges[0][pt].netID = (idx == 0) ? b : a;
_Globals.edges[0][pt].faceID = facet;
} else {
_Globals.edges[1][pt].netID = (idx == 0) ? b : a;
_Globals.edges[1][pt].faceID = facet;
}
}
result = a;
} else if (horizon != -1 && horizon != -2) {
result = horizon;
}
}
return result;
}
// -----------------------------------------------------------------
// Depth‑first search for surface area (exact copy of the C++ code)
// -----------------------------------------------------------------
void _dfsArea(int f) {
if (_Globals.facets[f].visitedTime == _Globals.timeStep) return;
_Globals.facets[f].visitedTime = _Globals.timeStep;
final Vector normal = _Globals.facets[f].plane.normal();
surfaceArea += normal.magnitude() / 2.0;
for (int i = 0; i < 3; ++i) {
_dfsArea(_Globals.facets[f].N[i]);
}
}
}
// ---------------------------------------------------------------------
// Helper that creates the four‑point seed tetrahedron
// ---------------------------------------------------------------------
ConvexHulls3d getInitialHull(List<Vector> points, int totalPoints) {
// -------------------------------------------------
// 1️⃣ Find the 6 extreme points (min / max on each axis)
// -------------------------------------------------
final List<Vector> extremes = List.filled(6, points[1]);
for (int i = 1; i <= totalPoints; ++i) {
final Vector p = points[i];
if (isGreaterThan(p.x, extremes[0].x)) extremes[0] = p;
if (isGreaterThan(extremes[1].x, p.x)) extremes[1] = p;
if (isGreaterThan(p.y, extremes[2].y)) extremes[2] = p;
if (isGreaterThan(extremes[3].y, p.y)) extremes[3] = p;
if (isGreaterThan(p.z, extremes[4].z)) extremes[4] = p;
if (isGreaterThan(extremes[5].z, p.z)) extremes[5] = p;
}
// -------------------------------------------------
// 2️⃣ Furthest pair among the six extremes → line
// -------------------------------------------------
Vector extreme0 = extremes[0];
Vector extreme1 = extremes[1];
for (int i = 0; i < 6; ++i) {
for (int j = i + 1; j < 6; ++j) {
final d = distancePointPoint(extremes[i], extremes[j]);
if (isGreaterThan(d, distancePointPoint(extreme0, extreme1))) {
extreme0 = extremes[i];
extreme1 = extremes[j];
}
}
}
// -------------------------------------------------
// 3️⃣ Furthest point from that line → third vertex
// -------------------------------------------------
final line = Line(extreme0, extreme1);
Vector extreme2 = extremes[0];
for (int i = 0; i < 6; ++i) {
if (isGreaterThan(
distancePointLine(extremes[i], line),
distancePointLine(extreme2, line))) {
extreme2 = extremes[i];
}
}
// -------------------------------------------------
// 4️⃣ Furthest point from the plane defined by the three above
// -------------------------------------------------
final basePlane = Plane(extreme0, extreme1, extreme2);
Vector extreme3 = points[1];
for (int i = 1; i <= totalPoints; ++i) {
final d1 = distancePointPlane(points[i], basePlane).abs();
final d2 = distancePointPlane(extreme3, basePlane).abs();
if (isGreaterThan(d1, d2)) extreme3 = points[i];
}
// If the point lies under the plane we have to flip two vertices
if (isGreaterThan(0.0, distancePointPlane(extreme3, basePlane))) {
final tmp = extreme1;
extreme1 = extreme2;
extreme2 = tmp;
}
// -------------------------------------------------
// 5️⃣ Create the first four facets (the seed tetrahedron)
// -------------------------------------------------
final List<int> f = List.filled(4, 0);
for (int i = 0; i < 4; ++i) {
_Globals.facets.add(Facet(_Globals.facets.length));
f[i] = _Globals.facets.length - 1;
}
_Globals.facets[f[0]].plane = Plane(extreme0, extreme2, extreme1);
_Globals.facets[f[1]].plane = Plane(extreme0, extreme1, extreme3);
_Globals.facets[f[2]].plane = Plane(extreme1, extreme2, extreme3);
_Globals.facets[f[3]].plane = Plane(extreme2, extreme0, extreme3);
_Globals.facets[f[0]].N = [f[3], f[2], f[1]];
_Globals.facets[f[1]].N = [f[0], f[2], f[3]];
_Globals.facets[f[2]].N = [f[0], f[3], f[1]];
_Globals.facets[f[3]].N = [f[0], f[1], f[2]];
// -------------------------------------------------
// 6️⃣ Allocate hull‑point containers (one per facet)
// -------------------------------------------------
for (int i = 0; i < 4; ++i) {
_Globals.hullPoints.add(<Vector>[]);
}
// -------------------------------------------------
// 7️⃣ Distribute the remaining points to the facet they lie “above”
// -------------------------------------------------
for (int i = 1; i <= totalPoints; ++i) {
final Vector p = points[i];
if (p.equals(extreme0) ||
p.equals(extreme1) ||
p.equals(extreme2) ||
p.equals(extreme3)) continue;
for (int j = 0; j < 4; ++j) {
if (isAbove(p, _Globals.facets[f[j]].plane)) {
_Globals.hullPoints[f[j]].add(p);
break;
}
}
}
return ConvexHulls3d(f[0]);
}
// ---------------------------------------------------------------------
// Main Quick‑Hull driver (exact port of the C++ `QuickHull3D` function)
// ---------------------------------------------------------------------
ConvexHulls3d quickHull3D(List<Vector> points, int totalPoints) {
ConvexHulls3d hull = getInitialHull(points, totalPoints);
// ---- initialise the processing queue --------------------------------
_Globals.queue.clear();
_Globals.queue.add(hull.index);
for (int i = 0; i < 3; ++i) {
_Globals.queue.add(_Globals.facets[hull.index].N[i]);
}
int newHorizon = 0;
while (_Globals.queue.isNotEmpty) {
final int nf = _Globals.queue.removeAt(0);
// facet already deleted or empty → skip (but remember the last non‑deleted)
if (_Globals.facets[nf].isDeleted || _Globals.hullPoints[nf].isEmpty) {
if (!_Globals.facets[nf].isDeleted) newHorizon = nf;
continue;
}
// -------------------------------------------------
// 1️⃣ Find the furthest point of the current facet
// -------------------------------------------------
Vector point = _Globals.hullPoints[nf][0];
for (final vec in _Globals.hullPoints[nf]) {
if (isGreaterThan(
distancePointPlane(vec, _Globals.facets[nf].plane),
distancePointPlane(point, _Globals.facets[nf].plane))) {
point = vec;
}
}
// -------------------------------------------------
// 2️⃣ Determine the horizon (all facets that will be removed)
// -------------------------------------------------
_Globals.timeStep++;
_Globals.resFDel.clear();
final int horizon = hull.getHorizon(nf, point, _Globals.resFDel);
// -------------------------------------------------
// 3️⃣ Build the new faces that replace the deleted ones
// -------------------------------------------------
_Globals.resFNew.clear();
_Globals.timeStep++;
int from = -1;
int lastF = -1;
int firstF = -1;
int cur = horizon;
while (_Globals.visitTime[cur] != _Globals.timeStep) {
// mark as visited in this sweep
_Globals.visitTime[cur] = _Globals.timeStep;
// ----- find the "next" edge on the horizon -----
int net, f;
if (_Globals.edges[0][cur].netID == from) {
net = _Globals.edges[1][cur].netID;
f = _Globals.edges[1][cur].faceID;
} else {
net = _Globals.edges[0][cur].netID;
f = _Globals.edges[0][cur].faceID;
}
// ----- locate the two vertices of the old facet that belong to this edge -----
int pt1 = 0, pt2 = 0;
for (int i = 0; i < 3; ++i) {
if (points[cur].equals(_Globals.facets[f].plane.vectorAtIndex(i))) pt1 = i;
if (points[net].equals(_Globals.facets[f].plane.vectorAtIndex(i))) pt2 = i;
}
// guarantee (pt1,pt2) are oriented consistently
if ((pt1 + 1) % 3 != pt2) {
final tmp = pt1;
pt1 = pt2;
pt2 = tmp;
}
// ----- create the new facet -----
final Plane newPlane = Plane(
_Globals.facets[f].plane.vectorAtIndex(pt2),
_Globals.facets[f].plane.vectorAtIndex(pt1),
point);
_Globals.facets.add(Facet(_Globals.facets.length, newPlane));
final int newF = _Globals.facets.length - 1;
_Globals.hullPoints.add(<Vector>[]);
_Globals.resFNew.add(newF);
// ----- link the new facet with the old neighbour -----
_Globals.facets[newF].N[0] = f;
_Globals.facets[f].N[pt1] = newF;
// ----- link consecutive new facets together (forming a fan) -----
if (lastF >= 0) {
// the orientation of the shared edge decides which N‑slot to fill
if (_Globals.facets[newF].plane.vectorAtIndex(1).equals(_Globals.facets[lastF].plane.u)) {
_Globals.facets[newF].N[1] = lastF;
_Globals.facets[lastF].N[2] = newF;
} else {
_Globals.facets[newF].N[2] = lastF;
_Globals.facets[lastF].N[1] = newF;
}
} else {
firstF = newF;
}
lastF = newF;
from = cur;
cur = net;
}
// ---- close the fan (first ↔ last) ----
if (_Globals.facets[firstF].plane.vectorAtIndex(1).equals(_Globals.facets[lastF].plane.u)) {
_Globals.facets[firstF].N[1] = lastF;
_Globals.facets[lastF].N[2] = firstF;
} else {
_Globals.facets[firstF].N[2] = lastF;
_Globals.facets[lastF].N[1] = firstF;
}
// -------------------------------------------------
// 4️⃣ Re‑assign points that were stored in the deleted facets
// -------------------------------------------------
_Globals.resPt.clear();
for (final int fId in _Globals.resFDel) {
_Globals.resPt.addAll(_Globals.hullPoints[fId]);
_Globals.hullPoints[fId].clear();
}
for (final Vector vec in _Globals.resPt) {
if (vec.equals(point)) continue;
for (final int fId in _Globals.resFNew) {
if (isAbove(vec, _Globals.facets[fId].plane)) {
_Globals.hullPoints[fId].add(vec);
break;
}
}
}
// -------------------------------------------------
// 5️⃣ Queue the new facets for further processing
// -------------------------------------------------
for (final int fId in _Globals.resFNew) {
_Globals.queue.add(fId);
}
}
// The last facet that still holds a point becomes the new “starting” facet.
hull.index = newHorizon;
return hull;
}
// ---------------------------------------------------------------------
// Helper that prepares the static containers (exact copy of C++ `prepareConvexHulls`)
// ---------------------------------------------------------------------
void prepareConvexHulls() {
// reserve index 0 (as the C++ code does)
_Globals.hullPoints.add(<Vector>[]);
_Globals.facets.add(Facet());
// initialise the edge matrix (already allocated, just reset the fields)
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < MAX_SIZE; ++j) {
_Globals.edges[i][j] = Edge();
}
}
}
// ---------------------------------------------------------------------
// Example driver (exactly the same points as the C++ main)
// ---------------------------------------------------------------------
void main() {
prepareConvexHulls();
const int n = 4;
// points[0] is a dummy placeholder – the original algorithm expects 1‑based indexing
final List<Vector> points = List.filled(n + 1, Vector());
points[1] = Vector(0, 0, 0, 1);
points[2] = Vector(1, 0, 0, 2);
points[3] = Vector(0, 1, 0, 3);
points[4] = Vector(0, 0, 1, 4);
final ConvexHulls3d hull = quickHull3D(points, n);
// three decimal places, same as the C++ `cout << fixed << setprecision(3)`
print(hull.getSurfaceArea().toStringAsFixed(3));
}
- Output:
2.366
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
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
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
// 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
# 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
/* QuickHull3D – Kotlin version (fixed + wrapper)
-------------------------------------------------
No algorithmic changes – only the two fixes explained above.
*/
import kotlin.math.abs
import kotlin.math.sqrt
private const val MAXN = 2500
private const val EPS = 1e-8
object QuickHull3D {
// -----------------------------------------------------------------
// Global structures (static fields in Java)
// -----------------------------------------------------------------
val FAC = mutableListOf<Facet>() // facets (index 0 is dummy)
val pts = mutableListOf<MutableList<Vect>>() // points belonging to each facet
var TIME = 0 // global timestamp
val e = Array(2) { Array(MAXN) { Edge() } } // adjacency edges (two per vertex)
val vistime = IntArray(MAXN) // used while walking the horizon
val queue = mutableListOf<Int>() // processing queue
val resfnew = mutableListOf<Int>() // new facets (temp)
val resfdel = mutableListOf<Int>() // deleted facets (temp)
val respt = mutableListOf<Vect>() // points to be reassigned (temp)
// -----------------------------------------------------------------
// Entry point – same test as the Java version
// -----------------------------------------------------------------
@JvmStatic
fun main(args: Array<String>) {
preConvexHulls()
// Example: unit tetrahedron (1‑based indexing, point[0] is a dummy)
val n = 4
val point = arrayOfNulls<Vect>(n + 1)
point[0] = Vect(0.0, 0.0, 0.0, 0) // dummy, never used by the algorithm
point[1] = Vect(0.0, 0.0, 0.0, 1)
point[2] = Vect(1.0, 0.0, 0.0, 2)
point[3] = Vect(0.0, 1.0, 0.0, 3)
point[4] = Vect(0.0, 0.0, 1.0, 4)
// “requireNoNulls” is now safe because point[0] is filled
val hull = quickHull3D(point.requireNoNulls(), n)
System.out.printf("%.3f%n", hull.surfaceArea())
}
// -----------------------------------------------------------------
private fun preConvexHulls() {
pts.add(mutableListOf()) // pts[0] – dummy list
FAC.add(Facet()) // FAC[0] – dummy facet
for (i in 0..1) {
for (j in 0 until MAXN) {
e[i][j] = Edge()
}
}
}
// -----------------------------------------------------------------
// Geometry primitives
// -----------------------------------------------------------------
data class Vect(
var X: Double,
var Y: Double,
var Z: Double,
var Id: Int
) {
fun sub(b: Vect) = Vect(X - b.X, Y - b.Y, Z - b.Z, 0)
fun cross(b: Vect) = Vect(
Y * b.Z - Z * b.Y,
Z * b.X - X * b.Z,
X * b.Y - Y * b.X,
0
)
fun dot(b: Vect) = X * b.X + Y * b.Y + Z * b.Z
fun mag() = sqrt(X * X + Y * Y + Z * Z)
fun eq(b: Vect) = eq(X, b.X) && eq(Y, b.Y) && eq(Z, b.Z)
}
data class Line(val U: Vect, val V: Vect)
class Plane() {
lateinit var U: Vect
lateinit var V: Vect
lateinit var W: Vect
constructor(u: Vect, v: Vect, w: Vect) : this() {
U = u; V = v; W = w
}
fun normal(): Vect = V.sub(U).cross(W.sub(U))
fun vecAt(i: Int): Vect = when (i) {
0 -> U
1 -> V
else -> W
}
fun vecId(i: Int) = vecAt(i).Id
}
class Facet {
var N = IntArray(3) // neighbour facet indices
var Id = 0
var visitTime = 0
var isDeleted = false
var P = Plane()
}
class Edge {
var netId = 0
var facetId = 0
}
// -----------------------------------------------------------------
// Convex hull container (same as Java inner class)
// -----------------------------------------------------------------
class ConvexHulls3d(var index: Int) {
var surfaceArea = 0.0
private fun dfsArea(f: Int) {
val F = FAC[f]
if (F.visitTime == TIME) return
F.visitTime = TIME
val nrm = F.P.normal()
surfaceArea += nrm.mag() / 2.0
for (i in 0..2) dfsArea(F.N[i])
}
fun surfaceArea(): Double {
if (gtr(surfaceArea, 0.0)) return surfaceArea
TIME++
dfsArea(index)
return surfaceArea
}
/** Horizon detection – returns one vertex of the horizon (or sentinel) */
fun getHorizon(f: Int, p: Vect, resDel: MutableList<Int>): Int {
val Ff = FAC[f]
if (!isAbove(p, Ff.P)) return 0
if (Ff.visitTime == TIME) return -1
Ff.visitTime = TIME
Ff.isDeleted = true
resDel.add(Ff.Id)
var ret = -2
for (i in 0..2) {
val ni = Ff.N[i]
val r = getHorizon(ni, p, resDel)
if (r == 0) {
val a = Ff.P.vecId(i)
val b = Ff.P.vecId((i + 1) % 3)
for (idx in 0..1) {
val pt = if (idx == 0) a else b
val facet = ni
if (vistime[pt] != TIME) {
vistime[pt] = TIME
e[0][pt].netId = if (idx == 0) b else a
e[0][pt].facetId = facet
} else {
e[1][pt].netId = if (idx == 0) b else a
e[1][pt].facetId = facet
}
}
ret = a
} else if (r != -1 && r != -2) {
ret = r
}
}
return ret
}
}
// -----------------------------------------------------------------
// Utility functions (same as Java static helpers)
// -----------------------------------------------------------------
private fun eq(a: Double, b: Double) = kotlin.math.abs(a - b) < EPS
private fun gtr(a: Double, b: Double) = a - b > EPS
private fun abs(x: Double) = kotlin.math.abs(x)
private fun distPointPlane(v: Vect, p: Plane): Double {
val num = v.sub(p.U).dot(p.normal())
val den = p.normal().mag()
return num / den
}
private fun distPointLine(v: Vect, l: Line): Double {
val d = v.sub(l.U).mag()
if (d == 0.0) return 0.0
return l.V.sub(l.U).cross(v.sub(l.U)).mag() / l.V.sub(l.U).mag()
}
private fun distPointPoint(a: Vect, b: Vect) = a.sub(b).mag()
private fun isAbove(v: Vect, p: Plane) = gtr(v.sub(p.U).dot(p.normal()), 0.0)
// -----------------------------------------------------------------
// Core algorithm (unchanged apart from the two fixes)
// -----------------------------------------------------------------
private fun getStart(point: Array<Vect>, totp: Int): ConvexHulls3d {
// ---- extremes in each coordinate direction ----
val extremes = Array(6) { point[0] }
for (i in 1..totp) {
val 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 among extremes ------------
var s0 = extremes[0]
var s1 = extremes[1]
for (i in 0 until 6) {
for (j in i + 1 until 6) {
val d = distPointPoint(extremes[i], extremes[j])
if (gtr(d, distPointPoint(s0, s1))) {
s0 = extremes[i]; s1 = extremes[j]
}
}
}
// ---- furthest point from that line ------------
val line = Line(s0, s1)
var s2 = extremes[0]
for (i in 0 until 6) {
if (gtr(distPointLine(extremes[i], line), distPointLine(s2, line))) {
s2 = extremes[i]
}
}
// ---- furthest point from the plane s0,s1,s2 ---
var s3 = point[1]
val basePlane = Plane(s0, s1, s2)
for (i in 1..totp) {
val d1 = abs(distPointPlane(point[i], basePlane))
val d2 = abs(distPointPlane(s3, basePlane))
if (gtr(d1, d2)) s3 = point[i]
}
// ---- fix orientation if needed ---------------
if (gtr(0.0, distPointPlane(s3, basePlane))) {
val tmp = s1; s1 = s2; s2 = tmp
}
// ---- create the 4 initial facets ---------------
val f = IntArray(4)
repeat(4) {
val F = Facet()
F.Id = FAC.size
FAC.add(F)
f[it] = FAC.size - 1
}
FAC[f[0]].P = Plane(s0, s2, s1)
FAC[f[1]].P = Plane(s0, s1, s3)
FAC[f[2]].P = Plane(s1, s2, s3)
FAC[f[3]].P = Plane(s2, s0, s3)
FAC[f[0]].N = intArrayOf(f[3], f[2], f[1])
FAC[f[1]].N = intArrayOf(f[0], f[2], f[3])
FAC[f[2]].N = intArrayOf(f[0], f[3], f[1])
FAC[f[3]].N = intArrayOf(f[0], f[1], f[2])
// ---- allocate point buckets for these facets ---
repeat(4) { pts.add(mutableListOf()) }
// ---- assign each non‑extreme point to a visible facet --
for (i in 1..totp) {
val v = point[i]
if (v.eq(s0) || v.eq(s1) || v.eq(s2) || v.eq(s3)) continue
for (j in 0..3) {
if (isAbove(v, FAC[f[j]].P)) {
pts[f[j]].add(v)
break
}
}
}
return ConvexHulls3d(f[0])
}
private fun quickHull3D(point: Array<Vect>, totp: Int): ConvexHulls3d {
val hull = getStart(point, totp)
// initialise queue with the starting facet and its three neighbours
queue.clear()
queue.add(hull.index)
repeat(3) { queue.add(FAC[hull.index].N[it]) }
var snew = 0
while (queue.isNotEmpty()) {
val nf = queue.removeAt(0)
val Fnf = FAC[nf]
// skip deleted facets or empty point sets
if (Fnf.isDeleted || pts[nf].isEmpty()) {
if (!Fnf.isDeleted) snew = nf
continue
}
// ---- farthest point of the current bucket ----
var p = pts[nf][0]
for (v in pts[nf]) {
if (gtr(distPointPlane(v, Fnf.P), distPointPlane(p, Fnf.P))) {
p = v
}
}
// ---- horizon detection -----------------------
TIME++
resfdel.clear()
val s = hull.getHorizon(nf, p, resfdel)
// ---- build new faces attached to the horizon ----
resfnew.clear()
TIME++
var from = -1
var lastf = -1
var fstf = -1
var cur = s
while (vistime[cur] != TIME) {
vistime[cur] = TIME
val (net, fIdx) = if (e[0][cur].netId == from) {
e[1][cur].netId to e[1][cur].facetId
} else {
e[0][cur].netId to e[0][cur].facetId
}
// locate the two vertices of facet fIdx that correspond to cur / net
val Ff = FAC[fIdx]
var pt1 = -1
var pt2 = -1
for (i in 0..2) {
if (point[cur].eq(Ff.P.vecAt(i))) pt1 = i
if (point[net].eq(Ff.P.vecAt(i))) pt2 = i
}
if ((pt1 + 1) % 3 != pt2) {
val tmp = pt1; pt1 = pt2; pt2 = tmp
}
// ---- create the new facet --------------------
val newF = Facet()
newF.Id = FAC.size
newF.P = Plane(Ff.P.vecAt(pt2), Ff.P.vecAt(pt1), p)
FAC.add(newF)
pts.add(mutableListOf())
val fnew = FAC.size - 1
resfnew.add(fnew)
// ---- neighbour links -------------------------
newF.N[0] = fIdx
Ff.N[pt1] = fnew
if (lastf >= 0) {
val Pf1 = newF.P
val Plast = FAC[lastf].P
if (Pf1.vecAt(1).eq(Plast.U)) {
newF.N[1] = lastf
FAC[lastf].N[2] = fnew
} else {
newF.N[2] = lastf
FAC[lastf].N[1] = fnew
}
} else {
fstf = fnew
}
// advance on the horizon loop
lastf = fnew
from = cur
cur = net
}
// ---- close the horizon loop (first ↔ last) ----
val Ffst = FAC[fstf]
val Flast = FAC[lastf]
// *** FIXED LINE *** – correct access to the vertex of the plane
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 points that disappeared with deleted facets ----
respt.clear()
for (fid in resfdel) {
respt.addAll(pts[fid])
pts[fid].clear()
}
// ---- re‑assign those points to the newly created facets ----
for (v in respt) {
if (v === p) continue
for (fid in resfnew) {
if (isAbove(v, FAC[fid].P)) {
pts[fid].add(v)
break
}
}
}
// ---- enqueue the new facets for later processing ------------
for (fid in resfnew) queue.add(fid)
}
hull.index = snew
return hull
}
}
/* -----------------------------------------------------------------
Tiny top‑level wrapper so that the JVM can find “main” in MainKt.
You can delete this function if you prefer to launch the program
with “java QuickHull3D” instead of “java MainKt”.
----------------------------------------------------------------- */
fun main(args: Array<String>) = QuickHull3D.main(args)
- Output:
2.366
-- 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
#!/usr/bin/env perl
use strict;
use warnings;
use List::Util qw(first);
# Constants
use constant MAX_SIZE => 2500;
use constant EPSILON => 0.00000001;
# Numerical utilities
sub is_greater_than {
my ($a, $b) = @_;
return ($a - $b) > EPSILON;
}
sub is_equal {
my ($a, $b) = @_;
return abs($a - $b) < EPSILON;
}
# Vector class
package Vector {
sub new {
my ($class, $x, $y, $z, $id) = @_;
$x //= 0;
$y //= 0;
$z //= 0;
$id //= 0;
return bless {
x => $x,
y => $y,
z => $z,
id => $id
}, $class;
}
sub subtract {
my ($self, $other) = @_;
return Vector->new(
$self->{x} - $other->{x},
$self->{y} - $other->{y},
$self->{z} - $other->{z}
);
}
sub vector_product {
my ($self, $other) = @_;
return Vector->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}
);
}
sub scalar_product {
my ($self, $other) = @_;
return $self->{x} * $other->{x} + $self->{y} * $other->{y} + $self->{z} * $other->{z};
}
sub magnitude {
my ($self) = @_;
return sqrt($self->{x} * $self->{x} + $self->{y} * $self->{y} + $self->{z} * $self->{z});
}
sub equals {
my ($self, $other) = @_;
return main::is_equal($self->{x}, $other->{x}) &&
main::is_equal($self->{y}, $other->{y}) &&
main::is_equal($self->{z}, $other->{z});
}
}
# Line class
package Line {
sub new {
my ($class, $u, $v) = @_;
return bless {
u => $u,
v => $v
}, $class;
}
}
# Plane class
package Plane {
sub new {
my ($class, $u, $v, $w) = @_;
return bless {
u => $u,
v => $v,
w => $w
}, $class;
}
sub normal {
my ($self) = @_;
return $self->{v}->subtract($self->{u})->vector_product($self->{w}->subtract($self->{u}));
}
sub vector_at_index {
my ($self, $i) = @_;
return $self->{u} if $i == 0;
return $self->{v} if $i == 1;
return $self->{w} if $i == 2;
return 0;
}
sub vector_id {
my ($self, $i) = @_;
return $self->vector_at_index($i)->{id};
}
}
# Facet class
package Facet {
sub new {
my ($class, $id, $plane) = @_;
return bless {
id => $id // 0,
plane => $plane,
N => [],
visited_time => 0,
is_deleted => 0
}, $class;
}
}
# Edge class
package Edge {
sub new {
my ($class) = @_;
return bless {
netID => 0,
faceID => 0
}, $class;
}
}
# ConvexHulls3d class
package ConvexHulls3d {
sub new {
my ($class, $index) = @_;
return bless {
index => $index,
surface_area => 0.0
}, $class;
}
sub get_surface_area {
my ($self) = @_;
if (main::is_greater_than($self->{surface_area}, 0.0)) {
return $self->{surface_area};
}
$main::time_step++;
$self->dfs_area($self->{index});
return $self->{surface_area};
}
sub dfs_area {
my ($self, $f) = @_;
if (($main::facets[$f]->{visited_time} // 0) == $main::time_step) {
return;
}
$main::facets[$f]->{visited_time} = $main::time_step;
my $normal = $main::facets[$f]->{plane}->normal();
$self->{surface_area} += $normal->magnitude() / 2.0;
for my $i (0..2) {
$self->dfs_area($main::facets[$f]->{N}->[$i]);
}
}
sub get_horizon {
my ($self, $f, $point, $resDel) = @_;
my $Ff = $main::facets[$f];
if (!main::is_above($point, $Ff->{plane})) {
return 0;
}
if (($Ff->{visited_time} // 0) == $main::time_step) {
return -1;
}
$Ff->{visited_time} = $main::time_step;
$Ff->{is_deleted} = 1;
push @$resDel, $Ff->{id};
my $result = -2;
for my $i (0..2) {
my $ni = $Ff->{N}->[$i];
my $horizon = $self->get_horizon($ni, $point, $resDel);
if ($horizon == 0) {
my $a = $main::facets[$f]->{plane}->vector_id($i);
my $b = $main::facets[$f]->{plane}->vector_id(($i + 1) % 3);
for my $idx (0..1) {
my $pt = ($idx == 0) ? $a : $b;
my $facet = $ni;
if (($main::visit_time[$pt] // 0) != $main::time_step) {
$main::visit_time[$pt] = $main::time_step;
$main::edges->[0]->[$pt]->{netID} = ($idx == 0) ? $b : $a;
$main::edges->[0]->[$pt]->{faceID} = $facet;
} else {
$main::edges->[1]->[$pt]->{netID} = ($idx == 0) ? $b : $a;
$main::edges->[1]->[$pt]->{faceID} = $facet;
}
}
$result = $a;
} elsif ($horizon != -1 && $horizon != -2) {
$result = $horizon;
}
}
return $result;
}
}
# Global variables
our @facets;
our @hull_points;
our $time_step = 0;
our $edges;
our @visit_time;
our @queue;
our @resfnew;
our @resfdel;
our @respt;
# Geometric utilities
sub distance_point_plane {
my ($vec, $plane) = @_;
my $num = $vec->subtract($plane->{u})->scalar_product($plane->normal());
my $den = $plane->normal()->magnitude();
return $num / $den;
}
sub distance_point_line {
my ($vec, $line) = @_;
my $length = $vec->subtract($line->{u})->magnitude();
return 0.0 if $length == 0.0;
return $line->{v}->subtract($line->{u})->vector_product($vec->subtract($line->{u}))->magnitude() /
$line->{v}->subtract($line->{u})->magnitude();
}
sub distance_point_point {
my ($a, $b) = @_;
return $a->subtract($b)->magnitude();
}
sub is_above {
my ($point, $plane) = @_;
return is_greater_than($point->subtract($plane->{u})->scalar_product($plane->normal()), 0.0);
}
sub prepare_convex_hulls {
# Reserve index 0
push @hull_points, [];
push @facets, Facet->new();
# Initialize edge vector
$edges = [];
for my $i (0..1) {
$edges->[$i] = [];
for my $j (0..MAX_SIZE-1) {
$edges->[$i]->[$j] = Edge->new();
}
}
# Initialize visit_time array
@visit_time = (0) x MAX_SIZE;
}
sub get_initial_hull {
my ($points, $total_points) = @_;
my @extremes = ($points->[1]) x 6;
for my $i (1..$total_points) {
my $point = $points->[$i];
$extremes[0] = $point if is_greater_than($point->{x}, $extremes[0]->{x});
$extremes[1] = $point if is_greater_than($extremes[1]->{x}, $point->{x});
$extremes[2] = $point if is_greater_than($point->{y}, $extremes[2]->{y});
$extremes[3] = $point if is_greater_than($extremes[3]->{y}, $point->{y});
$extremes[4] = $point if is_greater_than($point->{z}, $extremes[4]->{z});
$extremes[5] = $point if is_greater_than($extremes[5]->{z}, $point->{z});
}
# Furthest pair
my $extreme0 = $extremes[0];
my $extreme1 = $extremes[1];
for my $i (0..5) {
for my $j ($i+1..5) {
my $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
my $line = Line->new($extreme0, $extreme1);
my $extreme2 = $extremes[0];
for my $i (0..5) {
if (is_greater_than(distance_point_line($extremes[$i], $line), distance_point_line($extreme2, $line))) {
$extreme2 = $extremes[$i];
}
}
# Furthest from plane
my $extreme3 = $points->[1];
my $basePlane = Plane->new($extreme0, $extreme1, $extreme2);
for my $i (1..$total_points) {
my $distance1 = abs(distance_point_plane($points->[$i], $basePlane));
my $distance2 = abs(distance_point_plane($extreme3, $basePlane));
if (is_greater_than($distance1, $distance2)) {
$extreme3 = $points->[$i];
}
}
if (is_greater_than(0, distance_point_plane($extreme3, $basePlane))) {
($extreme1, $extreme2) = ($extreme2, $extreme1);
}
# Create 4 initial facets
my @f;
for my $i (0..3) {
push @facets, Facet->new(scalar @facets);
$f[$i] = $#facets;
}
$facets[$f[0]]->{plane} = Plane->new($extreme0, $extreme2, $extreme1);
$facets[$f[1]]->{plane} = Plane->new($extreme0, $extreme1, $extreme3);
$facets[$f[2]]->{plane} = Plane->new($extreme1, $extreme2, $extreme3);
$facets[$f[3]]->{plane} = Plane->new($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 array
for my $i (0..3) {
push @hull_points, [];
}
# Assign points
for my $i (1..$total_points) {
my $point = $points->[$i];
next if ($point->equals($extreme0) || $point->equals($extreme1) ||
$point->equals($extreme2) || $point->equals($extreme3));
for my $j (0..3) {
if (is_above($point, $facets[$f[$j]]->{plane})) {
push @{$hull_points[$f[$j]]}, $point;
last;
}
}
}
return ConvexHulls3d->new($f[0]);
}
sub QuickHull3D {
my ($points, $total_points) = @_;
my $hull = get_initial_hull($points, $total_points);
# Initialize queue
@queue = ();
push @queue, $hull->{index};
for my $i (0..2) {
push @queue, $facets[$hull->{index}]->{N}->[$i];
}
my $new_horizon = 0;
while (@queue) {
my $nf = shift @queue;
if ($facets[$nf]->{is_deleted} || !@{$hull_points[$nf]}) {
if (!$facets[$nf]->{is_deleted}) {
$new_horizon = $nf;
}
next;
}
# Farthest point
my $point = $hull_points[$nf]->[0];
for my $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 = ();
my $horizon = $hull->get_horizon($nf, $point, \@resfdel);
# Build new faces
@resfnew = ();
$time_step++;
my $from = -1;
my $last_f = -1;
my $first_f = -1;
while (($visit_time[$horizon] // 0) != $time_step) {
$visit_time[$horizon] = $time_step;
my ($net, $f, $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
my ($pt1, $pt2) = (0, 0);
for my $i (0..2) {
$pt1 = $i if $points->[$horizon]->equals($facets[$f]->{plane}->vector_at_index($i));
$pt2 = $i if $points->[$net]->equals($facets[$f]->{plane}->vector_at_index($i));
}
if (($pt1 + 1) % 3 != $pt2) {
($pt1, $pt2) = ($pt2, $pt1);
}
# Create new facet
push @facets, Facet->new(
scalar @facets,
Plane->new(
$facets[$f]->{plane}->vector_at_index($pt2),
$facets[$f]->{plane}->vector_at_index($pt1),
$point
)
);
$new_f = $#facets;
push @hull_points, [];
push @resfnew, $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 = ();
for my $f_id (@resfdel) {
push @respt, @{$hull_points[$f_id]};
$hull_points[$f_id] = [];
}
# Reassign
for my $vec (@respt) {
next if $vec->equals($point);
for my $f_id (@resfnew) {
if (is_above($vec, $facets[$f_id]->{plane})) {
push @{$hull_points[$f_id]}, $vec;
last;
}
}
}
# Enqueue new faces
push @queue, @resfnew;
}
$hull->{index} = $new_horizon;
return $hull;
}
# Main program
prepare_convex_hulls();
# Example: a tetrahedron
my $n = 4;
my @points = (undef); # Index 0 is unused
push @points, Vector->new(0, 0, 0, 1);
push @points, Vector->new(1, 0, 0, 2);
push @points, Vector->new(0, 1, 0, 3);
push @points, Vector->new(0, 0, 1, 4);
my $hull = QuickHull3D(\@points, $n);
printf "%.3f\n", $hull->get_surface_area();
- Output:
2.366
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
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
# 20250710 Raku programming solution
constant MAX_SIZE = 2500;
constant EPSILON = 1e-8;
sub is-equal($a, $b, :$eps = EPSILON) { abs($a - $b) < $eps }
sub is-greater-than($a, $b, :$eps = EPSILON) { ($a - $b) > $eps }
class Vector {
has Real ($.x, $.y, $.z) is rw;
has Int $.id is rw;
method subtract(Vector $other) {
Vector.new( :x($!x - $other.x), :y($!y - $other.y), :z($!z - $other.z) )
}
method cross-product(Vector $other) {
Vector.new(
x => $!y * $other.z - $!z * $other.y,
y => $!z * $other.x - $!x * $other.z,
z => $!x * $other.y - $!y * $other.x
)
}
method dot-product(Vector $other) {
$!x, $!y, $!z Z*+ $other.x, $other.y, $other.z # Seq
}
method magnitude { sqrt($!x² + $!y² + $!z²) }
method equals(Vector $other) {
all (($!x, $other.x),($!y,$other.y),($!z,$other.z)).map: { is-equal |$_ }
}
method gist { "Vector({$!x}, {$!y}, {$!z}, id={$!id})" }
}
class Line { has Vector ($.u, $.v) is rw }
class Plane {
has Vector ($.u, $.v, $.w) is rw;
method normal { $!v.subtract($!u).cross-product($!w.subtract($!u)) }
method vector-at(Int $i) { ($!u, $!v, $!w)[$i] // Vector.new }
method vector-id(Int $i) { self.vector-at($i).id }
}
class Facet {
has Int $.id is rw = 0;
has Int $.visited-time is rw = 0;
has Bool $.is-deleted is rw = False;
has Plane $.plane is rw;
has @.neighbors is rw;
}
class Edge { has Int ($.net-id, $.face-id) is rw }
# Global state variables
my ( @facets, @hull-points, @edges, @visit-time, @queue );
my ( @res-facets-new, @res-facets-deleted, @res-points, $time-step );
# Geometric utilities
sub distance-point-plane(Vector $vec, Plane $plane) {
$vec.subtract($plane.u).dot-product($plane.normal) / $plane.normal.magnitude
}
sub distance-point-line(Vector $vec, Line $line) {
return 0 if my $length = $vec.subtract($line.u).magnitude == 0;
$line.v.subtract($line.u).cross-product($vec.subtract($line.u)).magnitude /
$line.v.subtract($line.u).magnitude
}
sub distance-point-point(Vector $a, Vector $b) { $a.subtract($b).magnitude }
sub is-above(Vector $point, Plane $plane) {
is-greater-than($point.subtract($plane.u).dot-product($plane.normal), 0)
}
class ConvexHull3D {
has Int $.index is rw;
has Real $.surface-area is rw = 0;
method get-surface-area {
return $!surface-area if is-greater-than($!surface-area, 0);
$time-step++;
self!dfs-area($!index);
return $!surface-area
}
method !dfs-area(Int $facet-index) {
return if @facets[$facet-index].visited-time == $time-step;
@facets[$facet-index].visited-time = $time-step;
$!surface-area += ( @facets[$facet-index].plane.normal ).magnitude / 2;
do for ^3 -> $i { self!dfs-area(@facets[$facet-index].neighbors[$i]) }
}
method get-horizon(Int $facet-index, Vector $point, @res-deleted) {
return 0 unless is-above($point, @facets[$facet-index].plane);
return -1 if @facets[$facet-index].visited-time == $time-step;
@facets[$facet-index].visited-time = $time-step;
@facets[$facet-index].is-deleted = True;
@res-deleted.push(@facets[$facet-index].id);
my $result = -2;
for 0..2 -> $neighbor-index {
my $neighbor-facet-index = @facets[$facet-index].neighbors[$neighbor-index];
my $horizon = self.get-horizon($neighbor-facet-index, $point, @res-deleted);
if $horizon == 0 {
my $a = @facets[$facet-index].plane.vector-id($neighbor-index);
my $b = @facets[$facet-index].plane.vector-id(($neighbor-index + 1) % 3);
for (0, 1) -> $idx {
my $point-id = $idx == 0 ?? $a !! $b;
my $neighbor-facet = $neighbor-facet-index;
if @visit-time[$point-id] != $time-step {
@visit-time[$point-id] = $time-step;
@edges[0][$point-id].net-id = $idx == 0 ?? $b !! $a;
@edges[0][$point-id].face-id = $neighbor-facet;
} else {
@edges[1][$point-id].net-id = $idx == 0 ?? $b !! $a;
@edges[1][$point-id].face-id = $neighbor-facet;
}
}
$result = $a;
} elsif $horizon != -1 && $horizon != -2 {
$result = $horizon;
}
}
return $result
}
}
sub get-initial-hull(@points, Int $total-points) {
# Find extreme points in each dimension
my @extremes = @points[1] xx 6;
for @points[1..$total-points] -> $point {
@extremes[0] = $point if is-greater-than($point.x, @extremes[0].x);
@extremes[1] = $point if is-greater-than(@extremes[1].x, $point.x);
@extremes[2] = $point if is-greater-than($point.y, @extremes[2].y);
@extremes[3] = $point if is-greater-than(@extremes[3].y, $point.y);
@extremes[4] = $point if is-greater-than($point.z, @extremes[4].z);
@extremes[5] = $point if is-greater-than(@extremes[5].z, $point.z);
}
# Find furthest pair
my ($extreme0, $extreme1) = @extremes[0,1];
for ^6 -> $i {
for $i+1..^6 -> $j {
my $distance = distance-point-point(@extremes[$i], @extremes[$j]);
if is-greater-than($distance, distance-point-point($extreme0, $extreme1)) {
($extreme0, $extreme1) = @extremes[$i], @extremes[$j]
}
}
}
# Find furthest point from line
my $line = Line.new(u => $extreme0, v => $extreme1);
my $extreme2 = @extremes[0];
for @extremes -> $point {
if is-greater-than(distance-point-line($point, $line), distance-point-line($extreme2, $line)) {
$extreme2 = $point;
}
}
# Find furthest point from plane
my $extreme3 = @points[1];
my $base-plane = Plane.new(u => $extreme0, v => $extreme1, w => $extreme2);
for 1..$total-points -> $i {
my $dist1 = abs(distance-point-plane(@points[$i], $base-plane));
my $dist2 = abs(distance-point-plane($extreme3, $base-plane));
if is-greater-than($dist1, $dist2) { $extreme3 = @points[$i] }
}
# Ensure proper orientation
if is-greater-than(0, distance-point-plane($extreme3, $base-plane)) {
($extreme1, $extreme2) = $extreme2, $extreme1;
}
# Create 4 initial facets
my @facet-indices = ^4;
for ^4 -> $i {
@facets.push(Facet.new(id => @facets.elems));
@facet-indices[$i] = @facets.end;
}
@facets[@facet-indices[0]].plane = Plane.new(u => $extreme0, v => $extreme2, w => $extreme1);
@facets[@facet-indices[1]].plane = Plane.new(u => $extreme0, v => $extreme1, w => $extreme3);
@facets[@facet-indices[2]].plane = Plane.new(u => $extreme1, v => $extreme2, w => $extreme3);
@facets[@facet-indices[3]].plane = Plane.new(u => $extreme2, v => $extreme0, w => $extreme3);
@facets[@facet-indices[0]].neighbors = [@facet-indices[3], @facet-indices[2], @facet-indices[1]];
@facets[@facet-indices[1]].neighbors = [@facet-indices[0], @facet-indices[2], @facet-indices[3]];
@facets[@facet-indices[2]].neighbors = [@facet-indices[0], @facet-indices[3], @facet-indices[1]];
@facets[@facet-indices[3]].neighbors = [@facet-indices[0], @facet-indices[1], @facet-indices[2]];
# Assign remaining points to facets
for @points[1..$total-points] -> $point {
next if any(($extreme0,$extreme1,$extreme2,$extreme3)>>.&{$point.equals($_)});
for ^4 -> $j {
if is-above($point, @facets[@facet-indices[$j]].plane) {
@hull-points[@facet-indices[$j]].push($point);
last;
}
}
}
return ConvexHull3D.new(index => @facet-indices[0])
}
sub quick-hull-Td(@points, Int $total-points) {
my $hull = get-initial-hull(@points, $total-points);
# Initialize processing queue
@queue = [$hull.index];
@queue.append(@facets[$hull.index].neighbors[^3]);
my $new-horizon = 0;
while @queue {
my $next-facet-index = @queue.shift;
if @facets[$next-facet-index].is-deleted || !@hull-points[$next-facet-index] {
$new-horizon = $next-facet-index unless @facets[$next-facet-index].is-deleted;
next;
}
# Find furthest point from current facet
my $point = @hull-points[$next-facet-index][0];
for @hull-points[$next-facet-index] -> $vec {
if is-greater-than(distance-point-plane($vec, @facets[$next-facet-index].plane), distance-point-plane($point, @facets[$next-facet-index].plane)) {
$point = $vec;
}
}
# Find horizon
$time-step++;
@res-facets-deleted = ();
my $horizon = $hull.get-horizon($next-facet-index, $point, @res-facets-deleted);
# Build new faces around horizon
$time-step++;
my ($from, $last-facet-index, $first-facet-index) = -1 xx 3;
while @visit-time[$horizon] != $time-step {
@visit-time[$horizon] = $time-step;
my ($net, $facet-index, $new-facet-index);
if @edges[0][$horizon].net-id == $from {
$net = @edges[1][$horizon].net-id;
$facet-index = @edges[1][$horizon].face-id;
} else {
$net = @edges[0][$horizon].net-id;
$facet-index = @edges[0][$horizon].face-id;
}
# Find indices on facet f
my ($point-index1, $point-index2) = 0, 0;
for ^3 -> $i {
$point-index1 = $i if @points[$horizon].equals(@facets[$facet-index].plane.vector-at($i));
$point-index2 = $i if @points[$net].equals(@facets[$facet-index].plane.vector-at($i));
}
($point-index1, $point-index2) = $point-index2, $point-index1 if ($point-index1 + 1) % 3 != $point-index2;
# Create new facet
@facets.push(Facet.new(
id => @facets.elems,
plane => Plane.new(
u => @facets[$facet-index].plane.vector-at($point-index2),
v => @facets[$facet-index].plane.vector-at($point-index1),
w => $point
)
));
$new-facet-index = @facets.end;
@res-facets-new.push($new-facet-index);
@facets[$new-facet-index].neighbors[0] = $facet-index;
@facets[$facet-index].neighbors[$point-index1] = $new-facet-index;
if $last-facet-index >= 0 { # Link with previous new facet
if @facets[$new-facet-index].plane.vector-at(1).equals(@facets[$last-facet-index].plane.u) {
@facets[$new-facet-index].neighbors[1] = $last-facet-index;
@facets[$last-facet-index].neighbors[2] = $new-facet-index;
} else {
@facets[$new-facet-index].neighbors[2] = $last-facet-index;
@facets[$last-facet-index].neighbors[1] = $new-facet-index;
}
} else {
$first-facet-index = $new-facet-index;
}
($last-facet-index, $from, $horizon) = $new-facet-index, $horizon, $net
}
# Close the loop
if @facets[$first-facet-index].plane.vector-at(1).equals(@facets[$last-facet-index].plane.u) {
@facets[$first-facet-index].neighbors[1] = $last-facet-index;
@facets[$last-facet-index].neighbors[2] = $first-facet-index;
} else {
@facets[$first-facet-index].neighbors[2] = $last-facet-index;
@facets[$last-facet-index].neighbors[1] = $first-facet-index;
}
# Collect points from deleted facets
@res-points = gather for @res-facets-deleted -> $facet-id {
take @hull-points[$facet-id];
@hull-points[$facet-id] = [];
}
# Reassign points to new facets
for @res-points -> $vec {
next if $vec.equals($point);
for @res-facets-new -> $facet-id {
if is-above($vec, @facets[$facet-id].plane) {
@hull-points[$facet-id].push($vec) and last
}
}
}
@queue.append(@res-facets-new); # Enqueue new faces for processing
}
$hull.index = $new-horizon;
return $hull
}
# Example: a tetrahedron
(my @points = Vector.new).append( # placeholder at index 0
Vector.new(x => 0, y => 0, z => 0, id => 1),
Vector.new(x => 1, y => 0, z => 0, id => 2),
Vector.new(x => 0, y => 1, z => 0, id => 3),
Vector.new(x => 0, y => 0, z => 1, id => 4)
);
say quick-hull-Td(@points, 4).get-surface-area.fmt('%.3f');
You may Attempt This Online!
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
You may Attempt This Online!
import scala.collection.mutable
import scala.collection.mutable.{ArrayBuffer, ListBuffer}
import scala.math._
import scala.util.control.Breaks._
object QuickHull3D {
private val MAXN = 2500
private val EPS = 1e-8
// Mutable state (keeping similar structure to original)
private val facets = ArrayBuffer[Facet]()
private val pointsPerFacet = ArrayBuffer[ArrayBuffer[Vect]]()
private var TIME = 0
private val edges = Array.ofDim[Edge](2, MAXN)
private val vistime = Array.ofDim[Int](MAXN)
private val queue = mutable.Queue[Int]()
private val resfnew = ArrayBuffer[Int]()
private val resfdel = ArrayBuffer[Int]()
private val respt = ArrayBuffer[Vect]()
// Utility functions (moved up so they're available to case classes)
private def eq(a: Double, b: Double): Boolean = abs(a - b) < EPS
private def gtr(a: Double, b: Double): Boolean = a - b > EPS
private def abs(x: Double): Double = if (x < 0) -x else x
def main(args: Array[String]): Unit = {
preConvexHulls()
// Example: unit tetrahedron
val points = Array(
null, // index 0 unused
Vect(0, 0, 0, 1),
Vect(1, 0, 0, 2),
Vect(0, 1, 0, 3),
Vect(0, 0, 1, 4)
)
val hull = quickHull3D(points, 4)
println(f"${hull.getSurfaceArea}%.3f")
}
private def preConvexHulls(): Unit = {
// Reserve index 0
pointsPerFacet += ArrayBuffer[Vect]() // dummy point list[0]
facets += Facet() // dummy facet[0]
// Initialize edge array
for (i <- 0 until 2; j <- 0 until MAXN) {
edges(i)(j) = Edge()
}
}
case class Vect(x: Double, y: Double, z: Double, id: Int) {
def -(other: Vect): Vect = Vect(x - other.x, y - other.y, z - other.z, 0)
def cross(other: Vect): Vect = Vect(
y * other.z - z * other.y,
z * other.x - x * other.z,
x * other.y - y * other.x,
0
)
def dot(other: Vect): Double = x * other.x + y * other.y + z * other.z
def magnitude: Double = sqrt(x * x + y * y + z * z)
def ~=(other: Vect): Boolean =
QuickHull3D.eq(x, other.x) && QuickHull3D.eq(y, other.y) && QuickHull3D.eq(z, other.z)
}
case class Line(u: Vect, v: Vect)
case class Plane(u: Vect = Vect(0, 0, 0, 0), v: Vect = Vect(0, 0, 0, 0), w: Vect = Vect(0, 0, 0, 0)) {
def normal: Vect = (v - u).cross(w - u)
def vecAt(i: Int): Vect = i match {
case 0 => u
case 1 => v
case _ => w
}
def vecId(i: Int): Int = vecAt(i).id
}
case class Facet(
var neighbors: Array[Int] = Array(0, 0, 0),
var id: Int = 0,
var visitTime: Int = 0,
var isDeleted: Boolean = false,
var plane: Plane = Plane()
)
case class Edge(var netId: Int = 0, var facetId: Int = 0)
class ConvexHulls3d(var index: Int) {
private var surfaceArea: Double = 0.0
private def dfsArea(f: Int): Unit = {
val facet = facets(f)
if (facet.visitTime == TIME) return
facet.visitTime = TIME
val nrm = facet.plane.normal
surfaceArea += nrm.magnitude / 2.0
facet.neighbors.foreach(dfsArea)
}
def getSurfaceArea: Double = {
if (gtr(surfaceArea, 0.0)) return surfaceArea
TIME += 1
dfsArea(index)
surfaceArea
}
def getHorizon(f: Int, p: Vect, resDel: ArrayBuffer[Int]): Int = {
val facet = facets(f)
if (!isAbove(p, facet.plane)) return 0
if (facet.visitTime == TIME) return -1
facet.visitTime = TIME
facet.isDeleted = true
resDel += facet.id
var ret = -2
for (i <- facet.neighbors.indices) {
val ni = facet.neighbors(i)
val r = getHorizon(ni, p, resDel)
if (r == 0) {
val a = facet.plane.vecId(i)
val b = facet.plane.vecId((i + 1) % 3)
for (idx <- 0 until 2) {
val pt = if (idx == 0) a else b
val facetId = ni
if (vistime(pt) != TIME) {
vistime(pt) = TIME
edges(0)(pt).netId = if (idx == 0) b else a
edges(0)(pt).facetId = facetId
} else {
edges(1)(pt).netId = if (idx == 0) b else a
edges(1)(pt).facetId = facetId
}
}
ret = a
} else if (r != -1 && r != -2) {
ret = r
}
}
ret
}
}
private def distPointPlane(v: Vect, p: Plane): Double = {
val num = (v - p.u).dot(p.normal)
val den = p.normal.magnitude
num / den
}
private def distPointLine(v: Vect, l: Line): Double = {
val d = (v - l.u).magnitude
if (d == 0) return 0
(l.v - l.u).cross(v - l.u).magnitude / (l.v - l.u).magnitude
}
private def distPointPoint(a: Vect, b: Vect): Double = (a - b).magnitude
private def isAbove(v: Vect, p: Plane): Boolean =
gtr((v - p.u).dot(p.normal), 0)
private def getStart(points: Array[Vect], totalPoints: Int): ConvexHulls3d = {
// Find extreme points
val extremes = Array.fill(6)(points(1))
for (i <- 1 to totalPoints) {
val v = points(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
}
// Find furthest pair
var (s0, s1) = (extremes(0), extremes(1))
for (i <- extremes.indices; j <- i + 1 until extremes.length) {
val d = distPointPoint(extremes(i), extremes(j))
if (gtr(d, distPointPoint(s0, s1))) {
s0 = extremes(i)
s1 = extremes(j)
}
}
// Find point furthest from line
val line = Line(s0, s1)
var s2 = extremes(0)
for (extreme <- extremes) {
if (gtr(distPointLine(extreme, line), distPointLine(s2, line))) {
s2 = extreme
}
}
// Find point furthest from plane
var s3 = points(1)
val basePlane = Plane(s0, s1, s2)
for (i <- 1 to totalPoints) {
val d1 = abs(distPointPlane(points(i), basePlane))
val d2 = abs(distPointPlane(s3, basePlane))
if (gtr(d1, d2)) s3 = points(i)
}
if (gtr(0, distPointPlane(s3, basePlane))) {
val tmp = s1; s1 = s2; s2 = tmp
}
// Create 4 initial facets
val f = Array.ofDim[Int](4)
for (i <- 0 until 4) {
val facet = Facet()
facet.id = facets.length
facets += facet
f(i) = facets.length - 1
}
facets(f(0)).plane = Plane(s0, s2, s1)
facets(f(1)).plane = Plane(s0, s1, s3)
facets(f(2)).plane = Plane(s1, s2, s3)
facets(f(3)).plane = Plane(s2, s0, s3)
facets(f(0)).neighbors = Array(f(3), f(2), f(1))
facets(f(1)).neighbors = Array(f(0), f(2), f(3))
facets(f(2)).neighbors = Array(f(0), f(3), f(1))
facets(f(3)).neighbors = Array(f(0), f(1), f(2))
// Prepare point lists for the 4 facets
for (_ <- 0 until 4) {
pointsPerFacet += ArrayBuffer[Vect]()
}
// Assign points to facets
for (i <- 1 to totalPoints) {
val v = points(i)
if (!(v ~= s0) && !(v ~= s1) && !(v ~= s2) && !(v ~= s3)) {
breakable({
for (j <- 0 until 4) {
if (isAbove(v, facets(f(j)).plane)) {
pointsPerFacet(f(j)) += v
break
}
}
});
}
}
new ConvexHulls3d(f(0))
}
def quickHull3D(points: Array[Vect], totalPoints: Int): ConvexHulls3d = {
val hull = getStart(points, totalPoints)
queue.clear()
queue.enqueue(hull.index)
facets(hull.index).neighbors.foreach(queue.enqueue)
var snew = 0
while (queue.nonEmpty) {
breakable({
val nf = queue.dequeue()
val facet = facets(nf)
if (facet.isDeleted || pointsPerFacet(nf).isEmpty) {
if (!facet.isDeleted) snew = nf
break//continue
}
// Find farthest point from facet plane
var p = pointsPerFacet(nf).head
for (v <- pointsPerFacet(nf)) {
if (gtr(distPointPlane(v, facet.plane), distPointPlane(p, facet.plane))) {
p = v
}
}
// Find horizon
TIME += 1
resfdel.clear()
val s = hull.getHorizon(nf, p, resfdel)
// Build new faces
resfnew.clear()
TIME += 1
var from = -1
var lastf = -1
var fstf = -1
var currentS = s
while (vistime(currentS) != TIME) {
vistime(currentS) = TIME
val (net, fidx) = if (edges(0)(currentS).netId == from) {
(edges(1)(currentS).netId, edges(1)(currentS).facetId)
} else {
(edges(0)(currentS).netId, edges(0)(currentS).facetId)
}
// Find indices on facet fidx
val facetFidx = facets(fidx)
var (pt1, pt2) = (0, 0)
for (i <- 0 until 3) {
if (points(currentS) ~= facetFidx.plane.vecAt(i)) pt1 = i
if (points(net) ~= facetFidx.plane.vecAt(i)) pt2 = i
}
if ((pt1 + 1) % 3 != pt2) {
val tmp = pt1; pt1 = pt2; pt2 = tmp
}
// Create new facet
val newFacet = Facet()
newFacet.id = facets.length
newFacet.plane = Plane(
facetFidx.plane.vecAt(pt2),
facetFidx.plane.vecAt(pt1),
p
)
facets += newFacet
pointsPerFacet += ArrayBuffer[Vect]()
val fnew = facets.length - 1
resfnew += fnew
// Link neighborhoods
newFacet.neighbors(0) = fidx
facetFidx.neighbors(pt1) = fnew
if (lastf >= 0) {
val newPlane = newFacet.plane
val lastPlane = facets(lastf).plane
if (newPlane.vecAt(1) ~= lastPlane.u) {
newFacet.neighbors(1) = lastf
facets(lastf).neighbors(2) = fnew
} else {
newFacet.neighbors(2) = lastf
facets(lastf).neighbors(1) = fnew
}
} else {
fstf = fnew
}
lastf = fnew
from = currentS
currentS = net
}
// Close the loop
val firstFacet = facets(fstf)
val lastFacet = facets(lastf)
if (firstFacet.plane.vecAt(1) ~= lastFacet.plane.u) {
firstFacet.neighbors(1) = lastf
lastFacet.neighbors(2) = fstf
} else {
firstFacet.neighbors(2) = lastf
lastFacet.neighbors(1) = fstf
}
// Collect deleted points
respt.clear()
for (fid <- resfdel) {
respt ++= pointsPerFacet(fid)
pointsPerFacet(fid).clear()
}
// Reassign points
for (v <- respt if v != p) {
for (fid <- resfnew) {
if (isAbove(v, facets(fid).plane)) {
pointsPerFacet(fid) += v
break
}
}
}
// Enqueue new faces
resfnew.foreach(queue.enqueue)
});
}
hull.index = snew
hull
}
}
- Output:
2.366
import Foundation
// ------------------------------------------------------------
// MARK: - Constants & helpers
// ------------------------------------------------------------
let MAX_SIZE = 2_500
let EPSILON = 1e-8
@inline(__always) func isGreaterThan(_ a: Double, _ b: Double) -> Bool {
a - b > EPSILON
}
@inline(__always) func isEqual(_ a: Double, _ b: Double) -> Bool {
abs(a - b) < EPSILON
}
// ------------------------------------------------------------
// MARK: - Basic geometry
// ------------------------------------------------------------
struct Vector {
var x: Double
var y: Double
var z: Double
var id: Int // the original index of the point
init(_ x: Double = 0, _ y: Double = 0, _ z: Double = 0, _ id: Int = 0) {
self.x = x; self.y = y; self.z = z; self.id = id
}
func subtract(_ other: Vector) -> Vector {
Vector(x - other.x, y - other.y, z - other.z)
}
func vectorProduct(_ other: Vector) -> Vector {
Vector(
y * other.z - z * other.y,
z * other.x - x * other.z,
x * other.y - y * other.x
)
}
func scalarProduct(_ other: Vector) -> Double {
x * other.x + y * other.y + z * other.z
}
func magnitude() -> Double {
sqrt(x * x + y * y + z * z)
}
func equals(_ other: Vector) -> Bool {
isEqual(x, other.x) && isEqual(y, other.y) && isEqual(z, other.z)
}
}
struct Line {
var u: Vector
var v: Vector
init(_ u: Vector, _ v: Vector) {
self.u = u; self.v = v
}
}
struct Plane {
var u: Vector
var v: Vector
var w: Vector
init(_ u: Vector = Vector(),
_ v: Vector = Vector(),
_ w: Vector = Vector()) {
self.u = u; self.v = v; self.w = w
}
func normal() -> Vector {
v.subtract(u).vectorProduct(w.subtract(u))
}
func vector(at i: Int) -> Vector {
switch i {
case 0: return u
case 1: return v
case 2: return w
default: return Vector()
}
}
func vectorID(at i: Int) -> Int {
vector(at: i).id
}
}
// ------------------------------------------------------------
// MARK: - Facet, Edge and global containers
// ------------------------------------------------------------
class Facet {
var N: [Int] // neighbour facet indices (always size 3)
var id: Int
var visitedTime: Int = 0
var isDeleted: Bool = false
var plane: Plane
init(id: Int) {
self.id = id
self.N = [Int](repeating: 0, count: 3)
self.plane = Plane()
}
init(id: Int, plane: Plane) {
self.id = id
self.N = [Int](repeating: 0, count: 3)
self.plane = plane
}
init() {
self.id = 0
self.N = [Int](repeating: 0, count: 3)
self.plane = Plane()
}
}
struct Edge {
var netID: Int = 0
var faceID: Int = 0
}
// Global containers (mirroring the C++ version)
var facets: [Facet] = []
var hullPoints: [[Vector]] = [] // points that belong to each facet
var timeStep: Int = 0
var edges: [[Edge]] = Array(
repeating: Array(repeating: Edge(), count: MAX_SIZE),
count: 2
)
var visitTime: [Int] = Array(repeating: 0, count: MAX_SIZE)
var queue: [Int] = []
var resFNew: [Int] = []
var resFDel: [Int] = []
var resPt: [Vector] = []
// ------------------------------------------------------------
// MARK: - Distance / side tests
// ------------------------------------------------------------
func distancePointPlane(_ pt: Vector, _ plane: Plane) -> Double {
let num = pt.subtract(plane.u).scalarProduct(plane.normal())
let den = plane.normal().magnitude()
return num / den
}
func distancePointLine(_ pt: Vector, _ line: Line) -> Double {
let len = pt.subtract(line.u).magnitude()
if len == 0 { return 0.0 }
return line.v.subtract(line.u)
.vectorProduct(pt.subtract(line.u))
.magnitude() / line.v.subtract(line.u).magnitude()
}
func distancePointPoint(_ a: Vector, _ b: Vector) -> Double {
a.subtract(b).magnitude()
}
func isAbove(_ point: Vector, _ plane: Plane) -> Bool {
isGreaterThan(point.subtract(plane.u).scalarProduct(plane.normal()), 0.0)
}
// ------------------------------------------------------------
// MARK: - Convex hull object (the driver)
// ------------------------------------------------------------
final class ConvexHulls3D {
var index: Int // the “seed” facet for surface‑area computation
private var surfaceArea: Double = 0.0
init(index: Int) {
self.index = index
}
// ----------------------------------------------------------------
// public: surface area – cached after first computation
// ----------------------------------------------------------------
func getSurfaceArea() -> Double {
if isGreaterThan(surfaceArea, 0.0) { return surfaceArea }
timeStep += 1
dfsArea(from: index)
return surfaceArea
}
// ----------------------------------------------------------------
// private: depth‑first search of the facet adjacency graph
// ----------------------------------------------------------------
private func dfsArea(from f: Int) {
if facets[f].visitedTime == timeStep { return }
facets[f].visitedTime = timeStep
let normal = facets[f].plane.normal()
surfaceArea += normal.magnitude() / 2.0
for i in 0..<3 { dfsArea(from: facets[f].N[i]) }
}
// ----------------------------------------------------------------
// Horizon discovery – returns the first horizon edge (or 0 / -1 / -2)
// ----------------------------------------------------------------
func getHorizon(_ f: Int, _ point: Vector, _ resDel: inout [Int]) -> Int {
let facet = facets[f]
if !isAbove(point, facet.plane) { return 0 }
if facet.visitedTime == timeStep { return -1 }
facet.visitedTime = timeStep
facet.isDeleted = true
resDel.append(facet.id)
var result = -2
for i in 0..<3 {
let neigh = facet.N[i]
let horizon = getHorizon(neigh, point, &resDel)
if horizon == 0 {
// edge (a,b) of facet f opposite vertex i
let a = facet.plane.vectorID(at: i)
let b = facet.plane.vectorID(at: (i + 1) % 3)
for idx in 0..<2 {
let pt = (idx == 0) ? a : b
let fID = neigh
if visitTime[pt] != timeStep {
visitTime[pt] = timeStep
edges[0][pt].netID = (idx == 0) ? b : a
edges[0][pt].faceID = fID
} else {
edges[1][pt].netID = (idx == 0) ? b : a
edges[1][pt].faceID = fID
}
}
result = a
} else if horizon != -1 && horizon != -2 {
result = horizon
}
}
return result
}
}
// ------------------------------------------------------------
// MARK: - Preparations (allocate structures)
// ------------------------------------------------------------
func prepareConvexHulls() {
// reserve index 0 – the dummy entry used by the C++ program
hullPoints.append([])
facets.append(Facet())
// initialise edges matrix
for i in 0..<2 {
for j in 0..<MAX_SIZE {
edges[i][j] = Edge()
}
}
}
// ------------------------------------------------------------
// MARK: - Build the initial tetrahedron
// ------------------------------------------------------------
func getInitialHull(_ points: [Vector], _ totalPoints: Int) -> ConvexHulls3D {
// ---- 1. extreme points on each axis -----------------------------------------
var extremes = Array(repeating: points[1], count: 6)
for i in 1...totalPoints {
let p = points[i]
if isGreaterThan(p.x, extremes[0].x) { extremes[0] = p }
if isGreaterThan(extremes[1].x, p.x) { extremes[1] = p }
if isGreaterThan(p.y, extremes[2].y) { extremes[2] = p }
if isGreaterThan(extremes[3].y, p.y) { extremes[3] = p }
if isGreaterThan(p.z, extremes[4].z) { extremes[4] = p }
if isGreaterThan(extremes[5].z, p.z) { extremes[5] = p }
}
// ---- 2. farthest pair -------------------------------------------------------
var extreme0 = extremes[0]
var extreme1 = extremes[1]
for i in 0..<6 {
for j in i+1..<6 {
let d = distancePointPoint(extremes[i], extremes[j])
if isGreaterThan(d, distancePointPoint(extreme0, extreme1)) {
extreme0 = extremes[i]
extreme1 = extremes[j]
}
}
}
// ---- 3. farthest from that line ---------------------------------------------
let line = Line(extreme0, extreme1)
var extreme2 = extremes[0]
for i in 0..<6 {
if isGreaterThan(distancePointLine(extremes[i], line),
distancePointLine(extreme2, line)) {
extreme2 = extremes[i]
}
}
// ---- 4. farthest from the plane formed by the three points ------------------
var extreme3 = points[1] // any point that is not a dummy
let basePlane = Plane(extreme0, extreme1, extreme2)
for i in 1...totalPoints {
let d1 = abs(distancePointPlane(points[i], basePlane))
let d2 = abs(distancePointPlane(extreme3, basePlane))
if isGreaterThan(d1, d2) { extreme3 = points[i] }
}
// Ensure the orientation of the base plane is outward
if isGreaterThan(0.0, distancePointPlane(extreme3, basePlane)) {
swap(&extreme1, &extreme2)
}
// ---- 5. create the four seed facets -----------------------------------------
var f = [Int](repeating: 0, count: 4)
for i in 0..<4 {
facets.append(Facet(id: facets.count))
f[i] = facets.count - 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]]
// ---- 6. allocate a hull‑point list for each facet ---------------------------
for _ in 0..<4 { hullPoints.append([]) }
// ---- 7. distribute every remaining point to the facet that sees it ---------
for i in 1...totalPoints {
let p = points[i]
if p.equals(extreme0) || p.equals(extreme1) ||
p.equals(extreme2) || p.equals(extreme3) { continue }
for j in 0..<4 {
if isAbove(p, facets[f[j]].plane) {
hullPoints[f[j]].append(p)
break
}
}
}
return ConvexHulls3D(index: f[0])
}
// ------------------------------------------------------------
// MARK: - Main QuickHull driver
// ------------------------------------------------------------
func quickHull3D(_ points: [Vector], _ totalPoints: Int) -> ConvexHulls3D {
var hull = getInitialHull(points, totalPoints)
// initialise the processing queue
queue.removeAll()
queue.append(hull.index)
for i in 0..<3 { queue.append(facets[hull.index].N[i]) }
var newHorizon = 0
while !queue.isEmpty {
let nf = queue.removeFirst()
if facets[nf].isDeleted || hullPoints[nf].isEmpty {
if !facets[nf].isDeleted { newHorizon = nf }
continue
}
// ----- farthest point of the current facet ------------------------------
var farthest = hullPoints[nf][0]
for v in hullPoints[nf] {
if isGreaterThan(
distancePointPlane(v, facets[nf].plane),
distancePointPlane(farthest, facets[nf].plane)) {
farthest = v
}
}
// ----- horizon -----------------------------------------------------------
timeStep += 1
resFDel.removeAll()
let horizon = hull.getHorizon(nf, farthest, &resFDel)
// ----- build new facets --------------------------------------------------
resFNew.removeAll()
timeStep += 1
var from = -1
var lastF = -1
var firstF = -1
var curHorizon = horizon
while visitTime[curHorizon] != timeStep {
visitTime[curHorizon] = timeStep
// decide which edge instance (0/1) belongs to the current direction
let net: Int
let fID: Int
if edges[0][curHorizon].netID == from {
net = edges[1][curHorizon].netID
fID = edges[1][curHorizon].faceID
} else {
net = edges[0][curHorizon].netID
fID = edges[0][curHorizon].faceID
}
// locate the two vertices of facet f that correspond to
// - curHorizon (the “horizon” vertex)
// - net (the neighbour vertex)
var pt1 = 0
var pt2 = 0
for i in 0..<3 {
if points[curHorizon].equals(facets[fID].plane.vector(at: i)) {
pt1 = i
}
if points[net].equals(facets[fID].plane.vector(at: i)) {
pt2 = i
}
}
// make sure that (pt1,pt2) are ordered clockwise as in the C++ code
if (pt1 + 1) % 3 != pt2 {
swap(&pt1, &pt2)
}
// create the new facet
let newPlane = Plane(
facets[fID].plane.vector(at: pt2),
facets[fID].plane.vector(at: pt1),
farthest
)
facets.append(Facet(id: facets.count, plane: newPlane))
let newF = facets.count - 1
hullPoints.append([]) // allocate its point list
resFNew.append(newF)
// link the new facet with the old one
facets[newF].N[0] = fID
facets[fID].N[pt1] = newF
if lastF >= 0 {
// connect to the previous new facet (closing the fan)
if facets[newF].plane.vector(at: 1).equals(facets[lastF].plane.u) {
facets[newF].N[1] = lastF
facets[lastF].N[2] = newF
} else {
facets[newF].N[2] = lastF
facets[lastF].N[1] = newF
}
} else {
firstF = newF
}
lastF = newF
from = curHorizon
curHorizon = net
}
// ----- close the fan ----------------------------------------------------
if facets[firstF].plane.vector(at: 1).equals(facets[lastF].plane.u) {
facets[firstF].N[1] = lastF
facets[lastF].N[2] = firstF
} else {
facets[firstF].N[2] = lastF
facets[lastF].N[1] = firstF
}
// ----- collect all points that were attached to deleted facets ----------
resPt.removeAll()
for fID in resFDel {
resPt.append(contentsOf: hullPoints[fID])
hullPoints[fID].removeAll()
}
// ----- re‑assign those points to the newly created facets ---------------
for p in resPt where !p.equals(farthest) {
for fID in resFNew {
if isAbove(p, facets[fID].plane) {
hullPoints[fID].append(p)
break
}
}
}
// ----- enqueue newly created facets --------------------------------------
for fID in resFNew { queue.append(fID) }
}
hull.index = newHorizon
return hull
}
// ------------------------------------------------------------
// MARK: - Example / entry point
// ------------------------------------------------------------
func runExample() {
prepareConvexHulls()
// Example: a tetrahedron (the same data as the C++ main)
let n = 4
var points = [Vector](repeating: Vector(), count: n + 1) // points[0] unused
points[1] = Vector(0, 0, 0, 1)
points[2] = Vector(1, 0, 0, 2)
points[3] = Vector(0, 1, 0, 3)
points[4] = Vector(0, 0, 1, 4)
let hull = quickHull3D(points, n)
let area = hull.getSurfaceArea()
print(String(format: "%.3f", area))
}
// ------------------------------------------------------------
// MARK: - Run
// ------------------------------------------------------------
runExample()
- Output:
2.366
#!/usr/bin/env tclsh
# Constants
set MAX_SIZE 2500
set EPSILON 0.00000001
# Global variables
set facets [list]
set hull_points [list]
set time_step 0
set edges [list [list] [list]]
set visit_time [list]
set queue [list]
set resfnew [list]
set resfdel [list]
set respt [list]
# Numerical utilities
proc is_greater_than {a b} {
global EPSILON
return [expr {($a - $b) > $EPSILON}]
}
proc is_equal {a b} {
global EPSILON
return [expr {abs($a - $b) < $EPSILON}]
}
# Vector class implementation
proc vector_create {x y z id} {
return [dict create x $x y $y z $z id $id]
}
proc vector_subtract {vec1 vec2} {
set x1 [dict get $vec1 x]
set y1 [dict get $vec1 y]
set z1 [dict get $vec1 z]
set x2 [dict get $vec2 x]
set y2 [dict get $vec2 y]
set z2 [dict get $vec2 z]
return [vector_create [expr {$x1 - $x2}] [expr {$y1 - $y2}] [expr {$z1 - $z2}] 0]
}
proc vector_cross_product {vec1 vec2} {
set x1 [dict get $vec1 x]
set y1 [dict get $vec1 y]
set z1 [dict get $vec1 z]
set x2 [dict get $vec2 x]
set y2 [dict get $vec2 y]
set z2 [dict get $vec2 z]
set x [expr {$y1 * $z2 - $z1 * $y2}]
set y [expr {$z1 * $x2 - $x1 * $z2}]
set z [expr {$x1 * $y2 - $y1 * $x2}]
return [vector_create $x $y $z 0]
}
proc vector_dot_product {vec1 vec2} {
set x1 [dict get $vec1 x]
set y1 [dict get $vec1 y]
set z1 [dict get $vec1 z]
set x2 [dict get $vec2 x]
set y2 [dict get $vec2 y]
set z2 [dict get $vec2 z]
return [expr {$x1 * $x2 + $y1 * $y2 + $z1 * $z2}]
}
proc vector_magnitude {vec} {
set x [dict get $vec x]
set y [dict get $vec y]
set z [dict get $vec z]
return [expr {sqrt($x * $x + $y * $y + $z * $z)}]
}
proc vector_equals {vec1 vec2} {
set x1 [dict get $vec1 x]
set y1 [dict get $vec1 y]
set z1 [dict get $vec1 z]
set x2 [dict get $vec2 x]
set y2 [dict get $vec2 y]
set z2 [dict get $vec2 z]
return [expr {[is_equal $x1 $x2] && [is_equal $y1 $y2] && [is_equal $z1 $z2]}]
}
# Line class implementation
proc line_create {u v} {
return [dict create u $u v $v]
}
# Plane class implementation
proc plane_create {u v w} {
return [dict create u $u v $v w $w]
}
proc plane_normal {plane} {
set u [dict get $plane u]
set v [dict get $plane v]
set w [dict get $plane w]
set vec1 [vector_subtract $v $u]
set vec2 [vector_subtract $w $u]
return [vector_cross_product $vec1 $vec2]
}
proc plane_vector_at_index {plane i} {
switch $i {
0 { return [dict get $plane u] }
1 { return [dict get $plane v] }
2 { return [dict get $plane w] }
default { return [vector_create 0 0 0 0] }
}
}
proc plane_vector_id {plane i} {
set vec [plane_vector_at_index $plane $i]
return [dict get $vec id]
}
# Facet class implementation
proc facet_create {id plane} {
return [dict create id $id plane $plane N [list] visited_time 0 is_deleted false]
}
# Edge class implementation
proc edge_create {} {
return [dict create netID 0 faceID 0]
}
# Geometric utilities
proc distance_point_plane {vec plane} {
set u [dict get $plane u]
set normal [plane_normal $plane]
set diff [vector_subtract $vec $u]
set num [vector_dot_product $diff $normal]
set den [vector_magnitude $normal]
return [expr {$num / $den}]
}
proc distance_point_line {vec line} {
set u [dict get $line u]
set v [dict get $line v]
set diff1 [vector_subtract $vec $u]
set length [vector_magnitude $diff1]
if {$length == 0.0} {
return 0.0
}
set diff2 [vector_subtract $v $u]
set cross [vector_cross_product $diff2 $diff1]
set cross_mag [vector_magnitude $cross]
set line_mag [vector_magnitude $diff2]
return [expr {$cross_mag / $line_mag}]
}
proc distance_point_point {a b} {
set diff [vector_subtract $a $b]
return [vector_magnitude $diff]
}
proc is_above {point plane} {
set u [dict get $plane u]
set normal [plane_normal $plane]
set diff [vector_subtract $point $u]
set dot [vector_dot_product $diff $normal]
return [is_greater_than $dot 0.0]
}
# ConvexHull3d class implementation
proc convex_hull_create {index} {
return [dict create index $index surface_area 0.0]
}
proc get_surface_area {hull_ref} {
upvar $hull_ref hull
set surface_area [dict get $hull surface_area]
if {[is_greater_than $surface_area 0.0]} {
return $surface_area
}
global time_step
incr time_step
set index [dict get $hull index]
set area [dfs_area $index]
dict set hull surface_area $area
return $area
}
proc dfs_area {f} {
global facets time_step
set facet [lindex $facets $f]
if {[dict get $facet visited_time] == $time_step} {
return 0.0
}
dict set facet visited_time $time_step
lset facets $f $facet
set plane [dict get $facet plane]
set normal [plane_normal $plane]
set area [expr {[vector_magnitude $normal] / 2.0}]
set N [dict get $facet N]
foreach neighbor $N {
set area [expr {$area + [dfs_area $neighbor]}]
}
return $area
}
proc get_horizon {f point resDel_ref} {
upvar $resDel_ref resDel
global facets time_step
set facet [lindex $facets $f]
set plane [dict get $facet plane]
if {![is_above $point $plane]} {
return 0
}
if {[dict get $facet visited_time] == $time_step} {
return -1
}
dict set facet visited_time $time_step
dict set facet is_deleted true
lset facets $f $facet
lappend resDel [dict get $facet id]
set result -2
set N [dict get $facet N]
for {set i 0} {$i < 3} {incr i} {
set ni [lindex $N $i]
set horizon [get_horizon $ni $point resDel]
if {$horizon == 0} {
set a [plane_vector_id $plane $i]
set b [plane_vector_id $plane [expr {($i + 1) % 3}]]
# Handle edges
global edges visit_time
for {set idx 0} {$idx < 2} {incr idx} {
set pt [expr {$idx == 0 ? $a : $b}]
set facet_id $ni
if {[lindex $visit_time $pt] != $time_step} {
lset visit_time $pt $time_step
set edge [edge_create]
dict set edge netID [expr {$idx == 0 ? $b : $a}]
dict set edge faceID $facet_id
lset edges 0 $pt $edge
} else {
set edge [edge_create]
dict set edge netID [expr {$idx == 0 ? $b : $a}]
dict set edge faceID $facet_id
lset edges 1 $pt $edge
}
}
set result $a
} elseif {$horizon != -1 && $horizon != -2} {
set result $horizon
}
}
return $result
}
proc prepareConvexHulls {} {
global hull_points facets edges MAX_SIZE
# Reserve index 0
lappend hull_points [list]
lappend facets [facet_create 0 [plane_create [vector_create 0 0 0 0] [vector_create 0 0 0 0] [vector_create 0 0 0 0]]]
# Initialize edge vectors
for {set i 0} {$i < 2} {incr i} {
set edge_list [list]
for {set j 0} {$j < $MAX_SIZE} {incr j} {
lappend edge_list [edge_create]
}
lset edges $i $edge_list
}
# Initialize visit_time
global visit_time
for {set i 0} {$i < $MAX_SIZE} {incr i} {
lappend visit_time 0
}
}
proc get_initial_hull {points total_points} {
global facets hull_points
# Find extreme points
set extremes [list]
for {set i 0} {$i < 6} {incr i} {
lappend extremes [lindex $points 1]
}
for {set i 1} {$i <= $total_points} {incr i} {
set point [lindex $points $i]
set x [dict get $point x]
set y [dict get $point y]
set z [dict get $point z]
if {[is_greater_than $x [dict get [lindex $extremes 0] x]]} {
lset extremes 0 $point
}
if {[is_greater_than [dict get [lindex $extremes 1] x] $x]} {
lset extremes 1 $point
}
if {[is_greater_than $y [dict get [lindex $extremes 2] y]]} {
lset extremes 2 $point
}
if {[is_greater_than [dict get [lindex $extremes 3] y] $y]} {
lset extremes 3 $point
}
if {[is_greater_than $z [dict get [lindex $extremes 4] z]]} {
lset extremes 4 $point
}
if {[is_greater_than [dict get [lindex $extremes 5] z] $z]} {
lset extremes 5 $point
}
}
# Find furthest pair
set extreme0 [lindex $extremes 0]
set extreme1 [lindex $extremes 1]
for {set i 0} {$i < 6} {incr i} {
for {set j [expr {$i + 1}]} {$j < 6} {incr j} {
set distance [distance_point_point [lindex $extremes $i] [lindex $extremes $j]]
if {[is_greater_than $distance [distance_point_point $extreme0 $extreme1]]} {
set extreme0 [lindex $extremes $i]
set extreme1 [lindex $extremes $j]
}
}
}
# Find furthest from line
set line [line_create $extreme0 $extreme1]
set extreme2 [lindex $extremes 0]
for {set i 0} {$i < 6} {incr i} {
if {[is_greater_than [distance_point_line [lindex $extremes $i] $line] [distance_point_line $extreme2 $line]]} {
set extreme2 [lindex $extremes $i]
}
}
# Find furthest from plane
set extreme3 [lindex $points 1]
set basePlane [plane_create $extreme0 $extreme1 $extreme2]
for {set i 1} {$i <= $total_points} {incr i} {
set distance1 [expr {abs([distance_point_plane [lindex $points $i] $basePlane])}]
set distance2 [expr {abs([distance_point_plane $extreme3 $basePlane])}]
if {[is_greater_than $distance1 $distance2]} {
set extreme3 [lindex $points $i]
}
}
if {[is_greater_than 0 [distance_point_plane $extreme3 $basePlane]]} {
set temp $extreme1
set extreme1 $extreme2
set extreme2 $temp
}
# Create 4 initial facets
set f [list]
for {set i 0} {$i < 4} {incr i} {
set facet_id [llength $facets]
set plane [plane_create $extreme0 $extreme1 $extreme2]
lappend facets [facet_create $facet_id $plane]
lappend f $facet_id
}
# Set up planes
set f0 [lindex $f 0]
set f1 [lindex $f 1]
set f2 [lindex $f 2]
set f3 [lindex $f 3]
set facet0 [lindex $facets $f0]
dict set facet0 plane [plane_create $extreme0 $extreme2 $extreme1]
lset facets $f0 $facet0
set facet1 [lindex $facets $f1]
dict set facet1 plane [plane_create $extreme0 $extreme1 $extreme3]
lset facets $f1 $facet1
set facet2 [lindex $facets $f2]
dict set facet2 plane [plane_create $extreme1 $extreme2 $extreme3]
lset facets $f2 $facet2
set facet3 [lindex $facets $f3]
dict set facet3 plane [plane_create $extreme2 $extreme0 $extreme3]
lset facets $f3 $facet3
# Set neighbors
set facet0 [lindex $facets $f0]
dict set facet0 N [list $f3 $f2 $f1]
lset facets $f0 $facet0
set facet1 [lindex $facets $f1]
dict set facet1 N [list $f0 $f2 $f3]
lset facets $f1 $facet1
set facet2 [lindex $facets $f2]
dict set facet2 N [list $f0 $f3 $f1]
lset facets $f2 $facet2
set facet3 [lindex $facets $f3]
dict set facet3 N [list $f0 $f1 $f2]
lset facets $f3 $facet3
# Prepare hull_points
for {set i 0} {$i < 4} {incr i} {
lappend hull_points [list]
}
# Assign points
for {set i 1} {$i <= $total_points} {incr i} {
set point [lindex $points $i]
if {[vector_equals $point $extreme0] || [vector_equals $point $extreme1] ||
[vector_equals $point $extreme2] || [vector_equals $point $extreme3]} {
continue
}
for {set j 0} {$j < 4} {incr j} {
set fj [lindex $f $j]
set facet [lindex $facets $fj]
set plane [dict get $facet plane]
if {[is_above $point $plane]} {
set current_points [lindex $hull_points $fj]
lappend current_points $point
lset hull_points $fj $current_points
break
}
}
}
return [convex_hull_create $f0]
}
proc QuickHull3D {points total_points} {
global queue facets hull_points time_step resfdel resfnew respt edges visit_time
set hull [get_initial_hull $points $total_points]
# Initialize queue
set queue [list]
set hull_index [dict get $hull index]
lappend queue $hull_index
set facet [lindex $facets $hull_index]
set N [dict get $facet N]
foreach neighbor $N {
lappend queue $neighbor
}
set new_horizon 0
while {[llength $queue] > 0} {
set nf [lindex $queue 0]
set queue [lrange $queue 1 end]
set facet [lindex $facets $nf]
if {[dict get $facet is_deleted] || [llength [lindex $hull_points $nf]] == 0} {
if {![dict get $facet is_deleted]} {
set new_horizon $nf
}
continue
}
# Find farthest point
set point_list [lindex $hull_points $nf]
set point [lindex $point_list 0]
set plane [dict get $facet plane]
foreach vec $point_list {
if {[is_greater_than [distance_point_plane $vec $plane] [distance_point_plane $point $plane]]} {
set point $vec
}
}
# Find horizon
incr time_step
set resfdel [list]
set horizon [get_horizon $nf $point resfdel]
# Build new faces
set resfnew [list]
incr time_step
set from -1
set last_f -1
set first_f -1
while {[lindex $visit_time $horizon] != $time_step} {
lset visit_time $horizon $time_step
if {[dict get [lindex $edges 0 $horizon] netID] == $from} {
set net [dict get [lindex $edges 1 $horizon] netID]
set f [dict get [lindex $edges 1 $horizon] faceID]
} else {
set net [dict get [lindex $edges 0 $horizon] netID]
set f [dict get [lindex $edges 0 $horizon] faceID]
}
# Create new facet
set new_f [llength $facets]
set facet_f [lindex $facets $f]
set plane_f [dict get $facet_f plane]
# Find indices
set pt1 0
set pt2 0
for {set i 0} {$i < 3} {incr i} {
set vec_at_i [plane_vector_at_index $plane_f $i]
if {[vector_equals [lindex $points $horizon] $vec_at_i]} {
set pt1 $i
}
if {[vector_equals [lindex $points $net] $vec_at_i]} {
set pt2 $i
}
}
if {[expr {($pt1 + 1) % 3}] != $pt2} {
set temp $pt1
set pt1 $pt2
set pt2 $temp
}
set new_plane [plane_create [plane_vector_at_index $plane_f $pt2] [plane_vector_at_index $plane_f $pt1] $point]
lappend facets [facet_create $new_f $new_plane]
lappend hull_points [list]
lappend resfnew $new_f
# Set neighbors
set new_facet [lindex $facets $new_f]
set N [dict get $new_facet N]
lset N 0 $f
dict set new_facet N $N
lset facets $new_f $new_facet
set facet_f [lindex $facets $f]
set N_f [dict get $facet_f N]
lset N_f $pt1 $new_f
dict set facet_f N $N_f
lset facets $f $facet_f
if {$last_f >= 0} {
# Link with previous new facet
set new_facet [lindex $facets $new_f]
set last_facet [lindex $facets $last_f]
set new_plane [dict get $new_facet plane]
set last_plane [dict get $last_facet plane]
if {[vector_equals [plane_vector_at_index $new_plane 1] [dict get $last_plane u]]} {
set N [dict get $new_facet N]
lset N 1 $last_f
dict set new_facet N $N
lset facets $new_f $new_facet
set N [dict get $last_facet N]
lset N 2 $new_f
dict set last_facet N $N
lset facets $last_f $last_facet
} else {
set N [dict get $new_facet N]
lset N 2 $last_f
dict set new_facet N $N
lset facets $new_f $new_facet
set N [dict get $last_facet N]
lset N 1 $new_f
dict set last_facet N $N
lset facets $last_f $last_facet
}
} else {
set first_f $new_f
}
set last_f $new_f
set from $horizon
set horizon $net
}
# Close the loop
if {$first_f >= 0 && $last_f >= 0} {
set first_facet [lindex $facets $first_f]
set last_facet [lindex $facets $last_f]
set first_plane [dict get $first_facet plane]
set last_plane [dict get $last_facet plane]
if {[vector_equals [plane_vector_at_index $first_plane 1] [dict get $last_plane u]]} {
set N [dict get $first_facet N]
lset N 1 $last_f
dict set first_facet N $N
lset facets $first_f $first_facet
set N [dict get $last_facet N]
lset N 2 $first_f
dict set last_facet N $N
lset facets $last_f $last_facet
} else {
set N [dict get $first_facet N]
lset N 2 $last_f
dict set first_facet N $N
lset facets $first_f $first_facet
set N [dict get $last_facet N]
lset N 1 $first_f
dict set last_facet N $N
lset facets $last_f $last_facet
}
}
# Collect deleted points
set respt [list]
foreach f_id $resfdel {
set deleted_points [lindex $hull_points $f_id]
set respt [concat $respt $deleted_points]
lset hull_points $f_id [list]
}
# Reassign points
foreach vec $respt {
if {[vector_equals $vec $point]} {
continue
}
foreach f_id $resfnew {
set facet [lindex $facets $f_id]
set plane [dict get $facet plane]
if {[is_above $vec $plane]} {
set current_points [lindex $hull_points $f_id]
lappend current_points $vec
lset hull_points $f_id $current_points
break
}
}
}
# Enqueue new faces
foreach f_id $resfnew {
lappend queue $f_id
}
}
dict set hull index $new_horizon
return $hull
}
# Main execution
proc main {} {
prepareConvexHulls
# Example: a tetrahedron
set n 4
set points [list]
# points[0] is placeholder
lappend points [vector_create 0 0 0 0]
lappend points [vector_create 0 0 0 1]
lappend points [vector_create 1 0 0 2]
lappend points [vector_create 0 1 0 3]
lappend points [vector_create 0 0 1 4]
set hull [QuickHull3D $points $n]
set surface_area [get_surface_area hull]
puts [format "%.3f" $surface_area]
}
# Run the main function
main
- Output:
2.366
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
const std = @import("std");
const print = std.debug.print;
const ArrayList = std.ArrayList;
const Allocator = std.mem.Allocator;
const math = std.math;
const MAXN: usize = 2500;
const EPS: f64 = 1e-8;
const Vect = struct {
x: f64,
y: f64,
z: f64,
id: usize,
fn init(x: f64, y: f64, z: f64, id: usize) Vect {
return Vect{ .x = x, .y = y, .z = z, .id = id };
}
fn sub(self: Vect, other: Vect) Vect {
return Vect.init(self.x - other.x, self.y - other.y, self.z - other.z, 0);
}
fn cross(self: Vect, other: Vect) Vect {
return Vect.init(
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: Vect, other: Vect) f64 {
return self.x * other.x + self.y * other.y + self.z * other.z;
}
fn magnitude(self: Vect) f64 {
return math.sqrt(self.x * self.x + self.y * self.y + self.z * self.z);
}
fn eq(self: Vect, other: Vect) bool {
return floatEq(self.x, other.x) and floatEq(self.y, other.y) and floatEq(self.z, other.z);
}
};
const Line = struct {
u: Vect,
v: Vect,
fn init(u: Vect, v: Vect) Line {
return Line{ .u = u, .v = v };
}
};
const Plane = struct {
vec: [3]Vect,
fn init(my_u: Vect, my_v: Vect, my_w: Vect) Plane {
return Plane{ .vec = [3]Vect{ my_u, my_v, my_w } };
}
fn normal(self: Plane) Vect {
return self.vec[1].sub(self.vec[0]).cross(self.vec[2].sub(self.vec[0]));
}
fn u(self: Plane) Vect {
return self.vec[0];
}
};
fn gtr(a: f64, b: f64) bool {
return a - b > EPS;
}
fn floatEq(a: f64, b: f64) bool {
return -EPS < a - b and a - b < EPS;
}
fn absFloat(x: f64) f64 {
return if (gtr(0.0, x)) -x else x;
}
// Signed distance
fn distPointPlane(v: Vect, p: Plane) f64 {
const norm = p.normal();
return v.sub(p.u()).dot(norm) / norm.magnitude();
}
// Unsigned distance
fn distPointLine(v: Vect, f: Line) f64 {
if (v.sub(f.u).magnitude() > 0.0) {
return f.v.sub(f.u).cross(v.sub(f.u)).magnitude() / f.v.sub(f.u).magnitude();
} else {
return 0.0;
}
}
fn distPointPoint(u: Vect, v: Vect) f64 {
return u.sub(v).magnitude();
}
fn isAbove(v: Vect, p: Plane) bool {
return gtr(v.sub(p.u()).dot(p.normal()), 0.0);
}
const Facet = struct {
n: [3]usize, // neighbors, correspond to point (u->v, v->w, w->u)
id: usize,
vistime: usize, // access timestamp
isdel: bool,
p: Plane,
fn init(id: usize, p: Plane) Facet {
return Facet{
.n = [3]usize{ 0, 0, 0 },
.id = id,
.vistime = 0,
.isdel = false,
.p = p,
};
}
fn setNeighbors(self: *Facet, n1: usize, n2: usize, n3: usize) void {
self.n = [3]usize{ n1, n2, n3 };
}
};
const Edge = struct {
netid: usize,
facetid: usize,
fn init() Edge {
return Edge{ .netid = 0, .facetid = 0 };
}
};
const ConvexHulls3d = struct {
index: usize, // Index face
surfacearea: f64,
fn init(index: usize) ConvexHulls3d {
return ConvexHulls3d{
.index = index,
.surfacearea = 0.0,
};
}
fn dfsArea(self: *ConvexHulls3d, nf: usize, fac: []Facet, time: *usize) void {
// Already visited in current timestamp
if (fac[nf].vistime == time.*) {
return;
}
fac[nf].vistime = time.*;
self.surfacearea += fac[nf].p.normal().magnitude() / 2.0;
// Make a copy of neighbors to avoid borrow issues
const neighbors = fac[nf].n;
for (neighbors) |neighbor| {
self.dfsArea(neighbor, fac, time);
}
}
fn getSurfaceArea(self: *ConvexHulls3d, fac: []Facet, time: *usize) f64 {
if (gtr(self.surfacearea, 0.0)) {
return self.surfacearea;
}
time.* += 1;
self.dfsArea(self.index, fac, time);
return self.surfacearea;
}
fn getHorizon(
self: ConvexHulls3d,
f: usize,
p: Vect,
vistime: []usize,
e1: []Edge,
e2: []Edge,
resfdel: *ArrayList(usize),
fac: []Facet,
time: usize,
) !i32 {
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;
try resfdel.append(fac[f].id);
var ret: i32 = -2;
// Make a copy of neighbors to avoid borrow issues
const neighbors = fac[f].n;
for (neighbors, 0..) |neighbor, i| {
const res = try self.getHorizon(
neighbor,
p,
vistime,
e1,
e2,
resfdel,
fac,
time,
);
if (res == 0) {
const pt = [2]usize{
fac[f].p.vec[i].id,
fac[f].p.vec[(i + 1) % 3].id,
};
for (pt, 0..) |point_id, j| {
if (vistime[point_id] != time) {
vistime[point_id] = time;
e1[point_id].netid = pt[(j + 1) % 2];
e1[point_id].facetid = neighbor;
} else {
e2[point_id].netid = pt[(j + 1) % 2];
e2[point_id].facetid = neighbor;
}
}
ret = @intCast(pt[0]);
} else if (res != -1 and res != -2) {
// The face is enclosed in the middle
ret = res;
}
}
return ret;
}
};
// Construct initial simplex
fn getStart(allocator: Allocator, point: []Vect, totp: usize) !struct { ConvexHulls3d, ArrayList(Facet), ArrayList(ArrayList(Vect)) } {
var fac = ArrayList(Facet).init(allocator);
var pts = ArrayList(ArrayList(Vect)).init(allocator);
// Add a dummy facet at index 0
try fac.append(Facet.init(0, Plane.init(
Vect.init(0.0, 0.0, 0.0, 0),
Vect.init(0.0, 0.0, 0.0, 0),
Vect.init(0.0, 0.0, 0.0, 0),
)));
var pt = [_]Vect{point[1]} ** 6;
var s = [_]Vect{point[1]} ** 4;
// Find the maximum point of the coordinate axis
for (2..totp + 1) |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];
}
}
// Take the two points with the largest distance
for (0..6) |i| {
for (i + 1..6) |j| {
if (gtr(distPointPoint(pt[i], pt[j]), distPointPoint(s[0], s[1]))) {
s[0] = pt[i];
s[1] = pt[j];
}
}
}
// Take the point farthest from the line connecting the two points
for (0..6) |i| {
if (gtr(
distPointLine(pt[i], Line.init(s[0], s[1])),
distPointLine(s[2], Line.init(s[0], s[1])),
)) {
s[2] = pt[i];
}
}
// Take the point farthest from the face
for (1..totp + 1) |i| {
if (gtr(
absFloat(distPointPlane(point[i], Plane.init(s[0], s[1], s[2]))),
absFloat(distPointPlane(s[3], Plane.init(s[0], s[1], s[2]))),
)) {
s[3] = point[i];
}
}
// Ensure that the constructed face faces outwards
if (gtr(0.0, distPointPlane(s[3], Plane.init(s[0], s[1], s[2])))) {
// Swap s[1] and s[2]
const temp = s[1];
s[1] = s[2];
s[2] = temp;
}
// Construct simplex
var f = [4]usize{ 0, 0, 0, 0 };
for (0..4) |i| {
try fac.append(Facet.init(
fac.items.len,
Plane.init(
Vect.init(0.0, 0.0, 0.0, 0),
Vect.init(0.0, 0.0, 0.0, 0),
Vect.init(0.0, 0.0, 0.0, 0),
),
));
f[i] = fac.items.len - 1;
}
fac.items[f[0]].p = Plane.init(s[0], s[2], s[1]); // Bottom face
fac.items[f[1]].p = Plane.init(s[0], s[1], s[3]);
fac.items[f[2]].p = Plane.init(s[1], s[2], s[3]);
fac.items[f[3]].p = Plane.init(s[2], s[0], s[3]);
fac.items[f[0]].setNeighbors(f[3], f[2], f[1]);
fac.items[f[1]].setNeighbors(f[0], f[2], f[3]);
fac.items[f[2]].setNeighbors(f[0], f[3], f[1]);
fac.items[f[3]].setNeighbors(f[0], f[1], f[2]);
// Assign point set space
for (0..5) |_| {
try pts.append(ArrayList(Vect).init(allocator));
}
// Assign points to four faces
for (1..totp + 1) |i| {
if (point[i].eq(s[0]) or point[i].eq(s[1]) or point[i].eq(s[2]) or point[i].eq(s[3])) {
continue;
}
for (0..4) |j| {
if (isAbove(point[i], fac.items[f[j]].p)) {
try pts.items[f[j]].append(point[i]);
break;
}
}
}
// Return the initial simplex, using a face as index
return .{ ConvexHulls3d.init(f[0]), fac, pts };
}
fn quickHull3d(allocator: Allocator, point: []Vect, totp: usize) !struct { ConvexHulls3d, ArrayList(Facet) } {
const result = try getStart(allocator, point, totp);
var hull = result[0];
var fac = result[1];
var pts = result[2];
defer {
for (pts.items) |*pt_list| {
pt_list.deinit();
}
pts.deinit();
}
// Add the face of initial simplex to queue
var que = ArrayList(usize).init(allocator);
defer que.deinit();
try que.append(hull.index);
for (fac.items[hull.index].n) |neighbor| {
try que.append(neighbor);
}
// snew saves index face of the final convex hull
var snew: usize = 0;
var time: usize = 0;
// Border line graph information
var e1 = try allocator.alloc(Edge, MAXN);
defer allocator.free(e1);
var e2 = try allocator.alloc(Edge, MAXN);
defer allocator.free(e2);
for (0..MAXN) |i| {
e1[i] = Edge.init();
e2[i] = Edge.init();
}
// Timestamp of each point access
var vistime = try allocator.alloc(usize, MAXN);
defer allocator.free(vistime);
@memset(vistime, 0);
while (que.items.len > 0) {
const nf = que.orderedRemove(0);
// Skip if the current face has been deleted
if (fac.items[nf].isdel) {
continue;
}
// Skip if no vertices are allocated to the current face
if (pts.items[nf].items.len == 0) {
snew = nf;
continue;
}
// Find the farthest point from the face
var p = pts.items[nf].items[0];
for (1..pts.items[nf].items.len) |i| {
if (gtr(
distPointPlane(pts.items[nf].items[i], fac.items[nf].p),
distPointPlane(p, fac.items[nf].p),
)) {
p = pts.items[nf].items[i];
}
}
// Find the horizon
time += 1;
var resfdel = ArrayList(usize).init(allocator);
defer resfdel.deinit();
// The current face must be deleted, so start dfs from current face
const s_result = try hull.getHorizon(
nf,
p,
vistime,
e1,
e2,
&resfdel,
fac.items,
time,
);
// Iterate over horizon(go around a circle), construct new face
var resfnew = ArrayList(usize).init(allocator);
defer resfnew.deinit();
time += 1;
var from: usize = 0; // The previous visited point
var lastf: usize = 0; // The last created face
var fstf: usize = 0; // The first created face
var s: usize = @intCast(s_result);
while (vistime[s] != time) {
// Record whether the current point has been visited with timestamp
vistime[s] = time;
var net: usize = 0; // Next point
var f: usize = 0; // The unseen face connected to the current edge on horizon
var fnew: usize = 0; // 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
var pt1: i32 = -1;
var pt2: i32 = -1;
for (0..3) |i| {
if (point[s].eq(fac.items[f].p.vec[i])) {
pt1 = @intCast(i);
}
if (point[net].eq(fac.items[f].p.vec[i])) {
pt2 = @intCast(i);
}
}
// Make sure pt1->pt2 is arranged counterclockwise by adjacent face points
if ((@rem(pt1 + 1, 3)) != pt2) {
const temp = pt1;
pt1 = pt2;
pt2 = temp;
}
// The face constructed in this way faces outward
try fac.append(Facet.init(
fac.items.len,
Plane.init(
fac.items[f].p.vec[@intCast(pt2)],
fac.items[f].p.vec[@intCast(pt1)],
p,
),
));
fnew = fac.items.len - 1;
try pts.append(ArrayList(Vect).init(allocator));
try resfnew.append(fnew);
// Maintain adjacency information
fac.items[fnew].n[0] = f;
fac.items[f].n[@intCast(pt1)] = fnew;
if (lastf != 0) {
// Can't determine whether to traverse clockwise or counterclockwise in advance
// Maintain adjacency information between new faces
if (fac.items[fnew].p.vec[1].eq(fac.items[lastf].p.vec[0])) {
fac.items[fnew].n[1] = lastf;
fac.items[lastf].n[2] = fnew;
} else {
fac.items[fnew].n[2] = lastf;
fac.items[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.items[fstf].p.vec[1].eq(fac.items[lastf].p.vec[0])) {
fac.items[fstf].n[1] = lastf;
fac.items[lastf].n[2] = fstf;
} else {
fac.items[fstf].n[2] = lastf;
fac.items[lastf].n[1] = fstf;
}
// Get all the points to be assigned
var respt = ArrayList(Vect).init(allocator);
defer respt.deinit();
for (resfdel.items) |deleted_face| {
for (pts.items[deleted_face].items) |pt| {
try respt.append(pt);
}
pts.items[deleted_face].clearAndFree();
}
// Assign points
for (respt.items) |pt| {
if (pt.eq(p)) {
continue; // Skip the points used to create the new face
}
for (resfnew.items) |new_face| {
if (isAbove(pt, fac.items[new_face].p)) {
try pts.items[new_face].append(pt);
break; // Make sure the points are not reassigned
}
}
}
// Add the new face to queue
for (resfnew.items) |new_face| {
try que.append(new_face);
}
}
hull.index = snew;
return .{ hull, fac };
}
pub fn main() !void {
var gpa = std.heap.GeneralPurposeAllocator(.{}){};
defer _ = gpa.deinit();
const allocator = gpa.allocator();
const n = 4; // number of points
var point = try allocator.alloc(Vect, n + 1); // 0th element is a placeholder
defer allocator.free(point);
point[0] = Vect.init(0.0, 0.0, 0.0, 0); // placeholder
const my_input = [_][3]f64{
[_]f64{ 0.0, 0.0, 0.0 },
[_]f64{ 1.0, 0.0, 0.0 },
[_]f64{ 0.0, 1.0, 0.0 },
[_]f64{ 0.0, 0.0, 1.0 },
};
for (1..n + 1) |i| {
const x = my_input[i - 1][0];
const y = my_input[i - 1][1];
const z = my_input[i - 1][2];
point[i] = Vect.init(x, y, z, i);
}
const result = try quickHull3d(allocator, point, n);
var hull = result[0];
var fac = result[1];
defer fac.deinit();
var time: usize = 0;
print("{d:.3}\n", .{hull.getSurfaceArea(fac.items, &time)});
}
- Output:
2.366