Find new position for overlapping circles

Question

I am trying to write a code that, for a given list of circles (list1), it is able to find the positions for new circles (list2). list1 and list2 have the same length, because for each circle in list1 there must be a circle from list2.

1. Each pair of circles (let's say circle1 from list1 and circle2 from list2), must be as close together as possible,
2. circles from list2 must not overlap with circles from list1, while circles of the single lists can overlap each other.

list1 is fixed, so now I have to find the right position for circles from list2.

I wrote this simple function to recognize if 2 circles overlap:

def overlap(x1, y1, x2, y2, r1, r2):
    distSq = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)
    radSumSq = (r1 + r2) * (r1 + r2)
    if (distSq >= radSumSq):
        return False # no overlap
    else:
        return True  #overlap

and this is the list1:

with:

x=[14.11450195 14.14184093 14.15435028 14.16206741 14.16951752 14.17171097
 14.18569565 14.19700241 14.23129082 14.24083233 14.24290752 14.24968338
 14.2518959  14.26536751 14.27209759 14.27612877 14.2904377  14.29187012
 14.29409599 14.29618549 14.30615044 14.31624985 14.3206892  14.3228569
 14.36143875 14.36351967 14.36470699 14.36697292 14.37235737 14.41422081
 14.42583466 14.43226814 14.43319225 14.4437027  14.4557848  14.46592999
 14.47036076 14.47452068 14.47815609 14.52229309 14.53059006 14.53404236
 14.5411644 ] 
y=[-0.35319126 -0.44222349 -0.44763246 -0.35669261 -0.24366629 -0.3998799
 -0.38940558 -0.57744932 -0.45223859 -0.21021004 -0.44250247 -0.45866323
 -0.47203487 -0.51684451 -0.44884869 -0.2018993  -0.40296811 -0.23641759
 -0.18019417 -0.33391538 -0.53565156 -0.45215255 -0.40939832 -0.26936951
 -0.30894437 -0.55504167 -0.47177047 -0.45573688 -0.43100587 -0.5805912
 -0.21770373 -0.199422   -0.17372169 -0.38522363 -0.56950212 -0.56947368
 -0.48770753 -0.24940367 -0.31492445 -0.54263926 -0.53460872 -0.4053807
 -0.43733299]
radius = 0.014

---

Copy and pasteable...

x = [14.11450195,14.14184093,14.15435028,14.16206741,14.16951752,
     14.17171097,14.18569565,14.19700241,14.23129082,14.24083233,
     14.24290752,14.24968338,14.2518959,14.26536751,14.27209759,
     14.27612877,14.2904377,14.29187012,14.29409599,14.29618549,
     14.30615044,14.31624985,14.3206892,14.3228569,14.36143875,
     14.36351967,14.36470699,14.36697292,14.37235737,14.41422081,
     14.42583466,14.43226814,14.43319225,14.4437027,14.4557848,
     14.46592999,14.47036076,14.47452068,14.47815609,14.52229309,
     14.53059006,14.53404236,14.5411644]

y = [-0.35319126,-0.44222349,-0.44763246,-0.35669261,-0.24366629,
     -0.3998799,-0.38940558,-0.57744932,-0.45223859,-0.21021004,
     -0.44250247,-0.45866323,-0.47203487,-0.51684451,-0.44884869,
     -0.2018993,-0.40296811,-0.23641759,-0.18019417,-0.33391538,
     -0.53565156,-0.45215255,-0.40939832,-0.26936951,-0.30894437,
     -0.55504167,-0.47177047,-0.45573688,-0.43100587,-0.5805912,
     -0.21770373,-0.199422,-0.17372169,-0.38522363,-0.56950212,
     -0.56947368,-0.48770753,-0.24940367,-0.31492445,-0.54263926,
     -0.53460872,-0.4053807,-0.43733299]

---

Now I am not sure about what I have to do, my first idea is to draw circles of list2 taking x and y from list one and do something like x+c and y+c, where c is a fixed value. Then I can call my overlapping function and, if there is overlap I can increase the c value. In this way I have 2 for loops. Now, my questions are:

