From 0e6d67f63666c27a0fed7a741eb2653e20f84dc0 Mon Sep 17 00:00:00 2001 From: martinfouilleul Date: Mon, 3 Jul 2023 14:21:53 +0200 Subject: [PATCH] [wip, win32, canvas] Cubics segment setup --- examples/polygon/main.c | 10 +- src/glsl_shaders/common.glsl | 1 - src/glsl_shaders/segment_setup.glsl | 492 ++++++++++++++++++++++++++++ 3 files changed, 497 insertions(+), 6 deletions(-) diff --git a/examples/polygon/main.c b/examples/polygon/main.c index 35572b4..9401206 100644 --- a/examples/polygon/main.c +++ b/examples/polygon/main.c @@ -128,13 +128,13 @@ int main() mg_close_path(); mg_set_color_rgba(0, 0, 1, 1); mg_fill(); -/* - mg_move_to(2*400, 2*400); - mg_cubic_to(2*400, 2*200, 2*600, 2*500, 2*600, 2*400); + + mg_move_to(200, 450); + mg_cubic_to(200, 250, 400, 550, 400, 450); mg_close_path(); - mg_set_color_rgba(0, 0, 1, 1); + mg_set_color_rgba(1, 0.5, 0, 1); mg_fill(); -*/ + /* mg_set_joint(MG_JOINT_NONE); mg_set_max_joint_excursion(20); diff --git a/src/glsl_shaders/common.glsl b/src/glsl_shaders/common.glsl index a679600..914d894 100644 --- a/src/glsl_shaders/common.glsl +++ b/src/glsl_shaders/common.glsl @@ -44,7 +44,6 @@ struct mg_gl_segment int config; //TODO pack these int windingIncrement; vec4 box; - mat3 hullMatrix; mat3 implicitMatrix; float sign; vec2 hullVertex; diff --git a/src/glsl_shaders/segment_setup.glsl b/src/glsl_shaders/segment_setup.glsl index 90d3743..9fb505f 100644 --- a/src/glsl_shaders/segment_setup.glsl +++ b/src/glsl_shaders/segment_setup.glsl @@ -334,6 +334,492 @@ void quadratic_setup(vec2 p[4], int pathIndex) } } +int quadratic_roots_with_det(float a, float b, float c, float det, out float r[2]) +{ + int count = 0; + + if(a == 0) + { + if(b) + { + count = 1; + r[0] = -c/b; + } + } + else + { + b /= 2.0; + + if(det >= 0) + { + count = (det == 0) ? 1 : 2; + + if(b > 0) + { + float q = b + sqrt(det); + r[0] = -c/q; + r[1] = -q/a; + } + else if(b < 0) + { + float q = -b + sqrt(det); + r[0] = q/a; + r[1] = c/q; + } + else + { + float q = sqrt(-a*c); + if(abs(a) >= abs(c)) + { + r[0] = q/a; + r[1] = -q/a; + } + else + { + r[0] = -c/q; + r[1] = c/q; + } + } + } + } + if(count>1 && r[0] > r[1]) + { + float tmp = r[0]; + r[0] = r[1]; + r[1] = tmp; + } + return(count); +} + +int quadratic_roots(float a, float b, float c, out float r[2]) +{ + float det = square(b)/4. - a*c; + return(quadratic_roots_with_det(a, b, c, det, r)); +} + +vec2 cubic_blossom(vec2 p[4], float u, float v, float w) +{ + vec2 b10 = u*p[1] + (1-u)*p[0]; + vec2 b11 = u*p[2] + (1-u)*p[1]; + vec2 b12 = u*p[3] + (1-u)*p[2]; + vec2 b20 = v*b11 + (1-v)*b10; + vec2 b21 = v*b12 + (1-v)*b11; + vec2 b30 = w*b21 + (1-w)*b20; + return(b30); +} + +void cubic_slice(vec2 p[4], float s0, float s1, out vec2 sp[4]) +{ + /*NOTE: using blossoms to compute sub-curve control points ensure that the fourth point + of sub-curve (s0, s1) and the first point of sub-curve (s1, s3) match. + However, due to numerical errors, the evaluation of B(s=0) might not be equal to + p[0] (and likewise, B(s=1) might not equal p[3]). + We handle that case explicitly to ensure that we don't create gaps in the paths. + */ + sp[0] = (s0 == 0) ? p[0] : cubic_blossom(p, s0, s0, s0); + sp[1] = cubic_blossom(p, s0, s0, s1); + sp[2] = cubic_blossom(p, s0, s1, s1); + sp[3] = (s1 == 1) ? p[3] : cubic_blossom(p, s1, s1, s1); +} + +#define CUBIC_ERROR 0 +#define CUBIC_SERPENTINE 1 +#define CUBIC_CUSP 2 +#define CUBIC_CUSP_INFINITY 3 +#define CUBIC_LOOP 4 +#define CUBIC_DEGENERATE_QUADRATIC 5 +#define CUBIC_DEGENERATE_LINE 6 + +struct cubic_info +{ + int kind; + mat4 K; + vec2 ts[2]; + float d1; + float d2; + float d3; +}; + +cubic_info cubic_classify(vec2 c[4]) +{ + cubic_info result; + result.kind = CUBIC_ERROR; + mat4 F; + + /*NOTE(martin): + now, compute determinants d0, d1, d2, d3, which gives the coefficients of the + inflection points polynomial: + + I(t, s) = d0*t^3 - 3*d1*t^2*s + 3*d2*t*s^2 - d3*s^3 + + The roots of this polynomial are the inflection points of the parametric curve, in homogeneous + coordinates (ie we can have an inflection point at inifinity with s=0). + + |x3 y3 w3| |x3 y3 w3| |x3 y3 w3| |x2 y2 w2| + d0 = det |x2 y2 w2| d1 = -det |x2 y2 w2| d2 = det |x1 y1 w1| d3 = -det |x1 y1 w1| + |x1 y1 w1| |x0 y0 w0| |x0 y0 w0| |x0 y0 w0| + + In our case, the pi.w equal 1 (no point at infinity), so _in_the_power_basis_, w1 = w2 = w3 = 0 and w0 = 1 + (which also means d0 = 0) + + //WARN: there seems to be a mismatch between the signs of the d_i and the orientation test in the Loop-Blinn paper? + // flipping the sign of the d_i doesn't change the roots (and the implicit matrix), but it does change the orientation. + // Keeping the signs of the paper puts the interior on the left of parametric travel, unlike what's stated in the paper. + // this may very well be an error on my part that's cancelled by flipping the signs of the d_i though! + */ + + float d1 = -(c[3].y*c[2].x - c[3].x*c[2].y); + float d2 = -(c[3].x*c[1].y - c[3].y*c[1].x); + float d3 = -(c[2].y*c[1].x - c[2].x*c[1].y); + + result.d1 = d1; + result.d2 = d2; + result.d3 = d3; + + //NOTE(martin): compute the second factor of the discriminant discr(I) = d1^2*(3*d2^2 - 4*d3*d1) + float discrFactor2 = 3.0*square(d2) - 4.0*d3*d1; + + //NOTE(martin): each following case gives the number of roots, hence the category of the parametric curve + if(abs(d1) <= 1e-6 && abs(d2) <= 1e-6 && abs(d3) > 1e-6) + { + //NOTE(martin): quadratic degenerate case + //NOTE(martin): compute quadratic curve control point, which is at p0 + 1.5*(p1-p0) = 1.5*p1 - 0.5*p0 + result.kind = CUBIC_DEGENERATE_QUADRATIC; + } + else if( (discrFactor2 > 0 && abs(d1) > 1e-6) + ||(discrFactor2 == 0 && abs(d1) > 1e-6)) + { + //NOTE(martin): serpentine curve or cusp with inflection at infinity + // (these two cases are handled the same way). + //NOTE(martin): compute the solutions (tl, sl), (tm, sm), and (tn, sn) of the inflection point equation + float tmtl[2]; + quadratic_roots_with_det(1, -2*d2, (4./3.*d1*d3), (1./3.)*discrFactor2, tmtl); + + float tm = tmtl[0]; + float sm = 2*d1; + float tl = tmtl[1]; + float sl = 2*d1; + + float invNorm = 1/sqrt(square(tm) + square(sm)); + tm *= invNorm; + sm *= invNorm; + + invNorm = 1/sqrt(square(tl) + square(sl)); + tl *= invNorm; + sl *= invNorm; + + /*NOTE(martin): + the power basis coefficients of points k,l,m,n are collected into the rows of the 4x4 matrix F: + + | tl*tm tl^3 tm^3 1 | + | -sm*tl - sl*tm -3sl*tl^2 -3*sm*tm^2 0 | + | sl*sm 3*sl^2*tl 3*sm^2*tm 0 | + | 0 -sl^3 -sm^3 0 | + */ + result.kind = (discrFactor2 > 0 && d1 != 0) ? CUBIC_SERPENTINE : CUBIC_CUSP; + + F = mat4(tl*tm, -sm*tl-sl*tm, sl*sm, 0, + cube(tl), -3*sl*square(tl), 3*square(sl)*tl, -cube(sl), + cube(tm), -3*sm*square(tm), 3*square(sm)*tm, -cube(sm), + 1, 0, 0, 0); + + result.ts[0] = vec2(tm, sm); + result.ts[1] = vec2(tl, sl); + } + else if(discrFactor2 < 0 && abs(d1) > 1e-6) + { + //NOTE(martin): loop curve + result.kind = CUBIC_LOOP; + + float tetd[2]; + quadratic_roots_with_det(1, -2*d2, 4*(square(d2)-d1*d3), -discrFactor2, tetd); + + float td = tetd[1]; + float sd = 2*d1; + float te = tetd[0]; + float se = 2*d1; + + float invNorm = 1/sqrt(square(td) + square(sd)); + td *= invNorm; + sd *= invNorm; + + invNorm = 1/sqrt(square(te) + square(se)); + te *= invNorm; + se *= invNorm; + + /*NOTE(martin): + the power basis coefficients of points k,l,m,n are collected into the rows of the 4x4 matrix F: + + | td*te td^2*te td*te^2 1 | + | -se*td - sd*te -se*td^2 - 2sd*te*td -sd*te^2 - 2*se*td*te 0 | + | sd*se te*sd^2 + 2*se*td*sd td*se^2 + 2*sd*te*se 0 | + | 0 -sd^2*se -sd*se^2 0 | + */ + F = mat4(td*te, -se*td-sd*te, sd*se, 0, + square(td)*te, -se*square(td)-2*sd*te*td, te*square(sd)+2*se*td*sd, -square(sd)*se, + td*square(te), -sd*square(te)-2*se*td*te, td*square(se)+2*sd*te*se, -sd*square(se), + 1, 0, 0, 0); + + result.ts[0] = vec2(td, sd); + result.ts[1] = vec2(te, se); + } + else if(d2 != 0) + { + //NOTE(martin): cusp with cusp at infinity + float tl = d3; + float sl = 3*d2; + + float invNorm = 1/sqrt(square(tl)+square(sl)); + tl *= invNorm; + sl *= invNorm; + + /*NOTE(martin): + the power basis coefficients of points k,l,m,n are collected into the rows of the 4x4 matrix F: + + | tl tl^3 1 1 | + | -sl -3sl*tl^2 0 0 | + | 0 3*sl^2*tl 0 0 | + | 0 -sl^3 0 0 | + */ + result.kind = CUBIC_CUSP_INFINITY; + + F = mat4(tl, -sl, 0, 0, + cube(tl), -3*sl*square(tl), 3*square(sl)*tl, -cube(sl), + 1, 0, 0, 0, + 1, 0, 0, 0); + + result.ts[0] = vec2(tl, sl); + result.ts[1] = vec2(0, 0); + } + else + { + //NOTE(martin): line or point degenerate case + result.kind = CUBIC_DEGENERATE_LINE; + } + + /* + F is then multiplied by M3^(-1) on the left which yelds the bezier coefficients k, l, m, n + at the control points. + + | 1 0 0 0 | + M3^(-1) = | 1 1/3 0 0 | + | 1 2/3 1/3 0 | + | 1 1 1 1 | + */ + mat4 invM3 = mat4(1, 1, 1, 1, + 0, 1./3., 2./3., 1, + 0, 0, 1./3., 1, + 0, 0, 0, 1); + + result.K = transpose(invM3*F); + + return(result); +} + +vec2 select_hull_vertex(vec2 p0, vec2 p1, vec2 p2, vec2 p3) +{ + /*NOTE: check intersection of lines (p1-p0) and (p3-p2) + P = p0 + u(p1-p0) + P = p2 + w(p3-p2) + + control points are inside a right triangle so we should always find an intersection + */ + vec2 pm; + + float det = (p1.x - p0.x)*(p3.y - p2.y) - (p1.y - p0.y)*(p3.x - p2.x); + float sqrNorm0 = dot(p1-p0, p1-p0); + float sqrNorm1 = dot(p2-p3, p2-p3); + + if(abs(det) < 1e-3 || sqrNorm0 < 0.1 || sqrNorm1 < 0.1) + { + if(sqrNorm0 < sqrNorm1) + { + pm = p2; + } + else + { + pm = p1; + } + } + else + { + float u = ((p0.x - p2.x)*(p2.y - p3.y) - (p0.y - p2.y)*(p2.x - p3.x))/det; + pm = p0 + u*(p1-p0); + } + return(pm); +} + +void cubic_emit(cubic_info curve, vec2 p[4], float s0, float s1, vec2 sp[4], int pathIndex) +{ + int segIndex = push_segment(sp, MG_GL_CUBIC, pathIndex); + + vec2 v0 = p[0]; + vec2 v1 = p[3]; + vec2 v2; + mat3 K; + + //TODO: haul that up in caller + float sqrNorm0 = dot(p[1]-p[0], p[1]-p[0]); + float sqrNorm1 = dot(p[2]-p[3], p[2]-p[3]); + + if(dot(p[0]-p[3], p[0]-p[3]) > 1e-5) + { + if(sqrNorm0 >= sqrNorm1) + { + v2 = p[1]; + K = mat3(curve.K[0].xyz, curve.K[3].xyz, curve.K[1].xyz); + } + else + { + v2 = p[2]; + K = mat3(curve.K[0].xyz, curve.K[3].xyz, curve.K[2].xyz); + } + } + else + { + v1 = p[1]; + v2 = p[2]; + K = mat3(curve.K[0].xyz, curve.K[1].xyz, curve.K[2].xyz); + } + //NOTE: set matrices + + //TODO: should we compute matrix relative to a base point to avoid loss of precision + // when computing barycentric matrix? + + mat3 B = barycentric_matrix(v0, v1, v2); + + segmentBuffer.elements[segIndex].implicitMatrix = K*B; + segmentBuffer.elements[segIndex].hullVertex = select_hull_vertex(sp[0], sp[1], sp[2], sp[3]); + + //NOTE: compute sign flip + segmentBuffer.elements[segIndex].sign = 1; + + if( curve.kind == CUBIC_SERPENTINE + || curve.kind == CUBIC_CUSP) + { + segmentBuffer.elements[segIndex].sign = (curve.d1 < 0)? -1 : 1; + } + else if(curve.kind == CUBIC_LOOP) + { + float d1 = curve.d1; + float d2 = curve.d2; + float d3 = curve.d3; + + float H0 = d3*d1-square(d2) + d1*d2*s0 - square(d1)*square(s0); + float H1 = d3*d1-square(d2) + d1*d2*s1 - square(d1)*square(s1); + float H = (abs(H0) > abs(H1)) ? H0 : H1; + segmentBuffer.elements[segIndex].sign = (H*d1 > 0) ? -1 : 1; + } + + if(sp[3].y > sp[0].y) + { + segmentBuffer.elements[segIndex].sign *= -1; + } + + //NOTE: bin to tiles + bin_to_tiles(segIndex); +} + +void cubic_setup(vec2 p[4], int pathIndex) +{ + /*NOTE(martin): first convert the control points to power basis, multiplying by M3 + + | 1 0 0 0| |p0| |c0| + M3 = |-3 3 0 0|, B = |p1|, C = |c1| = M3*B + | 3 -6 3 0| |p2| |c2| + |-1 3 -3 1| |p3| |c3| + */ + vec2 c[4] = { + p[0], + 3.0*(p[1] - p[0]), + 3.0*(p[0] + p[2] - 2*p[1]), + 3.0*(p[1] - p[2]) + p[3] - p[0]}; + + //NOTE: get classification, implicit matrix, double points and inflection points + cubic_info curve = cubic_classify(c); + + if(curve.kind == CUBIC_DEGENERATE_LINE) + { + vec2 l[4] = {p[0], p[3], vec2(0), vec2(0)}; + line_setup(l, pathIndex); + return; + } + else if(curve.kind == CUBIC_DEGENERATE_QUADRATIC) + { + vec2 quadPoint = vec2(1.5*p[1].x - 0.5*p[0].x, 1.5*p[1].y - 0.5*p[0].y); + vec2 q[4] = {p[0], quadPoint, p[3], vec2(0)}; + quadratic_setup(q, pathIndex); + return; + } + + //NOTE: get the roots of B'(s) = 3.c3.s^2 + 2.c2.s + c1 + float rootsX[2]; + int rootCountX = quadratic_roots(3*c[3].x, 2*c[2].x, c[1].x, rootsX); + + float rootsY[2]; + int rootCountY = quadratic_roots(3*c[3].y, 2*c[2].y, c[1].y, rootsY); + + float roots[6]; + for(int i=0; i=0 && roots[j]>tmp) + { + roots[j+1] = roots[j]; + j--; + } + roots[j+1] = tmp; + } + + //NOTE: compute split points + float splits[8]; + int splitCount = 0; + splits[0] = 0; + splitCount++; + for(int i=0; i 0 && roots[i] < 1) + { + splits[splitCount] = roots[i]; + splitCount++; + } + } + splits[splitCount] = 1; + splitCount++; + + //NOTE: for each monotonic segment, compute hull matrix and sign, and emit segment + for(int sliceIndex=0; sliceIndex