Skip to content

Instantly share code, notes, and snippets.

@YclepticStudios
Created October 23, 2016 20:55
Show Gist options
  • Save YclepticStudios/c0d8cea56ee1e1d9714754ee9427085b to your computer and use it in GitHub Desktop.
Save YclepticStudios/c0d8cea56ee1e1d9714754ee9427085b to your computer and use it in GitHub Desktop.
Convex hulling algorithm for Unity; converts a point cloud into triangle indices for meshing.
/**
* ============================================================================
* MIT License
*
* Copyright (c) 2016 Eric Phillips
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
* ============================================================================
*
*
* This file implements a 3D generalization of the convex hulling algorithm
* known as the "QuickHull" algorithm. This algorithm generates a surface for a
* point cloud which contains all of the points and all of the surface features
* possible using only a convex mesh.
*
* The procedure followed here is based on the following example.
* http://thomasdiewald.com/blog/?p=1888
*
* Created by Eric Phillips on October 23, 2016.
*/
using System.Collections.Generic;
using UnityEngine;
public static class ConvexHull
{
/// <summary>
/// Class for holding faces and their data.
/// </summary>
private class FaceData
{
public Vector4 Plane; // Equation of a face's plane
public int[] FaceIndices; // Indices of the face's points
public List<int> VisibleIndices; // Indices of points visible to face
private Vector3[] _vertices; // Main vertex array
/// <summary>
/// Create a new data structure for a face.
/// </summary>
/// <param name="pt1Idx">Index of the first point on the face.</param>
/// <param name="pt2Idx">Index of the second point on the face.</param>
/// <param name="pt3Idx">Index of the third point on the face.</param>
/// <param name="vertices">Array of vertices.</param>
public FaceData(int pt1Idx, int pt2Idx, int pt3Idx, Vector3[] vertices)
{
Plane = Points2Plane(vertices[pt1Idx], vertices[pt2Idx],
vertices[pt3Idx]);
FaceIndices = new int[] { pt1Idx, pt2Idx, pt3Idx };
VisibleIndices = new List<int>();
_vertices = vertices;
}
/// <summary>
/// Get the point farthest from the plane in the positive direction.
/// </summary>
/// <returns>The index of the point.</returns>
public int GetFurthestPoint()
{
int maxIndex = 0;
float maxDistance = float.NegativeInfinity;
foreach (int index in VisibleIndices)
{
// Calculate something linearly equivalent to the plane distance
Vector3 v = _vertices[index] - _vertices[FaceIndices[0]];
float distance = Vector3.Dot(v, Plane);
if (distance > maxDistance)
{
maxIndex = index;
maxDistance = distance;
}
}
return maxIndex;
}
}
/// <summary>
/// Compute indices of triangles which form a convex hull around the given
/// array of vertices. This is formatted so that a set of mesh vertices can
/// be passed in, and the output can be directly assigned to the triangles
/// property of the mesh.
/// </summary>
/// <param name="vertices">Array of vertices passed in.</param>
/// <returns>Array of triangles as indices into the input array.</returns>
public static int[] Generate(Vector3[] vertices)
{
// Create initial simplex
int[] extremePoints = FindExtremePoints(vertices);
int[] initialSimplex = CreateInitialSimplex(extremePoints, vertices);
// Assign each vertex to the first face to which it is visible
List<FaceData> faceStack =
AssignInitialPointsToFaces(initialSimplex, vertices);
// Iterate until all faces have been processed
FaceData face;
while ((face = GetUnprocessedFace(faceStack)) != null)
{
// Get the point farthest from the plane in the positive direction
int maxPtIndex = face.GetFurthestPoint();
// Get all faces which are visible to this point and pop from stack
List<FaceData> visibleFaces = ExtractVisibleFaces(
vertices[maxPtIndex], faceStack);
// Extract the horizon edges of the visible faces
// All faces will be connected
Dictionary<string, int[]> horizonEdges =
ExtractHorizonEdges(visibleFaces);
// Create new points from horizon edges and max point
List<FaceData> newFaces = CreateNewFaces(maxPtIndex, visibleFaces,
vertices, horizonEdges);
// Add new faces to stack and repeat
faceStack.AddRange(newFaces);
}
// Compile one array of all the triangle indices
int[] indices = new int[faceStack.Count * 3];
for (int ii = 0; ii < faceStack.Count; ii++)
for (int jj = 0; jj < 3; jj++)
indices[ii * 3 + jj] = faceStack[ii].FaceIndices[jj];
return indices;
}
/// <summary>
/// Find the indices of the points in the vertices array which are at a
/// minimum or maximum on each of the three axis. So find the min and max
/// X, the min and max Y and the min and max Z points, and return their
/// indices.
/// </summary>
/// <param name="vertices">Array of vertices passed in.</param>
/// <returns>Array of six indices into the vertices array.</returns>
private static int[] FindExtremePoints(Vector3[] vertices)
{
// Setup variables
int[] extremePoints = new int[6];
float[] extremePointValues = new float[6];
for (int ii = 0; ii < extremePoints.Length - 1; ii += 2)
{
extremePoints[ii] = 0;
extremePoints[ii + 1] = 0;
extremePointValues[ii] = float.PositiveInfinity;
extremePointValues[ii + 1] = float.NegativeInfinity;
}
// Search point cloud
for (int ii = 0; ii < vertices.Length; ii++)
for (int jj = 0; jj < extremePoints.Length - 1; jj += 2)
{
float val = vertices[ii][jj / 2];
if (val < extremePointValues[jj])
{
extremePoints[jj] = ii;
extremePointValues[jj] = val;
}
else if (val > extremePointValues[jj + 1])
{
extremePoints[jj + 1] = ii;
extremePointValues[jj + 1] = val;
}
}
return extremePoints;
}
/// <summary>
/// Create a tetrahedron of maximum volume out of the extreme points.
/// This returns four indices in right-handed manner. The first three
/// define the base of the tetrahedron, and face away from the third point.
/// </summary>
/// <param name="extremePoints">The indices of the extreme points.</param>
/// <param name="vertices">Array of vertices passed in.</param>
/// <returns>The indices of the initial tetrahedron.</returns>
private static int[] CreateInitialSimplex(int[] extremePoints,
Vector3[] vertices)
{
int[] initialSimplex = new int[4];
// Find two most distent extreme points (base line of tetrahedron)
float maxDistance = float.NegativeInfinity;
for (int ii = 0; ii < extremePoints.Length; ii++)
for (int jj = ii + 1; jj < extremePoints.Length; jj++)
{
float distance = (vertices[extremePoints[ii]] -
vertices[extremePoints[jj]]).sqrMagnitude;
if (distance > maxDistance)
{
initialSimplex[0] = extremePoints[ii];
initialSimplex[1] = extremePoints[jj];
maxDistance = distance;
}
}
// Find the extreme point most distent from the line
maxDistance = float.NegativeInfinity;
Vector3 normal = vertices[initialSimplex[0]] -
vertices[initialSimplex[1]];
for (int ii = 0; ii < extremePoints.Length; ii++)
{
Vector3 v = vertices[extremePoints[ii]] -
vertices[initialSimplex[0]];
Vector3 rejection = Vector3.ProjectOnPlane(v, normal);
float distance = rejection.sqrMagnitude;
if (distance > maxDistance)
{
initialSimplex[2] = extremePoints[ii];
maxDistance = distance;
}
}
// Find the most distant of all the points from the plane of the
// triangle formed from the first three "initialSimplex" points
maxDistance = float.NegativeInfinity;
Vector3 v1 = vertices[initialSimplex[1]] - vertices[initialSimplex[0]];
Vector3 v2 = vertices[initialSimplex[2]] - vertices[initialSimplex[0]];
normal = Vector3.Cross(v1, v2);
for (int ii = 0; ii < vertices.Length; ii++)
{
Vector3 v = vertices[ii] - vertices[initialSimplex[0]];
float distance = Mathf.Abs(Vector3.Dot(v, normal));
if (distance > maxDistance)
{
initialSimplex[3] = ii;
maxDistance = distance;
}
}
// Swap the two first vertices if the final point is in front of the
// triangular base plane (this makes all faces of the tetrahedron point)
// outward
Vector4 baseFace = Points2Plane(vertices[initialSimplex[0]],
vertices[initialSimplex[1]],
vertices[initialSimplex[2]]);
if (PointAbovePlane(vertices[initialSimplex[3]], baseFace))
{
int t = initialSimplex[0];
initialSimplex[0] = initialSimplex[1];
initialSimplex[1] = t;
}
return initialSimplex;
}
/// <summary>
/// Create the first four faces and assign all the points to the first face
/// to which they are visible.
/// </summary>
/// <param name="initialSimplex">The indices of the initial simplex.</param>
/// <param name="vertices">Array of vertices.</param>
/// <returns>The initial list of faces with their data.</returns>
private static List<FaceData> AssignInitialPointsToFaces(
int[] initialSimplex, Vector3[] vertices)
{
// Create the first four faces
List<FaceData> faceData = new List<FaceData>();
faceData.Add(new FaceData(initialSimplex[0], initialSimplex[1],
initialSimplex[2], vertices));
faceData.Add(new FaceData(initialSimplex[0], initialSimplex[2],
initialSimplex[3], vertices));
faceData.Add(new FaceData(initialSimplex[1], initialSimplex[3],
initialSimplex[2], vertices));
faceData.Add(new FaceData(initialSimplex[1], initialSimplex[0],
initialSimplex[3], vertices));
// Assign each point to the first face to which it is visible
for (int ii = 0; ii < vertices.Length; ii++)
foreach (FaceData face in faceData)
if (PointAbovePlane(vertices[ii], face.Plane))
{
face.VisibleIndices.Add(ii);
break;
}
return faceData;
}
/// <summary>
/// Get a face from the stack which still needs to be processed.
/// </summary>
/// <param name="faceStack">The list of faces.</param>
/// <returns>The face to process next.</returns>
private static FaceData GetUnprocessedFace(List<FaceData> faceStack)
{
foreach (FaceData face in faceStack)
if (face.VisibleIndices.Count > 0)
return face;
return null;
}
/// <summary>
/// Find all the faces which can be seen from the given point. This is the
/// same as if the point were a light and we were looking for all faces
/// which the point illuminated. Faces pointing the opposite direction are
/// ignored. These faces are removed from the stack.
/// </summary>
/// <param name="pt">The point in question.</param>
/// <param name="faceStack">The list of faces to search.</param>
/// <returns>A new list of faces.</returns>
private static List<FaceData> ExtractVisibleFaces(Vector3 pt,
List<FaceData> faceStack)
{
List<FaceData> visibleFaces = new List<FaceData>();
foreach (FaceData face in faceStack)
if (PointAbovePlane(pt, face.Plane))
visibleFaces.Add(face);
foreach (FaceData face in visibleFaces)
faceStack.Remove(face);
return visibleFaces;
}
/// <summary>
/// Get pairings of all the edges around the set of visible faces.
/// This basically finds one outline around all the given faces.
/// This is returned as a dictionary to reduce data copying between
/// functions.
/// </summary>
/// <param name="visibleFaces">A list of the visible faces.</param>
/// <returns>A dictionary with the edge pairs.</returns>
private static Dictionary<string, int[]> ExtractHorizonEdges(
List<FaceData> visibleFaces)
{
Dictionary<string, int[]> edges = new Dictionary<string, int[]>();
foreach (FaceData face in visibleFaces)
for (int ii = 0; ii < face.FaceIndices.Length; ii++)
{
// Get indices of triangle vertices
int idx1 = face.FaceIndices[ii];
int idx2 = face.FaceIndices[(ii + 1) % face.FaceIndices.Length];
// Use keys to easily detect and remove duplicate edges
// If an edge appears twice, it is shared by two triangles,
// and not really an edge
string key = idx1 + "," + idx2;
string keyReversed = idx2 + "," + idx1;
if (edges.ContainsKey(key) || edges.ContainsKey(keyReversed))
{
edges.Remove(key);
edges.Remove(keyReversed);
}
else
edges.Add(key, new int[] { idx1, idx2 });
}
return edges;
}
/// <summary>
/// Create new faces between the horizon edges and the maximum point.
/// This is somewhat like extruding the edges out to the maximum point to
/// form new faces.
/// </summary>
/// <param name="maxPtIndex">Point farthest from current face.</param>
/// <param name="visibleFaces">List of current visible faces.</param>
/// <param name="vertices">Array of vertices.</param>
/// <param name="horizonEdges">Dictionary of pairs of indices.</param>
/// <returns>A new list of faces.</returns>
private static List<FaceData> CreateNewFaces(int maxPtIndex,
List<FaceData> visibleFaces, Vector3[] vertices,
Dictionary<string, int[]> horizonEdges)
{
List<FaceData> newFaces = new List<FaceData>();
// Create new faces from horizon edges and max point
foreach (int[] edge in horizonEdges.Values)
newFaces.Add(new FaceData(edge[0], edge[1], maxPtIndex, vertices));
// Assign all points off all visible faces to the first of the new
// faces to which the point is visible
foreach (FaceData face in visibleFaces)
foreach (int idx in face.VisibleIndices)
foreach (FaceData newFace in newFaces)
if (PointAbovePlane(vertices[idx], newFace.Plane))
{
newFace.VisibleIndices.Add(idx);
break;
}
return newFaces;
}
/// <summary>
/// Convert three points to a plane equation of the form
/// Ax + By + Cz + D = 0. The returned Vector4 contains A, B, C and D.
/// </summary>
/// <param name="pt1">The first point.</param>
/// <param name="pt2">The second point.</param>
/// <param name="pt3">The third point.</param>
/// <returns>A new Vector4.</returns>
private static Vector4 Points2Plane(Vector3 pt1, Vector3 pt2, Vector3 pt3)
{
Vector4 v = Vector3.Cross(pt2 - pt1, pt3 - pt1).normalized;
v.w = -Vector3.Dot(v, pt1);
return v;
}
/// <summary>
/// Checks if the point is on the front face of the plane.
/// This is the same as the point being able to "see" the front of the
/// plane.
/// </summary>
/// <param name="point">The point in question.</param>
/// <param name="plane">The plane in question.</param>
/// <returns>A boolean value.</returns>
private static bool PointAbovePlane(Vector3 point, Vector4 plane)
{
return Vector3.Dot(point, plane) + plane.w > 0.0000001f;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment