How can I plot the surface of a structure which is given by vectors in python?

Question

I would like to plot the surface of my data which is given by 3D vectors in cartesian coordinates x,y,z. The data can not be represented by a smooth function.

So first we generate some dummy data with the function eq_points(N_count, r) which returns an array points with the x,y,z coordinates of each point on the surface of our object. The quantity omega is the solid angle, and not of interest right now.

#credit to Markus Deserno from MPI
#https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
def eq_points(N_count, r):
    points = []
    a = 4*np.pi*r**2/N_count
    d = np.sqrt(a)
    M_theta = int(np.pi/d)
    d_theta = np.pi/M_theta
    d_phi = a/d_theta
    for m in range(M_theta):
        theta = np.pi*(m+0.5)/M_theta
        M_phi = int(2*np.pi*np.sin(theta)/d_phi)
        for n in range(M_phi):
            phi = 2*np.pi*n/M_phi

            points.append(np.array([r*np.sin(theta)*np.cos(phi),
                                    r*np.sin(theta)*np.sin(phi),
                                    r*np.cos(theta)]))

    omega = 4*np.pi/N_count

    return np.array(points), omega

#starting plotting sequence

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

points, omega = eq_points(400, 1.)

ax.scatter(points[:,0], points[:,1], points[:,2])
ax.scatter(0., 0., 0., c="r")
ax.set_xlabel(r'$x$ axis')
ax.set_ylabel(r'$y$ axis')
ax.set_zlabel(r'$Z$ axis')

plt.savefig("./sphere.png", format="png", dpi=300)
plt.clf()

The result is a sphere shown in the following figure. The blue points mark the data from the points array, while the red point is the origin.

I would like to get something like this taken from here. However the data in the mplot3d tutorial is always a result of a smooth function. Except to the ax.scatter() function which I used for my sphere plot.

So in the end my goal would be to plot some data showing only its surface. This data is produced by changing the radial distance to the origin of each blue point. Further more it would be necessary to make sure each point is in contact with the surface. How are the surfaces which are plotted here e.g. in plot_surface() constructed in detail? Some actual live data looks like this:

Answer 1

Solution to the question with the new specification that all points are touching the surface. Assuming that the angles are set by the user as shown in the example, it is easy to precompute the indices of the points forming the simplices making up the surface by computing the simplices of the hull formed by points on the unit sphere with the same angles as in the data set of interest. We can then use these indices to get the surface of interest.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from scipy.spatial import ConvexHull

def eq_points(N_count, r):
    points = []
    a = 4*np.pi*r**2/N_count
    d = np.sqrt(a)
    M_theta = int(np.pi/d)
    d_theta = np.pi/M_theta
    d_phi = a/d_theta
    for m in range(M_theta):
        theta = np.pi*(m+0.5)/M_theta
        M_phi = int(2*np.pi*np.sin(theta)/d_phi)
        for n in range(M_phi):
            phi = 2*np.pi*n/M_phi

            points.append(np.array([r*np.sin(theta)*np.cos(phi),
                                    r*np.sin(theta)*np.sin(phi),
                                    r*np.cos(theta)]))

    omega = 4*np.pi/N_count

    return np.array(points), omega

def eq_points_with_random_radius(N_count, r):
    points = []
    a = 4*np.pi*r**2/N_count
    d = np.sqrt(a)
    M_theta = int(np.pi/d)
    d_theta = np.pi/M_theta
    d_phi = a/d_theta
    for m in range(M_theta):
        theta = np.pi*(m+0.5)/M_theta
        M_phi = int(2*np.pi*np.sin(theta)/d_phi)
        for n in range(M_phi):
            phi = 2*np.pi*n/M_phi
            rr = r * np.random.rand()
            points.append(np.array([rr*np.sin(theta)*np.cos(phi),
                                    rr*np.sin(theta)*np.sin(phi),
                                    rr*np.cos(theta)]))

    omega = 4*np.pi/N_count

    return np.array(points), omega


N = 400
pts, _ = eq_points(N, 1.)
pts_rescaled, _ = eq_points_with_random_radius(N, 1.)
extremum = 2.

# plot points
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(pts_rescaled[:,0], pts_rescaled[:,1], pts_rescaled[:,2])
ax.set_xlim(-extremum, extremum)
ax.set_ylim(-extremum, extremum)
ax.set_zlim(-extremum, extremum)

# get indices of simplices making up the surface using points on unit sphere;
# index into rescaled points  
hull = ConvexHull(pts)
vertices = [pts_rescaled[s] for s in hull.simplices]

fig = plt.figure()
ax = Axes3D(fig)
triangles = Poly3DCollection(vertices, edgecolor='k')
ax.add_collection3d(triangles)
ax.set_xlim(-extremum, extremum)
ax.set_ylim(-extremum, extremum)
ax.set_zlim(-extremum, extremum)
plt.show()