只要找到三角形顶点的正确顺序,就可以使用plt.plot_trisurf(...)
函数在散点图上绘制3D曲面.SciPy有一个名为ConvexHull的函数,它可以找到数据集外部的点的简化值.这非常方便,但在这个示例中不会立即起作用,因为您的数据集不是凸的!
解决方案是通过将点扩展到远离中心直到它们形成球体来使数据凸起.有关这一点的可视化,请参见下面的内容.
在将头部变成球体后,现在可以调用ConvexHull(...)
来进行所需的三角测量.此三角剖分可以首先应用于球形头部(见下文).然后,头部可以缩小到原来的形式,三角剖分的顶点仍然有效!
这就是最终的产品!
代码
import numpy as np
import matplotlib.pyplot as plt
import csv
from scipy.spatial import KDTree
from scipy.spatial import ConvexHull
from matplotlib import cm
from matplotlib import animation
plt.style.use('dark_background')
# Data reader from a .csv file
def getData(file):
lstX = []
lstY = []
lstZ = []
with open(file, newline='\n') as f:
reader = csv.reader(f, quoting=csv.QUOTE_NONNUMERIC)
for row in reader:
lstX.append(row[0])
lstY.append(row[1])
lstZ.append(row[2])
return lstX, lstY, lstZ
# This function gets rid of the triangles at the base of the neck
# It just filters out any triangles which have at least one side longer than toler
def removeBigTriangs(points, inds, toler=35):
newInds = []
for ind in inds:
if ((np.sqrt(np.sum((points[ind[0]]-points[ind[1]])**2, axis=0))<toler) and
(np.sqrt(np.sum((points[ind[0]]-points[ind[2]])**2, axis=0))<toler) and
(np.sqrt(np.sum((points[ind[1]]-points[ind[2]])**2, axis=0))<toler)):
newInds.append(ind)
return np.array(newInds)
# this calculates the location of each point when it is expanded out to the sphere
def calcSpherePts(points, center):
kdtree = KDTree(points) # tree of nearest points
# d is an array of distances, i is array of indices
d, i = kdtree.query(center, points.shape[0])
spherePts = np.zeros(points.shape, dtype=float)
radius = np.amax(d)
for p in range(points.shape[0]):
spherePts[p] = points[i[p]] *radius /d[p]
return spherePts, i # points and the indices for where they were in the original lists
x,y,z = getData(".\coords3Ddetailed.csv")
pts = np.stack((x,y,z), axis=1)
# generating data
spherePts, sphereInd = calcSpherePts(pts, [0,0,0])
hull = ConvexHull(spherePts)
triangInds = hull.simplices # returns the list of indices for each triangle
triangInds = removeBigTriangs(pts[sphereInd], triangInds)
# plotting!
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter3D(pts[:,0], pts[:,1], pts[:,2], s=2, c='r', alpha=1.0)
ax.plot_trisurf(pts[sphereInd,0], pts[sphereInd,1], pts[sphereInd,2], triangles=triangInds, cmap=cm.Blues, alpha=1.0)
plt.show()