2016-04-28 6 views
0

表面曲率(平均、ガウス、および原理曲率)の値を表面にマッピングしようとしています。私は、人工的に生成された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 

すべての寸法が同じです。私の実装で何がうまくいかないのか誰にでも教えてくれますか?

+0

変数をカラーにマップすることができます(関数は3以上を返す)。このリンクを確認してください:http://stackoverflow.com/questions/36614456/how-do-i-color-individual-sections-of-a-3d-sphere-in-python/36616730#36616730 – armatita

+0

コメントありがとうございます。私はあなたが提案し、投稿を編集しようとしました。私はまだそれを働かせることはできません。もっと助けていただければ幸いです。 –

答えて

0

私はあなたが必要とするものを正確に把握する必要があると思います。あなたのコードを見て、私はあなたが使用していない変数を生成していることに気付きました。また、あなたは表面曲率を計算する機能を持っているようですが、目的を見ることができないnp.unique関数を使っていくつかの計算をしようとするよりも、そのエラーが表示される理由です。

それでは、これを想定してみましょう:

  • あなたは、各セルの曲率値を返す関数を持っています。
  • X、Y、Zのメッシュがそのサーフェスをプロットします。

    Value mapped to color in matpotlib surface

    import numpy 
    import matplotlib.pyplot as plt 
    from mpl_toolkits.mplot3d import Axes3D 
    from matplotlib import cm 
    import math 
    
    # Here would be the surface_curvature function 
    
    X = numpy.array(ptX) 
    Y = numpy.array(ptY) 
    Z = numpy.array(ptZ) 
    p,g,m= surface_curvature(X,Y,Z) 
    C = m.reshape(X.shape) 
    
    C -= C.min() 
    C /= C.max() 
    
    fig = plt.figure() 
    ax = fig.add_subplot(111,projection = '3d') 
    n = numpy.reshape(m,numpy.shape(X)) 
    ax.plot_surface(X,Y,Z,facecolors = cm.jet(C), rstride=1, cstride=1) 
    plt.show() 
    

    、私はこれを取得:

あなたのコードを使用して、そしてあなたにm変数を想定したが、曲率(再び、これはあなたのコードである)私がしなければ、これです

matplotlibサーフェスの色にマップされた値です。あなたが作成したCが実際の曲率でない場合は、それを置き換える必要があります。

関連する問題