Skip to content

Instantly share code, notes, and snippets.

@PardhavMaradani

PardhavMaradani/1_main.py

Last active September 29, 2021 15:53
Show Gist options
  • Save PardhavMaradani/5e86410e1cbe96fb96262745a9360a0c to your computer and use it in GitHub Desktop.
Save PardhavMaradani/5e86410e1cbe96fb96262745a9360a0c to your computer and use it in GitHub Desktop.
Download ZIP
Blog: Voronoi Whirls
Raw
1_main.py
# https://pardhav-m.blogspot.com/2020/07/voronoi-whirls.html
# https://trinket.io/embed/python/25fa0c9d1a
# Voronoi Whirls
# Basic Vornoi
import turtle
import random
import time
from Voronoi import Voronoi
from voronoi_helpers import get_voronoi_polygons
# voronoi variables to tweak
g_n_seeds = 20
g_show_seeds = True
g_fill_colors = True
# random
seed = random.randrange(65536)
random.seed(seed)
# screen setup
g_screen = turtle.Screen()
g_wh = int(g_screen.window_height())
g_ww = int(1.5 * g_wh)
g_screen.setup(g_ww, g_wh)
g_screen.tracer(0)
# turtle
t = turtle.Turtle()
t.ht()
t.speed(0)
# draw please wait
def draw_please_wait():
t.pu()
t.goto(-g_ww/2 + 20, g_wh/2 - 40)
t.pd()
t.write("Please wait...", font=("Arial", 16, "normal"))
t.update()
time.sleep(0.001)
t.clear()
# draw seeds
def draw_seeds(seeds):
t.pu()
for seed in seeds:
t.goto(seed)
t.dot(3)
# draw polygon
def draw_polygon(polygon, fill_color = False):
if fill_color:
R = random.randrange(0,255)
G = random.randrange(0,255)
B = random.randrange(0,255)
color = (R,G,B)
t.fillcolor(color)
for i, v in enumerate(polygon):
if i == 0:
t.pu()
t.goto(v)
# t.dot()
if i == 0:
t.pd()
t.begin_fill()
t.goto(polygon[0])
if fill_color:
t.end_fill()
# draw voronoi
def draw_voronoi():
# bounding box
maxx = g_ww/2
maxy = g_wh/2
minx = -g_ww/2
miny = -g_wh/2
# create seeds
seeds = []
for i in range(g_n_seeds):
x = random.randint(1, g_ww-1) - g_ww/2
y = random.randint(1, g_wh-1) - g_wh/2
seeds.append((x, y))
# get voronoi edges
voronoi = Voronoi(seeds)
voronoi.process()
voronoi_edges = voronoi.get_output()
# get voronoi polygons
polygons = get_voronoi_polygons(seeds, voronoi_edges, (minx, miny, maxx, maxy))
for polygon in polygons:
draw_polygon(polygon, g_fill_colors)
if g_show_seeds:
draw_seeds(seeds)
# update screen
g_screen.update()
# main
draw_please_wait()
draw_voronoi()
Raw
2_main.py
# https://pardhav-m.blogspot.com/2020/07/voronoi-whirls.html
# https://trinket.io/embed/python/e57355b2b3
# Voronoi Whirls
# Voronoi Whirls
import turtle
import random
import time
from Voronoi import Voronoi
from voronoi_helpers import get_voronoi_polygons
from polygon_helpers import draw_polygon_whirl, draw_polygon_concentric
# voronoi variables to tweak
g_n_seeds = 20
g_gradient = False
g_random_colors = False # when g_gradient = True
g_offset_factor = 0.1 # when g_gradient = False
g_whirl = True # else draw concentric lines
g_equidistant = False # when concentric lines
# random
seed = random.randrange(65536)
# seed = 30737 # samples
random.seed(seed)
# screen setup
g_screen = turtle.Screen()
g_wh = int(g_screen.window_height())
g_ww = int(1.5 * g_wh)
g_screen.setup(g_ww, g_wh)
g_screen.tracer(0)
# turtle
t = turtle.Turtle()
t.ht()
t.speed(0)
t.pensize(0.75)
# draw please wait
def draw_please_wait():
t.pu()
t.goto(-g_ww/2 + 20, g_wh/2 - 40)
t.pd()
t.write("Please wait...", font=("Arial", 16, "normal"))
t.update()
time.sleep(0.001)
t.clear()
# draw voronoi
def draw_voronoi():
# bounding box
maxx = g_ww/2
maxy = g_wh/2
minx = -g_ww/2
miny = -g_wh/2
# create seeds
seeds = []
for i in range(g_n_seeds):
x = random.randint(1, g_ww-1) - g_ww/2
y = random.randint(1, g_wh-1) - g_wh/2
seeds.append((x, y))
# get voronoi edges
voronoi = Voronoi(seeds)
voronoi.process()
voronoi_edges = voronoi.get_output()
# get voronoi polygons
polygons = get_voronoi_polygons(seeds, voronoi_edges, (minx, miny, maxx, maxy))
offset_factor = g_offset_factor
rgb_multiples = (1, 2, 3)
for polygon in polygons:
if g_gradient:
offset_factor = 0.075 if g_whirl else 0.003
if g_random_colors:
m1 = random.randint(0, 2)
m2 = random.randint(0, 2)
m3 = random.randint(0, 2)
rgb_multiples = (m1, m2, m3)
if g_whirl:
draw_polygon_whirl(t, polygon, offset_factor, g_gradient, rgb_multiples)
else:
if g_equidistant:
offset_factor = g_offset_factor / 2
draw_polygon_concentric(t, polygon, offset_factor, g_equidistant, g_gradient, rgb_multiples)
# update screen
g_screen.update()
# main
draw_please_wait()
draw_voronoi()
Raw
3_polygon_helpers.py
import math
# Source: https://algorithmtutor.com/Computational-Geometry/Area-of-a-polygon-given-a-set-of-points/
# Shoelace formula to calculate the area of a polygon
# the points must be sorted anticlockwise (or clockwise)
def polygon_area(vertices):
psum = 0
nsum = 0
for i in range(len(vertices)):
sindex = (i + 1) % len(vertices)
prod = vertices[i][0] * vertices[sindex][1]
psum += prod
for i in range(len(vertices)):
sindex = (i + 1) % len(vertices)
prod = vertices[sindex][0] * vertices[i][1]
nsum += prod
return abs(1/2*(psum - nsum))
# get centroid of polygon
# Source: https://bell0bytes.eu/centroid-convex/
def get_polygon_centroid(polygon):
centroidX = 0
centroidY = 0
det = 0
temDet = 0
j = 0
nVertices = len(polygon)
for i in range(nVertices):
if (i + 1 == nVertices):
j = 0
else:
j = i + 1
x1, y1 = polygon[i]
x2, y2 = polygon[j]
tempDet = x1 * y2 - x2 * y1
det += tempDet
centroidX += (x1 + x2) * tempDet
centroidY += (y1 + y2) * tempDet
centroidX = centroidX / (3 * det)
centroidY = centroidY / (3 * det)
return (centroidX, centroidY)
# draw polygon whirls recursively
def _draw_polygon_whirl(level, t, polygon, offset_factor, gradient, rgb_multiples):
max_levels = 500 if gradient else 250
# limit levels
if level > max_levels:
return
min_area = 50 if gradient else 100
# limit area
if polygon_area(polygon) < min_area:
return
gap = 10 if gradient else 50
# draw current polygon
for i, v in enumerate(polygon):
if i == 0:
t.pu()
if gradient:
r = level * 1.05
m1, m2, m3 = rgb_multiples
t.color(r*m1, r*m2, r*m3)
t.goto(v)
if i == 0:
t.pd()
# create new polygon verties
new_polygon = []
for i, p0 in enumerate(polygon):
if i == len(polygon) - 1:
i = 0
p1 = polygon[i + 1]
x0, y0 = p0
x1, y1 = p1
D = math.sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0))
d = gap * offset_factor
x2 = x0 + (d / D) * (x1 - x0)
y2 = y0 + (d / D) * (y1 - y0)
new_polygon.append((x2, y2))
# recursively call with new polygon
_draw_polygon_whirl(level + 1, t, new_polygon, offset_factor, gradient, rgb_multiples)
# draw polygon concentric recursively
def _draw_polygon_concentric(level, t, center, polygon, offset_factor, equidistant, gradient, rgb_multiples):
max_levels = 1000 if gradient else 250
# limit levels
if level > max_levels:
return
min_area = 50 if gradient else 100
# limit area
if polygon_area(polygon) < min_area:
return
m = 1
if equidistant:
m = level
# draw current polygon
for i, v in enumerate(polygon):
if i == 0:
t.pu()
if gradient:
r = 10 + (level)*0.25
m1, m2, m3 = rgb_multiples
t.color(r*m1, r*m2, r*m3)
t.goto(v)
if i == 0:
t.pd()
# create new polygon verties
new_polygon = []
for i, p0 in enumerate(polygon):
if i == len(polygon) - 1:
i = 0
p1 = center
x0, y0 = p0
x1, y1 = p1
x2 = x0 + offset_factor * m * (x1 - x0)
y2 = y0 + offset_factor * m * (y1 - y0)
new_polygon.append((x2, y2))
# recursively call with new polygon
_draw_polygon_concentric(level + 1, t, center, new_polygon, offset_factor, equidistant, gradient, rgb_multiples)
def draw_polygon_whirl(t, polygon, offset_factor, gradient = False, rgb_multiples = (1, 2, 3)):
_draw_polygon_whirl(0, t, polygon, offset_factor, gradient, rgb_multiples)
def draw_polygon_concentric(t, polygon, offset_factor, equidistant = False, gradient = False, rgb_multiples = (1, 2, 3)):
center = get_polygon_centroid(polygon)
if gradient:
equidistant = True
offset_factor = 0.00001
_draw_polygon_concentric(0, t, center, polygon, offset_factor, equidistant, gradient, rgb_multiples)
Raw
3_voronoi_helpers.py
# https://pardhav-m.blogspot.com/2020/07/voronoi-whirls.html
import math
from line_clip import cohensutherland
from DataType import Point
# check if point c is between points a and b
# Source: https://stackoverflow.com/questions/328107/how-can-you-determine-a-point-is-between-two-other-points-on-a-line-segment
def isBetween(a, b, c):
eps = 0.01
crossproduct = (c.y - a.y) * (b.x - a.x) - (c.x - a.x) * (b.y - a.y)
# compare versus epsilon for floating point values, or != 0 if using integers
if abs(crossproduct) > eps:
return False
dotproduct = (c.x - a.x) * (b.x - a.x) + (c.y - a.y)*(b.y - a.y)
if dotproduct < 0:
return False
squaredlengthba = (b.x - a.x)*(b.x - a.x) + (b.y - a.y)*(b.y - a.y)
if dotproduct > squaredlengthba:
return False
return True
# simplify polygon by removing multiple vertices on same line
def simplify_polygon(polygon):
simplified_polygon = []
end = len(polygon) - 1
for i, v in enumerate(polygon):
pi = i - 1 # previous index
ni = i + 1 # next index
if pi < 0:
pi = end
if ni > end:
ni = 0
pv = polygon[pi] # previous vertex
nv = polygon[ni] # next vertex
a = Point(pv[0], pv[1])
b = Point(nv[0], nv[1])
c = Point(v[0], v[1])
if not isBetween(a, b, c):
simplified_polygon.append(v)
simplified_polygon.append(simplified_polygon[0])
return simplified_polygon
def get_voronoi_polygons(seeds, edges, bounding_box):
(minx, miny, maxx, maxy) = bounding_box
polygons = []
# 1. clip edges
# save border vertices after clipping
border_vertices = [None] * 4
for i in range(4):
border_vertices[i] = []
# clip edges
clipped_edges = []
for edge in edges:
(x1, y1, x2, y2) = edge
(nx1, ny1, nx2, ny2) = cohensutherland(minx, maxy, maxx, miny, x1, y1, x2, y2)
if nx1 == None:
continue
clipped_edges.append((nx1, ny1, nx2, ny2))
# save edges on border for next step
if ny1 == miny or ny1 == maxy:
border_vertices[0 if ny1 == miny else 2].append(nx1)
if ny2 == miny or ny2 == maxy:
border_vertices[0 if ny2 == miny else 2].append(nx2)
if nx1 == minx or nx1 == maxx:
border_vertices[3 if nx1 == minx else 1].append(ny1)
if nx2 == minx or nx2 == maxx:
border_vertices[3 if nx2 == minx else 1].append(ny2)
# 2. connect borders
# bb 0, 2
for i in range(0, 4, 2):
y = miny if i == 0 else maxy
x1 = minx
for v in sorted(border_vertices[i]):
clipped_edges.append((x1, y, v, y))
x1 = v
clipped_edges.append((x1, y, maxx, y))
# bb 1, 3
for i in range(1, 4, 2):
x = minx if i == 3 else maxx
y1 = miny
for v in sorted(border_vertices[i]):
clipped_edges.append((x, y1, x, v))
y1 = v
clipped_edges.append((x, y1, x, maxy))
# 3. get polygons
# edges array for each seed
seed_edges = [None] * len(seeds)
for i, seed in enumerate(seeds):
seed_edges[i] = []
# Find the seed (cell) that the edge belongs to
for edge in clipped_edges:
(x1, y1, x2, y2) = edge
midx = (x1 + x2) / 2
midy = (y1 + y2) / 2
distances = []
i = 0
# calculate distance from edge to each seed
for seed in seeds:
sx, sy = seed
distance = math.sqrt((sx-midx)*(sx-midx) + (sy-midy)*(sy-midy))
distances.append((distance, i))
i += 1
# sort the distances to get closest seeds (cells)
sorted_distances = sorted(distances)
_, si1 = sorted_distances[0]
_, si2 = sorted_distances[1]
# for border edges, we only need one cell, for all others, two
only_one = False
if x1 == x2 and (x1 == minx or x1 == maxx):
only_one = True
if y1 == y2 and (y1 == miny or y1 == maxy):
only_one = True
seed_edges[si1].append(edge)
if not only_one:
seed_edges[si2].append(edge)
# 4. create polygons from edges
polygons = []
for i in range(len(seeds)):
# get vertices
vertices = []
for edge in seed_edges[i]:
(x1, y1, x2, y2) = edge
v1 = (x1, y1)
v2 = (x2, y2)
if v1 not in vertices:
vertices.append(v1)
if v2 not in vertices:
vertices.append(v2)
# find approx center of polygon from vertices
midx = 0
midy = 0
for v in vertices:
x, y = v
midx += x
midy += y
midx /= len(vertices)
midy /= len(vertices)
# calculate angle to approx center to sort vertices
sorted_vertices = []
for v in vertices:
x, y = v
angle = math.atan2((y - midy), (x - midx))
sorted_vertices.append((angle, v))
# create polygon from sorted vertices
polygon = []
for e in sorted(sorted_vertices):
_, v = e
polygon.append(v)
# 5. simplify polygons (remove multiple vertices on same line)
polygons.append(simplify_polygon(polygon))
# return polygons
return polygons
Raw
4_DataType.py
# Source: https://github.com/jansonh/Voronoi
from min_heapq import heappush, heappop
class Point:
x = 0.0
y = 0.0
def __init__(self, x, y):
self.x = x
self.y = y
class Event:
x = 0.0
p = None
a = None
valid = True
def __init__(self, x, p, a):
self.x = x
self.p = p
self.a = a
self.valid = True
class Arc:
p = None
pprev = None
pnext = None
e = None
s0 = None
s1 = None
def __init__(self, p, a=None, b=None):
self.p = p
self.pprev = a
self.pnext = b
self.e = None
self.s0 = None
self.s1 = None
class Segment:
start = None
end = None
done = False
def __init__(self, p):
self.start = p
self.end = None
self.done = False
def finish(self, p):
if self.done: return
self.end = p
self.done = True
class PriorityQueue:
def __init__(self):
self.pq = []
self.entry_finder = {}
self.counter = 0 #itertools.count()
def push(self, item):
# check for duplicate
if item in self.entry_finder: return
count = self.counter
self.counter += 1
# use x-coordinate as a primary key (heapq in python is min-heap)
entry = [item.x, count, item]
self.entry_finder[item] = entry
heappush(self.pq, entry)
def remove_entry(self, item):
entry = self.entry_finder.pop(item)
entry[-1] = 'Removed'
def pop(self):
while self.pq:
priority, count, item = heappop(self.pq)
if item is not 'Removed':
del self.entry_finder[item]
return item
raise KeyError('pop from an empty priority queue')
def top(self):
while self.pq:
priority, count, item = heappop(self.pq)
if item is not 'Removed':
del self.entry_finder[item]
self.push(item)
return item
raise KeyError('top from an empty priority queue')
def empty(self):
return not self.pq
Raw
4_line_clip.py
# Source: https://github.com/scivision/lineclipping-python-fortran/blob/master/pylineclip/__init__.py
'''
The MIT License (MIT)
Copyright (c) 2014 Michael Hirsch
reference: http://en.wikipedia.org/wiki/Cohen%E2%80%93Sutherland_algorithm
* The best way to Numba JIT this would probably be in the function calling this,
to include the loop itself inside the jit decoration.
'''
def cohensutherland(xmin, ymax, xmax, ymin, x1, y1, x2, y2):
"""Clips a line to a rectangular area.
This implements the Cohen-Sutherland line clipping algorithm. xmin,
ymax, xmax and ymin denote the clipping area, into which the line
defined by x1, y1 (start point) and x2, y2 (end point) will be
clipped.
If the line does not intersect with the rectangular clipping area,
four None values will be returned as tuple. Otherwise a tuple of the
clipped line points will be returned in the form (cx1, cy1, cx2, cy2).
"""
INSIDE, LEFT, RIGHT, LOWER, UPPER = 0, 1, 2, 4, 8
def _getclip(xa, ya):
# if dbglvl>1: print('point: '),; print(xa,ya)
p = INSIDE # default is inside
# consider x
if xa < xmin:
p |= LEFT
elif xa > xmax:
p |= RIGHT
# consider y
if ya < ymin:
p |= LOWER # bitwise OR
elif ya > ymax:
p |= UPPER # bitwise OR
return p
# check for trivially outside lines
k1 = _getclip(x1, y1)
k2 = _getclip(x2, y2)
# %% examine non-trivially outside points
# bitwise OR |
while (k1 | k2) != 0: # if both points are inside box (0000) , ACCEPT trivial whole line in box
# if line trivially outside window, REJECT
if (k1 & k2) != 0: # bitwise AND &
# if dbglvl>1: print(' REJECT trivially outside box')
# return nan, nan, nan, nan
return None, None, None, None
# non-trivial case, at least one point outside window
# this is not a bitwise or, it's the word "or"
opt = k1 or k2 # take first non-zero point, short circuit logic
if opt & UPPER: # these are bitwise ANDS
x = x1 + (x2 - x1) * (ymax - y1) / (y2 - y1)
y = ymax
elif opt & LOWER:
x = x1 + (x2 - x1) * (ymin - y1) / (y2 - y1)
y = ymin
elif opt & RIGHT:
y = y1 + (y2 - y1) * (xmax - x1) / (x2 - x1)
x = xmax
elif opt & LEFT:
y = y1 + (y2 - y1) * (xmin - x1) / (x2 - x1)
x = xmin
else:
raise RuntimeError('Undefined clipping state')
if opt == k1:
x1, y1 = x, y
k1 = _getclip(x1, y1)
elif opt == k2:
x2, y2 = x, y
k2 = _getclip(x2, y2)
return x1, y1, x2, y2
Raw
4_min_heapq.py
# Source: https://github.com/python/cpython/tree/3.8/Lib/heapq.py
def heappush(heap, item):
"""Push item onto heap, maintaining the heap invariant."""
heap.append(item)
_siftdown(heap, 0, len(heap)-1)
def heappop(heap):
"""Pop the smallest item off the heap, maintaining the heap invariant."""
lastelt = heap.pop() # raises appropriate IndexError if heap is empty
if heap:
returnitem = heap[0]
heap[0] = lastelt
_siftup(heap, 0)
return returnitem
return lastelt
def _siftdown(heap, startpos, pos):
newitem = heap[pos]
# Follow the path to the root, moving parents down until finding a place
# newitem fits.
while pos > startpos:
parentpos = (pos - 1) >> 1
parent = heap[parentpos]
if newitem < parent:
heap[pos] = parent
pos = parentpos
continue
break
heap[pos] = newitem
def _siftup(heap, pos):
endpos = len(heap)
startpos = pos
newitem = heap[pos]
# Bubble up the smaller child until hitting a leaf.
childpos = 2*pos + 1 # leftmost child position
while childpos < endpos:
# Set childpos to index of smaller child.
rightpos = childpos + 1
if rightpos < endpos and not heap[childpos] < heap[rightpos]:
childpos = rightpos
# Move the smaller child up.
heap[pos] = heap[childpos]
pos = childpos
childpos = 2*pos + 1
# The leaf at pos is empty now. Put newitem there, and bubble it up
# to its final resting place (by sifting its parents down).
heap[pos] = newitem
_siftdown(heap, startpos, pos)
Raw
4_Voronoi.py
# Source: https://github.com/jansonh/Voronoi
import math
from DataType import Point, Event, Arc, Segment, PriorityQueue
# Source: (C++) http://www.cs.hmc.edu/~mbrubeck/voronoi.html
class Voronoi:
def __init__(self, points):
self.output = [] # list of line segment
self.arc = None # binary tree for parabola arcs
self.points = PriorityQueue() # site events
self.event = PriorityQueue() # circle events
# bounding box
self.x0 = -0.0
self.x1 = -0.0
self.y0 = 0.0
self.y1 = 0.0
# insert points to site event
for pts in points:
point = Point(pts[0], pts[1])
self.points.push(point)
# keep track of bounding box size
if point.x < self.x0: self.x0 = point.x
if point.y < self.y0: self.y0 = point.y
if point.x > self.x1: self.x1 = point.x
if point.y > self.y1: self.y1 = point.y
# add margins to the bounding box
dx = (self.x1 - self.x0 + 1) / 5.0
dy = (self.y1 - self.y0 + 1) / 5.0
self.x0 = self.x0 - dx
self.x1 = self.x1 + dx
self.y0 = self.y0 - dy
self.y1 = self.y1 + dy
def process(self):
while not self.points.empty():
if not self.event.empty() and (self.event.top().x <= self.points.top().x):
self.process_event() # handle circle event
else:
self.process_point() # handle site event
# after all points, process remaining circle events
while not self.event.empty():
self.process_event()
self.finish_edges()
def process_point(self):
# get next event from site pq
p = self.points.pop()
# add new arc (parabola)
self.arc_insert(p)
def process_event(self):
# get next event from circle pq
e = self.event.pop()
if e.valid:
# start new edge
s = Segment(e.p)
self.output.append(s)
# remove associated arc (parabola)
a = e.a
if a.pprev is not None:
a.pprev.pnext = a.pnext
a.pprev.s1 = s
if a.pnext is not None:
a.pnext.pprev = a.pprev
a.pnext.s0 = s
# finish the edges before and after a
if a.s0 is not None: a.s0.finish(e.p)
if a.s1 is not None: a.s1.finish(e.p)
# recheck circle events on either side of p
if a.pprev is not None: self.check_circle_event(a.pprev, e.x)
if a.pnext is not None: self.check_circle_event(a.pnext, e.x)
def arc_insert(self, p):
if self.arc is None:
self.arc = Arc(p)
else:
# find the current arcs at p.y
i = self.arc
while i is not None:
flag, z = self.intersect(p, i)
if flag:
# new parabola intersects arc i
flag, zz = self.intersect(p, i.pnext)
if (i.pnext is not None) and (not flag):
i.pnext.pprev = Arc(i.p, i, i.pnext)
i.pnext = i.pnext.pprev
else:
i.pnext = Arc(i.p, i)
i.pnext.s1 = i.s1
# add p between i and i.pnext
i.pnext.pprev = Arc(p, i, i.pnext)
i.pnext = i.pnext.pprev
i = i.pnext # now i points to the new arc
# add new half-edges connected to i's endpoints
seg = Segment(z)
self.output.append(seg)
i.pprev.s1 = i.s0 = seg
seg = Segment(z)
self.output.append(seg)
i.pnext.s0 = i.s1 = seg
# check for new circle events around the new arc
self.check_circle_event(i, p.x)
self.check_circle_event(i.pprev, p.x)
self.check_circle_event(i.pnext, p.x)
return
i = i.pnext
# if p never intersects an arc, append it to the list
i = self.arc
while i.pnext is not None:
i = i.pnext
i.pnext = Arc(p, i)
# insert new segment between p and i
x = self.x0
y = (i.pnext.p.y + i.p.y) / 2.0;
start = Point(x, y)
seg = Segment(start)
i.s1 = i.pnext.s0 = seg
self.output.append(seg)
def check_circle_event(self, i, x0):
# look for a new circle event for arc i
if (i.e is not None) and (i.e.x <> self.x0):
i.e.valid = False
i.e = None
if (i.pprev is None) or (i.pnext is None): return
flag, x, o = self.circle(i.pprev.p, i.p, i.pnext.p)
if flag and (x > self.x0):
i.e = Event(x, o, i)
self.event.push(i.e)
def circle(self, a, b, c):
# check if bc is a "right turn" from ab
if ((b.x - a.x)*(c.y - a.y) - (c.x - a.x)*(b.y - a.y)) > 0: return False, None, None
# Joseph O'Rourke, Computational Geometry in C (2nd ed.) p.189
A = b.x - a.x
B = b.y - a.y
C = c.x - a.x
D = c.y - a.y
E = A*(a.x + b.x) + B*(a.y + b.y)
F = C*(a.x + c.x) + D*(a.y + c.y)
G = 2*(A*(c.y - b.y) - B*(c.x - b.x))
if (G == 0): return False, None, None # Points are co-linear
# point o is the center of the circle
ox = 1.0 * (D*E - B*F) / G
oy = 1.0 * (A*F - C*E) / G
# o.x plus radius equals max x coord
x = ox + math.sqrt((a.x-ox)**2 + (a.y-oy)**2)
o = Point(ox, oy)
return True, x, o
def intersect(self, p, i):
# check whether a new parabola at point p intersect with arc i
if (i is None): return False, None
if (i.p.x == p.x): return False, None
a = 0.0
b = 0.0
if i.pprev is not None:
a = (self.intersection(i.pprev.p, i.p, 1.0*p.x)).y
if i.pnext is not None:
b = (self.intersection(i.p, i.pnext.p, 1.0*p.x)).y
if (((i.pprev is None) or (a <= p.y)) and ((i.pnext is None) or (p.y <= b))):
py = p.y
px = 1.0 * ((i.p.x)**2 + (i.p.y-py)**2 - p.x**2) / (2*i.p.x - 2*p.x)
res = Point(px, py)
return True, res
return False, None
def intersection(self, p0, p1, l):
# get the intersection of two parabolas
p = p0
if (p0.x == p1.x):
py = (p0.y + p1.y) / 2.0
elif (p1.x == l):
py = p1.y
elif (p0.x == l):
py = p0.y
p = p1
else:
# use quadratic formula
z0 = 2.0 * (p0.x - l)
z1 = 2.0 * (p1.x - l)
a = 1.0/z0 - 1.0/z1;
b = -2.0 * (p0.y/z0 - p1.y/z1)
c = 1.0 * (p0.y**2 + p0.x**2 - l**2) / z0 - 1.0 * (p1.y**2 + p1.x**2 - l**2) / z1
py = 1.0 * (-b-math.sqrt(b*b - 4*a*c)) / (2*a)
px = 1.0 * (p.x**2 + (p.y-py)**2 - l**2) / (2*p.x-2*l)
res = Point(px, py)
return res
def finish_edges(self):
l = self.x1 + (self.x1 - self.x0) + (self.y1 - self.y0)
i = self.arc
while i.pnext is not None:
if i.s1 is not None:
p = self.intersection(i.p, i.pnext.p, l*2.0)
i.s1.finish(p)
i = i.pnext
def get_output(self):
res = []
for o in self.output:
p0 = o.start
p1 = o.end
if o.done:
res.append((p0.x, p0.y, p1.x, p1.y))
return res
@wenjun-dbt
Copy link

Thank you for sharing this code. After I run it on my computer it shows: TurtleGraphicsError: bad color sequence: (1.05, 2.1, 3.1500000000000004). Do you have idea to solve this problem? Thank you in advance.

@PardhavMaradani
Copy link
Author

The code works fine on Trinket as you can see from the blog page (https://pardhav-m.blogspot.com/2020/07/voronoi-whirls.html). Doing a quick search, I came across the following if you are not using Trinket: raspberrypilearning/turtley-amazing#3 - maybe you can try setting the colormode as suggested there?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment