5
\$\begingroup\$

Last version of library (performance has been improved drastically since posting).

I tried to implement the Quick Hull Algorithm for computing the convex hull of a finite set of D-dimensional points. Algorithm itself:

#pragma once

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <deque>
#include <set>
#include <map>
#include <list>
#include <iterator>
#include <algorithm>
#include <utility>
#include <numeric>
#include <exception>

#include <cassert>

struct bad_geometry
        : std::exception
{

    virtual
    ~bad_geometry() noexcept = default;

    bad_geometry() = default;

    explicit
    bad_geometry(const char * const _what)
        : what_(_what)
    { ; }

    explicit
    bad_geometry(std::string const & _what)
        : what_(_what)
    { ; }

    virtual
    const char *
    what() const noexcept override
    {
        return what_.c_str();
    }

private :

    std::string const what_ = "bad_geometry";

};

using boolean_type = bool;
using size_type = std::size_t;

template< typename F >
F
determinant(boost::numeric::ublas::matrix< F > _m)
{
    size_type const size_ = _m.size1();
    assert(_m.size2() == size_);
    boost::numeric::ublas::permutation_matrix< size_type > pm_(size_);
    if (0 < boost::numeric::ublas::lu_factorize(_m, pm_)) {
        return F(0.0L); // singular matrix
    } else {
        F determinant_(1.0L);
        for (size_type i = 0; i < size_; ++i) {
            if (i == pm_(i)) {
                determinant_ *= +_m(i, i);
            } else {
                determinant_ *= -_m(i, i);
            }
        }
        return determinant_;
    }
}

template< typename point_type >
struct convex_hull
{

    using G = typename point_type::value_type;
    using point_refs_type = std::deque< std::reference_wrapper< point_type const > >;
    using point_list = std::list< size_type >;
    using point_set = std::set< size_type >;
    using points_type = std::deque< size_type >;

    template< typename ForwardIterator >
    convex_hull(ForwardIterator _first, ForwardIterator _last)
        : dimension_(_first->size())
        , points_(_first, _last)
    {
        assert(0 < dimension_);
        for (point_type const & point_ : points_) {
            if (point_.size() != dimension_) {
                throw bad_geometry("dimensionalities does not match");
            }
        }
    }

    size_type dimension_;
    point_refs_type points_;
    point_list internal_set_;

    struct facet;

    using facet_set = std::set< size_type >;

    struct facet // (d - 1)-dimensional faces
    {

        template< typename ForwardIterator >
        facet(ForwardIterator first, ForwardIterator mid, ForwardIterator last,
              boolean_type const _outward)
            : vertices_(first, std::prev(mid))
            , points_(vertices_.cbegin(), vertices_.cend())
            , outward_(_outward)
        {
            auto const rest = vertices_.insert(vertices_.cend(), mid, last);
            points_.insert(points_.cend(), rest, vertices_.end());
            std::sort(points_.begin(), points_.end());
        }

        facet(points_type && _vertices,
              boolean_type const _outward,
              size_type const _neighbour)
            : vertices_(std::move(_vertices))
            , points_(vertices_.cbegin(), vertices_.cend())
            , outward_(_outward)
            , neighbours_({_neighbour})
        {
            std::sort(points_.begin(), points_.end());
        }

        points_type vertices_; // oriented
        points_type points_;   // non-oriented
        boolean_type outward_;
        facet_set neighbours_;
        points_type outside_set_; // if not empty, then first point is furthest for this facet

        boolean_type
        further(G const & _nearer, G const & _further) const
        {
            if (outward_) {
                return (_nearer < _further);
            } else {
                return (_further < _nearer);
            }
        }

    };

    using facets_map = std::map< size_type, facet >;
    using facets_type = std::deque< size_type >;

    facets_map facets_;

    boolean_type
    below(facet const & _facet, G const & _orientation) const
    {
        if (_facet.outward_) {
            return (G(0.0L) < _orientation);
        } else {
            return (_orientation < -G(0.0L));
        }
    }

    // http://math.stackexchange.com/questions/822741/
    template< typename vertices >
    G
    orientation(vertices const & _vertices, point_type const & _apex) const
    {
        size_type const size_ = _vertices.size(); // dimensionality of the subspace of interest
        assert(!(_apex.size() < size_));
        boost::numeric::ublas::matrix< G > m_(size_ + 1, size_ + 1);
        auto v_ = _vertices.cbegin();
        for (size_type i = 0; i < size_; ++i) {
            assert(v_ != _vertices.cend());
            point_type const & vertex_ = points_.at(*v_);
            ++v_;
            assert(!(vertex_.size() < size_));
            for (size_type j = 0; j < size_; ++j) {
                m_(i, j) = vertex_[j];
            }
            m_(i, size_) = G(1.0L);
        }
        for (size_type j = 0; j < size_; ++j) {
            m_(size_, j) = _apex[j];
        }
        m_(size_, size_) = G(1.0L);
        return determinant(std::move(m_));
    }

