From 39acb393201e9c0d29454509fa8a6ca2177f41db Mon Sep 17 00:00:00 2001
From: John Collins <collins@braincorporation.com>
Date: Wed, 29 Nov 2023 21:12:56 -0800
Subject: [PATCH] Delete quickhull, use builtin

---
 bindings/java/pom.xml                         |   11 +-
 .../src/main/cpp/manifold3d/QuickHull.hpp     |  995 --------
 .../main/cpp/manifold3d/QuickHullTomilov.hpp  | 1029 --------
 .../src/main/cpp/manifold3d/convex_hull.hpp   |  754 ------
 .../cpp/manifold3d/convex_hull_qhull_impl.hpp |  216 --
 bindings/java/src/main/cpp/manifold3d/vmath.h | 2201 -----------------
 .../src/main/java/manifold3d/ConvexHull.java  |   33 -
 .../main/java/manifold3d/LibraryPaths.java    |    2 +-
 .../src/main/java/manifold3d/Manifold.java    |   18 +-
 .../manifold3d/manifold/CrossSection.java     |   32 +-
 .../manifold/CrossSectionVector.java          |    6 +-
 11 files changed, 41 insertions(+), 5256 deletions(-)
 delete mode 100644 bindings/java/src/main/cpp/manifold3d/QuickHull.hpp
 delete mode 100644 bindings/java/src/main/cpp/manifold3d/QuickHullTomilov.hpp
 delete mode 100644 bindings/java/src/main/cpp/manifold3d/convex_hull.hpp
 delete mode 100644 bindings/java/src/main/cpp/manifold3d/convex_hull_qhull_impl.hpp
 delete mode 100644 bindings/java/src/main/cpp/manifold3d/vmath.h
 delete mode 100644 bindings/java/src/main/java/manifold3d/ConvexHull.java

diff --git a/bindings/java/pom.xml b/bindings/java/pom.xml
index db9b9009f..33b70833f 100644
--- a/bindings/java/pom.xml
+++ b/bindings/java/pom.xml
@@ -8,8 +8,8 @@
 
     <properties>
         <project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
-        <maven.compiler.source>17</maven.compiler.source>
-        <maven.compiler.target>17</maven.compiler.target>
+        <maven.compiler.source>11</maven.compiler.source>
+        <maven.compiler.target>11</maven.compiler.target>
     </properties>
 
     <dependencies>
@@ -60,6 +60,13 @@
                     <include>libmanifold.dylib</include>
                 </includes>
             </resource>
+            <resource>
+                <directory>../../build/src/third_party</directory>
+                <includes>
+                    <include>libquickhull.so</include>
+                    <include>libquickhull.dylib</include>
+                </includes>
+            </resource>
             <resource>
                 <directory>../../build/src/cross_section</directory>
                 <includes>
