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 >Fdeterminant(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:An example of output for rbox D3 n 33: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)?
Quick Hull Algorithm implementation for higher dimensions
c++;performance;algorithm;computational geometry
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 STUFFstruct 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 facetThis 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)).