pythonalgorithmmatplotlibgeometryconvex-hull
Polygon from set of points that lie on a boundary
I have following set of points that lie on a boundary and want to create the polygon that connects these points. For a person it is quite obvious what path to follow, but I am unable to find an algorithm that does the same and trying to solve it myself it all seems quite tricky and ambiguous occasionally. What is the best solution for this?
As a background.
This is the boundary for the julia set with constant = -0.624+0.435j with stable area defined after 100 iterations. I got these points by setting the stable points to 1 and all other to zero and then convolving with a 3x3 matrix [[1, 1, 1], [1, 1, 1], [1, 1, 1]] and select the points that have value 1. My experimenting code is as follows:
import numpy as np
from scipy.signal import convolve2d
import matplotlib.pyplot as plt
r_min, r_max = -1.5, 1.5
c_min, c_max = -2.0, 2.0
dpu = 50 # dots per unit - 50 dots per 1 units means 200 points per 4 units
max_iterations = 100
cmap='hot'
intval = 1 / dpu
r_range = np.arange(r_min, r_max + intval, intval)
c_range = np.arange(c_min, c_max + intval, intval)
constant = -0.624+0.435j
def z_func(point, constant):
z = point
stable = True
num_iterations = 1
while stable and num_iterations < max_iterations:
z = z**2 + constant
if abs(z) > max(abs(constant), 2):
stable = False
return (stable, num_iterations)
num_iterations += 1
return (stable, 0)
points = np.array([])
colors = np.array([])
stables = np.array([], dtype='bool')
progress = 0
for imag in c_range:
for real in r_range:
point = complex(real, imag)
points = np.append(points, point)
stable, color = z_func(point, constant)
stables = np.append(stables, stable)
colors = np.append(colors, color)
print(f'{100*progress/len(c_range)/len(r_range):3.2f}% completed\r', end='')
progress += len(r_range)
print(' \r', end='')
rows = len(r_range)
start = len(colors)
orig_field = []
for i_num in range(len(c_range)):
start -= rows
real_vals = [color for color in colors[start:start+rows]]
orig_field.append(real_vals)
orig_field = np.array(orig_field, dtype='int')
rows = len(r_range)
start = len(stables)
stable_field = []
for i_num in range(len(c_range)):
start -= rows
real_vals = [1 if val == True else 0 for val in stables[start:start+rows]]
stable_field.append(real_vals)
stable_field = np.array(stable_field, dtype='int')
kernel = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
stable_boundary = convolve2d(stable_field, kernel, mode='same')
boundary_points = []
cols, rows = stable_boundary.shape
assert cols == len(c_range), "check c_range and cols"
assert rows == len(r_range), "check r_range and rows"
zero_field = np.zeros((cols, rows))
for col in range(cols):
for row in range(rows):
if stable_boundary[col, row] in [1]:
real_val = r_range[row]
# invert cols as min imag value is highest col and vice versa
imag_val = c_range[cols-1 - col]
stable_boundary[col, row] = 1
boundary_points.append((real_val, imag_val))
else:
stable_boundary[col, row] = 0
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, figsize=(5, 5))
ax1.matshow(orig_field, cmap=cmap)
ax2.matshow(stable_field, cmap=cmap)
ax3.matshow(stable_boundary, cmap=cmap)
x = [point[0] for point in boundary_points]
y = [point[1] for point in boundary_points]
ax4.plot(x, y, 'o', c='r', markersize=0.5)
ax4.set_aspect(1)
plt.show()
Output with dpu = 200 and max_iterations = 100:
inspired by this Youtube video: What's so special about the Mandelbrot Set? - Numberphile
Solution
Thanks for the input. As it turned out this is indeed not as easy as it seems. In the end I have used the convex_hull and the alpha shape algorithms to deterimine boundary polygon(s) around the boundary points as shown the picture below. Top left is the juliaset where colors represent the number of iterations; top right black is unstable and white is stable; bottom left is a collection of points representing the boundary between unstable and stable; and bottom right is the collection of boundary polygons around the boundary points.
The code shows below:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
from matplotlib import patches as mpl_patches
from matplotlib.collections import PatchCollection
import shapely.geometry as geometry
from shapely.ops import cascaded_union, polygonize
from scipy.signal import convolve2d
from scipy.spatial import Delaunay # pylint: disable-msg=no-name-in-module
from descartes.patch import PolygonPatch
def juliaset_func(point, constant, max_iterations):
z = point
stable = True
num_iterations = 1
while stable and num_iterations < max_iterations:
z = z**2 + constant
if abs(z) > max(abs(constant), 2):
stable = False
return (stable, num_iterations)
num_iterations += 1
return (stable, num_iterations)
def create_juliaset(r_range, c_range, constant, max_iterations):
''' create a juliaset that returns two fields (matrices) - orig_field and
stable_field, where orig_field contains the number of iterations for
a point in the complex plane (r, c) and stable_field for each point
either whether the point is stable (True) or not stable (False)
'''
points = np.array([])
colors = np.array([])
stables = np.array([], dtype='bool')
progress = 0
for imag in c_range:
for real in r_range:
point = complex(real, imag)
points = np.append(points, point)
stable, color = juliaset_func(point, constant, max_iterations)
stables = np.append(stables, stable)
colors = np.append(colors, color)
print(f'{100*progress/len(c_range)/len(r_range):3.2f}% completed\r', end='')
progress += len(r_range)
print(' \r', end='')
rows = len(r_range)
start = len(colors)
orig_field = []
stable_field = []
for i_num in range(len(c_range)):
start -= rows
real_colors = [color for color in colors[start:start+rows]]
real_stables = [1 if val == True else 0 for val in stables[start:start+rows]]
orig_field.append(real_colors)
stable_field.append(real_stables)
orig_field = np.array(orig_field, dtype='int')
stable_field = np.array(stable_field, dtype='int')
return orig_field, stable_field
def find_boundary_points_of_stable_field(stable_field, r_range, c_range):
''' find the boundary points by convolving the stable_field with a 3x3
kernel of all ones and define the point on the boundary where the
convolution is 1.
'''
kernel = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]], dtype='int8')
stable_boundary = convolve2d(stable_field, kernel, mode='same')
rows = len(r_range)
cols = len(c_range)
boundary_points = []
for col in range(cols):
for row in range(rows):
# Note you can make the boundary 'thicker ' by
# expanding the range of possible values like [1, 2, 3]
if stable_boundary[col, row] in [1]:
real_val = r_range[row]
# invert cols as min imag value is highest col and vice versa
imag_val = c_range[cols-1 - col]
boundary_points.append((real_val, imag_val))
else:
pass
return [geometry.Point(val[0], val[1]) for val in boundary_points]
def alpha_shape(points, alpha):
''' determine the boundary of a cluster of points whereby 'sharpness' of
the boundary depends on alpha.
paramaters:
:points: list of shapely Point objects
:alpha: scalar
returns:
shapely Polygon object or MultiPolygon
edge_points: list of start and end point of each side of the polygons
'''
if len(points) < 4:
# When you have a triangle, there is no sense
# in computing an alpha shape.
return geometry.MultiPoint(list(points)).convex_hull
def add_edge(edges, edge_points, coords, i, j):
"""
Add a line between the i-th and j-th points,
if not in the list already
"""
if (i, j) in edges or (j, i) in edges:
# already added
return
edges.add((i, j))
edge_points.append((coords[[i, j]]))
coords = np.array([point.coords[0]
for point in points])
tri = Delaunay(coords)
edges = set()
edge_points = []
# loop over triangles:
# ia, ib, ic = indices of corner points of the
# triangle
for ia, ib, ic in tri.vertices:
pa = coords[ia]
pb = coords[ib]
pc = coords[ic]
# Lengths of sides of triangle
a = np.sqrt((pa[0]-pb[0])**2 + (pa[1]-pb[1])**2)
b = np.sqrt((pb[0]-pc[0])**2 + (pb[1]-pc[1])**2)
c = np.sqrt((pc[0]-pa[0])**2 + (pc[1]-pa[1])**2)
# Semiperimeter of triangle
s = (a + b + c)/2.0
# Area of triangle by Heron's formula
area = np.sqrt(s*(s-a)*(s-b)*(s-c))
circum_r = a*b*c/(4.0*area)
# Here's the radius filter.
if circum_r < alpha:
add_edge(edges, edge_points, coords, ia, ib)
add_edge(edges, edge_points, coords, ib, ic)
add_edge(edges, edge_points, coords, ic, ia)
m = geometry.MultiLineString(edge_points)
triangles = list(polygonize(m))
return cascaded_union(triangles), edge_points
def main():
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2, figsize=(5, 5))
# define limits, range and resolution in the complex plane
r_min, r_max = -1.5, 1.5
c_min, c_max = -1.1, 1.1
dpu = 100 # dots per unit - 50 dots per 1 units means 200 points per 4 units
intval = 1 / dpu
r_range = np.arange(r_min, r_max + intval, intval)
c_range = np.arange(c_min, c_max + intval, intval)
# create two matrixes (orig_field and stable_field) for the juliaset with
# constant
constant = -0.76 -0.10j
max_iterations = 50
orig_field, stable_field = create_juliaset(r_range, c_range,
constant,
max_iterations)
cmap='nipy_spectral'
ax1.matshow(orig_field, cmap=cmap, interpolation='bilinear')
ax2.matshow(stable_field, cmap=cmap)
# find points that are on the boundary of the stable field
boundary_points = find_boundary_points_of_stable_field(stable_field,
r_range, c_range)
x = [p.x for p in boundary_points]
y = [p.y for p in boundary_points]
ax3.plot(x, y, 'o', c='r', markersize=0.5)
ax3.set_xlim(r_min, r_max)
ax3.set_ylim(c_min, c_max)
ax3.set_aspect(1)
# find the boundary polygon using alpha_shape where 'sharpness' of the
# boundary is determined by the factor ALPHA
# a green boundary consists of multiple polygons, a red boundary on a single
# polygon
alpha = 0.03 # determines shape of the boundary polygon
bnd_polygon, _ = alpha_shape(boundary_points, alpha)
patches = []
if bnd_polygon.geom_type == 'Polygon':
patches.append(PolygonPatch(bnd_polygon))
ec = 'red'
else:
for poly in bnd_polygon:
patches.append(PolygonPatch(poly))
ec = 'green'
p = PatchCollection(patches, facecolor='none', edgecolor=ec, lw=1)
ax4.add_collection(p)
ax4.set_xlim(r_min, r_max)
ax4.set_ylim(c_min, c_max)
ax4.set_aspect(1)
plt.show()
if __name__ == "__main__":
main()