GteConvexHull3.h

Go to the documentation of this file.

// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2017
// Distributed under the Boost Software License, Version 1.0.
// http://www.boost.org/LICENSE_1_0.txt
// http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// File Version: 3.0.0 (2016/06/19)

#pragma once

// Compute the convex hull of 3D points using incremental insertion.  The only
// way to ensure a correct result for the input vertices (assumed to be exact)
// is to choose ComputeType for exact rational arithmetic.  You may use
// BSNumber.  No divisions are performed in this computation, so you do not
// have to use BSRational.
//
// The worst-case choices of N for Real of type BSNumber or BSRational with
// integer storage UIntegerFP32<N> are listed in the next table.  The numerical
// computations are encapsulated in PrimalQuery3<Real>::ToPlane.  We recommend
// using only BSNumber, because no divisions are performed in the convex-hull
// computations.
//
//    input type | compute type | N
//    -----------+--------------+------
//    float      | BSNumber     |    27
//    double     | BSNumber     |   197
//    float      | BSRational   |  2882
//    double     | BSRational   | 21688

#include <Mathematics/GteETManifoldMesh.h>
#include <Mathematics/GtePrimalQuery3.h>
#include <Mathematics/GteLine.h>
#include <Mathematics/GteHyperplane.h>
#include <functional>
#include <set>
#include <thread>
#include <vector>

namespace gte
{

template <typename InputType, typename ComputeType>
class ConvexHull3
{
public:
    // The class is a functor to support computing the convex hull of multiple
    // data sets using the same class object.  For multithreading in Update,
    // choose 'numThreads' subject to the constraints
    //     1 <= numThreads <= std::thread::hardware_concurrency().
    ConvexHull3(unsigned int numThreads = 1);

    // The input is the array of points whose convex hull is required.  The
    // epsilon value is used to determine the intrinsic dimensionality of the
    // vertices (d = 0, 1, 2, or 3).  When epsilon is positive, the
    // determination is fuzzy--points approximately the same point,
    // approximately on a line, approximately planar, or volumetric.
    bool operator()(int numPoints, Vector3<InputType> const* points,
        InputType epsilon);

    // Dimensional information.  If GetDimension() returns 1, the points lie
    // on a line P+t*D (fuzzy comparison when epsilon > 0).  You can sort
    // these if you need a polyline output by projecting onto the line each
    // vertex X = P+t*D, where t = Dot(D,X-P).  If GetDimension() returns 2,
    // the points line on a plane P+s*U+t*V (fuzzy comparison when
    // epsilon > 0).  You can project each point X = P+s*U+t*V, where
    // s = Dot(U,X-P) and t = Dot(V,X-P), then apply ConvexHull2 to the (s,t)
    // tuples.
    inline InputType GetEpsilon() const;
    inline int GetDimension() const;
    inline Line3<InputType> const& GetLine() const;
    inline Plane3<InputType> const& GetPlane() const;

    // Member access.
    inline int GetNumPoints() const;
    inline int GetNumUniquePoints() const;
    inline Vector3<InputType> const* GetPoints() const;
    inline PrimalQuery3<ComputeType> const& GetQuery() const;

    // The convex hull is a convex polyhedron with triangular faces.
    inline std::vector<TriangleKey<true>> const& GetHullUnordered() const;
    ETManifoldMesh const& GetHullMesh() const;

private:
    // Support for incremental insertion.
    void Update(int i);

    // The epsilon value is used for fuzzy determination of intrinsic
    // dimensionality.  If the dimension is 0, 1, or 2, the constructor
    // returns early.  The caller is responsible for retrieving the dimension
    // and taking an alternate path should the dimension be smaller than 3.
    // If the dimension is 0, the caller may as well treat all points[] as a
    // single point, say, points[0].  If the dimension is 1, the caller can
    // query for the approximating line and project points[] onto it for
    // further processing.  If the dimension is 2, the caller can query for
    // the approximating plane and project points[] onto it for further
    // processing.
    InputType mEpsilon;
    int mDimension;
    Line3<InputType> mLine;
    Plane3<InputType> mPlane;

    // The array of points used for geometric queries.  If you want to be
    // certain of a correct result, choose ComputeType to be BSNumber.
    std::vector<Vector3<ComputeType>> mComputePoints;
    PrimalQuery3<ComputeType> mQuery;

