Based on Paul Bourke's article.

In summary, an intuitive approach to tiling (filling a space with no overlaps) is to find a series of shapes with smaller areas that eventually add up to the area of the canvas and couple it with a collision detection system. One way to systematically find these areas is using the Riemann zeta function. Since the tiling is made up of circles, simply comparing the distance with the radii is sufficient for detecting collisions.

Written in Javascript with Rough.js

code