python で球面グリッド三角化

# coding: utf-8

# In[333]:

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import pyshtools as sh
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial import Delaunay
import math
import collections # 多値返却のため
get_ipython().magic(u'matplotlib notebook')


# In[334]:

# 経線、緯線数を与えて三角グリッドを作る
def spheregrid(nlat=10,nlon=10):
    phi = 2*math.pi/nlon * np.array([range(nlon)])
    theta = math.pi/(nlat-1) * np.array([range(nlat)])
    
    n = 2 + nlon * (nlat-2)
    
    ret = np.zeros((2,n))
    ret[0,n-1] = math.pi
    for i in range(nlat -2):
        for j in range(nlon):
            tmp = i * nlon + j + 1
            ret[0,tmp] = theta[0,i+1]
            ret[1,tmp] = phi[0,j]
    n_tri = nlon * 2 + nlon*2*(nlat-3)
    tri = np.zeros((n_tri,3),dtype=np.int)
    cnt = 0
    for i in range(nlon-1):
        tri[cnt,0] = 0
        tri[cnt,1] = i+1
        tri[cnt,2] = i+2
        cnt +=1
    tri[cnt,0] = 0
    tri[cnt,1] = nlon
    tri[cnt,2] = 1
    cnt +=1
    for i in range(nlat-3):
        for j in range(nlon-1):
            tri[cnt,0] = i*nlon + j + 1
            tri[cnt,1] = i*nlon + j + 2
            tri[cnt,2] = (i+1)*nlon + j + 2
            cnt += 1
            tri[cnt,0] = i*nlon + j + 1
            tri[cnt,1] = (i+1)*nlon + j + 1
            tri[cnt,2] = (i+1)*nlon + j + 2
            cnt += 1
        tri[cnt,0] = i*nlon + nlon
        tri[cnt,1] = i*nlon + 1
        tri[cnt,2] = (i+1)*nlon + 1
        cnt += 1
        tri[cnt,0] = i*nlon + nlon
        tri[cnt,1] = (i+1)*nlon + nlon
        tri[cnt,2] = (i+1)*nlon + 1
        cnt += 1
    for i in range(nlon-1):
        tri[cnt,0] = (nlat-3)*nlon + i + 1
        tri[cnt,1] = (nlat-3)*nlon + i + 2
        tri[cnt,2] = n-1
        cnt +=1
    tri[cnt,0] = (nlat-3)*nlon + nlon-1 + 1
    tri[cnt,1] = (nlat-3)*nlon + 1
    tri[cnt,2] = n-1
    
    X = np.ndarray((3,n))
    for i in range(ret[0,].size):
        X[0,i], X[1,i], X[2,i] = np.sin(ret[0,i])*np.cos(ret[1,i]), np.sin(ret[0,i])*np.sin(ret[1,i]), np.cos(ret[0,i])
    
    result = collections.namedtuple('result', 'X, angles, tri')
    return result(X=X, angles=ret, tri=tri)





# In[335]:

# 3d 正規乱点を単位球面に発生させる
def runitsphere(n,d=3):
    X = np.random.randn(n*d).reshape(n,d)
    r = np.sum(np.power(X,2), axis=1)
    X = X.T/np.power(r,0.5)
    return X


# In[336]:


# 球面上乱点を発生しその3角化をする関数を作る
def spheretri_random(n,s=0.001):
    X = runitsphere(n-1) # 点 (0,0,1)を加えて全部でn個の点を球面に置くためにn-1とする
    # (0,0,1)のごく近辺に3点を発生させる
    X_min = np.min(np.abs(X))
    xy_min = X_min * s
    three1x = xy_min
    three1y = 0
    three1z = math.sqrt(1-three1x**2-three1y**2)
    three2x = xy_min * (-1/2)
    three2y = xy_min * (math.sqrt(3)/2)
    three2z = math.sqrt(1-three2x**2-three2y**2)
    three3x = xy_min * (-1/2)
    three3y = xy_min * (-math.sqrt(3)/2)
    three3z = math.sqrt(1-three3x**2-three3y**2)
# 2通りの球面上点座標を作る
# (0,0,1)に近い3点をそれぞれ登録する場合
    x2 = np.append(X[0,],np.array([three1x,three2x,three3x]))
    y2 = np.append(X[1,],np.array([three1y,three2y,three3y]))
    z2 = np.append(X[2,],np.array([three1z,three2z,three3z]))
    # 3点を登録する代わりに(0,0,1)を登録する場合
    x = np.append(X[0,],np.array([0]))
    y = np.append(X[1,],np.array([0]))
    z = np.append(X[2,],np.array([1]))
