xref: /trunk/main/basegfx/source/curve/b2dcubicbezier.cxx (revision 42fb6e95e2f95f39ff1eee3336bd998cf8ae20a1)
1 /**************************************************************
2  *
3  * Licensed to the Apache Software Foundation (ASF) under one
4  * or more contributor license agreements.  See the NOTICE file
5  * distributed with this work for additional information
6  * regarding copyright ownership.  The ASF licenses this file
7  * to you under the Apache License, Version 2.0 (the
8  * "License"); you may not use this file except in compliance
9  * with the License.  You may obtain a copy of the License at
10  *
11  *   http://www.apache.org/licenses/LICENSE-2.0
12  *
13  * Unless required by applicable law or agreed to in writing,
14  * software distributed under the License is distributed on an
15  * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
16  * KIND, either express or implied.  See the License for the
17  * specific language governing permissions and limitations
18  * under the License.
19  *
20  *************************************************************/
21 
22 // MARKER(update_precomp.py): autogen include statement, do not remove
23 #include "precompiled_basegfx.hxx"
24 #include <basegfx/curve/b2dcubicbezier.hxx>
25 #include <basegfx/vector/b2dvector.hxx>
26 #include <basegfx/polygon/b2dpolygon.hxx>
27 #include <basegfx/numeric/ftools.hxx>
28 
29 #include <limits>
30 
31 // #i37443#
32 #define FACTOR_FOR_UNSHARPEN    (1.6)
33 #ifdef DBG_UTIL
34 static 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) > (F_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,bool bAllowUnsharpen)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             bool bAllowUnsharpen)           // #i37443# allow the criteria to get unsharp in recursions
119         {
120             sal_uInt16 nMaxRecursionDepth(8);
121             const B2DVector aLeft(rfEA - rfPA);
122             const B2DVector aRight(rfEB - rfPB);
123             bool bLeftEqualZero(aLeft.equalZero());
124             bool bRightEqualZero(aRight.equalZero());
125             bool bAllParallel(false);
126 
127             if(bLeftEqualZero && bRightEqualZero)
128             {
129                 nMaxRecursionDepth = 0;
130             }
131             else
132             {
133                 const B2DVector aBase(rfPB - rfPA);
134                 const bool bBaseEqualZero(aBase.equalZero()); // #i72104#
135 
136                 if(!bBaseEqualZero)
137                 {
138                     const bool bLeftParallel(bLeftEqualZero ? true : areParallel(aLeft, aBase));
139                     const bool bRightParallel(bRightEqualZero ? true : areParallel(aRight, aBase));
140 
141                     if(bLeftParallel && bRightParallel)
142                     {
143                         bAllParallel = true;
144 
145                         if(!bLeftEqualZero)
146                         {
147                             double fFactor;
148 
149                             if(fabs(aBase.getX()) > fabs(aBase.getY()))
150                             {
151                                 fFactor = aLeft.getX() / aBase.getX();
152                             }
153                             else
154                             {
155                                 fFactor = aLeft.getY() / aBase.getY();
156                             }
157 
158                             if(fFactor >= 0.0 && fFactor <= 1.0)
159                             {
160                                 bLeftEqualZero = true;
161                             }
162                         }
163 
164                         if(!bRightEqualZero)
165                         {
166                             double fFactor;
167 
168                             if(fabs(aBase.getX()) > fabs(aBase.getY()))
169                             {
170                                 fFactor = aRight.getX() / -aBase.getX();
171                             }
172                             else
173                             {
174                                 fFactor = aRight.getY() / -aBase.getY();
175                             }
176 
177                             if(fFactor >= 0.0 && fFactor <= 1.0)
178                             {
179                                 bRightEqualZero = true;
180                             }
181                         }
182 
183                         if(bLeftEqualZero && bRightEqualZero)
184                         {
185                             nMaxRecursionDepth = 0;
186                         }
187                     }
188                 }
189             }
190 
191             if(nMaxRecursionDepth)
192             {
193                 // divide at 0.5 ad test both edges for angle criteria
194                 const B2DPoint aS1L(average(rfPA, rfEA));
195                 const B2DPoint aS1C(average(rfEA, rfEB));
196                 const B2DPoint aS1R(average(rfEB, rfPB));
197                 const B2DPoint aS2L(average(aS1L, aS1C));
198                 const B2DPoint aS2R(average(aS1C, aS1R));
199                 const B2DPoint aS3C(average(aS2L, aS2R));
200 
201                 // test left
202                 bool bAngleIsSmallerLeft(bAllParallel && bLeftEqualZero);
203                 if(!bAngleIsSmallerLeft)
204                 {
205                     const B2DVector aLeftLeft(bLeftEqualZero ? aS2L - aS1L : aS1L - rfPA); // #i72104#
206                     const B2DVector aRightLeft(aS2L - aS3C);
207                     const double fCurrentAngleLeft(aLeftLeft.angle(aRightLeft));
208                     bAngleIsSmallerLeft = (fabs(fCurrentAngleLeft) > (F_PI - rfAngleBound));
209                 }
210 
211                 // test right
212                 bool bAngleIsSmallerRight(bAllParallel && bRightEqualZero);
213                 if(!bAngleIsSmallerRight)
214                 {
215                     const B2DVector aLeftRight(aS2R - aS3C);
216                     const B2DVector aRightRight(bRightEqualZero ? aS2R - aS1R : aS1R - rfPB); // #i72104#
217                     const double fCurrentAngleRight(aLeftRight.angle(aRightRight));
218                     bAngleIsSmallerRight = (fabs(fCurrentAngleRight) > (F_PI - rfAngleBound));
219                 }
220 
221                 if(bAngleIsSmallerLeft && bAngleIsSmallerRight)
222                 {
223                     // no recursion necessary at all
224                     nMaxRecursionDepth = 0;
225                 }
226                 else
227                 {
228                     // left
229                     if(bAngleIsSmallerLeft)
230                     {
231                         rTarget.append(aS3C);
232                     }
233                     else
234                     {
235                         ImpSubDivAngle(rfPA, aS1L, aS2L, aS3C, rTarget, rfAngleBound, bAllowUnsharpen, nMaxRecursionDepth);
236                     }
237 
238                     // right
239                     if(bAngleIsSmallerRight)
240                     {
241                         rTarget.append(rfPB);
242                     }
243                     else
244                     {
245                         ImpSubDivAngle(aS3C, aS2R, aS1R, rfPB, rTarget, rfAngleBound, bAllowUnsharpen, nMaxRecursionDepth);
246                     }
247                 }
248             }
249 
250             if(!nMaxRecursionDepth)
251             {
252                 rTarget.append(rfPB);
253             }
254         }
255 
ImpSubDivDistance(const B2DPoint & rfPA,const B2DPoint & rfEA,const B2DPoint & rfEB,const B2DPoint & rfPB,B2DPolygon & rTarget,double fDistanceBound2,double fLastDistanceError2,sal_uInt16 nMaxRecursionDepth)256         void ImpSubDivDistance(
257             const B2DPoint& rfPA,           // start point
258             const B2DPoint& rfEA,           // edge on A
259             const B2DPoint& rfEB,           // edge on B
260             const B2DPoint& rfPB,           // end point
261             B2DPolygon& rTarget,            // target polygon
262             double fDistanceBound2,         // quadratic distance criteria
263             double fLastDistanceError2,     // the last quadratic distance error
264             sal_uInt16 nMaxRecursionDepth)  // endless loop protection
265         {
266             if(nMaxRecursionDepth)
267             {
268                 // decide if another recursion is needed. If not, set
269                 // nMaxRecursionDepth to zero
270 
271                 // Perform bezier flatness test (lecture notes from R. Schaback,
272                 // Mathematics of Computer-Aided Design, Uni Goettingen, 2000)
273                 //
274                 // ||P(t) - L(t)|| <= max     ||b_j - b_0 - j/n(b_n - b_0)||
275                 //                    0<=j<=n
276                 //
277                 // What is calculated here is an upper bound to the distance from
278                 // a line through b_0 and b_3 (rfPA and P4 in our notation) and the
279                 // curve. We can drop 0 and n from the running indices, since the
280                 // argument of max becomes zero for those cases.
281                 const double fJ1x(rfEA.getX() - rfPA.getX() - 1.0/3.0*(rfPB.getX() - rfPA.getX()));
282                 const double fJ1y(rfEA.getY() - rfPA.getY() - 1.0/3.0*(rfPB.getY() - rfPA.getY()));
283                 const double fJ2x(rfEB.getX() - rfPA.getX() - 2.0/3.0*(rfPB.getX() - rfPA.getX()));
284                 const double fJ2y(rfEB.getY() - rfPA.getY() - 2.0/3.0*(rfPB.getY() - rfPA.getY()));
285                 const double fDistanceError2(::std::max(fJ1x*fJ1x + fJ1y*fJ1y, fJ2x*fJ2x + fJ2y*fJ2y));
286 
287                 // stop if error measure does not improve anymore. This is a
288                 // safety guard against floating point inaccuracies.
289                 // stop if distance from line is guaranteed to be bounded by d
290                 const bool bFurtherDivision(fLastDistanceError2 > fDistanceError2 && fDistanceError2 >= fDistanceBound2);
291 
292                 if(bFurtherDivision)
293                 {
294                     // remember last error value
295                     fLastDistanceError2 = fDistanceError2;
296                 }
297                 else
298                 {
299                     // stop recustion
300                     nMaxRecursionDepth = 0;
301                 }
302             }
303 
304             if(nMaxRecursionDepth)
305             {
306                 // divide at 0.5
307                 const B2DPoint aS1L(average(rfPA, rfEA));
308                 const B2DPoint aS1C(average(rfEA, rfEB));
309                 const B2DPoint aS1R(average(rfEB, rfPB));
310                 const B2DPoint aS2L(average(aS1L, aS1C));
311                 const B2DPoint aS2R(average(aS1C, aS1R));
312                 const B2DPoint aS3C(average(aS2L, aS2R));
313 
314                 // left recursion
315                 ImpSubDivDistance(rfPA, aS1L, aS2L, aS3C, rTarget, fDistanceBound2, fLastDistanceError2, nMaxRecursionDepth - 1);
316 
317                 // right recursion
318                 ImpSubDivDistance(aS3C, aS2R, aS1R, rfPB, rTarget, fDistanceBound2, fLastDistanceError2, nMaxRecursionDepth - 1);
319             }
320             else
321             {
322                 rTarget.append(rfPB);
323             }
324         }
325     } // end of anonymous namespace
326 } // end of namespace basegfx
327 
328 namespace basegfx
329 {
B2DCubicBezier(const B2DCubicBezier & rBezier)330     B2DCubicBezier::B2DCubicBezier(const B2DCubicBezier& rBezier)
331     :   maStartPoint(rBezier.maStartPoint),
332         maEndPoint(rBezier.maEndPoint),
333         maControlPointA(rBezier.maControlPointA),
334         maControlPointB(rBezier.maControlPointB)
335     {
336     }
337 
B2DCubicBezier()338     B2DCubicBezier::B2DCubicBezier()
339     {
340     }
341 
B2DCubicBezier(const B2DPoint & rStart,const B2DPoint & rEnd)342     B2DCubicBezier::B2DCubicBezier(const B2DPoint& rStart, const B2DPoint& rEnd)
343     :   maStartPoint(rStart),
344         maEndPoint(rEnd),
345         maControlPointA(rStart),
346         maControlPointB(rEnd)
347     {
348     }
349 
B2DCubicBezier(const B2DPoint & rStart,const B2DPoint & rControlPointA,const B2DPoint & rControlPointB,const B2DPoint & rEnd)350     B2DCubicBezier::B2DCubicBezier(const B2DPoint& rStart, const B2DPoint& rControlPointA, const B2DPoint& rControlPointB, const B2DPoint& rEnd)
351     :   maStartPoint(rStart),
352         maEndPoint(rEnd),
353         maControlPointA(rControlPointA),
354         maControlPointB(rControlPointB)
355     {
356     }
357 
~B2DCubicBezier()358     B2DCubicBezier::~B2DCubicBezier()
359     {
360     }
361 
362     // assignment operator
operator =(const B2DCubicBezier & rBezier)363     B2DCubicBezier& B2DCubicBezier::operator=(const B2DCubicBezier& rBezier)
364     {
365         maStartPoint = rBezier.maStartPoint;
366         maEndPoint = rBezier.maEndPoint;
367         maControlPointA = rBezier.maControlPointA;
368         maControlPointB = rBezier.maControlPointB;
369 
370         return *this;
371     }
372 
373     // compare operators
operator ==(const B2DCubicBezier & rBezier) const374     bool B2DCubicBezier::operator==(const B2DCubicBezier& rBezier) const
375     {
376         return (
377             maStartPoint == rBezier.maStartPoint
378             && maEndPoint == rBezier.maEndPoint
379             && maControlPointA == rBezier.maControlPointA
380             && maControlPointB == rBezier.maControlPointB
381         );
382     }
383 
operator !=(const B2DCubicBezier & rBezier) const384     bool B2DCubicBezier::operator!=(const B2DCubicBezier& rBezier) const
385     {
386         return (
387             maStartPoint != rBezier.maStartPoint
388             || maEndPoint != rBezier.maEndPoint
389             || maControlPointA != rBezier.maControlPointA
390             || maControlPointB != rBezier.maControlPointB
391         );
392     }
393 
equal(const B2DCubicBezier & rBezier) const394     bool B2DCubicBezier::equal(const B2DCubicBezier& rBezier) const
395     {
396         return (
397             maStartPoint.equal(rBezier.maStartPoint)
398             && maEndPoint.equal(rBezier.maEndPoint)
399             && maControlPointA.equal(rBezier.maControlPointA)
400             && maControlPointB.equal(rBezier.maControlPointB)
401         );
402     }
403 
404     // test if vectors are used
isBezier() const405     bool B2DCubicBezier::isBezier() const
406     {
407         if(maControlPointA != maStartPoint || maControlPointB != maEndPoint)
408         {
409             return true;
410         }
411 
412         return false;
413     }
414 
testAndSolveTrivialBezier()415     void B2DCubicBezier::testAndSolveTrivialBezier()
416     {
417         if(maControlPointA != maStartPoint || maControlPointB != maEndPoint)
418         {
419             const B2DVector aEdge(maEndPoint - maStartPoint);
420 
421             // controls parallel to edge can be trivial. No edge -> not parallel -> control can
422             // still not be trivial (e.g. balloon loop)
423             if(!aEdge.equalZero())
424             {
425                 // get control vectors
426                 const B2DVector aVecA(maControlPointA - maStartPoint);
427                 const B2DVector aVecB(maControlPointB - maEndPoint);
428 
429                 // check if trivial per se
430                 bool bAIsTrivial(aVecA.equalZero());
431                 bool bBIsTrivial(aVecB.equalZero());
432 
433                 // #i102241# prepare inverse edge length to normalize cross values;
434                 // else the small compare value used in fTools::equalZero
435                 // will be length dependent and this detection will work as less
436                 // precise as longer the edge is. In principle, the length of the control
437                 // vector would need to be used too, but to be trivial it is assumed to
438                 // be of roughly equal length to the edge, so edge length can be used
439                 // for both. Only needed when one of both is not trivial per se.
440                 const double fInverseEdgeLength(bAIsTrivial && bBIsTrivial
441                     ? 1.0
442                     : 1.0 / aEdge.getLength());
443 
444                 // if A is not zero, check if it could be
445                 if(!bAIsTrivial)
446                 {
447                     // #i102241# parallel to edge? Check aVecA, aEdge. Use cross() which does what
448                     // we need here with the precision we need
449                     const double fCross(aVecA.cross(aEdge) * fInverseEdgeLength);
450 
451                     if(fTools::equalZero(fCross))
452                     {
453                         // get scale to edge. Use bigger distance for numeric quality
454                         const double fScale(fabs(aEdge.getX()) > fabs(aEdge.getY())
455                             ? aVecA.getX() / aEdge.getX()
456                             : aVecA.getY() / aEdge.getY());
457 
458                         // relative end point of vector in edge range?
459                         if(fTools::moreOrEqual(fScale, 0.0) && fTools::lessOrEqual(fScale, 1.0))
460                         {
461                             bAIsTrivial = true;
462                         }
463                     }
464                 }
465 
466                 // if B is not zero, check if it could be, but only if A is already trivial;
467                 // else solve to trivial will not be possible for whole edge
468                 if(bAIsTrivial && !bBIsTrivial)
469                 {
470                     // parallel to edge? Check aVecB, aEdge
471                     const double fCross(aVecB.cross(aEdge) * fInverseEdgeLength);
472 
473                     if(fTools::equalZero(fCross))
474                     {
475                         // get scale to edge. Use bigger distance for numeric quality
476                         const double fScale(fabs(aEdge.getX()) > fabs(aEdge.getY())
477                             ? aVecB.getX() / aEdge.getX()
478                             : aVecB.getY() / aEdge.getY());
479 
480                         // end point of vector in edge range? Caution: controlB is directed AGAINST edge
481                         if(fTools::lessOrEqual(fScale, 0.0) && fTools::moreOrEqual(fScale, -1.0))
482                         {
483                             bBIsTrivial = true;
484                         }
485                     }
486                 }
487 
488                 // if both are/can be reduced, do it.
489                 // Not possible if only one is/can be reduced (!)
490                 if(bAIsTrivial && bBIsTrivial)
491                 {
492                     maControlPointA = maStartPoint;
493                     maControlPointB = maEndPoint;
494                 }
495             }
496         }
497     }
498 
499     namespace {
impGetLength(const B2DCubicBezier & rEdge,double fDeviation,sal_uInt32 nRecursionWatch)500         double impGetLength(const B2DCubicBezier& rEdge, double fDeviation, sal_uInt32 nRecursionWatch)
501         {
502             const double fEdgeLength(rEdge.getEdgeLength());
503             const double fControlPolygonLength(rEdge.getControlPolygonLength());
504             const double fCurrentDeviation(fTools::equalZero(fControlPolygonLength) ? 0.0 : 1.0 - (fEdgeLength / fControlPolygonLength));
505 
506             if(!nRecursionWatch || fTools:: lessOrEqual(fCurrentDeviation, fDeviation))
507             {
508                 return (fEdgeLength + fControlPolygonLength) * 0.5;
509             }
510             else
511             {
512                 B2DCubicBezier aLeft, aRight;
513                 const double fNewDeviation(fDeviation * 0.5);
514                 const sal_uInt32 nNewRecursionWatch(nRecursionWatch - 1);
515 
516                 rEdge.split(0.5, &aLeft, &aRight);
517 
518                 return impGetLength(aLeft, fNewDeviation, nNewRecursionWatch)
519                     + impGetLength(aRight, fNewDeviation, nNewRecursionWatch);
520             }
521         }
522     }
523 
getLength(double fDeviation) const524     double B2DCubicBezier::getLength(double fDeviation) const
525     {
526         if(isBezier())
527         {
528             if(fDeviation < 0.00000001)
529             {
530                 fDeviation = 0.00000001;
531             }
532 
533             return impGetLength(*this, fDeviation, 6);
534         }
535         else
536         {
537             return B2DVector(getEndPoint() - getStartPoint()).getLength();
538         }
539     }
540 
getEdgeLength() const541     double B2DCubicBezier::getEdgeLength() const
542     {
543         const B2DVector aEdge(maEndPoint - maStartPoint);
544         return aEdge.getLength();
545     }
546 
getControlPolygonLength() const547     double B2DCubicBezier::getControlPolygonLength() const
548     {
549         const B2DVector aVectorA(maControlPointA - maStartPoint);
550         const B2DVector aVectorB(maEndPoint - maControlPointB);
551 
552         if(!aVectorA.equalZero() || !aVectorB.equalZero())
553         {
554             const B2DVector aTop(maControlPointB - maControlPointA);
555             return (aVectorA.getLength() + aVectorB.getLength() + aTop.getLength());
556         }
557         else
558         {
559             return getEdgeLength();
560         }
561     }
562 
adaptiveSubdivideByAngle(B2DPolygon & rTarget,double fAngleBound,bool bAllowUnsharpen) const563     void B2DCubicBezier::adaptiveSubdivideByAngle(B2DPolygon& rTarget, double fAngleBound, bool bAllowUnsharpen) const
564     {
565         if(isBezier())
566         {
567             // use support method #i37443# and allow unsharpen the criteria
568             ImpSubDivAngleStart(maStartPoint, maControlPointA, maControlPointB, maEndPoint, rTarget, fAngleBound * F_PI180, bAllowUnsharpen);
569         }
570         else
571         {
572             rTarget.append(getEndPoint());
573         }
574     }
575 
getTangent(double t) const576     B2DVector B2DCubicBezier::getTangent(double t) const
577     {
578         if(fTools::lessOrEqual(t, 0.0))
579         {
580             // tangent in start point
581             B2DVector aTangent(getControlPointA() - getStartPoint());
582 
583             if(!aTangent.equalZero())
584             {
585                 return aTangent;
586             }
587 
588             // start point and control vector are the same, fallback
589             // to implicit start vector to control point B
590             aTangent = (getControlPointB() - getStartPoint()) * 0.3;
591 
592             if(!aTangent.equalZero())
593             {
594                 return aTangent;
595             }
596 
597             // not a bezier at all, return edge vector
598             return (getEndPoint() - getStartPoint()) * 0.3;
599         }
600         else if(fTools::moreOrEqual(t, 1.0))
601         {
602             // tangent in end point
603             B2DVector aTangent(getEndPoint() - getControlPointB());
604 
605             if(!aTangent.equalZero())
606             {
607                 return aTangent;
608             }
609 
610             // end point and control vector are the same, fallback
611             // to implicit start vector from control point A
612             aTangent = (getEndPoint() - getControlPointA()) * 0.3;
613 
614             if(!aTangent.equalZero())
615             {
616                 return aTangent;
617             }
618 
619             // not a bezier at all, return edge vector
620             return (getEndPoint() - getStartPoint()) * 0.3;
621         }
622         else
623         {
624             // t is in ]0.0 .. 1.0[. Split and extract
625             B2DCubicBezier aRight;
626             split(t, 0, &aRight);
627 
628             return aRight.getControlPointA() - aRight.getStartPoint();
629         }
630     }
631 
632     // #i37443# adaptive subdivide by nCount subdivisions
adaptiveSubdivideByCount(B2DPolygon & rTarget,sal_uInt32 nCount) const633     void B2DCubicBezier::adaptiveSubdivideByCount(B2DPolygon& rTarget, sal_uInt32 nCount) const
634     {
635         const double fLenFact(1.0 / static_cast< double >(nCount + 1));
636 
637         for(sal_uInt32 a(1); a <= nCount; a++)
638         {
639             const double fPos(static_cast< double >(a) * fLenFact);
640             rTarget.append(interpolatePoint(fPos));
641         }
642 
643         rTarget.append(getEndPoint());
644     }
645 
646     // adaptive subdivide by distance
adaptiveSubdivideByDistance(B2DPolygon & rTarget,double fDistanceBound) const647     void B2DCubicBezier::adaptiveSubdivideByDistance(B2DPolygon& rTarget, double fDistanceBound) const
648     {
649         if(isBezier())
650         {
651             ImpSubDivDistance(maStartPoint, maControlPointA, maControlPointB, maEndPoint, rTarget,
652                 fDistanceBound * fDistanceBound, ::std::numeric_limits<double>::max(), 30);
653         }
654         else
655         {
656             rTarget.append(getEndPoint());
657         }
658     }
659 
interpolatePoint(double t) const660     B2DPoint B2DCubicBezier::interpolatePoint(double t) const
661     {
662         OSL_ENSURE(t >= 0.0 && t <= 1.0, "B2DCubicBezier::interpolatePoint: Access out of range (!)");
663 
664         if(isBezier())
665         {
666             const B2DPoint aS1L(interpolate(maStartPoint, maControlPointA, t));
667             const B2DPoint aS1C(interpolate(maControlPointA, maControlPointB, t));
668             const B2DPoint aS1R(interpolate(maControlPointB, maEndPoint, t));
669             const B2DPoint aS2L(interpolate(aS1L, aS1C, t));
670             const B2DPoint aS2R(interpolate(aS1C, aS1R, t));
671 
672             return interpolate(aS2L, aS2R, t);
673         }
674         else
675         {
676             return interpolate(maStartPoint, maEndPoint, t);
677         }
678     }
679 
getSmallestDistancePointToBezierSegment(const B2DPoint & rTestPoint,double & rCut) const680     double B2DCubicBezier::getSmallestDistancePointToBezierSegment(const B2DPoint& rTestPoint, double& rCut) const
681     {
682         const sal_uInt32 nInitialDivisions(3L);
683         B2DPolygon aInitialPolygon;
684 
685         // as start make a fix division, creates nInitialDivisions + 2L points
686         aInitialPolygon.append(getStartPoint());
687         adaptiveSubdivideByCount(aInitialPolygon, nInitialDivisions);
688 
689         // now look for the closest point
690         const sal_uInt32 nPointCount(aInitialPolygon.count());
691         B2DVector aVector(rTestPoint - aInitialPolygon.getB2DPoint(0L));
692         double fQuadDist(aVector.getX() * aVector.getX() + aVector.getY() * aVector.getY());
693         double fNewQuadDist;
694         sal_uInt32 nSmallestIndex(0L);
695 
696         for(sal_uInt32 a(1L); a < nPointCount; a++)
697         {
698             aVector = B2DVector(rTestPoint - aInitialPolygon.getB2DPoint(a));
699             fNewQuadDist = aVector.getX() * aVector.getX() + aVector.getY() * aVector.getY();
700 
701             if(fNewQuadDist < fQuadDist)
702             {
703                 fQuadDist = fNewQuadDist;
704                 nSmallestIndex = a;
705             }
706         }
707 
708         // look right and left for even smaller distances
709         double fStepValue(1.0 / (double)((nPointCount - 1L) * 2L)); // half the edge step width
710         double fPosition((double)nSmallestIndex / (double)(nPointCount - 1L));
711         bool bDone(false);
712 
713         while(!bDone)
714         {
715             if(!bDone)
716             {
717                 // test left
718                 double fPosLeft(fPosition - fStepValue);
719 
720                 if(fPosLeft < 0.0)
721                 {
722                     fPosLeft = 0.0;
723                     aVector = B2DVector(rTestPoint - maStartPoint);
724                 }
725                 else
726                 {
727                     aVector = B2DVector(rTestPoint - interpolatePoint(fPosLeft));
728                 }
729 
730                 fNewQuadDist = aVector.getX() * aVector.getX() + aVector.getY() * aVector.getY();
731 
732                 if(fTools::less(fNewQuadDist, fQuadDist))
733                 {
734                     fQuadDist = fNewQuadDist;
735                     fPosition = fPosLeft;
736                 }
737                 else
738                 {
739                     // test right
740                     double fPosRight(fPosition + fStepValue);
741 
742                     if(fPosRight > 1.0)
743                     {
744                         fPosRight = 1.0;
745                         aVector = B2DVector(rTestPoint - maEndPoint);
746                     }
747                     else
748                     {
749                         aVector = B2DVector(rTestPoint - interpolatePoint(fPosRight));
750                     }
751 
752                     fNewQuadDist = aVector.getX() * aVector.getX() + aVector.getY() * aVector.getY();
753 
754                     if(fTools::less(fNewQuadDist, fQuadDist))
755                     {
756                         fQuadDist = fNewQuadDist;
757                         fPosition = fPosRight;
758                     }
759                     else
760                     {
761                         // not less left or right, done
762                         bDone = true;
763                     }
764                 }
765             }
766 
767             if(0.0 == fPosition || 1.0 == fPosition)
768             {
769                 // if we are completely left or right, we are done
770                 bDone = true;
771             }
772 
773             if(!bDone)
774             {
775                 // prepare next step value
776                 fStepValue /= 2.0;
777             }
778         }
779 
780         rCut = fPosition;
781         return sqrt(fQuadDist);
782     }
783 
split(double t,B2DCubicBezier * pBezierA,B2DCubicBezier * pBezierB) const784     void B2DCubicBezier::split(double t, B2DCubicBezier* pBezierA, B2DCubicBezier* pBezierB) const
785     {
786         OSL_ENSURE(t >= 0.0 && t <= 1.0, "B2DCubicBezier::split: Access out of range (!)");
787 
788         if(!pBezierA && !pBezierB)
789         {
790             return;
791         }
792 
793         if(isBezier())
794         {
795             const B2DPoint aS1L(interpolate(maStartPoint, maControlPointA, t));
796             const B2DPoint aS1C(interpolate(maControlPointA, maControlPointB, t));
797             const B2DPoint aS1R(interpolate(maControlPointB, maEndPoint, t));
798             const B2DPoint aS2L(interpolate(aS1L, aS1C, t));
799             const B2DPoint aS2R(interpolate(aS1C, aS1R, t));
800             const B2DPoint aS3C(interpolate(aS2L, aS2R, t));
801 
802             if(pBezierA)
803             {
804                 pBezierA->setStartPoint(maStartPoint);
805                 pBezierA->setEndPoint(aS3C);
806                 pBezierA->setControlPointA(aS1L);
807                 pBezierA->setControlPointB(aS2L);
808             }
809 
810             if(pBezierB)
811             {
812                 pBezierB->setStartPoint(aS3C);
813                 pBezierB->setEndPoint(maEndPoint);
814                 pBezierB->setControlPointA(aS2R);
815                 pBezierB->setControlPointB(aS1R);
816             }
817         }
818         else
819         {
820             const B2DPoint aSplit(interpolate(maStartPoint, maEndPoint, t));
821 
822             if(pBezierA)
823             {
824                 pBezierA->setStartPoint(maStartPoint);
825                 pBezierA->setEndPoint(aSplit);
826                 pBezierA->setControlPointA(maStartPoint);
827                 pBezierA->setControlPointB(aSplit);
828             }
829 
830             if(pBezierB)
831             {
832                 pBezierB->setStartPoint(aSplit);
833                 pBezierB->setEndPoint(maEndPoint);
834                 pBezierB->setControlPointA(aSplit);
835                 pBezierB->setControlPointB(maEndPoint);
836             }
837         }
838     }
839 
snippet(double fStart,double fEnd) const840     B2DCubicBezier B2DCubicBezier::snippet(double fStart, double fEnd) const
841     {
842         B2DCubicBezier aRetval;
843 
844         if(fTools::more(fStart, 1.0))
845         {
846             fStart = 1.0;
847         }
848         else if(fTools::less(fStart, 0.0))
849         {
850             fStart = 0.0;
851         }
852 
853         if(fTools::more(fEnd, 1.0))
854         {
855             fEnd = 1.0;
856         }
857         else if(fTools::less(fEnd, 0.0))
858         {
859             fEnd = 0.0;
860         }
861 
862         if(fEnd <= fStart)
863         {
864             // empty or NULL, create single point at center
865             const double fSplit((fEnd + fStart) * 0.5);
866             const B2DPoint aPoint(interpolate(getStartPoint(), getEndPoint(), fSplit));
867             aRetval.setStartPoint(aPoint);
868             aRetval.setEndPoint(aPoint);
869             aRetval.setControlPointA(aPoint);
870             aRetval.setControlPointB(aPoint);
871         }
872         else
873         {
874             if(isBezier())
875             {
876                 // copy bezier; cut off right, then cut off left. Do not forget to
877                 // adapt cut value when both cuts happen
878                 const bool bEndIsOne(fTools::equal(fEnd, 1.0));
879                 const bool bStartIsZero(fTools::equalZero(fStart));
880                 aRetval = *this;
881 
882                 if(!bEndIsOne)
883                 {
884                     aRetval.split(fEnd, &aRetval, 0);
885 
886                     if(!bStartIsZero)
887                     {
888                         fStart /= fEnd;
889                     }
890                 }
891 
892                 if(!bStartIsZero)
893                 {
894                     aRetval.split(fStart, 0, &aRetval);
895                 }
896             }
897             else
898             {
899                 // no bezier, create simple edge
900                 const B2DPoint aPointA(interpolate(getStartPoint(), getEndPoint(), fStart));
901                 const B2DPoint aPointB(interpolate(getStartPoint(), getEndPoint(), fEnd));
902                 aRetval.setStartPoint(aPointA);
903                 aRetval.setEndPoint(aPointB);
904                 aRetval.setControlPointA(aPointA);
905                 aRetval.setControlPointB(aPointB);
906             }
907         }
908 
909         return aRetval;
910     }
911 
getRange() const912     B2DRange B2DCubicBezier::getRange() const
913     {
914         B2DRange aRetval(maStartPoint, maEndPoint);
915 
916         aRetval.expand(maControlPointA);
917         aRetval.expand(maControlPointB);
918 
919         return aRetval;
920     }
921 
getMinimumExtremumPosition(double & rfResult) const922     bool B2DCubicBezier::getMinimumExtremumPosition(double& rfResult) const
923     {
924         ::std::vector< double > aAllResults;
925 
926         aAllResults.reserve(4);
927         getAllExtremumPositions(aAllResults);
928 
929         const sal_uInt32 nCount(aAllResults.size());
930 
931         if(!nCount)
932         {
933             return false;
934         }
935         else if(1 == nCount)
936         {
937             rfResult = aAllResults[0];
938             return true;
939         }
940         else
941         {
942             rfResult = *(::std::min_element(aAllResults.begin(), aAllResults.end()));
943             return true;
944         }
945     }
946 
947     namespace
948     {
impCheckExtremumResult(double fCandidate,::std::vector<double> & rResult)949         inline void impCheckExtremumResult(double fCandidate, ::std::vector< double >& rResult)
950         {
951             // check for range ]0.0 .. 1.0[ with excluding 1.0 and 0.0 clearly
952             // by using the equalZero test, NOT ::more or ::less which will use the
953             // ApproxEqual() which is too exact here
954             if(fCandidate > 0.0 && !fTools::equalZero(fCandidate))
955             {
956                 if(fCandidate < 1.0 && !fTools::equalZero(fCandidate - 1.0))
957                 {
958                     rResult.push_back(fCandidate);
959                 }
960             }
961         }
962     }
963 
getAllExtremumPositions(::std::vector<double> & rResults) const964     void B2DCubicBezier::getAllExtremumPositions(::std::vector< double >& rResults) const
965     {
966         rResults.clear();
967 
968         // calculate the x-extrema parameters by zeroing first x-derivative
969         // of the cubic bezier's parametric formula, which results in a
970         // quadratic equation: dBezier/dt = t*t*fAX - 2*t*fBX + fCX
971         const B2DPoint aControlDiff( maControlPointA - maControlPointB );
972         double fCX = maControlPointA.getX() - maStartPoint.getX();
973         const double fBX = fCX + aControlDiff.getX();
974         const double fAX = 3 * aControlDiff.getX() + (maEndPoint.getX() - maStartPoint.getX());
975 
976         if(fTools::equalZero(fCX))
977         {
978             // detect fCX equal zero and truncate to real zero value in that case
979             fCX = 0.0;
980         }
981 
982         if( !fTools::equalZero(fAX) )
983         {
984             // derivative is polynomial of order 2 => use binomial formula
985             const double fD = fBX*fBX - fAX*fCX;
986             if( fD >= 0.0 )
987             {
988                 const double fS = sqrt(fD);
989                 // calculate both roots (avoiding a numerically unstable subtraction)
990                 const double fQ = fBX + ((fBX >= 0) ? +fS : -fS);
991                 impCheckExtremumResult(fQ / fAX, rResults);
992                 if( !fTools::equalZero(fS) ) // ignore root multiplicity
993                     impCheckExtremumResult(fCX / fQ, rResults);
994             }
995         }
996         else if( !fTools::equalZero(fBX) )
997         {
998             // derivative is polynomial of order 1 => one extrema
999             impCheckExtremumResult(fCX / (2 * fBX), rResults);
1000         }
1001 
1002         // calculate the y-extrema parameters by zeroing first y-derivative
1003         double fCY = maControlPointA.getY() - maStartPoint.getY();
1004         const double fBY = fCY + aControlDiff.getY();
1005         const double fAY = 3 * aControlDiff.getY() + (maEndPoint.getY() - maStartPoint.getY());
1006 
1007         if(fTools::equalZero(fCY))
1008         {
1009             // detect fCY equal zero and truncate to real zero value in that case
1010             fCY = 0.0;
1011         }
1012 
1013         if( !fTools::equalZero(fAY) )
1014         {
1015             // derivative is polynomial of order 2 => use binomial formula
1016             const double fD = fBY*fBY - fAY*fCY;
1017             if( fD >= 0.0 )
1018             {
1019                 const double fS = sqrt(fD);
1020                 // calculate both roots (avoiding a numerically unstable subtraction)
1021                 const double fQ = fBY + ((fBY >= 0) ? +fS : -fS);
1022                 impCheckExtremumResult(fQ / fAY, rResults);
1023                 if( !fTools::equalZero(fS) ) // ignore root multiplicity
1024                     impCheckExtremumResult(fCY / fQ, rResults);
1025             }
1026         }
1027         else if( !fTools::equalZero(fBY) )
1028         {
1029             // derivative is polynomial of order 1 => one extrema
1030             impCheckExtremumResult(fCY / (2 * fBY), rResults);
1031         }
1032     }
1033 
getMaxDistancePositions(double pResult[2]) const1034     int B2DCubicBezier::getMaxDistancePositions( double pResult[2]) const
1035     {
1036         // the distance from the bezier to a line through start and end
1037         // is proportional to (ENDx-STARTx,ENDy-STARTy)*(+BEZIERy(t)-STARTy,-BEZIERx(t)-STARTx)
1038         // this distance becomes zero for at least t==0 and t==1
1039         // its extrema that are between 0..1 are interesting as split candidates
1040         // its derived function has the form dD/dt = fA*t^2 + 2*fB*t + fC
1041         const B2DPoint aRelativeEndPoint(maEndPoint-maStartPoint);
1042         const double fA = (3 * (maControlPointA.getX() - maControlPointB.getX()) + aRelativeEndPoint.getX()) * aRelativeEndPoint.getY()
1043                 - (3 * (maControlPointA.getY() - maControlPointB.getY()) + aRelativeEndPoint.getY()) * aRelativeEndPoint.getX();
1044         const double fB = (maControlPointB.getX() - 2 * maControlPointA.getX() + maStartPoint.getX()) * aRelativeEndPoint.getY()
1045                 - (maControlPointB.getY() - 2 * maControlPointA.getY() + maStartPoint.getY()) * aRelativeEndPoint.getX();
1046         const double fC = (maControlPointA.getX() - maStartPoint.getX()) * aRelativeEndPoint.getY()
1047                 - (maControlPointA.getY() - maStartPoint.getY()) * aRelativeEndPoint.getX();
1048 
1049         // test for degenerated case: order<2
1050         if( fTools::equalZero(fA) )
1051         {
1052             // test for degenerated case: order==0
1053             if( fTools::equalZero(fB) )
1054                 return 0;
1055 
1056             // solving the order==1 polynomial is trivial
1057             pResult[0] = -fC / (2*fB);
1058 
1059             // test root and ignore it when it is outside the curve
1060             int nCount = ((pResult[0] > 0) && (pResult[0] < 1));
1061             return nCount;
1062         }
1063 
1064         // derivative is polynomial of order 2
1065         // check if the polynomial has non-imaginary roots
1066         const double fD = fB*fB - fA*fC;
1067         if( fD >= 0.0 ) // TODO: is this test needed? geometrically not IMHO
1068         {
1069             // calculate first root (avoiding a numerically unstable subtraction)
1070             const double fS = sqrt(fD);
1071             const double fQ = -(fB + ((fB >= 0) ? +fS : -fS));
1072             pResult[0] = fQ / fA;
1073             // ignore root when it is outside the curve
1074             static const double fEps = 1e-9;
1075             int nCount = ((pResult[0] > fEps) && (pResult[0] < fEps));
1076 
1077             // ignore root multiplicity
1078             if( !fTools::equalZero(fD) )
1079             {
1080                 // calculate the other root
1081                 const double fRoot = fC / fQ;
1082                 // ignore root when it is outside the curve
1083                 if( (fRoot > fEps) && (fRoot < 1.0-fEps) )
1084                     pResult[ nCount++ ] = fRoot;
1085             }
1086 
1087             return nCount;
1088         }
1089 
1090         return 0;
1091     }
1092 
1093 } // end of namespace basegfx
1094 
1095 /* vim: set noet sw=4 ts=4: */
1096