Log-Gap Charts for Time Series

Sometimes you want to visualize your time-series data in a chart. Traditionally, this is done by slicing your time-scale into equal-sized bins and counting the events that fall into each one. The big players—Grafana, Datadog, CloudWatch—all do this, handing you a nice, neat histogram of events over time.

But binning the data destroys resolution.

Choose too large a bin, and you lose time resolution. You completely miss high-density burst events and long idle gaps; it all just averages out into a meaningless block. Choose too small a bin, and you lose aggregate resolution. You end up with hundreds of atomized bins containing zero, one, or two events each, flattening your chart into a needle-bed of noise. If you're lucky, you might find an intermediate bin size that still shows you something marginally interesting on both the count and time axes. But more often than not, there is no reasonable sweet spot. If you're particularly unlucky, a bin boundary will arbitrarily divide a burst right down the middle, and you'll miss the anomaly altogether.

Furthermore, binned data is highly sensitive to scale. If you zoom in on a bin, you don't get a closer look at the behavior; you just get a histogram with a single lonely bar in it. If you zoom out, you don't get a higher-elevation view of systemic trends; you simply add more bins.

Let's look at a live example. Here is a chart of the past 48 hours of HTTP requests hitting my web site, sliced into 3-hour bins:

We can see that some 3-hour blocks obviously have more requests than others, but the resolution is crude. We only have about fifteen or so blocks of data. It tells us almost nothing about how the traffic arrived. Was it steady across the three hour interval or did it come in bursts?

Here is the exact same 48 hours of live data, but rendered with 15-minute bins:

Now we can more clearly see the bursts, but for the steady state the y-axis is practically useless. Most of the bins have zero or one element in them. We get absolutely no feel for the macro, hourly rate.

Because we take the logarithm of the gap, the data become insensitive to scale. We can plot small intervals of a single second right next to large intervals of hours on the same axis without losing the shape of either. The steady state appears as a cloud of points high up in the chart, while burst traffic shows as vertical lines.

To read this kind of chart, you don't look at the individual points, but at the overall shape of the point cloud, the envelope, and the streaks. The highest points in the chart are the longest intervals between events, the lowest are the shortest. The middle of the point cloud indicates the median time between events. I have been using charts like this to plot time series events and they give a good feel of how events arrive.

Note: if you are getting broken links, that is probably the rate limiting on my site. Give it some time and reload later.

Motion without Integration

A game where things move around has equations of motion that describe how objects move. We need to solve these to produce the locations of the objects as a function of time, and we have to do it in lock step real time. Fortunately, we usually don’t need highly realistic models or highly accurate results. Game phyiscs, while inspired by real world physics, can be crude approximations.

If the velocity of an object changes over time, we need to integrate the velocity to determine the position. This is typically done with forward Euler integration: the current velocity is multiplied by a small Δt and this is added to the current position. In other words, we assume the object moves linearly over the amount of time represented by Δt. If Δt is small enough, the linear segments of motion will approximate the actual curve of motion.

In a game, we use the game loop, which runs every few milliseconds, as our source of time. We subtract the current time from the previous time the game loop was run to obtain our Δt. We update each object’s state using forward Euler integration over Δt.

Our game loop keeps track of the last time it was run and subtracts that from the current time. It runs entity-step! on each entity in the game, passing in dticks.

(let ((start-tick (sdl2:get-ticks)))
  (let ((last-tick start-tick))
    (loop
      (let* ((tick (- (sdl2:get-ticks) start-tick))
             (dticks (- tick last-tick)))
        (map nil (lambda (entity)
                   (entity-step! entity dticks))
                 entities)
        (setq last-tick tick)))))

entity-step! performs our forward Euler integration:

(defun entity-step! (entity dticks)
  (incf (position entity) (* (velocity entity) dticks)))

Assuming the game loop runs reasonably frequently, the position of the entity will be a reasonable approximation of the integral of the velocity. If the game loop gets delayed for some reason, e.g. garbage collection, the corresponding increase in the value of dticks will make up for it.

