我们从Python开源项目中，提取了以下14个代码示例，用于说明如何使用scipy.spatial.Delaunay()。
def concave_hull(points, alpha, delunay_args=None): """Computes the concave hull (alpha-shape) of a set of points. """ delunay_args = delunay_args or { 'furthest_site': False, 'incremental': False, 'qhull_options': None } triangulation = Delaunay(np.array(points)) alpha_complex = get_alpha_complex(alpha, points, triangulation.simplices) X, Y = [], [] for s in triangulation.simplices: X.append([points[s[k]][0] for k in [0, 1, 2, 0]]) Y.append([points[s[k]][1] for k in [0, 1, 2, 0]]) poly = Polygon(list(zip(X[0], Y[0]))) for i in range(1, len(X)): poly = poly.union(Polygon(list(zip(X[i], Y[i])))) return poly
def scaled_Delaunay(self, points): """ Return a scaled Delaunay mesh and scale factors """ scale_factors = [] points = np.array(points) for i in range(points.shape[1]): scale_factors.append(1.0/np.mean(points[:,i])) points[:,i] = points[:,i]*scale_factors[-1] mesh = Delaunay(points) for i in range(points.shape[1]): mesh.points[:,i] = mesh.points[:,i]/scale_factors[i] return mesh
def scaled_Delaunay(points): """ Return a scaled Delaunay mesh and scale factors """ scale_factors = [] points = np.array(points) for i in range(points.shape[1]): scale_factors.append(1.0/np.mean(points[:,i])) points[:,i] = points[:,i]*scale_factors[-1] mesh = Delaunay(points) for i in range(points.shape[1]): mesh.points[:,i] = mesh.points[:,i]/scale_factors[i] return mesh
def add_delaunay_triangulation_(self): """ Computes the Delaunay triangulation of the protein structure. This has been used in prior work. References: - Harrison, R. W., Yu, X. & Weber, I. T. Using triangulation to include target structure improves drug resistance prediction accuracy. in 1–1 (IEEE, 2013). doi:10.1109/ICCABS.2013.6629236 - Yu, X., Weber, I. T. & Harrison, R. W. Prediction of HIV drug resistance from genotype with encoded three-dimensional protein structure. BMC Genomics 15 Suppl 5, S1 (2014). Notes: 1. We do not use the add_interacting_resis function, because this interaction is computed on the CA atoms. Therefore, there is code duplication. For now, I have chosen to leave this code duplication in. """ ca_coords = self.dataframe[self.dataframe['atom'] == 'CA'] tri = Delaunay(ca_coords[['x', 'y', 'z']]) # this is the triangulation for simplex in tri.simplices: nodes = ca_coords.reset_index().ix[simplex]['node_id'] for n1, n2 in combinations(nodes, 2): if self.has_edge(n1, n2): self.edge[n1][n2]['kind'].add('delaunay') else: self.add_edge(n1, n2, kind={'delaunay'})
def get_alpha_complex(alpha, points, simplexes): """Obtain the alpha shape. Args: alpha (float): the paramter for the alpha shape points: data points simplexes: the list of indices that define 2-simplexes in the Delaunay triangulation """ return filter(lambda simplex: circumcircle(points, simplex)[1] < alpha, simplexes)
def TriAndPaint(img, points, outputIMG): tri = Delaunay(points) triList = points[tri.simplices] cMap = ListedColormap( np.array([ChooseColor(img, tr) for tr in triList])) # use core rgb # center = np.sum(points[tri.simplices], axis=1) / 3 # print(center) # cMap = ListedColormap( # np.array([img.getpixel((x, y)) for x, y in center]) / 256) color = np.array(range(len(triList))) # print(color) width, height = img.size plt.figure(figsize=(width, height), dpi=1) plt.tripcolor(points[:, 0], points[:, 1], tri.simplices.copy(), facecolors=color, cmap=cMap) # plt.tick_params(labelbottom='off', labelleft='off', # left='off', right='off', bottom='off', top='off') # plt.tight_layout(pad=0) plt.axis('off') plt.subplots_adjust(left=0, right=1, top=1, bottom=0) plt.xlim(0, width) plt.ylim(0, height) plt.gca().invert_yaxis() plt.savefig(outputIMG) # uncomment show() if you want to view when it's done # plt.show()
def generate_points(): global triangles,cellsize f = open('67250_5950_25.asc',mode='r') ncols = int(f.readline().split()[1]); nrows = int(f.readline().split()[1]); xllcenter = float(f.readline().split()[1]); yllcenter = float(f.readline().split()[1]); cellsize = float(f.readline().split()[1]); nodata_value = int(f.readline().split()[1]); # print "Columns in data file: ",ncols # print "Rows in data file: ",nrows # print "Coordinate X-wise center (SWEREF99): ",xllcenter # print "Coordinate Y-wise center (SWEREF99): ",yllcenter # print "Cell size in meters: ",cellsize # print "Value if no data available in point: ",nodata_value row = 0 col = 0 index = 0 arr = np.zeros((ncols*nrows,3)); arr = arr.astype(np.float32) for line in f: valarr = line.split() for val in valarr: row = index/ncols col = index%ncols val = (float(val)) arr[index] = [row,col,val] index += 1 #Delaunay triangulation of laser points tri = Delaunay(np.delete(arr,2,1)) triangles = arr[tri.simplices] triangles = np.vstack(triangles)
def __init__(self, nx, ny): """ :param nx: The number of cells in the x direction. :param ny: The number of cells in the y direction. """ points = list((x, y) for x in np.linspace(0, 1, nx + 1) for y in np.linspace(0, 1, ny + 1)) mesh = Delaunay(points) super(UnitSquareMesh, self).__init__(mesh.points, mesh.simplices)
def inner_points_mask(points): """Mask array into `points` where ``points[msk]`` are all "inner" points, i.e. `points` with one level of edge points removed. For 1D, this is simply points[1:-1,:] (assuming ordered points). For ND, we calculate and remove the convex hull. Parameters ---------- points : nd array (npoints, ndim) Returns ------- msk : (npoints, ndim) Bool array. """ msk = np.ones((points.shape[0],), dtype=bool) if points.shape[1] == 1: assert (np.diff(points[:,0]) >= 0.0).all(), ("points not monotonic") msk[0] = False msk[-1] = False else: from scipy.spatial import Delaunay tri = Delaunay(points) edge_idx = np.unique(tri.convex_hull) msk.put(edge_idx, False) return msk
def test_read_delaunay(): xs = np.linspace(0, 0.5, 3) ys = np.linspace(0, 0.2, 3) points = np.array(np.meshgrid(xs, ys)).T.reshape(-1, 2) tri = Delaunay(points) mesh = read_delaunay(points, tri) assert len(mesh.nodes) == 9 assert len(mesh.elements) == 8 assert len(mesh.edges) == 16
def delaunay_triangulation(points, plot=False): """ Extract a Delaunay's triangulation from the points """ tri = Delaunay(points) if plot: plt.triplot(points[:, 0], points[:, 1], tri.simplices.copy()) plt.plot(points[:, 0], points[:, 1], 'o') plt.show() return tri.simplices
def wet_circles(A, B, thetaA, thetaB): """Generates a mesh that wets the surface of circles A and B. Parameters ------------- A,B : Circle theta : list the number of radians that the wet covers and number of the points on the surface range """ vector = B.center - A.center if vector.x > 0: angleA = np.arctan(vector.y/vector.x) angleB = np.pi + angleA else: angleB = np.arctan(vector.y/vector.x) angleA = np.pi + angleB # print(vector) rA = A.radius rB = B.radius points = [] for t in ((np.arange(0, thetaA[1])/(thetaA[1]-1) - 0.5) * thetaA[0] + angleA): x = rA*np.cos(t) + A.center.x y = rA*np.sin(t) + A.center.y points.append([x, y]) mid = len(points) for t in ((np.arange(0, thetaB[1])/(thetaB[1]-1) - 0.5) * thetaB[0] + angleB): x = rB*np.cos(t) + B.center.x y = rB*np.sin(t) + B.center.y points.append([x, y]) points = np.array(points) # Triangulate the polygon tri = Delaunay(points) # Remove extra triangles # print(tri.simplices) mask = np.sum(tri.simplices < mid, 1) mask = np.logical_and(mask < 3, mask > 0) tri.simplices = tri.simplices[mask, :] # print(tri.simplices) m = Mesh() for t in tri.simplices: m.append(Triangle(Point([points[t[0], 0], points[t[0], 1]]), Point([points[t[1], 0], points[t[1], 1]]), Point([points[t[2], 0], points[t[2], 1]]))) return m
def delaunay3d(pyptlist, tolerance = 1e-06): """ This function constructs a mesh (list of triangle OCCfaces) from a list of points. Currently only works for two dimensional points. Parameters ---------- pyptlist : a list of tuples The list of points to be triangulated. A pypt is a tuple that documents the xyz coordinates of a pt e.g. (x,y,z), thus a pyptlist is a list of tuples e.g. [(x1,y1,z1), (x2,y2,z2), ...] tolerance : float, optional The minimal surface area of each triangulated face, Default = 1e-06.. Any faces smaller than the tolerance will be deleted. Returns ------- list of face : list of OCCfaces A list of meshed OCCfaces (triangles) constructed from the meshing. """ import numpy as np from scipy.spatial import Delaunay pyptlistx = [] pyptlisty = [] pyptlistz = [] for pypt in pyptlist: pyptlistx.append(pypt[0]) pyptlisty.append(pypt[1]) pyptlistz.append(pypt[2]) # u, v are parameterisation variables u = np.array(pyptlistx) v = np.array(pyptlisty) x = u y = v z = np.array(pyptlistz) # Triangulate parameter space to determine the triangles tri = Delaunay(np.array([u,v]).T) occtriangles = [] xyz = np.array([x,y,z]).T for verts in tri.simplices: pt1 = list(xyz[verts[0]]) pt2 = list(xyz[verts[1]]) pt3 = list(xyz[verts[2]]) occtriangle = make_polygon([pt1,pt2,pt3]) tri_area = calculate.face_area(occtriangle) if tri_area > tolerance: occtriangles.append(occtriangle) return occtriangles #======================================================================================================== #EDGE INPUTS #========================================================================================================
def warp_image(img, triangulation, base_points, coord): """ Realize the mesh warping phase triangulation is the Delaunay triangulation of the base points base_points are the coordinates of the landmark poitns of the reference image code inspired from http://www.learnopencv.com/warp-one-triangle-to-another-using-opencv-c-python/ """ all_points, coordinates = preprocess_image_before_triangulation(img) img_out = 255 * np.ones(img.shape, dtype=img.dtype) for t in triangulation: # triangles to map one another src_tri = np.array([[all_points[x][0], all_points[x][1]] for x in t]).astype(np.float32) dest_tri = np.array([[base_points[x][0], base_points[x][1]] for x in t]).astype(np.float32) # bounding boxes src_rect = cv2.boundingRect(np.array([src_tri])) dest_rect = cv2.boundingRect(np.array([dest_tri])) # crop images src_crop_tri = np.zeros((3, 2), dtype=np.float32) dest_crop_tri = np.zeros((3, 2)) for k in range(0, 3): for dim in range(0, 2): src_crop_tri[k][dim] = src_tri[k][dim] - src_rect[dim] dest_crop_tri[k][dim] = dest_tri[k][dim] - dest_rect[dim] src_crop_img = img[src_rect[1]:src_rect[1] + src_rect[3], src_rect[0]:src_rect[0] + src_rect[2]] # affine transformation estimation mat = cv2.getAffineTransform( np.float32(src_crop_tri), np.float32(dest_crop_tri) ) dest_crop_img = cv2.warpAffine( src_crop_img, mat, (dest_rect[2], dest_rect[3]), None, flags=cv2.INTER_LINEAR, borderMode=cv2.BORDER_REFLECT_101 ) # Use a mask to keep only the triangle pixels # Get mask by filling triangle mask = np.zeros((dest_rect[3], dest_rect[2], 3), dtype=np.float32) cv2.fillConvexPoly(mask, np.int32(dest_crop_tri), (1.0, 1.0, 1.0), 16, 0) # Apply mask to cropped region dest_crop_img = dest_crop_img * mask # Copy triangular region of the rectangular patch to the output image img_out[dest_rect[1]:dest_rect[1] + dest_rect[3], dest_rect[0]:dest_rect[0] + dest_rect[2]] = \ img_out[dest_rect[1]:dest_rect[1] + dest_rect[3], dest_rect[0]:dest_rect[0] + dest_rect[2]] * ( (1.0, 1.0, 1.0) - mask) img_out[dest_rect[1]:dest_rect[1] + dest_rect[3], dest_rect[0]:dest_rect[0] + dest_rect[2]] = \ img_out[dest_rect[1]:dest_rect[1] + dest_rect[3], dest_rect[0]:dest_rect[0] + dest_rect[2]] + dest_crop_img return img_out[coord[2]:coord[3], coord[0]:coord[1]]
