summaryrefslogtreecommitdiff
path: root/dviware/dvisvgm/src/Bezier.cpp
blob: 2d90122531253add65c27713457aee7336b3d5a0 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
/*************************************************************************
** Bezier.cpp                                                           **
**                                                                      **
** This file is part of dvisvgm -- a fast DVI to SVG converter          **
** Copyright (C) 2005-2023 Martin Gieseking <martin.gieseking@uos.de>   **
**                                                                      **
** This program is free software; you can redistribute it and/or        **
** modify it under the terms of the GNU General Public License as       **
** published by the Free Software Foundation; either version 3 of       **
** the License, or (at your option) any later version.                  **
**                                                                      **
** This program is distributed in the hope that it will be useful, but  **
** WITHOUT ANY WARRANTY; without even the implied warranty of           **
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the         **
** GNU General Public License for more details.                         **
**                                                                      **
** You should have received a copy of the GNU General Public License    **
** along with this program; if not, see <http://www.gnu.org/licenses/>. **
*************************************************************************/

#include <algorithm>
#include <utility>
#include "Bezier.hpp"
#include "Matrix.hpp"
#include "utility.hpp"

using namespace std;

QuadBezier::QuadBezier () {
	_points[0] = _points[1] = _points[2]  = DPair(0, 0);
}


QuadBezier::QuadBezier (const DPair &p0, const DPair &p1, const DPair &p2) {
	setPoints(p0, p1, p2);
}


void QuadBezier::setPoints(const DPair &p0, const DPair &p1, const DPair &p2) {
	_points[0] = p0;
	_points[1] = p1;
	_points[2] = p2;
}


/** Returns the value (curve point) at t. */
DPair QuadBezier::valueAt (double t) const {
	const double s = 1-t;
	return _points[0]*s*s + _points[1]*2.0*s*t + _points[2]*t*t;
}


/** Returns the value of the first derivative of the curve at t. */
DPair QuadBezier::derivativeAt (double t) const {
	return _points[0]*(2*t-2) + _points[1]*(2-4*t) + _points[2]*(2*t);
}


/** Returns the arc length of the curve from 0 to t. */
double QuadBezier::arclen (double t) const {
	return math::integral(0, t, 20, [this](double t) -> double {
		DPair deriv = derivativeAt(t);
		return sqrt(deriv.x()*deriv.x() + deriv.y()*deriv.y());
	});
}


////////////////////////////////////////////////////////////////////////////////


CubicBezier::CubicBezier () {
	_points[0] = _points[1] = _points[2] = _points[3] = DPair(0, 0);
}


CubicBezier::CubicBezier (const DPair &p0, const DPair &p1, const DPair &p2, const DPair &p3) {
	setPoints(p0, p1, p2, p3);
}


/** Creates a subcurve of a given Bézier curve.
 *  @param[in] source original curve to be clipped
 *  @param[in] t0 'time' parameter \f$\in[0,1]\f$ of source curve where the subcurve starts
 *  @param[in] t1 'time' parameter \f$\in[0,1]\f$ of source curve where the subcurve ends */
CubicBezier::CubicBezier (const CubicBezier &source, double t0, double t1) {
	if (t0 == t1)
		_points[0] = _points[1] = _points[2] = _points[3] = source.valueAt(t0);
	else {
		if (t0 > t1)
			swap(t0, t1);
		if (t0 == 0)
			source.subdivide(t1, this, nullptr);
		else if (t1 == 1)
			source.subdivide(t0, nullptr, this);
		else {
			CubicBezier subcurve;
			source.subdivide(t0, nullptr, &subcurve);
			subcurve.subdivide((t1-t0)/(1-t0), this, nullptr);
		}
	}
}


/** Creates a cubic Bézier from a quadratic one. */
CubicBezier::CubicBezier (const QuadBezier &qbezier) {
	const DPair &p0 = qbezier.point(0);
	const DPair &p1 = qbezier.point(1);
	const DPair &p2 = qbezier.point(2);
	setPoints(p0, p0+(p1-p0)*2.0/3.0, p2+(p1-p2)*2.0/3.0, p2);
}


void CubicBezier::setPoints(const DPair &p0, const DPair &p1, const DPair &p2, const DPair &p3) {
	_points[0] = p0;
	_points[1] = p1;
	_points[2] = p2;
	_points[3] = p3;
}


void CubicBezier::reverse() {
	swap(_points[0], _points[3]);
	swap(_points[1], _points[2]);
}


DPair CubicBezier::valueAt (double t) const {
	const double s = 1-t;
	return _points[0]*s*s*s + _points[1]*3.0*s*s*t + _points[2]*3.0*s*t*t + _points[3]*t*t*t;
}


/** Returns a value of the Bézier curve's blossom representation. */
DPair CubicBezier::blossomValue (double u, double v, double w) const {
	const double uv = u*v;
	const double uw = u*w;
	const double vw = v*w;
	const double uvw = u*v*w;
	return _points[0]*(1.0-u-v-w+uv+uw+vw-uvw)
		+_points[1]*(u+v+w-2.0*(uv+uw+vw)+3.0*uvw)
		+_points[2]*(uv+uw+vw-3.0*uvw)
		+_points[3]*uvw;
}


/** Splits the curve at t into two sub-curves. */
void CubicBezier::subdivide (double t, CubicBezier *bezier1, CubicBezier *bezier2) const {
	const double s = 1-t;
	DPair p01   = _points[0]*s + _points[1]*t;
	DPair p12   = _points[1]*s + _points[2]*t;
	DPair p23   = _points[2]*s + _points[3]*t;
	DPair p012  = p01*s + p12*t;
	DPair p123  = p12*s + p23*t;
	DPair p0123 = p012*s + p123*t;
	if (bezier1)
		bezier1->setPoints(_points[0], p01, p012, p0123);
	if (bezier2)
		bezier2->setPoints(p0123, p123, p23, _points[3]);
}


/** Approximates the current Bézier curve by a sequence of line segments.
 *  This is done by subdividing the curve several times using De Casteljau's algorithm.
 *  If a sub-curve is almost flat, i.e. \f$\sum\limits_{k=0}^2 |p_{k+1}-p_k| - |p_3-p_0| < \delta\f$,
 *  the curve is not further subdivided.
 *  @param[in] delta threshold where to stop further subdivisions (see description above)
 *  @param[out] p the resulting sequence of points defining the start/end points of the line segments
 *  @param[out] t corresponding curve parameters of the approximated points p: \f$ b(t_i)=p_i \f$
 *  @return number of points in vector p */
int CubicBezier::approximate (double delta, std::vector<DPair> &p, vector<double> *t) const {
	p.push_back(_points[0]);
	if (t)
		t->push_back(0);
	return approximate(delta, 0, 1, p, t);
}


int CubicBezier::approximate (double delta, double t0, double t1, vector<DPair> &p, vector<double> *t) const {
	// compute distance of adjacent control points
	const double l01 = (_points[1]-_points[0]).length();
	const double l12 = (_points[2]-_points[1]).length();
	const double l23 = (_points[3]-_points[2]).length();
	const double l03 = (_points[3]-_points[0]).length();
	if (l01+l12+l23-l03 < delta) { // is curve flat enough?
		p.push_back(_points[3]);    // => store endpoint
		if (t)
			t->push_back(t1);
	}
	else {
		// subdivide curve at b(0.5) and approximate the resulting parts separately
		CubicBezier b1, b2;
		subdivide(0.5, &b1, &b2);
		double tmid = (t0+t1)/2;
		b1.approximate(delta, t0, tmid, p, t);
		b2.approximate(delta, tmid, t1, p, t);
	}
	return static_cast<int>(p.size());
}


/** Returns the signed area of the triangle (p1, p2, p3). */
static inline double signed_area (const DPair &p1, const DPair &p2, const DPair &p3) {
	return ((p2.x()-p1.x())*(p3.y()-p1.y()) - (p3.x()-p1.x())*(p2.y()-p1.y()))/2.0;
}


static inline double dot_prod (const DPair &p1, const DPair &p2) {
	return p1.x()*p2.x() + p1.y()*p2.y();
}


/** Returns true if p3 is located between p1 and p2, i.e. p3 lays almost on the line
 *  between p1 and p2. */
static bool between (const DPair &p1, const DPair &p2, const DPair &p3, double delta) {
	double sqr_dist = dot_prod(p2-p1, p2-p1);
	double factor = sqr_dist == 0.0 ? 1.0 : sqr_dist;
	double area2 = abs(signed_area(p1, p2, p3));
	return area2*area2/factor < delta    // does p3 lay almost on the line through p1 and p2...
		&& min(p1.x(), p2.x()) <= p3.x()  // ...and on or inside the rectangle spanned by p1 and p2?
		&& max(p1.x(), p2.x()) >= p3.x()
		&& min(p1.y(), p2.y()) <= p3.y()
		&& max(p1.y(), p2.y()) >= p3.y();
}


static inline bool near (const DPair &p1, const DPair &p2, double delta) {
	DPair diff = p2-p1;
	return abs(diff.x()) < delta && abs(diff.y()) < delta;
}


/** Tries to reduce the degree of the Bézier curve. This only works if the number of
 *  control points can be reduces without changing the shape of the curve significantly.
 *  @param[in] delta deviation tolerance
 *  @param[in] p control points of the reduced curve
 *  @return degree of the reduced curve */
int CubicBezier::reduceDegree (double delta, vector<DPair> &p) const {
	p.clear();
	if (near(_points[0], _points[1], delta) && near(_points[0], _points[2], delta) && near(_points[0], _points[3], delta))
		p.push_back(_points[0]);
	else if (between(_points[0], _points[3], _points[1], delta) && between(_points[0], _points[3], _points[2], delta)) {
		p.push_back(_points[0]);
		p.push_back(_points[3]);
	}
	else if (near((_points[1]-_points[0])*1.5+_points[0], (_points[2]-_points[3])*1.5+_points[3], delta)) {
		p.push_back(_points[0]);
		p.push_back((_points[1]-_points[0])*1.5 + _points[0]);
		p.push_back(_points[3]);
	}
	else {
		p.resize(4);
		for (int i=0; i < 4; i++)
			p[i] = _points[i];
	}
	return static_cast<int>(p.size()-1);
}


/** Approximates the cubic Bézier curve by a sequence of quadratic ones.
 *  @param[in] precision specifies the precision of the approximation
 *  @return map containing the split parameters t_n together with the qudratic curves */
vector<QuadBezier> CubicBezier::toQuadBeziers (double precision, vector<double> *splitParams) const {
	vector<QuadBezier> qbeziers;
	toQuadBeziers(0, 1, precision, qbeziers, splitParams);
	return qbeziers;
}


/** Returns the "mid-point approximation" of this cubic Bézier. */
QuadBezier CubicBezier::midpointApproximation () const {
	 DPair p0 = (_points[1]*3.0 - _points[0])/2.0;
	 DPair p1 = (_points[2]*3.0 - _points[3])/2.0;
	 return QuadBezier(_points[0], (p0+p1)/2.0, _points[3]);
}


/** Approximates a segment of a cubic Bézier curve by a sequence of quadratic curves.
 *  The quadratic segments are computed by adaptive subdivision of the cubic curve
 *  as described at http://www.caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
 *  @param[in] t0 curve parameter of the segment's start point
 *  @param[in] t1 curve parameter of the segment's end point
 *  @param[in] precision maximum allowed distance between the curve points b(t) and the
 *  	corresponding approximated point c(t) on the quadratic curve
 *  @param[out] qbeziers the resulting quadratic curves sorted in ascending order by their split points
 *  @param[out] startParams the start parameters t[k] of the curves qbeziers[k] relative to *this. */
void CubicBezier::toQuadBeziers (double t0, double t1, double precision, vector<QuadBezier> &qbeziers, vector<double> *startParams) const {
	// If -p0+3p1-3p2+p3 = 0, the degree of the curve is <= 2 and it's not necessary do any approximation.
	// In this case, the control point of the quadratic Bézier curve is (-p0+3p1)/2 = (3p2-p3)/2.
	// Otherwise, the distance d between q1:=(-p0+3p1)/2 and q2:=(3p2-p3)/2 is != 0.
	// Now we compute the quadratic Bézier with start point p0, end point p3, and control point (q1+q2)/2,
	// the "mid-point approximation" (MPA) of b.
	// The maximal distance between the points of the original curve b and the corresponding ones
	// on the MPA is d*sqrt(3)/18. The same computations can be done for the cubic curves we get
	// when subdividing b at a parameter t. The maximal distance of these curves to their MPA
	// is t^3*d*sqrt(3)/18.
	// Based on the formula for this distance we compute the split point tmax for a given precision
	// and check 3 cases:
	// * tmax >= 1:   curve is quadratic, return MPA of b
	// * tmax >= 0.5: split b at t=0.5, return the MPAs of both segments
	// * tmax < 0.5:  split b at tmax and 1-tmax, return the MPAs of the first and third segment,
	//                recurse the algorithm for the middle segment
	DPair q1 = (_points[1]*3.0 - _points[0])/2.0;
	DPair q2 = (_points[2]*3.0 - _points[3])/2.0;
	double dist = (q2-q1).length();
	double tmax3 = 18.0/sqrt(3.0)*precision/dist;  // the cube of tmax
	if (tmax3 >= 1.0) {
		// curve is already quadratic, no subdivision necessary, return MPA
		qbeziers.emplace_back(QuadBezier(_points[0], (q1+q2)/2.0, _points[3]));
		if (startParams)
			startParams->push_back(t0);
	}
	else if (tmax3 >= 0.125) { // tmax >= 0.5
		// split the curve at t=0.5 and compute the MPA for both segments
		CubicBezier cbezier1, cbezier2;
		subdivide(0.5, &cbezier1, &cbezier2);
		qbeziers.emplace_back(cbezier1.midpointApproximation());
		qbeziers.emplace_back(cbezier2.midpointApproximation());
		if (startParams) {
			startParams->push_back(t0);
			startParams->push_back((t0+t1)/2);
		}
	}
	else { // tmax < 0.5
		double tmax = cbrt(tmax3);
		double smax = 1.0-tmax;
		double dt = t1-t0;
		// first segment can be approximated by its MPA
		qbeziers.emplace_back(CubicBezier(*this, 0, tmax).midpointApproximation());
		if (startParams)
			startParams->push_back(t0);
		// recurse for middle segment
		CubicBezier(*this, tmax, smax).toQuadBeziers(t0+tmax*dt, t0+smax*dt, precision, qbeziers, startParams);
		// third segment can be approximated by its MPA
		qbeziers.emplace_back(CubicBezier(*this, smax, 1).midpointApproximation());
		if (startParams)
			startParams->push_back(smax);
	}
}


/** Try to solve the quadratic equation ax^2 + bx + c = 0. */
static bool solve_quadratic_equation (double a, double b, double c, double &x1, double &x2) {
	if (a == 0) {
		if (b == 0)
			return false;
		x1 = x2 = -c/b;
	}
	else {
		double discr = b*b - 4*a*c;
		if (discr < 0)
			return false;
		double p = -b/a/2;
		double r = sqrt(discr)/a/2;
		x1 = p+r;
		x2 = p-r;
	}
	return true;
}


/** Returns a tight bounding box parallel to the x- and y-axis. */
BoundingBox CubicBezier::getBBox () const {
	BoundingBox bbox;
	// coefficients of the derivative
	DPair pa = _points[3] - _points[2]*3.0 + _points[1]*3.0 - _points[0];
	DPair pb = (_points[2]-_points[1]*2.0+_points[0])*2.0;
	DPair pc = _points[1]-_points[0];

	// compute extrema for t > 0 and t < 1
	double t1, t2;
	if (solve_quadratic_equation(pa.x(), pb.x(), pc.x(), t1, t2)) {
		if (t1 > 0.001 && t1 < 0.999)
			bbox.embed(valueAt(t1));
		if (t1 != t2 && t2 > 0.001 && t2 < 0.999)
			bbox.embed(valueAt(t2));
	}
	if (solve_quadratic_equation(pa.y(), pb.y(), pc.y(), t1, t2)) {
		if (t1 > 0.001 && t1 < 0.999)
			bbox.embed(valueAt(t1));
		if (t1 != t2 && t2 > 0.001 && t2 < 0.999)
			bbox.embed(valueAt(t2));
	}
	bbox.embed(_points[0]);
	bbox.embed(_points[3]);
	return bbox;
}


CubicBezier& CubicBezier::transform (const Matrix &matrix) {
	for (auto &point : _points)
		point = matrix * point;
	return *this;
}