delaunay2D_Delaunay_pythondelaunay_三角剖分_

# -*- coding: ascii -*- """ Simple structured Delaunay triangulation in 2D with Bowyer-Watson algorithm. Written by Jose M. Espadero ( http://github.com/jmespadero/pyDelaunay2D ) Based on code from Ayron Catteau. Published at http://github.com/ayron/delaunay Just pretend to be simple and didactic. The only requisite is numpy. Robust checks disabled by default. May not work in degenerate set of points. """ import numpy as np from math import sqrt class Delaunay2D: """ Class to compute a Delaunay triangulation in 2D ref: http://en.wikipedia.org/wiki/Bowyer-Watson_algorithm ref: http://www.geom.uiuc.edu/~samuelp/del_project.html """ def __init__(self, center=(0, 0), radius=9999): """ Init and create a new frame to contain the triangulation center -- Optional position for the center of the frame. Default (0,0) radius -- Optional distance from corners to the center. """ center = np.asarray(center) # Create coordinates for the corners of the frame self.coords = [center+radius*np.array((-1, -1)), center+radius*np.array((+1, -1)), center+radius*np.array((+1, +1)), center+radius*np.array((-1, +1))] # Create two dicts to store triangle neighbours and circumcircles. self.triangles = {} self.circles = {} # Create two CCW triangles for the frame T1 = (0, 1, 3) T2 = (2, 3, 1) self.triangles[T1] = [T2, None, None] self.triangles[T2] = [T1, None, None] # Compute circumcenters and circumradius for each triangle for t in self.triangles: self.circles[t] = self.circumcenter(t) def circumcenter(self, tri): """Compute circumcenter and circumradius of a triangle in 2D. Uses an extension of the method described here: http://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html """ pts = np.asarray([self.coords[v] for v in tri]) pts2 = np.dot(pts, pts.T) A = np.bmat([[2 * pts2, [[1], [1], [1]]], [[[1, 1, 1, 0]]]]) b = np.hstack((np.sum(pts * pts, axis=1), [1])) x = np.linalg.solve(A, b) bary_coords = x[:-1] center = np.dot(bary_coords, pts) # radius = np.linalg.norm(pts[0] - center) # euclidean distance radius = np.sum(np.square(pts[0] - center)) # squared distance return (center, radius) def inCircleFast(self, tri, p): """Check if point p is inside of precomputed circumcircle of tri. """ center, radius = self.circles[tri] return np.sum(np.square(center - p)) <= radius def inCircleRobust(self, tri, p): """Check if point p is inside of circumcircle around the triangle tri. This is a robust predicate, slower than compare distance to centers ref: http://www.cs.cmu.edu/~quake/robust.html """ m1 = np.asarray([self.coords[v] - p for v in tri]) m2 = np.sum(np.square(m1), axis=1).reshape((3, 1)) m = np.hstack((m1, m2)) # The 3x3 matrix to check return np.linalg.det(m) <= 0 def addPoint(self, p): """Add a point to the current DT, and refine it using Bowyer-Watson. """ p = np.asarray(p) idx = len(self.coords) # print("coords[", idx,"] ->",p) self.coords.append(p) # Search the triangle(s) whose circumcircle contains p bad_triangles = [] for T in self.triangles: # Choose one method: inCircleRobust(T, p) or inCircleFast(T, p) if self.inCircleFast(T, p): bad_triangles.append(T) # Find the CCW boundary (star shape) of the bad triangles, # expressed as a list of edges (point pairs) and the opposite # triangle to each edge. boundary = [] # Choose a "random" triangle and edge T = bad_triangles[0] edge = 0 # get the opposite triangle of this edge while True: # Check if edge of triangle T is on the boundary... # if opposite triangle of this edge is external to the list tri_op = self.triangles[T][edge] if tri_op not in bad_triangles: # Insert edge and external triangle into boundary list boundary.append((T[(edge+1) % 3], T[(edge-1) % 3], tri_op)) # Move to next CCW edge in this triangle edge = (edge + 1) % 3 # Check if boundary is a closed loop if boundary[0][0] == boundary[-1][1]: break else: # Move to next CCW edge in opposite triangle edge = (self.triangles[tri_op].index(T) + 1) % 3 T = tri_op # Remove triangles too near of point p of our solution for T in bad_triangles: del self.triangles[T] del self.circles[T] # Retriangle the hole left by bad_triangles new_triangles = [] for (e0, e1, tri_op) in boundary: # Create a new triangle using point p and edge extremes T = (idx, e0, e1) # Store circumcenter and circumradius of the triangle self.circles[T] = self.circumcenter(T) # Set opposite triangle of the edge as neighbour of T self.triangles[T] = [tri_op, None, None] # Try to set T as neighbour of the opposite triangle if tri_op: # search the neighbour of tri_op that use edge (e1, e0) for i, neigh in enumerate(self.triangles[tri_op]): if neigh: if e1 in neigh and e0 in neigh: # change link to use our new triangle self.triangles[tri_op][i] = T # Add triangle to a temporal list new_triangles.append(T) # Link the new triangles each another N = len(new_triangles) for i, T in enumerate(new_triangles): self.triangles[T][1] = new_triangles[(i+1) % N] # next self.triangles[T][2] = new_triangles[(i-1) % N] # previous def exportTriangles(self): """Export the current list of Delaunay triangles """ # Filter out triangles with any vertex in the extended BBox return [(a-4, b-4, c-4) for (a, b, c) in self.triangles if a > 3 and b > 3 and c > 3] def exportCircles(self): """Export the circumcircles as a list of (center, radius) """ # Remember to compute circumcircles if not done before # for t in self.triangles: # self.circles[t] = self.circumcenter(t) # Filter out triangles with any vertex in the extended BBox # Do sqrt of radius before of return return [(self.circles[(a, b, c)][0], sqrt(self.circles[(a, b, c)][1])) for (a, b, c) in self.triangles if a > 3 and b > 3 and c > 3] def exportDT(self): """Export the current set of Delaunay coordinates and triangles. """ # Filter out coordinates in the extended BBox coord = self.coords[4:] # Filter out triangles with any vertex in the extended BBox tris = [(a-4, b-4, c-4) for (a, b, c) in self.triangles if a > 3 and b > 3 and c > 3] return coord, tris def exportExtendedDT(self): """Export the Extended Delaunay Triangulation (with the frame vertex). """ return self.coords, list(self.triangles) def exportVoronoiRegions(self): """Export coordinates and regions of Voronoi diagram as indexed