Making Circles in Python and Julia

circles.jl

import ArgParse
import ImageSegmentation
import Images
import Luxor
import Statistics


int_round(f) = round(Int, f)


struct Circle{T}
    x::Int
    y::Int
    rgb::T
    r::Float64
end

# Is there a more normal way to do a keyword constructor for a type?
function Circle(;x::Int, y::Int, rgb::Tuple{Int64, Int64, Int64}, r::Float64)
    return Circle(x, y, rgb, r)
end


function parse_commandline()
    s = ArgParse.ArgParseSettings()
    s.exc_handler = ArgParse.debug_handler

    @ArgParse.add_arg_table s begin
        "--image", "-i"
            help = "Path to input image"
            arg_type = String
            required = true
        "--min_radius", "-r"
            help = "Minimum radius of generated circles"
            arg_type = Float64
            default = 2.0
    end

    return ArgParse.parse_args(s)
end


function scale(img, target_px)
    h, w = size(img)
    # Need sqrt because ratio gets applied in both x and y dimension.
    return Images.imresize(img, ratio=(target_px / (h*w))^.5)
end


function segment(img)
    # k controls segment size. 5 - 500 is a good range
    # min_size gets rid of small segments 
    return ImageSegmentation.felzenszwalb(img, 60, 100)
end


function make_mask(label_ix, imgh, imgw)
    # Create a mask array that's just large enough to contain
    # the segment we are working on. This will allow the distance transform
    # to operate the smallest array possible, for speed.
    pad = Dict(true => 0, false => 1)
    (top, left), (bottom, right) = map(ix -> ix.I, extrema(label_ix))

    # Allocate a pixel-wide boundary strip if we aren't at the image border.
    top_pad, bottom_pad = pad[top == 1], pad[bottom == imgh]
    left_pad, right_pad = pad[left == 1], pad[right == imgw]

    mask = ones(Bool,
        bottom - top + 1 + top_pad + bottom_pad,
        right - left + 1 + left_pad + right_pad
    )

    new_ix = map(x -> x - CartesianIndex(top - top_pad - 1, left - left_pad - 1), label_ix)
    mask[new_ix] .= 0

    return (left - left_pad, top - top_pad, mask)
end


function circle_indices(center, radius, height, width)
    # Code inspired by:
    # https://github.com/JuliaImages/ImageDraw.jl/blob/master/src/ellipse2d.jl

    indices = CartesianIndex{2}[]
    r = CartesianIndex(int_round(radius), int_round(radius))
    y, x = center.I
    # Examine (2*radius + 1)^2 spots, expecting array length ~ pi * radius^2.
    for ix in center - r: center + r
        row, col = ix.I
        if 1 <= row <= height && 1 <= col <= width
            if (row - y)^2 + (col - x)^2 <= radius^2
                push!(indices, ix)
            end
        end
    end

    return indices
end


function make_circles(img, seg, min_radius)
    float2int(f) = round(Int, 256*f)
    to_clr(l) = float2int.( (l.r, l.g, l.b) )

    h, w = size(img)

    # Add a big circle with the avg color for the background.
    circles = [Circle(
        x=int_round(w/2),
        y=int_round(h/2),
        rgb= to_clr(Statistics.mean(img)),
        r=(w*w + h*h)^.5 + 2
    )]

    for label in seg.segment_labels
        m = findall(x -> x == label, seg.image_indexmap)

        left, top, mask = make_mask(m, h, w)

        while true
            f = Images.feature_transform(mask)
            edt = Images.distance_transform(f)
            ix = argmax(edt)
            r = edt[ix]

            if r < min_radius || r == Inf
                break
            end

            circle_ix = circle_indices(ix, r, size(edt)...)
            mask[circle_ix] .= 1

            offset_ix = map(x -> x + CartesianIndex(top - 1, left - 1), circle_ix)
            clr = to_clr(Statistics.mean(img[offset_ix]))

            y, x = ix.I

            push!(circles, Circle(
                x=left + x,
                y=top + y,
                rgb=clr,
                r=r
            ))
        end
    end

    return circles
end


function draw(img, circles, filename)
    h, w = size(img)
    Luxor.Drawing(w, h, filename)
    # Paint from biggest to smallest
    for c in sort(circles, by=c -> c.r, rev=true)
        r, g, b = map(x -> x / 255, c.rgb)
        Luxor.setcolor(r, g, b)
        Luxor.circle(c.x, c.y, c.r - .5, :fill)
    end
    Luxor.finish()
