function to_bw(sequence image)
sequence color
integer lum
sequence hist = repeat(0,256)
integer l = 1, r = 256
integer ltot, rtot
    for i=1 to length(image) do
        for j=1 to length(image[i]) do
            color = sq_div(sq_and_bits(image[i][j], {#FF0000,#FF00,#FF}),
                                                    {#010000,#0100,#01})
            lum = floor(0.2126*color[1] + 0.7152*color[2] + 0.0722*color[3])
            image[i][j] = lum
            hist[lum+1] += 1
        end for
    end for
    ltot = hist[l]
    rtot = hist[r]
    while l!=r do
        if ltot<rtot then
            l += 1
            ltot += hist[l]
        else
            r -= 1
            rtot += hist[r]
        end if
    end while
    lum = l
    for i=1 to length(image) do
        for j=1 to length(image[i]) do
            image[i][j] = iff(image[i][j]<lum?black:white)
        end for
    end for
    return image
end function

sequence img = read_ppm("Lena.ppm")
    img = to_bw(img)
    write_ppm("LenaBW.ppm",img)