1. There is a way to avoid for loops?
2. Is there a smart solution to find a neighbor (circle from list2) for each circle from list1 (without overlaps with other circles from list2)?

Answer 1

Using numpy arrays, you can avoid for loops.

Setup from your example.

import numpy as np

#Using your x and y
c1 = np.array([x,y]).T

# random set of other centers within the same range as c1
c2 = np.random.random((10,2))
np.multiply(c2, c1.max(0)-c1.min(0),out = c2)
np.add(c2, c1.min(0), out=c2)

radius = 0.014
r = radius
min_d = (2*r)*(2*r)

plot_circles(c1,c2)    # see function at end

An array of distances from each center in c1 to each center in c2

def dist(c1,c2):
    dx = c1[:,0,None] - c2[:,0]
    dy = c1[:,1,None] - c2[:,1]
    return dx*dx + dy*dy

d = dist(c1,c2)

Or you could use scipy.spatial

from scipy.spatial import distance
d = distance.cdist(c1,c2,'sqeuclidean')

Create a 2d Boolean array for circles that intersect.

intersect = d <= min_d

Find the indices of overlapping circles from the two sets.

a,b = np.where(intersect)
plot_circles(c1[a],c2[b])

---

Using intersect or a and b to index c1,c2, and d you should be able to get groups of intersecting circles then figure out how to move the c2 centers - but I'll leave that for another question/answer. If a list2 circle intersects one list1 circle - find the line between the two and move along that line. If a list2 circle intersects more than one list1 circle - find the line between the two closestlist1circles and move thelitst2` circle along a line perpendicular to that. You didn't mention any constraints on moving the circles so maybe random movement then find the intersects again but that might be problematic. In the following image, it may be trivial to figure out how to move most of the red circles but the group circled in blue might require a different strategy.

Here are some examples for getting groups:

>>> for f,g,h in zip(c1[a],c2[b],d[a,b]):
    print(f,g,h)
>>> c1[intersect.any(1)],c2[intersect.any(0)]
>>> for (f,g) in zip(c2,intersect.T):
    if g.any():
            print(f.tolist(),c1[g].tolist())

---

import matplotlib as mpl
from matplotlib import pyplot as plt

def plot_circles(c1,c2):

    bounds = np.array([c1.min(0),c2.min(0),c1.max(0),c2.max(0)])
    xmin, ymin = bounds.min(0)
    xmax, ymax = bounds.max(0)

    circles1 = [mpl.patches.Circle(xy,radius=r,fill=False,edgecolor='g') for xy in c1]
    circles2 = [mpl.patches.Circle(xy,radius=r,fill=False,edgecolor='r') for xy in c2]
    fig = plt.figure()
    ax = fig.add_subplot(111)
    for c in circles2:
        ax.add_artist(c)
    for c in circles1:
        ax.add_artist(c)

    ax.set_xlim(xmin-r,xmax+r)
    ax.set_ylim(ymin-r,ymax+r)

    plt.show()
    plt.close()

Answer 2

This problem can very well be seen as an optimization problem. To be more precise, a nonlinear optimization problem with constraints. Since optimization strategies are not always so easy to understand, I will define the problem as simply as possible and also choose an approach that is as general as possible (but less efficient) and does not involve a lot of mathematics. As a spoiler: We are going to formulate the problem and the minimization process in less than 10 lines of code using the scipy library.

However, I will still provide hints on where you can get your hands even dirtier.

Formulating the problem

As a guide for a formulation of an NLP-class problem (Nonlinear Programming), you can go directly to the two requirements in the original post.

1. Each pair of circles must be as close together as possible -> Hint for a cost-function
2. Circles must not overlap with other (moved) circles -> Hint for a constraint

Cost function

Let's start with the formulation of the cost function to be minimized. Since the circles should be moved as little as possible (resulting in the closest possible neighborhood), a quadratic penalty term for the distances between the circles of the two lists can be chosen for the cost function:

