Im ersten Artikel hatten wir den
SMT-Solver Z3 kennen gelernt. Wir hatten in der zugehörigen Sprache
SMT-LIB-2 eine Spezifikation für die Rasterisierung von Geraden
verfasst und dann einen verbesserten Algorithmus entworfen, der mit
dem simplen in seiner Funktionsweise übereinstimmt. Während die
Korrespondenz dieser beiden Algorithmen manchen Lesern vielleicht noch
offensichtlich erschien, gehen wir in diesem Artikel einen signifikanten
Schritt weiter. Wir implementieren den
sog. Bresenham-Algorithmus,
bei dem auf den ersten Blick ganz und gar nicht mehr klar ist, warum
er überhaupt funktioniert.
Der nächste bitte
Der zweite Algorithmus aus unserem ersten Artikel ist zwar nicht mehr
ganz so verschwenderisch wie stures Ausmultiplizieren, aber er führt
im Algorithmenzustand immer noch eine reelle (rationale) Zahl y mit, welche in
jedem Schritt gerundet werden muss. Das ist immer noch recht aufwendig
und auch das ist recht unnötig. Wir haben uns am Anfang bewusst
eingeschränkt auf „flache“ (dx < dy), „aufsteigende“ (dy >= 0)
Linien. Bei der Rasterung solcher Linien geht‘s für jeden x-Schritt in
y-Richtung entweder einen Schritt nach oben, oder wir bleiben auf der
gleichen Höhe. Das klingt zum einen plausibel, ist andererseits aber
auch beweisbar:
(define-fun candidate-next-pixel-top ((st State-1))
Int
(+ 1 (state-1-y st)))
(define-fun candidate-next-pixel-bottom ((st State-1))
Int
(state-1-y st))
(push)
(declare-const st State-1)
(assert (state-1-valid st))
(assert (not (let ((next-st (step-1 st)))
(or
(= (candidate-next-pixel-bottom st)
(state-1-y next-st))
(= (candidate-next-pixel-top st)
(state-1-y next-st))))))
(check-sat)
(pop)
Wir benutzen hier und im folgenden die Ausdrücke (push) und (pop),
um Zwischenergebnisse zu prüfen, und zwar ohne, dass das einen
Einfluss auf das Endergebnis hat. Alles, was seit dem letzten
(push) deklariert oder assertet wurde, ist nach dem (pop) wieder
weg.
Die Entscheidung, welchen dieser beiden
Kandidatenpixel wir auswählen, treffen wir anhand der Abstände der
jeweiligen Mittelpunkte zur Ideallinie. Dieser Zusammenhang lässt sich
am besten anhand einer Grafik illustrieren:
Die Abstände zum oberen Pixel d-oben und zum unteren Pixel d-unten
werden wir später als expliziten Teil des Zustands mitführen. Jetzt –
in der Phase, wo wir uns noch klar werden müssen über alle
Zusammenhänge – können wir die beiden Abstände mit Funktionen
ableiten.
(define-fun state-1-d-oben ((st State-1))
Real
(-
(candidate-next-pixel-top st)
(state-1-exact-y (step-1 st))))
(define-fun state-1-d-unten ((st State-1))
Real
(-
(state-1-exact-y (step-1 st))
(candidate-next-pixel-bottom st)))
Hierbei ist wichtig, dass sich die beiden Abstände d-unten und
d-oben auf den nächsten Zustand beziehen. Deshalb müssen wir zur
Berechnung einmal mit step-1 weiterschalten.
Mit diesen Definitionen können wir prüfen, ob das Wissen über die beiden
Abstände ausreicht, um die Entscheidung zu treffen, ob wir den oberen
Kandidatenpixel auswählen müssen oder den unteren. Eine Entscheidung in
SMT-LIB-2 treffen wir mit dem ite-Ausdruck, was für „if-then-else“
steht:
(push)
(declare-const st State-1)
(assert (state-1-valid st))
(assert (not
(let ((next-y (state-1-y
(step-1 st)))
(dx (line-dx (state-1-line st)))
(dy (line-dy (state-1-line st)))
(x (state-1-x st))
(y (state-1-y st)))
(ite (< (state-1-d-oben st)
(state-1-d-unten st))
(= next-y
(candidate-next-pixel-top st))
(= next-y
(candidate-next-pixel-bottom st))))))
(check-sat)
(pop)
Z3 sagt hierzu sat. Das ist schlecht; unsere Aussage gilt nicht. Um
herauszufinden, was wir falsch gemacht haben, können wir einen
get-value-Aufruf hinter das (check-sat) schreiben. get-value
nimmt als Argument eine Liste mit Ausdrücken und druckt deren Werte
aus für das gefundene Modell (sprich: die gefundene
Variablebelegung). Wir lassen uns ein paar interessante Werte
ausdrucken.
(get-value
(st
(state-1-y st)
(state-1-y
(step-1 st))
(state-1-d-oben st)
(state-1-d-unten st)
(candidate-next-pixel-top st)
(candidate-next-pixel-bottom st)
))
Das Ergebnis davon ist:
sat
((st (mk-state-1 (mk-line 4 3) 1))
((state-1-y st) 1)
((state-1-y
(step-1 st)) 2)
((state-1-d-oben st) (/ 1.0 2.0))
((state-1-d-unten st) (/ 1.0 2.0))
((candidate-next-pixel-top st) 2)
((candidate-next-pixel-bottom st) 1))
Wir sehen, dass der Abstand zum oberen Kandidaten mit 1/2 genau
gleich groß ist wie der Abstand zum unteren Kandidaten. In der
Bedingung unseres ite steht ein <, also strikt kleiner, was in
diesem Fall nicht gegeben ist. Dann wählen wir den
unteren Pixel aus. Richtig gewesen wäre der obere, denn bei 1/2
rundet round-to-nearest auf statt ab. Ein naheliegender
Lösungsversuch wäre, das < in ein <= zu ändern. Und siehe da,
damit sagt Z3 unsat.
Die Bedingung für das ite ist jetzt (<= (state-1-d-oben st)
(state-1-d-unten st)), was immer noch die Berechnung des exakten
y-Werts erfordert. Davon wollen wir jetzt weg. Zu diesem Zwecke
formulieren wir die Bedingung nach und nach um. Diese Umformungen
können wir dabei jeweils von Z3 auf Korrektheit prüfen lassen.
Wir substituieren zunächst von Hand alle Funktionsaufrufe durch
deren jeweilige Implementierungen und formen nebenbei die Ungleichung
so um, dass links nur noch die Null steht.
(push)
(declare-const st State-1)
(assert (state-1-valid st))
(assert (not
(=
(<= (state-1-d-oben st)
(state-1-d-unten st))
(<= 0
(-
(* 2 (line-y (state-1-line (step-1 st))
(state-1-x (step-1 st))))
(* 2 (state-1-y st))
1)))))
(check-sat)
(pop)
Wir wissen außerdem, dass die Linie state-1-line sich nicht
ändert. Darüberhinaus ist das x nach einem step-1 einfach x +
1. Wenn wir zusätzlich noch die Definition von line-y und dann die
von line-slope einsetzen, erhalten wir:
(push)
(declare-const st State-1)
(assert (state-1-valid st))
(assert (not
(let ((next-y (state-1-y
(step-1 st)))
(dx (line-dx (state-1-line st)))
(dy (line-dy (state-1-line st)))
(x (state-1-x st))
(y (state-1-y st)))
(ite (<= 0
(-
(* 2 (* (/ dy dx)
(+ 1 x)))
(* 2 y)
1))
(= next-y
(candidate-next-pixel-top st))
(= next-y
(candidate-next-pixel-bottom st))))))
(check-sat)
(pop)
Die Ungleichung enthält (/ dy dx) als einzigen Bruch. Den können wir
beseitigen, indem wir die Ungleichung mit dx
durchmultiplizieren. Das ist zulässig, so lange dx ungleich 0
ist. Anders als auf Papier müssen wir die Zulässigkeit hier gar nicht
verargumentieren, sondern wir können die Umformulierung einfach machen
und Z3 sagt uns dann, ob das zulässig war. Die Ungleichung sieht jetzt so aus:
(<= 0
(-
(* 2 dy (+ 1 x))
(* 2 y dx)
dx))
Z3 sagt dazu immer noch unsat. Die Einschränkung (state-1-valid st)
enthält die Einschränkung, dass die im Zustand befindliche Linie
valide sein soll. In line-valid steht (> (line-dx l) 0). Wir
können mal probeweise diese Einschränkung löschen und Z3 fragen, ob
dann unsere Aussagen immer noch gelten. Erwartungsgemäß sagt Z3 jetzt
sat.
Der Ausdruck (- (* 2 (* dy (+ 1 x))) (* 2 y dx) dx) ist jetzt eine
Ganzzahl und zu deren Berechnung müssen wir auch nur
Ganzzahlarithmetik betreiben. Nur die Multiplikationen stören uns
noch. Aber um eine Multiplikation potenziell aufzulösen kennen wir ja
schon einen Trick: Wir merken uns den jeweiligen Wert aus dem
vorherigen Schritt und addieren für den nächsten Schritt was drauf
(oder ziehen was ab). Das sollte hier funktionieren, denn
das Produkt (* 2 dy (+ 1 x)) erhöht sich in jedem Schritt um genau
2 dy und das Produkt (* 2 y dx) erhöht sich je nachdem, ob wir den
unteren Kandidatenpixel oder den oberen auswählen entweder gar nicht
oder um ein konstantes (* 2 dx). Mit diesen formal verifizierten
Überlegungen sind wir gewappnet für den finalen Schlag: Wir bauen
einen Algorithmus, der eine Linie zeichnet und sich dabei nur
einfacher Ganzzahlarithmetik bedient.
I keep my ends out for the tie that binds
Wir machen wieder das übliche Schema. Zuerst definieren wir ein Zustandsobjekt. Das enthält neben dem üblichen Kram jetzt unser d, welches unsere Entscheidung für den oberen oder unteren nächsten Pixel treibt:
(declare-datatypes ()
((State-4 (mk-state-4 (state-4-dx Int)
(state-4-dy Int)
(state-4-x Int)
(state-4-y Int)
(state-4-d Int)
))))
Ein solcher Zustand ist valide, wenn die üblichen Einschränkungen gelten und zusätzlich unser d dem Wert entspricht, der uns vorschwebt:
(define-fun state-4-m ((st State-4))
Real
(/ (state-4-dy st)
(state-4-dx st)))
(define-fun state-4-valid ((st State-4))
Bool
(let ((dx (state-4-dx st))
(dy (state-4-dy st))
(m (state-4-m st))
(x (state-4-x st))
(y (state-4-y st))
(d (state-4-d st)))
(and
(< 0 dx)
(<= 0 dy dx)
(>= x 0)
;; y is actually the correct y for m and x
(= y (round-to-nearest (* m x)))
;; d is what we want it to describe
(= d
(- (* 2 dy (+ 1 x))
(* 2 y dx)
dx))
)))
Fehlt noch die Funktion für den Einstieg:
(define-fun into-state-4 ((st2 State-1))
State-4
(let ((dx (line-dx (state-1-line st2)))
(dy (line-dy (state-1-line st2)))
(x (state-1-x st2)))
(let ((y (round-to-nearest (* (/ dy dx) x))))
(mk-state-4
dx
dy
x
y
(- (* 2 dy (+ 1 x))
(* 2 y dx)
dx)))))
Und die Schrittfunktion:
(define-fun step-4 ((st State-4))
State-4
(let ((dx (state-4-dx st))
(dy (state-4-dy st))
(x (state-4-x st))
(y (state-4-y st))
(d (state-4-d st)))
(ite (<= 0 d)
;; Diagonalschritt
(mk-state-4
dx
dy
(+ x 1)
(+ y 1)
(+ d (- (* 2 dy)
(* 2 dx))))
;; Rechtsschritt
(mk-state-4
dx
dy
(+ x 1)
y
(+ d (* 2 dy))))))
Als Korrektheitsaussage würden wir gern wieder dem oben etablierten Schema entsprechen:
(declare-const st State-1)
(assert (state-1-valid st))
(assert
(not
(= (step-4
(into-state-4 st))
(into-state-4
(step-1 st)))))
Das terminiert leider nicht. Das ist leider gar nicht schön und zeigt
die Grenzen eines SMT-Solvers als Verifikationsinstrument. Wir können mit etwas Gefrickel aber eine Korrektheitsaussage hinbasteln, die äquivalent ist:
(declare-const before-1 State-1)
(assert (state-1-valid before-1))
(assert (not
(let ((before-4 (into-state-4 before-1)))
(let ((after-4 (step-4 before-4))
(after-1 (step-1 before-1)))
(and
(= (state-4-x after-4)
(state-1-x after-1))
(state-4-valid after-4)
(= (state-1-y after-1)
(state-4-y after-4))
)))))
(check-sat)
Die Bedingungen zu x und y bewirken, dass step-4 dasselbe
Ergebnis liefert wie step-1 und die Bedingung
(step-4-valid after-4) bewirkt, dass wir auch nicht den Pfad
der Tugend verlassen.
Because you‘re mine, I walk the line
Unsere finale Formulierung ist eine Instanz des
sog. Bresenham-Algorithmus. Wir haben diesen hier nicht nur korrekt
implementiert, sondern diese Korrektheit auch sichergestellt und
außerdem auf dem Weg dorthin einen Einblick erhalten, warum der
Algorithmus funktioniert.
Um die Korrektheit zu zeigen, formulierten wir Bedingungen mit
Rückgriff auf die reellen Zahlen – nicht Float, nicht Double, sondern
präzise reelle Zahlen. Das ist richtig, denn nur mit reellen Zahlen
lässt sich die Rasterisierung einer Geraden akkurat
beschreiben. Das bedeutet jedoch nicht, dass eine Implementierung
ebenfalls mit reellen Zahlen arbeiten muss. Das wäre für einen
Computer sehr aufwendig und mit dem Bresenham-Algorithmus haben wir
gezeigt, dass es ja auch ohne geht.
Im nächsten Artikel zeigen wir eine Implementierung des
Bresenham-Algorithmus in Haskell. Wir nutzen Liquid Haskell,
um sicherzustellen, dass auch diese Haskell-Implementierung korrekt im
Sinne unserer obigen Überlegungen ist.