end


function main()
    args = nothing
    try
        args = parse_commandline()
        println("args", args)
    catch err
        args = Dict("image" => "test.jpg", "min_radius" => 2)
        println("No args supplied, using defaults: ", args)
    end

    path = args["image"]
    name, ext = splitext(basename(path))


    println("Loading $(path) ...")
    @time img = Images.load(path)

    println("Resizing ...")
    @time img = scale(img, 1e6)  # Downsample to ~1M pixels
    Images.save("./$(name)-orig.png", img)

    println("Segmenting ... ")
    @time seg = segment(img)
    Images.save("./$(name)-seg.png", map(i-> seg.segment_means[i], seg.image_indexmap))

    println("Circling ...")
    @time circles = make_circles(img, seg, args["min_radius"])
    println("Made $(length(circles)) circles")

    println("Drawing ...")
    @time draw(img, circles, "./$(name)-circled.png")
end


if !isinteractive()
    main()
end

circles.py

from scipy import ndimage
from skimage import color, io, segmentation, transform
from skimage.draw import disk, circle_perimeter_aa, set_color
from PIL import Image

import argparse
import json
import numpy as np
import os
import time


int_round = lambda c: int(round(c))


def timed(function):
    def wrap(*args, **kwargs):
        start_time = time.time()
        result = function(*args, **kwargs)
        end_time  = time.time()
        duration = (end_time- start_time)*1000.0
        f_name = function.__name__
        print("{} took {:.3f} ms".format(f_name, duration))

        return result
    return wrap


@timed
def segment(img):
    return segmentation.felzenszwalb(img, scale=60, min_size=100)


def amax(a):
    """returns the max of an array, and the coordinates at which it occurs."""
    row, col = np.unravel_index(np.argmax(a), a.shape)
    return a[row][col], (row, col)


@timed
def scale(img, target_px):
    """Scales an input image to have roughly target_px pixels."""
    w, h, d = img.shape
    # Need sqrt below because scale gets applied in both x and y dimension.
    scale_factor = (float(target_px) / (w*h)) ** .5

    nw, nh = int(w * scale_factor), int(h * scale_factor)
    return transform.resize(img, (nw, nh))


def distance_transform(a, ft, ix):
    """
        a is a binary image of dtype np.int8
        ft = np.zeros((input.ndim,) + input.shape, dtype=np.int32)
        ix = np.indices(input.shape, dtype=ft.dtype)
    """
    ndimage._nd_image.euclidean_feature_transform(a, None, ft)

    # calculate the distance transform
    dt = ft - ix
    np.multiply(dt, dt, dt) # dt *= dt
    return np.sqrt(dt[0] + dt[1])


def make_mask(labeled_px, label):
    """ Create a mask array that's just large enough to contain
        the segment we are working on. This will make the distance transform
        (the slowest part of this code) much more efficient.
    """
    pad = {True:1, False:0}

    rows, cols = np.where(labeled_px == label)
    top, bottom = rows.min(), rows.max()
    left, right = cols.min(), cols.max()

    # We need to allocate a pixel-wide boundary strip around this segment
    # unless we're flush with one of the image sides
    h, w = labeled_px.shape
    top_pad, bottom_pad = pad[top != 0], pad[bottom != (h-1)]
    left_pad, right_pad = pad[left != 0], pad[right != (w-1)]

    mask = np.zeros((
        bottom - top + 1 + top_pad + bottom_pad,
        right - left + 1 + left_pad + right_pad
    ), dtype='int8')

    mask[rows - top + top_pad, cols - left + left_pad] = 1

    return (left - left_pad, top - top_pad, mask)


