1
/* cairo - a vector graphics library with display and print output
2
 *
3
 * Copyright © 2002 University of Southern California
4
 *
5
 * This library is free software; you can redistribute it and/or
6
 * modify it either under the terms of the GNU Lesser General Public
7
 * License version 2.1 as published by the Free Software Foundation
8
 * (the "LGPL") or, at your option, under the terms of the Mozilla
9
 * Public License Version 1.1 (the "MPL"). If you do not alter this
10
 * notice, a recipient may use your version of this file under either
11
 * the MPL or the LGPL.
12
 *
13
 * You should have received a copy of the LGPL along with this library
14
 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
15
 * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA
16
 * You should have received a copy of the MPL along with this library
17
 * in the file COPYING-MPL-1.1
18
 *
19
 * The contents of this file are subject to the Mozilla Public License
20
 * Version 1.1 (the "License"); you may not use this file except in
21
 * compliance with the License. You may obtain a copy of the License at
22
 * http://www.mozilla.org/MPL/
23
 *
24
 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
25
 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
26
 * the specific language governing rights and limitations.
27
 *
28
 * The Original Code is the cairo graphics library.
29
 *
30
 * The Initial Developer of the Original Code is University of Southern
31
 * California.
32
 *
33
 * Contributor(s):
34
 *	Carl D. Worth <cworth@cworth.org>
35
 */
36

            
37
#include "cairoint.h"
38

            
39
#include "cairo-box-inline.h"
40
#include "cairo-slope-private.h"
41

            
42
cairo_bool_t
43
86046
_cairo_spline_intersects (const cairo_point_t *a,
44
			  const cairo_point_t *b,
45
			  const cairo_point_t *c,
46
			  const cairo_point_t *d,
47
			  const cairo_box_t *box)
48
{
49
    cairo_box_t bounds;
50

            
51
154881
    if (_cairo_box_contains_point (box, a) ||
52
137256
	_cairo_box_contains_point (box, b) ||
53
136242
	_cairo_box_contains_point (box, c) ||
54
67821
	_cairo_box_contains_point (box, d))
55
    {
56
18804
	return TRUE;
57
    }
58

            
59
67242
    bounds.p2 = bounds.p1 = *a;
60
67242
    _cairo_box_add_point (&bounds, b);
61
67242
    _cairo_box_add_point (&bounds, c);
62
67242
    _cairo_box_add_point (&bounds, d);
63

            
64
67242
    if (bounds.p2.x <= box->p1.x || bounds.p1.x >= box->p2.x ||
65
34284
	bounds.p2.y <= box->p1.y || bounds.p1.y >= box->p2.y)
66
    {
67
50472
	return FALSE;
68
    }
69

            
70
#if 0 /* worth refining? */
71
    bounds.p2 = bounds.p1 = *a;
72
    _cairo_box_add_curve_to (&bounds, b, c, d);
73
    if (bounds.p2.x <= box->p1.x || bounds.p1.x >= box->p2.x ||
74
	bounds.p2.y <= box->p1.y || bounds.p1.y >= box->p2.y)
75
    {
76
	return FALSE;
77
    }
78
#endif
79

            
80
16770
    return TRUE;
81
}
82

            
83
cairo_bool_t
84
56462
_cairo_spline_init (cairo_spline_t *spline,
85
		    cairo_spline_add_point_func_t add_point_func,
86
		    void *closure,
87
		    const cairo_point_t *a, const cairo_point_t *b,
88
		    const cairo_point_t *c, const cairo_point_t *d)
89
{
90
    /* If both tangents are zero, this is just a straight line */
91
56462
    if (a->x == b->x && a->y == b->y && c->x == d->x && c->y == d->y)
92
	return FALSE;
93

            
94
56462
    spline->add_point_func = add_point_func;
95
56462
    spline->closure = closure;
96

            
97
56462
    spline->knots.a = *a;
98
56462
    spline->knots.b = *b;
99
56462
    spline->knots.c = *c;
100
56462
    spline->knots.d = *d;
101

            
102
56462
    if (a->x != b->x || a->y != b->y)
103
52547
	_cairo_slope_init (&spline->initial_slope, &spline->knots.a, &spline->knots.b);
104
3915
    else if (a->x != c->x || a->y != c->y)
105
3813
	_cairo_slope_init (&spline->initial_slope, &spline->knots.a, &spline->knots.c);
106
102
    else if (a->x != d->x || a->y != d->y)
107
102
	_cairo_slope_init (&spline->initial_slope, &spline->knots.a, &spline->knots.d);
108
    else
109
	return FALSE;
110

            
111
56462
    if (c->x != d->x || c->y != d->y)
112
54404
	_cairo_slope_init (&spline->final_slope, &spline->knots.c, &spline->knots.d);
113
2058
    else if (b->x != d->x || b->y != d->y)
114
1962
	_cairo_slope_init (&spline->final_slope, &spline->knots.b, &spline->knots.d);
115
    else
116
96
	return FALSE; /* just treat this as a straight-line from a -> d */
117

            
118
    /* XXX if the initial, final and vector are all equal, this is just a line */
119

            
120
56366
    return TRUE;
121
}
122

            
123
static cairo_status_t
124
769149
_cairo_spline_add_point (cairo_spline_t *spline,
125
			 const cairo_point_t *point,
126
			 const cairo_point_t *knot)
