Bézier curves/Intersections: Difference between revisions
Content added Content deleted
(Added Object Icon.) |
|||
Line 2,040: | Line 2,040: | ||
(-0.6810249, 2.68102500) |
(-0.6810249, 2.68102500) |
||
(-0.8549834, 1.34501657)</pre> |
(-0.8549834, 1.34501657)</pre> |
||
=={{header|ObjectIcon}}== |
|||
{{trans|Python}} |
|||
{{trans|Icon}} |
|||
See also: |
|||
* [[#Pascal|Pascal]] |
|||
This program is mostly a translation of the Python. |
|||
Note that that the Icon program could itself be run as an Object Icon program, after a few necessary modifications. ("import io" would have to be added, the implementations of "abs" and "max" would have to be removed, and maybe a few other minor changes would be needed.) |
|||
<syntaxhighlight lang="objecticon"> |
|||
# This is the algorithm of the Icon and Python implementations. |
|||
# Primarily it is a translation of the Python. Snippets are taken from |
|||
# the Icon. |
|||
# |
|||
# References on the symmetric 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. |
|||
# |
|||
import ipl.printf(printf) |
|||
procedure main () |
|||
local flatness_tolerance |
|||
local minimum_spacing |
|||
local p, q |
|||
local tbl, params, tp, tq, ppoint, qpoint |
|||
local keys, k |
|||
flatness_tolerance := 0.0001 |
|||
minimum_spacing := 0.000001 |
|||
p := Curve.from_controls (Point (-1.0, 0.0), |
|||
Point (0.0, 10.0), |
|||
Point (1.0, 0.0)) |
|||
q := Curve.from_controls (Point (2.0, 1.0), |
|||
Point (-8.0, 2.0), |
|||
Point ( 2.0, 3.0)) |
|||
tbl := table () |
|||
every params := find_intersections (p, q, flatness_tolerance) do |
|||
{ |
|||
tp := params[1]; ppoint := p.eval(tp) |
|||
tq := params[2]; qpoint := q.eval(tq) |
|||
not (length ((!tbl)[2].x - ppoint.x, |
|||
(!tbl)[2].y - ppoint.y) < minimum_spacing) & |
|||
not (length ((!tbl)[4].x - qpoint.x, |
|||
(!tbl)[4].y - qpoint.y) < minimum_spacing) & |
|||
tbl[tp] := [tp, ppoint, tq, qpoint] |
|||
} |
|||
every put (keys := [], key (tbl)) |
|||
keys := sort (keys) |
|||
printf ("\n") |
|||
printf (" convex up " || |
|||
" convex left\n") |
|||
every k := !keys do |
|||
{ |
|||
tp := tbl[k][1] |
|||
ppoint := tbl[k][2] |
|||
tq := tbl[k][3] |
|||
qpoint := tbl[k][4] |
|||
printf (" %11.8r (%11.8r, %11.8r) " || |
|||
"%11.8r (%11.8r, %11.8r)\n", |
|||
tp, ppoint.x, ppoint.y, tq, qpoint.x, qpoint.y) |
|||
} |
|||
printf ("\n") |
|||
end |
|||
# Generate t-parameter pairs detected as corresponding to intersection |
|||
# points of p and q. There may be duplicate detections. It is assumed |
|||
# those will be weeded out by later processing. The tol parameter |
|||
# specifies the "flatness tolerance" and is a relative tolerance. |
|||
procedure find_intersections (p, q, tol) |
|||
local i, j |
|||
local tp, tq, pportion, qportion |
|||
local workload, workpair, params |
|||
# The initial workload is the cartesian product of the monotonic |
|||
# portions of p and q, respectively. |
|||
tp := [0.0] ||| sort (p.critical_points()) ||| [1.0] |
|||
tq := [0.0] ||| sort (q.critical_points()) ||| [1.0] |
|||
workload := set () |
|||
every i := 1 to *tp - 1 do |
|||
every j := 1 to *tq - 1 do |
|||
{ |
|||
pportion := Portion (p, tp[i], tp[i + 1], |
|||
p.eval(tp[i]), p.eval(tp[i + 1])) |
|||
qportion := Portion (q, tq[j], tq[j + 1], |
|||
q.eval(tq[j]), q.eval(tq[j + 1])) |
|||
insert (workload, [pportion, qportion]) |
|||
} |
|||
while *workload ~= 0 do |
|||
{ |
|||
workpair := ?workload |
|||
delete (workload, workpair) |
|||
pportion := workpair[1] |
|||
qportion := workpair[2] |
|||
if rectangles_overlap (pportion.endpt0, pportion.endpt1, |
|||
qportion.endpt0, qportion.endpt1) then |
|||
{ |
|||
if pportion.flat_enough(tol) then |
|||
{ |
|||
if qportion.flat_enough(tol) then |
|||
{ |
|||
if params := segment_parameters(pportion.endpt0, |
|||
pportion.endpt1, |
|||
qportion.endpt0, |
|||
qportion.endpt1) then |
|||
{ |
|||
tp := params[1] |
|||
tq := params[2] |
|||
tp := (1 - tp) * pportion.t0 + tp * pportion.t1 |
|||
tq := (1 - tq) * qportion.t0 + tq * qportion.t1 |
|||
suspend [tp, tq] |
|||
} |
|||
} |
|||
else |
|||
every insert (workload, [pportion, qportion.split()]) |
|||
} |
|||
else |
|||
{ |
|||
if qportion.flat_enough(tol) then |
|||
every insert (workload, [pportion.split(), qportion]) |
|||
else |
|||
every insert (workload, [pportion.split(), |
|||
qportion.split()]) |
|||
} |
|||
} |
|||
} |
|||
end |
|||
class Point () |
|||
public const x, y |
|||
public new (x, y) |
|||
self.x := x |
|||
self.y := y |
|||
return |
|||
end |
|||
end |
|||
class SPower () # Non-parametric spline in s-power basis. |
|||
public const c0, c1, c2 |
|||
private critpoints |
|||
public new (c0, c1, c2) |
|||
self.c0 := c0 |
|||
self.c1 := c1 |
|||
self.c2 := c2 |
|||
return |
|||
end |
|||
# Evaluate at t. |
|||
public eval (t) |
|||
return (c0 + (c1 * t)) * (1.0 - t) + (c2 * t) |
|||
end |
|||
# Return the center coefficient for the [t0,t1] portion. (The other |
|||
# coefficients can be found with the eval method.) |
|||
public center_coef (t0, t1) |
|||
return c1 * ((t1 - t0 - t0) * t1 + (t0 * t0)) |
|||
end |
|||
# Return a set of independent variable values for the critical |
|||
# points that lie in (0,1). (This is a memoizing implementation.) |
|||
public critical_points () |
|||
local t |
|||
if /critpoints then |
|||
{ |
|||
critpoints := set () |
|||
if c1 ~= 0.0 then # If c1 is zero then the spline is linear. |
|||
{ |
|||
if c0 = c2 then |
|||
insert (critpoints, 0.5) # The spline is "pulse-like". |
|||
else |
|||
{ |
|||
# The root of the derivative is the critical point. |
|||
t = (0.5 * (c2 + c1 - c0)) / c1 |
|||
insert (critpoints, t) |
|||
} |
|||
} |
|||
} |
|||
return critpoints |
|||
end |
|||
end |
|||
class Curve () # Parametric spline in s-power basis. |
|||
public const x, y |
|||
public new (x, y) |
|||
self.x := x |
|||
self.y := y |
|||
return |
|||
end |
|||
public static from_controls (ctl0, ctl1, ctl2) |
|||
local c0x, c0y, c1x, c1y, c2x, c2y |
|||
c0x := ctl0.x |
|||
c0y := ctl0.y |
|||
c1x := (2.0 * ctl1.x) - ctl0.x - ctl2.x |
|||
c1y := (2.0 * ctl1.y) - ctl0.y - ctl2.y |
|||
c2x := ctl2.x |
|||
c2y := ctl2.y |
|||
return Curve (SPower (c0x, c1x, c2x), |
|||
SPower (c0y, c1y, c2y)) |
|||
end |
|||
# Evaluate at t. |
|||
public eval (t) |
|||
return Point (x.eval(t), y.eval(t)) |
|||
end |
|||
# Return a set of t-parameter values for the critical points that |
|||
# lie in (0,1). |
|||
public critical_points () |
|||
return (x.critical_points() ++ y.critical_points()) |
|||
end |
|||
end |
|||
class Portion () # The [t0,t1]-portion of a parametric spline. |
|||
public const curve, t0, t1, endpt0, endpt1 |
|||
public static default_num_pieces |
|||
private static init () |
|||
default_num_pieces := 2 |
|||
end |
|||
public new (curve, t0, t1, endpt0, endpt1) |
|||
self.curve := curve |
|||
self.t0 := t0 |
|||
self.t1 := t1 |
|||
self.endpt0 := endpt0 |
|||
self.endpt1 := endpt1 |
|||
return |
|||
end |
|||
# Is the Portion close enough to linear to be treated as a line |
|||
# segment? |
|||
public flat_enough (tol) |
|||
local xcentercoef, ycentercoef |
|||
local xlen, ylen |
|||
# 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. |
|||
xcentercoef := curve.x.center_coef(t0, t1) |
|||
ycentercoef := curve.y.center_coef(t0, t1) |
|||
xlen := endpt1.x - endpt0.x |
|||
ylen := endpt1.y - endpt0.y |
|||
return compare_lengths (0.25 * xcentercoef, |
|||
0.25 * ycentercoef, |
|||
tol * xlen, tol * ylen) <= 0 |
|||
end |
|||
# Generate num_pieces (or default_num_pieces) sections of the |
|||
# Portion. |
|||
public split (num_pieces) |
|||
local k1, k, ts, vals, i |
|||
/num_pieces := default_num_pieces |
|||
k1 := 1.0 / num_pieces |
|||
every put (ts := [], (k := k1 * (1 to num_pieces - 1) & |
|||
(1.0 - k) * t0 + k * t1)) |
|||
every put (vals := [], curve.eval(!ts)) |
|||
ts := [t0] ||| ts ||| [t1] |
|||
vals := [endpt0] ||| vals ||| [endpt1] |
|||
every i := 1 to *ts - 1 do |
|||
suspend Portion (curve, ts[i], ts[i + 1], vals[i], vals[i + 1]) |
|||
end |
|||
end |
|||
# Do the rectangles with corners at (ptA0,ptA1) and (ptB0,ptB1) |
|||
# overlap? |
|||
procedure rectangles_overlap (ptA0, ptA1, ptB0, ptB1) |
|||
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 |
|||
# 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. |
|||
procedure segment_parameters (ptA0, ptA1, ptB0, ptB1) |
|||
local ax0, ax1, ay0, ay1 |
|||
local bx0, bx1, by0, by1 |
|||
local axdiff, aydiff |
|||
local bxdiff, bydiff |
|||
local anumer, bnumer, denom |
|||
local tA, tB |
|||
ax0 := ptA0.x; ax1 := ptA1.x |
|||
ay0 := ptA0.y; ay1 := ptA1.y |
|||
bx0 := ptB0.x; bx1 := ptB1.x |
|||
by0 := ptB0.y; by1 := ptB1.y |
|||
axdiff := ax1 - ax0 |
|||
aydiff := ay1 - ay0 |
|||
bxdiff := bx1 - bx0 |
|||
bydiff := by1 - by0 |
|||
denom := (axdiff * bydiff) - (aydiff * bxdiff) |
|||
anumer := |
|||
(bxdiff * ay0) - (bydiff * ax0) + (bx0 * by1) - (bx1 * by0) |
|||
tA := anumer / denom |
|||
0 <= tA <= 1 | fail |
|||
bnumer := |
|||
-((axdiff * by0) - (aydiff * bx0) + (ax0 * ay1) - (ax1 * ay0)) |
|||
tB := bnumer / denom |
|||
0 <= tB <= 1 | fail |
|||
return [tA, tB] |
|||
end |
|||
# Length according to some norm, where (ax,ay) is a "measuring stick" |
|||
# vector. Here I use the max norm. |
|||
procedure length (a_x, a_y) |
|||
return max (abs (a_x), abs (a_y)) |
|||
end |
|||
# Having a compare_lengths function lets one compare lengths in the |
|||
# euclidean metric by comparing the squares of the lengths, and thus |
|||
# avoiding the square root. The following, however, is a general |
|||
# implementation. |
|||
procedure compare_lengths (a_x, a_y, b_x, b_y) |
|||
local len_a, len_b, cmpval |
|||
len_a := length (a_x, a_y) |
|||
len_b := length (b_x, b_y) |
|||
if len_a < len_b then |
|||
cmpval := -1 |
|||
else if len_a > len_b then |
|||
cmpval := 1 |
|||
else |
|||
cmpval := 0 |
|||
return cmpval |
|||
end |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
convex up convex left |
|||
0.07250828 (-0.85498344, 1.34501661) 0.17250830 (-0.85498373, 1.34501660) |
|||
0.15948753 (-0.68102494, 2.68102517) 0.84051247 (-0.68102517, 2.68102494) |
|||
0.82749170 ( 0.65498340, 2.85498373) 0.92749172 ( 0.65498339, 2.85498344) |
|||
0.94051247 ( 0.88102493, 1.11897533) 0.05948753 ( 0.88102467, 1.11897507) |
|||
</pre> |
|||
=={{header|Pascal}}== |
=={{header|Pascal}}== |