Tuesday, May 30, 2023

Raymarching in Lisp

It turns out there’s a nice functional variation of raytracing called raymarching. The algorithms involved are simple and elegant. The learning curve is shallow and you can generate great looking images without hairy trig or linear algebra.

We’ll follow the example of Georges Seurat and simply compute the color independently for each of myriads of pixels. This is efficiently done in parallel in real time on a GPU, but then you have to use shader language and I want to use Lisp. It is insanely inefficient to do this serially on the CPU in Lisp, but still fast enough to render an image in a couple of seconds.

Imagine you have some scene you want to render. There is a volume of 3-dimensional space the scene occupies. Now imagine we know for every point in 3-dimensional space how far away that point is from the nearest surface. This is a scalar value that can be assigned to any point. It is zero for every point that lies on a surface, positive for points above surfaces, and negative for points below. This is the SDF (Signed Distance Field). The SDF is all we need to know to generate a raytraced image of the scene.

We’ll use the SDF to feel our way through the scene. We’ll start at the tip of the ray we’re tracing. We don’t know where the surface is, but if we consult the SDF, we can determine a distance we can safely extend the ray without hitting any surface. From this new point, we can recur, again stepping along no further than the SDF at this new location permits. One of two things will happen: we either step forever or we converge on a surface.

(defun raymarch (sdf origin direction)
  (let iter ((distance 0)
             (count 0))
    (let* ((position (+ origin (* direction distance)))
           (free-path (funcall sdf position)))
      (if (< free-path +min-distance+)
          position  ;; a hit, a very palpable hit
          (unless (or (> count +max-raymarch-iterations+)
                      (> free-path +max-distance+))
            (iter (+ distance free-path) (+ count 1)))))))

To convert an SDF to a Seurat function, we trace an imaginary ray from our eye, through the screen, and into the scene. The ray origin is at your eye, and we’ll say that is about 3 units in front of the window. The ray will travel 3 units to the screen and hit the window at point (i,j), so the ray direction is (normalize (vector i j 3)). We march along the ray to find if we hit a surface. If we did, we compute the amount of light the camera sees using the Lambert shading model.

