Bresenham tasks in ATS
The following tasks are implemented here in ATS:
- Bitmap/Bresenham's_line_algorithm
- Bitmap/Midpoint_circle_algorithm
- Bitmap/Bézier_curves/Quadratic
- Bitmap/Bézier_curves/Cubic
You will also need the source files from Bitmap#ATS and Bitmap/Write_a_PPM_file#ATS.
The ATS static file
The following file should be named bitmap_bresenham_tasks.sats.
#define ATS_PACKNAME "Rosetta_Code.bresenham_tasks"
%{#
#include <math.h>
%}
staload "bitmap_task.sats"
fn {a : t@ype} {tk : tkind}
pixmap_draw_line_endpoints_bresenham
(pix : !pixmap a,
endpt0 : @(g0int tk, g0int tk),
endpt1 : @(g0int tk, g0int tk),
color : a) : void
fn {a : t@ype} {tk : tkind}
pixmap_draw_circle_midpoint_algorithm
(pix : !pixmap a,
midpt : @(g0int tk, g0int tk),
radius : g0int tk,
color : a) : void
fn {a : t@ype} {tk : tkind}
pixmap_draw_bezier_degree2_control_points
(pix : !pixmap a,
ctrlpt0 : @(g0float tk, g0float tk),
ctrlpt1 : @(g0float tk, g0float tk),
ctrlpt2 : @(g0float tk, g0float tk),
color : a) : void
fn {a : t@ype} {tk : tkind}
pixmap_draw_bezier_degree3_control_points
(pix : !pixmap a,
ctrlpt0 : @(g0float tk, g0float tk),
ctrlpt1 : @(g0float tk, g0float tk),
ctrlpt2 : @(g0float tk, g0float tk),
ctrlpt3 : @(g0float tk, g0float tk),
color : a) : void
overload pixmap_draw_line with pixmap_draw_line_endpoints_bresenham
overload pixmap_draw_circle with pixmap_draw_circle_midpoint_algorithm
overload pixmap_draw_bezier_degree2 with
pixmap_draw_bezier_degree2_control_points
overload pixmap_draw_bezier_degree3 with
pixmap_draw_bezier_degree3_control_points
overload pixmap_draw_bezier with pixmap_draw_bezier_degree2
overload pixmap_draw_bezier with pixmap_draw_bezier_degree3
overload draw_line with pixmap_draw_line
overload draw_circle with pixmap_draw_circle
overload draw_bezier with pixmap_draw_bezierThe ATS dynamic file
The following file should be named bitmap_bresenham_tasks.dats.
(*------------------------------------------------------------------*)
#define ATS_DYNLOADFLAG 0
#define ATS_PACKNAME "Rosetta_Code.bresenham_tasks"
#include "share/atspre_staload.hats"
staload "bitmap_task.sats"
staload _ = "bitmap_task.dats"
staload "bitmap_bresenham_tasks.sats"
(*------------------------------------------------------------------*)
(* The ats2-xprelude package implements these, but for this task we
will implement them ourselves: *)
extern castfn llint2size : {i : nat} llint i -<> size_t i
implement g1int2uint<llintknd,sizeknd> i = llint2size i
extern castfn size2llint : {i : int} size_t i -<> llint i
implement g1uint2int<sizeknd,llintknd> i =
let prval () = lemma_g1uint_param i in size2llint i end
(*------------------------------------------------------------------*)
(* The ats2-xprelude package has "sqrt" that selects the correct
function for a given floating point type. However, for this Rosetta
Code task, I am avoiding the dependency. We will define it for
double, only. *)
extern fn {tk : tkind} g0float_sqrt : g0float tk -<> g0float tk
implement g0float_sqrt<dblknd> x = $extfcall (double, "sqrt", x)
overload sqrt with g0float_sqrt
(*------------------------------------------------------------------*)
(* "Bisection iterators" -- a way to avoid using extra space when
doing any of many adaptive algorithms.
In ats2-xprelude, this facility is provided for exact rational
numbers, going either forwards or backwards. Here we will use
ullint and go only forwards. (WARNING: Using fixed size integers
causes some risk of overflow. This risk is similar to that, in a
recursive implementation, of using up all your stack.) *)
typedef bisect_iter = @{n = ullint, d = ullint}
fn {}
bisect_iter_make () :<> bisect_iter =
@{n = 0LLU, d = 1LLU}
fn {}
bisect_iter_done (iter : bisect_iter) :<> bool =
(iter.n = iter.d && iter.d = 1LLU)
fn {tk : tkind}
bisect_iter_interval (iter : bisect_iter)
:<> @(g0float tk, g0float tk) =
(* Get the endpoints of the current interval, as floating point
numbers between zero and one (inclusive). *)
let
(* The following works only for floating point types built into
C. For other types (such as mpfr), you may get a WRONG
result. *)
extern castfn ulli2f : ullint -<> g0float tk
in
@(ulli2f (iter.n) / ulli2f (iter.d),
ulli2f (succ (iter.n)) / ulli2f (iter.d))
end
fn {}
bisect_iter_bisect (iter : bisect_iter) : bisect_iter =
(* Bisect the current interval. *)
@{n = (iter.n << 1), d = (iter.d << 1)}
fn {}
bisect_iter_next (iter : bisect_iter) : bisect_iter =
(* Advance to the next interval. *)
let
extern castfn ull2ll : ullint -<> llint
val n = succ (iter.n)
(* WARNING: __builtin_ffsll is a GNU C extension, available also
with some other compilers people use with ATS. If you
need to, you can write a replacement for your system
(or take one from Gnulib). See
https://en.wikipedia.org/w/index.php?title=Find_first_set&oldid=1139170019
and
https://www.chessprogramming.org/BitScan#Bitscan_reverse
*)
val shift_amount : intGte 0 =
pred ($extfcall (intGte 1, "__builtin_ffsll", ull2ll n))
in
@{n = (n >> shift_amount), d = ((iter.d) >> shift_amount)}
end
(*------------------------------------------------------------------*)
(* Bresenham's line-drawing algorithm. See
https://en.wikipedia.org/w/index.php?title=Bresenham%27s_line_algorithm&oldid=1145813731
(This implementation may produce a fair bit of C code. I was not
parsimonious.) *)
extern fn {a : t@ype}
pixmap_draw_line$plot :
{w, h : int}
{x, y : nat | x < w; y < h}
(!pixmap (a, w, h), size_t x, size_t y, a) -> void
fn {a : t@ype} {tk : tkind}
pixmap_draw_line_shallow
{w, h : int}
{x0, y0, x1, y1 : int | x0 <= x1; y0 <= y1;
y1 - y0 <= x1 - x0}
(pix : !pixmap (a, w, h),
x0 : g1int (tk, x0),
y0 : g1int (tk, y0),
x1 : g1int (tk, x1),
y1 : g1int (tk, y1),
color : a)
: void =
let
val w : size_t w = width pix
and h : size_t h = height pix
val dx = x1 - x0
and dy = y1 - y0
val twice_dx = dx + dx
and twice_dy = dy + dy
val twice_dy_minus_dx = twice_dy - twice_dx
fun
loop {x, y : int | x <= x1; y <= y1}
{p : addr}
.<x1 - x>.
(pix : !pixmap (a, w, h, p),
x : g1int (tk, x),
y : g1int (tk, y),
D : g1int tk) : void =
begin
if (g1i2i 0 <= x) * (x < g1u2i w)
* (g1i2i 0 <= y) * (y < g1u2i h) then
pixmap_draw_line$plot<a> (pix, g1i2u x, g1i2u y, color);
if x <> x1 then
begin
if g1i2i 0 < D then
(* The call to "min" here is to provide a proof that the
value is no greater than y1. If you can come up with
an actual proof, then go ahead and get rid of the
"min". *)
loop (pix, succ x, min (succ y, y1),
D + twice_dy_minus_dx)
else
loop (pix, succ x, y, D + twice_dy)
end
end
in
loop (pix, x0, y0, twice_dy - dx)
end
fn {a : t@ype} {tk : tkind}
pixmap_draw_line_steep
{w, h : int}
{x0, y0, x1, y1 : int | x0 <= x1; y0 <= y1;
x1 - x0 <= y1 - y0}
(pix : !pixmap (a, w, h),
x0 : g1int (tk, x0),
y0 : g1int (tk, y0),
x1 : g1int (tk, x1),
y1 : g1int (tk, y1),
color : a)
: void =
let
val w : size_t w = width pix
and h : size_t h = height pix
val dx = x1 - x0
and dy = y1 - y0
val twice_dx = dx + dx
and twice_dy = dy + dy
val twice_dx_minus_dy = twice_dx - twice_dy
fun
loop {x, y : int | x <= x1; y <= y1}
{p : addr}
.<y1 - y>.
(pix : !pixmap (a, w, h, p),
x : g1int (tk, x),
y : g1int (tk, y),
D : g1int tk) : void =
begin
if (g1i2i 0 <= x) * (x < g1u2i w)
* (g1i2i 0 <= y) * (y < g1u2i h) then
pixmap_draw_line$plot<a> (pix, g1i2u x, g1i2u y, color);
if y <> y1 then
begin
if g1i2i 0 < D then
(* The call to "min" here is to provide a proof that the
value is no greater than x1. If you can come up with
an actual proof, then go ahead and get rid of the
"min". *)
loop (pix, min (succ x, x1), succ y,
D + twice_dx_minus_dy)
else
loop (pix, x, succ y, D + twice_dx)
end
end
in
loop (pix, x0, y0, twice_dx - dy)
end
implement {a} {tk}
pixmap_draw_line_endpoints_bresenham (pix, @(x0, y0), @(x1, y1),
color) =
let
macdef draw_shallow = pixmap_draw_line_shallow<a>
macdef draw_steep = pixmap_draw_line_steep<a>
val x0 = g1ofg0 x0 and y0 = g1ofg0 y0
and x1 = g1ofg0 x1 and y1 = g1ofg0 y1
in
if x0 <= x1 then
begin
if y0 <= y1 then
let
implement
pixmap_draw_line$plot<a> (pix, x, y, color) =
pix[x, y] := color
in
if y1 - y0 <= x1 - x0 then
draw_shallow (pix, x0, y0, x1, y1, color)
else
draw_steep (pix, x0, y0, x1, y1, color)
end
else
let
implement
pixmap_draw_line$plot<a> (pix, x, y, color) =
let
val h1 = pred (height pix)
in
pix[x, h1 - y] := color
end
val h1 : g1int tk = pred (g1u2i (height pix))
in
if y0 - y1 <= x1 - x0 then
draw_shallow (pix, x0, h1 - y0, x1, h1 - y1, color)
else
draw_steep (pix, x0, h1 - y0, x1, h1 - y1, color)
end
end
else
begin
if y0 <= y1 then
let
implement
pixmap_draw_line$plot<a> (pix, x, y, color) =
let
val w1 = pred (width pix)
in
pix[w1 - x, y] := color
end
val w1 : g1int tk = pred (g1u2i (width pix))
in
if y1 - y0 <= x0 - x1 then
draw_shallow (pix, w1 - x0, y0, w1 - x1, y1, color)
else
draw_steep (pix, w1 - x0, y0, w1 - x1, y1, color)
end
else
let
implement
pixmap_draw_line$plot<a> (pix, x, y, color) =
let
val w1 = pred (width pix)
and h1 = pred (height pix)
in
pix[w1 - x, h1 - y] := color
end
val w1 : g1int tk = pred (g1u2i (width pix))
and h1 : g1int tk = pred (g1u2i (height pix))
in
if y0 - y1 <= x0 - x1 then
draw_shallow (pix, w1 - x0, h1 - y0,
w1 - x1, h1 - y1, color)
else
draw_steep (pix, w1 - x0, h1 - y0,
w1 - x1, h1 - y1, color)
end
end
end
(*------------------------------------------------------------------*)
(* The midpoint algorithm for drawing a circle.
There are several closely related algorithms. The one I am using is
called "Midpoint Method" at https://archive.ph/g0Z9e (although that
presentation contains errors), and seems to be what most others
have been using.
The Wikipedia article is, at the time of this writing, featuring
"Jesko's Method". Presumably the page used to feature the "Midpoint
Method" proper.
I might mention that, when I use Wikipedia as a reference, I use
the constant link for a page. Pages not only change--they can
completely disappear. But constant links remain. Alternatively one
can save a page at archive.today or the Wayback Machine. *)
implement {a} {tk}
pixmap_draw_circle_midpoint_algorithm (pix, @(x0, y0), radius,
color) =
let
macdef zero = g1i2i 0
macdef one = g1i2i 1
val radius = g1ofg0 radius
val radius = abs radius
prval [radius : int] EQINT () = eqint_make_gint radius
prval () = prop_verify {0 <= radius} ()
val [x0 : int] x0 = g1ofg0 x0
and [y0 : int] y0 = g1ofg0 y0
fn {}
plot {x, y : int}
{w, h : int}
{p : addr}
(pix : !pixmap (a, w, h, p),
x : g1int (tk, x),
y : g1int (tk, y)) : void =
if (zero <= x) * (x < g1u2i (width pix))
* (zero <= y) * (y < g1u2i (height pix)) then
pix[x, y] := color
fun
loop {w, h : int}
{p : addr}
{x, y : int | y <= x + 2}
.<(x + 2) - y>.
(pix : !pixmap (a, w, h, p),
x : g1int (tk, x),
y : g1int (tk, y),
P : g1int tk) : void =
if y <= x then
begin
plot (pix, x0 - x, y0 - y);
plot (pix, x0 - x, y0 + y);
plot (pix, x0 + x, y0 - y);
plot (pix, x0 + x, y0 + y);
plot (pix, x0 - y, y0 - x);
plot (pix, x0 - y, y0 + x);
plot (pix, x0 + y, y0 - x);
plot (pix, x0 + y, y0 + x);
let
val y = succ y
in
if P <= zero then
loop (pix, x, y, succ (P + y + y))
else
let
val x = pred x
val diff = y - x
in
loop (pix, x, y, succ (P + diff + diff))
end
end
end
in
loop (pix, radius, zero, one - radius)
end
(*------------------------------------------------------------------*)
(* Symmetric power (Sánchez-Reyes) representation of Bézier
splines. S-power representation has the advantage that
coefficients, other than those representing the endpoints, tend
towards zero as the curve gets flatter.
References:
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.
*)
(* - - - - - - - - - - - - - - - - - - - - - - *)
(* Degree 2. *)
fn {tk : tkind}
bernstein2spower_degree2
@(c0 : g0float tk, c1 : g0float tk, c2 : g0float tk)
:<> @(g0float tk, g0float tk, g0float tk) =
(* Convert from Bernstein coefficients ("control points") to s-power
coefficients. *)
let
(* First convert to scaled Bernstein basis. *)
val c1 = c1 + c1
in
(* Then to s-power basis. *)
@(c0, c1 - c0 - c2, c2)
end
fn {tk : tkind}
spower_eval_degree2
(@(c0 : g0float tk, c1 : g0float tk, c2 : g0float tk),
t : g0float tk)
:<> g0float tk =
(* Evaluate the polynomial spower(c0, c1, c2) at t. (This function
is used in the test code.) *)
let
val one_minus_t = g0i2f 1 - t
val s = t * one_minus_t
val c1_s = c1 * s
val y0 = c0 + c1_s and y1 = c2 + c1_s
in
(y0 * one_minus_t) + (y1 * t)
end
fn {tk : tkind}
spower_portion_degree2
(@(c0 : g0float tk, c1 : g0float tk, c2 : g0float tk),
@(t0 : g0float tk, t1 : g0float tk))
: @(g0float tk, g0float tk, g0float tk) =
(* Compose spower(c0, c1, c2) with spower(t0, t1). This will map the
portion [t0,t1] onto [0,1]. (I used Maxima to compute these
coefficients symbolically.) *)
let
val t0_t0 = t0 * t0
and t0_t1 = t0 * t1
and t1_t1 = t1 * t1
and c2p1m0 = c2 + c1 - c0
val d0 = c0 + (c2p1m0 * t0) - (c1 * t0_t0)
and d1 = (c1 * t1_t1) - ((c1 + c1) * t0_t1) + (c1 * t0_t0)
and d2 = c0 + (c2p1m0 * t1) - (c1 * t1_t1)
in
@(d0, d1, d2)
end
fn {a : t@ype} {tk : tkind}
spower_draw_degree2
{w, h : int}
(pix : !pixmap (a, w, h),
cx : @(g0float tk, g0float tk, g0float tk),
cy : @(g0float tk, g0float tk, g0float tk),
color : a)
: void =
(* The algorithm is adaptive flattening by repeated bisection. This
is done in the s-power basis. (The algorithm of de Casteljau is
NOT used.) *)
let
extern castfn sz2lli : {i : int} size_t i -<> llint i
extern castfn dbl2lli : double -<> [i : int] llint i
extern castfn f2dbl : g0float tk -<> double
fn {}
round_to_lli (x : g0float tk) :<> [i : int] llint i =
let
val x = f2dbl x
val x = $extfcall (double, "rintl", x)
in
dbl2lli x
end
#define TOL 0.25 (* QUESTION: What is a good value for TOL? *)
fun
flatten_and_draw
{p : addr}
(pix : !pixmap (a, w, h, p),
iter : bisect_iter)
: void =
if ~bisect_iter_done iter then
let
val interv = bisect_iter_interval<tk> iter
val dx = spower_portion_degree2<tk> (cx, interv)
and dy = spower_portion_degree2<tk> (cy, interv)
in
if abs dx.1 < g0f2f TOL
&& abs dy.1 < g0f2f TOL then
let
val w = sz2lli (width pix) and h = sz2lli (height pix)
val dx0 = round_to_lli dx.0
and dy0 = round_to_lli dy.0
and dx2 = round_to_lli dx.2
and dy2 = round_to_lli dy.2
in
if (0LL <= dx0) * (dx0 < w)
* (0LL <= dy0) * (dy0 < h)
* (0LL <= dx2) * (dx2 < w)
* (0LL <= dy2) * (dy2 < h) then
draw_line<a><llintknd>
(pix, @(dx0, dy0), @(dx2, dy2), color);
flatten_and_draw (pix, bisect_iter_next iter)
end
else
flatten_and_draw (pix, bisect_iter_bisect iter)
end
in
flatten_and_draw (pix, bisect_iter_make ())
end
implement {a} {tk}
pixmap_draw_bezier_degree2_control_points (pix, pt0, pt1, pt2,
color) =
let
(* Bernstein coefficients. *)
val bx = @(pt0.0, pt1.0, pt2.0)
and by = @(pt0.1, pt1.1, pt2.1)
(* s-power coefficients. *)
val cx = bernstein2spower_degree2<tk> bx
and cy = bernstein2spower_degree2<tk> by
in
spower_draw_degree2<a><tk> (pix, cx, cy, color)
end
(* - - - - - - - - - - - - - - - - - - - - - - *)
(* Degree 3. *)
fn {tk : tkind}
bernstein2spower_degree3
@(c0 : g0float tk, c1 : g0float tk,
c2 : g0float tk, c3 : g0float tk)
:<> @(g0float tk, g0float tk, g0float tk, g0float tk) =
(* Convert from Bernstein coefficients ("control points") to s-power
coefficients. *)
let
(* First convert to scaled Bernstein basis. *)
val c1 = g0i2f 3 * c1
and c2 = g0i2f 3 * c2
in
(* Then to s-power basis. *)
@(c0, c1 - c0 - c0 - c3, c2 - c0 - c3 - c3, c3)
end
fn {tk : tkind}
spower_eval_degree3
(@(c0 : g0float tk, c1 : g0float tk,
c2 : g0float tk, c3 : g0float tk),
t : g0float tk)
:<> g0float tk =
(* Evaluate the polynomial spower(c0, c1, c2, c3) at t. (This
function is used in the test code.) *)
let
val one_minus_t = g0i2f 1 - t
val s = t * one_minus_t
val y0 = c0 + (c1 * s) and y1 = c3 + (c2 * s)
in
(y0 * one_minus_t) + (y1 * t)
end
fn {tk : tkind}
spower_portion_degree3
(@(c0 : g0float tk, c1 : g0float tk,
c2 : g0float tk, c3 : g0float tk),
@(t0 : g0float tk, t1 : g0float tk))
: @(g0float tk, g0float tk, g0float tk, g0float tk) =
(* Compose spower(c0, c1, c2, c3) with spower(t0, t1). This will map
the portion [t0,t1] onto [0,1]. (I used Maxima to compute these
coefficients symbolically.) *)
let
(* I do not try to find the most efficient Horner scheme. Let us
simply expand out the products of t0 and t1 that will be
needed, and try to group them a bit. *)
val t0_t0 = t0 * t0
and t0_t1 = t0 * t1
and t1_t1 = t1 * t1
val t0_t0_t0 = t0_t0 * t0
and t0_t0_t1 = t0_t0 * t1
and t0_t1_t1 = t0 * t1_t1
and t1_t1_t1 = t1 * t1_t1
and c3p1m0 = c3 + c1 - c0
and c2m1m1 = c2 - c1 - c1
and c1m2 = c1 - c2
val c2m1 = ~c1m2
and three_c1m2 = g0i2f 3 * c1m2
val two_c2m1 = c2m1 + c2m1
and two_c2m1m1 = c2m1m1 + c2m1m1
val d0 = c0 + (c3p1m0 * t0) + (c2m1m1 * t0_t0) + (c1m2 * t0_t0_t0)
and d1 = (c2m1 * t1_t1_t1) + (three_c1m2 * t0_t0_t1)
+ (two_c2m1m1 * t0_t1) + (two_c2m1 * t0_t0_t0)
- (c2m1m1 * (t0_t0 + t1_t1))
and d2 = (c2m1 * t0_t0_t0) + (three_c1m2 * t0_t1_t1)
+ (two_c2m1m1 * t0_t1) + (two_c2m1 * t1_t1_t1)
- (c2m1m1 * (t0_t0 + t1_t1))
and d3 = c0 + (c3p1m0 * t1) + (c2m1m1 * t1_t1) + (c1m2 * t1_t1_t1)
in
@(d0, d1, d2, d3)
end
fn {tk : tkind}
spower_inflections_degree3
(@(cx0 : g0float tk, cx1 : g0float tk,
cx2 : g0float tk, cx3 : g0float tk),
@(cy0 : g0float tk, cy1 : g0float tk,
cy2 : g0float tk, cy3 : g0float tk))
: @(g0float tk, g0float tk) =
(* Find any value of t at which there is an inflection or cusp. Note
that such values may lie outside [0,1]. *)
let
val nan : g0float tk = g0f2f ($extval (g0float tk, "NAN"))
(* Translate the curve so its first endpoint be at (0,0). *)
val cx3m0 = cx3 - cx0
and cy3m0 = cy3 - cy0
(* I used Maxima to compute this solution. It is for a curve
translated so its first endpoint should lie at (0,0). *)
val cx1cy2 = cx1 * cy2
and cx2cy1 = cx2 * cy1
and cx1cy3m0 = cx1 * cy3m0
and cx2cy3m0 = cx2 * cy3m0
and cx3m0cy1 = cx3m0 * cy1
and cx3m0cy2 = cx3m0 * cy2
val cx1cy2_m_cx2cy1 = cx1cy2 - cx2cy1
and cx1cy3m0_m_cx2cy3m0 = cx1cy3m0 - cx2cy3m0
and cx1cy3m0_m_cx3m0cy1 = cx1cy3m0 - cx3m0cy1
macdef compute_A = g0i2f 3 * cx1cy2_m_cx2cy1
macdef compute_B =
g0i2f 3 * (~cx1cy3m0_m_cx2cy3m0 + cx3m0cy1 - cx3m0cy2
- cx1cy2_m_cx2cy1)
macdef compute_C =
cx1cy3m0_m_cx3m0cy1 + cx1cy3m0_m_cx3m0cy1
+ cx3m0cy2 - cx2cy3m0 + cx1cy2_m_cx2cy1
val A = compute_A
in
if A = g0i2f 0 then
let
val B = compute_B
in
if B = g0i2f 0 then
@(nan, nan)
else
let
val C = compute_C
in
@(~C / B, nan)
end
end
else
let
val B = compute_B and C = compute_C
val discriminant = (B * B) - (g0i2f 4 * A * C)
in
if discriminant < g0i2f 0 then
@(nan, nan)
else
let
val sqrt_discr = sqrt<tk> discriminant
in
@((~B - sqrt_discr) / (A + A),
(~B + sqrt_discr) / (A + A))
end
end
end
fn {a : t@ype} {tk : tkind}
spower_draw_degree3
{w, h : int}
(pix : !pixmap (a, w, h),
cx : @(g0float tk, g0float tk, g0float tk, g0float tk),
cy : @(g0float tk, g0float tk, g0float tk, g0float tk),
color : a)
: void =
(* The algorithm is adaptive flattening by repeated bisection. This
is done in the s-power basis. (The algorithm of de Casteljau is
NOT used.) First, though, the curve is broken at its inflection
points, mainly in case such a point be a cusp. *)
let
extern castfn sz2lli : {i : int} size_t i -<> llint i
extern castfn dbl2lli : double -<> [i : int] llint i
extern castfn f2dbl : g0float tk -<> double
fn {}
round_to_lli (x : g0float tk) :<> [i : int] llint i =
let
val x = f2dbl x
val x = $extfcall (double, "rintl", x)
in
dbl2lli x
end
#define TOL 0.25 (* QUESTION: What is a good value for TOL? *)
fun
flatten_and_draw
{p : addr}
(pix : !pixmap (a, w, h, p),
cx : @(g0float tk, g0float tk,
g0float tk, g0float tk),
cy : @(g0float tk, g0float tk,
g0float tk, g0float tk),
iter : bisect_iter)
: void =
if ~bisect_iter_done iter then
let
val interv = bisect_iter_interval<tk> iter
val dx = spower_portion_degree3<tk> (cx, interv)
and dy = spower_portion_degree3<tk> (cy, interv)
in
if abs dx.1 + abs dx.2 < g0f2f TOL
&& abs dy.1 + abs dy.2 < g0f2f TOL then
let
val w = sz2lli (width pix) and h = sz2lli (height pix)
val dx0 = round_to_lli dx.0
and dy0 = round_to_lli dy.0
and dx3 = round_to_lli dx.3
and dy3 = round_to_lli dy.3
in
if (0LL <= dx0) * (dx0 < w)
* (0LL <= dy0) * (dy0 < h)
* (0LL <= dx3) * (dx3 < w)
* (0LL <= dy3) * (dy3 < h) then
draw_line<a><llintknd>
(pix, @(dx0, dy0), @(dx3, dy3), color);
flatten_and_draw (pix, cx, cy, bisect_iter_next iter)
end
else
flatten_and_draw (pix, cx, cy, bisect_iter_bisect iter)
end
val @(t1, t2) = spower_inflections_degree3 (cx, cy)
macdef zero = g0i2f 0
macdef one = g0i2f 1
in
if zero <= t1 && t1 <= one then
begin
if t1 = t2 then
let
val cx_0 = spower_portion_degree3<tk> (cx, @(zero, t1))
and cy_0 = spower_portion_degree3<tk> (cy, @(zero, t1))
and cx_1 = spower_portion_degree3<tk> (cx, @(t1, one))
and cy_1 = spower_portion_degree3<tk> (cy, @(t1, one))
in
flatten_and_draw (pix, cx_0, cy_0, bisect_iter_make ());
flatten_and_draw (pix, cx_1, cy_1, bisect_iter_make ())
end
else if zero <= t2 && t2 <= one then
let
val t1 = min (t1, t2)
and t2 = max (t1, t2)
val cx_0 = spower_portion_degree3<tk> (cx, @(zero, t1))
and cy_0 = spower_portion_degree3<tk> (cy, @(zero, t1))
and cx_1 = spower_portion_degree3<tk> (cx, @(t1, t2))
and cy_1 = spower_portion_degree3<tk> (cy, @(t1, t2))
and cx_2 = spower_portion_degree3<tk> (cx, @(t2, one))
and cy_2 = spower_portion_degree3<tk> (cy, @(t2, one))
in
flatten_and_draw (pix, cx_0, cy_0, bisect_iter_make ());
flatten_and_draw (pix, cx_1, cy_1, bisect_iter_make ());
flatten_and_draw (pix, cx_2, cy_2, bisect_iter_make ())
end
end
else if zero <= t2 && t2 <= one then
let
val cx_0 = spower_portion_degree3<tk> (cx, @(zero, t2))
and cy_0 = spower_portion_degree3<tk> (cy, @(zero, t2))
and cx_1 = spower_portion_degree3<tk> (cx, @(t2, one))
and cy_1 = spower_portion_degree3<tk> (cy, @(t2, one))
in
flatten_and_draw (pix, cx_0, cy_0, bisect_iter_make ());
flatten_and_draw (pix, cx_1, cy_1, bisect_iter_make ())
end
else
flatten_and_draw (pix, cx, cy, bisect_iter_make ())
end
implement {a} {tk}
pixmap_draw_bezier_degree3_control_points (pix, pt0, pt1, pt2, pt3,
color) =
let
(* Bernstein coefficients. *)
val bx = @(pt0.0, pt1.0, pt2.0, pt3.0)
and by = @(pt0.1, pt1.1, pt2.1, pt3.1)
(* s-power coefficients. *)
val cx = bernstein2spower_degree3<tk> bx
and cy = bernstein2spower_degree3<tk> by
in
spower_draw_degree3<a><tk> (pix, cx, cy, color)
end
(*------------------------------------------------------------------*)
#ifdef BITMAP_BRESENHAM_TASKS_TEST #then
staload "bitmap_write_ppm_task.sats"
staload _ = "bitmap_write_ppm_task.dats"
fn
draw_some_lines () : void =
let
val @(pfgc | pix) =
pixmap_make_elt<rgb24> (i2sz 200, i2sz 100,
rgb24_make (217, 217, 214))
val () = draw_line (pix, @(0, 0), @(199, 99),
rgb24_make (39, 93, 56))
val () = draw_line (pix, @(0, 99), @(199, 0),
rgb24_make (0, 51, 160))
val () = draw_line (pix, @(90, 0), @(110, 99),
rgb24_make (125, 85, 199))
val () = draw_line (pix, @(110, 0), @(90, 99),
rgb24_make (137, 60, 71))
val () = draw_line (pix, @(199, 89), @(0, 10),
rgb24_make (0, 134, 191))
val () = draw_line (pix, @(199, 20), @(0, 79),
rgb24_make (89, 203, 232))
val () = draw_line (pix, @(99, 79), @(119, 30),
rgb24_make (235, 111, 189))
val () = draw_line (pix, @(99, 79), @(79, 30),
rgb24_make (179, 155, 0))
val () = draw_line (pix, @(0, 10), @(199, 10),
rgb24_make (200, 130, 66))
val () = draw_line (pix, @(199, 90), @(0, 90),
rgb24_make (187, 166, 0))
val () = draw_line (pix, @(10, 0), @(10, 99),
rgb24_make (0, 119, 200))
val () = draw_line (pix, @(190, 99), @(190, 0),
rgb24_make (92, 78, 99))
val () = draw_line (pix, @(1199, ~10), @(~1000, 109),
rgb24_make (162, 0, 103))
val () = draw_line (pix, @(1199, 109), @(~1000, ~10),
rgb24_make (244, 54, 76))
val outf = fileref_open_exn ("some_lines.ppm", file_mode_w)
in
ignoret (pixmap_write_ppm<rgb24> (outf, pix));
fileref_close outf;
free (pfgc | pix)
end
fn
draw_some_circles () : void =
let
val @(pfgc | pix) =
pixmap_make_elt<rgb24> (i2sz 200, i2sz 100,
rgb24_make (215, 210, 203))
val () = draw_circle (pix, @(99, 49), 30, rgb24_make(98, 52, 18))
val () = draw_circle (pix, @(0, 0), 99, rgb24_make(34, 55, 43))
val () = draw_circle (pix, @(160, 55), 50, rgb24_make(0, 85, 140))
val outf = fileref_open_exn ("some_circles.ppm", file_mode_w)
in
ignoret (pixmap_write_ppm<rgb24> (outf, pix));
fileref_close outf;
free (pfgc | pix)
end
fn
test_portion_degree2 () : void =
let
val b = @(3.0, ~5.0, 17.0)
val c = bernstein2spower_degree2<dblknd> b
val t0 = 0.2 and t1 = 0.7
val d = spower_portion_degree2<dblknd> (c, @(t0, t1))
var t : double
in
for (t := 0.0; t <= 1.0; t := t + 0.025)
let
val f = spower_eval_degree2<dblknd> (c, t0 + ((t1 - t0) * t))
val g = spower_eval_degree2<dblknd> (d, t)
val- true = abs (f - g) < 0.00000001
in
end
end
fn
test_portion_degree3 () : void =
let
val b = @(3.0, ~5.0, 17.0, ~70.0)
val c = bernstein2spower_degree3<dblknd> b
val t0 = 0.2 and t1 = 0.7
val d = spower_portion_degree3<dblknd> (c, @(t0, t1))
var t : double
in
for (t := 0.0; t <= 1.0; t := t + 0.025)
let
val f = spower_eval_degree3<dblknd> (c, t0 + ((t1 - t0) * t))
val g = spower_eval_degree3<dblknd> (d, t)
val- true = abs (f - g) < 0.00000001
in
end
end
fn
test_inflections () : void =
begin
let
(* Some of the following test examples are from the paper
referenced at https://archive.ph/1mUxT *)
(* This should have a cusp at t = 0.5. *)
val bx = @(0.0, 1.0, 0.0, 1.0)
and by = @(0.0, 1.0, 1.0, 0.0)
val cx = bernstein2spower_degree3<dblknd> bx
and cy = bernstein2spower_degree3<dblknd> by
val @(t1, t2) = spower_inflections_degree3<dblknd> (cx, cy)
val- true = abs (t1 - 0.5) < 0.00000001
val- true = abs (t2 - 0.5) < 0.00000001
in
end;
let
(* This should have an inflection at t = 0.5 *)
val bx = @(0.0, 1.0, 2.0, 3.0)
and by = @(1.0, 1.0, ~1.0, ~1.0)
val cx = bernstein2spower_degree3<dblknd> bx
and cy = bernstein2spower_degree3<dblknd> by
val @(t1, t2) = spower_inflections_degree3<dblknd> (cx, cy)
val- true = abs (t1 - 0.5) < 0.00000001
in
end;
let
(* This should have one inflection point in [0,1]. *)
val bx = @(16.0, 185.0, 673.0, 810.0)
and by = @(467.0, 95.0, 545.0, 17.0)
val cx = bernstein2spower_degree3<dblknd> bx
and cy = bernstein2spower_degree3<dblknd> by
val @(t1, t2) = spower_inflections_degree3<dblknd> (cx, cy)
val- true = abs (t1 - 0.456590) < 0.00001
val- true = abs (t2 - ~24.047383) < 0.00001
in
end;
let
(* This should have two inflections in [0,1]. *)
val bx = @(859.0, 13.0, 781.0, 266.0)
and by = @(676.0, 422.0, 12.0, 425.0)
val cx = bernstein2spower_degree3<dblknd> bx
and cy = bernstein2spower_degree3<dblknd> by
val @(t1, t2) = spower_inflections_degree3<dblknd> (cx, cy)
val- true = abs (t1 - 0.681076) < 0.00001
val- true = abs (t2 - 0.705299) < 0.00001
in
end;
let
(* This should have two inflections in [0,1]. *)
val bx = @(872.0, 11.0, 779.0, 220.0)
and by = @(686.0, 423.0, 13.0, 376.0)
val cx = bernstein2spower_degree3<dblknd> bx
and cy = bernstein2spower_degree3<dblknd> by
val @(t1, t2) = spower_inflections_degree3<dblknd> (cx, cy)
val- true = abs (t1 - 0.588071) < 0.00001
val- true = abs (t2 - 0.886863) < 0.00001
in
end;
let
(* This should have two inflections in [0,1]. *)
val bx = @(819.0, 43.0, 826.0, 25.0)
and by = @(566.0, 18.0, 18.0, 533.0)
val cx = bernstein2spower_degree3<dblknd> bx
and cy = bernstein2spower_degree3<dblknd> by
val @(t1, t2) = spower_inflections_degree3<dblknd> (cx, cy)
val- true = abs (t1 - 0.476169) < 0.00001
val- true = abs (t2 - 0.539295) < 0.00001
in
end;
let
(* This should have two inflections in [0,1]. *)
val bx = @(819.0, 43.0, 826.0, 25.0)
and by = @(566.0, 18.0, 18.0, 533.0)
val cx = bernstein2spower_degree3<dblknd> bx
and cy = bernstein2spower_degree3<dblknd> by
val @(t1, t2) = spower_inflections_degree3<dblknd> (cx, cy)
val- true = abs (t1 - 0.476169) < 0.00001
val- true = abs (t2 - 0.539295) < 0.00001
in
end;
let
(* This should have two inflections in [0,1]. *)
val bx = @(884.0, 135.0, 678.0, 14.0)
and by = @(574.0, 14.0, 14.0, 566.0)
val cx = bernstein2spower_degree3<dblknd> bx
and cy = bernstein2spower_degree3<dblknd> by
val @(t1, t2) = spower_inflections_degree3<dblknd> (cx, cy)
val- true = abs (t1 - 0.320836) < 0.00001
val- true = abs (t2 - 0.682291) < 0.00001
in
end;
end
fn
draw_some_quadratics () : void =
let
val @(pfgc | pix) =
pixmap_make_elt<rgb24> (i2sz 200, i2sz 100,
rgb24_make (177, 228, 227))
val () = draw_bezier (pix, @(0.0, 49.0), @(0.0, 99.0),
@(99.0, 99.0), rgb24_make (166, 9, 61))
val () = draw_bezier (pix, @(99.0, 99.0), @(199.0, 99.0),
@(199.0, 49.0), rgb24_make (234, 39, 194))
val () = draw_bezier (pix, @(199.0, 49.0), @(199.0, 0.0),
@(99.0, 0.0), rgb24_make (113, 92, 42))
val () = draw_bezier (pix, @(99.0, 0.0), @(0.0, 0.0),
@(0.0, 49.0), rgb24_make (0, 83, 76))
val outf = fileref_open_exn ("some_quadratics.ppm", file_mode_w)
in
ignoret (pixmap_write_ppm<rgb24> (outf, pix));
fileref_close outf;
free (pfgc | pix)
end
fn
draw_some_cubics () : void =
let
val @(pfgc | pix) =
pixmap_make_elt<rgb24> (i2sz 200, i2sz 100,
rgb24_make (197, 207, 218))
(* This has a cusp. *)
val () = draw_bezier (pix, @(0.0, 0.0), @(99.0, 99.0),
@(0.0, 99.0), @(99.0, 0.0),
rgb24_make (166, 9, 61))
(* This self-intersects. *)
val () = draw_bezier (pix, @(100.0, 0.0), @(249.0, 99.0),
@(50.0, 99.0), @(199.0, 0.0),
rgb24_make (0, 81, 81))
(* This is a closed curve. *)
val () = draw_bezier (pix, @(50.0, 50.0), @(250.0, 20.0),
@(0.0, 90.0), @(50.0, 50.0),
rgb24_make (0, 83, 76))
val () = draw_bezier (pix, @(0.0, 0.0), @(49.0, 99.0),
@(149.0, 0.0), @(199.0, 99.0),
rgb24_make (113, 92, 42))
val () = draw_bezier (pix, @(199.0, 0.0), @(0.0, 0.0),
@(250.0, 200.0), @(0.0, 0.0),
rgb24_make (111, 44, 63))
(* This one is only partly inside the window. *)
val () = draw_bezier (pix, @(~10.0, 0.0), @(110.0, 0.0),
@(110.0, 99.0), @(~10.0, 99.0),
rgb24_make (111, 44, 63))
val outf = fileref_open_exn ("some_cubics.ppm", file_mode_w)
in
ignoret (pixmap_write_ppm<rgb24> (outf, pix));
fileref_close outf;
free (pfgc | pix)
end
implement
main0 () =
begin
draw_some_lines ();
draw_some_circles ();
test_portion_degree2 ();
test_portion_degree3 ();
test_inflections ();
draw_some_quadratics ();
draw_some_cubics ()
end
#endif
(*------------------------------------------------------------------*)- Output:
To compile and run the test program:
$ patscc -std=gnu2x -g -O2 -DATS_MEMALLOC_LIBC -DATS BITMAP_BRESENHAM_TASKS_TEST bitmap_{bresenham_tasks,{,write_ppm_}task}.{s,d}ats -lm
$ ./a.out
The program should create Portable Pixel Map files some_lines.ppm, some_circles.ppm, some_quadratics.ppm, and some_cubics.ppm:
A Maxima script I used
It might be useful if I posted the Maxima script I used to help me devise some of the Bézier curve code. The Maxima commands are not organized well, nor documented well, but nevertheless they might prove useful:
/* sbern <--> spow */
T1i : matrix(
[1, 0],
[0, 1]);
T1 : invert(T1i);
T2i : matrix(
[1, 0, 0],
[1, 1, 1],
[0, 0, 1]);
T2 : invert(T2i);
T3i : matrix(
[1, 0, 0, 0],
[2, 1, 0, 1],
[1, 0, 1, 2],
[0, 0, 0, 1]);
T3 : invert(T3i);
expand(T1.matrix([c0],[c1]));
expand(T2.matrix([c0],[c1],[c2]));
expand(T3.matrix([c0],[c1],[c2],[c3]));
/* sbern <--> mono */
S1i : transpose(matrix(
[1, -1],
[0, 1]));
S1 : invert(S1i);
S2i : transpose(matrix(
[1, -2, 1],
[0, 1, -1],
[0, 0, 1]));
S2 : invert(S2i);
S3i : transpose(matrix(
[1, -3, 3, -1],
[0, 1, -2, 1],
[0, 0, 1, -1],
[0, 0, 0, 1]));
S3 : invert(S3i);
F1 : expand(S1i.T1i.matrix([c0],[c1]));
F2 : expand(S2i.T2i.matrix([c0],[c1],[c2]));
F3 : expand(S3i.T3i.matrix([c0],[c1],[c2],[c3]));
/* portion t0..t1 */
g(x) := t0 + ((t1 - t0) * x);
f2(x) := expand(F2.matrix([1],[x],[x**2]));
expand(T2.S2.matrix(
[coeff(f2(g(t)), t, 0)],
[coeff(f2(g(t)), t, 1)],
[coeff(f2(g(t)), t, 2)]));
f3(x) := expand(F3.matrix([1],[x],[x**2],[x**3]));
Solution3 : expand(T3.S3.matrix(
[coeff(f3(g(t)), t, 0)],
[coeff(f3(g(t)), t, 1)],
[coeff(f3(g(t)), t, 2)],
[coeff(f3(g(t)), t, 3)]));
Solution3[1];
Solution3[2];
Solution3[3];
Solution3[4];
/* derivatives and critical points */
/* In s-power format: */
Mx : S3i.T3i.matrix([cx0],[cx1],[cx2],[cx3]);
My : S3i.T3i.matrix([cy0],[cy1],[cy2],[cy3]);
fx(t) := expand(transpose(Mx).matrix([1],[t],[t**2],[t**3]));
fy(t) := expand(transpose(My).matrix([1],[t],[t**2],[t**3]));
fxp(t) := diff(fx(t), t);
fxpp(t) := diff(fxp(t), t);
fyp(t) := diff(fy(t), t);
fypp(t) := diff(fyp(t), t);
cross(t) := expand((fxp(t) * fypp(t)) - (fyp(t) * fxpp(t)));
expand(fxp(t));
expand(fyp(t));
expand(cross(t) / 2); /* Notice that this is quadratic. */
/* Translating the first endpoint to (0,0), and the second to
cx3 - cx0. */
Mx : S3i.T3i.matrix([0],[cx1],[cx2],[cx3m0]);
My : S3i.T3i.matrix([0],[cy1],[cy2],[cy3m0]);
fx(t) := expand(transpose(Mx).matrix([1],[t],[t**2],[t**3]));
fy(t) := expand(transpose(My).matrix([1],[t],[t**2],[t**3]));
fxp(t) := diff(fx(t), t);
fxpp(t) := diff(fxp(t), t);
fyp(t) := diff(fy(t), t);
fypp(t) := diff(fyp(t), t);
cross(t) := expand((fxp(t) * fypp(t)) - (fyp(t) * fxpp(t)));
expand(fxp(t));
expand(fyp(t));
expand(cross(t) / 2); /* Notice that this is quadratic. */
Some of the calculations are done by assigning names to the coefficients of splines in the s-power basis, finding the corresponding coefficients in the ordinary monomial basis, then (using Maxima's inherent capabilities) performing an operation such as composition or differentiation. It is a sloppy script.