CODEWITHSUNDEEP
pythonmatplotlibplot
Dipole
Dipole

Reputation: 1910

Create surface grid from point cloud data in Python

Here is an example creating a point cloud which I then want to fit a grided surface to. The problem comes at the end when I try to pass in meshgrid arrays to a function which interpolated the data:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

# Create some point cloud data:
c = 1
a = 3
b = 4

slice = {}
t = np.linspace(0,2*np.pi,50)

for s in np.linspace(1,9,10):
    c = 5*s
    r = (-s**2+10.0*s)/10.0
    X = r*np.cos(t)
    Y = r*np.sin(t)
    Z = c*(Y**2/b**2 - X**2/a**2) + c
    slice[str(int(s))] = np.vstack([X,Y,Z])


# Visualize it:

fig = plt.figure()
ax = fig.gca(projection = '3d')

for k,v in slice.iteritems():
    print type(v)
    print np.shape(v)
    X = v[0,:]
    Y = v[1,:]
    Z = v[2,:]
    ax.scatter(X,Y,Z)
plt.show()

It looks like this:

enter image description here

I now need to create a surface mesh based on these points. There are multiple interpretations of surface in this case because I just have a point cloud rather than a function z = f(x,y) but the correct surface in this case should be the one that creates a hollow "warped cylinder". I thought of attacking the problem like this:

# stack all points from each slice into one vector for each coordinate:
Xs = []
Ys = []
Zs = []
for k,v in slice.iteritems():
    #ax.plot_surface(X,Y,Z)
    X = v[0,:]
    Y = v[1,:]
    Z = v[2,:]
    Xs = np.hstack((Xs,X))
    Ys = np.hstack((Ys,Y))
    Zs = np.hstack((Zs,Z))

XX, YY = np.meshgrid(Xs,Ys)

from scipy import interpolate
f = interpolate.interp2d(Xs,Ys,Zs, kind = 'cubic')
ZZ = f(XX,YY)

which I would then be able to plot using

fig = plt.figure()
ax = fig.gca(projection = '3d')

ax.plot_surface(XX, YY, ZZ)
plt.show()

However the interpolated function does not seem to accept arrays as inputs so this method might not work. Can anyone come up with a suggestion on how to do this properly?

Edit:

Actually the data is obviously not able to be represented as a function as it would not be one to one.

Upvotes: 5

Views: 4206

Answers (1)

zufall
zufall

Reputation: 174

I stumbled upon the same question and wondered why it has not been solved in the last 7 years. Here's my solution for any future reader based on plot_trisurf (and the corresponding code examples).

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib.tri as mtri

# Create some point cloud data:
a = 3
b = 4

# def grid of parametric variables
u = np.linspace(0,2*np.pi,50)
v = np.linspace(1,9,50)
U, V = np.meshgrid(u, v)
U, V = U.flatten(), V.flatten()

# Triangulate parameter space to determine the triangles
tri = mtri.Triangulation(U, V)

# get the transformed data as list
X,Y,Z = [],[],[]
for _u in u:
  for _v in v:
    r = (-_v**2+10.0*_v)/10.0
    x = r*np.cos(_u)
    y = r*np.sin(_u)
    z = 5*_v*(y**2/b**2 - x**2/a**2) + 5*_v
    X.append(x)
    Y.append(y)
    Z.append(z)

# Visualize it:
fig = plt.figure()
ax = fig.gca(projection = '3d')
ax.scatter(X,Y,Z, s=1, c='r')
ax.plot_trisurf(X, Y, Z, triangles=tri.triangles, alpha=.5)
plt.show()

This produces the following plot. enter image description here

Upvotes: 2

Related Questions