Bézier curves/Intersections: Difference between revisions
Added Object Icon.
(Added Object Icon.) |
|||
Line 2,040:
(-0.6810249, 2.68102500)
(-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}}==
| |||