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
python类Delaunay()的实例源码
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]]