# (0,0,1)近傍3点を区別する場合に関して
# (0,0,1) が無限遠となるように、球面を平面に開く
# そうすることでその二次元平面にて三角化する
# その上で(0,0,1)相当の3点は、最終的に、三角を作るものとする
    X2 = 2*x2/(2-(z2+1))
    Y2 = 2*y2/(2-(z2+1))
# ドロネー3角化(二次元平面にて行う)
# Triangulate parameter space to determine the triangles
    tri = Delaunay(np.array([X2,Y2]).T)
# 最後に追加した3点は、(0,0,1)のことなので、3角化の頂点座標のうち、追加点のノードIDを同じ値にする
    tri.simplices[tri.simplices > n-2]=n-1
# 返り値は、n個の点の座標とその三角形頂点リスト
    result = collections.namedtuple('result', 'X, tri')
    retX = np.vstack([np.vstack([x,y]),z])
    return result(X=retX, tri=tri.simplices)

# 三角形の面積を計算する
def triarea(v1,v2,v3):
    e12 = v2-v1
    e13 = v3-v1
    ip = np.dot(e12,e13)
    r1 = np.sqrt(np.sum(np.power(e12,2)))
    r2 = np.sqrt(np.sum(np.power(e13,2))) 
    costheta = ip/(r1*r2)
    sintheta = np.sqrt(1-np.power(costheta,2))
    area = 0.5 * r1*r2*sintheta
    return area
    
# 三角形の面積を三角化オブジェクトから作る
def triarea_ts(ts):
    ret = np.zeros(ts.tri.T[0,].size)
    for i in range(ts.tri.T[0,].size):
        id1 = ts.tri[i,0]
        id2 = ts.tri[i,1]
        id3 = ts.tri[i,2]
        if (id1-id2)*(id2-id3)*(id3-id1) == 0:
            ret[i] = 0
        else:
            v1 = np.array([ts.X[0,id1],ts.X[1,id1],ts.X[2,id1]])
            v2 = np.array([ts.X[0,id2],ts.X[1,id2],ts.X[2,id2]])
            v3 = np.array([ts.X[0,id3],ts.X[1,id3],ts.X[2,id3]])
            ret[i] = triarea(v1,v2,v3)
        
    return ret

# 三角形の重心座標を三角化オブジェクトから作る
def triweightedcenter(ts):
    ret = np.zeros((3,ts.tri.T[0,].size))
    for i in range(ts.tri.T[0,].size):
        id1 = ts.tri[i,0]
        id2 = ts.tri[i,1]
        id3 = ts.tri[i,2]

        v1 = np.array([ts.X[0,id1],ts.X[1,id1],ts.X[2,id1]])
        v2 = np.array([ts.X[0,id2],ts.X[1,id2],ts.X[2,id2]])
        v3 = np.array([ts.X[0,id3],ts.X[1,id3],ts.X[2,id3]])
        ret[0,i], ret[1,i], ret[2,i] = (v1 + v2 + v3)/3
        
    return ret

# 三角化プロット
def spheretriplot(ts):  
# 3Dplot
    fig = plt.figure()
    ax = Axes3D(fig)
# The triangles in parameter space determine which x, y, z points are
# connected by an edge
    ax.plot_trisurf(ts.X[0,],ts.X[1,],ts.X[2,], triangles=ts.tri, cmap=plt.cm.Spectral)
    #ax.plot_trisurf(ts.X[0,],ts.X[1,],ts.X[2,], triangles=ts.tri, color=col)
#ax.plot_trisurf(x, y, z, triangles=tri.simplices)
    plt.show()


# In[337]:

# 3Dデカルト座標を極座標に変える(r==0は考えない)
def decart2polar3d(x,y,z):
    r = math.sqrt(x**2+y**2+z**2)
    rxy = x**2+y**2
    if rxy==0:
        phi = 0
        theta = 0
        if z < 0:
            theta = math.pi
    else:
        if x==0:
            phi = math.pi/2
            if y < 0:
                phi = -math.pi/2
        else:
            phi = np.arctan(y/x)
        theta = np.arccos(z/r)
    return [phi,theta,r]
    
decart2polar3d(0,0,-1)

def decart2polar3d_ndarray(X):
    ret = np.zeros_like(X)
    for i in range(X[0,].size):
        tmp= decart2polar3d(X[0,i],X[1,i],X[2,i])
        ret[0,i],ret[1,i],ret[2,i] = tmp
    return ret


