Tuesday 11 February 2014

Parallel Delaunay Triangulation

Resource handling post will happen in a bit, as I got distracted after reading a fascinating paper entitled Localizing the Delaunay Triangulation and its Parallel Implementation (Renjie Chen/Craig Gotsman). The paper itself is quite googleable, although in some cases it hides behind a paywall.
After a few false starts, it translated to Haskell relatively cleanly. As ever, this is not terribly idiomatic but I thought it might be interesting to some people.

First, some elided code: I have some very, very simple linear algebra types kicking around of the form FloatN. They're designed to mimic the ones found in shading languages or in SIMD intrinsics rather than behave correctly, so you can multiply them together and add a scalar to them if you like.

Further, I cobbled together a quick k-d tree implementation for doing spatial queries. In this case, the only important one is fetching all the objects within a circle. For my purposes I'm only interested in periodic triangulations and corresponding Voronoi cells, and will be restricting the point distribution to be dense with a known minimum distance between points. Anticipating a move to a more efficient spatial structure at a later date, here's the interface I'm broadly coding against:

class SpatialStructure2D a where
  type Query a
  positionOf :: a -> Query a -> Float2
  queryCircle :: a -> Float2 -> Float -> [Query a]

In particular, the Query a type and positionOf function exist so some future structure can use indices into a flat array or vector of positional data, which will be useful when constructing vertex buffers later. Might be thinking too far ahead there.

On to the implementation! A quick record type holds the local Delaunay one-ring (vertices directly connected to the point under consideration in the triangulaiton) and the circumcircles described by this fan of triangles. The circumcenters obviously form our nice little Voronoi cell at the end.

data Delaunay = Delaunay {
  dlyPoint :: !Float2,
  dlyRing :: [Float2],
  vorCell :: [(Float2,Float)] } deriving (Eq, Show)

initDly :: Float2 -> Delaunay
initDly c = Delaunay {
  dlyPoint = c,
  dlyRing = p0,
  vorCell = zipWith (circumcircle c) p0 (rotateList p0)
  } where
    p0 = [float2 1.5 0.5, float2 0.5 1.5, float2 (-0.5) 0.5, float2 0.5 (-0.5)]

The utility function initDly just makes a massive quad around the test point so we can start from a known (but uselessly huge) triangulation. Some more tangential stuff:

circumcircle :: Float2 -> Float2 -> Float2 -> (Float2,Float)
circumcircle a b c = (cen,rad)
    cen = circumcenter a b c
    rad = circumradius a b c

rotateList :: [a] -> [a]
rotateList [] = []
rotateList (x:xs) = xs ++ [x]

rotateUntil :: (a -> Bool) -> [a] -> [a]
rotateUntil _ [] = []
rotateUntil _ [a] = [a]
rotateUntil p xs = go (take (length xs) (iterate rotateList xs))
    go [] = xs
    go (ys@(y:_) : yss)
      | p y = ys
      | otherwise = go yss

The functions to calculate the circumcenter and circumradius of a triangle are omitted because I shamelessly translated the implementations from here, and they're quite long and uninteresting. The list rotations are ugly but made things easy.

data Halfplane = Halfplane Float2 Float2 Float
halfplane :: Float2 -> Float2 -> Halfplane
halfplane v o = Halfplane o (normalise d) (0.5 * vecLen d) where d = v - o

halfplaneTest :: Halfplane -> Float2 -> Float
halfplaneTest (Halfplane o n d) x = d - (n `dot` (x - o))

inHalfplane :: Halfplane -> Float2 -> Bool
inHalfplane h x = halfplaneTest h x >= 0

Testing whether or not points are inside the halfplane defined by the bisector of two points is key to the incremental bit of the algorithm, where we test a point against the current local triangulation:

updateDT :: Delaunay -> Float2 -> Delaunay
updateDT locDT v = if inHvQ == q0 then locDT else locDT { dlyRing = p1, vorCell = q1 }
    c = dlyPoint locDT
    p0 = dlyRing locDT
    q0 = vorCell locDT
    hv = halfplane v (dlyPoint locDT)
    passHv = inHalfplane hv . fst . fst
    failHv = not . passHv
    --ensure the prefix of the one-ring/voronoi polygon is in halfplane
    qp0 = rotateUntil passHv (rotateUntil failHv (zip q0 p0))
    (inHv, outHv) = span passHv qp0
    (inHvQ, inHvP) = unzip inHv
    --avoid strictness, profit
    dp0 = snd (head outHv)
    dp1 = head inHvP
    newQ0 = circumcircle c dp0 v
    newQ1 = circumcircle c v dp1
    q1 = newQ0 : newQ1 : inHvQ
    p1 = dp0 : v : inHvP

