1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 /*
3  * This file is part of the LibreOffice project.
4  *
5  * This Source Code Form is subject to the terms of the Mozilla Public
6  * License, v. 2.0. If a copy of the MPL was not distributed with this
7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8  *
9  * This file incorporates work covered by the following license notice:
10  *
11  *   Licensed to the Apache Software Foundation (ASF) under one or more
12  *   contributor license agreements. See the NOTICE file distributed
13  *   with this work for additional information regarding copyright
14  *   ownership. The ASF licenses this file to you under the Apache
15  *   License, Version 2.0 (the "License"); you may not use this file
16  *   except in compliance with the License. You may obtain a copy of
17  *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
18  */
19 
20 #include <basegfx/curve/b2dcubicbezier.hxx>
21 #include <basegfx/vector/b2dvector.hxx>
22 #include <basegfx/polygon/b2dpolygon.hxx>
23 #include <basegfx/matrix/b2dhommatrix.hxx>
24 #include <basegfx/numeric/ftools.hxx>
25 
26 #include <osl/diagnose.h>
27 
28 #include <cmath>
29 #include <limits>
30 
31 // #i37443#
32 #define FACTOR_FOR_UNSHARPEN    (1.6)
33 #ifdef DBG_UTIL
34 const double fMultFactUnsharpen = FACTOR_FOR_UNSHARPEN;
35 #endif
36 
37 namespace basegfx
38 {
39     namespace
40     {
ImpSubDivAngle(const B2DPoint & rfPA,const B2DPoint & rfEA,const B2DPoint & rfEB,const B2DPoint & rfPB,B2DPolygon & rTarget,double fAngleBound,bool bAllowUnsharpen,sal_uInt16 nMaxRecursionDepth)41         void ImpSubDivAngle(
42             const B2DPoint& rfPA,           // start point
43             const B2DPoint& rfEA,           // edge on A
44             const B2DPoint& rfEB,           // edge on B
45             const B2DPoint& rfPB,           // end point
46             B2DPolygon& rTarget,            // target polygon
47             double fAngleBound,             // angle bound in [0.0 .. 2PI]
48             bool bAllowUnsharpen,           // #i37443# allow the criteria to get unsharp in recursions
49             sal_uInt16 nMaxRecursionDepth)  // endless loop protection
50         {
51             if(nMaxRecursionDepth)
52             {
53                 // do angle test
54                 B2DVector aLeft(rfEA - rfPA);
55                 B2DVector aRight(rfEB - rfPB);
56 
57                 // #i72104#
58                 if(aLeft.equalZero())
59                 {
60                     aLeft = rfEB - rfPA;
61                 }
62 
63                 if(aRight.equalZero())
64                 {
65                     aRight = rfEA - rfPB;
66                 }
67 
68                 const double fCurrentAngle(aLeft.angle(aRight));
69 
70                 if(fabs(fCurrentAngle) > (M_PI - fAngleBound))
71                 {
72                     // end recursion
73                     nMaxRecursionDepth = 0;
74                 }
75                 else
76                 {
77                     if(bAllowUnsharpen)
78                     {
79                         // #i37443# unsharpen criteria
80 #ifdef DBG_UTIL
81                         fAngleBound *= fMultFactUnsharpen;
82 #else
83                         fAngleBound *= FACTOR_FOR_UNSHARPEN;
84 #endif
85                     }
86                 }
87             }
88 
89             if(nMaxRecursionDepth)
90             {
91                 // divide at 0.5
92                 const B2DPoint aS1L(average(rfPA, rfEA));
93                 const B2DPoint aS1C(average(rfEA, rfEB));
94                 const B2DPoint aS1R(average(rfEB, rfPB));
95                 const B2DPoint aS2L(average(aS1L, aS1C));
96                 const B2DPoint aS2R(average(aS1C, aS1R));
97                 const B2DPoint aS3C(average(aS2L, aS2R));
98 
99                 // left recursion
100                 ImpSubDivAngle(rfPA, aS1L, aS2L, aS3C, rTarget, fAngleBound, bAllowUnsharpen, nMaxRecursionDepth - 1);
101 
102                 // right recursion
103                 ImpSubDivAngle(aS3C, aS2R, aS1R, rfPB, rTarget, fAngleBound, bAllowUnsharpen, nMaxRecursionDepth - 1);
104             }
105             else
106             {
107                 rTarget.append(rfPB);
108             }
109         }
110 
ImpSubDivAngleStart(const B2DPoint & rfPA,const B2DPoint & rfEA,const B2DPoint & rfEB,const B2DPoint & rfPB,B2DPolygon & rTarget,const double & rfAngleBound)111         void ImpSubDivAngleStart(
112             const B2DPoint& rfPA,           // start point
113             const B2DPoint& rfEA,           // edge on A
114             const B2DPoint& rfEB,           // edge on B
115             const B2DPoint& rfPB,           // end point
116             B2DPolygon& rTarget,            // target polygon
117             const double& rfAngleBound)     // angle bound in [0.0 .. 2PI]
118         {
119             sal_uInt16 nMaxRecursionDepth(8);
120             const B2DVector aLeft(rfEA - rfPA);
121             const B2DVector aRight(rfEB - rfPB);
122             bool bLeftEqualZero(aLeft.equalZero());
123             bool bRightEqualZero(aRight.equalZero());
124             bool bAllParallel(false);
125 
126             if(bLeftEqualZero && bRightEqualZero)
127             {
128                 nMaxRecursionDepth = 0;
129             }
130             else
131             {
132                 const B2DVector aBase(rfPB - rfPA);
133                 const bool bBaseEqualZero(aBase.equalZero()); // #i72104#
134 
135                 if(!bBaseEqualZero)
136                 {
137                     const bool bLeftParallel(bLeftEqualZero || areParallel(aLeft, aBase));
138                     const bool bRightParallel(bRightEqualZero || areParallel(aRight, aBase));
139 
140                     if(bLeftParallel && bRightParallel)
141                     {
142                         bAllParallel = true;
143 
144                         if(!bLeftEqualZero)
145                         {
146                             double fFactor;
147 
148                             if(fabs(aBase.getX()) > fabs(aBase.getY()))
149                             {
150                                 fFactor = aLeft.getX() / aBase.getX();
151                             }
152                             else
153                             {
154                                 fFactor = aLeft.getY() / aBase.getY();
155                             }
156 
157                             if(fFactor >= 0.0 && fFactor <= 1.0)
158                             {
159                                 bLeftEqualZero = true;
160                             }
161                         }
162 
163                         if(!bRightEqualZero)
164                         {
165                             double fFactor;
166 
167                             if(fabs(aBase.getX()) > fabs(aBase.getY()))
168                             {
169                                 fFactor = aRight.getX() / -aBase.getX();
170                             }
171                             else
172                             {
173                                 fFactor = aRight.getY() / -aBase.getY();
174                             }
175 
176                             if(fFactor >= 0.0 && fFactor <= 1.0)
177                             {
178                                 bRightEqualZero = true;
179                             }
180                         }
181 
182                         if(bLeftEqualZero && bRightEqualZero)
183                         {
184                             nMaxRecursionDepth = 0;
185                         }
186                     }
187                 }
188             }
189 
190             if(nMaxRecursionDepth)
191             {
192                 // divide at 0.5 ad test both edges for angle criteria
193                 const B2DPoint aS1L(average(rfPA, rfEA));
194                 const B2DPoint aS1C(average(rfEA, rfEB));
195                 const B2DPoint aS1R(average(rfEB, rfPB));
196                 const B2DPoint aS2L(average(aS1L, aS1C));
197                 const B2DPoint aS2R(average(aS1C, aS1R));
198                 const B2DPoint aS3C(average(aS2L, aS2R));
199 
200                 // test left
201                 bool bAngleIsSmallerLeft(bAllParallel && bLeftEqualZero);
202                 if(!bAngleIsSmallerLeft)
203                 {
204                     const B2DVector aLeftLeft(bLeftEqualZero ? aS2L - aS1L : aS1L - rfPA); // #i72104#
205                     const B2DVector aRightLeft(aS2L - aS3C);
206                     const double fCurrentAngleLeft(aLeftLeft.angle(aRightLeft));
207                     bAngleIsSmallerLeft = (fabs(fCurrentAngleLeft) > (M_PI - rfAngleBound));
208                 }
209 
210                 // test right
211                 bool bAngleIsSmallerRight(bAllParallel && bRightEqualZero);
212                 if(!bAngleIsSmallerRight)
213                 {
214                     const B2DVector aLeftRight(aS2R - aS3C);
215                     const B2DVector aRightRight(bRightEqualZero ? aS2R - aS1R : aS1R - rfPB); // #i72104#
216                     const double fCurrentAngleRight(aLeftRight.angle(aRightRight));
217                     bAngleIsSmallerRight = (fabs(fCurrentAngleRight) > (M_PI - rfAngleBound));
218                 }
219 
220                 if(bAngleIsSmallerLeft && bAngleIsSmallerRight)
221                 {
222                     // no recursion necessary at all
223                     nMaxRecursionDepth = 0;
224                 }
225                 else
226                 {
227                     // left
228                     if(bAngleIsSmallerLeft)
229                     {
230                         rTarget.append(aS3C);
231                     }
232                     else
233                     {
234                         ImpSubDivAngle(rfPA, aS1L, aS2L, aS3C, rTarget, rfAngleBound, true/*bAllowUnsharpen*/, nMaxRecursionDepth);
235                     }
236 
237                     // right
238                     if(bAngleIsSmallerRight)
239                     {
240                         rTarget.append(rfPB);
241                     }
242                     else
243                     {
244                         ImpSubDivAngle(aS3C, aS2R, aS1R, rfPB, rTarget, rfAngleBound, true/*bAllowUnsharpen*/, nMaxRecursionDepth);
245                     }
246                 }
247             }
248 
249             if(!nMaxRecursionDepth)
250             {
251                 rTarget.append(rfPB);
252             }
253         }
254 
ImpSubDivDistance(const B2DPoint & rfPA,const B2DPoint & rfEA,const B2DPoint & rfEB,const B2DPoint & rfPB,B2DPolygon & rTarget,double fDistanceBound2,double fLastDistanceError2,sal_uInt16 nMaxRecursionDepth)255         void ImpSubDivDistance(
256             const B2DPoint& rfPA,           // start point
257             const B2DPoint& rfEA,           // edge on A
258             const B2DPoint& rfEB,           // edge on B
259             const B2DPoint& rfPB,           // end point
260             B2DPolygon& rTarget,            // target polygon
261             double fDistanceBound2,         // quadratic distance criteria
262             double fLastDistanceError2,     // the last quadratic distance error
263             sal_uInt16 nMaxRecursionDepth)  // endless loop protection
264         {
265             if(nMaxRecursionDepth)
266             {
267                 // decide if another recursion is needed. If not, set
268                 // nMaxRecursionDepth to zero
269 
270                 // Perform bezier flatness test (lecture notes from R. Schaback,
271                 // Mathematics of Computer-Aided Design, Uni Goettingen, 2000)
272 
273                 // ||P(t) - L(t)|| <= max     ||b_j - b_0 - j/n(b_n - b_0)||
274                 //                    0<=j<=n
275 
276                 // What is calculated here is an upper bound to the distance from
277                 // a line through b_0 and b_3 (rfPA and P4 in our notation) and the
278                 // curve. We can drop 0 and n from the running indices, since the
279                 // argument of max becomes zero for those cases.
280                 const double fJ1x(rfEA.getX() - rfPA.getX() - 1.0/3.0*(rfPB.getX() - rfPA.getX()));
281                 const double fJ1y(rfEA.getY() - rfPA.getY() - 1.0/3.0*(rfPB.getY() - rfPA.getY()));
282                 const double fJ2x(rfEB.getX() - rfPA.getX() - 2.0/3.0*(rfPB.getX() - rfPA.getX()));
283                 const double fJ2y(rfEB.getY() - rfPA.getY() - 2.0/3.0*(rfPB.getY() - rfPA.getY()));
284                 const double fDistanceError2(std::max(fJ1x*fJ1x + fJ1y*fJ1y, fJ2x*fJ2x + fJ2y*fJ2y));
285 
286                 // stop if error measure does not improve anymore. This is a
287                 // safety guard against floating point inaccuracies.
288                 // stop if distance from line is guaranteed to be bounded by d
289                 const bool bFurtherDivision(fLastDistanceError2 > fDistanceError2 && fDistanceError2 >= fDistanceBound2);
290 
291                 if(bFurtherDivision)
292                 {
293                     // remember last error value
294                     fLastDistanceError2 = fDistanceError2;
295                 }
296                 else
297                 {
298                     // stop recursion
299                     nMaxRecursionDepth = 0;
300                 }
301             }
302 
303             if(nMaxRecursionDepth)
304             {
305                 // divide at 0.5
306                 const B2DPoint aS1L(average(rfPA, rfEA));
307                 const B2DPoint aS1C(average(rfEA, rfEB));
308                 const B2DPoint aS1R(average(rfEB, rfPB));
309                 const B2DPoint aS2L(average(aS1L, aS1C));
310                 const B2DPoint aS2R(average(aS1C, aS1R));
311                 const B2DPoint aS3C(average(aS2L, aS2R));
312 
313                 // left recursion
314                 ImpSubDivDistance(rfPA, aS1L, aS2L, aS3C, rTarget, fDistanceBound2, fLastDistanceError2, nMaxRecursionDepth - 1);
315 
316                 // right recursion
317                 ImpSubDivDistance(aS3C, aS2R, aS1R, rfPB, rTarget, fDistanceBound2, fLastDistanceError2, nMaxRecursionDepth - 1);
318             }
319             else
320             {
321                 rTarget.append(rfPB);
322             }
323         }
324     } // end of anonymous namespace
325 } // end of namespace basegfx
326 
327 namespace basegfx
328 {
329     B2DCubicBezier::B2DCubicBezier(const B2DCubicBezier&) = default;
330 
331     B2DCubicBezier::B2DCubicBezier() = default;
332 
B2DCubicBezier(const B2DPoint & rStart,const B2DPoint & rControlPointA,const B2DPoint & rControlPointB,const B2DPoint & rEnd)333     B2DCubicBezier::B2DCubicBezier(const B2DPoint& rStart, const B2DPoint& rControlPointA, const B2DPoint& rControlPointB, const B2DPoint& rEnd)
334     :   maStartPoint(rStart),
335         maEndPoint(rEnd),
336         maControlPointA(rControlPointA),
337         maControlPointB(rControlPointB)
338     {
339     }
340 
341     // assignment operator
342     B2DCubicBezier& B2DCubicBezier::operator=(const B2DCubicBezier&) = default;
343 
344     // compare operators
operator ==(const B2DCubicBezier & rBezier) const345     bool B2DCubicBezier::operator==(const B2DCubicBezier& rBezier) const
346     {
347         return (
348             maStartPoint == rBezier.maStartPoint
349             && maEndPoint == rBezier.maEndPoint
350             && maControlPointA == rBezier.maControlPointA
351             && maControlPointB == rBezier.maControlPointB
352         );
353     }
354 
equal(const B2DCubicBezier & rBezier) const355     bool B2DCubicBezier::equal(const B2DCubicBezier& rBezier) const
356     {
357         return (
358             maStartPoint.equal(rBezier.maStartPoint)
359             && maEndPoint.equal(rBezier.maEndPoint)
360             && maControlPointA.equal(rBezier.maControlPointA)
361             && maControlPointB.equal(rBezier.maControlPointB)
362         );
363     }
364 
365     // test if vectors are used
isBezier() const366     bool B2DCubicBezier::isBezier() const
367     {
368         return maControlPointA != maStartPoint || maControlPointB != maEndPoint;
369     }
370 
testAndSolveTrivialBezier()371     void B2DCubicBezier::testAndSolveTrivialBezier()
372     {
373         if(maControlPointA == maStartPoint && maControlPointB == maEndPoint)
374             return;
375 
376         const B2DVector aEdge(maEndPoint - maStartPoint);
377 
378         // controls parallel to edge can be trivial. No edge -> not parallel -> control can
379         // still not be trivial (e.g. ballon loop)
380         if(aEdge.equalZero())
381             return;
382 
383         // get control vectors
384         const B2DVector aVecA(maControlPointA - maStartPoint);
385         const B2DVector aVecB(maControlPointB - maEndPoint);
386 
387         // check if trivial per se
388         bool bAIsTrivial(aVecA.equalZero());
389         bool bBIsTrivial(aVecB.equalZero());
390 
391         // #i102241# prepare inverse edge length to normalize cross values;
392         // else the small compare value used in fTools::equalZero
393         // will be length dependent and this detection will work as less
394         // precise as longer the edge is. In principle, the length of the control
395         // vector would need to be used too, but to be trivial it is assumed to
396         // be of roughly equal length to the edge, so edge length can be used
397         // for both. Only needed when one of both is not trivial per se.
398         const double fInverseEdgeLength(bAIsTrivial && bBIsTrivial
399             ? 1.0
400             : 1.0 / aEdge.getLength());
401 
402         // if A is not zero, check if it could be
403         if(!bAIsTrivial)
404         {
405             // #i102241# parallel to edge? Check aVecA, aEdge. Use cross() which does what
406             // we need here with the precision we need
407             const double fCross(aVecA.cross(aEdge) * fInverseEdgeLength);
408 
409             if(fTools::equalZero(fCross))
410             {
411                 // get scale to edge. Use bigger distance for numeric quality
412                 const double fScale(fabs(aEdge.getX()) > fabs(aEdge.getY())
413                     ? aVecA.getX() / aEdge.getX()
414                     : aVecA.getY() / aEdge.getY());
415 
416                 // relative end point of vector in edge range?
417                 if (fTools::betweenOrEqualEither(fScale, 0.0, 1.0))
418                 {
419                     bAIsTrivial = true;
420                 }
421             }
422         }
423 
424         // if B is not zero, check if it could be, but only if A is already trivial;
425         // else solve to trivial will not be possible for whole edge
426         if(bAIsTrivial && !bBIsTrivial)
427         {
428             // parallel to edge? Check aVecB, aEdge
429             const double fCross(aVecB.cross(aEdge) * fInverseEdgeLength);
430 
431             if(fTools::equalZero(fCross))
432             {
433                 // get scale to edge. Use bigger distance for numeric quality
434                 const double fScale(fabs(aEdge.getX()) > fabs(aEdge.getY())
435                     ? aVecB.getX() / aEdge.getX()
436                     : aVecB.getY() / aEdge.getY());
437 
438                 // end point of vector in edge range? Caution: controlB is directed AGAINST edge
439                 if (fTools::betweenOrEqualEither(fScale, -1.0, 0.0))
440                 {
441                     bBIsTrivial = true;
442                 }
443             }
444         }
445 
446         // if both are/can be reduced, do it.
447         // Not possible if only one is/can be reduced (!)
448         if(bAIsTrivial && bBIsTrivial)
449         {
450             maControlPointA = maStartPoint;
451             maControlPointB = maEndPoint;
452         }
453     }
454 
455     namespace {
impGetLength(const B2DCubicBezier & rEdge,double fDeviation,sal_uInt32 nRecursionWatch)456         double impGetLength(const B2DCubicBezier& rEdge, double fDeviation, sal_uInt32 nRecursionWatch)
457         {
458             const double fEdgeLength(rEdge.getEdgeLength());
459             const double fControlPolygonLength(rEdge.getControlPolygonLength());
460             const double fCurrentDeviation(fTools::equalZero(fControlPolygonLength) ? 0.0 : 1.0 - (fEdgeLength / fControlPolygonLength));
461 
462             if(!nRecursionWatch || fTools:: lessOrEqual(fCurrentDeviation, fDeviation))
463             {
464                 return (fEdgeLength + fControlPolygonLength) * 0.5;
465             }
466             else
467             {
468                 B2DCubicBezier aLeft, aRight;
469                 const double fNewDeviation(fDeviation * 0.5);
470                 const sal_uInt32 nNewRecursionWatch(nRecursionWatch - 1);
471 
472                 rEdge.split(0.5, &aLeft, &aRight);
473 
474                 return impGetLength(aLeft, fNewDeviation, nNewRecursionWatch)
475                     + impGetLength(aRight, fNewDeviation, nNewRecursionWatch);
476             }
477         }
478     }
479 
getLength(double fDeviation) const480     double B2DCubicBezier::getLength(double fDeviation) const
481     {
482         if(isBezier())
483         {
484             if(fDeviation < 0.00000001)
485             {
486                 fDeviation = 0.00000001;
487             }
488 
489             return impGetLength(*this, fDeviation, 6);
490         }
491         else
492         {
493             return B2DVector(getEndPoint() - getStartPoint()).getLength();
494         }
495     }
496 
getEdgeLength() const497     double B2DCubicBezier::getEdgeLength() const
498     {
499         const B2DVector aEdge(maEndPoint - maStartPoint);
500         return aEdge.getLength();
501     }
502 
getControlPolygonLength() const503     double B2DCubicBezier::getControlPolygonLength() const
504     {
505         const B2DVector aVectorA(maControlPointA - maStartPoint);
506         const B2DVector aVectorB(maEndPoint - maControlPointB);
507 
508         if(!aVectorA.equalZero() || !aVectorB.equalZero())
509         {
510             const B2DVector aTop(maControlPointB - maControlPointA);
511             return (aVectorA.getLength() + aVectorB.getLength() + aTop.getLength());
512         }
513         else
514         {
515             return getEdgeLength();
516         }
517     }
518 
adaptiveSubdivideByAngle(B2DPolygon & rTarget,double fAngleBound) const519     void B2DCubicBezier::adaptiveSubdivideByAngle(B2DPolygon& rTarget, double fAngleBound) const
520     {
521         if(isBezier())
522         {
523             // use support method #i37443# and allow unsharpen the criteria
524             ImpSubDivAngleStart(maStartPoint, maControlPointA, maControlPointB, maEndPoint, rTarget,
525                                 deg2rad(fAngleBound));
526         }
527         else
528         {
529             rTarget.append(getEndPoint());
530         }
531     }
532 
getTangent(double t) const533     B2DVector B2DCubicBezier::getTangent(double t) const
534     {
535         if(t <= 0.0)
536         {
537             // tangent in start point
538             B2DVector aTangent(getControlPointA() - getStartPoint());
539 
540             if(!aTangent.equalZero())
541             {
542                 return aTangent;
543             }
544 
545             // start point and control vector are the same, fallback
546             // to implicit start vector to control point B
547             aTangent = (getControlPointB() - getStartPoint()) * 0.3;
548 
549             if(!aTangent.equalZero())
550             {
551                 return aTangent;
552             }
553 
554             // not a bezier at all, return edge vector
555             return (getEndPoint() - getStartPoint()) * 0.3;
556         }
557         else if(fTools::moreOrEqual(t, 1.0))
558         {
559             // tangent in end point
560             B2DVector aTangent(getEndPoint() - getControlPointB());
561 
562             if(!aTangent.equalZero())
563             {
564                 return aTangent;
565             }
566 
567             // end point and control vector are the same, fallback
568             // to implicit start vector from control point A
569             aTangent = (getEndPoint() - getControlPointA()) * 0.3;
570 
571             if(!aTangent.equalZero())
572             {
573                 return aTangent;
574             }
575 
576             // not a bezier at all, return edge vector
577             return (getEndPoint() - getStartPoint()) * 0.3;
578         }
579         else
580         {
581             // t is in ]0.0 .. 1.0[. Split and extract
582             B2DCubicBezier aRight;
583             split(t, nullptr, &aRight);
584 
585             return aRight.getControlPointA() - aRight.getStartPoint();
586         }
587     }
588 
589     // #i37443# adaptive subdivide by nCount subdivisions
adaptiveSubdivideByCount(B2DPolygon & rTarget,sal_uInt32 nCount) const590     void B2DCubicBezier::adaptiveSubdivideByCount(B2DPolygon& rTarget, sal_uInt32 nCount) const
591     {
592         const double fLenFact(1.0 / static_cast< double >(nCount + 1));
593 
594         for(sal_uInt32 a(1); a <= nCount; a++)
595         {
596             const double fPos(static_cast< double >(a) * fLenFact);
597             rTarget.append(interpolatePoint(fPos));
598         }
599 
600         rTarget.append(getEndPoint());
601     }
602 
603     // adaptive subdivide by distance
adaptiveSubdivideByDistance(B2DPolygon & rTarget,double fDistanceBound,int nRecurseLimit) const604     void B2DCubicBezier::adaptiveSubdivideByDistance(B2DPolygon& rTarget, double fDistanceBound, int nRecurseLimit) const
605     {
606         if(isBezier())
607         {
608             ImpSubDivDistance(maStartPoint, maControlPointA, maControlPointB, maEndPoint, rTarget,
609                 fDistanceBound * fDistanceBound, std::numeric_limits<double>::max(), nRecurseLimit);
610         }
611         else
612         {
613             rTarget.append(getEndPoint());
614         }
615     }
616 
interpolatePoint(double t) const617     B2DPoint B2DCubicBezier::interpolatePoint(double t) const
618     {
619         OSL_ENSURE(t >= 0.0 && t <= 1.0, "B2DCubicBezier::interpolatePoint: Access out of range (!)");
620 
621         if(isBezier())
622         {
623             const B2DPoint aS1L(interpolate(maStartPoint, maControlPointA, t));
624             const B2DPoint aS1C(interpolate(maControlPointA, maControlPointB, t));
625             const B2DPoint aS1R(interpolate(maControlPointB, maEndPoint, t));
626             const B2DPoint aS2L(interpolate(aS1L, aS1C, t));
627             const B2DPoint aS2R(interpolate(aS1C, aS1R, t));
628 
629             return interpolate(aS2L, aS2R, t);
630         }
631         else
632         {
633             return interpolate(maStartPoint, maEndPoint, t);
634         }
635     }
636 
getSmallestDistancePointToBezierSegment(const B2DPoint & rTestPoint,double & rCut) const637     double B2DCubicBezier::getSmallestDistancePointToBezierSegment(const B2DPoint& rTestPoint, double& rCut) const
638     {
639         const sal_uInt32 nInitialDivisions(3);
640         B2DPolygon aInitialPolygon;
641 
642         // as start make a fix division, creates nInitialDivisions + 2 points
643         aInitialPolygon.append(getStartPoint());
644         adaptiveSubdivideByCount(aInitialPolygon, nInitialDivisions);
645 
646         // now look for the closest point
647         const sal_uInt32 nPointCount(aInitialPolygon.count());
648         B2DVector aVector(rTestPoint - aInitialPolygon.getB2DPoint(0));
649         double pointDistance(std::hypot(aVector.getX(), aVector.getY()));
650         double newPointDistance;
651         sal_uInt32 nSmallestIndex(0);
652 
653         for(sal_uInt32 a(1); a < nPointCount; a++)
654         {
655             aVector = B2DVector(rTestPoint - aInitialPolygon.getB2DPoint(a));
656             newPointDistance = std::hypot(aVector.getX(), aVector.getY());
657             if(newPointDistance < pointDistance)
658             {
659                 pointDistance = newPointDistance;
660                 nSmallestIndex = a;
661             }
662         }
663 
664         // look right and left for even smaller distances
665         double fStepValue(1.0 / static_cast<double>((nPointCount - 1) * 2)); // half the edge step width
666         double fPosition(static_cast<double>(nSmallestIndex) / static_cast<double>(nPointCount - 1));
667 
668         while(true)
669         {
670             // test left
671             double fPosLeft(fPosition - fStepValue);
672 
673             if(fPosLeft < 0.0)
674             {
675                 fPosLeft = 0.0;
676                 aVector = B2DVector(rTestPoint - maStartPoint);
677             }
678             else
679             {
680                 aVector = B2DVector(rTestPoint - interpolatePoint(fPosLeft));
681             }
682 
683             newPointDistance = std::hypot(aVector.getX(), aVector.getY());
684 
685             if(fTools::less(newPointDistance, pointDistance))
686             {
687                 pointDistance = newPointDistance;
688                 fPosition = fPosLeft;
689             }
690             else
691             {
692                 // test right
693                 double fPosRight(fPosition + fStepValue);
694 
695                 if(fPosRight > 1.0)
696                 {
697                     fPosRight = 1.0;
698                     aVector = B2DVector(rTestPoint - maEndPoint);
699                 }
700                 else
701                 {
702                     aVector = B2DVector(rTestPoint - interpolatePoint(fPosRight));
703                 }
704 
705                 newPointDistance = std::hypot(aVector.getX(), aVector.getY());
706 
707                 if(fTools::less(newPointDistance, pointDistance))
708                 {
709                     pointDistance = newPointDistance;
710                     fPosition = fPosRight;
711                 }
712                 else
713                 {
714                     // not less left or right, done
715                     break;
716                 }
717             }
718 
719             if(fPosition == 0.0 || fPosition == 1.0)
720             {
721                 // if we are completely left or right, we are done
722                 break;
723             }
724 
725             // prepare next step value
726             fStepValue /= 2.0;
727         }
728 
729         rCut = fPosition;
730         return pointDistance;
731     }
732 
split(double t,B2DCubicBezier * pBezierA,B2DCubicBezier * pBezierB) const733     void B2DCubicBezier::split(double t, B2DCubicBezier* pBezierA, B2DCubicBezier* pBezierB) const
734     {
735         OSL_ENSURE(t >= 0.0 && t <= 1.0, "B2DCubicBezier::split: Access out of range (!)");
736 
737         if(!pBezierA && !pBezierB)
738         {
739             return;
740         }
741 
742         if(isBezier())
743         {
744             const B2DPoint aS1L(interpolate(maStartPoint, maControlPointA, t));
745             const B2DPoint aS1C(interpolate(maControlPointA, maControlPointB, t));
746             const B2DPoint aS1R(interpolate(maControlPointB, maEndPoint, t));
747             const B2DPoint aS2L(interpolate(aS1L, aS1C, t));
748             const B2DPoint aS2R(interpolate(aS1C, aS1R, t));
749             const B2DPoint aS3C(interpolate(aS2L, aS2R, t));
750 
751             if(pBezierA)
752             {
753                 pBezierA->setStartPoint(maStartPoint);
754                 pBezierA->setEndPoint(aS3C);
755                 pBezierA->setControlPointA(aS1L);
756                 pBezierA->setControlPointB(aS2L);
757             }
758 
759             if(pBezierB)
760             {
761                 pBezierB->setStartPoint(aS3C);
762                 pBezierB->setEndPoint(maEndPoint);
763                 pBezierB->setControlPointA(aS2R);
764                 pBezierB->setControlPointB(aS1R);
765             }
766         }
767         else
768         {
769             const B2DPoint aSplit(interpolate(maStartPoint, maEndPoint, t));
770 
771             if(pBezierA)
772             {
773                 pBezierA->setStartPoint(maStartPoint);
774                 pBezierA->setEndPoint(aSplit);
775                 pBezierA->setControlPointA(maStartPoint);
776                 pBezierA->setControlPointB(aSplit);
777             }
778 
779             if(pBezierB)
780             {
781                 pBezierB->setStartPoint(aSplit);
782                 pBezierB->setEndPoint(maEndPoint);
783                 pBezierB->setControlPointA(aSplit);
784                 pBezierB->setControlPointB(maEndPoint);
785             }
786         }
787     }
788 
snippet(double fStart,double fEnd) const789     B2DCubicBezier B2DCubicBezier::snippet(double fStart, double fEnd) const
790     {
791         B2DCubicBezier aRetval;
792 
793         fStart = std::clamp(fStart, 0.0, 1.0);
794         fEnd = std::clamp(fEnd, 0.0, 1.0);
795 
796         if(fEnd <= fStart)
797         {
798             // empty or NULL, create single point at center
799             const double fSplit((fEnd + fStart) * 0.5);
800             const B2DPoint aPoint(interpolate(getStartPoint(), getEndPoint(), fSplit));
801             aRetval.setStartPoint(aPoint);
802             aRetval.setEndPoint(aPoint);
803             aRetval.setControlPointA(aPoint);
804             aRetval.setControlPointB(aPoint);
805         }
806         else
807         {
808             if(isBezier())
809             {
810                 // copy bezier; cut off right, then cut off left. Do not forget to
811                 // adapt cut value when both cuts happen
812                 const bool bEndIsOne(fTools::equal(fEnd, 1.0));
813                 const bool bStartIsZero(fTools::equalZero(fStart));
814                 aRetval = *this;
815 
816                 if(!bEndIsOne)
817                 {
818                     aRetval.split(fEnd, &aRetval, nullptr);
819 
820                     if(!bStartIsZero)
821                     {
822                         fStart /= fEnd;
823                     }
824                 }
825 
826                 if(!bStartIsZero)
827                 {
828                     aRetval.split(fStart, nullptr, &aRetval);
829                 }
830             }
831             else
832             {
833                 // no bezier, create simple edge
834                 const B2DPoint aPointA(interpolate(getStartPoint(), getEndPoint(), fStart));
835                 const B2DPoint aPointB(interpolate(getStartPoint(), getEndPoint(), fEnd));
836                 aRetval.setStartPoint(aPointA);
837                 aRetval.setEndPoint(aPointB);
838                 aRetval.setControlPointA(aPointA);
839                 aRetval.setControlPointB(aPointB);
840             }
841         }
842 
843         return aRetval;
844     }
845 
getRange() const846     B2DRange B2DCubicBezier::getRange() const
847     {
848         B2DRange aRetval(maStartPoint, maEndPoint);
849 
850         aRetval.expand(maControlPointA);
851         aRetval.expand(maControlPointB);
852 
853         return aRetval;
854     }
855 
getMinimumExtremumPosition(double & rfResult) const856     bool B2DCubicBezier::getMinimumExtremumPosition(double& rfResult) const
857     {
858         std::vector< double > aAllResults;
859 
860         aAllResults.reserve(4);
861         getAllExtremumPositions(aAllResults);
862 
863         const sal_uInt32 nCount(aAllResults.size());
864 
865         if(!nCount)
866         {
867             return false;
868         }
869         else if(nCount == 1)
870         {
871             rfResult = aAllResults[0];
872             return true;
873         }
874         else
875         {
876             rfResult = *(std::min_element(aAllResults.begin(), aAllResults.end()));
877             return true;
878         }
879     }
880 
881     namespace
882     {
impCheckExtremumResult(double fCandidate,std::vector<double> & rResult)883         void impCheckExtremumResult(double fCandidate, std::vector< double >& rResult)
884         {
885             // check for range ]0.0 .. 1.0[ with excluding 1.0 and 0.0 clearly
886             // by using the equalZero test, NOT ::more or ::less which will use the
887             // ApproxEqual() which is too exact here
888             if(fCandidate > 0.0 && !fTools::equalZero(fCandidate))
889             {
890                 if(fCandidate < 1.0 && !fTools::equalZero(fCandidate - 1.0))
891                 {
892                     rResult.push_back(fCandidate);
893                 }
894             }
895         }
896     }
897 
getAllExtremumPositions(std::vector<double> & rResults) const898     void B2DCubicBezier::getAllExtremumPositions(std::vector< double >& rResults) const
899     {
900         rResults.clear();
901 
902         // calculate the x-extrema parameters by zeroing first x-derivative
903         // of the cubic bezier's parametric formula, which results in a
904         // quadratic equation: dBezier/dt = t*t*fAX - 2*t*fBX + fCX
905         const B2DPoint aControlDiff( maControlPointA - maControlPointB );
906         double fCX = maControlPointA.getX() - maStartPoint.getX();
907         const double fBX = fCX + aControlDiff.getX();
908         const double fAX = 3 * aControlDiff.getX() + (maEndPoint.getX() - maStartPoint.getX());
909 
910         if(fTools::equalZero(fCX))
911         {
912             // detect fCX equal zero and truncate to real zero value in that case
913             fCX = 0.0;
914         }
915 
916         if( !fTools::equalZero(fAX) )
917         {
918             // derivative is polynomial of order 2 => use binomial formula
919             const double fD = fBX*fBX - fAX*fCX;
920             if( fD >= 0.0 )
921             {
922                 const double fS = sqrt(fD);
923                 // calculate both roots (avoiding a numerically unstable subtraction)
924                 const double fQ = fBX + ((fBX >= 0) ? +fS : -fS);
925                 impCheckExtremumResult(fQ / fAX, rResults);
926                 if( !fTools::equalZero(fS) ) // ignore root multiplicity
927                     impCheckExtremumResult(fCX / fQ, rResults);
928             }
929         }
930         else if( !fTools::equalZero(fBX) )
931         {
932             // derivative is polynomial of order 1 => one extrema
933             impCheckExtremumResult(fCX / (2 * fBX), rResults);
934         }
935 
936         // calculate the y-extrema parameters by zeroing first y-derivative
937         double fCY = maControlPointA.getY() - maStartPoint.getY();
938         const double fBY = fCY + aControlDiff.getY();
939         const double fAY = 3 * aControlDiff.getY() + (maEndPoint.getY() - maStartPoint.getY());
940 
941         if(fTools::equalZero(fCY))
942         {
943             // detect fCY equal zero and truncate to real zero value in that case
944             fCY = 0.0;
945         }
946 
947         if( !fTools::equalZero(fAY) )
948         {
949             // derivative is polynomial of order 2 => use binomial formula
950             const double fD = fBY*fBY - fAY*fCY;
951             if( fD >= 0.0 )
952             {
953                 const double fS = sqrt(fD);
954                 // calculate both roots (avoiding a numerically unstable subtraction)
955                 const double fQ = fBY + ((fBY >= 0) ? +fS : -fS);
956                 impCheckExtremumResult(fQ / fAY, rResults);
957                 if( !fTools::equalZero(fS) ) // ignore root multiplicity
958                     impCheckExtremumResult(fCY / fQ, rResults);
959             }
960         }
961         else if( !fTools::equalZero(fBY) )
962         {
963             // derivative is polynomial of order 1 => one extrema
964             impCheckExtremumResult(fCY / (2 * fBY), rResults);
965         }
966     }
967 
transform(const basegfx::B2DHomMatrix & rMatrix)968     void B2DCubicBezier::transform(const basegfx::B2DHomMatrix& rMatrix)
969     {
970         if(rMatrix.isIdentity())
971             return;
972 
973         if(maControlPointA == maStartPoint)
974         {
975             maControlPointA = maStartPoint = rMatrix * maStartPoint;
976         }
977         else
978         {
979             maStartPoint *= rMatrix;
980             maControlPointA *= rMatrix;
981         }
982 
983         if(maControlPointB == maEndPoint)
984         {
985             maControlPointB = maEndPoint = rMatrix * maEndPoint;
986         }
987         else
988         {
989             maEndPoint *= rMatrix;
990             maControlPointB *= rMatrix;
991         }
992     }
993 
fround()994     void B2DCubicBezier::fround()
995     {
996         if(maControlPointA == maStartPoint)
997         {
998             maControlPointA = maStartPoint = basegfx::B2DPoint(
999                 std::round(maStartPoint.getX()),
1000                 std::round(maStartPoint.getY()));
1001         }
1002         else
1003         {
1004             maStartPoint = basegfx::B2DPoint(
1005                 std::round(maStartPoint.getX()),
1006                 std::round(maStartPoint.getY()));
1007             maControlPointA = basegfx::B2DPoint(
1008                 std::round(maControlPointA.getX()),
1009                 std::round(maControlPointA.getY()));
1010         }
1011 
1012         if(maControlPointB == maEndPoint)
1013         {
1014             maControlPointB = maEndPoint = basegfx::B2DPoint(
1015                 std::round(maEndPoint.getX()),
1016                 std::round(maEndPoint.getY()));
1017         }
1018         else
1019         {
1020             maEndPoint = basegfx::B2DPoint(
1021                 std::round(maEndPoint.getX()),
1022                 std::round(maEndPoint.getY()));
1023             maControlPointB = basegfx::B2DPoint(
1024                 std::round(maControlPointB.getX()),
1025                 std::round(maControlPointB.getY()));
1026         }
1027     }
1028 } // end of namespace basegfx
1029 
1030 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
1031