Contouring algorithm in Julia.

contour.jl

# contour(A, v)
# Implements the Marching Squares contouring algorithm (see:
# https://en.wikipedia.org/wiki/Marching_squares for details). The
# matrix `A` is an m x n scalar field and the scalar `v` is the
# value to be contoured. The return value is an (m-1) x (n-1) matrix
# of contour line type indices (see the Wikipedia article).
# TODO: Implement in parallel
function contour(A, v)
    rows, cols = size(A)

    B = Array(Int64, rows, cols)
    for i = 1:rows
        for j = 1:cols
            B[i,j] = (A[i,j] > v) ? 0 : 1
        end
    end

    C = Array(Int64, rows-1, cols-1)
    for i = 2:rows
        for j = 2:cols
            acc = 0
            acc = (acc | B[i-1,j-1]) << 1
            acc = (acc | B[i-1,j]) << 1
            acc = (acc | B[i,j]) << 1
            acc = acc | B[i,j-1]
            C[i-1,j-1] = acc
        end
    end

    C
end