GCC Code Coverage Report


Directory: ./
File: src/2geom/d2-sbasis.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 48 204 23.5%
Functions: 5 24 20.8%
Branches: 40 330 12.1%

Line Branch Exec Source
1 /**
2 * \file
3 * \brief Some two-dimensional SBasis operations
4 *//*
5 * Authors:
6 * MenTaLguy <mental@rydia.net>
7 * Jean-François Barraud <jf.barraud@gmail.com>
8 * Johan Engelen <j.b.c.engelen@alumnus.utwente.nl>
9 *
10 * Copyright 2007-2012 Authors
11 *
12 * This library is free software; you can redistribute it and/or
13 * modify it either under the terms of the GNU Lesser General Public
14 * License version 2.1 as published by the Free Software Foundation
15 * (the "LGPL") or, at your option, under the terms of the Mozilla
16 * Public License Version 1.1 (the "MPL"). If you do not alter this
17 * notice, a recipient may use your version of this file under either
18 * the MPL or the LGPL.
19 *
20 * You should have received a copy of the LGPL along with this library
21 * in the file COPYING-LGPL-2.1; if not, output to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 * You should have received a copy of the MPL along with this library
24 * in the file COPYING-MPL-1.1
25 *
26 * The contents of this file are subject to the Mozilla Public License
27 * Version 1.1 (the "License"); you may not use this file except in
28 * compliance with the License. You may obtain a copy of the License at
29 * http://www.mozilla.org/MPL/
30 *
31 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
32 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
33 * the specific language governing rights and limitations.
34 *
35 */
36
37 #include <2geom/d2.h>
38 #include <2geom/piecewise.h>
39
40 namespace Geom {
41
42 SBasis L2(D2<SBasis> const & a, unsigned k) { return sqrt(dot(a, a), k); }
43
44 D2<SBasis> multiply(Linear const & a, D2<SBasis> const & b) {
45 return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));
46 }
47
48 21066 D2<SBasis> multiply(SBasis const & a, D2<SBasis> const & b) {
49
3/6
✓ Branch 3 taken 21066 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 21066 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 21066 times.
✗ Branch 12 not taken.
42132 return D2<SBasis>(multiply(a, b[X]), multiply(a, b[Y]));
50 }
51
52 D2<SBasis> truncate(D2<SBasis> const & a, unsigned terms) {
53 return D2<SBasis>(truncate(a[X], terms), truncate(a[Y], terms));
54 }
55
56 unsigned sbasis_size(D2<SBasis> const & a) {
57 return std::max((unsigned) a[0].size(), (unsigned) a[1].size());
58 }
59
60 //TODO: Is this sensical? shouldn't it be like pythagorean or something?
61 double tail_error(D2<SBasis> const & a, unsigned tail) {
62 return std::max(a[0].tailError(tail), a[1].tailError(tail));
63 }
64
65 Piecewise<D2<SBasis> > sectionize(D2<Piecewise<SBasis> > const &a) {
66 Piecewise<SBasis> x = partition(a[0], a[1].cuts), y = partition(a[1], a[0].cuts);
67 assert(x.size() == y.size());
68 Piecewise<D2<SBasis> > ret;
69 for(unsigned i = 0; i < x.size(); i++)
70 ret.push_seg(D2<SBasis>(x[i], y[i]));
71 ret.cuts.insert(ret.cuts.end(), x.cuts.begin(), x.cuts.end());
72 return ret;
73 }
74
75 300 D2<Piecewise<SBasis> > make_cuts_independent(Piecewise<D2<SBasis> > const &a) {
76 300 D2<Piecewise<SBasis> > ret;
77
2/2
✓ Branch 0 taken 600 times.
✓ Branch 1 taken 300 times.
900 for(unsigned d = 0; d < 2; d++) {
78
2/2
✓ Branch 1 taken 24000 times.
✓ Branch 2 taken 600 times.
24600 for(unsigned i = 0; i < a.size(); i++)
79
1/2
✓ Branch 4 taken 24000 times.
✗ Branch 5 not taken.
24000 ret[d].push_seg(a[i][d]);
80
1/2
✓ Branch 9 taken 600 times.
✗ Branch 10 not taken.
600 ret[d].cuts.insert(ret[d].cuts.end(), a.cuts.begin(), a.cuts.end());
81 }
82 300 return ret;
83 }
84
85 300 Piecewise<D2<SBasis> > rot90(Piecewise<D2<SBasis> > const &M){
86 300 Piecewise<D2<SBasis> > result;
87
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 300 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
300 if (M.empty()) return M;
88
1/2
✓ Branch 2 taken 300 times.
✗ Branch 3 not taken.
300 result.push_cut(M.cuts[0]);
89
2/2
✓ Branch 1 taken 11700 times.
✓ Branch 2 taken 300 times.
12000 for (unsigned i=0; i<M.size(); i++){
90
2/4
✓ Branch 4 taken 11700 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 11700 times.
✗ Branch 8 not taken.
11700 result.push(rot90(M[i]),M.cuts[i+1]);
91 }
92 300 return result;
93 300 }
94
95 /** @brief Calculates the 'dot product' or 'inner product' of \c a and \c b
96 * @return \f[
97 * f(t) \rightarrow \left\{
98 * \begin{array}{c}
99 * a_1 \bullet b_1 \\
100 * a_2 \bullet b_2 \\
101 * \ldots \\
102 * a_n \bullet b_n \\
103 * \end{array}\right.
104 * \f]
105 * @relates Piecewise */
106 9600 Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Piecewise<D2<SBasis> > const &b)
107 {
108 9600 Piecewise<SBasis > result;
109
3/6
✓ Branch 1 taken 9600 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 9600 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 9600 times.
9600 if (a.empty() || b.empty()) return result;
110
1/2
✓ Branch 2 taken 9600 times.
✗ Branch 3 not taken.
9600 Piecewise<D2<SBasis> > aa = partition(a,b.cuts);
111
1/2
✓ Branch 2 taken 9600 times.
✗ Branch 3 not taken.
9600 Piecewise<D2<SBasis> > bb = partition(b,a.cuts);
112
113
1/2
✓ Branch 2 taken 9600 times.
✗ Branch 3 not taken.
9600 result.push_cut(aa.cuts.front());
114
2/2
✓ Branch 1 taken 9600 times.
✓ Branch 2 taken 9600 times.
19200 for (unsigned i=0; i<aa.size(); i++){
115
2/4
✓ Branch 5 taken 9600 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 9600 times.
✗ Branch 9 not taken.
9600 result.push(dot(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);
116 }
117 9600 return result;
118 9600 }
119
120 /** @brief Calculates the 'dot product' or 'inner product' of \c a and \c b
121 * @return \f[
122 * f(t) \rightarrow \left\{
123 * \begin{array}{c}
124 * a_1 \bullet b \\
125 * a_2 \bullet b \\
126 * \ldots \\
127 * a_n \bullet b \\
128 * \end{array}\right.
129 * \f]
130 * @relates Piecewise */
131 Piecewise<SBasis> dot(Piecewise<D2<SBasis> > const &a, Point const &b)
132 {
133 Piecewise<SBasis > result;
134 if (a.empty()) return result;
135
136 result.push_cut(a.cuts.front());
137 for (unsigned i = 0; i < a.size(); ++i){
138 result.push(dot(a.segs[i],b), a.cuts[i+1]);
139 }
140 return result;
141 }
142
143
144 Piecewise<SBasis> cross(Piecewise<D2<SBasis> > const &a,
145 Piecewise<D2<SBasis> > const &b){
146 Piecewise<SBasis > result;
147 if (a.empty() || b.empty()) return result;
148 Piecewise<D2<SBasis> > aa = partition(a,b.cuts);
149 Piecewise<D2<SBasis> > bb = partition(b,a.cuts);
150
151 result.push_cut(aa.cuts.front());
152 for (unsigned i=0; i<a.size(); i++){
153 result.push(cross(aa.segs[i],bb.segs[i]),aa.cuts[i+1]);
154 }
155 return result;
156 }
157
158 Piecewise<D2<SBasis> > operator*(Piecewise<D2<SBasis> > const &a, Affine const &m) {
159 Piecewise<D2<SBasis> > result;
160 if(a.empty()) return result;
161 result.push_cut(a.cuts[0]);
162 for (unsigned i = 0; i < a.size(); i++) {
163 result.push(a[i] * m, a.cuts[i+1]);
164 }
165 return result;
166 }
167
168 //if tol>0, only force continuity where the jump is smaller than tol.
169 300 Piecewise<D2<SBasis> > force_continuity(Piecewise<D2<SBasis> > const &f, double tol, bool closed)
170 {
171
1/4
✗ Branch 1 not taken.
✓ Branch 2 taken 300 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
300 if (f.size()==0) return f;
172
1/2
✓ Branch 2 taken 300 times.
✗ Branch 3 not taken.
300 Piecewise<D2<SBasis> > result=f;
173
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 300 times.
300 unsigned cur = (closed)? 0:1;
174
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 300 times.
300 unsigned prev = (closed)? f.size()-1:0;
175
2/2
✓ Branch 1 taken 11400 times.
✓ Branch 2 taken 300 times.
11700 while(cur<f.size()){
176 11400 Point pt0 = f.segs[prev].at1();
177 11400 Point pt1 = f.segs[cur ].at0();
178
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 11400 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 11400 times.
✓ Branch 10 taken 11400 times.
✗ Branch 11 not taken.
11400 if (tol<=0 || L2sq(pt0-pt1)<tol*tol){
179 11400 pt0 = (pt0+pt1)/2;
180
2/2
✓ Branch 0 taken 22800 times.
✓ Branch 1 taken 11400 times.
34200 for (unsigned dim=0; dim<2; dim++){
181 22800 SBasis &prev_sb=result.segs[prev][dim];
182 22800 SBasis &cur_sb =result.segs[cur][dim];
183 22800 Coord const c=pt0[dim];
184
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 22800 times.
22800 if (prev_sb.isZero(0)) {
185 prev_sb = SBasis(Linear(0.0, c));
186 } else {
187
1/2
✓ Branch 1 taken 22800 times.
✗ Branch 2 not taken.
22800 prev_sb[0][1] = c;
188 }
189
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 22800 times.
22800 if (cur_sb.isZero(0)) {
190 cur_sb = SBasis(Linear(c, 0.0));
191 } else {
192
1/2
✓ Branch 1 taken 22800 times.
✗ Branch 2 not taken.
22800 cur_sb[0][0] = c;
193 }
194 }
195 }
196 11400 prev = cur++;
197 }
198 300 return result;
199 300 }
200
201 std::vector<Geom::Piecewise<Geom::D2<Geom::SBasis> > >
202 split_at_discontinuities (Geom::Piecewise<Geom::D2<Geom::SBasis> > const & pwsbin, double tol)
203 {
204 using namespace Geom;
205 std::vector<Piecewise<D2<SBasis> > > ret;
206 unsigned piece_start = 0;
207 for (unsigned i=0; i<pwsbin.segs.size(); i++){
208 if (i==(pwsbin.segs.size()-1) || L2(pwsbin.segs[i].at1()- pwsbin.segs[i+1].at0()) > tol){
209 Piecewise<D2<SBasis> > piece;
210 piece.cuts.push_back(pwsbin.cuts[piece_start]);
211 for (unsigned j = piece_start; j<i+1; j++){
212 piece.segs.push_back(pwsbin.segs[j]);
213 piece.cuts.push_back(pwsbin.cuts[j+1]);
214 }
215 ret.push_back(piece);
216 piece_start = i+1;
217 }
218 }
219 return ret;
220 }
221
222 Point unitTangentAt(D2<SBasis> const & a, Coord t, unsigned n)
223 {
224 std::vector<Point> derivs = a.valueAndDerivatives(t, n);
225 for (unsigned deriv_n = 1; deriv_n < derivs.size(); deriv_n++) {
226 Coord length = derivs[deriv_n].length();
227 if ( ! are_near(length, 0) ) {
228 // length of derivative is non-zero, so return unit vector
229 return derivs[deriv_n] / length;
230 }
231 }
232 return Point (0,0);
233 }
234
235 static void set_first_point(Piecewise<D2<SBasis> > &f, Point const &a){
236 if ( f.empty() ){
237 f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(SBasis(Linear(a[X])), SBasis(Linear(a[Y])))));
238 return;
239 }
240 for (unsigned dim=0; dim<2; dim++){
241 f.segs.front()[dim][0][0] = a[dim];
242 }
243 }
244 static void set_last_point(Piecewise<D2<SBasis> > &f, Point const &a){
245 if ( f.empty() ){
246 f.concat(Piecewise<D2<SBasis> >(D2<SBasis>(SBasis(Linear(a[X])), SBasis(Linear(a[Y])))));
247 return;
248 }
249 for (unsigned dim=0; dim<2; dim++){
250 f.segs.back()[dim][0][1] = a[dim];
251 }
252 }
253
254 std::vector<Piecewise<D2<SBasis> > > fuse_nearby_ends(std::vector<Piecewise<D2<SBasis> > > const &f, double tol){
255
256 if ( f.empty()) return f;
257 std::vector<Piecewise<D2<SBasis> > > result;
258 std::vector<std::vector<unsigned> > pre_result;
259 for (unsigned i=0; i<f.size(); i++){
260 bool inserted = false;
261 Point a = f[i].firstValue();
262 Point b = f[i].lastValue();
263 for (auto & j : pre_result){
264 Point aj = f.at(j.back()).lastValue();
265 Point bj = f.at(j.front()).firstValue();
266 if ( L2(a-aj) < tol ) {
267 j.push_back(i);
268 inserted = true;
269 break;
270 }
271 if ( L2(b-bj) < tol ) {
272 j.insert(j.begin(),i);
273 inserted = true;
274 break;
275 }
276 }
277 if (!inserted) {
278 pre_result.emplace_back();
279 pre_result.back().push_back(i);
280 }
281 }
282 for (auto & i : pre_result){
283 Piecewise<D2<SBasis> > comp;
284 for (unsigned j=0; j<i.size(); j++){
285 Piecewise<D2<SBasis> > new_comp = f.at(i[j]);
286 if ( j>0 ){
287 set_first_point( new_comp, comp.segs.back().at1() );
288 }
289 comp.concat(new_comp);
290 }
291 if ( L2(comp.firstValue()-comp.lastValue()) < tol ){
292 //TODO: check sizes!!!
293 set_last_point( comp, comp.segs.front().at0() );
294 }
295 result.push_back(comp);
296 }
297 return result;
298 }
299
300 /*
301 * Computes the intersection of two sets given as (ordered) union of intervals.
302 */
303 static std::vector<Interval> intersect( std::vector<Interval> const &a, std::vector<Interval> const &b){
304 std::vector<Interval> result;
305 //TODO: use order!
306 for (auto i : a){
307 for (auto j : b){
308 OptInterval c( i );
309 c &= j;
310 if ( c ) {
311 result.push_back( *c );
312 }
313 }
314 }
315 return result;
316 }
317
318 std::vector<Interval> level_set( D2<SBasis> const &f, Rect region){
319 std::vector<Rect> regions( 1, region );
320 return level_sets( f, regions ).front();
321 }
322 std::vector<Interval> level_set( D2<SBasis> const &f, Point p, double tol){
323 Rect region(p, p);
324 region.expandBy( tol );
325 return level_set( f, region );
326 }
327 std::vector<std::vector<Interval> > level_sets( D2<SBasis> const &f, std::vector<Rect> regions){
328 std::vector<Interval> regsX (regions.size(), Interval() );
329 std::vector<Interval> regsY (regions.size(), Interval() );
330 for ( unsigned i=0; i < regions.size(); i++ ){
331 regsX[i] = regions[i][X];
332 regsY[i] = regions[i][Y];
333 }
334 std::vector<std::vector<Interval> > x_in_regs = level_sets( f[X], regsX );
335 std::vector<std::vector<Interval> > y_in_regs = level_sets( f[Y], regsY );
336 std::vector<std::vector<Interval> >result(regions.size(), std::vector<Interval>() );
337 for (unsigned i=0; i<regions.size(); i++){
338 result[i] = intersect ( x_in_regs[i], y_in_regs[i] );
339 }
340 return result;
341 }
342 std::vector<std::vector<Interval> > level_sets( D2<SBasis> const &f, std::vector<Point> pts, double tol){
343 std::vector<Rect> regions( pts.size(), Rect() );
344 for (unsigned i=0; i<pts.size(); i++){
345 regions[i] = Rect( pts[i], pts[i] );
346 regions[i].expandBy( tol );
347 }
348 return level_sets( f, regions );
349 }
350
351
352 } // namespace Geom
353
354
355 /*
356 Local Variables:
357 mode:c++
358 c-file-style:"stroustrup"
359 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
360 indent-tabs-mode:nil
361 fill-column:99
362 End:
363 */
364 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
365