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:
Advertisement
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:
