1// [Blend2D]
2// 2D Vector Graphics Powered by a JIT Compiler.
3//
4// [License]
5// Zlib - See LICENSE.md file in the package.
6
7#ifndef BLEND2D_BLGEOMETRY_P_H
8#define BLEND2D_BLGEOMETRY_P_H
9
10#include "./blgeometry.h"
11#include "./blmath_p.h"
12#include "./blsupport_p.h"
13
14//! \cond INTERNAL
15//! \addtogroup blend2d_internal
16//! \{
17
18// ============================================================================
19// [Math Extensions]
20// ============================================================================
21
22static BL_INLINE bool blIsNaN(const BLPoint& p) noexcept { return blIsNaN(p.x + p.y); }
23static BL_INLINE bool blIsFinite(const BLPoint& p) noexcept { return blIsFinite(p.x, p.y); }
24static BL_INLINE bool blIsFinite(const BLBox& b) noexcept { return blIsFinite(b.x0, b.y0, b.x1, b.y1); }
25static BL_INLINE bool blIsFinite(const BLRect& r) noexcept { return blIsFinite(r.x, r.y, r.w, r.h); }
26
27static BL_INLINE BLPoint blCopySign(const BLPoint& a, const BLPoint& b) noexcept {
28 return BLPoint(blCopySign(a.x, b.x), blCopySign(a.y, b.y));
29}
30
31static BL_INLINE BLPoint blSqrt(const BLPoint& p) noexcept {
32 return BLPoint(blSqrt(p.x), blSqrt(p.y));
33}
34
35static BL_INLINE size_t blSimplifiedQuadRoots(BLPoint dst[2], const BLPoint& a, const BLPoint& b, const BLPoint& c) noexcept {
36 BLPoint d = blMax(b * b - 4.0 * a * c, 0.0);
37 BLPoint s = blSqrt(d);
38 BLPoint q = -0.5 * (b + blCopySign(s, b));
39
40 dst[0] = q / a;
41 dst[1] = c / q;
42
43 return 2;
44}
45
46static BL_INLINE bool blIsZero(const BLPoint& p) noexcept {
47 return (p.x == 0) & (p.y == 0);
48}
49
50// ============================================================================
51// [IsValid]
52// ============================================================================
53
54static BL_INLINE bool blIsValid(const BLSizeI& size) noexcept { return (size.w > 0) & (size.h > 0); }
55static BL_INLINE bool blIsValid(const BLSize& size) noexcept { return (size.w > 0) & (size.h > 0); }
56
57static BL_INLINE bool blIsValid(const BLBoxI& box) noexcept { return (box.x0 < box.x1) & (box.y0 < box.y1); }
58static BL_INLINE bool blIsValid(const BLBox& box) noexcept { return (box.x0 < box.x1) & (box.y0 < box.y1); }
59
60static BL_INLINE bool blIsValid(const BLRectI& rect) noexcept {
61 BLOverflowFlag of = 0;
62 int x1 = blAddOverflow(rect.x, rect.w, &of);
63 int y1 = blAddOverflow(rect.y, rect.h, &of);
64 return (rect.x < x1) & (rect.y < y1) & (!of);
65}
66
67static BL_INLINE bool blIsValid(const BLRect& rect) noexcept {
68 double x1 = rect.x + rect.w;
69 double y1 = rect.y + rect.h;
70 return (rect.x < x1) & (rect.y < y1);
71}
72
73// ============================================================================
74// [Box/Rect Manipulation]
75// ============================================================================
76
77static BL_INLINE bool blIntersectBoxes(BLBoxI& dst, const BLBoxI& a, const BLBoxI& b) noexcept {
78 dst.reset(blMax(a.x0, b.x0), blMax(a.y0, b.y0),
79 blMin(a.x1, b.x1), blMin(a.y1, b.y1));
80 return (dst.x0 < dst.x1) & (dst.y0 < dst.y1);
81}
82
83static BL_INLINE bool blIntersectBoxes(BLBox& dst, const BLBox& a, const BLBox& b) noexcept {
84 dst.reset(blMax(a.x0, b.x0), blMax(a.y0, b.y0),
85 blMin(a.x1, b.x1), blMin(a.y1, b.y1));
86 return (dst.x0 < dst.x1) & (dst.y0 < dst.y1);
87}
88
89static BL_INLINE void blBoundBoxes(BLBox& box, const BLPoint& p) noexcept {
90 box.reset(blMin(box.x0, p.x), blMin(box.y0, p.y),
91 blMax(box.x1, p.x), blMax(box.y1, p.y));
92}
93
94static BL_INLINE void blBoundBoxes(BLBox& box, const BLBox& other) noexcept {
95 box.reset(blMin(box.x0, other.x0), blMin(box.y0, other.y0),
96 blMax(box.x1, other.x1), blMax(box.y1, other.y1));
97}
98
99static BL_INLINE void blBoundBoxes(BLBoxI& box, const BLBoxI& other) noexcept {
100 box.reset(blMin(box.x0, other.x0), blMin(box.y0, other.y0),
101 blMax(box.x1, other.x1), blMax(box.y1, other.y1));
102}
103
104static BL_INLINE bool blSubsumes(const BLBoxI& a, const BLBoxI& b) noexcept { return (a.x0 <= b.x0) & (a.y0 <= b.y0) & (a.x1 >= b.x1) & (a.y1 >= b.y1); }
105static BL_INLINE bool blSubsumes(const BLBox& a, const BLBox& b) noexcept { return (a.x0 <= b.x0) & (a.y0 <= b.y0) & (a.x1 >= b.x1) & (a.y1 >= b.y1); }
106
107static BL_INLINE bool blOverlaps(const BLBoxI& a, const BLBoxI& b) noexcept { return (a.x1 > b.x0) & (a.y1 > b.y0) & (a.x0 < b.x1) & (a.y0 < b.y1); }
108static BL_INLINE bool blOverlaps(const BLBox& a, const BLBox& b) noexcept { return (a.x1 > b.x0) & (a.y1 > b.y0) & (a.x0 < b.x1) & (a.y0 < b.y1); }
109
110// ============================================================================
111// [Point / Vector]
112// ============================================================================
113
114static BL_INLINE double blLengthSq(const BLPoint& v) noexcept { return v.x * v.x + v.y * v.y; }
115static BL_INLINE double blLengthSq(const BLPoint& a, const BLPoint& b) noexcept { return blLengthSq(b - a); }
116
117static BL_INLINE double blLength(const BLPoint& v) noexcept { return blSqrt(blLengthSq(v)); }
118static BL_INLINE double blLength(const BLPoint& a, const BLPoint& b) noexcept { return blSqrt(blLengthSq(a, b)); }
119
120static BL_INLINE BLPoint blNormal(const BLPoint& v) noexcept { return BLPoint(-v.y, v.x); }
121static BL_INLINE BLPoint blUnitVector(const BLPoint& v) noexcept { return v / blLength(v); }
122
123static BL_INLINE double blDotProduct(const BLPoint& a, const BLPoint& b) noexcept { return a.x * b.x + a.y * b.y; }
124static BL_INLINE double blCrossProduct(const BLPoint& a, const BLPoint& b) noexcept { return a.x * b.y - a.y * b.x; }
125
126// ============================================================================
127// [Line]
128// ============================================================================
129
130static BL_INLINE BLPoint blGetLineVectorIntersection(const BLPoint& p0, const BLPoint& v0, const BLPoint& p1, const BLPoint& v1) noexcept {
131 return p0 + blCrossProduct(p1 - p0, v1) / blCrossProduct(v0, v1) * v0;
132}
133
134// ============================================================================
135// [Quad]
136// ============================================================================
137
138// Quad Coefficients
139// -----------------
140//
141// A = p0 + 2*p1 + p2
142// B = -2*p0 + 2*p1
143// C = p0
144//
145// Quad Evaluation at `t`
146// ----------------------
147//
148// V = At^2 + Bt + C => t(At + B) + C
149
150static BL_INLINE void blSplitQuad(const BLPoint p[3], BLPoint aOut[3], BLPoint bOut[3]) noexcept {
151 BLPoint p01(blLerp(p[0], p[1]));
152 BLPoint p12(blLerp(p[1], p[2]));
153
154 aOut[0] = p[0];
155 aOut[1] = p01;
156 bOut[1] = p12;
157 bOut[2] = p[2];
158 aOut[2] = blLerp(p01, p12);
159 bOut[0] = aOut[2];
160}
161
162static BL_INLINE void blSplitQuad(const BLPoint p[3], BLPoint aOut[3], BLPoint bOut[3], double t) noexcept {
163 BLPoint p01(blLerp(p[0], p[1], t));
164 BLPoint p12(blLerp(p[1], p[2], t));
165
166 aOut[0] = p[0];
167 aOut[1] = p01;
168 bOut[1] = p12;
169 bOut[2] = p[2];
170 aOut[2] = blLerp(p01, p12, t);
171 bOut[0] = aOut[2];
172}
173
174static BL_INLINE void blSplitQuadBefore(const BLPoint p[3], BLPoint out[3], double t) noexcept {
175 BLPoint p01(blLerp(p[0], p[1], t));
176 BLPoint p12(blLerp(p[1], p[2], t));
177
178 out[0] = p[0];
179 out[1] = p01;
180 out[2] = blLerp(p01, p12, t);
181}
182
183static BL_INLINE void blSplitQuadAfter(const BLPoint p[3], BLPoint out[3], double t) noexcept {
184 BLPoint p01(blLerp(p[0], p[1], t));
185 BLPoint p12(blLerp(p[1], p[2], t));
186
187 out[0] = blLerp(p01, p12, t);
188 out[1] = p12;
189 out[2] = p[2];
190}
191
192static BL_INLINE void blSplitQuadBetween(const BLPoint p[3], BLPoint out[3], double t0, double t1) noexcept {
193 BLPoint t0p01 = blLerp(p[0], p[1], t0);
194 BLPoint t0p12 = blLerp(p[1], p[2], t0);
195
196 BLPoint t1p01 = blLerp(p[0], p[1], t1);
197 BLPoint t1p12 = blLerp(p[1], p[2], t1);
198
199 out[0] = blLerp(t0p01, t0p12, t0);
200 out[1] = blLerp(t0p01, t0p12, t1);
201 out[2] = blLerp(t1p01, t1p12, t1);
202}
203
204static BL_INLINE void blGetQuadCoefficients(const BLPoint p[3], BLPoint& a, BLPoint& b, BLPoint& c) noexcept {
205 BLPoint v1 = p[1] - p[0];
206 BLPoint v2 = p[2] - p[1];
207
208 a = v2 - v1;
209 b = v1 + v1;
210 c = p[0];
211}
212
213static BL_INLINE void blGetQuadDerivativeCoefficients(const BLPoint p[3], BLPoint& a, BLPoint& b) noexcept {
214 BLPoint v1 = p[1] - p[0];
215 BLPoint v2 = p[2] - p[1];
216
217 a = 2.0 * v2 - 2.0 * v1;
218 b = 2.0 * v1;
219}
220
221static BL_INLINE BLPoint blGetQuadValueAt(const BLPoint p[3], double t) noexcept {
222 BLPoint a, b, c;
223 blGetQuadCoefficients(p, a, b, c);
224 return (a * t + b) * t + c;
225}
226
227static BL_INLINE BLPoint blGetQuadValueAt(const BLPoint p[3], const BLPoint& t) noexcept {
228 BLPoint a, b, c;
229 blGetQuadCoefficients(p, a, b, c);
230 return (a * t + b) * t + c;
231}
232
233static BL_INLINE BLPoint blGetPreciseQuadValueAt(const BLPoint p[3], double t) noexcept {
234 return blLerp(blLerp(p[0], p[1], t), blLerp(p[1], p[2], t), t);
235}
236
237static BL_INLINE BLPoint blGetPreciseQuadValueAt(const BLPoint p[3], const BLPoint& t) noexcept {
238 return blLerp(blLerp(p[0], p[1], t), blLerp(p[1], p[2], t), t);
239}
240
241static BL_INLINE void blGetQuadExtremaPoint(const BLPoint p[3], BLPoint& out) noexcept {
242 BLPoint t = blClamp((p[0] - p[1]) / (p[0] - p[1] * 2.0 + p[2]), 0.0, 1.0);
243 out = blGetPreciseQuadValueAt(p, t);
244}
245
246static BL_INLINE double blGetQuadParameterAtAngle(const BLPoint p[3], double m) noexcept {
247 BLPoint qa, qb;
248 blGetQuadDerivativeCoefficients(p, qa, qb);
249
250 double aob = blDotProduct(qa, qb);
251 double axb = blCrossProduct(qa, qb);
252
253 if (aob == 0.0)
254 return 1.0;
255
256 // m * (bx * bx + by * by) / (|ax * by - ay * bx| - m * (ax * bx + ay * by));
257 return m * blLengthSq(qb) / (blAbs(axb) - m * aob);
258}
259
260static BL_INLINE double blGetQuadCurvatureMetric(const BLPoint p[3]) noexcept {
261 return blCrossProduct(p[2] - p[1], p[1] - p[0]);
262}
263
264static BL_INLINE size_t blGetQuadOffsetCuspTs(const BLPoint bez[3], double d, double tOut[2]) {
265 BLPoint qqa, qqb;
266 blGetQuadDerivativeCoefficients(bez, qqa, qqb);
267
268 double bxa = blCrossProduct(qqb, qqa);
269 double boa = blDotProduct(qqb, qqa);
270
271 if (bxa == 0)
272 return 0;
273
274 double alen2 = blLengthSq(qqa);
275 double blen2 = blLengthSq(qqb);
276
277 double fac = -1.0 / alen2;
278 double sqrt_ = blSqrt(boa * boa - alen2 * (blen2 - blCbrt(d * d * bxa * bxa)));
279
280 double t0 = fac * (boa + sqrt_);
281 double t1 = fac * (boa - sqrt_);
282
283 // We are only interested in (0, 1) range.
284 t0 = blMax(t0, 0.0);
285
286 size_t n = size_t(t0 > 0.0 && t0 < 1.0);
287 tOut[0] = t0;
288 tOut[n] = t1;
289 return n + size_t(t1 > t0 && t1 < 1.0);
290}
291
292//! Coverts quadratic curve to cubic curve.
293//!
294//! \code
295//! cubic[0] = q0
296//! cubic[1] = q0 + 2/3 * (q1 - q0)
297//! cubic[2] = q2 + 2/3 * (q1 - q2)
298//! cubic[3] = q2
299//! \endcode
300static BL_INLINE void blQuadToCubic(const BLPoint p[3], BLPoint cubicOut[4]) noexcept {
301 constexpr double k1Div3 = 1.0 / 3.0;
302 constexpr double k2Div3 = 2.0 / 3.0;
303
304 BLPoint tmp = p[1] * k2Div3;
305 cubicOut[0] = p[0];
306 cubicOut[3] = p[2];
307 cubicOut[1].reset(cubicOut[0] * k1Div3 + tmp);
308 cubicOut[2].reset(cubicOut[2] * k1Div3 + tmp);
309}
310
311class BLQuadCurveTsIter {
312public:
313 const double* ts;
314 const double* tsEnd;
315
316 BLPoint input[3];
317 BLPoint part[3];
318 BLPoint pTmp01;
319 BLPoint pTmp12;
320
321 BL_INLINE BLQuadCurveTsIter() noexcept
322 : ts(nullptr),
323 tsEnd(nullptr) {}
324
325 BL_INLINE BLQuadCurveTsIter(const BLPoint* input_, const double* ts_, size_t count_) noexcept {
326 reset(input_, ts_, count_);
327 }
328
329 BL_INLINE void reset(const BLPoint* input_, const double* ts_, size_t count_) noexcept {
330 // There must be always at least one T.
331 BL_ASSERT(count_ > 0);
332
333 input[0] = input_[0];
334 input[1] = input_[1];
335 input[2] = input_[2];
336 ts = ts_;
337 tsEnd = ts + count_;
338
339 // The first iterated curve is the same as if we split left side at `t`.
340 // This behaves identically to `blSplitQuadBefore()`, however, we cache
341 // `pTmp01` and `pTmp12` for reuse in `next()`.
342 double t = *ts++;
343 pTmp01 = blLerp(input[0], input[1], t);
344 pTmp12 = blLerp(input[1], input[2], t);
345
346 part[0] = input[0];
347 part[1] = pTmp01;
348 part[2] = blLerp(part[1], pTmp12, t);
349 }
350
351 BL_INLINE bool next() noexcept {
352 if (ts >= tsEnd)
353 return false;
354
355 double t = *ts++;
356 part[0] = part[2];
357 part[1] = blLerp(pTmp01, pTmp12, t);
358
359 pTmp01 = blLerp(input[0], input[1], t);
360 pTmp12 = blLerp(input[1], input[2], t);
361 part[2] = blLerp(pTmp01, pTmp12, t);
362 return true;
363 }
364};
365
366// ============================================================================
367// [Cubic]
368// ============================================================================
369
370// Cubic Coefficients
371// ------------------
372//
373// A = -p0 + 3*p1 - 3*p2 + p3 => 3*(p1 - p2) + p3 - p0
374// B = 3*p0 - 6*p1 + 3*p2 => 3*(p0 - 2*p2 + p3)
375// C = -3*p0 + 3*p1 => 3*(p1 - p0)
376// D = p0 => p0
377//
378// Cubic Evaluation at `t`
379// -----------------------
380//
381// V = At^3 + Bt^2 + Ct + D => t(t(At + B) + C) + D
382
383static BL_INLINE void blSplitCubic(const BLPoint p[4], BLPoint a[4], BLPoint b[4]) noexcept {
384 BLPoint p01(blLerp(p[0], p[1]));
385 BLPoint p12(blLerp(p[1], p[2]));
386 BLPoint p23(blLerp(p[2], p[3]));
387
388 a[0] = p[0];
389 a[1] = p01;
390 b[2] = p23;
391 b[3] = p[3];
392
393 a[2] = blLerp(p01, p12);
394 b[1] = blLerp(p12, p23);
395 a[3] = blLerp(a[2], b[1]);
396 b[0] = a[3];
397}
398
399static BL_INLINE void blSplitCubic(const BLPoint p[4], BLPoint before[4], BLPoint after[4], double t) noexcept {
400 BLPoint p01(blLerp(p[0], p[1], t));
401 BLPoint p12(blLerp(p[1], p[2], t));
402 BLPoint p23(blLerp(p[2], p[3], t));
403
404 before[0] = p[0];
405 before[1] = p01;
406 after[2] = p23;
407 after[3] = p[3];
408
409 before[2] = blLerp(p01, p12, t);
410 after[1] = blLerp(p12, p23, t);
411 before[3] = blLerp(before[2], after[1], t);
412 after[0] = before[3];
413}
414
415static BL_INLINE void blSplitCubicBefore(const BLPoint p[4], BLPoint a[4], double t) noexcept {
416 BLPoint p01(blLerp(p[0], p[1], t));
417 BLPoint p12(blLerp(p[1], p[2], t));
418 BLPoint p23(blLerp(p[2], p[3], t));
419
420 a[0] = p[0];
421 a[1] = p01;
422 a[2] = blLerp(p01, p12, t);
423 a[3] = blLerp(a[2], blLerp(p12, p23, t), t);
424}
425
426static BL_INLINE void blSplitCubicAfter(const BLPoint p[4], BLPoint b[4], double t) noexcept {
427 BLPoint p01(blLerp(p[0], p[1], t));
428 BLPoint p12(blLerp(p[1], p[2], t));
429 BLPoint p23(blLerp(p[2], p[3], t));
430
431 b[3] = p[3];
432 b[2] = p23;
433 b[1] = blLerp(p12, p23, t);
434 b[0] = blLerp(blLerp(p01, p12, t), b[1], t);
435}
436
437static BL_INLINE void blGetCubicCoefficients(const BLPoint p[4], BLPoint& a, BLPoint& b, BLPoint& c, BLPoint& d) noexcept {
438 BLPoint v1 = p[1] - p[0];
439 BLPoint v2 = p[2] - p[1];
440 BLPoint v3 = p[3] - p[2];
441
442 a = v3 - v2 - v2 + v1;
443 b = 3.0 * (v2 - v1);
444 c = 3.0 * v1;
445 d = p[0];
446}
447
448static BL_INLINE void blGetCubicDerivativeCoefficients(const BLPoint p[4], BLPoint& a, BLPoint& b, BLPoint& c) noexcept {
449 BLPoint v1 = p[1] - p[0];
450 BLPoint v2 = p[2] - p[1];
451 BLPoint v3 = p[3] - p[2];
452
453 a = 3.0 * (v3 - v2 - v2 + v1);
454 b = 6.0 * (v2 - v1);
455 c = 3.0 * v1;
456}
457
458static BL_INLINE BLPoint blGetCubicValueAt(const BLPoint p[4], double t) noexcept {
459 BLPoint a, b, c, d;
460 blGetCubicCoefficients(p, a, b, c, d);
461 return ((a * t + b) * t + c) * t + d;
462}
463
464static BL_INLINE BLPoint blGetCubicValueAt(const BLPoint p[4], const BLPoint& t) noexcept {
465 BLPoint a, b, c, d;
466 blGetCubicCoefficients(p, a, b, c, d);
467 return ((a * t + b) * t + c) * t + d;
468}
469
470static BL_INLINE BLPoint blGetPreciseCubicValueAt(const BLPoint p[4], double t) noexcept {
471 BLPoint p01(blLerp(p[0], p[1], t));
472 BLPoint p12(blLerp(p[1], p[2], t));
473 BLPoint p23(blLerp(p[2], p[3], t));
474
475 return blLerp(blLerp(p01, p12, t), blLerp(p12, p23, t), t);
476}
477
478static BL_INLINE BLPoint blGetPreciseCubicValueAt(const BLPoint p[4], const BLPoint& t) noexcept {
479 BLPoint p01(blLerp(p[0], p[1], t));
480 BLPoint p12(blLerp(p[1], p[2], t));
481 BLPoint p23(blLerp(p[2], p[3], t));
482 return blLerp(blLerp(p01, p12, t), blLerp(p12, p23, t), t);
483}
484
485static BL_INLINE BLPoint blGetCubicDerivativeAt(const BLPoint p[4], double t) noexcept {
486 BLPoint p01 = blLerp(p[0], p[1], t);
487 BLPoint p12 = blLerp(p[1], p[2], t);
488 BLPoint p23 = blLerp(p[2], p[3], t);
489
490 return 3.0 * (blLerp(p12, p23, t) - blLerp(p01, p12, t));
491}
492
493static BL_INLINE void blGetCubicExtremaPoints(const BLPoint p[4], BLPoint out[2]) noexcept {
494 BLPoint a, b, c;
495 blGetCubicDerivativeCoefficients(p, a, b, c);
496
497 BLPoint t[2];
498 blSimplifiedQuadRoots(t, a, b, c);
499
500 t[0] = blClamp(t[0], 0.0, 1.0);
501 t[1] = blClamp(t[1], 0.0, 1.0);
502
503 out[0] = blGetPreciseCubicValueAt(p, t[0]);
504 out[1] = blGetPreciseCubicValueAt(p, t[1]);
505}
506
507static BL_INLINE BLPoint blCubicMidPoint(const BLPoint p[4]) noexcept {
508 return (p[0] + p[3]) * 0.125 + (p[1] + p[2]) * 0.375;
509}
510
511static BL_INLINE BLPoint blGetCubicIdentity(const BLPoint p[4]) noexcept {
512 BLPoint v1 = p[1] - p[0];
513 BLPoint v2 = p[2] - p[1];
514 BLPoint v3 = p[3] - p[2];
515
516 return v3 - v2 - v2 + v1;
517}
518
519static BL_INLINE bool blIsCubicFlat(const BLPoint p[3], double f) {
520 if (p[3] == p[0]) {
521 BLPoint v = p[2] - p[1];
522 double a = blCrossProduct(v, p[1] - p[0]);
523 return 0.5625 * a * a <= f * f * blLengthSq(v);
524 }
525 else {
526 BLPoint v = p[3] - p[0];
527 double a1 = blCrossProduct(v, p[1] - p[0]);
528 double a2 = blCrossProduct(v, p[2] - p[0]);
529 return 0.5625 * blMax(a1 * a1, a2 * a2) <= f * f * blLengthSq(v);
530 }
531}
532
533static BL_INLINE void blGetCubicInflectionParameter(const BLPoint p[4], double& tc, double& tl) noexcept {
534 BLPoint a, b, c;
535 blGetCubicDerivativeCoefficients(p, a, b, c);
536
537 // To get the inflections C'(t) cross C''(t) = at^2 + bt + c = 0 needs to be solved for 't'.
538 // The first cooefficient of the quadratic formula is also the denominator.
539 double den = blCrossProduct(b, a);
540
541 if (den != 0) {
542 // Two roots might exist, solve with quadratic formula ('tl' is real).
543 tc = blCrossProduct(a, c) / den;
544 tl = tc * tc + blCrossProduct(b, c) / den;
545
546 // If 'tl < 0' there are two complex roots (no need to solve).
547 // If 'tl == 0' there is a real double root at tc (cusp case).
548 // If 'tl > 0' two real roots exist at 'tc - Sqrt(tl)' and 'tc + Sqrt(tl)'.
549 if (tl > 0)
550 tl = blSqrt(tl);
551 }
552 else {
553 // One real root might exist, solve linear case ('tl' is NaN).
554 tc = -0.5 * blCrossProduct(c, b) / blCrossProduct(c, a);
555 tl = blNaN<double>();
556 }
557}
558
559static BL_INLINE BLPoint blGetCubicStartTangent(const BLPoint p[4]) noexcept {
560 BLPoint out = p[1] - p[0];
561 BLPoint t20 = p[2] - p[0];
562 BLPoint t30 = p[3] - p[0];
563
564 if (blIsZero(out)) out = t20;
565 if (blIsZero(out)) out = t30;
566
567 return out;
568}
569
570static BL_INLINE BLPoint blGetCubicEndTangent(const BLPoint p[4]) noexcept {
571 BLPoint out = p[3] - p[2];
572 BLPoint t31 = p[3] - p[1];
573 BLPoint t30 = p[3] - p[0];
574
575 if (blIsZero(out)) out = t31;
576 if (blIsZero(out)) out = t30;
577
578 return out;
579}
580
581static BL_INLINE void blApproximateCubicWithTwoQuads(const BLPoint p[4], BLPoint quads[7]) noexcept {
582 BLPoint c1 = blLerp(p[0], p[1], 0.75);
583 BLPoint c2 = blLerp(p[3], p[2], 0.75);
584 BLPoint pm = blLerp(c1, c2);
585
586 if (c1 == p[0])
587 c1 = blGetLineVectorIntersection(p[0], blGetCubicStartTangent(p), pm, blGetCubicDerivativeAt(p, 0.5));
588
589 if (c2 == p[3])
590 c2 = blGetLineVectorIntersection(p[3], blGetCubicEndTangent(p), pm, blGetCubicDerivativeAt(p, 0.5));
591
592 quads[0] = p[0];
593 quads[1] = c1;
594 quads[2] = pm;
595 quads[3] = c2;
596 quads[4] = p[3];
597}
598
599template<typename Callback>
600static BL_INLINE BLResult blApproximateCubicWithQuads(const BLPoint p[4], double simplifyTolerance, const Callback& callback) noexcept {
601 // Tolerance consists of a prefactor (27/4 * 2^3) combined with `simplifyTolerance`.
602 double toleranceSq = blSquare(54.0 * simplifyTolerance);
603
604 // Smallest parameter step to satisfy tolerance condition.
605 double t = blPow(toleranceSq / blLengthSq(blGetCubicIdentity(p)), 1.0 / 6.0);
606
607 BLPoint cubic[7];
608 cubic[3] = p[0];
609 cubic[4] = p[1];
610 cubic[5] = p[2];
611 cubic[6] = p[3];
612
613 for (;;) {
614 BLPoint quads[5];
615 t = blMin(t, 1.0);
616
617 // Split the cubic:
618 // - cubic[0:3] contains the part before `t`.
619 // - cubic[3:7] contains the part after `t`.
620 blSplitCubic(cubic + 3, cubic, cubic + 3, t);
621 blApproximateCubicWithTwoQuads(cubic, quads);
622
623 for (size_t i = 0; i <= 2; i += 2)
624 BL_PROPAGATE(callback(quads + i));
625
626 if (t >= 1.0)
627 return BL_SUCCESS;
628
629 // Recalculate the parameter.
630 double oldT = t;
631 t = t / (1.0 - t);
632
633 if (oldT - t < 1e-3)
634 t += 0.01;
635 }
636}
637
638//! \}
639//! \endcond
640
641#endif // BLEND2D_BLGEOMETRY_P_H
642