plot 3D connected prism matplotlib based on vertices

Question

I have a set of coordinates x, y, z from 4-6 vectors. I want to plot the corresponding prism. But my lines are crossed and do not look like a prism at all.

I assume I have to sort my dataset, but I am not sure how or if this is the correct answer.

My plot

This is clearly erroneous:

Ideal prism

And this how it should look like:

Real prism for a 3 coordinates dataset

Script

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
import numpy as np

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
p = np.array([[-43.11150999, -118.14365791, -1100.99389988],
                [-27.97693445,-124.54828379, -1089.54038197],
                [-55.99892873, -120.42384095, -1084.32576297],
                [-40.75143664, -133.41566716, -1077.33745869],

              [-43.2165748, -34.770722, -1030.85272686],
              [-27.89568594, -43.06953117, -1021.03437003],
              [-56.072327, -44.66085799, -1019.15166512],
              [-40.75143814, -52.95966716, -1009.3333083]])

ax.scatter3D(p[:, 0], p[:, 1], p[:, 2])

x = np.array([[-43.11150999], [-27.97693445], [-55.99892873], [-40.75143664], [-43.2165748],[-27.89568594],[-56.072327],[-40.75143814]])
y = np.array([[-118.14365791], [-124.54828379], [-120.42384095], [-133.41566716], [-34.770722],[-43.06953117],[-44.66085799],[-52.95966716]])
z = np.array([[-1100.99389988], [-1089.54038197], [-1084.32576297], [-1077.33745869], [-1030.85272686],[-1021.03437003],[-1019.15166512],[-1009.3333083]])

labels = ['PT-EP-1n', 'PT-EP-2n', 'PT-EP-3n', 'PT-EP-4n', 'PT-TP-1n','PT-TP-2n','PT-TP-3n','PT-TP-4n']

x = x.flatten()
y = y.flatten()
z = z.flatten()

ax.scatter(x, y, z)
#give the labels to each point
for x, y, z, label in zip(x, y,z, labels):
    ax.text(x, y, z, label)

verts = [[p[0],p[1],p[2],p[3]],
          [p[1],p[2],p[6],p[5]],
          [p[2],p[3],p[7],p[6]],
          [p[3],p[0],p[4],p[7]],
          [p[0],p[1],p[5],p[4]],
          [p[4],p[5],p[6],p[7]]]

collection = Poly3DCollection(verts, linewidths=1, edgecolors='black', alpha=0.2, zsort='min')
face_color = "salmon"
collection.set_facecolor(face_color)
ax.add_collection3d(collection)

plt.show()

Answer 1

A very well-asked first question!

I think what you are looking for is the convex hull of your data points, which can be computed using scipy.spatial.ConvexHull. The problem with this approach is, however, that this function returns a set of triangles, which will not correspond to the set of faces of your prism in a 1-to-1 fashion. Rather, multiple co-planar triangles will often form a single face.

In a second step, you will hence need to simplify adjacent, co-planar triangles into a single face. If your vertices come from experimental data, then you might run into trouble in this step as the triangles may not actually be co-planar, and then the naive simplify procedure given below will fail.

Borrowing from @ImportanceOfBeingErnest's answer here (which is much simpler and faster than the accepted answer), you get:

#!/usr/bin/env python
# coding: utf-8

import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as a3

from scipy.spatial import ConvexHull

class Faces():
    def __init__(self,tri, sig_dig=12, method="convexhull"):
        self.method=method
        self.tri = np.around(np.array(tri), sig_dig)
        self.grpinx = list(range(len(tri)))
        norms = np.around([self.norm(s) for s in self.tri], sig_dig)
        _, self.inv = np.unique(norms,return_inverse=True, axis=0)

    def norm(self,sq):
        cr = np.cross(sq[2]-sq[0],sq[1]-sq[0])
        return np.abs(cr/np.linalg.norm(cr))

    def isneighbor(self, tr1,tr2):
        a = np.concatenate((tr1,tr2), axis=0)
        return len(a) == len(np.unique(a, axis=0))+2

    def order(self, v):
        if len(v) <= 3:
            return v
        v = np.unique(v, axis=0)
        n = self.norm(v[:3])
        y = np.cross(n,v[1]-v[0])
        y = y/np.linalg.norm(y)
        c = np.dot(v, np.c_[v[1]-v[0],y])
        if self.method == "convexhull":
            h = ConvexHull(c)
            return v[h.vertices]
        else:
            mean = np.mean(c,axis=0)
            d = c-mean
            s = np.arctan2(d[:,0], d[:,1])
            return v[np.argsort(s)]

    def simplify(self):
        for i, tri1 in enumerate(self.tri):
            for j,tri2 in enumerate(self.tri):
                if j > i:
                    if self.isneighbor(tri1,tri2) and \
                       self.inv[i]==self.inv[j]:
                        self.grpinx[j] = self.grpinx[i]
        groups = []
        for i in np.unique(self.grpinx):
            u = self.tri[self.grpinx == i]
            u = np.concatenate([d for d in u])
            u = self.order(u)
            groups.append(u)
        return groups


if __name__ == '__main__':

    x = np.array([[-43.11150999], [-27.97693445], [-55.99892873], [-40.75143664], [-43.2165748],[-27.89568594],[-56.072327],[-40.75143814]])
    y = np.array([[-118.14365791], [-124.54828379], [-120.42384095], [-133.41566716], [-34.770722],[-43.06953117],[-44.66085799],[-52.95966716]])
    z = np.array([[-1100.99389988], [-1089.54038197], [-1084.32576297], [-1077.33745869], [-1030.85272686],[-1021.03437003],[-1019.15166512],[-1009.3333083]])

    verts = np.c_[x, y, z]

    # compute the triangles that make up the convex hull of the data points
    hull = ConvexHull(verts)
    triangles = [verts[s] for s in hull.simplices]

    # combine co-planar triangles into a single face
    faces = Faces(triangles, sig_dig=1).simplify()

    # plot
    ax = a3.Axes3D(plt.figure())
    pc = a3.art3d.Poly3DCollection(faces,
                                   facecolor="salmon",
                                   edgecolor="k", alpha=0.9)
    ax.add_collection3d(pc)

    # define view
    ax.set_xlim(np.min(x), np.max(x))
    ax.set_ylim(np.min(y), np.max(y))
    ax.set_zlim(np.min(z), np.max(z))
    ax.dist=10
    ax.azim=30
    ax.elev=10

    plt.show()