I had some time on a flight to Boston the other week, and I decided I wanted to mess around a bit more with Löve2d's pixel shaders. I had the idea of creating some sort of GPU accelerated game of life. I figured the "canvas" feature of Löve would be a good fit. I could create a shader that takes the current iteration as a texture, and the output of the shader would be the pixels that should be "alive" in the next iteration. Then all I would have to do is constantly take this output and feed it back into the shader. It turned out to be simpler than I thought.

Here's the shader.

extern float grid_width;
extern float grid_height;

int cellOccupied(Image cells, vec2 cell_coord)  {
    vec2 coord = mod(cell_coord, 1.0);
    int x = int(coord * grid_width);
    int y = int(coord * grid_height);
    if(Texel(cells, coord) == vec4(1.0, 1.0, 1.0, 1.0)) {
        return 1;
    } else {
        return 0;
    }
}

int neighborCount(Image cells, vec2 c) {
    int neighborCount = 0;

    //we need to index the texture by coordinates [0, 1]
    float x = c.x;
    float y = c.y;
    float xInc = 1.0 / float(grid_width);
    float yInc = 1.0 / float(grid_height);

    //checks right 3 cells
    neighborCount += cellOccupied(cells, vec2(c.x - xInc, c.y + yInc));
    neighborCount += cellOccupied(cells, vec2(c.x - xInc, c.y));
    neighborCount += cellOccupied(cells, vec2(c.x - xInc, c.y - yInc));

    //checks left 3 cells
    neighborCount += cellOccupied(cells, vec2(c.x + xInc, c.y + yInc));
    neighborCount += cellOccupied(cells, vec2(c.x + xInc, c.y));
    neighborCount += cellOccupied(cells, vec2(c.x + xInc, c.y - yInc));

    //checks middle 2 cells
    neighborCount += cellOccupied(cells, vec2(c.x, c.y + yInc));
    neighborCount += cellOccupied(cells, vec2(c.x, c.y - yInc));

    return neighborCount;
}

vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 screen_coords) {
    bool alive = cellOccupied(texture, texture_coords) == 1;
    int neighborCount = neighborCount(texture, texture_coords);
    bool willBeAlive = false;

    //  if alive, lives if 2 or 3 surround it
    //  if dead, becomes alive is exactly 3 neighbors
    if ((alive && (neighborCount == 2 || neighborCount == 3)) || (!alive && neighborCount == 3)) {
            willBeAlive = true;
    }

    if (willBeAlive) {
        return vec4(1.0, 1.0, 1.0, 1.0);
    } else {
        return vec4(0.0, 0.0, 0.0, 1.0);
    }
}

Pixels (cells) that are alive are colored white by definition, and dead cells are black. I wrote a cellOccupied method to return a boolean if the pixel was white or black. This is complicated only by the fact that the texture coordinates are normalized on [0, 1]. I also made the grid wrap around the edges of the buffer. For each pixel, the fragment shader uses the neighborCount method in combination with the pixel's current state to return white or black for the next iteration, in accordance with the rules of the game of life.

The Lua script is very simplistic as well. The love.load() method creates a canvas to render the next iteration into using the shader. It also creates an image buffer to store the current iteration. The shader itself is also initalized.

function love.load(a)
    --  store these for later use (must be updated on window resize)
    windowWidth = lg.getWidth()
    windowHeight = lg.getHeight()

    --  how large the data buffers are (we are storing data in an image)
    dataWidth = math.ceil(windowWidth / sample)
    dataHeight = math.ceil(windowHeight / sample)

    --  create a canvas that an image can be drawn on, and then read from
    calculationCanvas = lg.newCanvas(dataWidth, dataHeight)

    --  create the shader, send the grid width and height
    calculationShader = lg.newShader("calculationShader.glsl")
    calculationShader:send("grid_width", dataWidth)
    calculationShader:send("grid_height", dataHeight)

    --  create an image to store the current iteration
    currentIteration = li.newImageData(windowWidth / sample, windowHeight / sample)
    currentIterationImage = lg.newImage(currentIteration)
    currentIterationImage:setFilter("linear", "nearest")
end

The love.draw() method simply draws the current iteration image into the window.

function love.draw()
    lg.draw(currentIterationImage, 0, 0, 0, sample)
end

Most of the interesting stuff happens inside the iterate() method. The current iteration is drawn into the offscreen "calculation" canvas, then the image data is taken out of the offscreen canvsas and stored to be drawn on the next call to love.draw(). I had some issues here with memory issues, Löve didn't seem to garbage collect frequently enough and the script would crash soon after launching after consuming a large amount of memory. It's bad practice, but calling the garabge collection method manually seems to fix the issue.

function iterate()
    --  draw the current iteration into a canvas, and let the pixel shader
    --  calculate the next iteration of the grid
    lg.setCanvas(calculationCanvas)
    lg.setShader(calculationShader)
    lg.draw(currentIterationImage)
    lg.setShader()
    lg.setCanvas()

    --  now, get the next iteration as image data out of the calculationCanvas
    currentIteration = calculationCanvas:newImageData()
    currentIterationImage = lg.newImage(currentIteration)
    currentIterationImage:setFilter("linear", "nearest")
    collectgarbage()
end

