See the other notebook for the grisly details and dead-ends.
Requirements:
numpyscipyscikit-learnpillowI recommend installing them with conda install.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
from scipy import signal
In [3]:
nx, ny = 100, 100
z = np.random.rand(nx, ny)
sizex, sizey = 30, 30
x, y = np.mgrid[-sizex:sizex+1, -sizey:sizey+1]
g = np.exp(-0.333*(x**2/float(sizex)+y**2/float(sizey)))
f = g/g.sum()
z = signal.convolve(z, f, mode='valid')
z = (z - z.min())/(z.max() - z.min())
In [8]:
cd ~/Dropbox/dev/rainbow/notebooks
In [10]:
# Note: interpolation introduces new colours.
plt.imshow(z, cmap="spectral")
#plt.imshow(z, cmap="spectral", interpolation='none')
plt.axis('off')
plt.savefig('data/cbar/test.jpg', bbox_inches='tight')
plt.show()
In [11]:
cmaps = [('Perceptually Uniform Sequential',
['viridis', 'inferno', 'plasma', 'magma']),
('Sequential', ['Blues', 'BuGn', 'BuPu',
'GnBu', 'Greens', 'Greys', 'Oranges', 'OrRd',
'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu',
'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd']),
('Sequential (2)', ['afmhot', 'autumn', 'bone', 'cool',
'copper', 'gist_heat', 'gray', 'hot',
'pink', 'spring', 'summer', 'winter']),
('Diverging', ['BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr',
'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral',
'seismic']),
('Qualitative', ['Accent', 'Dark2', 'Paired', 'Pastel1',
'Pastel2', 'Set1', 'Set2', 'Set3']),
('Miscellaneous', ['gist_earth', 'terrain', 'ocean', 'gist_stern',
'brg', 'CMRmap', 'cubehelix',
'gnuplot', 'gnuplot2', 'gist_ncar',
'nipy_spectral', 'jet', 'rainbow',
'gist_rainbow', 'hsv', 'flag', 'prism'])]
In [13]:
cd ~/Dropbox/dev/rainbow/notebooks
In [14]:
from PIL import Image
# img = Image.open('data/cbar/boxer.png')
# img = Image.open('data/cbar/fluid.png')
# img = Image.open('data/cbar/lisa.png')
# img = Image.open('data/cbar/redblu.png')
# img = Image.open('data/cbar/seismic.png')
# img = Image.open('data/cbar/drainage.jpg')
img = Image.open('data/cbar/test.png')
img
Out[14]:
In [15]:
img.size
Out[15]:
Instead of taking a random sample, let's take all pixels from a smaller version of the image.
In [16]:
# def resize_if_necessary(img, max_size=256):
# h, w = img.size
# if h * w > max_size**2:
# img = img.resize((max_size, max_size))
# return img
In [17]:
# img_sm = resize_if_necessary(img)
Cast as an array, and ignore the last channel (alpha), if there is one.
In [30]:
n_colours = 256
In [31]:
from sklearn.cluster import KMeans
from sklearn.utils import shuffle
In [32]:
im = np.asarray(img) / 255
h, w, d = im.shape
im_ = im.reshape((w * h, d))[:, :3]
sample = shuffle(im_, random_state=0)[:1000] # Defines training set size
Train:
In [33]:
kmeans = KMeans(n_clusters=n_colours, random_state=0).fit(sample)
Now I can make an RGB palette p — also known as a codebook in information theory terms:
In [34]:
p = kmeans.cluster_centers_
# I don't know why I need to do this, but I do. Floating point precision maybe.
p[p > 1] = 1
p[p < 0] = 0
The only problem with this p is that it is not in order — that it, there cluster centres are more or less randomly arranged. We will fix that in the next section.
The vector p is actually all we need, but if you want to see what the quantized image looks like, carry on:
In [35]:
# labels = kmeans.predict(im_)
# def recreate_image(palette, labels, h, w):
# image = np.zeros((h, w, palette.shape[1]))
# label_idx = 0
# for i in range(h):
# for j in range(w):
# image[i][j] = palette[labels[label_idx]]
# label_idx += 1
# return image
# q = recreate_image(p, labels, h, w)
# plt.imshow(q)
# plt.show()
In [36]:
from mpl_toolkits.mplot3d import Axes3D
# Set up the figure
fig = plt.figure(figsize=(8, 8))
# Result of TSP solver
ax = fig.add_subplot(111, projection='3d')
ax.scatter(*p.T, c=p, lw=0, s=40, alpha=1)
ax.plot(*p.T, color='k', alpha=0.4)
ax.set_title('Codebook')
plt.show()
To regularize this cloud, let's try wrapping it in a convex hull.
In [40]:
from scipy.spatial import ConvexHull
hull = ConvexHull(p)
In [51]:
# Set up the figure
fig = plt.figure(figsize=(8, 8))
# Result of TSP solver
ax = fig.add_subplot(111, projection='3d')
ax.scatter(*p.T, c=p, lw=0, s=40, alpha=1)
for simplex in hull.simplices:
ax.plot(p[simplex, 0], p[simplex, 1], p[simplex, 2], 'b', alpha=0.5)
I propose starting at the dark end (the end of the line nearest black) and crawling aling the line of points from there. This will make a nice organized sequence of codes — in our case, this will be the colourmap.
We can solve this problem as the travelling salesman problem. Find the shortest tour from black to 'the other end'.
To start with, we need the distances between all points. This is just a norm, but there's a convenient function scipy.spatial.pdist for finding distance pairs in n-space. Then squareform casts it into a square symmetric matrix, which is what we need for our TSP solver.
Other than creating a naive TSP solver in Python – let's face it, it'll be broken or slow or non-optimal — there are three good TSP solver options:
LKH and Concorde can be used via the TSP Python package (but note that it used to be called pyconcorde so you need to change the names of some functions — look at the source or use my fork.
Note that you need to add the Concorde and LKH libs to PATH as mentioned in the docs for pytsp.
In [25]:
from pytsp import run, dumps_matrix
We want to start the traversal of the locus at the end nearest black. This is a heuristic based on the observation that we tend to map darker colours to lower values. Adding this point allows us to start the TSP there, because it will move to the next nearest point — I think. We will remove it later.
Add black point to p:
In [26]:
p = np.vstack([[[0,0,0]], p])
In [27]:
p[:6]
Out[27]:
Remember that these points are essentially in random order. We are going to solve the find the shortest Hamiltonian path to organize them.
To do this, we need weights for the edges — and we'll use the distances between the points. There's a convenient function, scipy.spatial.pdist for finding distances in n-space. We will use the 2-norm, but the pdist function has a great many other metrics.
The convenience function squareform organizes the distances into a matrix.
In [28]:
from scipy.spatial.distance import pdist, squareform
# Make distance matrix.
dists = squareform(pdist(p, 'euclidean'))
# The values in `dists` are floats in the range 0 to sqrt(3).
# Normalize the values to int16s.
d = 32767 * dists / np.sqrt(3)
d = d.astype(np.int16)
# To use a TSP algo to solve the shortest Hamiltonian path problem,
# we need to add a point that is zero units from every other point.
row, col = d.shape
d = np.insert(d, row, 0, axis=0)
d = np.insert(d, col, 0, axis=1)
The zero-point trick is legit. Reference from E. L. Lawler, Jan Karel Lenstra, A. H. G. Rinnooy Kan, D. B. Shmoys (1985). The Traveling Salesman Problem: A Guided Tour of Combinatorial Optimization, 1st Edition. Wiley. 476 pp. ISBN 978-0471904137.
In [29]:
d
Out[29]:
Concorde algorithm — a little slower than LKH:
In [94]:
outf = "/tmp/myroute_concorde.tsp"
with open(outf, 'w') as f:
f.write(dumps_matrix(d, name="My Route"))
In [95]:
tour_concorde = run(outf, start=0, solver="Concorde")
LKH implementation:
In [96]:
outf = "/tmp/myroute_lkh.tsp"
with open(outf, 'w') as f:
f.write(dumps_matrix(d, name="My Route"))
In [97]:
tour_lkh = run(outf, start=0, solver="LKH")
In [98]:
#result = np.array(tour_concorde['tour'])
result = np.array(tour_lkh['tour'])
In [99]:
result
Out[99]:
Now result is the indices of points for the shortest path, shape (256,). And p is our quantized colormap, shape (256, 3). So we can select the points easily for an ordered colourmap.
The offsets are to account for the fact that we added a black point at the start and a zero point at the end.
In [100]:
c = p[result[1:-1]]
Ideally I'd like all the distances too, but it wouldn't be too hard to compute these.
Now let's look at it all.
In [101]:
from mpl_toolkits.mplot3d import Axes3D
# Set up the figure
fig = plt.figure(figsize=(8, 8))
# Result of TSP solver
ax = fig.add_subplot(111, projection='3d')
ax.scatter(*c.T, c=c, lw=0, s=40, alpha=1)
ax.plot(*c.T, color='k', alpha=0.4)
ax.set_title('TSP solver')
plt.show()
Check below an interactive version of the 3D plot. May help when there are complicated paths between points. You need to install plotly and colorlover (with pip) if you don't already have them.
In [103]:
import plotly.graph_objs as go
import colorlover as cl
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=True)
cb = cl.to_rgb(tuple(map(tuple, c*255)))
trace = go.Scatter3d(
name='TSP Sover',
x = c[:,0], y = c[:,1], z = c[:,2],
marker = dict(
size=4.,
color=cb
),
line=dict(
color='#000',
width=1,
),
)
data = [trace]
# Set the different layout properties of the figure:
layout = go.Layout(
autosize=False,
width=600,
height=600,
margin = dict(
t=0,b=0,l=0,r=0
),
scene = go.Scene(
xaxis=dict(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=True,
backgroundcolor='rgb(230, 230,230)'
),
yaxis=dict(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=True,
backgroundcolor='rgb(230, 230,230)'
),
zaxis=dict(
gridcolor='rgb(255, 255, 255)',
zerolinecolor='rgb(255, 255, 255)',
showbackground=True,
backgroundcolor='rgb(230, 230,230)'
),
aspectmode='cube',
camera=dict(
eye=dict(
x=1.7,
y=-1.7,
z=1,
)
),
)
)
fig = go.Figure(data=data, layout=layout)
iplot(fig, show_link=False)
In [104]:
from scipy.spatial import cKDTree
kdtree = cKDTree(c)
In [105]:
dx, ix = kdtree.query(im_)
In [106]:
plt.imshow(ix.reshape((h, w)), cmap='viridis')
plt.colorbar()
plt.show()
In [107]:
plt.imshow(dx.reshape((h, w)))
plt.colorbar()
plt.show()
In [108]:
fig = plt.figure(figsize=(18, 5))
ax0 = fig.add_subplot(131)
plt.imshow(im, interpolation='none')
ax1 = fig.add_subplot(132, projection='3d')
ax1.scatter(*c.T, c=c, lw=0, s=40, alpha=1)
ax1.plot(*c.T, color='k', alpha=0.5)
ax1.text(*c[0], ' start')
ax1.text(*c[-1], ' end')
ax2 = fig.add_subplot(133)
plt.imshow(ix.reshape((h, w)), cmap="jet", interpolation='none')
plt.colorbar(shrink=0.75)
plt.show()
In [109]:
cmaps = [('Perceptually Uniform Sequential',
['viridis', 'inferno', 'plasma', 'magma']),
('Sequential', ['Blues', 'BuGn', 'BuPu',
'GnBu', 'Greens', 'Greys', 'Oranges', 'OrRd',
'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu',
'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd']),
('Sequential (2)', ['afmhot', 'autumn', 'bone', 'cool',
'copper', 'gist_heat', 'gray', 'hot',
'pink', 'spring', 'summer', 'winter']),
('Diverging', ['BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr',
'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral',
'seismic']),
('Qualitative', ['Accent', 'Dark2', 'Paired', 'Pastel1',
'Pastel2', 'Set1', 'Set2', 'Set3']),
('Miscellaneous', ['gist_earth', 'terrain', 'ocean', 'gist_stern',
'brg', 'CMRmap', 'cubehelix',
'gnuplot', 'gnuplot2', 'gist_ncar',
'nipy_spectral', 'jet', 'rainbow',
'gist_rainbow', 'hsv', 'flag', 'prism'])]
In [ ]:
In [ ]: