bresenham.mlw 1.68 KB
Newer Older
Jean-Christophe Filliâtre's avatar
Jean-Christophe Filliâtre committed
1 2
(* Bresenham line drawing algorithm. *)

3
module M
4

5
  use import int.Int
6
  use import ref.Ref
7

8 9
  (*  Parameters.
      Without loss of generality, we can take [x1=0] and [y1=0].
10
      Thus the line to draw joins [(0,0)] to [(x2,y2)]
11
      and we have [deltax = x2] and [deltay = y2].
12 13 14 15
      Moreover we assume being in the first octant, i.e.
      [0 <= y2 <= x2]. The seven other cases can be easily
      deduced by symmetry. *)

16 17
  constant x2: int
  constant y2: int
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
18
  axiom first_octant: 0 <= y2 <= x2
19

20 21 22 23
  (* [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]. *)
24

25
  use import int.Abs
26

Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
27 28
  predicate best (x y: int) =
    forall y': int. abs (x2 * y - x * y2) <= abs (x2 * y' - x * y2)
29

30 31 32
  (** 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]. *)
33

34
  lemma closest :
35
    forall a b c: int.
36 37
    abs (2 * a * b - 2 * c) <= a ->
    forall b': int. abs (a * b - c) <= abs (a * b' - c)
38

39 40 41
  let bresenham () =
    let y = ref 0 in
    let e = ref (2 * y2 - x2) in
42 43 44 45 46 47
    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 };
48 49 50 51 52
      if !e < 0 then
        e := !e + 2 * y2
      else begin
        y := !y + 1;
        e := !e + 2 * (y2 - x2)
53
      end
54
    done
55

56
end