Delaunay triangulation in Python with opencv and igraph, strange behavior
So I have some code in Python (3.9.13) to obtain a Delaunay triangulation of a set of points in real time and analyze the graph properties. First I use OpenCV (opencv-python 4.6.0.66) Subdiv2D method to obtain the triangulation. Then I convert it in a graph I can analyze with igraph (igraph 0.10.3). But I am not sure why, once every few frames the graph produced by igraph is messed up such as shown in this image:
graph messed up (Left is OpenCV and right is igraph):
Else it is working properly.
good graph (Left is OpenCV and right is igraph):
Here is my demo code:
import time
import numpy
import cv2
import igraph as ig
# Draw a point
def draw_point(img, p, color):
cv2.circle(img, (int(p[0]), int(p[1])), 2, color, 0)
# Get a triangulation
def get_delaunay(subdiv):
return subdiv.getTriangleList()
# Draw delaunay triangles
def draw_delaunay(img, subdiv, delaunay_color):
triangleList = subdiv.getTriangleList()
size = img.shape
r = (0, 0, size[1], size[0])
for t in triangleList:
pt1 = (int(t[0]), int(t[1]))
pt2 = (int(t[2]), int(t[3]))
pt3 = (int(t[4]), int(t[5]))
cv2.line(img, pt1, pt2, delaunay_color, 1, cv2.LINE_AA, 0)
cv2.line(img, pt2, pt3, delaunay_color, 1, cv2.LINE_AA, 0)
cv2.line(img, pt3, pt1, delaunay_color, 1, cv2.LINE_AA, 0)
if __name__ == '__main__':
NUM_PART = 500
SIZE = 1000
REPEAT = 10
for iteration in range(REPEAT):
positions = numpy.random.randint(0, SIZE, size=(NUM_PART, 2))
print("There is {p} positions. And {up} unique position".format(p=len(positions), up=len(numpy.unique(positions, axis=1))))
# Create an instance of Subdiv2D
rect = (0, 0, SIZE, SIZE)
subdiv = cv2.Subdiv2D(rect)
timer = time.time()
# Insert points into subdiv
print("There is {} points in subdiv".format(len(positions)))
for p in positions:
p = p.astype("float32")
subdiv.insert(p)
# get triangulation
trilist = get_delaunay(subdiv)
print("Took {}s".format(round(time.time() - timer, 12)))
print("there is {} triangles in trilist".format(len(trilist)))
# create image
opencv_image = numpy.zeros((SIZE, SIZE, 3))
# Draw delaunay triangles
draw_delaunay(opencv_image, subdiv, (255, 255, 255))
# Draw points
for p in positions:
draw_point(opencv_image, p, (0, 0, 255))
timer = time.time()
n_vertices = NUM_PART
# create graph
g = ig.Graph(n=n_vertices, )
g.vs["name"] = range(NUM_PART)
print("graph name vector of length {l}:\n{v}".format(l=len(g.vs["name"]), v=g.vs["name"]))
# Inversion x positions
positionsx = [SIZE - pos for pos in positions[:, 0]]
g.vs["x"] = positions[:, 0]
print("graph x vector of length {l}:\n{v}".format(l=len(g.vs["x"]), v=g.vs["x"]))
g.vs["y"] = positions[:, 1]
print("graph y vector of length {l}:\n{v}".format(l=len(g.vs["y"]), v=g.vs["y"]))
print("Graph took {}s".format(round(time.time() - timer, 12)))
list_vtx = []
for tri in trilist:
vertex1, _ = subdiv.findNearest((tri[0], tri[1]))
vertex2, _ = subdiv.findNearest((tri[2], tri[3]))
vertex3, _ = subdiv.findNearest((tri[4], tri[5]))
list_vtx.extend([vertex3, vertex2, vertex1])
list_cleared = list(dict.fromkeys(list_vtx))
list_cleared.sort()
print("list cleared of length {len}: {lst}".format(len=len(list_cleared), lst=list_cleared))
for tri in trilist:
vertex1, _ = subdiv.findNearest((tri[0], tri[1]))
vertex2, _ = subdiv.findNearest((tri[2], tri[3]))
vertex3, _ = subdiv.findNearest((tri[4], tri[5]))
#print("vertex 1: {v} of position {p}".format(v=vertex1, p=(tri[0], tri[1])))
#print("vertex 2: {v} of position {p}".format(v=vertex2, p=(tri[2], tri[3])))
#print("vertex 3: {v} of position {p}".format(v=vertex3, p=(tri[4], tri[5])))
# -4 because https://stackoverflow.com/a/52377891/18493005
g.add_edges([
(vertex1 - 4, vertex2 - 4),
(vertex2 - 4, vertex3 - 4),
(vertex3 - 4, vertex1 - 4),
])
# simplify graph
g.simplify()
nodes = g.vs.indices
print(nodes)
print(subdiv)
# create image
igraph_image = numpy.zeros((SIZE, SIZE, 3))
for point in g.vs:
draw_point(igraph_image, (point["x"], point["y"]), (0, 0, 255))
for edge in g.es:
# print(edge.tuple)
# print(g.vs["x"][edge.tuple[0]])
cv2.line(igraph_image, (int(g.vs["x"][edge.tuple[0]]), int(g.vs["y"][edge.tuple[0]])),
(int(g.vs["x"][edge.tuple[1]]), int(g.vs["y"][edge.tuple[1]])), (255, 255, 255), 1, cv2.LINE_AA, 0)
numpy_horizontal = numpy.hstack((opencv_image, igraph_image))
# Show results
cv2.imshow('L: opencv || R: igraph', numpy_horizontal)
cv2.waitKey(0)
I try to have a repeatable result of my graph in igraph. But it is only working 80% of the time which is pretty strange behavior. Any idea of what are my mistakes here?
Edit: it seems to be a variation in the length of the list generated by:
trilist = get_delaunay(subdiv)
list_vtx = []
for tri in trilist:
vertex1, _ = subdiv.findNearest((tri[0], tri[1]))
vertex2, _ = subdiv.findNearest((tri[2], tri[3]))
vertex3, _ = subdiv.findNearest((tri[4], tri[5]))
list_vtx.extend([vertex3, vertex2, vertex1])
list_cleared = list(dict.fromkeys(list_vtx))
list_cleared.sort()
but I am not sure why.
Edit2:
After the modification sugested by Markus. I do not get a messed up graph anymore. But now the graph is missing some edges
x_pos = [0] * NUM_PART # create 0-filled array of x-positions
y_pos = [0] * NUM_PART # create 0-filled array of y-positions
edges = [] # create empty array of edges
# for each triangle add vertex positions and edges
for tri in trilist:
vertex1 = subdiv.findNearest((tri[0], tri[1]))[0] - 4
vertex2 = subdiv.findNearest((tri[2], tri[3]))[0] - 4
vertex3 = subdiv.findNearest((tri[4], tri[5]))[0] - 4
x_pos[vertex1] = tri[0]
y_pos[vertex1] = tri[1]
x_pos[vertex2] = tri[2]
y_pos[vertex2] = tri[3]
x_pos[vertex3] = tri[4]
y_pos[vertex3] = tri[5]
edges.append((vertex1, vertex2))
edges.append((vertex2, vertex3))
edges.append((vertex2, vertex3))
# create graph
g = ig.Graph(NUM_PART, edges)
g.vs["name"] = range(NUM_PART)
g.vs["x"] = x_pos
g.vs["y"] = y_pos
g.simplify()
The following image shows an overlay between 3 type of drawing (White=opencv , Red=Markus suggestion, Green + Red = previous method used)
Overlay of Markus solution in case of no mess up
Overlay of Markus solution in case of mess up
So Markus solution indeed remove the mess up, but also some edges, even in the case that was working previously.
So in fact my test code was working as expected. The issue was not from Subdiv2D or igraph but from the generation of my position.
I made a mistake verifying the uniqueness of my position with
len(numpy.unique(positions, axis=1))
but should have been using
len(numpy.unique(positions, axis=0)).
So when I used subdiv.findNearest()[0] or subdiv.locate()[2] I was in fact finding several points at the same position, and only the first index was thrown back by the function and so the graph was being messed up.
In order to generate unique position I uses the following code and the graph messing disappeared:
rng = numpy.random.default_rng()
xx = rng.choice(range(SIZE), NUM_PART, replace=False).astype("float32")
yy = rng.choice(range(SIZE), NUM_PART, replace=False).astype("float32")
pp= numpy.stack((xx,yy), axis=1)
The fact that numpy.random.randint(0, 1000, size=(500, 2)) was providing a similar position every 10 or so frame is pretty strange to me as the probability of getting two identical positions seems intuitively to be lower than 0.1