Voronoi mosaics

While looking for ideas to implement voronoi in sage, I stumbled upon a beautiful paper written by a group of japanese computer graphic professionals from the universities of Hokkaido and Tokyo: A Method for Creating Mosaic Images Using Voronoi Diagrams. The first step of their algorithm is simple yet brilliant: Start with any given image, and superimpose an hexagonal tiling of the plane. By a clever approximation scheme, modify the tiling to become a voronoi diagram that adaptively minimizes some approximation error. As a consequence, the resulting voronoi diagram is somehow adapted to the desired contours of the original image.

(Fig. 1)	(Fig. 2)	(Fig. 3)	(Fig. 4)

In a second step, they manually adjust the Voronoi image interactively by moving, adding, or deleting sites. They also take the liberty of adding visual effects by hand: emphasizing the outlines and color variations in each Voronoi region, so they look like actual pieces of stained glass (Fig. 4).

The following is a simple script in sage to initialize the mosaic. It loads the original image into the array img (Fig. 1), and creates the hexagonal voronoi image newimg (Fig. 2). Each cell on this image has a unique color, and this is obviously chosen so to minimize the error by constants:

from numpy import *
import scipy.ndimage
import scipy.misc

img=matplotlib.image.imread('/Users/blanco/Desktop/flower.png')
img=img[:,:,0:3]
newimg=zeros(img.shape)

# The parameter *scale* modifies the size of the initial hexagonal cells.
scale=sqrt(3.0)*4.0

# Creation of the initial hexagonal tiling.
# Centers are stored in *lattice*, cells in *voronoi*
lattice=[]
A=zeros(newimg[:,:,0].shape)
length=newimg[:,:,0].shape[0]
width=newimg[:,:,0].shape[1]
for k in range(floor(length/scale/sqrt(3)-0.5)+1):
    for j in range(floor(width/scale/3-0.5)+1):
        xcoord=Integer(max(0,min(floor(scale*sqrt(3)*(2*k+1)/2),length-1)))
        ycoord=Integer(max(0,min(floor(3*scale*(2*j+1)/2),width-1)))
        A[xcoord,ycoord]=1
        lattice.append([xcoord,ycoord,1.0])
for k in range(floor(length/scale/sqrt(3)+1)):
    for j in range(floor(width/scale/3+1)):
        xcoord=Integer(max(0,min(floor(scale*k*sqrt(3)),length-1)))
        ycoord=Integer(max(0,min(floor(3*scale*j),width-1)))
        A[xcoord,ycoord]=1
        lattice.append([xcoord,ycoord,1.0])
s=[[0,1,0],[1,1,1],[0,1,0]]
segmentation,segments=scipy.ndimage.label(A,s)
L1,L2=scipy.ndimage.distance_transform_edt(segmentation==0,
                              return_distances=False,
                              return_indices=True)
voronoi=segmentation[L1,L2]

# Resulting image is stored in the array *newimg*
rimg=img[:,:,0]
gimg=img[:,:,1]
bimg=img[:,:,2]
for k in range(1,voronoi.max()+1):
    base=sum(voronoi==k)
    rmean = sum(multiply(voronoi==k,rimg))/base
    gmean = sum(multiply(voronoi==k,gimg))/base
    bmean = sum(multiply(voronoi==k,bimg))/base
    newimg[voronoi==k,0]=rmean
    newimg[voronoi==k,1]=gmean
    newimg[voronoi==k,2]=bmean

In the second step, each cell of the voronoi diagram is allowed to change its center to the eight immediate neighborhood pixels. Among the nine possible cells arising (counting the original as well), we keep the one that offers the smallest error. We repeat this second step as many times as necessary, updating the diagram each time accordingly (see Fig. 3). A naïve script to realize this task in sage could go like this—note that the algorithm is designed with pedagogical purposes in mind, instead of seeking an optimal code that employes the usual trickery to gain speed and minimize memory allocation:

for k in range(len(lattice)):
    center  = lattice[k]
    xcenter = center[0]
    ycenter = center[1]
    ecenter = center[2]
    errors = zeros((3,3))
    for xoffset in range(-1,2):
        for yoffset in range(-1,2):
            testxcenter=xcenter+xoffset
            testycenter=ycenter+yoffset
            if testxcenter>-1
                     and testxcenter<length
                     and testxcenter>-1
                     and testycenter<width:
                testlattice=lattice
                testlattice[k] = [xcenter+xoffset,ycenter+yoffset,ecenter]
                testA=zeros(A.shape)
                for location in testlattice:
                    testA[location[0],location[1]]=1
                    testsegmentation,testsegments=scipy.ndimage.label(testA,s)
                    L1,L2=scipy.ndimage.distance_transform_edt(
                                                  testsegmentation==0,
                                                  return_distances=False,
                                                  return_indices=True)
                    testvoronoi=testsegmentation[L1,L2]
                    newimg=zeros(img.shape)
                    cell = testvoronoi[xcenter,ycenter]
                    base=sum(testvoronoi==cell)
                    rmean=sum(multiply(testvoronoi==cell,rimg))/base
                    gmean=sum(multiply(testvoronoi==cell,gimg))/base
                    bmean=sum(multiply(testvoronoi==cell,bimg))/base
                    newimg[testvoronoi==cell,0]=rmean
                    newimg[testvoronoi==cell,1]=gmean
                    newimg[testvoronoi==cell,2]=bmean
                    error = 0.0
                    for step in range(img.shape[2]):
                        error+=sum(img[testvoronoi==cell,step]
                                    -newimg[testvoronoi==cell,step])^2
                        errors[xoffset+1,yoffset+1] = sqrt(error)
    minerror = errors.min()
    minerrorposition = errors.argmin()
    xshift = minerrorposition//3 -1
    yshift = minerrorposition%3 -1
    newxcenter = xcenter+xshift
    newycenter = ycenter+yshift
    if newxcenter>-1
            and newxcenter<length
            and newycenter>-1
            and newycenter<width
            and minerror<ecenter:
        lattice[k]=[newxcenter,newycenter,minerror]

It is very desirable to design a third step that will actually accomplish the task of imposing the real contours of the image, and enhances them accordingly. In a future post I will present some ideas in this direction, based on two techniques: “constrained Voronoi diagrams”, and “optimal subdivision of cells by linear regression.”