Bézier curves/Intersections: Difference between revisions
Content added Content deleted
(Created Nim solution.) |
(Added Ada.) |
||
Line 12: | Line 12: | ||
* [[Bitmap/Bézier_curves/Quadratic|Bitmap/Bézier curves/Quadratic]] |
* [[Bitmap/Bézier_curves/Quadratic|Bitmap/Bézier curves/Quadratic]] |
||
* [[Steffensen's_method|Steffensen's method]] |
* [[Steffensen's_method|Steffensen's method]] |
||
=={{header|Ada}}== |
|||
This is a variant of the Bézier clipping method of [http://hdl.lib.byu.edu/1877/2822 Thomas W. Sederberg's lecture notes on "Computer Aided Geometric Design"], using symmetric power (s-power) polynomials instead of the usual "control points". I do not know that it is any better than the method of the s-power polynomial flattening, which is done on this page in several languages. The flattening method is considerably simpler. Also, I have doubts about the numerical stability of the clipping method. |
|||
On the other hand, the code below introduces a new way of doing Bézier clipping, by transformation of coordinates. And it also hints at a geometric interpretation of s-power curves, in which all the constituent "points" undergo rotations and scalings, but only the endpoints are subject to translations. |
|||
Ada is an Algol-like language that is case-insensitive. I have used an all-lowercase style with manual indentation. Function parameters may be passed either positionally or by name. Ada does have a few oddities, among which are that the term for a pointer is not "pointer" but "access". |
|||
<syntaxhighlight lang="ada"> |
|||
pragma wide_character_encoding (utf8); |
|||
pragma ada_2022; |
|||
pragma assertion_policy (check); |
|||
-- |
|||
-- The algorithm here is iterated "Bézier clipping" [1] of quadratic |
|||
-- Bézier curves in a plane. However, the curves will be represented |
|||
-- in the s-power basis [2,3], and there will be these other |
|||
-- differences from the algorithms suggested by [1]: |
|||
-- |
|||
-- * Before clipping, the problem will be moved to a coordinate system |
|||
-- in which the "skeleton" of the clipper curve—the line segment |
|||
-- between its endpoints—extends from (0,0) and (1,0). In this |
|||
-- coordinate system, the "distance function" is the y-coordinate. |
|||
-- |
|||
-- * The fat line will be the best possible parallel to the "skeleton" |
|||
-- of the clipper curve. This is easily done in the s-power basis, |
|||
-- once one has transformed coordinate systems as described above: |
|||
-- the center coefficient of the (y,t)-polynomial is exactly four |
|||
-- times the optimum signed fat line width. |
|||
-- |
|||
-- * Intersections with fat line boundaries will be found not by |
|||
-- extending the sides of control polygons, but instead by the |
|||
-- quadratic formula. |
|||
-- |
|||
-- * Mechanisms to ensure termination of the algorithm, and to improve |
|||
-- the chances of finding intersections, may differ. |
|||
-- |
|||
-- Inputs and outputs to the algorithms will be in single |
|||
-- precision. However, the internals will be in double precision. This |
|||
-- decision helps resolve some minor numerical quandaries. |
|||
-- |
|||
-- References: |
|||
-- |
|||
-- [1] Sederberg, Thomas W., "Computer Aided Geometric Design" |
|||
-- (2012). Faculty |
|||
-- Publications. 1. https://scholarsarchive.byu.edu/facpub/1 |
|||
-- http://hdl.lib.byu.edu/1877/2822 |
|||
-- |
|||
-- [2] J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial |
|||
-- power basis’, ACM Transactions on Graphics, vol 16 no 3, July |
|||
-- 1997, page 319. |
|||
-- |
|||
-- [3] 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. |
|||
-- |
|||
-- Suggestions for the future: |
|||
-- |
|||
-- * Combination of Bézier clipping and flatness testing. Remember |
|||
-- from other Rosetta Code examples, that in the s-power basis it is |
|||
-- simple to test the "flatness" of a parametric plane |
|||
-- curve. Suppose two curves are running nearly parallel to each |
|||
-- other, so they are difficult to reduce by clipping. Nevertheless, |
|||
-- they can be made flatter and flatter, so eventually they become |
|||
-- "segment-like" and can be tested by a line-intersection routine. |
|||
-- |
|||
with ada.text_io; |
|||
use ada.text_io; |
|||
with ada.float_text_io; |
|||
use ada.float_text_io; |
|||
with ada.numerics.elementary_functions; |
|||
use ada.numerics.elementary_functions; |
|||
with ada.numerics.long_elementary_functions; |
|||
use ada.numerics.long_elementary_functions; |
|||
with ada.unchecked_deallocation; |
|||
procedure bezier_intersections is |
|||
subtype real is float; |
|||
subtype lreal is long_float; |
|||
function between_0_and_1 (x : real) |
|||
return boolean is |
|||
begin |
|||
return 0.0 <= x and x <= 1.0; |
|||
end between_0_and_1; |
|||
function between_0_and_1 (x : lreal) |
|||
return boolean is |
|||
begin |
|||
return 0.0 <= x and x <= 1.0; |
|||
end between_0_and_1; |
|||
procedure quadratic_formula (a, b, c : lreal; |
|||
has_smaller : out boolean; |
|||
smaller : out lreal; |
|||
has_larger : out boolean; |
|||
larger : out lreal) |
|||
with post => has_smaller or not has_larger |
|||
is |
|||
-- |
|||
-- This routine is based on QUADPL from the Naval Surface Warfare |
|||
-- Center mathematical library nswc.f90 https://archive.ph/2IpWb |
|||
-- |
|||
-- We are interested only in real roots. Be aware that, if there |
|||
-- are two roots, they might be a double root. |
|||
-- |
|||
bh, e, d : lreal; |
|||
-- Avogadro’s number, by definition (since 2019): |
|||
arbitrary_value : constant lreal := 6.012214076e23; |
|||
begin |
|||
has_smaller := false; |
|||
has_larger := false; |
|||
smaller := arbitrary_value; |
|||
larger := arbitrary_value; |
|||
if a = 0.0 then |
|||
if b = 0.0 then |
|||
-- The equation is a constant equation. I do not know what to |
|||
-- do with that. There may be no solution, or there may be |
|||
-- infinitely many solutions. Let us treat it as having no |
|||
-- USEFUL solutions. |
|||
null; |
|||
else |
|||
-- The equation is linear. There truly is only one root. |
|||
has_smaller := true; |
|||
smaller := -c / b; |
|||
end if; |
|||
elsif c = 0.0 then |
|||
-- The equation is quadratic, but the constant term |
|||
-- vanishes. One of the roots is trivially zero. |
|||
has_smaller := true; |
|||
has_larger := true; |
|||
smaller := 0.0; |
|||
larger := -b / a; |
|||
else |
|||
-- Compute the discriminant, avoiding overflow. |
|||
bh := b / 2.0; |
|||
if abs (bh) >= abs (c) then |
|||
e := 1.0 - (a / bh) * (c / bh); |
|||
d := sqrt (abs (e)) * abs(bh); |
|||
else |
|||
e := (if c < 0.0 then -a else a); |
|||
e := bh * (bh / abs (c)) - e; |
|||
d := sqrt (abs (e)) * sqrt (abs (c)); |
|||
end if; |
|||
if e < 0.0 then |
|||
-- The solutions are complex conjugates, but we are interested |
|||
-- only in real solutions. |
|||
null; |
|||
else |
|||
has_smaller := true; |
|||
has_larger := true; |
|||
if bh >= 0.0 then |
|||
d := -d; |
|||
end if; |
|||
larger := (-bh + d) / a; |
|||
if larger = 0.0 then |
|||
smaller := 0.0; -- Obviously a double root. |
|||
else |
|||
smaller := (c / larger) / a; |
|||
end if; |
|||
end if; |
|||
end if; |
|||
end quadratic_formula; |
|||
type point is |
|||
record |
|||
x, y : real; |
|||
end record; |
|||
type control_point_curve is |
|||
record |
|||
pt0, pt1, pt2 : point; |
|||
end record; |
|||
type spower_poly is |
|||
record |
|||
c0, c1, c2 : real; |
|||
end record; |
|||
type spower_curve is |
|||
record |
|||
x, y : spower_poly; |
|||
end record; |
|||
function eval (poly : spower_poly; |
|||
t : real) |
|||
return real is |
|||
c0, c1, c2 : real; |
|||
begin |
|||
c0 := poly.c0; |
|||
c1 := poly.c1; |
|||
c2 := poly.c2; |
|||
return (if t <= 0.5 then |
|||
c0 + t * ((c2 - c0) + ((1.0 - t) * c1)) |
|||
else |
|||
c2 + (1.0 - t) * ((c0 - c2) + (t * c1))); |
|||
end eval; |
|||
function eval (curve : spower_curve; |
|||
t : real) |
|||
return point is |
|||
begin |
|||
return (x => eval (curve.x, t), |
|||
y => eval (curve.y, t)); |
|||
end eval; |
|||
function to_spower_curve (ctrl : control_point_curve) |
|||
return spower_curve is |
|||
curve : spower_curve; |
|||
p0x, p0y : long_float; |
|||
p1x, p1y : long_float; |
|||
p2x, p2y : long_float; |
|||
begin |
|||
curve.x.c0 := ctrl.pt0.x; |
|||
curve.y.c0 := ctrl.pt0.y; |
|||
curve.x.c2 := ctrl.pt2.x; |
|||
curve.y.c2 := ctrl.pt2.y; |
|||
p0x := long_float (ctrl.pt0.x); |
|||
p0y := long_float (ctrl.pt0.y); |
|||
p1x := long_float (ctrl.pt1.x); |
|||
p1y := long_float (ctrl.pt1.y); |
|||
p2x := long_float (ctrl.pt2.x); |
|||
p2y := long_float (ctrl.pt2.y); |
|||
curve.x.c1 := float (p1x + p1x - p0x - p2x); |
|||
curve.y.c1 := float (p1y + p1y - p0y - p2y); |
|||
return curve; |
|||
end to_spower_curve; |
|||
type opoint is |
|||
record |
|||
-- A point relative to the origin. |
|||
x, y : lreal; |
|||
end record; |
|||
type dpoint is |
|||
record |
|||
-- A point relative to some other point. |
|||
dx, dy : lreal; |
|||
end record; |
|||
function "-" (pt1, pt2 : opoint) |
|||
return dpoint is |
|||
begin |
|||
return (dx => pt1.x - pt2.x, |
|||
dy => pt1.y - pt2.y); |
|||
end "-"; |
|||
function "+" (opt : opoint; |
|||
dpt : dpoint) |
|||
return opoint is |
|||
begin |
|||
return (x => opt.x + dpt.dx, |
|||
y => opt.y + dpt.dy); |
|||
end "+"; |
|||
function "+" (pt1, pt2 : dpoint) |
|||
return dpoint is |
|||
begin |
|||
return (dx => pt1.dx + pt2.dx, |
|||
dy => pt1.dy + pt2.dy); |
|||
end "+"; |
|||
function "*" (scalar : lreal; |
|||
pt : dpoint) |
|||
return dpoint is |
|||
begin |
|||
return (dx => scalar * pt.dx, |
|||
dy => scalar * pt.dy); |
|||
end "*"; |
|||
type geocurve is |
|||
record |
|||
-- |
|||
-- A plane curve in s-power basis, represented geometrically. |
|||
-- |
|||
-- pt0 and pt2 form the curve’s ‘skeleton’. Their convex combination, |
|||
-- ((1.0-t)*pt0)+(t*pt1), form the linear backbone of the curve. |
|||
-- |
|||
-- (1.0-t)*t*pt1 forms (x,t) and (y,t) parabolas that are |
|||
-- convex-up/down, and which (because this is an s-power curve) |
|||
-- vanish as the curve becomes linear. |
|||
-- |
|||
-- The actual geocurve is the locus of (1.0-t)*t*pt1 relative to |
|||
-- the skeleton, respectively for each value of t in [0,1]. |
|||
-- |
|||
-- If [t0,t1] is not [0.0,1.0], then this is a portion of some |
|||
-- parent curve that is determined by context. |
|||
-- |
|||
-- (I feel as if all this can be crafted nicely in terms of a |
|||
-- geometric algebra, with special points not only at infinity |
|||
-- but also at the origin, but let us not worry about that for |
|||
-- now. Likely the calculations would not change much, but the |
|||
-- names of things might.) |
|||
-- |
|||
pt0 : opoint; |
|||
pt1 : dpoint; |
|||
pt2 : opoint; |
|||
t0, t1 : lreal; |
|||
end record; |
|||
function to_geocurve (scurve : spower_curve) |
|||
return geocurve is |
|||
gcurve : geocurve; |
|||
begin |
|||
gcurve.pt0.x := lreal (scurve.x.c0); |
|||
gcurve.pt0.y := lreal (scurve.y.c0); |
|||
gcurve.pt1.dx := lreal (scurve.x.c1); |
|||
gcurve.pt1.dy := lreal (scurve.y.c1); |
|||
gcurve.pt2.x := lreal (scurve.x.c2); |
|||
gcurve.pt2.y := lreal (scurve.y.c2); |
|||
gcurve.t0 := 0.0; |
|||
gcurve.t1 := 1.0; |
|||
return gcurve; |
|||
end to_geocurve; |
|||
function eval (gcurve : geocurve; |
|||
t : lreal) |
|||
return opoint is |
|||
pt0 : opoint; |
|||
pt1 : dpoint; |
|||
pt2 : opoint; |
|||
begin |
|||
pt0 := gcurve.pt0; |
|||
pt1 := gcurve.pt1; |
|||
pt2 := gcurve.pt2; |
|||
return (if t <= 0.5 then |
|||
pt0 + t * ((pt2 - pt0) + ((1.0 - t) * pt1)) |
|||
else |
|||
pt2 + (1.0 - t) * ((pt0 - pt2) + (t * pt1))); |
|||
end eval; |
|||
type affine_projection is |
|||
record |
|||
-- |
|||
-- An affine transformation used to project objects into the |
|||
-- coordinate system of what we are calling a "skeleton". |
|||
-- |
|||
-- If it is an opoint, first add transpose[dx,dy]. If it is a |
|||
-- dpoint, ignore dx,dy. Then, in either case, premultiply by |
|||
-- [[a,b],[c,d]]. (One can view dx,dy as a dpoint.) |
|||
-- |
|||
dx, dy : lreal; |
|||
a, b : lreal; |
|||
c, d : lreal; |
|||
end record; |
|||
function project (pt : opoint; |
|||
trans : affine_projection) |
|||
return opoint is |
|||
x, y : lreal; |
|||
begin |
|||
-- Apply an affine transformation to an opoint, projecting it from |
|||
-- one coordinate system to another. |
|||
x := pt.x + trans.dx; -- The origin may have moved. |
|||
y := pt.y + trans.dy; |
|||
return (x => trans.a * x + trans.b * y, -- Rotation, scaling, etc. |
|||
y => trans.c * x + trans.d * y); |
|||
end project; |
|||
function project (pt : dpoint; |
|||
trans : affine_projection) |
|||
return dpoint is |
|||
x, y : lreal; |
|||
begin |
|||
-- Apply an affine transformation to a dpoint, projecting it from |
|||
-- one coordinate system to another. |
|||
x := pt.dx; -- With dpoint, translation is not performed. |
|||
y := pt.dy; |
|||
return (dx => trans.a * x + trans.b * y, -- Rotation, scaling, &c. |
|||
dy => trans.c * x + trans.d * y); |
|||
end project; |
|||
function project (gcurve : geocurve; |
|||
trans : affine_projection) |
|||
return geocurve is |
|||
begin |
|||
-- Apply an affine transformation to a geocurve, projecting it |
|||
-- from one coordinate system to another. |
|||
return (pt0 => project (gcurve.pt0, trans), |
|||
pt1 => project (gcurve.pt1, trans), |
|||
pt2 => project (gcurve.pt2, trans), |
|||
t0 => 0.0, t1 => 1.0); |
|||
end project; |
|||
function segment_projection_to_0_1 (pt0, pt1 : opoint) |
|||
return affine_projection is |
|||
a, b : lreal; |
|||
denom : lreal; |
|||
begin |
|||
-- Return the transformation that projects pt0 to (0,0) and pt1 to |
|||
-- (1,0) without any distortions except translation, rotation, and |
|||
-- scaling. |
|||
a := pt1.x - pt0.x; |
|||
b := pt1.y - pt0.y; |
|||
denom := (a * a) + (b * b); |
|||
pragma assert (denom /= 0.0); |
|||
a := a / denom; |
|||
b := b / denom; |
|||
return (dx => -pt0.x, dy => -pt0.y, |
|||
a => a, b => b, |
|||
c => -b, d => a); |
|||
end segment_projection_to_0_1; |
|||
procedure test_segment_projection_to_0_1 is |
|||
trans : affine_projection; |
|||
pt : opoint; |
|||
begin |
|||
-- Some unit testing. |
|||
trans := segment_projection_to_0_1 ((1.0, 2.0), (-3.0, 5.0)); |
|||
pt := project ((1.0, 2.0), trans); |
|||
pragma assert (abs (pt.x) < 1.0e-10); |
|||
pragma assert (abs (pt.y) < 1.0e-10); |
|||
pt := project ((-3.0, 5.0), trans); |
|||
pragma assert (abs (pt.x - 1.0) < 1.0e-10); |
|||
pragma assert (abs (pt.y) < 1.0e-10); |
|||
trans := segment_projection_to_0_1 ((0.0, 2.0), (0.0, -50.0)); |
|||
pt := project ((0.0, 2.0), trans); |
|||
pragma assert (abs (pt.x) < 1.0e-10); |
|||
pragma assert (abs (pt.y) < 1.0e-10); |
|||
pt := project ((0.0, -50.0), trans); |
|||
pragma assert (abs (pt.x - 1.0) < 1.0e-10); |
|||
pragma assert (abs (pt.y) < 1.0e-10); |
|||
end test_segment_projection_to_0_1; |
|||
function skeleton_projection (clipper : geocurve) |
|||
return affine_projection is |
|||
begin |
|||
-- Return the transformation that projects the "skeleton" of the |
|||
-- clipper into a coordinate system where the skeleton goes from |
|||
-- the origin (0,0) to the point (1,0). |
|||
return segment_projection_to_0_1 (clipper.pt0, clipper.pt2); |
|||
end skeleton_projection; |
|||
procedure fat_line_bounds (projected_clipper : geocurve; |
|||
dmin, dmax : out lreal) is |
|||
d : lreal; |
|||
d1, d2 : lreal; |
|||
padding : constant lreal := lreal (real'model_epsilon) / 1000.0; |
|||
begin |
|||
-- |
|||
-- (1-t)*t is bounded on [0,1] by 0 and 1/4. Thus we have our fat |
|||
-- line boundaries. We leave some padding, however, for numerical |
|||
-- safety. The padding is much smaller than our tolerance. (We are |
|||
-- calculating in double precision, but to obtain results in |
|||
-- single precision.) |
|||
-- |
|||
-- The fat line bounds will be y=dmin and y=dmax. If you ignore |
|||
-- the padding, the fat line is the best possible among those |
|||
-- parallel to the skeleton. It is likely to be much tighter than |
|||
-- what you would infer from the Bernstein-basis control polygon. |
|||
-- |
|||
d := projected_clipper.pt1.dy / 4.0; |
|||
d1 := -padding * d; |
|||
d2 := (1.0 + padding) * d; |
|||
dmin := lreal'min (d1, d2); |
|||
dmax := lreal'max (d1, d2); |
|||
end fat_line_bounds; |
|||
type clipping_places is |
|||
(clip_everywhere, |
|||
clip_nowhere, |
|||
clip_left_and_right); |
|||
type clipping_plan is |
|||
record |
|||
where : clipping_places; |
|||
t_lft : lreal; |
|||
t_rgt : lreal; |
|||
end record; |
|||
function make_clipping_plan (projected_clippehend : geocurve; |
|||
dmin, dmax : lreal) |
|||
return clipping_plan is |
|||
ncrossings : integer range 0 .. 4; |
|||
tcrossings : array (1 .. 4) of lreal; |
|||
dcrossings : array (1 .. 4) of lreal; |
|||
has_smaller : boolean; |
|||
has_larger : boolean; |
|||
smaller : lreal; |
|||
larger : lreal; |
|||
c0, c1, c2 : lreal; |
|||
function outside_of_bounds |
|||
return boolean is |
|||
dy, ymin, ymax : lreal; |
|||
begin |
|||
-- (1-t)*t is bounded on [0,1] by 0 and 1/4. Use that fact to |
|||
-- make padding for a simple (although not especially tight) |
|||
-- bounds on y. |
|||
dy := abs (projected_clippehend.pt1.dy) / 4.0; |
|||
ymin := -dy + lreal'min (projected_clippehend.pt0.y, |
|||
projected_clippehend.pt2.y); |
|||
ymax := dy + lreal'max (projected_clippehend.pt0.y, |
|||
projected_clippehend.pt2.y); |
|||
return (ymax <= dmin or dmax <= ymin); |
|||
end outside_of_bounds; |
|||
procedure find_crossings (yvalue : lreal) is |
|||
begin |
|||
quadratic_formula (-c1, c1 + c2 - c0, c0 - yvalue, |
|||
has_smaller, smaller, has_larger, larger); |
|||
if has_smaller then |
|||
if between_0_and_1 (smaller) then |
|||
ncrossings := ncrossings + 1; |
|||
tcrossings(ncrossings) := smaller; |
|||
dcrossings(ncrossings) := yvalue; |
|||
end if; |
|||
end if; |
|||
if has_larger then |
|||
if between_0_and_1 (larger) then |
|||
ncrossings := ncrossings + 1; |
|||
tcrossings(ncrossings) := larger; |
|||
dcrossings(ncrossings) := yvalue; |
|||
end if; |
|||
end if; |
|||
end find_crossings; |
|||
function create_a_clipping_plan |
|||
return clipping_plan is |
|||
plan : clipping_plan; |
|||
begin |
|||
-- |
|||
-- For simplicity, there is just one situation in which we will |
|||
-- clip: when there is one crossing of dmin and one crossing of |
|||
-- dmax. |
|||
-- |
|||
if ncrossings = 2 and then dcrossings(1) /= dcrossings(2) then |
|||
plan.where := clip_left_and_right; |
|||
plan.t_lft := lreal'min (tcrossings(1), tcrossings(2)); |
|||
plan.t_rgt := lreal'max (tcrossings(1), tcrossings(2)); |
|||
else |
|||
plan.where := clip_nowhere; |
|||
-- Arbitrary values: |
|||
plan.t_lft := 0.0; |
|||
plan.t_rgt := 0.0; |
|||
end if; |
|||
return plan; |
|||
end create_a_clipping_plan; |
|||
function analyze_and_plan |
|||
return clipping_plan is |
|||
begin |
|||
c0 := projected_clippehend.pt0.y; |
|||
c1 := projected_clippehend.pt1.dy; |
|||
c2 := projected_clippehend.pt2.y; |
|||
ncrossings := 0; |
|||
find_crossings (dmin); |
|||
find_crossings (dmax); |
|||
return create_a_clipping_plan; |
|||
end analyze_and_plan; |
|||
begin |
|||
return (if outside_of_bounds then |
|||
(where => clip_everywhere, |
|||
-- Arbitrary values: |
|||
t_lft => 0.0, t_rgt => 0.0) |
|||
else |
|||
analyze_and_plan); |
|||
end make_clipping_plan; |
|||
function make_portion (gcurve : geocurve; |
|||
t0, t1 : lreal) |
|||
return geocurve is |
|||
pt0 : opoint; |
|||
pt1 : dpoint; |
|||
pt2 : opoint; |
|||
begin |
|||
-- |
|||
-- NOTE FOR THE FUTURE: This operation is a change of HOMOGENEOUS |
|||
-- coordinates. Here we are still using |
|||
-- regular euclidean coordinates, however. |
|||
-- |
|||
pt0 := eval (gcurve, t0); |
|||
pt1 := ((t1 - t0 - t0) * t1 + (t0 * t0)) * gcurve.pt1; |
|||
pt2 := eval (gcurve, t1); |
|||
return (pt0 => pt0, pt1 => pt1, pt2 => pt2, |
|||
t0 => t0, t1 => t1); |
|||
end make_portion; |
|||
procedure segment_parameters (a0, a1 : opoint; |
|||
b0, b1 : opoint; |
|||
ta, tb : out lreal) is |
|||
axdiff : lreal; |
|||
aydiff : lreal; |
|||
bxdiff : lreal; |
|||
bydiff : lreal; |
|||
denom : lreal; |
|||
ta1 : lreal; |
|||
tb1 : lreal; |
|||
begin |
|||
-- Do the line segments (a0,a1) and (b0,b1) intersect? If so, |
|||
-- return t-parameter values in ta,tb in [0,1] for the point of |
|||
-- intersection, treating them as parametric splines of degree |
|||
-- 1. Otherwise return ta=-1.0 and tb=-1.0. |
|||
ta := -1.0; |
|||
tb := -1.0; |
|||
axdiff := a1.x - a0.x; |
|||
aydiff := a1.y - a0.y; |
|||
bxdiff := b1.x - b0.x; |
|||
bydiff := b1.y - b0.y; |
|||
denom := (axdiff * bydiff) - (aydiff * bxdiff); |
|||
ta1 := ((bxdiff * a0.y) - (bydiff * a0.x) |
|||
+ (b0.x * b1.y) - (b1.x * b0.y)); |
|||
ta1 := ta1 / denom; |
|||
if between_0_and_1 (ta1) then |
|||
tb1 := -((axdiff * b0.y) - (aydiff * b0.x) |
|||
+ (a0.x * a1.y) - (a1.x * a0.y)); |
|||
tb1 := tb1 / denom; |
|||
if between_0_and_1 (tb1) then |
|||
ta := ta1; |
|||
tb := tb1; |
|||
end if; |
|||
end if; |
|||
end segment_parameters; |
|||
type t_pair; |
|||
type t_pair_access is access t_pair; |
|||
type t_pair is |
|||
record |
|||
tp, tq : real; -- Single precision results. |
|||
next : t_pair_access; |
|||
end record; |
|||
procedure t_pair_delete is |
|||
new ada.unchecked_deallocation (t_pair, t_pair_access); |
|||
function length (params : t_pair_access) |
|||
return natural is |
|||
n : natural; |
|||
p : t_pair_access; |
|||
begin |
|||
n := 0; |
|||
p := params; |
|||
while p /= null loop |
|||
n := n + 1; |
|||
p := p.next; |
|||
end loop; |
|||
return n; |
|||
end length; |
|||
procedure insert_t_pair (params : in out t_pair_access; |
|||
tp, tq : real) |
|||
with pre => between_0_and_1 (tp) and between_0_and_1 (tq) |
|||
is |
|||
pair : t_pair_access; |
|||
begin |
|||
-- Test for duplicates and also sort as we insert. This is a |
|||
-- recursive implementation. |
|||
if params = null then |
|||
pair := new t_pair; |
|||
params := pair; |
|||
pair.next := null; |
|||
pair.tp := tp; |
|||
pair.tq := tq; |
|||
elsif abs (params.tp - tp) <= real'model_epsilon then |
|||
null; -- A duplicate, to within epsilon. |
|||
elsif params.tp < tp then |
|||
insert_t_pair (params.next, tp, tq); |
|||
else |
|||
pair := new t_pair; |
|||
pair.next := params.next; |
|||
params.next := pair; |
|||
pair.tp := params.tp; |
|||
pair.tq := params.tq; |
|||
params.tp := tp; |
|||
params.tq := tq; |
|||
end if; |
|||
end insert_t_pair; |
|||
procedure insert_t_pair (params : in out t_pair_access; |
|||
tp, tq : lreal) |
|||
with pre => between_0_and_1 (tp) and between_0_and_1 (tq) |
|||
is |
|||
begin |
|||
insert_t_pair (params, real (tp), real (tq)); |
|||
end insert_t_pair; |
|||
type clipping_roles is |
|||
(pportion_clips_qportion, |
|||
qportion_clips_pportion); |
|||
type intersection_job; |
|||
type intersection_job_access is access intersection_job; |
|||
type intersection_job is |
|||
record |
|||
pportion : geocurve; |
|||
qportion : geocurve; |
|||
roles : clipping_roles; |
|||
next : intersection_job_access; |
|||
end record; |
|||
procedure intersection_job_delete is |
|||
new ada.unchecked_deallocation (intersection_job, |
|||
intersection_job_access); |
|||
function find_intersections (p, q : geocurve) |
|||
return t_pair_access is |
|||
params : t_pair_access; |
|||
workload : intersection_job_access; |
|||
pportion : geocurve; |
|||
qportion : geocurve; |
|||
roles : clipping_roles; |
|||
tp, tq : lreal; |
|||
tolerance : constant lreal := 0.4 * lreal (real'model_epsilon); |
|||
bisection_threshold : constant lreal := 0.2; |
|||
procedure insert_job (pportion : geocurve; |
|||
qportion : geocurve; |
|||
roles : clipping_roles) |
|||
with post => workload /= null |
|||
is |
|||
job : intersection_job_access; |
|||
begin |
|||
job := new intersection_job; |
|||
job.pportion := pportion; |
|||
job.qportion := qportion; |
|||
job.roles := roles; |
|||
job.next := workload; |
|||
workload := job; |
|||
end insert_job; |
|||
procedure insert_job (pportion : geocurve; |
|||
qportion : geocurve) |
|||
with post => workload /= null |
|||
is |
|||
roles : clipping_roles; |
|||
begin |
|||
if pportion.t1 - pportion.t0 <= qportion.t1 - qportion.t0 then |
|||
roles := pportion_clips_qportion; |
|||
else |
|||
roles := qportion_clips_pportion; |
|||
end if; |
|||
insert_job (pportion, qportion, roles); |
|||
end insert_job; |
|||
procedure choose_job |
|||
with pre => workload /= null |
|||
is |
|||
job : intersection_job_access; |
|||
begin |
|||
job := workload; |
|||
workload := workload.next; |
|||
pportion := job.pportion; |
|||
qportion := job.qportion; |
|||
roles := job.roles; |
|||
intersection_job_delete (job); |
|||
end choose_job; |
|||
function t_parameters_are_within_tolerance |
|||
return boolean is |
|||
begin |
|||
return (pportion.t1 - pportion.t0 <= tolerance |
|||
and qportion.t1 - qportion.t0 <= tolerance); |
|||
end t_parameters_are_within_tolerance; |
|||
procedure insert_intersection_parameters is |
|||
begin |
|||
segment_parameters (pportion.pt0, pportion.pt2, |
|||
qportion.pt0, qportion.pt2, tp, tq); |
|||
if between_0_and_1 (tp) and between_0_and_1 (tq) then |
|||
tp := (1.0 - tp) * pportion.t0 + tp * pportion.t1; |
|||
tq := (1.0 - tq) * qportion.t0 + tq * qportion.t1; |
|||
insert_t_pair (params, tp, tq); |
|||
end if; |
|||
end insert_intersection_parameters; |
|||
procedure no_intersections is |
|||
begin |
|||
-- Do nothing. The reason to have a procedure for this is so the |
|||
-- code will document itself. |
|||
null; |
|||
end no_intersections; |
|||
procedure bisect_them is |
|||
pportion1 : geocurve; |
|||
pportion2 : geocurve; |
|||
qportion1 : geocurve; |
|||
qportion2 : geocurve; |
|||
t0, t1, th : lreal; |
|||
begin |
|||
-- |
|||
-- Sederberg employs heuristics to decide what to bisect and who |
|||
-- should be the clipper. We will simply bisect both curves, and |
|||
-- use some heuristic coded elsewhere to choose the |
|||
-- clipper. |
|||
-- |
|||
-- (Bisection as it is normally done in the s-power basis is |
|||
-- equivalent to change of homogeneous coordinates. It can be |
|||
-- applied to p and q directly and, for quadratics, is a simple |
|||
-- operation.) |
|||
-- |
|||
t0 := pportion.t0; |
|||
t1 := pportion.t1; |
|||
th := 0.5 * t0 + 0.5 * t1; |
|||
pportion1 := make_portion (p, t0, th); |
|||
pportion2 := make_portion (p, th, t1); |
|||
t0 := qportion.t0; |
|||
t1 := qportion.t1; |
|||
th := 0.5 * t0 + 0.5 * t1; |
|||
qportion1 := make_portion (q, t0, th); |
|||
qportion2 := make_portion (q, th, t1); |
|||
insert_job (pportion1, qportion1); |
|||
insert_job (pportion1, qportion2); |
|||
insert_job (pportion2, qportion1); |
|||
insert_job (pportion2, qportion2); |
|||
end bisect_them; |
|||
procedure clip_and_maybe_bisect (plan : clipping_plan) is |
|||
clipping : geocurve; |
|||
t_lft : lreal; |
|||
t_rgt : lreal; |
|||
t0, t1 : lreal; |
|||
begin |
|||
t_lft := plan.t_lft; |
|||
t_rgt := plan.t_rgt; |
|||
if roles = pportion_clips_qportion then |
|||
t0 := (1.0 - t_lft) * qportion.t0 + t_lft * qportion.t1; |
|||
t1 := (1.0 - t_rgt) * qportion.t0 + t_rgt * qportion.t1; |
|||
clipping := make_portion (q, t0, t1); |
|||
if 1.0 - (t_rgt - t_lft) < bisection_threshold then |
|||
qportion := clipping; |
|||
bisect_them; |
|||
else |
|||
insert_job (pportion, clipping, qportion_clips_pportion); |
|||
end if; |
|||
else |
|||
t0 := (1.0 - t_lft) * pportion.t0 + t_lft * pportion.t1; |
|||
t1 := (1.0 - t_rgt) * pportion.t0 + t_rgt * pportion.t1; |
|||
clipping := make_portion (p, t0, t1); |
|||
if 1.0 - (t_rgt - t_lft) < bisection_threshold then |
|||
pportion := clipping; |
|||
bisect_them; |
|||
else |
|||
insert_job (clipping, qportion, pportion_clips_qportion); |
|||
end if; |
|||
end if; |
|||
end clip_and_maybe_bisect; |
|||
procedure perform_clipping is |
|||
clippehend, clipper : geocurve; |
|||
projected_clippehend : geocurve; |
|||
projected_clipper : geocurve; |
|||
trans : affine_projection; |
|||
dmin, dmax : lreal; |
|||
plan : clipping_plan; |
|||
begin |
|||
if roles = pportion_clips_qportion then |
|||
clipper := pportion; |
|||
clippehend := qportion; |
|||
else |
|||
clipper := qportion; |
|||
clippehend := pportion; |
|||
end if; |
|||
trans := skeleton_projection (clipper); |
|||
projected_clipper := project (clipper, trans); |
|||
fat_line_bounds (projected_clipper, dmin, dmax); |
|||
projected_clippehend := project (clippehend, trans); |
|||
plan := make_clipping_plan (projected_clippehend, dmin, dmax); |
|||
case plan.where is |
|||
when clip_everywhere => no_intersections; |
|||
when clip_nowhere => bisect_them; |
|||
when clip_left_and_right => clip_and_maybe_bisect (plan); |
|||
end case; |
|||
end perform_clipping; |
|||
begin |
|||
params := null; |
|||
workload := null; |
|||
insert_job (p, q); |
|||
while workload /= null loop |
|||
choose_job; |
|||
if t_parameters_are_within_tolerance then |
|||
insert_intersection_parameters; |
|||
else |
|||
perform_clipping; |
|||
end if; |
|||
end loop; |
|||
return params; |
|||
end find_intersections; |
|||
function find_intersections (p, q : spower_curve) |
|||
return t_pair_access is |
|||
begin |
|||
return find_intersections (to_geocurve (p), to_geocurve (q)); |
|||
end find_intersections; |
|||
pctrl, qctrl : control_point_curve; |
|||
p, q : spower_curve; |
|||
params, ptr : t_pair_access; |
|||
begin |
|||
test_segment_projection_to_0_1; |
|||
pctrl := ((x => -1.0, y => 0.0), |
|||
(x => 0.0, y => 10.0), |
|||
(x => 1.0, y => 0.0)); |
|||
qctrl := ((x => 2.0, y => 1.0), |
|||
(x => -8.0, y => 2.0), |
|||
(x => 2.0, y => 3.0)); |
|||
p := to_spower_curve (pctrl); |
|||
q := to_spower_curve (qctrl); |
|||
params := find_intersections (p, q); |
|||
set_col (to => 9); |
|||
put ("convex up"); |
|||
set_col (to => 44); |
|||
put_line ("convex left"); |
|||
ptr := params; |
|||
while ptr /= null loop |
|||
set_col (to => 3); |
|||
put (ptr.tp, fore => 1, aft => 6, exp => 0); |
|||
set_col (to => 13); |
|||
put ("("); |
|||
put (eval (p, ptr.tp).x, fore => 1, aft => 6, exp => 0); |
|||
put (", "); |
|||
put (eval (p, ptr.tp).y, fore => 1, aft => 6, exp => 0); |
|||
put (")"); |
|||
set_col (to => 38); |
|||
put (ptr.tq, fore => 1, aft => 6, exp => 0); |
|||
set_col (to => 48); |
|||
put ("("); |
|||
put (eval (q, ptr.tq).x, fore => 1, aft => 6, exp => 0); |
|||
put (", "); |
|||
put (eval (q, ptr.tq).y, fore => 1, aft => 6, exp => 0); |
|||
put_line (")"); |
|||
ptr := ptr.next; |
|||
end loop; |
|||
end bezier_intersections; |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> convex up convex left |
|||
0.072508 (-0.854983, 1.345017) 0.172508 (-0.854983, 1.345017) |
|||
0.159488 (-0.681025, 2.681026) 0.840512 (-0.681026, 2.681025) |
|||
0.827492 (0.654984, 2.854983) 0.927492 (0.654983, 2.854983) |
|||
0.940512 (0.881025, 1.118975) 0.059488 (0.881025, 1.118975) |
|||
</pre> |
|||
=={{header|ATS}}== |
=={{header|ATS}}== |