The love.update() method lets the user make pixels in the current iteration alive or dead using the two mouse buttons. If the "g" key is pressed, the script spawns a glider at the mouse location. The function also handles calling iterate() at the right time, since the script also allows the user to control down the iteration speed. A simple method handles creating the gliders. Ugly, but it works.

function love.update(dt)
    if love.mouse.isDown(1, 2) then
        local x = love.mouse.getX()
        local y = love.mouse.getY()

        local imageX = math.floor(x / sample);
        local imageY = math.floor(y / sample);

        if love.mouse.isDown(1) then
            if love.keyboard.isDown('g') then
                makeGlider(imageX, imageY)
            else
                currentIteration:setPixel(imageX, imageY, 255, 255, 255, 255)
            end
        else
            currentIteration:setPixel(imageX, imageY, 0, 0, 0, 255)
        end

        currentIterationImage:refresh()
    end

    if running then
        time = time + dt
        if time > (1 / iterationsPerSecond) then
            iterate()
            time = 0
        end
    else
        time = 0
    end
end

function makeGlider(x, y)
    currentIteration:setPixel(x, y - 1, 255, 255, 255, 255)
    currentIteration:setPixel(x + 1, y, 255, 255, 255, 255)
    currentIteration:setPixel(x - 1, y + 1, 255, 255, 255, 255)
    currentIteration:setPixel(x, y + 1, 255, 255, 255, 255)
    currentIteration:setPixel(x + 1, y + 1, 255, 255, 255, 255)
end

The last bit of the user interface was some keyboard shortcuts. The spacebar pauses or resumes the simulation. The "s" key allows single-stepping the iterations. The "up" and "down" arrow keys let the user adjust the simulation rate. The "r" key fills the entire window with random pixels and "c" clears the current window.

The most fun easter egg I created lets the user hit the "h" key to fill the entire screen with a grid of gliders. Extra fun is letting the simulation run and then disturbing the grid of gliders and watching the disturbance travel throughout the grid.

function love.keypressed(key, scancode, isrepeat)
    if key == 'space' then
        running = not running
        if running then
            lw.setTitle("game of life - running")
        else
            lw.setTitle("game of life - paused")
        end
    end

    if key == 's' then
        iterate()
    end

    if key == 'up' then
        iterationsPerSecond = iterationsPerSecond + 5
    end
  
    if key == 'down' then
        iterationsPerSecond = iterationsPerSecond - 5
    end

    if key == 'r' then
        for x = 0, (dataWidth - 1) do
            for y = 0, (dataHeight - 1) do
                local alive = love.math.random(0, 4)
                if alive < 1 then
                    currentIteration:setPixel(x, y, 255, 255, 255, 255)
                else
                    currentIteration:setPixel(x, y, 0, 0, 0, 255)
                end
            end
        end
        currentIterationImage:refresh()
    end

    if key == 'c' then
        clear()
    end

    if key == 'h' then
        for x = 1, dataWidth - 4, 5 do
            for y = 1, dataHeight - 4, 5 do
                makeGlider(x, y)
            end
        end
        currentIterationImage:refresh()
    end

end

That's about it! The speed is insane, at least other CPU-based game of life impelmentations I've written in the past. At a grid size that's native to my screen (1680x1050) I get over 200 iterations per second. That's 1,764,000 pixels updated at 200 times a second, which is something like 3 nanoseconds per-pixel. Pretty nuts.

This video doesn't really do it justice with all of the compression (watch out for aliasing if you watch it really small).

I've been trying to get my youngest brother to start making his own video games instead of spending all of his time playing them. I sat down with him over this winter break to show him how to write some Lua for use with Löve, a 2d game engine.

So the first thing I thought we'd write is an asteriods clone. The code really isn't worth posting, but here's a little screenshot.

We actually had a lot of fun implementing all of the basic asteriods stuff. The asteriods are randomly generated in size and shape. The ship can be piloted around with the classic controls, and we even added a little red triangle as the rocket that shows up with you're using the thruster.

I wanted to show my brother that the best part about writing your own software and games is being able to make the computer kind of do whatever you want. We started goofing around by making the missiles super long (or are they laser beams?), which gave us some funny effects, and then we found that you could get some cool results by imparting no extra velocity on the missiles. Typically the when you want to shoot a missile from the ship you add some missile velocity to the ship's velocity and that's the initial velocity of the missile. With no missile velocity, if the ship was moving along at constant speed and a missile was shot, then the missile would just continue floating along with the ship.

At this point, I was looking for some cooler stuff to show him, that that's when I stumbled upon Love's new pixel shaders feature. With this, we wrote a quick shader that covers whatever is drawn on the black blackground with a spinning rainbow wheel. Here's the shader code.

extern number screen_width;
extern number screen_height;
extern number time;

#define PI 3.14159

vec4 hsvToRgb(float hue, float saturation, float value, float alpha) {

    float chroma = value * saturation;
    float h = hue / 60.0;
    float x = chroma * (1 - abs(mod(h, 2) - 1));
    vec3 intermediate = vec3(0, 0, 0);
    if (0 <= h && h < 1) {
        intermediate = vec3(chroma, x, 0);
    } else if (1 <= h && h < 2) {
        intermediate = vec3(x, chroma, 0);
    } else if (2 <= h && h < 3) {
        intermediate = vec3(0, chroma, x);
    } else if (3 <= h && h < 4) {
        intermediate = vec3(0, x, chroma);
    } else if (4 <= h && h < 5) {
        intermediate = vec3(x, 0, chroma);
    } else if (5 <= h && h < 6) {
        intermediate = vec3(chroma, 0, x);
    }
    float m = value - chroma;
    return vec4(intermediate[0] + m, intermediate[1] + m, intermediate[2] + m, alpha);
}

vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 screen_coords ){
    float x = ((screen_coords.x / screen_width) * 2) - 1;
    float y = ((screen_coords.y / screen_height) * 2) - 1;

    if (color == vec4(1.0, 0, 0, 1.0)) {
        return color;
    }

    if (x != 0 && y != 0) {
        float angle = atan(y, x);
        if (angle < 0) {
            angle += 2 * PI;
        }
        angle = mod((angle * (180 / PI)) + time * 10, 360);
        float dist = sqrt(pow(x, 2) + pow(y, 2)) / sqrt(2);
        float saturation = (((sin(time) / 2) + 0.5) * 0.1) + 0.9;
        return hsvToRgb(angle, dist, 1.0, 1.0);
    } else {
        return vec4(1.0, 1.0, 1.0, sin(time));
    }
}

There's a small function to convert from the HSV space to RGB space. There has to be a better way of doing this but this is the algorithm I've always used (it's in the Wikipedia page for the HSV color space). The actual shader code is calculating the angle from the center of the screen, and then using this to set the hue of the pixel. The distance from the center sets the saturation, and then the whole thing is rotating in time, creating this spiral color wheel effect. If the original pixel was black then it simply returns black. This means that the only pixels that get colored are where things have been drawn. This was really great to get him involved, since normally dealing with OpenGL mainly involves mucking about with a library, or writing a few hundred lines of boilerplate code to compile the shaders and setup framebuffers and whatnot.

Here's the result of the weird missiles and the shader.

What's really neat is that the shader only acts wherever you have drawn objects on the screen. If you just want to see what the output of the shader is, a bounds-sized rectangle can be drawn, and all operations can be done in the shader. However, this allows for "masking" of the shader. So here's a mandlebrot shader acting on two moving circles. In this case, the mandlebrot shader uses the screenCords that's passed to the pixel shader.

A friend challenged me to make the mandlebrots inside individual circles that are moving independently of one another. For this the textureCoords shader variable can be used. In this case, you have to draw a texture. I used a 1x1 pixel image, and then limited limited the drawing in the shader to a circular shape by first testing the distance from the center of the texture before coloring any pixels. What's neat is that the texture rotation is taken into account (not 100% sure how), so if you draw the texture rotated using the love.graphics.draw method, whatever the shader does is inherently rotated. The code for this demonstration is below the video, if you're interested with playing with it.

extern number screen_width;
extern number screen_height;

vec4 effect(vec4 color, Image texture, vec2 texture_coords, vec2 screen_coords) {
  if(sqrt(pow(texture_coords.x * 2.0 - 1.0, 2) + pow(texture_coords.y * 2.0 - 1.0, 2)) > 1) {
    return vec4(0.0, 0.0, 0.0, 0.0);
  }

  float x0 = (texture_coords.x * 3.35) - 2.5;
  float y0 = (texture_coords.y * 2) - 1;

  float x = 0.0;
  float y = 0.0;
  int iteration = 0;
  int max_iterations = 100;

  while((pow(x, 2) + pow(y, 2) < 256) && iteration < max_iterations) {
    float xtemp = pow(x, 2) - pow(y, 2) + x0;
    y = 2 * x * y + y0;
    x = xtemp;
    iteration += 1;
  }

  if (iteration < 950) {
    return vec4(1.0, iteration / 30.0, iteration / 10.0, 1.0);
  } else {
    return vec4(0.0, 0.0, 0.0, 1.0);
  }

  return vec4(screen_coords[0] / screen_width, screen_coords[1] / screen_height, 0.2, 1.0);
}

shader = ""
screenWidth = 0
screenHeight = 0
imageWidth = 0
imageHeight = 0
image = ""

function love.load(a)
    shader = love.graphics.newShader("shader.frag")
    screenWidth = love.graphics.getWidth()
    screenHeight = love.graphics.getHeight()
    image = love.graphics.newImage("tex.png")
    imageWidth = image:getWidth()
    imageHeight = image:getHeight()
end

t = 0

function love.update(dt)
    t = t + dt
end

function love.draw(dt)
    love.graphics.setShader(shader)
    local prefs = {
        {3, 4, 8, 1, 0.2},
        {1, 5, 4, 1.2, 0.2},
        {4, 5, 3, 0.8, 0.6},
        {1, 2, 2, 2.1, 0.1},
        {1, 4, 10, 20, 20}
    }
    for i = 1, 5 do
        local pref = prefs[i]
        local x1 = (math.sin(pref[1] * t) * 0.5 + 0.5) * (screenWidth - imageWidth)
        local y1 = (math.cos(pref[2] * t) * 0.5 + 0.5) * (screenHeight - imageHeight)
        local scale1 = ((math.sin(pref[3] * t) * 0.5 + 0.5) * pref[5]) + 1
        love.graphics.draw(image, x1, y1, t * pref[4], 200)
    end
end

I wanted to try implementing one of the smooth coloring techniques, and I'm not really a master of complex analysis so, I pretty much just implemented the psuedocode in the Wikipedia article listed under the "continuous smooth coloring" section. It utilizes the magnitude of the iteration (which, for points in the plane that you've decided to "bail out", is equal to the bailout radius) and the iteration at which this bailout radius was achieved. This is then used to choose a color out of a nice list of colors.