diff --git a/bindings/java/src/main/cpp/manifold3d/QuickHull.hpp b/bindings/java/src/main/cpp/manifold3d/QuickHull.hpp
deleted file mode 100644
index 07dbedf35..000000000
--- a/bindings/java/src/main/cpp/manifold3d/QuickHull.hpp
+++ /dev/null
@@ -1,995 +0,0 @@
-/*
- * QuickHull.hpp
- *
- *  Created on: May 4, 2014
- *      Author: anttiku
- *
- *  Single header version of all header files from the QuickHull implementation
- *  at https://github.com/akuukka/quickhull by Antti Kuukka
- *
- *  Per the README.md file:
- *
- *  This implementation is 100% Public Domain.
- *
- *  Feel free to use.
- *
- *  C++11 is needed to compile it.
- */
-
-#ifndef QUICKHULL_HPP_
-#define QUICKHULL_HPP_
-
-#include "common.h"
-
-#include <cmath>
-#include <cstdint>
-#include <deque>
-#include <iostream>
-
-#include <array>
-#include <cassert>
-#include <fstream>
-#include <limits>
-#include <memory>
-#include <string>
-#include <unordered_map>
-#include <vector>
-
-extern "C" {
-#include "vmath.h"
-}
-
-/*
- * Implementation of the 3D QuickHull algorithm by Antti Kuukka
- *
- * No copyrights. What follows is 100% Public Domain.
- *
- *
- *
- * INPUT:  a list of points in 3D space (for example, vertices of a 3D mesh)
- *
- * OUTPUT: a ConvexHull object which provides vertex and index buffers of the
- * generated convex hull as a triangle mesh.
- *
- *
- *
- * The implementation is thread-safe if each thread is using its own QuickHull
- * object.
- *
- *
- * SUMMARY OF THE ALGORITHM:
- *         - Create initial simplex (tetrahedron) using extreme points. We
- *           have four faces now and they form a convex mesh M.
- *         - For each point, assign them to the first face for which they
- *           are on the positive side of (so each point is assigned to at most
- *           one face). Points inside the initial tetrahedron are left behind
- *           now and no longer affect the calculations.
- *         - Add all faces that have points assigned to them to Face Stack.
- *         - Iterate until Face Stack is empty:
- *              - Pop topmost face F from the stack
- *              - From the points assigned to F, pick the point P that is
- *                farthest away from the plane defined by F.
- *              - Find all faces of M that have P on their positive side. Let
- *                us call these the "visible faces".
- *              - Because of the way M is constructed, these faces are
- *                connected. Solve their horizon edge loop.
- *		- "Extrude to P": Create new faces by connecting P with the
- *		   points belonging to the horizon edge. Add the new faces to
- *		   M and remove the visible faces from M.
- *              - Each point that was assigned to visible faces is now assigned
- *                to at most one of the newly created faces.
- *              - Those new faces that have points assigned to them are added
- *                to the top of Face Stack.
- *          - M is now the convex hull.
- *
- * TO DO:
- *  - Implement a proper 2D QuickHull and use that to solve the degenerate 2D
- *    case (when all the points lie on the same plane in 3D space).
- */
-
-
-namespace quickhull {
-
-    typedef size_t size_t;
-
-    template <typename T>
-	class Vector3
-	{
-	    public:
-		Vector3() = default;
-
-		Vector3(T px, T py, T pz) : x(px), y(py), z(pz) {
-
-		}
-
-		T x,y,z;
-
-		T dotProduct(const Vector3& other) const {
-		    return x*other.x+y*other.y+z*other.z;
-		}
-
-		void normalize() {
-		    const T len = getLength();
-		    x/=len;
-		    y/=len;
-		    z/=len;
-		}
-
-		Vector3 getNormalized() const {
-		    const T len = getLength();
-		    return Vector3(x/len,y/len,z/len);
-		}
-
-		T getLength() const {
-		    return std::sqrt(x*x+y*y+z*z);
-		}
-
-		Vector3 operator-(const Vector3& other) const {
-		    return Vector3(x-other.x,y-other.y,z-other.z);
-		}
-
-		Vector3 operator+(const Vector3& other) const {
-		    return Vector3(x+other.x,y+other.y,z+other.z);
-		}
-
-		Vector3& operator+=(const Vector3& other) {
-		    x+=other.x;
-		    y+=other.y;
-		    z+=other.z;
-		    return *this;
-		}
-		Vector3& operator-=(const Vector3& other) {
-		    x-=other.x;
-		    y-=other.y;
-		    z-=other.z;
-		    return *this;
-		}
-		Vector3& operator*=(T c) {
-		    x*=c;
-		    y*=c;
-		    z*=c;
-		    return *this;
-		}
-
-		Vector3& operator/=(T c) {
-		    x/=c;
-		    y/=c;
-		    z/=c;
-		    return *this;
-		}
-
-		Vector3 operator-() const {
-		    return Vector3(-x,-y,-z);
-		}
-
-		template<typename S>
-		    Vector3 operator*(S c) const {
-			return Vector3(x*c,y*c,z*c);
-		    }
-
-		template<typename S>
-		    Vector3 operator/(S c) const {
-			return Vector3(x/c,y/c,z/c);
-		    }
-
-		T getLengthSquared() const {
-		    return x*x + y*y + z*z;
-		}
-
-		bool operator!=(const Vector3& o) const {
-		    return (!EQUAL(x, o.x) || !EQUAL(y, o.y) || !EQUAL(z, o.z));
-		}
-
-		// Projection onto another vector
-		Vector3 projection(const Vector3& o) const {
-		    T C = dotProduct(o)/o.getLengthSquared();
-		    return o*C;
-		}
-
-		Vector3 crossProduct (const Vector3& rhs ) {
-		    T a = y * rhs.z - z * rhs.y ;
-		    T b = z * rhs.x - x * rhs.z ;
-		    T c = x * rhs.y - y * rhs.x ;
-		    Vector3 product( a , b , c ) ;
-		    return product ;
-		}
-
-		T getDistanceTo(const Vector3& other) const {
-		    Vector3 diff = *this - other;
-		    return diff.getLength();
-		}
-
-		T getSquaredDistanceTo(const Vector3& other) const {
-		    const T dx = x-other.x;
-		    const T dy = y-other.y;
-		    const T dz = z-other.z;
-		    return dx*dx+dy*dy+dz*dz;
-		}
-
-	};
-
-    // Overload also << operator for easy printing of debug data
-    template <typename T>
-	inline std::ostream& operator<<(std::ostream& os, const Vector3<T>& vec) {
-	    os << "(" << vec.x << "," << vec.y << "," << vec.z << ")";
-	    return os;
-	}
-
-    template <typename T>
-	inline Vector3<T> operator*(T c, const Vector3<T>& v) {
-	    return Vector3<T>(v.x*c,v.y*c,v.z*c);
-	}
-
-
-    template <typename T>
-	struct Ray {
-	    const Vector3<T> m_S;
-	    const Vector3<T> m_V;
-	    const T m_VInvLengthSquared;
-
-	    Ray(const Vector3<T>& S,const Vector3<T>& V) : m_S(S), m_V(V), m_VInvLengthSquared(1/m_V.getLengthSquared()) {
-	    }
-	};
-
-    template<typename T>
-	class Pool {
-	    std::vector<std::shared_ptr<T>> m_data;
-	    public:
-	    void clear() {
-		m_data.clear();
-	    }
-
-	    void reclaim(std::shared_ptr<T>& ptr) {
-		m_data.push_back(std::move(ptr));
-	    }
-
-	    std::shared_ptr<T> get() {
-		if (m_data.size()==0) {
-		    return std::shared_ptr<T>(new T());
-		}
-		auto it = m_data.end()-1;
-		std::shared_ptr<T> r = std::move(*it);
-		m_data.erase(it);
-		return r;
-	    }
-
-	};
-
-    template<typename T>
-	class Plane {
-	    public:
-		Vector3<T> m_N;
-
-		// Signed distance (if normal is of length 1) to the plane from origin
-		T m_D;
-
-		// Normal length squared
-		T m_sqrNLength;
-
-		bool isPointOnPositiveSide(const Vector3<T>& Q) const {
-		    T d = m_N.dotProduct(Q)+m_D;
-		    if (d>=0) return true;
-		    return false;
-		}
-
-		Plane() = default;
-
-		// Construct a plane using normal N and any point P on the plane
-		Plane(const Vector3<T>& N, const Vector3<T>& P) : m_N(N), m_D(-N.dotProduct(P)), m_sqrNLength(m_N.x*m_N.x+m_N.y*m_N.y+m_N.z*m_N.z) {
-
-		}
-	};
-
-
-    template<typename T>
-	class VertexDataSource {
-	    const Vector3<T>* m_ptr;
-	    size_t m_count;
-
-	    public:
-	    VertexDataSource(const Vector3<T>* ptr, size_t count) : m_ptr(ptr), m_count(count) {
-
-	    }
-
-	    VertexDataSource(const std::vector<Vector3<T>>& vec) : m_ptr(&vec[0]), m_count(vec.size()) {
-
-	    }
-
-	    VertexDataSource() : m_ptr(nullptr), m_count(0) {
-
-	    }
-
-	    VertexDataSource(const VertexDataSource &) = default;
-
-	    VertexDataSource& operator=(const VertexDataSource& other) = default;
-
-	    size_t size() const {
-		return m_count;
-	    }
-
-	    const Vector3<T>& operator[](size_t index) const {
-		return m_ptr[index];
-	    }
-
-	    const Vector3<T>* begin() const {
-		return m_ptr;
-	    }
-
-	    const Vector3<T>* end() const {
-		return m_ptr + m_count;
-	    }
-	};
-
-    template <typename T>
-	class MeshBuilder {
-	    public:
-		struct HalfEdge {
-		    size_t m_endVertex;
-		    size_t m_opp;
-		    size_t m_face;
-		    size_t m_next;
-
-		    void disable() {
-			m_endVertex = std::numeric_limits<size_t>::max();
-		    }
-
-		    bool isDisabled() const {
-			return m_endVertex == std::numeric_limits<size_t>::max();
-		    }
-		};
-
-		struct Face {
-		    size_t m_he;
-		    Plane<T> m_P;
-		    T m_mostDistantPointDist;
-		    size_t m_mostDistantPoint;
-		    size_t m_visibilityCheckedOnIteration;
-		    std::uint8_t m_isVisibleFaceOnCurrentIteration : 1;
-		    std::uint8_t m_inFaceStack : 1;
-		    std::uint8_t m_horizonEdgesOnCurrentIteration : 3; // Bit for each half edge assigned to this face, each being 0 or 1 depending on whether the edge belongs to horizon edge
-		    std::shared_ptr<std::vector<size_t>> m_pointsOnPositiveSide;
-
-		    Face() : m_he(std::numeric_limits<size_t>::max()),
-		    m_mostDistantPointDist(0),
-		    m_mostDistantPoint(0),
-		    m_visibilityCheckedOnIteration(0),
-		    m_isVisibleFaceOnCurrentIteration(0),
-		    m_inFaceStack(0),
-		    m_horizonEdgesOnCurrentIteration(0)
-		    {
-
-		    }
-
-		    void disable() {
-			m_he = std::numeric_limits<size_t>::max();
-		    }
-
-		    bool isDisabled() const {
-			return m_he == std::numeric_limits<size_t>::max();
-		    }
-		};
-
-		// Mesh data
-		std::vector<Face> m_faces;
-		std::vector<HalfEdge> m_halfEdges;
-
-		// When the mesh is modified and faces and half edges are removed from it, we do not actually remove them from the container vectors.
-		// Insted, they are marked as disabled which means that the indices can be reused when we need to add new faces and half edges to the mesh.
-		// We store the free indices in the following vectors.
-		std::vector<size_t> m_disabledFaces,m_disabledHalfEdges;
-
-		size_t addFace() {
-		    if (m_disabledFaces.size()) {
-			size_t index = m_disabledFaces.back();
-			auto& f = m_faces[index];
-			assert(f.isDisabled());
-			assert(!f.m_pointsOnPositiveSide);
-			f.m_mostDistantPointDist = 0;
-			m_disabledFaces.pop_back();
-			return index;
-		    }
-		    m_faces.emplace_back();
-		    return m_faces.size()-1;
-		}
-
-		size_t addHalfEdge()	{
-		    if (m_disabledHalfEdges.size()) {
-			const size_t index = m_disabledHalfEdges.back();
-			m_disabledHalfEdges.pop_back();
-			return index;
-		    }
-		    m_halfEdges.emplace_back();
-		    return m_halfEdges.size()-1;
-		}
-
-		// Mark a face as disabled and return a pointer to the points that were on the positive of it.
-		std::shared_ptr<std::vector<size_t>> disableFace(size_t faceIndex) {
-		    auto& f = m_faces[faceIndex];
-		    f.disable();
-		    m_disabledFaces.push_back(faceIndex);
-		    return std::move(f.m_pointsOnPositiveSide);
-		}
-
-		void disableHalfEdge(size_t heIndex) {
-		    auto& he = m_halfEdges[heIndex];
-		    he.disable();
-		    m_disabledHalfEdges.push_back(heIndex);
-		}
-
-		MeshBuilder() = default;
-
-		// Create a mesh with initial tetrahedron ABCD. Dot product of AB with the normal of triangle ABC should be negative.
-		void setup(size_t a, size_t b, size_t c, size_t d) {
-		    // Create halfedges
-		    HalfEdge AB;
-		    AB.m_endVertex = b;
-		    AB.m_opp = 6;
-		    AB.m_face = 0;
-		    AB.m_next = 1;
-		    m_halfEdges.push_back(AB);
-
-		    HalfEdge BC;
-		    BC.m_endVertex = c;
-		    BC.m_opp = 9;
-		    BC.m_face = 0;
-		    BC.m_next = 2;
-		    m_halfEdges.push_back(BC);
-
-		    HalfEdge CA;
-		    CA.m_endVertex = a;
-		    CA.m_opp = 3;
-		    CA.m_face = 0;
-		    CA.m_next = 0;
-		    m_halfEdges.push_back(CA);
-
-		    HalfEdge AC;
-		    AC.m_endVertex = c;
-		    AC.m_opp = 2;
-		    AC.m_face = 1;
-		    AC.m_next = 4;
-		    m_halfEdges.push_back(AC);
-
-		    HalfEdge CD;
-		    CD.m_endVertex = d;
-		    CD.m_opp = 11;
-		    CD.m_face = 1;
-		    CD.m_next = 5;
-		    m_halfEdges.push_back(CD);
-
-		    HalfEdge DA;
-		    DA.m_endVertex = a;
-		    DA.m_opp = 7;
-		    DA.m_face = 1;
-		    DA.m_next = 3;
-		    m_halfEdges.push_back(DA);
-
-		    HalfEdge BA;
-		    BA.m_endVertex = a;
-		    BA.m_opp = 0;
-		    BA.m_face = 2;
-		    BA.m_next = 7;
-		    m_halfEdges.push_back(BA);
-
-		    HalfEdge AD;
-		    AD.m_endVertex = d;
-		    AD.m_opp = 5;
-		    AD.m_face = 2;
-		    AD.m_next = 8;
-		    m_halfEdges.push_back(AD);
-
-		    HalfEdge DB;
-		    DB.m_endVertex = b;
-		    DB.m_opp = 10;
-		    DB.m_face = 2;
-		    DB.m_next = 6;
-		    m_halfEdges.push_back(DB);
-
-		    HalfEdge CB;
-		    CB.m_endVertex = b;
-		    CB.m_opp = 1;
-		    CB.m_face = 3;
-		    CB.m_next = 10;
-		    m_halfEdges.push_back(CB);
-
-		    HalfEdge BD;
-		    BD.m_endVertex = d;
-		    BD.m_opp = 8;
-		    BD.m_face = 3;
-		    BD.m_next = 11;
-		    m_halfEdges.push_back(BD);
-
-		    HalfEdge DC;
-		    DC.m_endVertex = c;
-		    DC.m_opp = 4;
-		    DC.m_face = 3;
-		    DC.m_next = 9;
-		    m_halfEdges.push_back(DC);
-
-		    // Create faces
-		    Face ABC;
-		    ABC.m_he = 0;
-		    m_faces.push_back(std::move(ABC));
-
-		    Face ACD;
-		    ACD.m_he = 3;
-		    m_faces.push_back(std::move(ACD));
-
-		    Face BAD;
-		    BAD.m_he = 6;
-		    m_faces.push_back(std::move(BAD));
-
-		    Face CBD;
-		    CBD.m_he = 9;
-		    m_faces.push_back(std::move(CBD));
-		}
-
-		std::array<size_t,3> getVertexIndicesOfFace(const Face& f) const {
-		    std::array<size_t,3> v;
-		    const HalfEdge* he = &m_halfEdges[f.m_he];
-		    v[0] = he->m_endVertex;
-		    he = &m_halfEdges[he->m_next];
-		    v[1] = he->m_endVertex;
-		    he = &m_halfEdges[he->m_next];
-		    v[2] = he->m_endVertex;
-		    return v;
-		}
-
-		std::array<size_t,2> getVertexIndicesOfHalfEdge(const HalfEdge& he) const {
-		    return {{m_halfEdges[he.m_opp].m_endVertex,he.m_endVertex}};
-		}
-
-		std::array<size_t,3> getHalfEdgeIndicesOfFace(const Face& f) const {
-		    return {{f.m_he,m_halfEdges[f.m_he].m_next,m_halfEdges[m_halfEdges[f.m_he].m_next].m_next}};
-		}
-	};
-
-
-
-
-
-    template<typename FloatType, typename size_t>
-	class HalfEdgeMesh {
-	    public:
-
-		struct HalfEdge {
-		    size_t m_endVertex;
-		    size_t m_opp;
-		    size_t m_face;
-		    size_t m_next;
-		};
-
-		struct Face {
-		    size_t m_halfEdgeIndex; // Index of one of the half edges of this face
-		};
-
-		std::vector<Vector3<FloatType>> m_vertices;
-		std::vector<Face> m_faces;
-		std::vector<HalfEdge> m_halfEdges;
-
-		HalfEdgeMesh(const MeshBuilder<FloatType>& builderObject, const VertexDataSource<FloatType>& vertexData )
-		{
-		    std::unordered_map<size_t,size_t> faceMapping;
-		    std::unordered_map<size_t,size_t> halfEdgeMapping;
-		    std::unordered_map<size_t, size_t> vertexMapping;
-
-		    size_t i=0;
-		    for (const auto& face : builderObject.m_faces) {
-			if (!face.isDisabled()) {
-			    m_faces.push_back({static_cast<size_t>(face.m_he)});
-			    faceMapping[i] = m_faces.size()-1;
-
-			    const auto heIndices = builderObject.getHalfEdgeIndicesOfFace(face);
-			    for (const auto heIndex : heIndices) {
-				const size_t vertexIndex = builderObject.m_halfEdges[heIndex].m_endVertex;
-				if (vertexMapping.count(vertexIndex)==0) {
-				    m_vertices.push_back(vertexData[vertexIndex]);
-				    vertexMapping[vertexIndex] = m_vertices.size()-1;
-				}
-			    }
-			}
-			i++;
-		    }
-
-		    i=0;
-		    for (const auto& halfEdge : builderObject.m_halfEdges) {
-			if (!halfEdge.isDisabled()) {
-			    m_halfEdges.push_back({static_cast<size_t>(halfEdge.m_endVertex),static_cast<size_t>(halfEdge.m_opp),static_cast<size_t>(halfEdge.m_face),static_cast<size_t>(halfEdge.m_next)});
-			    halfEdgeMapping[i] = m_halfEdges.size()-1;
-			}
-			i++;
-		    }
-
-		    for (auto& face : m_faces) {
-			assert(halfEdgeMapping.count(face.m_halfEdgeIndex) == 1);
-			face.m_halfEdgeIndex = halfEdgeMapping[face.m_halfEdgeIndex];
-		    }
-
-		    for (auto& he : m_halfEdges) {
-			he.m_face = faceMapping[he.m_face];
-			he.m_opp = halfEdgeMapping[he.m_opp];
-			he.m_next = halfEdgeMapping[he.m_next];
-			he.m_endVertex = vertexMapping[he.m_endVertex];
-		    }
-		}
-
-	};
-
-
-    namespace mathutils {
-
-	template <typename T>
-	    inline T getSquaredDistanceBetweenPointAndRay(const Vector3<T>& p, const Ray<T>& r) {
-		const Vector3<T> s = p-r.m_S;
-		T t = s.dotProduct(r.m_V);
-		return s.getLengthSquared() - t*t*r.m_VInvLengthSquared;
-	    }
-
-	// Note that the unit of distance returned is relative to plane's normal's length (divide by N.getNormalized() if needed to get the "real" distance).
-	template <typename T>
-	    inline T getSignedDistanceToPlane(const Vector3<T>& v, const Plane<T>& p) {
-		return p.m_N.dotProduct(v) + p.m_D;
-	    }
-
-	template <typename T>
-	    inline Vector3<T> getTriangleNormal(const Vector3<T>& a,const Vector3<T>& b,const Vector3<T>& c) {
-		// We want to get (a-c).crossProduct(b-c) without constructing temp vectors
-		T x = a.x - c.x;
-		T y = a.y - c.y;
-		T z = a.z - c.z;
-		T rhsx = b.x - c.x;
-		T rhsy = b.y - c.y;
-		T rhsz = b.z - c.z;
-		T px = y * rhsz - z * rhsy ;
-		T py = z * rhsx - x * rhsz ;
-		T pz = x * rhsy - y * rhsx ;
-		return Vector3<T>(px,py,pz);
-	    }
-
-
-    }
-
-
-    template<typename T>
-	class ConvexHull {
-	    std::shared_ptr<std::vector<Vector3<T>>> m_optimizedVertexBuffer;
-	    VertexDataSource<T> m_vertices;
-	    std::vector<size_t> m_indices;
-	    public:
-	    ConvexHull() {}
-
-	    // Copy constructor
-	    ConvexHull(const ConvexHull& o) {
-		m_indices = o.m_indices;
-		if (o.m_optimizedVertexBuffer) {
-		    m_optimizedVertexBuffer.reset(new std::vector<Vector3<T>>(*o.m_optimizedVertexBuffer));
-		    m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
-		}
-		else {
-		    m_vertices = o.m_vertices;
-		}
-	    }
-
-	    ConvexHull& operator=(const ConvexHull& o) {
-		if (&o == this) {
-		    return *this;
-		}
-		m_indices = o.m_indices;
-		if (o.m_optimizedVertexBuffer) {
-		    m_optimizedVertexBuffer.reset(new std::vector<Vector3<T>>(*o.m_optimizedVertexBuffer));
-		    m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
-		}
-		else {
-		    m_vertices = o.m_vertices;
-		}
-		return *this;
-	    }
-
-	    ConvexHull(ConvexHull&& o) {
-		m_indices = std::move(o.m_indices);
-		if (o.m_optimizedVertexBuffer) {
-		    m_optimizedVertexBuffer = std::move(o.m_optimizedVertexBuffer);
-		    o.m_vertices = VertexDataSource<T>();
-		    m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
-		}
-		else {
-		    m_vertices = o.m_vertices;
-		}
-	    }
-
-	    ConvexHull& operator=(ConvexHull&& o) {
-		if (&o == this) {
-		    return *this;
-		}
-		m_indices = std::move(o.m_indices);
-		if (o.m_optimizedVertexBuffer) {
-		    m_optimizedVertexBuffer = std::move(o.m_optimizedVertexBuffer);
-		    o.m_vertices = VertexDataSource<T>();
-		    m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
-		}
-		else {
-		    m_vertices = o.m_vertices;
-		}
-		return *this;
-	    }
-
-	    // Construct vertex and index buffers from half edge mesh and pointcloud
-	    ConvexHull(const MeshBuilder<T>& mesh, const VertexDataSource<T>& pointCloud, bool CCW, bool useOriginalIndices) {
-		if (!useOriginalIndices) {
-		    m_optimizedVertexBuffer.reset(new std::vector<Vector3<T>>());
-		}
-
-		std::vector<bool> faceProcessed(mesh.m_faces.size(),false);
-		std::vector<size_t> faceStack;
-		std::unordered_map<size_t,size_t> vertexIndexMapping; // Map vertex indices from original point cloud to the new mesh vertex indices
-		for (size_t i = 0;i<mesh.m_faces.size();i++) {
-		    if (!mesh.m_faces[i].isDisabled()) {
-			faceStack.push_back(i);
-			break;
-		    }
-		}
-		if (faceStack.size()==0) {
-		    return;
-		}
-
-		const size_t iCCW = CCW ? 1 : 0;
-		const size_t finalMeshFaceCount = mesh.m_faces.size() - mesh.m_disabledFaces.size();
-		m_indices.reserve(finalMeshFaceCount*3);
-
-		while (faceStack.size()) {
-		    auto it = faceStack.end()-1;
-		    size_t top = *it;
-		    assert(!mesh.m_faces[top].isDisabled());
-		    faceStack.erase(it);
-		    if (faceProcessed[top]) {
-			continue;
-		    }
-		    else {
-			faceProcessed[top]=true;
-			auto halfEdges = mesh.getHalfEdgeIndicesOfFace(mesh.m_faces[top]);
-			size_t adjacent[] = {mesh.m_halfEdges[mesh.m_halfEdges[halfEdges[0]].m_opp].m_face,mesh.m_halfEdges[mesh.m_halfEdges[halfEdges[1]].m_opp].m_face,mesh.m_halfEdges[mesh.m_halfEdges[halfEdges[2]].m_opp].m_face};
-			for (auto a : adjacent) {
-			    if (!faceProcessed[a] && !mesh.m_faces[a].isDisabled()) {
-				faceStack.push_back(a);
-			    }
-			}
-			auto vertices = mesh.getVertexIndicesOfFace(mesh.m_faces[top]);
-			if (!useOriginalIndices) {
-			    for (auto& v : vertices) {
-				auto lit = vertexIndexMapping.find(v);
-				if (lit == vertexIndexMapping.end()) {
-				    m_optimizedVertexBuffer->push_back(pointCloud[v]);
-				    vertexIndexMapping[v] = m_optimizedVertexBuffer->size()-1;
-				    v = m_optimizedVertexBuffer->size()-1;
-				}
-				else {
-				    v = lit->second;
-				}
-			    }
-			}
-			m_indices.push_back(vertices[0]);
-			m_indices.push_back(vertices[1 + iCCW]);
-			m_indices.push_back(vertices[2 - iCCW]);
-		    }
-		}
-
-		if (!useOriginalIndices) {
-		    m_vertices = VertexDataSource<T>(*m_optimizedVertexBuffer);
-		}
-		else {
-		    m_vertices = pointCloud;
-		}
-	    }
-
-	    std::vector<size_t>& getIndexBuffer() {
-		return m_indices;
-	    }
-
-	    VertexDataSource<T>& getVertexBuffer() {
-		return m_vertices;
-	    }
-
-	};
-
-
-    template<typename FloatType>
-	class QuickHull {
-	    using vec3 = Vector3<FloatType>;
-
-	    static const FloatType Epsilon;
-
-	    FloatType m_epsilon, m_epsilonSquared, m_scale;
-	    bool m_planar;
-	    std::vector<vec3> m_planarPointCloudTemp;
-	    VertexDataSource<FloatType> m_vertexData;
-	    MeshBuilder<FloatType> m_mesh;
-	    std::array<size_t,6> m_extremeValues;
-
-	    // Temporary variables used during iteration process
-	    std::vector<size_t> m_newFaceIndices;
-	    std::vector<size_t> m_newHalfEdgeIndices;
-	    std::vector< std::shared_ptr<std::vector<size_t>> > m_disabledFacePointVectors;
-	    std::vector<size_t> m_visibleFaces;
-	    std::vector<size_t> m_horizonEdges;
-	    struct FaceData {
-		size_t m_faceIndex;
-		size_t m_enteredFromHalfEdge; // If the face turns out not to be visible, this half edge will be marked as horizon edge
-		FaceData(size_t fi, size_t he) : m_faceIndex(fi),m_enteredFromHalfEdge(he) {}
-	    };
-	    std::vector<FaceData> m_possiblyVisibleFaces;
-	    std::deque<size_t> m_faceList;
-
-	    // Create a half edge mesh representing the base tetrahedron from
-	    // which the QuickHull iteration proceeds. m_extremeValues must be
-	    // properly set up when this is called.
-	    void setupInitialTetrahedron();
-
-	    // Given a list of half edges, try to rearrange them so that they
-	    // form a loop. Return true on success.
-	    bool reorderHorizonEdges(std::vector<size_t>& horizonEdges);
-
-	    // Find indices of extreme values (max x, min x, max y, min y, max
-	    // z, min z) for the given point cloud
-	    std::array<size_t,6> getExtremeValues();
-
-	    // Compute scale of the vertex data.
-	    FloatType getScale(const std::array<size_t,6>& extremeValues);
-
-	    // Each face contains a unique pointer to a vector of indices.
-	    // However, many - often most - faces do not have any points on the
-	    // positive side of them especially at the the end of the
-	    // iteration. When a face is removed from the mesh, its associated
-	    // point vector, if such exists, is moved to the index vector pool,
-	    // and when we need to add new faces with points on the positive
-	    // side to the mesh, we reuse these vectors. This reduces the
-	    // amount of std::vectors we have to deal with, and impact on
-	    // performance is remarkable.
-	    Pool<std::vector<size_t>> m_indexVectorPool;
-	    inline std::shared_ptr<std::vector<size_t>> getIndexVectorFromPool();
-	    inline void reclaimToIndexVectorPool(std::shared_ptr<std::vector<size_t>>& ptr);
-
-	    // Associates a point with a face if the point resides on the
-	    // positive side of the plane. Returns true if the points was on
-	    // the positive side.
-	    inline bool addPointToFace(typename MeshBuilder<FloatType>::Face& f, size_t pointIndex);
-
-	    // This will update m_mesh from which we create the ConvexHull
-	    // object that getConvexHull function returns
-	    void createConvexHalfEdgeMesh();
-
-	    // Constructs the convex hull into a MeshBuilder object which can
-	    // be converted to a ConvexHull or Mesh object
-	    void buildMesh(const VertexDataSource<FloatType>& pointCloud, bool UNUSED(CCW), bool UNUSED(useOriginalIndices), FloatType eps);
-
-	    // The public getConvexHull functions will setup a VertexDataSource
-	    // object and call this
-	    ConvexHull<FloatType> getConvexHull(const VertexDataSource<FloatType>& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps);
-	    public:
-	    // Computes convex hull for a given point cloud.
-	    // Params:
-	    //   pointCloud: a vector of of 3D points
-	    //
-	    //   CCW: whether the output mesh triangles should have CCW
-	    //   orientation
-	    //
-	    //   useOriginalIndices: should the output mesh use same vertex
-	    //   indices as the original point cloud.  If this is false, then
-	    //   we generate a new vertex buffer which contains only the
-	    //   vertices that are part of the convex hull.
-	    //
-	    //   eps: minimum distance to a plane to consider a point being on
-	    //   positive of it (for a point cloud with scale 1)
-	    ConvexHull<FloatType> getConvexHull(const std::vector<Vector3<FloatType>>& pointCloud, bool CCW, bool useOriginalIndices, FloatType eps = Epsilon);
-
-	    // Computes convex hull for a given point cloud.
-	    // Params:
-	    //
-	    //   vertexData: pointer to the first 3D point of the point cloud
-	    //
-	    //   vertexCount: number of vertices in the point cloud
-	    //
-	    //   CCW: whether the output mesh triangles should have CCW
-	    //   orientation
-	    //
-	    //   useOriginalIndices: should the output mesh use same vertex
-	    //   indices as the original point cloud.  If this is false, then
-	    //   we generate a new vertex buffer which contains only the
-	    //   vertices that are part of the convex hull.
-	    //
-	    //   eps: minimum distance to a plane to consider a point being on
-	    //   positive of it (for a point cloud with scale 1)
-	    ConvexHull<FloatType> getConvexHull(const Vector3<FloatType>* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, FloatType eps = Epsilon);
-
-	    // Computes convex hull for a given point cloud. This
-	    // function assumes that the vertex data resides in
-	    // memory in the following format:
-	    // x_0,y_0,z_0,x_1,y_1,z_1,...
-	    // Params:
-	    //   vertexData: pointer to the X component of the first point of
-	    //   the point cloud.
-	    //
-	    //   vertexCount: number of vertices in the point cloud
-	    //
-	    //   CCW: whether the output mesh triangles should have CCW
-	    //   orientation
-	    //
-	    //   useOriginalIndices: should the output mesh use same vertex
-	    //   indices as the original point cloud.  If this is false, then
-	    //   we generate a new vertex buffer which contains only the
-	    //   vertices that are part of the convex hull.
-	    //
-	    //   eps: minimum distance to a plane to consider a point being on
-	    //   positive of it (for a point cloud with scale 1)
-	    ConvexHull<FloatType> getConvexHull(const FloatType* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, FloatType eps = Epsilon);
-
-	    // Computes convex hull for a given point cloud. This function
-	    // assumes that the vertex data resides in memory in the following
-	    // format: x_0,y_0,z_0,x_1,y_1,z_1,...
-	    // Params:
-	    //   vertexData: pointer to the X component of the first point of the point cloud.
-	    //
-	    //   vertexCount: number of vertices in the point cloud
-	    //
-	    //   CCW: whether the output mesh triangles should have CCW orientation
-	    //
-	    //   eps: minimum distance to a plane to consider a point being on
-	    //   positive of it (for a point cloud with scale 1)
-	    //
-	    // Returns:
-	    //   Convex hull of the point cloud as a mesh object with half edge structure.
-	    HalfEdgeMesh<FloatType, size_t> getConvexHullAsMesh(const FloatType* vertexData, size_t vertexCount, bool CCW, FloatType eps = Epsilon);
-	};
-
-    /*
-     * Inline function definitions
-     */
-
-    template<typename T>
-	std::shared_ptr<std::vector<size_t>> QuickHull<T>::getIndexVectorFromPool() {
-	    auto r = m_indexVectorPool.get();
-	    r->clear();
-	    return r;
-	}
-
-    template<typename T>
-	void QuickHull<T>::reclaimToIndexVectorPool(std::shared_ptr<std::vector<size_t>>& ptr) {
-	    const size_t oldSize = ptr->size();
-	    if ((oldSize+1)*128 < ptr->capacity()) {
-		// Reduce memory usage! Huge vectors are needed at the
-		// beginning of iteration when faces have many points on their
-		// positive side. Later on, smaller vectors will suffice.
-		ptr.reset();
-		return;
-	    }
-	    m_indexVectorPool.reclaim(ptr);
-	}
-
-    template<typename T>
-	bool QuickHull<T>::addPointToFace(typename MeshBuilder<T>::Face& f, size_t pointIndex) {
-	    const T D = mathutils::getSignedDistanceToPlane(m_vertexData[ pointIndex ],f.m_P);
-	    if (D>0 && D*D > m_epsilonSquared*f.m_P.m_sqrNLength) {
-		if (!f.m_pointsOnPositiveSide) {
-		    f.m_pointsOnPositiveSide = std::move(getIndexVectorFromPool());
-		}
-		f.m_pointsOnPositiveSide->push_back( pointIndex );
-		if (D > f.m_mostDistantPointDist) {
-		    f.m_mostDistantPointDist = D;
-		    f.m_mostDistantPoint = pointIndex;
-		}
-		return true;
-	    }
-	    return false;
-	}
-
-}
-
-
-#endif /* QUICKHULL_HPP_ */
-
-// Local Variables:
-// tab-width: 8
-// mode: C++
-// c-basic-offset: 4
-// indent-tabs-mode: t
-// c-file-style: "stroustrup"
-// End:
-// ex: shiftwidth=4 tabstop=8
diff --git a/bindings/java/src/main/cpp/manifold3d/QuickHullTomilov.hpp b/bindings/java/src/main/cpp/manifold3d/QuickHullTomilov.hpp
deleted file mode 100644
index 670352599..000000000
--- a/bindings/java/src/main/cpp/manifold3d/QuickHullTomilov.hpp
+++ /dev/null
@@ -1,1029 +0,0 @@
-/* Quickhull algorithm implementation
- *
- * Copyright (c) 2014-2015, Anatoliy V. Tomilov
- * All rights reserved.
- *
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following condition is met:
- * Redistributions of source code must retain the above copyright notice, this condition and the following disclaimer.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- */
-#pragma once
-
-#include <type_traits>
-#include <vector>
-#include <deque>
-#include <list>
-#include <map>
-#include <unordered_set>
-#include <unordered_map>
-#include <iterator>
-#include <memory>
-#include <algorithm>
-#include <numeric>
-#include <utility>
-#include <functional>
-
-#include <cstdint>
-#include <cmath>
-#include <cassert>
-
-namespace QuickHullTomilov {
-
-template< typename point_iterator,
-          typename value_type = std::decay_t< decltype(*std::cbegin(std::declval< typename std::iterator_traits< point_iterator >::value_type >())) > >
-struct quick_hull
-{
-
-    static_assert(std::is_base_of< std::forward_iterator_tag, typename std::iterator_traits< point_iterator >::iterator_category >::value,
-                  "multipass guarantee required");
-
-    using size_type = std::size_t;
-
-    size_type const dimension_;
-    value_type const & eps;
-
-    value_type const zero = value_type(0);
-    value_type const one = value_type(1);
-
-    using vector = std::vector< value_type >;
-
-private :
-
-    using vrow = value_type *;
-    using crow = value_type const *;
-    using matrix = std::vector< vrow >;
-
-    vector storage_;
-    vrow inner_point_;
-    matrix matrix_;
-    matrix det_matrix_;
-    matrix shadow_matrix_;
-
-public :
-
-    quick_hull(size_type, value_type const &&) = delete; // bind eps to lvalue only
-
-    quick_hull(size_type const _dimension,
-               value_type const & _eps)
-        : dimension_(_dimension)
-        , eps(_eps)
-        , storage_(dimension_ * dimension_ * 2 + dimension_)
-        , inner_point_(storage_.data())
-        , matrix_(dimension_)
-        , det_matrix_(dimension_)
-        , shadow_matrix_(dimension_)
-        , vertices_hashes_(dimension_)
-    {
-        assert(1 < dimension_);
-        assert(!(eps < zero));
-        for (vrow & row_ : matrix_) {
-            row_ = inner_point_;
-            inner_point_ += dimension_;
-        }
-        for (vrow & row_ : shadow_matrix_) {
-            row_ = inner_point_;
-            inner_point_ += dimension_;
-        }
-        assert(inner_point_ + dimension_ == &storage_.back() + 1);
-    }
-
-    using point_array = std::vector< point_iterator >;
-    using point_list  = std::list< point_iterator >;
-    using point_deque = std::deque< point_iterator >;
-    using facet_array = std::vector< size_type >;
-
-    struct facet // (d - 1)-dimensional face
-    {
-
-        // each neighbouring facet lies against corresponding vertex and vice versa
-        point_array vertices_; // dimension_ points (oriented)
-        facet_array neighbours_; // dimension_ neighbouring facets
-
-        point_list outside_; // if empty, then is convex hull's facet, else the first point (i.e. outside_.front()) is the furthest point from this facet
-        point_deque coplanar_; // containing coplanar points and vertices of coplanar facets as well
-
-        // equation of supporting hyperplane
-        vector normal_; // components of normalized normal vector
-        value_type D; // distance from the origin to the hyperplane
-
-        template< typename iterator >
-        value_type
-        distance(iterator const _point) const
-        {
-            using iterator_traits = std::iterator_traits< iterator >;
-            static_assert(std::is_base_of< std::input_iterator_tag, typename iterator_traits::iterator_category >::value);
-            return std::inner_product(std::cbegin(normal_), std::cend(normal_), _point, D);
-        }
-
-    };
-
-    using facets = std::deque< facet >;
-
-    facets facets_;
-
-    value_type
-    cos_of_dihedral_angle(facet const & _first, facet const & _second) const
-    {
-        return std::inner_product(std::cbegin(_first.normal_), std::cend(_first.normal_), std::cbegin(_second.normal_), zero);
-    }
-
-private :
-
-    void
-    make_facet(facet & _facet,
-               point_array const & _vertices,
-               size_type const _against,
-               point_iterator const _apex,
-               size_type const _neighbour)
-    {
-        assert(_vertices.size() == dimension_);
-        _facet.vertices_ = _vertices;
-        _facet.vertices_[_against] = _apex;
-        _facet.neighbours_.resize(dimension_);
-        _facet.neighbours_[_against] = _neighbour;
-        _facet.normal_.resize(dimension_);
-    }
-
-    template< typename iterator >
-    void
-    make_facet(facet & _facet,
-               iterator sbeg, // simplex
-               size_type const _vertex,
-               bool const _swap)
-    {
-        using iterator_traits = std::iterator_traits< iterator >;
-        static_assert(std::is_base_of< std::input_iterator_tag, typename iterator_traits::iterator_category >::value);
-        static_assert(std::is_constructible< point_iterator, typename iterator_traits::value_type >::value);
-        _facet.vertices_.reserve(dimension_);
-        _facet.neighbours_.reserve(dimension_);
-        for (size_type v = 0; v <= dimension_; ++v) {
-            if (v != _vertex) {
-                _facet.vertices_.push_back(*sbeg);
-                _facet.neighbours_.push_back(v);
-            }
-            ++sbeg;
-        }
-        if (_swap == (((dimension_ - _vertex) % 2) == 0)) {
-            using std::swap;
-            swap(_facet.vertices_.front(), _facet.vertices_.back());
-            swap(_facet.neighbours_.front(), _facet.neighbours_.back());
-        }
-        _facet.normal_.resize(dimension_);
-    }
-
-    void
-    reuse_facet(facet & _facet,
-                point_array const & _vertices,
-                size_type const _against,
-                point_iterator const _apex,
-                size_type const _neighbour)
-    {
-        assert(_vertices.size() == dimension_);
-        _facet.vertices_ = _vertices;
-        _facet.vertices_[_against] = _apex;
-        assert(_facet.neighbours_.size() == dimension_);
-        _facet.neighbours_[_against] = _neighbour;
-        assert(_facet.normal_.size() == dimension_);
-    }
-
-    void
-    copy_point(point_iterator const _from, vrow _to) const
-    {
-        std::copy_n(std::cbegin(*_from), dimension_, _to);
-    }
-
-    void
-    subtract(vrow _minuend, crow _subtrahend) const
-    {
-        for (size_type i = 0; i < dimension_; ++i) {
-            *_minuend++ -= *_subtrahend++;
-        }
-    }
-
-    void
-    gshift(vrow _augend, value_type const & _addend) const // shift Gaussian row
-    {
-        for (size_type i = 0; i < dimension_; ++i) {
-            *_augend++ += _addend;
-        }
-        *_augend += _addend;
-    }
-
-    void
-    divide(vrow _dividend, value_type const & _divisor) const
-    {
-        for (size_type i = 0; i < dimension_; ++i) {
-            *_dividend++ /= _divisor;
-        }
-    }
-
-    void
-    subtract_and_assign(vrow _assignee, vrow _minuend, crow _subtrahend) const
-    {
-        for (size_type i = 0; i < dimension_; ++i) {
-            *_assignee++ = (*_minuend++ -= *_subtrahend++);
-        }
-    }
-
-    void
-    multiply_and_add(vrow _assignee, crow _multiplicand, value_type const & _factor) const
-    {
-        for (size_type i = 0; i < dimension_; ++i) {
-            *_assignee++ += (*_multiplicand++ * _factor);
-        }
-    }
-
-    void
-    scale_and_shift(vrow _multiplicand, crow _direction, value_type const & _factor) const
-    {
-        for (size_type i = 0; i < dimension_; ++i) {
-            (*_multiplicand++ *= _factor) += *_direction++;
-        }
-    }
-
-    void
-    matrix_transpose_copy(point_array const & _vertices)
-    {
-        for (size_type r = 0; r < dimension_; ++r) {
-            auto v = std::cbegin(*_vertices[r]);
-            for (size_type c = 0; c < dimension_; ++c) {
-                shadow_matrix_[c][r] = *v;
-                ++v;
-            }
-        }
-    }
-
-    void
-    matrix_restore(size_type const _identity)
-    {
-        for (size_type c = 0; c < _identity; ++c) {
-            std::copy_n(shadow_matrix_[c], dimension_, matrix_[c]);
-        }
-        std::fill_n(matrix_[_identity], dimension_, one);
-        for (size_type c = _identity + 1; c < dimension_; ++c) {
-            std::copy_n(shadow_matrix_[c], dimension_, matrix_[c]);
-        }
-    }
-
-    void
-    matrix_restore()
-    {
-        for (size_type c = 0; c < dimension_; ++c) {
-            std::copy_n(shadow_matrix_[c], dimension_, matrix_[c]);
-        }
-    }
-
-    void
-    matrix_sqr(size_type const _size)
-    { // shadow_matrix_ = matrix_ * transposed matrix_
-        assert(_size < dimension_);
-        for (size_type r = 0; r < _size; ++r) {
-            vrow const lhs_ = shadow_matrix_[r];
-            crow const row_ = matrix_[r];
-            for (size_type c = 0; c < _size; ++c) {
-                lhs_[c] = std::inner_product(row_, row_ + _size, matrix_[c], zero);
-            }
-        }
-    }
-
-    // based on LUP decomposition (complexity is (d^3 / 3 + d^2 / 2 - 5 * d / 6) vs (2 * d^3 / 3 + d^2 + d / 3 - 2) for QR decomposition via Householder reflections)
-    value_type
-    det(matrix const & _matrix, size_type const _dimension) // hottest function (52% of runtime for D10)
-    { // det_matrix_ contains lower unit triangular matrix and upper triangular at return
-        assert(0 < _dimension);
-        value_type det_ = one;
-        std::copy_n(std::cbegin(_matrix), _dimension, std::begin(det_matrix_));
-        for (size_type i = 0; i < _dimension; ++i) {
-            vrow & mi_ = det_matrix_[i];
-            size_type pivot = i;
-            {
-                using std::abs;
-                value_type max_ = abs(mi_[i]);
-                size_type j = i;
-                while (++j < _dimension) {
-                    value_type y_ = abs(det_matrix_[j][i]);
-                    if (max_ < y_) {
-                        max_ = std::move(y_);
-                        pivot = j;
-                    }
-                }
-                if (!(eps < max_)) { // regular?
-                    det_ = zero; // singular
-                    break;
-                }
-            }
-            if (pivot != i) {
-                det_ = -det_; // each permutation flips sign of det
-                std::swap(mi_, det_matrix_[pivot]);
-            }
-            value_type const & dia_ = mi_[i];
-            det_ *= dia_; // det is multiple of diagonal elements
-            size_type j = i;
-            while (++j < _dimension) {
-                vrow const mj_ = det_matrix_[j];
-                value_type & mji_ = mj_[i];
-                mji_ /= dia_;
-                size_type k = i;
-                while (++k < _dimension) {
-                    mj_[k] -= mji_ * mi_[k];
-                }
-            }
-        }
-        return det_;
-    }
-
-    value_type
-    det()
-    {
-        return det(matrix_, dimension_);
-    }
-
-    void
-    set_hyperplane_equation(facet & _facet)
-    {
-        matrix_transpose_copy(_facet.vertices_);
-        matrix_restore();
-        _facet.D = -det();
-        value_type N = zero;
-        for (size_type i = 0; i < dimension_; ++i) {
-            matrix_restore(i);
-            value_type & n = _facet.normal_[i];
-            n = det();
-            N += n * n;
-        }
-        using std::sqrt;
-        N = sqrt(std::move(N));
-        divide(_facet.normal_.data(), N);
-        _facet.D /= std::move(N);
-        assert(_facet.distance(inner_point_) < zero);
-    }
-
-    bool
-    orthonormalize(point_list const & _affine_space,
-                   size_type const _rank,
-                   crow const _origin)
-    {
-        assert(!(dimension_ < _rank));
-        assert(!(_affine_space.size() < _rank));
-        auto vertex = std::begin(_affine_space);
-        for (size_type r = 0; r < _rank; ++r) { // affine space -> vector space
-            vrow const row_ = shadow_matrix_[r];
-            copy_point(*vertex, row_);
-            subtract(row_, _origin);
-            ++vertex;
-        }
-        for (size_type i = 0; i < _rank; ++i) { // Householder transformation
-            value_type sum_ = zero;
-            vrow const qri_ = shadow_matrix_[i];
-            for (size_type k = i; k < dimension_; ++k) {
-                value_type const & qrik_ = qri_[k];
-                sum_ += qrik_ * qrik_;
-            }
-            using std::sqrt;
-            value_type norm_ = sqrt(sum_);
-            if (!(eps < norm_)) {
-                return false;
-            }
-            value_type & qrii_ = qri_[i];
-            if (qrii_ < zero) {
-                norm_ = -norm_;
-            }
-            value_type factor_ = sqrt(std::move(sum_) + qrii_ * norm_);
-            if (!(eps < factor_)) {
-                return false;
-            }
-            qrii_ += std::move(norm_);
-            for (size_type k = i; k < dimension_; ++k) {
-                qri_[k] /= factor_;
-            }
-            size_type j = i;
-            while (++j < _rank) {
-                vrow const qrj_ = shadow_matrix_[j];
-                value_type s_ = zero;
-                for (size_type k = i; k < dimension_; ++k) {
-                    s_ += qri_[k] * qrj_[k];
-                }
-                for (size_type k = i; k < dimension_; ++k) {
-                    qrj_[k] -= qri_[k] * s_;
-                }
-            }
-        } // shadow_matrix_ is packed QR
-        return true;
-    }
-
-    void
-    forward_transformation(size_type const _rank) // calculation of Q
-    {
-        assert(!(dimension_ < _rank));
-        for (size_type i = 0; i < _rank; ++i) {
-            vrow const qi_ = matrix_[i];
-            std::fill_n(qi_, dimension_, zero);
-            qi_[i] = one;
-            size_type j = _rank;
-            while (0 < j) {
-                --j;
-                vrow const qrj_ = shadow_matrix_[j]; // containing packed QR
-                value_type s_ = zero;
-                for (size_type k = j; k < dimension_; ++k) {
-                    s_ += qrj_[k] * qi_[k];
-                }
-                for (size_type k = j; k < dimension_; ++k) {
-                    qi_[k] -= qrj_[k] * s_;
-                }
-            }
-        } // matrix_ is Q
-    }
-
-    bool
-    steal_best(point_list & _basis)
-    { // set moves a point which is furthest from affine subspace formed by points of "_basis" set from "outside_" set to "_basis"
-        assert(!_basis.empty());
-        size_type const rank_ = _basis.size() - 1;
-        assert(rank_ < dimension_);
-        vrow const origin_ = matrix_[rank_];
-        copy_point(_basis.back(), origin_);
-        if (!orthonormalize(_basis, rank_, origin_)) {
-            return false;
-        }
-        forward_transformation(rank_);
-        vrow const projection_ = shadow_matrix_.back();
-        vrow const apex_ = shadow_matrix_.front();
-        value_type distance_ = zero; // square of distance to the subspace
-        auto const oend = std::cend(outside_);
-        auto furthest = oend;
-        for (auto it = std::cbegin(outside_); it != oend; ++it) {
-            copy_point(*it, apex_);
-            subtract_and_assign(projection_, apex_, origin_); // turn translated space into vector space then project onto orthogonal subspace
-            for (size_type i = 0; i < rank_; ++i) {
-                crow const qi_ = matrix_[i];
-                multiply_and_add(projection_, qi_, -std::inner_product(qi_, qi_ + dimension_, apex_, zero));
-            }
-            value_type d_ = std::inner_product(projection_, projection_ + dimension_, projection_, zero);
-            if (distance_ < d_) {
-                distance_ = std::move(d_);
-                furthest = it;
-            }
-        }
-        if (furthest == oend) {
-            return false;
-        }
-        _basis.splice(std::cend(_basis), outside_, furthest);
-        return true;
-    }
-
-    facet_array removed_facets_;
-
-    std::pair< facet &, size_type const >
-    add_facet(point_array const & _vertices,
-              size_type const _against,
-              point_iterator const _apex,
-              size_type const _neighbour)
-    {
-        if (removed_facets_.empty()) {
-            size_type const f = facets_.size();
-            facets_.emplace_back();
-            facet & facet_ = facets_.back();
-            make_facet(facet_, _vertices, _against, _apex, _neighbour);
-            return {facet_, f};
-        } else {
-            size_type const f = removed_facets_.back();
-            removed_facets_.pop_back();
-            facet & facet_ = facets_[f];
-            reuse_facet(facet_, _vertices, _against, _apex, _neighbour);
-            return {facet_, f};
-        }
-    }
-
-    using ranking = std::multimap< value_type, size_type >;
-    using ranking_meta = std::unordered_map< size_type, typename ranking::iterator >;
-
-    ranking ranking_;
-    ranking_meta ranking_meta_;
-
-    void
-    rank(value_type && _orientation,
-         size_type const f)
-    {
-        if (eps < _orientation) {
-            ranking_meta_.emplace(f, ranking_.emplace(std::move(_orientation), f));
-        }
-    }
-
-    void
-    unrank(size_type const f)
-    {
-        auto const r = ranking_meta_.find(f);
-        if (r != std::end(ranking_meta_)) {
-            ranking_.erase(r->second);
-            ranking_meta_.erase(r);
-        }
-        removed_facets_.push_back(f);
-    }
-
-    point_list outside_;
-
-    value_type
-    partition(facet & _facet)
-    {
-        value_type distance_ = zero;
-        auto it = std::cbegin(outside_);
-        auto const oend = std::cend(outside_);
-        while (it != oend) {
-            auto const next = std::next(it);
-            value_type d_ = _facet.distance(std::cbegin(**it));
-            if (eps < d_) {
-                if (distance_ < d_) {
-                    distance_ = std::move(d_);
-                    _facet.outside_.splice(std::cbegin(_facet.outside_), outside_, it);
-                } else {
-                    _facet.outside_.splice(std::cend(_facet.outside_), outside_, it);
-                }
-            } else if (!(d_ < -eps)) {
-                _facet.coplanar_.push_back(*it);
-            }
-            it = next;
-        }
-        return distance_;
-    }
-
-    size_type
-    get_best_facet() const
-    {
-        assert(ranking_meta_.size() == ranking_.size());
-        return std::prev(std::cend(ranking_))->second;
-    }
-
-    void
-    replace_neighbour(size_type const f,
-                      size_type const _from,
-                      size_type const _to)
-    {
-        if (_from != _to) {
-            for (size_type & n : facets_[f].neighbours_) {
-                if (n == _from) {
-                    n = _to;
-                    return;
-                }
-            }
-        }
-    }
-
-    struct ridge
-    {
-
-        facet & facet_;
-        size_type const f;
-        size_type const v;
-        size_type const hash_;
-
-        bool
-        operator == (ridge const & _rhs) const noexcept
-        {
-            point_iterator const lskip = facet_.vertices_[v];
-            point_iterator const rskip = _rhs.facet_.vertices_[_rhs.v];
-            for (point_iterator const & l : facet_.vertices_) {
-                if (l != lskip) {
-                    bool found_ = false;
-                    for (point_iterator const & r : _rhs.facet_.vertices_) {
-                        if (r != rskip) {
-                            if (l == r) {
-                                found_ = true; // O(D^2) expensive
-                                break;
-                            }
-                        }
-                    }
-                    if (!found_) {
-                        return false;
-                    }
-                }
-            }
-            return true;
-        }
-
-    };
-
-    struct ridge_hash
-    {
-
-        size_type
-        operator () (ridge const & _ridge) const noexcept
-        {
-            return _ridge.hash_;
-        }
-
-    };
-
-    std::unordered_set< ridge, ridge_hash > unique_ridges_;
-    std::hash< typename std::iterator_traits< point_iterator >::value_type const * > point_hash_;
-    std::vector< size_type > vertices_hashes_;
-
-    void
-    find_adjacent_facets(facet & _facet,
-                         size_type const f,
-                         size_type const _skip)
-    {
-        size_type ridge_hash_ = 0;
-        for (size_type v = 0; v < dimension_; ++v) {
-            if (v != _skip) {
-                ridge_hash_ ^= (vertices_hashes_[v] = point_hash_(std::addressof(*_facet.vertices_[v])));
-            }
-        }
-        for (size_type v = 0; v < dimension_; ++v) {
-            if (v != _skip) { // neighbouring facet against apex (_skip-indexed) is known atm
-                auto const position = unique_ridges_.insert({_facet, f, v, (ridge_hash_ ^ vertices_hashes_[v])});
-                if (!position.second) {
-                    ridge const & ridge_ = *position.first;
-                    ridge_.facet_.neighbours_[ridge_.v] = f;
-                    _facet.neighbours_[v] = ridge_.f;
-                    unique_ridges_.erase(position.first);
-                }
-            }
-        }
-    }
-
-    using facet_unordered_set = std::unordered_set< size_type >;
-
-    facet_unordered_set visited_;
-    facet_unordered_set visible_;
-
-    bool
-    process_visibles(facet_array & _newfacets,
-                     size_type const f,
-                     point_iterator const _apex) // traverse the graph of visible facets
-    {
-        assert(!(visited_.size() < visible_.size()));
-        if (!visited_.insert(f).second) {
-            return (visible_.count(f) != 0);
-        }
-        facet & facet_ = facets_[f];
-        if (!(zero < facet_.distance(std::cbegin(*_apex)))) {
-            return false;
-        }
-        visible_.insert(f);
-        outside_.splice(std::cend(outside_), std::move(facet_.outside_));
-        facet_.coplanar_.clear();
-        for (size_type v = 0; v < dimension_; ++v) {
-            size_type const neighbour = facet_.neighbours_[v];
-            if (!process_visibles(_newfacets, neighbour, _apex)) { // recursive function
-                auto const newfacet = add_facet(facet_.vertices_, v, _apex, neighbour);
-                set_hyperplane_equation(newfacet.first);
-                _newfacets.push_back(newfacet.second);
-                replace_neighbour(neighbour, f, newfacet.second);
-                find_adjacent_facets(newfacet.first, newfacet.second, v);
-            }
-        }
-        unrank(f);
-        return true;
-    }
-
-    void
-    compactify()
-    {
-        size_type source = facets_.size();
-        assert(removed_facets_.size() < source);
-        assert(dimension_ < source - removed_facets_.size());
-        assert(ranking_.size() == ranking_meta_.size());
-        assert(!(source < ranking_.size()));
-        auto const rend = std::end(ranking_meta_);
-        std::sort(std::rbegin(removed_facets_), std::rend(removed_facets_));
-        for (size_type const destination : removed_facets_) {
-            assert(!(source < destination));
-            if (destination != --source) {
-                facet & facet_ = facets_[destination];
-                facet_ = std::move(facets_.back());
-                for (size_type const n : facet_.neighbours_) {
-                    replace_neighbour(n, source, destination);
-                }
-                auto const r = ranking_meta_.find(source);
-                if (r != rend) {
-                    r->second->second = destination;
-                    ranking_meta_.emplace(destination, std::move(r->second));
-                    ranking_meta_.erase(r);
-                }
-            }
-            facets_.pop_back();
-        }
-        removed_facets_.clear();
-    }
-
-    bool
-    check_local_convexity(facet const & facet_,
-                          size_type const f) const
-    {
-        assert(&facets_[f] == &facet_);
-        for (size_type const n : facet_.neighbours_) {
-            facet const & neighbour_ = facets_[n];
-            if (cos_of_dihedral_angle(facet_, neighbour_) < one) { // avoid roundoff error
-                for (size_type v = 0; v < dimension_; ++v) {
-                    if (neighbour_.neighbours_[v] == f) { // vertex v of neigbour_ facet is opposite to facet_
-                        value_type const distance_ = facet_.distance(std::cbegin(*neighbour_.vertices_[v]));
-                        if (eps < distance_) {
-                            return false; // facet is not locally convex at ridge, common for facet_ and neighbour_ facets
-                        } else {
-                            break;
-                        }
-                    }
-                }
-            }
-        }
-        return true;
-    }
-
-public :
-
-    template< typename iterator >
-    value_type
-    hypervolume(iterator first,
-                iterator const last) // hypervolume of parallelotope spanned on vectors from last vertex (vlast) to all the vertices lies in [vfirst, vlast)
-    {
-        using iterator_traits = std::iterator_traits< iterator >;
-        static_assert(std::is_base_of< std::input_iterator_tag, typename iterator_traits::iterator_category >::value);
-        static_assert(std::is_constructible< point_iterator, typename iterator_traits::value_type >::value);
-        if (first == last) {
-            return zero;
-        }
-        vrow const origin_ = shadow_matrix_.back();
-        copy_point(*last, origin_);
-        size_type rank_ = 0;
-        while (first != last) { // affine space -> vector space
-            assert(rank_ < dimension_);
-            vrow const row_ = matrix_[rank_];
-            copy_point(*first, row_);
-            subtract(row_, origin_);
-            ++rank_;
-            ++first;
-        }
-        if (rank_ == dimension_) {
-            return det(); // oriented hypervolume
-        } else {
-            matrix_sqr(rank_);
-            using std::sqrt;
-            return sqrt(det(shadow_matrix_, rank_)); // non-oriented rank_-dimensional measure
-        }
-    }
-
-    void
-    add_points(point_iterator beg,
-               point_iterator const end) // [beg; end)
-    {
-        while (beg != end) {
-            outside_.push_back(beg);
-            ++beg;
-        }
-    }
-
-    template< typename iterator >
-    void
-    add_points(iterator const beg,
-               iterator const end) // [beg; end)
-    {
-        using iterator_traits = std::iterator_traits< iterator >;
-        static_assert(std::is_base_of< std::input_iterator_tag, typename iterator_traits::iterator_category >::value);
-        static_assert(std::is_constructible< point_iterator, typename iterator_traits::value_type >::value);
-        std::copy(beg, end, std::back_inserter(outside_));
-    }
-
-    point_list
-    get_affine_basis()
-    {
-        assert(facets_.empty());
-        point_list basis_;
-        basis_.splice(std::cend(basis_), outside_, std::begin(outside_));
-        if (!steal_best(basis_)) {
-            return basis_; // can't find affinely independent second point
-        }
-        outside_.splice(std::cbegin(outside_), basis_, std::begin(basis_)); // reject first point to rejudge it
-        for (size_type i = 0; i < dimension_; ++i) {
-            if (!steal_best(basis_)) {
-                return basis_; // can't find (i + 2) affinely independent points
-            }
-        }
-        return basis_;
-    }
-
-    template< typename iterator >
-    value_type
-    create_initial_simplex(iterator const first,
-                           iterator const last) // [bfirst; blast]
-    {
-        using iterator_traits = std::iterator_traits< iterator >;
-        static_assert(std::is_base_of< std::forward_iterator_tag, typename iterator_traits::iterator_category >::value);
-        static_assert(std::is_constructible< point_iterator, typename iterator_traits::value_type >::value);
-        assert(static_cast< size_type >(std::distance(first, last)) == dimension_);
-        assert(facets_.empty());
-        {
-            copy_point(*last, inner_point_);
-            auto it = first;
-            while (it != last) {
-                auto x = std::cbegin(**it);
-                for (size_type i = 0; i < dimension_; ++i) {
-                    inner_point_[i] += *x;
-                    ++x;
-                }
-                ++it;
-            }
-            divide(inner_point_, value_type(dimension_ + 1));
-        }
-        value_type const volume_ = hypervolume(first, last);
-        bool const swap_ = (volume_ < zero);
-        for (size_type f = 0; f <= dimension_; ++f) {
-            facets_.emplace_back();
-            facet & facet_ = facets_.back();
-            make_facet(facet_, first, f, swap_);
-            set_hyperplane_equation(facet_);
-            rank(partition(facet_), f);
-        }
-        outside_.clear();
-        assert(check());
-        return volume_;
-    }
-
-    // Barber, C. B., D.P. Dobkin, and H.T. Huhdanpaa, 1995. "The Quickhull Algorithm for Convex Hulls", ACM Transactions on Mathematical Software.
-    void
-    create_convex_hull()
-    {
-        assert(facets_.size() == dimension_ + 1);
-        assert(removed_facets_.empty());
-        facet_array newfacets_;
-        while (!ranking_.empty()) {
-            size_type const f = get_best_facet();
-            point_list & o_ = facets_[f].outside_;
-            assert(!o_.empty());
-            point_iterator const apex = std::move(o_.front());
-            o_.pop_front();
-            if (!process_visibles(newfacets_, f, apex)) {
-                assert(false);
-            }
-            visited_.clear();
-            visible_.clear();
-            assert(unique_ridges_.empty());
-            for (size_type const n : newfacets_) {
-                facet & facet_ = facets_[n];
-                assert(check_local_convexity(facet_, n));
-                rank(partition(facet_), n);
-            }
-            newfacets_.clear();
-            outside_.clear();
-            //assert((compactify(), check()));
-        }
-        assert(ranking_meta_.empty());
-        compactify();
-    }
-
-    // Kurt Mehlhorn, Stefan Näher, Thomas Schilz, Stefan Schirra, Michael Seel, Raimund Seidel, and Christian Uhrig.
-    // Checking geometric programs or verification of geometric structures. In Proc. 12th Annu. ACM Sympos. Comput. Geom., pages 159–165, 1996.
-    bool
-    check() const
-    {
-        assert(dimension_ < facets_.size());
-        size_type facets_count_ = 0;
-        for (facet const & facet_ : facets_) { // check local convexity of all the facets
-            if (!check_local_convexity(facet_, facets_count_)) {
-                return false;
-            }
-            ++facets_count_;
-        }
-        facet const & first_ = facets_.front();
-        {
-            value_type const distance_ = first_.distance(inner_point_);
-            if (!(distance_ < zero)) {
-                return false; // inner point is not on negative side of the first facet, therefore structure is not convex
-            }
-        }
-        vector memory_(dimension_ * (4 + dimension_), zero);
-        vrow centroid_ = memory_.data();
-        vrow const ray_ = centroid_;
-        centroid_ += dimension_;
-        for (point_iterator const & v : first_.vertices_) {
-            auto x = std::cbegin(*v);
-            for (size_type i = 0; i < dimension_; ++i) {
-                ray_[i] += *x;
-                ++x;
-            }
-        }
-        divide(ray_, value_type(dimension_));
-        subtract(ray_, inner_point_);
-        {
-            value_type const dot_product_ = std::inner_product(ray_, ray_ + dimension_, first_.normal_.data(), zero);
-            if (!(zero < dot_product_)) { // ray is parallel to the plane or directed away from the plane
-                return false;
-            }
-        }
-        matrix g_{dimension_}; // storage (d * (d + 1)) for Gaussian elimination with partial pivoting
-        for (vrow & row_ : g_) {
-            row_ = centroid_;
-            centroid_ += (dimension_ + 1);
-        }
-        vrow const intersection_point_ = centroid_;
-        centroid_ += dimension_;
-        assert(centroid_ + dimension_ == &memory_.back() + 1);
-        for (size_type f = 1; f < facets_count_; ++f) {
-            using std::abs;
-            facet const & facet_ = facets_[f];
-            value_type const numerator_ = facet_.distance(inner_point_);
-            if (!(numerator_ < zero)) {
-                return false; // inner point is not on negative side of all the facets, i.e. structure is not convex
-            }
-            value_type const denominator_ = std::inner_product(ray_, ray_ + dimension_, facet_.normal_.data(), zero);
-            if (!(zero < denominator_)) { // ray is parallel to the plane or directed away from the plane
-                continue;
-            }
-            std::copy_n(ray_, dimension_, intersection_point_);
-            scale_and_shift(intersection_point_, inner_point_, -(numerator_ / denominator_));
-            for (size_type v = 0; v < dimension_; ++v) {
-                auto beg = std::cbegin(*facet_.vertices_[v]);
-                for (size_type r = 0; r < dimension_; ++r) {
-                    g_[r][v] = *beg;
-                    ++beg;
-                }
-            }
-            for (size_type r = 0; r < dimension_; ++r) {
-                vrow const gr_ = g_[r];
-                centroid_[r] = -std::accumulate(gr_, gr_ + dimension_, zero) / value_type(dimension_);
-                gr_[dimension_] = intersection_point_[r];
-            }
-            for (size_type r = 0; r < dimension_; ++r) {
-                vrow const gr_ = g_[r];
-                value_type & x_ = centroid_[r];
-                gshift(gr_, x_);
-                //assert(!(eps * value_type(dimension_) < std::accumulate(gr_, gr_ + dimension_, zero))); // now center of the facet coincides with the origin, but no one vertex does
-                auto const bounding_box = std::minmax_element(gr_, gr_ + dimension_);
-                x_ = *bounding_box.second - *bounding_box.first;
-                if (!(eps * value_type(dimension_) < x_)) {
-                    x_ = one;
-                }
-            }
-            for (size_type r = 0; r < dimension_; ++r) {
-                vrow const gr_ = g_[r];
-                gshift(gr_, centroid_[r]);
-            }
-            for (size_type i = 0; i < dimension_; ++i) { // Gaussian elimination
-                vrow & gi_ = g_[i];
-                value_type max_ = abs(gi_[i]);
-                size_type pivot = i;
-                {
-                    size_type p = i;
-                    while (++p < dimension_) {
-                        value_type y_ = abs(g_[p][i]);
-                        if (max_ < y_) {
-                            max_ = std::move(y_);
-                            pivot = p;
-                        }
-                    }
-                }
-                assert(eps < max_); // vertex must not match the origin after above transformations
-                if (pivot != i) {
-                    std::swap(gi_, g_[pivot]);
-                }
-                value_type & gii_ = gi_[i];
-                for (size_type j = i + 1; j < dimension_; ++j) {
-                    vrow const gj_ = g_[j];
-                    value_type & gji_ = gj_[i];
-                    gji_ /= gii_;
-                    for (size_type k = i + 1; k <= dimension_; ++k) {
-                        gj_[k] -= gji_ * gi_[k];
-                    }
-                    gji_ = zero;
-                }
-            } // g_ is upper triangular now
-            bool in_range_ = true;
-            {
-                size_type i = dimension_;
-                while (0 < i) {
-                    --i;
-                    vrow const gi_ = g_[i];
-                    value_type & xi_ = gi_[dimension_];
-                    for (size_type j = i + 1; j < dimension_; ++j) {
-                        xi_ -= gi_[j] * g_[j][dimension_];
-                    }
-                    value_type const & gii_ = gi_[i];
-                    assert(eps < abs(gii_)); // vertex must not match the origin
-                    xi_ /= gii_;
-                    if ((xi_ < zero) || (one < xi_)) {
-                        in_range_ = false; // barycentric coordinate does not lie in [0;1] interval => miss
-                        break;
-                    }
-                }
-            }
-            if (in_range_) {
-                return false; // hit
-            }
-        }
-        return true;
-    }
-
-};
-
-}
diff --git a/bindings/java/src/main/cpp/manifold3d/convex_hull.hpp b/bindings/java/src/main/cpp/manifold3d/convex_hull.hpp
deleted file mode 100644
index a2333a2bc..000000000
--- a/bindings/java/src/main/cpp/manifold3d/convex_hull.hpp
+++ /dev/null
@@ -1,754 +0,0 @@
-#pragma once
-
-#include <glm/glm.hpp>
-
-#include <QuickHullTomilov.hpp>
-#include <QuickHull.hpp>
-#include <stdexcept>
-#include <array>
-#include <iterator>
-#include <vector>
-#include <unordered_map>
-
-
-#include <limits>
-#include <random>
-#include <cstdlib>
-
-#include "manifold.h"
-#include "cross_section.h"
-
-#include "common.h"
-#include <cmath>
-#include <cassert>
-#include <iostream>
-#include <algorithm>
-
-namespace quickhull {
-
-//template<>
-//const double QuickHull<double>::Epsilon = 0.0000001;
-
-
-/*
- * Implementation of the algorithm
- */
-
-template<typename T>
-ConvexHull<T> QuickHull<T>::getConvexHull(const std::vector<Vector3<T>>& pointCloud, bool CCW, bool useOriginalIndices, T epsilon) {
-    VertexDataSource<T> vertexDataSource(pointCloud);
-    return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
-}
-
-template<typename T>
-ConvexHull<T> QuickHull<T>::getConvexHull(const Vector3<T>* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, T epsilon) {
-    VertexDataSource<T> vertexDataSource(vertexData,vertexCount);
-    return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
-}
-
-template<typename T>
-ConvexHull<T> QuickHull<T>::getConvexHull(const T* vertexData, size_t vertexCount, bool CCW, bool useOriginalIndices, T epsilon) {
-    VertexDataSource<T> vertexDataSource((const vec3*)vertexData,vertexCount);
-    return getConvexHull(vertexDataSource,CCW,useOriginalIndices,epsilon);
-}
-
-template<typename FloatType>
-HalfEdgeMesh<FloatType, size_t> QuickHull<FloatType>::getConvexHullAsMesh(const FloatType* vertexData, size_t vertexCount, bool CCW, FloatType epsilon) {
-    VertexDataSource<FloatType> vertexDataSource((const vec3*)vertexData,vertexCount);
-    buildMesh(vertexDataSource, CCW, false, epsilon);
-    return HalfEdgeMesh<FloatType, size_t>(m_mesh, m_vertexData);
-}
-
-template<typename T>
-void QuickHull<T>::buildMesh(const VertexDataSource<T>& pointCloud, bool UNUSED(CCW), bool UNUSED(useOriginalIndices), T epsilon) {
-    if (pointCloud.size()==0) {
-        m_mesh = MeshBuilder<T>();
-        return;
-    }
-    m_vertexData = pointCloud;
-
-    // Very first: find extreme values and use them to compute the scale of the point cloud.
-    m_extremeValues = getExtremeValues();
-    m_scale = getScale(m_extremeValues);
-
-    // Epsilon we use depends on the scale
-    m_epsilon = epsilon*m_scale;
-    m_epsilonSquared = m_epsilon*m_epsilon;
-
-    m_planar = false; // The planar case happens when all the points appear to lie on a two dimensional subspace of R^3.
-    createConvexHalfEdgeMesh();
-    if (m_planar) {
-        const size_t extraPointIndex = m_planarPointCloudTemp.size()-1;
-        for (auto& he : m_mesh.m_halfEdges) {
-            if (he.m_endVertex == extraPointIndex) {
-                he.m_endVertex = 0;
-            }
-        }
-        m_vertexData = pointCloud;
-        m_planarPointCloudTemp.clear();
-    }
-}
-
-template<typename T>
-ConvexHull<T> QuickHull<T>::getConvexHull(const VertexDataSource<T>& pointCloud, bool CCW, bool useOriginalIndices, T epsilon) {
-    buildMesh(pointCloud,CCW,useOriginalIndices,epsilon);
-    return ConvexHull<T>(m_mesh,m_vertexData, CCW, useOriginalIndices);
-}
-
-template<typename T>
-void QuickHull<T>::createConvexHalfEdgeMesh() {
-
-    m_visibleFaces.clear();
-    m_horizonEdges.clear();
-    m_possiblyVisibleFaces.clear();
-
-    // Compute base tetrahedron
-    setupInitialTetrahedron();
-
-    // Init face stack with those faces that have points assigned to them
-    m_faceList.clear();
-    for (size_t i=0;i < 4;i++) {
-        auto& f = m_mesh.m_faces[i];
-        if (f.m_pointsOnPositiveSide && f.m_pointsOnPositiveSide->size()>0) {
-            m_faceList.push_back(i);
-            f.m_inFaceStack = 1;
-        }
-    }
-
-    // Process faces until the face list is empty.
-    size_t iter = 0;
-    while (!m_faceList.empty()) {
-        iter++;
-        if (iter == std::numeric_limits<size_t>::max()) {
-            // Visible face traversal marks visited faces with iteration counter (to mark that the face has been visited on this iteration) and the max value represents unvisited faces. At this point we have to reset iteration counter. This shouldn't be an
-            // issue on 64 bit machines.
-            iter = 0;
-        }
-
-        const size_t topFaceIndex = m_faceList.front();
-        m_faceList.pop_front();
-
-        auto& tf = m_mesh.m_faces[topFaceIndex];
-        tf.m_inFaceStack = 0;
-
-        assert(!tf.m_pointsOnPositiveSide || tf.m_pointsOnPositiveSide->size()>0);
-        if (!tf.m_pointsOnPositiveSide || tf.isDisabled()) {
-            continue;
-        }
-
-        // Pick the most distant point to this triangle plane as the point to which we extrude
-        const vec3& activePoint = m_vertexData[tf.m_mostDistantPoint];
-        const size_t activePointIndex = tf.m_mostDistantPoint;
-
-        // Find out the faces that have our active point on their positive side (these are the "visible faces"). The face on top of the stack of course is one of them. At the same time, we create a list of horizon edges.
-        m_horizonEdges.clear();
-        m_possiblyVisibleFaces.clear();
-        m_visibleFaces.clear();
-        m_possiblyVisibleFaces.emplace_back(topFaceIndex,std::numeric_limits<size_t>::max());
-        while (m_possiblyVisibleFaces.size()) {
-            const auto faceData = m_possiblyVisibleFaces.back();
-            m_possiblyVisibleFaces.pop_back();
-            auto& pvf = m_mesh.m_faces[faceData.m_faceIndex];
-            assert(!pvf.isDisabled());
-
-            if (pvf.m_visibilityCheckedOnIteration == iter) {
-                if (pvf.m_isVisibleFaceOnCurrentIteration) {
-                    continue;
-                }
-            }
-            else {
-                const Plane<T>& P = pvf.m_P;
-                pvf.m_visibilityCheckedOnIteration = iter;
-                const T d = P.m_N.dotProduct(activePoint)+P.m_D;
-                if (d>0) {
-                    pvf.m_isVisibleFaceOnCurrentIteration = 1;
-                    pvf.m_horizonEdgesOnCurrentIteration = 0;
-                    m_visibleFaces.push_back(faceData.m_faceIndex);
-                    for (auto heIndex : m_mesh.getHalfEdgeIndicesOfFace(pvf)) {
-                        if (m_mesh.m_halfEdges[heIndex].m_opp != faceData.m_enteredFromHalfEdge) {
-                            m_possiblyVisibleFaces.emplace_back( m_mesh.m_halfEdges[m_mesh.m_halfEdges[heIndex].m_opp].m_face,heIndex );
-                        }
-                    }
-                    continue;
-                }
-                assert(faceData.m_faceIndex != topFaceIndex);
-            }
-
-            // The face is not visible. Therefore, the halfedge we came from is part of the horizon edge.
-            pvf.m_isVisibleFaceOnCurrentIteration = 0;
-            m_horizonEdges.push_back(faceData.m_enteredFromHalfEdge);
-            // Store which half edge is the horizon edge. The other half edges of the face will not be part of the final mesh so their data slots can by recycled.
-            const auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace(m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face]);
-            const std::int8_t ind = (halfEdges[0]==faceData.m_enteredFromHalfEdge) ? 0 : (halfEdges[1]==faceData.m_enteredFromHalfEdge ? 1 : 2);
-            m_mesh.m_faces[m_mesh.m_halfEdges[faceData.m_enteredFromHalfEdge].m_face].m_horizonEdgesOnCurrentIteration |= (1<<ind);
-        }
-        const size_t horizonEdgeCount = m_horizonEdges.size();
-
-        // Order horizon edges so that they form a loop. This may fail due to numerical instability in which case we give up trying to solve horizon edge for this point and accept a minor degeneration in the convex hull.
-        if (!reorderHorizonEdges(m_horizonEdges)) {
-            std::cerr << "Failed to solve horizon edge." << std::endl;
-            auto it = std::find(tf.m_pointsOnPositiveSide->begin(),tf.m_pointsOnPositiveSide->end(),activePointIndex);
-            tf.m_pointsOnPositiveSide->erase(it);
-            if (tf.m_pointsOnPositiveSide->size()==0) {
-                reclaimToIndexVectorPool(tf.m_pointsOnPositiveSide);
-            }
-            continue;
-        }
-
-        // Except for the horizon edges, all half edges of the visible faces can be marked as disabled. Their data slots will be reused.
-        // The faces will be disabled as well, but we need to remember the points that were on the positive side of them - therefore
-        // we save pointers to them.
-        m_newFaceIndices.clear();
-        m_newHalfEdgeIndices.clear();
-        m_disabledFacePointVectors.clear();
-        size_t disableCounter = 0;
-        for (auto faceIndex : m_visibleFaces) {
-            auto& disabledFace = m_mesh.m_faces[faceIndex];
-            auto halfEdges = m_mesh.getHalfEdgeIndicesOfFace(disabledFace);
-            for (size_t j=0;j<3;j++) {
-                if ((disabledFace.m_horizonEdgesOnCurrentIteration & (1<<j)) == 0) {
-                    if (disableCounter < horizonEdgeCount*2) {
-                        // Use on this iteration
-                        m_newHalfEdgeIndices.push_back(halfEdges[j]);
-                        disableCounter++;
-                    }
-                    else {
-                        // Mark for reusal on later iteration step
-                        m_mesh.disableHalfEdge(halfEdges[j]);
-                    }
-                }
-            }
-            // Disable the face, but retain pointer to the points that were on the positive side of it. We need to assign those points
-            // to the new faces we create shortly.
-            auto t = m_mesh.disableFace(faceIndex);
-            if (t) {
-                assert(t->size()); // Because we should not assign point vectors to faces unless needed...
-                m_disabledFacePointVectors.push_back(std::move(t));
-            }
-        }
-        if (disableCounter < horizonEdgeCount*2) {
-            const size_t newHalfEdgesNeeded = horizonEdgeCount*2-disableCounter;
-            for (size_t i=0;i<newHalfEdgesNeeded;i++) {
-                m_newHalfEdgeIndices.push_back(m_mesh.addHalfEdge());
-            }
-        }
-
-        // Create new faces using the edgeloop
-        for (size_t i = 0; i < horizonEdgeCount; i++) {
-            const size_t AB = m_horizonEdges[i];
-
-            auto horizonEdgeVertexIndices = m_mesh.getVertexIndicesOfHalfEdge(m_mesh.m_halfEdges[AB]);
-            size_t A,B,C;
-            A = horizonEdgeVertexIndices[0];
-            B = horizonEdgeVertexIndices[1];
-            C = activePointIndex;
-
-            const size_t newFaceIndex = m_mesh.addFace();
-            m_newFaceIndices.push_back(newFaceIndex);
-
-            const size_t CA = m_newHalfEdgeIndices[2*i+0];
-            const size_t BC = m_newHalfEdgeIndices[2*i+1];
-
-            m_mesh.m_halfEdges[AB].m_next = BC;
-            m_mesh.m_halfEdges[BC].m_next = CA;
-            m_mesh.m_halfEdges[CA].m_next = AB;
-
-            m_mesh.m_halfEdges[BC].m_face = newFaceIndex;
-            m_mesh.m_halfEdges[CA].m_face = newFaceIndex;
-            m_mesh.m_halfEdges[AB].m_face = newFaceIndex;
-
-            m_mesh.m_halfEdges[CA].m_endVertex = A;
-            m_mesh.m_halfEdges[BC].m_endVertex = C;
-
-            auto& newFace = m_mesh.m_faces[newFaceIndex];
-
-            const Vector3<T> planeNormal = mathutils::getTriangleNormal(m_vertexData[A],m_vertexData[B],activePoint);
-            newFace.m_P = Plane<T>(planeNormal,activePoint);
-            newFace.m_he = AB;
-
-            m_mesh.m_halfEdges[CA].m_opp = m_newHalfEdgeIndices[i>0 ? i*2-1 : 2*horizonEdgeCount-1];
-            m_mesh.m_halfEdges[BC].m_opp = m_newHalfEdgeIndices[((i+1)*2) % (horizonEdgeCount*2)];
-        }
-
-        // Assign points that were on the positive side of the disabled faces to the new faces.
-        for (auto& disabledPoints : m_disabledFacePointVectors) {
-            assert(disabledPoints);
-            for (const auto& point : *(disabledPoints)) {
-                if (point == activePointIndex) {
-                    continue;
-                }
-                for (size_t j=0;j<horizonEdgeCount;j++) {
-                    if (addPointToFace(m_mesh.m_faces[m_newFaceIndices[j]], point)) {
-                        break;
-                    }
-                }
-            }
-            // The points are no longer needed: we can move them to the vector pool for reuse.
-            reclaimToIndexVectorPool(disabledPoints);
-        }
-
-        // Increase face stack size if needed
-        for (const auto newFaceIndex : m_newFaceIndices) {
-            auto& newFace = m_mesh.m_faces[newFaceIndex];
-            if (newFace.m_pointsOnPositiveSide) {
-                assert(newFace.m_pointsOnPositiveSide->size()>0);
-                if (!newFace.m_inFaceStack) {
-                    m_faceList.push_back(newFaceIndex);
-                    newFace.m_inFaceStack = 1;
-                }
-            }
-        }
-    }
-
-    // Cleanup
-    m_indexVectorPool.clear();
-}
-
-/*
- * Private helper functions
- */
-
-template <typename T>
-std::array<size_t,6> QuickHull<T>::getExtremeValues() {
-    std::array<size_t,6> outIndices{{0,0,0,0,0,0}};
-    T extremeVals[6] = {m_vertexData[0].x,m_vertexData[0].x,m_vertexData[0].y,m_vertexData[0].y,m_vertexData[0].z,m_vertexData[0].z};
-    const size_t vCount = m_vertexData.size();
-    for (size_t i=1;i<vCount;i++) {
-        const Vector3<T>& pos = m_vertexData[i];
-        if (pos.x>extremeVals[0]) {
-            extremeVals[0]=pos.x;
-            outIndices[0]=(size_t)i;
-        }
-        else if (pos.x<extremeVals[1]) {
-            extremeVals[1]=pos.x;
-            outIndices[1]=(size_t)i;
-        }
-        if (pos.y>extremeVals[2]) {
-            extremeVals[2]=pos.y;
-            outIndices[2]=(size_t)i;
-        }
-        else if (pos.y<extremeVals[3]) {
-            extremeVals[3]=pos.y;
-            outIndices[3]=(size_t)i;
-        }
-        if (pos.z>extremeVals[4]) {
-            extremeVals[4]=pos.z;
-            outIndices[4]=(size_t)i;
-        }
-        else if (pos.z<extremeVals[5]) {
-            extremeVals[5]=pos.z;
-            outIndices[5]=(size_t)i;
-        }
-    }
-    return outIndices;
-}
-
-template<typename T>
-bool QuickHull<T>::reorderHorizonEdges(std::vector<size_t>& horizonEdges) {
-    const size_t horizonEdgeCount = horizonEdges.size();
-    for (size_t i=0;i<horizonEdgeCount-1;i++) {
-        const size_t endVertex = m_mesh.m_halfEdges[ horizonEdges[i] ].m_endVertex;
-        bool foundNext = false;
-        for (size_t j=i+1;j<horizonEdgeCount;j++) {
-            const size_t beginVertex = m_mesh.m_halfEdges[ m_mesh.m_halfEdges[horizonEdges[j]].m_opp ].m_endVertex;
-            if (beginVertex == endVertex) {
-                std::swap(horizonEdges[i+1],horizonEdges[j]);
-                foundNext = true;
-                break;
-            }
-        }
-        if (!foundNext) {
-            return false;
-        }
-    }
-    assert(m_mesh.m_halfEdges[ horizonEdges[horizonEdges.size()-1] ].m_endVertex == m_mesh.m_halfEdges[ m_mesh.m_halfEdges[horizonEdges[0]].m_opp ].m_endVertex);
-    return true;
-}
-
-template <typename T>
-T QuickHull<T>::getScale(const std::array<size_t,6>& extremeValues) {
-    T s = 0;
-    for (size_t i=0;i<6;i++) {
-        const T* v = (const T*)(&m_vertexData[extremeValues[i]]);
-        v += i/2;
-        auto a = std::abs(*v);
-        if (a>s) {
-            s = a;
-        }
-    }
-    return s;
-}
-
-template<typename T>
-void QuickHull<T>::setupInitialTetrahedron() {
-    const size_t vertexCount = m_vertexData.size();
-
-    // If we have at most 4 points, just return a degenerate tetrahedron:
-    if (vertexCount <= 4) {
-        size_t v[4] = {0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1)};
-        const Vector3<T> N = mathutils::getTriangleNormal(m_vertexData[v[0]],m_vertexData[v[1]],m_vertexData[v[2]]);
-        const Plane<T> trianglePlane(N,m_vertexData[v[0]]);
-        if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) {
-            std::swap(v[0],v[1]);
-        }
-        return m_mesh.setup(v[0],v[1],v[2],v[3]);
-    }
-
-    // Find two most distant extreme points.
-    T maxD = m_epsilonSquared;
-    std::pair<size_t,size_t> selectedPoints;
-    for (size_t i=0;i<6;i++) {
-        for (size_t j=i+1;j<6;j++) {
-            const T d = m_vertexData[ m_extremeValues[i] ].getSquaredDistanceTo( m_vertexData[ m_extremeValues[j] ] );
-            if (d > maxD) {
-                maxD=d;
-                selectedPoints={m_extremeValues[i],m_extremeValues[j]};
-            }
-        }
-    }
-    if (EQUAL(maxD, m_epsilonSquared)) {
-        // A degenerate case: the point cloud seems to consists of a single point
-        return m_mesh.setup(0,std::min((size_t)1,vertexCount-1),std::min((size_t)2,vertexCount-1),std::min((size_t)3,vertexCount-1));
-    }
-    assert(selectedPoints.first != selectedPoints.second);
-
-    // Find the most distant point to the line between the two chosen extreme points.
-    const Ray<T> r(m_vertexData[selectedPoints.first], (m_vertexData[selectedPoints.second] - m_vertexData[selectedPoints.first]));
-    maxD = m_epsilonSquared;
-    size_t maxI=std::numeric_limits<size_t>::max();
-    const size_t vCount = m_vertexData.size();
-    for (size_t i=0;i<vCount;i++) {
-        const T distToRay = mathutils::getSquaredDistanceBetweenPointAndRay(m_vertexData[i],r);
-        if (distToRay > maxD) {
-            maxD=distToRay;
-            maxI=i;
-        }
-    }
-    if (EQUAL(maxD, m_epsilonSquared)) {
-        // It appears that the point cloud belongs to a 1 dimensional
-        // subspace of R^3: convex hull has no volume => return a thin
-        // triangle Pick any point other than selectedPoints.first and
-        // selectedPoints.second as the third point of the triangle
-        auto it = std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) {
-            return ve != m_vertexData[selectedPoints.first] && ve != m_vertexData[selectedPoints.second];
-        });
-        const size_t thirdPoint = (it == m_vertexData.end()) ? selectedPoints.first : std::distance(m_vertexData.begin(),it);
-        it = std::find_if(m_vertexData.begin(),m_vertexData.end(),[&](const vec3& ve) {
-            return ve != m_vertexData[selectedPoints.first] && ve != m_vertexData[selectedPoints.second] && ve != m_vertexData[thirdPoint];
-        });
-        const size_t fourthPoint = (it == m_vertexData.end()) ? selectedPoints.first : std::distance(m_vertexData.begin(),it);
-        return m_mesh.setup(selectedPoints.first,selectedPoints.second,thirdPoint,fourthPoint);
-    }
-
-    // These three points form the base triangle for our tetrahedron.
-    assert(selectedPoints.first != maxI && selectedPoints.second != maxI);
-    std::array<size_t,3> baseTriangle{{selectedPoints.first, selectedPoints.second, maxI}};
-    const Vector3<T> baseTriangleVertices[]={ m_vertexData[baseTriangle[0]], m_vertexData[baseTriangle[1]],  m_vertexData[baseTriangle[2]] };
-
-    // Next step is to find the 4th vertex of the tetrahedron. We
-    // naturally choose the point farthest away from the triangle
-    // plane.
-    maxD=m_epsilon;
-    maxI=0;
-    const Vector3<T> N = mathutils::getTriangleNormal(baseTriangleVertices[0],baseTriangleVertices[1],baseTriangleVertices[2]);
-    Plane<T> trianglePlane(N,baseTriangleVertices[0]);
-    for (size_t i=0;i<vCount;i++) {
-        const T d = std::abs(mathutils::getSignedDistanceToPlane(m_vertexData[i],trianglePlane));
-        if (d>maxD) {
-            maxD=d;
-            maxI=i;
-        }
-    }
-    if (EQUAL(maxD, m_epsilon)) {
-        // All the points seem to lie on a 2D subspace of R^3. How to handle this? Well, let's add one extra point to the point cloud so that the convex hull will have volume.
-        m_planar = true;
-        const vec3 N1 = mathutils::getTriangleNormal(baseTriangleVertices[1],baseTriangleVertices[2],baseTriangleVertices[0]);
-        m_planarPointCloudTemp.clear();
-        m_planarPointCloudTemp.insert(m_planarPointCloudTemp.begin(),m_vertexData.begin(),m_vertexData.end());
-        const vec3 extraPoint = N1 + m_vertexData[0];
-        m_planarPointCloudTemp.push_back(extraPoint);
-        maxI = m_planarPointCloudTemp.size()-1;
-        m_vertexData = VertexDataSource<T>(m_planarPointCloudTemp);
-    }
-
-    // Enforce CCW orientation (if user prefers clockwise orientation, swap two vertices in each triangle when final mesh is created)
-    const Plane<T> triPlane(N,baseTriangleVertices[0]);
-    if (triPlane.isPointOnPositiveSide(m_vertexData[maxI])) {
-        std::swap(baseTriangle[0],baseTriangle[1]);
-    }
-
-    // Create a tetrahedron half edge mesh and compute planes defined by each triangle
-    m_mesh.setup(baseTriangle[0],baseTriangle[1],baseTriangle[2],maxI);
-    for (auto& f : m_mesh.m_faces) {
-        auto v = m_mesh.getVertexIndicesOfFace(f);
-        const Vector3<T>& va = m_vertexData[v[0]];
-        const Vector3<T>& vb = m_vertexData[v[1]];
-        const Vector3<T>& vc = m_vertexData[v[2]];
-        const Vector3<T> N1 = mathutils::getTriangleNormal(va, vb, vc);
-        const Plane<T> trianglePlane1(N1,va);
-        f.m_P = trianglePlane1;
-    }
-
-    // Finally we assign a face for each vertex outside the tetrahedron (vertices inside the tetrahedron have no role anymore)
-    for (size_t i=0;i<vCount;i++) {
-        for (auto& face : m_mesh.m_faces) {
-            if (addPointToFace(face, i)) {
-                break;
-            }
-        }
-    }
-}
-
-/*
- * Explicit template specifications for float and double
- */
-
-template class QuickHull<float>;
-template class QuickHull<double>;
-template<> const float QuickHull<float>::Epsilon = 0.0001f;
-template<> const double QuickHull<double>::Epsilon = 0.00001f;
-
-}
-
-// Local Variables:
-// tab-width: 8
-// mode: C++
-// c-basic-offset: 4
-// indent-tabs-mode: t
-// c-file-style: "stroustrup"
-// End:
-// ex: shiftwidth=4 tabstop=8
-
-
-namespace ConvexHull {
-
-struct vec3_hash {
-    size_t operator()(const std::array<float, 3>& vec) const {
-        std::hash<float> hash_fn;
-        return hash_fn(vec[0]) ^ hash_fn(vec[1]) ^ hash_fn(vec[2]);
-    }
-};
-
-struct vec2_hash {
-    size_t operator()(const std::array<float, 2>& vec) const {
-        std::hash<float> hash_fn;
-        return hash_fn(vec[0]) ^ hash_fn(vec[1]);
-    }
-};
-
-void QuickHull3DAntti(const std::vector<glm::vec3>& inputVerts, std::vector<glm::ivec3>& triVerts, std::vector<glm::vec3>& vertPos, const float precision) {
-    // Make sure inputVerts is not empty
-    if (inputVerts.empty()) {
-        return;
-    }
-
-    // Convert inputVerts to the format required by QuickHull
-    std::vector<quickhull::Vector3<float>> qhInputVerts;
-    for (const auto& v : inputVerts) {
-        qhInputVerts.emplace_back(v.x, v.y, v.z);
-    }
-
-    quickhull::QuickHull<float> qh;
-    auto hull = qh.getConvexHull(qhInputVerts, false, false);
-    auto indexBuffer = hull.getIndexBuffer();
-    auto vertexBuffer = hull.getVertexBuffer();
-
-    // Convert the result to the format required by triVerts and vertPos
-    for (const auto& v : vertexBuffer) {
-        vertPos.push_back(glm::vec3(v.x, v.y, v.z));
-    }
-
-    for (size_t i = 0; i < indexBuffer.size(); i += 3) {
-        triVerts.push_back(glm::ivec3(indexBuffer[i], indexBuffer[i + 1], indexBuffer[i + 2]));
-    }
-}
-
-void QuickHull3D(const std::vector<glm::vec3>& inputVerts, std::vector<glm::ivec3>& triVerts, std::vector<glm::vec3>& vertPos, const float precision) {
-
-    constexpr std::size_t dim = 3;
-    using PointType = std::array<float, dim>;
-    using Points = std::vector<PointType>;
-
-    Points points(inputVerts.size());
-
-    for (std::size_t i = 0; i < inputVerts.size(); ++i) {
-        auto pt = inputVerts[i];
-        points[i] = {pt.x, pt.y, pt.z};
-    }
-
-    const float eps = 0.0;
-    QuickHullTomilov::quick_hull<typename Points::const_iterator> quickhull{dim, eps};
-    quickhull.add_points(std::cbegin(points), std::cend(points));
-    auto initial_simplex = quickhull.get_affine_basis();
-    if (initial_simplex.size() < dim + 1) {
-        throw std::runtime_error("Convex Hull failled to find initial simplex!");
-    }
-
-    quickhull.create_initial_simplex(std::cbegin(initial_simplex), std::prev(std::cend(initial_simplex)));
-    quickhull.create_convex_hull();
-    if (!quickhull.check()) {
-        throw std::runtime_error("Convex hull check failed! (generally due to precision errors)");
-    }
-
-    std::unordered_map<std::array<float, 3>, int, vec3_hash> vertIndices;
-
-    for (const auto& facet : quickhull.facets_) {
-        for(const auto& vertex : facet.vertices_) {
-
-            auto vert = *vertex;
-            std::array<float, 3> arrVertex = {vert[0], vert[1], vert[2]};
-
-            if(vertIndices.count(arrVertex) == 0) {
-                vertIndices[arrVertex] = vertPos.size();
-                vertPos.push_back(glm::vec3(arrVertex[0], arrVertex[1], arrVertex[2]));
-            }
-        }
-    }
-
-    for (const auto& facet : quickhull.facets_) {
-
-        auto vert = *facet.vertices_[0];
-        int firstVertIndex = vertIndices[{vert[0], vert[1], vert[2]}];
-        for(std::size_t i = 1; i < facet.vertices_.size() - 1; ++i) {
-            auto currVert = *facet.vertices_[i];
-            auto nextVert = *facet.vertices_[i+1];
-            int secondVertIndex = vertIndices[{currVert[0], currVert[1], currVert[2]}];
-            int thirdVertIndex = vertIndices[{nextVert[0], nextVert[1], nextVert[2]}];
-            triVerts.push_back(glm::ivec3(firstVertIndex, secondVertIndex, thirdVertIndex));
-        }
-    }
-}
-
-int findMinYPointIndex(const std::vector<glm::vec2>& points) {
-    int minYPointIndex = 0;
-    for (size_t i = 1; i < points.size(); i++) {
-        if ((points[i].y < points[minYPointIndex].y) ||
-            (points[i].y == points[minYPointIndex].y && points[i].x > points[minYPointIndex].x)) {
-            minYPointIndex = i;
-        }
-    }
-    return minYPointIndex;
-}
-
-std::vector<glm::vec2> sortPointsCounterClockwise(const std::vector<glm::vec2>& points) {
-    std::vector<glm::vec2> sortedPoints(points);
-
-    // Find the bottom-most point (or one of them, if multiple)
-    int minYPointIndex = findMinYPointIndex(sortedPoints);
-
-    // Sort the points by angle from the line horizontal to minYPoint, counter-clockwise
-    glm::vec2 minYPoint = points[minYPointIndex];
-    std::sort(sortedPoints.begin(), sortedPoints.end(),
-        [minYPoint](const glm::vec2& p1, const glm::vec2& p2) -> bool {
-            double angle1 = atan2(p1.y - minYPoint.y, p1.x - minYPoint.x);
-            double angle2 = atan2(p2.y - minYPoint.y, p2.x - minYPoint.x);
-            if (angle1 < 0) angle1 += 2 * 3.141592653589;
-            if (angle2 < 0) angle2 += 2 * 3.141592653589;
-            return angle1 < angle2;
-        }
-    );
-
-    return sortedPoints;
-}
-
-manifold::SimplePolygon QuickHull2D(const manifold::SimplePolygon& inputVerts, const float precision) {
-
-    constexpr std::size_t dim = 2;
-    using PointType = std::array<float, dim>;
-    using Points = std::vector<PointType>;
-
-    Points points(inputVerts.size()); // input
-
-    for (std::size_t i = 0; i < inputVerts.size(); ++i) {
-        auto pt = inputVerts[i];
-        points[i] = {pt.x, pt.y};
-    }
-
-    QuickHullTomilov::quick_hull<typename Points::const_iterator> quickhull{dim, precision};
-    quickhull.add_points(std::cbegin(points), std::cend(points));
-    auto initial_simplex = quickhull.get_affine_basis();
-
-    quickhull.create_initial_simplex(std::cbegin(initial_simplex), std::prev(std::cend(initial_simplex)));
-    quickhull.create_convex_hull();
-
-    std::unordered_map<std::array<float, 2>, int, vec2_hash> vertIndices;
-    manifold::SimplePolygon ret;
-
-    for (const auto& facet : quickhull.facets_) {
-        for(const auto& vertex : facet.vertices_) {
-
-            auto vert = *vertex;
-            std::array<float, 2> arrVertex = {vert[0], vert[1]};
-
-            if(vertIndices.count(arrVertex) == 0) {
-                vertIndices[arrVertex] = ret.size();
-                ret.push_back(glm::vec3(arrVertex[0], arrVertex[1], arrVertex[2]));
-            }
-        }
-    }
-    return sortPointsCounterClockwise(ret);
-}
-
-manifold::Manifold ConvexHull(const manifold::Manifold manifold, const float precision = 0.0001) {
-    manifold::Mesh inputMesh = manifold.GetMesh();
-    manifold::Mesh outputMesh;
-    QuickHull3DAntti(inputMesh.vertPos, outputMesh.triVerts, outputMesh.vertPos, precision);
-    //orientMesh(outputMesh);
-    return manifold::Manifold(outputMesh);
-}
-
-manifold::Manifold ConvexHull(const manifold::Manifold manifold, const manifold::Manifold other, const float precision = 0.0001) {
-    manifold::Mesh inputMesh1 = manifold.GetMesh();
-    manifold::Mesh inputMesh2 = other.GetMesh();
-
-    // Combine vertices from input meshes
-    std::vector<glm::vec3> combinedVerts;
-    for (auto& vert: inputMesh1.vertPos) {
-        combinedVerts.push_back(vert);
-    }
-
-    for (auto& vert: inputMesh2.vertPos) {
-        combinedVerts.push_back(vert);
-    }
-
-    manifold::Mesh outputMesh;
-
-    QuickHull3DAntti(combinedVerts, outputMesh.triVerts, outputMesh.vertPos, precision);
-
-    //orientMesh(outputMesh);
-
-    return manifold::Manifold(outputMesh);
-}
-
-manifold::CrossSection ConvexHull(const manifold::CrossSection& cross_section, const float precision = 0.0001) {
-    manifold::SimplePolygon hullPoints;
-    for (auto& poly: cross_section.ToPolygons()) {
-        for (auto& pt: poly) {
-            hullPoints.push_back(glm::vec2(pt.x, pt.y));
-        }
-    }
-    manifold::SimplePolygon res = QuickHull2D(hullPoints, precision);
-    return manifold::CrossSection(res);
-}
-
-
-manifold::CrossSection ConvexHull(const manifold::CrossSection& cross_section, const manifold::CrossSection& other, const float precision = 0.0001) {
-    manifold::SimplePolygon hullPoints;
-
-    for (auto& poly: cross_section.ToPolygons()) {
-        for (auto& pt: poly) {
-            hullPoints.push_back(glm::vec2(pt.x, pt.y));
-        }
-    }
-
-    for (auto& poly: other.ToPolygons()) {
-        for (auto& pt: poly) {
-            hullPoints.push_back(glm::vec2(pt.x, pt.y));
-        }
-    }
-
-    manifold::SimplePolygon res = QuickHull2D(hullPoints, precision);
-    return manifold::CrossSection(res);
-}
-
-}
diff --git a/bindings/java/src/main/cpp/manifold3d/convex_hull_qhull_impl.hpp b/bindings/java/src/main/cpp/manifold3d/convex_hull_qhull_impl.hpp
deleted file mode 100644
index eeef4b6d5..000000000
--- a/bindings/java/src/main/cpp/manifold3d/convex_hull_qhull_impl.hpp
+++ /dev/null
@@ -1,216 +0,0 @@
-#pragma once
-
-extern "C" {
-#include "libqhull_r/libqhull_r.h"
-}
-
-#include <stdexcept>
-#include "public.h"
-#include <glm/glm.hpp>
-#include <vector>
-#include <queue>
-#include <map>
-#include <algorithm>
-
-namespace manifold {
-
-using namespace manifold;
-
-std::vector<std::vector<int>> buildAdjacency(const Mesh &mesh) {
-    std::map<std::pair<int, int>, std::vector<int>> edgeToFace;
-    std::vector<std::vector<int>> adjacency(mesh.triVerts.size());
-
-    for (int i = 0; i < mesh.triVerts.size(); ++i) {
-        glm::ivec3 face = mesh.triVerts[i];
-        for (int j = 0; j < 3; ++j) {
-            std::pair<int, int> edge(face[j], face[(j+1)%3]);
-            if (edge.first > edge.second) std::swap(edge.first, edge.second);
-
-            edgeToFace[edge].push_back(i);
-        }
-    }
-
-    for (const auto &pair : edgeToFace) {
-        for (int i = 0; i < pair.second.size(); ++i) {
-            for (int j = i+1; j < pair.second.size(); ++j) {
-                adjacency[pair.second[i]].push_back(pair.second[j]);
-                adjacency[pair.second[j]].push_back(pair.second[i]);
-            }
-        }
-    }
-
-    return adjacency;
-}
-
-void orientMesh(Mesh &mesh) {
-    std::vector<std::vector<int>> adjacency = buildAdjacency(mesh);
-    std::vector<bool> visited(mesh.triVerts.size(), false);
-    std::queue<int> q;
-
-    q.push(0);
-    visited[0] = true;
-
-    while (!q.empty()) {
-        int currentFaceIndex = q.front();
-        q.pop();
-        glm::ivec3 currentFace = mesh.triVerts[currentFaceIndex];
-
-        for (int neighbour : adjacency[currentFaceIndex]) {
-            if (visited[neighbour]) continue;
-
-            glm::ivec3 neighbourFace = mesh.triVerts[neighbour];
-
-            for (int i = 0; i < 3; ++i) {
-                if ((neighbourFace[i] == currentFace[0] && neighbourFace[(i+1)%3] == currentFace[1]) ||
-                    (neighbourFace[i] == currentFace[1] && neighbourFace[(i+1)%3] == currentFace[2]) ||
-                    (neighbourFace[i] == currentFace[2] && neighbourFace[(i+1)%3] == currentFace[0])) {
-                    std::swap(mesh.triVerts[neighbour][(i+1)%3], mesh.triVerts[neighbour][(i+2)%3]);
-                    break;
-                }
-            }
-
-            visited[neighbour] = true;
-            q.push(neighbour);
-        }
-    }
-}
-
-
-Mesh computeConvexHull3D(const std::vector<glm::vec3>& combinedVerts) {
-  // Convert vertices to coordinate array for Qhull
-  int dim = 3;
-  int n = combinedVerts.size();
-
-  coordT* points = new coordT[n*dim];
-  for(int i = 0; i < n; i++) {
-    points[i*dim] = combinedVerts[i].x;
-    points[i*dim+1] = combinedVerts[i].y;
-    points[i*dim+2] = combinedVerts[i].z;
-  }
-
-  boolT ismalloc = False;  // True if qhull should free points in qh_freeqhull() or reallocation
-  char flags[] = "qhull Qt";  // option flags for qhull, see qh_opt.htm
-
-  // create a new instance of Qhull
-  qhT qh_qh;
-  qh_zero(&qh_qh, stderr);
-
-  int exitcode = qh_new_qhull(&qh_qh, dim, n, points, ismalloc, flags, NULL, NULL);
-  if(exitcode != 0) {
-    //std::cout << "Convex Hull failled! Returning first mesh." << std::endl;
-    throw std::runtime_error("Convex Hull failled!");
-  }
-
-  // Create a new mesh for the convex hull
-  Mesh convexHull;
-
-  std::map<int, int> vertexIndexMap;
-  //// Iterate over the facets of the hull
-  qhT* qh = &qh_qh;
-  facetT *facet;
-  vertexT *vertex, **vertexp;
-  FORALLfacets {
-    glm::ivec3 tri;
-    int i = 0;
-
-    FOREACHvertex_(facet->vertices) {
-      int id = qh_pointid(&qh_qh, vertex->point);
-
-      // Check if the vertex is already added
-      if(vertexIndexMap.find(id) == vertexIndexMap.end()) {
-        convexHull.vertPos.push_back(combinedVerts[id]);
-        vertexIndexMap[id] = convexHull.vertPos.size() - 1;
-      }
-
-      // Add the index to the triangle
-      tri[i] = vertexIndexMap[id];
-      i++;
-    }
-
-    convexHull.triVerts.push_back(tri);
-  }
-
-  qh_freeqhull(&qh_qh, !qh_ALL);
-  delete[] points;
-
-  // Orient faces by right-hand rule.
-  orientMesh(convexHull);
-
-  return convexHull;
-}
-
-int findMinYPointIndex(const std::vector<glm::vec2>& points) {
-    int minYPointIndex = 0;
-    for (size_t i = 1; i < points.size(); i++) {
-        if ((points[i].y < points[minYPointIndex].y) ||
-            (points[i].y == points[minYPointIndex].y && points[i].x > points[minYPointIndex].x)) {
-            minYPointIndex = i;
-        }
-    }
-    return minYPointIndex;
-}
-
-std::vector<glm::vec2> sortPointsCounterClockwise(const std::vector<glm::vec2>& points) {
-    std::vector<glm::vec2> sortedPoints(points);
-
-    // Find the bottom-most point (or one of them, if multiple)
-    int minYPointIndex = findMinYPointIndex(sortedPoints);
-
-    // Sort the points by angle from the line horizontal to minYPoint, counter-clockwise
-    glm::vec2 minYPoint = points[minYPointIndex];
-    std::sort(sortedPoints.begin(), sortedPoints.end(),
-        [minYPoint](const glm::vec2& p1, const glm::vec2& p2) -> bool {
-            double angle1 = atan2(p1.y - minYPoint.y, p1.x - minYPoint.x);
-            double angle2 = atan2(p2.y - minYPoint.y, p2.x - minYPoint.x);
-            if (angle1 < 0) angle1 += 2 * 3.141592653589;
-            if (angle2 < 0) angle2 += 2 * 3.141592653589;
-            return angle1 < angle2;
-        }
-    );
-
-    return sortedPoints;
-}
-
-SimplePolygon computeConvexHull2D(const SimplePolygon& allPoints) {
-  // Convert points to coordinate array for Qhull
-  int dim = 2;  // We're now in 2D
-  int n = allPoints.size();
-
-  coordT* points = new coordT[n*dim];
-  for(int i = 0; i < n; i++) {
-    points[i*dim] = allPoints[i].x;
-    points[i*dim+1] = allPoints[i].y;
-  }
-
-  boolT ismalloc = False;
-  char flags[] = "qhull Qt";
-
-  // create a new instance of Qhull
-  qhT qh_qh;
-  qh_zero(&qh_qh, stderr);
-
-  int exitcode = qh_new_qhull(&qh_qh, dim, n, points, ismalloc, flags, NULL, NULL);
-  if(exitcode != 0) {
-    //std::cout << "Convex Hull failed! Returning first polygon." << std::endl;
-    return allPoints;
-  }
-
-  // Create a new polygon for the convex hull
-  SimplePolygon convexHull;
-
-  qhT* qh = &qh_qh;
-  vertexT *vertex, **vertexp;
-  FORALLvertices {
-    int id = qh_pointid(&qh_qh, vertex->point);
-    convexHull.push_back(allPoints[id]);
-  }
-
-  qh_freeqhull(&qh_qh, !qh_ALL);
-  delete[] points;
-
-  // It's not clear why qhull does not sanely order vertices. For now we sort them ourselves.
-  SimplePolygon sortedPoints = sortPointsCounterClockwise(convexHull);
-  return sortedPoints;
-}
-
-}
diff --git a/bindings/java/src/main/cpp/manifold3d/vmath.h b/bindings/java/src/main/cpp/manifold3d/vmath.h
deleted file mode 100644
index c73c78cf6..000000000
--- a/bindings/java/src/main/cpp/manifold3d/vmath.h
+++ /dev/null
@@ -1,2201 +0,0 @@
-/*                         V M A T H . H
- * BRL-CAD
- *
- * Copyright (c) 2004-2023 United States Government as represented by
- * the U.S. Army Research Laboratory.
- *
- * This library is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * version 2.1 as published by the Free Software Foundation.
- *
- * This library is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with this file; see the file named COPYING for more
- * information.
- */
-/** @addtogroup vmath */
-/** @{ */
-/** @file vmath.h
- *
- * @brief fundamental vector, matrix, quaternion math macros
- *
- * VMATH defines commonly needed macros for 2D/3D/4D math involving:
- *
- *   points (point2d_t, point_t, and hpoint_t),
- *   vectors (vect2d_t, vect_t, and hvect_t),
- *   quaternions (quat_t),
- *   planes (plane_t), and
- *   4x4 matrices (mat_t).
- *
- * By default, all floating point numbers are stored in arrays using
- * the 'fastf_t' type definition.  It should be manually typedef'd to
- * the "fastest" 64-bit floating point type available on the current
- * hardware with at least 64 bits of precision.  On 16 and 32 bit
- * machines, this is typically "double", but on 64 bit machines, it
- * could be "float".
- *
- * Matrix array elements have the following positions in the matrix:
- * @code
- *			|  0  1  2  3 |		| 0 |
- *	[ 0 1 2 3 ]	|  4  5  6  7 |		| 1 |
- *			|  8  9 10 11 |		| 2 |
- *			| 12 13 14 15 |		| 3 |
- *
- * preVector (vect_t) Matrix (mat_t) postVector (vect_t)
- * @endcode
- *
- * Note that while many people in the computer graphics field use
- * post-multiplication with row vectors (i.e., vector * matrix *
- * matrix ...) VMATH uses the more traditional representation of
- * column vectors (i.e., ... matrix * matrix * vector).  (The matrices
- * in these two representations are the transposes of each
- * other). Therefore, when transforming a vector by a matrix,
- * pre-multiplication is used, i.e.:
- *
- * view_vec = model2view_mat * model_vec
- *
- * Furthermore, additional transformations are multiplied on the left,
- * i.e.:
- *
- * @code
- * vec'  =  T1 * vec
- * vec'' =  T2 * T1 * vec  =  T2 * vec'
- * @endcode
- *
- * The most notable implication of this is the location of the "delta"
- * (translation) values in the matrix, i.e.:
- *
- * @code
- * x'   (R0  R1  R2 Dx) x
- * y' = (R4  R5  R6 Dy) * y
- * z'   (R8  R9 R10 Dz) z
- * w'   (0   0   0  1/s) w
- * @endcode
- *
- * Note -
- *@n  vect_t objects are 3-tuples
- *@n hvect_t objects are 4-tuples
- *
- * Most of these macros require that the result be in separate
- * storage, distinct from the input parameters, except where noted.
- *
- * IMPLEMENTOR NOTES
- *
- * When writing macros like this, it is very important that any
- * variables declared within a macro code blocks start with an
- * underscore in order to (hopefully) minimize any name conflicts with
- * user-provided parameters, such as _f in the following example:
- *
- * @code
- * #define ABC() do { double _f; do stuff; } while (0)
- * @endcode
- *
- * All of the macros that introduce a scope like the preceding
- * example are written as do { } while (0) loops in order to require
- * callers provide a trailing semicolon (e.g., ABC();).  This helps
- * preserve source code formatting.
- */
-
-#ifndef VMATH_H
-#define VMATH_H
-
-#include "common.h"
-
-/* needed for additional math defines on Windows when including math.h */
-#define _USE_MATH_DEFINES 1
-
-/* for sqrt(), sin(), cos(), rint(), M_PI, INFINITY (HUGE_VAL), and more */
-#include <math.h>
-
-/* for floating point tolerances and other math constants */
-#include <float.h>
-
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-
-#ifndef M_
-#  define M_		XXX /**< all with 36-digits of precision */
-#endif
-
-#ifndef M_1_2PI
-#  define M_1_2PI	0.159154943091895335768883763372514362  /**< 1/(2*pi) */
-#endif
-#ifndef M_1_PI
-#  define M_1_PI	0.318309886183790671537767526745028724  /**< 1/pi */
-#endif
-#ifndef M_2_PI
-#  define M_2_PI	0.636619772367581343075535053490057448  /**< 2/pi */
-#endif
-#ifndef M_2_SQRTPI
-#  define M_2_SQRTPI	1.12837916709551257389615890312154517   /**< 2/sqrt(pi) */
-#endif
-#ifndef M_E
-#  define M_E		2.71828182845904523536028747135266250   /**< e */
-#endif
-#ifndef M_EULER
-#  define M_EULER	0.577215664901532860606512090082402431  /**< Euler's constant */
-#endif
-#ifndef M_LOG2E
-#  define M_LOG2E	1.44269504088896340735992468100189214   /**< log_2(e) */
-#endif
-#ifndef M_LOG10E
-#  define M_LOG10E	0.434294481903251827651128918916605082  /**< log_10(e) */
-#endif
-#ifndef M_LN2
-#  define M_LN2		0.693147180559945309417232121458176568  /**< log_e(2) */
-#endif
-#ifndef M_LN10
-#  define M_LN10	2.30258509299404568401799145468436421   /**< log_e(10) */
-#endif
-#ifndef M_LNPI
-#  define M_LNPI	1.14472988584940017414342735135305871   /** log_e(pi) */
-#endif
-#ifndef M_PI
-#  define M_PI		3.14159265358979323846264338327950288   /**< pi */
-#endif
-#ifndef M_2PI
-#  define M_2PI		6.28318530717958647692528676655900576   /**< 2*pi */
-#endif
-#ifndef M_PI_2
-#  define M_PI_2	1.57079632679489661923132169163975144   /**< pi/2 */
-#endif
-#ifndef M_PI_3
-#  define M_PI_3	1.04719755119659774615421446109316763   /**< pi/3 */
-#endif
-#ifndef M_PI_4
-#  define M_PI_4	0.785398163397448309615660845819875721  /**< pi/4 */
-#endif
-#ifndef M_SQRT1_2
-#  define M_SQRT1_2	0.707106781186547524400844362104849039  /**< sqrt(1/2) */
-#endif
-#ifndef M_SQRT2
-#  define M_SQRT2	1.41421356237309504880168872420969808   /**< sqrt(2) */
-#endif
-#ifndef M_SQRT3
-#  define M_SQRT3	1.73205080756887729352744634150587237   /**< sqrt(3) */
-#endif
-#ifndef M_SQRTPI
-#  define M_SQRTPI	1.77245385090551602729816748334114518   /**< sqrt(pi) */
-#endif
-
-#ifndef DEG2RAD
-#  define DEG2RAD	0.0174532925199432957692369076848861271 /**< pi/180 */
-#endif
-#ifndef RAD2DEG
-#  define RAD2DEG	57.2957795130823208767981548141051703   /**< 180/pi */
-#endif
-
-
-/**
- * Definitions about limits of floating point representation
- * Eventually, should be tied to type of hardware (IEEE, IBM, Cray)
- * used to implement the fastf_t type.
- *
- * MAX_FASTF - Very close to the largest value that can be held by a
- * fastf_t without overflow.  Typically specified as an integer power
- * of ten, to make the value easy to spot when printed.  TODO: macro
- * function syntax instead of constant (DEPRECATED)
- *
- * SQRT_MAX_FASTF - sqrt(MAX_FASTF), or slightly smaller.  Any number
- * larger than this, if squared, can be expected to * produce an
- * overflow.  TODO: macro function syntax instead of constant
- * (DEPRECATED)
- *
- * SMALL_FASTF - Very close to the smallest value that can be
- * represented while still being greater than zero.  Any number
- * smaller than this (and non-negative) can be considered to be
- * zero; dividing by such a number can be expected to produce a
- * divide-by-zero error.  All divisors should be checked against
- * this value before actual division is performed.  TODO: macro
- * function syntax instead of constant (DEPRECATED)
- *
- * SQRT_SMALL_FASTF - sqrt(SMALL_FASTF), or slightly larger.  The
- * value of this is quite a lot larger than that of SMALL_FASTF.  Any
- * number smaller than this, when squared, can be expected to produce
- * a zero result.  TODO: macro function syntax instead of constant
- * (DEPRECATED)
- *
- */
-#if defined(vax)
-/* DEC VAX "D" format, the most restrictive */
-#  define MAX_FASTF		1.0e37	/* Very close to the largest number */
-#  define SQRT_MAX_FASTF	1.0e18	/* This squared just avoids overflow */
-#  define SMALL_FASTF		1.0e-37	/* Anything smaller is zero */
-#  define SQRT_SMALL_FASTF	1.0e-18	/* This squared gives zero */
-#else
-/* IBM format, being the next most restrictive format */
-#  define MAX_FASTF		1.0e73	/* Very close to the largest number */
-#  define SQRT_MAX_FASTF	1.0e36	/* This squared just avoids overflow */
-#  define SMALL_FASTF		1.0e-77	/* Anything smaller is zero */
-#  if defined(aux)
-#    define SQRT_SMALL_FASTF	1.0e-40 /* _doprnt error in libc */
-#  else
-#    define SQRT_SMALL_FASTF	1.0e-39	/* This squared gives zero */
-#  endif
-#endif
-
-/** DEPRECATED, do not use */
-#define SMALL SQRT_SMALL_FASTF
-
-/**
- * It is necessary to have a representation of 1.0/0.0 or log(0),
- * i.e., "infinity" that fits within the dynamic range of the machine
- * being used.  This constant places an upper bound on the size object
- * which can be represented in the model.  With IEEE 754 floating
- * point, this may print as 'inf' and is represented with all 1 bits
- * in the biased-exponent field and all 0 bits in the fraction with
- * the sign indicating positive (0) or negative (1) infinity.
- * However, we do not assume or rely on IEEE 754 floating point.
- */
-#ifndef INFINITY
-#  if defined(DBL_MAX)
-#    define INFINITY ((fastf_t)DBL_MAX)
-#  elif defined(HUGE_VAL)
-#    define INFINITY ((fastf_t)HUGE_VAL)
-#  elif defined(MAXDOUBLE)
-#    define INFINITY ((fastf_t)MAXDOUBLE)
-#  elif defined(HUGE)
-#    define INFINITY ((fastf_t)HUGE)
-/* fall back to a single-precision limit */
-#  elif defined(FLT_MAX)
-#    define INFINITY ((fastf_t)FLT_MAX)
-#  elif defined(HUGE_VALF)
-#    define INFINITY ((fastf_t)HUGE_VALF)
-#  elif defined(MAXFLOAT)
-#    define INFINITY ((fastf_t)MAXFLOAT)
-#  else
-     /* all else fails, just pick something big slightly under the
-      * 32-bit single-precision floating point limit for IEEE 754.
-      */
-#    define INFINITY ((fastf_t)1.0e38)
-#  endif
-#endif
-
-
-/* minimum computation tolerances */
-#ifdef vax
-#  define VDIVIDE_TOL		(1.0e-10)
-#  define VUNITIZE_TOL		(1.0e-7)
-#else
-#  ifdef DBL_EPSILON
-#    define VDIVIDE_TOL		(DBL_EPSILON)
-#  else
-#    define VDIVIDE_TOL		(1.0e-20)
-#  endif
-#  ifdef FLT_EPSILON
-#    define VUNITIZE_TOL	(FLT_EPSILON)
-#  else
-#    define VUNITIZE_TOL	(1.0e-15)
-#  endif
-#endif
-
-
-/** @brief number of fastf_t's per vect2d_t */
-#define ELEMENTS_PER_VECT2D	2
-
-/** @brief number of fastf_t's per point2d_t */
-#define ELEMENTS_PER_POINT2D	2
-
-/** @brief number of fastf_t's per vect_t */
-#define ELEMENTS_PER_VECT	3
-
-/** @brief number of fastf_t's per point_t */
-#define ELEMENTS_PER_POINT	3
-
-/** @brief number of fastf_t's per hvect_t (homogeneous vector) */
-#define ELEMENTS_PER_HVECT	4
-
-/** @brief number of fastf_t's per hpt_t (homogeneous point) */
-#define ELEMENTS_PER_HPOINT	4
-
-/** @brief number of fastf_t's per plane_t */
-#define ELEMENTS_PER_PLANE	4
-
-/** @brief number of fastf_t's per mat_t */
-#define ELEMENTS_PER_MAT	(ELEMENTS_PER_PLANE*ELEMENTS_PER_PLANE)
-
-
-/*
- * Fundamental types
- */
-
-/** @brief fastest 64-bit (or larger) floating point type */
-typedef double fastf_t;
-
-/** @brief 2-tuple vector */
-typedef fastf_t vect2d_t[ELEMENTS_PER_VECT2D];
-
-/** @brief pointer to a 2-tuple vector */
-typedef fastf_t *vect2dp_t;
-
-/** @brief 2-tuple point */
-typedef fastf_t point2d_t[ELEMENTS_PER_POINT2D];
-
-/** @brief pointer to a 2-tuple point */
-typedef fastf_t *point2dp_t;
-
-/** @brief 3-tuple vector */
-typedef fastf_t vect_t[ELEMENTS_PER_VECT];
-
-/** @brief pointer to a 3-tuple vector */
-typedef fastf_t *vectp_t;
-
-/** @brief 3-tuple point */
-typedef fastf_t point_t[ELEMENTS_PER_POINT];
-
-/** @brief pointer to a 3-tuple point */
-typedef fastf_t *pointp_t;
-
-/** @brief 4-tuple vector */
-typedef fastf_t hvect_t[ELEMENTS_PER_HVECT];
-
-/** @brief 4-element quaternion */
-typedef hvect_t quat_t;
-
-/** @brief 4-tuple point */
-typedef fastf_t hpoint_t[ELEMENTS_PER_HPOINT];
-
-/** @brief 4x4 matrix */
-typedef fastf_t mat_t[ELEMENTS_PER_MAT];
-
-/** @brief pointer to a 4x4 matrix */
-typedef fastf_t *matp_t;
-
-/**
- * @brief Definition of a plane equation
- *
- * A plane is defined by a unit-length outward pointing normal vector
- * (N), and the perpendicular (shortest) distance from the origin to
- * the plane (in element N[W]).
- *
- * The plane consists of all points P=(x, y, z) such that
- *@n VDOT(P, N) - N[W] == 0
- *@n that is,
- *@n N[X]*x + N[Y]*y + N[Z]*z - N[W] == 0
- *
- * The inside of the halfspace bounded by the plane
- * consists of all points P such that
- *@n VDOT(P, N) - N[W] <= 0
- *
- * A ray with direction D is classified w.r.t. the plane by
- *
- *@n VDOT(D, N) < 0 ray enters halfspace defined by plane
- *@n VDOT(D, N) == 0 ray is parallel to plane
- *@n VDOT(D, N) > 0 ray exits halfspace defined by plane
- */
-typedef fastf_t plane_t[ELEMENTS_PER_PLANE];
-
-/** Vector component names for homogeneous (4-tuple) vectors */
-typedef enum vmath_vector_component_ {
-    X = 0,
-    Y = 1,
-    Z = 2,
-    W = 3,
-    H = W
-} vmath_vector_component;
-
-/**
- * Locations of deltas (MD*) and scaling values (MS*) in a 4x4
- * Homogeneous Transform matrix
- */
-typedef enum vmath_matrix_component_ {
-    MSX = 0,
-    MDX = 3,
-    MSY = 5,
-    MDY = 7,
-    MSZ = 10,
-    MDZ = 11,
-    MSA = 15
-} vmath_matrix_component;
-
-/**
- * Evaluates truthfully whether a number is not within valid range of
- * INFINITY to -INFINITY exclusive (open set).
- */
-#define INVALID(n) (!((n) > -INFINITY && (n) < INFINITY))
-
-/**
- * Evaluates truthfully whether any components of a vector are not
- * within a valid range.
- */
-#define VINVALID(v) (INVALID((v)[X]) || INVALID((v)[Y]) || INVALID((v)[Z]))
-
-/**
- * Evaluates truthfully whether any components of a 2D vector are not
- * within a valid range.
- */
-#define V2INVALID(v) (INVALID((v)[X]) || INVALID((v)[Y]))
-
-/**
- * Evaluates truthfully whether any components of a 4D vector are not
- * within a valid range.
- */
-#define HINVALID(v) (INVALID((v)[X]) || INVALID((v)[Y]) || INVALID((v)[Z]) || INVALID((v)[W]))
-
-/**
- * Return truthfully whether a value is within a specified epsilon
- * distance from zero.
- */
-#ifdef KEITH_WANTS_THIS
-/* this is a proposed change to equality/zero testing.  prior behavior
- * evaluated as an open set.  this would change the behavior to that
- * of a closed set so that you can perform exact comparisons against
- * the tolerance and get a match.  examples that fail with the current
- * macro: tol=0.1; 1.1 == 1.0 or tol=0; 1==1
- *
- * these need to be tested carefully to make sure we pass ALL
- * regression and integration tests, which will require some
- * concerted effort to coordinate prior to a release.  first step is
- * to evaluate impact on performance and behavior of our tests.
- */
-#  define NEAR_ZERO(val, epsilon)	(!(((val) < -epsilon) || ((val) > epsilon)))
-#  define NEAR_ZERO(val, epsilon)	(!(((val) < -epsilon)) && !(((val) > epsilon)))
-#else
-#  define NEAR_ZERO(val, epsilon)	(((val) > -epsilon) && ((val) < epsilon))
-#endif
-
-/**
- * Return truthfully whether all elements of a given vector are within
- * a specified epsilon distance from zero.
- */
-#define VNEAR_ZERO(v, tol) \
-	(NEAR_ZERO(v[X], tol) \
-	 && NEAR_ZERO(v[Y], tol) \
-	 && NEAR_ZERO(v[Z], tol))
-
-/**
- * Test for all elements of `v' being smaller than `tol'.
- * Version for degree 2 vectors.
- */
-#define V2NEAR_ZERO(v, tol) (NEAR_ZERO(v[X], tol) && NEAR_ZERO(v[Y], tol))
-
-/**
- * Test for all elements of `v' being smaller than `tol'.
- * Version for degree 2 vectors.
- */
-#define HNEAR_ZERO(v, tol) \
-    (NEAR_ZERO(v[X], tol) \
-     && NEAR_ZERO(v[Y], tol) \
-     && NEAR_ZERO(v[Z], tol) \
-     && NEAR_ZERO(h[W], tol))
-
-
-/**
- * Return truthfully whether a value is within a minimum
- * representation tolerance from zero.
- */
-#define ZERO(_a) NEAR_ZERO((_a), SMALL_FASTF)
-
-/**
- * Return truthfully whether a vector is within a minimum
- * representation tolerance from zero.
- */
-#define VZERO(_a) VNEAR_ZERO((_a), SMALL_FASTF)
-
-/**
- * Return truthfully whether a 2d vector is within a minimum
- * representation tolerance from zero.
- */
-#define V2ZERO(_a) V2NEAR_ZERO((_a), SMALL_FASTF)
-
-/**
- * Return truthfully whether a homogenized 4-element vector is within
- * a minimum representation tolerance from zero.
- */
-#define HZERO(_a) HNEAR_ZERO((_a), SMALL_FASTF)
-
-
-/**
- * Return truthfully whether two values are within a specified epsilon
- * distance from each other.
- */
-#define NEAR_EQUAL(_a, _b, _tol) NEAR_ZERO((_a) - (_b), (_tol))
-
-/**
- * Return truthfully whether two 3D vectors are approximately equal,
- * within a specified absolute tolerance.
- */
-#define VNEAR_EQUAL(_a, _b, _tol) \
-    (NEAR_EQUAL((_a)[X], (_b)[X], (_tol)) \
-     && NEAR_EQUAL((_a)[Y], (_b)[Y], (_tol)) \
-     && NEAR_EQUAL((_a)[Z], (_b)[Z], (_tol)))
-
-/**
- * Return truthfully whether two 2D vectors are approximately equal,
- * within a specified absolute tolerance.
- */
-#define V2NEAR_EQUAL(a, b, tol) \
-    (NEAR_EQUAL((a)[X], (b)[X], tol) \
-     && NEAR_EQUAL((a)[Y], (b)[Y], tol))
-
-/**
- * Return truthfully whether two 4D vectors are approximately equal,
- * within a specified absolute tolerance.
- */
-#define HNEAR_EQUAL(_a, _b, _tol) \
-    (NEAR_EQUAL((_a)[X], (_b)[X], (_tol)) \
-     && NEAR_EQUAL((_a)[Y], (_b)[Y], (_tol)) \
-     && NEAR_EQUAL((_a)[Z], (_b)[Z], (_tol)) \
-     && NEAR_EQUAL((_a)[W], (_b)[W], (_tol)))
-
-/**
- * Return truthfully whether two values are within a minimum
- * representation tolerance from each other.
- */
-#define EQUAL(_a, _b) NEAR_EQUAL((_a), (_b), SMALL_FASTF)
-
-
-/**
- * Return truthfully whether two vectors are equal within a minimum
- * representation tolerance.
- */
-#define VEQUAL(_a, _b) VNEAR_EQUAL((_a), (_b), SMALL_FASTF)
-
-/**
- * @brief Return truthfully whether two 2D vectors are equal within
- * a minimum representation tolerance.
- */
-#define V2EQUAL(_a, _b) V2NEAR_EQUAL((_a), (_b), SMALL_FASTF)
-
-/**
- * @brief Return truthfully whether two higher degree vectors are
- * equal within a minimum representation tolerance.
- */
-#define HEQUAL(_a, _b)  HNEAR_EQUAL((_a), (_b), SMALL_FASTF)
-
-
-/** @brief Compute distance from a point to a plane. */
-#define DIST_PNT_PLANE(_pt, _pl) (VDOT(_pt, _pl) - (_pl)[W])
-
-/** @brief Compute distance between two points. */
-#define DIST_PNT_PNT_SQ(_a, _b) \
-	((_a)[X]-(_b)[X])*((_a)[X]-(_b)[X]) + \
-	((_a)[Y]-(_b)[Y])*((_a)[Y]-(_b)[Y]) + \
-	((_a)[Z]-(_b)[Z])*((_a)[Z]-(_b)[Z])
-#define DIST_PNT_PNT(_a, _b) sqrt(DIST_PNT_PNT_SQ(_a, _b))
-
-/** @brief Compute distance between two 2D points. */
-#define DIST_PNT2_PNT2_SQ(_a, _b) \
-	((_a)[X]-(_b)[X])*((_a)[X]-(_b)[X]) + \
-	((_a)[Y]-(_b)[Y])*((_a)[Y]-(_b)[Y])
-#define DIST_PNT2_PNT2(_a, _b) sqrt(DIST_PNT2_PNT2_SQ(_a, _b))
-
-/** @brief set translation values of 4x4 matrix with x, y, z values. */
-#define MAT_DELTAS(_m, _x, _y, _z) do { \
-	(_m)[MDX] = (_x); \
-	(_m)[MDY] = (_y); \
-	(_m)[MDZ] = (_z); \
-    } while (0)
-
-/** @brief set translation values of 4x4 matrix from a vector. */
-#define MAT_DELTAS_VEC(_m, _v)	\
-	MAT_DELTAS(_m, (_v)[X], (_v)[Y], (_v)[Z])
-
-/**
- * @brief set translation values of 4x4 matrix from a reversed
- * vector.
- */
-#define MAT_DELTAS_VEC_NEG(_m, _v)	\
-	MAT_DELTAS(_m, -(_v)[X], -(_v)[Y], -(_v)[Z])
-
-/** @brief get translation values of 4x4 matrix to a vector. */
-#define MAT_DELTAS_GET(_v, _m) do { \
-	(_v)[X] = (_m)[MDX]; \
-	(_v)[Y] = (_m)[MDY]; \
-	(_v)[Z] = (_m)[MDZ]; \
-    } while (0)
-
-/**
- * @brief get translation values of 4x4 matrix to a vector,
- * reversed.
- */
-#define MAT_DELTAS_GET_NEG(_v, _m) do { \
-	(_v)[X] = -(_m)[MDX]; \
-	(_v)[Y] = -(_m)[MDY]; \
-	(_v)[Z] = -(_m)[MDZ]; \
-    } while (0)
-
-/**
- * @brief increment translation elements in a 4x4 matrix with x, y, z
- * values.
- */
-#define MAT_DELTAS_ADD(_m, _x, _y, _z) do { \
-	(_m)[MDX] += (_x); \
-	(_m)[MDY] += (_y); \
-	(_m)[MDZ] += (_z); \
-    } while (0)
-
-/**
- * @brief increment translation elements in a 4x4 matrix from a
- * vector.
- */
-#define MAT_DELTAS_ADD_VEC(_m, _v) do { \
-	(_m)[MDX] += (_v)[X]; \
-	(_m)[MDY] += (_v)[Y]; \
-	(_m)[MDZ] += (_v)[Z]; \
-    } while (0)
-
-/**
- * @brief decrement translation elements in a 4x4 matrix with x, y, z
- * values.
- */
-#define MAT_DELTAS_SUB(_m, _x, _y, _z) do { \
-	(_m)[MDX] -= (_x); \
-	(_m)[MDY] -= (_y); \
-	(_m)[MDZ] -= (_z); \
-    } while (0)
-
-/**
- * @brief decrement translation elements in a 4x4 matrix from a
- * vector.
- */
-#define MAT_DELTAS_SUB_VEC(_m, _v) do { \
-	(_m)[MDX] -= (_v)[X]; \
-	(_m)[MDY] -= (_v)[Y]; \
-	(_m)[MDZ] -= (_v)[Z]; \
-    } while (0)
-
-/**
- * @brief decrement translation elements in a 4x4 matrix with x, y, z
- * values.
- */
-#define MAT_DELTAS_MUL(_m, _x, _y, _z) do { \
-	(_m)[MDX] *= (_x); \
-	(_m)[MDY] *= (_y); \
-	(_m)[MDZ] *= (_z); \
-    } while (0)
-
-/**
- * @brief decrement translation elements in a 4x4 matrix from a
- * vector.
- */
-#define MAT_DELTAS_MUL_VEC(_m, _v) do { \
-	(_m)[MDX] *= (_v)[X]; \
-	(_m)[MDY] *= (_v)[Y]; \
-	(_m)[MDZ] *= (_v)[Z]; \
-    } while (0)
-
-/** @brief set scale of 4x4 matrix from xyz. */
-#define MAT_SCALE(_m, _x, _y, _z) do { \
-	(_m)[MSX] = _x; \
-	(_m)[MSY] = _y; \
-	(_m)[MSZ] = _z; \
-    } while (0)
-
-/** @brief set scale of 4x4 matrix from vector. */
-#define MAT_SCALE_VEC(_m, _v) do { \
-	(_m)[MSX] = (_v)[X]; \
-	(_m)[MSY] = (_v)[Y]; \
-	(_m)[MSZ] = (_v)[Z]; \
-    } while (0)
-
-/** @brief set uniform scale of 4x4 matrix from scalar. */
-#define MAT_SCALE_ALL(_m, _s) (_m)[MSA] = (_s)
-
-/** @brief add to scaling elements in a 4x4 matrix from xyz. */
-#define MAT_SCALE_ADD(_m, _x, _y, _z) do { \
-	(_m)[MSX] += _x; \
-	(_m)[MSY] += _y; \
-	(_m)[MSZ] += _z; \
-    } while (0)
-
-/** @brief add to scaling elements in a 4x4 matrix from vector. */
-#define MAT_SCALE_ADD_VEC(_m, _v) do { \
-	(_m)[MSX] += (_v)[X]; \
-	(_m)[MSY] += (_v)[Y]; \
-	(_m)[MSZ] += (_v)[Z]; \
-    } while (0)
-
-/** @brief subtract from scaling elements in a 4x4 matrix from xyz. */
-#define MAT_SCALE_SUB(_m, _x, _y, _z) do { \
-	(_m)[MSX] -= _x; \
-	(_m)[MSY] -= _y; \
-	(_m)[MSZ] -= _z; \
-    } while (0)
-
-/**
- * @brief subtract from scaling elements in a 4x4 matrix from
- * vector.
- */
-#define MAT_SCALE_SUB_VEC(_m, _v) do { \
-	(_m)[MSX] -= (_v)[X]; \
-	(_m)[MSY] -= (_v)[Y]; \
-	(_m)[MSZ] -= (_v)[Z]; \
-    } while (0)
-
-/** @brief multiply scaling elements in a 4x4 matrix from xyz. */
-#define MAT_SCALE_MUL(_m, _x, _y, _z) do { \
-	(_m)[MSX] *= _x; \
-	(_m)[MSY] *= _y; \
-	(_m)[MSZ] *= _z; \
-    } while (0)
-
-/** @brief multiply scaling elements in a 4x4 matrix from vector. */
-#define MAT_SCALE_MUL_VEC(_m, _v) do { \
-	(_m)[MSX] *= (_v)[X]; \
-	(_m)[MSY] *= (_v)[Y]; \
-	(_m)[MSZ] *= (_v)[Z]; \
-    } while (0)
-
-
-/**
- * In following are macro versions of librt/mat.c functions for when
- * speed really matters.
- */
-
-
-/** @brief Zero a matrix. */
-#define MAT_ZERO(m) do { \
-	(m)[0] = (m)[1] = (m)[2] = (m)[3] = \
-	(m)[4] = (m)[5] = (m)[6] = (m)[7] = \
-	(m)[8] = (m)[9] = (m)[10] = (m)[11] = \
-	(m)[12] = (m)[13] = (m)[14] = (m)[15] = 0.0; \
-    } while (0)
-
-/** @brief Set matrix to identity. */
-#define MAT_IDN(m) do { \
-	(m)[1] = (m)[2] = (m)[3] = (m)[4] = \
-	(m)[6] = (m)[7] = (m)[8] = (m)[9] = \
-	(m)[11] = (m)[12] = (m)[13] = (m)[14] = 0.0; \
-	(m)[0] = (m)[5] = (m)[10] = (m)[15] = 1.0; \
-    } while (0)
-
-/**
- * @brief set t to the transpose of matrix m
- *
- * NOTE: This implementation will not transpose in-place or
- * overlapping matrices (e.g., MAT_TRANSPOSE(m, m) will be wrong).
- */
-#define MAT_TRANSPOSE(t, m) do { \
-	(t)[0] = (m)[0]; \
-	(t)[4] = (m)[1]; \
-	(t)[8] = (m)[2]; \
-	(t)[12] = (m)[3]; \
-	(t)[1] = (m)[4]; \
-	(t)[5] = (m)[5]; \
-	(t)[9] = (m)[6]; \
-	(t)[13] = (m)[7]; \
-	(t)[2] = (m)[8]; \
-	(t)[6] = (m)[9]; \
-	(t)[10] = (m)[10]; \
-	(t)[14] = (m)[11]; \
-	(t)[3] = (m)[12]; \
-	(t)[7] = (m)[13]; \
-	(t)[11] = (m)[14]; \
-	(t)[15] = (m)[15]; \
-    } while (0)
-
-/** @brief Copy a matrix `m' into `c'. */
-#define MAT_COPY(c, m) do { \
-	(c)[0] = (m)[0]; \
-	(c)[1] = (m)[1]; \
-	(c)[2] = (m)[2]; \
-	(c)[3] = (m)[3]; \
-	(c)[4] = (m)[4]; \
-	(c)[5] = (m)[5]; \
-	(c)[6] = (m)[6]; \
-	(c)[7] = (m)[7]; \
-	(c)[8] = (m)[8]; \
-	(c)[9] = (m)[9]; \
-	(c)[10] = (m)[10]; \
-	(c)[11] = (m)[11]; \
-	(c)[12] = (m)[12]; \
-	(c)[13] = (m)[13]; \
-	(c)[14] = (m)[14]; \
-	(c)[15] = (m)[15]; \
-    } while (0)
-
-/** @brief Set 3D vector at `o' to have coordinates `a', `b', and `c'. */
-#define VSET(o, a, b, c) do { \
-	(o)[X] = (a); \
-	(o)[Y] = (b); \
-	(o)[Z] = (c); \
-    } while (0)
-
-/** @brief Set 2D vector at `o' to have coordinates `a' and `b'. */
-#define V2SET(o, a, b) do { \
-	(o)[X] = (a); \
-	(o)[Y] = (b); \
-    } while (0)
-
-/** @brief Set 4D vector at `o' to homogeneous coordinates `a', `b', `c', and `d'. */
-#define HSET(o, a, b, c, d) do { \
-	(o)[X] = (a); \
-	(o)[Y] = (b); \
-	(o)[Z] = (c); \
-	(o)[W] = (d); \
-    } while (0)
-
-
-/** @brief Set all elements of 3D vector to same scalar value. */
-#define VSETALL(v, s) do { \
-	(v)[X] = (v)[Y] = (v)[Z] = (s); \
-    } while (0)
-
-/** @brief Set 2D vector elements to same scalar value. */
-#define V2SETALL(v, s) do { \
-	(v)[X] = (v)[Y] = (s); \
-    } while (0)
-
-/** @brief Set 4D vector elements to same scalar value. */
-#define HSETALL(v, s) do { \
-	(v)[X] = (v)[Y] = (v)[Z] = (v)[W] = (s); \
-    } while (0)
-
-
-/** @brief Set all elements of N-vector to same scalar value. */
-#define VSETALLN(v, s, n) do { \
-	size_t _j; \
-	for (_j=0; _j < (size_t)(n); _j++) v[_j]=(s); \
-    } while (0)
-
-
-/** @brief Transfer 3D vector at `v' to vector at `o'. */
-#define VMOVE(o, v) do { \
-	(o)[X] = (v)[X]; \
-	(o)[Y] = (v)[Y]; \
-	(o)[Z] = (v)[Z]; \
-    } while (0)
-
-/** @brief Move a 2D vector at `v' to vector at `o'. */
-#define V2MOVE(o, v) do { \
-	(o)[X] = (v)[X]; \
-	(o)[Y] = (v)[Y]; \
-    } while (0)
-
-/** @brief Move a homogeneous 4-tuple at `v' to `o'. */
-#define HMOVE(o, v) do { \
-	(o)[X] = (v)[X]; \
-	(o)[Y] = (v)[Y]; \
-	(o)[Z] = (v)[Z]; \
-	(o)[W] = (v)[W]; \
-    } while (0)
-
-/** @brief Transfer vector of length `n' at `v' to vector at `o'. */
-#define VMOVEN(o, v, n) do { \
-	size_t _vmove; \
-	for (_vmove = 0; _vmove < (size_t)(n); _vmove++) { \
-	    (o)[_vmove] = (v)[_vmove]; \
-	} \
-    } while (0)
-
-
-/**
- * @brief Reverse the direction of 3D vector `v' and store it in `o'.
- *
- * NOTE: Reversing in place works (i.e., VREVERSE(v, v))
- */
-#define VREVERSE(o, v) do { \
-	(o)[X] = -(v)[X]; \
-	(o)[Y] = -(v)[Y]; \
-	(o)[Z] = -(v)[Z]; \
-    } while (0)
-
-/**
- * @brief Reverse the direction of 2D vector `v' and store it in `o'.
- *
- * NOTE: Reversing in place works (i.e., V2REVERSE(v, v))
- */
-#define V2REVERSE(o, v) do { \
-	(o)[X] = -(v)[X]; \
-	(o)[Y] = -(v)[Y]; \
-    } while (0)
-
-/**
- * @brief Same as VREVERSE, but for a 4-tuple.  Also useful on plane_t
- * objects.
- *
- * NOTE: Reversing in place works (i.e., HREVERSE(v, v))
- */
-#define HREVERSE(o, v) do { \
-	(o)[X] = -(v)[X]; \
-	(o)[Y] = -(v)[Y]; \
-	(o)[Z] = -(v)[Z]; \
-	(o)[W] = -(v)[W]; \
-    } while (0)
-
-/** @brief Add 3D vectors at `a' and `b', store result at `o'. */
-#define VADD2(o, a, b) do { \
-	(o)[X] = (a)[X] + (b)[X]; \
-	(o)[Y] = (a)[Y] + (b)[Y]; \
-	(o)[Z] = (a)[Z] + (b)[Z]; \
-    } while (0)
-
-/** @brief Add 2D vectors at `a' and `b', store result at `o'. */
-#define V2ADD2(o, a, b) do { \
-	(o)[X] = (a)[X] + (b)[X]; \
-	(o)[Y] = (a)[Y] + (b)[Y]; \
-    } while (0)
-
-/** @brief Add 4D vectors at `a' and `b', store result at `o'. */
-#define HADD2(o, a, b) do { \
-	(o)[X] = (a)[X] + (b)[X]; \
-	(o)[Y] = (a)[Y] + (b)[Y]; \
-	(o)[Z] = (a)[Z] + (b)[Z]; \
-	(o)[W] = (a)[W] + (b)[W]; \
-    } while (0)
-
-/**
- * @brief Add vectors of length `n' at `a' and `b', store result at
- * `o'.
- */
-#define VADD2N(o, a, b, n) do { \
-	size_t _vadd2; \
-	for (_vadd2 = 0; _vadd2 < (size_t)(n); _vadd2++) { \
-		(o)[_vadd2] = (a)[_vadd2] + (b)[_vadd2]; \
-	} \
-    } while (0)
-
-
-/**
- * @brief Subtract 3D vector at `b' from vector at `a', store result at
- * `o'.
- */
-#define VSUB2(o, a, b) do { \
-	(o)[X] = (a)[X] - (b)[X]; \
-	(o)[Y] = (a)[Y] - (b)[Y]; \
-	(o)[Z] = (a)[Z] - (b)[Z]; \
-    } while (0)
-
-/**
- * @brief Subtract 2D vector at `b' from vector at `a', store result at
- * `o'.
- */
-#define V2SUB2(o, a, b) do { \
-	(o)[X] = (a)[X] - (b)[X]; \
-	(o)[Y] = (a)[Y] - (b)[Y]; \
-    } while (0)
-
-/**
- * @brief Subtract 4D vector at `b' from vector at `a', store result at
- * `o'.
- */
-#define HSUB2(o, a, b) do { \
-	(o)[X] = (a)[X] - (b)[X]; \
-	(o)[Y] = (a)[Y] - (b)[Y]; \
-	(o)[Z] = (a)[Z] - (b)[Z]; \
-	(o)[W] = (a)[W] - (b)[W]; \
-    } while (0)
-
-/**
- * @brief Subtract `n' length vector at `b' from `n' length vector at
- * `a', store result at `o'.
- */
-#define VSUB2N(o, a, b, n) do { \
-	size_t _vsub2; \
-	for (_vsub2 = 0; _vsub2 < (size_t)(n); _vsub2++) { \
-		(o)[_vsub2] = (a)[_vsub2] - (b)[_vsub2]; \
-	} \
-    } while (0)
-
-
-/** @brief 3D Vectors:  O = A - B - C */
-#define VSUB3(o, a, b, c) do { \
-	(o)[X] = (a)[X] - (b)[X] - (c)[X]; \
-	(o)[Y] = (a)[Y] - (b)[Y] - (c)[Y]; \
-	(o)[Z] = (a)[Z] - (b)[Z] - (c)[Z]; \
-    } while (0)
-
-/** @brief 2D Vectors:  O = A - B - C */
-#define V2SUB3(o, a, b, c) do { \
-	(o)[X] = (a)[X] - (b)[X] - (c)[X]; \
-	(o)[Y] = (a)[Y] - (b)[Y] - (c)[Y]; \
-    } while (0)
-
-/** @brief 4D Vectors:  O = A - B - C */
-#define HSUB3(o, a, b, c) do { \
-	(o)[X] = (a)[X] - (b)[X] - (c)[X]; \
-	(o)[Y] = (a)[Y] - (b)[Y] - (c)[Y]; \
-	(o)[Z] = (a)[Z] - (b)[Z] - (c)[Z]; \
-	(o)[W] = (a)[W] - (b)[W] - (c)[W]; \
-    } while (0)
-
-/** @brief Vectors:  O = A - B - C for vectors of length `n'. */
-#define VSUB3N(o, a, b, c, n) do { \
-	size_t _vsub3; \
-	for (_vsub3 = 0; _vsub3 < (size_t)(n); _vsub3++) { \
-		(o)[_vsub3] = (a)[_vsub3] - (b)[_vsub3] - (c)[_vsub3]; \
-	} \
-    } while (0)
-
-
-/** @brief Add 3 3D vectors at `a', `b', and `c', store result at `o'. */
-#define VADD3(o, a, b, c) do { \
-	(o)[X] = (a)[X] + (b)[X] + (c)[X]; \
-	(o)[Y] = (a)[Y] + (b)[Y] + (c)[Y]; \
-	(o)[Z] = (a)[Z] + (b)[Z] + (c)[Z]; \
-    } while (0)
-
-/** @brief Add 3 2D vectors at `a', `b', and `c', store result at `o'. */
-#define V2ADD3(o, a, b, c) do { \
-	(o)[X] = (a)[X] + (b)[X] + (c)[X]; \
-	(o)[Y] = (a)[Y] + (b)[Y] + (c)[Y]; \
-    } while (0)
-
-/** @brief Add 3 4D vectors at `a', `b', and `c', store result at `o'. */
-#define HADD3(o, a, b, c) do { \
-	(o)[X] = (a)[X] + (b)[X] + (c)[X]; \
-	(o)[Y] = (a)[Y] + (b)[Y] + (c)[Y]; \
-	(o)[Z] = (a)[Z] + (b)[Z] + (c)[Z]; \
-	(o)[W] = (a)[W] + (b)[W] + (c)[W]; \
-    } while (0)
-
-/**
- * @brief Add 3 vectors of length `n' at `a', `b', and `c', store
- * result at `o'.
- */
-#define VADD3N(o, a, b, c, n) do { \
-	size_t _vadd3; \
-	for (_vadd3 = 0; _vadd3 < (size_t)(n); _vadd3++) { \
-		(o)[_vadd3] = (a)[_vadd3] + (b)[_vadd3] + (c)[_vadd3]; \
-	} \
-    } while (0)
-
-
-/**
- * @brief Add 4 vectors at `a', `b', `c', and `d', store result at
- * `o'.
- */
-#define VADD4(o, a, b, c, d) do { \
-	(o)[X] = (a)[X] + (b)[X] + (c)[X] + (d)[X]; \
-	(o)[Y] = (a)[Y] + (b)[Y] + (c)[Y] + (d)[Y]; \
-	(o)[Z] = (a)[Z] + (b)[Z] + (c)[Z] + (d)[Z]; \
-    } while (0)
-
-/**
- * @brief Add 4 2D vectors at `a', `b', `c', and `d', store result at
- * `o'.
- */
-#define V2ADD4(o, a, b, c, d) do { \
-	(o)[X] = (a)[X] + (b)[X] + (c)[X] + (d)[X]; \
-	(o)[Y] = (a)[Y] + (b)[Y] + (c)[Y] + (d)[Y]; \
-    } while (0)
-
-/**
- * @brief Add 4 4D vectors at `a', `b', `c', and `d', store result at
- * `o'.
- */
-#define HADD4(o, a, b, c, d) do { \
-	(o)[X] = (a)[X] + (b)[X] + (c)[X] + (d)[X]; \
-	(o)[Y] = (a)[Y] + (b)[Y] + (c)[Y] + (d)[Y]; \
-	(o)[Z] = (a)[Z] + (b)[Z] + (c)[Z] + (d)[Z]; \
-	(o)[W] = (a)[W] + (b)[W] + (c)[W] + (d)[W]; \
-    } while (0)
-
-/**
- * @brief Add 4 `n' length vectors at `a', `b', `c', and `d', store
- * result at `o'.
- */
-#define VADD4N(o, a, b, c, d, n) do { \
-	size_t _vadd4;		   \
-	for (_vadd4 = 0; _vadd4 < (size_t)(n); _vadd4++) { \
-		(o)[_vadd4] = (a)[_vadd4] + (b)[_vadd4] + (c)[_vadd4] + (d)[_vadd4]; \
-	} \
-    } while (0)
-
-
-/** @brief Scale 3D vector at `v' by scalar `s', store result at `o'. */
-#define VSCALE(o, v, s) do { \
-	(o)[X] = (v)[X] * (s); \
-	(o)[Y] = (v)[Y] * (s); \
-	(o)[Z] = (v)[Z] * (s); \
-    } while (0)
-
-/** @brief Scale 2D vector at `v' by scalar `s', store result at `o'. */
-#define V2SCALE(o, v, s) do { \
-	(o)[X] = (v)[X] * (s); \
-	(o)[Y] = (v)[Y] * (s); \
-    } while (0)
-
-/** @brief Scale 4D vector at `v' by scalar `s', store result at `o'. */
-#define HSCALE(o, v, s) do { \
-	(o)[X] = (v)[X] * (s); \
-	(o)[Y] = (v)[Y] * (s); \
-	(o)[Z] = (v)[Z] * (s); \
-	(o)[W] = (v)[W] * (s); \
-    } while (0)
-
-/**
- * @brief Scale vector of length `n' at `v' by scalar `s', store
- * result at `o'
- */
-#define VSCALEN(o, v, s, n) do { \
-	size_t _vscale; \
-	for (_vscale = 0; _vscale < (size_t)(n); _vscale++) { \
-		(o)[_vscale] = (v)[_vscale] * (s); \
-	} \
-    } while (0)
-
-/** @brief Normalize vector `v' to be a unit vector. */
-#define VUNITIZE(v) do { \
-	double _f = MAGSQ(v); \
-	if (! NEAR_EQUAL(_f, 1.0, VUNITIZE_TOL)) { \
-		_f = sqrt(_f); \
-		if (_f < VDIVIDE_TOL) { \
-			VSETALL((v), 0.0); \
-		} else { \
-			_f = 1.0/_f; \
-			(v)[X] *= _f; (v)[Y] *= _f; (v)[Z] *= _f; \
-		} \
-	} \
-    } while (0)
-
-/** @brief Normalize 2D vector `v' to be a unit vector. */
-#define V2UNITIZE(v) do { \
-	double _f = MAG2SQ(v); \
-	if (! NEAR_EQUAL(_f, 1.0, VUNITIZE_TOL)) { \
-		_f = sqrt(_f); \
-		if (_f < VDIVIDE_TOL) { \
-			V2SETALL((v), 0.0); \
-		} else { \
-			_f = 1.0/_f; \
-			(v)[X] *= _f; (v)[Y] *= _f; \
-		} \
-	} \
-    } while (0)
-
-/**
- * @brief Find the sum of two points, and scale the result.  Often
- * used to find the midpoint.
- */
-#define VADD2SCALE(o, a, b, s) do { \
-			(o)[X] = ((a)[X] + (b)[X]) * (s); \
-			(o)[Y] = ((a)[Y] + (b)[Y]) * (s); \
-			(o)[Z] = ((a)[Z] + (b)[Z]) * (s); \
-    } while (0)
-
-/**
- * @brief Find the sum of two vectors of length `n', and scale the
- * result by `s'.  Often used to find the midpoint.
- */
-#define VADD2SCALEN(o, a, b, s, n) do { \
-	size_t _vadd2scale; \
-	for (_vadd2scale = 0; \
-	     _vadd2scale < (size_t)(n); \
-	     _vadd2scale++) { \
-	    (o)[_vadd2scale] = ((a)[_vadd2scale] + (b)[_vadd2scale]) * (s); \
-	} \
-    } while (0)
-
-/**
- * @brief Find the difference between two points, and scale result.
- * Often used to compute bounding sphere radius given rpp points.
- */
-#define VSUB2SCALE(o, a, b, s) do { \
-			(o)[X] = ((a)[X] - (b)[X]) * (s); \
-			(o)[Y] = ((a)[Y] - (b)[Y]) * (s); \
-			(o)[Z] = ((a)[Z] - (b)[Z]) * (s); \
-    } while (0)
-
-/**
- * @brief Find the difference between two vectors of length `n', and
- * scale result by `s'.
- */
-#define VSUB2SCALEN(o, a, b, s, n) do { \
-	size_t _vsub2scale; \
-	for (_vsub2scale = 0; \
-	     _vsub2scale < (size_t)(n); \
-	     _vsub2scale++) { \
-	    (o)[_vsub2scale] = ((a)[_vsub2scale] - (b)[_vsub2scale]) * (s); \
-	} \
-    } while (0)
-
-
-/**
- * @brief Combine together three vectors, all scaled by scalars.
- *
- * DEPRECATED: API inconsistent, use combo of other macros.
- */
-#define VCOMB3(o, a, b, c, d, e, f) do { \
-	(o)[X] = (a) * (b)[X] + (c) * (d)[X] + (e) * (f)[X]; \
-	(o)[Y] = (a) * (b)[Y] + (c) * (d)[Y] + (e) * (f)[Y]; \
-	(o)[Z] = (a) * (b)[Z] + (c) * (d)[Z] + (e) * (f)[Z]; \
-    } while (0)
-
-/**
- * @brief Combine together three vectors of length `n', all scaled by
- * scalars.
- *
- * DEPRECATED: API inconsistent, use combo of other macros.
- */
-#define VCOMB3N(o, a, b, c, d, e, f, n) do { \
-	size_t _vcomb3; \
-	for (_vcomb3 = 0; \
-	     _vcomb3 < (size_t)(n); \
-	     _vcomb3++) { \
-	    (o)[_vcomb3] = (a) * (b)[_vcomb3] + (c) * (d)[_vcomb3] + (e) * (f)[_vcomb3]; \
-	} \
-    } while (0)
-
-/**
- * @brief Combine together 2 vectors, both scaled by scalars.
- */
-#define VCOMB2(o, sa, va, sb, vb) do { \
-	(o)[X] = (sa) * (va)[X] + (sb) * (vb)[X]; \
-	(o)[Y] = (sa) * (va)[Y] + (sb) * (vb)[Y]; \
-	(o)[Z] = (sa) * (va)[Z] + (sb) * (vb)[Z]; \
-    } while (0)
-
-/**
- * @brief Combine together 2 vectors of length `n', both scaled by
- * scalars.
- */
-#define VCOMB2N(o, sa, a, sb, b, n) do { \
-	size_t _vcomb2; \
-	for (_vcomb2 = 0; \
-	     _vcomb2 < (size_t)(n); \
-	     _vcomb2++) { \
-	    (o)[_vcomb2] = (sa) * (va)[_vcomb2] + (sb) * (vb)[_vcomb2]; \
-	} \
-    } while (0)
-
-/**
- * DEPRECATED.
- */
-#define VJOIN4(a, b, c, d, e, f, g, h, i, j) do { \
-	(a)[X] = (b)[X] + (c)*(d)[X] + (e)*(f)[X] + (g)*(h)[X] + (i)*(j)[X]; \
-	(a)[Y] = (b)[Y] + (c)*(d)[Y] + (e)*(f)[Y] + (g)*(h)[Y] + (i)*(j)[Y]; \
-	(a)[Z] = (b)[Z] + (c)*(d)[Z] + (e)*(f)[Z] + (g)*(h)[Z] + (i)*(j)[Z]; \
-    } while (0)
-
-/**
- * Join three scaled vectors to a base `a', storing the result in `o'.
- */
-#define VJOIN3(o, a, sb, b, sc, c, sd, d) do { \
-	(o)[X] = (a)[X] + (sb)*(b)[X] + (sc)*(c)[X] + (sd)*(d)[X]; \
-	(o)[Y] = (a)[Y] + (sb)*(b)[Y] + (sc)*(c)[Y] + (sd)*(d)[Y]; \
-	(o)[Z] = (a)[Z] + (sb)*(b)[Z] + (sc)*(c)[Z] + (sd)*(d)[Z]; \
-    } while (0)
-
-
-/**
- * @brief Compose 3D vector at `o' of:
- * Vector at `a' plus
- * scalar `sb' times vector at `b' plus
- * scalar `sc' times vector at `c'
- */
-#define VJOIN2(o, a, sb, b, sc, c) do { \
-	(o)[X] = (a)[X] + (sb) * (b)[X] + (sc) * (c)[X]; \
-	(o)[Y] = (a)[Y] + (sb) * (b)[Y] + (sc) * (c)[Y]; \
-	(o)[Z] = (a)[Z] + (sb) * (b)[Z] + (sc) * (c)[Z]; \
-    } while (0)
-
-/**
- * @brief Compose 2D vector at `o' of:
- * Vector at `a' plus
- * scalar `sb' times vector at `b' plus
- * scalar `sc' times vector at `c'
- */
-#define V2JOIN2(o, a, sb, b, sc, c) do { \
-	(o)[X] = (a)[X] + (sb) * (b)[X] + (sc) * (c)[X]; \
-	(o)[Y] = (a)[Y] + (sb) * (b)[Y] + (sc) * (c)[Y]; \
-    } while (0)
-
-/**
- * @brief Compose 4D vector at `o' of:
- * Vector at `a' plus
- * scalar `sb' times vector at `b' plus
- * scalar `sc' times vector at `c'
- */
-#define HJOIN2(o, a, sb, b, sc, c) do { \
-	(o)[X] = (a)[X] + (sb) * (b)[X] + (sc) * (c)[X]; \
-	(o)[Y] = (a)[Y] + (sb) * (b)[Y] + (sc) * (c)[Y]; \
-	(o)[Z] = (a)[Z] + (sb) * (b)[Z] + (sc) * (c)[Z]; \
-	(o)[W] = (a)[W] + (sb) * (b)[W] + (sc) * (c)[W]; \
-    } while (0)
-
-#define VJOIN2N(o, a, sb, b, sc, c, n) do { \
-	size_t _vjoin2; \
-	for (_vjoin2 = 0; \
-	     _vjoin2 < (size_t)(n); \
-	     _vjoin2++) { \
-	    (o)[_vjoin2] = (a)[_vjoin2] + (sb) * (b)[_vjoin2] + (sc) * (c)[_vjoin2]; \
-	} \
-    } while (0)
-
-
-/**
- * Compose 3D vector at `o' of:
- * vector at `a' plus
- * scalar `sb' times vector at `b'
- *
- * This is basically a shorthand for VSCALE();VADD2();.
- */
-#define VJOIN1(o, a, sb, b) do { \
-	(o)[X] = (a)[X] + (sb) * (b)[X]; \
-	(o)[Y] = (a)[Y] + (sb) * (b)[Y]; \
-	(o)[Z] = (a)[Z] + (sb) * (b)[Z]; \
-    } while (0)
-
-/**
- * Compose 2D vector at `o' of:
- * vector at `a' plus
- * scalar `sb' times vector at `b'
- *
- * This is basically a shorthand for V2SCALE();V2ADD2();.
- */
-#define V2JOIN1(o, a, sb, b) do { \
-	(o)[X] = (a)[X] + (sb) * (b)[X]; \
-	(o)[Y] = (a)[Y] + (sb) * (b)[Y]; \
-    } while (0)
-
-/**
- * Compose 4D vector at `o' of:
- * vector at `a' plus
- * scalar `sb' times vector at `b'
- *
- * This is basically a shorthand for HSCALE();HADD2();.
- */
-#define HJOIN1(o, a, sb, b) do { \
-	(o)[X] = (a)[X] + (sb) * (b)[X]; \
-	(o)[Y] = (a)[Y] + (sb) * (b)[Y]; \
-	(o)[Z] = (a)[Z] + (sb) * (b)[Z]; \
-	(o)[W] = (a)[W] + (sb) * (b)[W]; \
-    } while (0)
-
-/**
- * Compose `n'-D vector at `o' of:
- * vector at `a' plus
- * scalar `sb' times vector at `b'
- *
- * This is basically a shorthand for VSCALEN();VADD2N();.
- */
-#define VJOIN1N(o, a, sb, b, n) do { \
-	size_t _vjoin1; \
-	for (_vjoin1 = 0; \
-	     _vjoin1 < (size_t)(n); \
-	     _vjoin1++) { \
-	    (o)[_vjoin1] = (a)[_vjoin1] + (sb) * (b)[_vjoin1]; \
-	} \
-    } while (0)
-
-
-/**
- * @brief Blend into vector `o'
- * scalar `sa' times vector at `a' plus
- * scalar `sb' times vector at `b'
- */
-#define VBLEND2(o, sa, a, sb, b) do { \
-	(o)[X] = (sa) * (a)[X] + (sb) * (b)[X]; \
-	(o)[Y] = (sa) * (a)[Y] + (sb) * (b)[Y]; \
-	(o)[Z] = (sa) * (a)[Z] + (sb) * (b)[Z]; \
-    } while (0)
-
-/**
- * @brief Blend into vector `o'
- * scalar `sa' times vector at `a' plus
- * scalar `sb' times vector at `b'
- */
-#define VBLEND2N(o, sa, a, sb, b, n) do { \
-	size_t _vblend2; \
-	for (_vblend2 = 0; \
-	     _vblend2 < (size_t)(n); \
-	     _vblend2++) { \
-	    (b)[_vblend2] = (sa) * (a)[_vblend2] + (sb) * (b)[_vblend2]; \
-	} \
-    } while (0)
-
-
-/**
- * @brief Project vector `a' onto `b'
- * vector `c' is the component of `a' parallel to `b'
- *     "    `d' "   "     "      "   "  orthogonal "   "
- *
- * FIXME: consistency, the result should come first
- */
-#define VPROJECT(a, b, c, d) do { \
-    VSCALE(c, b, VDOT(a, b) / VDOT(b, b)); \
-    VSUB2(d, a, c); \
-    } while (0)
-
-/** @brief Return scalar magnitude squared of vector at `v' */
-#define MAGSQ(v)	((v)[X]*(v)[X] + (v)[Y]*(v)[Y] + (v)[Z]*(v)[Z])
-#define MAG2SQ(v)	((v)[X]*(v)[X] + (v)[Y]*(v)[Y])
-
-
-/**
- * @brief Return scalar magnitude of the 3D vector `a'.  This is
- * otherwise known as the Euclidean norm of the provided vector..
- */
-#define MAGNITUDE(v) sqrt(MAGSQ(v))
-
-/**
- * @brief Return scalar magnitude of the 2D vector at `a'.  This is
- * otherwise known as the Euclidean norm of the provided vector..
- */
-#define MAGNITUDE2(v) sqrt(MAG2SQ(v))
-
-/**
- * @brief Store cross product of 3D vectors at `a' and `b' in vector at `o'.
- *
- * NOTE: The "right hand rule" applies. If closing your right hand
- * goes from `a' to `b', then your thumb points in the direction of
- * the cross product.
- *
- * If the angle from `a' to `b' goes clockwise, then the result vector
- * points "into" the plane (inward normal).  Example: a=(0, 1, 0),
- * b=(1, 0, 0), then aXb=(0, 0, -1).
- *
- * If the angle from `a' to `b' goes counter-clockwise, then the
- * result vector points "out" of the plane.  This outward pointing
- * normal is the BRL-CAD convention.
- */
-#define VCROSS(o, a, b) do { \
-	(o)[X] = (a)[Y] * (b)[Z] - (a)[Z] * (b)[Y]; \
-	(o)[Y] = (a)[Z] * (b)[X] - (a)[X] * (b)[Z]; \
-	(o)[Z] = (a)[X] * (b)[Y] - (a)[Y] * (b)[X]; \
-    } while (0)
-
-/**
- * Return the analog of a cross product for 2D vectors `a' and `b'
- * as a scalar value.  If a = (ax, ay) and b = (bx, by), then the analog
- * of a x b is det(a*b) = ax*by - ay*bx
- */
-#define V2CROSS(a, b) ((a)[X] * (b)[Y] - (a)[Y] * (b)[X])
-
-/**
- * TODO: implement me
- */
-#define HCROSS(a, b, c)
-
-
-/** @brief Compute dot product of vectors at `a' and `b'. */
-#define VDOT(a, b)	((a)[X]*(b)[X] + (a)[Y]*(b)[Y] + (a)[Z]*(b)[Z])
-
-#define V2DOT(a, b)	((a)[X]*(b)[X] + (a)[Y]*(b)[Y])
-
-#define HDOT(a, b)	((a)[X]*(b)[X] + (a)[Y]*(b)[Y] + (a)[Z]*(b)[Z] + (a)[W]*(b)[W])
-
-
-/**
- * @brief Linearly interpolate between two 3D vectors `a' and `b' by
- * interpolant `t', expected in the range [0,1], with result in `o'.
- *
- * NOTE: We intentionally use the form "o = a*(1-t) + b*t" which is
- * mathematically equivalent to "o = a + (b-a)*t".  The latter might
- * result in fewer math operations but cannot guarantee o==v1 when
- * t==1 due to floating-point arithmetic error.
- */
-#define VLERP(o, a, b, t) do { \
-	(o)[X] = (a)[X] * (1 - (t)) + (b)[X] * (t); \
-	(o)[Y] = (a)[Y] * (1 - (t)) + (b)[Y] * (t); \
-	(o)[Z] = (a)[Z] * (1 - (t)) + (b)[Z] * (t); \
-    } while (0)
-
-/**
- * @brief Linearly interpolate between two 2D vectors `a' and `b' by
- * interpolant `t', expected in the range [0,1], with result in `o'.
- *
- * NOTE: We intentionally use the form "o = a*(1-t) + b*t" which is
- * mathematically equivalent to "o = a + (b-a)*t".  The latter might
- * result in fewer math operations but cannot guarantee o==v1 when
- * t==1 due to floating-point arithmetic error.
- */
-#define V2LERP(o, a, b, t) do { \
-	(o)[X] = (a)[X] * (1 - (t)) + (b)[X] * (t); \
-	(o)[Y] = (a)[Y] * (1 - (t)) + (b)[Y] * (t); \
-    } while (0)
-
-/**
- * @brief Linearly interpolate between two 4D vectors `a' and `b' by
- * interpolant `t', expected in the range [0,1], with result in `o'.
- *
- * NOTE: We intentionally use the form "o = a*(1-t) + b*t" which is
- * mathematically equivalent to "o = a + (b-a)*t".  The latter might
- * result in fewer math operations but cannot guarantee o==v1 when
- * t==1 due to floating-point arithmetic error.
- */
-#define HLERP(o, a, b, t) do { \
-	(o)[X] = (a)[X] * (1 - (t)) + (b)[X] * (t); \
-	(o)[Y] = (a)[Y] * (1 - (t)) + (b)[Y] * (t); \
-	(o)[Z] = (a)[Z] * (1 - (t)) + (b)[Z] * (t); \
-	(o)[W] = (a)[W] * (1 - (t)) + (b)[W] * (t); \
-    } while (0)
-
-
-/**
- * @brief Subtract two points to make a vector, dot with another
- * vector.  Returns the dot product scalar value.
- */
-#define VSUB2DOT(_pt2, _pt, _vec)	(\
-	((_pt2)[X] - (_pt)[X]) * (_vec)[X] + \
-	((_pt2)[Y] - (_pt)[Y]) * (_vec)[Y] + \
-	((_pt2)[Z] - (_pt)[Z]) * (_vec)[Z])
-
-/**
- * @brief Turn a vector into comma-separated list of elements, for
- * variable argument subroutines (e.g. printf()).
- */
-#define V2ARGS(a)	(a)[X], (a)[Y]
-#define V3ARGS(a)	(a)[X], (a)[Y], (a)[Z]
-#define V4ARGS(a)	(a)[X], (a)[Y], (a)[Z], (a)[W]
-
-/**
- * Clamp values within tolerance of an integer to that value.
- *
- * For example, INTCLAMP(10.0000123123) evaluates to 10.0
- *
- * NOTE: should use VDIVIDE_TOL here, but cannot yet.  we use
- * VUINITIZE_TOL until floats are replaced universally with fastf_t's
- * since their epsilon is considerably less than that of a double.
- */
-#define INTCLAMP(_a) (NEAR_EQUAL((_a), rint(_a), VUNITIZE_TOL) ? rint(_a) : (_a))
-
-/** Clamp a 3D vector's elements to nearby integer values. */
-#define VINTCLAMP(_v) do { \
-	(_v)[X] = INTCLAMP((_v)[X]); \
-	(_v)[Y] = INTCLAMP((_v)[Y]); \
-	(_v)[Z] = INTCLAMP((_v)[Z]); \
-    } while (0)
-
-/** Clamp a 2D vector's elements to nearby integer values. */
-#define V2INTCLAMP(_v) do { \
-	(_v)[X] = INTCLAMP((_v)[X]); \
-	(_v)[Y] = INTCLAMP((_v)[Y]); \
-    } while (0)
-
-/** Clamp a 4D vector's elements to nearby integer values. */
-#define HINTCLAMP(_v) do { \
-	VINTCLAMP(_v); \
-	(_v)[W] = INTCLAMP((_v)[W]); \
-    } while (0)
-
-
-/** @brief integer clamped versions of the previous arg macros. */
-#define V2INTCLAMPARGS(a) INTCLAMP((a)[X]), INTCLAMP((a)[Y])
-/** @brief integer clamped versions of the previous arg macros. */
-#define V3INTCLAMPARGS(a) INTCLAMP((a)[X]), INTCLAMP((a)[Y]), INTCLAMP((a)[Z])
-/** @brief integer clamped versions of the previous arg macros. */
-#define V4INTCLAMPARGS(a) INTCLAMP((a)[X]), INTCLAMP((a)[Y]), INTCLAMP((a)[Z]), INTCLAMP((a)[W])
-
-/** @brief Print vector name and components on stderr. */
-#define V2PRINT(a, b)	\
-	fprintf(stderr, "%s (%.6f, %.6g)\n", a, V2ARGS(b));
-#define VPRINT(a, b)	\
-	fprintf(stderr, "%s (%.6f, %.6f, %.6f)\n", a, V3ARGS(b));
-#define HPRINT(a, b)	\
-	fprintf(stderr, "%s (%.6f, %.6f, %.6f, %.6f)\n", a, V4ARGS(b));
-
-/**
- * @brief Included below are integer clamped versions of the previous
- * print macros.
- */
-
-#define V2INTCLAMPPRINT(a, b)	\
-	fprintf(stderr, "%s (%g, %g)\n", a, V2INTCLAMPARGS(b));
-#define VINTCLAMPPRINT(a, b)	\
-	fprintf(stderr, "%s (%g, %g, %g)\n", a, V3INTCLAMPARGS(b));
-#define HINTCLAMPPRINT(a, b)	\
-	fprintf(stderr, "%s (%g, %g, %g, %g)\n", a, V4INTCLAMPARGS(b));
-
-
-/** @brief Vector element multiplication.  Really: diagonal matrix X vect. */
-#define VELMUL(o, a, b) do { \
-	(o)[X] = (a)[X] * (b)[X]; \
-	(o)[Y] = (a)[Y] * (b)[Y]; \
-	(o)[Z] = (a)[Z] * (b)[Z]; \
-    } while (0)
-
-#define VELMUL3(o, a, b, c) do { \
-	(o)[X] = (a)[X] * (b)[X] * (c)[X]; \
-	(o)[Y] = (a)[Y] * (b)[Y] * (c)[Y]; \
-	(o)[Z] = (a)[Z] * (b)[Z] * (c)[Z]; \
-    } while (0)
-
-/** @brief Similar to VELMUL. */
-#define VELDIV(o, a, b) do { \
-	(o)[X] = (a)[X] / (b)[X]; \
-	(o)[Y] = (a)[Y] / (b)[Y]; \
-	(o)[Z] = (a)[Z] / (b)[Z]; \
-    } while (0)
-
-/**
- * @brief Given a direction vector, compute the inverses of each element.
- * When division by zero would have occurred, mark inverse as INFINITY.
- */
-#define VINVDIR(_inv, _dir) do { \
-	if ((_dir)[X] < -SQRT_SMALL_FASTF || (_dir)[X] > SQRT_SMALL_FASTF) { \
-		(_inv)[X]=1.0/(_dir)[X]; \
-	} else { \
-		(_dir)[X] = 0.0; \
-		(_inv)[X] = INFINITY; \
-	} \
-	if ((_dir)[Y] < -SQRT_SMALL_FASTF || (_dir)[Y] > SQRT_SMALL_FASTF) { \
-		(_inv)[Y]=1.0/(_dir)[Y]; \
-	} else { \
-		(_dir)[Y] = 0.0; \
-		(_inv)[Y] = INFINITY; \
-	} \
-	if ((_dir)[Z] < -SQRT_SMALL_FASTF || (_dir)[Z] > SQRT_SMALL_FASTF) { \
-		(_inv)[Z]=1.0/(_dir)[Z]; \
-	} else { \
-		(_dir)[Z] = 0.0; \
-		(_inv)[Z] = INFINITY; \
-	} \
-    } while (0)
-
-/**
- * @brief Apply the 3x3 part of a mat_t to a 3-tuple.  This rotates a
- * vector without scaling it (changing its length).
- */
-#define MAT3X3VEC(o, mat, vec) do { \
-	(o)[X] = (mat)[X]*(vec)[X]+(mat)[Y]*(vec)[Y] + (mat)[ 2]*(vec)[Z]; \
-	(o)[Y] = (mat)[4]*(vec)[X]+(mat)[5]*(vec)[Y] + (mat)[ 6]*(vec)[Z]; \
-	(o)[Z] = (mat)[8]*(vec)[X]+(mat)[9]*(vec)[Y] + (mat)[10]*(vec)[Z]; \
-    } while (0)
-
-/** @brief Multiply a 3-tuple by the 3x3 part of a mat_t. */
-#define VEC3X3MAT(o, i, m) do { \
-	(o)[X] = (i)[X]*(m)[X] + (i)[Y]*(m)[4] + (i)[Z]*(m)[8]; \
-	(o)[Y] = (i)[X]*(m)[1] + (i)[Y]*(m)[5] + (i)[Z]*(m)[9]; \
-	(o)[Z] = (i)[X]*(m)[2] + (i)[Y]*(m)[6] + (i)[Z]*(m)[10]; \
-    } while (0)
-
-/** @brief Apply the 3x3 part of a mat_t to a 2-tuple (Z part=0). */
-#define MAT3X2VEC(o, mat, vec) do { \
-	(o)[X] = (mat)[0]*(vec)[X] + (mat)[Y]*(vec)[Y]; \
-	(o)[Y] = (mat)[4]*(vec)[X] + (mat)[5]*(vec)[Y]; \
-	(o)[Z] = (mat)[8]*(vec)[X] + (mat)[9]*(vec)[Y]; \
-    } while (0)
-
-/** @brief Multiply a 2-tuple (Z=0) by the 3x3 part of a mat_t. */
-#define VEC2X3MAT(o, i, m) do { \
-	(o)[X] = (i)[X]*(m)[0] + (i)[Y]*(m)[4]; \
-	(o)[Y] = (i)[X]*(m)[1] + (i)[Y]*(m)[5]; \
-	(o)[Z] = (i)[X]*(m)[2] + (i)[Y]*(m)[6]; \
-    } while (0)
-
-/**
- * @brief Apply a 4x4 matrix to a 3-tuple which is an absolute Point
- * in space.  Output and input points should be separate arrays.
- */
-#define MAT4X3PNT(o, m, i) do { \
-	double _f; \
-	_f = 1.0/((m)[12]*(i)[X] + (m)[13]*(i)[Y] + (m)[14]*(i)[Z] + (m)[15]); \
-	(o)[X]=((m)[0]*(i)[X] + (m)[1]*(i)[Y] + (m)[ 2]*(i)[Z] + (m)[3]) * _f; \
-	(o)[Y]=((m)[4]*(i)[X] + (m)[5]*(i)[Y] + (m)[ 6]*(i)[Z] + (m)[7]) * _f; \
-	(o)[Z]=((m)[8]*(i)[X] + (m)[9]*(i)[Y] + (m)[10]*(i)[Z] + (m)[11])* _f; \
-    } while (0)
-
-/**
- * @brief Multiply an Absolute 3-Point by a full 4x4 matrix.  Output
- * and input points should be separate arrays.
- */
-#define PNT3X4MAT(o, i, m) do { \
-	double _f; \
-	_f = 1.0/((i)[X]*(m)[3] + (i)[Y]*(m)[7] + (i)[Z]*(m)[11] + (m)[15]); \
-	(o)[X]=((i)[X]*(m)[0] + (i)[Y]*(m)[4] + (i)[Z]*(m)[8] + (m)[12]) * _f; \
-	(o)[Y]=((i)[X]*(m)[1] + (i)[Y]*(m)[5] + (i)[Z]*(m)[9] + (m)[13]) * _f; \
-	(o)[Z]=((i)[X]*(m)[2] + (i)[Y]*(m)[6] + (i)[Z]*(m)[10] + (m)[14])* _f; \
-    } while (0)
-
-/**
- * @brief Multiply an Absolute hvect_t 4-Point by a full 4x4 matrix.
- * Output and input points should be separate arrays.
- */
-#define MAT4X4PNT(o, m, i) do { \
-	(o)[X]=(m)[ 0]*(i)[X] + (m)[ 1]*(i)[Y] + (m)[ 2]*(i)[Z] + (m)[ 3]*(i)[W]; \
-	(o)[Y]=(m)[ 4]*(i)[X] + (m)[ 5]*(i)[Y] + (m)[ 6]*(i)[Z] + (m)[ 7]*(i)[W]; \
-	(o)[Z]=(m)[ 8]*(i)[X] + (m)[ 9]*(i)[Y] + (m)[10]*(i)[Z] + (m)[11]*(i)[W]; \
-	(o)[W]=(m)[12]*(i)[X] + (m)[13]*(i)[Y] + (m)[14]*(i)[Z] + (m)[15]*(i)[W]; \
-    } while (0)
-
-/**
- * @brief Apply a 4x4 matrix to a 3-tuple which is a relative Vector
- * in space. This macro can scale the length of the vector if [15] !=
- * 1.0.  Output and input vectors should be separate arrays.
- */
-#define MAT4X3VEC(o, m, i) do { \
-	double _f; \
-	_f = 1.0/((m)[15]); \
-	(o)[X] = ((m)[0]*(i)[X] + (m)[1]*(i)[Y] + (m)[ 2]*(i)[Z]) * _f; \
-	(o)[Y] = ((m)[4]*(i)[X] + (m)[5]*(i)[Y] + (m)[ 6]*(i)[Z]) * _f; \
-	(o)[Z] = ((m)[8]*(i)[X] + (m)[9]*(i)[Y] + (m)[10]*(i)[Z]) * _f; \
-    } while (0)
-
-#define MAT4XSCALOR(o, m, i) do { \
-	(o) = (i) / (m)[15]; \
-    } while (0)
-
-/**
- * @brief Multiply a Relative 3-Vector by most of a 4x4 matrix.
- * Output and input vectors should be separate arrays.
- */
-#define VEC3X4MAT(o, i, m) do { \
-	double _f; \
-	_f = 1.0/((m)[15]); \
-	(o)[X] = ((i)[X]*(m)[0] + (i)[Y]*(m)[4] + (i)[Z]*(m)[8]) * _f; \
-	(o)[Y] = ((i)[X]*(m)[1] + (i)[Y]*(m)[5] + (i)[Z]*(m)[9]) * _f; \
-	(o)[Z] = ((i)[X]*(m)[2] + (i)[Y]*(m)[6] + (i)[Z]*(m)[10]) * _f; \
-    } while (0)
-
-/** @brief Multiply a Relative 2-Vector by most of a 4x4 matrix. */
-#define VEC2X4MAT(o, i, m) do { \
-	double _f; \
-	_f = 1.0/((m)[15]); \
-	(o)[X] = ((i)[X]*(m)[0] + (i)[Y]*(m)[4]) * _f; \
-	(o)[Y] = ((i)[X]*(m)[1] + (i)[Y]*(m)[5]) * _f; \
-	(o)[Z] = ((i)[X]*(m)[2] + (i)[Y]*(m)[6]) * _f; \
-    } while (0)
-
-/**
- * @brief Included below are macros to update min and max X, Y, Z
- * values to contain a point
- */
-
-#define V_MIN(r, s) if ((r) > (s)) r = (s)
-
-#define V_MAX(r, s) if ((r) < (s)) r = (s)
-
-#ifdef VMIN
-#  undef VMIN
-#endif
-#define VMIN(r, s) do { \
-	V_MIN((r)[X], (s)[X]); V_MIN((r)[Y], (s)[Y]); V_MIN((r)[Z], (s)[Z]); \
-    } while (0)
-
-#ifdef VMAX
-#  undef VMAX
-#endif
-#define VMAX(r, s) do { \
-	V_MAX((r)[X], (s)[X]); V_MAX((r)[Y], (s)[Y]); V_MAX((r)[Z], (s)[Z]); \
-    } while (0)
-
-#ifdef VMINMAX
-#  undef VMINMAX
-#endif
-#define VMINMAX(min, max, pt) do { \
-	VMIN((min), (pt)); VMAX((max), (pt)); \
-    } while (0)
-
-/**
- * @brief Included below are macros to update min and max X, Y
- * values to contain a point
- */
-
-#define V2MIN(r, s) do { \
-	V_MIN((r)[X], (s)[X]); V_MIN((r)[Y], (s)[Y]); \
-    } while (0)
-
-#define V2MAX(r, s) do { \
-	V_MAX((r)[X], (s)[X]); V_MAX((r)[Y], (s)[Y]); \
-    } while (0)
-
-#define V2MINMAX(min, max, pt) do { \
-	V2MIN((min), (pt)); V2MAX((max), (pt)); \
-    } while (0)
-
-/**
- * clamp a value to a low/high number.
- */
-#define CLAMP(_v, _l, _h) V_MAX((_v), (_l)); else V_MIN((_v), (_h))
-
-
-/**
- * @brief Divide out homogeneous parameter from hvect_t, creating
- * vect_t.
- */
-#define HDIVIDE(o, v) do { \
-	(o)[X] = (v)[X] / (v)[W]; \
-	(o)[Y] = (v)[Y] / (v)[W]; \
-	(o)[Z] = (v)[Z] / (v)[W]; \
-    } while (0)
-
-/**
- * @brief Quaternion math definitions.
- *
- * Note that the [W] component will be put in the last (i.e., third)
- * place rather than the first [X] (i.e., [0]) place, so that the X,
- * Y, and Z elements will be compatible with vectors.  Only
- * QUAT_FROM_VROT macros depend on component locations, however.
- */
-
-/**
- * @brief Create Quaternion from Vector and Rotation about vector.
- *
- * To produce a quaternion representing a rotation by PI radians about
- * X-axis:
- *
- * VSET(axis, 1, 0, 0);
- * QUAT_FROM_VROT(quat, M_PI, axis);
- * or
- * QUAT_FROM_ROT(quat, M_PI, 1.0, 0.0, 0.0, 0.0);
- *
- * Alternatively, in degrees:
- * QUAT_FROM_ROT_DEG(quat, 180.0, 1.0, 0.0, 0.0, 0.0);
- */
-#define QUAT_FROM_ROT(q, r, x, y, z) do { \
-	fastf_t _rot = (r) * 0.5; \
-	QSET(q, x, y, z, cos(_rot)); \
-	VUNITIZE(q); \
-	_rot = sin(_rot); /* _rot is really just a temp variable now */ \
-	VSCALE(q, q, _rot); \
-    } while (0)
-
-#define QUAT_FROM_VROT(q, r, v) do { \
-	fastf_t _rot = (r) * 0.5; \
-	VMOVE(q, v); \
-	VUNITIZE(q); \
-	(q)[W] = cos(_rot); \
-	_rot = sin(_rot); /* _rot is really just a temp variable now */ \
-	VSCALE(q, q, _rot); \
-    } while (0)
-
-#define QUAT_FROM_VROT_DEG(q, r, v) \
-	QUAT_FROM_VROT(q, ((r)*DEG2RAD), v)
-
-#define QUAT_FROM_ROT_DEG(q, r, x, y, z) \
-	QUAT_FROM_ROT(q, ((r)*DEG2RAD), x, y, z)
-
-
-/**
- * @brief Set quaternion at `a' to have coordinates `b', `c', `d', and
- * `e'.
- */
-#define QSET(a, b, c, d, e) do { \
-	(a)[X] = (b); \
-	(a)[Y] = (c); \
-	(a)[Z] = (d); \
-	(a)[W] = (e); \
-    } while (0)
-
-/** @brief Transfer quaternion at `b' to quaternion at `a'. */
-#define QMOVE(a, b) do { \
-	(a)[X] = (b)[X]; \
-	(a)[Y] = (b)[Y]; \
-	(a)[Z] = (b)[Z]; \
-	(a)[W] = (b)[W]; \
-    } while (0)
-
-/** @brief Add quaternions at `b' and `c', store result at `a'. */
-#define QADD2(a, b, c) do { \
-	(a)[X] = (b)[X] + (c)[X]; \
-	(a)[Y] = (b)[Y] + (c)[Y]; \
-	(a)[Z] = (b)[Z] + (c)[Z]; \
-	(a)[W] = (b)[W] + (c)[W]; \
-    } while (0)
-
-/**
- * @brief Subtract quaternion at `c' from quaternion at `b', store
- * result at `a'.
- */
-#define QSUB2(a, b, c) do { \
-	(a)[X] = (b)[X] - (c)[X]; \
-	(a)[Y] = (b)[Y] - (c)[Y]; \
-	(a)[Z] = (b)[Z] - (c)[Z]; \
-	(a)[W] = (b)[W] - (c)[W]; \
-    } while (0)
-
-/**
- * @brief Scale quaternion at `b' by scalar `c', store result at
- * `a'.
- */
-#define QSCALE(a, b, c) do { \
-	(a)[X] = (b)[X] * (c); \
-	(a)[Y] = (b)[Y] * (c); \
-	(a)[Z] = (b)[Z] * (c); \
-	(a)[W] = (b)[W] * (c); \
-    } while (0)
-
-/** @brief Normalize quaternion 'a' to be a unit quaternion. */
-#define QUNITIZE(a) do { \
-	double _f; \
-	_f = QMAGNITUDE(a); \
-	if (_f < VDIVIDE_TOL) _f = 0.0; else _f = 1.0/_f; \
-	(a)[X] *= _f; (a)[Y] *= _f; (a)[Z] *= _f; (a)[W] *= _f; \
-    } while (0)
-
-/** @brief Return scalar magnitude squared of quaternion at `a'. */
-#define QMAGSQ(a) \
-	((a)[X]*(a)[X] + (a)[Y]*(a)[Y] \
-	+ (a)[Z]*(a)[Z] + (a)[W]*(a)[W])
-
-/** @brief Return scalar magnitude of quaternion at `a'. */
-#define QMAGNITUDE(a) sqrt(QMAGSQ(a))
-
-/** @brief Compute dot product of quaternions at `a' and `b'. */
-#define QDOT(a, b) \
-	((a)[X]*(b)[X] + (a)[Y]*(b)[Y] \
-	+ (a)[Z]*(b)[Z] + (a)[W]*(b)[W])
-
-/**
- * @brief Compute quaternion product a = b * c
- *
- * a[W] = b[W]*c[W] - VDOT(b, c);
- * VCROSS(temp, b, c);
- * VJOIN2(a, temp, b[W], c, c[W], b);
- */
-#define QMUL(a, b, c) do { \
-	(a)[W] = (b)[W]*(c)[W] - (b)[X]*(c)[X] - (b)[Y]*(c)[Y] - (b)[Z]*(c)[Z]; \
-	(a)[X] = (b)[W]*(c)[X] + (b)[X]*(c)[W] + (b)[Y]*(c)[Z] - (b)[Z]*(c)[Y]; \
-	(a)[Y] = (b)[W]*(c)[Y] + (b)[Y]*(c)[W] + (b)[Z]*(c)[X] - (b)[X]*(c)[Z]; \
-	(a)[Z] = (b)[W]*(c)[Z] + (b)[Z]*(c)[W] + (b)[X]*(c)[Y] - (b)[Y]*(c)[X]; \
-    } while (0)
-
-/** @brief Conjugate quaternion */
-#define QCONJUGATE(a, b) do { \
-	(a)[X] = -(b)[X]; \
-	(a)[Y] = -(b)[Y]; \
-	(a)[Z] = -(b)[Z]; \
-	(a)[W] =  (b)[W]; \
-    } while (0)
-
-/** @brief Multiplicative inverse quaternion */
-#define QINVERSE(a, b) do { \
-	double _f = QMAGSQ(b); \
-	if (_f < VDIVIDE_TOL) _f = 0.0; else _f = 1.0/_f; \
-	(a)[X] = -(b)[X] * _f; \
-	(a)[Y] = -(b)[Y] * _f; \
-	(a)[Z] = -(b)[Z] * _f; \
-	(a)[W] =  (b)[W] * _f; \
-    } while (0)
-
-/**
- * @brief Blend into quaternion `a'
- *
- * scalar `b' times quaternion at `c' plus
- * scalar `d' times quaternion at `e'
- */
-#define QBLEND2(a, b, c, d, e) do { \
-	(a)[X] = (b) * (c)[X] + (d) * (e)[X]; \
-	(a)[Y] = (b) * (c)[Y] + (d) * (e)[Y]; \
-	(a)[Z] = (b) * (c)[Z] + (d) * (e)[Z]; \
-	(a)[W] = (b) * (c)[W] + (d) * (e)[W]; \
-    } while (0)
-
-/**
- * Macros for dealing with 3-D "extents", aka bounding boxes, that are
- * represented as axis-aligned right parallelepipeds (RPPs).  This is
- * stored as two points: a min point, and a max point.  RPP 1 is
- * defined by lo1, hi1, RPP 2 by lo2, hi2.
- */
-
-/**
- * Compare two bounding boxes and return true if they are disjoint.
- */
-#define V3RPP_DISJOINT(_l1, _h1, _l2, _h2) \
-      ((_l1)[X] > (_h2)[X] || (_l1)[Y] > (_h2)[Y] || (_l1)[Z] > (_h2)[Z] || \
-	(_l2)[X] > (_h1)[X] || (_l2)[Y] > (_h1)[Y] || (_l2)[Z] > (_h1)[Z])
-
-/**
- * Compare two bounding boxes and return true if they are disjoint
- * by at least distance tolerance.
- */
-#define V3RPP_DISJOINT_TOL(_l1, _h1, _l2, _h2, _t) \
-       ((_l1)[X] > (_h2)[X] + (_t) || \
-	(_l1)[Y] > (_h2)[Y] + (_t) || \
-	(_l1)[Z] > (_h2)[Z] + (_t) || \
-	(_l2)[X] > (_h1)[X] + (_t) || \
-	(_l2)[Y] > (_h1)[Y] + (_t) || \
-	(_l2)[Z] > (_h1)[Z] + (_t))
-
-/** Compare two bounding boxes and return true If they overlap. */
-#define V3RPP_OVERLAP(_l1, _h1, _l2, _h2) \
-    (! ((_l1)[X] > (_h2)[X] || (_l1)[Y] > (_h2)[Y] || (_l1)[Z] > (_h2)[Z] || \
-	(_l2)[X] > (_h1)[X] || (_l2)[Y] > (_h1)[Y] || (_l2)[Z] > (_h1)[Z]))
-
-/**
- * @brief If two extents overlap within distance tolerance, return
- * true.
- */
-#define V3RPP_OVERLAP_TOL(_l1, _h1, _l2, _h2, _t) \
-    (! ((_l1)[X] > (_h2)[X] + (_t) || \
-	(_l1)[Y] > (_h2)[Y] + (_t) || \
-	(_l1)[Z] > (_h2)[Z] + (_t) || \
-	(_l2)[X] > (_h1)[X] + (_t) || \
-	(_l2)[Y] > (_h1)[Y] + (_t) || \
-	(_l2)[Z] > (_h1)[Z] + (_t)))
-
-/**
- * @brief Is the point within or on the boundary of the RPP?
- *
- * FIXME: should not be using >= <=, '=' case is unreliable
- */
-#define V3PNT_IN_RPP(_pt, _lo, _hi)	(\
-	(_pt)[X] >= (_lo)[X] && (_pt)[X] <= (_hi)[X] && \
-	(_pt)[Y] >= (_lo)[Y] && (_pt)[Y] <= (_hi)[Y] && \
-	(_pt)[Z] >= (_lo)[Z] && (_pt)[Z] <= (_hi)[Z])
-
-/**
- * @brief Within the distance tolerance, is the point within the RPP?
- *
- * FIXME: should not be using >= <=, '=' case is unreliable
- */
-#define V3PNT_IN_RPP_TOL(_pt, _lo, _hi, _t)	(\
-	(_pt)[X] >= (_lo)[X]-(_t) && (_pt)[X] <= (_hi)[X]+(_t) && \
-	(_pt)[Y] >= (_lo)[Y]-(_t) && (_pt)[Y] <= (_hi)[Y]+(_t) && \
-	(_pt)[Z] >= (_lo)[Z]-(_t) && (_pt)[Z] <= (_hi)[Z]+(_t))
-
-/**
- * @brief Is the point outside the RPP by at least the distance tolerance?
- * This will not return true if the point is on the RPP.
- */
-#define V3PNT_OUT_RPP_TOL(_pt, _lo, _hi, _t)      (\
-	(_pt)[X] < (_lo)[X]-(_t) || (_pt)[X] > (_hi)[X]+(_t) || \
-	(_pt)[Y] < (_lo)[Y]-(_t) || (_pt)[Y] > (_hi)[Y]+(_t) || \
-	(_pt)[Z] < (_lo)[Z]-(_t) || (_pt)[Z] > (_hi)[Z]+(_t))
-
-/**
- * @brief Determine if one bounding box is within another.  Also
- * returns true if the boxes are the same.
- *
- * FIXME: should not be using >= <=, '=' case is unreliable
- */
-#define V3RPP1_IN_RPP2(_lo1, _hi1, _lo2, _hi2)	(\
-	(_lo1)[X] >= (_lo2)[X] && (_hi1)[X] <= (_hi2)[X] && \
-	(_lo1)[Y] >= (_lo2)[Y] && (_hi1)[Y] <= (_hi2)[Y] && \
-	(_lo1)[Z] >= (_lo2)[Z] && (_hi1)[Z] <= (_hi2)[Z])
-
-/** Convert an azimuth/elevation to a direction vector. */
-#define V3DIR_FROM_AZEL(_d, _a, _e) do { \
-	fastf_t _c_e = cos(_e); \
-	(_d)[X] = cos(_a) * _c_e; \
-	(_d)[Y] = sin(_a) * _c_e; \
-	(_d)[Z] = sin(_e); \
-    } while (0)
-
-/** Convert a direction vector to azimuth/elevation (in radians). */
-#define AZEL_FROM_V3DIR(_a, _e, _d) do { \
-	(_a) = ((NEAR_ZERO((_d)[X], SMALL_FASTF)) && (NEAR_ZERO((_d)[Y], SMALL_FASTF))) ? 0.0 : atan2(-((_d)[Y]), -((_d)[X])) * -RAD2DEG; \
-	(_e) = atan2(-((_d)[Z]), sqrt((_d)[X]*(_d)[X] + (_d)[Y]*(_d)[Y])) * -RAD2DEG; \
-    } while (0)
-
-
-/** Swap two 3D vectors */
-#define VSWAP(_a, _b) do { \
-	fastf_t _t; \
-	_t = (_a)[X]; \
-	(_a)[X] = (_b)[X]; \
-	(_b)[X] = _t; \
-	_t = (_a)[Y]; \
-	(_a)[Y] = (_b)[Y]; \
-	(_b)[Y] = _t; \
-	_t = (_a)[Z]; \
-	(_a)[Z] = (_b)[Z]; \
-	(_b)[Z] = _t; \
-    } while (0)
-
-/** Swap two 2D vectors */
-#define V2SWAP(_a, _b) do { \
-	fastf_t _t; \
-	_t = (_a)[X]; \
-	(_a)[X] = (_b)[X]; \
-	(_b)[X] = _t; \
-	_t = (_a)[Y]; \
-	(_a)[Y] = (_b)[Y]; \
-	(_b)[Y] = _t; \
-    } while (0)
-
-/** Swap two 4D vectors */
-#define HSWAP(_a, _b) do { \
-	fastf_t _t; \
-	_t = (_a)[X]; \
-	(_a)[X] = (_b)[X]; \
-	(_b)[X] = _t; \
-	_t = (_a)[Y]; \
-	(_a)[Y] = (_b)[Y]; \
-	(_b)[Y] = _t; \
-	_t = (_a)[Z]; \
-	(_a)[Z] = (_b)[Z]; \
-	(_b)[Z] = _t; \
-	_t = (_a)[W]; \
-	(_a)[W] = (_b)[W]; \
-	(_b)[W] = _t; \
-    } while (0)
-
-/** Swap two 4x4 matrices */
-#define MAT_SWAP(_a, _b) do { \
-	mat_t _t; \
-	MAT_COPY(_t, (_a)); \
-	MAT_COPY((_a), (_b)); \
-	MAT_COPY((_b), _t); \
-    } while (0)
-
-/*** Macros suitable for declaration statement initialization. ***/
-
-/**
- * 3D vector macro suitable for declaration statement initialization.
- * this sets all vector elements to the specified value similar to
- * VSETALL() but as an initializer array declaration instead of as a
- * statement.
- */
-#define VINITALL(_v) {(_v), (_v), (_v)}
-
-/**
- * 2D vector macro suitable for declaration statement initialization.
- * this sets all vector elements to the specified value similar to
- * VSETALLN(hvect_t,val,2) but as an initializer array declaration
- * instead of as a statement.
- */
-#define V2INITALL(_v) {(_v), (_v), (_v)}
-
-/**
- * 4D homogeneous vector macro suitable for declaration statement
- * initialization.  this sets all vector elements to the specified
- * value similar to VSETALLN(hvect_t,val,4) but as an initializer
- * array declaration instead of as a statement.
- */
-#define HINITALL(_v) {(_v), (_v), (_v), (_v)}
-
-/**
- * 3D vector macro suitable for declaration statement initialization.
- * this sets all vector elements to zero similar to calling
- * VSETALL(0.0) but as an initializer array declaration instead of as
- * a statement.
- */
-#define VINIT_ZERO {0.0, 0.0, 0.0}
-
-/**
- * 2D vector macro suitable for declaration statement initialization.
- * this sets all vector elements to zero similar to calling
- * V2SETALL(0.0) but as an initializer array declaration instead of as
- * a statement.
- */
-#define V2INIT_ZERO {0.0, 0.0}
-
-/**
- * 4D homogeneous vector macro suitable for declaration statement
- * initialization.  this sets all vector elements to zero similar to
- * calling VSETALLN(hvect_t,0.0,4) but as an initializer array
- * declaration instead of as a statement.
- */
-#define HINIT_ZERO {0.0, 0.0, 0.0, 0.0}
-
-/**
- * matrix macro suitable for declaration statement initialization.
- * this sets up an identity matrix similar to calling MAT_IDN but as
- * an initializer array declaration instead of as a statement.
- */
-#define MAT_INIT_IDN {1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0}
-
-/**
- * matrix macro suitable for declaration statement initialization.
- * this sets up a zero matrix similar to calling MAT_ZERO but as an
- * initializer array declaration instead of as a statement.
- */
-#define MAT_INIT_ZERO {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}
-
-
-#ifdef __cplusplus
-} /* end extern "C" */
-#endif
-
-#endif /* VMATH_H */
-
-/** @} */
-/*
- * Local Variables:
- * mode: C
- * tab-width: 8
- * indent-tabs-mode: t
- * c-file-style: "stroustrup"
- * End:
- * ex: shiftwidth=4 tabstop=8
- */
diff --git a/bindings/java/src/main/java/manifold3d/ConvexHull.java b/bindings/java/src/main/java/manifold3d/ConvexHull.java
deleted file mode 100644
index b9a557ce7..000000000
--- a/bindings/java/src/main/java/manifold3d/ConvexHull.java
+++ /dev/null
@@ -1,33 +0,0 @@
-package manifold3d;
-
-import org.bytedeco.javacpp.*;
-import org.bytedeco.javacpp.annotation.*;
-
-import manifold3d.LibraryPaths;
-import manifold3d.Manifold;
-import manifold3d.manifold.CrossSection;
-import manifold3d.glm.DoubleVec3Vector;
-import manifold3d.glm.DoubleMat4x3Vector;
-import manifold3d.glm.DoubleMat4x3Vector;
-import manifold3d.UIntVecVector;
-
-import manifold3d.Manifold;
-import java.nio.DoubleBuffer;
-import java.nio.IntBuffer;
-
-@Platform(compiler = "cpp17", include = {"convex_hull.hpp"}, linkpath = { LibraryPaths.MANIFOLD_LIB_DIR, LibraryPaths.MANIFOLD_LIB_DIR_WINDOWS, LibraryPaths.QHULL}, link = {"manifold"})
-public class ConvexHull extends Pointer {
-    static { Loader.load(); }
-
-    public ConvexHull() { }
-
-    public static native @ByVal Manifold ConvexHull(@Const @ByRef Manifold manifold);
-    public static native @ByVal Manifold ConvexHull(@Const @ByRef Manifold manifold, @Const float precision);
-    public static native @ByVal Manifold ConvexHull(@Const @ByRef Manifold manifold, @Const @ByRef Manifold other);
-    public static native @ByVal Manifold ConvexHull(@Const @ByRef Manifold manifold, @Const @ByRef Manifold other, float precision);
-
-    public static native @ByVal CrossSection ConvexHull(@Const @ByRef CrossSection crossSection);
-    public static native @ByVal CrossSection ConvexHull(@Const @ByRef CrossSection crossSection, @Const float precision);
-    public static native @ByVal CrossSection ConvexHull(@Const @ByRef CrossSection crossSection, @Const @ByRef CrossSection other);
-    public static native @ByVal CrossSection ConvexHull(@Const @ByRef CrossSection crossSection, @Const @ByRef CrossSection other, float precision);
-}
diff --git a/bindings/java/src/main/java/manifold3d/LibraryPaths.java b/bindings/java/src/main/java/manifold3d/LibraryPaths.java
index 6c9374047..66c71ee99 100644
--- a/bindings/java/src/main/java/manifold3d/LibraryPaths.java
+++ b/bindings/java/src/main/java/manifold3d/LibraryPaths.java
@@ -1,11 +1,11 @@
 package manifold3d;
 
 public class LibraryPaths {
-    public static final String QHULL = "src/third_party/qhull/lib";
     public static final String MANIFOLD_BUILD_DIR = "../../build";
 
     public static final String MANIFOLD_LIB_DIR = MANIFOLD_BUILD_DIR + "/src/manifold";
     public static final String CROSS_SECTION_LIB_DIR = MANIFOLD_BUILD_DIR + "/src/cross_section";
+    public static final String QUICKHULL_LIB_DIR = MANIFOLD_BUILD_DIR + "/src/third_party";
     public static final String MANIFOLD_LIB_DIR_WINDOWS = MANIFOLD_LIB_DIR + "/Release";
     public static final String MESHIO_LIB_DIR = MANIFOLD_BUILD_DIR + "/meshIO";
     public static final String MESHIO_LIB_DIR_WINDOWS = MESHIO_LIB_DIR + "/Release";
diff --git a/bindings/java/src/main/java/manifold3d/Manifold.java b/bindings/java/src/main/java/manifold3d/Manifold.java
index 929ff68ae..d742cd895 100644
--- a/bindings/java/src/main/java/manifold3d/Manifold.java
+++ b/bindings/java/src/main/java/manifold3d/Manifold.java
@@ -19,7 +19,6 @@
 import manifold3d.manifold.CrossSection;
 
 import manifold3d.pub.DoubleMesh;
-import manifold3d.ConvexHull;
 import manifold3d.pub.Box;
 import manifold3d.pub.Properties;
 import manifold3d.pub.SmoothnessVector;
@@ -36,8 +35,7 @@
 @Platform(compiler = "cpp17", include = {"manifold.h", "meshIO.h"},
           linkpath = { LibraryPaths.MANIFOLD_LIB_DIR,
                        LibraryPaths.CROSS_SECTION_LIB_DIR,
-                       LibraryPaths.MANIFOLD_LIB_DIR_WINDOWS,
-                       LibraryPaths.QHULL },
+                       LibraryPaths.MANIFOLD_LIB_DIR_WINDOWS },
           link = { "manifold", "cross_section" })
 @Namespace("manifold")
 public class Manifold extends Pointer {
@@ -50,6 +48,7 @@ public class Manifold extends Pointer {
                 System.load(Loader.extractResource("/libClipper2.so.1.2.1", null, "libClipper2", ".so").getAbsolutePath());
                 System.load(Loader.extractResource("/libcross_section.so", null, "libcross_section", ".so").getAbsolutePath());
                 System.load(Loader.extractResource("/libmanifold.so", null, "libmanifold", ".so").getAbsolutePath());
+                System.load(Loader.extractResource("/libquickhull.so", null, "libquickhull", ".so").getAbsolutePath());
             } catch (IOException e) {
                 throw new RuntimeException(e);
             }
@@ -63,6 +62,8 @@ public class Manifold extends Pointer {
                 System.load(Loader.extractResource("/cross_section.dll", null, "cross_section", ".dll").getAbsolutePath());
                 System.out.println("Loading manifold");
                 System.load(Loader.extractResource("/manifold.dll", null, "manifold", ".dll").getAbsolutePath());
+                System.out.println("Loading QuickHull");
+                System.load(Loader.extractResource("/quickhull.dll", null, "quickhull", ".dll").getAbsolutePath());
                 System.out.println("Finished Loading.");
             } catch (IOException e) {
                 throw new RuntimeException(e);
@@ -77,6 +78,8 @@ public class Manifold extends Pointer {
                 System.load(Loader.extractResource("/libcross_section.dylib", null, "libcross_section", ".dylib").getAbsolutePath());
                 System.out.println("Loading Manifold");
                 System.load(Loader.extractResource("/libmanifold.dylib", null, "libmanifold", ".dylib").getAbsolutePath());
+                System.out.println("Loading QuickHull");
+                System.load(Loader.extractResource("/libquickhull.dylib", null, "libquickhull", ".dylib").getAbsolutePath());
                 System.out.println("Finished Loading.");
 
             } catch (IOException e) {
@@ -122,13 +125,12 @@ public class Manifold extends Pointer {
     @Name("OriginalID") public native int originalID();
     @Name("AsOriginal") public native @ByVal Manifold asOriginal();
 
-    //@Name("ConvexHull") public native @ByVal Manifold convexHull(@ByRef Manifold other);
+    @Name("Hull") public native @ByVal Manifold convexHull();
+    @Name("Hull") public native @ByVal Manifold convexHull(@ByRef ManifoldVector others);
 
-    public @ByVal Manifold convexHull() {
-        return ConvexHull.ConvexHull(this);
-    }
     public @ByVal Manifold convexHull(@ByRef Manifold other) {
-        return ConvexHull.ConvexHull(this, other);
+        Manifold[] vec  = new Manifold[]{other};
+        return this.convexHull(new ManifoldVector(vec));
     }
 
     //// Modifiers
diff --git a/bindings/java/src/main/java/manifold3d/manifold/CrossSection.java b/bindings/java/src/main/java/manifold3d/manifold/CrossSection.java
index ea2b0eb85..0254ecade 100644
--- a/bindings/java/src/main/java/manifold3d/manifold/CrossSection.java
+++ b/bindings/java/src/main/java/manifold3d/manifold/CrossSection.java
@@ -4,7 +4,6 @@
 import org.bytedeco.javacpp.annotation.*;
 
 import manifold3d.LibraryPaths;
-import manifold3d.ConvexHull;
 import manifold3d.glm.DoubleVec2;
 import manifold3d.glm.DoubleMat3x2;
 import manifold3d.manifold.Rect;
@@ -13,7 +12,12 @@
 import manifold3d.pub.SimplePolygon;
 import manifold3d.pub.Polygons;
 
-@Platform(compiler = "cpp17", include = { "cross_section.h" }, linkpath = { LibraryPaths.CROSS_SECTION_LIB_DIR, LibraryPaths.MANIFOLD_LIB_DIR_WINDOWS }, link = { "cross_section" })
+@Platform(compiler = "cpp17",
+          include = { "cross_section.h" },
+          linkpath = { LibraryPaths.CROSS_SECTION_LIB_DIR,
+                       LibraryPaths.MANIFOLD_LIB_DIR,
+                       LibraryPaths.QUICKHULL_LIB_DIR },
+          link = { "cross_section", "quickhull", "manifold" })
 @Namespace("manifold")
 public class CrossSection extends Pointer {
     static { Loader.load(); }
@@ -21,8 +25,11 @@ public class CrossSection extends Pointer {
     public CrossSection() { allocate(); }
     private native void allocate();
 
+    public CrossSection(@Const @ByRef CrossSection other) { allocate(other); }
+    private native void allocate(@Const @ByRef CrossSection other);
+
     @Name("operator=")
-    public native @ByRef CrossSection put(@ByRef CrossSection other);
+    public native @ByRef CrossSection put(@Const @ByRef CrossSection other);
 
     public enum FillRule {
         EvenOdd,   ///< Only odd numbered sub-regions are filled.
@@ -74,20 +81,13 @@ public CrossSection translateY(double y) {
     @Name("Boolean") public native @ByVal CrossSection booleanOp(@ByRef CrossSection second, @Cast("manifold::OpType") int op);
     public static native @ByVal CrossSection BatchBoolean(@ByRef CrossSectionVector sections, @Cast("manifold::OpType") int op);
 
-    public @ByVal CrossSection convexHull() {
-        return ConvexHull.ConvexHull(this);
-    }
-
-    public @ByVal CrossSection convexHull(float precision) {
-        return ConvexHull.ConvexHull(this, precision);
-    }
-
-    public @ByVal CrossSection convexHull(@Const @ByRef CrossSection other) {
-       return ConvexHull.ConvexHull(this, other);
-    }
+    @Name("Hull") public native @ByVal CrossSection convexHull();
+    @Name("Hull") public native @ByVal CrossSection convexHull(@ByRef SimplePolygon pts);
+    @Name("Hull") public native @ByVal CrossSection convexHull(@ByRef Polygons pts);
+    @Name("Hull") public native @ByVal CrossSection convexHull(@ByRef CrossSectionVector sections);
 
-    public @ByVal CrossSection convexHull(@Const @ByRef CrossSection other, float precision) {
-       return ConvexHull.ConvexHull(this, other, precision);
+    public @ByVal CrossSection convexHull(@ByRef CrossSection other) {
+        return this.convexHull(new CrossSectionVector(new CrossSection[]{other}));
     }
 
     @Name("operator+") public native @ByVal CrossSection add(@ByRef CrossSection rhs);
diff --git a/bindings/java/src/main/java/manifold3d/manifold/CrossSectionVector.java b/bindings/java/src/main/java/manifold3d/manifold/CrossSectionVector.java
index b73ee733c..e5d073422 100644
--- a/bindings/java/src/main/java/manifold3d/manifold/CrossSectionVector.java
+++ b/bindings/java/src/main/java/manifold3d/manifold/CrossSectionVector.java
@@ -11,7 +11,11 @@
 import java.lang.Iterable;
 import java.util.NoSuchElementException;
 
-@Platform(compiler = "cpp17", include = {"manifold.h", "<vector>"}, linkpath = { LibraryPaths.MANIFOLD_LIB_DIR, LibraryPaths.MANIFOLD_LIB_DIR_WINDOWS }, link = { "manifold" })
+@Platform(compiler = "cpp17",
+          include = {"manifold.h", "<vector>"},
+          linkpath = { LibraryPaths.CROSS_SECTION_LIB_DIR,
+                       LibraryPaths.MANIFOLD_LIB_DIR_WINDOWS },
+          link = { "cross_section" })
 @Name("std::vector<manifold::CrossSection>")
 public class CrossSectionVector extends Pointer implements Iterable<CrossSection> {
     static { Loader.load(); }