# Drawing a cloud of points

I decided that it would be nice to write some code to draw a cloud of points.

Of course my first thought was to draw some random points but the result is not cloud-ish, it is well... random.

In the end I decided to use the Mitchell random points algoritmh generating pretty good results:

• Blue = Random

• Red = Cloud (Michelin)

from __future__ import division
import matplotlib.pyplot as plt
import random

# The more the nearest the points
RANDOM_POINTS_PER_GENERATION = 6

def random_point():
return [random.random(), random.random()]

def average(lst):
return sum(lst) / len(lst)

def distance(a,b):
return ((a[0]-b[0])**2 + (a[1]-b[1])**2)

def nearest(random_points, seen_points):
def distance_to_points(point):
return average([distance(point, seen_point)
for seen_point in seen_points])
return min(random_points, key=distance_to_points)

def next_point(points):
random_points = [random_point() for _ in range(RANDOM_POINTS_PER_GENERATION)]
return nearest(random_points, points)

def mitchell(number_of_points):
"""
Generates almost random points, but at each iteration
k points are generated and the nearest is taken, in
order to have a better looking points-cloud.
"""
points = [random_point()]
while len(points) < number_of_points:
points.append(next_point(points))
return points

def plot_points(points, color):
plt.plot(*zip(*points), color=color, linestyle='none', marker='.',
markerfacecolor='blue', markersize=4)

if __name__ == "__main__":
plot_points(mitchell(250), "red")
plot_points([random_point() for _ in range(250)], "blue")
plt.show()


1. mitchell is the only function with a docstring. What do all the other functions do? How do I call them? What do they return?

2. A function's name should describe what it does, or what it returns, rather than the name of algorithm it uses (which is typically an implementation detail). So I would use a name like point_cloud instead of mitchell.

3. Your implementation picks the random sample that is closest to all the other points, but the algorithm described in Mitchell's 1991 paper "Spectrally optimal sampling for distribution ray tracing" picks the random sample that is furthest away from any other point (see section 4, "Sequential Poisson-Disk Sampling"). This is so different that I think it would be better not to use the name "Mitchell" to describe the algorithm.

4. It would make the code easier to use if RANDOM_POINTS_PER_GENERATION were a keyword argument to mitchell, rather than a global variable. Consider a program that tries to determine a good value for this parameter, or a program with two pieces of code that need different values for the parameter.

5. mitchell(0) returns a list of one point, but surely this should return a list of no points?

while len(points) < number_of_points:


write:

for _ in range(number_of_points - 1):


since you know that each iteration adds exactly one point.

7. The function average is built into Python under the name statistics.mean.

8. The function distance is misleadly named: it actually calculates the squared distance between points. Similarly, distance_to_points returns the mean square distance to the points.

9. The function nearest doesn't need to take the average: since seen_point has the same length in every case, the sum of the squared distances would be just as good.

10. Instead of creating a temporary list that gets thrown away immediately:

average([distance(point, seen_point) for seen_point in seen_points])


use a generator expression:

sum(distance(point, seen_point) for seen_point in seen_points)

11. At each stage nearest finds the average squared distance from a candidate point to each point in the cloud, making the overall runtime $O(kn^2)$. However, it's possible to do better than that. For each candidate point $q$, the value that needs to be computed is the sum of squared distances, $$\sum_j \left|p_j - q\right|^2.$$ Using the dot product we can expand the square, getting $$\sum_j (p_j - q)·(p_j - q).$$ Multiplying out, we get $$\sum_j p_j·p_j - 2q·\sum_j p_j + q.q \sum_j 1.$$ Since we taking the minimum of this value over all the candidate points $q$, we don't need the first term $\sum_j p_j·p_j$ as this is the same for all $q$ and so cannot affect the result.

What this means is that if we keep a running sum $\sum_j p_j$ then we can find the candidate point with the minimum sum of squared distances to all the points in the cloud in $O(k)$, bringing the runtime down to $O(kn)$.

To implement this, we need a representation of vectors that supports sums and dot products. Here I'm using NumPy arrays:

import numpy as np

def point_cloud(n, k=6):
"""Generate a point cloud with n points. At each step k (default 6)
candidate points are generated and the one whose sum of squared
distances to the other points is smallest is used.

"""
if n == 0: return
p = np.random.random(size=2)
yield p
sum_p = p.copy()
for i in range(1, n):
# Generate k candidate points.
q = np.random.random(size=(k, 2))

# Sum of squared distances from each candidate point to all
# the points generated so far, less the sum of squares of the
# points generated so far (as this doesn't affect the minimum,
# we don't need it here).
s = i * (q * q).sum(axis=1) - 2 * (q * sum_p).sum(axis=1)

# Pick the candidate that's closest to the points generated so
# far.
p = q[np.argmin(s)]
yield p
sum_p += p