I generated lists of smoothly varying colors using this HCL picker, which lets you drag a line between two colors on a nice HCL color plot. I clicked the + a bunch of times to generate a few hundred points along the line, and then hit the copy button. This copies the list as hex values. I wrote a small Lua script to convert the hex values to a static array in Lua syntax, that could be easily handed off to the shader.

A strange feature of Löve here is the love.filesystem.load method. Looking into it more it's probably using the Lua dofile builtin. It allowed me to create an almost JSON-like file of data. The colorData.lua file simply defines a global colors list that is accessable in the global scope after the file is executed.

--  load and execute the color data, which creates a "colors" global
local colorData = love.filesystem.load("colorData.lua")
colorData()

I wrote a small user interface that enables the user to change the zoom speed using the scroll wheel, as well as cycle through a few different color pallets. It also allows adjustable oversampling to achieve a configurable amount of anti-aliasing. The fastest rendering occurs when this is disabled. I also added the ability to export the current viewport at a very large size to an image file.

The trickiest part was figuring out how to zoom in or out at the point the mouse was hovering in the complex plane. It ended up being a simple calculation. The first step was to take the mouse coordinates in the window, convert them to a coordinate in the complex plane, shrink the domain, and then do the inverse calculation to ensure the point in the complex plane is still at the same window coordinate after the domain has shrunk. Zooming out is an identical calculation, but the domain becomes larger instead of smaller. The relevant code is below.

if love.mouse.isDown(1, 2) then

    --  find the x and y point in the current fractal domain
    local x = love.mouse.getX()
    local y = love.mouse.getY()
    local range = domain / aspectRatio
    local xClick = xCenter + ((x / windowWidth) * domain) - (domain / 2)
    local yClick = yCenter + ((y / windowHeight) * range) - (range / 2)

    --  are we zooming in or out?
    if love.mouse.isDown(1) then
        domain = domain * zoomSpeed
    elseif love.mouse.isDown(2) then
        domain = domain / zoomSpeed
    end

    --  need to set a new center so that this point stays at the
    --  location of the click, instead of zooming in at the center, this is
    --  essentially the inverse of the previous calculation, just using
    --  the new domain after zooming in or out
    local range = domain / aspectRatio
    xCenter = xClick - ((x / windowWidth) * domain) + (domain / 2)
    yCenter = yClick - ((y / windowHeight) * range) + (range / 2)

    --  update the fractal in the offscreen buffer
    drawFractal(fractalCanvasWidth, fractalCanvasHeight, fractalCanvas)
end

To avoid re-drawing the factal every frame, the fractal is rendered into an offscreen buffer and only updated when the user modifies the display in any way (zooming in, panning, changing color pallette). The drawFractal() function does the drawing into the offscreen buffer.

function drawFractal(width, height, drawCanvas, sendSize)

    --  send all of the information to the shader
    print(width, height, xCenter, yCenter, domain)
    print(drawCanvas)
    print(fractalCanvas)
    shader:send("window_width", width)
    shader:send("window_height", height)
    shader:send("x_center", xCenter)
    shader:send("y_center", yCenter)
    shader:send("domain", domain)

    --  draw the fractal to the drawCanvas using the shader by filling the
    --  canvas with a filled rectangle
    love.graphics.setCanvas(drawCanvas)
    love.graphics.setShader(shader)
    love.graphics.rectangle("fill", 0, 0, width, height)
    love.graphics.setShader()
    love.graphics.setCanvas()

end

This keeps the love.draw() method short and sweet. If the fractal canvas was oversized to accomplish the oversampling, the love.graphics.draw() method downsizes the buffer while drawing it to the screen.

function love.draw(dt)

    --  draw the fractal canvas, rescaling it to the window's size
    love.graphics.draw(fractalCanvas, 0, 0, 0, 1 / oversample)

    --  print the current zoom speed to the screen
    if displayZoomSpeed then
        local x = 5
        local y = windowHeight - 18
        love.graphics.print("zoom speed : " .. zoomSpeed, x, y, 0, 0.8)
    end

end

Here are some images and a video of the little program. Much prettier (and faster) than my previous attempt with OpenCL.

The limitation now is the finite precicion of the floating point math done inside the shader. I can't think of an easy workaround, my first idea would be to do arbitrary precision floating point math in software, somehow in an OpenGL shader. This might be a brute-force solution though, there might be more mathematical way around it.

As always, the full code listing is on github. I also made a post on the Löve2D forums here. The members there help fix some of the shader code to make it more platform independent. I also got it running on my iPhone, which is neat.

I've always kind of wondered how 3D renderers do their magic. A few years ago I found a book in a relative's library from the 1990's on some basic 3D computer graphics techniques, and since then I've been itching to write a simple software rendering engine. The book is 3D Computer Graphics by Alan Watt (2nd Edition), there's a 3rd edition out now that you might be able to get your hands on. The book was first printed in 1989, and thus doesn't include anything about hardware rendering or giant 3D graphics libraries, which I found really nice. I've done some stuff with OpenGL before, but since most of the tutorials and guides I've read are simply describing how OpenGL works, they don't usually spend the time to explain most of the the mathematics that are involved.

