# Copyright (C) Frank Zaverl, Jr.
# See file LICENSE for license information.
from s3dlib import __version__
import math
import warnings
import copy
from functools import reduce
#from time import time
import numpy as np
from scipy import interpolate, spatial
from skimage import measure
from matplotlib import cm, colors, image, tri
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
#.. simple warning string
warnings.formatwarning = lambda m,c,n,l,line=0 : str(m)+'\n'
_MAXREZ = 8 # in subclass surfaces, maximum recursive triangulation.
_MAXPLREZ = 10 # in ParametricLine lines, maximum recursive segmentation.
_MINRAD = 0.01 # in subclassed polar and spherical surface split basetypes.
_SPLITSIZE = 0.0001 # for split grid geometries, phi offset from 2*pi
_COOR_KWARGS = { # used for clipping & contour lines and line-filled surfaces
'coor': [ [0], [0, 'planar','xyz','p','P'],
[1, 'c', 'C', 'cylinder', 'cylindrical'],
[2, 's', 'S', 'sphere', 'spherical'],
[3, 'x', 'X'],
[4, 'y', 'Y'],
[5, 'z', 'Z'] ]
}
_COORSYS = { "XYZ": 0, "POLAR": 1, "CYLINDRICAL": 2, "SPHERICAL": 3, "LONGLAT": 4 }
_ILLUM = (1,0,1) # default light direction, (cmap_normals, shade, hilite & line direction).
_DFT_VLIM = [ -1.0, 1.0 ] # default face color value bounds
_DFT_VERTVLIM = [ 0.0, 1.0 ] # default vector magnitude value bounds
_DFT_ALR = 0.25 # default axis length ratio, head size to vector magnitude.
_DFT_VIEW = [ 30, -60 ] # default, axes elev and azim
# +----------------------------------------------------------+
# | S#Dlib v_1.2.0 was developed using Matplotlib v_3.8.0 |
# +----------------------------------------------------------+
# FutDev: Future Development notes.
# MROI: Minimum Return On coding effort Investment.
# (the juice isn't worth the squeeze for current release)
[docs]class Surface3DCollection(Poly3DCollection):
"""
Base class for 3D surfaces.
"""
_geomType = _COORSYS['XYZ']
[docs] @staticmethod
def chull(points, **kargs) :
"""
Surface3DCollection Convex Hull for a set of 3D points.
Parameters
----------
points : N x 3 float array-like
Returns
-------
Surface3DCollection object
"""
points = np.array(points)
hull = spatial.ConvexHull(points)
verts = points[hull.vertices]
cPnt = verts.mean(axis=0)
# set the faceIndices referencing the verts, not the points....
revlkup = { val:i for i,val in enumerate(hull.vertices) }
fIndices = [None]*len( hull.simplices )
for i,indx in enumerate(hull.simplices):
a,b,c = indx
fIndices[i] = [ revlkup[a] , revlkup[b] , revlkup[c] ]
# set the face indices using the RHR for outward normals from the center.
temp = [None]*len(fIndices)
for i,findx in enumerate(fIndices) :
a,b,c = findx
fcenter = (verts[a] + verts[b] + verts[c])/3
nordir = np.cross(verts[b]-verts[a],verts[c]-verts[b])
dotprod = np.dot(fcenter-cPnt,nordir)
temp[i] = [a,b,c] if dotprod>0 else [c,b,a]
fIndices = temp
surface = Surface3DCollection(verts,fIndices, **kargs)
return surface
[docs] @staticmethod
def implsurf(operation, drez=2.0, domain=1, name=None, **kwargs) :
"""
Surface3DCollection defined by an implicit function f(x,y,z) = 0.
Parameters
----------
operation : function object
Implicit function that takes one
argument, a 3xN Numpy array of x,y,z coordinates.
The function returns a single scalar value.
The value of zero defines the implicit surface.
drez : Float, optional, default: 2.0
Multiplier for the number of subdivisions within the domain.
Frez values range from 1 to 10. Function evaluation will be
the order of frez cubed.
domain : a number, list or array, default: 1
The domain of the function evaluation. For a number, n,
the x,y,z axes domains will be [-n,n]. For a 1-dimensional 2-element list,
[a,b] will be assigned the domain for all 3 axes. Using a list of list (array),
as [ [a,b],[c,d],[e,f] ], the domain is assigned individually for each of the
three coordinate axes.
Returns
-------
Surface3DCollection object
"""
#.....................................
def get_VF(op,N,dm) :
XD,YD,ZD = dm[0],dm[1],dm[2]
dmT = dm.T
A = dmT[1]-dmT[0]
B=np.array( [ XD[0],YD[0],ZD[0] ])
sp,level = 1/N, 0.001 # FutDev: level should not be a constant.
Nj = (N+1)*1j
xyz = np.mgrid[XD[0]:XD[1]:Nj,YD[0]:YD[1]:Nj, ZD[0]:ZD[1]:Nj ]
verts, fvIndices, _, _ = measure.marching_cubes(op(xyz), level, spacing=(sp, sp, sp))
C = np.empty([len(verts),3])
C[:] = A
verts = C*verts + B
return verts,fvIndices
#.....................................
dftDm = Surface3DCollection._interpret_domain(domain)
if not isinstance(drez, (int,float)) :
raise ValueError("Incorrect drez argument datatype passed to impf:", drez)
if drez<1 or drez>10 :
raise ValueError("Error: incorrect value drez passed to impf:", drez)
N = int( 10*drez )
verts, faces = get_VF(operation,N,dftDm)
surface = Surface3DCollection( verts, faces, **kwargs)
if name is None: name = _getFunctionName(operation,name)
if name is not None: surface._geomName = name
return surface
[docs] @staticmethod
def triangulateBase(rez,baseVcoor,baseFaceVertexIndices, midVectFunc) :
"""
Recursively subdivide triangles.
Each recursion subdivides each triangle by four.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
baseVcoor : V x 3 float array
An array of 'v' number of xyz vertex coordinates.
baseFaceVertextIndices : F x 3 int array
An array of 'F' number of face vertex indices.
midVectFunc : function object
A function that takes two xyz coordinate (list of 3)
representing a surface face edge. Returns one xyz coordinate
at the bisection of the edge, mapped on the surface.
Returns
-------
indexObj : a dictionary of vertex indices, for a surface
of F faces, E edges and V vertices.
'face' : F x 3 int array
'edge' : E x 2 int array
vertCoor : V x 3 float array
An array of V number of xyz vertex coordinates.
"""
vCoor = baseVcoor.copy()
fvIndices = []
evIndices = []
#.....................................
def getVerticesLine(degree, leftIndex, rightIndex) :
def recurs_edgeCoor(degree, indexOrigin, isRight, leftIndex, rightIndex, Eindex) :
m = degree - 1
delta = int(np.power(2,m))
if isRight : delta = -delta
i = indexOrigin - delta
mid_coor = midVectFunc( vCoor[leftIndex], vCoor[rightIndex] )
currentCoorIndex = len(vCoor)
vCoor.append(mid_coor)
Eindex[i] = currentCoorIndex
if m <= 0 :
evIndices.append( [leftIndex, currentCoorIndex] )
evIndices.append( [currentCoorIndex, rightIndex] )
return
recurs_edgeCoor(m,i,False,leftIndex,currentCoorIndex, Eindex)
recurs_edgeCoor(m,i,True,currentCoorIndex,rightIndex, Eindex)
return
#.....................................
iOrigin = int(np.power(2,degree))
Eindex = [None]*(iOrigin+1)
Eindex[0] = leftIndex
Eindex[iOrigin] = rightIndex
if degree==0 :
evIndices.append( [leftIndex, rightIndex] )
return Eindex
recurs_edgeCoor(degree, iOrigin, False, leftIndex, rightIndex, Eindex)
return Eindex
#.....................................
def recurs_faceIndices(A,B,C,rez,atCenter=True):
if rez==0 :
abc = [ A[0], B[0], C[0] ]
fvIndices.append( abc )
if atCenter : evIndices.append( [ B[0], C[0] ] )
return
# -----------------------------------------------------
J = int( (len(A)-1)/2 )
K = J + 1
# construct the center sub triangle edge vertices....
if rez==1 :
xz = [ A[1], C[1] ]
yx = [ B[1], A[1] ]
zy = [ C[1], B[1] ]
else:
xz = getVerticesLine(rez-1, A[J], C[J] )
yx = getVerticesLine(rez-1, B[J], A[J] )
zy = getVerticesLine(rez-1, C[J], B[J] )
# construct the 4 sub triangles...
recurs_faceIndices(A[:K],xz,C[J:],rez-1)
recurs_faceIndices(B[:K],yx,A[J:],rez-1)
recurs_faceIndices(C[:K],zy,B[J:],rez-1)
recurs_faceIndices(yx[::-1],zy[::-1],xz[::-1],rez-1,False)
return
#.....................................
# generate the base edge indices for each face
baseFaceEdgeIndices = []
edgeArrayIndexMatrix = [[None for x in range(len(vCoor))] for y in range(len(vCoor))]
for face in baseFaceVertexIndices :
A = [ face[0],face[1]]
B = [ face[1],face[2]]
C = [ face[2],face[0]]
edgeIndices = [A,B,C]
baseFaceEdgeIndices.append(edgeIndices)
for vertexIndex in edgeIndices :
i = vertexIndex[0]
j = vertexIndex[1]
if edgeArrayIndexMatrix[i][j] is None:
edgeArrayIndexMatrix[i][j] = getVerticesLine(rez, i, j )
edgeArrayIndexMatrix[j][i] = edgeArrayIndexMatrix[i][j].copy()[::-1]
# trianglulate the base triangles
for face in baseFaceEdgeIndices :
A = edgeArrayIndexMatrix[face[0][0]][face[0][1]]
B = edgeArrayIndexMatrix[face[1][0]][face[1][1]]
C = edgeArrayIndexMatrix[face[2][0]][face[2][1]]
x = (rez > 0)
recurs_faceIndices(A,B,C,rez, x)
vertexCoor = np.array(vCoor)
indexObj = { 'face': np.array(fvIndices) , 'edge': np.array(evIndices) }
return indexObj ,vertexCoor
[docs] @staticmethod
def rectangulateBase(rez,baseVcoor,baseFaceVertexIndices, midVectFunc) :
"""
Recursively subdivide quadrilaterals.
Each recursion subdivides each quadrilateral by four.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of the rectangle base faces.
Rez values range from 0 to 7.
baseVcoor : V x 3 float array
An array of 'v' number of xyz vertex coordinates.
baseFaceVertextIndices : F x 4 int array
An array of 'F' number of face vertex indices.
midVectFunc : function object
A function that takes two xyz coordinate (list of 4)
representing a surface face edge. Returns one xyz coordinate
at the bisection of the edge, mapped on the surface.
Returns
-------
indexObj : a dictionary of vertex indices, for a surface
of F faces, E edges and V vertices.
'face' : F x 4 int array
'edge' : E x 2 int array
vertCoor : V x 3 float array
An array of V number of xyz vertex coordinates.
"""
vCoor = baseVcoor.copy()
fvIndices = []
evIndices = []
#.....................................
def getVerticesLine(degree, leftIndex, rightIndex) :
def recurs_edgeCoor(degree, indexOrigin, isRight, leftIndex, rightIndex, Eindex) :
m = degree - 1
delta = int(np.power(2,m))
if isRight : delta = -delta
i = indexOrigin - delta
mid_coor = midVectFunc( vCoor[leftIndex], vCoor[rightIndex] )
currentCoorIndex = len(vCoor)
vCoor.append(mid_coor)
Eindex[i] = currentCoorIndex
if m <= 0 :
evIndices.append( [leftIndex, currentCoorIndex] )
evIndices.append( [currentCoorIndex, rightIndex] )
return
recurs_edgeCoor(m,i,False,leftIndex,currentCoorIndex, Eindex)
recurs_edgeCoor(m,i,True,currentCoorIndex,rightIndex, Eindex)
return
#.....................................
iOrigin = int(np.power(2,degree))
Eindex = [None]*(iOrigin+1)
Eindex[0] = leftIndex
Eindex[iOrigin] = rightIndex
if degree==0 :
evIndices.append( [leftIndex, rightIndex] )
return Eindex
recurs_edgeCoor(degree, iOrigin, False, leftIndex, rightIndex, Eindex)
return Eindex
#.....................................
def recurs_faceIndices(A,B,C,D,rez):
# similar to inner function of triangulateBase
if rez==0 :
abc = [ A[0], B[0], C[0], D[0] ]
fvIndices.append( abc )
return
# -----------------------------------------------------
J = int( (len(A)-1)/2 )
K = J + 1
# midpoint indices ...............
mid_coor = midVectFunc( vCoor[ D[J] ], vCoor[ B[J] ] )
i_E = len(vCoor)
vCoor.append(mid_coor)
# arrays of the four inner quadrilaterals ..............
pE = getVerticesLine(rez-1, A[J], i_E )
qE = getVerticesLine(rez-1, B[J], i_E )
rE = getVerticesLine(rez-1, C[J], i_E )
sE = getVerticesLine(rez-1, D[J], i_E )
# sudivide the four inner quadraterals .................
recurs_faceIndices( A[:K], pE, sE[::-1], D[J:], rez-1)
recurs_faceIndices( B[:K], qE, pE[::-1], A[J:], rez-1)
recurs_faceIndices( C[:K], rE, qE[::-1], B[J:], rez-1)
recurs_faceIndices( D[:K], sE, rE[::-1], C[J:], rez-1)
return
#.....................................
# generate the base edge indices for each face
baseFaceEdgeIndices = []
edgeArrayIndexMatrix = [[None for x in range(len(vCoor))] for y in range(len(vCoor))]
for face in baseFaceVertexIndices :
A = [ face[0],face[1]]
B = [ face[1],face[2]]
C = [ face[2],face[3]]
D = [ face[3],face[0]]
edgeIndices = [A,B,C,D]
baseFaceEdgeIndices.append(edgeIndices)
for vertexIndex in edgeIndices :
i = vertexIndex[0]
j = vertexIndex[1]
if edgeArrayIndexMatrix[i][j] is None:
edgeArrayIndexMatrix[i][j] = getVerticesLine(rez, i, j )
edgeArrayIndexMatrix[j][i] = edgeArrayIndexMatrix[i][j].copy()[::-1]
# rectangulate the base quadrilaterals
for face in baseFaceEdgeIndices :
A = edgeArrayIndexMatrix[face[0][0]][face[0][1]]
B = edgeArrayIndexMatrix[face[1][0]][face[1][1]]
C = edgeArrayIndexMatrix[face[2][0]][face[2][1]]
D = edgeArrayIndexMatrix[face[3][0]][face[3][1]]
recurs_faceIndices(A,B,C,D,rez)
vertexCoor = np.array(vCoor)
indexObj = { 'face': np.array(fvIndices) , 'edge': np.array(evIndices) }
return indexObj ,vertexCoor
[docs] @staticmethod
def coor_convert(xyz, tocart=False) :
"""Coordinate transformation.
To be overriddden by any subclass not in Cartesian coordinates.
"""
xyz = np.array(xyz)
return xyz
@staticmethod
def _interpret_domain(domain) :
"""
Process domain argument for implsurf and CubicSurface.domain to
the domain for x,y,z domains.
"""
errCode = 0
if isinstance(domain, (int,float)) :
dftDm = np.array( [ [-domain,domain ], [-domain,domain ], [-domain,domain ] ] )
else :
correctDomain = True
try :
npDm = np.array(domain, dtype=float)
shDm = npDm.shape
nbDm = npDm.ndim
correctDomain = nbDm==1 or nbDm==2
errCode = 1
except :
correctDomain = False
errCode = 2
if correctDomain :
if nbDm == 1 :
dftDm = np.array( [ domain, domain, domain] )
correctDomain = shDm[0] == 2
errCode=3
elif nbDm == 2 :
dftDm = np.array( domain )
correctDomain = shDm[0]==3 and shDm[1]==2
errCode=4
else :
correctDomain = False
errCode=5
if not correctDomain :
raise ValueError("Incorrect domain argument passed to impf.", errCode)
return dftDm
@staticmethod
def _map3Dto2Dsquare(xyz) :
""" Map surface to rectangular unit coordinates (0,1). """
x,y,_ = xyz
a = (x+1)/2
b = (y+1)/2
return [a,b]
@staticmethod
def _sfev(rez,fevArray) :
""" Calculate number of surface faces, edges, vertices. """
# vestigial
Fo, phi, Vx = fevArray
eta = np.power(2,rez)
F = eta*eta*Fo
E = np.int( 3*F/2 + eta*phi )
V = np.int( F/2 + eta*phi + Vx )
return [F,E,V]
@classmethod
def _grid_face_indices(cls,U, V, isSplit=False, isTria=False) :
""" method only called by the classmethod _square_net or
the CylindricalSurface classmethod _cyl_vol_square_net
"""
if isSplit :
t = np.linspace(0,V*(U+1)-1,V*(U+1),dtype=int )
ts = t.reshape( (V,U+1) )
a = ts[:,0:U].flatten()
findx = np.array( [ a, a+1, a+U+2, a+U+1 ] )
else :
a = np.linspace(0,V*U-1,V*U,dtype=int )
bs = a.reshape( (V,U) ).T
order = np.append(np.linspace(1,U-1,U-1,dtype=int),0)
b = bs[order].T.flatten()
findx = np.array( [ a, b, b+U, a+U ] )
if isTria :
a,b,c,d = findx
findx = np.concatenate( ( [a,b,c],[c,d,a] ), axis=1 )
return findx.T
@classmethod
def _square_net(cls, numV, numU, basetype,minRad, splitSize, name, **kwargs) :
""" method only called by inherited class classmethods 'grid'. """
divLimits = [1,3]
if cls._geomType is _COORSYS['XYZ'] : divLimits = [1,1]
if cls._geomType is _COORSYS['SPHERICAL'] : divLimits = [2,3]
if numV < divLimits[0] :
raise ValueError('Grid division {} must be greater or equal to {}'.format(str(numV),divLimits[0] ))
if numU < divLimits[1] :
raise ValueError('Grid division {} must be greater or equal to {}'.format(str(numU),divLimits[1] ))
if basetype == None : basetype = 'd'
if minRad is None : minRad = _MINRAD
if splitSize is None : splitSize = _SPLITSIZE
basetype = basetype.lower()
baselist = ['d','s','q','w','r','x']
if basetype not in baselist :
raise ValueError("Grid basetype '{}' is not recognized. Possible values are: {}.".format(basetype,list(baselist)))
isSplit = (basetype == 's') or (basetype == 'w') or (basetype == 'x')
needsTriag = not ( (basetype == 'r') or (basetype == 'x') )
hasCenter = (basetype == 'd') or (basetype == 's')
# hasCenter does not apply for cylindrical grids of this type.
if cls._geomType is _COORSYS['CYLINDRICAL'] : hasCenter=False
# Note: numX and num_X are the number of divisions and the number of coor positions.
if cls._geomType is _COORSYS['XYZ'] :
isSplit = True
hasCenter = False
num_U,min_U, max_U = numU+1,-1.0,1.0
num_V = numV + 1
coor_U = np.linspace(min_U,max_U,num_U)
U_coor = np.tile(coor_U,num_V)
else :
# T angle ranges, used for disk,cylinder and sphere.....
num_U = numU+1 if isSplit else numU
if hasCenter :
num_V = numV if cls._geomType is _COORSYS['POLAR'] else numV-1
else :
num_V = numV+1
min_U = 0.0
max_U = (1-1/num_U)*2*np.pi
if isSplit : max_U = (2-splitSize)*np.pi
coor_U = np.linspace(min_U,max_U,num_U)
U_coor = np.tile(coor_U,num_V)
totVerts = num_U*num_V
if cls._geomType is _COORSYS['XYZ'] :
val_W = 0.0 # z is constant at 0.
coor_V = np.linspace(-1.0,1.0,num_V)
order_abc = lambda u,v,w :[u,v,w]
if cls._geomType is _COORSYS['POLAR'] :
val_W = 0.0 # z is constant at 0.
end_V = 1/numV if hasCenter else minRad
coor_V = np.linspace(1.0,end_V,num_V)
order_abc = lambda u,v,w : [v,u,w]
if cls._geomType is _COORSYS['CYLINDRICAL'] :
val_W = 1.0 # r is constant at 1.
coor_V = np.linspace(-1.0,1.0,num_V)
order_abc = lambda u,v,w : [w,u,v]
if cls._geomType is _COORSYS['SPHERICAL'] :
val_W = 1.0 # r is constant at 1.
st_angle = np.pi/numV if hasCenter else np.arcsin(minRad)
stt_V = np.pi - st_angle
end_V = st_angle
coor_V = np.linspace(stt_V,end_V,num_V)
order_abc = lambda u,v,w : [w,u,v]
W_coor = np.tile(val_W,totVerts)
V_coor = np.tile(coor_V,(num_U,1))
V_coor = np.ravel(V_coor,order='F')
abc = order_abc(U_coor,V_coor,W_coor)
abc = np.array(abc)
xyz = cls.coor_convert(abc,True)
verts = np.array(xyz).T
numVt = num_V-1 if hasCenter else numV
faceIndices = cls._grid_face_indices(numU,numVt,isSplit,needsTriag)
if hasCenter :
# add top triangular faces for polar & spherical surfaces ...
topVert = np.array([[0.0,0.0,val_W]])
topIndex = totVerts
startIndex = topIndex - num_U
endIndex = topIndex-2 if isSplit else topIndex-1
a = np.linspace(startIndex,endIndex,numU,dtype=int)
b = a+1 if isSplit else np.append( a[1:], [startIndex] )
c = np.full(numU,topIndex)
abc = np.array([a,b,c]).T
verts = np.append(verts,topVert, axis=0)
faceIndices = np.append(faceIndices,abc, axis=0)
if cls._geomType is _COORSYS['SPHERICAL'] :
# add bottom triangular faces for spherical surfaces ...
btmVert = np.array([[0.0,0.0,-val_W]])
btmIndex = totVerts+1
startIndex,endIndex = 0,numU-1
a = np.linspace(startIndex,endIndex,numU,dtype=int)
b = a+1 if isSplit else np.append( a[1:], [startIndex] )
c = np.full(numU,btmIndex)
bac = np.array([b,a,c]).T # bac instead for abc for outward normals.
verts = np.append(verts,btmVert, axis=0)
faceIndices = np.append(faceIndices,bac, axis=0)
surface = cls(**kwargs) # set to the class default object
# reset verts, faceIndices and edges (retain class coordinate system).
surface.vertexCoor = verts
surface.fvIndices = faceIndices
surface.evIndices = surface._get_edges_from_faces(faceIndices)
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( verts[faceIndices] )
surface._set_geometric_bounds()
surface.name = name
surface._basetype = 'grid_'+basetype
return surface
@classmethod
def _fev(cls,rez,basetype,dercls) :
""" Sets up _sfev method, based on calls from the subclass dercls. """
# vestigial
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be an int, 0 <= rez <= {}'.format(rez,_MAXREZ))
if basetype is None: basetype = dercls._default_base
basesurf = dercls._base_surfaces
if basetype not in basesurf :
raise ValueError('Basetype {} is not recognized. Possible values are: {}'.format(basetype,list(basesurf)))
return cls._sfev(rez, basesurf[basetype]['fevParam'] )
@classmethod
def _uvw_to_xyz(cls,rst,abc) :
""" Rotational transformation at a coordinate.
To be overridden by any subclass not in Cartesian coordinates.
"""
return np.transpose(rst)
@classmethod
def _disp_to_xyz( cls, uvw, theta=None, phi=None) :
""" Rotational transformation from rotation angles. """
if theta is None : return uvw
# each set of uvw has a different rotation matrix
# that must be evaluated one at a time. MROI
if phi is None : phi = np.zeros( len(theta) )
uvwT = np.transpose(uvw)
dispVect = []
for i in range(len(theta)) :
rotMx = eulerRot(theta[i],phi[i], inrad=True)
vec = np.dot( uvwT[i], rotMx )
dispVect.append(vec)
return dispVect
def _postProc_surfaceColors(self,onlyFaces=False) :
""" Matplotlib will 'fillin' colors during rendering using the
facecolor array. (this is NOT the assigned color to faces).
However, S3Dlib will use the inherited _facecolors
for the assigned face colors.
Whenever color is 'externally' assigned or reassigned,
need to directly assign colors for each face for any
method that uses facecolors. This occurs initially
in the method before colors are assigned
(set_color, +, clipping, shading, edges, etc.)
Note that when initialization using **kargs,
onlyFaces should be True so don't reset edgecolors
if they are assigned.
"""
N = len(self.fvIndices) # number of faces
initFC = self._dfacecolors # update for v_1.2.0
inLen = len(initFC) # number of facecolors
if inLen == N : return
if isinstance(initFC,np.ndarray) : initFC = initFC.tolist()
n = math.ceil(N/inLen)
colorcoll = initFC*n
if onlyFaces : self.set_facecolor(colorcoll[0:N])
#if onlyFaces : self._dfacecolors = colorcoll[0:N] # update for v_1.2.0
else : self.set_color(colorcoll[0:N])
return
def _set_index_relations(self) :
""" from the vertex and faceIndices, set the indices:
evIndices: Array, for each edge, 2 vertex indicies.
vfIndicesList: List, for each vertex, face indices.
efIndicesList: List, for each edge, 0 to 2 face indices.
"""
edgeIndices = self._get_edges_from_faces(self.fvIndices)
self.evIndices = np.array(edgeIndices)
self.vfIndicesList = self._vertexCommonFaceIndices()
self.efIndicesList = self._edgeCommonFaceIndices()
return
def __init__(self, vertexCoor, faceIndices, name=None, vcolor=None, **kwargs) :
"""
3D surface of connected F faces, each face with N number of vertices.
The surface has V number of vertices and E number of edges.
Parameters
----------
vertexCoor : V x 3 float array-like
An array of V number of xyz vertex coordinates.
faceIndices : F x N int list or array
A list or array of F number of faces with N vertex indices.
For lists, N may be different for each face.
name : string, optional, default: None.
Descriptive identifier for the geometry.
vcolor : V x 3 float array-like, default: None.
An array of V number of vertex colors.
(not implemented, reserved for future development)
Other Parameters
----------------
**kwargs
All other arguments are passed on to mpl_toolkits.mplot3d.art3d.Poly3DCollection.
Valid argument names include: color, edgecolors, facecolors, linewidths, cmap.
"""
self.vertexCoor = np.array(vertexCoor,dtype=float)
self.baseVertexCoor = np.array(vertexCoor,dtype=float)
self.fvIndices = np.array(faceIndices, dtype=object) # update for v_1.2.0
if self.fvIndices.ndim == 1 :
# assume input is a list of faces with different number of edges.
self._initedges = self._get_edges_from_faces(self.fvIndices)
self._initverts = np.array(vertexCoor,dtype=float)
fvIndNew = self._triangulate_faceIndices(faceIndices)
fcolor = None
# note: if edgecolors are not defined, default to defined facecolor list,
# not the defined color ( continue with different 'magic' ) MROI.
# following 'may' be vestigial to be removed ?
if 'color' in kwargs : fcolor = kwargs['color']
if 'facecolor' in kwargs : fcolor = kwargs['facecolor']
if colors.is_color_like(fcolor) :
# single color assigned
fvIndNew = self._triangulate_faceIndices(faceIndices)
else :
if fcolor is None :
# no color assigned
fvIndNew = self._triangulate_faceIndices(faceIndices)
else : # list of colors assigned
fvIndNew, faceColorsNew = self._triangulate_faceIndices_faceColors(faceIndices,fcolor)
kwargs['facecolor'] = faceColorsNew
kwargs['color'] = faceColorsNew
self.fvIndices = np.array(fvIndNew)
else :
self.fvIndices = np.array(faceIndices) # update for v_1.2.0
self._set_index_relations()
f = self.fvIndices
v = self.vertexCoor
super().__init__(v[f], **kwargs)
# note: in following, only set the facecolors since Matplotlib auto-setting
# the edges, if not explicitely passed with edgecolor kwargs.
self._postProc_surfaceColors(True)
#self._geomName = None # Redundant from next line, internal tracking to determine if set...
self.name = name
self.valuesName = None
# ========================================================================
self.coorType = _COORSYS["XYZ"]
self.baseSurfaceColor = np.array(self._dfacecolors[0]).flatten().tolist() # update for v_1.2.0
self.vertexColor = np.array(self._dfacecolors[0]).copy()[np.newaxis,:] # update for v_1.2.0
# following assigned when surface color applied using an operation
self.vertexValues = None
self.faceValues = None
if len (self._dedgecolors) == 0 : # update for v_1.2.0
self.baseEdgeColor = self.baseSurfaceColor
else :
self.baseEdgeColor = np.array(self._dedgecolors[0]).flatten().tolist() # update for v_1.2.0
self.isSurfaceColored = False # flag set True for color array export.
self.onlyColoredFaces = False # flag set True for cmap normals.
self._surfaceShapeChanged = False # control flag for image operations
self._isCompositeSurface = False # control flag for a composite object.
self.disp_Vector = None # overridden in subclass.
self._normal_scale = 1.0 # overridden in subclass default
self.restoredClip = None # dictionary to restore from clip operation.
self.vector_field = None # for vector fields when transformed.
self._svd_dict = None # results from map_geom_from_svd
self.scale = [1,1,1] # for axis when transformed.
self.translation = [0,0,0] # for axis when transformed.
self.rotation = np.identity(3) # for axis when transformed.
self._rez = 0 # used by child classes.
self._bounds = {}
self._set_geometric_bounds()
self._bounds['vlim'] = _DFT_VLIM
self._bounds['vertvlim'] = _DFT_VERTVLIM
[docs] def __add__(self, othr) :
"""
Combine two surface objects into a single surface object.
Parameters
----------
othr : Surface3DCollection object
Returns
-------
Surface3DCollection object
"""
if not isinstance(othr, Surface3DCollection) :
raise ValueError('Add operations can only apply between Surface3DCollection objects.')
# MROI: should check if 2 surface faces are both triangular or quadrilateral
topVerCoor = self.vertexCoor
botVerCoor = othr.vertexCoor
totVerCoor = np.append(topVerCoor,botVerCoor,axis=0)
nextIndex = len(topVerCoor)
topFVindices = self.fvIndices
botFVindices = othr.fvIndices
FVtail = np.add(botFVindices,nextIndex)
totFVindices = np.append(topFVindices,FVtail,axis=0)
obj = Surface3DCollection(totVerCoor,totFVindices)
obj._isCompositeSurface = True
top_onearr = np.ones( len(topFVindices) )[:, np.newaxis]
top_orig_colors = self._dfacecolors # update for v_1.2.0
if len(top_orig_colors) == 1 : top_orig_colors = top_orig_colors*top_onearr
bot_onearr = np.ones( len(botFVindices) )[:, np.newaxis]
bot_orig_colors = othr._dfacecolors # update for v_1.2.0
if len(bot_orig_colors) == 1 : bot_orig_colors = bot_orig_colors*bot_onearr
total_colors = np.append(top_orig_colors,bot_orig_colors,axis=0)
obj.set_color(total_colors)
obj._set_geometric_bounds()
return obj
def __str__(self) :
# parent class doesn't have specific properties...
# Note: self._rez may be a string or int. probably not good but MROI
try: bs = ' ( {} , {} )'.format(self._rez,self._basetype)
except: bs = ' '
numFaces = len(self.fvIndices)
numVerts = len(self.vertexCoor)
name = self.__class__.__name__
sz = ': faces: {}, vertices: {}'.format(numFaces,numVerts)
val = name+bs+sz
return val
def _calc_rez(self,N0) :
# string val for equivalent grid rez
numFaces = len(self.fvIndices)
cRez = np.log2(numFaces/N0)/2
return '{:.1f}'.format(cRez)
def _vertexCommonFaceIndices(self) :
""" list of face indices at each vertex. """
# note: a vertex may have NO faces (empty faceIndex list)
# due to clipping or a vertex not used for a face, so:
vfIndicesList = [ [] for i in range(len(self.vertexCoor))]
# note: usually vertices may have 1 to 6 common faces
for i,face in enumerate(self.fvIndices) :
for vInx in face :
if vfIndicesList[vInx] is None : vfIndicesList[vInx] = [i]
else : vfIndicesList[vInx].append(i)
return vfIndicesList
def _edgeCommonFaceIndices(self) :
""" list of face indices adjacent to each edge. """
efIndicesList = [None]*len(self.evIndices)
# note: edges may have 1 or 2 common faces
for i,eIndex in enumerate(self.evIndices) :
head_vfi = self.vfIndicesList[ eIndex[0] ]
tail_vfi = self.vfIndicesList[ eIndex[1] ]
for fi in head_vfi :
if (fi in tail_vfi) :
if efIndicesList[i] is None: efIndicesList[i] = [fi]
else: efIndicesList[i].append(fi)
return efIndicesList
def _contourLineCollection( self, dotprod, projVect, *args ) :
"""
ColorLine3DCollection of segments on faces interecting a plane.
dotprod - dot product between each vertex and plane normal.
projVect - interection plane normal
*args - list of distances of the plane to the origin.
"""
# FutDev: There are 'vestigal' loops during early development
# debugging which need to be cleaned out. (MROI for release)
# ---------------------------------------------------------------
def getLineVals( dotprod,aLen,projVect) :
# ................................................................
def order_value( face,eIndices, conCoor, planeNormal) :
A = eIndices[0]
B = eIndices[1]
fvi = self.fvIndices[face]
verts = self.vertexCoor[ fvi ]
faceNormal = self._get_face_normals([verts])
edgeVert_A = conCoor[A]
edgeVert_B = conCoor[B]
A_to_B = np.subtract(edgeVert_B, edgeVert_A)
edgeCross = np.cross(A_to_B, planeNormal )
edgeCross= np.divide( edgeCross, np.linalg.norm(edgeCross) )
test = np.dot(faceNormal, edgeCross)
return [A,B] if test>0 else [B,A]
# ................................................................
# .. get intersect coor for each edge........
evi = self.evIndices
verts = self.vertexCoor
dp_a = aLen - dotprod
eElev = dp_a[evi]
# eElev is the distance of the two edge vertices above/below the plane.
# edges cross the plane if one vertex is above, the other below.
# hence the product of the two distances is negative if edge crosses the plane.
edChk = eElev[:,0]*eElev[:,1]
xCoor = [] # intersect coordinates.
evix = [] # the edge index of the intersect coordinate.
for i in range( len(evi) ) :
if edChk[i] < 0.0 :
p1 = verts[ evi[i,0] ]
p2 = verts[ evi[i,1] ]
dp1 = dotprod [ evi[i,0] ]
dp2 = dotprod [ evi[i,1] ]
a_mdp1 = dp_a [ evi[i,0] ]
lmb = a_mdp1/( dp2 - dp1 )
px = lmb*np.subtract(p2,p1) + p1
xCoor.append(px)
evix.append(i)
# .. determine faces having edges which have an intersect.....
if len(xCoor) == 0 : return None,None,None
efi = self.efIndicesList
numbXedges = len(evix)
findexDict = {} # faces(key=face index) which have two edges (value=xCoor index)
for i in range(numbXedges) :
ei = evix[i]
efiArr = efi[ei]
for fi in efiArr :
if fi in findexDict :
curVal = findexDict[fi]
if isinstance(curVal,list) :
# DevNote: this modifications is a consequence of quadrilateral
# faces, which may have more than one intersectioon to a plane.
# This only eliminates runtime error. Solution is to use a surface
# with lrez > 0, then generated contours.
newVal = [curVal[0],i]
else :
newVal = [curVal,i]
findexDict[fi] = newVal
else :
findexDict[fi] = i
# remove face interset points if face only has one intersect edge
# ( intersect may be on a vertex )
#
test = [ x for x in findexDict.values() if type(x) is list ]
if len(test) == 0 : return None,None,None
#
face, coorInx = [], []
for key,value in findexDict.items() :
if type(value) is list :
valueX = order_value(key,value,xCoor,projVect)
coorInx.append(valueX)
face.append(key)
if len(coorInx) == 0 : return None,None,None
return np.array(xCoor), coorInx, face
# ---------------------------------------------------------------
def getColors(colors, indices) :
con_colors = [None]*len(indices)
for i in range(len(indices)) :
con_colors[i] = colors[indices[i]]
return con_colors
# ---------------------------------------------------------------
# get the surface colors which will be applied to the default contour colors
self._postProc_surfaceColors()
surf_colors = self._dfacecolors # update for v_1.2.0
if type(args[0]) is tuple : args=args[0]
line = None
for aLen in args :
vertexCoor,segmIndices,faceIndices = getLineVals( dotprod,aLen,projVect )
if vertexCoor is None :
warnings.warn( 'Surface contourLines not found for dist value of {}'.format(aLen) )
continue
tline = ColorLine3DCollection(vertexCoor,segmIndices)
if tline is None : continue
tline.set_color(getColors(surf_colors,faceIndices)) # default surface color
if line is None: line = tline
else: line += tline
return line
def _transformVector(self, orgCoor, operation, returnxyz) :
""" Functional transformation of coordinates. """
xyz = np.transpose(orgCoor)
abc = self.coor_convert(xyz)
rst = np.array(operation(abc))
if returnxyz : XYZ = rst
else : XYZ = self.coor_convert( rst , tocart=True )
return np.transpose(XYZ)
def _get_vectorfield_Line3DCollection(self, vector_field, color, width, alr) :
"""Vector3DCollection at surface vertices."""
self.vector_field = vector_field
lcol = Vector3DCollection(self.vertexCoor,vector_field,alr,colors=color, linewidths=width)
# --- 1.1 mod --------------------
lcol.coorType = self.coorType
lcol.uvwOrientationType = self.coorType
return lcol
def _get_vertex_normals(self) :
"""Unit normals of at vertex coordinates. """
# The normalized unweighted average of the surface
# normals of the faces that contain that vertex.
verts = self.vertexCoor[ self.fvIndices ]
faceNormalCoor = self._get_face_normals(verts)
# note: following list are the face indices that each vertex is a member.
# Since the number of faces a vertex is attached ranges from 1 - 6,
# needs a list, not an np.array.
vfIndicesList = [ None ] *len(self.vertexCoor)
for i in range(len(vfIndicesList)) : vfIndicesList[i] = []
for fInx in range( len(self.fvIndices) ) :
face = self.fvIndices[fInx]
for vInx in face :
vfIndicesList[vInx].append(fInx)
# DevNote: probably a 'cleaner' method? MROI
normals = [ None ] *len(self.vertexCoor)
for i in range(len(normals)) :
vert = vfIndicesList[i]
vertList = faceNormalCoor[ vert ]
vertSum = np.sum(vertList,axis=0)/len(vertList)
unitVector = np.divide( vertSum, np.linalg.norm(vertSum) )
normals[i] = unitVector
normals = np.array(normals)
return normals
def _get_face_normals(self,faceCoor,ausf=None) :
""" Unit normals of triangular face coordinates."""
vfT = np.transpose(faceCoor, (1,2,0) )
vAB = np.subtract( vfT[1], vfT[0] )
vAC = np.subtract( vfT[2], vfT[0] )
vABt = np.transpose(vAB)
vACt = np.transpose(vAC)
if ausf is not None:
vABt = vABt*ausf
vACt = vACt*ausf
cross = np.cross(vABt,vACt)
sz = np.linalg.norm(cross,axis=1)[:,np.newaxis]
return np.divide(cross,sz)
def _get_face_centers(self,faceCoor) :
""" Face centers from face vertex coordinates. """
vfT = np.transpose(faceCoor, (1,2,0) )
numVerts = vfT.shape[0]
sumV = np.sum(vfT, axis=0)
return np.transpose(sumV)/numVerts
def _set_geometric_bounds(self) :
""" set the bounds dictionary from surface vertex coordinates. """
_set_geomBounds(self.vertexCoor,self._bounds)
return
def _viewportCoor(self,xyzCoor,viewport=None) :
""" UV coordinates for image texture mapping """
# better way to do this? MROI
xyz = np.transpose(xyzCoor)
ab = self._map3Dto2Dsquare( xyz )
a,b = ab
if viewport is None : return a,b,np.full(len(a),True)
As,Bs,Ae,Be = viewport
Aref = 1.0-As
Bref = 1.0-Bs
Adelta = np.mod(Ae+Aref,1.0)
Bdelta = np.mod(Be+Bref,1.0)
if Adelta <=0.0 : Adelta = 1.0
if Bdelta <=0.0 : Bdelta = 1.0
Ap = np.mod(a+Aref,1.0)
Bp = np.mod(b+Bref,1.0)
Av = np.divide(Ap,Adelta)
Bv = np.divide(Bp,Bdelta)
inViewport = Av<=1
inViewport = np.where(Bv>1,np.full(len(inViewport),False),inViewport)
Av = np.clip(Av,0.0,1.0)
Bv = np.clip(Bv,0.0,1.0)
return Av,Bv,inViewport
def _triangulate_faceIndices(self,fvIndx) :
fvIndices, _ = self._triangulate_faceIndices_faceColors(fvIndx)
return fvIndices
def _triangulate_faceIndices_faceColors(self,fvIndx,faceColors=None) :
"""
Subdivide all faces into triangles. Two uses with number face edges is 3,4,5 or 6:
1. initial triangulate collection of faces having a differnt number of edges
called from __init__ with faceIndices passed as a list.
2. triangulate collection of faces with identical number of edges
called from triangulate where faceIndicess is a numpy array
"""
# .............................................................
def divideFace4(face) :
dist_a = np.linalg.norm( self.vertexCoor[ face[0]] - self.vertexCoor[ face[2]] )
dist_b = np.linalg.norm( self.vertexCoor[ face[1]] - self.vertexCoor[ face[3]] )
if dist_a < dist_b :
face1 = [ face[0], face[1], face[2] ]
face2 = [ face[2], face[3], face[0] ]
else :
face1 = [ face[1], face[2], face[3] ]
face2 = [ face[3], face[0], face[1] ]
return face1, face2
def divideFace5(face) :
diag = lambda i : self.vertexCoor[ face[i]] - self.vertexCoor[ face[(i+2)%5] ]
fc3 = lambda i : [ face[i], face[(i+1)%5], face[(i+2)%5] ]
fc4 = lambda i : [ face[(i+2)%5], face[(i+3)%5], face[(i+4)%5], face[(i+5)%5] ]
dist = [None]*5
for i in range(5) : dist[i] = np.linalg.norm( diag(i) )
minDiag = 0
for i in range(1,5) :
if dist[i] < dist[minDiag] : minDiag = i
face1 = fc3(minDiag)
otherFace = fc4(minDiag)
face2, face3 = divideFace4(otherFace)
return face1, face2, face3
def divideFaces6(face) :
# .. not the best, but easiest (could calc midpoint). MROI
dist = [None]*3
dist[0] = np.linalg.norm( self.vertexCoor[ face[0]] - self.vertexCoor[ face[3]] )
dist[1] = np.linalg.norm( self.vertexCoor[ face[1]] - self.vertexCoor[ face[4]] )
dist[2] = np.linalg.norm( self.vertexCoor[ face[2]] - self.vertexCoor[ face[5]] )
iMin = 0
if dist[1] < dist[0] : iMin = 1
if dist[2] < dist[iMin] : iMin = 2
faceA = [ face[iMin], face[ iMin+1 ], face[ iMin+2 ], face[ iMin+3 ] ]
faceB = [ face[iMin+3], face[(iMin+4)%6], face[(iMin+5)%6], face[(iMin+6)%6] ]
face1, face2 = divideFace4(faceA)
face3, face4 = divideFace4(faceB)
return face1, face2, face3, face4
# .............................................................
# FutDev: easy to understand but need to cleanup, MROI
subFaceIndices = []
fColor = []
setColors = False
if faceColors is not None:
setColors = True
if len(fvIndx) != len(faceColors) :
# if faceColors, then this should not happen, but report if so....
warnings.warn('Internal ERROR: tfc001')
return None
for i, face in enumerate(fvIndx) :
numVerts = len(face)
if numVerts == 3 :
subFaceIndices.append(face)
if setColors : fColor += [faceColors[i]]*1
elif numVerts == 4 :
f1, f2 = divideFace4(face)
subFaceIndices.append(f1)
subFaceIndices.append(f2)
if setColors : fColor += [faceColors[i]]*2
elif numVerts == 5 :
f1, f2, f3 = divideFace5(face)
subFaceIndices.append(f1)
subFaceIndices.append(f2)
subFaceIndices.append(f3)
if setColors : fColor += [faceColors[i]]*3
elif numVerts == 6 :
f1, f2, f3, f4 = divideFaces6(face)
subFaceIndices.append(f1)
subFaceIndices.append(f2)
subFaceIndices.append(f3)
subFaceIndices.append(f4)
if setColors : fColor += [faceColors[i]]*4
return np.array(subFaceIndices), fColor
def _get_edges_from_faces(self,fvIndx) :
"""
Used in __init__ if edgeIndices=None to create edgeIndices.
"""
# .............................................................
def edges_from_faceList( fvIndx ) :
# NOTE: very poor execution time. (~2 orders-of-magnitude slower)
# face index list of lists can't be used as a Numpy array
# since edge list of multi-type faces.
warnings.warn('multi-type polyhedron faces: use of non-optimized edge calc.')
keyList, edgeList = [],[]
strkey = lambda a,b : str(max(a,b)) + '_' + str(min(a,b))
for face in fvIndx :
numEdges = len(face)
for i in range(numEdges) :
a = face[i]
b = face[ (i+1)%numEdges]
key = strkey(a,b)
if key not in keyList :
keyList.append(key)
edgeList.append([a,b])
return np.array(edgeList)
# .............................................................
# following line only needed when called for mixed polygon faces
if fvIndx.ndim == 1 : return edges_from_faceList(fvIndx)
numEdges = fvIndx.shape[1]
indicies = [ (i+1)%numEdges for i in range(numEdges) ]
ofs_fvIndx = fvIndx[:,indicies]
stepA = np.stack( (fvIndx,ofs_fvIndx), axis=-1 )
stepB = np.reshape( stepA, (-1,2 ) )
stepC = np.sort(stepB,axis=1)
u, indices = np.unique(stepC, axis=0, return_index=True)
edgeList = stepC[indices]
return edgeList
@property
def _dfacecolors(self) :
"""
Direct access to the private property _facecolor.
Appears that the self.get_facecolor() needs a call to the
self.do_3d_projection() which will 'sort the 2D version by view depth'.
for render construction. Need the original order of colors.
BANDAID fix: update for v_1.2.0
FutDev: define a separate private property to hold these values
and not use the Matplotlib inherited value?
"""
return self._facecolors
@property
def _dedgecolors(self) :
"""
Direct access to the private property _edgecolor.
"""
return self._edgecolors
@property
def vlim(self) :
'''Range of values associated with color'''
return self._bounds['vlim']
@vlim.setter
def vlim(self,val) :
self._bounds['vlim'] = val
return
@property
def name(self) :
'''Descriptive identifier for the surface geometry.'''
if self._geomName is None : return ''
return self._geomName
@name.setter
def name(self, val) :
self._geomName = val
self.set_label( '' if val is None else val )
return
@property
def cname(self) :
'''Descriptive identifier for values indicated by color.'''
if self.valuesName is None : return ''
return self.valuesName
@cname.setter
def cname(self, val) :
self.valuesName = val
return
@property
def bounds(self) :
"""
Dictionary of surface geometric and value ranges.
Each dictionary value is a 2 float array of minimum
and maximum values of the surface. Keys are:
'xlim' : x-coordinate
'ylim' : y-coordinate
'zlim' : z-coordinate
'r_xy' : radial distance from the z axis
'rorg' : radial distance from the origin
'vlim' : value functional assignments.
Values are assigned from the geometry and color mapping methods,
including clipping.
"""
return self._bounds
@property
def cBar_ScalarMappable(self) :
"""matplotlib.cm.ScalarMappable object for surface values."""
vmin, vmax = self._bounds['vlim']
norm = colors.Normalize(vmin=vmin,vmax=vmax)
objMap = self.get_cmap()
sm = cm.ScalarMappable(cmap=objMap, norm=norm)
sm.set_array([])
return sm
@property
def edges(self) :
"""ColorLine3DCollection of the surface edges from facecolors."""
# ......................................................
def get_edge_color() :
efcList = copy.copy(self.efIndicesList)
# an edge may have 1 or 2 adjacent faces....
for i,cList in enumerate(self.efIndicesList) :
if len(cList)==1 : efcList[i].append(efcList[i][0])
efcArray = np.array(efcList)
self._postProc_surfaceColors(True)
colorA = self.facecolors[efcArray.T[0]]
colorB = self.facecolors[efcArray.T[1]]
# average of adjacent face RGB colors
return ( colorA+colorB)/2
# ......................................................
v = self.vertexCoor
e = self._get_edges_from_faces(self.fvIndices) # accounts for clipping...
lcol = ColorLine3DCollection(v,e,'edges')
lcol._ledgecolors = get_edge_color() # update for v_1.2.0
return lcol
@property
def initedges(self) :
''' ColorLine3DCollection of the initial surface edges.
ONLY available for surfaces where the number of
face edges differ among faces.
'''
if self._initedges is None :
raise ValueError('initedges property not available, all initial faces have equal number of edges (use edges property)')
e = self._initedges
v = self._initverts
lcol = ColorLine3DCollection(v,e,'edges')
return lcol
@property
def vertices(self) :
"""A 3 x N array of surface vertices."""
v = self.vertexCoor
m = np.transpose(v)
return m
@property
def facecenters(self) :
"""A 3 x N array of surface face center coordinates."""
verts = self.vertexCoor[ self.fvIndices ]
faceCenterCoor = self._get_face_centers(verts)
m = np.transpose(faceCenterCoor)
return m
@property
def facecolors(self) :
"""A N x 4 array of surface face colors."""
self._postProc_surfaceColors()
return self._dfacecolors # update for v_1.2.0
@property
def svd_dict(self) :
"""
A dictionary of results from a Singular Value Decomposition of the data.
Keys are:
'disarr' : array of N floats, normalized to the surface size.
'sigma' : standard deviation
'trans' : [rotation matrix, scaling, translation]
Values are assigned using the map_geom_from_svd method.
The dictionary is None if the method has not been called by the surface object.
"""
return self._svd_dict
@property
def area_h2b(self) :
""" A 2 x N array of normalized N face areas and shapes."""
f,v = self.fvIndices, self.vertexCoor
if f.shape[1] == 3 : # triangular
A = v[ f[:,1] ] - v[ f[:,0] ]
B = v[ f[:,2] ] - v[ f[:,0] ]
C = v[ f[:,2] ] - v[ f[:,1] ]
area = 0.5*np.linalg.norm( np.cross(A,B),axis=1 )
aveArea = np.average(area)
normArea = area/aveArea
edges = np.linalg.norm( np.array( [A,B,C]),axis=2)
maxsizes = np.amax(edges,axis=0)
f = 4*np.sqrt(3)/3 # equilateral triangle
skew = f*area/(maxsizes*maxsizes)
else : # quadrateral is two triangles.
b02 = v[ f[:,2] ] - v[ f[:,0] ]
b13 = v[ f[:,3] ] - v[ f[:,1] ]
base = np.array( [b02,b13] )
diag = np.linalg.norm( base ,axis=2)
i = np.where(diag[0]>diag[1], np.zeros(len(f),int), np.ones(len(f),int) )
maxsizes = diag[i,1]
A = b02
# top triangle ....
B = v[ f[:,3] ] - v[ f[:,0] ]
areaTop = 0.5*np.linalg.norm( np.cross(A,B),axis=1 )
# bottom triangle ....
B = v[ f[:,1] ] - v[ f[:,2] ]
areaBtm = 0.5*np.linalg.norm( np.cross(-A,B),axis=1 )
area = np.add( areaTop,areaBtm)
aveArea = np.average(area)
normArea = area/aveArea
f = 1 # square
skew = f*area/(maxsizes*maxsizes)
# ....................
return [normArea,skew]
[docs] def triangulate(self, rez=0) :
"""
PLANAR subdivision of each face into triangular faces.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of the base faces
into triangular faces, four faces per rez.
Rez values range from 0 to 7.
For surfaces with 4 vertices per face, each face is first
subdivided into two faces, then recursion proceeds.
For surfaces with 5 vertices per face, each face is first
subdivided into three faces, then recursion proceeds.
Returns
-------
self : surface object
"""
# .............................................................
def midVFun(vectA, vectB) :
mid = np.add(vectA,vectB)
mid = np.multiply(0.5,mid)
return mid
# .............................................................
self._postProc_surfaceColors(True)
if self.fvIndices.shape[1] > 3 :
# if needed, break faces into triangles before proceeding....
subDiv = self.fvIndices.shape[1] - 2
self.fvIndices = self._triangulate_faceIndices(self.fvIndices)
newFaceColors = []
for fColor in self._dfacecolors : # update for v_1.2.0
newFaceColors += [fColor]*subDiv
self.set_facecolor( np.array(newFaceColors) )
indexObj ,vertexCoor = self.triangulateBase(rez,list(self.vertexCoor),self.fvIndices, midVFun)
self.fvIndices = indexObj['face']
self.vertexCoor = vertexCoor
self.evIndices = self._get_edges_from_faces(self.fvIndices)
self.vfIndicesList = self._vertexCommonFaceIndices()
self.efIndicesList = self._edgeCommonFaceIndices()
v = self.vertexCoor
verts = v[self.fvIndices]
self.set_verts( verts )
self._set_geometric_bounds()
if rez != 0 :
newFaceColors = []
subDiv = 4**rez
for fColor in self._dfacecolors : # update for v_1.2.0
newFaceColors += [fColor]*subDiv
self.set_facecolor( np.array(newFaceColors) )
return self
[docs] def normalize_scale(self) :
"""Scaling normalization array and reciprocal."""
bounds = self.bounds
scX = ( bounds['xlim'][1] - bounds['xlim'][0] )
scY = ( bounds['ylim'][1] - bounds['ylim'][0] )
scZ = ( bounds['zlim'][1] - bounds['zlim'][0] )
recip = [scX,scY,scZ]
scale = np.reciprocal(recip)
return scale, recip
[docs] def vertexnormals(self, **kargs) :
"""
Unit directions or Vector3DCollection of vertex normals at vertices.
Parameters
----------
scale: number, optional
If not spectified, scaled proportional to the mean edge
length of surface base faces.
v3d: boolean, optional, default: True
If True, Vector3DCollection is returned, else a N x 3 array
of unit vector vertex normals is returned if False.
color : str or float list of length 3 or 4.
RGB or RGBA color, either a Matplotlib format string or a
list of color values in range [0,1].
width : number, optional, default: 1
Line width of the vector.
alr : scalar, optional, default: 0.25
Axis length ratio, head size to vector magnitude.
Returns
-------
Vector3DCollection object or array of unit direction vectors
"""
kargDft = { 'scale':self._normal_scale, 'v3d':True, 'color':None, 'width':1, 'alr':_DFT_ALR }
KW = KWprocessor(kargDft,'vectorfield_from_op',**kargs)
scale = KW.getVal('scale')
isV3D = KW.getVal('v3d')
color = KW.getVal('color')
width = KW.getVal('width')
alr = KW.getVal('alr')
verts = self.vertexCoor
norms = self._get_vertex_normals()
if not isV3D : return norms
name = 'vertex_normals'
lcol = Vector3DCollection(verts,norms*scale,alr,name,colors=color, linewidths=width)
lcol.coorType = self.coorType
lcol.uvwOrientationType = _COORSYS["XYZ"]
if color is None : lcol._set_vectColor( self.vertexColor )
return lcol
[docs] def facenormals(self, **kargs) :
"""
Unit directions or Vector3DCollection of face normals at face centers.
Parameters
----------
scale: number, optional
If not spectified, scaled proportional to the mean edge
length of surface base faces.
v3d: boolean, optional, default: True
If True, Vector3DCollection is returned, else a N x 3 array
of unit vector face normals is returned if False.
color : str or float list of length 3 or 4.
RGB or RGBA color, either a Matplotlib format string or a
list of color values in range [0,1].
width : number, optional, default: 1
Line width of the vector.
alr : scalar, optional, default: 0.25
Axis length ratio, head size to vector magnitude.
Returns
-------
Vector3DCollection object or array of unit direction vectors
"""
#Note: default color:None so that can use as a flag to use face colors if not defined.
kargDft = { 'scale':self._normal_scale, 'v3d':True, 'color':None, 'width':1, 'alr':_DFT_ALR }
KW = KWprocessor(kargDft,'vectorfield_from_op',**kargs)
scale = KW.getVal('scale')
isV3D = KW.getVal('v3d')
color = KW.getVal('color')
width = KW.getVal('width')
alr = KW.getVal('alr')
verts = self.vertexCoor[ self.fvIndices ]
faceCenterCoor = self._get_face_centers(verts)
faceNormalCoor = self._get_face_normals(verts)
if not isV3D : return faceNormalCoor
name = 'face_normals'
lcol = Vector3DCollection(faceCenterCoor,faceNormalCoor*scale,alr,name,colors=color, linewidths=width)
# --- 1.1 mod --------------------
lcol.coorType = self.coorType
lcol.uvwOrientationType = _COORSYS["XYZ"]
if color is None : lcol._set_vectColor(self._dfacecolors) # update for v_1.2.0
return lcol
[docs] def dispfield_from_op(self, operation, **kargs) :
"""
Vector3DCollection of displacement of vertices to a different position.
Parameters
----------
operation : function object
Function that takes one coordinate
argument, a 3xN Numpy array.
The function returns a 3xN array of vectors.
returnxyz : bool { True, False }, optional, default: False
By default, native coordinates are returned by the
operation function. If set True, the operation
returns xyz Cartesian coordinates.
useBase : bool { True, False }, optional, default: False
When set False, vertices of the surface, prior to
any geometric mapping or transforms, are passed to
the operation function. Otherwise, when True, the
current surface vertices are passed.
scale: number, optional, default: 1.
color : str or float list of length 3 or 4.
RGB or RGBA color, either a Matplotlib format string or a
list of color values in range [0,1].
width : number, optional, default: 1
Line width of the vector.
alr : scalar, optional, default: 0.25
Axis length ratio, head size to vector magnitude.
Returns
-------
Vector3DCollection object
"""
kargDft = { 'returnxyz':False, 'useBase':False,
'scale':1, 'color':'black', 'width':1, 'alr':_DFT_ALR }
KW = KWprocessor(kargDft,'dispfield_from_op',**kargs)
returnxyz = KW.getVal('returnxyz')
useBase = KW.getVal('useBase')
scale = KW.getVal('scale')
color = KW.getVal('color')
width = KW.getVal('width')
alr = KW.getVal('alr')
if self._isCompositeSurface :
warnings.warn('Vector operation not available for combined shapes.')
return
start_coor = self.vertexCoor
if useBase : start_coor = self.baseVertexCoor
v = self._transformVector(start_coor, operation, returnxyz)
delta = v - self.vertexCoor
delta = scale*delta
# --- 1.1 mod --------------------
vField = self._get_vectorfield_Line3DCollection(delta, color, width, alr)
if returnxyz : vField.uvwOrientationType = _COORSYS["XYZ"]
return vField
[docs] def vectorfield_from_op(self, operation, **kargs) :
"""
Vector3DCollection of vectors in u,v,w coordinates at surface vertices.
Parameters
----------
operation : function object
Function that takes one coordinate
argument, a 3xN Numpy array of native coordinates.
The function returns a 3xN array of vectors.
scale: number, optional, default: 1.
color : str or float list of length 3 or 4.
RGB or RGBA color, either a Matplotlib format string or a
list of color values in range [0,1].
width : number, optional, default: 1
Line width of the vector.
alr : scalar, optional, default: 0.25
Axis length ratio, head size to vector magnitude.
Returns
-------
Vector3DCollection object
"""
kargDft = { 'scale':1, 'color':'black', 'width':1, 'alr':_DFT_ALR }
KW = KWprocessor(kargDft,'vectorfield_from_op',**kargs)
scale = KW.getVal('scale')
color = KW.getVal('color')
width = KW.getVal('width')
alr = KW.getVal('alr')
if self._isCompositeSurface :
warnings.warn('Vector operation not available for combined shapes.')
return
xyz = np.transpose(self.vertexCoor)
abc = self.coor_convert(xyz)
abc = np.array(abc)
rst = np.array(operation(abc))
rst = scale*rst
delta = self._uvw_to_xyz(rst,abc)
# --- 1.1 mod --------------------
delta = np.array(delta)
return self._get_vectorfield_Line3DCollection(delta, color, width, alr)
[docs] def vectorfield_to_surface(self, surface, **kargs) :
"""
Vector3DCollection of vectors from surface to surface vertex coordinates.
Parameters
----------
surface : surface object
Surface that matches the calling surface vectors.
Should be the same basetype and rez for surface subclasses.
scale: number, optional, default: 1.
color : str or float list of length 3 or 4.
RGB or RGBA color, either a Matplotlib format string or a
list of color values in range [0,1].
width : number, optional, default: 1
Line width of the vector.
alr : scalar, optional, default: 0.25
Axis length ratio, head size to vector magnitude.
Raises
------
ValueError
Mismatched surfaces based on the number of vertices.
Returns
-------
Vector3DCollection object
"""
kargDft = { 'scale':1, 'color':'black', 'width':1, 'alr':_DFT_ALR }
KW = KWprocessor(kargDft,'vectorfield_to_surface',**kargs)
scale = KW.getVal('scale')
color = KW.getVal('color')
width = KW.getVal('width')
alr = KW.getVal('alr')
surLen, selfLen = len(surface.vertexCoor) , len(self.vertexCoor)
#if surLen is not selfLen : # 1.1 err. corrected
if surLen != selfLen :
raise ValueError('Surfaces have unequal number of vertices, {} != {}'.format(surLen, selfLen))
delta = surface.vertexCoor - self.vertexCoor
delta = scale*delta
# --- 1.1 mod --------------------
vField = self._get_vectorfield_Line3DCollection(delta, color, width, alr)
vField.uvwOrientationType = _COORSYS["XYZ"]
return vField
def _map_color_from_image_vert(self, img, viewport ) :
a,b,inViewport = self._viewportCoor(self.vertexCoor,viewport)
height, width = np.subtract(img.shape[:2],1)
M_index = ( (1-b)*height ).astype(int)
N_index = ( a*width ).astype(int)
#'''
orig_colors = self.vertexColor
if len(orig_colors) == 1 :
onearr = np.ones( len(self.vertexCoor) )[:, np.newaxis]
orig_colors = orig_colors*onearr
colorMap = []
oneAlp = np.array([1.0])
if viewport is None : colorMap = img[M_index,N_index]
else :
# FutDev: use numpy methods to optimize efficiency?
for i in range( len(orig_colors) ) :
colorMap.append(orig_colors[i])
if inViewport[i] :
temp = img[M_index[i],N_index[i]]
#.. note: img may or may not have an alpha channel..
temp2 = temp[0:3]
colorMap[i] = np.concatenate((temp2,oneAlp))
self.vertexColor = np.array(colorMap)
return
[docs] def map_color_from_image(self, fname, viewport=None) :
"""
Assign image color values to surface face colors.
This method is not available for composite surfaces.
The entire image domain is applied to the coloring operation.
Parameters
----------
fname : str or file-like
The image file to read: a filename, a URL or a file-like object opened
in read-binary mode. Currently restricted to PNG format.
viewport : 4D array-like, optional, default: None
Viewport defines the subdomain of the surface onto which
the image is mapped. The subdomain is set by normalized
native surface coordinates in a 4D array.
The entire surface is colored for the default of None.
Returns
-------
self : surface object
"""
if self._isCompositeSurface :
warnings.warn('Image operation not available for combined shapes.')
return self
if self._surfaceShapeChanged :
warnings.warn('Image operation may be anomalous after shape modification.')
if viewport is not None :
viewport = np.array(viewport)
if np.any( (viewport<0) | (viewport>1) ) :
raise ValueError('viewport list values, {}, must be between 0 and 1'.format(viewport))
if len(viewport) != 4 :
raise ValueError('viewport list must be length 4, found {}'.format(len(viewport)))
verts = self.vertexCoor[ self.fvIndices ]
faceCenterCoor = self._get_face_centers(verts)
a,b,inViewport = self._viewportCoor(faceCenterCoor,viewport)
img = image.imread(fname)
height, width = np.subtract(img.shape[:2],1)
M_index = ( (1-b)*height ).astype(int)
N_index = ( a*width ).astype(int)
orig_colors = self._dfacecolors # update for v_1.2.0
if len(orig_colors) == 1 :
onearr = np.ones( len(faceCenterCoor) )[:, np.newaxis]
orig_colors = orig_colors*onearr
colorMap = []
if viewport is None : colorMap = img[M_index,N_index]
else :
# FutDev: use numpy methods to optimize efficiency?
for i in range( len(orig_colors) ) :
colorMap.append(orig_colors[i])
if inViewport[i] : colorMap[i] = img[M_index[i],N_index[i]]
self.set_color(colorMap)
self.isSurfaceColored = True
self.vertexValues = None
# ========================================================================
self._map_color_from_image_vert(img,viewport )
# ========================================================================
return self
def _map_color_from_op_vert(self, operation, rgb=True) :
xyz = np.transpose(self.vertexCoor)
abc = self.coor_convert(xyz)
colors = np.array(operation(abc))
colors = np.transpose(colors)
colors = np.clip(colors,0,1)
if rgb : RGB = colors
else : RGB = cm.colors.hsv_to_rgb(colors)
if RGB.shape[1] == 3 : # add alpha channel.
ones = np.array( [np.ones( RGB.shape[0] )] ).T
RGB = np.concatenate( (RGB, ones ), axis=1 )
self.vertexColor = RGB
self._bounds['vertvlim'] = _DFT_VERTVLIM
return
[docs] def map_color_from_op(self, operation, rgb=True, cname=None) :
"""
Assignment of face color from a function.
Face colors are assigned from a function
of face coordinates.
Parameters
----------
operation : function object
Function that takes one argument,
a 3xN Numpy array of native coordinates.
The function returns a 3xN color value.
rgb : bool {True, False}, optional, default: True
By default, RGB color values are returned by the
operation function. If set False, the operation
returns HSV color values.
cname : string, optional, default: None or op function name.
Descriptive identifier for values indicated by color.
Returns
-------
self : surface object
"""
verts = self.vertexCoor[ self.fvIndices ]
faceCenterCoor = self._get_face_centers(verts)
xyz = np.transpose(faceCenterCoor)
abc = self.coor_convert(xyz)
colors = np.array(operation(abc))
colors = np.transpose(colors)
colors = np.clip(colors,0,1)
if rgb : RGB = colors
else : RGB = cm.colors.hsv_to_rgb(colors)
self.set_color(RGB)
self.isSurfaceColored = True
cname = _getFunctionName(operation,cname)
if cname is not None : self.valuesName = cname
self._bounds['vlim'] = _DFT_VLIM
self._map_color_from_op_vert(operation, rgb)
return self
def _map_cmap_from_datagrid_vert(self, datagrid, cmap=None, viewport=None) :
a,b,inViewport = self._viewportCoor(self.vertexCoor,viewport)
data = datagrid
dmax = np.amax(datagrid)
dmin = np.amin(datagrid)
xd = np.linspace(0, 1, data.shape[0] )
yd = np.linspace(0, 1, data.shape[1] )
g = interpolate.interp2d(yd, xd, data, kind='cubic')
d = []
# FutDev: use numpy methods to optimize efficiency.
for ab in np.transpose([a,b]) :
d.append( g(ab[0],ab[1])[0] )
d = np.where(inViewport,d,np.full(len(d),dmin))
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cmap = self.get_cmap()
orig_colors = self.vertexColor
if len(orig_colors) == 1 :
onearr = np.ones( len(self.vertexCoor) )[:, np.newaxis]
orig_colors = orig_colors*onearr
norm = colors.Normalize(dmin,dmax)
temp = cmap(norm(d))
if viewport is None :
colorMap = temp
else :
colorMap = np.where(inViewport[:, np.newaxis], temp, orig_colors )
# FutDev: set_color used instead of set_facecolor to eliminate 'gaps' between faces, <<<<<
# however, set edgecolor required after illumination if edges are to be displayed. <<<<<
# HOWEVER !!!! if colorMap has a alpha<1, these don't register with the edges. <<<<<
#self.vertexValues = d <<<< NO.. interpolated values only in viewport,
self.vertexColor = np.array(colorMap)
return
[docs] def map_cmap_from_datagrid(self, datagrid, cmap=None, viewport=None, cname=None) :
"""
Face color assignment using a 2D datagrid surface.
Datagrid values are normalized in the range 0 to 1.
This method is not available for composite surfaces.
The entire datagrid domain is applied to the geometric operation.
Parameters
----------
datagrid : 2D float array
cmap : str or Colormap, optional
A Colormap instance or registered colormap name.
If not assigned, the surface Colormap is used.
The colormap maps the datagrid values to colors.
viewport : 4D array-like, optional, default: None
Viewport defines the subdomain of the surface onto which
the datagrid is mapped. The subdomain is set by normalized
native surface coordinates in a 4D array.
The entire surface is colored for the default of None.
cname : string, optional, default: None.
Descriptive identifier for values indicated by color.
Returns
-------
self : surface object
"""
if self._isCompositeSurface :
warnings.warn('Datagrid operation not available for combined shapes.')
return self
if self._surfaceShapeChanged :
warnings.warn('Datagrid operation may be anomalous after shape modification.')
verts = self.vertexCoor[ self.fvIndices ]
faceCenterCoor = self._get_face_centers(verts)
a,b,inViewport = self._viewportCoor(faceCenterCoor,viewport)
data = datagrid
dmax = np.amax(datagrid)
dmin = np.amin(datagrid)
xd = np.linspace(0, 1, data.shape[0] )
yd = np.linspace(0, 1, data.shape[1] )
g = interpolate.interp2d(yd, xd, data, kind='cubic')
d = []
# FutDev: use numpy methods to optimize efficiency.
for ab in np.transpose([a,b]) :
d.append( g(ab[0],ab[1])[0] )
d = np.where(inViewport,d,np.full(len(d),dmin))
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cmap = self.get_cmap()
orig_colors = self._dfacecolors # update for v_1.2.0
if len(orig_colors) == 1 :
onearr = np.ones( len(faceCenterCoor) )[:, np.newaxis]
orig_colors = orig_colors*onearr
norm = colors.Normalize(dmin,dmax)
temp = cmap(norm(d))
if viewport is None :
colorMap = temp
else :
colorMap = np.where(inViewport[:, np.newaxis], temp, orig_colors )
# FutDev: set_color used instead of set_facecolor to eliminate 'gaps' between faces,
# however, set edgecolor required after illumination if edges are to be displayed.
# HOWEVER !!!! if colorMap has a alpha<1, these don't register with the edges.
self.set_color(colorMap)
self.isSurfaceColored = True
self._bounds['vlim'] = [ dmin, dmax ]
self.valuesName = cname
# ========================================================================
self._map_cmap_from_datagrid_vert(datagrid, cmap, viewport )
# ========================================================================
return self
def _getColorMap_fromNormals(self,fvo,cmap, direction, isAbs) :
""" used by map_cmap_from_normals() & _map_cmap_from_normals_vert() """
unitVector = lambda v : np.divide( v, np.linalg.norm(v) )
incidentLight = unitVector( direction )
d = np.dot(fvo,incidentLight)
if isAbs : v = 1 - np.abs(d)
else :
v = ( d + 1 )/2
v = np.abs(v) # to insure positive values near zero
colorMap = cmap(v)
return colorMap
def _map_cmap_from_normals_vert(self, cmap, direction, isAbs) :
vobj = self._get_vertex_normals()
cMap = self._getColorMap_fromNormals(vobj, cmap, direction, isAbs)
self.vertexColor = np.array(cMap)
return
[docs] def map_cmap_from_normals(self, cmap=None, direction=None, cname=None, isAbs=False) :
"""
Face color assignment using normals relative to a direction.
The dot product of face normals with the direction
is used to assign face colors from a colormap.
Parameters
----------
cmap : str or Colormap, optional
A Colormap instance or registered colormap name.
If not assigned, the surface Colormap is used.
The colormap maps the dot product values to colors.
direction : list of size 3, optional, default: [1,0,1]
A 3D vector in xyz Cartesian coordinates designating
the direction of the illumination source.
If assigned the Axes3D, will use the view direction.
cname : string, optional, default: None.
Descriptive identifier for values indicated by color.
isAbs : boolean, optional, default: False
If set True, the absolute value of the dot product
between face normals and illumination direction is
used.
Returns
-------
self : surface object
"""
if direction is not None:
temp = direction
try : # check if direction holds the elev, azim (ie. 3D axis)
elev,azim = temp.elev, temp.azim
direction = elev_azim_2vector(elev,azim)
except :
direction = temp
if len(direction) != 3 :
raise ValueError('Cmap from normals direction must be a array of length 3, or a Matplotlib Axes3D object.')
else : direction = _ILLUM
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cm_colorMap = self.get_cmap()
verts = self.vertexCoor[ self.fvIndices ]
fvobj = self._get_face_normals(verts)
cMap = self._getColorMap_fromNormals(fvobj, cm_colorMap, direction, isAbs)
self.set_cmap(cm_colorMap)
self.set_color(cMap)
self.isSurfaceColored = True
self.onlyColoredFaces = True
self.valuesName = cname
# ========================================================================
self._map_cmap_from_normals_vert( cm_colorMap, direction, isAbs )
# ========================================================================
return self
def _map_cmap_from_op_vert(self, operation, cmap) :
xyz = np.transpose(self.vertexCoor)
# .....
abc = self.coor_convert(xyz)
v = np.array(operation(abc))
self.vertexValues = v
norm = colors.Normalize(vmin=v.min(),vmax=v.max())
colorMap = cmap(norm(v))
self.vertexColor = colorMap
colorMap = np.array(colorMap)
self._bounds['vertvlim'] = [ v.min(), v.max() ]
return
[docs] def map_cmap_from_op(self, operation=None, cmap=None, cname=None) :
"""
Functional assignment of a color from a color map.
Face coordinates are used to calculate a scalar
which is then used to assign face colors from a
colormap.
Parameters
----------
operation : function object, default : None
Function that takes one argument,
a 3xN Numpy array of native coordinates.
The function returns a Numpy array of scalar values.
If function is None, function will map in the z-direction
or r-direction, dependent on the native coordinates.
cmap : str or Colormap, optional
A Colormap instance or registered colormap name.
If not assigned, the surface Colormap is used.
The colormap maps the function return values to colors.
cname : string, optional, default: None or op function name.
Descriptive identifier for values indicated by color.
Returns
-------
self : surface object
"""
dftDirName = None
if operation is None :
opDir = 2
dftDirName = 'Z-direction'
if self.coorType == _COORSYS["SPHERICAL"] : opDir = 0
if self.coorType == _COORSYS["CYLINDRICAL"] : opDir = 0
if opDir == 0 : dftDirName = 'R-direction'
operation = lambda c : c[opDir]
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cmap = self.get_cmap()
verts = self.vertexCoor[ self.fvIndices ]
faceCenterCoor = self._get_face_centers(verts)
xyz = np.transpose(faceCenterCoor)
# .....
abc = self.coor_convert(xyz)
v = np.array(operation(abc))
self.faceValues = v
vmin, vmax = v.min(), v.max()
norm = colors.Normalize(vmin=v.min(),vmax=v.max())
colorMap = cmap(norm(v))
# FutDev: set_color used instead of set_facecolor to eliminate 'gaps' between faces,
# however, set edgecolor required after illumination if edges are to be displayed.
# HOWEVER !!!! if colorMap has a alpha<1, these don't register with the edges.
self.set_color(colorMap)
self.isSurfaceColored = True
self._bounds['vlim'] = [ vmin, vmax ]
#cname = _getFunctionName(operation,cname)
trialName = None
if cname is None :
if dftDirName is not None : trialName = dftDirName
else : trialName = _getFunctionName(operation,None)
else: trialName = cname
self.valuesName = trialName
# ========================================================================
self._map_cmap_from_op_vert(operation, cmap)
# ========================================================================
return self
[docs] def map_geom_from_datagrid(self,datagrid, scale=1.0, viewport=None, name=None ) :
"""
Append surface vertices proportional to 2D datagrid surface.
Datagrid values are normalized in the range 0 to 1.
This method is not available for composite surfaces.
The entire datagrid domain is applied to the geometric operation.
Parameters
----------
datagrid : 2D array
scale : number, optional, default: 1
viewport : 4D array-like, optional, default: None
Viewport defines the subdomain of the surface onto which
the datagrid is mapped. The subdomain is set by normalized
native surface coordinates in a 4D array.
The entire surface is geometrically mapped for the default of None.
name : string, optional, default: None
Descriptive identifier for the geometry.
Returns
-------
self : surface object
"""
if self._isCompositeSurface :
warnings.warn('Datagrid operation not available for combined shapes.')
return self
if self._surfaceShapeChanged :
#warnings.warn('Datagrid operation may be anomalous after shape modification.')
warnings.warn('Datagrid operation not available after shape modification.')
return self
a,b,inViewport = self._viewportCoor(self.vertexCoor,viewport)
data = datagrid
dmax = np.amax(datagrid)
dmin = np.amin(datagrid)
delta = dmax-dmin
data = (datagrid - dmin)/delta
xd = np.linspace(0, 1, data.shape[0] )
yd = np.linspace(0, 1, data.shape[1] )
g = interpolate.interp2d(yd, xd, data, kind='cubic')
# FutDev: optimize efficiency (?).
d = []
for ab in np.transpose([a,b]) :
d.append( g(ab[0],ab[1])[0] )
d = np.where(inViewport,d,np.zeros(len(d)))
d = d[:, np.newaxis]
# FutDev: 'fVal' only applicable for initial subclassed
# surfaces to get surface normals at a vertex.
# Future work: to make this generic for any surface.
fVal = np.array([0.,0.,1.])
if self.disp_Vector is not None :
fVal = self.vertexCoor*self.disp_Vector
v = self.vertexCoor + scale*d*fVal
verts = v[self.fvIndices]
self.vertexCoor = v
self._set_geometric_bounds()
self.set_verts( verts )
self._surfaceShapeChanged = True
if name is not None: self._geomName = name
return self
[docs] def map_geom_from_image(self, fname, scale=1.0, viewport=None, cref='v', hzero=0, name=None ) :
"""
Append surface vertices proportional to image color component values.
For planar and polar coordinates, the z-coordinate is appended.
For cylindrical and spherical coordinates, the r-coordinate is appended.
This method is not available for composite surfaces.
The entire image domain is applied to the geometric operation.
Parameters
----------
fname : str or file-like
The image file to read: a filename, a URL or a file-like object opened
in read-binary mode. Currently restricted to PNG format.
scale : number
cref : {'r','g','b','h','s','v'}, optional, default: 'v'
Sets which color tuple value to use for geometry mapping.
Note: only the first character of the string is evaluated.
hzero : number or None, optional, default: 0
Argument is used if cref is set to 'h'.
The hzero magnitude indicates the Hue value for a zero coordinate
diplacement. The hzero sign (positive or negative),
indicates the direction of increasing value in the
appended coordinate with hue value. Range of hzero is [-1,1].
viewport : 4D array-like, optional, default: None
Viewport defines the subdomain of the surface onto which
the image is mapped. The subdomain is set by normalized
native surface coordinates in a 4D array.
The entire surface is geometrically mapped for the default of None.
name : string, optional, default: None
Descriptive identifier for the geometry.
Returns
-------
self : surface object
"""
# ..............................................................
def getIndex(cref) :
cref = cref[0].lower()
valIndex = 'rgbhsv'.find(cref)
if valIndex < 0 :
raise ValueError('Invalid cref argument: ' + str(cref))
isHSV = False
if valIndex > 2 :
valIndex -= 3
isHSV = True
return valIndex, isHSV
# ..............................................................
if self._isCompositeSurface :
warnings.warn('Image operation not available for combined shapes.')
return self
if self._surfaceShapeChanged :
warnings.warn('Image operation may be anomalous after shape modification.')
if (hzero<-1) or (hzero>1) :
raise ValueError('hzero must be between -1 and 1, found {}'.format(hzero))
ho = (1.0 -hzero)
hsgn = np.sign(hzero)
h_ref = lambda h : np.mod( ho + hsgn*h , 1.0)
# ...................
a,b,inViewport = self._viewportCoor(self.vertexCoor,viewport)
img = image.imread(fname)
height, width = np.subtract(img.shape[:2],1)
x_index = ( (1-b)*height ).astype(int)
y_index = ( a*width ).astype(int)
valIndex, isHSV = getIndex(cref)
rgbcolor = img[x_index,y_index]
if isHSV :
hsvcolor = colors.rgb_to_hsv(rgbcolor)
hsvcolor = np.transpose(hsvcolor)
if valIndex==0 : d= h_ref(hsvcolor[0])
else : d = hsvcolor[valIndex]
else :
rgbcolor = np.transpose(rgbcolor)
d = rgbcolor[valIndex]
d = np.where(inViewport,d,np.zeros(len(d)))
d = d[:, np.newaxis]
# FutDev: 'fVal' only applicable for initial subclassed surfaces
# surfaces to get surface normals at a vertex.
# Future work: to make this generic for any surface.
fVal = np.array([0.,0.,1.])
if self.disp_Vector is not None :
fVal = self.vertexCoor*self.disp_Vector
v = self.vertexCoor + scale*d*fVal
verts = v[self.fvIndices]
self.vertexCoor = v
self._set_geometric_bounds()
self.set_verts( verts )
self._surfaceShapeChanged = True
if name is not None: self._geomName = name
return self
[docs] def map_geom_from_op(self, operation, returnxyz=False, name=None) :
"""
Functional transformation of surface vertex coordinates.
This method is not available for composite surfaces.
Parameters
----------
operation : function object
Coordinate mapping function that takes one
argument, a 3xN Numpy array of native coordinates.
The function returns a 3xN array of coordinates.
returnxyz : bool { True, False }, optional, default: False
By default, native coordinates are returned by the
operation function. If set True, the operation
returns xyz Cartesian coordinates.
name : string, optional, default: None or op function name.
Descriptive identifier for the geometry.
Returns
-------
self : surface object
"""
return self._map_geom_from_op( operation, returnxyz, name, check=True)
def _map_geom_from_op(self, operation, returnxyz, name, check) :
# used by domain with check=False to allow composites scaling.
# otherwise, not allowed by map_geom_from_op method.
if self._isCompositeSurface and check :
warnings.warn('Map operation MUST use xyz coordinates for combined shapes.')
#warnings.warn('Map operation not available for combined shapes.')
#return self
v = self._transformVector(self.vertexCoor, operation, returnxyz)
verts = v[self.fvIndices]
self.vertexCoor = v
self._set_geometric_bounds()
self.set_verts( verts )
self._surfaceShapeChanged = True
if name is None : name = self._geomName
name = _getFunctionName(operation,name)
if name is not None: self._geomName = name
return self
[docs] def clip(self, operation, usexyz=False) :
"""
Remove faces from the surface based on position.
Parameters
----------
operation : function object
Function that takes one argument, a 3xN Numpy array of coordinates.
The function returns a bool { True, False } indicating if
the face, at the face centered coordinate, is to be
retained (True), otherwise, the face is removed from the
surface (False).
usexyz : bool { True, False }, default: False
If True, face centers are passed to the operation function
in xyz coordinates. If False, face centers are passed
in native coordinates.
Returns
-------
self : surface object
"""
# Clip operation MUST use xyz coordinates for combined shapes.
# since native coordinate may be undefined.
if self._isCompositeSurface : usexyz = True
verts = self.vertexCoor[ self.fvIndices ]
faceCenterCoor = self._get_face_centers(verts)
self._postProc_surfaceColors()
orig_colors = self._dfacecolors # update for v_1.2.0
xyz = np.transpose(faceCenterCoor)
if usexyz is False : xyz = self.coor_convert(xyz)
abc = xyz
shouldKeep = np.array(operation(abc))
if np.any(shouldKeep) is None :
warnings.warn('WARNING: Clipping resulted in no faces found, surface return unclipped.')
return self
fvi = self.fvIndices[shouldKeep]
fcol = orig_colors[shouldKeep]
# this permits multiple clips without changing original, so restore is available.
if self.restoredClip is None :
self.restoredClip = { 'faces' : self.fvIndices, 'colors' : orig_colors }
fvi = np.array(fvi)
self.fvIndices = fvi
self.set_color(fcol)
v = self.vertexCoor
vfvi = v[fvi]
self.set_verts( vfvi )
clipVerts = np.reshape(vfvi,(-1,3))
_set_geomBounds( clipVerts,self.bounds )
self._set_index_relations()
return self
[docs] def clip_alpha(self,alphaCut,useval=False) :
"""
Remove faces from the surface based on alpha transparency.
Parameters
----------
alphaCut : scalar
If the face color alpha is greater than alphaCut, the face is
retained, otherwise, the face is removed from the surface.
useval : bool { True, False }, default: False
If True, clipping based on color HSV value. Otherwise,
clipping is based on the color alpha value.
Returns
-------
self : surface object
"""
self._postProc_surfaceColors()
orig_colors = self._dfacecolors # update for v_1.2.0
if useval :
orig_hsv = colors.rgb_to_hsv(orig_colors[:,:3])
alphas = np.ravel( orig_hsv[:,2:3] )
else :
alphas = np.ravel( orig_colors[:,3:4] )
cuts = np.full(len(alphas),alphaCut)
shouldKeep = alphas>cuts
if np.any(shouldKeep) is None :
warnings.warn('WARNING: Alpha clipping resulted in no faces found, surface return unclipped.')
return self
fvi = self.fvIndices[shouldKeep]
fcol = orig_colors[shouldKeep]
# this permits multiple clips without changing original, so restore is available.
if self.restoredClip is None :
self.restoredClip = { 'faces' : self.fvIndices, 'colors' : orig_colors }
fvi = np.array(fvi)
self.fvIndices = fvi
self.set_color(fcol)
v = self.vertexCoor
vfvi = v[fvi]
self.set_verts( vfvi )
clipVerts = np.reshape(vfvi,(-1,3))
_set_geomBounds( clipVerts,self.bounds )
self._set_index_relations()
return self
[docs] def clip_plane(self,dist,**kargs) :
"""
Remove faces from the surface based on a clip surface.
Parameters
----------
dist : float
Distance from the origin to the clip surface,
along the direction vector.
direction : array-like, optional, default: [0,0,1]
A xyz vector normal to the intersection plane
for a planar clip surface or axial direction of
a cylinder for a cylindrical clip surface.
coor : integer or string indicating the type of clip surface:
0, p, P, xyz,planar - planar (default)
1, c, C, cylinder,pplar,cylindrical - cylinder
2, s, S, sphere,spherical - sphere
3, x, X - y-z plane
4, y, Y - x-z plane
5, z, Z - x-y plane
Returns
-------
self : surface object
"""
# ...........................................................
def clip_op(xyz,dist,coor,direction) :
abc = np.transpose(xyz)
dotprod = _coor_dotprod(direction,coor,abc)
if coor > 0 :
if dist < 0 :
shouldclip = np.greater(dotprod,np.full(len(dotprod),-dist))
else :
shouldclip = np.less(dotprod,np.full(len(dotprod),dist))
else :
shouldclip = np.less(dotprod,np.full(len(dotprod),dist))
return shouldclip
# ...........................................................
dirDft = { 'direction': [0,0,1.0] }
kargDft = { **_COOR_KWARGS, **dirDft }
KW = KWprocessor( kargDft, 'clip_plane', **kargs)
coor = KW.getVal('coor')
direction = KW.getVal('direction')
if coor==3 : coor, direction = 0, [ 1,0,0 ]
if coor==4 : coor, direction = 0, [ 0,1,0 ]
if coor==5 : coor, direction = 0, [ 0,0,1 ]
return self.clip( lambda xyz : clip_op(xyz,dist,coor,direction),True)
[docs] def clip_normals(self,direction=None) :
"""
Remove faces from the surface based on a direction relative to face normals.
Parameters
----------
direction : array-like or 3Daxis, optional, default: (30,-60)
The vector direction for viewing (elev,azim) or 3Daxis.
Returns
-------
self : surface object
"""
if direction is not None:
temp = direction
try : # check if direction holds the elev, azim (ie. 3D axis)
elev,azim = temp.elev, temp.azim
direction = np.array( elev_azim_2vector(elev,azim) )
except :
direction = temp
if len(direction) != 3 :
raise ValueError('Clip from normals direction must be a array of length 3, or a Matplotlib Axes3D object.')
else : direction = np.array( elev_azim_2vector( *_DFT_VIEW ) )
verts = self.vertexCoor[ self.fvIndices ]
faceNormalCoor = self._get_face_normals(verts)
self._postProc_surfaceColors()
orig_colors = self._dfacecolors # update for v_1.2.0
shouldKeep = np.dot(faceNormalCoor,direction) > 0.0
if np.any(shouldKeep) is None :
warnings.warn('WARNING: Normal clipping resulted in no faces found, surface return unclipped.')
return self
fvi = self.fvIndices[shouldKeep]
fcol = orig_colors[shouldKeep]
# this permits multiple clips without changing original, so restore is available.
if self.restoredClip is None :
self.restoredClip = { 'faces' : self.fvIndices, 'colors' : orig_colors }
fvi = np.array(fvi)
self.fvIndices = fvi
self.set_color(fcol)
v = self.vertexCoor
vfvi = v[fvi]
self.set_verts( vfvi )
clipVerts = np.reshape(vfvi,(-1,3))
_set_geomBounds( clipVerts,self.bounds )
self._set_index_relations()
return self
def _restore_from_clip(self) :
"""
Reset face and colors prior to any clip operation.
(in development)
"""
# FutDev: this method needs thorough testing.
# probably doesn't work as is.
if self.restoredClip is None :
warnings.warn('There is no clipped surface to restore.')
return self
fvi = self.restoredClip['faces']
fcol = self.restoredClip['colors']
self.restoredClip = None
self.fvIndices = fvi
self.set_color(fcol)
v = self.vertexCoor
self.set_verts(v[ np.array(fvi) ])
self._set_geometric_bounds() # reset to original bounds...
return self
[docs] def set_surface_alpha(self, alpha, constant=False, adjustlw=True) :
"""
Adjust the face color alpha values and linewidth of the surface.
Parameters
----------
alpha : scalar
Alpha is in the range 0 to 1.
constant : bool { True, False }, optional, False
If False, face color values are multiplied by
alpha. If True, all face colors alpha channels
are assigned to alpha.
adjustlw : bool { True, False }, optional, True
If True, the surface linewidth will be set
to make the edges visually transparent with
the faces for Matplotlib rendering.
Returns
-------
self : surface object
"""
if (alpha<0.0) or (alpha>1.0) :
raise ValueError('surface alpha must be between 0 and 1, found {}'.format(alpha))
self._postProc_surfaceColors()
orig_colors = self._dfacecolors # update for v_1.2.0
if constant :
np.put_along_axis(orig_colors, np.array([[3]]), alpha, axis=1)
colorMap = orig_colors
else:
colorMap = np.multiply(orig_colors,np.array([1.0,1.0,1.0,alpha]))
self.set_color(colorMap)
if adjustlw :
# heuristic approach to hiding face edges. BANDAID
lw = 0.0063*np.exp(4.2*alpha)
self.set_linewidth(lw)
return self
[docs] def shade(self, depth=0, direction=None, contrast=None, isAbs=False, ax=None, rview=False) :
"""
Reduce surface HSV color Value based on face normals.
The dot product of face normals with the illumination direction
is used to adjust face HSV color value.
Parameters
----------
depth : scalar, optional, default: 0
Minimum color value of shaded surface faces.
Depth value ranges from 0 to 1.
direction : array-like, optional, default: (1,0,1)
A xyz vector pointing to the illumination source.
contrast : scalar, optional, default: 1
Shading contrast adjustment from low to high with a value
of 1 for linear variations with the normal direction.
Contrast value ranges from 0.1 to 3.
isAbs : bool, optional, default: False
If True, the absolute value of the dot product is
is used to determine color value.
ax : Matplotlib 3D axes.
rview : boolean, default: False
If True, direction is relative to the view,
otherwise, relative to the axes.
Returns
-------
self : surface object
"""
if direction is None : direction = _ILLUM
viewdir,ausf = None,None
if ax is not None:
ausf = _axisUSF(ax)
viewdir = elev_azim_2vector(ax.elev,ax.azim)
if rview : direction = rtv(direction,ax.elev,ax.azim)
if np.any( (depth<0) | (depth>1) ) :
raise ValueError('depth values, {}, must be between 0 and 1'.format(depth))
if contrast is not None:
if (contrast<0.1 or contrast>3) :
raise ValueError('contrast must be between 0.1 and 3. , found {}'.format(contrast))
verts = self.vertexCoor[ np.array(self.fvIndices) ]
norms = self._get_face_normals(verts,ausf)
unitDirection = np.divide( direction, np.linalg.norm(direction) )
dprod = np.dot( norms, unitDirection )
if viewdir is not None :
vdot = np.dot(norms,viewdir)
dprod = np.where(np.less(vdot,[0.0]),-dprod,dprod)
# 0<d<1, heuristic function for
# normalized domain of effective dot product, d=f(dprod)
if isAbs : d = 1 - np.abs(dprod)
else :
d = ( dprod + 1 )/2
d = np.abs(d) # to insure positive values near zero
if contrast is not None :
d = 0.5*( 1 + np.sign(dprod)*np.power( np.abs(dprod) , 1/contrast) )
# extract HSV, leaving alpha unchanged.
self._postProc_surfaceColors()
orig_colors = self._dfacecolors # update for v_1.2.0
alphas = orig_colors[:,3:4]
fc_less_alpha = orig_colors[:,:3]
hsv_vals = cm.colors.rgb_to_hsv(fc_less_alpha)
# adjust color value with multiplier.
cvm = ((1-depth)*d + depth*np.ones(len(d)))[:, np.newaxis]
hsv_vals[:,2] = (cvm*hsv_vals[:,2:3])[:,0]
rgb_vals = cm.colors.hsv_to_rgb(hsv_vals)
shade_colors = np.concatenate((rgb_vals, alphas), axis=1)
self.set_color(shade_colors) # NOTE: this will also set appropriate edge colors.
return self
[docs] def hilite(self, height=1, direction=None, focus=None, ax=None, rview=False) :
"""
Increase surface HSV color Value and reduce Saturation, based on face normals.
The dot product of face normals with the illumination direction
is used to adjust face HSV color value and saturation.
Faces having normals with a negative component in the direction of
illumination are not affected. The reduction of color saturation is
proportional to the increase in color value.
Parameters
----------
height : scalar, optional, default: 1
Maximum color value of highlighted surface faces.
Height value ranges from 0 to 1.
direction : array-like, optional, default: (1,0,1)
A xyz vector pointing to the illumination source.
focus : scalar, optional, default: 1
Highlighting focus adjustment from low to high with a value
of 1 for quadratic variations with the normal direction.
Focus value ranges from 0.1 to 3.
ax : Matplotlib 3D axes.
rview : boolean, default: False
If True, direction is relative to the view,
otherwise, relative to the axes.
Returns
-------
self : surface object
"""
if direction is None : direction = _ILLUM
viewdir,ausf = None,None
if ax is not None:
ausf = _axisUSF(ax)
viewdir = elev_azim_2vector(ax.elev,ax.azim)
if rview : direction = rtv(direction,ax.elev,ax.azim)
if np.any( (height<0) | (height>1) ) :
raise ValueError('height values, {}, must be between 0 and 1'.format(height))
if focus is not None:
if (focus<0.1 or focus>3) :
raise ValueError('focus must be between 0.1 and 3. , found {}'.format(focus))
verts = self.vertexCoor[ self.fvIndices ]
norms = self._get_face_normals(verts,ausf)
unitDirection = np.divide( direction, np.linalg.norm(direction) )
dprod = np.dot( norms, unitDirection )
if viewdir is not None :
vdot = np.dot(norms,viewdir)
dprod = np.where(np.less(vdot,[0.0]),-dprod,dprod)
# 0<d<1, heuristic function for
# normalized domain of effective dot product, d=f(dprod)
d = np.where(dprod<0 , 0, dprod)
if focus is None : focus = 1
d0 = np.power(focus,1/2)/2.0
y = (d-d0)/(1-d0)
d = np.where(y<0 , 0, y)
d = np.power(d,1.0+2*focus)
# extract HSV, leaving alpha unchanged.
self._postProc_surfaceColors()
orig_colors = self._dfacecolors # update for v_1.2.0
alphas = orig_colors[:,3:4]
fc_less_alpha = orig_colors[:,:3]
hsv_vals = cm.colors.rgb_to_hsv(fc_less_alpha)
# adjust color saturation and value, linear with normalized domain.
hd = height*d
omhd = (np.ones( len(norms) ) - hd)
omhd = omhd[:,np.newaxis]
hsv_vals[:,1] = (omhd*hsv_vals[:,1:2])[:,0]
hsv_vals[:,2] = (omhd*hsv_vals[:,2:3])[:,0]
hsv_vals[:,2] += hd
rgb_vals = cm.colors.hsv_to_rgb(hsv_vals)
hilite_colors = np.concatenate((rgb_vals, alphas), axis=1)
self.set_color(hilite_colors)
return self
[docs] def fade(self,depth=0,elev=None,azim=None,ax=None) :
"""
Reduce surface face opacity based on face position relative
to the view orientation.
Parameters
----------
depth : scalar, optional, default: 0
Minimum opacity to 1 for face opacity from back
to front face center position.
Depth value ranges from 0 to 1.
elev : scalar, optional, default: 30
Elevation of the axes view.
azim : scalar, optional, default: -60
Azimuth of the axes view.
ax: Matplotlib 3D axes.
If not None, elev and azim are taken from the ax.
Returns
-------
self : surface object
"""
if ax is None :
if elev is None: elev = _DFT_VIEW[0]
if azim is None: azim = _DFT_VIEW[1]
else :
elev = ax.elev
azim = ax.azim
direction = elev_azim_2vector(elev, azim)
unitDirection = np.divide( direction, np.linalg.norm(direction) )
if np.any( (depth<0) or (depth>1) ) :
raise ValueError('depth values, {}, must be between 0 and 1'.format(depth))
verts = self.vertexCoor[ np.array(self.fvIndices) ]
faceCenterCoor = self._get_face_centers(verts)
self._postProc_surfaceColors()
colors = self._dfacecolors # update for v_1.2.0
if len(faceCenterCoor) == 1 : return self # fade not available for single face line
dtprod = np.dot(faceCenterCoor,unitDirection)
dprange = np.amax(dtprod)-np.amin(dtprod)
norm_dp = (dtprod - np.amin(dtprod))/dprange
fadex = (1-depth)*norm_dp + depth
orig_colors = colors
fc_less_alpha = orig_colors[:,:3]
fadex = fadex*orig_colors[:,3] # update for v_1.2.0
faded_colors = np.concatenate((fc_less_alpha, fadex[:,np.newaxis] ), axis=1)
self.set_color(faded_colors)
# heuristic approach to hiding face edges. BANDAID update for v_1.2.0
alpha = (depth+1)/2 # take the midpoint of the gradient.
lw = 0.0063*np.exp(4.2*alpha)
self.set_linewidth(lw)
return self
[docs] def map_geom_from_svd(self, data, pc=None ) :
"""
Transform surface geometry based on a PCA analysis of a data set.
Result values contained in the surface property svd_dict.
Parameters
----------
data : N x 3 array of x,y,z data coordinates.
pc : scalar, optional, default: None
Sets the minimum size of the surface based on the
relative number of data. Value range from 0 to 1,
with 1 for 100%. If None, one standard deviation
of the data sets the surface size.
Returns
-------
self : surface object
"""
if self._isCompositeSurface :
warnings.warn('Data operation not available for combined shapes.')
return self
data_mean = data.mean(axis=0)
data_centered = data - data_mean
Ua, sa, Vta = np.linalg.svd(data_centered)
dataAve = np.mean(sa)
scale = sa/dataAve
tranData = np.divide(np.dot(data_centered,Vta.T),scale )
disArr=np.linalg.norm(tranData,axis=1)
disArrSum = np.linalg.norm(disArr)
lendisArr = len(disArr)
refLen = disArrSum/np.sqrt(lendisArr)
s=refLen
disArr = np.divide(disArr,refLen)
if pc is not None:
dlen = len(disArr)
dindex = int(dlen*pc)
if dindex >= dlen : dindex = dlen - 1
if dindex < 0 : dindex = 0
s = np.sort(disArr)[dindex]
disArr = disArr/s
s = s*refLen
self.transform(Vta,s*scale,data_mean)
self._svd_dict = {
'disarr': disArr, 'sigma': refLen, 'trans': [Vta,s*scale,data_mean]
}
return self
[docs] def contourLines(self,*dist,**kargs) :
"""
ColorLine3DCollection contour lines on the surface.
Parameters
----------
dist : floats
Distances from the origin to the surface of intersection,
along the direction vector.
direction : array-like, optional, default: [0,0,1]
A xyz vector normal to the intersection planes
for planar contours or axial direction of
cylinders for cylindrical contours.
name : string identifier (default None)
color : color (default is the surface colors)
coor : integer or string indicating the type of contour surface:
0, p, P, xyz,planar - planar (default, planar & polar)
1, c, C, cylinder,pplar,cylindrical - cylinder (default, cylindrical)
2, s, S, sphere,spherical - sphere (default, spherical)
3, x, X - y-z plane
4, y, Y - x-z plane
5, z, Z - x-y plane
Returns
-------
self : line object
"""
extDft = { 'direction':[0,0,1.0], 'name': None, 'color': None }
kargDft = { **_COOR_KWARGS, **extDft }
# reset default 'coor' for cylindrical and spherical coor surfaces.
coorVal = 0
if self.coorType == _COORSYS["CYLINDRICAL"] : coorVal = 1
if self.coorType == _COORSYS["SPHERICAL"] : coorVal = 2
if coorVal != 0 :
temp = copy.deepcopy(_COOR_KWARGS)
temp['coor'][0] = [coorVal]
kargDft = { **temp, **extDft }
KW = KWprocessor(kargDft,'contourLines',**kargs)
cIndex = KW.getVal('coor')
direction = KW.getVal('direction')
name = KW.getVal('name')
color = KW.getVal('color')
if direction is None : direction = extDft['direction']
if cIndex==3 : cIndex, direction = 0, [ 1,0,0 ]
if cIndex==4 : cIndex, direction = 0, [ 0,1,0 ]
if cIndex==5 : cIndex, direction = 0, [ 0,0,1 ]
dotprod = _coor_dotprod(direction,cIndex,self.vertexCoor)
conLines = self._contourLineCollection( dotprod,direction,*dist )
if conLines is None :
raise ValueError('No surface.contourLines were found.')
if cIndex == 0 : conLines._planeNormal = direction
if name is not None : conLines.name = kargs['name']
if color is not None : conLines.set_color(kargs['color'])
return conLines
[docs] def contourLineSet(self,numb=2,**kargs) :
"""
ColorLine3DCollection of evenly spaced contour lines on the surface.
Parameters
----------
numb : integer
Number of contour lines intersecting a surface.
Numb value must be getter than zero.
direction : array-like, optional, default: [0,0,1]
A xyz vector normal to the intersection planes
for planar contours or axial direction of
cylinders for cylindrical contours.
name : string identifier (default None)
color : color (default is the surface colors)
coor : integer or string indicating the type of contour surface:
0, p, P, xyz,planar - planar (default, planar & polar)
1, c, C, cylinder,pplar,cylindrical - cylinder (default, cylindrical)
2, s, S, sphere,spherical - sphere (default, spherical)
3, x, X - y-z plane
4, y, Y - x-z plane
5, z, Z - x-y plane
Returns
-------
self : line object
"""
extDft = { 'direction':[0,0,1.0], 'name': None, 'color': None }
kargDft = { **_COOR_KWARGS, **extDft }
# reset default 'coor' for cylindrical and spherical coor surfaces.
coorVal = 0
if self.coorType == _COORSYS["CYLINDRICAL"] : coorVal = 1
if self.coorType == _COORSYS["SPHERICAL"] : coorVal = 2
if coorVal != 0 :
temp = copy.deepcopy(_COOR_KWARGS)
temp['coor'][0] = [coorVal]
kargDft = { **temp, **extDft }
KW = KWprocessor(kargDft,'contourLineSet',**kargs)
cIndex = KW.getVal('coor')
direction = KW.getVal('direction')
name = KW.getVal('name')
color = KW.getVal('color')
dotprod = _coor_dotprod(direction,cIndex,self.vertexCoor)
maxd = np.amax(dotprod)
mind = np.amin(dotprod)
convals = np.linspace(mind,maxd,numb+2)[1:numb+1]
convals = tuple(convals)
conLines = self._contourLineCollection( dotprod,direction,convals )
if cIndex == 0 : conLines._planeNormal = direction
if name is not None : conLines.name = kargs['name']
if color is not None : conLines.set_color(kargs['color'])
return conLines
[docs] def evert(self) :
"""
Reverse direction of face normals.
"""
self.fvIndices = np.flip(self.fvIndices, axis=1)
return self
[docs] def domain(self, xlim=None, ylim=None, zlim=None) :
"""
Set the domain of the surface.
Used for setting the domain of the base.
Parameters
----------
xlim, ylim, zlim : 2d arrays
Minimum and maximum values of the arrays are used to scale
and translate the surface from an intial domains of
[ -1, 1 ]
If set to None, the domains are unchanged.
Returns
-------
self : surface object
"""
def setDomain(xyz, mx,bx, my,by, mz,bz) :
x,y,z = xyz
X = mx*x + bx
Y = my*y + by
Z = mz*z + bz
return X,Y,Z
def getCont( lims ) :
m, b = 1.0, 0.0
if lims is not None :
# allow single bndry to be input as a float
if isinstance(lims, (int,float)) : lims = [-lims,lims]
m, b = 0.5*(lims[1]-lims[0]) , 0.5*(lims[1]+lims[0])
return m, b
mx, bx = getCont(xlim)
my, by = getCont(ylim)
mz, bz = getCont(zlim)
self._map_geom_from_op(lambda A : setDomain(A, mx,bx, my,by, mz,bz),False,None,False)
return self
[docs]class CubicSurface(Surface3DCollection) :
"""
Cubic 3D surface in Cartesian coordinates with rectangular faces.
Methods are inherited from the Surface3DCollection.
"""
@staticmethod
def _get_cubic_surface_dictionary() :
""" Constructs the dictionary of base planar surface geometries. """
surfacesDictionary = {}
baseVert = []
baseVert = [ [ 1,-1, -1], [ 1, 1, -1], [-1, 1, -1], [-1,-1, -1],
[ 1,-1, 1], [ 1, 1, 1], [-1, 1, 1], [-1,-1, 1] ]
cubeFaces=(
[0,1,5,4], [1,2,6,5], [2,3,7,6], [3,0,4,7],
[4,5,6,7], [0,3,2,1]
)
Fc = len(cubeFaces)
cube = { 'baseFaceIndices' : cubeFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fc,0,2] }
surfacesDictionary['cubic'] = cube
return surfacesDictionary
_base_surfaces = _get_cubic_surface_dictionary.__func__()
_default_base = 'cubic' # default assignment is indicated in __init Docstring.
@classmethod
def _Xfev(cls,rez, basetype=None) :
"""
Calculates the number of faces, edges, and vertices of a CubicSurface.
Parameters
----------
rez : integer
Number of recursive subdivisions of a triangulated base faces.
basetype : None
Placeholder argument has no effect.
Returns
-------
list: the number of faces, edges, and vertices.
"""
F = 6*(4**rez)
return [ F, 2*F, F+2 ]
def __init__(self, rez=0, name=None, **kwargs):
"""
Create a cubic surface of 2 x 2 x 2 units.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of the rectangulated base faces.
Rez values range from 0 to 7.
name : string, optional, default: None.
Descriptive identifier for the geometry.
Raises
------
ValueError
If rez is not an integer in range 0 to 7.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
"""
# ................................................................
def midVFun(vectA, vectB) :
"""Mid-point between two points in the xy plane."""
mid = np.add(vectA,vectB)
mid = np.multiply(0.5,mid)
return mid
# ................................................................
if not isinstance(rez, int) :
raise ValueError('Incorrect rez type, must ba of type int')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be an int, 0 <= rez <= {}'.format(rez,_MAXREZ))
basetype = self._default_base
baseSurfObj = self._base_surfaces[basetype]
baseVcoor = baseSurfObj['baseVerticesCoor']
baseFaceVertexIndices =baseSurfObj['baseFaceIndices']
indexObj,vertexCoor = self.rectangulateBase(rez,baseVcoor,baseFaceVertexIndices, midVFun)
faceIndices = indexObj['face']
super().__init__(vertexCoor, faceIndices, name, **kwargs)
self.coorType = _COORSYS["XYZ"]
self._normal_scale = _normLength( len(faceIndices), 0.61, 1.50 )
self._rez = rez
self._basetype = basetype
return
[docs] def domain(self, xlim, ylim=None, zlim=None) :
"""
Set the domain of the cubic surface.
Used for setting the min,max values of the base cube.
Note: this CubicSurface method uses an alternative list of arguments
from the base class domain method. Using the xlim arg abbreviates setting
identical domains for all three axes.
Parameters
----------
xlim : a number, list or array, default: 1
The domain of the function evaluation. For a number, n,
the x,y,z axes domains will be [-n,n]. For a 1-dimensional 2-element list,
[a,b] will be assigned the domain for all 3 axes. Using a list of list (array),
as [ [a,b],[c,d],[e,f] ], the domain is assigned individually for each of the
three coordinate axes.
Must be a 2D array if ylim or zlim is set.
ylim, zlim : 2D array, default: None
Set the ylim and/or zlim.
Returns
-------
self : CubicSurface object
"""
# the following allows alternative 'arg list' used in super.
test = ylim is None and zlim is None
if not test:
super().domain(xlim,ylim,zlim)
return self
domain = Surface3DCollection._interpret_domain(xlim)
scl = lambda min,max : (max-min)/2
trn = lambda min,max : (max+min)/2
scale = [ scl(d[0],d[1]) for d in domain ]
translate = [ trn(d[0],d[1]) for d in domain ]
self.transform(None,scale,translate)
return self
[docs]class PlanarSurface(Surface3DCollection) :
"""
Flat square 3D surface in Cartesian coordinates.
Methods are inherited from the Surface3DCollection.
Planar surface geometries are created in this subclass.
"""
[docs] @staticmethod
def rand(rez=0,seed=None,name=None,kind=None,**kargs) :
"""
PlanarSurface with random triangular faces.
Parameters
----------
rez : number, optional, default: 0
Number of 'recursive' subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
seed : integer, optional, default: None
An initialization seed for random generation.
name : string, optional, default: None.
Descriptive identifier for the geometry.
kind : string, optional, default: None
(reserved, not implemented)
Raises
------
ValueError
If rez is not an integer in range 0 to 7.
If basetype is not recognized for this surface constructor.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
Returns
-------
surface : PlanarSurface object.
"""
if not isinstance(rez, (int,float)) :
raise ValueError('Incorrect rez type, must be a number')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be an int, 0 <= rez <= {}'.format(rez,_MAXREZ))
if seed is not None :
if not isinstance(seed, int) :
raise ValueError('Incorrect seed type, must ba of type int')
np.random.seed(seed)
N = int( 8*(4**rez)/2 + 4*(2**rez) + 1 )
x = 2*np.random.rand( N ) - 1
y = 2*np.random.rand( N ) - 1
xydata = np.array([x,y]).T
tess = spatial.Delaunay(xydata)
if name is None : name = 'random mesh'
surface = PlanarSurface(name=name,**kargs) # set to the class default object
# reset verts, faceIndices and edges (retain class coordinate system).
verts = np.zeros( (N,3))
verts[:,:2] = tess.points
faceIndices = tess.simplices
surface.vertexCoor = verts
surface.fvIndices = faceIndices
surface.evIndices = surface._get_edges_from_faces(faceIndices)
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( verts[faceIndices] )
surface._set_geometric_bounds()
surface._rez = str(rez) if isinstance(rez,int) else '{:.1f}'.format(rez)
surface._basetype = 'rand_'+str(seed)
return surface
@staticmethod
def _get_planar_surfaces_dictionary() :
""" Constructs the dictionary of base planar surface geometries. """
surfacesDictionary = {}
baseVert = [ [0,0,0],
[1,1,0], [-1,1,0], [-1,-1,0], [1,-1,0],
[1,0,0], [0,1,0], [-1,0,0], [0,-1,0] ]
qVert = [ baseVert[i] for i in range(0,5) ]
quadFaces=(
[0,1,2],[0,2,3],[0,3,4],[0,4,1] )
octFaces1=(
[0,5,6],[0,6,7],[0,7,8],[0,8,5],
[5,1,6],[6,2,7],[7,3,8],[8,4,5] )
octFaces2=(
[0,5,1],[0,1,6],[0,6,2],[0,2,7],
[0,7,3],[0,3,8],[0,8,4],[0,4,5] )
gridFaces=(
[0,5,1,6], [0,6,2,7], [0,7,3,8], [0,8,4,5]
)
Fq = len(quadFaces)
F1 = len(octFaces1)
F2 = len(octFaces2)
Fc = len(gridFaces)
quad = { 'baseFaceIndices' : quadFaces , 'baseVerticesCoor' : qVert, 'fevParam' : [Fq,2,1] }
oct1 = { 'baseFaceIndices' : octFaces1 , 'baseVerticesCoor' : baseVert, 'fevParam' : [F1,4,1] }
oct2 = { 'baseFaceIndices' : octFaces2 , 'baseVerticesCoor' : baseVert, 'fevParam' : [F2,4,1] }
grid = { 'baseFaceIndices' : gridFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fc,0,2] }
surfacesDictionary['quad'] = quad
surfacesDictionary['oct1'] = oct1
surfacesDictionary['oct2'] = oct2
surfacesDictionary['squ'] = grid
# ----------------------------------------------------------
mr = 0.01
MR = mr*np.sqrt(2)
qVertS=copy.copy( surfacesDictionary['quad']['baseVerticesCoor'] )
qVertS[0] = [ MR, 0, 0]
qVertS.append([ 0, MR, 0])
qVertS.append([-MR, 0, 0])
qVertS.append([0 ,-MR, 0])
baseVertS = copy.copy(surfacesDictionary['oct1']['baseVerticesCoor'])
baseVertS[0] = [ mr, mr, 0]
baseVertS.append([-mr, mr, 0])
baseVertS.append([-mr,-mr, 0])
baseVertS.append([ mr,-mr, 0])
quadFacesS=(
[5,1,2],[6,2,3],[7,3,4],[0,4,1],
[1,5,0],[2,6,5],[3,7,6],[4,0,7] )
octFaces1S=(
[ 0, 5, 6],[ 9, 6, 7],[10, 7, 8],[11, 8, 5],
[ 5, 1, 6],[ 6, 2, 7],[ 7, 3, 8],[ 8, 4, 5],
[ 5, 0,11],[ 6, 9, 0],[ 7,10, 9],[ 8,11,10] )
octFaces2S=(
[ 0, 5, 1],[ 0, 1, 6],[ 9, 6, 2],[ 9, 2, 7],
[10, 7, 3],[10, 3, 8],[11, 8, 4],[11, 4, 5],
[ 5, 0,11],[ 6, 9, 0],[ 7,10, 9],[ 8,11,10] )
FqS = len(quadFacesS)
F1S = len(octFaces1S)
F2S = len(octFaces2S)
quadS = { 'baseFaceIndices' : quadFacesS , 'baseVerticesCoor' : qVertS, 'fevParam' : [FqS,2,1] }
oct1S = { 'baseFaceIndices' : octFaces1S , 'baseVerticesCoor' : baseVertS, 'fevParam' : [F1S,4,1] }
oct2S = { 'baseFaceIndices' : octFaces2S , 'baseVerticesCoor' : baseVertS, 'fevParam' : [F2S,4,1] }
surfacesDictionary['quad_c'] = quadS
surfacesDictionary['oct1_c'] = oct1S
surfacesDictionary['oct2_c'] = oct2S
return surfacesDictionary
@staticmethod
def _midVectorFun(vectA, vectB) :
"""Mid-point between two points in the xy plane."""
mid = np.add(vectA,vectB)
mid = np.multiply(0.5,mid)
return mid
_base_surfaces = _get_planar_surfaces_dictionary.__func__()
_default_base = 'quad' # default assignment is indicated in __init Docstring.
[docs] @classmethod
def grid(cls,xdiv, ydiv, basetype='r', minrad=None, name=None, **kwargs) :
"""
PlanarSurface with X and Y rectangular faces.
Parameters
----------
xdiv : integer, required.
Number of X divisions.
ydiv : integer, required.
Number of Y divisions.
basetype : {'d','r'}, optional, default: 'r'
Starting surface geometries from which the surface is
constructed indicating face shape.
minrad : placeholder, not used.
name : string, optional, default: None.
Descriptive identifier.
Raises
------
ValueError
If xdiv is not an integer greater or equal to 1.
If ydiv is not an integer greater or equal to 1.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
Returns
-------
surface : PlanarSurface object.
"""
splitSize = None
minrad = None
# note: reversal of nang, nrad argsuments compared to other inherited class grid methods.
nrad, nang = ydiv, xdiv
surface = cls._square_net(nrad, nang, basetype,minrad, splitSize, name, **kwargs)
surface._postProc_surfaceColors(True)
surface._rez = surface._calc_rez( 4 )
return surface
[docs] @classmethod
def pntsurf(cls, xyz_points, name=None, **kargs) :
"""
PlanarSurface from points in Cartesian coordinates.
Parameters
----------
xyz_points : array, required.
An N X 3 array of N Cartesian coordinate points.
name : string, optional, default: None.
Descriptive identifier for the point surface.
Returns
-------
surface : PlanarSurface object
"""
dataPts = np.array(xyz_points)
xyzPts = dataPts.T
fvIndices = tri.Triangulation(xyzPts[0],xyzPts[1]).triangles
base_surface = Surface3DCollection(dataPts, fvIndices, color='tan').shade()
faceIndices = base_surface.fvIndices
verts = dataPts
if name is None : name = 'point surface'
surface = PlanarSurface(name=name,**kargs)
surface.vertexCoor = verts
surface.fvIndices = faceIndices
surface.evIndices = surface._get_edges_from_faces(faceIndices)
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( verts[faceIndices] )
surface._set_geometric_bounds()
surface._rez = surface._calc_rez(4)
surface._basetype = 'pntsurf'
return surface
@classmethod
def _Xfev(cls,rez, basetype=None) :
"""
Calculates the number of faces, edges, and vertices of a PlanarSurface.
Parameters
----------
rez : integer
Number of recursive subdivisions of a triangulated base faces.
basetype : string
Starting surface geometries from which the surface would be
constructed using recursive subdivisions of the triangular faces.
Returns
-------
list: the number of faces, edges, and vertices.
"""
return super()._fev(rez,basetype,cls)
[docs] @staticmethod
def meshgrid(X,Y,Z,revnormal=False,grid=False, name=None) :
"""
PlanarSurface object based on a meshgrid.
Parameters
----------
X,Y,Z : N x M arrays
Cartesian x,y,z coordinate values.
revnormal : { True, False }, boolean, default: False
Reverse the face normal directions if True.
grid : { True, False }, boolean, default: False
If True, polygon faces are rectangular. If
False, polygon faces are triangular.
name : string, optional, default: None.
Descriptive identifier for the mesh surface.
Returns
-------
surface : PlanarSurface object
"""
vertCoor = np.reshape(np.stack((X,Y,Z), axis=2 ), (-1,3) )
lenY, lenX = Z.shape
mtr = np.reshape( np.arange(lenX*lenY) , (lenY,lenX))
A = mtr[:,0:lenX-1][0:lenY-1,:].flatten()
B = A + 1
C = A + lenX +1
D = A + lenX
if revnormal :
t = B
B = D
D = t
faceIndices = np.stack( (A,B,C,D) ).T
if grid : obj = PlanarSurface( basetype='squ')
else : obj = PlanarSurface( )
obj.vertexCoor = vertCoor
obj.fvIndices = faceIndices
obj._normal_scale = _normLength( len(faceIndices), 0, 0, 4 )
obj._set_geometric_bounds()
obj._surfaceShapeChanged = True
v = obj.vertexCoor
fvi = obj.fvIndices
if grid :
obj.vfIndicesList = obj._vertexCommonFaceIndices()
edgeIndices = obj._get_edges_from_faces(obj.fvIndices)
obj.evIndices = np.array(edgeIndices)
obj.efIndicesList = obj._edgeCommonFaceIndices()
obj.set_verts(v[fvi])
if not grid : obj.triangulate()
if name is None : name = 'meshgrid'
obj.name = name
return obj
def __init__(self, rez=0, basetype=None, name=None, **kwargs):
"""
Create flat square surface of 2 x 2 units.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
basetype : {'quad','oct1','oct2','squ'}, optional, default: 'quad'
Starting surface geometries from which the surface is
constructed using recursive subdivisions of the faces.
name : string, optional, default: None.
Descriptive identifier for the geometry.
Raises
------
ValueError
If rez is not an integer in range 0 to 7.
If basetype is not recognized for this surface constructor.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
"""
if not isinstance(rez, int) :
raise ValueError('Incorrect rez type, must ba of type int')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be an int, 0 <= rez <= {}'.format(rez,_MAXREZ))
if basetype is None : basetype = self._default_base
if basetype not in self._base_surfaces :
raise ValueError('Basetype {} is not recognized. Possible values are: {}'.format(basetype,list(self._base_surfaces)))
baseSurfObj = self._base_surfaces[basetype]
baseVcoor = baseSurfObj['baseVerticesCoor']
baseFaceVertexIndices =baseSurfObj['baseFaceIndices']
if basetype == 'squ' :
indexObj,vertexCoor = self.rectangulateBase(rez,baseVcoor,baseFaceVertexIndices, self._midVectorFun)
else :
indexObj,vertexCoor = self.triangulateBase(rez,baseVcoor,baseFaceVertexIndices, self._midVectorFun)
faceIndices = indexObj['face']
super().__init__(vertexCoor, faceIndices, name, **kwargs)
self.coorType = _COORSYS["XYZ"]
self._normal_scale = _normLength( len(faceIndices), 0, 0, 4 )
self._rez = rez
self._basetype = basetype
return
[docs] def scale_dataframe(self,X,Y,Z) :
"""
Scaling the surface geometry based on a datagrid.
Parameters
----------
X, Y, Z : N x M arrays
Minimum and maximum values of the arrays are used to scale
and translate the surface from an intial domain of
[ (-1,1), (-1,1), (0,1) ]
Returns
-------
self : PlanarSurface object
"""
return _scale_dataframe(self,X,Y,Z)
[docs] def domain(self, xlim=None, ylim=None, zcoor=None) :
"""
Set the domain of the planar surface.
Used for setting the limits of the base plane.
Parameters
----------
xlim, ylim : 2d arrays
Minimum and maximum values of the arrays are used to scale
and translate the surface from an intial domains of
[ -1, 1 ]
If set to None, the domains are unchanged.
zcoor: number, default : 0
Z offset of the plane.
If set to None, the Z offset is unchanged.
Returns
-------
self : PlanarSurface object
"""
def setDomain(xyz, mx,bx,my,by,mz) :
x,y,z = xyz
X = mx*x + bx
Y = my*y + by
Z = z if mz is None else np.full(len(z),mz)
return X,Y,Z
def getCont( lims ) :
m, b = 1.0, 0.0
if lims is not None :
# allow single bndry to be input as a float
if isinstance(lims, (int,float)) : lims = [-lims,lims]
m, b = 0.5*(lims[1]-lims[0]) , 0.5*(lims[1]+lims[0])
return m, b
mx, bx = getCont(xlim)
my, by = getCont(ylim)
self.map_geom_from_op(lambda A : setDomain(A, mx,bx,my,by,zcoor))
return self
[docs]class PolarSurface(Surface3DCollection) :
"""
Flat disk 3D surface in polar coordinates.
Methods are inherited from the Surface3DCollection.
Polar surface geometries are created in this subclass.
"""
[docs] @staticmethod
def rand(rez=0,seed=None,name=None,kind=None,**kargs) :
"""
PolarSurface with random triangular faces.
Parameters
----------
rez : number, optional, default: 0
Number of 'recursive' subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
seed : integer, optional, default: None
An initialization seed for random generation.
name : string, optional, default: None.
Descriptive identifier for the geometry.
kind : string, optional, default: None
(reserved, not implemented)
Raises
------
ValueError
If rez is not an integer in range 0 to 7.
If basetype is not recognized for this surface constructor.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
Returns
-------
surface : PolarSurface object.
"""
#.........................................................
def random_linDist(size) :
nLoop,nControl,samples = 0, 10**6, []
while len(samples)<size and nLoop<nControl:
x=np.random.random_sample()
prob=x/2
assert prob>=0 and prob<=1
if np.random.random_sample() <=prob: samples += [x]
nLoop+=1
return np.array(samples)
#.........................................................
if not isinstance(rez, (int,float)) :
raise ValueError('Incorrect rez type, must be a number.')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be a number, 0 <= rez <= {}'.format(rez,_MAXREZ))
if seed is not None :
if not isinstance(seed, int) :
raise ValueError('Incorrect seed type, must ba of type int')
np.random.seed(seed)
N = int( 6*(4**rez)/2 + 3*(2**rez) + 1 )
r = random_linDist( N )
t = 2*np.pi*np.random.rand( N )
z = np.zeros( N )
data = PolarSurface.coor_convert([r,t,z],True).T
rtdata = data[:,:2]
tess = spatial.Delaunay(rtdata)
if name is None : name = 'random mesh'
surface = PolarSurface(name=name,**kargs) # set to the class default object
# reset verts, faceIndices and edges (retain class coordinate system).
verts = np.zeros( (N,3))
verts[:,:2] = tess.points
faceIndices = tess.simplices
surface.vertexCoor = verts
surface.fvIndices = faceIndices
surface.evIndices = surface._get_edges_from_faces(faceIndices)
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( verts[faceIndices] )
surface._set_geometric_bounds()
surface._rez = str(rez) if isinstance(rez,int) else '{:.1f}'.format(rez)
surface._basetype = 'rand_'+str(seed)
return surface
@staticmethod
def _get_polar_surfaces_dictionary() :
"""Constructs the dictionary of base polar surface geometries."""
surfacesDictionary = {}
soff = 0.0001
x0 = 0.5
y0 = np.sqrt(3)/2
baseVert = [ [0,0,0],
[1,0,0], [x0,y0,0], [-x0,y0,0],
[-1,0,0], [-x0,-y0,0], [x0,-y0,0] ]
qVert = [ [0,0,0],
[1,0,0], [0,1,0], [-1,0,0], [0,-1,0] ]
x0S, y0S = soff*y0, soff*x0
baseVertS = [ [1,-y0S,0], [1, y0S,0],
[x0,y0,0], [-x0,y0,0], [-1,0,0],
[-x0,-y0,0], [x0,-y0,0], [ x0S, y0S,0],
[ 0, soff,0], [-x0S, y0S,0],
[-x0S, -y0S,0], [ 0,-soff,0], [ x0S, -y0S,0] ]
qVertS = [ [1,-soff,0],
[1, soff,0], [0,1,0], [-1,0,0], [0,-1,0],
[ soff, soff,0], [-soff ,soff,0], [-soff,-soff,0], [soff,-soff,0] ]
baseFaces=([0,1,2],[0,2,3],[0,3,4],[0,4,5],[0,5,6],[0,6,1])
quadFaces=([0,1,2],[0,2,3],[0,3,4],[0,4,1])
baseFacesS=( [7,1,2], [8,2,3], [9,3,4], [10,4,5], [11,5,6], [12,6,0],
[2,8,7], [3,9,8], [4,10,9], [5,11,10], [6,12,11] )
quadFacesS=( [5,1,2],[6,2,3],[7,3,4],[8,4,0],
[2,6,5],[3,7,6],[4,8,7] )
Fs = len(quadFaces)
Fh = len(baseFaces)
FsS = len(quadFacesS)
FhS = len(baseFacesS)
square = { 'baseFaceIndices' : quadFaces , 'baseVerticesCoor' : qVert, 'fevParam' : [Fs,2,1] }
hexag = { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fh,3,1] }
squareS = { 'baseFaceIndices' : quadFacesS , 'baseVerticesCoor' : qVertS, 'fevParam' : [FsS,4.5,1] }
hexagS = { 'baseFaceIndices' : baseFacesS , 'baseVerticesCoor' : baseVertS, 'fevParam' : [FhS,6.5,1] }
surfacesDictionary['squ'] = square
surfacesDictionary['hex'] = hexag
surfacesDictionary['squ_s'] = squareS
surfacesDictionary['hex_s'] = hexagS
# ... 'special case' only for fev method for polar surfaces.
surfacesDictionary['squ_c'] = { 'fevParam' : [16,5,1] }
surfacesDictionary['hex_c'] = { 'fevParam' : [12,4,1] }
return surfacesDictionary
@staticmethod
def _midVectorFun(vectA, vectB) :
"""Mid-point at the mid radial position between two points in a plane."""
mid = np.add(vectA,vectB)
mid = np.multiply(0.5,mid)
a = (np.linalg.norm(vectA)+np.linalg.norm(vectB))/2
b = np.linalg.norm(mid)
mid = np.multiply(mid,a/b)
return mid
@staticmethod
def _map3Dto2Dsquare(xyz) :
"""Override map surface to rectagular unit coordinates (0,1)."""
x,y,_ = xyz
theta = np.arctan2(y,x)
theta = np.where( theta<0, 2.0*np.pi + theta ,theta)
a = theta/(2*np.pi)
r = np.sqrt(x*x+y*y)
b = 1-r
a = np.clip(a,0.0,1.0)
b = np.clip(b,0.0,1.0)
return [a,b]
[docs] @staticmethod
def coor_convert(xyz, tocart=False) :
"""
Tranformation between polar and Cartesian coordinates.
The angular coordinate is in radians with domain 0 to 2pi.
Parameters
----------
xyz : 3xN array
N number of vectors in either polar or cartesian coordinates.
tocart : { True, False }, default: False
If True, input is polar and output is cartesian.
If False, input is cartesian and output is polar.
Returns
-------
list: 3xN array
Transformed coordinates.
"""
if tocart :
r, theta, z = xyz
x = r * np.cos(theta)
y = r * np.sin(theta)
abc = [x,y,z]
else:
x,y,z = xyz
R = np.sqrt( x*x + y*y )
theta = np.arctan2(y,x)
theta = np.where(theta<0,theta+2*np.pi,theta)
abc = [R,theta,z]
abc = np.array(abc)
return abc
_base_surfaces = _get_polar_surfaces_dictionary.__func__()
_default_base = 'hex' # default assignment is indicated in __init Docstring.
_geomType = _COORSYS['POLAR']
@classmethod
def _uvw_to_xyz( cls,uvw, rtz) :
"""Override rotational transformation at a polar coordinate."""
return super()._disp_to_xyz( uvw, rtz[1] )
[docs] @classmethod
def grid(cls,nrad, nang, basetype='d', minrad=None, name=None, **kwargs) :
"""
PolarSurface with axial and radial faces.
Parameters
----------
nrad : integer, required.
Number of radial subdivisions. Minimum value is 1
nang : integer, required.
Number of angular subdivisions. Minimum value is 3.
basetype : {'d','s','q','w','r','x'}, optional, default: 'd'
Starting surface geometries from which the surface is
constructed indicating face shape and if split.
minrad : scalar, optional, default: 0.01
The minimum distance from the z-axis of any vertex coordinate.
Use is dependent on basetype.
name : string, optional, default: None.
Descriptive identifier.
Raises
------
ValueError
If nrad is not an integer greater or equal to 1.
If nang is not an integer greater or equal to 3.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
Returns
-------
surface : PolarSurface object.
"""
Nfaces = { 'd': 3, 's': 3, 'q':6, 'w':6, 'r':3, 'x':3 } # 1,3
splitSize = None
surface = cls._square_net(nrad, nang, basetype,minrad, splitSize, name, **kwargs)
surface._postProc_surfaceColors(True)
surface._rez = surface._calc_rez( 6 )
return surface
[docs] @classmethod
def pntsurf(cls, rtz_points, name=None, **kargs) :
"""
PolarSurface from points in polar coordinates.
Parameters
----------
rtz_points : array, required.
An N X 3 array of N polar coordinate points.
name : string, optional, default: None.
Descriptive identifier for the point surface.
Returns
-------
surface : PolarSurface object
"""
dataPts = np.array(rtz_points)
xyzPts = PolarSurface.coor_convert(dataPts.T,True)
fvIndices = tri.Triangulation(xyzPts[0],xyzPts[1]).triangles
base_surface = Surface3DCollection(dataPts, fvIndices, color='tan').shade()
faceIndices = base_surface.fvIndices
verts = xyzPts.T
if name is None : name = 'point surface'
surface = PolarSurface(name=name,**kargs)
surface.vertexCoor = verts
surface.fvIndices = faceIndices
surface.evIndices = surface._get_edges_from_faces(faceIndices)
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( verts[faceIndices] )
surface._set_geometric_bounds()
surface._rez = surface._calc_rez(4)
surface._basetype = 'pntsurf'
return surface
@classmethod
def _Xfev(cls,rez, basetype=None, minrad=None) :
"""
Calculates number of faces, edges, and vertices of a PolarSurface.
Parameters
----------
rez : integer
Number of recursive subdivisions of a triangulated base faces.
basetype : string
Starting surface geometries from which the surface would be
constructed using recursive subdivisions of the triangular faces.
Returns
-------
list: the number of faces, edges, and vertices.
"""
if isinstance(rez,(list,tuple)) : # check for disk...
nrad = rez[0]
nang = rez[1]
# only check is ending in s
isSplit = basetype[ len(basetype)-1] == 's'
if minrad is None :
f = 2*nrad*nang - nang
e = 3*nrad*nang - nang
v = nrad*nang + 1
if isSplit :
e += nrad
v += nrad
else : # any value will trigger this condition.
f = 2*nrad*nang
e = 3*nrad*nang + nang
v = (nrad+1)*nang
if isSplit :
e += nrad
v += nrad + 1
return [f,e,v]
return super()._fev(rez,basetype,cls)
@staticmethod
def _flat_split_cyclinder(rez,basetype,minrad, **kwargs) :
"""
Transform cylindrical base surface to polar coordinates.
note : This method provides an alternative to the
poler basetypes of 'hex_s' and 'squ_s'.
The grids constructed here may provide smoother
geometries depending on the mapping functions.
"""
# .................................................
def cyl2pol(rtz):
r,t,z = rtz
b = (1.0+minrad)/2.0
m = -(1.0-minrad)/2.0
R = m*z + b
Z = 0*z
return R,t,Z
# .................................................
if basetype == 'squ_c' : cylbase = 'squ_s'
if basetype == 'hex_c' : cylbase = 'tri_s'
cyl_obj = CylindricalSurface(rez,cylbase, **kwargs).map_geom_from_op(cyl2pol)
baseFaceVertexIndices = cyl_obj._base_surfaces[cylbase] ['baseFaceIndices']
fvIndices = cyl_obj.fvIndices
vertexCoor = cyl_obj.vertexCoor
evIndices = cyl_obj.evIndices
return vertexCoor, fvIndices, evIndices, baseFaceVertexIndices,
def __init__(self, rez=0, basetype=None, minrad=None, name=None, **kwargs):
"""
Create flat disk surface of unit radius.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
basetype : {'squ','hex','squ_s','hex_s','squ_c','hex_c'}, optional, default: 'hex'
Starting surface geometries from which the surface is
constructed using recursive subdivisions of the triangular faces.
minrad : scalar, optional, default: 0.01
For basetypes 'squ_c' and 'hex_c, the minimum
distance from the z-axis of any vertex coordinate.
name : string, optional, default: None.
Descriptive identifier for the geometry.
Raises
------
ValueError
If rez is not an integer in range 0 to 7.
If basetype is not recognized for this surface constructor.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
"""
if minrad is None : minrad = _MINRAD
if not isinstance(rez, int) :
raise ValueError('Incorrect rez type, must ba of type int')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be an int, 0 <= rez <= {}'.format(rez,_MAXREZ))
if basetype is None : basetype = self._default_base
# ... check for 'special case' (ugly code but needed)
if basetype == 'squ_c' or basetype == 'hex_c' :
vertexCoor, faceIndices, edgeIndices, bfvi = self._flat_split_cyclinder(rez,basetype,minrad, **kwargs)
baseFaceVertexIndices = bfvi
else :
if basetype not in self._base_surfaces :
raise ValueError('Basetype {} is not recognized. Possible values are: {}'.format(basetype,list(self._base_surfaces)))
baseSurfObj = self._base_surfaces[basetype]
baseVcoor = baseSurfObj['baseVerticesCoor']
baseFaceVertexIndices =baseSurfObj['baseFaceIndices']
indexObj,vertexCoor = self.triangulateBase(rez,baseVcoor,baseFaceVertexIndices, self._midVectorFun)
faceIndices = indexObj['face']
edgeIndices = indexObj['edge']
super().__init__(vertexCoor, faceIndices, name, **kwargs)
self.coorType = _COORSYS["POLAR"]
self._normal_scale = _normLength( len(faceIndices), 1.2, 1.55, np.pi )
self._rez = rez
self._basetype = basetype
return
[docs] def domain(self, radius=None, zcoor=None) :
"""
Set the domain of the polar surface.
Used for setting the limits of the base polar plane.
Parameters
----------
radius: number, default : 1
Radius of the polar surface.
If set to None, the radius is unchanged.
zcoor: number, default : 0
Z offset of the polar plane.
If set to None, the Z offset is unchanged.
Returns
-------
self : PolarSurface object
"""
def setDomain(rtz, mr, mz) :
r,t,z = rtz
R = mr*r
Z = z if mz is None else np.full(len(z),mz)
return R,t,Z
def getCont(lims) :
m = 1.0
if lims is not None : m = lims
return m
mr = getCont(radius)
self.map_geom_from_op(lambda A : setDomain(A, mr,zcoor ))
return self
[docs]class CylindricalSurface(Surface3DCollection) :
"""
Cylindrical 3D surface in cylindrical coordinates.
Methods are inherited from the Surface3DCollection.
Cylindrical surface geometries are created in this subclass.
"""
[docs] @staticmethod
def rand(rez=0,seed=None,name=None,kind=None,**kargs) :
"""
CylindricalSurface with random triangular faces.
Parameters
----------
rez : number, optional, default: 0
Number of 'recursive' subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
seed : integer, optional, default: None
An initialization seed for random generation.
name : string, optional, default: None.
Descriptive identifier for the geometry.
kind : string, optional, default: None
(reserved, not implemented)
Raises
------
ValueError
If rez is not an integer in range 0 to 7.
If basetype is not recognized for this surface constructor.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
Returns
-------
surface : CylindricalSurface object.
"""
#.........................................................
def random_sinDist(size) :
nLoop,nControl,samples = 0, 10**6, []
while len(samples)<size and nLoop<nControl:
x=np.random.random_sample()
prob=np.sin(x*np.pi)/2
assert prob>=0 and prob<=1
if np.random.random_sample() <=prob: samples += [x]
nLoop+=1
return np.array(samples)
#.........................................................
if not isinstance(rez, (int, float)) :
raise ValueError('Incorrect rez type, must be a number.')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be a number, 0 <= rez <= {}'.format(rez,_MAXREZ))
if seed is not None :
if not isinstance(seed, int) :
raise ValueError('Incorrect seed type, must ba of type int')
np.random.seed(seed)
N = int(12*(4**rez)/2 + 3*(2**rez))
r = np.ones( N )
t = np.random.uniform(0,2*np.pi, N)
z = np.random.uniform(-1,1, N)
dataPts = np.array( [r,t,z] ).T
if name is None : name = 'random mesh'
surface = CylindricalSurface.pntsurf(dataPts, name=None, **kargs)
surface._rez = str(rez) if isinstance(rez,int) else '{:.1f}'.format(rez)
surface._basetype = 'rand_'+str(seed)
return surface
@staticmethod
def _get_cylindrical_surfaces_dictionary() :
"""Constructs the dictionary of base cylindrical surface geometries. """
surfacesDictionary = {}
def construct_tricylinder() :
y0 = np.sqrt(3.0)/2.0
x0 = 0.5
baseVert = [ [1,0,0], [-x0,-y0,0], [-x0,y0,0],
[-1,0,1], [x0,-y0,1], [x0,y0,1],
[-1,0,-1], [x0,-y0,-1], [x0,y0,-1] ]
baseFaces=(
[0,2,5],[2,1,3],[1,0,4],
[0,5,4],[2,3,5],[1,4,3],
[0,8,2],[2,6,1],[1,7,0],
[0,7,8],[2,8,6],[1,6,7])
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,3,0] }
def construct_tricylinder_v() :
vfi_d = construct_tricylinder()
baseVert = vfi_d['baseVerticesCoor']
baseFaces = vfi_d['baseFaceIndices']
baseVert = baseVert + [[0,0,1],[0,0,-1]]
baseFaces = list(baseFaces)
addFaces = [
[3,4,9],[4,5,9],[5,3,9],
[6,8,10],[8,7,10],[7,6,10] ]
baseFaces = baseFaces + addFaces
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
def construct_tri2cylinder() :
y0 = np.sqrt(3.0)/2.0
x0 = 0.5
baseVert = [ [-1,0,0], [x0,-y0,0], [x0,y0,0],
[1,0,1], [-x0,-y0,1], [-x0,y0,1],
[1,0,-1], [-x0,-y0,-1], [-x0,y0,-1] ]
baseFaces=(
[0,4,5],[1,3,4],[2,5,3],
[0,8,7],[1,7,6],[2,6,8],
[0,7,4],[1,6,3],[2,8,5],
[0,5,8],[1,4,7],[2,3,6])
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,3,0] }
def construct_tri2cylinder_v() :
vfi_d = construct_tri2cylinder()
baseVert = vfi_d['baseVerticesCoor']
baseFaces = vfi_d['baseFaceIndices']
baseVert = baseVert + [[0,0,1],[0,0,-1]]
baseFaces = list(baseFaces)
addFaces = [
[3,9,4],[4,9,5],[5,9,3],
[6,10,8],[7,10,6],[8,10,7] ]
baseFaces = baseFaces + addFaces
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
def construct_tri2cylinder_sliced() :
y0 = np.sqrt(3.0)/2.0
x0 = 0.5
yof = 0.0001
baseVert = [ [-1,0,0], [x0,-y0,0], [x0,y0,0],
[1,yof,1], [-x0,-y0,1], [-x0,y0,1],
[1,yof,-1], [-x0,-y0,-1], [-x0,y0,-1],
[1,-yof,1], [1,-yof,-1] ]
baseFaces=(
[0,4,5],[1,9,4],[2,5,3],
[0,8,7],[1,7,10],[2,6,8],
[0,7,4],[1,10,9],[2,8,5],
[0,5,8],[1,4,7],[2,3,6])
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,4,1] }
def construct_tri2cylinder_sliced_v() :
vfi_d = construct_tri2cylinder_sliced()
baseVert = vfi_d['baseVerticesCoor']
baseFaces = vfi_d['baseFaceIndices']
baseVert = baseVert + [[0,0,1],[0,0,-1]]
baseFaces = list(baseFaces)
addFaces = [
[3,5,11],[5,4,11],[4,9,11],
[8,6,12],[7,8,12],[10,7,12] ]
baseFaces = baseFaces + addFaces
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
def construct_quadcylinder() :
xy0 = 1.0/np.sqrt(2.0)
baseVert = [ [1,0,0], [0,-1,0], [-1,0,0], [0,1,0],
[xy0,xy0,1], [-xy0,xy0,1], [-xy0,-xy0,1], [xy0,-xy0,1],
[xy0,xy0,-1], [-xy0,xy0,-1], [-xy0,-xy0,-1], [xy0,-xy0,-1] ]
baseFaces=(
[1,0,7],[2,1,6],[3,2,5],[0,3,4],
[4,3,5],[5,2,6],[6,1,7],[7,0,4],
[0,8,3],[3,9,2],[2,10,1],[1,11,0],
[0,11,8],[3,8,9],[2,9,10],[1,10,11] )
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,4,0] }
def construct_quadcylinder_v() :
vfi_d = construct_quadcylinder()
baseVert = vfi_d['baseVerticesCoor']
baseFaces = vfi_d['baseFaceIndices']
baseVert = baseVert + [[0,0,1],[0,0,-1]]
baseFaces = list(baseFaces)
addFaces = [
[4,5,12],[5,6,12],[6,7,12],[7,4,12],
[9,8,13],[10,9,13],[11,10,13],[8,11,13] ]
baseFaces = baseFaces + addFaces
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
def construct_quad2cylinder() :
xy0 = 1.0/np.sqrt(2.0)
baseVert = [ [xy0,xy0,0], [-xy0,xy0,0], [-xy0,-xy0,0], [xy0,-xy0,0],
[1,0,1], [0,-1,1], [-1,0,1], [0,1,1],
[1,0,-1], [0,-1,-1], [-1,0,-1], [0,1,-1] ]
baseFaces=(
[0,7,4],[1,6,7],[2,5,6],[3,4,5],
[0,8,11],[1,11,10],[2,10,9],[3,9,8],
[0,4,8],[1,7,11],[2,6,10],[3,5,9],
[0,11,7],[1,10,6],[2,9,5],[3,8,4] )
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,4,0] }
def construct_quad2cylinder_v() :
vfi_d = construct_quad2cylinder()
baseVert = vfi_d['baseVerticesCoor']
baseFaces = vfi_d['baseFaceIndices']
baseVert = baseVert + [[0,0,1],[0,0,-1]]
baseFaces = list(baseFaces)
addFaces = [
[5,4,12],[6,5,12],[7,6,12],[4,7,12],
[8,9,13],[9,10,13],[10,11,13],[11,8,13] ]
baseFaces = baseFaces + addFaces
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
def construct_quad2cylinder_sliced() :
xy0 = 1.0/np.sqrt(2.0)
yof = 0.0001
baseVert = [ [xy0,xy0,0], [-xy0,xy0,0], [-xy0,-xy0,0], [xy0,-xy0,0],
[1,yof,1], [0,-1,1], [-1,0,1], [0,1,1],
[1,yof,-1], [0,-1,-1], [-1,0,-1], [0,1,-1],
[1,-yof,1], [1,-yof,-1] ]
baseFaces=(
[0,7,4],[1,6,7],[2,5,6],[3,12,5],
[0,8,11],[1,11,10],[2,10,9],[3,9,13],
[0,4,8],[1,7,11],[2,6,10],[3,5,9],
[0,11,7],[1,10,6],[2,9,5],[3,13,12] )
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,5,1] }
def construct_quad2cylinder_sliced_v() :
vfi_d = construct_quad2cylinder_sliced()
baseVert = vfi_d['baseVerticesCoor']
baseFaces = vfi_d['baseFaceIndices']
baseVert = baseVert + [[0,0,1],[0,0,-1]]
baseFaces = list(baseFaces)
addFaces = [
[4,7,14],[7,6,14],[6,5,14],[5,12,14],
[11,8,15],[10,11,15],[9,10,15],[13,9,15] ]
baseFaces = baseFaces + addFaces
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
surfacesDictionary['tri'] = construct_tricylinder()
surfacesDictionary['squ'] = construct_quadcylinder()
surfacesDictionary['tri2'] = construct_tri2cylinder()
surfacesDictionary['squ2'] = construct_quad2cylinder()
surfacesDictionary['tri_s'] = construct_tri2cylinder_sliced()
surfacesDictionary['squ_s'] = construct_quad2cylinder_sliced()
surfacesDictionary['tri_v'] = construct_tricylinder_v()
surfacesDictionary['squ_v'] = construct_quadcylinder_v()
surfacesDictionary['tri2_v'] = construct_tri2cylinder_v()
surfacesDictionary['squ2_v'] = construct_quad2cylinder_v()
surfacesDictionary['tri_sv'] = construct_tri2cylinder_sliced_v()
surfacesDictionary['squ_sv'] = construct_quad2cylinder_sliced_v()
return surfacesDictionary
@staticmethod
def _midVectorFun(vectA, vectB) :
"""Mid-point between two points on a unit cylinder."""
mid = np.add(vectA,vectB)
mid = np.multiply(0.5,mid)
z = np.dot(mid,[0,0,1])
v = np.subtract(mid,[0,0,z])
xy = v/np.linalg.norm(v)
# xy is unit vector at cylinder radius, but
# is average of radius for top & bottom faces.
zA = np.dot(vectA,[0,0,1])
vA = np.subtract(vectA,[0,0,zA])
zB = np.dot(vectB,[0,0,1])
vB = np.subtract(vectB,[0,0,zB])
a = ( np.linalg.norm(vA) + np.linalg.norm(vB) )/2.0
xy = a*xy
return np.add(xy,[0,0,z])
@staticmethod
def _map3Dto2Dsquare(xyz) :
"""Override map surface to rectagular unit coordinates (0,1)."""
x,y,z = xyz
theta = np.arctan2(y,x)
theta = np.where( theta<0, 2.0*np.pi + theta ,theta)
a = theta/(2*np.pi)
b = z/2 + 0.5
a = np.clip(a,0.0,1.0)
b = np.clip(b,0.0,1.0)
return [a,b]
[docs] @staticmethod
def coor_convert(xyz, tocart=False) :
"""
Transformation between cylindrical and Cartesian coordinates.
The angular coordinate is in radians with domain 0 to 2pi.
Parameters
----------
xyz : 3xN
N number of vectors in either cylindrical or cartesian coordinates,
including the z coordinate( which remains unchanged).
tocart : { True, False }, default: False
If True, input is cylindrical and output is cartesian.
If False, input is cartesian and output is cylindrical.
Returns
-------
list: 3xN array
Transformed coordinates.
"""
if tocart :
r, theta, z = xyz
x = r * np.cos(theta)
y = r * np.sin(theta)
abc = [x,y,z]
else:
x,y,z = xyz
R = np.sqrt( x*x + y*y )
theta = np.arctan2(y,x)
theta = np.where(theta<0,theta+2*np.pi,theta)
abc = [R,theta,z]
abc = np.array(abc)
return abc
_base_surfaces = _get_cylindrical_surfaces_dictionary.__func__()
_default_base = 'tri' # default assignment is indicated in __init Docstring.
_geomType = _COORSYS['CYLINDRICAL']
@classmethod
def _cyl_vol_square_net(cls, numZ, numU, basetype,minRad, splitSize, numR, name=None, **kwargs) :
"""
method only called by the CylinderSurface classmethod grid ONLY for volume cylinder objects.
"""
# DevNote: this method uses the same algorithm as the base class _square_net classmethod.
# This was done to simplify the code for this special case rather than
# add additional 'convoluted' logic to the _square_net method.
# For this case, the algorithm combines one cylindrical and two polar grids.
# Note: numX and num_X are the number of divisions and the number of coor positions.
if numR is None : numR = int((numZ+1)/2)
if basetype is None : basetype = 'd'
if minRad is None : minRad = _MINRAD
if splitSize is None : splitSize = _SPLITSIZE
basetype = basetype.lower()
baselist = ['d','s','q','w','r','x']
if basetype not in baselist :
raise ValueError("Grid basetype '{}' is not recognized. Possible values are: {}.".format(basetype,list(baselist)))
isSplit = (basetype == 's') or (basetype == 'w') or (basetype == 'x')
needsTriag = not ( (basetype == 'r') or (basetype == 'x') )
hasCenter = (basetype == 'd') or (basetype == 's')
# T angle ranges .....
num_U = numU+1 if isSplit else numU
min_U = 0.0
max_U = (1-1/num_U)*2*np.pi
if isSplit : max_U = (2-splitSize)*np.pi
coor_U = np.linspace(min_U,max_U,num_U)
# R disk ranges ........
if hasCenter :
num_R,inn_R = numR-1, 1.0/numR
out_R = 1.0-(1.0/numR)
else :
num_R,inn_R = numR, minRad
out_R = 1.0-(1.0/numR) + minRad/numR
coor_R_bot = np.linspace(inn_R,out_R, num_R)
coor_R_mid = np.full(numZ+1,1.0)
coor_R_top = np.linspace(out_R,inn_R, num_R)
num_Z = numZ+1
coor_Z_bot = np.full(num_R,-1.0)
coor_Z_mid = np.linspace(-1.0,1.0,num_Z)
coor_Z_top = np.full(num_R, 1.0)
coor_R = np.concatenate( (coor_R_bot,coor_R_mid,coor_R_top),axis=0)
coor_Z = np.concatenate( (coor_Z_bot,coor_Z_mid,coor_Z_top),axis=0)
numV = 2*numR + numZ
num_V = 2*num_R + num_Z
totVerts = num_U*num_V
U_coor = np.tile(coor_U,num_V)
R_coor = np.tile(coor_R,(num_U,1))
R_coor = np.ravel(R_coor,order='F')
Z_coor = np.tile(coor_Z,(num_U,1))
Z_coor = np.ravel(Z_coor,order='F')
abc = np.array([R_coor,U_coor,Z_coor])
xyz = CylindricalSurface.coor_convert(abc,True)
verts = np.array(xyz).T
numVt = num_V-1 if hasCenter else numV
faceIndices = cls._grid_face_indices(numU,numVt,isSplit,needsTriag)
if hasCenter :
# add top triangular faces ...
topVert = np.array([[0.0,0.0,1]])
topIndex = totVerts
startIndex = topIndex - num_U
endIndex = topIndex-2 if isSplit else topIndex-1
a = np.linspace(startIndex,endIndex,numU,dtype=int)
b = a+1 if isSplit else np.append( a[1:], [startIndex] )
c = np.full(numU,topIndex)
abc = np.array([a,b,c]).T
verts = np.append(verts,topVert, axis=0)
faceIndices = np.append(faceIndices,abc, axis=0)
# add bottom triangular facess ...
btmVert = np.array([[0.0,0.0,-1]])
btmIndex = totVerts+1
startIndex,endIndex = 0,numU-1
a = np.linspace(startIndex,endIndex,numU,dtype=int)
b = a+1 if isSplit else np.append( a[1:], [startIndex] )
c = np.full(numU,btmIndex)
bac = np.array([b,a,c]).T # bac instead for abc for outward normals.
verts = np.append(verts,btmVert, axis=0)
faceIndices = np.append(faceIndices,bac, axis=0)
surface = cls(**kwargs) # set to the class default object
# reset verts, faceIndices and edges (retain class coordinate system).
surface.vertexCoor = verts
surface.fvIndices = faceIndices
surface.evIndices = surface._get_edges_from_faces(faceIndices)
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( verts[faceIndices] )
surface._set_geometric_bounds()
surface.name = name
surface._basetype = 'grid_'+basetype+'v'
return surface
@classmethod
def _uvw_to_xyz(cls, uvw, rtz) :
"""Override rotational tranformation at a cylindrical coordinate."""
return super()._disp_to_xyz( uvw, rtz[1], None )
[docs] @classmethod
def grid(cls, nver ,nang, basetype='d',minrad=None, numR=None, name=None, **kwargs) :
"""
CylindricalSurface with axial and vertical faces.
Parameters
----------
nver : integer, required.
Number of vertical subdivisions. Minimum value is 1
nang : integer, required.
Number of angular subdivisions. Minimum value is 3.
basetype : {'d','q','r','s','w','x', 'dv','qv','rv','sv','wv','xv'}, optional, default: 'd'
Starting surface geometries from which the surface is
constructed indicating face shape and if split.
minrad : scalar, optional, default: 0.01
The minimum distance from the z-axis of any vertex coordinate.
Use is dependent on basetype.
numR : integer, optional, default: None
Only applicable for volume cylinder geometries.
Number of radial subdivisions for volume cylinders.
When None, default is half the number of vertical subdivisions.
name : string, optional, default: None.
Descriptive identifier.
Raises
------
ValueError
If nver is not an integer greater or equal to 1.
If nang is not an integer greater or equal to 3.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
Returns
-------
surface : CylindricalSurface object.
"""
Nfaces = { 'd': 6, 's': 6, 'q':6, 'w':6, 'r':3, 'x':3,
'dv': 12, 'sv': 12, 'qv':18, 'wv':18, 'rv':9, 'xv':9 } # 1,3
splitSize = None # FutDev: currently, not user accessable.
dftNfaces = 12
basetype = basetype.lower()
if basetype.endswith('v') :
dftNfaces = 18
subtype = basetype[0]
surface = CylindricalSurface._cyl_vol_square_net( nver ,nang, subtype, minrad, splitSize, numR, name, **kwargs)
else :
surface = cls._square_net(nver ,nang,basetype,minrad,splitSize, name, **kwargs)
surface._postProc_surfaceColors(True)
surface._rez = surface._calc_rez( dftNfaces )
return surface
[docs] @classmethod
def pntsurf(cls, rtz_points, name=None, **kargs) :
"""
CylindricalSurface from points in cylindrical coordinates.
Parameters
----------
rtz_points : array, required.
An N X 3 array of N cylindrical coordinate points.
name : string, optional, default: None.
Descriptive identifier for the point surface.
Returns
-------
surface : CylindricalSurface object
"""
dataPts = np.array(rtz_points)
cylPts = np.array(dataPts.T)
cylPts[0] = np.amax(cylPts[0])
xyz_cylPts = CylindricalSurface.coor_convert(cylPts,True)
minZ,maxZ = np.amin(cylPts[2]), np.amax(cylPts[2])
offsetZ = (maxZ+minZ)/2
xyz_cylPts[2] = xyz_cylPts[2] - offsetZ
sphPts = SphericalSurface.coor_convert(xyz_cylPts,False)
rmax = np.amax(sphPts[0])
sphPts[0] = rmax
xyz_sphPts = SphericalSurface.coor_convert(sphPts,True)
ptTop,ptBtm = [0,0,rmax], [0,0,-rmax]
xyz_sphPts = np.append(xyz_sphPts.T,[ptTop,ptBtm],axis=0).T
# note: the following surface method chull corrects face normal
# orientations so spatial.ConvexHull was not directly used.
hull_surface = Surface3DCollection.chull(xyz_sphPts.T)
verts = CylindricalSurface.coor_convert(dataPts.T,True).T
# remove faces attached to the two top & bottom vertices.
N = len(dataPts)
fvIndices = hull_surface.fvIndices
maxIndices = np.max(fvIndices,axis=1)
shouldkeep = np.where(maxIndices < N)
faceIndices = fvIndices[shouldkeep]
if name is None : name = 'point surface'
surface = CylindricalSurface(name=name,**kargs)
surface.vertexCoor = verts
surface.fvIndices = faceIndices
surface.evIndices = surface._get_edges_from_faces(faceIndices)
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( verts[faceIndices] )
surface._set_geometric_bounds()
surface._rez = surface._calc_rez(12)
surface._basetype = 'pntsurf'
return surface
@classmethod
def _Xfev(cls,rez, basetype=None, minrad=None) :
"""
Calculates number of faces, edges, and vertices of a CylindricalSurface.
Parameters
----------
rez : integer
Number of recursive subdivisions of a triangulated base faces.
basetype : string
Starting surface geometries from which the surface would be
constructed using recursive subdivisions of the triangular faces.
Returns
-------
list: the number of faces, edges, and vertices.
"""
if isinstance(rez,(list,tuple)) : # check for ring...
nver = rez[0]
nang = rez[1]
# only check is ending in s
isSplit = basetype[ len(basetype)-1] == 's'
f = 2*nver*nang
e = 3*nver*nang + nang
v = (nver + 1)*nang
if isSplit :
e += nver
v += nver + 1
return [f,e,v]
return super()._fev(rez,basetype,cls)
def __init__(self, rez=0, basetype=None, name=None, **kwargs):
"""
Create cylindrical surface of unit radius and length 2.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
basetype : {'tri','squ','tri2','squ2','tri_s','squ_s', 'tri_v','squ_v','tri2_v','squ2_v','tri_sv','squ_sv'}, optional, default: 'tri'
Starting surface geometries from which the surface is
constructed using recursive subdivisions of the triangular faces.
name : string, optional, default: None.
Descriptive identifier for the geometry.
Raises
------
ValueError
If rez is not an integer in range 0 to 7.
If basetype is not recognized for this surface constructor.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
"""
if not isinstance(rez, int) :
raise ValueError('Incorrect rez type, must ba of type int')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be an int, 0 <= rez <= {}'.format(rez,_MAXREZ))
if basetype is None : basetype = self._default_base
if basetype not in self._base_surfaces :
raise ValueError('Basetype {} is not recognized. Possible values are: {}'.format(basetype,list(self._base_surfaces)))
baseSurfObj = self._base_surfaces[basetype]
baseVcoor = baseSurfObj['baseVerticesCoor']
baseFaceVertexIndices =baseSurfObj['baseFaceIndices']
indexObj,vertexCoor = self.triangulateBase(rez,baseVcoor,baseFaceVertexIndices, self._midVectorFun)
faceIndices = indexObj['face']
super().__init__(vertexCoor, faceIndices, name, **kwargs)
self.coorType = _COORSYS["CYLINDRICAL"]
self._normal_scale = _normLength( len(faceIndices), 1.14, 1.31)
self._rez = rez
self._basetype = basetype
self.disp_Vector = np.array( [1,1,0] )
return None
[docs] def domain(self, radius=None, zlim=None) :
"""
Set the domain of the cylindrical surface.
Used for setting the limits of the base cylinder.
Parameters
----------
radius: number, default : 1
Radius of the cylindrical surface.
If set to None, the domain is unchanged.
zlim : 2d array
Minimum and maximum values of the arrays are used to scale
and translate the surface from an intial domain of
[ -1, 1 ]
If set to None, the Z-domain is unchanged.
Returns
-------
self : CylindricalSurface object
"""
def setDomain(rtz, mr, mz,bz) :
r,t,z = rtz
R = mr*r
Z = mz*z + bz
return R,t,Z
def getCont_r(lims) :
m = 1.0
if lims is not None : m = lims
return m
def getCont( lims ) :
m, b = 1.0, 0.0
if lims is not None :
# allow single bndry to be input as a float
if isinstance(lims, (int,float)) : lims = [-lims,lims]
m, b = 0.5*(lims[1]-lims[0]) , 0.5*(lims[1]+lims[0])
return m, b
mr = getCont_r(radius)
mz, bz = getCont(zlim)
self.map_geom_from_op(lambda A : setDomain(A, mr, mz,bz))
return self
[docs]class SphericalSurface(Surface3DCollection) :
"""
Spherical 3D surface in spherical coordinates.
Methods are inherited from the Surface3DCollection.
Spherical surface geometries are created in this subclass.
"""
[docs] @staticmethod
def rand(rez=0,seed=None,name=None,kind=None,**kargs) :
"""
SphericalSurface with random triangular faces.
Parameters
----------
rez : number, optional, default: 0
Number of 'recursive' subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
seed : integer, optional, default: None
An initialization seed for random generation.
name : string, optional, default: None.
Descriptive identifier for the geometry.
kind : string, optional, default: None
(reserved, not implemented)
Raises
------
ValueError
If rez is not an integer in range 0 to 7.
If basetype is not recognized for this surface constructor.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
Returns
-------
surface : SphericalSurface object.
"""
#.........................................................
def random_sinDist(size) :
nLoop,nControl,samples = 0, 10**6, []
while len(samples)<size and nLoop<nControl:
x=np.random.random_sample()
prob=np.sin(x*np.pi)/2
assert prob>=0 and prob<=1
if np.random.random_sample() <=prob: samples += [x]
nLoop+=1
return np.array(samples)
#.........................................................
if not isinstance(rez, (int, float)) :
raise ValueError('Incorrect rez type, must be a number.')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be a number, 0 <= rez <= {}'.format(rez,_MAXREZ))
if seed is not None :
if not isinstance(seed, int) :
raise ValueError('Incorrect seed type, must ba of type int')
np.random.seed(seed)
N = int(20*(4**rez)/2 + 2)
r = np.ones( N )
t = 2*np.pi*np.random.rand( N )
p = np.pi*random_sinDist( N )
data = SphericalSurface.coor_convert([r,t,p],True).T
hullsurf = Surface3DCollection.chull(data)
if name is None : name = 'random mesh'
surface = SphericalSurface(name=name,**kargs) # set to the class default object
# reset verts, faceIndices and edges (retain class coordinate system).
verts = hullsurf.vertexCoor
faceIndices = hullsurf.fvIndices
surface.vertexCoor = verts
surface.fvIndices = faceIndices
surface.evIndices = surface._get_edges_from_faces(faceIndices)
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( verts[faceIndices] )
surface._set_geometric_bounds()
surface._rez = str(rez) if isinstance(rez,int) else '{:.1f}'.format(rez)
surface._basetype = 'rand_'+str(seed)
return surface
[docs] @staticmethod
def platonic(rez=0,basetype=None,name=None, **kwargs) :
"""
Platonic Solid surfaces.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
basetype : {'tetra','octa','icosa','cube','dodeca'}, optional, default: 'icosa'
Starting surface geometries from which the surface is
constructed using recursive subdivisions of the triangular faces.
name : string, optional, default: None.
Descriptive identifier for the geometry.
Returns
-------
SphericalSurface object
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.SphericalSurface'.
"""
#.....................................
def fev_cube_a() :
v = [ [-1,-1,-1],[-1, 1,-1],[ 1, 1,-1],[ 1,-1,-1],
[-1,-1, 1],[-1, 1, 1],[ 1, 1, 1],[ 1,-1, 1] ]
f = [ [0,1,2,3],[4,7,6,5],[0,3,7,4],[3,2,6,7],[2,1,5,6],[1,0,4,5]]
e = [ [3,2],[2,1],[1,0],[0,3], [7,6],[6,5],[5,4],[4,7], [0,4],[3,7],[2,6],[1,5]]
return np.array(v),np.array(f),np.array(e)
def _preTriangulateDodoc(v,f,e,fc) :
# (special Case)
# Only used for rez=1 of a platonic dodecahedron. Each pentagon
# is broken into 5 isoceles triangles from a centered vertex.
subFaceIndices,subEdgeIndices = [],[]
newCoors,newColor = [],[]
nvInx = len(v)
for i, face in enumerate(f) :
for n0 in range(5) :
n1 = (n0+1)%5
subFaceIndices.append( [face[n0],face[n1],nvInx+i] )
subEdgeIndices.append( [face[n0],nvInx+i] )
newCoors.append( np.sum(v[face], axis=0)/5 )
newColor += [fc[i]]*5
newverts = np.append(v,newCoors,axis=0)
newfvInx = np.array(subFaceIndices)
newEvInx = np.append(e,subEdgeIndices,axis=0)
newColor = np.array(newColor)
return newverts,newfvInx,newEvInx,newColor
#.....................................
if not isinstance(rez, int) :
raise ValueError('Incorrect rez type, must ba of type int')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be an int, 0 <= rez <= {}'.format(rez,_MAXREZ))
if basetype is None : basetype = SphericalSurface._default_base
baseDict = { 'tetra':"tetrahedron",'octa':"octahedron",
'icosa':"icosahedron",'cube':"cube",'dodeca':"dodecahedron",
'cube_a':"cube_align" }
baseCheck = baseDict.keys()
if basetype not in baseCheck :
raise ValueError('Basetype {} is not recognized. Possible values are: {}'.format(basetype,baseCheck) )
if basetype == 'cube_a' :
surface = SphericalSurface(basetype='cube',name=name, **kwargs)
else :
surface = SphericalSurface(basetype=basetype,name=name, **kwargs)
surface._rez = rez
if basetype == 'cube_a' :
v,f,e = fev_cube_a()
surface.vertexCoor = v
surface.fvIndices = f
surface.evIndices = e
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( v[f] )
if basetype == 'cube' :
v,f,e = SphericalSurface.get_cube()
surface.vertexCoor = v
surface.fvIndices = f
surface.evIndices = e
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( v[f] )
if basetype == 'dodeca' :
v,f,e = SphericalSurface.get_dodecahedron()
if rez>0 :
fc = surface._dfacecolors # update for v_1.2.0
v,f,e,fc = _preTriangulateDodoc(v,f,e,fc)
surface.set_facecolor(fc)
rez -= 1
surface.vertexCoor = v
surface.fvIndices = f
surface.evIndices = e
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( v[f] )
surface._postProc_surfaceColors(True)
if rez > 0 :
surface.triangulate(rez)
surface._set_geometric_bounds()
surface._basetype = basetype + '*'
surface.name = baseDict[basetype] if name is None else name
return surface
[docs] @staticmethod
def get_dodecahedron() :
"""Numpy arrays of dodecahedron vertices (20,3) and face vertex indices (12,5)
and edge vertex indices (30,2)
"""
v,f,e = SphericalSurface._get_platonic_solids_dictionary('dodeca')
return np.array(v), np.array(f), np.array(e)
[docs] @staticmethod
def get_cube() :
"""Numpy arrays of cube vertices (8,3) and face vertex indices 6,5)
and edge vertex indices (12,2)
"""
v,f,e = SphericalSurface._get_platonic_solids_dictionary('cube')
return np.array(v), np.array(f), np.array(e)
@staticmethod
def _get_platonic_solids_dictionary(getArrays = None) :
"""Constructs the dictionary of base spherical surface geometries."""
# ....................................................
def get_cube() :
zeta = 1/3
x0 = np.sqrt(2)/3
y0 = x0*np.sqrt(3)
chi = 2*x0
cubeBaseVert = [ [0,0,1],
[chi,0,zeta], [-x0,y0,zeta], [-x0,-y0,zeta],
[-chi,0,-zeta], [x0,-y0,-zeta], [x0,y0,-zeta],
[0,0,-1] ]
cubeBaseFaces=[
[0,1,6,2],[0,2,4,3],[0,3,5,1],
[7,6,1,5],[7,4,2,6],[7,5,3,4] ]
cubeBaseEdges = [ [0,1], [0,2], [0,3], [7,4], [7,5], [7,6],
[1,6], [6,2], [2,4], [4,3], [3,5], [5,1] ]
return cubeBaseVert, cubeBaseFaces, cubeBaseEdges
def get_dodecahedron() :
zeta = 1/3
ksi = np.sqrt(5)/3
lam_squ = 8/9
hEdge_squ = (1-ksi)/2
Ax = 2/3
Apx = 1/3
Bx = np.sqrt(lam_squ-hEdge_squ)
Cx = ksi
Dx = 1-Bx
Apy = np.sqrt(3)/3
By = np.sqrt(hEdge_squ)
Cy = Apy
Dy = By + Cy
# -----------------------
docTop = [ [0,0,1],
[Ax, 0, ksi], [-Apx, Apy, ksi], [-Apx, -Apy, ksi],
[Cx, Cy, zeta], [Dx, Dy, zeta], [-Bx, By, zeta],
[-Bx, -By, zeta], [Dx, -Dy, zeta], [Cx, -Cy, zeta]
]
docBottom = []
for ver in docTop : docBottom.append( np.multiply(ver,[-1,1,-1]).tolist() )
docbaseVert = np.concatenate((docTop,docBottom)).tolist()
docbaseFaces = [
[0,1,4,5,2], [0,2,6,7,3] , [0,3,8,9,1],
[1,9,17,16,4], [2,5,15,14,6], [3,7,19,18,8],
[10,12,16,17,13], [10,11,14,15,12], [10,13,18,19,11],
[11,19,7,6,14],[13,17,9,8,18],[12,15,5,4,16]
]
docbaseEdges = [
[ 0, 1], [ 0, 2], [ 0, 3], [ 4, 5], [ 6, 7], [ 8, 9],
[ 1, 9], [ 1, 4], [ 2, 5], [ 2, 6], [ 3, 7], [ 3, 8],
[10,11], [10,12], [10,13], [14,15], [16,17], [18,19],
[11,19], [11,14], [12,15], [12,16], [13,17], [13,18],
[ 6,14], [ 4,16], [ 8,18], [ 5,15], [ 9,17], [ 7,19]
]
return docbaseVert, docbaseFaces, docbaseEdges
def construct_tetrahedron() :
cSqr = (np.tan(np.pi/3))*(np.tan(np.pi/3))
zeta = (cSqr-1)/(2*cSqr)
chi = np.sqrt(3*cSqr-1)*np.sqrt(cSqr+1)/(2*cSqr)
x0 = chi*np.cos(np.pi/3)
y0 = chi*np.sin(np.pi/3)
baseVert = [ [0,0,1],
[chi,0,-zeta], [-x0,y0,-zeta], [-x0,-y0,-zeta] ]
baseFaces=([0,1,2],[0,2,3],[0,3,1],[1,3,2])
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
def construct_cube() :
cubeBaseVert, cubeBaseFaces, cubeBaseEdges = get_cube()
baseFaces = []
baseVert = cubeBaseVert.copy()
for face in cubeBaseFaces :
vcSum = [0,0,0]
# get square face center
for vI in face : vcSum = np.add(vcSum,cubeBaseVert[vI])
vcSum = np.divide(vcSum,4)
vcCenter = vcSum/np.linalg.norm(vcSum)
vIndex = len(baseVert)
baseVert.append(vcCenter)
# subdivide cube face into 4 triangles
for i in range(0,4) :
newFace = []
newFace.append(vIndex)
newFace.append(face[i])
newFace.append(face[ (i+1)%4])
baseFaces.append(newFace)
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
def construct_cubeS() :
zeta = 1/3
x0 = np.sqrt(2)/3
y0 = x0*np.sqrt(3)
chi = 2*x0
soff = 0.01
x0S, y0S = soff/2, soff*np.sqrt(3.0)/2.0
w = y0S/2.
cubeBaseVert = [ [chi,-w,zeta], [chi, w,zeta],
[-x0,y0,zeta], [-x0,-y0,zeta], [-chi,0,-zeta], [x0,-y0,-zeta], [x0,y0,-zeta],
[ soff, w, 1], [ x0S, y0S, 1], [ -x0S, y0S, 1], [ -soff, 0, 1], [ -x0S,-y0S, 1],
[ x0S,-y0S, 1], [ soff, -w, 1],
[ soff, w,-1], [ x0S, y0S,-1], [ -x0S, y0S,-1], [ -soff, 0,-1], [ -x0S,-y0S,-1],
[ x0S,-y0S,-1], [ soff, -w,-1] ]
cubeBaseFaces=(
[ 7, 1, 6, 8], [ 8, 6, 2, 9], [ 9, 2, 4,10], [10, 4, 3,11], [11, 3, 5,12], [12, 5, 0,13],
[14,15, 6, 1], [15,16, 2, 6], [16,17, 4, 2], [17,18, 3, 4], [18,19, 5, 3], [19,20, 0, 5] )
baseFaces = []
baseVert = cubeBaseVert.copy()
for face in cubeBaseFaces :
vcSum = [0,0,0]
# get square face center
for vI in face : vcSum = np.add(vcSum,cubeBaseVert[vI])
vcSum = np.divide(vcSum,4)
vcCenter = vcSum/np.linalg.norm(vcSum)
vIndex = len(baseVert)
baseVert.append(vcCenter)
# subdivide cube face into 4 triangles
for i in range(0,4) :
newFace = []
newFace.append(vIndex)
newFace.append(face[i])
newFace.append(face[ (i+1)%4])
baseFaces.append(newFace)
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,8,1] }
def construct_octahedron():
baseVert = [ [0,0,1],
[1,0,0], [0,1,0], [-1,0,0], [0,-1,0],
[0,0,-1] ]
baseFaces=([0,1,2],[0,2,3],[0,3,4],[0,4,1],
[5,2,1],[5,3,2],[5,4,3],[5,1,4])
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
def construct_octahedronS():
soff = 0.01
baseVertS = []
x0S, y0S = soff/np.sqrt(2.0), soff/np.sqrt(2.0)
baseVertS = [ [1,-y0S,0],
[1, y0S,0], [0,1,0], [-1,0,0], [0,-1,0],
[ x0S, y0S, 1], [-x0S, y0S, 1], [-x0S, -y0S, 1], [ x0S, -y0S, 1],
[ x0S, y0S,-1], [-x0S, y0S,-1], [-x0S, -y0S,-1], [ x0S, -y0S,-1] ]
baseFacesS=( [5,1,2], [6,2,3], [7,3,4], [8,4,0],
[2,6,5], [3,7,6], [4,8,7],
[9,2,1], [10,3,2], [11,4,3], [12,0,4],
[2,9,10], [3,10,11], [4,11,12] )
FoS = len(baseFacesS)
return { 'baseFaceIndices' : baseFacesS , 'baseVerticesCoor' : baseVertS, 'fevParam' : [FoS,5,1] }
def construct_dodecahedron() :
docbaseVert, docbaseFaces, docbaseEdges = get_dodecahedron()
baseFaces = []
baseVert = docbaseVert.copy()
for face in docbaseFaces :
vcSum = [0,0,0]
# get pentagon face center
for vI in face : vcSum = np.add(vcSum,docbaseVert[vI])
vcSum = np.divide(vcSum,4)
vcCenter = vcSum/np.linalg.norm(vcSum)
vIndex = len(baseVert)
baseVert.append(vcCenter)
# subdivide pentagon face into 5 triangles
for i in range(0,5) :
newFace = []
newFace.append(vIndex)
newFace.append(face[i])
newFace.append(face[ (i+1)%5])
baseFaces.append(newFace)
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
def construct_icosahedron():
cSqr = (np.tan(np.pi/5))*(np.tan(np.pi/5))
zeta = (1-cSqr)/(2*cSqr)
chi = np.sqrt(3*cSqr-1)*np.sqrt(cSqr+1)/(2*cSqr)
thetaInc = 2*np.pi/5
thetaOff = thetaInc/2
baseVert = []
baseVert.append([0,0,1])
for i in range(0,5) :
theta = i*thetaInc
baseVert.append( [ chi*np.cos(theta), chi*np.sin(theta), zeta] )
for i in range(0,5) :
theta = thetaOff+i*thetaInc
baseVert.append( [ chi*np.cos(theta), chi*np.sin(theta), -zeta] )
baseVert.append([0,0,-1])
baseFaces=([0,1,2],[0,2,3],[0,3,4],[0,4,5],[0,5,1],
[6,2,1],[7,3,2],[8,4,3],[9,5,4],[10,1,5],
[2,6,7],[3,7,8],[4,8,9],[5,9,10],[1,10,6],
[11,7,6],[11,8,7],[11,9,8],[11,10,9],[11,6,10])
Fo = len(baseFaces)
return { 'baseFaceIndices' : baseFaces , 'baseVerticesCoor' : baseVert, 'fevParam' : [Fo,0,2] }
# ....................................................
if getArrays == 'cube' : return get_cube()
if getArrays == 'dodeca' : return get_dodecahedron()
surfacesDictionary = {}
surfacesDictionary['tetra'] = construct_tetrahedron()
surfacesDictionary['cube'] = construct_cube()
surfacesDictionary['cube_s'] = construct_cubeS()
surfacesDictionary['octa'] = construct_octahedron()
surfacesDictionary['octa_s'] = construct_octahedronS()
surfacesDictionary['dodeca'] = construct_dodecahedron()
surfacesDictionary['icosa'] = construct_icosahedron()
# ... 'special case' only for fev method for spherical surfaces.
surfacesDictionary['octa_c'] = { 'fevParam' : [16,5,1] }
surfacesDictionary['cube_c'] = { 'fevParam' : [12,4,1] }
# ... 'special case' only for fev method for platonic surfaces.
surfacesDictionary['cube_a'] = { 'fevParam' : [12,0,2] }
surfacesDictionary['cube*'] = { 'fevParam' : [12,0,2] }
surfacesDictionary['dodeca*'] = { 'fevParam' : [36,0,2] }
return surfacesDictionary
@staticmethod
def _midVectorFun(vectA, vectB) :
"""Mid-point between two points on a unit sphere."""
mid = np.add(vectA,vectB)
return mid/np.linalg.norm(mid)
@staticmethod
def _map3Dto2Dsquare(xyz) :
"""Override map surface to rectagular unit coordinates (0,1)."""
x,y,z = xyz
t = np.arctan2(y,x)
theta = np.where( t<0, 2.0*np.pi + t ,t)
phi = np.arcsin(z)
a = theta/(2*np.pi)
b = phi/np.pi + 0.5
a = np.clip(a,0.0,1.0)
b = np.clip(b,0.0,1.0)
return [a,b]
[docs] @staticmethod
def coor_convert(xyz, tocart=False) :
"""
Transformation between spherical and Cartesian coordinates.
The azimuthal coordinate is in radians with domain 0 to 2pi.
The polar angle is in radians with domain 0 to pi.
Parameters
----------
xyz : 3xN array
N number of vectors in either spherical or cartesian coordinates,
tocart : { True, False }, default: False
If True, input is spherical and output is cartesian.
If False, input is cartesian and output is spherical.
Returns
-------
list: 3xN array
Transformed coordinates.
"""
if tocart :
r, theta, phi = xyz
x = r*np.sin(phi)*np.cos(theta)
y = r*np.sin(phi)*np.sin(theta)
z = r*np.cos(phi)
abc = [x,y,z]
else:
x,y,z = xyz
R = np.sqrt( x*x + y*y + z*z )
phi = np.nan_to_num(np.arccos(z/R))
theta = np.arctan2(y,x)
theta = np.where(theta<0,theta+2*np.pi,theta)
abc = [R,theta,phi]
abc = np.array(abc)
return abc
_base_surfaces = _get_platonic_solids_dictionary.__func__()
_default_base = 'icosa' # default assignment is indicated in __init Docstring.
_geomType = _COORSYS['SPHERICAL']
@classmethod
def _uvw_to_xyz(cls, uvw, rtp) :
"""Override rotational transformation at a spherical coordinate."""
return super()._disp_to_xyz( uvw, rtp[1], rtp[2] )
[docs] @classmethod
def grid(cls,nlat, nlng, basetype='d', minrad=None, name=None, **kwargs) :
"""
SphericalSurface with latitude and longitude faces.
Parameters
----------
nlat : integer, required.
Number of latitude subdivisions. Minimum value is 2
nlng : integer, required.
Number of longitude subdivisions. Minimum value is 3.
basetype : {'d','s','q','w','r','x'}, optional, default: 'd'
Starting surface geometries from which the surface is
constructed indicating face shape and if split.
minrad : scalar, optional, default: 0.01
The minimum distance from the z-axis of any vertex coordinate.
Use is dependent on basetype.
name : string, optional, default: None.
Descriptive identifier.
Raises
------
ValueError
If nlat is not an integer greater or equal to 2.
If nlng is not an integer greater or equal to 3.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
Returns
-------
surface : SphericalSurface object.
"""
Nfaces = { 'd': 6, 's': 6, 'q':12, 'w':12, 'r':6, 'x':6 } # 2,3
splitSize = None
nrad, nang = nlat, nlng
surface = cls._square_net(nrad, nang, basetype,minrad, splitSize, name, **kwargs)
surface._postProc_surfaceColors(True)
surface._rez = surface._calc_rez( 20 )
return surface
[docs] @classmethod
def pntsurf(cls, rtp_points, name=None, **kargs) :
"""
SphericalSurface from points in spherical coordinates.
Parameters
----------
rtp_points : array, required.
An N X 3 array of N spherical coordinate points.
name : string, optional, default: None.
Descriptive identifier for the point surface.
Returns
-------
surface : SphericalSurface object
"""
dataPts = np.array(rtp_points)
sphPts = dataPts.T
xyzPts = SphericalSurface.coor_convert(dataPts.T,True)
rmax = np.amax(sphPts[0])
sphPts[0] = rmax
xyz_sphPts = SphericalSurface.coor_convert(sphPts,True)
hull_surface = Surface3DCollection.chull(xyz_sphPts.T)
verts = xyzPts.T
faceIndices = hull_surface.fvIndices
if name is None : name = 'point surface'
surface = SphericalSurface(name=name,**kargs)
surface.vertexCoor = verts
surface.fvIndices = faceIndices
surface.evIndices = surface._get_edges_from_faces(faceIndices)
surface.vfIndicesList = surface._vertexCommonFaceIndices()
surface.efIndicesList = surface._edgeCommonFaceIndices()
surface.set_verts( verts[faceIndices] )
surface._set_geometric_bounds()
surface._rez = surface._calc_rez(20)
surface._basetype = 'pntsurf'
return surface
@classmethod
def _Xfev(cls, rez, basetype=None, minrad=None) :
"""
Calculates number of faces, edges, and vertices of a SphericalSurface.
Parameters
----------
rez : integer
Number of recursive subdivisions of a triangulated base faces.
basetype : string
Starting surface geometries from which the surface would be
constructed using recursive subdivisions of the triangular faces.
Returns
-------
list: the number of faces, edges, and vertices.
"""
if isinstance(rez,(list,tuple)) : # check for globe...
nlat = rez[0]
nlng = rez[1]
# only check is ending in s
isSplit = basetype[ len(basetype)-1] == 's'
if minrad is None :
f = 2*nlat*nlng - 2*nlng
e = 3*nlat*nlng - 3*nlng
v = nlat*nlng - nlng + 2
if isSplit :
e += nlat
v += nlat - 1
else : # any value will trigger this condition.
f = 2*nlat*nlng
e = 3*nlat*nlng + nlng
v = nlat*nlng + nlng
if isSplit :
e += nlat
v += nlat + 1
return [f,e,v]
return super()._fev(rez,basetype,cls)
@staticmethod
def _spherical_split_cylinder(rez,basetype,minrad, **kwargs) :
"""
Transform cylindrical base surface to spherical coordinates.
note : This method provides an alternative to the
spherical basetypes of 'hex_s' and 'square_s'.
The grids constructed here may provide smoother
geometries depending on the mapping functions.
"""
# .................................................
def cyl2sph(rtz):
r,t,z = rtz
minphi = np.arccos(minrad)
phi = z*minphi
Z = np.sin(phi)
R = np.cos(phi)
return R,t,Z
# .................................................
if basetype == 'octa_c' : cylbase = 'squ_s'
if basetype == 'cube_c' : cylbase = 'tri_s'
cyl_obj = CylindricalSurface(rez,cylbase, **kwargs).map_geom_from_op(cyl2sph)
baseFaceVertexIndices = cyl_obj._base_surfaces[cylbase] ['baseFaceIndices']
fvIndices = cyl_obj.fvIndices
vertexCoor = cyl_obj.vertexCoor
evIndices = cyl_obj.evIndices
return vertexCoor, fvIndices, evIndices, baseFaceVertexIndices,
def __init__(self, rez=0, basetype=None, minrad=None, name=None, **kwargs):
"""
Create spherical surface of unit radius.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of the triangulated base faces.
Rez values range from 0 to 7.
basetype : {'tetra','octa','icosa','cube','dodeca', \
'octa_s','cube_s','octa_c','cube_c'}, optional, default: 'icosa'
Starting surface geometries from which the surface is
constructed using recursive subdivisions of the triangular faces.
All basetypes are based on the five platonic solids.
minrad : scalar, optional, default: 0.01
For basetypes 'octa_c' and 'cube_c, the minimum
distance from the z-axis of any vertex coordinate.
name : string, optional, default: None.
Descriptive identifier for the geometry.
Raises
------
ValueError
If rez is not an integer in range 0 to 7.
If basetype is not recognized for this surface constructor.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.Surface3DCollection'.
"""
if minrad is None : minrad = _MINRAD
if not isinstance(rez, int) :
raise ValueError('Incorrect rez type, must ba of type int')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be an int, 0 <= rez <= {}'.format(rez,_MAXREZ))
if basetype is None : basetype = self._default_base
# ... check for 'special case' (ugly code but needed, MROI)
if basetype == 'octa_c' or basetype == 'cube_c' :
vertexCoor, faceIndices, edgeIndices, bfvi = self._spherical_split_cylinder(rez,basetype,minrad, **kwargs)
baseFaceVertexIndices = bfvi
else :
if basetype not in self._base_surfaces :
raise ValueError('Basetype {} is not recognized. Possible values are: {}'.format(basetype,list(self._base_surfaces)))
baseSurfObj = self._base_surfaces[basetype]
baseVcoor = baseSurfObj['baseVerticesCoor']
baseFaceVertexIndices =baseSurfObj['baseFaceIndices']
indexObj,vertexCoor = self.triangulateBase(rez,baseVcoor,baseFaceVertexIndices, self._midVectorFun)
faceIndices = indexObj['face']
edgeIndices = indexObj['edge']
super().__init__(vertexCoor, faceIndices, name, **kwargs)
self.coorType = _COORSYS["SPHERICAL"]
self._normal_scale = _normLength( len(faceIndices), 0.61, 1.52)
self._rez = rez
self._basetype = basetype
self.disp_Vector = np.array( [1,1,1] )
return None
[docs] def domain(self, radius=None) :
"""
Set the domain of the spherical surface.
Used for setting the radius of the base sphere.
Parameters
----------
radius: number, default : 1
Radius of the spherical surface.
If set to None, the radius is unchanged.
Returns
-------
self : SphericalSurface object
"""
def setDomain(rtp, mr ) :
r,t,p = rtp
R = mr*r
return R,t,p
def getCont(lims) :
m = 1.0
if lims is not None : m = lims
return m
mr = getCont(radius)
self.map_geom_from_op(lambda A : setDomain(A, mr ))
return self
[docs]class Vector3DCollection(Line3DCollection) :
"""
Collection of 3D vectors represented as arrows.
"""
# FutDev: additional methods to add functionality.
def __init__(self, location, vect , alr=None, name=None, **kwargs) :
"""
Create a collection of 3D vectors.
Parameters
----------
location : N x 3 float array
Cartesian coordinate location (tails) for N number of vectors.
vect : N x 3 float array
N number of vectors in Cartesian coordinates.
alr : scalar, optional, default: 0.25
Axis length ratio, head size to vector magnitude.
name : string identifier
Other Parameters
----------------
**kwargs
All other parameters are passed on to mpl_toolkits.mplot3d.art3d.Line3DCollection.
Valid keywords include: colors, linewidths.
"""
# note: self coor types primarially used for export info.
self.coorType = _COORSYS["XYZ"]
self.uvwOrientationType = _COORSYS["XYZ"]
if alr is None : alr = _DFT_ALR
self._alr = alr
self.baseLocation = np.array(location, dtype=float)
self.baseVect = np.array(vect, dtype=float)
lines,XYZ,UVW = self._quiverLines(location,vect,alr)
# note: _quiverLines may remove zero length vect
self.location = XYZ
self.vect = UVW
if 'color' in kwargs :
try:
test = colors.to_rgba( kwargs['color'] )
except:
del kwargs['color']
warnings.warn("Only a SINGLE color value may be used for vector instantiation")
super().__init__(lines, **kwargs)
#self._geomName = name # Redunant from next line, internal tracking to determine if set...
self.name = name
self.valuesName = None
self.baseColor = np.array(self.get_color()[0]).flatten().tolist()
self._bounds = {}
_set_geomBounds(self.location,self._bounds)
self._bounds['vlim'] = _DFT_VLIM #.. from operation
self._bounds['vertvlim'] = _DFT_VERTVLIM #.. from direction or magnitude
def _set_vectColor(self, color) :
# note: for a list of four colors [ A,B,C,D ], result is
# [ A,B,C,D, A,A, B,B, C,C, D,D ].
# line segment colors, then a pair for each arrow head (2 lines).
head = np.arange(len(color))
tail = np.ravel( [ [head],[head] ] , order='F')
index = np.append(head,tail)
#self._edgecolors = color[index]
self.set_color(color[index]) # update for v_1.2.0
return
def __str__(self) :
# note: each vector representation is composed of 3 lines
numbVect = int(len(self._segments3d)/3)
name = self.__class__.__name__
val = name + ': vectors: {}'.format(numbVect)
return val
def _quiverLines( self, vC, vN, alr) :
# ..................................................
def quiver( *args, length=1, arrow_length_ratio=alr, pivot='tail', normalize=False):
# taken directly from mpl_toolkits.mplot3d.axes3d source
# and removed axis reference. Acknowledgement:
# Copyright (c) 2012-2013 Matplotlib Development Team; All Rights Reserved
# ...............................................
def calc_arrow(uvw, angle=15):
"""
To calculate the arrow head. uvw should be a unit vector.
We normalize it here:
"""
# get unit direction vector perpendicular to (u,v,w)
norm = np.linalg.norm(uvw[:2])
if norm > 0:
x = uvw[1] / norm
y = -uvw[0] / norm
else:
x, y = 0, 1
# compute the two arrowhead direction unit vectors
ra = np.radians(angle)
c = np.cos(ra)
s = np.sin(ra)
# construct the rotation matrices
Rpos = np.array([[c+(x**2)*(1-c), x*y*(1-c), y*s],
[y*x*(1-c), c+(y**2)*(1-c), -x*s],
[-y*s, x*s, c]])
# opposite rotation negates all the sin terms
Rneg = Rpos.copy()
Rneg[[0, 1, 2, 2], [2, 2, 0, 1]] = \
-Rneg[[0, 1, 2, 2], [2, 2, 0, 1]]
# multiply them to get the rotated vector
return Rpos.dot(uvw), Rneg.dot(uvw)
# ...............................................
argi = 6
if len(args) < argi:
raise ValueError('Wrong number of arguments. Expected %d got %d' %
(argi, len(args)))
input_args = args[:argi]
masks = [k.mask for k in input_args
if isinstance(k, np.ma.MaskedArray)]
bcast = np.broadcast_arrays(*input_args, *masks)
input_args = bcast[:argi]
masks = bcast[argi:]
if masks:
mask = reduce(np.logical_or, masks)
input_args = [np.ma.array(k, mask=mask).compressed()
for k in input_args]
else:
input_args = [np.ravel(k) for k in input_args]
if any(len(v) == 0 for v in input_args):
return []
shaft_dt = np.array([0, length])
arrow_dt = shaft_dt * arrow_length_ratio
if pivot == 'tail':
shaft_dt -= length
elif pivot == 'middle':
shaft_dt -= length/2.
elif pivot != 'tip':
raise ValueError('Invalid pivot argument: ' + str(pivot))
XYZ = np.column_stack(input_args[:3])
UVW = np.column_stack(input_args[3:argi]).astype(float)
norm = np.linalg.norm(UVW, axis=1)
mask = norm > 0
XYZ = XYZ[mask]
if normalize:
UVW = UVW[mask] / norm[mask].reshape((-1, 1))
else:
UVW = UVW[mask]
if len(XYZ) > 0:
shafts = (XYZ - np.multiply.outer(shaft_dt, UVW)).swapaxes(0, 1)
head_dirs = np.array([calc_arrow(d) for d in UVW])
heads = shafts[:, :1] - np.multiply.outer(arrow_dt, head_dirs)
heads.shape = (len(arrow_dt), -1, 3)
heads = heads.swapaxes(0, 1)
lines = [*shafts, *heads]
else:
lines = []
return lines , XYZ,UVW # <--- added return of XYZ,UVW
# ..................................................
x,y,z = np.transpose( vC )
u,v,w = np.transpose( vN )
lines, XYZ, UVW = quiver(x,y,z,u,v,w)
return lines, XYZ, UVW
@property
def alr(self) :
'''Arrow head to length ratio.'''
return self._alr
@alr.setter
def alr(self, val) :
location = self.baseLocation
vect = self.baseVect
if val is None : val = _DFT_ALR
self._alr = val
lines,_,_ = self._quiverLines(location,vect,self._alr)
self.set_segments(lines)
return
@property
def vlim(self) :
'''Range of values associated with color'''
return self._bounds['vlim']
@vlim.setter
def vlim(self,val) :
self._bounds['vlim'] = val
return
@property
def name(self) :
'''Descriptive identifier for the vector geometry.'''
if self._geomName is None : return ''
return self._geomName
@name.setter
def name(self, val) :
self._geomName = val
self.set_label( '' if val is None else val)
return
@property
def cname(self) :
'''Descriptive identifier for values indicated by color.'''
if self.valuesName is None : return ''
return self.valuesName
@cname.setter
def cname(self, val) :
self.valuesName = val
return
@property
def bounds(self) :
"""
Dictionary of vector geometric and value ranges.
Each dictionary value is a 2 float array of minimum
and maximum values of the vector location. Keys are:
'xlim' : x-coordinate
'ylim' : y-coordinate
'zlim' : z-coordinate
'r_xy' : radial distance from the z axis
'rorg' : radial distance from the origin
'vlim' : value functional assignments.
'vertvlim: magnitude
Values are assigned from the geometry and color mapping methods.
"""
return self._bounds
[docs] def map_color_from_op( self, operation, rgb=True, cname=None ) :
"""
Assignment of vector color from a function.
Vector colors are assigned from a function
of direction and location coordinates.
Parameters
----------
operation : function object
Function that takes two arguments, both
a 3xN Numpy array of xyz coordinates.
The first and second arguments are the location
and direction, respectively.
The function returns a 3xN color value.
rgb : bool {True, False}, optional, default: True
By default, RGB color values are returned by the
operation function. If set False, the operation
returns HSV color values.
Returns
-------
self : Vector3DCollection object
"""
xyz = np.transpose(self.location)
uvw = np.transpose(self.vect)
# .....
colors = np.array(operation(xyz,uvw))
colors = np.transpose(colors)
colors = np.clip(colors,0,1)
if rgb : RGB = colors
else : RGB = cm.colors.hsv_to_rgb(colors)
if RGB.shape[1] == 3 :
ones = np.ones(RGB.shape[0])[:,np.newaxis]
RGB = np.concatenate((RGB,ones),axis=1)
self._set_vectColor(RGB)
cname = _getFunctionName(operation,cname)
if cname is not None : self.valuesName = cname
return self
[docs] def map_cmap_from_op( self, operation, cmap=None, cname=None ) :
"""
Functional assignment of a vector color from a color map.
Location and direction coordinates are used to calculate a scalar
which is then used to assign face colors from a
colormap.
Parameters
----------
operation : function object
Function that takes two arguments, both
a 3xN Numpy array of xyz coordinates.
The first and second arguments are the location
and direction, respectively.
The function returns a Numpy array of scalar values.
cmap : str or Colormap, optional
A Colormap instance or registered colormap name.
If not assigned, the surface Colormap is used.
The colormap maps the function return values to colors.
Returns
-------
self : Vector3DCollection object
"""
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cm_colorMap = self.get_cmap()
xyz = np.transpose(self.location)
uvw = np.transpose(self.vect)
# .....
v = np.array(operation(xyz,uvw))
norm = colors.Normalize(vmin=v.min(),vmax=v.max())
color = cm_colorMap(norm(v))
self._set_vectColor(color)
cname = _getFunctionName(operation,cname)
if cname is not None: self.valuesName = cname
self._bounds['vlim'] = [ v.min(), v.max() ]
return self
[docs] def map_cmap_from_magnitude( self, cmap=None, cname=None ) :
"""
Vector color assignment using vector magnitude.
Parameters
----------
cmap : str or Colormap, optional, default: 'viridis'
A Colormap instance or registered colormap name.
If not assigned, the surface Colormap is used.
The colormap maps the dot product values to colors.
Returns
-------
self : Vector3DCollection object
"""
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cm_colorMap = self.get_cmap()
mag = np.linalg.norm(self.vect,axis=1)
norm = colors.Normalize(vmin=mag.min(),vmax=mag.max())
color = cm_colorMap(norm(mag))
self._set_vectColor(color)
self.cname = cname
if cname is None: self.cname = 'Magnitude'
self._bounds['vertvlim'] = [ mag.min(), mag.max() ]
return self
[docs] def map_cmap_from_direction( self, cmap=None, direction=[1,1,1], cname=None ) :
"""
Vector color assignment using vector direction relative to direction argument.
The dot product of vector direction with the argument direction
is used to assign vector colors from a colormap.
Parameters
----------
cmap : str or Colormap, optional, default: 'viridis'
A Colormap instance or registered colormap name.
If not assigned, the surface Colormap is used.
The colormap maps the dot product values to colors.
direction : list of size 3, optional, default: [1,1,1]
A 3D vector in xyz Cartesian coordinates designating
the reference direction.
refCoor : string, optional, default: "XYZ"
Direction coordinate system for the evaluation.
(Not implimented)
Returns
-------
self : Vector3DCollection object
"""
# FutDev: allow input direction to be in other coor systems (refCoor)
#.....................................
def getColorMap(fvo,cmap, direction) :
unitVector = lambda v : np.divide( v, np.linalg.norm(v) )
incidentLight = unitVector( direction )
d = np.dot(fvo,incidentLight)
v = np.add(d,1)
v = np.divide(v,2)
v = np.abs(v)
norm = colors.Normalize(vmin=v.min(),vmax=v.max())
colorMap = cmap(norm(v))
self._bounds['vertvlim'] = [ v.min(), v.max() ]
return colorMap
#.....................................
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cm_colorMap = self.get_cmap()
color = getColorMap(self.vect, cm_colorMap, direction)
self._set_vectColor(color)
self.cname = cname
if cname is None: self.cname = 'Direction'
return self
[docs]class ColorLine3DCollection(Line3DCollection) :
"""
Base class for non-uniformly colored 3D lines.
"""
# FutDev: This class leaves a lot of room for improvement.
# Need better use of Numpy and better approach to
# differentiating between collections, lines & segments
# to accomodate Matplot rendering and future export. MROI
[docs] @staticmethod
def meshgrid(X,Y,Z, name=None) :
"""
ColorLine3DCollection, a set of xy slice lines based on a meshgrid.
Parameters
----------
X,Y,Z : N x M arrays
Cartesian x,y,z coordinate values.
name : string, optional, default: None.
Descriptive identifier for the mesh surface.
Returns
-------
surface : ColorLine3DCollection object
"""
# Needs cleanup MROI
x_lineVertsPts = np.array(X)[0]
y_lineVertsPts = np.array(Y)[:,0]
z_lines = np.array(Z)
numb_X_lines = len(x_lineVertsPts)
numb_Y_lines = len(y_lineVertsPts)
indices = [None]*( numb_X_lines + numb_Y_lines )
index = 0
verts = []
# ..................................................
numVerts = numb_Y_lines
lineVertsPts = y_lineVertsPts
z_X_lines = z_lines
for iX in range(numb_X_lines) :
for iY in range(numVerts) :
verts.append( [ x_lineVertsPts[iX] , lineVertsPts[iY], z_X_lines[iY,iX] ])
for i in range(0,numb_X_lines) :
indices[i] = [None]*(numVerts)
for j in range(numVerts) :
indices[i][j] = index
index += 1
# ..................................................
numVerts = numb_X_lines
lineVertsPts = x_lineVertsPts
z_Y_lines = z_lines
for iY in range(numb_Y_lines) :
for iX in range(numVerts) :
verts.append( [ lineVertsPts[iX], y_lineVertsPts[iY] ,z_Y_lines[iY,iX] ])
for i in range(numb_X_lines,numb_X_lines+numb_Y_lines) :
indices[i] = [None]*(numVerts)
for j in range(numVerts) :
indices[i][j] = index
index += 1
verts = np.array(verts)
return ColorLine3DCollection(verts,indices,name=name)
def __init__(self, vertexCoor, segmIndices, name=None, **kwargs) :
"""
A set 3D lines, each line with a varying number of sequential vertices.
Parameters
----------
vertexCoor : V x 3 float array-like
An array of V number of xyz vertex coordinates.
segmIndices : L x Nl int list
An list of L number of Nl line vertex indices, where
Nl is the number of vertices for the Lth line.
name : string, optional, default: None.
Descriptive identifier.
Other Parameters
----------------
**kwargs
All other parameters are passed on to mpl_toolkits.mplot3d.art3d.Line3DCollection.
Valid keywords include: colors, linewidths.
"""
self.vertexCoor = np.array(vertexCoor,dtype=float)
if segmIndices is None :
self.lineIndices = [ [*range(len(vertexCoor))] ]
linesegs = self.lineIndices
indices = []
for line in linesegs :
numSegs = len(line) - 1
for seg in range (numSegs) :
t = [ line[seg],line[seg+1] ]
indices.append(t)
self.segmIndices = indices
else :
# all indices must be lists, not numpy arrays.
segmIndices = list(segmIndices)
# segmIndices will be shredded for any color operation
# forming multiple single segment lines.
self.segmIndices = segmIndices.copy()
# lineIndices are appended for an add operation,
# where lineIndices preserve individual 'lines' even though
# all lines are broken into single segment lines omce shredded.
self.lineIndices = segmIndices.copy()
self.isShredded = False
S = self.segmIndices
v = self.vertexCoor
segments = [None]*len(S)
for verlist in range(len(S)) :
verArr = [None]*len( S[verlist] )
for verInx in range( len( S[verlist] ) ) :
verArr[verInx] = v[ S[verlist][verInx] ]
segments[verlist] = verArr
super().__init__(segments, **kwargs)
#self._geomName = None # Redundant from next line, internal tracking to determine if set...
self.name = name
self.valuesName = None
# ...............................................................
self.coorType = _COORSYS["XYZ"]
self.baseColor = np.array(self.get_color()[0]).flatten().tolist()
self.vertexColor = np.array( [ self.baseColor ] )
# following assigned when surface color applied using an operation
self.vertexValues = None
self.segValues = None
# ...............................................................
self.set_joinstyle('round')
self.set_capstyle('round')
self.scale = [1,1,1] # for axis when transformed.
self.translation = [0,0,0] # for axis when transformed.
self.rotation = np.identity(3).tolist() # for axis when transformed.
self.restoredClip = None
self.rez = None
self._bounds = {}
self._set_geometric_bounds()
self._bounds['vlim'] = _DFT_VLIM
self._bounds['vertvlim'] = _DFT_VERTVLIM
self._planeNormal = None #.. set for planar contours, for FutDev
@staticmethod
def _midVectorFun(vectA, vectB) :
"""Mid-point between two points in the xy plane."""
mid = np.add(vectA,vectB)
mid = np.multiply(0.5,mid)
return mid
@staticmethod
def _extractLines(lineSegments) :
# Edge array sorting to line arrays, which may be loops. Edge arrays
# must have two indices in consistent order.
# This is an 'unshred' function.
# ..................................................................
def extractLine(edgeDict) :
# use a dictionary to build a 'sort-of' single linked list
# note: index-list would conain many emply values.
edgeDict_fl = edgeDict
k = list(edgeDict.keys())[0]
v = list(edgeDict.values())[0]
line = [k,v]
edgeDict_fl.pop(k)
# add segments to the tail :
continueBuild = True
while continueBuild :
continueBuild = False
end = line[len(line)-1]
if end in edgeDict_fl :
line.append(edgeDict_fl[end])
edgeDict_fl.pop(end)
continueBuild = True
if line[len(line)-1] == line[0] :
return line,edgeDict_fl
# then line is not a loop so continue, adding segments to head :
edgeDict_lf = { value:key for key,value in edgeDict_fl.items()}
continueBuild = True
while continueBuild :
continueBuild = False
start = line[0]
if start in edgeDict_lf :
line.insert(0,edgeDict_lf[start])
edgeDict_lf.pop(start)
continueBuild = True
edgeDict_nxt = { value:key for key,value in edgeDict_lf.items() }
return line, edgeDict_nxt
# ..................................................................
edgeDict = { litem[0]:litem[1] for litem in lineSegments }
lines = []
while(len(edgeDict)>0) :
line,edgeDict = extractLine(edgeDict)
lines.append(line)
return lines
[docs] def __add__(self, othr) :
"""
Combine two line objects into a single line object.
Parameters
----------
othr : ColorLine3DCollection object
Returns
-------
ColorLine3DCollection object
"""
# ............................................................................
def get_lineColors( colors, lines) :
colors = np.array(colors)
nlines = len(lines)
if len(colors) != nlines :
colors = np.tile(colors,(nlines,1))
return colors
# ............................................................................
if not isinstance(othr, ColorLine3DCollection) :
raise ValueError('Add operations can only apply between ColorLine3DCollection objects.')
# vertexCoor are array-like.
topVerCoor = self.vertexCoor
botVerCoor = othr.vertexCoor
totVerCoor = np.append(topVerCoor,botVerCoor,axis=0)
nextIndex = len(topVerCoor)
# lineIndices are int Lists
toplineIndices = self.lineIndices
botlineIndices = copy.copy(othr.lineIndices)
for i in range(len(botlineIndices)) :
botlineIndices[i] = np.add( botlineIndices[i] , nextIndex ).tolist()
totLineIndices = copy.copy(toplineIndices)
for line in botlineIndices :
totLineIndices.append(line)
# segmIndices are int Lists
topSegmIndices = self.segmIndices
botSegmIndices = copy.copy(othr.segmIndices)
for i in range(len(botSegmIndices)) :
botSegmIndices[i] = np.add( botSegmIndices[i] , nextIndex ).tolist()
totSegmIndices = copy.copy(topSegmIndices)
for line in botSegmIndices :
totSegmIndices.append(line)
# edgecolors are array-like
topEdgeColors = get_lineColors(self._ledgecolors,topSegmIndices) # update for v_1.2.0
botEdgeColors = get_lineColors(othr._ledgecolors,botSegmIndices) # update for v_1.2.0
totEdgeColors = np.append(topEdgeColors,botEdgeColors,axis=0)
# FutDev: 'add' will not retain line collections. Needed for export
# but sufficient for rendering. MROI
newLine = ColorLine3DCollection(totVerCoor,totSegmIndices,color=totEdgeColors)
newLine.isShredded = self.isShredded and othr.isShredded
newLine._planeNormal = None #.. could check here if same normals.
newLine._set_geometric_bounds()
return newLine
def __str__(self) :
# note: each line representation is composed of multiple segment lines.
numVerts = len(self.vertexCoor)
numLines = len(self.lineIndices)
numSeg = 0
for line in self._segments3d :
numSeg += len(line) -1
name = self.__class__.__name__
if self.rez is not None : name = name + ' ({}):'.format(self.rez)
val = name + ': lines: {}, segs: {}, verts: {}'.format(numLines,numSeg,numVerts)
return val
def _postProc_segmentColors(self) :
""" Matplotlib will 'fillin' colors during rendering using the
facecolor array. (this is NOT the assigned color to edges).
However, S3Dlib will use the inherited _edgecolors
for the assigned edeg colors.
Whenever color is 'externally' assigned or reassigned,
need to directly assign colors for each edge for any
method that uses edgecolors. This occurs initially
in the method before colors are assigned
(set_color, +, clipping, shading, edges, etc.)
"""
N = len(self.segmIndices) # number of segments
initEC = self._ledgecolors # update for v_1.2.0
inLen = len(initEC) # number of edgecolors
if inLen == N : return
if isinstance(initEC,np.ndarray) : initEC = initEC.tolist()
n = math.ceil(N/inLen)
colorcoll = initEC*n
self.set_color(colorcoll[0:N])
return
def _shred(self) :
# convert multi-segment lines to multiple one-segment lines.
# Required for application of non-uniform line color.
linesegs = self.lineIndices
indices = []
for line in linesegs :
numSegs = len(line) - 1
for seg in range (numSegs) :
t = [ line[seg],line[seg+1] ]
indices.append(t)
segments = self.vertexCoor[np.array(indices)]
self.segmIndices = indices
self.set_segments(segments)
self.isShredded = True
return segments
def _shred_segm_color(self,numDiv) :
initColors = self._ledgecolors # colors before line shredding # update for v_1.2.0
numColors = len(initColors)
linesegs = self.segmIndices
colors = []
colorIndex = 0
for line in linesegs :
currentColor = initColors[colorIndex]
for seg in range (numDiv) :
colors.append(currentColor)
colorIndex = (colorIndex+1)%numColors
colors = np.array(colors) # colors after line shredding
return colors
def _shred_color(self) :
# Since line colors are in a SINGLE List of colors
# for mulitple lines with multiple segment colors,
# colors for a multi-segment line of a single color
# is reconstructed into multiple colors of the
# same color. ( see _shred )
# Individual segemnt colors are needed for line
# shading, clipping, and alpha setting.
initColors = self._ledgecolors # colors before line shredding # update for v_1.2.0
numColors = len(initColors)
linesegs = self.lineIndices
colors = []
colorIndex = 0
for line in linesegs :
numSegs = len(line) - 1
currentColor = initColors[colorIndex]
for seg in range (numSegs) :
colors.append(currentColor)
colorIndex = (colorIndex+1)%numColors
colors = np.array(colors) # colors after line shredding
return colors
def _get_segment_centers(self,segments) :
# segment centers from segment point coordinates.
# Assumed segments of a shredded line.
vfT = np.transpose(segments, (1,2,0) )
numVerts = vfT.shape[0]
sumV = np.sum(vfT, axis=0)
return np.transpose(sumV)/numVerts
def _get_segment_directions(self,segments) :
# segment unit directions from segment point coordinates.
# Assumed segments of a shredded line.
vfT = np.transpose(segments, (1,0,2) )
direction = vfT[1] - vfT[0]
sz = np.linalg.norm(direction, axis=1)[:,np.newaxis]
unitVector = np.divide( direction, sz )
return unitVector
def _set_geometric_bounds(self) :
_set_geomBounds(self.vertexCoor,self._bounds)
return
def _get_segments( self, v, S ) :
# Since line segments are in a List of multiple lines,
# line vertices are set per line in the line collection
# to reconstruct the segment coordinates ( self._segments3d )
# DevNote: got to be a better way of doing this, but line
# sets are lists with varying number of segments
# until shredded, then all the same number of segs.
if self.isShredded : return v[np.array(S)]
segments = [None]*len(S)
for verlist in range(len(S)) :
verArr = [None]*len( S[verlist] )
for verInx in range( len( S[verlist] ) ) :
verArr[verInx] = v[ S[verlist][verInx] ]
segments[verlist] = verArr
return segments
# ---------------------------------------------------------------------+
# Note: |
# the following five _c...._vert methods assign color and values for |
# each vertex which are only used for FutDev export of vertex colors |
# (segment center colors are used for the matplotlib rendering) |
# ---------------------------------------------------------------------+
def _map_cmap_from_sequence_vert(self,cmap) :
i = np.linspace(0.0,1.0, len(self.vertexCoor) )
norm = colors.Normalize(vmin=0.0,vmax=1.0)
color = cmap(norm(i))
self.vertexColor = color
self.vertexValues = list( range(len(self.vertexCoor)) )
self._bounds['vertvlim'] = [ 0.0, 1.0 ]
return
def _map_cmap_from_direction_vert( self, edgecolors, edgevalues, vBounds ) :
# since the direction is dependent on two vertices,
# the segment colors are used for the vertexColors.
# The last vertex color is assigned the last segment
# segment color as an approximation.
numEdges = len(edgecolors)
lastColor = edgecolors[numEdges-1]
vColor = copy.copy(edgecolors)
vColor = np.append(vColor,[lastColor],axis=0)
lastValue = edgevalues[numEdges-1]
vValue = copy.copy(edgevalues)
vValue = np.append(vValue,[lastValue],axis=0)
self.vertexColor = vColor
self.vertexValues = vValue
self._bounds['vertvlim'] = vBounds
return
def _map_cmap_from_datagrid_vert(self,g, delta, dmin, cmap) :
x,y,_ = np.transpose(self.vertexCoor)
# normalize xy plane -1 to 1
x = np.interp( x, (x.min(), x.max()), (-1.0, 1.0) )
y = np.interp( y, (y.min(), y.max()), (-1.0, 1.0) )
trial = [ g( x[i],y[i])[0] for i in range(len(x)) ]
v = np.array(trial)*delta + dmin
vmin, vmax = v.min(), v.max()
norm = colors.Normalize(vmin=v.min(),vmax=v.max())
color = cmap(norm(v))
self.vertexColor = color
self.vertexValues = v
self._bounds['vertvlim'] = [ vmin, vmax ]
return
def _map_cmap_from_op_vert(self, operation, cmap) :
xyz = np.transpose(self.vertexCoor)
v = np.array(operation(xyz))
vmin, vmax = v.min(), v.max()
norm = colors.Normalize(vmin=v.min(),vmax=v.max())
color = cmap(norm(v))
self.vertexColor = color
self.vertexValues = v
self._bounds['vertvlim'] = [ vmin, vmax ]
return
def _map_color_from_op_vert(self, operation, rgb) :
xyz = np.transpose(self.vertexCoor)
colors = np.array(operation(xyz))
colors = np.transpose(colors)
colors = np.clip(colors,0,1)
if rgb : RGB = colors
else : RGB = cm.colors.hsv_to_rgb(colors)
if RGB.shape[1] == 3 :
ones = np.ones(RGB.shape[0])[:,np.newaxis]
RGB = np.concatenate((RGB,ones),axis=1)
self.vertexColor = RGB
self.vertexValues = None
self._bounds['vertvlim'] = [ 0.0, 1.0 ]
return
# ---------------------------------------------------------------------
@property
def _ledgecolors(self) :
return self._edgecolors
@_ledgecolors.setter
def _ledgecolors(self,val) :
self._edgecolors = val
return
@property
def vlim(self) :
'''Range of scalar values which may be associated with color'''
return self._bounds['vlim']
@vlim.setter
def vlim(self,val) :
self._bounds['vlim'] = val
return
@property
def name(self) :
'''Descriptive identifier for the line geometry.'''
if self._geomName is None : return ''
return self._geomName
@name.setter
def name(self, val) :
self._geomName = val
self.set_label( '' if val is None else val)
return
@property
def cname(self) :
'''Descriptive identifier for values indicated by color.'''
if self.valuesName is None : return ''
return self.valuesName
@cname.setter
def cname(self, val) :
self.valuesName = val
return
@property
def vertices(self) :
"""A 3 x N array of line vertices."""
v = self.vertexCoor
m = np.transpose(v)
return m
@property
def segmentcenters(self) :
"""A 3 x N array of line segment centers."""
segments = self._segments3d
totSegs = []
for line in segments :
npline = np.array(line)
for segIx in range(npline.shape[0]-1) :
center = np.add(npline[segIx],npline[segIx+1])/2
totSegs.append(center)
m = np.transpose(np.array(totSegs))
return m
@property
def segmentdirections(self) :
"""A N x 3 array of line segment direction vectors."""
# -----------------------------------------
segments = self._segments3d
if not self.isShredded :
segments = self._shred()
# -----------------------------------------
vfT = np.transpose(segments, (1,0,2) )
direction = vfT[1] - vfT[0]
return direction
@property
def segmentcolors(self) :
"""A N x 4 array of segment colors."""
return self._ledgecolors # update for v_1.2.0
@property
def bounds(self) :
"""
Dictionary of line geometric and value ranges.
Each dictionary value is a 2 float array of minimum
and maximum values of the line. Keys are:
'xlim' : x-coordinate
'ylim' : y-coordinate
'zlim' : z-coordinate
'r_xy' : radial distance from the z axis
'rorg' : radial distance from the origin
'vlim' : value functional assignments.
Values are assigned from the geometry and color mapping methods,
including cllipping.
"""
return self._bounds
@property
def cBar_ScalarMappable(self) :
"""matplotlib.cm.ScalarMappable object for line values. """
return self.stcBar_ScalarMappable()
[docs] def stcBar_ScalarMappable(self, stripe=None, bgcolor='k') :
"""
ScalarMappable object with striped color values.
Primarily useful for display of contour colorbars.
Parameters
----------
stripe : integer, optional, default: None
Number of color stripes. If None, no
stripes are formed.
bgcolor : color, optional, default: 'black'
Color between the stripes.
Returns
-------
self : Matplotlib.cm.ScalarMappable object
"""
# .......................................
def stripeCmap( cmap, numb, background='k', cname=None) :
if isinstance(cmap,str) :
cmap = cm.get_cmap(cmap)
numbSegs = 256
maxIndx = numbSegs-1
bgrnd = cm.colors.to_rgba_array(background)
cList = np.tile( bgrnd, (numbSegs,1))
fIndex = np.linspace(0.0,1.0,numb)
for i in fIndex :
index = int(i*(maxIndx))
cList[index] = cmap(i)
if index < maxIndx : cList[index+1] = cmap(i)
else : cList[index-2] = cmap(i)
if index > 0 : cList[index-1] = cmap(i)
else : cList[index+2] = cmap(i)
return colors.ListedColormap(cList)
# .......................................
vmin, vmax = self._bounds['vlim']
norm = colors.Normalize(vmin=vmin,vmax=vmax)
objMap = self.get_cmap()
if stripe is not None:
objMap = stripeCmap( objMap, stripe, bgcolor)
sm = cm.ScalarMappable(cmap=objMap, norm=norm)
sm.set_array([])
return sm
[docs] def shred(self, rez=0) :
"""
Recursively halve each line segment into separate lines.
Parameters
----------
rez : integer, optional, default: 0
Number of recursive subdivisions of a line segment.
Values range from 0 to 7.
Returns
-------
self : line object
"""
# .............................................................
def midSegmentVertices(numDiv,tailCoor, headCoor) :
# return added subdivision vertices, in order from
# tail to head.
minLen = 1.0/numDiv
divArr = np.linspace(minLen,1.0-minLen, numDiv-1)
delta = np.array(headCoor) - np.array(tailCoor)
vSet = [ tailCoor + np.multiply(delta,mult) for mult in divArr ]
vSet = np.array(vSet)
return vSet
# .............................................................
if not isinstance(rez, int) :
raise ValueError('Incorrect rez type, must ba of type int')
if (rez<0) or (rez>_MAXREZ) :
raise ValueError('Incorrest rez of {}, must be an int, 0 <= rez <= {}'.format(rez,_MAXREZ))
numDiv = int(2**rez)
# FutDev: following seems to work in cases checked but still not 100% confident.
if not self.isShredded :
colors = self._shred_color()
self._shred()
self._ledgecolors = colors # update for v_1.2.0
# DEVNOTE: colors now refer to segment colors, not line colors.
if rez==0: return self
else :
colors = self._shred_segm_color(numDiv)
numSeg = 0
for line in self.lineIndices : numSeg += len(line)-1
new_vertexCoor = np.copy(self.vertexCoor)
new_lineIndices = []
new_start_vert_index = len(self.vertexCoor)
for line in self.lineIndices :
single_lineIndex = []
numLineSegs = len(line)-1
lastLineVertIndex = numLineSegs
for i in range(numLineSegs ) :
tail_vert = self.vertexCoor[ line[i] ]
head_vert = self.vertexCoor[ line[i+1] ]
single_lineIndex.append( line[i] )
new_verts = midSegmentVertices(numDiv,tail_vert, head_vert)
new_vertexCoor = np.concatenate( (new_vertexCoor,new_verts),axis=0)
next_start_index = new_start_vert_index + len(new_verts)
temp = [*range(new_start_vert_index,next_start_index)]
single_lineIndex.extend( temp )
new_start_vert_index = next_start_index
if (i+1) == lastLineVertIndex : single_lineIndex.append( line[i+1] )
new_lineIndices.append(single_lineIndex)
self.vertexCoor = new_vertexCoor
self.lineIndices = new_lineIndices
self._shred()
self._ledgecolors = colors # update for v_1.2.0
return self
[docs] def map_cmap_from_sequence( self, cmap=None, sop=None, cname=None ) :
"""
Sequential color assignment of each line segment.
Parameters
----------
cmap : str or Colormap, optional, default: 'viridis'
A Colormap instance or registered colormap name.
sop : function object, optional, default: None
Function that takes one argument,
a N Numpy array of scalar values ranging 0 to 1.
The function returns a Numpy array of scalar values.
cname : str, optional, default: function object name
if not a lambda function, otherwise ''.
Identifier for the color values.
Returns
-------
self : line object
"""
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cm_colorMap = self.get_cmap()
self.valuesName = 'Sequence'
# -----------------------------------------
segments = self._segments3d
if not self.isShredded :
segments = self._shred()
# -----------------------------------------
i = np.linspace(0.0,1.0, len(segments) )
v = i
vmin, vmax = 0.0, 1.0
if sop is not None :
temp = _getFunctionName(sop,cname)
if temp is not None : self.valuesName = temp
v = np.array(sop(i))
vmin, vmax = v.min(), v.max()
norm = colors.Normalize(vmin=vmin,vmax=vmax)
self._ledgecolors = cm_colorMap(norm(v)) # update for v_1.2.0
RGB = cm_colorMap(norm(v))
self.set_color(RGB) # <<<<<<<<<<<<<<<<<<<<<<< debug
self.segValues = v
self._bounds['vlim'] = [ vmin, vmax ]
# .............................................
self._map_cmap_from_sequence_vert(cm_colorMap)
# .............................................
return self
[docs] def map_cmap_from_direction( self, direction=None, cmap=None, isAbs=False, cname=None ) :
"""
Assignment of line color based on line segment direction.
The dot product of line segment directions with the direction paramenter
is used to assign line segment colors from a colormap.
Parameters
----------
direction : array-like, optional, default: (1,0,1)
A xyz vector.
cmap : str or Colormap, optional, default: 'viridis'
A Colormap instance or registered colormap name.
isAbs : boolean, optional, default: False
If True, the absolute value of the dot product
will control the color.
cname : str, optional, default is None
Identifier for the color values.
Returns
-------
self : line object
"""
if direction is None : direction = _ILLUM
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cm_colorMap = self.get_cmap()
self.valuesName = cname
# -----------------------------------------
segments = self._segments3d
if not self.isShredded :
segments = self._shred()
# -----------------------------------------
segDir = self._get_segment_directions(segments)
unitDirection = np.divide( direction, np.linalg.norm(direction) )
v = np.dot(segDir,unitDirection)
if isAbs : v = -np.abs(v)
vmin, vmax = v.min(), v.max()
norm = colors.Normalize(vmin=v.min(),vmax=v.max())
self._ledgecolors = cm_colorMap(norm(v)) # update for v_1.2.0
RGB = cm_colorMap(norm(v))
self.set_color(RGB) # <<<<<<<<<<<<<<<<<<<<<<< debug
self.segValues = v
self._bounds['vlim'] = [ vmin, vmax ]
# .............................................
self._map_cmap_from_direction_vert( self._ledgecolors, self.segValues, self._bounds['vlim'] ) # update for v_1.2.0
# .............................................
return self
[docs] def map_cmap_from_datagrid(self, datagrid, cmap=None, cname=None ) :
"""
Line color assignment using a 2D datagrid.
Datagrid values are normalized in the range 0 to 1.
Parameters
----------
datagrid : 2D float array
cmap : str or Colormap, optional
A Colormap instance or registered colormap name.
If not assigned, the surface Colormap is used.
The colormap maps the datagrid values to colors.
cname : string, optional, default: None.
Descriptive identifier for values indicated by color.
Returns
-------
self : line object
"""
data = datagrid
dmax = np.amax(datagrid)
dmin = np.amin(datagrid)
delta = dmax-dmin
data = (datagrid - dmin)/delta
xd = np.linspace(-1, 1, data.shape[0] )
yd = np.linspace(-1, 1, data.shape[1] )
g = interpolate.interp2d(yd, xd, data, kind='cubic')
# -----------------------------------------
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cm_colorMap = self.get_cmap()
self.valuesName = cname
# -----------------------------------------
segments = self._segments3d
if not self.isShredded :
segments = self._shred()
# -----------------------------------------
segCoor = self._get_segment_centers(segments)
x,y,_ = np.transpose(segCoor)
# normalize xy plane -1 to 1
x = np.interp( x, (x.min(), x.max()), (-1.0, 1.0) )
y = np.interp( y, (y.min(), y.max()), (-1.0, 1.0) )
trial = [ g( x[i],y[i])[0] for i in range(len(x)) ]
v = np.array(trial)*delta + dmin
vmin, vmax = v.min(), v.max()
norm = colors.Normalize(vmin=v.min(),vmax=v.max())
self._ledgecolors = cm_colorMap(norm(v)) # update for v_1.2.0
RGB = cm_colorMap(norm(v))
self.set_color(RGB) # <<<<<<<<<<<<<<<<<<<<<<< debug
self.segValues = v
self._bounds['vlim'] = [ vmin, vmax ]
# .............................................
self._map_cmap_from_datagrid_vert(g, delta, dmin, cm_colorMap)
# .............................................
return self
[docs] def map_cmap_from_op( self, operation, cmap=None, cname=None ) :
"""
Functional assignment of a color from a color map.
Line segment center coordinates are used to calculate a scalar which
is then used to assign segment colors from a colormap.
Parameters
----------
operation : function object
Function that takes one argument,
a 3xN Numpy array of native coordinates.
The function returns a Numpy array of scalar values.
cmap : str or Colormap, optional, default: 'viridis'
A Colormap instance or registered colormap name.
cname : str, optional, default: function object name
if not a lambda function, otherwise ''.
Identifier for the color values.
Returns
-------
self : line object
"""
cname = _getFunctionName(operation,cname)
if cname is not None : self.valuesName = cname
if cmap is not None :
if isinstance(cmap,str) :
self.set_cmap(cm.get_cmap(cmap))
else :
self.set_cmap(cmap)
cm_colorMap = self.get_cmap()
# -----------------------------------------
segments = self._segments3d
if not self.isShredded :
segments = self._shred()
# -----------------------------------------
segCoor = self._get_segment_centers(segments)
# .....
xyz = np.transpose(segCoor)
v = np.array(operation(xyz))
vmin, vmax = v.min(), v.max()
norm = colors.Normalize(vmin=v.min(),vmax=v.max())
self._ledgecolors = cm_colorMap(norm(v)) # update for v_1.2.0
RGB = cm_colorMap(norm(v))
self.set_color(RGB) # <<<<<<<<<<<<<<<<<<<<<<< debug
self.segValues = v
self._bounds['vlim'] = [ vmin, vmax ]
# ========================================================================
self._map_cmap_from_op_vert(operation, cm_colorMap)
# ========================================================================
return self
[docs] def map_color_from_op( self, operation, rgb=True, cname=None ) :
"""
Assignment of line color from a function.
Line segment colors are assigned from a function of line
segment center coordinates.
Parameters
----------
operation : function object
Function that takes one argument,
a 3xN Numpy array of native coordinates.
The function returns a 3xN color value.
rgb : bool {True, False}, optional, default: True
By default, RGB color values are returned by the
operation function. If set False, the operation
returns HSV color values.
cname : str, optional, default: function object name
if not a lambda function, otherwise ''.
Identifier for the color values.
Returns
-------
self : line object
"""
cname = _getFunctionName(operation,cname)
if cname is not None : self.valuesName = cname
# -----------------------------------------
segments = self._segments3d
if not self.isShredded :
segments =self._shred()
# -----------------------------------------
segCoor = self._get_segment_centers(segments)
xyz = np.transpose(segCoor)
colors = np.array(operation(xyz))
colors = np.transpose(colors)
colors = np.clip(colors,0,1)
if rgb : RGB = colors
else : RGB = cm.colors.hsv_to_rgb(colors)
if RGB.shape[1] == 3 :
ones = np.ones(RGB.shape[0])[:,np.newaxis]
RGB = np.concatenate((RGB,ones),axis=1)
self._ledgecolors = RGB # update for v_1.2.0
self.set_color(RGB) # <<<<<<<<<<<<<<<<<<<<<<< debug
self.segValues = None
self._bounds['vlim'] = [ 0.0, 1.0 ]
# ========================================================================
self._map_color_from_op_vert( operation, rgb )
# ========================================================================
return self
[docs] def map_geom_from_op(self, operation, name=None) :
"""
Functional transformation of line vertex coordinates.
Parameters
----------
operation : function object
Coordinate mapping function that takes one
argument, a 3xN Numpy array of native coordinates.
The function returns a 3xN array of coordinates.
name : string, optional, default: None or op function name.
Descriptive identifier for the geometry.
Returns
-------
self : line object
"""
orgCoor = self.vertexCoor
xyz = np.transpose(orgCoor)
XYZ = np.array(operation(xyz))
v = np.transpose(XYZ)
self.vertexCoor = v
segments = self._get_segments( self.vertexCoor, self.segmIndices )
self.set_segments(segments)
self._set_geometric_bounds()
if name is None : name = self._geomName
name = _getFunctionName(operation,name)
if name is not None : self._geomName = name
return self
[docs] def map_to_plane(self, dist, **kargs) :
"""
Project line onto a mapping surface.
Parameters
----------
dist : number, optional, default : 0.0
distance from the origin to the surface,
along the direction vector.
direction : array-like, optional, default: [0,0,1]
A xyz vector normal to the intersection planes
for planar mapping or axial direction of
a cylinder for cylindrical mapping.
coor : integer or string indicating the type of mapping surface:
0, p, P, xyz,planar - planar (default)
1, c, C, cylinder,pplar,cylindrical - cylinder
2, s, S, sphere,spherical - sphere
3, x, X - y-z plane
4, y, Y - x-z plane
5, z, Z - x-y plane
Returns
-------
self : line object
"""
dirDft = { 'direction': [0,0,1.0] }
kargDft = { **_COOR_KWARGS, **dirDft }
KW = KWprocessor( kargDft, 'map_to_plane', **kargs)
coor = KW.getVal('coor')
direction = KW.getVal('direction')
if coor==3 : coor, direction = 0, [ 1,0,0 ]
if coor==4 : coor, direction = 0, [ 0,1,0 ]
if coor==5 : coor, direction = 0, [ 0,0,1 ]
# .......................................................
def project_plane(xyz) :
verts = xyz.T
unitDirection = np.divide( direction, np.linalg.norm(direction) )
dprod = np.dot( verts, unitDirection )
vdot = np.array([dprod]).T
dUnit = np.tile(unitDirection, (len(verts),1) )
vUnit = vdot*dUnit
dUnit = dist*dUnit
projVerts = verts + (dUnit - vUnit)
return projVerts.T
def project_cyl(xyz) :
verts = xyz.T
unitDirection = np.divide( direction, np.linalg.norm(direction) )
dmag = np.dot(verts,unitDirection)
uVect = np.tile(unitDirection,(len(dmag),1) )
dVect = np.array([dmag]).T
A = uVect*dVect
zeta = verts - A
zetanorm = np.linalg.norm(zeta,axis=1)
znorm = np.tile(zetanorm, (3,1))
unitZeta = np.divide( zeta, znorm.T)
# above may report nan warning, but following corrects these values.
unitZeta[ zetanorm <= 0.0, : ] = [0.0,0.0,0.0]
projVerts = A + dist*unitZeta
return projVerts.T
def project_sph(xyz) :
verts = xyz.T
dprod = np.linalg.norm(verts,axis=1)
tnorm = np.tile(dprod, (3,1) )
unitVerts = np.divide( verts, tnorm.T )
projVerts = dist * unitVerts
return projVerts.T
# .......................................................
project = project_plane
if coor == 1 : project = project_cyl
if coor == 2 : project = project_sph
return self.map_geom_from_op( project )
[docs] def clip(self, operation) :
"""
Remove segments from the line.
NOTE: all operations are in xyz coordinates.
Parameters
----------
operation : function object
Function that takes one argument, a 3xN Numpy array of coordinates.
The function returns N array of bool { True, False } indicating if
the segment, at the segment center coordinate, is to be
retained (True). Otherwise, the segment is removed from the
line (False).
Returns
-------
self : line object
"""
# -----------------------------------------
colors = self._ledgecolors # update for v_1.2.0
segments = self._segments3d
if self.restoredClip is None :
self.restoredClip = {
'indices' : self.segmIndices,
'colors' : colors,
'isShredded' : self.isShredded }
if not self.isShredded :
# .... note:colors modified, not created
colors = self._shred_color()
segments = self._shred()
self._ledgecolors = colors # update for v_1.2.0
else :
lenCol = len(self._ledgecolors) # update for v_1.2.0
lenSeg = len(segments)
if lenCol != lenSeg : colors = self._shred_color()
# -----------------------------------------
segCoor = self._get_segment_centers(segments)
xyz = np.transpose(segCoor)
shouldKeep = np.array(operation(xyz))
if np.any(shouldKeep) is None :
warnings.warn('WARNING: Clipping resulted in no lines found, line return unclipped.')
return self
clipIndices = np.array(self.segmIndices)[shouldKeep]
clipColors = colors[shouldKeep]
clipColors = np.array(clipColors)
clipCoor = self.vertexCoor[ np.array(clipIndices) ] # the lineset is shredded so can do this.
self.set_segments(clipCoor)
self.set_color(clipColors)
# FutDev: 'correct' definition should use nparrays for segmIndices not lists
# self.segmIndices = clipIndices ... bug fix, just a bandaid.
self.segmIndices = clipIndices.tolist()
return self
[docs] def clip_plane(self,dist,**kargs) :
"""
Remove segments from the line based on a clip surface.
Parameters
----------
dist : number, optional, default : 0.0
Distance from the origin to the clip surface,
along the direction vector.
direction : array-like, optional, default: [0,0,1]
A xyz vector normal to the intersection plane
for a planar clip surface or axial direction of
a cylinder for a cylindrical clip surface.
coor : integer or string indicating the type of clip surface:
0, p, P, xyz,planar - planar (default)
1, c, C, cylinder,pplar,cylindrical - cylinder
2, s, S, sphere,spherical - sphere
3, x, X - y-z plane
4, y, Y - x-z plane
5, z, Z - x-y plane
Returns
-------
self : line object
"""
# ...........................................................
def clip_op(xyz,dist,coor,direction) :
abc = np.transpose(xyz)
dotprod = _coor_dotprod(direction,coor,abc)
shouldclip = np.less(dotprod,np.full(len(dotprod),dist))
return shouldclip
# ...........................................................
dirDft = { 'direction': [0,0,1.0] }
kargDft = { **_COOR_KWARGS, **dirDft }
KW = KWprocessor( kargDft, 'clip_plane', **kargs)
coor = KW.getVal('coor')
direction = KW.getVal('direction')
if coor==3 : coor, direction = 0, [ 1,0,0 ]
if coor==4 : coor, direction = 0, [ 0,1,0 ]
if coor==5 : coor, direction = 0, [ 0,0,1 ]
return self.clip(lambda xyz : clip_op(xyz,dist,coor,direction) )
def _restore_from_clip(self) :
"""
Reset segment and colors prior to any clip operation.
(in development)
"""
if self.restoredClip is None :
warnings.warn('There are no clipped lines to restore.')
return self
self.isShredded = self.restoredClip['isShredded']
self.set_color( self.restoredClip['colors'] )
self.segmIndices = self.restoredClip['indices']
segmIndices = self.segmIndices
segments = self._get_segments( self.vertexCoor, segmIndices )
self.set_segments(segments)
self.restoredClip = None
return self
[docs] def set_line_alpha(self, alpha, constant=False) :
"""
Adjust the alpha value of the line segment colors.
Parameters
----------
alpha : scalar
Alpha is in the domain 0 to 1.
constant : bool { True, False }, optional, False
If False, segment color values are multiplied by
alpha. If True, all segment colors alpha channels
are assigned to alpha.
Returns
-------
self : line object
"""
if (alpha<0.0) or (alpha>1.0) :
raise ValueError('line alpha must be between 0 and 1, found {}'.format(alpha))
orig_colors = self._ledgecolors # update for v_1.2.0
if constant :
np.put_along_axis(orig_colors, np.array([[3]]), alpha, axis=1)
colorMap = orig_colors
else:
colorMap = np.multiply(orig_colors,np.array([1.0,1.0,1.0,alpha]))
self._ledgecolors = colorMap # update for v_1.2.0
return self
[docs] def shade(self, depth=0, direction=None, contrast=None) :
"""
Reduce line HSV color Value based on line segment direction.
The dot product of line segment with the illumination direction
is used to adjust face HSV color value.
Parameters
----------
depth : scalar, optional, default: 0
Minimum color value of shaded line segments.
Depth value ranges from 0 to 1.
direction : array-like, optional, default: (1,0,1)
A xyz vector pointing to the illumination source.
contrast : scalar, optional, default: 1
Shading contrast adjustment from low to high with a value
of 1 for linear variations with the line segment direction.
Contrast value ranges from 0.1 to 3.
Returns
-------
self : line object
"""
if direction is None : direction = _ILLUM
if np.any( (depth<0) | (depth>1) ) :
raise ValueError('depth values, {}, must be between 0 and 1'.format(depth))
if contrast is not None:
if (contrast<0.1 or contrast>3) :
raise ValueError('contrast must be between 0.1 and 3. , found {}'.format(contrast))
# -----------------------------------------
colors = self._ledgecolors # update for v_1.2.0
segments = self._segments3d
if not self.isShredded :
# .... note:colors modified, not created
colors = self._shred_color()
segments = self._shred()
self._ledgecolors = colors # update for v_1.2.0
else :
lenCol = len(self._ledgecolors) # update for v_1.2.0
lenSeg = len(segments)
if lenCol != lenSeg : colors = self._shred_color()
# -----------------------------------------
segDir = self._get_segment_directions(segments)
unitDirection = np.divide( direction, np.linalg.norm(direction) )
dtprod = np.dot(segDir,unitDirection)
# for Lines, effective shade is normal to line
d = 1- np.abs(dtprod)
if contrast is not None :
d = 0.5*( 1 + np.sign(dtprod)*np.power( np.abs(dtprod) , 1/contrast) )
orig_colors = colors
alphas = orig_colors[:,3:4]
fc_less_alpha = orig_colors[:,:3]
hsv_vals = cm.colors.rgb_to_hsv(fc_less_alpha)
# adjust color value with multiplier.
cvm = ((1-depth)*d + depth*np.ones(len(d)))[:, np.newaxis]
hsv_vals[:,2] = (cvm*hsv_vals[:,2:3])[:,0]
rgb_vals = cm.colors.hsv_to_rgb(hsv_vals)
shade_colors = np.concatenate((rgb_vals, alphas), axis=1)
self._ledgecolors = shade_colors # update for v_1.2.0
self.set_color(shade_colors)
return self
[docs] def fade(self,depth=0,elev=None,azim=None,ax=None) :
"""
Reduce line opacity based on line segment position relative
to the view orientation.
Parameters
----------
depth : scalar, optional, default: 0
Minimum opacity to 1 for linear line opacity from back
to front line segment position.
Depth value ranges from 0 to 1.
elev : scalar, optional, default: 30
Elevation of the axes view.
azim : scalar, optional, default: -60
Azimuth of the axes view.
ax: Matplotlib 3D axes.
If not None, elev and azim are taken from the ax.
Returns
-------
self : line object
"""
if ax is None :
if elev is None: elev = _DFT_VIEW[0]
if azim is None: azim = _DFT_VIEW[1]
else :
elev = ax.elev
azim = ax.azim
direction = elev_azim_2vector(elev, azim)
unitDirection = np.divide( direction, np.linalg.norm(direction) )
if np.any( (depth<0) or (depth>1) ) :
raise ValueError('depth values, {}, must be between 0 and 1'.format(depth))
# -----------------------------------------
#colors = self._edgecolors
colors = self.get_edgecolor() # update for v_1.2.0
segments = self._segments3d
if not self.isShredded :
# .... note:colors modified, not created
colors = self._shred_color()
segments = self._shred()
#self._edgecolors = colors
self.set_edgecolor(colors) # update for v_1.2.0
else :
#lenCol = len(self._edgecolors)
lenCol = len(colors) # update for v_1.2.0
lenSeg = len(segments)
if lenCol != lenSeg : colors = self._shred_color()
# -----------------------------------------
segCenters = self._get_segment_centers(segments)
if len(segCenters) == 1 : return self # fade not available for single segment line
dtprod = np.dot(segCenters,unitDirection)
dprange = np.amax(dtprod)-np.amin(dtprod)
norm_dp = (dtprod - np.amin(dtprod))/dprange
fadex = (1-depth)*norm_dp + depth
orig_colors = colors
fc_less_alpha = orig_colors[:,:3]
fadex = fadex*orig_colors[:,3] # update for v_1.2.0
faded_colors = np.concatenate((fc_less_alpha, fadex[:,np.newaxis] ), axis=1)
#self._edgecolors = faded_colors
self.set_edgecolor(faded_colors) # update for v_1.2.0
self.set_color(faded_colors) # update for v_1.2.0
return self
def _rectangulateFilledSurFace(self,rez,baseVcoor,baseFaceVertexIndices, midVectFunc) :
vCoor = baseVcoor.copy()
fvIndices = []
evIndices = []
#.....................................
def getVerticesLine(degree, leftIndex, rightIndex) :
def recurs_edgeCoor(degree, indexOrigin, isRight, leftIndex, rightIndex, Eindex) :
m = degree - 1
delta = int(np.power(2,m))
if isRight : delta = -delta
i = indexOrigin - delta
mid_coor = midVectFunc( vCoor[leftIndex], vCoor[rightIndex] )
currentCoorIndex = len(vCoor)
vCoor.append(mid_coor)
Eindex[i] = currentCoorIndex
if m <= 0 :
evIndices.append( [leftIndex, currentCoorIndex] )
evIndices.append( [currentCoorIndex, rightIndex] )
return
recurs_edgeCoor(m,i,False,leftIndex,currentCoorIndex, Eindex)
recurs_edgeCoor(m,i,True,currentCoorIndex,rightIndex, Eindex)
return
#.....................................
iOrigin = int(np.power(2,degree))
Eindex = [None]*(iOrigin+1)
Eindex[0] = leftIndex
Eindex[iOrigin] = rightIndex
if degree==0 :
evIndices.append( [leftIndex, rightIndex] )
return Eindex
recurs_edgeCoor(degree, iOrigin, False, leftIndex, rightIndex, Eindex)
return Eindex
#.....................................
def recurs_faceIndices(A,C,rez):
# similar to inner function of triangulateBase
if rez==0 :
abc = [ A[0], A[1], C[0], C[1] ]
fvIndices.append( abc )
return
# -----------------------------------------------------
J = int( (len(A)-1)/2 )
K = J + 1
# sudivide the two inner quadraterals .................
recurs_faceIndices( A[:K], C[J:], rez-1)
recurs_faceIndices( C[:K], A[J:], rez-1)
return
#.....................................
i_prev_2 = -1 # check for first face only.
for face in baseFaceVertexIndices :
if i_prev_2 == face[3] :
B = D[::-1]
else :
B = getVerticesLine(rez, face[3],face[0] )
i_prev_2 = face[2]
D = getVerticesLine(rez, face[1],face[2] )
recurs_faceIndices(B,D,rez)
vertexCoor = np.array(vCoor)
indexObj = { 'face': np.array(fvIndices) , 'edge': np.array(evIndices) }
return indexObj ,vertexCoor
def _get_filled_surface(self, addedverts, lrez=0, sclr=None, alpha=None, name=None) :
# Dual purpose algorithm for lines projected to a plane (True) or
# Parametric lines projected to another Parametric line.
#.....................................
def prepSegs() :
# devNote: even if lines are not shredded, each
# line could possibly have different colors.
# Hence, must construct surfaces from 'shredded'
# lines without actually shredding the line.
# The color array is then for each polygon.
# The following is actually somewhat a duplication
# of the _shred.
if not self.isShredded :
color = self._shred_color()
linesegs = self.lineIndices
indices = []
for line in linesegs :
numSegs = len(line) - 1
for seg in range (numSegs) :
t = [ line[seg],line[seg+1] ]
indices.append(t)
vertIx = np.array(indices)
else :
color = self._edgecolors # update for v_1.2.0
vertIx= np.array(self.segmIndices)
return vertIx, color
#.....................................
topSegInx,color = prepSegs()
stInx = len(self.vertexCoor)
botSegInx = topSegInx + stInx
topSegInx = np.flip(topSegInx,axis=1)
fcIndx = np.concatenate( (botSegInx,topSegInx), axis=1)
sverts = np.concatenate((self.vertexCoor,addedverts),axis=0)
scolor = np.array(color)
if lrez > 0 :
sverts = list(sverts)
indexObj ,vertexCoor = self._rectangulateFilledSurFace(lrez,sverts,fcIndx, self._midVectorFun)
sverts = vertexCoor
fcIndx = indexObj['face']
scolor = np.repeat(scolor,2**lrez,axis=0)
if sclr is not None : scolor = sclr
if alpha is not None:
surface = Surface3DCollection(sverts, fcIndx, name=name, color=scolor, alpha=alpha)
else :
surface = Surface3DCollection(sverts, fcIndx, name=name,color=scolor)
return surface
[docs] def get_filled_surface(self,**kargs) :
"""
Surface3DCollection from a line projected to a surface.
Parameters
----------
dist : number, optional, default : 0.0
distance from the origin to the projected surface,
along the direction vector.
direction : array-like, optional, default: [0,0,1]
A xyz vector normal to the intersection plane
for planar projections or axial direction of
a cylinder for cylindrical projections.
coor : integer or string indicating the type of projection surface:
0, p, P, xyz,planar - planar (default)
1, c, C, cylinder,pplar,cylindrical - cylinder
2, s, S, sphere,spherical - sphere
3, x, X - y-z plane
4, y, Y - x-z plane
5, z, Z - x-y plane
lrez : positive integer, optional, default : 0
projection recursive subdivisions.
color : Color, optional, default : line color
alpha : scalar, optional, default : 1
Face color alpha value, in the range 0 to 1.
name : string, optional, default: None.
Descriptive identifier for the geometry.
Returns
-------
Surface3DCollection object
"""
extDft = { 'direction':[0,0,1], 'dist':0.0, 'lrez':0, 'alpha':None, 'color':None, 'name':None }
kargDft = { **_COOR_KWARGS, **extDft }
KW = KWprocessor(kargDft,'get_filled_surface',**kargs)
drctn = KW.getVal('direction')
axis = KW.getVal('coor')
base = KW.getVal('dist')
lrez = KW.getVal('lrez')
alpha = KW.getVal('alpha')
color = KW.getVal('color')
name = KW.getVal('name')
if base is None : base = 0.0
if axis==3 : axis, drctn = 0, [ 1,0,0 ]
if axis==4 : axis, drctn = 0, [ 0,1,0 ]
if axis==5 : axis, drctn = 0, [ 0,0,1 ]
selfcopy = copy.copy(self)
mappedself = selfcopy.map_to_plane( base, direction=drctn, coor=axis)
surf = self.get_surface_to_line(mappedself, lrez=lrez, alpha=alpha, color=color, name=name)
return surf
[docs] def get_surface_to_line(self, line, **kargs) :
"""
Surface3DCollection from a line projected to another line.
Parameters
----------
line : line object with the same number of vertices
as the calling object.
lrez : positive integer, optional, default : 0
projection recursive subdivisions.
color : Color, optional, optional, default : line color
alpha : scalar, optional, default : None (1)
Face color alpha value, in the range 0 to 1.
name : string, optional, default: None.
Descriptive identifier for the geometry.
Returns
-------
Surface3DCollection object
"""
kargDft = { 'lrez':0, 'alpha':None, 'color':None, 'name':None }
KW = KWprocessor(kargDft,'get_surface_to_line',**kargs)
lrez = KW.getVal('lrez')
alpha = KW.getVal('alpha')
name = KW.getVal('name')
color= KW.getVal('color')
len1, len2 = len(self.vertexCoor), len(line.vertexCoor)
if len1 != len2 :
raise ValueError('Mismaatched line rez for Parametric Lines: {} != {}'.format(len1,len2))
addedverts = line.vertexCoor
return self._get_filled_surface(addedverts, lrez=lrez, sclr=color, alpha=alpha, name=name)
[docs]class ParametricLine(ColorLine3DCollection) :
def __init__( self, rez , operation=None, name=None, **kwargs ) :
"""
Create single line from a parametric function with domain [0,1]
Parameters
----------
rez : integer, optional, default: 0
Number of recursive bisection of line segments.
Rez values range from 0 to 7. If rez is less than zero,
negative value will be the number of segments.
operation : function object
Function that takes one argument, a Numpy float
array. The function domain is [0.0, 1,0] and returns
a 3 X N array of xyz coordinates.
name : str, optional, default: function object name
if not a lambda function, otherwise string ''.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.ColorLine3DCollection'.
"""
# -----------------------------------------
def default_line(numVerts) :
# xyz axes lines [-1,1]
lineVals = np.linspace(-1.0,1.0, int(numVerts) )
zeros = np.zeros( int(numVerts) )
Xaxis = np.asarray( [lineVals, zeros, zeros ] ).T
Yaxis = np.asarray( [zeros, lineVals, zeros ] ).T
Zaxis = np.asarray( [zeros, zeros, lineVals ] ).T
verts = np.vstack( (Xaxis,Yaxis,Zaxis ))
indices = [None]*3
index = 0
for i in range(0,3) :
indices[i] = [None]*(numVerts)
for j in range(numVerts) :
indices[i][j] = index
index += 1
return verts, indices
# -----------------------------------------
if not isinstance(rez,int) :
raise ValueError('ParametricLine: Incorrest rez of {}, must be an int'.format(rez))
if (rez>_MAXPLREZ) :
raise ValueError('ParametricLine:Incorrest rez of {}, must less than or equal to {}'.format(rez,_MAXPLREZ))
if rez<0 :
self.numSeg = -rez
# FutDev: should check if > 2**(_MAXPLREZ+1)
else :
self.numSeg = int(2**(rez+1))
numVerts = self.numSeg + 1
if operation is None :
verts, indices = default_line(numVerts)
else :
name = _getFunctionName(operation,name)
lineVals = np.linspace(0.0,1.0, self.numSeg+1)
verts = np.array(operation(lineVals)).T
verts = verts.tolist()
indices = [ [*range(len(verts))] ]
super().__init__(verts, indices, name, **kwargs)
self.rez = rez
def _getSlice(self, g, numVerts, domain, xVal=None, yVal=None ) :
if xVal is None and yVal is None : return None, None
numLines = 0
if xVal is not None : numLines += 1
if yVal is not None : numLines += 1
indices = [None]*( numLines )
numb_X_lines = 0
index = 0
verts = []
if xVal is not None :
numb_X_lines = 1
lineVertsPts = np.linspace(domain['ymin'], domain['ymax'], numVerts )
x_lineVertsPts = np.array([xVal] )
z_X_lines = g(x_lineVertsPts,lineVertsPts)
for iY in range(numVerts) :
verts.append( [ xVal , lineVertsPts[iY], z_X_lines[iY,0] ])
indices[0] = [None]*(numVerts)
for j in range(numVerts) :
indices[0][j] = index
index += 1
if yVal is not None :
lineVertsPts = np.linspace(domain['xmin'], domain['xmax'], numVerts )
y_lineVertsPts = np.array([yVal] )
z_Y_lines = g(lineVertsPts, y_lineVertsPts)
# the following, compensates for interp2d for datagrid
# which drops the first dimension, see map_xySlice_from_op ???
z_Y_lines = np.array([g(lineVertsPts, y_lineVertsPts)])
for iX in range(numVerts) :
verts.append( [ lineVertsPts[iX], yVal, z_Y_lines[0,iX] ])
indices[numb_X_lines] = [None]*(numVerts)
for j in range(numVerts) :
indices[numb_X_lines][j] = index
index += 1
verts = np.array(verts)
return verts, indices
def _getSliceSet(self, g, numVerts, domain, xsamples, ysamples ) :
if xsamples is None and ysamples is None : return None, None
numb_X_lines, numb_Y_lines = 0,0
if xsamples is not None : numb_X_lines = xsamples
if ysamples is not None : numb_Y_lines = ysamples
indices = [None]*( numb_X_lines + numb_Y_lines )
index = 0
verts = []
if xsamples is not None :
lineVertsPts = np.linspace(domain['ymin'], domain['ymax'], numVerts )
x_lineVertsPts = np.linspace(domain['xmin'], domain['xmax'], numb_X_lines )
if xsamples == 1 : x_lineVertsPts = np.asarray([0.0])
z_X_lines = g(x_lineVertsPts,lineVertsPts)
for iX in range(numb_X_lines) :
for iY in range(numVerts) :
verts.append( [ x_lineVertsPts[iX] , lineVertsPts[iY], z_X_lines[iY,iX] ])
for i in range(0,numb_X_lines) :
indices[i] = [None]*(numVerts)
for j in range(numVerts) :
indices[i][j] = index
index += 1
if ysamples is not None :
lineVertsPts = np.linspace(domain['xmin'], domain['xmax'], numVerts )
y_lineVertsPts = np.linspace(domain['ymin'], domain['ymax'], numb_Y_lines )
if ysamples == 1 : y_lineVertsPts = np.asarray([0.0])
z_Y_lines = g(lineVertsPts,y_lineVertsPts)
for iY in range(numb_Y_lines) :
for iX in range(numVerts) :
verts.append( [ lineVertsPts[iX], y_lineVertsPts[iY] ,z_Y_lines[iY,iX] ])
for i in range(numb_X_lines,numb_X_lines+numb_Y_lines) :
indices[i] = [None]*(numVerts)
for j in range(numVerts) :
indices[i][j] = index
index += 1
verts = np.array(verts)
return verts, indices
def _map_xySlice(self, g, **kargs ) :
xplane, yplane, xsamples, ysamples = None,None,None,None
if 'xplane' in kargs : xplane = kargs['xplane']
if 'yplane' in kargs : yplane = kargs['yplane']
if 'xset' in kargs : xsamples = kargs['xset']
if 'yset' in kargs : ysamples = kargs['yset']
planeNormal = None # both x and y slice planes are defined.
if (xplane is None) and (xsamples is None) :
planeNormal = [0,1,0] # must be in the x-z plane
if (yplane is None) and (ysamples is None) :
planeNormal = [1,0,0] # must be in the y-z plane
domain = {'xmin':-1, 'xmax':1, 'ymin':-1, 'ymax':1}
def checkArr( A ) :
if isinstance(A,(list,tuple)) : return A[0], A[1]
else: return -A,A
if 'xlim' in kargs :
domain['xmin'], domain['xmax'] = checkArr(kargs['xlim'])
if 'ylim' in kargs :
domain['ymin'], domain['ymax'] = checkArr(kargs['ylim'])
# -----------------------------------------
if ( (xsamples is None and ysamples is None) and
(xplane is None and yplane is None) ) : return self
numVerts = self.numSeg + 1
vertSingle,indicesSingle = self._getSlice(g,numVerts,domain,xplane, yplane)
vertsSet,indicesSet = self._getSliceSet(g,numVerts,domain,xsamples, ysamples)
if vertSingle is not None :
verts,indices = vertSingle,indicesSingle
if vertsSet is not None :
startIndex = len(verts)
verts = np.concatenate( (verts,vertsSet), axis=0)
shiftedIndicesSet = [ [ i+startIndex for i in line ] for line in indicesSet]
numbIndiceSingle = len(indicesSingle)
numbIndicesSet = len(shiftedIndicesSet)
temp = [None]*(numbIndiceSingle + numbIndicesSet)
for i in range(numbIndiceSingle) : temp[i] = indicesSingle[i]
for i in range(numbIndicesSet) : temp[i+numbIndiceSingle] = shiftedIndicesSet[i]
indices = temp
else :
verts,indices = vertsSet,indicesSet
self.lineIndices = indices
self.vertexCoor = verts
self.segmIndices = indices
segs = self._get_segments( verts, indices )
self.set_segments(segs)
if len(self.get_color() ) == 1 :
c = self.get_color()[0]
self.set_color( [c]*len(indices) )
self._set_geometric_bounds()
self._planeNormal = planeNormal
return self
[docs] def map_xySlice_from_op(self, operation, **kargs ) :
"""
Using a functional operation, lines of z=f(x) or z=f(y) in domain -1 to 1.
Parameters
----------
operation : function object
Function that takes one argument, a Numpy float
3 X N array of xyz coordinates in the domain.
The return value is a 3 X N array of xyz coordinates.
xplane : contour at a constant value of x=xplane
yplane : contour at a constant value of y=xplane
xset : number of evenly spaced contours, y-axis normal plane
yset : number of evenly spaced contours, x-axis normal plane
xlim : X domain of operation ( default: [-1,1] )
ylim : Y domain of operation ( default: [-1,1] )
Returns
-------
self : ParametricLine object
"""
def zValue_generator(xArr, yArr, op) :
xlen = len(xArr)
ylen = len(yArr)
xcol = np.tile( [1.0,0.0,0.0], (xlen,1))
xVal = np.multiply(xcol.T,xArr).T
xVal3 = np.tile(xVal, (ylen,1))
ycol = np.tile( [0.0,1.0,0.0] , (ylen,1))
yVal = np.multiply(ycol.T,yArr).T
yVal3 = np.repeat(yVal,xlen,axis=0)
pts = xVal3 + yVal3
pts = np.reshape(pts,(xlen*ylen,3))
x,y,z = op(pts.T)
z = np.reshape(z,(ylen,xlen))
# the following, compensates for interp2d for datagrid
# which drops the first dimension, see _getSlice
if ylen == 1 : z = z[0]
return z
def g(xArr, yArr) :
return zValue_generator(xArr, yArr, operation)
# -----------------------------------------
name = _getFunctionName(operation,self._geomName)
if name is not None : self._geomName = name
return self._map_xySlice( g, **kargs )
[docs] def map_xySlice_from_datagrid(self, datagrid, scale=None, **kargs ) :
"""
Using a datagrid, lines of z=f(x) or z=f(y) in domain [-1, 1]
Parameters
----------
datagrid : 2D float array
scale: multiplier for the Z direction.
xplane : contour at a constant value of x=xplane
yplane : contour at a constant value of y=xplane
xset : number of evenly spaced contours, y-axis normal plane
yset : number of evenly spaced contours, x-axis normal plane
Returns
-------
self : ParametricLine object
"""
if 'xlim' in kargs or 'ylim' in kargs :
warnings.warn('xlim and ylim args for Datagrid xySlice not available.')
return self
data = datagrid
dmax = np.amax(datagrid)
dmin = np.amin(datagrid)
delta = dmax-dmin
data = (datagrid - dmin)/delta
xd = np.linspace(-1, 1, data.shape[0] )
yd = np.linspace(-1, 1, data.shape[1] )
g = interpolate.interp2d(yd, xd, data, kind='cubic')
# -----------------------------------------
self._map_xySlice( g, **kargs )
if scale is not None :
self.vertexCoor = np.multiply(self.vertexCoor,[1,1,scale])
indices = self.segmIndices
verts= self.vertexCoor
segs = self._get_segments( verts, indices )
self.set_segments(segs)
return self
[docs] def scale_dataframe(self,X,Y,Z) :
"""
Scaling the line geometry based on a datagrid.
Parameters
----------
X, Y, Z : N x M arrays
Minimum and maximum values of the arrays are used to scale
and translate the surface from an intial domain of
[ (-1,1), (-1,1), (0,1) ]
Returns
-------
self : ParametricLine object
"""
return _scale_dataframe(self,X,Y,Z)
[docs]class SegmentLine(ColorLine3DCollection) :
def __init__( self, vertexCoor , name=None, **kwargs ) :
"""
A single line of squential segments from an array of vertices.
Parameters
----------
vertexCoor : V x 3 numerical array-like.
An array of V number of xyz vertex coordinates.
name : string, optional, default: None.
Descriptive identifier.
Other Parameters
----------------
**kwargs
All other parameters are passed on to 's3dlib.surface.ColorLine3DCollection'.
"""
line = vertexCoor # <--- super will make this to a numpy float array
indices = [ [*range(len(line))] ]
super().__init__(line, indices, name, **kwargs)
# =================================================================================+
# DevNote : Although Surface, Vector and Line collection classes have similar |
# methods (eg. map..., name, clip ..), differences among these classes |
# are sufficiently different to not warrant a class that these three |
# classes could use with multiple inheritence. Instead, the following |
# methods are common among these three classes. |
# FutDev may be useful to define such a class, but MROI. |
# |
# Overall ColorLine3DCollection class needs a architectural change. |
# Possibly a list of ColorSegment3DCollection objects which are |
# contained within each ColorLine3DCollectio object. Each list |
# is a basically a line. Analogous to Surface collections which |
# which are a collection of polyhedrons. This would be an advantage |
# for writing and reading to files. |
# |
# Becoming a priority: Need to have separate object properties for |
# face colors and edge colors, and don't use the Matplotlib properties |
# for these in Poly3DCollection, Line3DCollection. Too much 'magic' |
# happening which is out of control. Then use set_color |
# when needed. This applies to both surface and line objects. |
# This also would allow constructing classes independent of |
# Poly3DCollection, Line3DCollection for future developement, |
# allowing an interface to matplotlib to be constructed instead |
# inheritance. |
# =================================================================================+
class KWprocessor():
""" Access values from defaults or kwargs. """
def __init__(self,valDict, methodStr=None, checkInputKeys=True, **kargs) :
"""
Validate input of a method's kwargs. Internal use ONLY.
Parameters
----------
valDict : dictionary of argument names, defaults and valid values.
key - keyword argument name
value - single default value, or list of 2D lists to be validated:
If a list of list, the structure is:
value[0] - ['default value'] <- single value
value[1] - [ a, b, c, ..... ] <- list of a list of valid values
methodStr : string, optional, default: None
Identifier of the method called where a validation error occured
which is shown in the error message.
checkInputKeys : bool, default: True
If True, check if all kargs are in the valDict.
Unknown keys will produce a ValueError.
If False, ignore kargs not in valDict.
kargs :
The dictionary object containing the input keyword:values to be extracted.
"""
# FutDev: could use a dictionary for value[1] to validate ranges and types.
self.valDict = valDict
self.kargs = kargs
self.method = methodStr
if checkInputKeys :
if len(kargs) == 0 : return
if len(valDict) == 0 : return
errMsg = "Keyword, {}, is not recognized. \n List of possible keys are: {}"
if methodStr is not None : errMsg = '['+methodStr+'] - ' + errMsg
for key in kargs.keys() :
if not key in valDict.keys() :
raise ValueError( errMsg.format(key, list(valDict.keys()) ) )
return
def _defVal(self,v) :
value = self.valDict[v]
if isinstance(value,list) :
firstVal = value[0]
return firstVal[0] if isinstance(firstVal,(list,tuple)) else value
return value
def _valueKeywordError(self,keyStr,usrVal) :
dicVal = self.valDict[keyStr]
errList = dicVal[1]
for i in range(2,len(dicVal) ) : errList += dicVal[i]
errA = "Keyword, {}, value of {} is not recognized.\n List of possible values are: {}"
errB = "Keyword, {}, value of '{}' is not recognized.\n List of possible values are: {}"
errStr = errB if isinstance(usrVal,str) else errA
if self.method is not None : errStr = '['+self.method+'] - ' + errStr
return errStr.format(keyStr,str(usrVal),errList)
def getVal(self,keyStr,checkValidity=True) :
"""
Return the named argument value.
Parameters
----------
keyStr : string. Keyword
checkValidity : bool, optional, default: True
If True, karg value is checked to be valid
among the list of valid values.
Raises
------
ValueError
if invalid keyword assignment.
Returns
-------
value of keyword, or the default value if the keyStr is not found in kargs.
"""
def listOfList(v) :
if not isinstance(v,(list,tuple)) : return False
return isinstance(v[0],(list,tuple))
kargs = self.kargs
if len(kargs) == 0 : return self._defVal(keyStr)
if keyStr not in self.valDict : # only for development purposes.
print("################ CODE ERROR KW001 #####################")
return None
if keyStr in kargs :
dicVal = self.valDict[keyStr]
usrVal = kargs[keyStr]
if listOfList(dicVal) :
for vLst in dicVal :
if usrVal in set(vLst) : return vLst[0] #<< value is the first.
if not checkValidity : return usrVal
else :
raise ValueError( self._valueKeywordError(keyStr,usrVal))
return usrVal
return self._defVal(keyStr)
def filter(self):
"""
Removes defaults from kargs.
Use implies that more than default kargs are to be used,
thus 'checkInputKeys' should be set to False at __init__.
"""
for key in self.valDict:
try : del self.kargs[ key ]
except: pass
return self.kargs
# FutDev: include transform
def _axisUSF(axes) :
# axis unit scale factor
xlim = axes.get_xlim()
ylim = axes.get_ylim()
zlim = axes.get_zlim()
xR = xlim[1]-xlim[0]
yR = ylim[1]-ylim[0]
zR = zlim[1]-zlim[0]
sf = np.array([ yR*zR, xR*zR, xR*yR ])
return sf/np.linalg.norm(sf)
def _coor_dotprod(direction,cIndex,refCoor) :
# refCoor may be vertex(contours) or face(clip) coordinates.
# determines the distance from the refCoor points to either:
# 0. a plane, (planar coordinates)
# 1. a line, (cylindrical coordinates)
# 2. a point, (spherical coordinates)
# used for surface contours (vertices), surface clipping (faces)
# and line clipping (vertices).
verts = refCoor
if cIndex == 1 : # ..................... cylindrical
if direction is None :
dotprod = np.linalg.norm(verts[:,:2],axis=1) # z-axis default
else :
unitDirection = np.divide( direction, np.linalg.norm(direction) )
dmag = np.dot(verts,unitDirection)
uVect = np.tile(unitDirection,(len(dmag),1) )
dVect = np.array([dmag]).T
zeta = verts - uVect*dVect
dotprod = np.linalg.norm(zeta,axis=1)
elif cIndex == 2 : # ..................... spherical
dotprod = np.linalg.norm(verts,axis=1)
else : # ..................... planar
if direction is None : direction = [0,0,1.0] # xy-plane default
unitDirection = np.divide( direction, np.linalg.norm(direction) )
dotprod = np.dot(verts,unitDirection)
return dotprod
def _getFunctionName(op,xName=None) :
if xName is not None : return xName
# assign function name if not a lambda function
fobar = lambda x : x
if op.__name__ != fobar.__name__ :
name = op.__name__
else :
name = None
return name
def _set_geomBounds(v,bounds) :
"""Set _bounds dictionary values from vertex coordinates."""
# v = self.vertexCoor
xlim,ylim,zlim = np.transpose( [np.amin(v ,axis=0), np.amax(v ,axis=0) ] )
bounds['xlim'] = xlim
bounds['ylim'] = ylim
bounds['zlim'] = zlim
rmax_xy = np.amax( np.linalg.norm(v[:,0:2], axis=1) )
rmin_xy = np.amin( np.linalg.norm(v[:,0:2], axis=1) )
bounds['r_xy'] = [rmin_xy, rmax_xy ]
rmax_xyz = np.amax( np.linalg.norm(v, axis=1) )
rmin_xyz = np.amin( np.linalg.norm(v, axis=1) )
bounds['rorg'] = [rmin_xyz, rmax_xyz ]
bounds['rlim'] = [rmax_xy, rmax_xyz ] # .. vestigal
return
def _scale_dataframe(obj,X,Y,Z) :
"""
Linear scale and translate the geometry.
Used for scaling a geometry based on a datagrid.
Parameters
----------
X, Y, Z : 2d arrays
Minimum and maximum values of the arrays are used to scale
and translate the surface from an intial domain of
[ (-1,1), (-1,1), (0,1) ]
Returns
-------
obj : first calling argument, object
"""
Xmin, Xmax = np.amin(X), np.amax(X)
Ymin, Ymax = np.amin(Y), np.amax(Y)
Zmin, Zmax = np.amin(Z), np.amax(Z)
Xscale = (Xmax-Xmin)/2
X0 = (Xmax+Xmin)/2
Yscale = (Ymax-Ymin)/2
Y0 = (Ymax+Ymin)/2
Zscale = Zmax-Zmin
Z0 = Zmin
obj.transform(scale=[Xscale,Yscale,Zscale],translate=[X0,Y0,Z0])
return obj
def _transformAxis(selfProp, lenmult=1.0, width=1.5, colors=None, negaxis=True) :
"""
Line3DCollection of a transformed coordinate axis.
just take method from surface version 1.0 and allow use for all
other transformable objects.
Called from the interface method get_transformAxis
Returns
-------
Line3DCollection
"""
self_scale = selfProp['scale']
self_rotation = selfProp['rotation']
self_translation = selfProp['translation']
scale = np.multiply(lenmult,self_scale)
dDelta = np.multiply(np.identity(3).tolist(),scale).tolist()
npdot = np.dot(dDelta,self_rotation)
npdot2 = np.expand_dims(npdot,axis=1)
if negaxis :
lines = np.concatenate( (npdot2,-npdot2),axis=1)
else :
lines = np.insert(npdot2,0,[0,0,0],axis=1)
lines = np.add(lines,self_translation)
if colors is None :
colorMap = ['r','g','b']
else :
colorMap = colors
axisCol = Line3DCollection(lines, colors=colorMap, linewidths=width)
return axisCol
def _normLength( N, exp, fac, limitarea=None) :
"""
Default face normal vector length.
Length proportional to the mean triangulated face edge length.
"""
if limitarea is None : limitarea = 4*np.pi
# aratio - area ratio as a function of the number of faces, N.
# heuristic values, exp and fac, are adjustment for small rez surfaces,
# with aratio limit to 1 for large N. (continuous smooth surface)
aratio = 1.0 - fac*np.power(N,-exp)
sc = 0.5*(1.520) # scale of normalized edge length
nLen = sc*np.sqrt(aratio*limitarea/N)
return nLen
# =================================================================================+
# Utility Functions: |
# Useful for construction of 3D Matplotlib figures composed of S3Dlib objects. |
# =================================================================================+
[docs]def setupAxis(axes, **kargs ) :
"""
Add XYZ origin coordinate axis Vector3DCollection to the axes.
Parameters
----------
axes : Matplotlib 3D axes object.
length : coordinate axis length (single or 3-value list), default: 1.5
color : axis color (single or 3-value list), default: 'black'
offset : axis vector offset from the origin (single or 3-value list), default: 0
labels : axis label (single or 3-value list), default: ['X','Y','Z']
pad : label padding from the vector axis head (single or 3-value list), default: 0.15
width : axis line widths, default: 2
negaxis : show negative axis line, default: False
alr : axis length ratio, head size to axis length, default: 0.2
Returns
-------
axes
"""
kargDft = { 'length': 1.5, 'color':'k', 'labels':['X','Y','Z'],
'offset':0.0, 'pad': 0.15, 'width':2, 'alr':0.2, 'negaxis':False}
KW = KWprocessor(kargDft,'setupAxis',**kargs)
length = KW.getVal('length')
color = KW.getVal('color')
offset = KW.getVal('offset')
labels = KW.getVal('labels')
pad = KW.getVal('pad')
width = KW.getVal('width')
negaxis = KW.getVal('negaxis')
alr = KW.getVal('alr')
# Note: not checking array length==3, or a valid color.
if color is None : color='k'
if not colors.is_color_like(color) : # check if arg is a color array, not an array of colors
color = np.array ( [colors.to_rgba(color[0]), colors.to_rgba(color[1]), colors.to_rgba(color[2]) ] )
else :
color = np.array ( [colors.to_rgba(color), colors.to_rgba(color), colors.to_rgba(color)] )
if not isinstance(length,(list,tuple,np.ndarray)) : length = [length,length,length]
length=np.array(length)
if not isinstance(offset,(list,tuple,np.ndarray)) : offset = [offset,offset,offset]
offset = np.array(offset)
if labels is None : labels = ['X','Y','Z']
if not isinstance(labels,(list,tuple,np.ndarray)) : labels = [labels,labels,labels]
if not isinstance(pad,(list,tuple,np.ndarray)) : pad = [pad,pad,pad]
pad = np.array(pad)
"""Set XYZ origin coordinate axis arrows."""
seglength = length-offset
endCoor = np.multiply(np.identity(3),seglength)
strCoor = np.multiply(np.identity(3),offset)
xyz = np.transpose( strCoor )
uvw = np.transpose( endCoor )
vcf = Vector3DCollection(xyz,uvw,alr=alr, color=color[0], linewidth=width)
vcf._set_vectColor(color)
axes.add_collection3d( vcf )
if negaxis :
endCoor = np.multiply(np.identity(3), -1.0*length)
strCoor = np.multiply(np.identity(3), -1.0*offset)
xyz = np.transpose( strCoor )
uvw = np.transpose( endCoor )
vertexCoor = np.concatenate( (endCoor,strCoor),axis=0 )
segmIndices = [ [0,3],[1,4],[2,5] ]
line = ColorLine3DCollection(vertexCoor,segmIndices, color=color, linewidth=width)
axes.add_collection3d( line )
pos = length + pad
axes.text(pos[0],0,0,labels[0], color = color[0], fontweight='bold', fontsize='large',
horizontalalignment='center', verticalalignment='center')
axes.text(0,pos[1],0,labels[1], color = color[1], fontweight='bold', fontsize='large',
horizontalalignment='center', verticalalignment='center')
axes.text(0,0,pos[2],labels[2], color = color[2], fontweight='bold', fontsize='large',
horizontalalignment='center', verticalalignment='center')
return axes
[docs]def standardAxis( axes, **kargs ) :
"""
Add XYZ origin coordinate axis Vector3DCollection and
set the axis viewing angle, projection type and clear planes.
Parameters
----------
axes : Matplotlib 3D axes object.
kargs: identical to those for setupAxis function.
Returns
-------
axes
"""
setupAxis(axes, **kargs)
minmax = (-1,1)
axes.set(xlim=minmax, ylim=minmax, zlim=minmax)
axes.view_init(30, 30)
axes.set_proj_type('ortho')
axes.set_axis_off()
return axes
[docs]def auto_scale(axes,*obj3d,**kargs) :
"""
Scale Axes3d using object bounds.
Parameters
----------
axes : Matplotlib 3D axes.
obj3d : positional arguments of s3dlib objects
uscale : scaling in the range (0.5,2.5)
if neither uscale or rscale are defined -
. Each axes ranges are independently
. set by sizes of input objects.
if both uscale and rscale are defined -
. uscale will be ignored.
if defined using an int or float -
. Each axes ranges have the same scalar
. length, set by the values of the objects.
. Out of range values will set the scaling
. to the minimum or maximum.
else - scaling is one.
rscale : scaling in the range (0.5,2.5)
if neither uscale or rscale are defined -
. Each axes ranges are independently
. set by sizes of input objects.
if both uscale and rscale are defined -
. uscale will be ignored.
if defined using an int or float -
. Each axes ranges have the same minimum
. and maximum, set by the values of the
. objects radial extent from the origin.
. Out of range values will set the scaling
. to the maximum.
else - scaling is one.
Returns
-------
axes, which has been scaled
"""
X,Y,Z,R = [],[],[],[]
for obj in obj3d :
X.append(obj.bounds['xlim'][0])
X.append(obj.bounds['xlim'][1])
Y.append(obj.bounds['ylim'][0])
Y.append(obj.bounds['ylim'][1])
Z.append(obj.bounds['zlim'][0])
Z.append(obj.bounds['zlim'][1])
R.append(obj.bounds['rorg'][1])
if 'uscale' in kargs :
sc = kargs['uscale']
if not isinstance(sc,(int,float)): sc=1.0
if sc > 2.0 : sc = 2.0
if sc < 0.5 : sc = 0.5
AxLen = lambda a : max(a)-min(a)
AxMpt = lambda a : 0.5*( max(a) + min(a) )
maxAlen = sc*max( [AxLen(X), AxLen(Y), AxLen(Z)] )
X = [ (AxMpt(X)-0.5*maxAlen), (AxMpt(X)+0.5*maxAlen) ]
Y = [ (AxMpt(Y)-0.5*maxAlen), (AxMpt(Y)+0.5*maxAlen) ]
Z = [ (AxMpt(Z)-0.5*maxAlen), (AxMpt(Z)+0.5*maxAlen) ]
if 'rscale' in kargs :
sc = kargs['rscale']
if not isinstance(sc,(int,float)): sc=1.0
if sc > 2.0 : sc = 2.0
if sc < 0.5 : sc = 0.5
rmax = sc*max(R)
minmax = (-rmax,rmax)
axes.set(xlim=minmax, ylim=minmax, zlim=minmax )
return axes
axes.auto_scale_xyz(X,Y,Z)
return axes
[docs]def eulerRot(theta, phi, psi=0, useXconv=True, inrad=False) :
"""Rotation matrix from Euler angles, X or Y convention."""
# X-convention (zxz') and Y-convention (zyz') rotations.
tht_Rad = theta
phi_Rad = phi
psi_Rad = psi
if not inrad :
tht_Rad = theta*np.pi /180.0
phi_Rad = phi*np.pi /180.0
psi_Rad = psi*np.pi /180.0
cT,sT = np.cos(tht_Rad), np.sin(tht_Rad)
cP,sP = np.cos(phi_Rad), np.sin(phi_Rad)
cS,sS = np.cos(psi_Rad), np.sin(psi_Rad)
T_t = [ [ cT, sT, 0 ], [-sT, cT, 0 ], [ 0 , 0 , 1 ] ]
T_X = [ [ 1 , 0 , 0 ], [ 0 , cP, sP], [ 0 ,-sP, cP] ]
T_Y = [ [ cP, 0 ,-sP], [ 0 , 1 , 0 ], [ sP, 0 , cP] ]
T_s = [ [ cS, sS, 0 ], [-sS, cS, 0 ], [ 0 , 0 , 1 ] ]
if useXconv : T_p = T_X
else: T_p = T_Y
M = np.matmul( T_p, T_t )
if psi != 0 : M = np.matmul( T_s, M )
return M
[docs]def axisRot(alpha, direction, inrad=False) :
"""Rotation matrix from rotation about a directional axis"""
# see Wikipedia: Quaternions and spatial rotation
alp_Rad = -alpha
if not inrad :
alp_Rad = -alpha*np.pi /180.0
u = np.divide( direction, np.linalg.norm(direction) )
sT = np.sin(alp_Rad/2.0)
qR = np.cos(alp_Rad/2.0)
qI, qJ, qK = sT*u[0], sT*u[1], sT*u[2]
M = [
[ (1-2*(qJ*qJ + qK*qK)), 2*(qI*qJ - qK*qR) , 2*(qI*qK + qJ*qR) ] ,
[ 2*(qI*qJ + qK*qR) , (1-2*(qI*qI + qK*qK)), 2*(qJ*qK - qI*qR) ] ,
[ 2*(qI*qK - qJ*qR) , 2*(qJ*qK + qI*qR) , (1-2*(qI*qI + qJ*qJ)) ]
]
return M
[docs]def vectRot(xDirection, in_planeDirection) :
"""Rotations matrix from two vectors """
ipDir = np.array(in_planeDirection)
xDir = np.array(xDirection)
zDir = np.cross(xDir,ipDir)
yDir = np.cross(zDir,xDir)
vect = np.array( [xDir,yDir,zDir] )
vnorm = np.linalg.norm(vect,axis=1)
unitVects = np.divide(vect.T,vnorm).T
return unitVects
[docs]def rtv(direction,elev=None,azim=None) :
"""Transform direction vector 'relative to view' axis viewing angles."""
if elev is None : elev = _DFT_VIEW[0]
if azim is None : azim = _DFT_VIEW[1]
azim_Rad = azim*np.pi /180.0
elev_Rad = elev*np.pi /180.0
cA,sA = np.cos(azim_Rad), np.sin(azim_Rad)
cE,sE = np.cos(elev_Rad), np.sin(elev_Rad)
TM = [ [ cE*cA, cE*sA, sE ], [ -sA, cA, 0.0 ], [ -sE*cA, -sE*sA, cE ] ]
reldir = np.dot(direction,TM)
return reldir
[docs]def elev_azim_2vector(elev, azim, inrad=False) :
"""Unit vector in the terms of angular position."""
tht_Rad = azim
phi_Rad = np.pi/2 - elev
if not inrad :
tht_Rad = azim*np.pi /180.0
phi_Rad = (90.0-elev)*np.pi /180.0
k = np.cos(phi_Rad)
r = np.sin(phi_Rad)
i = r*np.cos(tht_Rad)
j = r*np.sin(tht_Rad)
return np.array( [i,j,k] )
[docs]def frame_to_value(x, *valList, **kargs) :
"""
Y = f(x) from a collection of linear segment functions f(x)
Parameters
----------
x : float value between 0 and 1.
valList : list of values, evenly spaced along x. For multiple lists, all must be the same length.
stepped : boolean, default is False.
If True, function evaluate to a constant of the first value between
valList values. If False, function is linear between valList values.
For True, number of linear segments is the valList length.
For False, number of linear segments is one less than valList length.
Returns
-------
List of values of size valList arguments number.
"""
# 0 < x < 1
def vFromList(f,A) :
nSeg = len(A) - 1
i = math.floor(f*nSeg)
if i>=nSeg : return A[nSeg]
x = (f*nSeg) % 1
val = A[i] + x*( A[i+1]-A[i] )
return val
def vFromStep(f,A) :
nSeg = len(A)
i = math.floor(f*nSeg)
if i>=nSeg : return A[nSeg-1]
if i<0 : return A[0]
val = A[i]
return val
stepped = False
if 'step' in kargs : stepped = kargs['step']
vals = []
for vList in valList :
if stepped : vals.append( vFromStep(x,vList) )
else : vals.append( vFromList(x,vList) )
return vals
[docs]def density_function(vals, bins=(10,10), scale=True, xbnd=True, ybnd=True, kind='linear') :
"""
2D histogram density, Data_density = f(x,y) from a x,y dataset.
Parameters
----------
vals : 2 X N array of x,y data values.
bins : [int,int], the number of bins for each dimension.
scale : bool or number, default: True
If True, normalize values to the average.
If False, scale is 1, the density value of data per
bin area. When numerical, multiplier for
the density value.
xbnd : bool, default : True
If True, x values are extended for density of zero.
For cyclic mapping, this should be set to False.
ybnd : bool, default : True
If True, y values are extended for density of zero.
For cyclic mapping, this should be set to False.
kind : {‘linear’, ‘cubic’, ‘quintic’}, default: 'linear'
The kind of spline interpolation to use.
Returns
-------
A function f(x,y)
"""
vals = np.array(vals)
bins = np.array(bins)
Zvals,xedges,yedges = np.histogram2d(*vals,bins=bins)
X,Y = 0,1
minval = np.amin(vals,axis=1)
maxval = np.amax(vals,axis=1)
datarng = (maxval-minval)
area = datarng[0]*datarng[1]
delta = datarng/bins
endval = maxval + delta/2
if xbnd :
zeros = np.zeros( (Zvals.shape[X]+2,Zvals.shape[Y]) )
zeros[1:bins[X]+1,:] = Zvals
Zvals = zeros
xedges = np.array([*(xedges-delta[X]/2),endval[X]])
else :
xedges = xedges+delta[X]/2
xedges = xedges[:len(xedges)-1]
if ybnd :
zeros = np.zeros( (Zvals.shape[X],Zvals.shape[Y]+2) )
zeros[:,1:bins[Y]+1] = Zvals
Zvals = zeros
yedges = np.array([*(yedges-delta[Y]/2),endval[Y]])
else :
yedges = yedges+delta[Y]/2
yedges = yedges[:len(yedges)-1]
if isinstance(scale,bool) :
scale = area*(bins[0]*bins[1])/len(vals[0]) if scale else 1.0
Zvals = scale*Zvals/area
f = interpolate.interp2d(yedges,xedges,Zvals,kind=kind)
def density(xa,ya):
size = min( len(xa),len(ya) )
vals = [ f(ya[i],xa[i])[0] for i in range(size) ]
return np.array(vals)
return density
# =================================================================================+
# Only vertices and face indices are exported/imported. |
# FutDev, include : |
# 1. metadata |
# 2. face/vertex color |
# 3. vertex normals for external use for Gouraud shading, etc. |
# 4. line objects |
# 5. code in separate eximport module |
# 6. file types (ascii,binary) : stl, ply, X3D, etc. |
# =================================================================================+
[docs]def get_surfgeom_from_obj(file, **kargs) :
""" Create surface from obj formatted 'file'.
kargs passed to instantiate surface in constructor.
'Proof of concept' for future versions.
"""
# .....................................................................
def file_records(in_file):
try :
f = open(in_file)
buf = f.read()
f.close()
except :
raise ValueError('ERROR: OBJ file not accessed !!!')
for b in buf.split('\n'):
b = ' '.join(b.split()) # remove extra spaces.
if b.startswith('v '):
yield ['v', [float(x) for x in b.split(" ")[1:]]]
elif b.startswith('f '):
datalist = [x for x in b.split(" ")[1:] ]
# only get face indices at the first index...
yield [ 'f', [ int(y.split("/")[0])-1 for y in datalist] ]
else:
yield ['', ""]
def get_v_f(in_file):
vertices,faces = [],[]
for k, v in file_records(in_file):
if k == 'v': vertices.append(v)
elif k == 'f': faces.append(v)
if not len(faces) or not len(vertices): return None, None
return vertices, faces
# .....................................................................
v,f = get_v_f(file)
v = np.array(v)
v = v[:,[0,2,1]]
v = np.multiply(v,[1,-1,1])
if v is not None: surface = Surface3DCollection(v,f,**kargs)
else: surface = SphericalSurface(0,'tetra',color='r')
return surface
[docs]def save_surfgeom_to_obj(file,surface) :
""" Export 'surface' geometry to obj formatted file.
'Proof of concept' for future versions.
"""
v = surface.vertexCoor
v = v[:,[0,2,1]]
v = np.multiply(v,[1,1,-1])
f = surface.fvIndices+1 # obj file indices start at 1
vs = [ 'v '+np.array2string(x,prefix='',suffix='')[1:-1] for x in v ]
fs = [ 'f '+np.array2string(x,prefix='',suffix='')[1:-1] for x in f ]
geomLines = vs + fs
with open(file,"w" ) as fileout :
for line in geomLines :
fileout.write(line+'\n')
fileout.close()
return
