# Image creation for every possible color

I'm trying to create an image with every possible color. It would start with a seed pixel and then place randomly-generated RGB pixels around it. Future placements would be based on whichever open spot had the average of the pixels surrounding it closest to the new color to be placed.

from PIL import Image
import numpy as np
from random import randint
import sys
import random
import itertools

sys.setcheckinterval(10000)

def moddistance3(x1,y1,z1,x2,y2,z2):  #get relative distance between two 3D points
x = abs(x1 - x2)
y = abs(y1 - y2)
z = abs(z1 - z2)
return (x + y + z)

def genColor(unused): #generate random color (not used anymore)
test = 0
while test == 0:
red = randint(0,255)
green = randint(0,255)
blue = randint(0,255)
if unused[red,green,blue] == 1:
test = 1
return (red,green,blue)

def surroundAvg(points,unfilled):
surrounding = {}
count = len(points)
for inc in xrange(count):
neighbors = filledNeighbors(points[inc][0],points[inc][1],unfilled)
nearcount = len(neighbors)
pixred = 0
pixgreen = 0
pixblue = 0
for num in xrange(nearcount):
(temp_red,temp_green,temp_blue) = pixels[neighbors[num][0],neighbors[num][1]]
pixred = pixred + temp_red
pixgreen = pixgreen + temp_green
pixblue = pixblue + temp_blue
pixred = pixred / nearcount
pixgreen = pixgreen / nearcount
pixblue = pixblue / nearcount
surrounding[(points[inc][0],points[inc][1])] = (pixred,pixgreen,pixblue)
return surrounding

def genPoint(perim,unfilled,averages,red,green,blue):
num_test = len(perim)
test = 0
least_diff = 9999
nearby = []
for point in xrange(num_test):
i = perim[point][0]
j = perim[point][1]
pixred = averages[(i,j)][0]
pixgreen = averages[(i,j)][1]
pixblue = averages[(i,j)][2]
diff = abs(red - pixred) + abs(green - pixgreen) + abs(blue - pixblue)
if diff < least_diff or test == 0:
least_diff = diff
newx = i
newy = j
test = 1
return newx,newy

def cubegen():  #create the cube of colors with each color having its own number
cube = np.zeros(16777216,dtype=np.object)
num = 0
for red in xrange(0,256):
for green in xrange(0,256):
for blue in xrange(0,256):
cube[num] = [red,green,blue]
num += 1
return cube

def getNeighbors(x,y,unfilled):
Prod = itertools.product
toremove = []
neighbors = list(Prod(range(x-1,x+2),range(y-1,y+2)))
for num in xrange(len(neighbors)):
i,j = neighbors[num]
if j > 4095 or i > 4095 or unfilled[(i,j)] == 0 or j < 0 or i < 0:
toremove.append((i,j))
map(neighbors.remove,toremove)
return neighbors

def filledNeighbors(x,y,unfilled):
Prod = itertools.product
toremove = []
neighbors = list(Prod(range(x-1,x+2),range(y-1,y+2)))
#neighbors = filter(lambda i,j: j < 4096 and i < 4096 and unfilled[i,j] == 0 and j > -1 and i > -1,allneighbors)
for num in xrange(len(neighbors)):
i,j = neighbors[num]
if j > 4095 or i > 4095 or unfilled[(i,j)] == 1 or j < 0 or i < 0:
toremove.append((i,j))
map(neighbors.remove,toremove)
return neighbors

img = Image.new('RGB', (4096,4096)) # create a new black image
pixels = img.load() # create the pixel map

colorList = range(16777216)
colorCube = cubegen()
print("Color cube created successfully")
unfilled = {}
for x in xrange(4096):
for y in xrange(4096):
unfilled[(x,y)] = 1
startx = 2048
starty = 2048
random.shuffle(colorList)
print("Color list shuffled successfully")
color = colorList[0]
(red,green,blue) = colorCube[color]
pixels[startx,starty] = (red,green,blue)
unfilled[(startx,starty)] = 0
perim_empty = getNeighbors(startx,starty,unfilled)
edge = []
#edge.append((startx,starty))
avg = surroundAvg(perim_empty,unfilled)
print("First point placed successfully.")
#appendEdge = edge.append
#removeEdge = edge.remove
appendPerim = perim_empty.append
removePerim = perim_empty.remove
updateAvg = avg.update

for iteration in xrange(1,16777216):
temp = {}
color = colorList[iteration]
(red,green,blue) = colorCube[color]
(i,j) = genPoint(perim_empty,unfilled,avg,red,green,blue)
unfilled[(i,j)] = 0
pixels[i,j] = (red,green,blue)
new_neighbors = getNeighbors(i,j,unfilled)
map(appendPerim,new_neighbors)
temp = surroundAvg(new_neighbors,unfilled)
updateAvg(temp)
removePerim((i,j))
#appendEdge((i,j))