When I got the book for the first time, I made this little experiment, simply trying to figure out where the eight points in a cube end up after a projection from 3D space into screen space, but something wasn't quite right (if you can notice in the video below), and I kind of gave up on it then and hadn't tried again until now.

This entry is not really intended to be a tutorial of sorts, since there are plenty of them, but more of just a log of some of the fun I had in learning this stuff. If you're interested in a tutorial, there's a great series here by David Rousset that has a lot of the fundamentals.

matlab experiementation

Because my previous experiment had some issues I could never track down and I was never sure where in the pipeline the issues were, I wanted to make sure I had some of the main concepts down before writing a lot of my own matrix math code. I figured Matlab has a nice line function built in, and it has all of the matrix operations built into the langugage, so I figured it'd be a good starting point (even though a lot of my peers tend to bemoan Matlab's existance). I suppose I could have done it in any other langugage, but I didn't want to muck around with matrix libraries and such and was itching to get my feet wet.

I first wrote some functions to give back some standard transformation matricies for quick use. These are ugly, but I was really just looking to get something working to prove the math behind it all.

function [matrix] = perspectiveMatrix(r, l, t, b, f, n)
    matrix = [
        2/(r-l) 0 0 -(r+l)/(r-l) ; 
        0 2/(t-b) 0 -(t+b)/(t-b) ;
        0 0 1/(f-n) -n/(f-n) ;
        0 0 0 1
        ];
end

function [matrix] = xrotationMatrix(deg)
    rad = deg * pi / 180;
    matrix = [
        1 0 0 0 ;
        0 cos(rad) sin(rad) 0 ;
        0 -sin(rad) cos(rad) 0 ;
        0 0 0 1
        ];
end

function [matrix] = yrotationMatrix(deg)
    rad = deg * pi / 180;
    matrix = [
        cos(rad) 0 -sin(rad) 0 ;
        0 1 0 0 ;
        sin(rad) 0 cos(rad) 0 ;
        0 0 0 1
        ];
end

function [matrix] = zrotationMatrix(deg)
    rad = deg * pi / 180;
    matrix = [
        cos(rad) sin(rad) 0 0 ;
        -sin(rad) cos(rad) 0 0 ;
        0 0 1 0 ;
        0 0 0 1
        ];
end

function [matrix] = translationMatrix(x, y, z)
    matrix = [
        1 0 0 0 ;
        0 1 0 0 ;
        0 0 1 0 ;
        x y z 1
        ];
end

The meat of the code is below.

