GCC Code Coverage Report

Directory: Bembel/src/
File: Bembel/src/Spline/DeBoor.hpp
Date: 2024-03-19 14:38:05
Exec Total Coverage
Lines: 50 50 100.0%
Functions: 2 2 100.0%
Branches: 37 48 77.1%
Line Branch Exec Source
1 // This file is part of Bembel, the higher order C++ boundary element library.
2 //
3 // Copyright (C) 2022 see <http://www.bembel.eu>
4 //
5 // It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
6 // M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
7 // Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
8 // source code is subject to the GNU General Public License version 3 and
9 // provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
10 // information.
11 #ifndef BEMBEL_SRC_SPLINE_DEBOOR_HPP_
12 #define BEMBEL_SRC_SPLINE_DEBOOR_HPP_
13
14 namespace Bembel {
15 namespace Spl {
16 /**
17 * \ingroup Spline
18 * \brief "By the book" implementations of the Cox-DeBoor formula. Inefficient,
19 * do not use at bottlenecks.
20 **/
21 template <typename T>
22 64254 Eigen::Matrix<T, -1, -1> DeBoor(
23 Eigen::Matrix<T, -1, -1> const &control_points,
24 const std::vector<double> &knot_vector,
25 const std::vector<double> &evaluation_points) noexcept {
26 64254 const int cols_control_points = control_points.cols();
27 64254 const int rows_control_points = control_points.rows();
28 64254 const int polynomial_degree = knot_vector.size() - cols_control_points - 1;
29 64254 const int size_evaluation_points = evaluation_points.size();
30
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 64254 times.
 64254 assert(polynomial_degree >= 0);
31 64254 Eigen::Matrix<T, -1, -1> out(rows_control_points, size_evaluation_points);
32
33
2/2
✓ Branch 0 taken 64740 times.
✓ Branch 1 taken 64254 times.
 128994 for (int j = size_evaluation_points - 1; j >= 0; j--) {
34 64740 const double &x = evaluation_points[j];
35
36
2/4
✓ Branch 1 taken 64740 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64740 times.
✗ Branch 5 not taken.
 64740 assert(x >= knot_vector[polynomial_degree] &&
37 x <= knot_vector[cols_control_points + polynomial_degree]);
38 64740 int l = 0;
39 {
40
6/6
✓ Branch 1 taken 923394 times.
✓ Branch 2 taken 59108 times.
✓ Branch 3 taken 917762 times.
✓ Branch 4 taken 5632 times.
✓ Branch 5 taken 917762 times.
✓ Branch 6 taken 64740 times.
 982502 while (knot_vector[l] <= x && l != cols_control_points) l++;
41 64740 l = l - 1;
42 }
43
44
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 64740 times.
 64740 assert(l >= 0);
45 64740 Eigen::Matrix<T, -1, -1> temp = control_points.block(
46 64740 0, l - polynomial_degree, rows_control_points, polynomial_degree + 1);
47
48 /// ISO C++ forbids variable length array
49
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 64740 times.
 64740 assert(polynomial_degree <= Bembel::Constants::MaxP);
50 double ws[Bembel::Constants::MaxP];
51 // Iterators remain the same size, order is correct
52
2/2
✓ Branch 0 taken 852302 times.
✓ Branch 1 taken 64740 times.
 917042 for (int k = polynomial_degree; 0 != k; k--) {
53
2/2
✓ Branch 0 taken 6894742 times.
✓ Branch 1 taken 852302 times.
 7747044 for (int i = 0; i < k; i++) {
54 6894742 ws[i] = (x - knot_vector[l - k + 1 + i]) /
55 6894742 (knot_vector[l + 1 + i] - knot_vector[l - k + 1 + i]);
56 }
57
2/2
✓ Branch 0 taken 6894742 times.
✓ Branch 1 taken 852302 times.
 7747044 for (int i = 0; i < k; i++) {
58 6894742 temp.col(i) =
59 13789484 ((ws[i] * temp.col(1 + i)) + ((1 - ws[i]) * temp.col(i))).eval();
60 }
61 }
62 64740 out.col(j) = temp.col(0);
63 }
64 64254 return out;
65 }
66
67 template <typename T>
68 474292 std::vector<T> DeBoor(std::vector<T> const &control_points,
69 const std::vector<double> &knot_vector,
70 const std::vector<double> &evaluation_points) noexcept {
71 474292 const int cols_control_points = control_points.size();
72 474292 const int polynomial_degree = knot_vector.size() - cols_control_points - 1;
73 474292 const int size_evaluation_points = evaluation_points.size();
74
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 474292 times.
 474292 assert(polynomial_degree >= 0);
75 474292 std::vector<T> out(size_evaluation_points);
76
77
2/2
✓ Branch 0 taken 1963180 times.
✓ Branch 1 taken 474292 times.
 2437472 for (int j = size_evaluation_points - 1; j >= 0; j--) {
78 1963180 const double &x = evaluation_points[j];
79
80
2/4
✓ Branch 1 taken 1963180 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1963180 times.
✗ Branch 5 not taken.
 1963180 assert(x >= knot_vector[polynomial_degree] &&
81 x <= knot_vector[cols_control_points + polynomial_degree]);
82 1963180 int l = 0;
83 {
84
5/6
✓ Branch 1 taken 14012456 times.
✓ Branch 2 taken 1963180 times.
✓ Branch 3 taken 14012456 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 14012456 times.
✓ Branch 6 taken 1963180 times.
 15975636 while (knot_vector[l] <= x && l != cols_control_points) l++;
85 1963180 l = l - 1;
86 }
87
88
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1963180 times.
 1963180 assert(l >= 0);
89
90 3926360 std::vector<T> temp(control_points.begin() + l - polynomial_degree,
91 3926360 control_points.begin() + l + 1);
92
93 // ISO C++ forbids variable length array
94
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1963180 times.
 1963180 assert(polynomial_degree <= 19 &&
95 "Polynomial degree not supported in DeBoor.");
96 double ws[20];
97 // Iterators remain the same size, order is correct
98
2/2
✓ Branch 0 taken 6064604 times.
✓ Branch 1 taken 1963180 times.
 8027784 for (int k = polynomial_degree; 0 != k; k--) {
99
2/2
✓ Branch 0 taken 13261220 times.
✓ Branch 1 taken 6064604 times.
 19325824 for (int i = 0; i < k; i++) {
100 13261220 ws[i] = (x - knot_vector[l - k + 1 + i]) /
101 13261220 (knot_vector[l + 1 + i] - knot_vector[l - k + 1 + i]);
102 }
103
2/2
✓ Branch 0 taken 13261220 times.
✓ Branch 1 taken 6064604 times.
 19325824 for (int i = 0; i < k; i++) {
104 13261220 temp[i] = ((ws[i] * temp[1 + i]) + ((1 - ws[i]) * temp[i]));
105 }
106 }
107 1963180 out[j] = temp[0];
108 }
109 474292 return out;
110 }
111
112 template <typename T>
113 std::vector<Eigen::Matrix<T, -1, -1>> DeBoor(
114 std::vector<Eigen::Matrix<T, -1, -1>> const &control_points,
115 std::vector<double> const &knot_vector,
116 std::vector<double> const &evaluation_points) noexcept {
117 const int dimension = control_points.size();
118 const int cols_control_points = control_points[0].cols();
119 const int rows_control_points = control_points[0].rows();
120 const int polynomial_degree = knot_vector.size() - cols_control_points - 1;
121 const int size_evaluation_points = evaluation_points.size();
122 assert(polynomial_degree >= 0);
123 std::vector<Eigen::Matrix<T, -1, -1>> out(dimension);
124
125 for (int l = 0; l < dimension; l++)
126 out[l].resize(rows_control_points, size_evaluation_points);
127
128 for (int j = size_evaluation_points - 1; j >= 0; j--) {
129 const double &x = evaluation_points[j];
130
131 assert(x >= knot_vector[polynomial_degree] &&
132 x <= knot_vector[cols_control_points + polynomial_degree]);
133 int l = 0;
134 {
135 while (knot_vector[l] <= x && l != cols_control_points) l++;
136 l = l - 1;
137 }
138
139 assert(l >= 0);
140 std::vector<Eigen::Matrix<T, -1, -1>> temp(dimension);
141 for (int h = 0; h < dimension; h++)
142 temp[h] = control_points[h].block(
143 0, l - polynomial_degree, rows_control_points, polynomial_degree + 1);
144
145 // ISO C++ forbids variable length array
146 assert(polynomial_degree <= Bembel::Constants::MaxP);
147 double ws[Bembel::Constants::MaxP];
148 // Iterators remain the same size, order is correct
149 for (int k = polynomial_degree; 0 != k; k--) {
150 for (int i = 0; i < k; i++) {
151 ws[i] = (x - knot_vector[l - k + 1 + i]) /
152 (knot_vector[l + 1 + i] - knot_vector[l - k + 1 + i]);
153 }
154
155 for (int h = 0; h < dimension; h++) {
156 for (int i = 0; i < k; i++) {
157 temp[h].col(i) =
158 ((ws[i] * temp[h].col(1 + i)) + ((1 - ws[i]) * temp[h].col(i)))
159 .eval();
160 }
161 }
162 }
163 for (int h = 0; h < dimension; h++) out[h].col(j) = temp[h].col(0);
164 }
165 return out;
166 }
167 } // namespace Spl
168 } // namespace Bembel
169 #endif // BEMBEL_SRC_SPLINE_DEBOOR_HPP_
170