From 0fcd73fb4abd31c466b411d5951ab0a8a34afee5 Mon Sep 17 00:00:00 2001 From: Anuj Verma Date: Tue, 18 Aug 2020 10:17:46 +0530 Subject: [PATCH] [sdf] Added essential math functions. * src/sdf/ftsdf.c (cube_root, arc_cos): Added functions to compute cube root and cosine inverse of a 16.16 fixed point number. * src/sdf/ftsdf.c (solve_quadratic_equation, solve_cubic_equation): Added functions to find roots of quadratic and cubic polynomial equations. --- src/sdf/ftsdf.c | 224 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 224 insertions(+) diff --git a/src/sdf/ftsdf.c b/src/sdf/ftsdf.c index 61611f279..7fcd12a3c 100644 --- a/src/sdf/ftsdf.c +++ b/src/sdf/ftsdf.c @@ -1251,4 +1251,228 @@ return error; } + /************************************************************************** + * + * math functions + * + */ + +#if !USE_NEWTON_FOR_CONIC + + /* [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 + cube_root( FT_16D16 val ) + { + /* [IMPORTANT]: This function is not good as it may */ + /* not break, so use a lookup table instead. Or we */ + /* can use algorithm similar to `square_root'. */ + + FT_Int v, g, c; + + + if ( val == 0 || + val == -FT_INT_16D16( 1 ) || + val == FT_INT_16D16( 1 ) ) + return val; + + v = val < 0 ? -val : val; + g = square_root( v ); + c = 0; + + while ( 1 ) + { + c = FT_MulFix( FT_MulFix( g, g ), g ) - v; + c = FT_DivFix( c, 3 * FT_MulFix( g, g ) ); + + g -= c; + + if ( ( c < 0 ? -c : c ) < 30 ) + break; + } + + return val < 0 ? -g : g; + } + + /* The function calculate the perpendicular */ + /* using 1 - ( base ^ 2 ) and then use arc */ + /* tan to compute the angle. */ + static FT_16D16 + arc_cos( FT_16D16 val ) + { + FT_16D16 p, b = val; + FT_16D16 one = FT_INT_16D16( 1 ); + + + if ( b > one ) b = one; + if ( b < -one ) b = -one; + + p = one - FT_MulFix( b, b ); + p = square_root( p ); + + 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 number of real */ + /* roots of the equation. */ + /* The procedure can be found at: */ + /* https://mathworld.wolfram.com/CubicFormula.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 otherwise solve quadratic*/ + if ( a == 0 || FT_ABS( a ) < 16 ) + 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; + } + } + +#endif + /* END */