Again, ugly. Apart from the noise of packing and unpacking tuples, this code is essentially rotating both the one-ring and Voronoi cell until the circumcenters outside the halfplane under consideration comprise the suffix of the list. If no points fall outside the halfplane, the triangulation remains unchanged and we just return it. Otherwise we introduce the new point in place of those we removed and stitch everything together. Note that although we furiously and pointlessly rotate both the one-ring and Voronoi polygon, they never change winding order.

Finally, to actually calculate the local triangulation from scratch:

localDT :: (SpatialStructure2D a, Query a ~ Float2) => a -> Float -> Float2 -> Delaunay
localDT tree r c = updateFrom dly0
    dly0 = foldl' updateDT (initDly c) (filter (/= c) v0)
    v0 = queryCircle tree c r
    updateFrom dly = if null vs then dly else updateFrom (updateDT dly (head $ head vs))
        fl = filter (`notElem` (c: (dlyRing dly)))
        vss :: [[Float2]]
        vss = map (fl . uncurry (queryCircle tree)) (vorCell dly)
        vs = dropWhile null vss

Technically I should rewrite this to use the positionOf and relax the restriction on Query a. Oops. Anyway!
Because we know the minimum distance between points and that the point set is relatively dense, passing in an initial query radius to quickly cut the initial working set down to size is a huge win. Then all that is left is to try each circumcircle in the current triangulation - if it might add a new point, update the local triangulation and recurse. When every circumcircle is empty (excepting those points already in the triangulation), we're done.

Now all that we need to do is something like map (localDT tree (r*2.0)) points `using` parListChunk 32 rseq, where tree is a spatial structure containing our random points, r is the minimum distance between those random points... and we get parallel computation of all the interesting things.

Now, there are lots and lots of bits of code that are ugly and somewhat hacky in here, and it could certainly be optimised hugely (profiling suggests most of the time is spent traversing my k-d tree), this still manages to spit the periodic Voronoi diagram for ~64k points to a .obj file on disc in around ten seconds (on an 8-core i7). Should absolutely be sufficient for my purpose as it stands.

And courtesy of Blender and the subdivision modifier, here are those little cells with a bit of smoothing applied:

Tuesday 4 February 2014

Rendering With Haskell

I've been on a bit of a Haskell kick lately.

The reason? I mostly blame Light Table - I've been testing the alpha version with various Clojure toy projects and enjoying myself mightily, when I somewhat randomly opened a .hs file and discovered not only syntax highlighting but also relatively clean, not-overly-smart indentation. Haskell indentation is (IMO) much nicer than that of Python in that the layout is more permissive and you can always fall back on braces and semicolons in a pinch, but its nice when the editor doesn't try to be too clever (and shoots you in the foot with its cleverness). I hate feeling like I'm fighting a text editing tool, and have nearly come to blows with Visual Assist in the course of the day job.

So, I spent a bunch of time wrapping up the basics of the OpenGLRaw package - HOpenGL sadly lacks support for matrix uniforms where shaders are concerned, and exposes/uses the Compatibility profile. I'm trying to restrict myself to the relatively small surface area of OpenGL 3.1 - 3.2 Core. GLFW-b provides window/input management; JuicyPixels seems reasonable for loading various shapes of image file and doesn't involve linking against DevIL. Some really, really ugly Parsec code for dealing with .obj files exported by Blender completes the foundations for a cube-with-moveable-camera prototype.

Now the fiddly bits. Resource management, tying the game's logic to the rendering thread, input handler stacks, GUI, all that fun stuff. Trying to keep as much of it as possible in idiomatic, functional code and outside the IO monad. Much of it would be simplified by sticking to a single thread, but there's no point having all these lovely fast cores lying around if we don't use them.

Resource management is the first priority, because having a bunch of explicit paths lying about your main source file along with various attendant bits of code to load and render textures/meshes/shaders is... horrible. At this point the topic starts to veer into the venerable Wavefront .obj format and it's auxilliary material definitions, which thankfully form a very simple chain of dependencies for what files to grab and turn into OpenGL resources before you can render something. I'll go into that next time I remember to blog, I think.

Friday 18 January 2013

Sweden, Ho!

I moved again! Now living in Uppsala, working in Stockholm. The past few months have been pretty exciting and we just moved to a new flat, but things are settling down. I've had a little time to poke at my personal projects, although I'd always like more.

My recent obsession has been with managing state. For many things, I seem to be tending towards a pipeline approach:

(defn foo [state x y z]
  (let [{:keys [a b c]} state]
    (assoc state 
           :a (bar x a)
           :b (baz b z y)
           :c (inc c))))

(defn quux [state]
   (-> state
       (foo "x" [] 0)
       (foo "y" [:r :s] 42)))

(defn outer-loop [state] (recur (quux state))

And so on and so forth. Passing state around explicitly as maps, writing functions that take the current state as their first argument and return that state modified with assoc, dissoc, assoc-in, update-in, merge etc.

Honestly, I'm not sure that this approach is the best way to do what I want. I have spent a little time looking at Clojure's state monad, but I suspect this is overkill (and despite being simple in practice, the whys and wherefores of all monads ever continue to befuddle me). The piped-stateblob approach translates nicely into ClojureScript with WebGL too; this is another thing I have been investigating of late.

Other advantages include very quick inspection or serialisation of a system's state (pprint to the rescue), and a single obvious point to recurse around in a thread's main loop. The stateblobs can also be useful carriers for communicating between threads.

On the other hand there's nothing particularly elegant or nice about carrying around the entire world with you wherever you go. And it works well enough that I'm starting to feel it is my hammer with which to treat every problem as a nail.

Using stateblobs everwhere: an ugly turn-based for prototype graphical RL #734.

Wednesday 14 March 2012

Clojure-flavoured Rendering Framework: State

My rendering framework, such as it is, sits on clojars.org and does... not very much. That should change. I can't guess how useful it'd be to anyone else, but it could certainly be a lot more friendly.

With that in mind, I started taking a critical look at a few emergent antipatterns in the very core of the framework: the OpenGL wrapper.

This example, I am happy with. Mostly.

Last time, I mentioned how the rendering functions have state very explicitly threaded through them in the form of an argument. A bunch of my wrapper functions do very simple things to mutate GL state. Some of them reset it afterwards, or wrap those that don't in a friendly way. For example, matrix stuff is done something like this:

  (ogl/translate 1 0 0)
  (ogl/scale 5.0)
    (ogl/vertex 0 0 0)
    (ogl/vertex 1 0 0)
    (ogl/vertex 1 1 0)))

Hopefully that's fairly obvious to someone familiar with OpenGL. The with-matrix wraps a glPush/PopMatrix pair of calls around all the expresions within it. Likewise triangles calls glBegin/glEnd with the appropriate primitive type. I debated not including this ancient method of drawing stuff, but it's so handy for debugging and quick tests... anyway, the point is that even though translate and scale mutate some OpenGL state, they're intended to be used within this kind of protected block. Maybe they should all have exclaimation marks, I'm a little bit hazy of when it's polite to start exclaiming when you're not mutating your arguments.

Obviously the triangles call and all the stuff therein is side effecting like crazy, but I'm not sure how much value attaches to trying to make those calls look different. If anyone has an opinion on the aesthetics of such things, please yell.

A lot of the rendering stuff behaves similarly to with-matrix. I have equivalent calls to wrap GL attribute fiddling, wireframe rendering (mostly debug friendliness), capturing the results of rendering expressions to render targets, assigning vertex buffers, and binding shaders. Some of these are could do with a refactor to make them more consistent, but they generally behave.

But then there's this...

So that's all well and good. Then I came across the with-matrix-mode macro.


This one looks so similar to the with-matrix example, but instead of pushing and popping from the current matrix stack to allow safe manipulation of a transformation, it changes which matrix stack is currently active. OpenGL has a few of these, but the only ones that are really relevant are your "world" transform (GL_MODELVIEW) and your view/projection transform (GL_PROJECTION). The former moves vertices around, the latter manipulates your notional camera, as it were.

