hull calculation

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

math/geo2.h

+/**
+ * advanced geometry helpers
+ * @author Tobias Weber <tobias.weber@tum.de>
+ * @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 <memory>
+
+	#include <Qhull.h>
+	#include <QhullFacetList.h>
+	#include <QhullVertexSet.h>
+#endif
+
+
+namespace tl {
+
+namespace ublas = boost::numeric::ublas;
+namespace math = boost::math;
+
+
+
+#ifndef NO_QHULL
+
+using t_real_qh = double;
+
+template<class t_vec = ublas::vector<double>,
+template<class...> class t_cont = std::vector,
+class T = typename t_vec::value_type>
+t_cont<t_cont<t_vec>> get_convexhull_qh(const t_cont<t_vec>& vecVerts)
+{
+	t_cont<t_cont<t_vec>> 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<t_real_qh[]> mem{new t_real_qh[len]};
+
+	std::size_t i=0;
+	for(const t_vec& vert : vecVerts)
+	{
+		for(int d=0; d<dim; ++d)
+		{
+			mem[i] = t_real_qh(vert[d]);
+			++i;
+		}
+	}
+
+
+	orgQhull::Qhull qhull{"tlibs_geo2", dim, int(vecVerts.size()), mem.get(), ""};
+
+	orgQhull::QhullFacetList facets = qhull.facetList();
+	//std::cout << "Convex hull has " << facets.size() << " facets.";
+
+	for(orgQhull::QhullLinkedList<orgQhull::QhullFacet>::iterator iter=facets.begin();
+		iter!=facets.end(); ++iter)
+	{
+		// triangulable?
+		if(iter->isUpperDelaunay())
+			continue;
+
+		t_cont<t_vec> vecPoly;
+		orgQhull::QhullVertexSet vertices = iter->vertices();
+		for(orgQhull::QhullSet<orgQhull::QhullVertex>::iterator iterVertex=vertices.begin();
+			iterVertex!=vertices.end(); ++iterVertex)
+		{
+			orgQhull::QhullPoint point = (*iterVertex).point();
+			t_vec vecPoint(dim);
+			for(int i=0; i<dim; ++i)
+				vecPoint[i] = T(point[i]);
+
+			vecPoly.emplace_back(std::move(vecPoint));
+		}
+
+		sort_poly_verts<t_vec, t_cont, T>(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<class t_vec = ublas::vector<double>,
+template<class...> class t_cont = std::vector,
+class T = typename t_vec::value_type>
+t_cont<t_cont<t_vec>> get_convexhull(
+	const t_cont<t_vec>& vecVerts, T eps = tl::get_epsilon<T>())
+{
+#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<t_vec, t_cont, T>(vecVerts, eps);
+#else
+	return get_convexhull_qh<t_vec, t_cont, T>(vecVerts);
+#endif
+}
+
+
+}
+
+#endif

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 <iostream>
 #include <fstream>
-#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<T>;
 
-int main()
+
+void SaveX3d(const std::vector<std::vector<t_vec>>& vecPolys,
+	const std::vector<t_vec>& vecVertices, const char* pcFile)
 {
-	std::ifstream ifstr("/home/tweber/tmp/coord.dat");
-	std::vector<t_vec> vecVertices;
 	t_vec vecCentral = tl::zero_v<t_vec>(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<T> x3d;
 	for(const t_vec& vec : vecVertices)
 	{
-		tl::X3dTrafo *pTrafo = new tl::X3dTrafo();
+		tl::X3dTrafo<T> *pTrafo = new tl::X3dTrafo<T>();
 		pTrafo->SetTrans(vec - vecCentral);
 
-		tl::X3dSphere *pSphere = new tl::X3dSphere(0.025);
+		tl::X3dSphere<T> *pSphere = new tl::X3dSphere<T>(0.025);
 		pSphere->SetColor(tl::make_vec({1., 0., 0.}));
 		pTrafo->AddChild(pSphere);
 
 		x3d.GetScene().AddChild(pTrafo);
 	}
 
-
-	std::vector<std::vector<t_vec>> vecPolys = tl::verts_to_polyhedron<t_vec, std::vector, T>(vecVertices);
-
 	for(const auto& vecPoly : vecPolys)
 	{
-		tl::X3dPolygon *pPoly = new tl::X3dPolygon();
+		tl::X3dPolygon<T> *pPoly = new tl::X3dPolygon<T>();
 		//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<t_vec> vecVertices;
+
+	while(1)
+	{
+		t_vec vec(3);
+		ifstr >> vec[0] >> vec[1] >> vec[2];
+		if(!ifstr)
+			break;
+
+		vecVertices.push_back(vec);
+	}
 
+	std::vector<std::vector<t_vec>> vecPolys = tl::verts_to_polyhedron<t_vec, std::vector, T>(vecVertices, 1e-4);
+	std::vector<std::vector<t_vec>> vecPolys_qh = tl::get_convexhull_qh<t_vec, std::vector, T>(vecVertices);
 
-	x3d.Save("/home/tweber/tmp/coord.x3d");
+	SaveX3d(vecPolys, vecVertices, "/tmp/coord.x3d");
+	SaveX3d(vecPolys_qh, vecVertices, "/tmp/coord_qh.x3d");
 
 	return 0;
 }