From 2c3b8d0c57c5b952e11831a052ab2408dfdb327f Mon Sep 17 00:00:00 2001 From: Thorsten Behrens Date: Mon, 8 Dec 2003 12:24:09 +0000 Subject: Merged to AW's latest changes, added preliminary bezier testcases --- basegfx/source/curve/b2dbeziertools.cxx | 505 +++++++++++++++++--------------- 1 file changed, 263 insertions(+), 242 deletions(-) (limited to 'basegfx/source/curve') diff --git a/basegfx/source/curve/b2dbeziertools.cxx b/basegfx/source/curve/b2dbeziertools.cxx index 67c25aacb17a..cd6a764a08ec 100644 --- a/basegfx/source/curve/b2dbeziertools.cxx +++ b/basegfx/source/curve/b2dbeziertools.cxx @@ -2,9 +2,9 @@ * * $RCSfile: b2dbeziertools.cxx,v $ * - * $Revision: 1.4 $ + * $Revision: 1.5 $ * - * last change: $Author: aw $ $Date: 2003-11-28 11:18:00 $ + * last change: $Author: thb $ $Date: 2003-12-08 13:24:06 $ * * The Contents of this file are made available subject to the terms of * either of the following licenses @@ -63,11 +63,9 @@ #include #ifndef _SOLAR_H -#include +#include // F_PI definition, no link dependency #endif -#include - #ifndef _BGFX_CURVE_B2DCUBICBEZIER_HXX #include #endif @@ -92,247 +90,270 @@ #include #endif +#include -namespace basegfx + +namespace { - namespace + class DistanceErrorFunctor { - class DistanceErrorFunctor + public: + DistanceErrorFunctor( const double& distance ) : + mfDistance2( distance*distance ), + mfLastDistanceError2( ::std::numeric_limits::max() ) { - public: - DistanceErrorFunctor( const double& distance ) : - mfDistance2( distance*distance ), - mfLastDistanceError2( ::std::numeric_limits::max() ) - { - } + } - bool subdivideFurther( const double& P1x, const double& P1y, - const double& P2x, const double& P2y, - const double& P3x, const double& P3y, - const double& P4x, const double& P4y, - const double&, const double& ) // last two values not used here - { - // Perform bezier flatness test (lecture notes from R. Schaback, - // Mathematics of Computer-Aided Design, Uni Goettingen, 2000) - // - // ||P(t) - L(t)|| <= max ||b_j - b_0 - j/n(b_n - b_0)|| - // 0<=j<=n - // - // What is calculated here is an upper bound to the distance from - // a line through b_0 and b_3 (P1 and P4 in our notation) and the - // curve. We can drop 0 and n from the running indices, since the - // argument of max becomes zero for those cases. - const double fJ1x( P2x - P1x - 1.0/3.0*(P4x - P1x) ); - const double fJ1y( P2y - P1y - 1.0/3.0*(P4y - P1y) ); - const double fJ2x( P3x - P1x - 2.0/3.0*(P4x - P1x) ); - const double fJ2y( P3y - P1y - 2.0/3.0*(P4y - P1y) ); - const double distanceError2( ::std::max( fJ1x*fJ1x + fJ1y*fJ1y, - fJ2x*fJ2x + fJ2y*fJ2y) ); - - // stop if error measure does not improve anymore. This is a - // safety guard against floating point inaccuracies. - // stop if distance from line is guaranteed to be bounded by d - bool bRet( mfLastDistanceError2 > distanceError2 && - distanceError2 >= mfDistance2 ); - - mfLastDistanceError2 = distanceError2; - - return bRet; - } + bool needsFurtherSubdivision( const double& P1x, const double& P1y, + const double& P2x, const double& P2y, + const double& P3x, const double& P3y, + const double& P4x, const double& P4y, + const double&, const double& ) // last two values not used here + { + // Perform bezier flatness test (lecture notes from R. Schaback, + // Mathematics of Computer-Aided Design, Uni Goettingen, 2000) + // + // ||P(t) - L(t)|| <= max ||b_j - b_0 - j/n(b_n - b_0)|| + // 0<=j<=n + // + // What is calculated here is an upper bound to the distance from + // a line through b_0 and b_3 (P1 and P4 in our notation) and the + // curve. We can drop 0 and n from the running indices, since the + // argument of max becomes zero for those cases. + const double fJ1x( P2x - P1x - 1.0/3.0*(P4x - P1x) ); + const double fJ1y( P2y - P1y - 1.0/3.0*(P4y - P1y) ); + const double fJ2x( P3x - P1x - 2.0/3.0*(P4x - P1x) ); + const double fJ2y( P3y - P1y - 2.0/3.0*(P4y - P1y) ); + const double distanceError2( ::std::max( fJ1x*fJ1x + fJ1y*fJ1y, + fJ2x*fJ2x + fJ2y*fJ2y) ); + + // stop if error measure does not improve anymore. This is a + // safety guard against floating point inaccuracies. + // stop if distance from line is guaranteed to be bounded by d + bool bRet( mfLastDistanceError2 > distanceError2 && + distanceError2 >= mfDistance2 ); + + mfLastDistanceError2 = distanceError2; + + return bRet; + } - private: - double mfDistance2; - double mfLastDistanceError2; - }; + private: + double mfDistance2; + double mfLastDistanceError2; + }; - class AngleErrorFunctor + class AngleErrorFunctor + { + public: + AngleErrorFunctor( const double& angleBounds ) : + // take positive branch, convert from degree to radiant + mfTanAngle( tan( fabs(angleBounds)/180.0*F_PI ) ), + mfLastTanAngle( ::std::numeric_limits::max() ) { - public: - AngleErrorFunctor( const double& angleBounds ) : - mfTanAngle( angleBounds * F_PI180 ), - mfLastTanAngle( ::std::numeric_limits::max() ) - { - } + } - bool subdivideFurther( const double P1x, const double P1y, - const double P2x, const double P2y, - const double P3x, const double P3y, - const double P4x, const double P4y, - const double Pdx, const double Pdy ) + bool needsFurtherSubdivision( const double P1x, const double P1y, + const double P2x, const double P2y, + const double P3x, const double P3y, + const double P4x, const double P4y, + const double Pdx, const double Pdy ) + { + // Test angle differences between two lines (ad + // and bd), meeting in the t=0.5 division point + // (d), and the angle from the other ends of those + // lines (b and a, resp.) to the tangents to the + // curve at this points: + // + // *__________ + // ......*b + // ... + // .. + // . + // * *d + // | . + // | . + // | . + // | . + // |. + // |. + // * + // a + // + // The idea of this test is based on the fact that + // small angle deviations between straight lines + // go by mostly unnoticed by the human eye. When + // we stop subdivision according to this test, the + // resulting line polygon is guaranteed to have + // angle differences between adjacent lines lower + // than the given threshold. + // + // When using half of the angle bound for the + // difference to the tangents at a or b, resp., + // this procedure guarantees that no angle in the + // resulting line polygon is larger than the + // specified angle bound. This is because during + // subdivision, adjacent curve segments will have + // collinear tangent vectors, thus, when each + // side's line segments differs by at most angle/2 + // from that tangent, the summed difference will + // be at most angle (this was modeled after an + // idea from Armin Weiss). + + // To stay within the notation above, a equals P1, + // the other end point of the tangent starting at + // a is P2, d is Pd, and so forth. + const ::basegfx::B2DVector vecAD( Pdx - P1x, Pdy - P1y ); + const ::basegfx::B2DVector vecDB( P4x - Pdx, P4y - Pdy ); + + const double scalarVecADDB( vecAD.scalar( vecDB ) ); + const double crossVecADDB( vecAD.cross( vecDB ) ); + + const ::basegfx::B2DVector vecStartTangent( P2x - P1x, P2y - P1y ); + const ::basegfx::B2DVector vecEndTangent( P4x - P3x, P4y - P3y ); + + const double scalarVecStartTangentAD( vecStartTangent.scalar( vecAD ) ); + const double crossVecStartTangentAD( vecStartTangent.cross( vecAD ) ); + + const double scalarVecDBEndTangent( vecDB.scalar( vecEndTangent ) ); + const double crossVecDBEndTangent( vecDB.cross( vecEndTangent ) ); + + double fCurrAngle( ::std::numeric_limits::max() ); + + // anyone has zero denominator? then we're at + // +infinity, anyway, and should better keep on + // subdividing + if( !::basegfx::fTools::equalZero( scalarVecADDB ) && + !::basegfx::fTools::equalZero( scalarVecStartTangentAD ) && + !::basegfx::fTools::equalZero( scalarVecDBEndTangent ) ) { - // Test angle differences between two lines (ad - // and bd), meeting in the t=0.5 division point - // (d), and the angle from the other ends of those - // lines (b and a, resp.) to the tangents to the - // curve at this points: - // - // *__________ - // ......*b - // ... - // .. - // . - // * *d - // | . - // | . - // | . - // | . - // |. - // |. - // * - // a - // - // When using half of the angle bound for the - // difference to the tangents at a or b, resp., - // this procedure guarantees that no angle in the - // resulting line polygon is larger than the - // specified angle bound. This is because during - // subdivision, adjacent curve segments will have - // collinear tangent vectors, thus, when each - // side's line segments differs by at most angle/2 - // from that tangent, the summed difference will - // be at most angle (this was modeled after an - // idea from Armin Weiss). - - // To stay within the notation above, a equals P1, - // the other end point of the tangent starting at - // a is P2, d is Pd, and so forth. The - const B2DVector vecAD( Pdx - P1x, Pdy - P1y ); - const B2DVector vecDB( P4x - Pdx, P4y - Pdy ); - - const double scalarVecADDB( vecAD.scalar( vecDB ) ); - const double crossVecADDB( vecAD.cross( vecDB ) ); - - const B2DVector vecStartTangent( P2x - P1x, P2y - P1y ); - const B2DVector vecEndTangent( P4x - P3x, P4y - P3y ); - - const double scalarVecStartTangentAD( vecStartTangent.scalar( vecAD ) ); - const double crossVecStartTangentAD( vecStartTangent.cross( vecAD ) ); - - const double scalarVecDBEndTangent( vecDB.scalar( vecEndTangent ) ); - const double crossVecDBEndTangent( vecDB.cross( vecEndTangent ) ); - - - double fCurrAngle( ::std::numeric_limits::max() ); - - // anyone has zero denominator? then we're at - // +infinity, anyway - if( !fTools::equalZero( scalarVecADDB ) && - !fTools::equalZero( scalarVecStartTangentAD ) && - !fTools::equalZero( scalarVecDBEndTangent ) ) + // now, sieve out quadrants which have + // deviating angles of more than 90 + // degrees. By inspection, this is everything + // with a negative scalar product. If we + // encounter such a negative scalar product, + // we can simply keep on subdividing, since at + // least one angle is then >90 degrees. + if( !::basegfx::fTools::less( scalarVecADDB, 0.0 ) && + !::basegfx::fTools::less( scalarVecStartTangentAD, 0.0 ) && + !::basegfx::fTools::less( scalarVecDBEndTangent, 0.0 ) ) { - if( scalarVecADDB > 0.0 && - scalarVecStartTangentAD > 0.0 && - scalarVecDBEndTangent > 0.0 ) - { - fCurrAngle = ::std::max( fabs( atan2( crossVecADDB, scalarVecADDB ) ), - ::std::max( fabs( atan2( crossVecStartTangentAD, scalarVecStartTangentAD ) ), - fabs( atan2( crossVecDBEndTangent, scalarVecDBEndTangent ) ) ) ); - } + // take the maximum tangens of angle + // deviation, to compare against the threshold + // below + fCurrAngle = ::std::max( fabs( crossVecADDB / scalarVecADDB ), + ::std::max( fabs( crossVecStartTangentAD / scalarVecStartTangentAD ), + fabs( crossVecDBEndTangent / scalarVecDBEndTangent ) ) ); } - - // stop if error measure does not improve anymore. This is a - // safety guard against floating point inaccuracies. - // stop if angle difference is guaranteed to be bounded by mfTanAngle - bool bRet( mfLastTanAngle > fCurrAngle && - fCurrAngle >= mfTanAngle ); - - mfLastTanAngle = fCurrAngle; - - return bRet; } - private: - double mfTanAngle; - double mfLastTanAngle; - }; - - - /* Recursively subdivide cubic bezier curve via deCasteljau. + // stop if error measure does not improve anymore. This is a + // safety guard against floating point inaccuracies. + // stop if angle difference is guaranteed to be bounded by mfTanAngle + bool bRet( mfLastTanAngle > fCurrAngle && + fCurrAngle >= mfTanAngle ); - @param rPoly - Polygon to append generated points to + mfLastTanAngle = fCurrAngle; - @param d2 - Maximal squared difference of curve to a straight line - - @param P* - Exactly four points, interpreted as support and control points of - a cubic bezier curve. - - @param old_distance2 - Last squared distance to line for this recursion - path. Used as an end condition, if it is no longer - improving. + return bRet; + } - @param recursionDepth - Depth of recursion. Used as a termination criterion, to - prevent endless looping. - */ - template < class ErrorFunctor > int ImplAdaptiveSubdivide( B2DPolygon& rPoly, - ErrorFunctor rErrorFunctor, - const double P1x, const double P1y, - const double P2x, const double P2y, - const double P3x, const double P3y, - const double P4x, const double P4y, - int recursionDepth ) + private: + double mfTanAngle; + double mfLastTanAngle; + }; + + /* Recursively subdivide cubic bezier curve via deCasteljau. + + @param rPoly + Polygon to append generated points to + + @param d2 + Maximal squared difference of curve to a straight line + + @param P* + Exactly four points, interpreted as support and control points of + a cubic bezier curve. + + @param old_distance2 + Last squared distance to line for this recursion + path. Used as an end condition, if it is no longer + improving. + + @param recursionDepth + Depth of recursion. Used as a termination criterion, to + prevent endless looping. + */ + template < class ErrorFunctor > int ImplAdaptiveSubdivide( ::basegfx::B2DPolygon& rPoly, + ErrorFunctor rErrorFunctor, + const double P1x, const double P1y, + const double P2x, const double P2y, + const double P3x, const double P3y, + const double P4x, const double P4y, + int recursionDepth ) + { + // Hard limit on recursion depth, empiric number. + enum {maxRecursionDepth=128}; + + // deCasteljau bezier arc, split at t=0.5 + // Foley/vanDam, p. 508 + + // Note that for the pure distance error method, this + // subdivision could be moved into the if-branch. But + // since this accounts for saved work only for the + // very last subdivision step, and we need the + // subdivided curve for the angle criterium, I think + // it's justified here. + const double L1x( P1x ), L1y( P1y ); + const double L2x( (P1x + P2x)*0.5 ), L2y( (P1y + P2y)*0.5 ); + const double Hx ( (P2x + P3x)*0.5 ), Hy ( (P2y + P3y)*0.5 ); + const double L3x( (L2x + Hx)*0.5 ), L3y( (L2y + Hy)*0.5 ); + const double R4x( P4x ), R4y( P4y ); + const double R3x( (P3x + P4x)*0.5 ), R3y( (P3y + P4y)*0.5 ); + const double R2x( (Hx + R3x)*0.5 ), R2y( (Hy + R3y)*0.5 ); + const double R1x( (L3x + R2x)*0.5 ), R1y( (L3y + R2y)*0.5 ); + const double L4x( R1x ), L4y( R1y ); + + // stop at recursion level 128. This is a safety guard against + // floating point inaccuracies. + if( recursionDepth < maxRecursionDepth && + rErrorFunctor.needsFurtherSubdivision( P1x, P1y, + P2x, P2y, + P3x, P3y, + P4x, P4y, + R1x, R1y ) ) { - // Hard limit on recursion depth, empiric number. - enum {maxRecursionDepth=128}; - - // deCasteljau bezier arc, split at t=0.5 - // Foley/vanDam, p. 508 - - // Note that for the pure distance error method, this - // subdivision could be moved into the if-branch. But - // since this accounts for saved work only for the - // very last subdivision step, and we need the - // subdivided curve for the angle criterium, I think - // it's justified here. - const double L1x( P1x ), L1y( P1y ); - const double L2x( (P1x + P2x)*0.5 ), L2y( (P1y + P2y)*0.5 ); - const double Hx ( (P2x + P3x)*0.5 ), Hy ( (P2y + P3y)*0.5 ); - const double L3x( (L2x + Hx)*0.5 ), L3y( (L2y + Hy)*0.5 ); - const double R4x( P4x ), R4y( P4y ); - const double R3x( (P3x + P4x)*0.5 ), R3y( (P3y + P4y)*0.5 ); - const double R2x( (Hx + R3x)*0.5 ), R2y( (Hy + R3y)*0.5 ); - const double R1x( (L3x + R2x)*0.5 ), R1y( (L3y + R2y)*0.5 ); - const double L4x( R1x ), L4y( R1y ); - - // stop at recursion level 128. This is a safety guard against - // floating point inaccuracies. - if( recursionDepth < maxRecursionDepth && - rErrorFunctor.subdivideFurther( P1x, P1y, - P2x, P2y, - P3x, P3y, - P4x, P4y, - R1x, R1y ) ) - { - // subdivide further - ++recursionDepth; + // subdivide further + ++recursionDepth; - int nGeneratedPoints(0); + int nGeneratedPoints(0); - nGeneratedPoints += ImplAdaptiveSubdivide(rPoly, rErrorFunctor, L1x, L1y, L2x, L2y, L3x, L3y, L4x, L4y, recursionDepth); - nGeneratedPoints += ImplAdaptiveSubdivide(rPoly, rErrorFunctor, R1x, R1y, R2x, R2y, R3x, R3y, R4x, R4y, recursionDepth); + nGeneratedPoints += ImplAdaptiveSubdivide(rPoly, rErrorFunctor, L1x, L1y, L2x, L2y, L3x, L3y, L4x, L4y, recursionDepth); + nGeneratedPoints += ImplAdaptiveSubdivide(rPoly, rErrorFunctor, R1x, R1y, R2x, R2y, R3x, R3y, R4x, R4y, recursionDepth); - // return number of points generated in this - // recursion branch - return nGeneratedPoints; - } - else - { - // requested resolution reached. Add end points to - // output iterator. order is preserved, since - // this is so to say depth first traversal. - rPoly.append( B2DPoint( P1x, P1y ) ); - - // return number of points generated in this - // recursion branch - return 1; - } + // return number of points generated in this + // recursion branch + return nGeneratedPoints; } + else + { + // requested resolution reached. Add end points to + // output iterator. order is preserved, since + // this is so to say depth first traversal. + rPoly.append( ::basegfx::B2DPoint( P1x, P1y ) ); + + // return number of points generated in this + // recursion branch + return 1; + } + } +} + +namespace basegfx +{ +// TODO: For bezier length calculations, use Jens Gravesens algorithm (e.g. Graphic Gems V, p. 199) // LATER #if 0 /* Approximate given cubic bezier curve by quadratic bezier segments */ @@ -482,12 +503,12 @@ namespace basegfx } } } -#endif } +#endif - sal_Int32 adaptiveSubdivideByDistance( B2DPolygon& rPoly, - const B2DCubicBezier& rCurve, - double distanceBounds ) + sal_Int32 adaptiveSubdivideByDistance( B2DPolygon& rPoly, + const B2DCubicBezier& rCurve, + double distanceBounds ) { const B2DPoint start( rCurve.getStartPoint() ); const B2DPoint control1( rCurve.getControlPointA() ); @@ -495,15 +516,15 @@ namespace basegfx const B2DPoint end( rCurve.getEndPoint() ); return ImplAdaptiveSubdivide( rPoly, - DistanceErrorFunctor( distanceBounds ), - start.getX(), start.getY(), - control1.getX(), control1.getY(), - control2.getX(), control2.getY(), - end.getX(), end.getY(), - 0 ); + DistanceErrorFunctor( distanceBounds ), + start.getX(), start.getY(), + control1.getX(), control1.getY(), + control2.getX(), control2.getY(), + end.getX(), end.getY(), + 0 ); } - sal_Int32 adaptiveSubdivideByAngle( B2DPolygon& rPoly, + sal_Int32 adaptiveSubdivideByAngle( B2DPolygon& rPoly, const B2DCubicBezier& rCurve, double angleBounds ) { @@ -513,17 +534,17 @@ namespace basegfx const B2DPoint end( rCurve.getEndPoint() ); return ImplAdaptiveSubdivide( rPoly, - AngleErrorFunctor( angleBounds ), - start.getX(), start.getY(), - control1.getX(), control1.getY(), - control2.getX(), control2.getY(), - end.getX(), end.getY(), - 0 ); + AngleErrorFunctor( angleBounds ), + start.getX(), start.getY(), + control1.getX(), control1.getY(), + control2.getX(), control2.getY(), + end.getX(), end.getY(), + 0 ); } - sal_Int32 adaptiveSubdivideByDistance( B2DPolygon& rPoly, - const B2DQuadraticBezier& rCurve, - double distanceBounds ) + sal_Int32 adaptiveSubdivideByDistance( B2DPolygon& rPoly, + const B2DQuadraticBezier& rCurve, + double distanceBounds ) { // TODO return 0; -- cgit v1.2.3