hmn
/
orca
5
Fork 2
Code Issues 21 Pull Requests 1 Packages Projects Releases Wiki Activity

[mtl renderer] Fixed curve slicing which used matrix operation and re-parameterization, which could create gaps in path. Now use blossoms, which ensure endpoints of subcurves match

This commit is contained in:
Martin Fouilleul 2023-04-07 17:17:55 +02:00
parent d1fab449bc
commit 65b5a4b52a
3 changed files with 352 additions and 339 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/mtl_renderer.h
View File

@ -54,7 +54,8 @@ typedef struct mg_mtl_segment
vector_float4 box; vector_float4 box;
matrix_float3x3 hullMatrix; matrix_float3x3 hullMatrix;
matrix_float3x3 implicitMatrix; matrix_float3x3 implicitMatrix;
float sign;
vector_float2 hullVertex;
int debugID; int debugID;
} mg_mtl_segment; } mg_mtl_segment;

3
src/mtl_renderer.m
View File

@ -632,8 +632,7 @@ void mg_mtl_canvas_render(mg_canvas_backend* interface,
if(primitive->cmd == MG_CMD_STROKE) if(primitive->cmd == MG_CMD_STROKE)
{ {
continue; mg_mtl_render_stroke(&context, pathElements + primitive->path.startIndex, &primitive->path);
// mg_mtl_render_stroke(&context, pathElements + primitive->path.startIndex, &primitive->path);
} }
else else
{ {

657
src/mtl_renderer.metal
View File

@ -258,6 +258,11 @@ kernel void mtl_path_setup(constant int* pathCount [[buffer(0)]],
} }
} }
float ccw(float2 a, float2 b, float2 c)
{
return((b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x));
}
int mtl_side_of_segment(float2 p, const device mg_mtl_segment* seg, mtl_log_context log = {.enabled = false}) int mtl_side_of_segment(float2 p, const device mg_mtl_segment* seg, mtl_log_context log = {.enabled = false})
{ {
int side = 0; int side = 0;
@ -318,17 +323,55 @@ int mtl_side_of_segment(float2 p, const device mg_mtl_segment* seg, mtl_log_cont
case MG_MTL_CUBIC: case MG_MTL_CUBIC:
{ {
float3 ph = {p.x, p.y, 1}; /*
float3 hullCoords = seg->hullMatrix * ph; float2 a, b, c;
if(all(hullCoords > 0)) switch(seg->config)
{ {
case MG_MTL_TL:
a = seg->box.xy;
b = seg->box.zw;
break;
case MG_MTL_BR:
a = seg->box.zw;
b = seg->box.xy;
break;
case MG_MTL_TR:
a = seg->box.xw;
b = seg->box.zy;
break;
case MG_MTL_BL:
a = seg->box.zy;
b = seg->box.xw;
break;
}
c = seg->hullVertex;
if(ccw(b, c, p) > 0 && ccw(c, a, p) > 0)
{
float3 ph = {p.x, p.y, 1};
float3 klm = seg->implicitMatrix * ph; float3 klm = seg->implicitMatrix * ph;
side = (klm.x*klm.x*klm.x - klm.y*klm.z < 0)? -1 : 1; side = (seg->sign*(klm.x*klm.x*klm.x - klm.y*klm.z) < 0)? -1 : 1;
} }
else else
{ {
side = (seg->config == MG_MTL_BL || seg->config == MG_MTL_TL) ? -1 : 1; side = (seg->config == MG_MTL_BL || seg->config == MG_MTL_TL) ? -1 : 1;
} }
*/
float3 ph = {p.x, p.y, 1};
float3 hullCoords = seg->hullMatrix * ph;
if(all(hullCoords > 0))
{
float3 klm = seg->implicitMatrix * ph;
side = (seg->sign*(klm.x*klm.x*klm.x - klm.y*klm.z) < 0)? -1 : 1;
}
else
{
side = (seg->config == MG_MTL_BL || seg->config == MG_MTL_TL) ? -1 : 1;
}
} break; } break;
} }
} }
@ -336,9 +379,9 @@ int mtl_side_of_segment(float2 p, const device mg_mtl_segment* seg, mtl_log_cont
return(side); return(side);
} }
typedef struct mtl_segment_setup_context typedef struct mtl_segment_setup_context
{ {
int pathIndex;
device atomic_int* segmentCount; device atomic_int* segmentCount;
device mg_mtl_segment* segmentBuffer; device mg_mtl_segment* segmentBuffer;
const device mg_mtl_path_queue* pathQueue; const device mg_mtl_path_queue* pathQueue;
@ -347,6 +390,9 @@ typedef struct mtl_segment_setup_context
device atomic_int* tileOpCount; device atomic_int* tileOpCount;
int tileSize; int tileSize;
mtl_log_context log; mtl_log_context log;
int pathIndex;
} mtl_segment_setup_context; } mtl_segment_setup_context;
void mtl_segment_bin_to_tiles(thread mtl_segment_setup_context* context, device mg_mtl_segment* seg) void mtl_segment_bin_to_tiles(thread mtl_segment_setup_context* context, device mg_mtl_segment* seg)
@ -428,6 +474,9 @@ void mtl_segment_bin_to_tiles(thread mtl_segment_setup_context* context, device
//NOTE: if the segment crosses the tile's bottom boundary, update the tile's winding offset //NOTE: if the segment crosses the tile's bottom boundary, update the tile's winding offset
if(crossB) if(crossB)
{ {
mtl_log(context->log, "cross bottom boundary, increment ");
mtl_log_f32(context->log, seg->windingIncrement);
mtl_log(context->log, "\n");
atomic_fetch_add_explicit(&tile->windingOffset, seg->windingIncrement, memory_order_relaxed); atomic_fetch_add_explicit(&tile->windingOffset, seg->windingIncrement, memory_order_relaxed);
} }
@ -538,22 +587,25 @@ void mtl_line_setup(thread mtl_segment_setup_context* context, float2 p[2])
mtl_segment_bin_to_tiles(context, seg); mtl_segment_bin_to_tiles(context, seg);
} }
float2 mtl_quadratic_blossom(float2 p[3], float u, float v)
{
float2 b10 = u*p[1] + (1-u)*p[0];
float2 b11 = u*p[2] + (1-u)*p[1];
float2 b20 = v*b11 + (1-v)*b10;
return(b20);
}
void mtl_quadratic_slice(float2 p[3], float s0, float s1, float2 sp[3]) void mtl_quadratic_slice(float2 p[3], float s0, float s1, float2 sp[3])
{ {
//NOTE cut curve between splitPoint[i] and splitPoint[i+1] /*NOTE: using blossoms to compute sub-curve control points ensure that the fourth point
float sr = (s1-s0)/(1-s0); 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
sp[0] = (s0-1)*(s0-1)*p[0] p[0] (and likewise, B(s=1) might not equal p[3]).
- 2*(s0-1)*s0*p[1] We handle that case explicitly to ensure that we don't create gaps in the paths.
+ s0*s0*p[2]; */
sp[0] = (s0 == 0) ? p[0] : mtl_quadratic_blossom(p, s0, s0);
sp[1] = (s0-1)*(s0-1)*(1-sr)*p[0] sp[1] = mtl_quadratic_blossom(p, s0, s1);
+ ((1-s0)*sr - 2*(s0-1)*(1-sr)*s0)*p[1] sp[2] = (s1 == 1) ? p[2] : mtl_quadratic_blossom(p, s1, s1);
+ (s0*s0*(1-sr) + s0*sr)*p[2];
sp[2] = (s0-1)*(s0-1)*(1-sr)*(1-sr)*p[0]
- 2*((s0-1)*s0*(sr-1)*(sr-1)+ (1-s0)*(sr-1)*sr)*p[1]
+ (s0*s0*(sr-1)*(sr-1) - 2*s0*(sr-1)*sr + sr*sr)*p[2];
} }
int mtl_quadratic_monotonize(float2 p[3], float splits[4]) int mtl_quadratic_monotonize(float2 p[3], float splits[4])
@ -714,152 +766,58 @@ int mtl_quadratic_roots(float a, float b, float c, thread float* r, mtl_log_cont
return(mtl_quadratic_roots_with_det(a, b, c, det, r, log)); return(mtl_quadratic_roots_with_det(a, b, c, det, r, log));
} }
float2 mtl_cubic_blossom(float2 p[4], float u, float v, float w)
{
float2 b10 = u*p[1] + (1-u)*p[0];
float2 b11 = u*p[2] + (1-u)*p[1];
float2 b12 = u*p[3] + (1-u)*p[2];
float2 b20 = v*b11 + (1-v)*b10;
float2 b21 = v*b12 + (1-v)*b11;
float2 b30 = w*b21 + (1-w)*b20;
return(b30);
}
void mtl_cubic_slice(float2 p[4], float s0, float s1, float2 sp[4]) void mtl_cubic_slice(float2 p[4], float s0, float s1, float2 sp[4])
{ {
float sr = (s1 - s0)/(1-s0); /*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.
matrix_float4x4 rightCut = {{-cube(s0-1), 0, 0, 0}, However, due to numerical errors, the evaluation of B(s=0) might not be equal to
{3*square(s0-1)*s0, square(s0-1), 0, 0}, p[0] (and likewise, B(s=1) might not equal p[3]).
{-3*(s0-1)*square(s0), -2*(s0-1)*s0, 1-s0, 0}, We handle that case explicitly to ensure that we don't create gaps in the paths.
{cube(s0), square(s0), s0, 1}}; */
sp[0] = (s0 == 0) ? p[0] : mtl_cubic_blossom(p, s0, s0, s0);
matrix_float4x4 leftCut = {{1, 1-sr, square(sr-1), -cube(sr-1)}, sp[1] = mtl_cubic_blossom(p, s0, s0, s1);
{0, sr, -2*(sr-1)*sr, 3*square(sr-1)*sr}, sp[2] = mtl_cubic_blossom(p, s0, s1, s1);
{0, 0, square(sr), -3*(sr-1)*square(sr)}, sp[3] = (s1 == 1) ? p[3] : mtl_cubic_blossom(p, s1, s1, s1);
{0, 0, 0, cube(sr)}};
float4 px = {p[0].x, p[1].x, p[2].x, p[3].x};
float4 py = {p[0].y, p[1].y, p[2].y, p[3].y};
float4 qx = leftCut*rightCut*px;
float4 qy = leftCut*rightCut*py;
sp[0] = float2(qx.x, qy.x);
sp[1] = float2(qx.y, qy.y);
sp[2] = float2(qx.z, qy.z);
sp[3] = float2(qx.w, qy.w);
} }
int mtl_cubic_monotonize(float2 p[4], float splits[8], mtl_log_context log) typedef enum {
{
//NOTE(martin): first convert the control points to power basis
float2 c[4];
c[0] = p[0];
c[1] = 3*(p[1]-p[0]);
c[2] = 3*(p[0] - 2*p[1] + p[2]);
c[3] = 3*(p[1] - p[2]) + p[3] - p[0];
//NOTE: compute the roots of the derivative
float roots[6];
int rootCount = mtl_quadratic_roots(3*c[3].x, 2*c[2].x, c[1].x, roots);
rootCount += mtl_quadratic_roots(3*c[3].y, 2*c[2].y, c[1].y, roots+rootCount);
//NOTE: compute inflection points
rootCount += mtl_quadratic_roots(3*(c[2].x*c[3].y - c[3].x*c[2].y),
3*(c[1].x*c[3].y - c[1].y*c[3].x),
(c[1].x*c[2].y - c[1].y*c[2].x),
roots + rootCount,
log);
/*
mtl_log(log, "bezier basis: ");
log_cubic_bezier(p, log);
mtl_log(log, "power basis: ");
log_cubic_bezier(c, log);
mtl_log(log, "inflection equation: ");
mtl_log_f32(log, 3*(c[2].x*c[3].y-c[3].x*c[2].y));
mtl_log(log, ", ");
mtl_log_f32(log, 3*(c[1].x*c[3].y-c[1].y*c[3].x));
mtl_log(log, ", ");
mtl_log_f32(log, (c[1].x*c[2].y-c[1].y*c[2].x));
mtl_log(log, "\n");
mtl_log(log, "inflection split count: ");
mtl_log_i32(log, rootCount-tmp);
mtl_log(log, "\n");
*/
//NOTE: sort roots
for(int i=1; i<rootCount; i++)
{
float tmp = roots[i];
int j = i-1;
while(j>=0 && roots[j]>tmp)
{
roots[j+1] = roots[j];
j--;
}
roots[j+1] = tmp;
}
//NOTE: compute split points
int splitCount = 0;
splits[0] = 0;
splitCount++;
for(int i=0; i<rootCount; i++)
{
/*
mtl_log(log, "root: ");
mtl_log_f32(log, roots[i]);
mtl_log(log, "\n");
*/
if(roots[i] > 0 && roots[i] < 1)
{
splits[splitCount] = roots[i];
splitCount++;
}
}
splits[splitCount] = 1;
splitCount++;
//NOTE: return number of split points
return(splitCount);
}
typedef enum
{
MTL_CUBIC_ERROR, MTL_CUBIC_ERROR,
MTL_CUBIC_DEGENERATE_LINE, MTL_CUBIC_SERPENTINE,
MTL_CUBIC_DEGENERATE_QUADRATIC,
MTL_CUBIC_LOOP_SPLIT,
MTL_CUBIC_LOOP_OK,
MTL_CUBIC_CUSP, MTL_CUBIC_CUSP,
MTL_CUBIC_SERPENTINE MTL_CUBIC_CUSP_INFINITY,
MTL_CUBIC_LOOP,
MTL_CUBIC_DEGENERATE_QUADRATIC,
MTL_CUBIC_DEGENERATE_LINE,
} mtl_cubic_kind; } mtl_cubic_kind;
typedef struct mtl_cubic_info typedef struct mtl_cubic_info
{ {
mtl_cubic_kind kind; mtl_cubic_kind kind;
matrix_float4x4 K; matrix_float4x4 K;
float2 quadPoint; float2 ts[2];
float split; float d1;
float d2;
float d3;
} mtl_cubic_info; } mtl_cubic_info;
mtl_cubic_info mtl_cubic_classify(thread float2* p, mtl_log_context log = {.enabled = false}) mtl_cubic_info mtl_cubic_classify(thread float2* c, mtl_log_context log = {.enabled = false})
{ {
mtl_cubic_info result = {MTL_CUBIC_ERROR}; mtl_cubic_info result = {MTL_CUBIC_ERROR};
matrix_float4x4 F; matrix_float4x4 F;
/*NOTE(martin): first convert the control points to power basis, multiplying by M3
| 1 0 0 0|
M3 = |-3 3 0 0|
| 3 -6 3 0|
|-1 3 -3 1|
ie:
c0 = p0
c1 = -3*p0 + 3*p1
c2 = 3*p0 - 6*p1 + 3*p2
c3 = -p0 + 3*p1 - 3*p2 + p3
*/
float2 c1 = 3.0*(p[1] - p[0]);
float2 c2 = 3.0*(p[0] + p[2] - 2*p[1]);
float2 c3 = 3.0*(p[1] - p[2]) + p[3] - p[0];
/*NOTE(martin): /*NOTE(martin):
now, compute determinants d0, d1, d2, d3, which gives the coefficients of the now, compute determinants d0, d1, d2, d3, which gives the coefficients of the
inflection points polynomial: inflection points polynomial:
@ -882,17 +840,13 @@ mtl_cubic_info mtl_cubic_classify(thread float2* p, mtl_log_context log = {.enab
// this may very well be an error on my part that's cancelled by flipping the signs of the d_i though! // 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);
mtl_log(log, "bezier basis: "); float d2 = -(c[3].x*c[1].y - c[3].y*c[1].x);
log_cubic_bezier(p, log); float d3 = -(c[2].y*c[1].x - c[2].x*c[1].y);
float2 c[4] = {p[0], c1, c2, c3}; result.d1 = d1;
mtl_log(log, "power basis: "); result.d2 = d2;
log_cubic_bezier(c, log); result.d3 = d3;
*/
float d1 = -(c3.y*c2.x - c3.x*c2.y);
float d2 = -(c3.x*c1.y - c3.y*c1.x);
float d3 = -(c2.y*c1.x - c2.x*c1.y);
// mtl_log(log, "d1 = "); // mtl_log(log, "d1 = ");
/* mtl_log_f32(log, d1); /* mtl_log_f32(log, d1);
@ -906,16 +860,17 @@ mtl_cubic_info mtl_cubic_classify(thread float2* p, mtl_log_context log = {.enab
float discrFactor2 = 3.0*square(d2) - 4.0*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 //NOTE(martin): each following case gives the number of roots, hence the category of the parametric curve
if(fabs(d1) < 1e-6 && fabs(d2) < 1e-6 && fabs(d3) > 1e-6) if(fabs(d1) <= 1e-6 && fabs(d2) <= 1e-6 && fabs(d3) > 1e-6)
{ {
//NOTE(martin): quadratic degenerate case //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 //NOTE(martin): compute quadratic curve control point, which is at p0 + 1.5*(p1-p0) = 1.5*p1 - 0.5*p0
result.kind = MTL_CUBIC_DEGENERATE_QUADRATIC; result.kind = MTL_CUBIC_DEGENERATE_QUADRATIC;
result.quadPoint = float2(1.5*p[1].x - 0.5*p[0].x, 1.5*p[1].y - 0.5*p[0].y);
} }
else if( (discrFactor2 > 0 && fabs(d1) > 1e-6) else if( (discrFactor2 > 0 && fabs(d1) > 1e-6)
||(discrFactor2 == 0 && fabs(d1) > 1e-6)) ||(discrFactor2 == 0 && fabs(d1) > 1e-6))
{ {
//mtl_log(log, "cusp or serpentine\n");
//NOTE(martin): serpentine curve or cusp with inflection at infinity //NOTE(martin): serpentine curve or cusp with inflection at infinity
// (these two cases are handled the same way). // (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 //NOTE(martin): compute the solutions (tl, sl), (tm, sm), and (tn, sn) of the inflection point equation
@ -950,20 +905,16 @@ mtl_cubic_info mtl_cubic_classify(thread float2* p, mtl_log_context log = {.enab
{cube(tm), -3*sm*square(tm), 3*square(sm)*tm, -cube(sm)}, {cube(tm), -3*sm*square(tm), 3*square(sm)*tm, -cube(sm)},
{1, 0, 0, 0}}; {1, 0, 0, 0}};
//NOTE: if necessary, flip sign of k and l to ensure the interior is west from the curve result.ts[0] = (float2){tm, sm};
float flip = (d1 < 0)? -1 : 1; result.ts[1] = (float2){tl, sl};
if(p[3].y > p[0].y)
{
flip *= -1;
}
F[0] *= flip;
F[1] *= flip;
} }
else if(discrFactor2 < 0 && fabs(d1) > 1e-6) else if(discrFactor2 < 0 && fabs(d1) > 1e-6)
{ {
// mtl_log(log, "loop\n");
//NOTE(martin): loop curve //NOTE(martin): loop curve
result.kind = MTL_CUBIC_LOOP;
float tetd[2]; float tetd[2];
mtl_quadratic_roots_with_det(1, -2*d2, 4*(square(d2)-d1*d3), -discrFactor2, tetd, log); mtl_quadratic_roots_with_det(1, -2*d2, 4*(square(d2)-d1*d3), -discrFactor2, tetd, log);
@ -983,7 +934,7 @@ mtl_cubic_info mtl_cubic_classify(thread float2* p, mtl_log_context log = {.enab
//NOTE(martin): if one of the parameters (td/sd) or (te/se) is in the interval [0,1], the double point //NOTE(martin): if one of the parameters (td/sd) or (te/se) is in the interval [0,1], the double point
// is inside the control points convex hull and would cause a shading anomaly. If this is // is inside the control points convex hull and would cause a shading anomaly. If this is
// the case, subdivide the curve at that point // the case, subdivide the curve at that point
//* /*
mtl_log(log, "td = "); mtl_log(log, "td = ");
mtl_log_f32(log, td); mtl_log_f32(log, td);
mtl_log(log, ", sd = "); mtl_log(log, ", sd = ");
@ -998,20 +949,6 @@ mtl_cubic_info mtl_cubic_classify(thread float2* p, mtl_log_context log = {.enab
mtl_log_f32(log, te/se); mtl_log_f32(log, te/se);
mtl_log(log, "\n"); mtl_log(log, "\n");
//*/ //*/
//TODO: investigate better margins here. The problem is that if we have a double point around 0 or 1,
// splitting the curve might also produce a root in [0, 1] due to numerical errors.
if(sd != 0 && td/sd < 0.99 && td/sd > 0.01)
{
result.kind = MTL_CUBIC_LOOP_SPLIT;
result.split = td/sd;
}
else if(se != 0 && te/se < 0.99 && te/se > 0.01)
{
result.kind = MTL_CUBIC_LOOP_SPLIT;
result.split = te/se;
}
else
{
/*NOTE(martin): /*NOTE(martin):
the power basis coefficients of points k,l,m,n are collected into the rows of the 4x4 matrix F: the power basis coefficients of points k,l,m,n are collected into the rows of the 4x4 matrix F:
@ -1020,44 +957,20 @@ mtl_cubic_info mtl_cubic_classify(thread float2* p, mtl_log_context log = {.enab
| sd*se te*sd^2 + 2*se*td*sd td*se^2 + 2*sd*te*se 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 | | 0 -sd^2*se -sd*se^2 0 |
*/ */
result.kind = MTL_CUBIC_LOOP_OK;
F = (matrix_float4x4){{td*te, -se*td-sd*te, sd*se, 0}, F = (matrix_float4x4){{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}, {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)}, {td*square(te), -sd*square(te)-2*se*td*te, td*square(se)+2*sd*te*se, -sd*square(se)},
{1, 0, 0, 0}}; {1, 0, 0, 0}};
//NOTE: if necessary, flip sign of k and l to ensure the interior is west from the curve result.ts[0] = (float2){td, sd};
float H0 = d3*d1-square(d2); result.ts[1] = (float2){te, se};
float H1 = d3*d1-square(d2) + d1*d2 - square(d1);
float H = (abs(H0) > abs(H1)) ? H0 : H1;
float flip = (H*d1 > 0) ? -1 : 1;
/*
mtl_log(log, "H0 = ");
mtl_log_f32(log, H0);
mtl_log(log, ", H1 = ");
mtl_log_f32(log, H1);
mtl_log(log, ", flip = ");
mtl_log_f32(log, flip);
mtl_log(log, "\n");
*/
if(p[3].y > p[0].y)
{
/* mtl_log(log, "fixed flip = ");
mtl_log_f32(log, flip);
mtl_log(log, "\n");
*/
flip *= -1;
}
F[0] *= flip;
F[1] *= flip;
}
} }
else if(d2 != 0) else if(d2 != 0)
{ {
//NOTE(martin): cusp with cusp at infinity //NOTE(martin): cusp with cusp at infinity
// mtl_log(log, "cusp at infinity\n");
float tl = d3; float tl = d3;
float sl = 3*d2; float sl = 3*d2;
@ -1073,17 +986,15 @@ mtl_cubic_info mtl_cubic_classify(thread float2* p, mtl_log_context log = {.enab
| 0 3*sl^2*tl 0 0 | | 0 3*sl^2*tl 0 0 |
| 0 -sl^3 0 0 | | 0 -sl^3 0 0 |
*/ */
result.kind = MTL_CUBIC_CUSP; result.kind = MTL_CUBIC_CUSP_INFINITY;
F = (matrix_float4x4){{tl, -sl, 0, 0}, F = (matrix_float4x4){{tl, -sl, 0, 0},
{cube(tl), -3*sl*square(tl), 3*square(sl)*tl, -cube(sl)}, {cube(tl), -3*sl*square(tl), 3*square(sl)*tl, -cube(sl)},
{1, 0, 0, 0}, {1, 0, 0, 0},
{1, 0, 0, 0}}; {1, 0, 0, 0}};
//NOTE: if necessary, flip sign of k and l to ensure the interior is west from the curve result.ts[0] = (float2){tl, sl};
float flip = (p[3].y > p[0].y) ? -1 : 1; result.ts[1] = (float2){0, 0};
F[0] *= flip;
F[1] *= flip;
} }
else else
{ {
@ -1110,6 +1021,42 @@ mtl_cubic_info mtl_cubic_classify(thread float2* p, mtl_log_context log = {.enab
return(result); return(result);
} }
float2 mtl_select_hull_vertex(float2 p0, float2 p1, float2 p2, float2 p3, mtl_log_context log)
{
/*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
*/
float2 pm;
float det = (p1.x - p0.x)*(p3.y - p2.y) - (p1.y - p0.y)*(p3.x - p2.x);
float sqrNorm0 = length_squared(p1-p0);
float sqrNorm1 = length_squared(p2-p3);
if(fabs(det) < 1e-3 || sqrNorm0 < 0.1 || sqrNorm1 < 0.1)
{
float sqrNorm0 = length_squared(p1-p0);
float sqrNorm1 = length_squared(p2-p3);
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);
}
matrix_float3x3 mtl_hull_matrix(float2 p0, float2 p1, float2 p2, float2 p3, mtl_log_context log) matrix_float3x3 mtl_hull_matrix(float2 p0, float2 p1, float2 p2, float2 p3, mtl_log_context log)
{ {
/*NOTE: check intersection of lines (p1-p0) and (p3-p2) /*NOTE: check intersection of lines (p1-p0) and (p3-p2)
@ -1148,9 +1095,9 @@ matrix_float3x3 mtl_hull_matrix(float2 p0, float2 p1, float2 p2, float2 p3, mtl_
return(m); return(m);
} }
void mtl_cubic_emit(thread mtl_segment_setup_context* context, float2 p[4], mtl_cubic_info info) void mtl_cubic_emit(thread mtl_segment_setup_context* context, mtl_cubic_info curve, float2 p[4], float s0, float s1, float2 sp[4])
{ {
device mg_mtl_segment* seg = mtl_segment_push(context, p, MG_MTL_CUBIC); device mg_mtl_segment* seg = mtl_segment_push(context, sp, MG_MTL_CUBIC);
float2 v0 = p[0]; float2 v0 = p[0];
float2 v1 = p[3]; float2 v1 = p[3];
@ -1160,122 +1107,174 @@ void mtl_cubic_emit(thread mtl_segment_setup_context* context, float2 p[4], mtl_
float sqrNorm0 = length_squared(p[1]-p[0]); float sqrNorm0 = length_squared(p[1]-p[0]);
float sqrNorm1 = length_squared(p[2]-p[3]); float sqrNorm1 = length_squared(p[2]-p[3]);
//TODO: should not be the local sub-curve, but the global curve!!!
if(length_squared(p[0]-p[3]) > 1e-5)
{
if(sqrNorm0 >= sqrNorm1) if(sqrNorm0 >= sqrNorm1)
{ {
v2 = p[1]; v2 = p[1];
K = {info.K[0].xyz, info.K[3].xyz, info.K[1].xyz}; K = {curve.K[0].xyz, curve.K[3].xyz, curve.K[1].xyz};
} }
else else
{ {
v2 = p[2]; v2 = p[2];
K = {info.K[0].xyz, info.K[3].xyz, info.K[2].xyz}; K = {curve.K[0].xyz, curve.K[3].xyz, curve.K[2].xyz};
} }
//NOTE: set matrices and bin segment }
else
{
v1 = p[1];
v2 = p[2];
K = {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?
matrix_float3x3 B = mtl_barycentric_matrix(v0, v1, v2); matrix_float3x3 B = mtl_barycentric_matrix(v0, v1, v2);
seg->implicitMatrix = K*B; seg->implicitMatrix = K*B;
seg->hullMatrix = mtl_hull_matrix(p[0], p[1], p[2], p[3], context->log); seg->hullMatrix = mtl_hull_matrix(sp[0], sp[1], sp[2], sp[3], context->log);
seg->hullVertex = mtl_select_hull_vertex(sp[0], sp[1], sp[2], sp[3], context->log);
//NOTE: compute sign flip
seg->sign = 1;
if(curve.kind == MTL_CUBIC_SERPENTINE
|| curve.kind == MTL_CUBIC_CUSP)
{
seg->sign = (curve.d1 < 0)? -1 : 1;
}
else if(curve.kind == MTL_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;
seg->sign = (H*d1 > 0) ? -1 : 1;
}
if(sp[3].y > sp[0].y)
{
seg->sign *= -1;
}
//NOTE: bin to tiles
mtl_segment_bin_to_tiles(context, seg); mtl_segment_bin_to_tiles(context, seg);
} }
void mtl_cubic_setup(thread mtl_segment_setup_context* context, float2 p[4]) void mtl_cubic_setup(thread mtl_segment_setup_context* context, float2 p[4])
{ {
float splits[8]; /*NOTE(martin): first convert the control points to power basis, multiplying by M3
int splitCount = mtl_cubic_monotonize(p, splits, context->log);
mtl_log(context->log, "curve = "); | 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|
*/
float2 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]};
/*
mtl_log(context->log, "bezier basis: ");
log_cubic_bezier(p, context->log); log_cubic_bezier(p, context->log);
mtl_log(context->log, "split count = "); mtl_log(context->log, "power basis: ");
log_cubic_bezier(c, context->log);
*/
//NOTE: get classification, implicit matrix, double points and inflection points
mtl_cubic_info curve = mtl_cubic_classify(c, context->log);
if(curve.kind == MTL_CUBIC_DEGENERATE_LINE)
{
float2 l[2] = {p[0], p[3]};
mtl_line_setup(context, l);
return;
}
else if(curve.kind == MTL_CUBIC_DEGENERATE_QUADRATIC)
{
float2 quadPoint = float2(1.5*p[1].x - 0.5*p[0].x, 1.5*p[1].y - 0.5*p[0].y);
float2 q[3] = {p[0], quadPoint, p[3]};
mtl_quadratic_setup(context, q);
return;
}
//NOTE: get the roots of B'(s) = 3.c3.s^2 + 2.c2.s + c1
float roots[6];
int rootCount = mtl_quadratic_roots(3*c[3].x, 2*c[2].x, c[1].x, roots);
rootCount += mtl_quadratic_roots(3*c[3].y, 2*c[2].y, c[1].y, roots + rootCount);
//NOTE: add double points and inflection points to roots if finite
for(int i=0; i<2; i++)
{
if(curve.ts[i].y)
{
roots[rootCount] = curve.ts[i].x / curve.ts[i].y;
rootCount++;
}
}
//NOTE: sort roots
for(int i=1; i<rootCount; i++)
{
float tmp = roots[i];
int j = i-1;
while(j>=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<rootCount; i++)
{
if(roots[i] > 0 && roots[i] < 1)
{
splits[splitCount] = roots[i];
splitCount++;
}
}
splits[splitCount] = 1;
splitCount++;
mtl_log(context->log, "monotonic segment count = ");
mtl_log_i32(context->log, splitCount-1); mtl_log_i32(context->log, splitCount-1);
mtl_log(context->log, "\n"); mtl_log(context->log, "\n");
//NOTE: produce bézier curve for each consecutive pair of roots //NOTE: for each monotonic segment, compute hull matrix and sign, and emit segment
for(int sliceIndex=0; sliceIndex<splitCount-1; sliceIndex++) for(int sliceIndex=0; sliceIndex<splitCount-1; sliceIndex++)
{ {
/////////////////////////////////////DEBUG float s0 = splits[sliceIndex];
/* if(sliceIndex != 0) float s1 = splits[sliceIndex+1];
{
continue;
}
*/
float2 sp[4]; float2 sp[4];
mtl_cubic_slice(p, splits[sliceIndex], splits[sliceIndex+1], sp); mtl_cubic_slice(p, s0, s1, sp);
mtl_log(context->log, "slice = "); /////////////////////// we should ensure that these have the same endpoints as the original curve and all endpoints are joining
mtl_log(context->log, "monotonic slice ");
mtl_log_i32(context->log, sliceIndex);
mtl_log(context->log, " ( ");
mtl_log_f32(context->log, s0);
mtl_log(context->log, " <= s <= ");
mtl_log_f32(context->log, s1);
mtl_log(context->log," ): ");
log_cubic_bezier(sp, context->log); log_cubic_bezier(sp, context->log);
mtl_cubic_info curve = mtl_cubic_classify(sp, context->log);
switch(curve.kind)
{
case MTL_CUBIC_ERROR:
mtl_log(context->log, "cubic curve classification error\n");
break;
case MTL_CUBIC_DEGENERATE_LINE: mtl_cubic_emit(context, curve, p, s0, s1, sp);
{
float2 l[2] = {sp[0], sp[1]};
mtl_line_setup(context, l);
} break;
case MTL_CUBIC_DEGENERATE_QUADRATIC:
{
float2 q[3] = {sp[0], curve.quadPoint, sp[3]};
mtl_quadratic_setup(context, q);
} break;
case MTL_CUBIC_LOOP_SPLIT:
{
mtl_log(context->log, "loop split: \n");
mtl_log_f32(context->log, curve.split);
mtl_log(context->log, "\n");
//NOTE: split and reclassify, check that we have a valid loop and emit
float2 ssp[8];
mtl_cubic_slice(sp, 0, curve.split, ssp);
mtl_cubic_slice(sp, curve.split, 1, ssp+4);
for(int i=0; i<2; i++)
{
mtl_cubic_info splitCurve = mtl_cubic_classify(ssp + 4*i, context->log);
mtl_log(context->log, "loop slice \n");
mtl_log_i32(context->log, i);
mtl_log(context->log, ": ");
log_cubic_bezier(ssp+i*4, context->log);
mtl_log_i32(context->log, splitCurve.kind);
mtl_log(context->log, "\n");
////////////////////////////////////////////////////////////////////////////////////
//TODO: here the result of mtl_cubic_classify seems to be changed if we print something
// inside it...
// Anyway, we shouldn't reclassify split curves, just find the new hull matrix?
////////////////////////////////////////////////////////////////////////////////////
CONTINUE_HERE;
if(splitCurve.kind == MTL_CUBIC_LOOP_SPLIT)
{
mtl_log(context->log, "loop split error (");
mtl_log_f32(context->log, splitCurve.split);
mtl_log(context->log, ") ****************************************\n");
}
else
{
mtl_cubic_emit(context, ssp + 4*i, splitCurve);
}
}
} break;
case MTL_CUBIC_LOOP_OK:
case MTL_CUBIC_CUSP:
case MTL_CUBIC_SERPENTINE:
{
mtl_cubic_emit(context, sp, curve);
} break;
}
} }
} }
@ -1296,20 +1295,20 @@ kernel void mtl_segment_setup(constant int* elementCount [[buffer(0)]],
{ {
const device mg_mtl_path_elt* elt = &elementBuffer[eltIndex]; const device mg_mtl_path_elt* elt = &elementBuffer[eltIndex];
//28 //28
// 125 // 125
// 112 // 112
if(elt->pathIndex != 124) // if(elt->pathIndex != 96)
{ {
return; // return;
} }
/*
if(elt->localEltIndex != 4)// && elt->localEltIndex != 3) if(elt->localEltIndex != 21)// && elt->localEltIndex != 22)
{ {
return; return;
} }
*/
const device mg_mtl_path_queue* pathQueue = &pathQueueBuffer[elt->pathIndex]; const device mg_mtl_path_queue* pathQueue = &pathQueueBuffer[elt->pathIndex];
device mg_mtl_tile_queue* tileQueues = &tileQueueBuffer[pathQueue->tileQueues]; device mg_mtl_tile_queue* tileQueues = &tileQueueBuffer[pathQueue->tileQueues];
@ -1324,33 +1323,31 @@ kernel void mtl_segment_setup(constant int* elementCount [[buffer(0)]],
.tileSize = tileSize[0], .tileSize = tileSize[0],
.log.buffer = logBuffer, .log.buffer = logBuffer,
.log.offset = logOffsetBuffer, .log.offset = logOffsetBuffer,
.log.enabled = true}; .log.enabled = false};
switch(elt->kind) switch(elt->kind)
{ {
case MG_MTL_LINE: case MG_MTL_LINE:
{ {
float2 p[2] = {elt->p[0]*scale[0], elt->p[1]*scale[0]}; float2 p[2] = {elt->p[0]*scale[0], elt->p[1]*scale[0]};
mtl_log(setupCtx.log, "line: "); mtl_log(setupCtx.log, "line: ");
log_line(p, setupCtx.log); log_line(p, setupCtx.log);
mtl_line_setup(&setupCtx, p); mtl_line_setup(&setupCtx, p);
} break; } break;
case MG_MTL_QUADRATIC: case MG_MTL_QUADRATIC:
{ {
float2 p[3] = {elt->p[0]*scale[0], elt->p[1]*scale[0], elt->p[2]*scale[0]}; float2 p[3] = {elt->p[0]*scale[0], elt->p[1]*scale[0], elt->p[2]*scale[0]};
mtl_log(setupCtx.log, "quadratic: "); mtl_log(setupCtx.log, "quadratic: ");
log_quadratic_bezier(p, setupCtx.log); log_quadratic_bezier(p, setupCtx.log);
mtl_quadratic_setup(&setupCtx, p); mtl_quadratic_setup(&setupCtx, p);
} break; } break;
case MG_MTL_CUBIC: case MG_MTL_CUBIC:
{ {
float2 p[4] = {elt->p[0]*scale[0], elt->p[1]*scale[0], elt->p[2]*scale[0], elt->p[3]*scale[0]}; float2 p[4] = {elt->p[0]*scale[0], elt->p[1]*scale[0], elt->p[2]*scale[0], elt->p[3]*scale[0]};
mtl_log(setupCtx.log, "cubic: ");
log_cubic_bezier(p, setupCtx.log);
mtl_cubic_setup(&setupCtx, p); mtl_cubic_setup(&setupCtx, p);
} break; } break;
@ -1518,16 +1515,30 @@ kernel void mtl_raster(const device int* screenTilesBuffer [[buffer(0)]],
pathIndex = op->index; pathIndex = op->index;
winding = op->windingOffset; winding = op->windingOffset;
if(op->next != -1) /*
if(winding < 0)
{
color = float4(1, 0, 0, 1);
}
else if(winding > 0)
{ {
color = float4(0, 1, 0, 1); color = float4(0, 1, 0, 1);
} }
if(winding && ((pixelCoord.x % 16) == 0))
{
outTexture.write(color, uint2(pixelCoord));
return;
}
if(op->next != -1)
{
color = float4(0, 0.5, 1, 1);
}
*/
} }
else if(op->kind == MG_MTL_OP_SEGMENT) else if(op->kind == MG_MTL_OP_SEGMENT)
{ {
// outTexture.write(float4(1, 0, 0, 1), uint2(pixelCoord));
// return;
const device mg_mtl_segment* seg = &segmentBuffer[op->index]; const device mg_mtl_segment* seg = &segmentBuffer[op->index];
if( (pixelCoord.y > seg->box.y) if( (pixelCoord.y > seg->box.y)
@ -1539,7 +1550,7 @@ kernel void mtl_raster(const device int* screenTilesBuffer [[buffer(0)]],
if(op->crossRight) if(op->crossRight)
{ {
color = float4(0, 1, 1, 1); //color = float4(0, 1, 1, 1);
if( (seg->config == MG_MTL_BR || seg->config == MG_MTL_TL) if( (seg->config == MG_MTL_BR || seg->config == MG_MTL_TL)
&&(pixelCoord.y > seg->box.w)) &&(pixelCoord.y > seg->box.w))
@ -1565,12 +1576,14 @@ kernel void mtl_raster(const device int* screenTilesBuffer [[buffer(0)]],
color = color*(1-pathColor.a) + pathColor; color = color*(1-pathColor.a) + pathColor;
} }
/*
if( (pixelCoord.x % tileSize[0] == 0) if( (pixelCoord.x % tileSize[0] == 0)
||(pixelCoord.y % tileSize[0] == 0)) ||(pixelCoord.y % tileSize[0] == 0))
{ {
outTexture.write(float4(0, 0, 0, 1), uint2(pixelCoord)); outTexture.write(float4(0, 0, 0, 1), uint2(pixelCoord));
return; return;
} }
*/
outTexture.write(color, uint2(pixelCoord)); outTexture.write(color, uint2(pixelCoord));
} }