Controlled circle packing with Python

Question

I'm trying to port to Python the "Controlled Circle Packing with Processing" algorithm that I found here:

http://www.codeplastic.com/2017/09/09/controlled-circle-packing-with-processing/?replytocom=22#respond

For now my goal is just to make it work, before I tweak it for my own needs. This question is not about the best way to do circle packing.

So far here is what I have:

#!/usr/bin/python
# coding: utf-8

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


class Ball:

    def __init__(self, x, y, radius):

        self.r = radius

        self.acceleration = np.array([0, 0])

        self.velocity = np.array([uniform(0, 1),
                                  uniform(0, 1)])

        self.position = np.array([x, y])


    @property
    def x(self):
        return self.position[0]

    @property
    def y(self):
        return self.position[1]


    def applyForce(self, force):

        self.acceleration = np.add(self.acceleration, force)


    def update(self):

        self.velocity = np.add(self.velocity, self.acceleration)
        self.position = np.add(self.position, self.velocity)
        self.acceleration *= 0


class Pack:

    def __init__(self, radius, list_balls):

        self.list_balls = list_balls
        self.r = radius
        self.list_separate_forces = [np.array([0, 0])] * len(self.list_balls)
        self.list_near_balls = [0] * len(self.list_balls)


    def _normalize(self, v):

        norm = np.linalg.norm(v)
        if norm == 0:
            return v
        return v / norm


    def run(self):

        for i in range(300):
            print(i)
            for ball in self.list_balls:
                self.checkBorders(ball)
                self.checkBallPositions(ball)
                self.applySeparationForcesToBall(ball)

    def checkBorders(self, ball):

        if (ball.x - ball.r) < - self.r or (ball.x + ball.r) > self.r:
            ball.velocity[0] *= -1
            ball.update()
        if (ball.y - ball.r) < -self.r or (ball.y + ball.r) > self.r:
            ball.velocity[1] *= -1
            ball.update()


    def checkBallPositions(self, ball):

        list_neighbours = [e for e in self.list_balls if e is not ball]

        for neighbour in list_neighbours:

            d = self._distanceBalls(ball, neighbour)

            if d < (ball.r + neighbour.r):
                return

        ball.velocity[0] = 0
        ball.velocity[1] = 0


    def getSeparationForce(self, c1, c2):

        steer = np.array([0, 0])

        d = self._distanceBalls(c1, c2)

        if d > 0 and d < (c1.r + c2.r):
            diff = np.subtract(c1.position, c2.position)
            diff = self._normalize(diff)
            diff = np.divide(diff, d)
            steer = np.add(steer, diff)

        return steer


    def _distanceBalls(self, c1, c2):

        x1, y1 = c1.x, c1.y
        x2, y2 = c2.x, c2.y

        dist = np.sqrt((x2 - x1)**2 + (y2 - y1)**2)

        return dist


    def applySeparationForcesToBall(self, ball):

        i = self.list_balls.index(ball)

        list_neighbours = [e for e in self.list_balls if e is not ball]

        for neighbour in list_neighbours:
            j = self.list_balls.index(neighbour)
            forceij = self.getSeparationForce(ball, neighbour)

            if np.linalg.norm(forceij) > 0:
                self.list_separate_forces[i] = np.add(self.list_separate_forces[i], forceij)
                self.list_separate_forces[j] = np.subtract(self.list_separate_forces[j], forceij)
                self.list_near_balls[i] += 1
                self.list_near_balls[j] += 1

        if self.list_near_balls[i] > 0:
            self.list_separate_forces[i] = np.divide(self.list_separate_forces[i], self.list_near_balls[i])


        if np.linalg.norm(self.list_separate_forces[i]) > 0:
            self.list_separate_forces[i] = self._normalize(self.list_separate_forces[i])
            self.list_separate_forces[i] = np.subtract(self.list_separate_forces[i], ball.velocity)
            self.list_separate_forces[i] = np.clip(self.list_separate_forces[i], a_min=0, a_max=np.array([1]))

        separation = self.list_separate_forces[i]
        ball.applyForce(separation)
        ball.update()


list_balls = list()

for i in range(10):
    b = Ball(0, 0, 7)
    list_balls.append(b)


p = Pack(30, list_balls)
p.run()

plt.axes()

# Big container
circle = plt.Circle((0, 0), radius=30, fc='none', ec='k')
plt.gca().add_patch(circle)

