Raytracing on a grid, in Haskell

I attempt to write something in Haskell every now and then. Mainly this is because I think we're going to have to move beyond the imperative model in order to effectively use our increasingly parallel computers. Languages like C are built around the assumption that there is a single locus of control which is reading data from memory, modifying it, and writing it back. Multi-threaded programming is possible in C but it's devilishly difficult.

Haskell's a good-looking language, in my opinion; it is parenthesized like Scheme, but due to its equational and algebraic syntax needs far less of them. It uses layout like Python to indicate block structure, which eliminates even more syntactic verbiage. It's got static typing with type inference, so the compiler catches my stupid mistakes with a minimum of additional typing. In general, when the program compiles, it works, which is a very different experience from coding in, say, Scheme or Python.

For fun, I ported my grid raytracing test program from C++ to Haskell. Here's the original C++ code:

#include <cmath>
#include <windows.h>
#include <gl/glut.h>

static void draw_frame();
static void window_was_resized(int width, int height);
static void key_pressed(unsigned char key, int x, int y);
static void mouse_activity(int button, int state, int x, int y);

static const initial_window_size_x = 800;
static const initial_window_size_y = 800;

static const pixels_per_unit = 64;

static int g_window_width  = 0;
static int g_window_height = 0;

static double g_line_x0 = 4;
static double g_line_y0 = 1;

static double g_line_x1 = 1;
static double g_line_y1 = 1;

static int g_nr_points_defined = 2;

int main(int argc, char * argv[])
{
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
    glutInitWindowPosition
    (
        (GetSystemMetrics (SM_CXSCREEN) - initial_window_size_x) / 2,
        (GetSystemMetrics (SM_CYSCREEN) - initial_window_size_y) / 2
    );
    glutInitWindowSize(initial_window_size_x, initial_window_size_y);
    glutCreateWindow("Line Draw Test");

    glutReshapeFunc(window_was_resized);
    glutDisplayFunc(draw_frame);
    glutKeyboardFunc(key_pressed);
    glutMouseFunc(mouse_activity);

    glEnable(GL_COLOR_MATERIAL);

    glutMainLoop();

    return 0;
}

static void plot_squares(double x0, double y0, double x1, double y1)
{
    double dx = fabs(x1 - x0);
    double dy = fabs(y1 - y0);
    int x = int(floor(x0));
    int y = int(floor(y0));

    int n = 1;
    int x_inc, y_inc;
    double error;

    if (x1 > x0)
    {
        x_inc = 1;
        n += int(floor(x1)) - x;
        error = (ceil(x0) - x0) * dy;
    }
    else
    {
        x_inc = -1;
        n += x - int(floor(x1));
        error = (x0 - floor(x0)) * dy;
    }

    if (y1 > y0)
    {
        y_inc = 1;
        n += int(floor(y1)) - y;
        error -= (ceil(y0) - y0) * dx;
    }
    else
    {
        y_inc = -1;
        n += y - int(floor(y1));
        error -= (y0 - floor(y0)) * dx;
    }

    glColor3d(0.25, 0.25, 0.5);
    glBegin(GL_QUADS);

    for (; n > 0; --n)
    {
        glVertex2i(x,   y  );
        glVertex2i(x+1, y  );
        glVertex2i(x+1, y+1);
        glVertex2i(x,   y+1);

        if (error > 0)
        {
            y += y_inc;
            error -= dx;
        }
        else
        {
            x += x_inc;
            error += dy;
        }
    }

    glEnd();
}

static void draw_frame()
{
    glClear(GL_COLOR_BUFFER_BIT);

    // Set up the projection matrix.

    double window_world_size_x = double(g_window_width ) / pixels_per_unit;
    double window_world_size_y = double(g_window_height) / pixels_per_unit;

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluOrtho2D(0, window_world_size_x, 0, window_world_size_y);

    // Plot the line's squares if the line is defined.

    if (g_nr_points_defined >= 2)
    {
        plot_squares(g_line_x0, g_line_y0, g_line_x1, g_line_y1);
    }

    // Draw a background grid.

    int number_of_lines_x = int(ceil(window_world_size_x));
    int number_of_lines_y = int(ceil(window_world_size_y));

    glColor3d(0.5, 0.5, 0.5);
    glBegin(GL_LINES);

    int i;
    for (i = 0; i < number_of_lines_x; ++ i)
    {
        glVertex2d(i, 0);
        glVertex2d(i, window_world_size_y);
    }

    for (i = 0; i < number_of_lines_y; ++ i)
    {
        glVertex2d(0, i);
        glVertex2d(window_world_size_x, i);
    }

    glEnd();

    // Draw the actual line, if defined.

    if (g_nr_points_defined == 1)
    {
        glColor3d(1, 1, 1);
        glBegin(GL_POINTS);
        glVertex2d(g_line_x0, g_line_y0);
        glEnd();
    }
    else if (g_nr_points_defined >= 2)
    {
        glColor3d(1, 1, 1);
        glBegin(GL_LINES);
        glVertex2d(g_line_x0, g_line_y0);
        glVertex2d(g_line_x1, g_line_y1);
        glEnd();
    }

    glutSwapBuffers();
}

