Store cubic bezier control points, instead of storing them as quadratics.
Signed-off-by: Connor McAdams conmanx360@gmail.com --- dlls/d2d1/geometry.c | 474 +++++++++++++++++++++++++++++++++++++------ 1 file changed, 407 insertions(+), 67 deletions(-)
diff --git a/dlls/d2d1/geometry.c b/dlls/d2d1/geometry.c index c18b648aef..cd21ef21c5 100644 --- a/dlls/d2d1/geometry.c +++ b/dlls/d2d1/geometry.c @@ -48,8 +48,10 @@ enum d2d_vertex_type { D2D_VERTEX_TYPE_NONE, D2D_VERTEX_TYPE_LINE, - D2D_VERTEX_TYPE_BEZIER, - D2D_VERTEX_TYPE_SPLIT_BEZIER, + D2D_VERTEX_TYPE_QUADRATIC_BEZIER, + D2D_VERTEX_TYPE_CUBIC_BEZIER, + D2D_VERTEX_TYPE_SPLIT_QUAD_BEZIER, + D2D_VERTEX_TYPE_SPLIT_CUBIC_BEZIER, };
struct d2d_segment_idx @@ -409,6 +411,68 @@ static void d2d_point_normalise(D2D1_POINT_2F *p) d2d_point_scale(p, 1.0f / l); }
+static void d2d_bezier_cubic_to_quad(const D2D1_POINT_2F *p0, const D2D1_POINT_2F *p1, + const D2D1_POINT_2F *p2, const D2D1_POINT_2F *p3, D2D1_POINT_2F *c0) +{ + c0->x = (p1->x + p2->x) * 0.75f; + c0->y = (p1->y + p2->y) * 0.75f; + c0->x -= (p0->x + p3->x) * 0.25f; + c0->y -= (p0->y + p3->y) * 0.25f; +} + +static void d2d_bezier_split_cubic(const D2D1_POINT_2F *p0, const D2D1_POINT_2F *p1, + const D2D1_POINT_2F *p2, const D2D1_POINT_2F *p3, float t, D2D1_BEZIER_SEGMENT *left, + D2D1_BEZIER_SEGMENT *right, D2D1_POINT_2F *center) +{ + D2D1_POINT_2F q[4], r[3], mid; + + d2d_point_lerp(&q[0], p0, p1, t); + d2d_point_lerp(&q[1], p1, p2, t); + d2d_point_lerp(&q[2], p2, p3, t); + + d2d_point_lerp(&r[0], &q[0], &q[1], t); + d2d_point_lerp(&r[1], &q[1], &q[2], t); + d2d_point_lerp(&mid, &r[0], &r[1], t); + + if (center) + *center = mid; + + if (left) + { + left->point1 = q[0]; + left->point2 = r[0]; + left->point3 = mid; + } + + if (right) + { + right->point1 = r[1]; + right->point2 = q[2]; + right->point3 = *p3; + } +} + +static BOOL d2d_vertex_type_is_bezier(enum d2d_vertex_type t) +{ + return (t == D2D_VERTEX_TYPE_QUADRATIC_BEZIER || t == D2D_VERTEX_TYPE_CUBIC_BEZIER || + t == D2D_VERTEX_TYPE_SPLIT_QUAD_BEZIER || t == D2D_VERTEX_TYPE_SPLIT_CUBIC_BEZIER); +} + +static BOOL d2d_vertex_type_is_cubic_bezier(enum d2d_vertex_type t) +{ + return (t == D2D_VERTEX_TYPE_CUBIC_BEZIER || t == D2D_VERTEX_TYPE_SPLIT_CUBIC_BEZIER); +} + +static BOOL d2d_vertex_type_is_split_bezier(enum d2d_vertex_type t) +{ + return (t == D2D_VERTEX_TYPE_SPLIT_QUAD_BEZIER || t == D2D_VERTEX_TYPE_SPLIT_CUBIC_BEZIER); +} + +static BOOL d2d_vertex_type_is_unsplit_bezier(enum d2d_vertex_type t) +{ + return (t == D2D_VERTEX_TYPE_QUADRATIC_BEZIER || t == D2D_VERTEX_TYPE_CUBIC_BEZIER); +} + /* This implementation is based on the paper "Adaptive Precision * Floating-Point Arithmetic and Fast Robust Geometric Predicates" and * associated (Public Domain) code by Jonathan Richard Shewchuk. */ @@ -613,37 +677,71 @@ static BOOL d2d_figure_add_vertex(struct d2d_figure *figure, D2D1_POINT_2F verte return TRUE; }
-static BOOL d2d_figure_insert_bezier_control(struct d2d_figure *figure, size_t idx, const D2D1_POINT_2F *p) +static BOOL d2d_figure_insert_bezier_controls(struct d2d_figure *figure, size_t idx, size_t count, + const D2D1_POINT_2F **p) { + unsigned int i; + if (!d2d_array_reserve((void **)&figure->bezier_controls, &figure->bezier_controls_size, - figure->bezier_control_count + 1, sizeof(*figure->bezier_controls))) + figure->bezier_control_count + count, sizeof(*figure->bezier_controls))) { ERR("Failed to grow bezier controls array.\n"); return FALSE; }
- memmove(&figure->bezier_controls[idx + 1], &figure->bezier_controls[idx], + memmove(&figure->bezier_controls[idx + count], &figure->bezier_controls[idx], (figure->bezier_control_count - idx) * sizeof(*figure->bezier_controls)); - figure->bezier_controls[idx] = *p; - ++figure->bezier_control_count; + + for (i = 0; i < count; ++i) + figure->bezier_controls[idx + i] = *p[i]; + + figure->bezier_control_count += count;
return TRUE; }
-static BOOL d2d_figure_add_bezier_control(struct d2d_figure *figure, const D2D1_POINT_2F *p) +static BOOL d2d_figure_add_bezier_controls(struct d2d_figure *figure, size_t count, const D2D1_POINT_2F **p) { + unsigned int i; + if (!d2d_array_reserve((void **)&figure->bezier_controls, &figure->bezier_controls_size, - figure->bezier_control_count + 1, sizeof(*figure->bezier_controls))) + figure->bezier_control_count + count, sizeof(*figure->bezier_controls))) { ERR("Failed to grow bezier controls array.\n"); return FALSE; }
- figure->bezier_controls[figure->bezier_control_count++] = *p; + for (i = 0; i < count; ++i) + figure->bezier_controls[figure->bezier_control_count + i] = *p[i]; + + figure->bezier_control_count += count;
return TRUE; }
+static BOOL d2d_figure_insert_bezier_control(struct d2d_figure *figure, size_t idx, const D2D1_POINT_2F *p) +{ + return d2d_figure_insert_bezier_controls(figure, idx, 1, &p); +} + +static BOOL d2d_figure_add_bezier_control(struct d2d_figure *figure, const D2D1_POINT_2F *p) +{ + return d2d_figure_add_bezier_controls(figure, 1, &p); +} + +static unsigned int d2d_figure_get_bezier_control_count(struct d2d_figure *figure) +{ + unsigned int i, control_count; + + for (i = control_count = 0; i < figure->vertex_count; ++i) + { + if (d2d_vertex_type_is_bezier(figure->vertex_types[i])) + ++control_count; + } + + return control_count; +} + static void d2d_cdt_edge_rot(struct d2d_cdt_edge_ref *dst, const struct d2d_cdt_edge_ref *src) { dst->idx = src->idx; @@ -1723,15 +1821,28 @@ static BOOL d2d_geometry_intersect_bezier_line(struct d2d_geometry *geometry, const D2D1_POINT_2F *p[3], *q[2]; const struct d2d_figure *figure; float y[3], root, theta, d, e; + enum d2d_vertex_type type; + D2D1_POINT_2F tmp; size_t next;
figure = &geometry->u.path.figures[idx_p->figure_idx]; + type = figure->vertex_types[idx_p->vertex_idx]; p[0] = &figure->vertices[idx_p->vertex_idx]; - p[1] = &figure->bezier_controls[idx_p->control_idx]; next = idx_p->vertex_idx + 1; if (next == figure->vertex_count) next = 0; p[2] = &figure->vertices[next]; + if (d2d_vertex_type_is_cubic_bezier(type)) + { + d2d_bezier_cubic_to_quad(p[0], + &figure->bezier_controls[idx_p->control_idx], + &figure->bezier_controls[idx_p->control_idx + 1], + p[2], &tmp); + + p[1] = &tmp; + } + else + p[1] = &figure->bezier_controls[idx_p->control_idx];
figure = &geometry->u.path.figures[idx_q->figure_idx]; q[0] = &figure->vertices[idx_q->vertex_idx]; @@ -1799,25 +1910,48 @@ static BOOL d2d_geometry_intersect_bezier_bezier(struct d2d_geometry *geometry, const D2D1_POINT_2F *p[3], *q[3]; const struct d2d_figure *figure; D2D_RECT_F p_bounds, q_bounds; - D2D1_POINT_2F intersection; + D2D1_POINT_2F intersection, tmp_p, tmp_q; + enum d2d_vertex_type type_p, type_q; float centre_p, centre_q; size_t next;
figure = &geometry->u.path.figures[idx_p->figure_idx]; + type_p = figure->vertex_types[idx_p->vertex_idx]; p[0] = &figure->vertices[idx_p->vertex_idx]; - p[1] = &figure->bezier_controls[idx_p->control_idx]; next = idx_p->vertex_idx + 1; if (next == figure->vertex_count) next = 0; p[2] = &figure->vertices[next]; + if (d2d_vertex_type_is_cubic_bezier(type_p)) + { + d2d_bezier_cubic_to_quad(p[0], + &figure->bezier_controls[idx_p->control_idx], + &figure->bezier_controls[idx_p->control_idx + 1], + p[2], &tmp_p); + + p[1] = &tmp_p; + } + else + p[1] = &figure->bezier_controls[idx_p->control_idx];
figure = &geometry->u.path.figures[idx_q->figure_idx]; + type_q = figure->vertex_types[idx_q->vertex_idx]; q[0] = &figure->vertices[idx_q->vertex_idx]; - q[1] = &figure->bezier_controls[idx_q->control_idx]; next = idx_q->vertex_idx + 1; if (next == figure->vertex_count) next = 0; q[2] = &figure->vertices[next]; + if (d2d_vertex_type_is_cubic_bezier(type_q)) + { + d2d_bezier_cubic_to_quad(q[0], + &figure->bezier_controls[idx_q->control_idx], + &figure->bezier_controls[idx_q->control_idx + 1], + q[2], &tmp_q); + + q[1] = &tmp_q; + } + else + q[1] = &figure->bezier_controls[idx_q->control_idx];
d2d_rect_get_bezier_segment_bounds(&p_bounds, p[0], p[1], p[2], start_p, end_p); d2d_rect_get_bezier_segment_bounds(&q_bounds, q[0], q[1], q[2], start_q, end_q); @@ -1861,9 +1995,10 @@ static BOOL d2d_geometry_apply_intersections(struct d2d_geometry *geometry, { size_t vertex_offset, control_offset, next, i; struct d2d_geometry_intersection *inter; - enum d2d_vertex_type vertex_type; - const D2D1_POINT_2F *p[3]; + enum d2d_vertex_type vertex_type, split_type; + const D2D1_POINT_2F *p[4], *c[2]; struct d2d_figure *figure; + D2D1_BEZIER_SEGMENT b[2]; D2D1_POINT_2F q[2]; float t, t_prev;
@@ -1875,7 +2010,7 @@ static BOOL d2d_geometry_apply_intersections(struct d2d_geometry *geometry,
figure = &geometry->u.path.figures[inter->figure_idx]; vertex_type = figure->vertex_types[inter->vertex_idx + vertex_offset]; - if (vertex_type != D2D_VERTEX_TYPE_BEZIER && vertex_type != D2D_VERTEX_TYPE_SPLIT_BEZIER) + if (!d2d_vertex_type_is_bezier(vertex_type)) { if (!d2d_figure_insert_vertex(&geometry->u.path.figures[inter->figure_idx], inter->vertex_idx + vertex_offset + 1, inter->p)) @@ -1902,19 +2037,43 @@ static BOOL d2d_geometry_apply_intersections(struct d2d_geometry *geometry, next = inter->vertex_idx + vertex_offset + 1; if (next == figure->vertex_count) next = 0; - p[2] = &figure->vertices[next];
- d2d_point_lerp(&q[0], p[0], p[1], t); - d2d_point_lerp(&q[1], p[1], p[2], t); + if (d2d_vertex_type_is_cubic_bezier(vertex_type)) + { + p[2] = &figure->bezier_controls[inter->control_idx + control_offset + 1]; + p[3] = &figure->vertices[next]; + + d2d_bezier_split_cubic(p[0], p[1], p[2], p[3], t, &b[0], &b[1], NULL);
- figure->bezier_controls[inter->control_idx + control_offset] = q[0]; - if (!(d2d_figure_insert_bezier_control(figure, inter->control_idx + control_offset + 1, &q[1]))) - return FALSE; - ++control_offset; + figure->bezier_controls[inter->control_idx + control_offset] = b[0].point1; + figure->bezier_controls[inter->control_idx + control_offset + 1] = b[0].point2; + c[0] = &b[1].point1; + c[1] = &b[1].point2; + + if (!(d2d_figure_insert_bezier_controls(figure, inter->control_idx + control_offset + 2, 2, c))) + return FALSE; + control_offset += 2; + + split_type = D2D_VERTEX_TYPE_SPLIT_CUBIC_BEZIER; + } + else + { + p[2] = &figure->vertices[next]; + + d2d_point_lerp(&q[0], p[0], p[1], t); + d2d_point_lerp(&q[1], p[1], p[2], t); + + figure->bezier_controls[inter->control_idx + control_offset] = q[0]; + if (!(d2d_figure_insert_bezier_control(figure, inter->control_idx + control_offset + 1, &q[1]))) + return FALSE; + ++control_offset; + + split_type = D2D_VERTEX_TYPE_SPLIT_QUAD_BEZIER; + }
if (!(d2d_figure_insert_vertex(figure, inter->vertex_idx + vertex_offset + 1, inter->p))) return FALSE; - figure->vertex_types[inter->vertex_idx + vertex_offset + 1] = D2D_VERTEX_TYPE_SPLIT_BEZIER; + figure->vertex_types[inter->vertex_idx + vertex_offset + 1] = split_type; ++vertex_offset; }
@@ -1961,9 +2120,9 @@ static BOOL d2d_geometry_intersect_self(struct d2d_geometry *geometry) for (idx_q.vertex_idx = 0; idx_q.vertex_idx < max_q; ++idx_q.vertex_idx) { type_q = figure_q->vertex_types[idx_q.vertex_idx]; - if (type_q == D2D_VERTEX_TYPE_BEZIER) + if (d2d_vertex_type_is_unsplit_bezier(type_q)) { - if (type_p == D2D_VERTEX_TYPE_BEZIER) + if (d2d_vertex_type_is_unsplit_bezier(type_p)) { if (!d2d_geometry_intersect_bezier_bezier(geometry, &intersections, &idx_p, 0.0f, 1.0f, &idx_q, 0.0f, 1.0f)) @@ -1974,11 +2133,14 @@ static BOOL d2d_geometry_intersect_self(struct d2d_geometry *geometry) if (!d2d_geometry_intersect_bezier_line(geometry, &intersections, &idx_q, &idx_p)) goto done; } - ++idx_q.control_idx; + if (type_q == D2D_VERTEX_TYPE_QUADRATIC_BEZIER) + ++idx_q.control_idx; + else if (type_q == D2D_VERTEX_TYPE_CUBIC_BEZIER) + idx_q.control_idx += 2; } else { - if (type_p == D2D_VERTEX_TYPE_BEZIER) + if (d2d_vertex_type_is_unsplit_bezier(type_p)) { if (!d2d_geometry_intersect_bezier_line(geometry, &intersections, &idx_p, &idx_q)) goto done; @@ -1991,8 +2153,10 @@ static BOOL d2d_geometry_intersect_self(struct d2d_geometry *geometry) } } } - if (type_p == D2D_VERTEX_TYPE_BEZIER) + if (type_p == D2D_VERTEX_TYPE_QUADRATIC_BEZIER) ++idx_p.control_idx; + else if (type_p == D2D_VERTEX_TYPE_CUBIC_BEZIER) + idx_p.control_idx += 2; } }
@@ -2328,7 +2492,7 @@ static BOOL d2d_geometry_add_figure_outline(struct d2d_geometry *geometry, if (!i) { prev_type = figure->vertex_types[figure->vertex_count - 1]; - if (prev_type == D2D_VERTEX_TYPE_BEZIER) + if (d2d_vertex_type_is_unsplit_bezier(prev_type)) prev = &figure->bezier_controls[figure->bezier_control_count - 1]; else prev = &figure->vertices[figure->vertex_count - 1]; @@ -2336,19 +2500,22 @@ static BOOL d2d_geometry_add_figure_outline(struct d2d_geometry *geometry, else { prev_type = figure->vertex_types[i - 1]; - if (prev_type == D2D_VERTEX_TYPE_BEZIER) + if (d2d_vertex_type_is_unsplit_bezier(prev_type)) prev = &figure->bezier_controls[bezier_idx - 1]; else prev = &figure->vertices[i - 1]; }
- if (type == D2D_VERTEX_TYPE_BEZIER) + if (d2d_vertex_type_is_unsplit_bezier(type)) next = &figure->bezier_controls[bezier_idx++]; else if (i == figure->vertex_count - 1) next = &figure->vertices[0]; else next = &figure->vertices[i + 1];
+ if (type == D2D_VERTEX_TYPE_CUBIC_BEZIER) + bezier_idx++; + if (figure_end == D2D1_FIGURE_END_CLOSED || (i && i < figure->vertex_count - 1)) { D2D1_POINT_2F q_next, q_prev; @@ -2372,15 +2539,23 @@ static BOOL d2d_geometry_add_figure_outline(struct d2d_geometry *geometry, ERR("Failed to add line segment.\n"); return FALSE; } - else if (type == D2D_VERTEX_TYPE_BEZIER) + else if (d2d_vertex_type_is_unsplit_bezier(type)) { const D2D1_POINT_2F *p2; + D2D1_POINT_2F tmp;
if (i == figure->vertex_count - 1) p2 = &figure->vertices[0]; else p2 = &figure->vertices[i + 1];
+ if (type == D2D_VERTEX_TYPE_CUBIC_BEZIER) + { + d2d_bezier_cubic_to_quad(p0, &figure->bezier_controls[bezier_idx - 2], + &figure->bezier_controls[bezier_idx - 1], p2, &tmp); + next = &tmp; + } + if (!d2d_geometry_outline_add_bezier_segment(geometry, p0, next, p2)) { ERR("Failed to add bezier segment.\n"); @@ -2561,6 +2736,7 @@ static void STDMETHODCALLTYPE d2d_geometry_sink_AddBeziers(ID2D1GeometrySink *if { struct d2d_geometry *geometry = impl_from_ID2D1GeometrySink(iface); struct d2d_figure *figure = &geometry->u.path.figures[geometry->u.path.figure_count - 1]; + const D2D1_POINT_2F *c[2]; D2D1_POINT_2F p; unsigned int i;
@@ -2581,12 +2757,15 @@ static void STDMETHODCALLTYPE d2d_geometry_sink_AddBeziers(ID2D1GeometrySink *if p.y = (beziers[i].point1.y + beziers[i].point2.y) * 0.75f; p.x -= (figure->vertices[figure->vertex_count - 1].x + beziers[i].point3.x) * 0.25f; p.y -= (figure->vertices[figure->vertex_count - 1].y + beziers[i].point3.y) * 0.25f; - figure->vertex_types[figure->vertex_count - 1] = D2D_VERTEX_TYPE_BEZIER; + figure->vertex_types[figure->vertex_count - 1] = D2D_VERTEX_TYPE_CUBIC_BEZIER;
d2d_rect_get_bezier_bounds(&bezier_bounds, &figure->vertices[figure->vertex_count - 1], &p, &beziers[i].point3);
- if (!d2d_figure_add_bezier_control(figure, &p)) + c[0] = &beziers[i].point1; + c[1] = &beziers[i].point2; + + if (!d2d_figure_add_bezier_controls(figure, 2, c)) { ERR("Failed to add bezier control.\n"); geometry->u.path.state = D2D_GEOMETRY_STATE_ERROR; @@ -2659,23 +2838,29 @@ static void d2d_path_geometry_free_figures(struct d2d_geometry *geometry)
static BOOL d2d_geometry_get_bezier_segment_idx(struct d2d_geometry *geometry, struct d2d_segment_idx *idx, BOOL next) { + struct d2d_figure *figure = &geometry->u.path.figures[idx->figure_idx]; + enum d2d_vertex_type type = figure->vertex_types[idx->vertex_idx]; + if (next) { + if (d2d_vertex_type_is_cubic_bezier(type)) + ++idx->control_idx; + ++idx->vertex_idx; ++idx->control_idx; }
for (; idx->figure_idx < geometry->u.path.figure_count; ++idx->figure_idx, idx->vertex_idx = idx->control_idx = 0) { - struct d2d_figure *figure = &geometry->u.path.figures[idx->figure_idx]; + figure = &geometry->u.path.figures[idx->figure_idx];
if (!figure->bezier_control_count || figure->flags & D2D_FIGURE_FLAG_HOLLOW) continue;
for (; idx->vertex_idx < figure->vertex_count; ++idx->vertex_idx) { - if (figure->vertex_types[idx->vertex_idx] == D2D_VERTEX_TYPE_BEZIER - || figure->vertex_types[idx->vertex_idx] == D2D_VERTEX_TYPE_SPLIT_BEZIER) + type = figure->vertex_types[idx->vertex_idx]; + if (d2d_vertex_type_is_bezier(type)) return TRUE; } } @@ -2700,25 +2885,48 @@ static BOOL d2d_geometry_check_bezier_overlap(struct d2d_geometry *geometry, { const D2D1_POINT_2F *a[3], *b[3], *p[2], *q; const struct d2d_figure *figure; - D2D1_POINT_2F v_q[3], v_p, v_qp; + D2D1_POINT_2F v_q[3], v_p, v_qp, tmp_p, tmp_q; + enum d2d_vertex_type type; unsigned int i, j, score; float det, t;
figure = &geometry->u.path.figures[idx_p->figure_idx]; + type = figure->vertex_types[idx_p->vertex_idx]; a[0] = &figure->vertices[idx_p->vertex_idx]; - a[1] = &figure->bezier_controls[idx_p->control_idx]; if (idx_p->vertex_idx == figure->vertex_count - 1) a[2] = &figure->vertices[0]; else a[2] = &figure->vertices[idx_p->vertex_idx + 1]; + if (d2d_vertex_type_is_cubic_bezier(type)) + { + d2d_bezier_cubic_to_quad(a[0], + &figure->bezier_controls[idx_p->control_idx], + &figure->bezier_controls[idx_p->control_idx + 1], + a[2], &tmp_p); + + a[1] = &tmp_p; + } + else + a[1] = &figure->bezier_controls[idx_p->control_idx];
figure = &geometry->u.path.figures[idx_q->figure_idx]; + type = figure->vertex_types[idx_q->vertex_idx]; b[0] = &figure->vertices[idx_q->vertex_idx]; - b[1] = &figure->bezier_controls[idx_q->control_idx]; if (idx_q->vertex_idx == figure->vertex_count - 1) b[2] = &figure->vertices[0]; else b[2] = &figure->vertices[idx_q->vertex_idx + 1]; + if (d2d_vertex_type_is_cubic_bezier(type)) + { + d2d_bezier_cubic_to_quad(b[0], + &figure->bezier_controls[idx_q->control_idx], + &figure->bezier_controls[idx_q->control_idx + 1], + b[2], &tmp_q); + + b[1] = &tmp_q; + } + else + b[1] = &figure->bezier_controls[idx_q->control_idx];
if (d2d_point_ccw(a[0], a[1], a[2]) == 0.0f || d2d_point_ccw(b[0], b[1], b[2]) == 0.0f) return FALSE; @@ -2786,40 +2994,81 @@ static BOOL d2d_geometry_check_bezier_overlap(struct d2d_geometry *geometry, static float d2d_geometry_bezier_ccw(struct d2d_geometry *geometry, const struct d2d_segment_idx *idx) { const struct d2d_figure *figure = &geometry->u.path.figures[idx->figure_idx]; + enum d2d_vertex_type type = figure->vertex_types[idx->vertex_idx]; size_t next = idx->vertex_idx + 1; + const D2D1_POINT_2F *p; + D2D1_POINT_2F tmp;
if (next == figure->vertex_count) next = 0;
- return d2d_point_ccw(&figure->vertices[idx->vertex_idx], - &figure->bezier_controls[idx->control_idx], &figure->vertices[next]); + if (d2d_vertex_type_is_cubic_bezier(type)) + { + d2d_bezier_cubic_to_quad(&figure->vertices[idx->vertex_idx], + &figure->bezier_controls[idx->control_idx], + &figure->bezier_controls[idx->control_idx + 1], + &figure->vertices[next], &tmp); + + p = &tmp; + } + else + p = &figure->bezier_controls[idx->control_idx]; + + return d2d_point_ccw(&figure->vertices[idx->vertex_idx], p, &figure->vertices[next]); }
static BOOL d2d_geometry_split_bezier(struct d2d_geometry *geometry, const struct d2d_segment_idx *idx) { - const D2D1_POINT_2F *p[3]; + const D2D1_POINT_2F *p[4], *c[2]; struct d2d_figure *figure; - D2D1_POINT_2F q[3]; + enum d2d_vertex_type type; + D2D1_BEZIER_SEGMENT b[2]; + D2D1_POINT_2F q[2], mid; size_t next;
figure = &geometry->u.path.figures[idx->figure_idx]; + type = figure->vertex_types[idx->vertex_idx]; p[0] = &figure->vertices[idx->vertex_idx]; p[1] = &figure->bezier_controls[idx->control_idx]; next = idx->vertex_idx + 1; if (next == figure->vertex_count) next = 0; - p[2] = &figure->vertices[next];
- d2d_point_lerp(&q[0], p[0], p[1], 0.5f); - d2d_point_lerp(&q[1], p[1], p[2], 0.5f); - d2d_point_lerp(&q[2], &q[0], &q[1], 0.5f); + if (d2d_vertex_type_is_cubic_bezier(type)) + { + p[2] = &figure->bezier_controls[idx->control_idx + 1]; + p[3] = &figure->vertices[next];
- figure->bezier_controls[idx->control_idx] = q[0]; - if (!(d2d_figure_insert_bezier_control(figure, idx->control_idx + 1, &q[1]))) - return FALSE; - if (!(d2d_figure_insert_vertex(figure, idx->vertex_idx + 1, q[2]))) + d2d_bezier_split_cubic(p[0], p[1], p[2], p[3], 0.5f, &b[0], &b[1], &mid); + + figure->bezier_controls[idx->control_idx] = b[0].point1; + figure->bezier_controls[idx->control_idx + 1] = b[0].point2; + c[0] = &b[1].point1; + c[1] = &b[1].point2; + + if (!(d2d_figure_insert_bezier_controls(figure, idx->control_idx + 2, 2, c))) + return FALSE; + + type = D2D_VERTEX_TYPE_SPLIT_CUBIC_BEZIER; + } + else + { + p[2] = &figure->vertices[next]; + + d2d_point_lerp(&q[0], p[0], p[1], 0.5f); + d2d_point_lerp(&q[1], p[1], p[2], 0.5f); + d2d_point_lerp(&mid, &q[0], &q[1], 0.5f); + + figure->bezier_controls[idx->control_idx] = q[0]; + if (!(d2d_figure_insert_bezier_control(figure, idx->control_idx + 1, &q[1]))) + return FALSE; + + type = D2D_VERTEX_TYPE_SPLIT_QUAD_BEZIER; + } + + if (!(d2d_figure_insert_vertex(figure, idx->vertex_idx + 1, mid))) return FALSE; - figure->vertex_types[idx->vertex_idx + 1] = D2D_VERTEX_TYPE_SPLIT_BEZIER; + figure->vertex_types[idx->vertex_idx + 1] = type;
return TRUE; } @@ -2831,6 +3080,7 @@ static HRESULT d2d_geometry_resolve_beziers(struct d2d_geometry *geometry) const D2D1_POINT_2F *p[3]; struct d2d_figure *figure; size_t bezier_idx, i; + D2D1_POINT_2F tmp;
if (!d2d_geometry_get_first_bezier_segment_idx(geometry, &idx_p)) return S_OK; @@ -2838,6 +3088,8 @@ static HRESULT d2d_geometry_resolve_beziers(struct d2d_geometry *geometry) /* Split overlapping bezier control triangles. */ while (d2d_geometry_get_next_bezier_segment_idx(geometry, &idx_p)) { + figure = &geometry->u.path.figures[idx_p.figure_idx]; + d2d_geometry_get_first_bezier_segment_idx(geometry, &idx_q); while (idx_q.figure_idx < idx_p.figure_idx || idx_q.vertex_idx < idx_p.vertex_idx) { @@ -2849,6 +3101,9 @@ static HRESULT d2d_geometry_resolve_beziers(struct d2d_geometry *geometry) return E_OUTOFMEMORY; if (idx_p.figure_idx == idx_q.figure_idx) { + if (d2d_vertex_type_is_cubic_bezier(figure->vertex_types[idx_p.vertex_idx])) + ++idx_p.control_idx; + ++idx_p.vertex_idx; ++idx_p.control_idx; } @@ -2865,9 +3120,11 @@ static HRESULT d2d_geometry_resolve_beziers(struct d2d_geometry *geometry)
for (i = 0; i < geometry->u.path.figure_count; ++i) { - if (geometry->u.path.figures[i].flags & D2D_FIGURE_FLAG_HOLLOW) + figure = &geometry->u.path.figures[i]; + if (figure->flags & D2D_FIGURE_FLAG_HOLLOW) continue; - geometry->fill.bezier_vertex_count += 3 * geometry->u.path.figures[i].bezier_control_count; + + geometry->fill.bezier_vertex_count += 3 * d2d_figure_get_bezier_control_count(figure); }
if (!(geometry->fill.bezier_vertices = heap_calloc(geometry->fill.bezier_vertex_count, @@ -2886,9 +3143,23 @@ static HRESULT d2d_geometry_resolve_beziers(struct d2d_geometry *geometry)
figure = &geometry->u.path.figures[idx_p.figure_idx]; p[0] = &figure->vertices[idx_p.vertex_idx]; - p[1] = &figure->bezier_controls[idx_p.control_idx];
i = idx_p.vertex_idx + 1; + if (d2d_vertex_type_is_cubic_bezier(figure->vertex_types[idx_p.vertex_idx])) + { + if (i == figure->vertex_count) + p[2] = &figure->vertices[0]; + else + p[2] = &figure->vertices[i]; + + d2d_bezier_cubic_to_quad(p[0], &figure->bezier_controls[idx_p.control_idx], + &figure->bezier_controls[idx_p.control_idx + 1], p[2], &tmp); + + p[1] = &tmp; + } + else + p[1] = &figure->bezier_controls[idx_p.control_idx]; + if (d2d_path_geometry_point_inside(geometry, p[1], FALSE)) { sign = 1.0f; @@ -3002,7 +3273,7 @@ static void STDMETHODCALLTYPE d2d_geometry_sink_AddQuadraticBeziers(ID2D1Geometr d2d_rect_get_bezier_bounds(&bezier_bounds, &figure->vertices[figure->vertex_count - 1], &beziers[i].point1, &beziers[i].point2);
- figure->vertex_types[figure->vertex_count - 1] = D2D_VERTEX_TYPE_BEZIER; + figure->vertex_types[figure->vertex_count - 1] = D2D_VERTEX_TYPE_QUADRATIC_BEZIER; if (!d2d_figure_add_bezier_control(figure, &beziers[i].point1)) { ERR("Failed to add bezier.\n"); @@ -3198,7 +3469,7 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_GetBounds(ID2D1PathGeometry * for (bezier_idx = 0, ++j; j < figure->vertex_count; ++j) { if (figure->vertex_types[j] == D2D_VERTEX_TYPE_NONE - || figure->vertex_types[j] == D2D_VERTEX_TYPE_SPLIT_BEZIER) + || d2d_vertex_type_is_split_bezier(figure->vertex_types[j])) continue;
switch (type) @@ -3209,7 +3480,7 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_GetBounds(ID2D1PathGeometry * d2d_rect_expand(bounds, &p); break;
- case D2D_VERTEX_TYPE_BEZIER: + case D2D_VERTEX_TYPE_QUADRATIC_BEZIER: p1 = figure->original_bezier_controls[bezier_idx++]; d2d_point_transform(&p1, transform, p1.x, p1.y); p2 = figure->vertices[j]; @@ -3219,6 +3490,20 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_GetBounds(ID2D1PathGeometry * p = p2; break;
+ case D2D_VERTEX_TYPE_CUBIC_BEZIER: + d2d_bezier_cubic_to_quad(&figure->vertices[j - 1], + &figure->original_bezier_controls[bezier_idx], + &figure->original_bezier_controls[bezier_idx + 1], + &figure->vertices[j], &p1); + bezier_idx += 2; + d2d_point_transform(&p1, transform, p1.x, p1.y); + p2 = figure->vertices[j]; + d2d_point_transform(&p2, transform, p2.x, p2.y); + d2d_rect_get_bezier_bounds(&bezier_bounds, &p, &p1, &p2); + d2d_rect_union(bounds, &bezier_bounds); + p = p2; + break; + default: FIXME("Unhandled vertex type %#x.\n", type); p = figure->vertices[j]; @@ -3230,9 +3515,20 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_GetBounds(ID2D1PathGeometry * type = figure->vertex_types[j]; }
- if (type == D2D_VERTEX_TYPE_BEZIER) + if (d2d_vertex_type_is_unsplit_bezier(type)) { - p1 = figure->original_bezier_controls[bezier_idx++]; + if (type == D2D_VERTEX_TYPE_CUBIC_BEZIER) + { + d2d_bezier_cubic_to_quad(&figure->vertices[j - 1], + &figure->original_bezier_controls[bezier_idx], + &figure->original_bezier_controls[bezier_idx + 1], + &figure->vertices[0], &p1); + + bezier_idx += 2; + } + else + p1 = figure->original_bezier_controls[bezier_idx++]; + d2d_point_transform(&p1, transform, p1.x, p1.y); p2 = figure->vertices[0]; d2d_point_transform(&p2, transform, p2.x, p2.y); @@ -3365,6 +3661,23 @@ static void d2d_geometry_simplify_quadratic(ID2D1SimplifiedGeometrySink *sink, ID2D1SimplifiedGeometrySink_AddBeziers(sink, &b, 1); }
+static void d2d_geometry_simplify_cubic(ID2D1SimplifiedGeometrySink *sink, + D2D1_GEOMETRY_SIMPLIFICATION_OPTION option, const D2D1_POINT_2F *p0, + const D2D1_POINT_2F *p1, const D2D1_POINT_2F *p2, const D2D1_POINT_2F *p3, + float tolerance) +{ + D2D1_BEZIER_SEGMENT b; + + b.point1 = *p1; + b.point2 = *p2; + b.point3 = *p3; + + if (option == D2D1_GEOMETRY_SIMPLIFICATION_OPTION_LINES) + d2d_geometry_flatten_cubic(sink, p0, &b, tolerance); + else + ID2D1SimplifiedGeometrySink_AddBeziers(sink, &b, 1); +} + static HRESULT STDMETHODCALLTYPE d2d_path_geometry_Simplify(ID2D1PathGeometry *iface, D2D1_GEOMETRY_SIMPLIFICATION_OPTION option, const D2D1_MATRIX_3X2_F *transform, float tolerance, ID2D1SimplifiedGeometrySink *sink) @@ -3373,7 +3686,7 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_Simplify(ID2D1PathGeometry *i enum d2d_vertex_type type = D2D_VERTEX_TYPE_NONE; unsigned int i, j, bezier_idx; D2D1_FIGURE_BEGIN begin; - D2D1_POINT_2F p, p1, p2; + D2D1_POINT_2F p, p1, p2, p3; D2D1_FIGURE_END end;
TRACE("iface %p, option %#x, transform %p, tolerance %.8e, sink %p.\n", @@ -3401,7 +3714,7 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_Simplify(ID2D1PathGeometry *i for (bezier_idx = 0, ++j; j < figure->vertex_count; ++j) { if (figure->vertex_types[j] == D2D_VERTEX_TYPE_NONE - || figure->vertex_types[j] == D2D_VERTEX_TYPE_SPLIT_BEZIER) + || d2d_vertex_type_is_split_bezier(figure->vertex_types[j])) continue;
switch (type) @@ -3413,7 +3726,7 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_Simplify(ID2D1PathGeometry *i ID2D1SimplifiedGeometrySink_AddLines(sink, &p, 1); break;
- case D2D_VERTEX_TYPE_BEZIER: + case D2D_VERTEX_TYPE_QUADRATIC_BEZIER: p1 = figure->original_bezier_controls[bezier_idx++]; if (transform) d2d_point_transform(&p1, transform, p1.x, p1.y); @@ -3424,6 +3737,20 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_Simplify(ID2D1PathGeometry *i p = p2; break;
+ case D2D_VERTEX_TYPE_CUBIC_BEZIER: + p1 = figure->original_bezier_controls[bezier_idx++]; + if (transform) + d2d_point_transform(&p1, transform, p1.x, p1.y); + p2 = figure->original_bezier_controls[bezier_idx++]; + if (transform) + d2d_point_transform(&p2, transform, p2.x, p2.y); + p3 = figure->vertices[j]; + if (transform) + d2d_point_transform(&p3, transform, p3.x, p3.y); + d2d_geometry_simplify_cubic(sink, option, &p, &p1, &p2, &p3, tolerance); + p = p3; + break; + default: FIXME("Unhandled vertex type %#x.\n", type); p = figure->vertices[j]; @@ -3436,7 +3763,7 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_Simplify(ID2D1PathGeometry *i type = figure->vertex_types[j]; }
- if (type == D2D_VERTEX_TYPE_BEZIER) + if (type == D2D_VERTEX_TYPE_QUADRATIC_BEZIER) { p1 = figure->original_bezier_controls[bezier_idx++]; if (transform) @@ -3446,6 +3773,19 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_Simplify(ID2D1PathGeometry *i d2d_point_transform(&p2, transform, p2.x, p2.y); d2d_geometry_simplify_quadratic(sink, option, &p, &p1, &p2, tolerance); } + else if (type == D2D_VERTEX_TYPE_CUBIC_BEZIER) + { + p1 = figure->original_bezier_controls[bezier_idx++]; + if (transform) + d2d_point_transform(&p1, transform, p1.x, p1.y); + p2 = figure->original_bezier_controls[bezier_idx++]; + if (transform) + d2d_point_transform(&p2, transform, p2.x, p2.y); + p3 = figure->vertices[0]; + if (transform) + d2d_point_transform(&p3, transform, p3.x, p3.y); + d2d_geometry_simplify_cubic(sink, option, &p, &p1, &p2, &p3, tolerance); + }
end = figure->flags & D2D_FIGURE_FLAG_CLOSED ? D2D1_FIGURE_END_CLOSED : D2D1_FIGURE_END_OPEN; ID2D1SimplifiedGeometrySink_EndFigure(sink, end);
Add functions for calculating the bounds of a cubic bezier.
Signed-off-by: Connor McAdams conmanx360@gmail.com --- dlls/d2d1/geometry.c | 189 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 159 insertions(+), 30 deletions(-)
diff --git a/dlls/d2d1/geometry.c b/dlls/d2d1/geometry.c index cd21ef21c5..b1d22d7264 100644 --- a/dlls/d2d1/geometry.c +++ b/dlls/d2d1/geometry.c @@ -394,7 +394,7 @@ static void d2d_point_lerp(D2D1_POINT_2F *out, out->y = a->y * (1.0f - t) + b->y * t; }
-static void d2d_point_calculate_bezier(D2D1_POINT_2F *out, const D2D1_POINT_2F *p0, +static void d2d_point_calculate_quadratic_bezier(D2D1_POINT_2F *out, const D2D1_POINT_2F *p0, const D2D1_POINT_2F *p1, const D2D1_POINT_2F *p2, float t) { float t_c = 1.0f - t; @@ -403,6 +403,20 @@ static void d2d_point_calculate_bezier(D2D1_POINT_2F *out, const D2D1_POINT_2F * out->y = t_c * (t_c * p0->y + t * p1->y) + t * (t_c * p1->y + t * p2->y); }
+static void d2d_point_calculate_cubic_bezier(D2D1_POINT_2F *out, const D2D1_POINT_2F *p0, + const D2D1_POINT_2F *p1, const D2D1_POINT_2F *p2, const D2D1_POINT_2F *p3, float t) +{ + D2D1_POINT_2F q[5]; + + d2d_point_lerp(&q[0], p0, p1, t); + d2d_point_lerp(&q[1], p1, p2, t); + d2d_point_lerp(&q[2], p2, p3, t); + + d2d_point_lerp(&q[3], &q[0], &q[1], t); + d2d_point_lerp(&q[4], &q[1], &q[2], t); + d2d_point_lerp(out, &q[3], &q[4], t); +} + static void d2d_point_normalise(D2D1_POINT_2F *p) { float l; @@ -557,6 +571,37 @@ static float d2d_point_ccw(const D2D1_POINT_2F *a, const D2D1_POINT_2F *b, const return det_d[det_d_len - 1]; }
+static unsigned int d2d_solve_quadratic_roots(float a, float b, float c, float *roots) +{ + float sq_root, root, tmp; + unsigned int root_count = 0; + + tmp = b * b - 4.0f * a * c; + + if (tmp >= 0.0f) + { + sq_root = sqrtf(tmp); + + if (b >= 0.0f) + root = (-b - sq_root) / (2.0f * a); + else + root = (2.0f * c) / (-b + sq_root); + + if (root > 0.0f && root < 1.0f) + roots[root_count++] = root; + + if (b >= 0.0f) + root = (2.0f * c) / (-b - sq_root); + else + root = (-b + sq_root) / (2.0f * a); + + if (root > 0.0f && root < 1.0f) + roots[root_count++] = root; + } + + return root_count; +} + static void d2d_rect_union(D2D1_RECT_F *l, const D2D1_RECT_F *r) { l->left = min(l->left, r->left); @@ -570,7 +615,7 @@ static BOOL d2d_rect_check_overlap(const D2D_RECT_F *p, const D2D_RECT_F *q) return p->left < q->right && p->top < q->bottom && p->right > q->left && p->bottom > q->top; }
-static void d2d_rect_get_bezier_bounds(D2D_RECT_F *bounds, const D2D1_POINT_2F *p0, +static void d2d_rect_get_quadratic_bezier_bounds(D2D_RECT_F *bounds, const D2D1_POINT_2F *p0, const D2D1_POINT_2F *p1, const D2D1_POINT_2F *p2) { D2D1_POINT_2F p; @@ -592,18 +637,101 @@ static void d2d_rect_get_bezier_bounds(D2D_RECT_F *bounds, const D2D1_POINT_2F * root = (p0->x - p1->x) / (p2->x - 2.0f * p1->x + p0->x); if (root > 0.0f && root < 1.0f) { - d2d_point_calculate_bezier(&p, p0, p1, p2, root); + d2d_point_calculate_quadratic_bezier(&p, p0, p1, p2, root); d2d_rect_expand(bounds, &p); }
root = (p0->y - p1->y) / (p2->y - 2.0f * p1->y + p0->y); if (root > 0.0f && root < 1.0f) { - d2d_point_calculate_bezier(&p, p0, p1, p2, root); + d2d_point_calculate_quadratic_bezier(&p, p0, p1, p2, root); d2d_rect_expand(bounds, &p); } }
+static BOOL d2d_check_if_cubic_is_quad(const D2D1_POINT_2F *p0, const D2D1_POINT_2F *p1, + const D2D1_POINT_2F *p2, const D2D1_POINT_2F *p3) +{ + D2D1_POINT_2F tmp; + + tmp.x = -p0->x + 3.0f * p1->x - 3.0f * p2->x + p3->x; + tmp.y = -p0->y + 3.0f * p1->y - 3.0f * p2->y + p3->y; + + if ((tmp.x < 0.0005f && tmp.x > -0.0005f) && (tmp.y < 0.0005f && tmp.y > -0.0005f)) + return TRUE; + else + return FALSE; +} + +static void d2d_rect_get_cubic_bezier_bounds(D2D_RECT_F *bounds, const D2D1_POINT_2F *p0, + const D2D1_POINT_2F *p1, const D2D1_POINT_2F *p2, const D2D1_POINT_2F *p3) +{ + D2D1_POINT_2F p, v1, v2, v3; + float roots[2], a[2], b[2], c[2]; + unsigned int i, y, root_count; + + if (d2d_check_if_cubic_is_quad(p0, p1, p2, p3)) + { + d2d_bezier_cubic_to_quad(p0, p1, p2, p3, &p); + d2d_rect_get_quadratic_bezier_bounds(bounds, p0, &p, p3); + return; + } + + bounds->left = p0->x; + bounds->top = p0->y; + bounds->right = p0->x; + bounds->bottom = p0->y; + + d2d_rect_expand(bounds, p3); + + /* + * f(t) = (1 - t)³P₀ + 3(1 - t)²tP₁ + 3(1 - t)t²P₂ + t³P₃ + * f'(t) = 3(1 - t)²(P₁ - P₀) + 6(1 - t)t(P₂ - P₁) + 3t²(P₃ - P₂) + * + * Or, from https://pomax.github.io/bezierinfo/#extremities + * V₁ = 3(P₁ - P₀) + * V₂ = 3(P₂ - P₁) + * V₃ = 3(P₃ - P₂) + * f'(t) = V₁(1 - t)² + 2V₂(1 - t)t + V₃t² + * = (V₁ - 2V₂ + V₃)t² + 2(V₂ - V₁)t + V₁ + * + * And new quadratic coefficients a, b, and c are: + * a = V₁ - 2V₂ + V₃ + * b = 2(V₂ - V₁) + * c = V₁ + * + * f'(t) = 0 + * t = (-b ± √(b² - 4ac)) / 2a + */ + d2d_point_subtract(&v1, p1, p0); + d2d_point_scale(&v1, 3.0f); + + d2d_point_subtract(&v2, p2, p1); + d2d_point_scale(&v2, 3.0f); + + d2d_point_subtract(&v3, p3, p2); + d2d_point_scale(&v3, 3.0f); + + a[0] = v1.x - 2.0f * v2.x + v3.x; + a[1] = v1.y - 2.0f * v2.y + v3.y; + + b[0] = 2.0f * (v2.x - v1.x); + b[1] = 2.0f * (v2.y - v1.y); + + c[0] = v1.x; + c[1] = v1.y; + + for (i = 0; i < 2; ++i) + { + root_count = d2d_solve_quadratic_roots(a[i], b[i], c[i], roots); + for (y = 0; y < root_count; ++y) + { + d2d_point_calculate_cubic_bezier(&p, p0, p1, p2, p3, roots[y]); + d2d_rect_expand(bounds, &p); + } + } +} + static void d2d_rect_get_bezier_segment_bounds(D2D_RECT_F *bounds, const D2D1_POINT_2F *p0, const D2D1_POINT_2F *p1, const D2D1_POINT_2F *p2, float start, float end) { @@ -618,7 +746,7 @@ static void d2d_rect_get_bezier_segment_bounds(D2D_RECT_F *bounds, const D2D1_PO d2d_point_lerp(&r[0], &r[1], p2, end); d2d_point_lerp(&q[2], &q[1], &r[0], end);
- d2d_rect_get_bezier_bounds(bounds, &q[0], &q[1], &q[2]); + d2d_rect_get_quadratic_bezier_bounds(bounds, &q[0], &q[1], &q[2]); }
static BOOL d2d_figure_insert_vertex(struct d2d_figure *figure, size_t idx, D2D1_POINT_2F vertex) @@ -1795,7 +1923,7 @@ static BOOL d2d_geometry_add_bezier_line_intersections(struct d2d_geometry *geom D2D1_POINT_2F intersection; float t;
- d2d_point_calculate_bezier(&intersection, p[0], p[1], p[2], s); + d2d_point_calculate_quadratic_bezier(&intersection, p[0], p[1], p[2], s); if (fabsf(q[1]->x - q[0]->x) > fabsf(q[1]->y - q[0]->y)) t = (intersection.x - q[0]->x) / (q[1]->x - q[0]->x); else @@ -1964,7 +2092,7 @@ static BOOL d2d_geometry_intersect_bezier_bezier(struct d2d_geometry *geometry,
if (end_p - start_p < 1e-3f) { - d2d_point_calculate_bezier(&intersection, p[0], p[1], p[2], centre_p); + d2d_point_calculate_quadratic_bezier(&intersection, p[0], p[1], p[2], centre_p); if (start_p > 0.0f && end_p < 1.0f && !d2d_geometry_intersections_add(intersections, idx_p, centre_p, intersection)) return FALSE; @@ -2759,8 +2887,8 @@ static void STDMETHODCALLTYPE d2d_geometry_sink_AddBeziers(ID2D1GeometrySink *if p.y -= (figure->vertices[figure->vertex_count - 1].y + beziers[i].point3.y) * 0.25f; figure->vertex_types[figure->vertex_count - 1] = D2D_VERTEX_TYPE_CUBIC_BEZIER;
- d2d_rect_get_bezier_bounds(&bezier_bounds, &figure->vertices[figure->vertex_count - 1], - &p, &beziers[i].point3); + d2d_rect_get_cubic_bezier_bounds(&bezier_bounds, &figure->vertices[figure->vertex_count - 1], + &beziers[i].point1, &beziers[i].point2, &beziers[i].point3);
c[0] = &beziers[i].point1; c[1] = &beziers[i].point2; @@ -3270,7 +3398,7 @@ static void STDMETHODCALLTYPE d2d_geometry_sink_AddQuadraticBeziers(ID2D1Geometr { D2D1_RECT_F bezier_bounds;
- d2d_rect_get_bezier_bounds(&bezier_bounds, &figure->vertices[figure->vertex_count - 1], + d2d_rect_get_quadratic_bezier_bounds(&bezier_bounds, &figure->vertices[figure->vertex_count - 1], &beziers[i].point1, &beziers[i].point2);
figure->vertex_types[figure->vertex_count - 1] = D2D_VERTEX_TYPE_QUADRATIC_BEZIER; @@ -3440,7 +3568,7 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_GetBounds(ID2D1PathGeometry * const struct d2d_figure *figure = &geometry->u.path.figures[i]; enum d2d_vertex_type type = D2D_VERTEX_TYPE_NONE; D2D1_RECT_F bezier_bounds; - D2D1_POINT_2F p, p1, p2; + D2D1_POINT_2F p, p1, p2, p3; size_t j, bezier_idx;
if (figure->flags & D2D_FIGURE_FLAG_HOLLOW) @@ -3485,23 +3613,21 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_GetBounds(ID2D1PathGeometry * d2d_point_transform(&p1, transform, p1.x, p1.y); p2 = figure->vertices[j]; d2d_point_transform(&p2, transform, p2.x, p2.y); - d2d_rect_get_bezier_bounds(&bezier_bounds, &p, &p1, &p2); + d2d_rect_get_quadratic_bezier_bounds(&bezier_bounds, &p, &p1, &p2); d2d_rect_union(bounds, &bezier_bounds); p = p2; break;
case D2D_VERTEX_TYPE_CUBIC_BEZIER: - d2d_bezier_cubic_to_quad(&figure->vertices[j - 1], - &figure->original_bezier_controls[bezier_idx], - &figure->original_bezier_controls[bezier_idx + 1], - &figure->vertices[j], &p1); - bezier_idx += 2; + p1 = figure->original_bezier_controls[bezier_idx++]; d2d_point_transform(&p1, transform, p1.x, p1.y); - p2 = figure->vertices[j]; + p2 = figure->original_bezier_controls[bezier_idx++]; d2d_point_transform(&p2, transform, p2.x, p2.y); - d2d_rect_get_bezier_bounds(&bezier_bounds, &p, &p1, &p2); + p3 = figure->vertices[j]; + d2d_point_transform(&p3, transform, p3.x, p3.y); + d2d_rect_get_cubic_bezier_bounds(&bezier_bounds, &p, &p1, &p2, &p3); d2d_rect_union(bounds, &bezier_bounds); - p = p2; + p = p3; break;
default: @@ -3519,20 +3645,23 @@ static HRESULT STDMETHODCALLTYPE d2d_path_geometry_GetBounds(ID2D1PathGeometry * { if (type == D2D_VERTEX_TYPE_CUBIC_BEZIER) { - d2d_bezier_cubic_to_quad(&figure->vertices[j - 1], - &figure->original_bezier_controls[bezier_idx], - &figure->original_bezier_controls[bezier_idx + 1], - &figure->vertices[0], &p1); - - bezier_idx += 2; + p1 = figure->original_bezier_controls[bezier_idx++]; + d2d_point_transform(&p1, transform, p1.x, p1.y); + p2 = figure->original_bezier_controls[bezier_idx++]; + d2d_point_transform(&p2, transform, p2.x, p2.y); + p3 = figure->vertices[0]; + d2d_point_transform(&p3, transform, p3.x, p3.y); + d2d_rect_get_cubic_bezier_bounds(&bezier_bounds, &p, &p1, &p2, &p3); } else + { p1 = figure->original_bezier_controls[bezier_idx++]; + d2d_point_transform(&p1, transform, p1.x, p1.y); + p2 = figure->vertices[0]; + d2d_point_transform(&p2, transform, p2.x, p2.y); + d2d_rect_get_quadratic_bezier_bounds(&bezier_bounds, &p, &p1, &p2); + }
- d2d_point_transform(&p1, transform, p1.x, p1.y); - p2 = figure->vertices[0]; - d2d_point_transform(&p2, transform, p2.x, p2.y); - d2d_rect_get_bezier_bounds(&bezier_bounds, &p, &p1, &p2); d2d_rect_union(bounds, &bezier_bounds); } }
Signed-off-by: Connor McAdams conmanx360@gmail.com --- dlls/d2d1/geometry.c | 69 +++++++++++++++++++++++++------------------- 1 file changed, 40 insertions(+), 29 deletions(-)
diff --git a/dlls/d2d1/geometry.c b/dlls/d2d1/geometry.c index b1d22d7264..5fe33246b2 100644 --- a/dlls/d2d1/geometry.c +++ b/dlls/d2d1/geometry.c @@ -732,8 +732,9 @@ static void d2d_rect_get_cubic_bezier_bounds(D2D_RECT_F *bounds, const D2D1_POIN } }
-static void d2d_rect_get_bezier_segment_bounds(D2D_RECT_F *bounds, const D2D1_POINT_2F *p0, - const D2D1_POINT_2F *p1, const D2D1_POINT_2F *p2, float start, float end) +static void d2d_rect_get_quadratic_bezier_segment_bounds(D2D_RECT_F *bounds, + const D2D1_POINT_2F *p0, const D2D1_POINT_2F *p1, const D2D1_POINT_2F *p2, + float start, float end) { D2D1_POINT_2F q[3], r[2];
@@ -749,6 +750,21 @@ static void d2d_rect_get_bezier_segment_bounds(D2D_RECT_F *bounds, const D2D1_PO d2d_rect_get_quadratic_bezier_bounds(bounds, &q[0], &q[1], &q[2]); }
+static void d2d_rect_get_cubic_bezier_segment_bounds(D2D_RECT_F *bounds, + const D2D1_POINT_2F *p0, const D2D1_POINT_2F *p1, const D2D1_POINT_2F *p2, + const D2D1_POINT_2F *p3, float start, float end) +{ + D2D1_BEZIER_SEGMENT left, right; + D2D1_POINT_2F mid; + + d2d_bezier_split_cubic(p0, p1, p2, p3, start, NULL, &right, &mid); + + end = (end - start) / (1.0f - start); + d2d_bezier_split_cubic(&mid, &right.point1, &right.point2, p3, end, &left, NULL, NULL); + + d2d_rect_get_cubic_bezier_bounds(bounds, &mid, &left.point1, &left.point2, &left.point3); +} + static BOOL d2d_figure_insert_vertex(struct d2d_figure *figure, size_t idx, D2D1_POINT_2F vertex) { if (!d2d_array_reserve((void **)&figure->vertices, &figure->vertices_size, @@ -2035,10 +2051,10 @@ static BOOL d2d_geometry_intersect_bezier_bezier(struct d2d_geometry *geometry, const struct d2d_segment_idx *idx_p, float start_p, float end_p, const struct d2d_segment_idx *idx_q, float start_q, float end_q) { - const D2D1_POINT_2F *p[3], *q[3]; + const D2D1_POINT_2F *p[4], *q[4]; const struct d2d_figure *figure; D2D_RECT_F p_bounds, q_bounds; - D2D1_POINT_2F intersection, tmp_p, tmp_q; + D2D1_POINT_2F intersection; enum d2d_vertex_type type_p, type_q; float centre_p, centre_q; size_t next; @@ -2046,43 +2062,34 @@ static BOOL d2d_geometry_intersect_bezier_bezier(struct d2d_geometry *geometry, figure = &geometry->u.path.figures[idx_p->figure_idx]; type_p = figure->vertex_types[idx_p->vertex_idx]; p[0] = &figure->vertices[idx_p->vertex_idx]; + p[1] = &figure->bezier_controls[idx_p->control_idx]; + if (type_p == D2D_VERTEX_TYPE_CUBIC_BEZIER) + p[2] = &figure->bezier_controls[idx_p->control_idx + 1]; next = idx_p->vertex_idx + 1; if (next == figure->vertex_count) next = 0; - p[2] = &figure->vertices[next]; - if (d2d_vertex_type_is_cubic_bezier(type_p)) - { - d2d_bezier_cubic_to_quad(p[0], - &figure->bezier_controls[idx_p->control_idx], - &figure->bezier_controls[idx_p->control_idx + 1], - p[2], &tmp_p); - - p[1] = &tmp_p; - } - else - p[1] = &figure->bezier_controls[idx_p->control_idx]; + p[3] = &figure->vertices[next];
figure = &geometry->u.path.figures[idx_q->figure_idx]; type_q = figure->vertex_types[idx_q->vertex_idx]; q[0] = &figure->vertices[idx_q->vertex_idx]; + q[1] = &figure->bezier_controls[idx_q->control_idx]; + if (type_q == D2D_VERTEX_TYPE_CUBIC_BEZIER) + q[2] = &figure->bezier_controls[idx_q->control_idx + 1]; next = idx_q->vertex_idx + 1; if (next == figure->vertex_count) next = 0; - q[2] = &figure->vertices[next]; - if (d2d_vertex_type_is_cubic_bezier(type_q)) - { - d2d_bezier_cubic_to_quad(q[0], - &figure->bezier_controls[idx_q->control_idx], - &figure->bezier_controls[idx_q->control_idx + 1], - q[2], &tmp_q); + q[3] = &figure->vertices[next];
- q[1] = &tmp_q; - } + if (type_p == D2D_VERTEX_TYPE_CUBIC_BEZIER) + d2d_rect_get_cubic_bezier_segment_bounds(&p_bounds, p[0], p[1], p[2], p[3], start_p, end_p); else - q[1] = &figure->bezier_controls[idx_q->control_idx]; + d2d_rect_get_quadratic_bezier_segment_bounds(&p_bounds, p[0], p[1], p[3], start_p, end_p);
- d2d_rect_get_bezier_segment_bounds(&p_bounds, p[0], p[1], p[2], start_p, end_p); - d2d_rect_get_bezier_segment_bounds(&q_bounds, q[0], q[1], q[2], start_q, end_q); + if (type_q == D2D_VERTEX_TYPE_CUBIC_BEZIER) + d2d_rect_get_cubic_bezier_segment_bounds(&q_bounds, q[0], q[1], q[2], q[3], start_q, end_q); + else + d2d_rect_get_quadratic_bezier_segment_bounds(&q_bounds, q[0], q[1], q[3], start_q, end_q);
if (!d2d_rect_check_overlap(&p_bounds, &q_bounds)) return TRUE; @@ -2092,7 +2099,11 @@ static BOOL d2d_geometry_intersect_bezier_bezier(struct d2d_geometry *geometry,
if (end_p - start_p < 1e-3f) { - d2d_point_calculate_quadratic_bezier(&intersection, p[0], p[1], p[2], centre_p); + if (type_p == D2D_VERTEX_TYPE_CUBIC_BEZIER) + d2d_point_calculate_cubic_bezier(&intersection, p[0], p[1], p[2], p[3], centre_p); + else + d2d_point_calculate_quadratic_bezier(&intersection, p[0], p[1], p[3], centre_p); + if (start_p > 0.0f && end_p < 1.0f && !d2d_geometry_intersections_add(intersections, idx_p, centre_p, intersection)) return FALSE;
Signed-off-by: Connor McAdams conmanx360@gmail.com --- dlls/d2d1/geometry.c | 216 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 181 insertions(+), 35 deletions(-)
diff --git a/dlls/d2d1/geometry.c b/dlls/d2d1/geometry.c index 5fe33246b2..6d7a904898 100644 --- a/dlls/d2d1/geometry.c +++ b/dlls/d2d1/geometry.c @@ -1936,10 +1936,18 @@ static BOOL d2d_geometry_add_bezier_line_intersections(struct d2d_geometry *geom struct d2d_geometry_intersections *intersections, const struct d2d_segment_idx *idx_p, const D2D1_POINT_2F **p, const struct d2d_segment_idx *idx_q, const D2D1_POINT_2F **q, float s) { + const struct d2d_figure *figure; + enum d2d_vertex_type type; D2D1_POINT_2F intersection; float t;
- d2d_point_calculate_quadratic_bezier(&intersection, p[0], p[1], p[2], s); + figure = &geometry->u.path.figures[idx_p->figure_idx]; + type = figure->vertex_types[idx_p->vertex_idx]; + if (type == D2D_VERTEX_TYPE_CUBIC_BEZIER) + d2d_point_calculate_cubic_bezier(&intersection, p[0], p[1], p[2], p[3], s); + else + d2d_point_calculate_quadratic_bezier(&intersection, p[0], p[1], p[3], s); + if (fabsf(q[1]->x - q[0]->x) > fabsf(q[1]->y - q[0]->y)) t = (intersection.x - q[0]->x) / (q[1]->x - q[0]->x); else @@ -1958,48 +1966,18 @@ static BOOL d2d_geometry_add_bezier_line_intersections(struct d2d_geometry *geom return TRUE; }
-static BOOL d2d_geometry_intersect_bezier_line(struct d2d_geometry *geometry, +static BOOL d2d_geometry_get_quadratic_bezier_roots(struct d2d_geometry *geometry, struct d2d_geometry_intersections *intersections, - const struct d2d_segment_idx *idx_p, const struct d2d_segment_idx *idx_q) + const struct d2d_segment_idx *idx_p, const D2D1_POINT_2F **p, + const struct d2d_segment_idx *idx_q, const D2D1_POINT_2F **q) { - const D2D1_POINT_2F *p[3], *q[2]; - const struct d2d_figure *figure; float y[3], root, theta, d, e; - enum d2d_vertex_type type; - D2D1_POINT_2F tmp; - size_t next; - - figure = &geometry->u.path.figures[idx_p->figure_idx]; - type = figure->vertex_types[idx_p->vertex_idx]; - p[0] = &figure->vertices[idx_p->vertex_idx]; - next = idx_p->vertex_idx + 1; - if (next == figure->vertex_count) - next = 0; - p[2] = &figure->vertices[next]; - if (d2d_vertex_type_is_cubic_bezier(type)) - { - d2d_bezier_cubic_to_quad(p[0], - &figure->bezier_controls[idx_p->control_idx], - &figure->bezier_controls[idx_p->control_idx + 1], - p[2], &tmp); - - p[1] = &tmp; - } - else - p[1] = &figure->bezier_controls[idx_p->control_idx]; - - figure = &geometry->u.path.figures[idx_q->figure_idx]; - q[0] = &figure->vertices[idx_q->vertex_idx]; - next = idx_q->vertex_idx + 1; - if (next == figure->vertex_count) - next = 0; - q[1] = &figure->vertices[next];
/* Align the line with x-axis. */ theta = -atan2f(q[1]->y - q[0]->y, q[1]->x - q[0]->x); y[0] = (p[0]->x - q[0]->x) * sinf(theta) + (p[0]->y - q[0]->y) * cosf(theta); y[1] = (p[1]->x - q[0]->x) * sinf(theta) + (p[1]->y - q[0]->y) * cosf(theta); - y[2] = (p[2]->x - q[0]->x) * sinf(theta) + (p[2]->y - q[0]->y) * cosf(theta); + y[2] = (p[3]->x - q[0]->x) * sinf(theta) + (p[3]->y - q[0]->y) * cosf(theta);
/* Intersect the transformed curve with the x-axis. * @@ -2046,6 +2024,174 @@ static BOOL d2d_geometry_intersect_bezier_line(struct d2d_geometry *geometry, return TRUE; }
+/* Cubic root finding method adapted from code at + * https://pomax.github.io/bezierinfo/#extremities */ +static float d2d_cubic_bezier_cuberoot(float a) +{ + if (a < 0.0f) + return -powf(-a, 1.0f / 3.0f); + else + return powf(a, 1.0f / 3.0f); +} + +static int d2d_cubic_bezier_get_roots(float *y, float *roots) +{ + float mp3, r, t, cosphi, phi, crtr, t1, u1, v1; + float p, p_3, q, q2, disc, sd, tmp, a, b, c, d; + unsigned int root_count; + + /* First, we need to convert the bezier coefficients to 'power basis' + * coefficients so that we can use a generic cubic root solving equation. */ + a = (3.0f * y[0] - 6.0f * y[1] + 3.0f * y[2]); + b = (-3.0f * y[0] + 3.0f * y[1]); + c = y[0]; + d = (-y[0] + 3.0f * y[1] - 3.0f * y[2] + y[3]); + + /* Check if the curve is actually a quadratic. */ + if ((d < 0.0005f && d > -0.0005f)) + { + /* Check if it's actually a line. */ + if ((a < 0.0005f && a > -0.0005f)) + { + /* Check if it's just a point. If it is, no roots. */ + if ((b < 0.0005f && b > -0.0005f)) + return 0; + + roots[0] = -c / b; + return 1; + } + + return d2d_solve_quadratic_roots(a, b, c, roots); + } + + a /= d; + b /= d; + c /= d; + + p = (3.0f * b - a * a) / 3.0f; + p_3 = p / 3.0f; + q = (2.0f * a * a * a - 9.0f * a * b + 27.0f * c) / 27.0f; + q2 = q / 2.0f; + disc = q2 * q2 + p_3 * p_3 * p_3; + + if ((disc < 0.0005f && disc > -0.0005f)) + disc = 0.0f; + + if (disc < 0.0f) + { + /* Three real roots. */ + mp3 = -p / 3.0f; + r = sqrtf(mp3 * mp3 * mp3); + t = -q / (2.0f * r); + + if (t < -1.0f) + cosphi = -1.0f; + else if (t > 1.0f) + cosphi = 1.0f; + else + cosphi = t; + + phi = acosf(cosphi); + crtr = d2d_cubic_bezier_cuberoot(r); + t1 = 2.0f * crtr; + + roots[0] = t1 * cosf(phi / 3) - a / 3.0f; + roots[1] = t1 * cosf((phi + 2 * M_PI) / 3) - a / 3.0f; + roots[2] = t1 * cosf((phi + 4 * M_PI) / 3) - a / 3.0f; + + root_count = 3; + } + else if (disc == 0.0f) + { + /* Three real roots, but two are equal. */ + if (q2 < 0.0f) + tmp = d2d_cubic_bezier_cuberoot(-q2); + else + tmp = -d2d_cubic_bezier_cuberoot(q2); + + roots[0] = 2.0f * tmp - (a / 3.0f); + roots[1] = -tmp - (a / 3.0f); + + root_count = 2; + } + else + { + /* One real root, and two complex roots. */ + sd = sqrtf(disc); + u1 = d2d_cubic_bezier_cuberoot(sd - q2); + v1 = d2d_cubic_bezier_cuberoot(sd + q2); + roots[0] = u1 - v1 - a / 3.0f; + + root_count = 1; + } + + return root_count; +} + +static BOOL d2d_geometry_get_cubic_bezier_roots(struct d2d_geometry *geometry, + struct d2d_geometry_intersections *intersections, + const struct d2d_segment_idx *idx_p, const D2D1_POINT_2F **p, + const struct d2d_segment_idx *idx_q, const D2D1_POINT_2F **q) +{ + float y[4], roots[3], theta; + unsigned int num_roots, i; + + /* Align the line with x-axis. */ + theta = -atan2f(q[1]->y - q[0]->y, q[1]->x - q[0]->x); + /* Rotate the Y coordinates of the cubic bezier. */ + y[0] = (p[0]->x - q[0]->x) * sinf(theta) + (p[0]->y - q[0]->y) * cosf(theta); + y[1] = (p[1]->x - q[0]->x) * sinf(theta) + (p[1]->y - q[0]->y) * cosf(theta); + y[2] = (p[2]->x - q[0]->x) * sinf(theta) + (p[2]->y - q[0]->y) * cosf(theta); + y[3] = (p[3]->x - q[0]->x) * sinf(theta) + (p[3]->y - q[0]->y) * cosf(theta); + + num_roots = d2d_cubic_bezier_get_roots(y, roots); + + for (i = 0; i < num_roots; ++i) + { + if (roots[i] >= 0.0f && roots[i] <= 1.0f) + { + if (!d2d_geometry_add_bezier_line_intersections(geometry, + intersections, idx_p, p, idx_q, q, roots[i])) + return FALSE; + } + } + + return TRUE; +} + +static BOOL d2d_geometry_intersect_bezier_line(struct d2d_geometry *geometry, + struct d2d_geometry_intersections *intersections, + const struct d2d_segment_idx *idx_p, const struct d2d_segment_idx *idx_q) +{ + const D2D1_POINT_2F *p[4], *q[2]; + const struct d2d_figure *figure; + enum d2d_vertex_type type; + size_t next; + + figure = &geometry->u.path.figures[idx_p->figure_idx]; + type = figure->vertex_types[idx_p->vertex_idx]; + p[0] = &figure->vertices[idx_p->vertex_idx]; + p[1] = &figure->bezier_controls[idx_p->control_idx]; + if (type == D2D_VERTEX_TYPE_CUBIC_BEZIER) + p[2] = &figure->bezier_controls[idx_p->control_idx + 1]; + next = idx_p->vertex_idx + 1; + if (next == figure->vertex_count) + next = 0; + p[3] = &figure->vertices[next]; + + figure = &geometry->u.path.figures[idx_q->figure_idx]; + q[0] = &figure->vertices[idx_q->vertex_idx]; + next = idx_q->vertex_idx + 1; + if (next == figure->vertex_count) + next = 0; + q[1] = &figure->vertices[next]; + + if (type == D2D_VERTEX_TYPE_CUBIC_BEZIER) + return d2d_geometry_get_cubic_bezier_roots(geometry, intersections, idx_p, p, idx_q, q); + else + return d2d_geometry_get_quadratic_bezier_roots(geometry, intersections, idx_p, p, idx_q, q); +} + static BOOL d2d_geometry_intersect_bezier_bezier(struct d2d_geometry *geometry, struct d2d_geometry_intersections *intersections, const struct d2d_segment_idx *idx_p, float start_p, float end_p,
Hi,
Il 31/03/20 22:11, Connor McAdams ha scritto:
- p = (3.0f * b - a * a) / 3.0f;
- p_3 = p / 3.0f;
- q = (2.0f * a * a * a - 9.0f * a * b + 27.0f * c) / 27.0f;
- q2 = q / 2.0f;
- disc = q2 * q2 + p_3 * p_3 * p_3;
I don't think it is advisable to use Cardano's formula to solve a cubic, because it can go numerically bad in a lot of ways. After a quick chat with a friend into numerical computing, it seems that it is probably better to do a few iterations of a numerical method, like the Aberth method (see the Wikipedia page), which should converge rather quickly and more or less always if the initial values are chosen not to be symmetrical around the real axis. However, it requires to use complex numbers.
A good criterion for halting the iterative search is when abs(p/p') is smaller than epsilon, because that condition implies that a root is in a neighbourhood of size approximately epsilon.
HTH, Giovanni.
On 4/1/20 19:04, Giovanni Mascellani wrote:
Hi,
Il 31/03/20 22:11, Connor McAdams ha scritto:
- p = (3.0f * b - a * a) / 3.0f;
- p_3 = p / 3.0f;
- q = (2.0f * a * a * a - 9.0f * a * b + 27.0f * c) / 27.0f;
- q2 = q / 2.0f;
- disc = q2 * q2 + p_3 * p_3 * p_3;
I don't think it is advisable to use Cardano's formula to solve a cubic, because it can go numerically bad in a lot of ways. After a quick chat with a friend into numerical computing, it seems that it is probably better to do a few iterations of a numerical method, like the Aberth method (see the Wikipedia page), which should converge rather quickly and more or less always if the initial values are chosen not to be symmetrical around the real axis. However, it requires to use complex numbers.
Given the complex roots are not needed here and the polynomial is always cubic, is this generic method really beneficial? It would probably be simpler and quicker to find one root x1 with simple bisection, then divide the polynomial into (x - x1) and deal with remaining quadratic equation.
Il 01/04/20 18:46, Paul Gofman ha scritto:
Given the complex roots are not needed here and the polynomial is always cubic, is this generic method really beneficial? It would probably be simpler and quicker to find one root x1 with simple bisection, then divide the polynomial into (x - x1) and deal with remaining quadratic equation.
This kind of division is typically numerically unstable. It might be that for cubic polynomials the problem is not very apparent, but given that Aberth method doesn't look that difficult to me to implement, I would go for that one. Also, I think that the Aberth method, being a variation on Newton's method, is much quicker than bisection.
That said, I don't have much time to write code right now, so I won't get in the way of those who are. Just a suggestion.
Giovanni.
On 4/1/20 20:03, Giovanni Mascellani wrote:
Il 01/04/20 18:46, Paul Gofman ha scritto:
Given the complex roots are not needed here and the polynomial is always cubic, is this generic method really beneficial? It would probably be simpler and quicker to find one root x1 with simple bisection, then divide the polynomial into (x - x1) and deal with remaining quadratic equation.
This kind of division is typically numerically unstable. It might be that for cubic polynomials the problem is not very apparent,
Yes, factoring out the roots from a high degree polynomial can accumulate the error, but how's that a problem for just one root?
Also, I think just using double precision in analytical solution will avoid any practical stability problems in this case.
I will see if I can get an Alberth method function written without too much difficulty and test it out, maybe run some speed tests between the two.
Also, I guess I should switch over to doubles then as well. I'll see if that makes any differences.
On Wed, Apr 1, 2020 at 1:14 PM Paul Gofman gofmanp@gmail.com wrote:
On 4/1/20 20:03, Giovanni Mascellani wrote:
Il 01/04/20 18:46, Paul Gofman ha scritto:
Given the complex roots are not needed here and the polynomial is always cubic, is this generic method really beneficial? It would probably be simpler and quicker to find one root x1 with simple bisection, then divide the polynomial into (x - x1) and deal with remaining quadratic equation.
This kind of division is typically numerically unstable. It might be that for cubic polynomials the problem is not very apparent,
Yes, factoring out the roots from a high degree polynomial can accumulate the error, but how's that a problem for just one root?
Also, I think just using double precision in analytical solution will avoid any practical stability problems in this case.
On 4/2/20 19:42, Connor McAdams wrote:
I will see if I can get an Alberth method function written without too much difficulty and test it out, maybe run some speed tests between the two.
I still can't understand the motivation why would you pursue a more complicated and general method which is not required for the problem, has more work to do by finding all the complex roots and having special requirements for initial value estimation. Analytical solution is probably good enough, even using doubles there is probably an overshoot. If for any reason you don't like it one root can be found by bisection. That is especially handy given you are interested in roots in 0.0 - 1.0 interval only, so you immediately know how to bound your search instead of caring about correct initial estimation for a generic iterative method). Then, given the polynomial is P_3(x) = a*x^3 + b*x^2 + c*x + d, and x1 is root (that is, P_3(x1) = 0), factoring this root out gives:
P_2(x) = a*x^2 + (b + a*x1)*x + (c + b*x1 + a*x1^2)*x.
What do you expect to go wrong with these quadratic polynomial coefficients accuracy? Especially given you are interested in a root in [0-1] interval only? If you could expect some special polynomial coefficients there that would mean that you would have troubles just evaluating such a polynomial without searching for any roots and would have to reformulate the problem in some way, but I doubt this is a real concern.
I suppose it would be more practical to make some tests for this, from higher level (testing how the intersection works, including degraded and corner cases like vertical lines, horizontal lines etc.) to the inner root finding solution function (given it is very easy to test by substituting the found roots in the initial polynomial).
Also things like this:
+ theta = -atan2f(q[1]->y - q[0]->y, q[1]->x - q[0]->x); + /* Rotate the Y coordinates of the cubic bezier. */ + y[0] = (p[0]->x - q[0]->x) * sinf(theta) + (p[0]->y - q[0]->y) * cosf(theta); + y[1] = (p[1]->x - q[0]->x) * sinf(theta) + (p[1]->y - q[0]->y) * cosf(theta); + y[2] = (p[2]->x - q[0]->x) * sinf(theta) + (p[2]->y - q[0]->y) * cosf(theta); + y[3] = (p[3]->x - q[0]->x) * sinf(theta) + (p[3]->y - q[0]->y) * cosf(theta);
Apart from computing sin and cos 4 times each, you can avoid using any transcendental function by normalizing (q[1]->y - q[0]->y, q[1]->x - q[0]->x) vector and just using its components instead of sin, cos.
Also, I guess I should switch over to doubles then as well. I'll see if that makes any differences.
On Wed, Apr 1, 2020 at 1:14 PM Paul Gofman gofmanp@gmail.com wrote:
On 4/1/20 20:03, Giovanni Mascellani wrote:
Il 01/04/20 18:46, Paul Gofman ha scritto:
Given the complex roots are not needed here and the polynomial is always cubic, is this generic method really beneficial? It would probably be simpler and quicker to find one root x1 with simple bisection, then divide the polynomial into (x - x1) and deal with remaining quadratic equation.
This kind of division is typically numerically unstable. It might be that for cubic polynomials the problem is not very apparent,
Yes, factoring out the roots from a high degree polynomial can accumulate the error, but how's that a problem for just one root?
Also, I think just using double precision in analytical solution will avoid any practical stability problems in this case.
-
Connor McAdams -
Giovanni Mascellani -
Paul Gofman