[FIXED] How can one plot a 3D surface in matplotlib by points' coordinates?

Issue

After awhole day of searching, in desperation, I address to you, my dear fellows. I want to draw a 3D surface of human head, for which I have found nice of 3D coordinates (you can download it from my Google drive here). Using 3D scatter plot, everything looks beautiful:

For my further purposes, I'd like to plot it as a 3D surface as it looks like in real life. As far as I could conclude, the best way to do so in matplotlib is plot_trisurf function, where one can explicitly put ther points' coordinates and enjoy the result, however, for the same coordinates as in the scatter plot above, I recieve the following ugly result:

The problem seems to be that triangles are plotted by random triplets of points rather than the closest ones.

I tried to sort the points minimizing the distance between neighbouring ones, which had no effect. I also tried to reproduce the meshgrid-based surface plotting that is basically used in all of the matplotlib examples. None of that worked. I also considered using other libraries, but most of them are either matplotlib-based or interactive. The latter is an overkill, since I'm going to use it in a realtime calculations, so sticking to the matplotlib API is prioritized. The idea seems pretty basic: plot a surface by 3D coordinates by dyeing the triangles of three closest points, although I didn't manage to find such a function. Hopefully, any of you know what can be done to overcome this issue. Thanks in advance.

Solution

It is possible to plot the 3D surface over your scatter plot using the plt.plot_trisurf(...) function as long as you find the right ordering of vertices for the triangles. There is a function from SciPy called ConvexHull which finds the simplices of the points on the outside of the data set. This is very handy, but does not immediately work on this example because your data set is not convex!

The solution is to make the data convex by expanding the points away from the center until they form a sphere. See below for a visualization of this.

After turning the head into a sphere it is now possible to call ConvexHull(...) to get the desired triangulation. This triangulation can be applied to the spherical head first (see below). Then, the head can be shrunk back into its original form, and the triangulation's vertices are still valid!

This is the final product!

Code

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()

Answered By - trent