static void window_was_resized(int width, int height)
{
    g_window_width  = width;
    g_window_height = height;

    glViewport(0, 0, width, height);
}

static void key_pressed(unsigned char key, int x, int y)
{
    if (key == 27)
    {
        exit (0);
    }
}

static void mouse_activity(int button, int state, int x, int y)
{
    if (state == GLUT_DOWN)
    {
        if (button == GLUT_LEFT_BUTTON)
        {
            // Get the coordinates of the mouse click in world coordinates.

            GLdouble modelview_matrix[16];
            GLdouble projection_matrix[16];
            GLint viewport[4];

            glGetDoublev(GL_MODELVIEW_MATRIX, modelview_matrix);
            glGetDoublev(GL_PROJECTION_MATRIX, projection_matrix);
            glGetIntegerv(GL_VIEWPORT, viewport);

            double wx, wy, wz;

            gluUnProject(x, (g_window_height - 1) - y, 0, modelview_matrix,
                projection_matrix, viewport, &wx, &wy, &wz);

            if (g_nr_points_defined == 1)
            {
                g_line_x1 = wx;
                g_line_y1 = wy;
                g_nr_points_defined = 2;
            }
            else
            {
                g_line_x0 = wx;
                g_line_y0 = wy;
                g_nr_points_defined = 1;
            }

            glutPostRedisplay();
        }
    }
}

Here is the same program written in Haskell:

-- Raytracing on a grid

import Data.IORef (IORef, newIORef, readIORef, modifyIORef)
import System.Exit (exitWith, ExitCode(ExitSuccess))
import Graphics.Rendering.OpenGL
import Graphics.UI.GLUT

initialWindowSizeX = 800
initialWindowSizeY = 800

pixelsPerUnit = 64

main = do
    getArgsAndInitialize
    initialDisplayMode $= [RGBAMode, DoubleBuffered]
    initialWindowSize $= Size initialWindowSizeX initialWindowSizeY
    (Size screenSizeX screenSizeY) <- get screenSize
    let initialPos = Position x y where
        x = (screenSizeX - initialWindowSizeX) `div` 2
        y = (screenSizeY - initialWindowSizeY) `div` 2
    initialWindowPosition $= initialPos
    createWindow "Grid Raytracing Demo"

    endpoints <- newIORef []

    displayCallback $= display endpoints
    reshapeCallback $= Just reshape
    keyboardMouseCallback $= Just (keyboard endpoints)
    matrixMode $= Projection
    
    mainLoop

display endpoints = do
    clear [ColorBuffer]
    points <- readIORef endpoints
    windowPixelSize <- get windowSize
    drawSquares points
    drawGrid windowPixelSize
    drawLine points
    swapBuffers

reshape size@(Size w h) = do
    viewport $= (Position 0 0, size)
    loadIdentity
    ortho2D 0 ((fromIntegral w) / pixelsPerUnit) 0 ((fromIntegral h) / pixelsPerUnit)

keyboard _ (Char '\27') Down _ _ = exitWith ExitSuccess
keyboard endpoints (MouseButton LeftButton) Down _ pos = do
    (wx, wy) <- worldFromScreen pos
    modifyIORef endpoints (addEndpoint (wx, wy))
    postRedisplay Nothing
keyboard _ _ _ _ _ = return ()

addEndpoint pos (_:_:_) = [pos]
addEndpoint pos endPoints = pos : endPoints

drawSquares (end1 : end0 : _) = do
    currentColor $= Color4 0.25 0.25 0.5 1
    renderPrimitive Quads $ mapM_ drawQuad (squaresOnLine end0 end1)
    where
        drawQuad (x, y) = do
            vertex $ Vertex2 x0 y0
            vertex $ Vertex2 x1 y0
            vertex $ Vertex2 x1 y1
            vertex $ Vertex2 x0 y1
            where
                x0 :: GLint
                x0 = fromIntegral x
                x1 = fromIntegral x + 1
                y0 = fromIntegral y
                y1 = fromIntegral y + 1
drawSquares _ = return ()

drawGrid (Size sizeX sizeY) = do
    currentColor $= Color4 0.5 0.5 0.5 1
    renderPrimitive Lines $ mapM_ (\(x, y) -> vertex $ Vertex2 x y) points
    where
        points = (interleave minYs maxYs) ++ (interleave minXs maxXs)
        minXs = zip (repeat 0) ys
        maxXs = zip (repeat maxX) ys
        minYs = zip xs (repeat 0)
        maxYs = zip xs (repeat maxY)
        xs = take (ceiling maxX) [0..]
        ys = take (ceiling maxY) [0..]
        maxX = (fromIntegral sizeX) / pixelsPerUnit
        maxY = (fromIntegral sizeY) / pixelsPerUnit

interleave (x:xs) (y:ys) = x : y : interleave xs ys
interleave _ _ = []

