diff --git a/[GSoC]ChangeLog b/[GSoC]ChangeLog index 47509fb65..74fa792be 100644 --- a/[GSoC]ChangeLog +++ b/[GSoC]ChangeLog @@ -1,3 +1,12 @@ +2020-06-30 Anuj Verma + + [sdf] Added functions to solve polynomial equations. + + * src/sdf/ftsdf.c (solve_quadratic_equation, + solve_cubic_equation): Added functions to solve + quadratic and cubic equations. These will be used for + conic bezier curves only. + 2020-06-30 Anuj Verma * src/sdf/ftsdf.c (square_root, cube_root, arc_cos): diff --git a/src/sdf/ftsdf.c b/src/sdf/ftsdf.c index ed20df963..5691643f6 100644 --- a/src/sdf/ftsdf.c +++ b/src/sdf/ftsdf.c @@ -625,6 +625,10 @@ return q; } + /* [NOTE]: All the functions below down until rasterizer */ + /* can be avoided if we decide to subdivide the */ + /* curve into lines. */ + /* This function uses newton's iteration to find */ /* cube root of a fixed point integer. */ static FT_16D16 @@ -678,6 +682,173 @@ return FT_Atan2( b, p ); } + /* The function compute the roots of a quadratic */ + /* polynomial, assigns it to `out' and returns the */ + /* number of real roots of the equation. */ + /* The procedure can be found at: */ + /* https://mathworld.wolfram.com/QuadraticFormula.html */ + static FT_UShort + solve_quadratic_equation( FT_26D6 a, + FT_26D6 b, + FT_26D6 c, + FT_16D16 out[2] ) + { + FT_16D16 discriminant = 0; + + + a = FT_26D6_16D16( a ); + b = FT_26D6_16D16( b ); + c = FT_26D6_16D16( c ); + + if ( a == 0 ) + { + if ( b == 0 ) + { + return 0; + } + else + { + out[0] = FT_DivFix( -c, b ); + return 1; + } + } + + discriminant = FT_MulFix( b, b ) - 4 * FT_MulFix( a, c ); + + if ( discriminant < 0 ) + { + return 0; + } + else if ( discriminant == 0 ) + { + out[0] = FT_DivFix( -b, 2 * a ); + + return 1; + } + else + { + discriminant = square_root( discriminant ); + out[0] = FT_DivFix( -b + discriminant, 2 * a ); + out[1] = FT_DivFix( -b - discriminant, 2 * a ); + + return 2; + } + } + + /* The function compute the roots of a cubic polynomial */ + /* assigns it to `out' and returns the umber of real */ + /* roots of the equation. */ + /* The procedure can be found at: */ + /* https://mathworld.wolfram.com/QuadraticFormula.html */ + static FT_UShort + solve_cubic_equation( FT_26D6 a, + FT_26D6 b, + FT_26D6 c, + FT_26D6 d, + FT_16D16 out[3] ) + { + FT_16D16 q = 0; /* intermediate */ + FT_16D16 r = 0; /* intermediate */ + + FT_16D16 a2 = b; /* x^2 coefficients */ + FT_16D16 a1 = c; /* x coefficients */ + FT_16D16 a0 = d; /* constant */ + + FT_16D16 q3 = 0; + FT_16D16 r2 = 0; + FT_16D16 a23 = 0; + FT_16D16 a22 = 0; + FT_16D16 a1x2 = 0; + + + /* cutoff value for `a' to be a cubic */ + if ( a == 0 || FT_ABS( a ) < 16 ) + { + /* quadratic equation */ + return solve_quadratic_equation( b, c, d, out ); + } + if ( d == 0 ) + { + out[0] = 0; + return solve_quadratic_equation( a, b, c, out + 1 ) + 1; + } + + /* normalize the coefficients, this also makes them 16.16 */ + a2 = FT_DivFix( a2, a ); + a1 = FT_DivFix( a1, a ); + a0 = FT_DivFix( a0, a ); + + /* compute intermediates */ + a1x2 = FT_MulFix( a1, a2 ); + a22 = FT_MulFix( a2, a2 ); + a23 = FT_MulFix( a22, a2 ); + + q = ( 3 * a1 - a22 ) / 9; + r = ( 9 * a1x2 - 27 * a0 - 2 * a23 ) / 54; + + /* [BUG]: `q3' and `r2' still causes underflow. */ + + q3 = FT_MulFix( q, q ); + q3 = FT_MulFix( q3, q ); + + r2 = FT_MulFix( r, r ); + + if ( q3 < 0 && r2 < -q3 ) + { + FT_16D16 t = 0; + + + q3 = square_root( -q3 ); + t = FT_DivFix( r, q3 ); + if ( t > ( 1 << 16 ) ) t = ( 1 << 16 ); + if ( t < -( 1 << 16 ) ) t = -( 1 << 16 ); + + t = arc_cos( t ); + a2 /= 3; + q = 2 * square_root( -q ); + out[0] = FT_MulFix( q, FT_Cos( t / 3 ) ) - a2; + out[1] = FT_MulFix( q, FT_Cos( ( t + FT_ANGLE_PI * 2 ) / 3 ) ) - a2; + out[2] = FT_MulFix( q, FT_Cos( ( t + FT_ANGLE_PI * 4 ) / 3 ) ) - a2; + + return 3; + } + else if ( r2 == -q3 ) + { + FT_16D16 s = 0; + + + s = cube_root( r ); + a2 /= -3; + out[0] = a2 + ( 2 * s ); + out[1] = a2 - s; + + return 2; + } + else + { + FT_16D16 s = 0; + FT_16D16 t = 0; + FT_16D16 dis = 0; + + + if ( q3 == 0 ) + { + dis = FT_ABS( r ); + } + else + { + dis = square_root( q3 + r2 ); + } + + s = cube_root( r + dis ); + t = cube_root( r - dis ); + a2 /= -3; + out[0] = ( a2 + ( s + t ) ); + + return 1; + } + } + /*************************************************************************/ /*************************************************************************/ /** **/