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:

The link to the MatLab tutorial is diffuse. Can you specify a section within the tutorial? Do you mean https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#d-plots-in-3d, for instance? — brannerchinese, Dec 07 '17 at 13:33
Thanks for the clarification. I suggest changing the link at "here" to https://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html#tri-surface-plots. — brannerchinese, Dec 07 '17 at 13:38
I would like to add pictures and my results, but I cant due to the lack of reputation so far — andre.klptk, Dec 07 '17 at 13:42
Reputation accrues before you know it, and you can return to older questions and improve them with (for instance) images and detail. Patience, colleague! — brannerchinese, Dec 07 '17 at 13:43
I have a problem understanding how that surface should be created for the general case. Even in two dimension there would be a lot of different options to draw a path (see e.g. the image of [this answer](https://stackoverflow.com/questions/44025403/how-to-use-matplotlib-path-to-draw-polygon)). Which one would you choose here? I guess once you have a clear idea of that, a solution with matplotlib will be possible. — ImportanceOfBeingErnest, Dec 07 '17 at 21:09
Can you guarantee that no two points have the same angle wrt to the origin? — Paul Brodersen, Dec 08 '17 at 09:35
@Paul yes, I can guarantee this. The data is produced by altering the radii of the `points` array yielded by `eq_points()` function with positive floats. So there occurs no flip (rotation around 180° around an arbitrary axis). I am really sorry for underspecfiying this. — andre.klptk, Dec 08 '17 at 09:52
Ok, in this case, you are probably looking for the solution with the minimum surface, which is easy enough. Pseudocode: for each point, find the 6 closest nearest neighbours (in angles), and get the 6 triangles that make up the hexagon formed by those neighbours. Plot the set of triangles as outlined in my solution. — Paul Brodersen, Dec 08 '17 at 10:22
Also, if you are setting those points/angles, you can just precompute the hull as I do below (which gives you the **indices** of the points forming each simplex), and then apply your scaling to the points, then plot the simplices by indexing into your rescaled points. If post a set of points on the unit sphere and a corresponding set of points with rescaled radius then I can demonstrate how to achieve this. — Paul Brodersen, Dec 08 '17 at 10:37

Paul Brodersen · Answer 1 · 2017-12-07T17:49:55.243

I would suggest finding the hull, and then plotting the simplices (i.e. the triangles forming the hull). Make sure to update the x,y,z-limits appropriately.

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

N = 1000
pts = np.random.randn(N, 3)

# exclude outliers
# obviously, this is data dependent
cutoff = 3.
is_outlier = np.any(np.abs(pts) > cutoff, axis=1)
pts = pts[~is_outlier]

# plot points
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(pts[:,0], pts[:,1], pts[:,2])

ax.set_xlim(-(cutoff +1), cutoff+1)
ax.set_ylim(-(cutoff +1), cutoff+1)
ax.set_zlim(-(cutoff +1), cutoff+1)

# get and plot hull
hull = ConvexHull(pts)

fig = plt.figure()
ax = Axes3D(fig)
vertices = [pts[s] for s in hull.simplices]
triangles = Poly3DCollection(vertices, edgecolor='k')
ax.add_collection3d(triangles)

ax.set_xlim(-(cutoff +1), cutoff+1)
ax.set_ylim(-(cutoff +1), cutoff+1)
ax.set_zlim(-(cutoff +1), cutoff+1)
plt.show()

Thanks for the answer. How do I make sure every point of the data is in contact with the surface? This is not solved by a convex hull surface. — andre.klptk, Dec 07 '17 at 19:51
The original question did not specify that each point should be in contact with the surface. If you think you need to change the question substantially after someone has given an answer, then you should really raise a new one instead. Also, if you don't want the hull, I agree with Ernest that the problem is then underspecified. — Paul Brodersen, Dec 08 '17 at 09:33

score 1 · Accepted Answer · answered Dec 08 '17 at 10:57

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

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

2 Answers2