drawLine ((x0, y0) : (x1, y1) : _) = do
    currentColor $= Color4 1 1 1 1
    renderPrimitive Lines $ do
        vertex $ Vertex2 x0 y0
        vertex $ Vertex2 x1 y1
drawLine ((x0, y0) : _) = do
    currentColor $= Color4 1 1 1 1
    renderPrimitive Points $ do
        vertex $ Vertex2 x0 y0
drawLine _ = return ()

worldFromScreen (Position sx sy) = do
    viewport@(_, Size _ viewSizeY) <- get viewport
    projectionMatrix <- get (matrix $ Just Projection) :: IO (GLmatrix GLdouble)
    modelviewMatrix <- get (matrix $ Just $ Modelview 0) :: IO (GLmatrix GLdouble)
    let screenPos = Vertex3 (fromIntegral sx) (fromIntegral ((viewSizeY - 1) - sy)) 0
    (Vertex3 wx wy wz) <- unProject screenPos projectionMatrix modelviewMatrix viewport
    return (wx, wy)

squaresOnLine (x0, y0) (x1, y1) = take n (genSquares (floor x0) (floor y0) error) where
    n = 1 + abs((floor x1) - (floor x0)) + abs((floor y1) - (floor y0))
    dx = abs (x1 - x0)
    dy = abs (y1 - y0)
    xInc = if x1 > x0 then 1 else -1
    yInc = if y1 > y0 then 1 else -1
    error = dy * tNextX - dx * tNextY
    tNextX
        | x1 > x0   = (fromIntegral (ceiling x0)) - x0
        | otherwise = x0 - (fromIntegral (floor x0))
    tNextY
        | y1 > y0   = (fromIntegral (ceiling y0)) - y0
        | otherwise = y0 - (fromIntegral (floor y0))
    genSquares x y error
        | error > 0 = (x, y) : genSquares x (y + yInc) (error - dx)
        | otherwise = (x, y) : genSquares (x + xInc) y (error + dy)

Both versions use OpenGL with the GLUT library for handling window creation/resizing and keyboard/mouse input. The Haskell version clocks in at around 126 lines, while the C++ version has around 230 lines. I have much more experience writing C++ than Haskell, so the Haskell version may not be great style.

Obvious differences:

Haskell, like Python, uses indentation to indicate the block structure of the code, instead of explicit markers. In C++ I prefer to keep my block markers lined up vertically to aid in matching them visually; this adds lines to the code.

The C++ code contains much more type annotation than the Haskell code. (Haskell type annotations are the phrases beginning with :: such as :: GLint.) Because types are largely inferred, this Haskell program contains only three type annotations. These are needed because Haskell's OpenGL binding is parameterized on input types, so for instance the Vertex2 constructor can take integer, 32-, or 64-bit floating point values as input. Because the input numbers don't offer any clues, the compiler needs a hint to know which one to use. If the constructors for the different types were named differently (Vertex2i, Vertex2f, and Vertex2d, say, as is done in the C API) then these type annotations would not be necessary.

C++ inherits C's single-pass compilation model to some extent, so forward declarations are necessary to be able to call functions before the lines where they are defined. Haskell allows you to define and use functions in any order.

Because OpenGL is implemented as a state machine, and because of the nature of input and output in general, most of the Haskell code has an imperative structure. In fact, the only functions that are purely functional (meaning all inputs and outputs are explicitly listed) are addEndpoint, interleave, and squaresOnLine.

Conditionals tend to be represented without if-statements in the Haskell code. Some if statements transform into multiple function definitions, with each definition covering a portion of the possible input. For instance, drawLine draws nothing if no endpoints are defined; draws a single point if only one endpoint has been defined; and draws a line between endpoints otherwise. Other if statements transform into patterns on a single function definition. For instance, tNextX in squaresOnLine is defined differently depending on the relative magnitudes of x0 and x1.

Loops, likewise, are less obvious in the Haskell version. The first line of squaresOnLine calls take to take the first n items of a list, for instance. It gets an infinite list from the genSquares subfunction which uses recursion to produce the list. Recursion is the heart of how loops are handled in Haskell, but often it's possible to call a function to handle the recursion internally. The mapM_ in drawGrid is an example of this.

Some less obvious differences:

The Haskell program separates concerns better. For instance, the squaresOnLine function generates all the squares intersected by the line segment, but drawing those squares is done in drawSquares. If we need to find the squares to do something else it is easy to reuse squaresOnLine. This could (and probably should) be done in the C++ version, but would likely add a fair amount of additional code.

Internal functions are extremely handy. Note that genSquares inside squaresOnLine can see the variables dx, dy, xInc, and yInc defined in the outer function. This simplifies the interface to genSquares so that it only contains the values that change during generation.

Haskell has no implicit conversion from integer to floating-point types, so lots of fromIntegral calls are needed.

Managing mutable state is kind of a pain. Fortunately for this program, the only state we have is the endpoints that the user has chosen. The program creates a mutable variable (an “IORef”) at startup, modifies it when the user clicks the mouse, and reads it to draw the view.