From 324ae0c7f1a392133f23d11cca86c7568ab02641 Mon Sep 17 00:00:00 2001 From: jmespadero Date: Tue, 13 Dec 2016 20:37:25 +0100 Subject: [PATCH] Initial commit --- .gitignore | 48 ++++++++++++ delaunay2D.py | 190 +++++++++++++++++++++++++++++++++++++++++++++ delaunay2D_test.py | 55 +++++++++++++ 3 files changed, 293 insertions(+) create mode 100644 .gitignore create mode 100644 delaunay2D.py create mode 100755 delaunay2D_test.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..8892c8d --- /dev/null +++ b/.gitignore @@ -0,0 +1,48 @@ +# Compiled source # +################### +*.com +*.class +*.dll +*.exe +*.o +*.so + +lib + +# Packages # +############ +# it's better to unpack these files and commit the raw source +# git has its own built in compression methods +*.7z +*.dmg +*.gz +*.iso +*.jar +*.rar +*.tar +*.zip + +# Logs and databases # +###################### +*.log +*.sql +*.sqlite + +# OS generated files # +###################### +.DS_Store +.DS_Store? +._* +.Spotlight-V100 +.Trashes +ehthumbs.db +Thumbs.db + +# Backups files +*.pyo +*.pyc +*~ +*.blend1 +tmp* + +output*blend diff --git a/delaunay2D.py b/delaunay2D.py new file mode 100644 index 0000000..abd71f4 --- /dev/null +++ b/delaunay2D.py @@ -0,0 +1,190 @@ +""" +Simple structured Delaunay triangulation in 2D + +Written by Jose M. Espadero ( http://github.com/jmespadero/pyDelaunay2D ) +Based on code from Ayron Catteau. Published at http://github.com/ayron/delaunay + +Pretend to be simple and didactic. The only requisite is numpy. +Lack of robustness checks, so do not spect to work on degenerate set of points. +""" + +import numpy as np + +class Triangle: + """Simple class to store an indexed triangle with neighbour information.""" + + def __init__(self, a, b, c): + """Constructor from 3 indexes to vertex""" + self.v = [a, b, c] + self.neighbour = [None]*3 # Adjacent triangles + #Avoid computing circumcenter here, so the triangle does + #not need to know its coordinates + self.center = None + self.radius = None + + def __repr__(self): + """Dump indexes as a text string""" + return 'tri< ' + str(self.v) + ' >' + +class Delaunay2D: + """ + Class to compute a Delaunay triangulation in 2D + ref: http://en.wikipedia.org/wiki/Bowyer%E2%80%93Watson_algorithm + ref: http://www.geom.uiuc.edu/~samuelp/del_project.html + """ + + def __init__(self, seeds, margin=None): + """Create a new set of triangles from a 2D vertex array + seeds -- Array of 2D coordinates for the vertex + margin -- Optional parameter to set the distance of extended BBox + """ + + seeds = np.asarray(seeds) + + # Compute BBox + smin= np.amin(seeds, axis=0) + smax= np.amax(seeds, axis=0) + #If no margin is given, compute one based on the diagonal + if not margin: + margin = 50*np.linalg.norm(smax- smin) + + # Extend the BBox coordinates by the margin. Store it. + smin -= margin + smax += margin + self.margin = margin + + # Create a two triangles 'frame' big enough to contains all seeds + self.coords = [smin, np.array((smax[0], smin[1])), + smax, np.array((smin[0], smax[1]))] + + T1 = Triangle(0, 3, 1) + T2 = Triangle(2, 1, 3) + T1.neighbour[0] = T2 + T2.neighbour[0] = T1 + T1.center, T1.radius = self.Circumcenter(T1) + T2.center, T2.radius = self.Circumcenter(T2) + self.triangles = [T1, T2] + + #Insert all seeds one by one + for s in seeds: + self.AddPoint(s) + + def Circumcenter(self, tri): + """Compute Circumcenter and circumradius of a triangle. + 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.v]) + A = np.bmat( [[ 2*np.dot(pts,pts.T), np.ones((3,1)) ], + [ np.ones((1,3)) , np.zeros((1,1)) ]] ) + + b = np.hstack((np.sum(pts * pts, axis=1), np.ones((1)))) + x = np.linalg.solve(A,b) + bary_coords = x[:-1] + + center = np.dot(bary_coords,pts) + radius = np.linalg.norm(pts[0] - center) + #print("center:", center, "radius:", radius) + return (center,radius) + + def AddPoint(self, p): + """Add a new point to the current triangulation. + """ + p = np.asarray(p) + idx = len(self.coords) + #print("coords[", idx,"] ->",p) + self.coords.append(p) + + # Search for the triangle(s) too near of p + bad_triangles = [] + for T in self.triangles: + if np.linalg.norm(T.center - p) <= T.radius: + bad_triangles.append(T) + + # Find the convex hull of the bad triangles. + # Expressed as a list of edges (point pairs) in ccw order + boundary = self.Boundary(bad_triangles) + + #Remove triangles too near of point p + for T in bad_triangles: + self.triangles.remove(T) + + # Retriangle the hole left by those triangles + new_triangles = [] + for edge in boundary: + #Create a new Triangle using point p and the edge + T = Triangle(idx, edge[0], edge[1]) + T.center, T.radius = self.Circumcenter(T) + + #Set triangle stored at edge[2] as neighbour of T + T.neighbour[0] = edge[2] + + #Set T as neighbour of triangle stored at edge[2] + if T.neighbour[0]: + #search the reversed edge in T.neighbour[0].v + tmp_v = T.neighbour[0].v + for i in range(3): + if edge[1] == tmp_v[i] and edge[0] == tmp_v[(i+1)%3]: + T.neighbour[0].neighbour[(i+2)%3] = T + + #Add triangle to a temporal list + new_triangles.append(T) + + # Link the new triangles + N = len(new_triangles) + for i, T in enumerate(new_triangles): + T.neighbour[2] = new_triangles[(i-1) % N] # back + T.neighbour[1] = new_triangles[(i+1) % N] # forward + + #Add triangles to our triangulation + self.triangles.extend(new_triangles) + + + def Boundary(self, bad_triangles): + """Build the CCW boundary of a set of triangles. + """ + + boundary = [] + + #Backtracking on each triangle searching edges in the boundary + T = bad_triangles[0] + edge = 0 + while True: + if T.neighbour[edge] in bad_triangles: + #Move to neighbour triangle + last = T + T = T.neighbour[edge] + edge = (T.neighbour.index(last) + 1) % 3 + + else: # Found an edge that is on the boundary + # Add edge to boundary list + boundary.append((T.v[(edge+1)%3], T.v[(edge+2)%3], T.neighbour[edge])) + #Move to next edge in this triangle + edge = (edge + 1) % 3 + + #Check if boundary is a closed loop + if boundary[0][0] == boundary[-1][1]: + break + + return boundary + + def exportDT(self): + """Export the current Delaunay Triangulation. + """ + #Filter out coordinates in the extended BBox + xs = [p[0] for p in self.coords[4:] ] + ys = [p[1] for p in self.coords[4:] ] + + #Filter out triangles with any vertex in the extended BBox + ETS = [t.v for t in self.triangles ] + tris = [(a-4,b-4,c-4) for (a,b,c) in ETS if a > 3 and b > 3 and c > 3] + return xs, ys, tris + + def exportExtendedDT(self): + """Export the Extended Delaunay Triangulation (include the frame vertex). + """ + xs = [p[0] for p in self.coords ] + ys = [p[1] for p in self.coords ] + tris = [t.v for t in self.triangles ] + return xs, ys, tris + diff --git a/delaunay2D_test.py b/delaunay2D_test.py new file mode 100755 index 0000000..5aa2d67 --- /dev/null +++ b/delaunay2D_test.py @@ -0,0 +1,55 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Simple delaunay2D test + +Written by Jose M. Espadero +""" +import numpy as np +import delaunay2D as D + +if __name__ == '__main__': + + ########################################################### + # Generate 'numSeeds' random seeds in a square of size 'radius' + numSeeds = 24 + radius = 100 + seeds = radius * np.random.random((numSeeds, 2)) + print("seeds:\n", seeds) + print("BBox Min:", np.amin(seeds, axis=0), "Bbox Max: ", np.amax(seeds, axis=0) ) + + #Compute our Delaunay triangulation of seeds. + #It is recommended to use the radius to create a noticeable margin + DT = D.Delaunay2D(seeds, 100 * radius) + + #Extract our result (without extended BBox coordinates) + DT_x, DT_y, DT_tris = DT.exportDT() + print("DT_tris:", len(DT_tris), "triangles") + print("DT_tris:\n", DT_tris) + + """ + How to plot triangular grids. (just for debug) + """ + import matplotlib.pyplot as plt + import matplotlib.tri + + # Create a plot with matplotlib.pyplot + fig, ax = plt.subplots() + ax.margins(0.1) + ax.set_aspect('equal') + + #plot our Delaunay triangulation (plot in blue) + ax.triplot(matplotlib.tri.Triangulation(DT_x, DT_y, DT_tris), 'bo-.') + + #DEBUG: Use matplotlib method to create a Delaunay triangulation (plot in green) + #DEBUG: It should be equal to our result in DT_tris (plot in blue) + #DEBUG: If boundary is diferent, try to increase the value of your margin + #ax.triplot(matplotlib.tri.Triangulation(DT_x, DT_y), 'g-.') + + #DEBUG: plot our extended triangulation (plot in red) + #EDT_x, EDT_y, EDT_tris = DT.exportExtendedDT() + #print("Extended DT_tris:", len(EDT_tris), "triangles") + #ax.triplot(matplotlib.tri.Triangulation(EDT_x, EDT_y, EDT_tris), 'ro-.') + + plt.show() +