for c in list_balls:
    ball = plt.Circle((c.x, c.y), radius=c.r, picker=True, fc='none', ec='k')
    plt.gca().add_patch(ball)

plt.axis('scaled')
plt.show()

The code was originally written with Processing, I did my best to use numpy instead. I'm not quite sure of my checkBallPosition, the original author uses a count variable that looks useless to me. I also wonder why the steer vector in the original code has a dimension of 3.

So far, here is what my code yields:

The circles (I had to rename them balls to not conflict with Circle from matplotlib) overlap and don't seem to get away from each other. I don't think I'm really far but I would need a bit of help to find what's wrong with my code. Could you help me please ?

EDIT: I realize that I probably need to do several passes. Maybe the Processing package (language ?) runs the run function several times. It actually makes sense to me, this problem is very similar to molecular mechanics optimization and it's an iterative process.

My question can now be a bit more specific: it seems the checkBorders function doesn't do its job properly and doesn't rebound the circles properly. But given its simplicity, I'd say the bug is in applySeparationForcesToBall, I probably don't apply the forces correctly.

Please ask a more concrete question. I see that you took some time to write a good post, but the actual question in the end is just asking for help debugging your code. Try to identify which part does not work as you expect it and make this a concrete questions (e.g. "is it correct when I use the 2-norm in X to calculate Y") which can have a concrete answer. — allo, Mar 26 '18 at 15:19
I haven't tried to understand your code, but it looks strange to me that you divide `r` (which is a radius) by 2 in some places. — Arndt Jonasson, Mar 26 '18 at 15:55
Besides, if two circles are created overlapping, wouldn't they just keep bouncing to and fro? — Arndt Jonasson, Mar 26 '18 at 15:55
the ```r/2``` seems strange to me as well. That's one thing I don't understand. As for your 2nd comment, I don't think so, if they overlap, ```checkBallPositions``` doesn't zero their velocities (which are randomly initialized), so will un-overlap at some point. — JPFrancoia, Mar 26 '18 at 16:00
Have you tried debugging your code? Which line of code is behaving differently from what you expected? — Kevin Workman, Mar 26 '18 at 16:22
Yes, of course I have tried. But I can't really find a line that behaves "weirdly", even if I'm sure there is at least one. I updated my question to make it a bit clearer. I realized I need to make several iterations. The problem probably comes from ```applySeparationForcesToBall```, now the circles get ejected from inside the container. — JPFrancoia, Mar 26 '18 at 16:42
The Processing code you're talking about is an animation that updates 60 times per second, and doesn't contain any logic for staying inside a rectangle. You'll have better luck if you just implement the algorithm from scratch instead of trying to hammer an existing program to do something else. I'm going to remove the [tag:processing] tag because this does not contain a question about Processing. — Kevin Workman, Mar 26 '18 at 22:18

score 9 · Accepted Answer · answered Mar 28 '18 at 10:17

Ok after days of fiddling, I managed to do it:

Here is the complete code:

#!/usr/bin/python
# coding: utf-8

"""
http://www.codeplastic.com/2017/09/09/controlled-circle-packing-with-processing/
https://stackoverflow.com/questions/573084/how-to-calculate-bounce-angle/573206#573206
https://stackoverflow.com/questions/4613345/python-pygame-ball-collision-with-interior-of-circle
"""

import numpy as np
import matplotlib.pyplot as plt
from random import randint
from random import uniform
from matplotlib import animation


class Ball:

    def __init__(self, x, y, radius):

        self.r = radius

        self.acceleration = np.array([0, 0])

        self.velocity = np.array([uniform(0, 1),
                                  uniform(0, 1)])

        self.position = np.array([x, y])


    @property
    def x(self):
        return self.position[0]

    @property
    def y(self):
        return self.position[1]


    def applyForce(self, force):

        self.acceleration = np.add(self.acceleration, force)

    def _normalize(self, v):

        norm = np.linalg.norm(v)
        if norm == 0:
            return v
        return v / norm


    def update(self):

        self.velocity = np.add(self.velocity, self.acceleration)
        self.position = np.add(self.position, self.velocity)
        self.acceleration *= 0


