Goose

22_image_filters.goose

// Image kernels over row slices. A grayscale image is generated, blurred,
// and run through a Sobel edge detector, and the result is written as a
// PGM file. The kernels take the three source rows and the destination row
// as slices whose lengths are asserted once per row, so inside the pixel
// loop every index is a loop counter plus a constant against a known
// length: the compiler proves the pixel loops' bounds checks away and the
// C backend sees the same loop it would see from C.
//
// Build and run (from this directory):
//   goose -o image_filters.c 22_image_filters.goose && cl image_filters.c && image_filters

import std;
import math;
import os;

let W = 512;
let H = 384;

// A test picture: a gradient with rings, a bright disc and some noise.
fn generate(img: u8[>..]&) {
    var r = rng { 99 };
    for y in H {
        for x in W {
            let dx = (x - W / 2) as f64;
            let dy = (y - H / 2) as f64;
            let d = sqrt(dx * dx + dy * dy);
            var v = 96.0 + 64.0 * sin(d / 12.0) + (x as f64) / W * 60.0;
            if d < 60.0 { v = 230.0; }
            if (x / 64 + y / 64) % 2 == 0 && y > H * 3 / 4 { v = 40.0; }
            v += r.rand_int(21) - 10;
            img.push(clamp(v, 0.0, 255.0) as u8);
        }
    }
}

// A 3x3 box blur of the interior; the border is copied. The rows are
// sliced once per output row, and the inner loop indexes those.
fn blur(src: u8[:], dst: u8[:]) {
    assert(src.len == W * H && dst.len == W * H);
    copy_into(dst[..W], src[..W]);
    copy_into(dst[(H - 1) * W..], src[(H - 1) * W..]);
    for y in 1..H - 1 {
        let a = src[(y - 1) * W..y * W];
        let b = src[y * W..(y + 1) * W];
        let c = src[(y + 1) * W..(y + 2) * W];
        let d = dst[y * W..(y + 1) * W];
        assert(a.len == W && b.len == W && c.len == W && d.len == W);
        d[0] = b[0];
        d[W - 1] = b[W - 1];
        for x in 1..W - 1 {
            let s = a[x - 1] as u16 + a[x] + a[x + 1]
                  + b[x - 1] + b[x] + b[x + 1]
                  + c[x - 1] + c[x] + c[x + 1];
            d[x] = (s / 9) as! u8;
        }
    }
}

// Sobel edge magnitude, same shape: two 3x3 kernels, one per axis.
fn sobel(src: u8[:], dst: u8[:]) {
    assert(src.len == W * H && dst.len == W * H);
    fill(dst[..W], 0);
    fill(dst[(H - 1) * W..], 0);
    for y in 1..H - 1 {
        let a = src[(y - 1) * W..y * W];
        let b = src[y * W..(y + 1) * W];
        let c = src[(y + 1) * W..(y + 2) * W];
        let d = dst[y * W..(y + 1) * W];
        assert(a.len == W && b.len == W && c.len == W && d.len == W);
        d[0] = 0;
        d[W - 1] = 0;
        for x in 1..W - 1 {
            let gx = (a[x + 1] as i32 - a[x - 1]) + 2 * (b[x + 1] as i32 - b[x - 1])
                   + (c[x + 1] as i32 - c[x - 1]);
            let gy = (c[x - 1] as i32 - a[x - 1]) + 2 * (c[x] as i32 - a[x])
                   + (c[x + 1] as i32 - a[x + 1]);
            let g = abs(gx) + abs(gy);
            d[x] = if g > 255 { 255 } else { g as! u8 };
        }
    }
}

fn histogram(img: u8[:]) -> i64[8] {
    var h: i64[8] = [0; 8];
    for v in img { h[v / 32] += 1; }
    return h;
}

fn write_pgm(path: u8[:], img: u8[:]) {
    var out: u8[>..] = [];
    format(out, "P5 ", W, " ", H, " 255\n");
    out.append(img);
    print(path, ": ", write_file(path, out), " (", out.len, " bytes)");
}

fn main() {
    var img: u8[>..] = [];
    generate(img);
    var tmp: u8[>..] = [];
    tmp.resize(W * H, 0);
    var edges: u8[>..] = [];
    edges.resize(W * H, 0);

    let t0 = clock();
    for i in 4 {                         // four blur passes, ping-ponging the buffers
        blur(img, tmp);
        blur(tmp, img);
    }
    let t1 = clock();
    sobel(img, edges);
    let t2 = clock();

    var note: u8[>..] = [];
    format(note, "8 blur passes: ");
    format_flt(note, (t1 - t0) * 1000.0, 2);
    format(note, " ms, sobel: ");
    format_flt(note, (t2 - t1) * 1000.0, 2);
    format(note, " ms\n");
    write_stderr(note);

    print("histogram of the blurred image: ", histogram(img));
    print("histogram of the edges: ", histogram(edges));
    var strong = 0;
    for v in edges { if v > 128 { strong++; } }
    print(strong, " strong edge pixels of ", W * H);
    write_pgm("blurred.pgm", img);
    write_pgm("edges.pgm", edges);
}