function a = draw(h, event)
    ClearLinesFromAxes()
    lines = [
        [[-50 -50 -100], [-50 50 -50]],       
        [[-50 50 -50], [50 50 -50]],
        [[50 50 -50], [50 -50 -50]],
        [[50 -50 -50], [-50 -50 -50]],
        %
        % more line pairs here, but removed for conciseness
        %

    ];

    lines = reshape(lines, [length(lines), 3, 2]);
    length(lines);
    t = perspectiveMatrix(-100, 100, -100, 100, -100, 100);
    data = guidata(gcf);
    rx = yrotationMatrix(data(1).Value);
    ry = xrotationMatrix(data(2).Value);
    rz = zrotationMatrix(data(3).Value);
    move = translationMatrix(data(4).Value, data(5).Value, data(6).Value);
    for m = [1: length(lines)]
        pointa = [lines(m, :, 1), 1]';
        pointb = [lines(m, :, 2), 1]';
        if (length(lines) - m) <= 44
            transforma = (pointa' * (rz * rx * ry * t))';
            transformb = (pointb' * (rz * rx * ry * t))';
        else
            transforma = ((pointa' * move) * (rz * rx * ry) * t)';
            transformb = ((pointb' * move) * (rz * rx * ry) * t)';
        end
        xya = transforma(1:2);
        xyb = transformb(1:2);
        x = [xya(1) xyb(1)];
        y = [xya(2) xyb(2)];
        l = line(x, y);
        i = length(lines) - m;
        if i == 0
            set(l, 'Color', 'red');
        elseif i == 1
            set(l, 'Color', 'magenta');
        elseif i == 2
            set(l, 'Color', 'green');
        elseif i <= 44
            set(l, 'Color', [.9, .9, .9]);
        end
    end

    xlim([-2, 2])
    ylim([-2, 2])
end

There's some GUI setup boilerplate not shown. There are six sliders shown on top of the plot. Three shift all of the objects around in space, and the other three adjust the transformation matrix to "move" the camera around. The camera tranformation is summed up into the translationMatrix function. The rest of the code is just going through all of the sets of points, transforming them into 2D screen space, and then using the Matlab line(x, y) function to draw then onto a plot. The draw function is called every time a slider value changes. The only weird thing is coloring the first three lines in the lines list different colors. This let me have three different colored axes to make sure everything made sense.

Here are two demos of this script. The first is a simple pyramid shape, and the second is after I added a few more structures and the axes to the output.

writing a matrix in swift (eg. getting sidetracked)

I've been wanting to get deeper into Swift too. So I figured this would be a good candidate project. I looked into a few Swift matrix libaries and really didn't see anyting I liked, so I rolled my own, and what started an "okay, I just want to be able to create and multiply some matricies" project turned into "this thing is like 500 lines now" undertaking.

I don't want to get too detailed on it here, so you can check Github for the full code listing (it's just one file) and the project itself for some more information. Any source code posted here will probably have little repetitive tidbits removed.

All elements in the matrix are doubles, I decided not to go with a generic implementation to simplify matters a bit, and are stored column-major ordered inside a Swift array. For the most part, I copied a lot of the nice functions that Matlab has built in, and tried to maximize the amount of code reuse in the library itself. For example, instead of having two funcitons ones and zeros (which return a prefilled MxN matrix with ones and zeros respectively) that both handle creating a new matrix and prefilling it, I wrote one function prefilled which returns a matrix with MxN elements that are all equal to one number, which is called by ones and zeros. Elementary I know, but I was proud of how clean some of the code came out to be.

//  returns a rows*cols matrix with all elements equal to 'prefill'
class func prefilled(rows: Int, _ cols: Int, pad: Double) -> Matrix {
    let elements = [Double](count: rows * cols, repeatedValue: pad)
    return Matrix.init(nr: rows, nc: cols, elements: elements)
}

//  returns a square prefilled matrix
class func prefilled(size: Int, pad: Double) -> Matrix {
    return Matrix.prefilled(size, size, pad: pad)
}

//  returns a rows*cols matrix filled with zeros
class func zeros(rows: Int, _ cols: Int) -> Matrix {
    return Matrix.prefilled(rows, cols, pad: 0.0)
}

//  returns a square zeros matrix
class func zeros(size: Int) -> Matrix {
    return Matrix.zeros(size, size)
}

//  returns a rows*cols matrix filled with ones
class func ones(rows: Int, _ cols: Int) -> Matrix {
    return Matrix.prefilled(rows, cols, pad: 1.0)
}

//  returns a square ones matrix
class func ones(size: Int) -> Matrix {
    return Matrix.ones(size, size)
}

Here are three more methods that I wrote in a similar way, one top-level method diagonal allows the creation of a diagonally filled matrix with the non-diagonal elements set to any arbitray number. it is then overridden to allow a default value of 0. The identity matrix function then uses the square-matrix diagonal method.

//  returns a fill.count*fill.count matrix with the diagonal elements equal to the elements in
//  fill, from top-left to bottom-right, and all other elements equal to 'pad'
class func diagonal(fill: [Double], pad: Double) -> Matrix {
    let result = Matrix.prefilled(fill.count, fill.count, pad: pad)
    for (row, element) in fill.enumerate() {
        result.elements[(row * fill.count) + row] = element
    }
    return result
}

//  same as diagonal(fill:, pad:), but uses a default pad of zero
class func diagonal(fill: [Double]) -> Matrix {
    return Matrix.diagonal(fill, pad: 0.0)
}

//  returns an size*size sized identity matrix
class func I(size: Int) -> Matrix {
    let ones = [Double](count: size, repeatedValue: 1.0)
    return Matrix.diagonal(ones)
}

Again, in the same vein, instead of writing four functions to perform the four elementary operations on an element-wise basis, I wrote one function that can "combine" two equally sized matricies element-by-element, that takes any function that accepts two doubles and returns a double. Then, the four elementary operations use this base function instead. Similarly, I wrote a map function that does the same, except it maps one matrix to a new matrix through a uinary operation.

//  apply a function to every element
func map(function: (Double) -> Double) -> Matrix {
    let result = Matrix.zeros(nr, nc)
    for (index, element) in elements.enumerate() {
        result.elements[index] = function(element)
    }
    return result
}

//  combine two equally sized matrices into one using a function
func combine(matrix: Matrix, operation: (Double, Double) -> Double) -> Matrix {
    assert(Matrix.dimensionsEqual(self, matrix), "Dimension Mismatch: Dimensions Must Equal")
    let result = Matrix.zeros(nr, nc)
    for (index, element) in elements.enumerate() {
        result.elements[index] = operation(element, matrix.elements[index])
    }
    return result
}

//  element-wise addition, subtraction, multiplication and division
func add(matrix: Matrix) -> Matrix {
    return self.combine(matrix, operation: +)
}

//  add, subtract, multiply and divide by a scalar
func add(scalar: Double) -> Matrix {
    return self.map{$0 + scalar}
}

Using Swift's operator overloading features was really nice as well. I defined the * infix operator on two matricies as matrix multiplication (much like Matlab) and then defined ∆* and ∆/ as element-by-element multiplication and division. However, I have not yet implemented matrix "division", because I have not implemented matrix inversion.

//  add two equivalently sized matrices, element by element
func + (left: Matrix, right: Matrix) -> Matrix {
    return left.add(right)
}

//  element-wise multiplication
infix operator ∆* { associativity left precedence 150 }
func ∆* (left: Matrix, right: Matrix) -> Matrix {
    return left.multiplyElements(right)
}

//  scalar addition
func + (m: Matrix, s: Double) -> Matrix {
    return m.add(s)
}

I also overloaded (and this is probably the coolest bit) the subscript accessor methods. I wanted to be able to index the matricies by any permutation of numbers, arrays, ranges and StrideThroughs. This allows for getting (and setting portions of) new matricies by writing mat(1..5, [1, 2, 9]) or mat(5, stride(from: 1, to: 9, by: 3)) and such. However, I didn't want to rewrite all of these things a billion time for every combination of the two datatypes, so instead every overloaded instance of the subscript accessors simply cast whatever arguments they have into arrays, and then call one main function.

// this array subscript accessor is the basis for the rest of the other accessors, it allows for
// getting and setting groupings of rows and columns in any order, when setting, the fill matrix
// must either be equal to the size of the submatrix that is being set, or a singular element
// that is to be filled in all spaces defined in 'rows' and 'cols'
subscript(rows: Array<Int>, cols: Array<Int>) -> Matrix {
    get {
        let rStart = rows.minElement()!
        let rEnd = rows.maxElement()!
        let cStart = cols.minElement()!
        let cEnd = rows.maxElement()!
        assert(rangeInBounds(rStart...rEnd, cStart...cEnd), "Getting [Rows, Cols] Out Of Bounds")
        let result = Matrix.zeros(rows.count, cols.count)
        for (placeRow, getRow) in rows.enumerate() {
            for (placeCol, getCol) in cols.enumerate() {
                result.elements[(placeRow * cols.count) + placeCol] = elements[(getRow * nc) + getCol]
            }
        }
        return result
    }
    set (fill) {
        let rStart = rows.minElement()!
        let rEnd = rows.maxElement()!
        let cStart = cols.minElement()!
        let cEnd = rows.maxElement()!
        assert(rangeInBounds(rStart...rEnd, cStart...cEnd), "Setting [Rows, Cols] Out Of Bounds")
        let singular = fill.elements.count == 1
        if !singular {
            assert(fill.rangeInBounds(0...(rows.count - 1), 0...(cols.count - 1)), "Supplied Matrix Is Insufficiently Sized")
        }
        for (getRow, placeRow) in rows.enumerate() {
            for (getCol, placeCol) in cols.enumerate() {
                elements[(placeRow * nc) + placeCol] = singular ? fill.elements[0] : fill.elements[(getRow * fill.nc) + getCol]
            }
        }
    }
}

//  rows and cols can be indexed by ranges or by integer strides
subscript(rows: Range<Int>, cols: Range<Int>) -> Matrix {
    get {
        return self[[Int](rows), [Int](cols)]
    }
    set(fill) {
        self[[Int](rows), [Int](cols)] = fill
    }
}

I think the coolest "functional" code I wrote though was in the matrix multiplication method. I wrote a sum method and a transpose method first, and then wrote this. Note, allRows and allCols are methods that return a range that specifies the minimum and maximum row or column indicies for the matrix. The "summing" operation is one line which uses the subscript accessors to get the correct LHS column matrix, and then multiplies its elements by the transposed (turns a row matrix into a column matrix) version of the correct RHS row, and then takes the sum of the result. I thought it was quite beautiful, at least compared to an even deeper nested loop.

// multiply two matrices
func multiply(matrix: Matrix) -> Matrix {
    assert(nc == matrix.nr, "Multiplication Dimension Mismatch")
    let result = Matrix.zeros(nr, matrix.nc)
    for row in result.allRows {
        for col in result.allCols {
            let value = self[row, allCols].multiplyElements(matrix[matrix.allRows, col].transpose()).sum
            result[row, col] = value
        }
    }
    return result
}

actually rendering some stuff now

I really hate mucking with CoreGraphics and NSView's drawRect(), so I wrote a quick little "render buffer" which has the ability to set individual pixel values to colors, and buffer updates (so all pixels get pushed at once to the screen). I wanted this to be cross platform (well, to the iPhone at least) so I wrote a little interface that could potentially be implemented in UIKit as well. This is hopefully the only UIKit/Cocoa specific code in the renderer, so we'll see how it goes and if I end up porting to iOS in the future.

///////////////////////////////////////////////////////////////////////////////////////////////////
//  MARK: RenderSurfaceProtocol
///////////////////////////////////////////////////////////////////////////////////////////////////
//  In order for the renderer to be cross platform, this protocol is used to represent something
//  that can render pixels to a display, it should allow for one state to be saved and
//  continuously drawn, while new pixel data is given to it, when draw() is called, the new pixel
//  data should be shown on the display, but the buffer should not be cleared until clear is called

protocol RenderSurfaceProtocol {
    //  returns the size of the buffer in pixels
    var size: CGSize { get }
    
    init(bufferSize: CGSize)

    //  sets a pixel in the buffer to an RGB color
    func setPixel(x: Int, _ y: Int, _ r: Int, _ g: Int, _ b: Int)

    //  clears all pixels to black
    func clear()

    //  clears all pixels to a given color
    func clearToColor(r: Int, _ g: Int, _ b: Int)

    //  puts the current buffer on the screen, does not clear the buffer
    func draw()
}

I don't think I'll bore you with the actual implementation. I did have some weird issues though when doing this. Honestly I can't remember what it was that caused only some of the screen to update, but I'm sure it was something dumb. Here's a video so you can laugh at the hilarity. I was drawing random circles in the buffer continuously, and for some reason it only updates in a specific region. However, when the window is being resized, it seemed to draw in the right place, until the mouse was released and it went back to acting strange. I fixed it though, so it's all good.

From here, I kind of just went about copying what I had done in Matlab. I did implement a very awful version of Bresenham's line algorithm because I only had a buffer that would draw individual pixels. Note the hardcoding of the viewport bounds.

func drawLine(x1: Int, _ y1: Int, _ x2: Int, _ y2: Int, _ r: Int, _ g: Int, _ b: Int) {
    if(x2 < x1) {
        drawLine(x2, y2, x1, y1, r, g, b)
    } else {
        let dx = x2 - x1
        let dy = y2 - y1
        if dx > 0 {
            for x in x1...x2 {
                let slope: Double = Double(dy) / Double(dx)
                let y: Double = (slope * Double(x - x1)) + Double(y1)
                if(x > 640 || x < 0 || y > 480 || y < 0) {
                    
                } else {
                    renderSurface.setPixel(x, Int(y), r, g, b)
                }
            }
        } else {
            if y2 > y1 {
                for y in y1...y2 {
                    renderSurface.setPixel(x1, y, r, g, b)
                }
            } else {
                for y in y2...y1 {
                    renderSurface.setPixel(x1, y, r, g, b)
                }
            }
        }

    }
}

To make things a little more extensible, I created a Renderable interface, which simply described something that had a collection of triangular faces. In order for it to be dynamic, I added a construct() method which would have the item fill its list of faces.

//  Describes an object that is renderable, renderable objects should have a list of verticies and
//  faces, a boolean which represents if the verticies and faces lists have been filled 
//  ('constructed') and a construct function which calculates the verticies and faces, filling 
//  the respective lists

protocol Renderable {
    
    //  the verticies that comprise the object, should be row vectors (in an array), with
    //  three items, x, y, and z respectively
    var vertices: [[Double]] { get }
    
    //  the faces of the object, should be a list of 3-item integer lists, refering to indicies
    //  in the verticies list
    var faces: [[Int]] { get }
    
    //  the origin of the object as a row matrix [[x, y, z]]
    var origin: Matrix { get }
    
    //  describes if the object has had it's verticis and faces lists filled, if an object is not
    //  constructed at init, this should be false until construct is called, and falsified if
    //  a property is modified that would modify the verticies or faces
    var constructed: Bool { get }
    
    //  fills the verticies and faces matricies
    func construct()
    
}

Here's an example of a pyramid object that implements the interface. In this case it's hand coded and really simple.

class Pyramid: BaseRenderableObject {
    
    override func construct() {
        construct(50.0)
    }
    
    func construct(size: Double) {
        let d: Double = size / 2.0
        //fill the verticies array
        vertices = [
            [-d,  0, -d],
            [-d,  0,  d],
            [ d,  0,  d],
            [ d,  0, -d],
            [ 0,  d,  0],
        ]
        
        faces = [
            //  bottom
            [0, 2, 1],
            [0, 2, 3],
            
            //  sides
            [0, 4, 1],
            [1, 4, 2],
            [2, 4, 3],
            [3, 4, 0]
        ]
        
        constructed = true
    }
    
}

Here's an example of a plane (like the flat surface, not an airplane) object that has an adjustable subdivision property, meaning it can either consist of just two triangles, or subdivide it's size into smaller "grid" subsections for more detail. Note I was kind of naughty here and didn't implement the construct() method and instead opted to call a different construct method manually. I should fix this later, but it works for now.

class Plane: BaseRenderableObject {
   
   override func construct() {
       
   }
   
   func construct(size: Double, subdivisions: Int, primaryAxis: Int, secondaryAxis: Int) {
       let d: Double = size / 2.0
       //fill the verticies array
       for primaryIndex in 0...subdivisions {
           let i: Double = -d + ((size/Double(subdivisions)) * Double(primaryIndex))
           for secondaryIndex in 0...subdivisions {
               let j: Double = -d + ((size/Double(subdivisions)) * Double(secondaryIndex))
               var vertex: [Double] = [0.0, 0.0, 0.0]
               vertex[primaryAxis] = i
               vertex[secondaryAxis] = j
               vertices.append(vertex)
               if primaryIndex != subdivisions && secondaryIndex != subdivisions {
                   let tl = (primaryIndex * (subdivisions + 1)) + secondaryIndex
                   let bl = ((primaryIndex + 1) * (subdivisions + 1)) + secondaryIndex
                   print([tl + 1, bl, tl])
                   print([tl + 1, bl, bl + 1])
                   faces.append([tl + 1, bl, tl])
                   faces.append([tl + 1, bl, bl + 1])
               }
           }
       }   
       constructed = true
   }    
}

The main rendering code is largely the same as the Matlab code, just even uglier. It caries out the same functions. I may make it public later after cleaning it up a bit.

Here's some videos from the process. Hilariously, the first result I got in Swift was the same thing I got when I tried doing this three years ago, where my cube would mysteriously "squish" down into a square. I tried forever thinking about how you could look at a cube and have it appear as just a square, and I think it's imposisble. I'm not 100% in the video the dots are even visible with the compression.

Just for grins I decided to connect these incorrectly located dots with lines, for humor's sake.

I reviewed the code some more, and fixed the bug.

This is when I added the planes. I also added some sliders to mess with the camera location.

Finally I added perspective and some more sliders. This took me a bit of time to wrap my head around how it works, but I finally got it (except for the bugs).

So this is it for now. I'd like to actually fill the triangles now to make the wireframe views into something more substantial, but I think to get there I'm going to have to do some major refactoring and optimization. It'd be fun to implement a z-buffer and then some backface culling. I'd also like to implement some basic lighting and a fly around camera. If nothing else, this has been a really fun introduction into 3D computer graphics (and some more advanced Swift) and now I think I'd feel comfortable applying some of this math in OpenGL shaders.

If you'd like the code send me an email, it's not really polished and really in no state to be put on Github.