GCC Code Coverage Report

Line	Branch	Exec	Source
1
2
3
4			#include <2geom/basic-intersection.h>
5			#include <2geom/sbasis-to-bezier.h>
6			#include <2geom/exception.h>
7
8
9			#include <gsl/gsl_vector.h>
10			#include <gsl/gsl_multiroots.h>
11
12
13			unsigned intersect_steps = 0;
14
15			using std::vector;
16
17			namespace Geom {
18
19			class OldBezier {
20			public:
21			std::vector<Geom::Point> p;
22		✗	OldBezier() {
23		✗	}
24			void split(double t, OldBezier &a, OldBezier &b) const;
25			Point operator()(double t) const;
26
27		✗	~OldBezier() {}
28
29		✗	void bounds(double &minax, double &maxax,
30			double &minay, double &maxay) {
31			// Compute bounding box for a
32		✗	minax = p[0][X]; // These are the most likely to be extremal
33		✗	maxax = p.back()[X];
34		✗	if( minax > maxax )
35		✗	std::swap(minax, maxax);
36		✗	for(unsigned i = 1; i < p.size()-1; i++) {
37		✗	if( p[i][X] < minax )
38		✗	minax = p[i][X];
39		✗	else if( p[i][X] > maxax )
40		✗	maxax = p[i][X];
41			}
42
43		✗	minay = p[0][Y]; // These are the most likely to be extremal
44		✗	maxay = p.back()[Y];
45		✗	if( minay > maxay )
46		✗	std::swap(minay, maxay);
47		✗	for(unsigned i = 1; i < p.size()-1; i++) {
48		✗	if( p[i][Y] < minay )
49		✗	minay = p[i][Y];
50		✗	else if( p[i][Y] > maxay )
51		✗	maxay = p[i][Y];
52			}
53
54		✗	}
55
56			};
57
58			static void
59			find_intersections_bezier_recursive(std::vector<std::pair<double, double> > & xs,
60			OldBezier a,
61			OldBezier b);
62
63			void
64		✗	find_intersections_bezier_recursive( std::vector<std::pair<double, double> > &xs,
65			vector<Geom::Point> const & A,
66			vector<Geom::Point> const & B,
67			double /*precision*/) {
68		✗	OldBezier a, b;
69		✗	a.p = A;
70		✗	b.p = B;
71		✗	return find_intersections_bezier_recursive(xs, a,b);
72		✗	}
73
74
75			/*
76			* split the curve at the midpoint, returning an array with the two parts
77			* Temporary storage is minimized by using part of the storage for the result
78			* to hold an intermediate value until it is no longer needed.
79			*/
80		✗	void OldBezier::split(double t, OldBezier &left, OldBezier &right) const {
81		✗	const unsigned sz = p.size();
82			//Geom::Point Vtemp[sz][sz];
83		✗	std::vector< std::vector< Geom::Point > > Vtemp;
84		✗	for (size_t i = 0; i < sz; ++i )
85		✗	Vtemp[i].reserve(sz);
86
87			/* Copy control points */
88		✗	std::copy(p.begin(), p.end(), Vtemp[0].begin());
89
90			/* Triangle computation */
91		✗	for (unsigned i = 1; i < sz; i++) {
92		✗	for (unsigned j = 0; j < sz - i; j++) {
93		✗	Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
94			}
95			}
96
97		✗	left.p.resize(sz);
98		✗	right.p.resize(sz);
99		✗	for (unsigned j = 0; j < sz; j++)
100		✗	left.p[j] = Vtemp[j][0];
101		✗	for (unsigned j = 0; j < sz; j++)
102		✗	right.p[j] = Vtemp[sz-1-j][j];
103		✗	}
104
105			#if 0
106			/*
107			* split the curve at the midpoint, returning an array with the two parts
108			* Temporary storage is minimized by using part of the storage for the result
109			* to hold an intermediate value until it is no longer needed.
110			*/
111			Point OldBezier::operator()(double t) const {
112			const unsigned sz = p.size();
113			Geom::Point Vtemp[sz][sz];
114
115			/* Copy control points */
116			std::copy(p.begin(), p.end(), Vtemp[0]);
117
118			/* Triangle computation */
119			for (unsigned i = 1; i < sz; i++) {
120			for (unsigned j = 0; j < sz - i; j++) {
121			Vtemp[i][j] = lerp(t, Vtemp[i-1][j], Vtemp[i-1][j+1]);
122			}
123			}
124			return Vtemp[sz-1][0];
125			}
126			#endif
127
128			// suggested by Sederberg.
129		✗	Point OldBezier::operator()(double const t) const {
130		✗	size_t const n = p.size()-1;
131		✗	Point r;
132		✗	for(int dim = 0; dim < 2; dim++) {
133		✗	double const u = 1.0 - t;
134		✗	double bc = 1;
135		✗	double tn = 1;
136		✗	double tmp = p[0][dim]*u;
137		✗	for(size_t i=1; i<n; i++){
138		✗	tn = tn*t;
139		✗	bc = bc*(n-i+1)/i;
140		✗	tmp = (tmp + tn*bc*p[i][dim])*u;
141			}
142		✗	r[dim] = (tmp + tn*t*p[n][dim]);
143			}
144		✗	return r;
145			}
146
147
148			/*
149			* Test the bounding boxes of two OldBezier curves for interference.
150			* Several observations:
151			* First, it is cheaper to compute the bounding box of the second curve
152			* and test its bounding box for interference than to use a more direct
153			* approach of comparing all control points of the second curve with
154			* the various edges of the bounding box of the first curve to test
155			* for interference.
156			* Second, after a few subdivisions it is highly probable that two corners
157			* of the bounding box of a given Bezier curve are the first and last
158			* control point. Once this happens once, it happens for all subsequent
159			* subcurves. It might be worth putting in a test and then short-circuit
160			* code for further subdivision levels.
161			* Third, in the final comparison (the interference test) the comparisons
162			* should both permit equality. We want to find intersections even if they
163			* occur at the ends of segments.
164			* Finally, there are tighter bounding boxes that can be derived. It isn't
165			* clear whether the higher probability of rejection (and hence fewer
166			* subdivisions and tests) is worth the extra work.
167			*/
168
169		✗	bool intersect_BB( OldBezier a, OldBezier b ) {
170		✗	double minax, maxax, minay, maxay;
171		✗	a.bounds(minax, maxax, minay, maxay);
172		✗	double minbx, maxbx, minby, maxby;
173		✗	b.bounds(minbx, maxbx, minby, maxby);
174			// Test bounding box of b against bounding box of a
175			// Not >= : need boundary case
176		✗	return !( ( minax > maxbx ) || ( minay > maxby )
177		✗	|| ( minbx > maxax ) || ( minby > maxay ) );
178			}
179
180			/*
181			* Recursively intersect two curves keeping track of their real parameters
182			* and depths of intersection.
183			* The results are returned in a 2-D array of doubles indicating the parameters
184			* for which intersections are found. The parameters are in the order the
185			* intersections were found, which is probably not in sorted order.
186			* When an intersection is found, the parameter value for each of the two
187			* is stored in the index elements array, and the index is incremented.
188			*
189			* If either of the curves has subdivisions left before it is straight
190			* (depth > 0)
191			* that curve (possibly both) is (are) subdivided at its (their) midpoint(s).
192			* the depth(s) is (are) decremented, and the parameter value(s) corresponding
193			* to the midpoints(s) is (are) computed.
194			* Then each of the subcurves of one curve is intersected with each of the
195			* subcurves of the other curve, first by testing the bounding boxes for
196			* interference. If there is any bounding box interference, the corresponding
197			* subcurves are recursively intersected.
198			*
199			* If neither curve has subdivisions left, the line segments from the first
200			* to last control point of each segment are intersected. (Actually the
201			* only the parameter value corresponding to the intersection point is found).
202			*
203			* The apriori flatness test is probably more efficient than testing at each
204			* level of recursion, although a test after three or four levels would
205			* probably be worthwhile, since many curves become flat faster than their
206			* asymptotic rate for the first few levels of recursion.
207			*
208			* The bounding box test fails much more frequently than it succeeds, providing
209			* substantial pruning of the search space.
210			*
211			* Each (sub)curve is subdivided only once, hence it is not possible that for
212			* one final line intersection test the subdivision was at one level, while
213			* for another final line intersection test the subdivision (of the same curve)
214			* was at another. Since the line segments share endpoints, the intersection
215			* is robust: a near-tangential intersection will yield zero or two
216			* intersections.
217			*/
218		✗	void recursively_intersect( OldBezier a, double t0, double t1, int deptha,
219			OldBezier b, double u0, double u1, int depthb,
220			std::vector<std::pair<double, double> > &parameters)
221			{
222		✗	intersect_steps ++;
223			//std::cout << deptha << std::endl;
224		✗	if( deptha > 0 )
225			{
226		✗	OldBezier A[2];
227		✗	a.split(0.5, A[0], A[1]);
228		✗	double tmid = (t0+t1)*0.5;
229		✗	deptha--;
230		✗	if( depthb > 0 )
231			{
232		✗	OldBezier B[2];
233		✗	b.split(0.5, B[0], B[1]);
234		✗	double umid = (u0+u1)*0.5;
235		✗	depthb--;
236		✗	if( intersect_BB( A[0], B[0] ) )
237		✗	recursively_intersect( A[0], t0, tmid, deptha,
238			B[0], u0, umid, depthb,
239			parameters );
240		✗	if( intersect_BB( A[1], B[0] ) )
241		✗	recursively_intersect( A[1], tmid, t1, deptha,
242			B[0], u0, umid, depthb,
243			parameters );
244		✗	if( intersect_BB( A[0], B[1] ) )
245		✗	recursively_intersect( A[0], t0, tmid, deptha,
246			B[1], umid, u1, depthb,
247			parameters );
248		✗	if( intersect_BB( A[1], B[1] ) )
249		✗	recursively_intersect( A[1], tmid, t1, deptha,
250			B[1], umid, u1, depthb,
251			parameters );
252		✗	}
253			else
254			{
255		✗	if( intersect_BB( A[0], b ) )
256		✗	recursively_intersect( A[0], t0, tmid, deptha,
257			b, u0, u1, depthb,
258			parameters );
259		✗	if( intersect_BB( A[1], b ) )
260		✗	recursively_intersect( A[1], tmid, t1, deptha,
261			b, u0, u1, depthb,
262			parameters );
263			}
264		✗	}
265			else
266		✗	if( depthb > 0 )
267			{
268		✗	OldBezier B[2];
269		✗	b.split(0.5, B[0], B[1]);
270		✗	double umid = (u0 + u1)*0.5;
271		✗	depthb--;
272		✗	if( intersect_BB( a, B[0] ) )
273		✗	recursively_intersect( a, t0, t1, deptha,
274			B[0], u0, umid, depthb,
275			parameters );
276		✗	if( intersect_BB( a, B[1] ) )
277		✗	recursively_intersect( a, t0, t1, deptha,
278			B[0], umid, u1, depthb,
279			parameters );
280		✗	}
281			else // Both segments are fully subdivided; now do line segments
282			{
283		✗	double xlk = a.p.back()[X] - a.p[0][X];
284		✗	double ylk = a.p.back()[Y] - a.p[0][Y];
285		✗	double xnm = b.p.back()[X] - b.p[0][X];
286		✗	double ynm = b.p.back()[Y] - b.p[0][Y];
287		✗	double xmk = b.p[0][X] - a.p[0][X];
288		✗	double ymk = b.p[0][Y] - a.p[0][Y];
289		✗	double det = xnm * ylk - ynm * xlk;
290		✗	if( 1.0 + det == 1.0 )
291		✗	return;
292			else
293			{
294		✗	double detinv = 1.0 / det;
295		✗	double s = ( xnm * ymk - ynm *xmk ) * detinv;
296		✗	double t = ( xlk * ymk - ylk * xmk ) * detinv;
297		✗	if( ( s < 0.0 ) || ( s > 1.0 ) || ( t < 0.0 ) || ( t > 1.0 ) )
298		✗	return;
299		✗	parameters.emplace_back(t0 + s * ( t1 - t0 ),
300		✗	u0 + t * ( u1 - u0 ));
301			}
302			}
303			}
304
305			inline double log4( double x ) { return log(x)/log(4.); }
306
307			/*
308			* Wang's theorem is used to estimate the level of subdivision required,
309			* but only if the bounding boxes interfere at the top level.
310			* Assuming there is a possible intersection, recursively_intersect is
311			* used to find all the parameters corresponding to intersection points.
312			* these are then sorted and returned in an array.
313			*/
314
315		✗	double Lmax(Point p) {
316		✗	return std::max(fabs(p[X]), fabs(p[Y]));
317			}
318
319
320		✗	unsigned wangs_theorem(OldBezier /*a*/) {
321		✗	return 6; // seems a good approximation!
322
323			/*
324			const double INV_EPS = (1L<<14); // The value of 1.0 / (1L<<14) is enough for most applications
325
326			double la1 = Lmax( ( a.p[2] - a.p[1] ) - (a.p[1] - a.p[0]) );
327			double la2 = Lmax( ( a.p[3] - a.p[2] ) - (a.p[2] - a.p[1]) );
328			double l0 = std::max(la1, la2);
329			unsigned ra;
330			if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 )
331			ra = 0;
332			else
333			ra = (unsigned)ceil( log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) );
334			//std::cout << ra << std::endl;
335			return ra;*/
336			}
337
338			struct rparams
339			{
340			OldBezier &A;
341			OldBezier &B;
342			};
343
344			/*static int
345			intersect_polish_f (const gsl_vector * x, void *params,
346			gsl_vector * f)
347			{
348			const double x0 = gsl_vector_get (x, 0);
349			const double x1 = gsl_vector_get (x, 1);
350
351			Geom::Point dx = ((struct rparams *) params)->A(x0) -
352			((struct rparams *) params)->B(x1);
353
354			gsl_vector_set (f, 0, dx[0]);
355			gsl_vector_set (f, 1, dx[1]);
356
357			return GSL_SUCCESS;
358			}*/
359
360			/*union dbl_64{
361			long long i64;
362			double d64;
363			};*/
364
365			/*static double EpsilonBy(double value, int eps)
366			{
367			dbl_64 s;
368			s.d64 = value;
369			s.i64 += eps;
370			return s.d64;
371			}*/
372
373			/*
374			static void intersect_polish_root (OldBezier &A, double &s,
375			OldBezier &B, double &t) {
376			const gsl_multiroot_fsolver_type *T;
377			gsl_multiroot_fsolver *sol;
378
379			int status;
380			size_t iter = 0;
381
382			const size_t n = 2;
383			struct rparams p = {A, B};
384			gsl_multiroot_function f = {&intersect_polish_f, n, &p};
385
386			double x_init[2] = {s, t};
387			gsl_vector *x = gsl_vector_alloc (n);
388
389			gsl_vector_set (x, 0, x_init[0]);
390			gsl_vector_set (x, 1, x_init[1]);
391
392			T = gsl_multiroot_fsolver_hybrids;
393			sol = gsl_multiroot_fsolver_alloc (T, 2);
394			gsl_multiroot_fsolver_set (sol, &f, x);
395
396			do
397			{
398			iter++;
399			status = gsl_multiroot_fsolver_iterate (sol);
400
401			if (status) // check if solver is stuck
402			break;
403
404			status =
405			gsl_multiroot_test_residual (sol->f, 1e-12);
406			}
407			while (status == GSL_CONTINUE && iter < 1000);
408
409			s = gsl_vector_get (sol->x, 0);
410			t = gsl_vector_get (sol->x, 1);
411
412			gsl_multiroot_fsolver_free (sol);
413			gsl_vector_free (x);
414
415			// This code does a neighbourhood search for minor improvements.
416			double best_v = L1(A(s) - B(t));
417			//std::cout << "------\n" << best_v << std::endl;
418			Point best(s,t);
419			while (true) {
420			Point trial = best;
421			double trial_v = best_v;
422			for(int nsi = -1; nsi < 2; nsi++) {
423			for(int nti = -1; nti < 2; nti++) {
424			Point n(EpsilonBy(best[0], nsi),
425			EpsilonBy(best[1], nti));
426			double c = L1(A(n[0]) - B(n[1]));
427			//std::cout << c << "; ";
428			if (c < trial_v) {
429			trial = n;
430			trial_v = c;
431			}
432			}
433			}
434			if(trial == best) {
435			//std::cout << "\n" << s << " -> " << s - best[0] << std::endl;
436			//std::cout << t << " -> " << t - best[1] << std::endl;
437			//std::cout << best_v << std::endl;
438			s = best[0];
439			t = best[1];
440			return;
441			} else {
442			best = trial;
443			best_v = trial_v;
444			}
445			}
446			}*/
447
448
449		✗	void find_intersections_bezier_recursive( std::vector<std::pair<double, double> > &xs,
450			OldBezier a, OldBezier b)
451			{
452		✗	if( intersect_BB( a, b ) )
453			{
454		✗	recursively_intersect( a, 0., 1., wangs_theorem(a),
455		✗	b, 0., 1., wangs_theorem(b),
456			xs);
457			}
458			/*for(unsigned i = 0; i < xs.size(); i++)
459			intersect_polish_root(a, xs[i].first,
460			b, xs[i].second);*/
461		✗	std::sort(xs.begin(), xs.end());
462		✗	}
463
464
465			};
466
467			/*
468			Local Variables:
469			mode:c++
470			c-file-style:"stroustrup"
471			c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
472			indent-tabs-mode:nil
473			fill-column:99
474			End:
475			*/
476			// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
477