(defun sdf->seurat (sdf)
  (let ((eye-position (vector 0 0 -4))
        (light-direction (normalize (vector 20 40 -30))))
    (lambda (i j)
      (let* ((ray-direction (normalize (vector i j 3)))
             (hit (raymarch sdf eye-position ray-direction)))
        (if hit
            (* #(0 1 0) (lambert sdf hit light-direction))
            (vector 0 0 0))))))

Lambert shading is proportional to the angle between the surface and the light falling on it, so we take the dot product of the light direction with the normal to the surface at the point the light hits it. If we know the SDF, we can approximate the normal vector at a point by probing the SDF nearby the point and seeing how it changes.

(defun lambert (sdf hit light-direction)
  (dot (pseudonormal sdf hit) light-direction))

(defun pseudonormal (sdf position)
  (let ((o (funcall sdf position))
        (dsdx (funcall sdf (+ #(0.001 0 0) position)))
        (dsdy (funcall sdf (+ #(0 0.001 0) position)))
        (dsdz (funcall sdf (+ #(0 0 0.001) position))))
      (normalize (vector (- dsdx o) (- dsdy o) (- dsdz o)))))

These are all you need to generate good looking 3-d images from a SDF. Now the SDFs for primitive geometric shapes are pretty simple. Here is the SDF for a sphere.

(defun sdf-sphere (position radius)
  (lambda (vector)
    (- (length (- vector position)) radius)))

and the SDF for the ground plane

(defun sdf-ground (h)
  (lambda (vector)
    (+ (svref vector 1) h)))

Given the SDF for two objects, you can use higher order functions to compose them into a scene. Taking the minimum of two SDFs will give you the union of the shapes. Taking the maximum will give you the intersection of two shapes. Other higher order functions on SDFs can blend two SDFs. This has the effect of morphing the shapes together in the image.

I like this approach to raytracing because the arthimetic is straightforward and obvious. You only need the simplest of vector arithmetic, and you don’t need linear algebra or matrix math to get started (although you’ll want project matrixes later on when you want to move your camera around). I’m more comfortable with recursive functions than 3x3 matrices.

This approach to raytracing is best done on a graphics card. These algorithms are pretty straightforward to code up in shader language, but shader language is fairly primitive and doesn’t have higher order functions or closures. Code written in shader language has to be converted to not use closures and HOFs.


Friday, October 21, 2022

Lisp: Second impression

My first impressions of Lisp were not good. I didn't see how navigating list structure was of any use. It seemed to be just a more cumbersome way of getting at the data.

In fact, my first impressions of computer science were not very positive. I enjoyed hobbyist coding on my TRS-80, but “real” programming was tedious and the proscriptions of “doing it the correct way” took the joy out of it. I explored other options for my major. Fate intervened. Over the next year I realized my calling was EECS, so in my sophomore year I took all the intro courses.

I had heard that the introductory computer science course used Lisp. That was disappointing, but I started hearing things about Lisp that made me think I should take a second look. I learned that Lisp was considered the premier language of MIT's Artificial Intelligence Laboratory. It was invented by hackers and designed to be a programmable programming language that was infinitely customizable. The lab had developed special computers that ran Lisp on the hardware. The OS was even written in Lisp. I wasn't looking forward to car and cdr'ing my way through endless cons cells, but I figured that there had to be more going on.

6.001 was unlike the prior computer courses I had taken. The course was not about how to instruct a computer to perform a task — the course was about expressing ideas as computation. To me, this seemed a much better way to approach computers. Professor Hal Abelson was co-lecturing the course. He said that he chose Lisp as the teaching language because it was easier to express ideas clearly.

Two things stood out to me in the first lecture. Professor Abelson showed the recursive and iterative versions of factorial. Of course I had seen recursive factorial from the earlier course and I knew how it worked. Clearly the iterative version must work the same way. (One of my early hangups about Lisp was all the recursion.) I was suprised to find out that the Lisp system would automatically detect tail recursive cases and turn them into iteration. Evidentally, the makers of Lisp had put some thought into this.

Professor Abelson also demonstrated first class functions. He wrote a procedure that numerically approximates the derivative of a function. He then used that in a generic Newton's method solver. This is all straightforward stuff, but to a newbie like me, I thought it was amazing. In just a few lines of code we were doing simple calculus.

It was a mystery to me how first class functions were implemented, but I could see how they were used in the Newton's method solver. The code Professor Abelson wrote was clear and obvious. It captured the concept of derivatives and iterative improvement concisely, and it effectively computed answers to boot. I had to try it. Right after the lecture I went to lab and started typing examples at the REPL. Sure enough, they worked as advertised. A tail-recursive loop really didn't push any stack. It didn't leak even the tiniest bit of memory, no matter how long the loop. I tried the Newton's method solver to take cube roots. I passed the cube function to the derivative function and the result was a function that was numerically close to the derivative.

Now I was a bit more impressed with Lisp than I was earlier. I wasn't completely sold, but I could see some potential here. I wanted to learn a bit more before I dismissed it entirely. It took me several months to become a Lisp fan. The parenthesis were a small hurdle — it took me a couple of weeks to get the hang of let forms. There was a week or two of navigating cons cells to wade through. But I eventually came to love the language.

My first impression of Lisp was poor. The uselessness of traversing random list structure was unmotivating. My second impression was better. Professor Abelson teaching directly from preprints of S&ICP might have had something to do with it.


Wednesday, October 19, 2022

Lisp: First Impressions

My first exposure to Lisp was in the summer of 1981. I was taking a summer school intro to computers. The course was taught on a PDP-11, and for the first few weeks we programmed in Macro-11 assembly language. For the last couple of weeks they introduced Lisp.

Frankly, I wasn't impressed.

The course started by talking about linked lists and how you could navigate them with car and cdr. We then went on to build more complicated structures like alists and plists. This was an old-fashioned lisp, so we used things like getprop and putprop to set symbol properties.

The subject matter wasn't difficult to understand (though chasing pointers around list structure is error prone). Since we had just been learning Macro-11, it was natural to play with linked list structure in assembly code. We wrote assembly code to look things up in a plist.

My impression was that Lisp was centered around manipulating these rather cumbersome data structures called cons cells. Linked lists of cons cells have obvious disadvantages when compared to arrays. This makes the language tedious to work with.

The summer school course was my first “real” college course in computers. I was put off. “Real” computing wasn't as much fun as I had hoped it would be. I definitely wouldn't be considering it as a major, let alone a career. I wasn't interested in Lisp at all.

to be continued

Wednesday, September 28, 2022

Observationally Functional

A couple of years back I wrote a Java microservice that talks to Jenkins to find the state of a series of builds. The code was structured to be “observationally functional” — there were plenty of side effects, but the main data abstractions behaved as if they were immutable as far as the abstract API was concerned. This allows us to treat code that uses these objects as if it were pure functional code.

If a data structure is observationally functional, then regardless of what the implementation does, there is no way to observe side effects at the abstract level. Primarily, this means that if you call a function twice with the same arguments, you always get the same answer. (This implies, but it isn't obvious, that calling a function should not mutate anything that would cause a different function to change.) This restriction has a lot of wiggle room. You can certainly side effect anything local to the abstraction that doesn't get returned to the caller. You can side effect data until the point it is returned to the caller.

The main data abstraction my microservice works with is a representation of the build metadata tree on the Jenkins server. The higher level code walks this tree looking for builds and metadata. The code maintains the illusion that the tree is a local data structure, but the implementation of the tree contains URL references to data that is stored on the Jenkins server. As the higher level code walks the tree, the lower level code fetches the data from the Jenkins server on demand and caches it.

Writing the code this way allows me to separate the data transfer and marshaling parts from the data traversal and analysis part. The tree, though it is mutated as it is traversed, is immutable in the parts that have already been visited. The caching code, which actually mutates the tree, needs to be synchronized across multiple threads, but the traversal code does not. Nodes in the tree that have already been visited are never mutated, so no synchronization is needed.

Once the caching tree abstraction was written, the higher level code simply walks the tree, selecting and filtering nodes, then reading the field values in the nodes. But the higher level code can be treated as if it were pure functional because there are no observable side effects. An advantage of pure functional code is that it is trivially thread safe, so my microservice can run hundreds of threads in parallel, each walking separate parts of the Jenkins tree and none interfering with the other. The only part of the code that uses synchronization is the tree caching code.

This implementation approach was quite fruitful. Once the code was tested with a single thread, it was obvious that multiple threads ought to work (because they couldn't observe each other's side effects) and when I turned the thread count up, no debugging was necessary. The code has been running continuously with dozens of threads for the past couple of years with no timing, synchronization, or race condition bugs.


Wednesday, September 7, 2022

Playing with raycasting

I wanted to learn about raycasting. Raycasting is like a simplified version of ray tracing. As in ray tracing, you examine the environment by projecting a ray from your current location out in some direction to determine what is visible in that direction. But we simplify the problem by only considering two dimensions.

I was also interested in making the graphics more functional and less dependent upon side effects. Now obviously rendering an image to the screen is going to involve side effects, but we can refactor the rendering problem into two subproblems, a pure function that maps the world to an image and the procedure that displays the image.

I'll put the code below. The run procedure implements the event loop state machine. It keeps track of the world and calls next-world on the current world to update the world as time passes. next-world just maps next-state over the objects in the world. next-state does not mutate an object, rather it returns a new object in the new state. Every 13 milliseconds, run calls render-world!, which calls render! on each object in the world.

We're going to use raycasting to fake up a first-person view of a two-dimensional maze. From a position within the maze, we'll cast a ray in a direction and see how far away the wall is. If we peer at the wall through a narrow slit in just that direction, it will appear as a vertical line with height inversely proportional to its distance. If we sweep the ray direction and stack the vertical lines next to each other, it will create three dimensional effect.

The render! method for a fp-view will side effect the screen, but we'll compute the contents functionally. We'll go through each column on the screen and call (vraster fp-view column) to compute a color and a height and we'll draw a vertical line of that height in that color in that column.

(defmethod render! (renderer (fp-view fp-view))
  (dotimes (column +window-width+)
    (multiple-value-bind (r g b height)
        (vraster fp-view column)
      (sdl2:set-render-draw-color renderer r g b #xFF)
      (sdl2:render-draw-line renderer
                             column (round (+ (/ +window-height+ 2) (/ height 2)))
                             column (round (- (/ +window-height+ 2) (/ height 2)))))))

vraster is a function that returns the color and height of the wall on a particular column on the screen. It figures out the angle at which to cast to a ray and calls range to find the distance to the nearest wall at that angle. This is sufficient to determine the wall height for that column, but the first person effect is enhanced significantly if you tint the color according to the distance and the direction of the wall. Knowing the distance, we compute the exact point px, py that the ray hit. It's a wall in the x direction if the y coordinate is an integer and vice versa.

(defparameter +field-of-view+ (/ pi 4))

(defun column->theta (column)
  (- (* (/ column +window-width+) +field-of-view+) (/ +field-of-view+ 2)))

(defun vraster (fp-view column)
  (let* ((location (get-location fp-view))
         (theta (column->theta column))
         (distance (range location theta))
         (px (+ (get-x location) (* (sin (+ (get-theta location) theta)) distance)))
         (py (+ (get-y location) (* (cos (+ (get-theta location) theta)) distance)))
         (wx (< (abs (- py (round py))) 0.05))
         (wy (< (abs (- px (round px))) 0.05)))
    (values
     (min #xFF (floor (/ (if wx #xff #x00) distance)))
     (min #xFF (floor (/ #xFF distance)))
     (min #xFF (floor (/ (if wy #xff #x00) distance)))
     (min +window-height+ (/ (* +window-height+ 2) distance)))))

So we've factored the rendering of a frame into a procedure that draws on the screen and a function that returns what to draw. The direct advantage of this is that we can determine what we should draw without actually drawing it. As an example, suppose we wanted to generate a stereo pair of images. The only thing we need to change is the render! method. It will now compute the view from two slightly different locations and put one set of columns on the left and the other on the right.

(defmethod render! (renderer (fp-view fp-view))
  (dotimes (column (/ +window-width+ 2))
    (multiple-value-bind (r g b height)
        (vraster (left-eye (get-location fp-view)) (* column 2))
      (sdl2:set-render-draw-color renderer r g b #xFF)
      (sdl2:render-draw-line renderer
                             column (round (+ (/ +window-height+ 2) (/ height 2)))
                             column (round (- (/ +window-height+ 2) (/ height 2)))))
    (multiple-value-bind (r g b height)
        (vraster (right-eye (get-location fp-view)) (* column 2))
      (sdl2:set-render-draw-color renderer r g b #xFF)
      (sdl2:render-draw-line renderer
                             (+ column (/ +window-width+ 2)) (round (+ (/ +window-height+ 2) (/ height 2)))
                             (+ column (/ +window-width+ 2)) (round (- (/ +window-height+ 2) (/ height 2)))))))

In this and in a previous post I've gone through the effort of writing some graphics code while avoiding unnecessary side effects. Typical graphics examples and tutorials are stuffed to the brim with global variables, state, and side effects. I wanted to see which side effects were intrinsic to graphics and which are simply incidental to how the examples are coded. It appears that large amounts of the global state and side effects are unnecessary and a more functional approach is reasonable.

As promised, here is the code.

;;; -*- Lisp -*-

(defpackage "RAYCAST"
  (:shadowing-import-from "NAMED-LET" "LET")
  (:use "COMMON-LISP" "NAMED-LET""))

(in-package "RAYCAST")

(defparameter +window-height+ 480)
(defparameter +window-width+ 640)

(defgeneric next-state (object dt)
  (:method ((object t) dt) object))
(defgeneric render! (renderer thing))

(defun next-world (previous-world dt)
  (map 'list (lambda (object) (next-state object dt)) previous-world))

(defun render-world! (renderer world)
  (sdl2:set-render-draw-color renderer #x00 #x00 #x00 #xFF)
  (sdl2:render-clear renderer)
  (mapc (lambda (object) (render! renderer object)) world)
  (sdl2:render-present renderer))

(defun run (initial-world)
  (sdl2:with-init (:video)
    (sdl2:with-window (window
                       :h +window-height+
                       :w +window-width+
                       :flags '(:shown))
      (sdl2:with-renderer (renderer window :index -1 :flags '(:accelerated :presentvsync))

        (let ((last-ticks 0)
              (render-ticker 0)
              (title-ticker 0)
              (sim-count 0)
              (frame-count 0)
              (world initial-world))

          (flet ((title-tick! (dticks)
                   (incf title-ticker dticks)
                   (when (>= title-ticker 1000)
                     (decf title-ticker 1000)
                     (sdl2:set-window-title window
                                            (format nil "Sim rate: ~d, Frame rate: ~d"
                                                    sim-count frame-count))
                     (setq sim-count 0)
                     (setq frame-count 0)))

                 (world-tick! (dticks)
                   (incf sim-count)
                   (setq world (next-world world (/ dticks 1000))))

                 (render-tick! (dticks)
                   (incf render-ticker dticks)
                   (when (>= render-ticker 13)
                     (incf frame-count)
                     (decf render-ticker 13)
                     (render-world! renderer world))))

            (sdl2:with-event-loop (:method :poll)

              (:idle ()
                     (let ((this-ticks (sdl2:get-ticks)))
                       (if (= this-ticks last-ticks)
                           (sdl2:delay 1)
                           (let ((dticks (- this-ticks last-ticks)))
                             (setq last-ticks this-ticks)
                             (title-tick! dticks)
                             (world-tick! dticks)
                             (render-tick! dticks)))))

              (:keydown (:keysym keysym)
                        (case (sdl2:scancode keysym)
                          ((:scancode-x :scancode-escape) (sdl2:push-quit-event))
                          ((:scancode-left :scancode-right
                            :scancode-up :scancode-down
                            :scancode-pageup :scancode-pagedown)
                           nil)
                          (t (format *trace-output* "~&Keydown: ~s" (sdl2:scancode keysym))
                           (force-output *trace-output*))))

              (:quit () t)
              )))))))

(defparameter +maze+
  #2a((1 1 1 1 1 1 1 1 1 1 1 1)
      (1 0 0 1 0 0 0 0 0 0 0 1)
      (1 0 0 1 0 0 0 0 0 0 0 1)
      (1 0 1 1 0 0 0 1 0 1 0 1)
      (1 0 0 0 0 0 0 0 0 0 0 1)
      (1 0 0 0 0 0 0 1 0 1 0 1)
      (1 0 0 0 0 0 0 0 0 0 0 1)
      (1 0 0 0 0 0 0 0 0 0 0 1)
      (1 0 0 0 0 1 0 0 0 1 0 1)
      (1 0 0 0 0 0 0 0 0 0 0 1)
      (1 1 1 1 1 1 1 1 1 1 1 1)))

(defclass location ()
  ((maze :initarg :maze
         :initform +maze+
         :reader get-maze)
   (x :initarg :x
      :initform 2.5
      :reader get-x)
   (y :initarg :y
      :initform 2.5
      :reader get-y)
   (theta :initarg :theta
          :initform 0
          :reader get-theta)))

(defun ud-input ()
  (- (if (sdl2:keyboard-state-p :scancode-up) 1 0)
     (if (sdl2:keyboard-state-p :scancode-down) 1 0)))

(defun lr-input ()
  (- (if (sdl2:keyboard-state-p :scancode-right) 1 0)
     (if (sdl2:keyboard-state-p :scancode-left) 1 0)))

(defun pg-input ()
  (- (if (sdl2:keyboard-state-p :scancode-pageup) 1 0)
     (if (sdl2:keyboard-state-p :scancode-pagedown) 1 0)))

(defun canonicalize-angle (angle)
  (cond ((>= angle pi) (canonicalize-angle (- angle (* pi 2))))
        ((>= angle (- pi)) angle)
        (t (canonicalize-angle (+ angle (* pi 2))))))

(defparameter +translation-rate+ 3.0) ;; tiles per second
(defparameter +rotation-rate+ pi) ;; radians per second

(defmethod next-state ((location location) dt)
  (let ((fbstep (* (ud-input) +translation-rate+ dt))
        (lrstep (* (pg-input) +translation-rate+ dt))
        (thstep (* (lr-input) +rotation-rate+ dt))

        (old-x (get-x location))
        (old-y (get-y location))
        (cos-theta (cos (get-theta location)))
        (sin-theta (sin (get-theta location))))

    (let ((new-x (+ old-x (* sin-theta fbstep) (- (* cos-theta lrstep))))
          (new-y (+ old-y (* cos-theta fbstep) (+ (* sin-theta lrstep))))
          (new-theta (canonicalize-angle (+ (get-theta location) thstep))))
      (cond ((zerop (aref (get-maze location) (floor new-x) (floor new-y)))
             (make-instance 'location :x new-x :y new-y :theta new-theta))
            ((zerop (aref (get-maze location) (floor old-x) (floor new-y)))
             (make-instance 'location :x old-x :y new-y :theta new-theta))
            ((zerop (aref (get-maze location) (floor new-x) (floor old-y)))
             (make-instance 'location :x new-x :y old-y :theta new-theta))
            (t
             (make-instance 'location :x old-x :y old-y :theta new-theta))))))

(defclass fp-view ()
  ((location :initarg :location
             :reader get-location)))

(defmethod next-state ((fp-view fp-view) dt)
  (make-instance 'fp-view :location (next-state (get-location fp-view) dt)))

(defun range (location relative-theta)
  (let* ((angle (+ (get-theta location) relative-theta))

         (dx/dS (sin angle))
         (dy/dS (cos angle))

         (x-step (if (< dx/dS 0) -1 1))
         (y-step (if (< dy/dS 0) -1 1))

         (dS/dx (abs (/ 1 (if (zerop dx/dS) 1e-30 dx/dS))))
         (dS/dy (abs (/ 1 (if (zerop dy/dS) 1e-30 dy/dS)))))

    (let dda ((next-x (* dS/dx
                         (if (< dx/dS 0)
                             (- (get-x location) (floor (get-x location)))
                             (- (+ 1.0 (floor (get-x location))) (get-x location)))))
              (mapx (floor (get-x location)))
              (next-y (* dS/dy
                         (if (< dy/dS 0)
                             (- (get-y location) (floor (get-y location)))
                             (- (+ 1.0 (floor (get-y location))) (get-y location)))))
              (mapy (floor (get-y location)))
              (distance 0))
      (cond ((not (zerop (aref (get-maze location) mapx mapy))) distance)
            ((< next-x next-y)
             (dda (+ next-x dS/dx) (+ mapx x-step)
                  next-y mapy
                  next-x))
            (t
             (dda next-x mapx
                  (+ next-y dS/dy) (+ mapy y-step)
                  next-y))))))

(defparameter +field-of-view+ (/ pi 4))

(defun column->theta (column)
  (- (* (/ column +window-width+) +field-of-view+) (/ +field-of-view+ 2)))

(defun vraster (location column)
  (let* ((theta (column->theta column))
         (distance (range location theta))
         (px (+ (get-x location) (* (sin (+ (get-theta location) theta)) distance)))
         (py (+ (get-y location) (* (cos (+ (get-theta location) theta)) distance)))
         (wx (< (abs (- py (round py))) 0.05))
         (wy (< (abs (- px (round px))) 0.05)))
    (values
     (min #xFF (floor (/ (if wx #xff #x00) distance)))
     (min #xFF (floor (/ #xfF distance)))
     (min #xFF (floor (/ (if wy #xff #x00) distance)))
     (min +window-height+ (/ (* +window-height+ 2) distance)))))

(defmethod render! (renderer (fp-view fp-view))
  (dotimes (column +window-width+)
    (multiple-value-bind (r g b height)
        (vraster (get-location fp-view) column)
      (sdl2:set-render-draw-color renderer r g b #xFF)
      (sdl2:render-draw-line renderer
                             column (round (+ (/ +window-height+ 2) (/ height 2)))
                             column (round (- (/ +window-height+ 2) (/ height 2)))))))

;; (run (list (make-instance 'fp-view :location (make-instance 'location))))

Monday, September 5, 2022

Drawing a circle

SDL is your bare bones graphics interface. It gives you primitives like draw-point, draw-line, and draw-rectangle, but you're on your own if you want to draw a circle. Naturally, I cribbed the code from stackoverflow.

Small circles didn't look right — they were a little squarish. The code worked by walking along pixels in one direction and accumulating an error term. When the error term got large enough, it would be reset and the code would advance a step along the other pixel axis. The error term was computed using integer math so that the circle was drawn quickly. The problem is that the integer math has rounding error and the rounding error is noticable with small circles.

For reference, it is straightforward to draw an exact circle. A circle isn't a function, but circle segment between 0 and 45 degrees is a function. If we mirror that segment eight ways horizontally, vertically, and at 90 degrees, we'll get a full circle.

(defun eightfold-point (renderer center-x center-y x y)
  (sdl2:render-draw-point renderer (+ center-x x) (+ center-y y))
  (sdl2:render-draw-point renderer (+ center-x x) (- center-y y))
  (sdl2:render-draw-point renderer (- center-x x) (+ center-y y))
  (sdl2:render-draw-point renderer (- center-x x) (- center-y y))

  (sdl2:render-draw-point renderer (+ center-x y) (- center-y x))
  (sdl2:render-draw-point renderer (+ center-x y) (+ center-y x))
  (sdl2:render-draw-point renderer (- center-x y) (- center-y x))
  (sdl2:render-draw-point renderer (- center-x y) (+ center-y x)))

(defun draw-circle (renderer center-x center-y radius)
  (let ((r-squared (* radius radius)))
    (dotimes (x (1+ (ceiling (/ radius (sqrt 2)))))
      (eightfold-point renderer
                       center-x center-y
                       x (round (sqrt (- r-squared (* x x))))))))
This gives much rounder looking circles than the code I cribbed from stackoverflow.

The problem, of course, is that this code computes a square root on each iteration. These days, computers are fast and that's not a big issue, but let's try to improve things.

On each iteration, we are computing the square root of r2-x2. This quantity changes between iterations, but not by much. You can compute a square root pretty quickly using Newton's method. The square root computed last iteration isn't that far from the new square root, so we can use the previous square root as the initial guess for Newton's method. Since we started pretty close to the right answer, we only need a single pass of Newton's method to get close enough to the square root for the current iteration.

(defun average (a b) (/ (+ a b) 2))

(defun draw-circle1 (renderer center-x center-y radius)
  (let ((x-limit (ceiling (/ radius (sqrt 2)))))
    (do ((x 0 (1+ x))
         (o 1 (+ o 2))
         (r2-x2 (* radius radius) (- r2-x2 o))
         (y radius (round (average y (/ (- r2-x2 o) y)))))
        ((> x x-limit))
      (eightfold-point renderer center-x center-y x y))))

This gives us pretty round circles — rounder than the integer methods, but not as round as the sqrt method. It requires less arithmetic than the sqrt method, but more than the integer method. What is more annoying, squarish circles or the amount of time it takes to draw round ones?


Wednesday, August 17, 2022

Playing with graphics

I wanted to play with some graphics. I don't know much about graphics, so I wanted to start with the basics. I played around with a couple of demos and I found that easiest to get reliably working was SDL2.

After downloading the SDL binary library and installing the FFI bindings with Quicklisp, I was off and running. You can find numerous SDL demos and tutorials on line and I tried a number of them. After I felt confident I decided to try something simple.

One thing I've noticed about graphics programs is the ubiquity of mutable state. Everything seems mutable and is freely modified and global variables abound. As a mostly functional programmer, I am alarmed by this. I wanted to see where we'd get if we tried to be more functional in our approach and avoid mutable data structures where practical.

Now the pixels on the screen had best be mutable, and I'm not trying to put a functional abstraction over the drawing primitives. We'll encapsulate the rest of the state in a state machine that is driven by the SDL event loop. The state machine will keep track of time and the current world. The current world is simply an immutable list of immutable objects. The state machine can transition through a render! phase, where it renders all the objects in the current world to a fresh frame. It attempts to do this about 75 times a second. The state machine can also transition through a next-world phase, where the current world and a delta-t are used to compute a new version of the world.

Our run program will take the initial list of objects. We'll start by initializing SDL, creating a window, and allocating a renderer for that window:

(defun run (initial-world)
  (sdl2:with-init (:video)
    (sdl2:with-window (window
                       :h +window-height+
                       :w +window-width+
                       :flags '(:shown))
      (sdl2:with-renderer (renderer window :index -1 :flags '())
        ... )))

Now we need the event loop state. last-ticks records the value from sdl2:get-ticks from the last time we processed the :idle event. This will be used to compute the elapsed time in ticks. render-ticker will record how many ticks have elapsed since the last time we rendered a frame to the screen. When render-ticker exceeds a certain amount, we'll call (render! current-world) and reset the ticker to zero. title-ticker will record how many ticks have occurred since the last time the window title was updated. When title-ticker exceeds a certain amount, we'll call sdl2:set-window-title to update the window title with some stats. sim-count is simply the number of times we've iterated next-world and frame-count is the number of times we've called render!. These are reset to zero every time we refresh the window title, so we'll have the frames per second and the world steps per second in the window title.

        (let ((last-ticks 0)
              (render-ticker 0)
              (title-ticker 0)
              (sim-count 0)
              (frame-count 0)
              (world initial-world))

          (flet ((title-tick! (dticks)
                   (incf title-ticker dticks)
                   (when (>= title-ticker 1000)
                     (decf title-ticker 1000)
                     (sdl2:set-window-title window (format nil "Sim rate: ~d, Frame rate: ~d" sim-count frame-count))
                     (setq sim-count 0)
                     (setq frame-count 0)))

                 (world-tick! (dticks)
                   (incf sim-count)
                   (setq world (next-world world (/ dticks 1000))))

                 (render-tick! (dticks)
                   (incf render-ticker dticks)
                   (when (>= render-ticker 13)
                     (incf frame-count)
                     (decf render-ticker 13)
                     (render-world! renderer world))))

Now we can run the event loop. The idle event is where the action happens:

          (sdl2:with-event-loop (:method :poll)

              (:idle ()
                     (let ((this-ticks (sdl2:get-ticks)))
                       (if (= this-ticks last-ticks)
                           (sdl2:delay 1)
                           (let ((dticks (- this-ticks last-ticks)))
                             (setq last-ticks this-ticks)
                             (title-tick! dticks)
                             (world-tick! dticks)
                             (render-tick! dticks)))))

              (:keydown (:keysym keysym)
                        (case (sdl2:scancode keysym)
                          (:scancode-escape (sdl2:push-quit-event))
                          (:scancode-x      (sdl2:push-quit-event))))

              (:quit () t))

Now that's a bunch of state, but it's more or less under control because what we have is a state machine and the state variables aren't accessible to anything.

render-world! is straightforward. It clears the renderer, calls render! on every object in the world, and presents the renderer for display.

(defun render-world! (renderer world)
  (sdl2:set-render-draw-color renderer #x00 #x00 #x00 #xFF)
  (sdl2:render-clear renderer)
  (mapc (lambda (object) (render! renderer object)) world)
  (sdl2:render-present renderer)
  )

next-world is a function that maps the current world to the next. It basically calls next on each object in the world and accumulate the results. We want objects to be able to go away, so if (next object) returns nil, we don't accumulate anything in the new world. If next returns the object unchanged, it will be accumulated unchanged in the next world. (next object) returns a new version of an object to simulate an update to the object. We want to be able to increase the amount of objects, so we allow (next object) to return a list of objects to be accumulated.

(defun next-world (previous-world dt)
  (fold-left
   (lambda (items item)
     (let ((more (next item dt)))
       (cond ((null more) items)
             ((consp more) (append more items))
             (t (cons more items)))))
   '()
   previous-world))

We'll start with a user-controlled player.

(defclass player ()
  ((x :initarg :x
      :reader get-x)
   (y :initarg :y
      :reader get-y)))

Everything that is to be displayed needs a render! method. This one just draws a little green triangle facing up.

(defmethod render! (renderer (player player))
  (let ((x (floor (get-x player)))
        (y (floor (get-y player))))
    (sdl2:set-render-draw-color renderer #x00 #xFF #x00 #xFF)
    (sdl2:render-draw-line renderer (- x 8) (+ y 8) (+ x 8) (+ y 8))
    (sdl2:render-draw-line renderer (- x 8) (+ y 8) (- x 1) (- y 16))
    (sdl2:render-draw-line renderer (+ x 8) (+ y 8) x (- y 16))
    (sdl2:render-draw-point renderer x y)
    ))

The next method computes the player in the next world:


(defun x-input ()
  (- (if (sdl2:keyboard-state-p :scancode-right)
         1
         0)
     (if (sdl2:keyboard-state-p :scancode-left)
         1
         0)))

(defun y-input ()
  (- (if (sdl2:keyboard-state-p :scancode-down)
         1
         0)
     (if (sdl2:keyboard-state-p :scancode-up)
         1
         0)))

(defparameter +player-speed+ 200.0) ;; pixels per second

(defmethod next ((player player) dt)
  (let ((new-x (max 8  (min (- +window-width+ 8)
                            (+ (get-x player)
                               (* (x-input) +player-speed+ dt)))))
        (new-y (max 16 (min (- +window-height+ 8)
                            (+ (get-y player)
                               (* (y-input) +player-speed+ dt))))))
    (make-instance 'player :x new-x :y new-y)))

Once we've defined a render! method and a next method, we're ready to go. If we call run on a list containing a player object, we'll have our little player on the screen controllable with the arrow keys.

An enemy ship can be defined.

(defclass enemy ()
  ((x :initarg :x :reader get-x)
   (y :initarg :y :reader get-y)
   (dx :initarg :dx :reader get-dx)
   (dy :initarg :dy :reader get-dy)))

(defmethod next ((enemy enemy) dt)
  (let ((new-x (+ (get-x enemy) (* (get-dx enemy) dt)))
        (new-y (+ (get-y enemy) (* (get-dy enemy) dt))))
    (when (and (>= new-x 8)
               (< new-x (+ +window-width+ 8))
               (>= new-y 8)
               (< new-y (- +window-height+ 16)))
      (make-instance 'enemy
                     :x new-x
                     :y new-y
                     :dx (get-dx enemy)
                     :dy (get-dy enemy)))))

;;; Render method omitted

As given, enemy ships will drift at constant speed until they run off the screen. We'd like to replenish the supply, so we'll make an enemy spawner:

(defclass enemy-spawner ()
  ((timer :initarg :timer :initform 0 :reader get-timer)))

(defmethod next ((spawner enemy-spawner) dt)
  (let ((new-time (- (get-timer spawner) dt)))
    (if (> new-time 0)
        (make-instance 'enemy-spawner :timer new-time)
        (list (make-instance 'enemy-spawner :timer (+ 1 (random 4)))
              (make-instance 'enemy :x (random (+ (- +window-width+ 32) 16))
                                    :y 16
                                    :dx (- 25 (random 50))
                                    :dy (+ (random 100) 50))))))

(defmethod render! (renderer (enemy-spawner enemy-spawner))
  nil)
The render! method doesn't do anything so a spawner doesn't have an image. It simply has a timer. To compute the next spawner, we subtract dt and create a new spawner with the reduced amount of time. If that's not a positive amount of time, though, we create two objects: a new spawner with somewhere between 1 and 5 seconds time and a new enemy ship.

We'll modify our player to allow him to shoot at the enemy:

(defclass player ()
  ((x :initarg :x
      :reader get-x)
   (y :initarg :y
      :reader get-y)
   (fire-cycle :initarg :fire-cycle
               :initform 0
               :reader get-fire-cycle)))

(defmethod next ((player player) dt)
  (let ((new-x (limit 8 (- +window-width+ 8) (+ (get-x player) (* (x-input) +player-speed+ dt))))
        (new-y (limit 16 (- +window-height+ 8) (+ (get-y player) (* (y-input) +player-speed+ dt))))
        (next-fire-cycle (- (get-fire-cycle player) dt)))
    (if (and (sdl2:keyboard-state-p :scancode-space)
             (< next-fire-cycle 0))
        (list
         (make-instance 'player
                        :x new-x
                        :y new-y
                        :fire-cycle .1)
         (make-instance 'bullet
                               :x (- (get-x player) 8)
                               :y (- (get-y player) 16)
                               :dx 0
                               :dy (- +bullet-speed+))
         (make-instance 'bullet
                               :x (+ (get-x player) 8)
                               :y (- (get-y player) 16)
                               :dx 0
                               :dy (- +bullet-speed+)))
        (make-instance 'player
                       :x new-x
                       :y new-y
                       :fire-cycle next-fire-cycle))))

A bullet is a simple moving object:

(defclass bullet ()
  ((x :initarg :x :reader get-x)
   (y :initarg :y :reader get-y)
   (dx :initarg :dx :reader get-dx)
   (dy :initarg :dy :reader get-dy)))

(defmethod next ((bullet bullet) dt)
  (let ((new-x (+ (get-x bullet) (* (get-dx bullet) dt)))
        (new-y (+ (get-y bullet) (* (get-dy bullet) dt))))
    (when (and (>= new-x 0)
               (< new-x +window-width+)
               (>= new-y 0)
               (< new-y +window-height+))
      (make-instance 'bullet
                     :x new-x
                     :y new-y
                     :dx (get-dx bullet)
                     :dy (get-dy bullet)))))

At this point we can move around the screen and shoot at enemies that spawn periodically. The problem is that the bullets go right through the enemy. We need to handle object collisions. We'll modify the next-world function. As it loops over the objects in the world, it will perform an inner loop that checks for collisions with other objects. If two objects collide, a function is called to get the collision results and those results are added to the list of objects in the world. If an object doesn't collide with anything, the next method is called to get the next version of the object.

(defun next-world (previous-world dt)
  (let outer ((tail previous-world)
              (next-world '()))
    (cond ((consp tail)
           (let ((this (car tail)))
             (let inner ((those (cdr tail)))
               (cond ((consp those)
                      (let ((that (car those))
                            (others (cdr those)))
                        (if (collides? this that)
                            (outer (append (collide this that) (delete that (cdr tail)))
                                   next-world)
                            (inner others))))
                     ((null those)
                      (outer (cdr tail)
                             (let ((more (next this dt)))
                               (cond ((consp more) (append more next-world))
                                     ((null more) next-world)
                                     (t (cons more next-world))))))
                     (t (error "Bad list."))))))
          ((null tail) next-world)
          (t (error "Bad list.")))))

We define collides? as a generic function that returns nil by default

(defgeneric collides? (this that)
  (:method ((this t) (that t)) nil)
  )
so that most objects don't collide. In the case where something does collide, we'll define collide as a generic function that returns nil by default
(defgeneric collide (this that)
  (:method ((this t) (that t)) nil)
  )
so when two objects collide, they simply disappear.

collides? will be called on pairs of objects in no particular order, so method pairs will be needed to handle both orders. We'll define collides? methods on bullets and enemies that checks if the bullet is within the bounding box of the enemy:

(defmethod collides? ((this bullet) (that enemy))
  (and (> (get-x this) (- (get-x that) 8))
       (< (get-x this) (+ (get-x that) 8))
       (> (get-y this) (- (get-y that) 8))
       (< (get-y this) (+ (get-y that) 8))))

(defmethod collides? ((this enemy) (that bullet))
  (collides? that this))

At this point, we can shoot enemy ships. The default method for collide between an enemy and a bullet returns nil so the enemy and the bullet simply disappear. If we were fancy, we could arrange for it to return an explosion object or several debris objects.

It would be nice to keep a tally of the number of enemy ships we have shot. We don't have to add any extra machinery for this. We create a score class and a point class:

(defclass score ()
  ((value :initarg :value
          :reader get-value)))

(defmethod next ((score score) dt) score)

;;; Render method prints score on screen.

(defclass point ()
  ((value :initarg :value
          :reader get-value)))

(defmethod next ((point point) dt) point)

(defmethod render! (renderer (point point)) nil)
Scores and points are immutable objects without positions, but we'll define methods so that when a score and a point collide, the result is a higher score.
(defmethod collides? ((this point) (that score)) t)
(defmethod collides? ((this score) (that point)) t)

(defmethod collide ((this point) (that score))
  (list (make-instance 'score
                       :font (get-font that)
                       :value (+ (get-value this) (get-value that)))))

(defmethod collide ((this score) (that point))
  (collide that this))
Now we'll define a bullet colliding with an enemy to produce a point:
(defmethod collide ((this bullet) (that enemy))
  (list (make-instance 'point :value 1)
        ;; add explosion object here
        ))

(defmethod collide ((this enemy) (that bullet))
  (collide that this))
So when you shoot an enemy, the bullet and enemy disappear to be replaced by a point. On the next update, the point will collide with the score to be replaced with an updated score.

At this point we've got a little demo game where we can fly a ship around and shoot enemies and a running score is kept. The world model is immutable and worlds are functions of previous worlds. I'll call it a successful proof of concept.

But did this buy us anything? We don't have mutable state per se, but we've kind of cheated. When we create new versions of an object, each version is immutable, but the sequence of versions taken as a whole seem to be evolving over time. For example, consider the score. At each time step, there is an immutable score object, but over time what is considered the current score changes. We've eliminated the direct problems of mutation, but we've introduced the problem of keeping track of what series of immutable versions correspond to a single evolving instance.

In this small example, we're not keeping track of the evolving objects. For instance, each bullet, as it is updated from step to step, is actually created anew at its new position on each step. The old bullet instance is dropped and the new instance really has no idea how it got there. Bullets are such simple objects that this doesn't matter, but the current score is different. It makes sense for there to be a singleton score object that increases over time, but we haven't built that in to our model. Instead, we've designed a set of collision interactions that drive the score.

We've eliminated the direct mutable state in our objects and our world, but sometimes we want to model stateful objects. We therefore create objects that represent state transitions (e.g. points) and then use the collision mechanism to combine the transition objects with the objects that represent physical entities. That seems a bit convoluted, and I don't think it will scale.

On the other hand, we do gain the immediate benefits of the world and the objects being immutable. Saving and restoring a world is trivial. Reasoning about objects at each update is easy because the objects don't change, but we now have to reason about how objects appear to change in the long run.

The tradeoff of immutable objects is increased allocation. But although a lot more consing is happening, most of it can be quickly reclaimed by the generational collector, so noticable GC pauses are infrequent. I haven't measured the performance, but it is certainly adequate for the little example I wrote. If you had a mechanism to reason about the objects as linear types (Sufficiently Smart Compiler™), you could determine when you can update objects in place and avoid reallocating.

The world model, simply a list of objects, is flexible, but not well structured. For instance, the current score and the bullets are among the elements mixed together in this list. You'd have to search this list to find specific elements or filter this list to find elements of a certain type.

The simple collision model is O(n2), so it can't handle a ton of objects. A more sophisticated world model would be needed to keep track of different classes of collidable objects to avoid the O(n2) search. For example, if bullets were kept separately, we could avoid checking if they collide with each other.

The point of this exercise was to play with graphics and see what you can do without the mutable state that is alarmingly ubiquitous. It turns out that you can go pretty far, but it's kind of strange.