Bézier curves/Intersections: Difference between revisions
Content added Content deleted
(→C implementation 2: Added gnu::const attributes to try to help an optimizer.) |
(Added Icon.) |
||
Line 1,215: | Line 1,215: | ||
(-0.681025, 2.681025) |
(-0.681025, 2.681025) |
||
(-0.854983, 1.345017) |
(-0.854983, 1.345017) |
||
</pre> |
|||
=={{header|Icon}}== |
|||
This implementation combines rectangle overlap (as in the Modula-2) and curve flattening (as in one of the C implementations), and tries to arrange the calculations efficiently. I think it may be a new algorithm worth serious consideration. Curve flattening alone requires too many line-segment-intersection checks, but rectangle overlap "prunes the tree". |
|||
Furthermore, the algorithm returns <math>t</math>-parameter pairs, which is ideal. |
|||
(Icon is a very good language in which to express this algorithm. The problem is few people will know how to read it. :D ) |
|||
<syntaxhighlight lang="icon"> |
|||
# This program combines the methods of the 2nd C implementation (which |
|||
# by itself is inefficient) with those of the Modula-2 implementation, |
|||
# and then rearranges the computations to try to achieve greater |
|||
# efficiency. |
|||
# |
|||
# The algorithm actually returns t-parameters for two curves, as a |
|||
# pair for each intersection point. This is exactly what one might |
|||
# want: for instance, to break a font glyph, made from two or more |
|||
# other glyphs, into pieces at the points of intersection of all the |
|||
# outlines. |
|||
# |
|||
# The code below is written to illustrate the algorithm rather than to |
|||
# squeeze performance out of Icon. For instance, I use a "set" to |
|||
# store the workload, and, when choosing the next workitem-pair to |
|||
# work on, do so by random selection. It would be faster, certainly, |
|||
# to use an Icon "list", as either a stack or a queue or both. |
|||
# |
|||
# It is also possible, of course, to cast the algorithm as a recursive |
|||
# procedure. |
|||
# |
|||
# |
|||
# References on the s-power basis: |
|||
# |
|||
# J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power |
|||
# basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997, |
|||
# page 319. |
|||
# |
|||
# J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis |
|||
# in geometry processing’, ACM Transactions on Graphics, vol 19 |
|||
# no 1, January 2000, page 35. |
|||
# |
|||
record point (x, y) |
|||
record spower (c0, c1, c2) |
|||
record curve (x, y) |
|||
record workitem (P, t0, t1, pt0, pt1) |
|||
$define P_controls [point (-1, 0), point ( 0, 10), point ( 1, 0)] |
|||
$define Q_controls [point ( 2, 1), point (-8, 2), point ( 2, 3)] |
|||
$define DEFAULT_NUMPIECES 2 # Bisection. |
|||
# Tolerance of the ratio of a bound on the non-linear component to the |
|||
# length of the segment. I use a max norm but you can use your |
|||
# favorite norm. |
|||
$define DEFAULT_FLATNESS_TOLERANCE 0.0001 |
|||
# For demonstration, I choose a minimum spacing between intersection |
|||
# points equal to several times single precision machine epsilon. I |
|||
# measure distance using a max norm, but you can use your favorite |
|||
# norm. |
|||
$define DEFAULT_MINIMUM_SPACING 1e-6 |
|||
procedure main () |
|||
local P, Q |
|||
local intersections, xsect |
|||
P := controls_to_curve ! P_controls |
|||
Q := controls_to_curve ! Q_controls |
|||
intersections := find_intersections (P, Q) |
|||
write () |
|||
write_tabbed_line (" convex up\t" || |
|||
" convex left\t") |
|||
every xsect := !intersections do |
|||
write_tabbed_line (xsect[1] || "\t(" || |
|||
xsect[2].x || ", " || xsect[2].y || ")\t" || |
|||
xsect[3] || "\t(" || |
|||
xsect[4].x || ", " || xsect[4].y || ")") |
|||
write () |
|||
end |
|||
procedure write_tabbed_line (line) |
|||
write (detab (line, 18, 56, 74)) |
|||
end |
|||
procedure find_intersections (P, Q, tol, spacing) |
|||
# Return a list of [tP, pointP, tQ, pointQ] for the intersections, |
|||
# sorted by tP values. |
|||
local workload, ts |
|||
local tP, ptP |
|||
local tQ, ptQ |
|||
local tbl, intersections |
|||
/tol := DEFAULT_FLATNESS_TOLERANCE |
|||
/spacing := DEFAULT_MINIMUM_SPACING |
|||
workload := initial_workload (P, Q) |
|||
tbl := table () |
|||
every ts := process_workload (tol, workload) do |
|||
{ |
|||
tP := ts[1]; ptP := curve_eval (P, tP) |
|||
tQ := ts[2]; ptQ := curve_eval (Q, tQ) |
|||
not (max (abs ((!tbl)[2].x - ptP.x), |
|||
abs ((!tbl)[2].y - ptP.y)) < spacing) & |
|||
not (max (abs ((!tbl)[4].x - ptQ.x), |
|||
abs ((!tbl)[4].y - ptQ.y)) < spacing) & |
|||
tbl[tP] := [tP, ptP, tQ, ptQ] |
|||
} |
|||
tbl := sort (tbl, 1) |
|||
every put (intersections := [], (!tbl)[2]) |
|||
return intersections |
|||
end |
|||
procedure process_workload (tol, workload) |
|||
# Generate pairs of t-parameters. |
|||
local pair, ts |
|||
while *workload ~= 0 do |
|||
{ |
|||
pair := ?workload |
|||
delete (workload, pair) |
|||
if rectangles_overlap (pair[1].pt0, pair[1].pt1, |
|||
pair[2].pt0, pair[2].pt1) then |
|||
{ |
|||
if flat_enough (tol, pair[1]) then |
|||
{ |
|||
if flat_enough (tol, pair[2]) then |
|||
{ |
|||
if ts := segment_parameters (pair[1].pt0, pair[1].pt1, |
|||
pair[2].pt0, pair[2].pt1) then |
|||
suspend [(1 - ts[1]) * pair[1].t0 + ts[1] * pair[1].t1, |
|||
(1 - ts[2]) * pair[2].t0 + ts[2] * pair[2].t1] |
|||
} |
|||
else |
|||
every insert (workload, [pair[1], |
|||
split_workitem (pair[2])]) |
|||
} |
|||
else |
|||
{ |
|||
if flat_enough (tol, pair[2]) then |
|||
every insert (workload, [split_workitem (pair[1]), |
|||
pair[2]]) |
|||
else |
|||
every insert (workload, [split_workitem (pair[1]), |
|||
split_workitem (pair[2])]) |
|||
} |
|||
} |
|||
} |
|||
end |
|||
procedure split_workitem (W, num_pieces) |
|||
# Split a workitem in pieces and generate the pieces. |
|||
local fraction, t1_t0, ts, pts, i |
|||
/num_pieces := DEFAULT_NUMPIECES |
|||
fraction := 1.0 / num_pieces |
|||
t1_t0 := W.t1 - W.t0 |
|||
every put (ts := [], |
|||
W.t0 + (1 to num_pieces - 1) * fraction * t1_t0) |
|||
every put (pts := [], curve_eval (W.P, !ts)) |
|||
ts := [W.t0] ||| ts ||| [W.t1] |
|||
pts := [W.pt0] ||| pts ||| [W.pt1] |
|||
every i := 1 to *pts - 1 do |
|||
suspend (workitem (W.P, ts[i], ts[i + 1], pts[i], pts[i + 1])) |
|||
end |
|||
procedure initial_workload (P, Q) |
|||
# Create an initial workload set, by breaking P and Q at any extreme |
|||
# points. |
|||
local workload |
|||
every insert (workload := set (), [break_at_extreme_points (P), |
|||
break_at_extreme_points (Q)]) |
|||
return workload |
|||
end |
|||
procedure break_at_extreme_points (P) |
|||
# Generate workitems for the curve P, after breaking it at any |
|||
# extreme points. |
|||
local ts, pts, i |
|||
ts := [0] ||| sort (curve_extreme_points (P)) ||| [1] |
|||
every put (pts := [], curve_eval (P, !ts)) |
|||
every i := 1 to *pts - 1 do |
|||
suspend (workitem (P, ts[i], ts[i + 1], pts[i], pts[i + 1])) |
|||
end |
|||
procedure flat_enough (tol, P, t0, t1, pt0, pt1) |
|||
# Is the [t0,t1] portion of the curve P flat enough to be treated as |
|||
# a straight line between pt0 and pt1, where pt0 and pt1 are the |
|||
# endpoints of the portion? |
|||
local error, length |
|||
# Let flat_enough be called this way, where W is a workitem: |
|||
# flat_enough(tol,W) |
|||
if /t0 then |
|||
{ |
|||
pt1 := P.pt1 |
|||
pt0 := P.pt0 |
|||
t1 := P.t1 |
|||
t0 := P.t0 |
|||
P := P.P |
|||
} |
|||
# pt0 and pt1 probably have been computed before and saved, but if |
|||
# necessary they could be computed now: |
|||
/pt0 := curve_eval (P, t0) |
|||
/pt1 := curve_eval (P, t1) |
|||
# The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to |
|||
# remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached |
|||
# at t=1/2. That accounts for the 1/4=0.25 in the following, which |
|||
# uses "max norm" length measurements. (Substitute your favorite |
|||
# norm.) |
|||
error := 0.25 * max (abs (spower_center_coef (P.x, t0, t1)), |
|||
abs (spower_center_coef (P.y, t0, t1))) |
|||
length := max (abs (pt1.x - pt0.x), abs (pt1.y - pt0.y)) |
|||
((error <= length * tol) & return) | fail |
|||
end |
|||
procedure curve_eval (P, t) |
|||
# Return the point that lies on the curve P at parameter value t. |
|||
return point (spower_eval (P.x, t), spower_eval (P.y, t)) |
|||
end |
|||
procedure curve_extreme_points (P) |
|||
# Return a set containing parameter values (values of t) for the |
|||
# extreme points of curve P. |
|||
local ts |
|||
ts := set () |
|||
insert (ts, spower_extreme_point (P.x)) |
|||
insert (ts, spower_extreme_point (P.y)) |
|||
return ts |
|||
end |
|||
procedure spower_eval (p, t) |
|||
# Evaluate the s-power spline p at t. |
|||
return (p.c0 + (p.c1 * t)) * (1 - t) + (p.c2 * t) |
|||
end |
|||
procedure spower_center_coef (p, t0, t1) |
|||
# Return the center coefficient for the [t0,t1] portion of the |
|||
# s-power spline p. |
|||
if /t1 then { t1 := t0[2]; t0 := t0[1] } # Allow a list as t0. |
|||
return p.c1 * ((t1 - t0 - t0) * t1 + (t0 * t0)) |
|||
end |
|||
procedure spower_extreme_point (p) |
|||
# Return t in [0,1] where p is extreme, else fail. |
|||
local t |
|||
p.c1 = 0 & fail # The spline is linear |
|||
p.c0 = p.c2 & return 0.5 # The spline is "pulse-like". |
|||
t := (0.5 * (p.c2 + p.c1 - p.c0)) / p.c1 # Root of the derivative. |
|||
0 <= t <= 1 & return t |
|||
fail |
|||
end |
|||
procedure rectangles_overlap (ptA0, ptA1, ptB0, ptB1) |
|||
# Do the rectangles with corners at (ptA0,ptA1) and (ptB0,ptB1) |
|||
# overlap? |
|||
local ax0, ay0, ax1, ay1 |
|||
local bx0, by0, bx1, by1 |
|||
ax0 := ptA0.x; ax1 := ptA1.x |
|||
bx0 := ptB0.x; bx1 := ptB1.x |
|||
if ax1 < ax0 then ax0 :=: ax1 |
|||
if bx1 < bx0 then bx0 :=: bx1 |
|||
bx1 < ax0 & fail |
|||
ax1 < bx0 & fail |
|||
ay0 := ptA0.y; ay1 := ptA1.y |
|||
by0 := ptB0.y; by1 := ptB1.y |
|||
if ay1 < ay0 then ay0 :=: ay1 |
|||
if by1 < by0 then by0 :=: by1 |
|||
by1 < ay0 & fail |
|||
ay1 < by0 & fail |
|||
return |
|||
end |
|||
procedure segment_parameters (ptA0, ptA1, ptB0, ptB1) |
|||
# Return the respective [0,1] parameters of line segments |
|||
# (ptA0,ptA1) and (ptB0,ptB1), for their intersection point. Fail if |
|||
# there are not such parameters. |
|||
local ax0, ax1, ay0, ay1 |
|||
local bx0, bx1, by0, by1 |
|||
local ax1_ax0, ay1_ay0 |
|||
local bx1_bx0, by1_by0 |
|||
local anumer, bnumer, denom |
|||
local tA, tB, x, y |
|||
ax0 := ptA0.x; ax1 := ptA1.x |
|||
ay0 := ptA0.y; ay1 := ptA1.y |
|||
bx0 := ptB0.x; bx1 := ptB1.x |
|||
by0 := ptB0.y; by1 := ptB1.y |
|||
ax1_ax0 := ax1 - ax0 |
|||
ay1_ay0 := ay1 - ay0 |
|||
bx1_bx0 := bx1 - bx0 |
|||
by1_by0 := by1 - by0 |
|||
denom := (ax1_ax0 * by1_by0) - (ay1_ay0 * bx1_bx0) |
|||
anumer := |
|||
(bx1_bx0 * ay0) - (by1_by0 * ax0) + (bx0 * by1) - (bx1 * by0) |
|||
tA := anumer / denom |
|||
0 <= tA <= 1 | fail |
|||
bnumer := |
|||
-((ax1_ax0 * by0) - (ay1_ay0 * bx0) + (ax0 * ay1) - (ax1 * ay0)) |
|||
tB := bnumer / denom |
|||
0 <= tB <= 1 | fail |
|||
return [tA, tB] |
|||
end |
|||
procedure controls_to_curve (p0, p1, p2) |
|||
# Convert control points to a curve in s-power basis. |
|||
return curve (spower (p0.x, (2 * p1.x) - p0.x - p2.x, p2.x), |
|||
spower (p0.y, (2 * p1.y) - p0.y - p2.y, p2.y)) |
|||
end |
|||
procedure abs (x) |
|||
return (if x < 0 then -x else x) |
|||
end |
|||
procedure max (x, y) |
|||
return (if x < y then y else x) |
|||
end |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
For each estimated intersection point, the program prints out the <math>t</math>-parameters and the corresponding values of the curves. |
|||
<pre> |
|||
convex up convex left |
|||
0.07250828117 (-0.8549834377, 1.345016607) 0.1725082997 (-0.8549837251, 1.345016599) |
|||
0.1594875309 (-0.6810249382, 2.681025168) 0.8405124691 (-0.681025168, 2.681024938) |
|||
0.8274917003 (0.6549834005, 2.854983725) 0.9274917188 (0.6549833933, 2.854983438) |
|||
0.9405124667 (0.8810249334, 1.118975334) 0.05948753331 (0.8810246662, 1.118975067) |
|||
</pre> |
</pre> |
||