    template< typename vertices >
    G
    orientation(vertices const & _vertices, size_type const _apex) const
    {
        return orientation(_vertices, points_.at(_apex));
    }

    G
    abs(G const & _x) const
    {
        return (_x < G(0.0L)) ? -_x : _x;
    }

    G
    abs(G && _x) const
    {
        return (_x < G(0.0L)) ? -std::move(_x) : std::move(_x);
    }

    G
    steal_best(point_list & _from, point_list & _to) const
    {
        auto it = _from.begin();
        auto const end = _from.end();
        G orientation_ = orientation(_to, *it);
        auto furthest = it;
        while (++it != end) {
            G const o_ = orientation(_to, *it);
            if (abs(orientation_) < abs(o_)) {
                orientation_ = o_;
                furthest = it;
            }
        }
        if (!(G(0.0L) < abs(orientation_))) {
            throw bad_geometry("can't find linearly independent point");
        }
        _to.splice(_to.end(), _from, furthest);
        return orientation_;
    }

    using ranking_type = std::multimap< G, size_type >;
    using ranking_meta_type = std::map< size_type, typename ranking_type::iterator >;
    ranking_type ranking_;
    ranking_meta_type ranking_meta_;

    void
    rank(G const _orientation, size_type const _facet)
    {
        if (G(0.0L) < _orientation) {
            auto const r = ranking_.emplace(_orientation, _facet);
            ranking_meta_.emplace(_facet, r);
        }
    }

    void
    unrank(size_type const _facet)
    {
        auto const r = ranking_meta_.find(_facet);
        if (r != ranking_meta_.end()) {
            ranking_.erase(r->second);
            ranking_meta_.erase(r);
        }
    }

    size_type
    get_furthest(size_type const _bad_value) const
    {
        if (ranking_.empty()) {
            return _bad_value;
        } else {
            auto const r = std::prev(ranking_.cend());
            return r->second;
        }
    }

    G
    partition(facet & _facet, point_list & _points)
    {
        auto it = _points.begin();
        auto const end = _points.end();
        points_type & outside_set_ = _facet.outside_set_;
        G orientation_(0.0L);
        while (it != end) {
            auto const next = std::next(it);
            G const o_ = orientation(_facet.vertices_, *it);
            if (below(_facet, o_)) {
                if (outside_set_.empty() || _facet.further(orientation_, o_)) {
                    orientation_ = o_;
                    outside_set_.push_front(*it);
                } else {
                    outside_set_.push_back(*it);
                }
                _points.erase(it);
            }
            it = next;
        }
        return abs(orientation_);
    }

    struct counter
            : std::iterator< std::output_iterator_tag, void, void, void, void >
    {

        counter(size_type & _counter)
            : counter_(_counter)
        { ; }

        counter &
        operator ++ ()
        {
            ++counter_;
            return *this;
        }

        counter
        operator ++ (int)
        {
            return operator ++ ();
        }

        counter &
        operator * ()
        {
            return *this;
        }

        template< typename T >
        counter &
        operator = (T &&)
        {
            return *this;
        }

    private :

        size_type & counter_;

    };

    void
    create_simplex()
    {
        {
            size_type const size_ = points_.size();
            assert(dimension_ < size_);
            for (size_type i = 0; i < size_; ++i) {
                internal_set_.push_back(i);
            }
        }
        point_list vertices_;
        vertices_.splice(vertices_.end(), internal_set_, internal_set_.begin());
        for (size_type i = 0; i < dimension_; ++i) {
            steal_best(internal_set_, vertices_);
        }
        assert(vertices_.size() == 1 + dimension_); // (N + 1) vertices defining a simplex
        internal_set_.splice(internal_set_.end(), vertices_, vertices_.begin());
        assert(vertices_.size() == dimension_); // N vertices defining a facet
        boolean_type outward_ = !(G(0.0L) < steal_best(internal_set_, vertices_)); // is top oriented?
        auto const vbeg = vertices_.cbegin();
        auto const vend = vertices_.cend();
#ifndef NDEBUG
        point_type inner_point_ = points_.at(*vbeg);
        {
            auto it = vbeg;
            while (++it != vend) {
                inner_point_ += points_.at(*it);
            }
            inner_point_ /= G(1 + dimension_);
        }
#endif
        auto const fend = facets_.end();
        for (auto exclusive = vend; exclusive != vbeg; --exclusive) {
            size_type const n = facets_.size();
            auto const f = facets_.emplace_hint(fend, n, facet(vbeg, exclusive, vend, outward_));
            facet & facet_ = f->second;
            rank(partition(facet_, internal_set_), n);
            assert(outward_ == !(G(0.0L) < orientation(facet_.vertices_, inner_point_)));
            outward_ = !outward_;
        }
        assert(dimension_ + 1 == facets_.size());
        {
            auto const fbeg = facets_.begin();
            for (auto i = fbeg; i != fend; ++i) {
                facet_set & neighbours_ = i->second.neighbours_;
                for (auto j = fbeg; j != fend; ++j) {
                    if (j != i) {
                        neighbours_.emplace_hint(neighbours_.end(), j->first);
                    }
                }
            }
        }
    }

