表面曲率(平均、ガウス、および原理曲率)の値を表面にマッピングしようとしています。私は、人工的に生成された3Dサーフェス(例:円筒)の曲率値を計算しました。私が得ようとしている3D表面は、mean curvature mapped to surfaceのようなものです。誰かがこれを得る方法で私を導くことができますか?表面の曲率をマッピングするか?
私が作成しています表面のためのコードは次のとおり
import math
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
xindex = []
yindex = []
zindex = []
x = []
y = []
z = []
count = 1
surfaceSt = []
import numpy
numpy.set_printoptions(threshold=numpy.nan)
#surfaceStX = numpy.empty((10,36))
#surfaceStY = numpy.empty((10,36))
#surfaceStZ = numpy.empty((10,36))
surfaceStZ = []
surfaceStX = []
surfaceStY = []
for i in range(1,21):
if i < 11:
x = []
y = []
z = []
pt = []
ptX = []
ptY = []
ptZ = []
for t in range(0,360,10):
x = i*math.sin(math.radians(t))
y = i*math.cos(math.radians(t))
z = i-1
ptX.append(x)
ptY.append(y)
ptZ.append(z)
pt.append([x,y,z])
ptX.append(ptX[0])
ptY.append(ptY[0])
ptZ.append(ptZ[0])
surfaceStX.append(ptX)
surfaceStY.append(ptY)
surfaceStZ.append(ptZ)
# numpy.append(surfaceStX,ptX)
# numpy.append(surfaceStY,ptY)
# numpy.append(surfaceStZ,ptZ)
#ax.scatter(x,y,z)
elif i >= 11:
x = []
y = []
z = []
pt = []
ptX = []
ptY = []
ptZ = []
for t in range(0,360,10):
x = (i-count)*math.sin(math.radians(t))
y = (i-count)*math.cos(math.radians(t))
z = i-1
ptX.append(x)
ptY.append(y)
ptZ.append(z)
pt.append([x,y,z])
ptX.append(ptX[0])
ptY.append(ptY[0])
ptZ.append(ptZ[0])
surfaceStX.append(ptX)
surfaceStY.append(ptY)
surfaceStZ.append(ptZ)
count +=2
X = numpy.array(surfaceStX)
Y = numpy.array(surfaceStY)
Z = numpy.array(surfaceStZ)
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, rstride=1, cstride=1,shade = 'True')
from surfaceCurvature import surface_curvature
Pcurvature,Gcurvature,Mcurvature = surface_curvature(X,Y,Z)
plt.show()
マイ表面曲率計算を(礼儀:https://github.com/sujithTSR/surface-curvature):以下に示す。
def surface_curvature(X,Y,Z):
\t (lr,lb)=X.shape
\t #print lr
\t #print "awfshss-------------"
\t #print lb
#First Derivatives
\t Xv,Xu=np.gradient(X)
\t Yv,Yu=np.gradient(Y)
\t Zv,Zu=np.gradient(Z)
#Second Derivatives
\t Xuv,Xuu=np.gradient(Xu)
\t Yuv,Yuu=np.gradient(Yu)
\t Zuv,Zuu=np.gradient(Zu)
\t Xvv,Xuv=np.gradient(Xv)
\t Yvv,Yuv=np.gradient(Yv)
\t Zvv,Zuv=np.gradient(Zv)
#2D to 1D conversion
#Reshape to 1D vectors
\t Xu=np.reshape(Xu,lr*lb)
\t Yu=np.reshape(Yu,lr*lb)
\t Zu=np.reshape(Zu,lr*lb)
\t Xv=np.reshape(Xv,lr*lb)
\t Yv=np.reshape(Yv,lr*lb)
\t Zv=np.reshape(Zv,lr*lb)
\t Xuu=np.reshape(Xuu,lr*lb)
\t Yuu=np.reshape(Yuu,lr*lb)
\t Zuu=np.reshape(Zuu,lr*lb)
\t Xuv=np.reshape(Xuv,lr*lb)
\t Yuv=np.reshape(Yuv,lr*lb)
\t Zuv=np.reshape(Zuv,lr*lb)
\t Xvv=np.reshape(Xvv,lr*lb)
\t Yvv=np.reshape(Yvv,lr*lb)
\t Zvv=np.reshape(Zvv,lr*lb)
\t Xu=np.c_[Xu, Yu, Zu]
\t Xv=np.c_[Xv, Yv, Zv]
\t Xuu=np.c_[Xuu, Yuu, Zuu]
\t Xuv=np.c_[Xuv, Yuv, Zuv]
\t Xvv=np.c_[Xvv, Yvv, Zvv]
# First fundamental Coeffecients of the surface (E,F,G)
\t
\t E=np.einsum('ij,ij->i', Xu, Xu)
\t F=np.einsum('ij,ij->i', Xu, Xv)
\t G=np.einsum('ij,ij->i', Xv, Xv)
\t m=np.cross(Xu,Xv,axisa=1, axisb=1)
\t p=np.sqrt(np.einsum('ij,ij->i', m, m))
\t n=m/np.c_[p,p,p]
# Second fundamental Coeffecients of the surface (L,M,N), (e,f,g)
\t L= np.einsum('ij,ij->i', Xuu, n) #e
\t M= np.einsum('ij,ij->i', Xuv, n) #f
\t N= np.einsum('ij,ij->i', Xvv, n) #g
# Gaussian Curvature
\t K=(L*N-M**2)/(E*G-F**2)
\t K=np.reshape(K,lr*lb)
# Mean Curvature
\t H = (E*N + G*L - 2*F*M)/((E*G - F**2))
\t H = np.reshape(H,lr*lb)
# Principle Curvatures
\t Pmax = H + np.sqrt(H**2 - K)
\t Pmin = H - np.sqrt(H**2 - K)
#[Pmax, Pmin]
\t Principle = [Pmax,Pmin]
\t return Principle,K,H
EDIT 1:
私はarmatitaによって提供されたリンクに基づいていくつか試しました。以下の私のコードです:
'''
Creat half cylinder
'''
import numpy
import matplotlib.pyplot as plt
import math
ptX= []
ptY = []
ptZ = []
ptX1 = []
ptY1 = []
ptZ1 = []
for i in range(0,10):
x = []
y = []
z = []
for t in range(0,200,20):
x.append(10*math.cos(math.radians(t)))
y.append(10*math.sin(math.radians(t)))
z.append(i)
x1= 5*math.cos(math.radians(t))
y1 = 5*math.sin(math.radians(t))
z1 = i
ptX1.append(x1)
ptY1.append(y1)
ptZ1.append(z1)
ptX.append(x)
ptY.append(y)
ptZ.append(z)
X = numpy.array(ptX)
Y = numpy.array(ptY)
Z = numpy.array(ptZ)
fig = plt.figure()
ax = fig.add_subplot(111,projection = '3d')
from surfaceCurvature import surface_curvature
p,g,m= surface_curvature(X,Y,Z)
n = numpy.reshape(m,numpy.shape(X))
ax.plot_surface(X,Y,Z, rstride=1, cstride=1)
plt.show()
'''
Map mean curvature to color
'''
import numpy as np
X1 = X.ravel()
Y1 = Y.ravel()
Z1 = Z.ravel()
from scipy.interpolate import RectBivariateSpline
# Define the points at the centers of the faces:
y_coords, x_coords = np.unique(Y1), np.unique(X1)
y_centers, x_centers = [ arr[:-1] + np.diff(arr)/2 for arr in (y_coords, x_coords)]
# Convert back to a 2D grid, required for plot_surface:
#Y1 = Y.reshape(y_coords.size, -1)
#X1 = X.reshape(-1, x_coords.size)
#Z1 = Z.reshape(X.shape)
C = m.reshape(X.shape)
C -= C.min()
C /= C.max()
interp_func = RectBivariateSpline(x_coords, y_coords, C.T, kx=1, ky=1)
I get the following error:
raise TypeError('y dimension of z must have same number of y')
TypeError: y dimension of z must have same number of elements as y
すべての寸法が同じです。私の実装で何がうまくいかないのか誰にでも教えてくれますか?
変数をカラーにマップすることができます(関数は3以上を返す)。このリンクを確認してください:http://stackoverflow.com/questions/36614456/how-do-i-color-individual-sections-of-a-3d-sphere-in-python/36616730#36616730 – armatita
コメントありがとうございます。私はあなたが提案し、投稿を編集しようとしました。私はまだそれを働かせることはできません。もっと助けていただければ幸いです。 –