Wiki Agenda Contact Version française

Bresenham line drawing algorithm

We prove optimality: for each column x, the point (x,y) which is drawn is the closest to the rational line between (0,0) and (x2,y2).

Bresenham algorithm is interesting because its code only uses linear arithmetic but its proof of optimality requires non-linear arithmetic.


Authors: Jean-Christophe Filliâtre

Topics: Non-linear Arithmetic

Tools: Why3

see also the index (by topic, by tool, by reference, by year)


(* Bresenham line drawing algorithm. *)

module M

  use int.Int
  use 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. *)

  (* `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 int.Abs

  predicate best (x2 y2 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 (x2 y2:int)
    requires { 0 <= y2 <= x2 }
  = 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 x2 y2 x !y };
      if !e < 0 then
        e := !e + 2 * y2
      else begin
        y := !y + 1;
        e := !e + 2 * (y2 - x2)
      end
    done

end

download ZIP archive