2

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.

Sandstone Microstructure

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

Circular Shapes

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

Christoph Rackwitz
  • 11,317
  • 4
  • 27
  • 36
Alpha
  • 399
  • 3
  • 9
  • Can you please share the code that you're using / trying to combine into one script. This will help to answer your question. – gincard Jun 12 '22 at 09:04
  • Naive approach could be: 1. Create ransom bezier patch. 2. Test whether it doesnt intersect with any other created patch. 3. Repeat. This will be much less efficient than circles, because you have to test after generation, while for circles you can define a working radius after testing the distances. – Micka Jun 12 '22 at 09:58
  • throw shapes on canvas, perform **collision detection**, and shove shapes around. your shapes are convex so this seems reasonably well-behaved. you could even do that in raster space: consider centroid of intersection relative to each shape's center, which gives you a vector to try moving. – Christoph Rackwitz Jun 12 '22 at 14:33
  • Can someone provide a code for throwing these generated images onto the canvas in python? I am able to do that but it is coming out to be quite uniform and not random – Alpha Jun 12 '22 at 19:15
  • Note that your post doesn't actually contain a question, though, so... what [on-topic thing](/help/on-topic) do you need _help_ with? There's a bunch of off-the-shelf physics engines you can use here, so that you don't have to implement the parts you're going to fail at anyway (through no fault of your own: physics engines already account for all the edge cases you're going to run into again and again and again if you try to roll your own). – Mike 'Pomax' Kamermans Jun 13 '22 at 00:15

1 Answers1

1

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: enter image description here

Mike B
  • 2,136
  • 2
  • 12
  • 31