    void
    create_convex_hull()
    {
        create_simplex();
        size_type facet_key = facets_.size(); // unique key for facets_
        assert(facet_key == dimension_ + 1);
        auto const fend = facets_.end();
        for (size_type furthest = get_furthest(facet_key); furthest != facet_key; furthest = get_furthest(facet_key)) {
            facet & facet_ = facets_.at(furthest);
            facet_set visible_facets_{furthest};
            size_type const apex = facet_.outside_set_.front();
            facet_.outside_set_.pop_front();
            {
                facet_set pool_ = facet_.neighbours_;
                facet_set visited_{furthest};
                while (!pool_.empty()) {
                    auto const first = pool_.begin();
                    size_type const f = *first;
                    facet const & facet_ = facets_.at(f);
                    if (below(facet_, orientation(facet_.vertices_, apex))) {
                        visible_facets_.insert(f);
                        std::set_difference(facet_.neighbours_.cbegin(), facet_.neighbours_.cend(),
                                            visited_.cbegin(), visited_.cend(),
                                            std::inserter(pool_, pool_.end()));
                    }
                    visited_.insert(f);
                    pool_.erase(first);
                }
            }
            // the boundary of visible facets is the set of horizon ridges
            // Each ridge signifies the adjacency of two facets.
            facets_type newfacets_;
            auto const vfend = visible_facets_.end();
            for (size_type const v : visible_facets_) {
                facet const & visible_facet_ = facets_.at(v);
                points_type const & vertices_ = visible_facet_.vertices_;
                for (size_type const n : visible_facet_.neighbours_) {
                    if (visible_facets_.find(n) == vfend) { // facets intersection with keeping of points order
                        facet & horizon_facet_ = facets_.at(n);
                        point_set horizon_(horizon_facet_.vertices_.cbegin(),
                                           horizon_facet_.vertices_.cend()); // n * log(n) +
                        auto const hend = horizon_.end();
                        points_type ridge_; // horizon ridge and furthest point
                        for (size_type const p : vertices_) {
                            auto const h = horizon_.find(p);
                            if (h == hend) {
                                ridge_.push_back(apex);
                            } else {
                                ridge_.push_back(p);
                                horizon_.erase(h);
                            }
                        }
                        assert(horizon_.size() == 1); // horizon_ contains the only invisible point beyond the horizon
                        assert(ridge_.size() == dimension_); // ridge_ contains newfacet vertices (ridge + current furthest point)
                        { // replace visible facet became internal with newly created facet in adjacency
                            horizon_facet_.neighbours_.erase(v);
                            horizon_facet_.neighbours_.insert(horizon_facet_.neighbours_.cend(), facet_key);
                        }
                        newfacets_.push_back(facet_key);
                        facets_.emplace_hint(fend, facet_key, facet(std::move(ridge_), visible_facet_.outward_, n));
                        ++facet_key;
                    }
                }
            }
            {
                auto const nend = newfacets_.end();
                for (auto first = newfacets_.begin(); first != nend; ++first) {
                    size_type const f = *first;
                    facet & first_ = facets_.at(f);
                    points_type const & ofirst_ = first_.points_;
                    for (auto second = std::next(first); second != nend; ++second) {
                        size_type const s = *second;
                        facet & second_ = facets_.at(s);
                        points_type const & osecond_ = second_.points_;
                        size_type count_ = 0;
                        std::set_difference(ofirst_.cbegin(), ofirst_.cend(),
                                            osecond_.cbegin(), osecond_.cend(),
                                            counter{count_});
                        if (count_ == 1) {
                            first_.neighbours_.insert(s);
                            second_.neighbours_.insert(f);
                        }
                    }
                }
            }
            point_list outside_set_;
            for (size_type const v : visible_facets_) {
                auto const visible_facet = facets_.find(v);
                assert(visible_facet != fend);
                facet const & visible_facet_ = visible_facet->second;
                outside_set_.insert(outside_set_.cend(),
                                    visible_facet_.outside_set_.cbegin(),
                                    visible_facet_.outside_set_.cend());
                facets_.erase(visible_facet);
                unrank(v);
            }
            for (size_type const n : newfacets_) {
                rank(partition(facets_.at(n), outside_set_), n);
            }
            internal_set_.splice(internal_set_.cend(), outside_set_);
        }
    }

};

And a test for generating an output for gnuplot (D1, D2 or D3):

#include "quickhull.hpp"

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <valarray>

#include <cstdlib>
#include <cstdio>

int main()
{
    using G = double;
    using point_type = std::valarray< G >;
    using H = convex_hull< point_type >;

    std::ifstream ifs_;
    ifs_.open("points.txt"); // rbox n D3 s 100 > points.txt
    if (!ifs_.is_open()) {
        std::cerr << "file is not open" << std::endl;
        return EXIT_FAILURE;
    }
    std::string line_;
    if (!std::getline(ifs_, line_)) {
        std::cerr << "no dim at first line" << std::endl;
        return EXIT_FAILURE;
    }
    size_type const dim_ = std::stoll(line_);
    if (!std::getline(ifs_, line_)) {
        std::cerr << "no count at second line" << std::endl;
        return EXIT_FAILURE;
    }
    size_type const count_ = std::stoll(line_);
    std::deque< point_type > points_;
    while (std::getline(ifs_, line_)) {
        points_.emplace_back(dim_);
        point_type & point_ = points_.back();
        std::istringstream iss(line_);
        for (size_type i = 0; i < dim_; ++i) {
            if (!(iss >> point_[i])) {
                std::cerr << "bad value faced at line " << points_.size() << " of data" << std::endl;
                return EXIT_FAILURE;
            }
        }
    }
    assert(points_.size() == count_);
    if (count_ != points_.size()) {
        std::cerr << "input file format error" << std::endl;
        return EXIT_FAILURE;
    }
    std::ofstream ofs_;
    ofs_.open("script.txt"); // gnuplot> load 'script.txt'
    if (!ofs_.is_open()) {
        std::cerr << "output file cannot be truncated" << std::endl;
        return EXIT_FAILURE;
    }
    ofs_ << "reset" << std::endl;
    ofs_ << "set autoscale" << std::endl;
    switch (dim_) {
    case 1 : {
        ofs_ << "set xrange [-0.5:0.5];" << std::endl;
        ofs_ << "plot";
        break;
    }
    case 2 : {
        ofs_ << "set size square; set xrange [-0.5:0.5]; set yrange [-0.5:0.5];" << std::endl;
        ofs_ << "plot";
        break;
    }
    case 3 : {
        ofs_ << "set view equal xyz; set view 0,0; set xrange [-0.5:0.5]; set yrange [-0.5:0.5]; set zrange [-0.5:0.5]; set xyplane at 0.0" << std::endl;
        ofs_ << "splot";
        break;
    }
    default : {
        std::cerr << "dimensionality value (" << dim_ << ") is out of supported range" << std::endl;
        return EXIT_FAILURE;
    }
    }
    std::cout.rdbuf()->pubsetbuf(nullptr, 0);
    std::cout << "D = " << dim_ << std::endl;
    std::cout << "N = " << count_ << std::endl;
    H convex_hull_(points_.cbegin(), points_.cend());
    convex_hull_.create_convex_hull();
    //convex_hull_.create_simplex();
    //return EXIT_SUCCESS;
    ofs_ << " '-' with points, '-' with labels offset 0,char 1";
    for (std::size_t i = 0; i < convex_hull_.facets_.size(); ++i) {
        ofs_ << ", '-' with lines notitle";
    }
    ofs_ << std::endl;
    for (point_type const & point_ : points_) {
        for (G const & coordinate_ : point_) {
            ofs_ << coordinate_ << ' ';
        }
        ofs_ << std::endl;
    }
    ofs_ << 'e' << std::endl;
    std::size_t i = 0;
    for (point_type const & point_ : points_) {
        for (G const & coordinate_ : point_) {
            ofs_ << coordinate_ << ' ';
        }
        ofs_ << i << std::endl;
        ++i;
    }
    ofs_ << 'e' << std::endl;
    for (auto const & facet_ : convex_hull_.facets_) {
        auto const & vertices_ = facet_.second.vertices_;
        std::cout << "facet #" << facet_.first << std::endl;
        for (size_type const vertex_ : vertices_) {
            for (G const & coordinate_ : points_.at(vertex_)) {
                ofs_ << coordinate_ << ' ';
            }
            ofs_ << std::endl;
        }
        point_type const & first_vertex_ = points_.at(vertices_.front());
        for (G const & coordinate_ : first_vertex_) {
            ofs_ << coordinate_ << ' ';
        }
        ofs_ << std::endl;
        ofs_ << 'e' << std::endl;
    }
    return EXIT_SUCCESS;
}

The correctness tested visually (for not too large set of points and for vaious combinations of rbox utility specific keys (D3: cube, diamond, cospherical points, helix, simplex, cone, etc. in all meaningfull combinations).

An example of output for rbox D2 n 33:

rbox D2 n 33 > points.txt && quickhull

An example of output for rbox D3 n 33:

not too visual

I tried to compare the time for large set of points against original qhull utility, and I faced significant comparative retardation (milliseconds vs seconds). I can't find rapidly: where my implementation is appreciably different from qhull (and from paper C. Bradford Barber et al. "The Quickhull Algorithm for Convex Hulls")? Maybe the degradation is only due to dynamic character of memory management (rather than qhull's one)?

Improve this question
\$\endgroup\$
1
  • \$\begingroup\$ My greater mistake is to calculate of point orientation relative to hyperplane (det() calculation 2/3*d<sup>3</sup> for each point) instead of simple calculation of hyperplane equation (d + 1 matrices of d * d size) and further only 2 * d floating-point operations to inner product calculation for each point. \$\endgroup\$
    – Tomilov Anatoliy
    Jun 15, 2014 at 14:28

1 Answer 1

Reset to default
3
\$\begingroup\$

None of this exception is needed:

struct bad_geometry
        : std::exception
{

    virtual
    ~bad_geometry() noexcept = default;

    bad_geometry() = default;

    explicit
    bad_geometry(const char * const _what)
        : what_(_what)
    { ; }

    explicit
    bad_geometry(std::string const & _what)
        : what_(_what)
    { ; }

    virtual
    const char *
    what() const noexcept override
    {
        return what_.c_str();
    }

private :

    std::string const what_ = "bad_geometry";

};

Replace with:

struct bad_geometry
        : std::runtime_error
{
    explicit
    bad_geometry(char const* what)
        : std::rutime_error(what)
    {}

    explicit
    bad_geometry(std::string const& what)
        : std::runtime_error(what.c_str())
    {}
};

Pass by value:

template< typename F >
F determinant(boost::numeric::ublas::matrix< F > _m)

Here _m is being passed by value. Thus it is probably generating a copy of the matrix. Passs by const reference to avoid the copy.

template< typename F >
F determinant(boost::numeric::ublas::matrix<F> const& m)

Lets stick all the using clauses together:

using points_type = std::deque< size_type >;

/// LOTS OF STUFF

struct facet;

using facet_set = std::set< size_type >;  // Stick this with the others.

If your code is inlined (like this); Then its easier to read a class if you put the variables first. Then I can read the constructors with a context of what values need to be set in what order:

// These defined way after the constructors use them.

    points_type vertices_; // oriented
    points_type points_;   // non-oriented
    boolean_type outward_;
    facet_set neighbours_;
    points_type outside_set_; // if not empty, then first point is furthest for this facet

This one is a take it or leave it (not a big deal).

   // This is fine.
        if (outward_) {
            return (_nearer < _further);
        } else {
            return (_further < _nearer);
        }

   // But I would write it like this:

        return (outward_)
             ? (_nearer < _further)
             : (_further < _nearer);

Your code is very dense (and thus hard to read). A bit more white space is not going to kill you and some explanation. Its hard to tell if there is any more unnecessary copying because of your use of auto (I hope you only use it for things where the type does not matter (iterators and the like)).

Improve this answer
\$\endgroup\$
1
  • \$\begingroup\$ Thank you for the reply. Especially first clause is usefull. Many others is matter of dispute. Namely, WRT matrix copying: LU-decomposition is modifying operation, so I should provide a copy, or, actually, non-const reference to it, because matrix is not used hereinafter in algorithm). WRT using of the ternary operator: in my opinion my version is more expressive. auto keyword used only for (I assume that they are "light-weight") iterators. I will correct the defects. Very grateful. \$\endgroup\$
    – Tomilov Anatoliy
    Jun 14, 2014 at 4:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.