diff --git a/LICENSES b/LICENSES index 55d229be9649baee82d276953cf082da4c5252b3..33e0ef904d2ec34dc3f80ce8027c9892e42a5f8f 100644 --- a/LICENSES +++ b/LICENSES @@ -148,6 +148,18 @@ License: GNU GPL version 2 (see below) or FreeType License +-------------------------------------------------------------------------------- +Qhull +URL: http://www.qhull.org/ +URL: https://github.com/qhull/qhull + +Copyright (c) 1993-2019 by C.B. Barber + +License URL: https://github.com/qhull/qhull/blob/master/COPYING.txt +-------------------------------------------------------------------------------- + + + -------------------------------------------------------------------------------- Gnuplot URL: http://www.gnuplot.info diff --git a/math/geo2.h b/math/geo2.h new file mode 100644 index 0000000000000000000000000000000000000000..ac060d0993cb8eb583d32d740820716aaf14f3a3 --- /dev/null +++ b/math/geo2.h @@ -0,0 +1,121 @@ +/** + * advanced geometry helpers + * @author Tobias Weber + * @date 25-apr-2020 + * @license GPLv2 or GPLv3 + * + * @see: https://github.com/t-weber/misc/blob/master/geo/qhulltst.cpp + */ + +#ifndef __TLIBS_GEO2_H__ +#define __TLIBS_GEO2_H__ + +#include "geo.h" + +#ifndef NO_QHULL + #include + + #include + #include + #include +#endif + + +namespace tl { + +namespace ublas = boost::numeric::ublas; +namespace math = boost::math; + + + +#ifndef NO_QHULL + +using t_real_qh = double; + +template, +template class t_cont = std::vector, +class T = typename t_vec::value_type> +t_cont> get_convexhull_qh(const t_cont& vecVerts) +{ + t_cont> vecPolys; + const t_vec vecCentre = mean_value(vecVerts); + + // copy vertices + int dim = vecVerts[0].size(); + std::size_t len = vecVerts.size()*dim; + std::unique_ptr mem{new t_real_qh[len]}; + + std::size_t i=0; + for(const t_vec& vert : vecVerts) + { + for(int d=0; d::iterator iter=facets.begin(); + iter!=facets.end(); ++iter) + { + // triangulable? + if(iter->isUpperDelaunay()) + continue; + + t_cont vecPoly; + orgQhull::QhullVertexSet vertices = iter->vertices(); + for(orgQhull::QhullSet::iterator iterVertex=vertices.begin(); + iterVertex!=vertices.end(); ++iterVertex) + { + orgQhull::QhullPoint point = (*iterVertex).point(); + t_vec vecPoint(dim); + for(int i=0; i(vecPoly, vecCentre); + vecPolys.emplace_back(std::move(vecPoly)); + } + + // too few polygons => remove polyhedron + if(vecPolys.size() < 3) + vecPolys = decltype(vecPolys){}; + + return vecPolys; +} + +#endif + + + +/** + * creates polyhedron faces out of vertices + */ +template, +template class t_cont = std::vector, +class T = typename t_vec::value_type> +t_cont> get_convexhull( + const t_cont& vecVerts, T eps = tl::get_epsilon()) +{ +#ifdef NO_QHULL + // note that his function is more restrictive than a convex hull: + // it searches for polyhedra and does not permit additional vertices + // inside the polyhedron + return verts_to_polyhedron(vecVerts, eps); +#else + return get_convexhull_qh(vecVerts); +#endif +} + + +} + +#endif diff --git a/tools/test/coordpoly.cpp b/tools/test/coordpoly.cpp index 3cc8cd0875e7e3bc0385aaa9c82ad69bb69ee8ce..a9ac522224d9fffc345e58ac5269f5a1c14b4731 100644 --- a/tools/test/coordpoly.cpp +++ b/tools/test/coordpoly.cpp @@ -4,52 +4,38 @@ * @license GPLv2 or GPLv3 */ -// gcc -DNO_LAPACK -o coordpoly coordpoly.cpp ../log/log.cpp ../file/x3d.cpp -std=c++11 -lstdc++ -lm +// g++ -std=c++11 -I. -I/usr/include/libqhullcpp -DNO_LAPACK -o coordpoly tools/test/coordpoly.cpp log/log.cpp file/x3d.cpp -lqhull_r -lqhullcpp #include #include -#include "../math/geo.h" -#include "../file/x3d.h" +#include "math/geo2.h" +#include "file/x3d.h" using T = double; using t_vec = tl::ublas::vector; -int main() + +void SaveX3d(const std::vector>& vecPolys, + const std::vector& vecVertices, const char* pcFile) { - std::ifstream ifstr("/home/tweber/tmp/coord.dat"); - std::vector vecVertices; t_vec vecCentral = tl::zero_v(3); - while(1) - { - t_vec vec(3); - ifstr >> vec[0] >> vec[1] >> vec[2]; - if(!ifstr) - break; - - vecVertices.push_back(vec); - } - - - tl::X3d x3d; + tl::X3d x3d; for(const t_vec& vec : vecVertices) { - tl::X3dTrafo *pTrafo = new tl::X3dTrafo(); + tl::X3dTrafo *pTrafo = new tl::X3dTrafo(); pTrafo->SetTrans(vec - vecCentral); - tl::X3dSphere *pSphere = new tl::X3dSphere(0.025); + tl::X3dSphere *pSphere = new tl::X3dSphere(0.025); pSphere->SetColor(tl::make_vec({1., 0., 0.})); pTrafo->AddChild(pSphere); x3d.GetScene().AddChild(pTrafo); } - - std::vector> vecPolys = tl::verts_to_polyhedron(vecVertices); - for(const auto& vecPoly : vecPolys) { - tl::X3dPolygon *pPoly = new tl::X3dPolygon(); + tl::X3dPolygon *pPoly = new tl::X3dPolygon(); //tl::X3dLines *pPoly = new tl::X3dLines(); pPoly->SetColor(tl::make_vec({1., 1., 0.})); @@ -59,9 +45,30 @@ int main() x3d.GetScene().AddChild(pPoly); } + x3d.Save(pcFile); +} + + +int main() +{ + std::ifstream ifstr("/tmp/coord.dat"); + std::vector vecVertices; + + while(1) + { + t_vec vec(3); + ifstr >> vec[0] >> vec[1] >> vec[2]; + if(!ifstr) + break; + + vecVertices.push_back(vec); + } + std::vector> vecPolys = tl::verts_to_polyhedron(vecVertices, 1e-4); + std::vector> vecPolys_qh = tl::get_convexhull_qh(vecVertices); - x3d.Save("/home/tweber/tmp/coord.x3d"); + SaveX3d(vecPolys, vecVertices, "/tmp/coord.x3d"); + SaveX3d(vecPolys_qh, vecVertices, "/tmp/coord_qh.x3d"); return 0; }