But sometimes we have the fortune of having a closed form solution to the time integral. If this is the case, we can compute the entity’s state directly and we don’t need to integrate the state changes on every iteration of the game loop. Let me give two examples.

In the platformer game, there are potions that the player can take to restore his health. A potion appears as a vial that floats above the ground. A “bobbing” effect is achieved by changing the height above the ground in a time varying manner. Now we could implement this by having the game loop modify the y position of the potion on each iteration, but there is a closed form solution. We can compute the y position at any time by adding a constant base y position to a factor of the sine of the current time. This is easily done by specializing the get-y method on potions:

;;; Make potions bob up and down.
(defmethod get-y ((object potion))
  (+ (call-next-method)
     (* 5 (sin (* 2 pi (/ (sdl2:get-ticks) 1000))))))

call-next-method invokes the next most specific method — in this case, the get-y method of the entity class, which returns the base y position. We add a sine wave based on the current time.

This is only thing we need to modify to make a potion bob up and down. Values dependent on the y position, such as the attackbox, will bob up and down with the potion because it uses the get-y method to determine what the potion is.

A second example is cannonball motion. Cannonballs travel linearly away from the cannon at a constant velocity. Rather than stepping the cannonball by adding its velocity to its position on each update, we specialize the get-x method on cannonballs

(defmethod get-x ((cannonball cannonball))
  (+ (call-next-method)
     (* (get-speed cannonball)
        (- (sdl2:get-ticks) (get-start-tick cannonball)))))

If we specialize both the get-x and get-y methods we can easily get objects to trace out any desired parametric curve, like a Lissajous figure.

Animation: Two Methods

Animated sprites are much cooler that static ones and really make your game come to life. They are pretty easy to implement, but there are a couple of ways to implement them, and I've seen several projects use the less optimal strategy.

The less optimal strategy is straightforward. We keep a time counter in the animated entity and increment it every time we update the entity. When the counter exceeds the amount of time for an animation frame, we advance to the next frame and reset the ticker.

(defclass entity ()
  ((animation-ticker :initform 0 :accessor animation-ticker)))

(defun entity-step! (entity dticks)
  (incf (animation-ticker entity) dticks)
  (when (> (animation-ticker entity) ticks-per-frame)
    (advance-to-next-frame! entity)
    (decf (animation-ticker entity) ticks-per-frame)))

Now this works, but the reason it is less optimal is that it involves maintaining time varying state in the entity. Recall that the game has two primary threads running: a render thread which runs 60 times a second and draws the game on the screen, and an game loop thread which runs maybe 200 times a second and updates the state of the entities in the game. The point of this is to separate and abstract the drawing of the game from the running of the mechanics of the game. By using this less optimal method, we involve the game loop thread in maintaining state that is used only by the render thead. To anthromorphize, the entity “knows” that it is animated and is taking an active role in keeping things moving.

A better way is to put the responsibility of animation solely within the render thread by having it select which animation frame to display when it draws the sprite. The sprite contains a vector of animation frames and a timestamp at which the animation was started. When redrawing the sprite, the elapsed time is computed by subtracting the start timestamp from the current time, and the frame to display is computed by dividing the elapsed time by the time per frame modulo the length of the frame vector.

(defclass animation ()
  ((start-tick :initform (get-ticks) :reader start-tick)
   (frames :initarg frames :reader frames)))

(defun display-animation! (animation)
   (let* ((elapsed-time (- (get-ticks) (start-tick animation)))
          (absolute-frame (floor elapsed-time ticks-per-frame))
          (frame (mod absolute-frame (length (frames animation)))))
     (display-frame! (svref (frames animation) frame))))

This way of doing it no longer needs time varying state, and we’ve abstracted the display of the animation from the handling of the entity’s state. The entity no longer has to “know” that it is animated or do anything about it. There no longer has to be communication between the game loop thread and the render thread on every animation frame. This could be important if the render thread is running on a different machine than the game loop thread. If the game loop thread lags (maybe because of garbage collection or entering the error handler), we still get smooth animation.

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.

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))))

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 r²-x². 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?