Source code for s3dlib.surface

# 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

#import matplotlib as mpl
from matplotlib import cm, colors, image
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


# FutDev: Future Development notes.
# MROI:   Minimum Return On coding effort Investment.
#         (the juice isn't worth the squeeze for current release)

# +------------------------------------------------------------------------------+
# |  The current verion of s3dlib was developed using Matplotlib version 3.0.2   |
# |  As a result, the following Matplotlib 3DCollection 'private' properties     |
# |  were needed for development:                                                |
# |     _facecolors, _edgecolors, _facecolors3d, _edgecolors3D                   |
# |  Note: further versions of Matplotlib include methods:                       |
# |     get_facecolor() ,       get_edgecolor()                                  |
# |  Private properties to be removed in future S3Dlib versions, MROI.           |
# +------------------------------------------------------------------------------+


[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 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=True) : """Coordinate transformation. To be overriddden by any subclass not in Cartesian coordinates. """ xyz = np.array(xyz) return xyz
@staticmethod def _map3Dto2Dsquare(xyz) : """ Map surface to rectangular unit coordinates (0,1). """ x,y,z = 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._facecolors 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]) 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, **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. 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) 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) 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._facecolors3d[0]).flatten().tolist() self.vertexColor = np.array(self._facecolors3d[0]).copy()[np.newaxis,:] # following assigned when surface color applied using an operation self.vertexValues = None self.faceValues = None if len (self._edgecolors3d) is 0 : self.baseEdgeColor = self.baseSurfaceColor else : self.baseEdgeColor = np.array(self._edgecolors3d[0]).flatten().tolist() 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._bounds = {} self._set_geometric_bounds() self._bounds['vlim'] = _DFT_VLIM self._bounds['vertvlim'] = _DFT_VERTVLIM
[docs] def __add__(self, othr) : """Surface3DCollection object from vertices, faces and colors of two surface objects. """ 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._facecolors 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._facecolors 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._facecolors 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, fColors = 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 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._edgecolors = get_edge_color() 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._facecolors @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] is 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._facecolors : 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._facecolors : 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._facecolors) 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._facecolors 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] is 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._facecolors 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. 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 """ if operation is None : operation = lambda c : c[2] 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) if cname is not None : self.valuesName = cname # ======================================================================== 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._facecolors 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._facecolors 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._facecolors 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._facecolors 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 hidding face edges. 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._facecolors 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._facecolors 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._facecolors 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] faded_colors = np.concatenate((fc_less_alpha, fadex[:,np.newaxis] ), axis=1) self.set_color(faded_colors) return self
[docs] def transform(self, rotate=None, scale=1.0, translate=[0,0,0] ) : """ Linear transformation of the surface object. Parameters ---------- rotate : array-like, optional, default: None A 3 by 3 rotation matrix. scale : 3D array or scalar, optional, default: 1.0 Multipliers of the object coordinates, in xyz coordinates. If a single scalar, three multipliers are the same value. translate : array-like, optional, default: [0,0,0] A xyz translation vector of object origin. Returns ------- self : surface object """ if rotate is None : rotate = np.identity(3) if isinstance(scale,(int,float)) : scale = np.full(3,scale).astype(float) #..................................... def isOrth(M) : M = np.array(M) tol = 0.00001 prod = np.dot(M,M.T) np.fill_diagonal(prod,0.0) return not (np.abs(prod)>tol).any() #..................................... trans = lambda V,S,T : np.add( np.dot( np.multiply(V,S),rotate) ,T) #..................................... if not isOrth(rotate) : warnings.warn('Transform rotate matrix is NOT Orthoginal.') v = trans( self.vertexCoor, scale, translate ) if self.vector_field is not None : nScal = np.divide(scale, np.linalg.norm(scale) ) self.vector_field = trans( self.vector_field, nScal, [0,0,0] ) self.vertexCoor = v self._set_geometric_bounds() self.set_verts( v[ np.array(self.fvIndices) ] ) self._surfaceShapeChanged = True self.scale = scale self.translation = translate self.rotation = rotate return self
[docs] def get_transformAxis(self, **kargs) : """ Line3DCollection of the 'last' transformed surface coordinate axis. Parameters ---------- lenmult : scalar or 3D array, optional, default: 1 Scalar multiplier of the three coordinate axis. width : scalar, optional, default: 1.5 Line width of the coordinate axis. color : a Color or 3D array of Colors, optional, default: ['r','g','b'] negaxis : boolean, optional, default: False If True, include the negative axis, otherwise axis start at origin. Returns ------- Line3DCollection object """ kargDft = { 'lenmult':1.0, 'width':1.5, 'color':['r','g','b'], 'negaxis':False } KW = KWprocessor(kargDft,'get_transformAxis',**kargs) lenmult = KW.getVal('lenmult') width = KW.getVal('width') color = KW.getVal('color') negaxis = KW.getVal('negaxis') selfProp = { 'scale' : self.scale, 'rotation' : self.rotation, 'translation' : self.translation } return _transformAxis(selfProp,lenmult, width, color, negaxis)
[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) 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 """ extDft = { 'direction':[0,0,1.0], 'name': None, 'color': None } kargDft = { **_COOR_KWARGS, **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) 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 """ extDft = { 'direction':[0,0,1.0], 'name': None, 'color': None } kargDft = { **_COOR_KWARGS, **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]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) 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
@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) : """PlanarSurface object based on a meshgrid.""" 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() 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 is '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) : """ Linear scale and translate the surface. Used for scaling a surface 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 ------- 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) 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,z = 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) : """ Override coordinate transformation. 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, including the z coordinate( which remains unchanged). 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
@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] is '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 is 'squ_c' : cylbase = 'squ_s' if basetype is '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 is 'squ_c' or basetype is '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) : """(reserved, not implemented)""" surface = CylindricalSurface(rez,name=name,**kargs) warnings.warn('Default CylindricalSurface object, not rand, returned') 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) : """ Override coordinate transformation. 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
@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] is '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) 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.""" #..................................... 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 is 'cube_a' : surface = SphericalSurface(basetype='cube',name=name, **kwargs) else : surface = SphericalSurface(basetype=basetype,name=name, **kwargs) surface._rez = rez if basetype is '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 is '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 is 'dodeca' : v,f,e = SphericalSurface.get_dodecahedron() if rez>0 : fc = surface._facecolors 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 None 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 is 'cube' : return get_cube() if getArrays is '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) : """ Override coordinate transformation. 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
@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] is '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 is 'octa_c' : cylbase = 'squ_s' if basetype is '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 is 'octa_c' or basetype is '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] 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,X,Y = 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] is 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.""" # 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) : """ColorLine3DCollection object from vertices, segments and colors of two line objects. """ # ............................................................................ 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._edgecolors,topSegmIndices) botEdgeColors = get_lineColors(othr._edgecolors,botSegmIndices) 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 _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._edgecolors # colors before line shredding 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._edgecolors # colors before line shredding numColors = len(initColors) linesegs = self.lineIndices colors