From a9874ef5ffaa2db8ed4792b91560f67cd03a34ec Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 23 Jan 2018 18:10:31 +0100
Subject: [PATCH] [geometry] Add convex hull of point cloud + triangulation and
test
---
dumux/common/geometry/grahamconvexhull.hh | 163 ++++++++++++++++++
test/common/geometry/CMakeLists.txt | 5 +-
.../geometry/test_graham_convex_hull.cc | 125 ++++++++++++++
test/common/geometry/writetriangulation.hh | 84 +++++++++
4 files changed, 376 insertions(+), 1 deletion(-)
create mode 100644 dumux/common/geometry/grahamconvexhull.hh
create mode 100644 test/common/geometry/test_graham_convex_hull.cc
create mode 100644 test/common/geometry/writetriangulation.hh
diff --git a/dumux/common/geometry/grahamconvexhull.hh b/dumux/common/geometry/grahamconvexhull.hh
new file mode 100644
index 0000000000..67ef578a83
--- /dev/null
+++ b/dumux/common/geometry/grahamconvexhull.hh
@@ -0,0 +1,163 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief A function to compute the convex hull of a point cloud
+ * and a function to triangulate the polygon area spanned by the convex hull
+ */
+#ifndef DUMUX_GRAHAM_CONVEX_HULL_HH
+#define DUMUX_GRAHAM_CONVEX_HULL_HH
+
+#include <vector>
+#include <array>
+#include <algorithm>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+
+#include <dumux/common/math.hh>
+
+namespace Dumux {
+
+/*!
+ * \brief Returns the orientation of a sequence a-->b-->c in one plane (defined by normal vector)
+ * \return -1 if a-->b-->c forms a counter-clockwise turn (given the normal vector)
+ * +1 for a clockwise turn,
+ * 0 if they are on one line (colinear)
+ */
+int getOrientation(const Dune::FieldVector<double, 3>& a,
+ const Dune::FieldVector<double, 3>& b,
+ const Dune::FieldVector<double, 3>& c,
+ const Dune::FieldVector<double, 3>& normal)
+{
+ const auto d = b-a;
+ const auto e = c-b;
+ const auto f = Dumux::crossProduct(d, e);
+ const auto area = f*normal;
+ return Dumux::sign(-area);
+}
+
+/*!
+ * \brief Compute the points making up the convex hull around the given set of unordered points
+ * \note We assume that all points are coplanar
+ * \note this algorithm changes the order of the given points a bit
+ * as they are unordered anyway this shouldn't matter too much
+ */
+std::vector<Dune::FieldVector<double, 3>>
+grahamConvexHull2d3d(std::vector<Dune::FieldVector<double, 3>>& points)
+{
+ using Point = Dune::FieldVector<double, 3>;
+ std::vector<Point> convexHull; convexHull.reserve(50);
+
+ // return empty convex hull
+ if (points.size() < 3)
+ return convexHull;
+
+ // return the points (already just one triangle)
+ if (points.size() == 3)
+ return points;
+
+ // get the normal vector of the plane
+ const auto a = points[1] - points[0];
+ const auto b = points[2] - points[0];
+ auto normal = Dumux::crossProduct(a, b);
+ normal /= normal.two_norm();
+
+ // find the element with the smallest y coordinate (if y is the same, smallest x coordinate)
+ auto minIt = std::min_element(points.begin(), points.end(),
+ [](const auto& a, const auto& b)
+ { return a[1] != b[1] ? a[1] < b[1] : a[0] < b[0]; });
+
+ // swap the smallest element to the front
+ std::iter_swap(minIt, points.begin());
+
+ // choose the first (min element) as the pivot point
+ // sort in counter-clockwise order around pivot point
+ const auto pivot = points[0];
+ std::sort(points.begin()+1, points.end(), [&](const auto& a, const auto& b)
+ {
+ int order = getOrientation(pivot, a, b, normal);
+ if (order == 0)
+ return (a-pivot).two_norm() < (b-pivot).two_norm();
+ else
+ return (order == -1);
+ });
+
+ // push the first three points
+ convexHull.push_back(points[0]);
+ convexHull.push_back(points[1]);
+ convexHull.push_back(points[2]);
+
+ // This is the heart of the algorithm
+ // Pop_back until the last point in the queue forms a counter-clockwise oriented line
+ // with the new vertices. Then add new points to the queue.
+ for (std::size_t i = 3; i < points.size(); ++i)
+ {
+ Point p = convexHull.back();
+ convexHull.pop_back();
+ // keep popping until the orientation a->b->currentp is counter-clockwise
+ while (getOrientation(convexHull.back(), p, points[i], normal) != -1)
+ {
+ p = convexHull.back();
+ convexHull.pop_back();
+ }
+
+ // add back the last popped point and this point
+ convexHull.emplace_back(std::move(p));
+ convexHull.push_back(points[i]);
+ }
+
+ return convexHull;
+}
+
+/*!
+ * \brief Triangulate area given points of the convex hull
+ * \note Assumes all points of the convex hull are coplanar
+ * \note This inserts a mid point and connects all corners with that point to triangles
+ */
+std::vector<std::array<Dune::FieldVector<double, 3>, 3> >
+triangulateConvexHull(const std::vector<Dune::FieldVector<double, 3>>& convexHull)
+{
+ using Point = Dune::FieldVector<double, 3>;
+ using Triangle = std::array<Point, 3>;
+
+ if (convexHull.size() < 3)
+ DUNE_THROW(Dune::InvalidStateException, "Try to triangulate point cloud with less than 3 points!");
+
+ if (convexHull.size() == 3)
+ return std::vector<Triangle>(1, {convexHull[0], convexHull[1], convexHull[2]});
+
+ Point midPoint(0.0);
+ for (const auto p : convexHull)
+ midPoint += p;
+ midPoint /= convexHull.size();
+
+ std::vector<Triangle> triangulation;
+ triangulation.reserve(convexHull.size());
+
+ for (std::size_t i = 0; i < convexHull.size()-1; ++i)
+ triangulation.emplace_back(Triangle{midPoint, convexHull[i], convexHull[i+1]});
+
+ triangulation.emplace_back(Triangle{midPoint, convexHull[convexHull.size()-1], convexHull[0]});
+
+ return triangulation;
+}
+
+} // end namespace Dumux
+
+# endif
diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
index 41465b6623..2d10a00d2e 100644
--- a/test/common/geometry/CMakeLists.txt
+++ b/test/common/geometry/CMakeLists.txt
@@ -1,8 +1,11 @@
dune_add_test(SOURCES test_1d3d_intersection.cc)
+dune_add_test(SOURCES test_graham_convex_hull.cc)
#install sources
install(FILES
test_1d3d_intersection.cc
+test_graham_convex_hull.cc
+writetriangulation.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/common/geometry)
-set(CMAKE_BUILD_TYPE Release)
+set(CMAKE_BUILD_TYPE RelWithDebInfo)
diff --git a/test/common/geometry/test_graham_convex_hull.cc b/test/common/geometry/test_graham_convex_hull.cc
new file mode 100644
index 0000000000..d08cf7f859
--- /dev/null
+++ b/test/common/geometry/test_graham_convex_hull.cc
@@ -0,0 +1,125 @@
+// graham convex hull test + triangulation
+#include <config.h>
+
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <random>
+#include <vector>
+#include <algorithm>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/timer.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/common/geometry/grahamconvexhull.hh>
+#include <test/common/geometry/writetriangulation.hh>
+
+template<class Scalar>
+class UniformDistributedRandomNumber
+{
+public:
+ UniformDistributedRandomNumber(Scalar min, Scalar max)
+ : gen_(r_()), dist_(min, max) {}
+
+ Scalar operator()()
+ { return dist_(gen_); }
+
+private:
+ std::random_device r_;
+ std::default_random_engine gen_;
+ std::uniform_real_distribution<Scalar> dist_;
+};
+
+Dune::FieldVector<double, 3>
+randomPointOnPlane(const Dune::FieldVector<double, 3>& origin,
+ const Dune::FieldVector<double, 3>& normal,
+ double min = -1.0, double max = 1.0)
+{
+ using Point = Dune::FieldVector<double, 3>;
+
+ UniformDistributedRandomNumber<double> dice(min, max);
+ Point p{dice(), dice(), dice()};
+ p = Dumux::crossProduct(normal, p);
+ p += origin;
+ return p;
+}
+
+void writeCSV(const std::vector<Dune::FieldVector<double, 3>>& p,
+ const std::string& filename)
+{
+ std::ofstream fout(filename + ".csv");
+ fout << "X, Y, Z\n";
+ for (const auto& point : p)
+ fout << point[0] << ", " << point[1] << ", " << point[2] << '\n';
+ fout << std::endl;
+}
+
+void writeVTKPolyData(const std::vector<Dune::FieldVector<double, 3>>& p,
+ const std::string& filename)
+{
+ std::ofstream fout(filename + ".vtp");
+ fout << "<?xml version=\"1.0\"?>\n"
+ << "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\">\n"
+ << " <PolyData>\n"
+ << " <Piece NumberOfPoints=\"" << p.size() << "\" NumberOfLines=\"" << p.size() << "\">\n"
+ << " <Points>\n"
+ << " <DataArray type=\"Float32\" Name=\"Coordinates\" NumberOfComponents=\"3\" format=\"ascii\">\n";
+
+ for (const auto& point : p)
+ fout << point << " ";
+ fout << '\n';
+
+ fout << " </DataArray>\n"
+ << " </Points>\n"
+ << " <Lines>\n"
+ << " <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+
+ for (std::size_t i = 0; i < p.size()-1; ++i)
+ fout << i << " " << i+1 << " ";
+ fout << p.size()-1 << " 0\n";
+
+ fout << " </DataArray>\n";
+ fout << " <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+
+ for (std::size_t i = 1; i <= p.size(); ++i)
+ fout << i*2 << " ";
+ fout << '\n';
+
+ fout << " </DataArray>\n"
+ << " </Lines>\n"
+ << " </Piece>\n"
+ << "</PolyData>\n"
+ << "</VTKFile>\n";
+}
+
+int main(int argc, char* argv[])
+{
+ using Point = Dune::FieldVector<double, 3>;
+ // create 100 random points
+ std::vector<Point> points(50);
+ // normal vector of the plane
+ const Point normal({0.0, 1.0, 1.0});
+ const Point origin({3.0, 0.0, 0.0});
+ for (auto&& p : points)
+ p = randomPointOnPlane(origin, normal, -10.0, 10.0);
+
+ writeCSV(points, "points");
+
+ Dune::Timer timer;
+ auto convexHullPoints = Dumux::grahamConvexHull2d3d(points);
+ std::cout << "Computed convex hull of " << points.size() << " points in "
+ << timer.elapsed() << " seconds." << std::endl;
+
+ writeVTKPolyData(convexHullPoints, "convexhull");
+
+ Dune::Timer timer2;
+ auto triangles = Dumux::triangulateConvexHull(convexHullPoints);
+ std::cout << "Computed triangulation of convex hull with " << convexHullPoints.size()
+ << " points in " << timer2.elapsed() << " seconds." << std::endl;
+
+ Dumux::writeVTKPolyDataTriangle(triangles, "triangulation");
+
+ return 0;
+}
diff --git a/test/common/geometry/writetriangulation.hh b/test/common/geometry/writetriangulation.hh
new file mode 100644
index 0000000000..4aa8e184a0
--- /dev/null
+++ b/test/common/geometry/writetriangulation.hh
@@ -0,0 +1,84 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief A function to write a triangulation to vtk
+ */
+#ifndef DUMUX_TEST_WRITE_TRIANGULATION_HH
+#define DUMUX_TEST_WRITE_TRIANGULATION_HH
+
+#include <fstream>
+#include <string>
+#include <iterator>
+
+namespace Dumux {
+
+//! TriangleVector has to be a nested vector of triangles in 3d
+template<class TriangleVector>
+void writeVTKPolyDataTriangle(const TriangleVector& triangles,
+ const std::string& filename)
+{
+ std::ofstream fout(filename + ".vtp");
+ fout << "<?xml version=\"1.0\"?>\n"
+ << "<VTKFile type=\"PolyData\" version=\"0.1\" byte_order=\"LittleEndian\">\n"
+ << " <PolyData>\n"
+ << " <Piece NumberOfPoints=\"" << triangles.size()*3 << "\" NumberOfLines=\"" << triangles.size()*3 << "\">\n"
+ << " <Points>\n"
+ << " <DataArray type=\"Float32\" Name=\"Coordinates\" NumberOfComponents=\"3\" format=\"ascii\">\n";
+
+ for (const auto& t : triangles)
+ for (const auto& p : t)
+ fout << p << " ";
+ fout << '\n';
+
+ fout << " </DataArray>\n"
+ << " </Points>\n"
+ << " <Lines>\n"
+ << " <DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+
+ int offset = 0;
+ for (const auto& t : triangles)
+ {
+ for (std::size_t i = 0; i < t.size()-1; ++i)
+ fout << i + offset*3 << " " << i+1 + offset*3 << " ";
+ fout << t.size()-1 + offset*3 << " " << offset*3 << " ";
+ ++offset;
+ }
+
+ fout << " </DataArray>\n";
+ fout << " <DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+
+ offset = 0;
+ for (const auto& t : triangles)
+ {
+ for (std::size_t i = 1; i <= t.size(); ++i)
+ fout << i*2 + offset*6 << " ";
+ ++offset;
+ }
+ fout << '\n';
+
+ fout << " </DataArray>\n"
+ << " </Lines>\n"
+ << " </Piece>\n"
+ << "</PolyData>\n"
+ << "</VTKFile>\n";
+}
+
+} // end namespace Dumux
+
+# endif
--
GitLab