(* Bresenham line drawing algorithm. *)
module M
use import int.Int
use import ref.Ref
(* Parameters.
Without loss of generality, we can take [x1=0] and [y1=0].
Thus the line to draw joins [(0,0)] to [(x2,y2)]
and we have [deltax = x2] and [deltay = y2].
Moreover we assume being in the first octant, i.e.
[0 <= y2 <= x2]. The seven other cases can be easily
deduced by symmetry. *)
constant x2: int
constant y2: int
axiom first_octant: 0 <= y2 <= x2
(* [best x y] expresses that the point [(x,y)] is the best
possible point i.e. the closest to the real line
i.e. for all y', we have |y - x*y2/x2| <= |y' - x*y2/x2|
We stay in type [int] by multiplying everything by [x2]. *)
use import int.Abs
predicate best (x y: int) =
forall y': int. abs (x2 * y - x * y2) <= abs (x2 * y' - x * y2)
(** Key lemma for Bresenham's proof: if [b] is at distance less or equal
than [1/2] from the rational [c/a], then it is the closest such integer.
We express this property using integers by multiplying everything by [2a]. *)
lemma closest :
forall a b c: int.
abs (2 * a * b - 2 * c) <= a ->
forall b': int. abs (a * b - c) <= abs (a * b' - c)
let bresenham () =
let y = ref 0 in
let e = ref (2 * y2 - x2) in
for x = 0 to x2 do
invariant { !e = 2 * (x + 1) * y2 - (2 * !y + 1) * x2 }
invariant { 2 * (y2 - x2) <= !e <= 2 * y2 }
(* here we would plot (x, y),
so we assert this is the best possible row y for column x *)
assert { best x !y };
if !e < 0 then
e := !e + 2 * y2
else begin
y := !y + 1;
e := !e + 2 * (y2 - x2)
end
done
end