Non overlapping random shapes on 2D plane

Objective: I want to create random non-overlapping irregular rounded shapes (contours) on a 2D plane, similar to a sandstone microstructure as shown below, to perform an oil flooding experiment using computer vision.

Approach: I have previously done a similar thing but with circular shapes instead of random shapes. The result is as below.

Also, I am aware of the programming way of creating closed Bezier curves.

Help: But I am not able to combine both these steps into a single script. Also if there are any other alternatives it is more than welcome.

Code:

1. Circular Shapes inside 2D plane

2. Bezier Curves

Answer

I took most of the code from ImportanceOfBeingErnest’s answer here: Create random shape/contour using matplotlib. Some fine-tuning is needed to tighten the space between the random shapes but the basic idea is there

import numpy as np
import random
import math
import matplotlib.pyplot as plt
from scipy.special import binom

def dist(x1, y1, x2, y2):
  return np.sqrt((x1-x2)**2 + (y1-y2)**2)

bernstein = lambda n, k, t: binom(n,k)* t**k * (1.-t)**(n-k)

def bezier(points, num=100):
    N = len(points)
    t = np.linspace(0, 1, num=num)
    curve = np.zeros((num, 2))
    for i in range(N):
        curve += np.outer(bernstein(N - 1, i, t), points[i])
    return curve

class Segment():
    def __init__(self, p1, p2, angle1, angle2, **kw):
        self.p1 = p1; self.p2 = p2
        self.angle1 = angle1; self.angle2 = angle2
        self.numpoints = kw.get("numpoints", 100)
        r = kw.get("r", 4)
        d = np.sqrt(np.sum((self.p2-self.p1)**2))
        self.r = r*d
        self.p = np.zeros((4,2))
        self.p[0,:] = self.p1[:]
        self.p[3,:] = self.p2[:]
        self.calc_intermediate_points(self.r)

    def calc_intermediate_points(self,r):
        self.p[1,:] = self.p1 + np.array([self.r*np.cos(self.angle1),
                                    self.r*np.sin(self.angle1)])
        self.p[2,:] = self.p2 + np.array([self.r*np.cos(self.angle2+np.pi),
                                    self.r*np.sin(self.angle2+np.pi)])
        self.curve = bezier(self.p,self.numpoints)


def get_curve(points, **kw):
    segments = []
    for i in range(len(points)-1):
        seg = Segment(points[i,:2], points[i+1,:2], points[i,2],points[i+1,2],**kw)
        segments.append(seg)
    curve = np.concatenate([s.curve for s in segments])
    return segments, curve

def ccw_sort(p):
    d = p-np.mean(p,axis=0)
    s = np.arctan2(d[:,0], d[:,1])
    return p[np.argsort(s),:]

def get_bezier_curve(a, rad=2, edgy=0):
    """ given an array of points *a*, create a curve through
    those points. 
    *rad* is a number between 0 and 1 to steer the distance of
          control points.
    *edgy* is a parameter which controls how "edgy" the curve is,
           edgy=0 is smoothest."""
    p = np.arctan(edgy)/np.pi+.5
    a = ccw_sort(a)
    a = np.append(a, np.atleast_2d(a[0,:]), axis=0)
    d = np.diff(a, axis=0)
    ang = np.arctan2(d[:,1],d[:,0])
    f = lambda ang : (ang>=0)*ang + (ang<0)*(ang+2*np.pi)
    ang = f(ang)
    ang1 = ang
    ang2 = np.roll(ang,1)
    ang = p*ang1 + (1-p)*ang2 + (np.abs(ang2-ang1) > np.pi )*np.pi
    ang = np.append(ang, [ang[0]])
    a = np.append(a, np.atleast_2d(ang).T, axis=1)
    s, c = get_curve(a, r=rad, method="var")
    x,y = c.T
    return x,y, a


def get_random_points(n=3, scale=2, mindst=None, rec=0):
    """ create n random points in the unit square, which are *mindst*
    apart, then scale them."""
    mindst = mindst or .7/n
    a = np.random.rand(n,2)
    d = np.sqrt(np.sum(np.diff(ccw_sort(a), axis=0), axis=1)**2)
    if np.all(d >= mindst) or rec>=200:
        return a*scale
    else:
        return get_random_points(n=n, scale=scale, mindst=mindst, rec=rec+1)
    
x_list = [] #list of x coordinates of circular inclusions
y_list = [] #list of y coordinates of circular inclusions
r_list = [] #list of radii of circular inclusions
n = 0       #number of circular inclusions
rad = 0.3
edgy = 0.05
fig, ax = plt.subplots(figsize=(8, 4))
ax.set(xlim=(-20, 20), ylim = (-10, 10)); #size of rectangular space (40 X 20) can be varied

while n < 80: #choosing number of inclusions (50), choose such a number so it fits inside the rectangular space
  x = round(random.uniform(-20, 20), 2) #generating random x coordinate
  y = round(random.uniform(-10, 10), 2) #generating random y coordinate

  overlap = False #checking and removing overlapping inclusions
  for j in range(n):
    d = dist(x, y, x_list[j], y_list[j])
    if d < 1.5:
      overlap = True

  if overlap == False:    
    x_list.append(x)
    y_list.append(y)
    n += 1
for (x,y) in zip(x_list, y_list):
  a = get_random_points(n=4, scale=2) + [x, y]
  x,y, _ = get_bezier_curve(a,rad=rad, edgy=edgy)
  plt.plot(x,y)
  
plt.show()

A random generated result: