While looking for ideas to implement voronoi in sage, we 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 $latex L_2$ 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 width=newimg[:,:,0].shape 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 ycenter = center ecenter = center 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,location]=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): 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 we will present some ideas in this direction, based on two techniques: "constrained Voronoi diagrams", and "optimal subdivision of cells by linear regression."