## beta2-spline.pkg
# Compiled by:
# src/lib/x-kit/draw/xkit-draw.sublib### "When the blind beetle crawls over the surface of a globe,
### he doesn't realize that the track he has covered is curved.
### I was lucky enough to have spotted it."
###
### -- Albert Einstein
stipulate
package g2d = geometry2d; # geometry2d is from src/lib/std/2d/geometry2d.pkgherein
package beta2_spline
: (weak) Beta2_Spline # Beta2_Spline is from src/lib/x-kit/draw/beta2-spline.api {
include package geometry2d; # geometry2d is from src/lib/std/2d/geometry2d.pkg fun round x
=
if (x > 0.0) floor (x+0.5);
else -1*floor(-x+0.5);
fi;
fun add_seg ([], x0, y0, x1, y1)
=>
[ { col => round x0, row => round y0 },
{ col => round x1, row => round y1 }
];
add_seg (l, x0, y0, _, _)
=>
{ col => round x0, row => round y0 } ! l;
end;
# is_flat:
# Returns TRUE if the polygon determined by the four points
# is flat enough. Flatness is measured by the maximum distance
# of (x1, y1) and (x2, y2) from the line determined by (x0, y0)
# and (x3, y3). In addition, check that p1, p2 are close to the
# line segment. To do this, make sure they are roughly within
# the circle with center (p0+p3)/2 and radius = |p3-p0|/2+flatness
#
fun is_flat { x0, y0, x1, y1, x2, y2, x3, y3 }
=
{ fun sqr x = x * x;
#
dx = x3 - x0;
dy = y3 - y0;
midx = 0.5*dx;
midy = 0.5*dy;
dist2 = sqr dy + sqr dx;
flatness2 = sqr 1.0 * dist2;
halfd2 = 0.25*dist2;
fun in_flat_range (x, y)
=
sqr (dy * x - dx * y) <= flatness2
and
{ d = sqr (x-midx) + sqr (y-midy);
d <= halfd2
or
sqr (d-halfd2) <= flatness2;
};
in_flat_range (x1-x0, y1-y0)
and
in_flat_range (x2-x0, y2-y0);
};
# bezier:
# Recursively compute a Bezier cubic section. If the points
# determine a polygon flat enough to be represented as a line
# segment, the segment is added to the list. Otherwise, the
# the curve is bisected and each part is recursively computed,
# with the lists concatenated.
#
# From "The Beta2-split: A special case of the Beta-spline Curve and
# Surface Representation." B. A. Barsky and A. D. DeRose. IEEE, 1985,
# http://www.eecs.berkeley.edu/Pubs/TechRpts/1983/CSD-83-152.pdf
# as adapted by Crispin Goswell for xps.
#
fun bezier (arg as { x0, y0, x1, y1, x2, y2, x3, y3 }, l)
=
if (is_flat arg)
#
add_seg (l, x0, y0, x3, y3);
else
mid_x = (x0 + x3) / 8.0 + 3.0 * (x1 + x2) / 8.0;
mid_y = (y0 + y3) / 8.0 + 3.0 * (y1 + y2) / 8.0;
l' = bezier (
{ x0 => mid_x, y0 => mid_y,
x1 => (x1+x3) / 4.0 + x2 / 2.0, y1 => (y1+y3) / 4.0 + y2 / 2.0,
x2 => (x2+x3) / 2.0, y2 => (y2+y3) / 2.0,
x3, y3
},
l
);
bezier (
{ x0, y0,
x1 => (x0+x1) / 2.0, y1 => (y0+y1) / 2.0,
x2 => (x0+x2) / 4.0 + x1 / 2.0, y2 => (y0+y2) / 4.0 + y1 / 2.0,
x3 => mid_x, y3 => mid_y
},
l'
);
fi;
# curve:
# Given four points [p0, p1, p2, p3], return a list of points corresponding to
# to a Bezier cubic section, starting at p0, ending at p3, with p1, p2 as
# control points.
#
fun curve
( { col=>x0, row=>y0 }, { col=>x1, row=>y1 },
{ col=>x2, row=>y2 }, { col=>x3, row=>y3 }
)
=
bezier (
{ x0 => float(x0), y0 => float(y0),
x1 => float(x1), y1 => float(y1),
x2 => float(x2), y2 => float(y2),
x3 => float(x3), y3 => float(y3)
},
[]
);
# doSpline:
# Given four points (p0, p1, p2, p3), return a list of points corresponding to
# to the B-spline curve section, accumulating the results on the argument list.
# We compute the curve by determining the corresponding Bezier control points,
# and then use the Bezier routines above.
#
fun do_spline
( { col=>x0, row=>y0 },
{ col=>x1, row=>y1 },
{ col=>x2, row=>y2 },
{ col=>x3, row=>y3 },
l
)
=
{ x0 = float(x0); y0 = float(y0);
x1 = float(x1); y1 = float(y1);
x2 = float(x2); y2 = float(y2);
x3 = float(x3); y3 = float(y3);
bezier (
{ x0 => (x0 + 4.0*x1 + x2)/6.0, y0 => (y0 + 4.0*y1 + y2)/6.0,
x1 => (2.0*x1 + x2)/3.0, y1 => (2.0*y1 + y2)/3.0,
x2 => (x1 + 2.0*x2)/3.0, y2 => (y1 + 2.0*y2)/3.0,
x3 => (x1 + 4.0*x2 + x3)/6.0, y3 => (y1 + 4.0*y2 + y3)/6.0
},
l
);
};
# loopSpline:
# Given a list of spline control points, generate the corresponding
# spline. Since we accumulate on the head of the list, we assume
# the calling function has reversed the list of control points.
# The loop continues as long as there are 4 control points left.
#
fun loop_spline (p3, p2, p1, p0 ! tl, l)
=>
loop_spline (p2, p1, p0, tl, do_spline (p0, p1, p2, p3, l));
loop_spline (_, _, _, [], l)
=>
l;
end;
# simpleBSpline:
# Compute a simple B-spline with the given control points.
#
fun simple_bspline arg
=
case (reverse arg)
#
(p3 ! p2 ! p1 ! tl) => loop_spline (p3, p2, p1, tl,[]);
_ => arg;
esac;
# bSpline:
# Compute a B-spline using the given control points.
# In addition, we constrain the resultant spline to connect the
# first and last points by adding copies of these points.
#
fun b_spline (arg as (p0 ! _ ! _ ! _))
=>
loop_spline (pn, pn, pn, tl, [])
where
my (pn, tl)
=
case (reverse (p0 ! p0 ! arg)) (pn ! tl) => (pn, tl);
/* */ _ => raise exception DIE "Bug: Unsupported case in b_spline";
esac;
end;
b_spline l
=>
l;
end;
# closed_bSpline:
# Compute a closed B-spline. This is done by repeating the first
# three points at the end of the list.
# Note that the first and last points of the result are the same.
#
fun closed_bspline (arg as (p0 ! p1 ! p2 ! _))
=>
loop_spline (p2, p1, p0, reverse arg,[]);
closed_bspline l
=>
l;
end;
}; # package spline
end;