#if iteration % 20 == 0:
#   toremove = []
#   appendToRemove = toremove.append
#   for num in xrange(len(edge)):
#       nearby = getNeighbors(edge[num][0],edge[num][1],unfilled)
#       if len(nearby) == 0:
#           appendToRemove(edge[num])
#for num in xrange(len(toremove)):
#   edge.remove(toremove[num])
#   map(removeEdge,toremove)

if iteration % 500 == 0:
print("Iteration %d complete" %iteration)
if iteration == 100000 or iteration == 500000 or iteration ==1000000 or iteration == 5000000 or iteration == 10000000 or iteration == 15000000:
img.save("Perimeter Averaging -- %d iterations.bmp" %iteration)
img.save("Perimeter Averaging Final.bmp")
img.show()


The problem is that when I try to run this, it takes days to even go through 1,000,000 of the colors, and slows down considerably as it goes. I can't figure out how to make it take less time, and I know there must be a way to do this that doesn't take months. I'm new to code and am teaching myself, so please forgive any obvious fixes I've totally overlooked.

• I don't fully understand your goal, but there are basically two approaches you could take to finding the problem: 1) Analyse the runtime complexity. Are you doing something here that is O(N^4) or O(2^N)? 2) Profile the code and see where it is spending most of the time. – bsa Jun 3 '17 at 23:58
• It spends most of its time in the genPoint function. I can't figure out how to make it take less time to go through its processes. – wolf53135 Jun 4 '17 at 18:26

### 1. Review

Before addressing the performance, the code needs some work:

1. Instead of the mysterious number 4096, make a named constant and show how you compute it:

IMG_SIZE = 2**12  # Side of square image with 2**24 pixels.

2. In cubegen you set dtype=np.object. But this is wasteful: the result of the cubegen is 1,207,959,648 bytes in size. Since all the numbers in this array are in the range 0 to 255, you could use np.uint8, and this would reduce the memory usage to 50,331,792 bytes (one-24th of the size).

3. cubegen can be implemented using numpy.arange, numpy.meshgrid, and numpy.stack:

def cubegen():
"""Return array of colour values with shape (256**3, 3)."""
c = np.arange(256, dtype=np.uint8)
r, g, b = np.meshgrid(c, c, c, indexing='ij', copy=False)
return np.stack((r, g, b), axis=-1).reshape(-1, 3)


This is about ten times faster than the code in the post.

4. Shuffle the array of colours directly, using numpy.random.shuffle, and avoid the need for colorList.

5. The dictionary of unfilled pixels is 1,476,395,080 bytes in size. To save a lot of space, make a NumPy array of Booleans:

unfilled = np.ones((IMG_SIZE,) * 2, dtype=np.bool)


This is just 16,777,328 bytes long (one-88th of the size).

6. Instead of iterating over a range, iterate directly over the colour cube.

7. getNeighbors and filledNeighbors are identical except for the test against the unfilled array. So these two functions can be combined into one function that takes the condition to test. (In fact we'll see later that we don't need filledNeighbors.)

8. Now that unfilled is a NumPy array, getNeighbours can be simplified:

# Offset of the 8 neighbours of a point.
DELTA_I = np.array([-1, -1, -1,  0, 0,  1, 1, 1])
DELTA_J = np.array([-1,  0,  1, -1, 1, -1, 0, 1])

def neighbors(point, unfilled):
"""Return arrays of coordinates of unfilled neighbors of point."""
i, j = point
ni = DELTA_I + i
nj = DELTA_J + j


(Note that this only considers the 8 neighbors of the point, whereas the code in the post considers the point too.)

9. The code in the post recomputes the surrounding average for every point in the perimeter. But this involves a lot of wasted work, because for most pixels in the perimeter, the average does not change from one step to the next. Only the pixels that are neighbors of the pixel that was just placed need to change. So if we maintain a running average of the filled-in neighbours of every pixel in the image, we can avoid nearly all this work.

### 2. Revised code

This answer grew quite long, so there are a few other tricks in the revised code that are not mentioned above. See if you can spot them and figure out how they work.

import numpy as np
import scipy.misc

def cubegen():
"""Return array of colour values with shape (256**3, 3)."""
c = np.arange(256, dtype=np.uint8)
r, g, b = np.meshgrid(c, c, c, indexing='ij', copy=False)
return np.stack((r, g, b), axis=-1).reshape(-1, 3)

# Side of square image with 2**24 pixels.
IMG_SIZE = 2 ** 12

# Offset of the 8 neighbours of a point.
DELTA_I = np.array([-1, -1, -1,  0, 0,  1, 1, 1])
DELTA_J = np.array([-1,  0,  1, -1, 1, -1, 0, 1])

# States of a pixel: empty, on the perimeter, or full.
EMPTY, PERIM, FULL = range(3)

