GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/geometry/grahamconvexhull.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 68 69 98.6%
Functions: 5 5 100.0%
Branches: 44 66 66.7%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Geometry
10 * \brief A function to compute the convex hull of a point cloud
11 */
12 #ifndef DUMUX_GEOMETRY_GRAHAM_CONVEX_HULL_HH
13 #define DUMUX_GEOMETRY_GRAHAM_CONVEX_HULL_HH
14
15 #include <vector>
16 #include <array>
17 #include <algorithm>
18 #include <iterator>
19
20 #include <dune/common/exceptions.hh>
21 #include <dune/common/fvector.hh>
22
23 #include <dumux/common/math.hh>
24
25 namespace Dumux {
26
27 /*!
28 * \ingroup Geometry
29 * \brief Returns the orientation of a sequence a-->b-->c in one plane (defined by normal vector)
30 * \return -1 if a-->b-->c forms a counter-clockwise turn (given the normal vector)
31 * +1 for a clockwise turn,
32 * 0 if they are on one line (colinear)
33 */
34 template<class ctype>
35 648127 int getOrientation(const Dune::FieldVector<ctype, 3>& a,
36 const Dune::FieldVector<ctype, 3>& b,
37 const Dune::FieldVector<ctype, 3>& c,
38 const Dune::FieldVector<ctype, 3>& normal)
39 {
40 648127 const auto d = b-a;
41 648127 const auto e = c-b;
42 648127 const auto f = Dumux::crossProduct(d, e);
43 648127 const auto area = f*normal;
44 648127 return Dumux::sign(-area);
45 }
46
47 /*!
48 * \ingroup Geometry
49 * \brief Compute the points making up the convex hull around the given set of unordered points
50 * \note We assume that all points are coplanar and there are no identical points in the list
51 * \note This algorithm changes the order of the given points a bit
52 * as they are unordered anyway this shouldn't matter too much
53 */
54 template<int dim, class ctype,
55 std::enable_if_t<(dim==2), int> = 0>
56 std::vector<Dune::FieldVector<ctype, 3>>
57 115956 grahamConvexHullImpl(std::vector<Dune::FieldVector<ctype, 3>>& points)
58 {
59 using Point = Dune::FieldVector<ctype, 3>;
60
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 115956 times.
 115956 std::vector<Point> convexHull;
61
62 // return empty convex hull
63
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 115956 times.
 115956 if (points.size() < 3)
64 ✗ return convexHull;
65
66 // return the points (already just one triangle)
67
2/2
✓ Branch 0 taken 3950 times.
✓ Branch 1 taken 112006 times.
 115956 if (points.size() == 3)
68
1/2
✓ Branch 1 taken 3950 times.
✗ Branch 2 not taken.
 3950 return points;
69
70 // try to compute the normal vector of the plane
71 112006 const auto a = points[1] - points[0];
72 112006 auto b = points[2] - points[0];
73 112006 auto normal = Dumux::crossProduct(a, b);
74
75 // make sure the normal vector has non-zero length
76 112006 std::size_t k = 2;
77 112006 auto norm = normal.two_norm();
78
4/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 112005 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 1 times.
 112010 while (norm == 0.0 && k < points.size()-1)
79 {
80 4 b = points[++k];
81 4 normal = Dumux::crossProduct(a, b);
82 4 norm = normal.two_norm();
83 }
84
85 // if all given points are colinear -> return empty convex hull
86
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 112005 times.
 112006 if (norm == 0.0)
87 1 return convexHull;
88
89 using std::sqrt;
90 112005 const auto eps = 1e-7*sqrt(norm);
91 112005 normal /= norm;
92
93 // find the element with the smallest x coordinate (if x is the same, smallest y coordinate, and so on...)
94
2/2
✓ Branch 0 taken 213591 times.
✓ Branch 1 taken 125293 times.
 450889 auto minIt = std::min_element(points.begin(), points.end(), [&eps](const auto& a, const auto& b)
95 {
96 using std::abs;
97
4/4
✓ Branch 0 taken 213591 times.
✓ Branch 1 taken 125293 times.
✓ Branch 2 taken 111261 times.
✓ Branch 3 taken 14032 times.
 338884 return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : (abs(a[1]-b[1]) > eps ? a[1] < b[1] : (a[2] < b[2])));
98 });
99
100 // swap the smallest element to the front
101 112005 std::iter_swap(minIt, points.begin());
102
103 // choose the first (min element) as the pivot point
104 // sort in counter-clockwise order around pivot point
105 112005 const auto pivot = points[0];
106 112005 std::sort(points.begin()+1, points.end(), [&](const auto& a, const auto& b)
107 {
108 522959 const int order = getOrientation(pivot, a, b, normal);
109
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 522948 times.
 522959 if (order == 0)
110 55 return (a-pivot).two_norm() < (b-pivot).two_norm();
111 else
112 522948 return (order == -1);
113 });
114
115 // push the first three points
116
1/2
✓ Branch 1 taken 112005 times.
✗ Branch 2 not taken.
 112005 convexHull.reserve(50);
117
1/2
✓ Branch 1 taken 112005 times.
✗ Branch 2 not taken.
 112005 convexHull.push_back(points[0]);
118
1/2
✓ Branch 1 taken 112005 times.
✗ Branch 2 not taken.
 112005 convexHull.push_back(points[1]);
119
1/2
✓ Branch 1 taken 112005 times.
✗ Branch 2 not taken.
 112005 convexHull.push_back(points[2]);
120
121 // This is the heart of the algorithm
122 // pop_back until the last point in the queue forms a counter-clockwise oriented line
123 // with the next two vertices. Then add these points to the queue.
124
2/2
✓ Branch 0 taken 114872 times.
✓ Branch 1 taken 112005 times.
 226877 for (std::size_t i = 3; i < points.size(); ++i)
125 {
126 114872 Point p = convexHull.back();
127 114872 convexHull.pop_back();
128 // keep popping until the orientation a->b->currentp is counter-clockwise
129
2/2
✓ Branch 1 taken 96 times.
✓ Branch 2 taken 114872 times.
 114968 while (getOrientation(convexHull.back(), p, points[i], normal) != -1)
130 {
131 // make sure the queue doesn't get empty
132
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 94 times.
 96 if (convexHull.size() == 1)
133 {
134 // before we reach i=size-1 there has to be a good candidate
135 // as not all points are colinear (a non-zero plane normal exists)
136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 assert(i < points.size()-1);
137 2 p = points[i++];
138 }
139 else
140 {
141 94 p = convexHull.back();
142 94 convexHull.pop_back();
143 }
144 }
145
146 // add back the last popped point and this point
147
1/2
✓ Branch 1 taken 114872 times.
✗ Branch 2 not taken.
 114872 convexHull.emplace_back(std::move(p));
148
1/2
✓ Branch 1 taken 114872 times.
✗ Branch 2 not taken.
 114872 convexHull.push_back(points[i]);
149 }
150
151 112005 return convexHull;
152 115956 }
153
154 /*!
155 * \ingroup Geometry
156 * \brief Compute the points making up the convex hull around the given set of unordered points
157 * \note This is the specialization for 2d space. Here, we make use of the generic implementation
158 * for the case of coplanar points in 3d space (a more efficient implementation could be provided).
159 */
160 template<int dim, class ctype,
161 std::enable_if_t<(dim==2), int> = 0>
162 std::vector<Dune::FieldVector<ctype, 2>>
163 78952 grahamConvexHullImpl(const std::vector<Dune::FieldVector<ctype, 2>>& points)
164 {
165
1/2
✓ Branch 1 taken 78952 times.
✗ Branch 2 not taken.
 78952 std::vector<Dune::FieldVector<ctype, 3>> points3D;
166
1/2
✓ Branch 1 taken 78952 times.
✗ Branch 2 not taken.
 78952 points3D.reserve(points.size());
167
1/2
✓ Branch 1 taken 78952 times.
✗ Branch 2 not taken.
 78952 std::transform(points.begin(), points.end(), std::back_inserter(points3D),
168 315772 [](const auto& p) { return Dune::FieldVector<ctype, 3>({p[0], p[1], 0.0}); });
169
170
1/2
✓ Branch 1 taken 78952 times.
✗ Branch 2 not taken.
 78952 const auto result3D = grahamConvexHullImpl<2>(points3D);
171
172
1/2
✓ Branch 1 taken 78952 times.
✗ Branch 2 not taken.
 78952 std::vector<Dune::FieldVector<ctype, 2>> result2D;
173
1/2
✓ Branch 1 taken 78952 times.
✗ Branch 2 not taken.
 78952 result2D.reserve(result3D.size());
174
1/2
✓ Branch 1 taken 78952 times.
✗ Branch 2 not taken.
 78952 std::transform(result3D.begin(), result3D.end(), std::back_inserter(result2D),
175 315772 [](const auto& p) { return Dune::FieldVector<ctype, 2>({p[0], p[1]}); });
176
177 78952 return result2D;
178 157904 }
179
180 /*!
181 * \ingroup Geometry
182 * \brief Compute the points making up the convex hull around the given set of unordered points
183 * \note We assume that all points are coplanar and there are no identical points in the list
184 */
185 template<int dim, class ctype, int dimWorld>
186 92473 std::vector<Dune::FieldVector<ctype, dimWorld>> grahamConvexHull(std::vector<Dune::FieldVector<ctype, dimWorld>>& points)
187 {
188
4/8
✓ Branch 1 taken 89689 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 102 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2245 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 437 times.
✗ Branch 11 not taken.
 92473 return grahamConvexHullImpl<dim>(points);
189 }
190
191 /*!
192 * \ingroup Geometry
193 * \brief Compute the points making up the convex hull around the given set of unordered points
194 * \note We assume that all points are coplanar and there are no identical points in the list
195 * \note This is the overload if we are not allowed to write into the given points vector
196 */
197 template<int dim, class ctype, int dimWorld>
198 23483 std::vector<Dune::FieldVector<ctype, dimWorld>> grahamConvexHull(const std::vector<Dune::FieldVector<ctype, dimWorld>>& points)
199 {
200 23483 auto copyPoints = points;
201
1/2
✓ Branch 1 taken 23483 times.
✗ Branch 2 not taken.
 23483 return grahamConvexHullImpl<dim>(copyPoints);
202 23483 }
203
204 } // end namespace Dumux
205
206 #endif
207