class Pack:

    def __init__(self, radius, list_balls):

        self.iter = 0
        self.list_balls = list_balls
        self.r = radius
        self.list_separate_forces = [np.array([0, 0])] * len(self.list_balls)
        self.list_near_balls = [0] * len(self.list_balls)
        self.wait = True


    def _normalize(self, v):

        norm = np.linalg.norm(v)
        if norm == 0:
            return v
        return v / norm


    def run(self):

        self.iter += 1
        for ball in self.list_balls:
            self.checkBorders(ball)
            self.checkBallPositions(ball)
            self.applySeparationForcesToBall(ball)
            print(ball.position)

        print("\n")


    def checkBorders(self, ball):

        d = np.sqrt(ball.x**2 + ball.y**2)

        if d >= self.r - ball.r:

            vr = self._normalize(ball.velocity) * ball.r

            # P1 is collision point between circle and container
            P1x = ball.x + vr[0]
            P1y = ball.y + vr[1]
            P1 = np.array([P1x, P1y])

            # Normal vector
            n_v = -1 * self._normalize(P1)

            u = np.dot(ball.velocity, n_v) * n_v
            w = np.subtract(ball.velocity, u)

            ball.velocity = np.subtract(w, u)

            ball.update()



    def checkBallPositions(self, ball):

        i = self.list_balls.index(ball)

        # for neighbour in list_neighbours:
        # ot a full loop; if we had two full loops, we'd compare every
        # particle to every other particle twice over (and compare each
        # particle to itself)
        for neighbour in self.list_balls[i + 1:]:

            d = self._distanceBalls(ball, neighbour)

            if d < (ball.r + neighbour.r):
                return

        ball.velocity[0] = 0
        ball.velocity[1] = 0


    def getSeparationForce(self, c1, c2):

        steer = np.array([0, 0])

        d = self._distanceBalls(c1, c2)

        if d > 0 and d < (c1.r + c2.r):
            diff = np.subtract(c1.position, c2.position)
            diff = self._normalize(diff)
            diff = np.divide(diff, 1 / d**2)
            steer = np.add(steer, diff)

        return steer


    def _distanceBalls(self, c1, c2):

        x1, y1 = c1.x, c1.y
        x2, y2 = c2.x, c2.y

        dist = np.sqrt((x2 - x1)**2 + (y2 - y1)**2)

        return dist


    def applySeparationForcesToBall(self, ball):

        i = self.list_balls.index(ball)

        for neighbour in self.list_balls[i + 1:]:
            j = self.list_balls.index(neighbour)
            forceij = self.getSeparationForce(ball, neighbour)

            if np.linalg.norm(forceij) > 0:
                self.list_separate_forces[i] = np.add(self.list_separate_forces[i], forceij)
                self.list_separate_forces[j] = np.subtract(self.list_separate_forces[j], forceij)
                self.list_near_balls[i] += 1
                self.list_near_balls[j] += 1

        if np.linalg.norm(self.list_separate_forces[i]) > 0:
            self.list_separate_forces[i] = np.subtract(self.list_separate_forces[i], ball.velocity)

        if self.list_near_balls[i] > 0:
            self.list_separate_forces[i] = np.divide(self.list_separate_forces[i], self.list_near_balls[i])


        separation = self.list_separate_forces[i]
        ball.applyForce(separation)
        ball.update()


list_balls = list()

for i in range(25):
    # b = Ball(randint(-15, 15), randint(-15, 15), 5)
    b = Ball(0, 0, 5)
    list_balls.append(b)


p = Pack(30, list_balls)

fig = plt.figure()

circle = plt.Circle((0, 0), radius=30, fc='none', ec='k')
plt.gca().add_patch(circle)
plt.axis('scaled')
plt.axes().set_xlim(-50, 50)
plt.axes().set_ylim(-50, 50)


def draw(i):

    patches = []

    p.run()
    fig.clf()
    circle = plt.Circle((0, 0), radius=30, fc='none', ec='k')
    plt.gca().add_patch(circle)
    plt.axis('scaled')
    plt.axes().set_xlim(-50, 50)
    plt.axes().set_ylim(-50, 50)

    for c in list_balls:
        ball = plt.Circle((c.x, c.y), radius=c.r, picker=True, fc='none', ec='k')
        patches.append(plt.gca().add_patch(ball))

    return patches


co = False
anim = animation.FuncAnimation(fig, draw,
                               frames=500, interval=2, blit=True)


# plt.show()


anim.save('line2.gif', dpi=80, writer='imagemagick')

From the original code, I modified the checkBorder function to bounce the circles properly from the edge, and changed the separation force between circles, it was too low. I know my question was a bit too vague from start, but I would have appreciated a more constructive feedback.

Controlled circle packing with Python

1 Answers1