import scipy.spatial.distance as sd

def cost_function(new_positions, old_positions):
    new_positions = np.reshape(new_positions, (-1, 2))
    return np.trace(sd.cdist(new_positions, old_positions, metric='sqeuclidean'))

Why quadratic? Partly because of differentiability and for stochastic reasons (think of the circles as normally distributed measurement errors -> least squares is then a maximum likelihood estimator). By exploiting the structure of the cost function, the efficiency of the optimization can be increased (elimination of sqrt). By the way, this problem is related to nonlinear regression, where (nonlinear) least squares are also used.

Now that we have a cost function at hand, we also have a good way to evaluate our optimization. To be able to compare solutions of different optimization strategies, we simply pass the newly calculated positions to the cost function.

Let's give it a try: For example, let us use the calculated positions from the Voronoi approach (by Paul Brodersen).

print(cost_function(new_positions, old_positions))
# prints 0.007999244511697411

That's a pretty good value if you ask me. Considering that the cost function spits out zero when there is no displacement at all, this cost is pretty close. We can now try to outperform this value by using classical optimization!

Non-linear constraint

We know that circles must not overlap with other circles in the new set. If we translate this into a constraint, we find that the lower bound for the distance is 2 times the radius and the upper bound is simply infinity.

import scipy.spatial.distance as sd
from scipy.optimize import NonlinearConstraint

def cons_f(x):
    x = np.reshape(x, (-1, 2))
    return sd.pdist(x)

nonlinear_constraint = NonlinearConstraint(cons_f, 2*radius, np.inf, jac='2-point')

Here we make life easy by approximating the Jacobi matrix via finite differences (see parameter jac='2-point'). At this point it should be said that we can increase the efficiency here, by formulating the derivatives of the first and second order ourselves instead of using approximations. But this is left to the interested reader. (It is not that hard, because we use quite simple mathematical expressions for distance calculation here.)

One additional note: You can also set a boundary constraint for the positions themselves not to exceed a specified region. This can then be used as another parameter. (See scipy.optimize.Bounds)

Minimizing the cost function under constraints

Now we have both ingredients, the cost function and the constraint, in place. Now let's minimize the whole thing!

from scipy.optimize import minimize

res = minimize(lambda x: cost_function(x, positions), positions.flatten(), method='SLSQP',
               jac="2-point", constraints=[nonlinear_constraint])

As you can see, we approximate the first derivatives here as well. You can also go deeper here and set up the derivatives yourself (analytically).

Also note that we must always pass the parameters (an nx2 vector specifying the positions of the new layout for n circles) as a flat vector. For this reason, reshaping can be found several times in the code.

Evaluation, summary and visualization

Let's see how the optimization result performs in our cost function:

new_positions = np.reshape(res.x, (-1,2))
print(cost_function(new_positions, old_positions))
# prints 0.0010314079483565686

Starting from the Voronoi approach, we actually reduced the cost by another 87%! Thanks to the power of modern optimization strategies, we can solve a lot of problems in no time.

Of course, it would be interesting to see how the shifted circles look now: Circles after Optimization

Performance: 77.1 ms ± 1.17 ms

The entire code:

from scipy.optimize import minimize
import scipy.spatial.distance as sd
from scipy.optimize import NonlinearConstraint

# Given by original post
positions = np.array([x, y]).T

def cost_function(new_positions, old_positions):
    new_positions = np.reshape(new_positions, (-1, 2))
    return np.trace(sd.cdist(new_positions, old_positions, metric='sqeuclidean'))

def cons_f(x):
    x = np.reshape(x, (-1, 2))
    return sd.pdist(x)

nonlinear_constraint = NonlinearConstraint(cons_f, 2*radius, np.inf, jac='2-point')
res = minimize(lambda x: cost_function(x, positions), positions.flatten(), method='SLSQP',
               jac="2-point", constraints=[nonlinear_constraint])