127
{
128
    cairo_point_t *prev;
129
    cairo_slope_t slope;
130

            
131
769149
    prev = &spline->last_point;
132
769149
    if (prev->x == point->x && prev->y == point->y)
133
56366
	return CAIRO_STATUS_SUCCESS;
134

            
135
712783
    _cairo_slope_init (&slope, point, knot);
136

            
137
712783
    spline->last_point = *point;
138
712783
    return spline->add_point_func (spline->closure, point, &slope);
139
}
140

            
141
static void
142
4276698
_lerp_half (const cairo_point_t *a, const cairo_point_t *b, cairo_point_t *result)
143
{
144
4276698
    result->x = a->x + ((b->x - a->x) >> 1);
145
4276698
    result->y = a->y + ((b->y - a->y) >> 1);
146
4276698
}
147

            
148
static void
149
712783
_de_casteljau (cairo_spline_knots_t *s1, cairo_spline_knots_t *s2)
150
{
151
    cairo_point_t ab, bc, cd;
152
    cairo_point_t abbc, bccd;
153
    cairo_point_t final;
154

            
155
712783
    _lerp_half (&s1->a, &s1->b, &ab);
156
712783
    _lerp_half (&s1->b, &s1->c, &bc);
157
712783
    _lerp_half (&s1->c, &s1->d, &cd);
158
712783
    _lerp_half (&ab, &bc, &abbc);
159
712783
    _lerp_half (&bc, &cd, &bccd);
160
712783
    _lerp_half (&abbc, &bccd, &final);
161

            
162
712783
    s2->a = final;
163
712783
    s2->b = bccd;
164
712783
    s2->c = cd;
165
712783
    s2->d = s1->d;
166

            
167
712783
    s1->b = ab;
168
712783
    s1->c = abbc;
169
712783
    s1->d = final;
170
712783
}
171

            
172
/* Return an upper bound on the error (squared) that could result from
173
 * approximating a spline as a line segment connecting the two endpoints. */
174
static double
175
1481932
_cairo_spline_error_squared (const cairo_spline_knots_t *knots)
176
{
177
    double bdx, bdy, berr;
178
    double cdx, cdy, cerr;
179

            
180
    /* We are going to compute the distance (squared) between each of the b
181
     * and c control points and the segment a-b. The maximum of these two
182
     * distances will be our approximation error. */
183

            
184
1481932
    bdx = _cairo_fixed_to_double (knots->b.x - knots->a.x);
185
1481932
    bdy = _cairo_fixed_to_double (knots->b.y - knots->a.y);
186

            
187
1481932
    cdx = _cairo_fixed_to_double (knots->c.x - knots->a.x);
188
1481932
    cdy = _cairo_fixed_to_double (knots->c.y - knots->a.y);
189

            
190
1481932
    if (knots->a.x != knots->d.x || knots->a.y != knots->d.y) {
191
	/* Intersection point (px):
192
	 *     px = p1 + u(p2 - p1)
193
	 *     (p - px) ∙ (p2 - p1) = 0
194
	 * Thus:
195
	 *     u = ((p - p1) ∙ (p2 - p1)) / ∥p2 - p1∥²;
196
	 */
197

            
198
	double dx, dy, u, v;
199

            
200
1481896
	dx = _cairo_fixed_to_double (knots->d.x - knots->a.x);
201
1481896
	dy = _cairo_fixed_to_double (knots->d.y - knots->a.y);
202
1481896
	 v = dx * dx + dy * dy;
203

            
204
1481896
	u = bdx * dx + bdy * dy;
205
1481896
	if (u <= 0) {
206
	    /* bdx -= 0;
207
	     * bdy -= 0;
208
	     */
209
1465809
	} else if (u >= v) {
210
1218
	    bdx -= dx;
211
1218
	    bdy -= dy;
212
	} else {
213
1464591
	    bdx -= u/v * dx;
214
1464591
	    bdy -= u/v * dy;
215
	}
216

            
217
1481896
	u = cdx * dx + cdy * dy;
218
1481896
	if (u <= 0) {
219
	    /* cdx -= 0;
220
	     * cdy -= 0;
221
	     */
222
1480321
	} else if (u >= v) {
223
10951
	    cdx -= dx;
224
10951
	    cdy -= dy;
225
	} else {
226
1469370
	    cdx -= u/v * dx;
227
1469370
	    cdy -= u/v * dy;
228
	}
229
    }
230

            
231
1481932
    berr = bdx * bdx + bdy * bdy;
232
1481932
    cerr = cdx * cdx + cdy * cdy;
233
1481932
    if (berr > cerr)
234
700999
	return berr;
235
    else
236
780933
	return cerr;
237
}
238

            
239
static cairo_status_t
240
1481932
_cairo_spline_decompose_into (cairo_spline_knots_t *s1,
241
			      double tolerance_squared,
242
			      cairo_spline_t *result)
243
{
244
    cairo_spline_knots_t s2;
245
    cairo_status_t status;
246

            
247
1481932
    if (_cairo_spline_error_squared (s1) < tolerance_squared)
248
769149
	return _cairo_spline_add_point (result, &s1->a, &s1->b);
249

            
250
712783
    _de_casteljau (s1, &s2);
251

            
252
712783
    status = _cairo_spline_decompose_into (s1, tolerance_squared, result);
253
712783
    if (unlikely (status))
254
	return status;
255

            
256
712783
    return _cairo_spline_decompose_into (&s2, tolerance_squared, result);
257
}
258

            
259
cairo_status_t
260
56366
_cairo_spline_decompose (cairo_spline_t *spline, double tolerance)
261
{
262
    cairo_spline_knots_t s1;
263
    cairo_status_t status;
264

            
265
56366
    s1 = spline->knots;
266
56366
    spline->last_point = s1.a;
267
56366
    status = _cairo_spline_decompose_into (&s1, tolerance * tolerance, spline);
268
56366
    if (unlikely (status))
269
	return status;
270

            
271
56366
    return spline->add_point_func (spline->closure,
272
56366
				   &spline->knots.d, &spline->final_slope);
273
}
274

            
275
/* Note: this function is only good for computing bounds in device space. */
276
cairo_status_t
277
303960
_cairo_spline_bound (cairo_spline_add_point_func_t add_point_func,
278
		     void *closure,
279
		     const cairo_point_t *p0, const cairo_point_t *p1,
280
		     const cairo_point_t *p2, const cairo_point_t *p3)
