Adding one or more attractors to a set of random 2D points
<p>Given a set of pseudorandom 2D points generated like so:</p>

<pre><code>points = random.sample([[x, y] for x in xrange(width) for y in yrange(height)], 100)
</code></pre>

<p>I would like to be able to add one or more non random attractor points around which the other points will be drawn. I don't plan to animate this so it doesn't need to be very efficient but I would like to be able to specify how the random points are drawn to each of the attractor points (for instance based on the square of the distance and a given gravitational constant) and then feed my points to a function that returns a modified version of the original list:</p>

<pre><code>points = random.sample([[x, y] for x in xrange(width) for y in xrange(height)], 100)

attractors = [(25, 102), (456, 300), (102, 562)]

def attract(random_points_list, attractor_points_list):
    (...)
    return modified_points_list

new_points_list = attract(points, attractors)
</code></pre>

<p>This new points list would then be used to seed a Voronoi diagram (not part of this question).</p>

<p>This turned out to be a more challenging problem than my initial estimate of it as a pure implementation-effort task. The difficulty is in defining a good attraction model.</p>

<ol>
<li><p>Simulating plain free-fall on attractors (as if in a real gravity field created by multiple point masses) is problematic since you must specify the duration of this process. If the duration is short enough, the displacements will be small and clustering around attractors will not be noticeable. If the duration is long enough then all points will fall on the attractors or be too close to them.</p></li>
<li><p>Computing the new position of each point in one shot (without doing time-based simulation) is simpler, but the question is whether the final position of every point must be affected by <em>all</em> attractors or <em>only the closest one</em> of them. The latter approach (attract-to-closest-one-only) proved to produce more visually appealing results. I could not achieve good results with the former approach (note however, that I tried only relatively simple attraction functions).</p></li>
</ol>

<p>Python 3.4 code with visualization using <code>matplotlib</code> follows:</p>

<pre><code>#!/usr/bin/env python3

import random
import numpy as np
import matplotlib.pyplot as plt

def dist(p1, p2):
    return np.linalg.norm(np.asfarray(p1) - np.asfarray(p2))


def closest_neighbor_index(p, attractors):
    min_d = float('inf')
    closest = None
    for i,a in enumerate(attractors):
        d = dist(p, a)
        if d &lt; min_d:
            closest, min_d = i, d
    return closest


def group_by_closest_neighbor(points, attractors):
    g = []
    for a in attractors:
        g.append([])
    for p in points:
        g[closest_neighbor_index(p, attractors)].append(p)
    return g


def attracted_point(p, a, f):
    a = np.asfarray(a)
    p = np.asfarray(p)
    r = p - a
    d = np.linalg.norm(r)
    new_d = f(d)
    assert(new_d &lt;= d)
    return a + r * new_d/d


def attracted_point_list(points, attractor, f):
    result=[]
    for p in points:
        result.append(attracted_point(p, attractor, f))
    return result


# Each point is attracted only to its closest attractor (as if the other
# attractors don't exist).
def attract_to_closest(points, attractors, f):
    redistributed_points = []
    grouped_points = group_by_closest_neighbor(points, attractors)
    for a,g in zip(attractors, grouped_points):
        redistributed_points.extend(attracted_point_list(g,a,f))
    return redistributed_points


def attraction_translation(p, a, f):
    return attracted_point(p, a, f) - p


# Each point is attracted by multiple attracters.
# The resulting point is the average of the would-be positions
# computed for each attractor as if the other attractors didn't exist.
def multiattract(points, attractors, f):
    redistributed_points = []
    n = float(len(attractors))
    for p in points:
        p = np.asfarray(p)
        t = np.zeros_like(p)
        for a in attractors:
            t += attraction_translation(p,a,f)
        redistributed_points.append(p+t/n)
    return redistributed_points


def attract(points, attractors, f):
    """ Draw points toward attractors

    points and attractors must be lists of points (2-tuples of the form (x, y)).

    f maps distance of the point from an attractor to the new distance value,
    i.e. for a single point P and attractor A, f(distance(P, A)) defines the
    distance of P from A in its new (attracted) location.

    0 &lt;= f(x) &lt;= x must hold for all non-negative values of x.
    """

    # multiattract() doesn't work well with simple attraction functions
    # return multiattract(points, attractors, f);
    return attract_to_closest(points, attractors, f);

if __name__ == '__main__':
    width=400
    height=300
    points = random.sample([[x, y] for x in range(width) for y in range(height)], 100)
    attractors = [(25, 102), (256, 256), (302, 62)]

    new_points = attract(points, attractors, lambda d: d*d/(d+100))
    #plt.scatter(*zip(*points), marker='+', s=32)
    plt.scatter(*zip(*new_points))
    plt.scatter(*zip(*attractors), color='red', marker='x', s=64, linewidths=2)

    plt.show()
</code></pre>

python game-physics voronoi