GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/geometry/grahamconvexhull.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 62 65 95.4%
Functions: 5 8 62.5%
Branches: 77 150 51.3%
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-FileCopyrightInfo: 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 599181 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 599181 const auto d = b-a;
41 599181 const auto e = c-b;
42 599181 const auto f = Dumux::crossProduct(d, e);
43 599181 const auto area = f*normal;
44 1198362 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 104028 grahamConvexHullImpl(std::vector<Dune::FieldVector<ctype, 3>>& points)
58 {
59 using Point = Dune::FieldVector<ctype, 3>;
60
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 104028 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 104028 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 208056 std::vector<Point> convexHull;
61
62 // return empty convex hull
63
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 104028 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 104028 times.
 208056 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 100078 times.
 104028 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 200156 const auto a = points[1] - points[0];
72 200156 auto b = points[2] - points[0];
73 100078 auto normal = Dumux::crossProduct(a, b);
74
75 // make sure the normal vector has non-zero length
76 100078 std::size_t k = 2;
77 200156 auto norm = normal.two_norm();
78
6/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 100077 times.
✓ Branch 2 taken 4 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 4 times.
✓ Branch 5 taken 1 times.
 100082 while (norm == 0.0 && k < points.size()-1)
79 {
80 4 b = points[++k];
81 4 normal = Dumux::crossProduct(a, b);
82 8 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 100077 times.
 100078 if (norm == 0.0)
87 1 return convexHull;
88
89 using std::sqrt;
90 100077 const auto eps = 1e-7*sqrt(norm);
91 100077 normal /= norm;
92
93 // find the element with the smallest x coordinate (if x is the same, smallest y coordinate, and so on...)
94 300231 auto minIt = std::min_element(points.begin(), points.end(), [&eps](const auto& a, const auto& b)
95 {
96 using std::abs;
97
16/32
✓ Branch 0 taken 198239 times.
✓ Branch 1 taken 104861 times.
✓ Branch 2 taken 198239 times.
✓ Branch 3 taken 104861 times.
✓ Branch 4 taken 198239 times.
✓ Branch 5 taken 104861 times.
✓ Branch 6 taken 198239 times.
✓ Branch 7 taken 104861 times.
✓ Branch 8 taken 98853 times.
✓ Branch 9 taken 6008 times.
✓ Branch 10 taken 98853 times.
✓ Branch 11 taken 6008 times.
✓ Branch 12 taken 98853 times.
✓ Branch 13 taken 6008 times.
✓ Branch 14 taken 98853 times.
✓ Branch 15 taken 6008 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
 1212400 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 200154 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 100077 const auto pivot = points[0];
106 400308 std::sort(points.begin()+1, points.end(), [&](const auto& a, const auto& b)
107 {
108 485943 const int order = getOrientation(pivot, a, b, normal);
109
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 485932 times.
 485943 if (order == 0)
110 33 return (a-pivot).two_norm() < (b-pivot).two_norm();
111 else
112 485932 return (order == -1);
113 });
114
115 // push the first three points
116
1/2
✓ Branch 1 taken 100077 times.
✗ Branch 2 not taken.
 100077 convexHull.reserve(50);
117
1/2
✓ Branch 1 taken 100077 times.
✗ Branch 2 not taken.
 100077 convexHull.push_back(points[0]);
118
2/4
✓ Branch 1 taken 100077 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100077 times.
✗ Branch 5 not taken.
 200154 convexHull.push_back(points[1]);
119
2/4
✓ Branch 1 taken 100077 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 100077 times.
✗ Branch 5 not taken.
 200154 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
4/4
✓ Branch 0 taken 102944 times.
✓ Branch 1 taken 100077 times.
✓ Branch 2 taken 102944 times.
✓ Branch 3 taken 100077 times.
 303098 for (std::size_t i = 3; i < points.size(); ++i)
125 {
126 102944 Point p = convexHull.back();
127 102944 convexHull.pop_back();
128 // keep popping until the orientation a->b->currentp is counter-clockwise
129
2/2
✓ Branch 3 taken 94 times.
✓ Branch 4 taken 102944 times.
 309114 while (getOrientation(convexHull.back(), p, points[i], normal) != -1)
130 {
131 // make sure the queue doesn't get empty
132
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 92 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 92 times.
 188 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
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 4 assert(i < points.size()-1);
137 4 p = points[i++];
138 }
139 else
140 {
141 92 p = convexHull.back();
142 92 convexHull.pop_back();
143 }
144 }
145
146 // add back the last popped point and this point
147
1/2
✓ Branch 1 taken 102944 times.
✗ Branch 2 not taken.
 102944 convexHull.emplace_back(std::move(p));
148
2/4
✓ Branch 1 taken 102944 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 102944 times.
✗ Branch 5 not taken.
 205888 convexHull.push_back(points[i]);
149 }
150
151 100077 return convexHull;
152 }
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 78820 grahamConvexHullImpl(const std::vector<Dune::FieldVector<ctype, 2>>& points)
164 {
165
1/4
✓ Branch 1 taken 78820 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 78820 std::vector<Dune::FieldVector<ctype, 3>> points3D;
166
2/4
✓ Branch 1 taken 78820 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 78820 times.
✗ Branch 5 not taken.
 157640 points3D.reserve(points.size());
167
4/8
✓ Branch 1 taken 78820 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 78820 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 78820 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 78820 times.
✗ Branch 11 not taken.
 315280 std::transform(points.begin(), points.end(), std::back_inserter(points3D),
168 ✗ [](const auto& p) { return Dune::FieldVector<ctype, 3>({p[0], p[1], 0.0}); });
169
170
3/8
✓ Branch 1 taken 78820 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 78820 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 78820 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 236460 const auto result3D = grahamConvexHullImpl<2>(points3D);
171
172
1/2
✓ Branch 1 taken 78820 times.
✗ Branch 2 not taken.
 78820 std::vector<Dune::FieldVector<ctype, 2>> result2D;
173
2/4
✓ Branch 1 taken 78820 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 78820 times.
✗ Branch 5 not taken.
 157640 result2D.reserve(result3D.size());
174
4/8
✓ Branch 1 taken 78820 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 78820 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 78820 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 78820 times.
✗ Branch 11 not taken.
 315280 std::transform(result3D.begin(), result3D.end(), std::back_inserter(result2D),
175 ✗ [](const auto& p) { return Dune::FieldVector<ctype, 2>({p[0], p[1]}); });
176
177 78820 return result2D;
178 }
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 std::vector<Dune::FieldVector<ctype, dimWorld>> grahamConvexHull(std::vector<Dune::FieldVector<ctype, dimWorld>>& points)
187 {
188
4/8
✓ Branch 1 taken 89077 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.
 91861 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 12167 std::vector<Dune::FieldVector<ctype, dimWorld>> grahamConvexHull(const std::vector<Dune::FieldVector<ctype, dimWorld>>& points)
199 {
200
1/4
✓ Branch 1 taken 12167 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 24334 auto copyPoints = points;
201
1/2
✓ Branch 1 taken 12167 times.
✗ Branch 2 not taken.
 24334 return grahamConvexHullImpl<dim>(copyPoints);
202 }
203
204 } // end namespace Dumux
205
206 #endif
207