[Question]
  [
I've compiled a mosaic of 2025 headshots from the avatars of the [top Stack Overflow users](https://stackoverflow.com/users?page=1&tab=reputation&filter=all).  
*(Click the image to view it full-size.)*
[![StackOverflow headshots mosaic](https://i.stack.imgur.com/GY1es.jpg)](http://www.pictureshack.us/images/25745_avatars.png)
Your task is to write an algorithm that will create an accurate photomosaic of another image using the 48×48 pixel avatars from this 45×45 grid of them.
## Test Images
Here are the test images. The first is, of course, a light bulb!  
*(They are not full-sized here. Click an image to view it at full-size. Half-sized versions are available for [The Kiss](http://www.pictureshack.us/images/65020_kiss_half.jpg), [A Sunday Afternoon...](http://www.pictureshack.us/images/75361_sunday_half.jpg), [Steve Jobs](http://www.pictureshack.us/images/73806_steve_half.jpg), and [the spheres](http://www.pictureshack.us/images/99408_spheres_half.png).)*
[![lightbulb](https://i.stack.imgur.com/GysVo.png)](http://www.pictureshack.us/images/77388_lightbulb.png)
[![The Kiss](https://i.stack.imgur.com/Ts649.jpg)](http://www.pictureshack.us/images/61555_kiss.jpg)
[![A Sunday Afternoon on the Island of La Grande Jatte](https://i.stack.imgur.com/11PKc.jpg)](http://www.pictureshack.us/images/45148_sunday.jpg)
[![Steve Jobs](https://i.stack.imgur.com/LsL3f.jpg)](http://www.pictureshack.us/images/47824_steve.jpg)
[![spheres](https://i.stack.imgur.com/wLzFs.png)](http://www.pictureshack.us/images/21104_spheres.png)
Thanks to Wikipedia for all but the raytraced spheres.
At full-size these images all have dimensions divisible by 48. The larger ones had to be JPEGs so they could be compressed enough to upload.
## Scoring
This is a popularity contest. The submission with mosaics that most accurately depict the original images should be voted up. I will accept the highest voted answer in a week or two.
## Rules
* Your photomosaics must be entirely composed of **unaltered** 48×48 pixel avatars taken from the mosaic above, arranged in a grid.
* You may reuse an avatar in a mosaic. (Indeed for the larger test images you'll have to.)
* Show your output, but keep in mind that the test images are very large, and **StackExchange only allows posting of images up to 2MB**. So compress your images or host them somewhere else and put smaller versions here.
* To be confirmed the winner you must provide PNG versions of your light bulb or spheres mosaics. This is so I can validate them (see below) to ensure you are not adding extra colors to the avatars to make the mosaics look better.
## Validator
This Python script can be used to check whether a completed mosaic really uses unaltered avatars. Just set `toValidate` and `allTiles`. It is unlikely to work for JPEGs or other lossy formats since it compares things exactly, pixel-for-pixel.
```
from PIL import Image, ImageChops
toValidate = 'test.png' #test.png is the mosaic to validate
allTiles = 'avatars.png' #avatars.png is the grid of 2025 48x48 avatars
def equal(img1, img2):
    return ImageChops.difference(img1, img2).getbbox() is None
def getTiles(mosaic, (w, h)):
    tiles = {}
    for i in range(mosaic.size[0] / w):
        for j in range(mosaic.size[1] / h):
            x, y = i * w, j * h
            tiles[(i, j)] = mosaic.crop((x, y, x + w, y + h))
    return tiles
def validateMosaic(mosaic, allTiles, tileSize):
    w, h = tileSize
    if mosaic.size[0] % w != 0 or mosaic.size[1] % h != 0:
        print 'Tiles do not fit mosaic.'
    elif allTiles.size[0] % w != 0 or allTiles.size[1] % h != 0:
        print 'Tiles do not fit allTiles.'
    else:
        pool = getTiles(allTiles, tileSize)
        tiles = getTiles(mosaic, tileSize)
        matches = lambda tile: equal(tiles[pos], tile)
        success = True
        for pos in tiles:
            if not any(map(matches, pool.values())):
                print 'Tile in row %s, column %s was not found in allTiles.' % (pos[1] + 1, pos[0] + 1)
                success = False
        if success:
            print 'Mosaic is valid.'
            return
    print 'MOSAIC IS INVALID!'
validateMosaic(Image.open(toValidate).convert('RGB'), Image.open(allTiles).convert('RGB'), (48, 48))
```
Good luck everyone! I cant wait to see the results.
**Note:** I know photomosaic algorithms are easy to find online, but they aren't yet on this site. I'm really hoping we see something more interesting than the usual *"average each tile and each grid space and match them up"* algorithm.
  
[Answer]
# Mathematica, with control for granularity
This uses the 48 x 48 pixel photos, as required. By default, it will swap those pixels for a corresponding 48x48 pixel square from the image to be approximated.
However, the size of the destination squares can be set to be smaller than 48 x 48, allowing for greater fidelity to detail. (see the examples below).
## Preprocessing the palette
`collage` is the image containing the photos to serve as the palette.
`picsColors`is a list of individual photos paired with their mean red, mean green, and mean blue values.
targetColorToPhoto[]` takes the average color of the target swath and finds the photo from the palette that best matches it.
```
parts=Flatten[ImagePartition[collage,48],1];
picsColors={#,c=Mean[Flatten[ImageData[#],1]]}&/@parts;
nf=Nearest[picsColors[[All,2]]];
targetColorToPhoto[p_]:=Cases[picsColors,{pic_,nf[p,1][[1]]}:>pic][[1]]
```
**Example**
Let's find the photo that best matches RGBColor[0.640, 0.134, 0.249]:
![example1](https://i.stack.imgur.com/QSap7.png)
---
## photoMosaic
```
photoMosaic[rawPic_, targetSwathSize_: 48] :=
 Module[{targetPic, data, dims, tiles, tileReplacements, gallery},
  targetPic = Image[data = ImageData[rawPic] /. {r_, g_, b_, _} :> {r, g, b}];
  dims = Dimensions[data];
  tiles = ImagePartition[targetPic, targetSwathSize];
  tileReplacements = targetColorToPhoto /@ (Mean[Flatten[ImageData[#], 1]] & /@Flatten[tiles, 1]);
  gallery = ArrayReshape[tileReplacements, {dims[[1]]/targetSwathSize,dims[[2]]/targetSwathSize}];
  ImageAssemble[gallery]]
```
`photoMosaic takes as input the raw picture we will make a photo mosaic of.
`targetPic` will remove a fourth parameter (of PNGs and some JPGs), leaving only R, G, B.
`dims` are the dimensions of `targetPic`.
`tiles` are the little squares that together comprise the target picture.
`targetSwathSize is the granularity parameter; it defaults at 48 (x48).`
`tileReplacements` are the photos that match each tile, in proper order.
`gallery` is the set of tile-replacements (photos) with the proper dimensionality (i.e. the number of rows and columns that match the tiles).
`ImageAssembly` joins up the mosaic into a continuous output image.
---
**Examples**
This replaces each 12x12 square from the image, Sunday, with a corresponding 48 x 48 pixel photograph that best matches it for average color.
```
photoMosaic[sunday, 12]
```
![sunday2](https://i.stack.imgur.com/a2fn8.jpg)
---
Sunday (detail)
![top hat](https://i.stack.imgur.com/GOW6T.jpg)
---
```
photoMosaic[lightbulb, 6]
```
![lightbulb 6](https://i.stack.imgur.com/k3K4t.jpg)
---
```
photoMosaic[stevejobs, 24]
```
![steve jobs 24](https://i.stack.imgur.com/eweC7.jpg)
---
Detail, stevejobs.
![jobs detail](https://i.stack.imgur.com/U7Uf7.jpg)
---
```
photoMosaic[kiss, 24]
```
![kiss](https://i.stack.imgur.com/acPZS.jpg)
---
Detail of the Kiss:
![detail kiss](https://i.stack.imgur.com/GaSBA.png)
---
```
photoMosaic[spheres, 24]
```
![spheres](https://i.stack.imgur.com/9v8FN.jpg)
[Answer]
### Java, average distance
```
package photomosaic;
import java.awt.image.BufferedImage;
import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import javax.imageio.ImageIO;
public class MainClass {
    static final String FILE_IMAGE = "45148_sunday.jpg";
    static final String FILE_AVATARS = "25745_avatars.png";
    static final String FILE_RESULT = "mosaic.png";
    static final int BLOCKSIZE = 48;
    static final int SCALING = 4;
    static final int RADIUS = 3;
    public static void main(String[] args) throws IOException {
        BufferedImage imageAvatars = ImageIO.read(new File(FILE_AVATARS));
        int[] avatars = deBlock(imageAvatars, BLOCKSIZE);
        BufferedImage image = ImageIO.read(new File(FILE_IMAGE));
        int[] data = deBlock(image, BLOCKSIZE);
        // perform the mosaic search on a downscaled version
        int[] avatarsScaled = scaleDown(avatars, BLOCKSIZE, SCALING);
        int[] dataScaled = scaleDown(data, BLOCKSIZE, SCALING);
        int[] bests = mosaicize(dataScaled, avatarsScaled, (BLOCKSIZE / SCALING) * (BLOCKSIZE / SCALING), image.getWidth() / BLOCKSIZE);
        // rebuild the image from the mosaic tiles
        reBuild(bests, data, avatars, BLOCKSIZE);
        reBlock(image, data, BLOCKSIZE);
        ImageIO.write(image, "png", new File(FILE_RESULT));
    }
    // a simple downscale function using simple averaging
    private static int[] scaleDown(int[] data, int size, int scale) {
        int newsize = size / scale;
        int newpixels = newsize * newsize;
        int[] result = new int[data.length / scale / scale];
        for (int r = 0; r < result.length; r += newpixels) {
            for (int y = 0; y < newsize; y++) {
                for (int x = 0; x < newsize; x++) {
                    int avgR = 0;
                    int avgG = 0;
                    int avgB = 0;
                    for (int sy = 0; sy < scale; sy++) {
                        for (int sx = 0; sx < scale; sx++) {
                            int dt = data[r * scale * scale + (y * scale + sy) * size + (x * scale + sx)];
                            avgR += (dt & 0xFF0000) >> 16;
                            avgG += (dt & 0xFF00) >> 8;
                            avgB += (dt & 0xFF) >> 0;
                        }
                    }
                    avgR /= scale * scale;
                    avgG /= scale * scale;
                    avgB /= scale * scale;
                    result[r + y * newsize + x] = 0xFF000000 + (avgR << 16) + (avgG << 8) + (avgB << 0);
                }
            }
        }
        return result;
    }
    // the mosaicize algorithm: take the avatar with least pixel-wise distance
    private static int[] mosaicize(int[] data, int[] avatars, int pixels, int tilesPerLine) {
        int tiles = data.length / pixels;
        // use random order for tile search
        List<Integer> ts = new ArrayList<Integer>();
        for (int t = 0; t < tiles; t++) {
            ts.add(t);
        }
        Collections.shuffle(ts);
        // result array
        int[] bests = new int[tiles];
        Arrays.fill(bests, -1);
        // make list of offsets to be ignored
        List<Integer> ignores = new ArrayList<Integer>();
        for (int sy = -RADIUS; sy <= RADIUS; sy++) {
            for (int sx = -RADIUS; sx <= RADIUS; sx++) {
                if (sx * sx + sy * sy <= RADIUS * RADIUS) {
                    ignores.add(sy * tilesPerLine + sx);
                }
            }
        }
        for (int t : ts) {
            int b = t * pixels;
            int bestsum = Integer.MAX_VALUE;
            for (int at = 0; at < avatars.length / pixels; at++) {
                int a = at * pixels;
                int sum = 0;
                for (int i = 0; i < pixels; i++) {
                    int r1 = (avatars[a + i] & 0xFF0000) >> 16;
                    int g1 = (avatars[a + i] & 0xFF00) >> 8;
                    int b1 = (avatars[a + i] & 0xFF) >> 0;
                    int r2 = (data[b + i] & 0xFF0000) >> 16;
                    int g2 = (data[b + i] & 0xFF00) >> 8;
                    int b2 = (data[b + i] & 0xFF) >> 0;
                    int dr = (r1 - r2) * 30;
                    int dg = (g1 - g2) * 59;
                    int db = (b1 - b2) * 11;
                    sum += Math.sqrt(dr * dr + dg * dg + db * db);
                }
                if (sum < bestsum) {
                    boolean ignore = false;
                    for (int io : ignores) {
                        if (t + io >= 0 && t + io < bests.length && bests[t + io] == at) {
                            ignore = true;
                            break;
                        }
                    }
                    if (!ignore) {
                        bestsum = sum;
                        bests[t] = at;
                    }
                }
            }
        }
        return bests;
    }
    // build image from avatar tiles
    private static void reBuild(int[] bests, int[] data, int[] avatars, int size) {
        for (int r = 0; r < bests.length; r++) {
            System.arraycopy(avatars, bests[r] * size * size, data, r * size * size, size * size);
        }
    }
    // splits the image into blocks and concatenates all the blocks
    private static int[] deBlock(BufferedImage image, int size) {
        int[] result = new int[image.getWidth() * image.getHeight()];
        int r = 0;
        for (int fy = 0; fy < image.getHeight(); fy += size) {
            for (int fx = 0; fx < image.getWidth(); fx += size) {
                for (int l = 0; l < size; l++) {
                    image.getRGB(fx, fy + l, size, 1, result, r * size * size + l * size, size);
                }
                r++;
            }
        }
        return result;
    }
    // unsplits the block version into the original image format
    private static void reBlock(BufferedImage image, int[] data, int size) {
        int r = 0;
        for (int fy = 0; fy < image.getHeight(); fy += size) {
            for (int fx = 0; fx < image.getWidth(); fx += size) {
                for (int l = 0; l < size; l++) {
                    image.setRGB(fx, fy + l, size, 1, data, r * size * size + l * size, size);
                }
                r++;
            }
        }
    }
}
```
The algorithm performs a search through all avatar tiles for each grid space separately.
Due to the small sizes I didn't implement any sophisticated data structurs or search algorithms but simply brute-force the whole space.
This code does not do any modifications to the tiles (e.g. no adaption to destination colors).
### Results
*Click for full-size image.*
[![light buld](https://i.stack.imgur.com/xwoxk.png)](http://www.pictureshack.us/images/42778_mosaic_04_lightbulb.png)
[![spheres](https://i.stack.imgur.com/b4Psv.png)](http://www.pictureshack.us/images/57865_mosaic_04_spheres.png)  
[![sunday](https://i.stack.imgur.com/ftTVx.png)](http://www.pictureshack.us/images/48448_mosaic_04_sunday.png)
### Effect of radius
Using `radius` you can reduce the repetitiveness of the tiles in the result. Setting `radius=0` there is no effect. E.g. `radius=3` supresses the same tile within a radius of 3 tiles.
[![light buld](https://i.stack.imgur.com/xwoxk.png)](http://www.pictureshack.us/images/42778_mosaic_04_lightbulb.png)
[![sunday](https://i.stack.imgur.com/ftTVx.png)](http://www.pictureshack.us/images/48448_mosaic_04_sunday.png)
radius=0
[![light buld](https://i.stack.imgur.com/sNhls.png)](http://www.pictureshack.us/images/82587_mosaic_04_r03_lightbulb.png)  
[![light buld](https://i.stack.imgur.com/WxxH2.png)](http://www.pictureshack.us/images/22259_mosaic_04_r03_sunday.png)  
### Effect of scaling factor
Using the `scaling` factor we can determine how the matching tile is searched. `scaling=1` means searching for a pixel-perfect match while `scaling=48` does a average-tile search.
[![scaling 48](https://i.stack.imgur.com/UtHKr.png)](http://www.pictureshack.us/images/45469_mosaic_48_sunday.png)  
scaling=48
[![scaling 16](https://i.stack.imgur.com/7ao0C.png)](http://www.pictureshack.us/images/72507_mosaic_16_sunday.png)  
scaling=16
[![scaling 4](https://i.stack.imgur.com/kk4oP.png)](http://www.pictureshack.us/images/48448_mosaic_04_sunday.png)  
scaling=4
[![scaling 1](https://i.stack.imgur.com/5bN7t.png)](http://www.pictureshack.us/images/10567_mosaic_01_sunday.png)  
scaling=1
[Answer]
**JS**
Same as in previous golf: <http://jsfiddle.net/eithe/J7jEk/> :D
(this time called with `unique: false, {pixel_2: {width: 48, height: 48}, pixel_1: {width: 48, height: 48}}`) (don't treat palette to use one pixel once, palette pixels are 48x48 swatches, shape pixels are 48x48 swatches).
Currently it searches through the avatars list to find the nearest matching by weight of selected algorithm, however it doesn't perform any color uniformity matching (something I need to take a look at.
* [![balanced](https://i.stack.imgur.com/1RFS0.png)](https://i.stack.imgur.com/Mt8dZ.png)
* [![lab](https://i.stack.imgur.com/Rr3Tp.png)](https://i.stack.imgur.com/B4J1z.jpg)
Unfortunately I'm not able to play around with larger images, because my RAM runs out :D If possible I'd appreciate smaller output images. If using 1/2 of image size provided, here's Sunday Afternoon:
* ![balanced](https://i.stack.imgur.com/V1R9w.png)
* ![lab](https://i.stack.imgur.com/594iG.png)
[Answer]
## GLSL
The difference between this challenge and the one at [American Gothic in the palette of Mona Lisa: Rearrange the pixels](https://codegolf.stackexchange.com/questions/33172/american-gothic-in-the-palette-of-mona-lisa-rearrange-the-pixels) interested me, because the mosaic tiles can be re-used, whereas the pixels couldn't. This means that it is possible to easily parallelize the algorithm, so I decided to attempt a massively parallel version. By "massively" I mean using the 1344 shader cores on my desktop's GTX670 all simultaneously, via GLSL.
## Method
The actual tile-matching is simple: I calculate the RGB distance between each pixel in a target area and the mosaic tile area, and choose the tile with the lowest difference (weighted by brightness values). The tile index is written out in the red and green color attributes of the fragment, then after all the fragments have been rendered I read the values back out of the framebuffer and build the output image from those indices. The actual implementation is quite a hack; instead of creating an FBO I just opened a window and rendered into it, but GLFW can't open windows at arbitrarily small resolutions, so I create the window bigger than required, then draw a small rectangle that is the correct size so that it has one fragment per tile that maps to the source image. The entire MSVC2013 solution is available at <https://bitbucket.org/Gibgezr/mosaicmaker>
It requires GLFW/FreeImage/GLEW/GLM to compile, and OpenGL 3.3 or better drivers/video card to run.
## Fragment Shader Source
```
#version 330 core
uniform sampler2D sourceTexture;
uniform sampler2D mosaicTexture;
in vec2 v_texcoord;
out vec4 out_Color;
void main(void)
{   
    ivec2 sourceSize = textureSize(sourceTexture, 0);
    ivec2 mosaicSize = textureSize(mosaicTexture, 0);
    float num_pixels = mosaicSize.x/45.f;
    vec4 myTexel;
    float difference = 0.f;
    //initialize smallest difference to a large value
    float smallest_difference = 255.0f*255.0f*255.0f;
    int closest_x = 0, closest_y = 0;
    vec2 pixel_size_src = vec2( 1.0f/sourceSize.x, 1.0f/sourceSize.y);
    vec2 pixel_size_mosaic = vec2( 1.0f/mosaicSize.x , 1.0f/mosaicSize.y);
    vec2 incoming_texcoord;
    //adjust incoming uv to bottom corner of the tile space
    incoming_texcoord.x =  v_texcoord.x - 1.0f/(sourceSize.x / num_pixels * 2.0f);
    incoming_texcoord.y =  v_texcoord.y - 1.0f/(sourceSize.y / num_pixels * 2.0f);
    vec2 texcoord_mosaic;
    vec2 pixelcoord_src, pixelcoord_mosaic;
    vec4 pixel_src, pixel_mosaic;
    //loop over all of the mosaic tiles
    for(int i = 0; i < 45; ++i)
    {
        for(int j = 0; j < 45; ++j)
        {
            difference = 0.f;
            texcoord_mosaic = vec2(j * pixel_size_mosaic.x * num_pixels, i * pixel_size_mosaic.y * num_pixels);
            //loop over all of the pixels in the images, adding to the difference
            for(int y = 0; y < num_pixels; ++y)
            {
                for(int x = 0; x < num_pixels; ++x)
                {
                    pixelcoord_src = vec2(incoming_texcoord.x + x * pixel_size_src.x, incoming_texcoord.y + y * pixel_size_src.y);                  
                    pixelcoord_mosaic = vec2(texcoord_mosaic.x + x * pixel_size_mosaic.x, texcoord_mosaic.y + y * pixel_size_mosaic.y); 
                    pixel_src = texture(sourceTexture, pixelcoord_src);
                    pixel_mosaic = texture(mosaicTexture, pixelcoord_mosaic);
                    pixel_src *= 255.0f;
                    pixel_mosaic *= 255.0f;
                    difference += (pixel_src.x - pixel_mosaic.x) * (pixel_src.x - pixel_mosaic.x) * 0.5f+
                        (pixel_src.y - pixel_mosaic.y) * (pixel_src.y - pixel_mosaic.y) +
                        (pixel_src.z - pixel_mosaic.z) * (pixel_src.z - pixel_mosaic.z) * 0.17f;
                }
            }
            if(difference < smallest_difference)
            {
                smallest_difference = difference;
                closest_x = j;
                closest_y = i;
            }               
        }
    }
    myTexel.x = float(closest_x)/255.f;
    myTexel.y = float(closest_y)/255.f;
    myTexel.z = 0.f;
    myTexel.w = 0.f;    
    out_Color = myTexel;
}
```
## Results
The images render almost instantly, so the parallelization was a success. The downside is that I can't make the individual fragments rely on the output of any other fragments, so there is no way to get the significant quality increase you can get by not picking the same tile twice within a certain range. So, fast results, but limited quality due to massive repetitions of tiles. All in all, it was fun.
<https://i.stack.imgur.com/MO7NU.jpg> for full-size versions.
![enter image description here](https://i.stack.imgur.com/9tVlw.jpg)
[Answer]
# Python
Here goes for the first Python solution, using a mean approach. We can evolve from here. The rest of the images are [here](https://i.stack.imgur.com/6pDtQ.jpg).
![sunday](https://i.stack.imgur.com/wsVf2.jpg)
![steve](https://i.stack.imgur.com/Qx85w.jpg)
```
from PIL import Image
import numpy as np
def calcmean(b):
    meansum = 0
    for k in range(len(b)):
        meansum = meansum + (k+1)*b[k]
    return meansum/sum(b)    
        
def gettiles(imageh,w,h):
    k = 0 
    tiles = {}
    for x in range(0,imageh.size[0],w):
        for y in range(0,imageh.size[1],h):
            a=imageh.crop((x, y, x + w, y + h))
            b=a.resize([1,1], Image.ANTIALIAS)
            tiles[k] = [a,x,y,calcmean(b.histogram()[0:256]) \
                             ,calcmean(b.histogram()[256:256*2]) \
                             ,calcmean(b.histogram()[256*2:256*3])]
            k = k + 1
    return tiles
w = 48 
h = 48
srch = Image.open('25745_avatars.png').convert('RGB')
files = ['21104_spheres.png', '45148_sunday.jpg', '47824_steve.jpg', '61555_kiss.jpg', '77388_lightbulb.png']
for f in range(len(files)):
    desh = Image.open(files[f]).convert('RGB')
    
    srctiles = gettiles(srch,w,h)
    destiles = gettiles(desh,w,h)
    
    #build proximity matrix 
    pm = np.zeros((len(destiles),len(srctiles)))
    for d in range(len(destiles)):
        for s in range(len(srctiles)):
            pm[d,s] = (srctiles[s][3]-destiles[d][3])**2 + \
                      (srctiles[s][4]-destiles[d][4])**2 + \
                      (srctiles[s][5]-destiles[d][5])**2
    
    for k in range(len(destiles)):
        j = list(pm[k,:]).index(min(pm[k,:]))
        desh.paste(srctiles[j][0], (destiles[k][1], destiles[k][2]))
    
    desh.save(files[f].replace('.','_m'+'.'))
   
```
[Answer]
# Yet Another Python Solution - Average Based (RGB vs L*a*b\*)
## Results (There are some slightly differences)
**Bulb - RGB**
[full view](http://www.pictureshack.us/images/3669_bulb_photomosaic.png)
![bulb_rgb](https://i.stack.imgur.com/OH0V4.png)
**Bulb - Lab**
[full view](http://www.pictureshack.us/images/62570_bulb_photomosaic.png)
![bulb_lab](https://i.stack.imgur.com/HUoC5.png)
**Steve - RGB**
[full view](http://www.pictureshack.us/images/3496_steve_photomosaic.png)
![steve_rgb](https://i.stack.imgur.com/vDlj4.png)
**Steve - Lab**
[full view](http://www.pictureshack.us/images/95937_steve_photomosaic.png)
![steve_lab](https://i.stack.imgur.com/c6fgU.png)
**Spheres - RGB**
[full view](http://www.pictureshack.us/images/9886_spheres_photomosaic.png)
![spheres_rgb](https://i.stack.imgur.com/sFZa6.png)
**Spheres - Lab**
[full view](http://www.pictureshack.us/images/2992_spheres_photomosaic.png)
![spheres_lab](https://i.stack.imgur.com/EvCM6.png)
**Sunday - RGB**
[full view](http://www.pictureshack.us/images/83218_sunday_photomosaic.png)
![sunday_rgb](https://i.stack.imgur.com/y3I6w.png)
**Sunday - Lab**
[full view](http://www.pictureshack.us/images/14201_sunday_photomosaic.png)
![sunday_lab](https://i.stack.imgur.com/rMmGC.png)
**Kiss - RGB**
[full view](https://i.stack.imgur.com/ULEpd.jpg)
![kiss_rgb](https://i.stack.imgur.com/jf9VA.jpg)
**Kiss - Lab**
[full view](https://i.stack.imgur.com/b0wfK.jpg)
![kiss_lab](https://i.stack.imgur.com/QqA8B.jpg)
## Code
requires [python-colormath](http://python-colormath.readthedocs.org/en/latest/index.html) for Lab
```
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from PIL import Image
from colormath.color_objects import LabColor,sRGBColor
from colormath.color_conversions import convert_color
from colormath.color_diff import delta_e_cie1976
def build_photomosaic_ils(mosaic_im,target_im,block_width,block_height,colordiff,new_filename):
    mosaic_width=mosaic_im.size[0]              #dimensions of the target image
    mosaic_height=mosaic_im.size[1]
    target_width=target_im.size[0]              #dimensions of the target image
    target_height=target_im.size[1]
    target_grid_width,target_grid_height=get_grid_dimensions(target_width,target_height,block_width,block_height)       #dimensions of the target grid
    mosaic_grid_width,mosaic_grid_height=get_grid_dimensions(mosaic_width,mosaic_height,block_width,block_height)       #dimensions of the mosaic grid
    target_nboxes=target_grid_width*target_grid_height
    mosaic_nboxes=mosaic_grid_width*mosaic_grid_height
    print "Computing the average color of each photo in the mosaic..."
    mosaic_color_averages=compute_block_avg(mosaic_im,block_width,block_height)
    print "Computing the average color of each block in the target photo ..."
    target_color_averages=compute_block_avg(target_im,block_width,block_height)
    print "Computing photomosaic ..."
    photomosaic=[0]*target_nboxes
    for n in xrange(target_nboxes):
        print "%.2f " % (n/float(target_nboxes)*100)+"%"
        for z in xrange(mosaic_nboxes):
            current_diff=colordiff(target_color_averages[n],mosaic_color_averages[photomosaic[n]])
            candidate_diff=colordiff(target_color_averages[n],mosaic_color_averages[z])
            if(candidate_diff<current_diff):
                photomosaic[n]=z
    print "Building final image ..."
    build_final_solution(photomosaic,mosaic_im,target_nboxes,target_im,target_grid_width,block_height,block_width,new_filename)
def build_initial_solution(target_nboxes,mosaic_nboxes):
    candidate=[0]*target_nboxes
    for n in xrange(target_nboxes):
        candidate[n]=random.randint(0,mosaic_nboxes-1)
    return candidate
def build_final_solution(best,mosaic_im,target_nboxes,target_im,target_grid_width,block_height,block_width,new_filename):
    for n in xrange(target_nboxes):
        i=(n%target_grid_width)*block_width             #i,j -> upper left point of the target image
        j=(n/target_grid_width)*block_height
        box = (i,j,i+block_width,j+block_height)        
        #get the best photo from the mosaic
        best_photo_im=get_block(mosaic_im,best[n],block_width,block_height)
        #paste the best photo found back into the image
        target_im.paste(best_photo_im,box)
    target_im.save(new_filename);
#get dimensions of the image grid
def get_grid_dimensions(im_width,im_height,block_width,block_height):
    grid_width=im_width/block_width     #dimensions of the target image grid
    grid_height=im_height/block_height
    return grid_width,grid_height
#compute the fitness of given candidate solution
def fitness(candidate,mosaic_color_averages,mosaic_nboxes,target_color_averages,target_nboxes):
    error=0.0
    for i in xrange(target_nboxes):
        error+=colordiff_rgb(mosaic_color_averages[candidate[i]],target_color_averages[i])
    return error
#get a list of color averages, i.e, the average color of each block in the given image
def compute_block_avg(im,block_height,block_width):
    width=im.size[0]
    height=im.size[1]
    grid_width_dim=width/block_width                    #dimension of the grid
    grid_height_dim=height/block_height
    nblocks=grid_width_dim*grid_height_dim              #number of blocks
    avg_colors=[]
    for i in xrange(nblocks):
        avg_colors+=[avg_color(get_block(im,i,block_width,block_height))]
    return avg_colors
#returns the average RGB color of a given image
def avg_color(im):
    avg_r=avg_g=avg_b=0.0
    pixels=im.getdata()
    size=len(pixels)
    for p in pixels:
        avg_r+=p[0]/float(size)
        avg_g+=p[1]/float(size)
        avg_b+=p[2]/float(size)
    return (avg_r,avg_g,avg_b)
#get the nth block of the image
def get_block(im,n,block_width,block_height):
    width=im.size[0]
    grid_width_dim=width/block_width                        #dimension of the grid
    i=(n%grid_width_dim)*block_width                        #i,j -> upper left point of the target block
    j=(n/grid_width_dim)*block_height
    box = (i,j,i+block_width,j+block_height)
    block_im = im.crop(box)
    return block_im
#calculate color difference of two pixels in the RGB space
#less is better
def colordiff_rgb(pixel1,pixel2):
    delta_red=pixel1[0]-pixel2[0]
    delta_green=pixel1[1]-pixel2[1]
    delta_blue=pixel1[2]-pixel2[2]
    fit=delta_red**2+delta_green**2+delta_blue**2
    return fit
#http://python-colormath.readthedocs.org/en/latest/index.html
#calculate color difference of two pixels in the L*ab space
#less is better
def colordiff_lab(pixel1,pixel2):
    #convert rgb values to L*ab values
    rgb_pixel_1=sRGBColor(pixel1[0],pixel1[1],pixel1[2],True)
    lab_1= convert_color(rgb_pixel_1, LabColor)
    rgb_pixel_2=sRGBColor(pixel2[0],pixel2[1],pixel2[2],True)
    lab_2= convert_color(rgb_pixel_2, LabColor)
    #calculate delta e
    delta_e = delta_e_cie1976(lab_1, lab_2)
    return delta_e
if __name__ == '__main__':
    mosaic="images/25745_avatars.png"
    targets=["images/lightbulb.png","images/steve.jpg","images/sunday.jpg","images/spheres.png","images/kiss.jpg"]
    target=targets[0]
    mosaic_im=Image.open(mosaic)
    target_im=Image.open(target)
    new_filename=target.split(".")[0]+"_photomosaic.png"
    colordiff=colordiff_rgb
    build_photomosaic_ils(mosaic_im,target_im,48,48,colordiff,new_filename)
```
]