@timed
def make_circles(img, labeled_px, min_radius):
    # Add a big circle with the avg color for the background.
    h, w, _ = img.shape
    circles = [dict(
        x=int_round(w/2.0),
        y=int_round(h/2.0),
        rgb=list(map(int_round, 255*np.average(img, (0, 1)))),
        r=int_round(2 + (w*w + h*h)**.5) # 2 to account for rounding, being off-center.
    )]

    # OK, now make some circles
    labels = sorted(set(labeled_px.flatten()))
    for label in labels:
        left, top, mask = make_mask(labeled_px, label)
        # preallocate arrays, to avoid a bunch of reallocs
        indices = np.zeros((mask.ndim,) + mask.shape, dtype='int32')
        ix = np.indices(mask.shape, dtype='float64')

        # 1) Compute Euclidean Distance Transform (edt)
        # 2) For each max in the edt, draw a circle, then zero it out.
        # 3) repeat until the edt doesn't have entries larger than min_radius
        while True:
            edt = distance_transform(mask, indices, ix)
            r, (row, col) = amax(edt)

            if r < min_radius:
                break

            cr, cc = disk((row, col), r, shape=mask.shape)

            # Mask out the circle
            mask[cr, cc] = 0

            # Compute average color. Does this need to be done here?
            clr = list(map(int_round, 255 * np.average(img[
                top  + cr,
                left + cc
            ], 0)))

            # Store painting instructions
            circles.append(dict(x=int(left + col), y=int(top + row), r=r, rgb=clr))

    return circles


@timed
def draw(img, circles, filename):
    avg = np.array(list(map(int_round, np.average(img, (0, 1)))), dtype='uint8')
    out = avg * np.ones(img.shape, dtype='uint8')
    for c in sorted(circles, key=lambda d: d['r'], reverse=True):
        x, y, r = c['x'], c['y'], int_round(c['r'] - 1)
        clr = c['rgb']

        # Compute the circle's border, antialiased, then the interior
        rr, cc, val = circle_perimeter_aa(y, x, r, out.shape)
        set_color(out, (rr, cc), clr, val)

        # Interior
        out[disk((y, x), r, shape=out.shape)] = clr

    io.imsave(filename, out)


@timed
def main(filename, min_radius=3):
    def imsave(s, img):
        # save a float image without annoying warnings
        io.imsave(s, (255 * img).astype(np.uint8))

    print("Loading %s ..." % filename)
    img = io.imread(filename)
    name, ext = os.path.splitext(os.path.basename(filename))

    print("Resizing ...")
    img = scale(img, 1e6)
    imsave(filename.replace(ext, '-orig.png'), img)

    print("Segmenting ...")
    seg = segment(img)
    labeled_px = color.label2rgb(seg, img, kind='avg', bg_label=-1)
    imsave(filename.replace(ext, '-seg.png'), labeled_px)

    print("Circling ...")
    circles = make_circles(img, seg, min_radius)

    print("Drawing ...")
    draw(img, circles, filename.replace(ext, '-circled.png'))

    print("Saving ...")
    with open('./%s-circles.json' % name, 'w') as f:
        json.dump(circles, f, indent=4)


if __name__ == '__main__':
    ap = argparse.ArgumentParser()
    ap.add_argument("-i", "--image", required=False, default="./test.jpg",
        help="Path to input image")
    ap.add_argument("-r", "--min_radius", required=False, default=2, type=float,
        help="Minimum radius of generated circles")

    args = vars(ap.parse_args())
    main(args['image'], args['min_radius'])

Here is the test image, test.jpg:

Total runtime is 2:41.73, broken down as follows:

 $ time julia circles.jl
No args supplied, using defaults: Dict{String,Any}("min_radius" => 2,"image" => "test.jpg")
Loading test.jpg ...
  1.659704 seconds (4.87 M allocations: 276.811 MiB, 6.11% gc time)
Resizing ...
  0.463902 seconds (1.48 M allocations: 93.129 MiB, 3.75% gc time)
Segmenting ...
  1.178172 seconds (2.64 M allocations: 242.965 MiB, 11.19% gc time)
Circling ...
143.547013 seconds (1.95 M allocations: 66.293 GiB, 3.88% gc time)
Made 7843 circles
Drawing ...
  0.683826 seconds (189.30 k allocations: 12.377 MiB)
julia circles.jl  156.09s user 5.58s system 99% cpu 2:41.73 total

Python timing is 17.7 seconds faster on the same machine:

$ time python circles.py
Loading ./test.jpg ...
Resizing ...
scale took 350.042 ms
Segmenting ...
segment took 1878.778 ms
Circling ...
make_circles took 135736.593 ms
Drawing ...
draw took 1844.562 ms
Saving ...
main took 142669.479 ms
python circles.py  138.25s user 5.17s system 99% cpu 2:24.03 total

posted about this on the Julia discourse

Output should look like this:

Nice script!