281
{
282
    double x0, x1, x2, x3;
283
    double y0, y1, y2, y3;
284
    double a, b, c;
285
    double t[4];
286
303960
    int t_num = 0, i;
287
    cairo_status_t status;
288

            
289
303960
    x0 = _cairo_fixed_to_double (p0->x);
290
303960
    y0 = _cairo_fixed_to_double (p0->y);
291
303960
    x1 = _cairo_fixed_to_double (p1->x);
292
303960
    y1 = _cairo_fixed_to_double (p1->y);
293
303960
    x2 = _cairo_fixed_to_double (p2->x);
294
303960
    y2 = _cairo_fixed_to_double (p2->y);
295
303960
    x3 = _cairo_fixed_to_double (p3->x);
296
303960
    y3 = _cairo_fixed_to_double (p3->y);
297

            
298
    /* The spline can be written as a polynomial of the four points:
299
     *
300
     *   (1-t)³p0 + 3t(1-t)²p1 + 3t²(1-t)p2 + t³p3
301
     *
302
     * for 0≤t≤1.  Now, the X and Y components of the spline follow the
303
     * same polynomial but with x and y replaced for p.  To find the
304
     * bounds of the spline, we just need to find the X and Y bounds.
305
     * To find the bound, we take the derivative and equal it to zero,
306
     * and solve to find the t's that give the extreme points.
307
     *
308
     * Here is the derivative of the curve, sorted on t:
309
     *
310
     *   3t²(-p0+3p1-3p2+p3) + 2t(3p0-6p1+3p2) -3p0+3p1
311
     *
312
     * Let:
313
     *
314
     *   a = -p0+3p1-3p2+p3
315
     *   b =  p0-2p1+p2
316
     *   c = -p0+p1
317
     *
318
     * Gives:
319
     *
320
     *   a.t² + 2b.t + c = 0
321
     *
322
     * With:
323
     *
324
     *   delta = b*b - a*c
325
     *
326
     * the extreme points are at -c/2b if a is zero, at (-b±√delta)/a if
327
     * delta is positive, and at -b/a if delta is zero.
328
     */
329

            
330
#define ADD(t0) \
331
    { \
332
	double _t0 = (t0); \
333
	if (0 < _t0 && _t0 < 1) \
334
	    t[t_num++] = _t0; \
335
    }
336

            
337
#define FIND_EXTREMES(a,b,c) \
338
    { \
339
	if (a == 0) { \
340
	    if (b != 0) \
341
		ADD (-c / (2*b)); \
342
	} else { \
343
	    double b2 = b * b; \
344
	    double delta = b2 - a * c; \
345
	    if (delta > 0) { \
346
		cairo_bool_t feasible; \
347
		double _2ab = 2 * a * b; \
348
		/* We are only interested in solutions t that satisfy 0<t<1 \
349
		 * here.  We do some checks to avoid sqrt if the solutions \
350
		 * are not in that range.  The checks can be derived from: \
351
		 * \
352
		 *   0 < (-b±√delta)/a < 1 \
353
		 */ \
354
		if (_2ab >= 0) \
355
		    feasible = delta > b2 && delta < a*a + b2 + _2ab; \
356
		else if (-b / a >= 1) \
357
		    feasible = delta < b2 && delta > a*a + b2 + _2ab; \
358
		else \
359
		    feasible = delta < b2 || delta < a*a + b2 + _2ab; \
360
	        \
361
		if (unlikely (feasible)) { \
362
		    double sqrt_delta = sqrt (delta); \
363
		    ADD ((-b - sqrt_delta) / a); \
364
		    ADD ((-b + sqrt_delta) / a); \
365
		} \
366
	    } else if (delta == 0) { \
367
		ADD (-b / a); \
368
	    } \
369
	} \
370
    }
371

            
372
    /* Find X extremes */
373
303960
    a = -x0 + 3*x1 - 3*x2 + x3;
374
303960
    b =  x0 - 2*x1 + x2;
375
303960
    c = -x0 + x1;
376
303960
    FIND_EXTREMES (a, b, c);
377

            
378
    /* Find Y extremes */
379
303960
    a = -y0 + 3*y1 - 3*y2 + y3;
380
303960
    b =  y0 - 2*y1 + y2;
381
303960
    c = -y0 + y1;
382
303960
    FIND_EXTREMES (a, b, c);
383

            
384
303960
    status = add_point_func (closure, p0, NULL);
385
303960
    if (unlikely (status))
386
	return status;
387

            
388
610440
    for (i = 0; i < t_num; i++) {
389
	cairo_point_t p;
390
	double x, y;
391
        double t_1_0, t_0_1;
392
        double t_2_0, t_0_2;
393
        double t_3_0, t_2_1_3, t_1_2_3, t_0_3;
394

            
395
306480
        t_1_0 = t[i];          /*      t  */
396
306480
        t_0_1 = 1 - t_1_0;     /* (1 - t) */
397

            
398
306480
        t_2_0 = t_1_0 * t_1_0; /*      t  *      t  */
399
306480
        t_0_2 = t_0_1 * t_0_1; /* (1 - t) * (1 - t) */
400

            
401
306480
        t_3_0   = t_2_0 * t_1_0;     /*      t  *      t  *      t      */
402
306480
        t_2_1_3 = t_2_0 * t_0_1 * 3; /*      t  *      t  * (1 - t) * 3 */
403
306480
        t_1_2_3 = t_1_0 * t_0_2 * 3; /*      t  * (1 - t) * (1 - t) * 3 */
404
306480
        t_0_3   = t_0_1 * t_0_2;     /* (1 - t) * (1 - t) * (1 - t)     */
405

            
406
        /* Bezier polynomial */
407
306480
        x = x0 * t_0_3
408
306480
          + x1 * t_1_2_3
409
306480
          + x2 * t_2_1_3
410
306480
          + x3 * t_3_0;
411
306480
        y = y0 * t_0_3
412
306480
          + y1 * t_1_2_3
413
306480
          + y2 * t_2_1_3
414
306480
          + y3 * t_3_0;
415

            
416
306480
	p.x = _cairo_fixed_from_double (x);
417
306480
	p.y = _cairo_fixed_from_double (y);
418
306480
	status = add_point_func (closure, &p, NULL);
419
306480
	if (unlikely (status))
420
	    return status;
421
    }
422

            
423
303960
    return add_point_func (closure, p3, NULL);
424
}