    int mNumPoints;
    int mNumUniquePoints;
    Vector3<InputType> const* mPoints;
    std::vector<TriangleKey<true>> mHullUnordered;
    mutable ETManifoldMesh mHullMesh;
    unsigned int mNumThreads;
};


template <typename InputType, typename ComputeType>
ConvexHull3<InputType, ComputeType>::ConvexHull3(unsigned int numThreads)
    :
    mEpsilon((InputType)0),
    mDimension(0),
    mLine(Vector3<InputType>::Zero(), Vector3<InputType>::Zero()),
    mPlane(Vector3<InputType>::Zero(), (InputType)0),
    mNumPoints(0),
    mNumUniquePoints(0),
    mPoints(nullptr),
    mNumThreads(numThreads)
{
}

template <typename InputType, typename ComputeType>
bool ConvexHull3<InputType, ComputeType>::operator()(int numPoints,
    Vector3<InputType> const* points, InputType epsilon)
{
    mEpsilon = std::max(epsilon, (InputType)0);
    mDimension = 0;
    mLine.origin = Vector3<InputType>::Zero();
    mLine.direction = Vector3<InputType>::Zero();
    mPlane.normal = Vector3<InputType>::Zero();
    mPlane.constant = (InputType)0;
    mNumPoints = numPoints;
    mNumUniquePoints = 0;
    mPoints = points;
    mHullUnordered.clear();
    mHullMesh.Clear();

    int i, j;
    if (mNumPoints < 4)
    {
        // ConvexHull3 should be called with at least four points.
        return false;
    }

    IntrinsicsVector3<InputType> info(mNumPoints, mPoints, mEpsilon);
    if (info.dimension == 0)
    {
        // The set is (nearly) a point.
        return false;
    }

    if (info.dimension == 1)
    {
        // The set is (nearly) collinear.
        mDimension = 1;
        mLine = Line3<InputType>(info.origin, info.direction[0]);
        return false;
    }

    if (info.dimension == 2)
    {
        // The set is (nearly) coplanar.
        mDimension = 2;
        mPlane = Plane3<InputType>(UnitCross(info.direction[0],
            info.direction[1]), info.origin);
        return false;
    }

    mDimension = 3;

    // Compute the vertices for the queries.
    mComputePoints.resize(mNumPoints);
    mQuery.Set(mNumPoints, &mComputePoints[0]);
    for (i = 0; i < mNumPoints; ++i)
    {
        for (j = 0; j < 3; ++j)
        {
            mComputePoints[i][j] = points[i][j];
        }
    }

    // Insert the faces of the (nondegenerate) tetrahedron constructed by the
    // call to GetInformation.
    if (!info.extremeCCW)
    {
        std::swap(info.extreme[2], info.extreme[3]);
    }

    mHullUnordered.push_back(TriangleKey<true>(info.extreme[1],
        info.extreme[2], info.extreme[3]));
    mHullUnordered.push_back(TriangleKey<true>(info.extreme[0],
        info.extreme[3], info.extreme[2]));
    mHullUnordered.push_back(TriangleKey<true>(info.extreme[0],
        info.extreme[1], info.extreme[3]));
    mHullUnordered.push_back(TriangleKey<true>(info.extreme[0],
        info.extreme[2], info.extreme[1]));

    // Incrementally update the hull.  The set of processed points is
    // maintained to eliminate duplicates, either in the original input
    // points or in the points obtained by snap rounding.
    std::set<Vector3<InputType>> processed;
    for (i = 0; i < 4; ++i)
    {
        processed.insert(points[info.extreme[i]]);
    }
    for (i = 0; i < mNumPoints; ++i)
    {
        if (processed.find(points[i]) == processed.end())
        {
            Update(i);
            processed.insert(points[i]);
        }
    }
    mNumUniquePoints = static_cast<int>(processed.size());
    return true;
}

template <typename InputType, typename ComputeType> inline
InputType ConvexHull3<InputType, ComputeType>::GetEpsilon() const
{
    return mEpsilon;
}

template <typename InputType, typename ComputeType> inline
int ConvexHull3<InputType, ComputeType>::GetDimension() const
{
    return mDimension;
}

template <typename InputType, typename ComputeType> inline
Line3<InputType> const& ConvexHull3<InputType, ComputeType>::GetLine() const
{
    return mLine;
}

template <typename InputType, typename ComputeType> inline
Plane3<InputType> const& ConvexHull3<InputType, ComputeType>::GetPlane() const
{
    return mPlane;
}

template <typename InputType, typename ComputeType> inline
int ConvexHull3<InputType, ComputeType>::GetNumPoints() const
{
    return mNumPoints;
}

template <typename InputType, typename ComputeType> inline
int ConvexHull3<InputType, ComputeType>::GetNumUniquePoints() const
{
    return mNumUniquePoints;
}

template <typename InputType, typename ComputeType> inline
Vector3<InputType> const* ConvexHull3<InputType, ComputeType>::GetPoints() const
{
    return mPoints;
}

template <typename InputType, typename ComputeType> inline
PrimalQuery3<ComputeType> const& ConvexHull3<InputType, ComputeType>::GetQuery() const
{
    return mQuery;
}

template <typename InputType, typename ComputeType> inline
std::vector<TriangleKey<true>> const& ConvexHull3<InputType, ComputeType>::GetHullUnordered() const
{
    return mHullUnordered;
}

template <typename InputType, typename ComputeType>
ETManifoldMesh const& ConvexHull3<InputType, ComputeType>::GetHullMesh() const
{
    // Create the mesh only on demand.
    if (mHullMesh.GetTriangles().size() == 0)
    {
        for (auto const& tri : mHullUnordered)
        {
            mHullMesh.Insert(tri.V[0], tri.V[1], tri.V[2]);
        }
    }

    return mHullMesh;
}

template <typename InputType, typename ComputeType>
void ConvexHull3<InputType, ComputeType>::Update(int i)
{
    // The terminator that separates visible faces from nonvisible faces is
    // constructed by this code.  Visible faces for the incoming hull are
    // removed, and the boundary of that set of triangles is the terminator.
    // New visible faces are added using the incoming point and the edges of
    // the terminator.
    //
    // A simple algorithm for computing terminator edges is the following.
    // Back-facing triangles are located and the three edges are processed.
    // The first time an edge is visited, insert it into the terminator.  If
    // it is visited a second time, the edge is removed because it is shared
    // by another back-facing triangle and, therefore, cannot be a terminator
    // edge.  After visiting all back-facing triangles, the only remaining
    // edges in the map are the terminator edges.
    //
    // The order of vertices of an edge is important for adding new faces with
    // the correct vertex winding.  However, the edge "toggle" (insert edge,
    // remove edge) should use edges with unordered vertices, because the
    // edge shared by one triangle has opposite ordering relative to that of
    // the other triangle.  The map uses unordered edges as the keys but
    // stores the ordered edge as the value.  This avoids having to look up
    // an edge twice in a map with ordered edge keys.

    unsigned int numFaces = static_cast<unsigned int>(mHullUnordered.size());
    std::vector<int> queryResult(numFaces);
    if (mNumThreads > 1 && numFaces >= mNumThreads)
    {
        // Partition the data for multiple threads.
        unsigned int numFacesPerThread = numFaces / mNumThreads;
        std::vector<unsigned int> jmin(mNumThreads), jmax(mNumThreads);
        for (unsigned int t = 0; t < mNumThreads; ++t)
        {
            jmin[t] = t * numFacesPerThread;
            jmax[t] = jmin[t] + numFacesPerThread - 1;
        }
        jmax[mNumThreads - 1] = numFaces - 1;

        // Execute the point-plane queries in multiple threads.
        std::vector<std::thread> process(mNumThreads);
        for (unsigned int t = 0; t < mNumThreads; ++t)
        {
            process[t] = std::thread([this, i, t, &jmin, &jmax, 
                &queryResult]()
            {
                for (unsigned int j = jmin[t]; j <= jmax[t]; ++j)
                {
                    TriangleKey<true> const& tri = mHullUnordered[j];
                    queryResult[j] = mQuery.ToPlane(i, tri.V[0], tri.V[1], tri.V[2]);
                }
            });
        }

        // Wait for all threads to finish.
        for (unsigned int t = 0; t < mNumThreads; ++t)
        {
            process[t].join();
        }
    }
    else
    {
        unsigned int j = 0;
        for (auto const& tri : mHullUnordered)
        {
            queryResult[j++] = mQuery.ToPlane(i, tri.V[0], tri.V[1], tri.V[2]);
        }
    }

    std::map<EdgeKey<false>, std::pair<int, int>> terminator;
    std::vector<TriangleKey<true>> backFaces;
    bool existsFrontFacingTriangle = false;
    unsigned int j = 0;
    for (auto const& tri : mHullUnordered)
    {
        if (queryResult[j++] <= 0)
        {
            // The triangle is back facing.  These include triangles that
            // are coplanar with the incoming point.
            backFaces.push_back(tri);

            // The current hull is a 2-manifold watertight mesh.  The
            // terminator edges are those shared with a front-facing triangle.
            // The first time an edge of a back-facing triangle is visited,
            // insert it into the terminator.  If it is visited a second time,
            // the edge is removed because it is shared by another back-facing
            // triangle.  After all back-facing triangles are visited, the
            // only remaining edges are shared by a single back-facing
            // triangle, which makes them terminator edges.
            for (int j0 = 2, j1 = 0; j1 < 3; j0 = j1++)
            {
                int v0 = tri.V[j0], v1 = tri.V[j1];
                EdgeKey<false> edge(v0, v1);
                auto iter = terminator.find(edge);
                if (iter == terminator.end())
                {
                    // The edge is visited for the first time.
                    terminator.insert(std::make_pair(edge, std::make_pair(v0, v1)));
                }
                else
                {
                    // The edge is visited for the second time.
                    terminator.erase(edge);
                }
            }
        }
        else
        {
            // If there are no strictly front-facing triangles, then the
            // incoming point is inside or on the convex hull.  If we get
            // to this code, then the point is truly outside and we can
            // update the hull.
            existsFrontFacingTriangle = true;
        }
    }

    if (!existsFrontFacingTriangle)
    {
        // The incoming point is inside or on the current hull, so no
        // update of the hull is necessary.
        return;
    }

    // The updated hull contains the triangles not visible to the incoming
    // point.
    mHullUnordered = backFaces;

    // Insert the triangles formed by the incoming point and the terminator
    // edges.
    for (auto const& edge : terminator)
    {
        mHullUnordered.push_back(TriangleKey<true>(i, edge.second.second, edge.second.first));
    }
}


}