The reason this makes me twitch is that this doesn't in any way revert the changes you make within it. For example, here's a little function that sets up an orthographic view in place of the conventional camera-like perspective view:

(defn orthographic
  "Loads an orthographic projection matrix to the top
of the projection matrix stack."
  ([left right bottom top near far]
     (with-matrix-mode :projection
       (GL11/glOrtho left right bottom top near far)))
  ([left right bottom top]
     (orthographic left right bottom top -1.0 1.0)))

Sure, at the end of the function the matrix mode has been restored to whatever it was at the start, but there's no attempt made to preserve the state of the projection matrix. That's rather the point, in fact - it gets murdered and replaced. Everything after this will be rendered using the orthographic projection.

I suppose this function is named in a misleading fashion. Anyone assuming the state of the GL will be unchanged after executing it will be disappointed, whereas elsewhere I'm trying to leave things as I found them. On the other hand, changing matrix mode is rather rare compared to rendering stuff, you generally do it once to set up your camera and then leave it. This kind of mutate-and-forget approach makes it relatively easy to split up the rest of your rendering into small chunks, as long as you keep things in the correct order.

Elsewhere similar compromises might sensibly be made with the excuse of efficiency; there's no sense constantly setting and then resetting a bit of state if the next expression simply sets it again, or worse, if it sets it back to the same value. All told, I probably just want to find a consistent way to signal that some functions do not behave with respect to render state, whilst the default remains that all functions clean up their toys after playing with them.

Wednesday 7 March 2012


Minor tweaks of late.

Restructured the rendering quite a bit. Keeping things very, very simple still: the rendering thread owns a list of functions that execute each frame, in the order in which they were added. Dumb as a pile of bricks, and assumes that each such render-function will actually do a bunch of stuff including culling.

Originally the rendering functions were all nullary and their results discarded. I swapped this over to each taking a single state argument, a hash map of arbitrary rendering guff, and each returns a (potentially modified) state blob. So the core of the rendering loop is a nice, cuddly-ish*
(reduce (fn [s f] (f s)) initial-state render-funcs).

Right now the uses for this blob are limited. I'm just bunging a few things which were previously updated from different threads in there, such as camera data (ad-hoc derefing of atoms and such gives inconsistent rendering and general ickyness). So it's more like a way to carry around an immutable copy of the game state as of the start of the frame than anything truy useful. Suspect it'll be handy later to mirror the GL state and avoid heavy state-changing operations where they're not necessary, do pre- and post-invariant checks and other contract-y things, accumulate frame stats such as objects rendered... stuff.

Image is just some randomly coloured lights scattered around, this time with a bit less gloom. Still quite gloomy. I find myself wondering if it might be better to work with entirely solid hexagonal prisms as a starting point, rather than more elaborate geometry. Hmm.

* Cuddly-ish, but not very cuddly. I would really rather use something like (flip apply) instead of the embedded anonymous function, but cannot quite seem to conjure the right way to do this without doing something slightly horrible. Will look at it later with fresh eyes and doubtless slap my forehead at my own stupidity.

Monday 27 February 2012


Decided it was time to fling some deferred rendering at the project, because who doesn't like taking huge chunks of VRAM to do simple things? Thanks to Jotaf's comments last week, I have a much easier method for making and rendering walls, so it only took a few minutes of Blender-torture to generate new ones.

The question of style is one I haven't resolved yet, alas. Making tech is far, far easier than trying to hang a scene together as a coherent whole. I know the cultural roots the game world will draw from, but not whether to go for a cell-shaded or painterly or somewhat realistic approach. The latter is probably most tedious in terms of rendering tech and modelling, but is almost certainly the easier to at least rough out textures based on photographic sources.

Sunday 19 February 2012

2D Is Easier

Annoyed with my inability to find a nice way to create 3D tiles, I made a little tiny tile-editor for the 2D version.

Everything's so much easier in 2D. Ah well.


I added some water and grass tiles, but not enough for a new post. Trying to do something in the spirit of pixel art but that is amenable to rotation around 60 degree increments (and requiring less in the way of raw talent). I think I might need to try flat colour polygons and a finer subidivision, but it's an interesting clean look.