def main(n=2**24, filename='cr164870.png'):
"""Place the first n pixels and save the result to filename."""
color_cube = cubegen()
np.random.shuffle(color_cube)

# Shape of image, plus 1-pixel-wide border all around.
img_shape = (IMG_SIZE + 2,) * 2

# RGB image under construction.
pixels = np.zeros(img_shape + (3,), dtype=np.uint8)

# Average color of the filled-in neighbors of each point.
average = np.zeros_like(pixels, dtype=np.float32)

# Number of filled-in neighbors of each point.
count_neighbors = np.zeros(img_shape, dtype=np.uint8)

# State of each pixel.
state = np.zeros(img_shape, dtype=np.uint8)

# Fill in the 1-pixel-wide border (this prevents the neighbors
# function from going out of bounds).
state[0,:] = state[-1,:] = state[:,0] = state[:,-1] = FULL

def neighbors(point):
"""Return arrays of coordinates of unfilled neighbors of point."""
i, j = point
ni = DELTA_I + i
nj = DELTA_J + j
mask = state[ni, nj] != FULL

# Start in the center of the image.
mid = IMG_SIZE // 2 + 1

# Arrays of coordinates of points on the perimeter.
perim_i = np.array([mid])
perim_j = np.array([mid])

for _, color in zip(range(n), (color_cube)):
# Manhattan distance in color space between color and average
# neighborhood of point.
dist = np.abs(color - average[perim_i, perim_j]).sum(axis=-1)

# Find the closest location on perimeter by color distance.
closest = np.argmin(dist)
point = perim_i[closest], perim_j[closest]

# Fill it in.
state[point] = FULL
pixels[point] = color
perim_i = np.delete(perim_i, closest)
perim_j = np.delete(perim_j, closest)

# Get the coordinates of the unfilled neighbors of point.
ni, nj = neighbors(point)

# Update the average neighborhood color of each neighbor of point.
count_neighbors[ni, nj] += 1
m = count_neighbors[ni, nj].reshape(-1, 1)
average[ni, nj] = (average[ni, nj] * (m - 1) + color) / m

# Update the perimeter.
mask = state[ni, nj] == EMPTY
state[ni, nj] = PERIM
perim_i = np.append(perim_i, ni)
perim_j = np.append(perim_j, nj)

scipy.misc.imsave(filename, pixels[1:-1,1:-1])


### 3. Performance analysis

For every pixel in the image, the algorithm has to find the best point in the perimeter. This takes time proportional to the length of perimeter. How long is that? Well, let's take a look at the output. Here's an image showing the position after a million pixels have been placed:

Its clear that the filled-in pixels form a roughly circular region, and so the length of the perimeter is proportional to the square root of the number of pixels placed so far. This means that we expect the runtime for placing $n$ pixels to be roughly proportional to $n^{3\over2}$.

The table below has some measurements of the code in §2. Here $t(n)$ is the time in seconds taken to place the first $n$ pixels on the image (not including the time taken to initialize the color_cube array at the beginning, or the time to save the image at the end).

      n       t(n)
------- ----------
3000       0.14
10000       1.14
30000       3.98
100000      17.86
300000      86.60
1000000     430.49


Let's plot these and fit to $t(n) = an^{3\over2}$:

Based on this, I estimate that it will take the code in §2 about 30,000 seconds (about 8 hours) to plot all 16,777,216 pixels. This is a lot better than the original code, but it seems as if we ought to be able to improve it further.

### 4. Improved algorithm

To get runtime that's better than $Θ(n^{3/2})$ we must avoid doing work at each step that's proportional to the length of the perimeter. But then how can we find the best point?

Well, finding the best point is a nearest-neighbour search in color space, and this can be solved efficiently using a space-partitioning data structure like a $k$-d tree. With careful attention to book-keeping, I think it should be possible to get the runtime down to $O(n \log n)$.

Unfortunately, SciPy's built-in scipy.spatial.KDTree is not suitable for this use case, because it is not modifiable in-place (after changing the underlying data you have to rebuild the whole tree).

So you're stuck with having to write your mutable $k$-d tree, I'm afraid. Good luck!

### 5. Update

8 hours didn't seem so long, so I ran it overnight. The total runtime was about 28,000 seconds, so the rough estimate above was pretty good. Unfortunately the whole image (at over 30 MB) is too big to upload to Stack Exchange, but here's a version that's scaled down to 768×768:

• Wow! Thank you so much for your help! This is so much more than I expected to receive lol. I definitely need to work on my knowledge of modules and Python commands. – wolf53135 Jun 5 '17 at 11:44
• If I could vote for this answer more than once, I would. Wow. – bsa Jun 5 '17 at 13:16
• @Gareth Rees isn´t p the parameter of KDTree.query which you can set to 1 for Manhattan distance? – juvian Jun 6 '17 at 18:02
• @juvian: So it is! I'll update the post accordingly. – Gareth Rees Jun 6 '17 at 18:41