# In[338]:

# 三角の頂点・重心座標の場の値を球面調和関数係数numpy.ndarrayから作る
def SHval(Xpolar,coef):
    valreal = np.zeros(Xpolar[0,].size)
    valcomplex = np.zeros(Xpolar[0,].size,'complex')
    for i in range(Xpolar[0,].size):
    #valreal += [ coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])]
    #valcomplex += [coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])]
        valreal[i] = coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])
        valcomplex[i] = coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])
        
    result = collections.namedtuple('result', 'real, clx')
    return result(real = valreal,clx=valcomplex)


# In[339]:

# 球面になめらかな場を与える
# 球面調和関数の係数のセットとして与える
# x.coeffsは[0,,]にm>=0[1,,]にm<0が格納される
# [*,0,k]はabs(m)=kの値を格納する
def randSHcoef(lmax):
    num_coef = (lmax+1)**2
    value_re = np.random.randn(num_coef)
#print value_re
    x = sh.SHCoeffs.from_zeros(lmax)

    for l in range(0, lmax+1):
        for m in range(-l, l+1):
          x.set_coeffs(value_re[l*(l+1)+m], l, m)
    return x


# In[340]:

# scipy.special.sph_harmを使ってSHcoeffから、球面座標の値を計算
def coeff2valreal(cf,phi,theta):
    ret = 0
    #lmax = cf
    for l in range(0, cf[0,0,].size):
        for m in range(-l, l+1):
            if m >= 0:
                tmp = cf[0,l,m]
                ret += sph_harm(m,l,phi,theta).real
            else:
                tmp = cf[1,l,np.abs(m)]
                ret += sph_harm(m,l,phi,theta).imag
    return ret

def coeff2valcomplex(cf,phi,theta):
    ret = 0+1j
    #lmax = cf
    for l in range(0, cf[0,0,].size):
        for m in range(-l, l+1):
            if m >= 0:
                tmp = cf[0,l,m]
            else:
                tmp = cf[1,l,np.abs(m)]
            ret += sph_harm(m,l,phi,theta)
    return ret


# # 球面上一様乱点とその三角化オブジェクト tsを作る
# n = 200
# ts2 = spheretri_random(n)
# #print ts2
# ts = spheregrid(nlat=10,nlon=10)
# #print ts
# #print(ts.tri)
# #print ts.tri.T[0,].size
# # 三角形の面積
# areas = triarea_ts(ts)
# # 三角形の重心座標
# weightedcenter = triweightedcenter(ts)
# #print weightedcenter
# # 頂点と三角形重心の極座標
# Xpolar = decart2polar3d_ndarray(ts.X)
# #Xpolar = ts.angles
# WCpolar = decart2polar3d_ndarray(weightedcenter)
# #print weightedcenter.size
# #print WCpolar.size
# spheretriplot(ts)

# In[343]:

# 球面場モデルを作る
f = randSHcoef(5)


# In[344]:

# 場の値を計算
fvalV = SHval(Xpolar,f.coeffs)
fvalF = SHval(WCpolar,f.coeffs)
#print WCpolar[0,].size
#print WCpolar
#print fvalV
#print fvalF.clx.size


# In[ ]:




# In[345]:

x = randSHcoef(3)
print x.coeffs
# 係数はそのまま回転できる
x_rot = x.rotate(1,2,3, degrees=False)
print x_rot.coeffs

xexp = x.expand()


# In[346]:

coef = randSHcoef(3).coeffs
#coef[0,0,].size
print coeff2valreal(coef,Xpolar[0,0],Xpolar[1,0])
print coeff2valcomplex(coef,Xpolar[0,0],Xpolar[1,0])
#Xpolar[0,0]

valreal = np.zeros(Xpolar[0,].size)
valcomplex = np.zeros(Xpolar[0,].size,'complex')
for i in range(Xpolar[0,].size):
    #valreal += [ coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])]
    #valcomplex += [coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])]
    valreal[i] = coeff2valreal(coef,Xpolar[0,i],Xpolar[1,i])
    valcomplex[i] = coeff2valcomplex(coef,Xpolar[0,i],Xpolar[1,i])



# In[348]:

plt.figure()
plt.plot(valreal,valcomplex.real)
plt.show()
plt.figure()
plt.plot(valcomplex.real,valcomplex.imag)
plt.show()


# In[ ]: