<s> from collections import OrderedDict <EOL> import numpy as np <EOL> from scipy . spatial import distance , Delaunay <EOL> from scipy import sparse <EOL> import scipy . sparse . linalg <EOL> import functools <EOL> import numexpr as ne <EOL> def _memo ( fn ) : <EOL> """<STR_LIT>""" <EOL> @ functools . wraps ( fn ) <EOL> def memofn ( self , * args , ** kwargs ) : <EOL> if id ( fn ) not in self . _cache : <EOL> self . _cache [ id ( fn ) ] = fn ( self ) <EOL> return self . _cache [ id ( fn ) ] <EOL> return memofn <EOL> class Surface ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , pts , polys ) : <EOL> """<STR_LIT>""" <EOL> self . pts = pts . astype ( np . double ) <EOL> self . polys = polys <EOL> self . _cache = dict ( ) <EOL> self . _rlfac_solvers = dict ( ) <EOL> self . _nLC_solvers = dict ( ) <EOL> @ property <EOL> @ _memo <EOL> def ppts ( self ) : <EOL> """<STR_LIT>""" <EOL> return self . pts [ self . polys ] <EOL> @ property <EOL> @ _memo <EOL> def connected ( self ) : <EOL> """<STR_LIT>""" <EOL> npt = len ( self . pts ) <EOL> npoly = len ( self . polys ) <EOL> return sparse . coo_matrix ( ( np . ones ( ( <NUM_LIT:3> * npoly , ) ) , <EOL> ( np . hstack ( self . polys . T ) , <EOL> np . tile ( range ( npoly ) , ( <NUM_LIT:1> , <NUM_LIT:3> ) ) . squeeze ( ) ) ) , <EOL> ( npt , npoly ) ) . tocsr ( ) <EOL> @ property <EOL> @ _memo <EOL> def adj ( self ) : <EOL> """<STR_LIT>""" <EOL> npt = len ( self . pts ) <EOL> npoly = len ( self . polys ) <EOL> adj1 = sparse . coo_matrix ( ( np . ones ( ( npoly , ) ) , <EOL> ( self . polys [ : , <NUM_LIT:0> ] , self . polys [ : , <NUM_LIT:1> ] ) ) , ( npt , npt ) ) <EOL> adj2 = sparse . coo_matrix ( ( np . ones ( ( npoly , ) ) , <EOL> ( self . polys [ : , <NUM_LIT:0> ] , self . polys [ : , <NUM_LIT:2> ] ) ) , ( npt , npt ) ) <EOL> adj3 = sparse . coo_matrix ( ( np . ones ( ( npoly , ) ) , <EOL> ( self . polys [ : , <NUM_LIT:1> ] , self . polys [ : , <NUM_LIT:2> ] ) ) , ( npt , npt ) ) <EOL> alladj = ( adj1 + adj2 + adj3 ) . tocsr ( ) <EOL> return alladj + alladj . T <EOL> @ property <EOL> @ _memo <EOL> def face_normals ( self ) : <EOL> """<STR_LIT>""" <EOL> nnfnorms = np . cross ( self . ppts [ : , <NUM_LIT:1> ] - self . ppts [ : , <NUM_LIT:0> ] , <EOL> self . ppts [ : , <NUM_LIT:2> ] - self . ppts [ : , <NUM_LIT:0> ] ) <EOL> nfnorms = nnfnorms / np . sqrt ( ( nnfnorms ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) [ : , np . newaxis ] <EOL> return np . nan_to_num ( nfnorms ) <EOL> @ property <EOL> @ _memo <EOL> def vertex_normals ( self ) : <EOL> """<STR_LIT>""" <EOL> nnvnorms = np . nan_to_num ( self . connected . dot ( self . face_normals ) / self . connected . sum ( <NUM_LIT:1> ) ) . A <EOL> return nnvnorms / np . sqrt ( ( nnvnorms ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) [ : , np . newaxis ] <EOL> @ property <EOL> @ _memo <EOL> def face_areas ( self ) : <EOL> """<STR_LIT>""" <EOL> nnfnorms = np . cross ( self . ppts [ : , <NUM_LIT:1> ] - self . ppts [ : , <NUM_LIT:0> ] , <EOL> self . ppts [ : , <NUM_LIT:2> ] - self . ppts [ : , <NUM_LIT:0> ] ) <EOL> return np . sqrt ( ( nnfnorms ** <NUM_LIT:2> ) . sum ( - <NUM_LIT:1> ) ) / <NUM_LIT:2> <EOL> @ property <EOL> @ _memo <EOL> def cotangent_weights ( self ) : <EOL> """<STR_LIT>""" <EOL> ppts = self . ppts <EOL> cots1 = ( ( ppts [ : , <NUM_LIT:1> ] - ppts [ : , <NUM_LIT:0> ] ) * <EOL> ( ppts [ : , <NUM_LIT:2> ] - ppts [ : , <NUM_LIT:0> ] ) ) . sum ( <NUM_LIT:1> ) / np . sqrt ( ( np . cross ( ppts [ : , <NUM_LIT:1> ] - ppts [ : , <NUM_LIT:0> ] , <EOL> ppts [ : , <NUM_LIT:2> ] - ppts [ : , <NUM_LIT:0> ] ) ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) <EOL> cots2 = ( ( ppts [ : , <NUM_LIT:2> ] - ppts [ : , <NUM_LIT:1> ] ) * <EOL> ( ppts [ : , <NUM_LIT:0> ] - ppts [ : , <NUM_LIT:1> ] ) ) . sum ( <NUM_LIT:1> ) / np . sqrt ( ( np . cross ( ppts [ : , <NUM_LIT:2> ] - ppts [ : , <NUM_LIT:1> ] , <EOL> ppts [ : , <NUM_LIT:0> ] - ppts [ : , <NUM_LIT:1> ] ) ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) <EOL> cots3 = ( ( ppts [ : , <NUM_LIT:0> ] - ppts [ : , <NUM_LIT:2> ] ) * <EOL> ( ppts [ : , <NUM_LIT:1> ] - ppts [ : , <NUM_LIT:2> ] ) ) . sum ( <NUM_LIT:1> ) / np . sqrt ( ( np . cross ( ppts [ : , <NUM_LIT:0> ] - ppts [ : , <NUM_LIT:2> ] , <EOL> ppts [ : , <NUM_LIT:1> ] - ppts [ : , <NUM_LIT:2> ] ) ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) <EOL> cots = np . vstack ( [ cots1 , cots2 , cots3 ] ) <EOL> cots [ np . isinf ( cots ) ] = <NUM_LIT:0> <EOL> cots [ np . isnan ( cots ) ] = <NUM_LIT:0> <EOL> return cots <EOL> @ property <EOL> @ _memo <EOL> def laplace_operator ( self ) : <EOL> """<STR_LIT>""" <EOL> D = self . connected . dot ( self . face_areas ) / <NUM_LIT> <EOL> npt = len ( self . pts ) <EOL> cots1 , cots2 , cots3 = self . cotangent_weights <EOL> W1 = sparse . coo_matrix ( ( cots1 , ( self . polys [ : , <NUM_LIT:1> ] , self . polys [ : , <NUM_LIT:2> ] ) ) , ( npt , npt ) ) <EOL> W2 = sparse . coo_matrix ( ( cots2 , ( self . polys [ : , <NUM_LIT:2> ] , self . polys [ : , <NUM_LIT:0> ] ) ) , ( npt , npt ) ) <EOL> W3 = sparse . coo_matrix ( ( cots3 , ( self . polys [ : , <NUM_LIT:0> ] , self . polys [ : , <NUM_LIT:1> ] ) ) , ( npt , npt ) ) <EOL> W = ( W1 + W1 . T + W2 + W2 . T + W3 + W3 . T ) . tocsr ( ) / <NUM_LIT> <EOL> V = sparse . dia_matrix ( ( np . array ( W . sum ( <NUM_LIT:0> ) ) . ravel ( ) , [ <NUM_LIT:0> ] ) , ( npt , npt ) ) <EOL> Be1 = sparse . coo_matrix ( ( self . face_areas , ( self . polys [ : , <NUM_LIT:1> ] , self . polys [ : , <NUM_LIT:2> ] ) ) , ( npt , npt ) ) <EOL> Be2 = sparse . coo_matrix ( ( self . face_areas , ( self . polys [ : , <NUM_LIT:2> ] , self . polys [ : , <NUM_LIT:0> ] ) ) , ( npt , npt ) ) <EOL> Be3 = sparse . coo_matrix ( ( self . face_areas , ( self . polys [ : , <NUM_LIT:0> ] , self . polys [ : , <NUM_LIT:1> ] ) ) , ( npt , npt ) ) <EOL> Bd = self . connected . dot ( self . face_areas ) / <NUM_LIT:6> <EOL> dBd = scipy . sparse . dia_matrix ( ( Bd , [ <NUM_LIT:0> ] ) , ( len ( D ) , len ( D ) ) ) <EOL> B = ( Be1 + Be1 . T + Be2 + Be2 . T + Be3 + Be3 . T ) / <NUM_LIT:12> + dBd <EOL> return B , D , W , V <EOL> def mean_curvature ( self ) : <EOL> """<STR_LIT>""" <EOL> B , D , W , V = self . laplace_operator <EOL> npt = len ( D ) <EOL> Dinv = sparse . dia_matrix ( ( D ** - <NUM_LIT:1> , [ <NUM_LIT:0> ] ) , ( npt , npt ) ) . tocsr ( ) <EOL> L = Dinv . dot ( ( V - W ) ) <EOL> curv = ( L . dot ( self . pts ) * self . vertex_normals ) . sum ( <NUM_LIT:1> ) <EOL> return curv <EOL> def smooth ( self , scalars , factor = <NUM_LIT:1.0> , iterations = <NUM_LIT:1> ) : <EOL> """<STR_LIT>""" <EOL> if factor == <NUM_LIT:0.0> : <EOL> return scalars <EOL> B , D , W , V = self . laplace_operator <EOL> npt = len ( D ) <EOL> lfac = sparse . dia_matrix ( ( D , [ <NUM_LIT:0> ] ) , ( npt , npt ) ) - factor * ( W - V ) <EOL> goodrows = np . nonzero ( ~ np . array ( lfac . sum ( <NUM_LIT:0> ) == <NUM_LIT:0> ) . ravel ( ) ) [ <NUM_LIT:0> ] <EOL> lfac_solver = sparse . linalg . dsolve . factorized ( lfac [ goodrows ] [ : , goodrows ] ) <EOL> to_smooth = scalars . copy ( ) <EOL> for _ in range ( iterations ) : <EOL> from_smooth = lfac_solver ( ( D * to_smooth ) [ goodrows ] ) <EOL> to_smooth [ goodrows ] = from_smooth <EOL> smscalars = np . zeros ( scalars . shape ) <EOL> smscalars [ goodrows ] = from_smooth <EOL> return smscalars <EOL> @ property <EOL> @ _memo <EOL> def avg_edge_length ( self ) : <EOL> """<STR_LIT>""" <EOL> adj = self . adj <EOL> tadj = sparse . triu ( adj , <NUM_LIT:1> ) <EOL> edgelens = np . sqrt ( ( ( self . pts [ tadj . row ] - self . pts [ tadj . col ] ) ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) <EOL> return edgelens . mean ( ) <EOL> def surface_gradient ( self , scalars , at_verts = True ) : <EOL> """<STR_LIT>""" <EOL> pu = scalars [ self . polys ] <EOL> fe12 , fe23 , fe31 = [ f . T for f in self . _facenorm_cross_edge ] <EOL> pu1 , pu2 , pu3 = pu . T <EOL> fa = self . face_areas <EOL> gradu = np . nan_to_num ( ne . evaluate ( "<STR_LIT>" ) . T ) <EOL> if at_verts : <EOL> return ( self . connected . dot ( gradu ) . T / self . connected . sum ( <NUM_LIT:1> ) . A . squeeze ( ) ) . T <EOL> return gradu <EOL> def _create_biharmonic_solver ( self , boundary_verts , clip_D = <NUM_LIT:0.1> ) : <EOL> """<STR_LIT>""" <EOL> try : <EOL> from sksparse . sparse . cholmod import cholesky <EOL> except ImportError : <EOL> from scikits . sparse . cholmod import cholesky <EOL> B , D , W , V = self . laplace_operator <EOL> npt = len ( D ) <EOL> g = np . nonzero ( D > <NUM_LIT:0> ) [ <NUM_LIT:0> ] <EOL> notboundary = np . setdiff1d ( np . arange ( npt ) [ g ] , boundary_verts ) <EOL> D = np . clip ( D , clip_D , D . max ( ) ) <EOL> Dinv = sparse . dia_matrix ( ( D ** - <NUM_LIT:1> , [ <NUM_LIT:0> ] ) , ( npt , npt ) ) . tocsr ( ) <EOL> L = Dinv . dot ( ( V - W ) ) <EOL> lhs = ( V - W ) . dot ( L ) <EOL> lhsfac = cholesky ( lhs [ notboundary ] [ : , notboundary ] ) <EOL> return lhs , D , Dinv , lhsfac , notboundary <EOL> def _create_interp ( self , verts , bhsolver = None ) : <EOL> """<STR_LIT>""" <EOL> if bhsolver is None : <EOL> lhs , D , Dinv , lhsfac , notb = self . _create_biharmonic_solver ( verts ) <EOL> else : <EOL> lhs , D , Dinv , lhsfac , notb = bhsolver <EOL> npt = len ( D ) <EOL> def _interp ( vals ) : <EOL> """<STR_LIT>""" <EOL> v2 = np . atleast_2d ( vals ) <EOL> nd , nv = v2 . shape <EOL> ij = np . zeros ( ( <NUM_LIT:2> , nv * nd ) ) <EOL> ij [ <NUM_LIT:0> ] = np . array ( verts ) [ np . repeat ( np . arange ( nv ) , nd ) ] <EOL> ij [ <NUM_LIT:1> ] = np . tile ( np . arange ( nd ) , nv ) <EOL> r = sparse . csr_matrix ( ( vals . T . ravel ( ) , ij ) , shape = ( npt , nd ) ) <EOL> vr = lhs . dot ( r ) <EOL> phi = lhsfac . solve_A ( - vr [ notb ] . todense ( ) ) <EOL> tphi = np . zeros ( ( npt , nd ) ) <EOL> tphi [ notb ] = phi <EOL> tphi [ verts ] = v2 . T <EOL> return tphi <EOL> return _interp <EOL> def interp ( self , verts , vals ) : <EOL> """<STR_LIT>""" <EOL> return self . _create_interp ( verts ) ( vals ) <EOL> @ property <EOL> @ _memo <EOL> def _facenorm_cross_edge ( self ) : <EOL> ppts = self . ppts <EOL> fnorms = self . face_normals <EOL> fe12 = np . cross ( fnorms , ppts [ : , <NUM_LIT:1> ] - ppts [ : , <NUM_LIT:0> ] ) <EOL> fe23 = np . cross ( fnorms , ppts [ : , <NUM_LIT:2> ] - ppts [ : , <NUM_LIT:1> ] ) <EOL> fe31 = np . cross ( fnorms , ppts [ : , <NUM_LIT:0> ] - ppts [ : , <NUM_LIT:2> ] ) <EOL> return fe12 , fe23 , fe31 <EOL> def approx_geodesic_distance ( self , verts , m = <NUM_LIT:0.1> ) : <EOL> npt = len ( self . pts ) <EOL> t = m * self . avg_edge_length ** <NUM_LIT:2> <EOL> if m not in self . _rlfac_solvers : <EOL> B , D , W , V = self . laplace_operator <EOL> nLC = W - V <EOL> spD = sparse . dia_matrix ( ( D , [ <NUM_LIT:0> ] ) , ( npt , npt ) ) . tocsr ( ) <EOL> lfac = spD - t * nLC <EOL> goodrows = np . nonzero ( ~ np . array ( lfac . sum ( <NUM_LIT:0> ) == <NUM_LIT:0> ) . ravel ( ) ) [ <NUM_LIT:0> ] <EOL> self . _goodrows = goodrows <EOL> self . _rlfac_solvers [ m ] = sparse . linalg . dsolve . factorized ( lfac [ goodrows ] [ : , goodrows ] ) <EOL> u0 = np . zeros ( ( npt , ) ) <EOL> u0 [ verts ] = <NUM_LIT:1.0> <EOL> goodu = self . _rlfac_solvers [ m ] ( u0 [ self . _goodrows ] ) <EOL> u = np . zeros ( ( npt , ) ) <EOL> u [ self . _goodrows ] = goodu <EOL> return - <NUM_LIT:4> * t * np . log ( u ) <EOL> def geodesic_distance ( self , verts , m = <NUM_LIT:1.0> , fem = False ) : <EOL> """<STR_LIT>""" <EOL> npt = len ( self . pts ) <EOL> if m not in self . _rlfac_solvers or m not in self . _nLC_solvers : <EOL> B , D , W , V = self . laplace_operator <EOL> nLC = W - V <EOL> if not fem : <EOL> spD = sparse . dia_matrix ( ( D , [ <NUM_LIT:0> ] ) , ( npt , npt ) ) . tocsr ( ) <EOL> else : <EOL> spD = B <EOL> t = m * self . avg_edge_length ** <NUM_LIT:2> <EOL> lfac = spD - t * nLC <EOL> goodrows = np . nonzero ( ~ np . array ( lfac . sum ( <NUM_LIT:0> ) == <NUM_LIT:0> ) . ravel ( ) ) [ <NUM_LIT:0> ] <EOL> self . _goodrows = goodrows <EOL> self . _rlfac_solvers [ m ] = sparse . linalg . dsolve . factorized ( lfac [ goodrows ] [ : , goodrows ] ) <EOL> self . _nLC_solvers [ m ] = sparse . linalg . dsolve . factorized ( nLC [ goodrows ] [ : , goodrows ] ) <EOL> u0 = np . zeros ( ( npt , ) ) <EOL> u0 [ verts ] = <NUM_LIT:1.0> <EOL> goodu = self . _rlfac_solvers [ m ] ( u0 [ self . _goodrows ] ) <EOL> u = np . zeros ( ( npt , ) ) <EOL> u [ self . _goodrows ] = goodu <EOL> gradu = self . surface_gradient ( u , at_verts = False ) <EOL> graduT = gradu . T <EOL> gusum = ne . evaluate ( "<STR_LIT>" ) <EOL> X = np . nan_to_num ( ne . evaluate ( "<STR_LIT>" ) . T ) <EOL> c32 , c13 , c21 = self . _cot_edge <EOL> x1 = <NUM_LIT:0.5> * ( c32 * X ) . sum ( <NUM_LIT:1> ) <EOL> x2 = <NUM_LIT:0.5> * ( c13 * X ) . sum ( <NUM_LIT:1> ) <EOL> x3 = <NUM_LIT:0.5> * ( c21 * X ) . sum ( <NUM_LIT:1> ) <EOL> conn1 , conn2 , conn3 = self . _polyconn <EOL> divx = conn1 . dot ( x1 ) + conn2 . dot ( x2 ) + conn3 . dot ( x3 ) <EOL> goodphi = self . _nLC_solvers [ m ] ( divx [ self . _goodrows ] ) <EOL> phi = np . zeros ( ( npt , ) ) <EOL> phi [ self . _goodrows ] = goodphi - goodphi . min ( ) <EOL> phi [ verts ] = <NUM_LIT:0.0> <EOL> return phi <EOL> @ property <EOL> @ _memo <EOL> def _cot_edge ( self ) : <EOL> ppts = self . ppts <EOL> cots1 , cots2 , cots3 = self . cotangent_weights <EOL> c3 = cots3 [ : , np . newaxis ] * ( ppts [ : , <NUM_LIT:1> ] - ppts [ : , <NUM_LIT:0> ] ) <EOL> c2 = cots2 [ : , np . newaxis ] * ( ppts [ : , <NUM_LIT:0> ] - ppts [ : , <NUM_LIT:2> ] ) <EOL> c1 = cots1 [ : , np . newaxis ] * ( ppts [ : , <NUM_LIT:2> ] - ppts [ : , <NUM_LIT:1> ] ) <EOL> c32 = c3 - c2 <EOL> c13 = c1 - c3 <EOL> c21 = c2 - c1 <EOL> return c32 , c13 , c21 <EOL> @ property <EOL> @ _memo <EOL> def _polyconn ( self ) : <EOL> npt = len ( self . pts ) <EOL> npoly = len ( self . polys ) <EOL> o = np . ones ( ( npoly , ) ) <EOL> c1 = sparse . coo_matrix ( ( o , ( self . polys [ : , <NUM_LIT:0> ] , range ( npoly ) ) ) , ( npt , npoly ) ) . tocsr ( ) <EOL> c2 = sparse . coo_matrix ( ( o , ( self . polys [ : , <NUM_LIT:1> ] , range ( npoly ) ) ) , ( npt , npoly ) ) . tocsr ( ) <EOL> c3 = sparse . coo_matrix ( ( o , ( self . polys [ : , <NUM_LIT:2> ] , range ( npoly ) ) ) , ( npt , npoly ) ) . tocsr ( ) <EOL> return c1 , c2 , c3 <EOL> @ property <EOL> @ _memo <EOL> def graph ( self ) : <EOL> """<STR_LIT>""" <EOL> import networkx as nx <EOL> def iter_surfedges ( tris ) : <EOL> for a , b , c in tris : <EOL> yield a , b <EOL> yield b , c <EOL> yield a , c <EOL> def make_surface_graph ( tris ) : <EOL> graph = nx . Graph ( ) <EOL> graph . add_edges_from ( iter_surfedges ( tris ) ) <EOL> return graph <EOL> return make_surface_graph ( self . polys ) <EOL> def get_graph ( self ) : <EOL> return self . graph <EOL> def extract_chunk ( self , nfaces = <NUM_LIT:100> , seed = None , auxpts = None ) : <EOL> '''<STR_LIT>''' <EOL> node = seed <EOL> if seed is None : <EOL> node = np . random . randint ( len ( self . pts ) ) <EOL> ptmap = dict ( ) <EOL> queue = [ node ] <EOL> faces = set ( ) <EOL> visited = set ( [ node ] ) <EOL> while len ( faces ) < nfaces and len ( queue ) > <NUM_LIT:0> : <EOL> node = queue . pop ( <NUM_LIT:0> ) <EOL> for face in self . connected [ node ] . indices : <EOL> if face not in faces : <EOL> faces . add ( face ) <EOL> for pt in self . polys [ face ] : <EOL> if pt not in visited : <EOL> visited . add ( pt ) <EOL> queue . append ( pt ) <EOL> pts , aux , polys = [ ] , [ ] , [ ] <EOL> for face in faces : <EOL> for pt in self . polys [ face ] : <EOL> if pt not in ptmap : <EOL> ptmap [ pt ] = len ( pts ) <EOL> pts . append ( self . pts [ pt ] ) <EOL> if auxpts is not None : <EOL> aux . append ( auxpts [ pt ] ) <EOL> polys . append ( [ ptmap [ p ] for p in self . polys [ face ] ] ) <EOL> if auxpts is not None : <EOL> return np . array ( pts ) , np . array ( aux ) , np . array ( polys ) <EOL> return np . array ( pts ) , np . array ( polys ) <EOL> def polyhedra ( self , wm ) : <EOL> '''<STR_LIT>''' <EOL> for p , facerow in enumerate ( self . connected ) : <EOL> faces = facerow . indices <EOL> pts , polys = _ptset ( ) , _quadset ( ) <EOL> if len ( faces ) > <NUM_LIT:0> : <EOL> poly = np . roll ( self . polys [ faces [ <NUM_LIT:0> ] ] , - np . nonzero ( self . polys [ faces [ <NUM_LIT:0> ] ] == p ) [ <NUM_LIT:0> ] [ <NUM_LIT:0> ] ) <EOL> assert pts [ wm [ p ] ] == <NUM_LIT:0> <EOL> assert pts [ self . pts [ p ] ] == <NUM_LIT:1> <EOL> pts [ wm [ poly [ [ <NUM_LIT:0> , <NUM_LIT:1> ] ] ] . mean ( <NUM_LIT:0> ) ] <EOL> pts [ self . pts [ poly [ [ <NUM_LIT:0> , <NUM_LIT:1> ] ] ] . mean ( <NUM_LIT:0> ) ] <EOL> for face in faces : <EOL> poly = np . roll ( self . polys [ face ] , - np . nonzero ( self . polys [ face ] == p ) [ <NUM_LIT:0> ] [ <NUM_LIT:0> ] ) <EOL> a = pts [ wm [ poly ] . mean ( <NUM_LIT:0> ) ] <EOL> b = pts [ self . pts [ poly ] . mean ( <NUM_LIT:0> ) ] <EOL> c = pts [ wm [ poly [ [ <NUM_LIT:0> , <NUM_LIT:2> ] ] ] . mean ( <NUM_LIT:0> ) ] <EOL> d = pts [ self . pts [ poly [ [ <NUM_LIT:0> , <NUM_LIT:2> ] ] ] . mean ( <NUM_LIT:0> ) ] <EOL> e = pts [ wm [ poly [ [ <NUM_LIT:0> , <NUM_LIT:1> ] ] ] . mean ( <NUM_LIT:0> ) ] <EOL> f = pts [ self . pts [ poly [ [ <NUM_LIT:0> , <NUM_LIT:1> ] ] ] . mean ( <NUM_LIT:0> ) ] <EOL> polys ( ( <NUM_LIT:0> , c , a , e ) ) <EOL> polys ( ( <NUM_LIT:1> , f , b , d ) ) <EOL> polys ( ( <NUM_LIT:1> , d , c , <NUM_LIT:0> ) ) <EOL> polys ( ( <NUM_LIT:1> , <NUM_LIT:0> , e , f ) ) <EOL> polys ( ( f , e , a , b ) ) <EOL> polys ( ( d , b , a , c ) ) <EOL> yield pts . points , np . array ( list ( polys . triangles ) ) <EOL> def patches ( self , auxpts = None , n = <NUM_LIT:1> ) : <EOL> def align_polys ( p , polys ) : <EOL> x , y = np . nonzero ( polys == p ) <EOL> y = np . vstack ( [ y , ( y + <NUM_LIT:1> ) % <NUM_LIT:3> , ( y + <NUM_LIT:2> ) % <NUM_LIT:3> ] ) . T <EOL> return polys [ np . tile ( x , [ <NUM_LIT:3> , <NUM_LIT:1> ] ) . T , y ] <EOL> def half_edge_align ( p , pts , polys ) : <EOL> poly = align_polys ( p , polys ) <EOL> mid = pts [ poly ] . mean ( <NUM_LIT:1> ) <EOL> left = pts [ poly [ : , [ <NUM_LIT:0> , <NUM_LIT:2> ] ] ] . mean ( <NUM_LIT:1> ) <EOL> right = pts [ poly [ : , [ <NUM_LIT:0> , <NUM_LIT:1> ] ] ] . mean ( <NUM_LIT:1> ) <EOL> s1 = np . array ( np . broadcast_arrays ( pts [ p ] , mid , left ) ) . swapaxes ( <NUM_LIT:0> , <NUM_LIT:1> ) <EOL> s2 = np . array ( np . broadcast_arrays ( pts [ p ] , mid , right ) ) . swapaxes ( <NUM_LIT:0> , <NUM_LIT:1> ) <EOL> return np . vstack ( [ s1 , s2 ] ) <EOL> def half_edge ( p , pts , polys ) : <EOL> poly = align_polys ( p , polys ) <EOL> mid = pts [ poly ] . mean ( <NUM_LIT:1> ) <EOL> left = pts [ poly [ : , [ <NUM_LIT:0> , <NUM_LIT:2> ] ] ] . mean ( <NUM_LIT:1> ) <EOL> right = pts [ poly [ : , [ <NUM_LIT:0> , <NUM_LIT:1> ] ] ] . mean ( <NUM_LIT:1> ) <EOL> stack = np . vstack ( [ mid , left , right , pts [ p ] ] ) <EOL> return stack [ ( distance . cdist ( stack , stack ) == <NUM_LIT:0> ) . sum ( <NUM_LIT:0> ) == <NUM_LIT:1> ] <EOL> for p , facerow in enumerate ( self . connected ) : <EOL> faces = facerow . indices <EOL> if len ( faces ) > <NUM_LIT:0> : <EOL> if n == <NUM_LIT:1> : <EOL> if auxpts is not None : <EOL> pidx = np . unique ( self . polys [ faces ] ) <EOL> yield np . vstack ( [ self . pts [ pidx ] , auxpts [ pidx ] ] ) <EOL> else : <EOL> yield self . pts [ self . polys [ faces ] ] <EOL> elif n == <NUM_LIT:0.5> : <EOL> if auxpts is not None : <EOL> pts = half_edge ( p , self . pts , self . polys [ faces ] ) <EOL> aux = half_edge ( p , auxpts , self . polys [ faces ] ) <EOL> yield np . vstack ( [ pts , aux ] ) <EOL> else : <EOL> yield half_edge_align ( p , self . pts , self . polys [ faces ] ) <EOL> else : <EOL> raise ValueError <EOL> else : <EOL> yield None <EOL> def edge_collapse ( self , p1 , p2 , target ) : <EOL> raise NotImplementedError <EOL> face1 = self . connected [ p1 ] <EOL> face2 = self . connected [ p2 ] <EOL> class _ptset ( object ) : <EOL> def __init__ ( self ) : <EOL> self . idx = OrderedDict ( ) <EOL> def __getitem__ ( self , idx ) : <EOL> idx = tuple ( idx ) <EOL> if idx not in self . idx : <EOL> self . idx [ idx ] = len ( self . idx ) <EOL> return self . idx [ idx ] <EOL> @ property <EOL> def points ( self ) : <EOL> return np . array ( list ( self . idx . keys ( ) ) ) <EOL> class _quadset ( object ) : <EOL> def __init__ ( self ) : <EOL> self . polys = dict ( ) <EOL> def __call__ ( self , quad ) : <EOL> idx = tuple ( sorted ( quad ) ) <EOL> if idx in self . polys : <EOL> del self . polys [ idx ] <EOL> else : <EOL> self . polys [ idx ] = quad <EOL> @ property <EOL> def triangles ( self ) : <EOL> for quad in list ( self . polys . values ( ) ) : <EOL> yield quad [ : <NUM_LIT:3> ] <EOL> yield [ quad [ <NUM_LIT:0> ] , quad [ <NUM_LIT:2> ] , quad [ <NUM_LIT:3> ] ] <EOL> class Distortion ( object ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self , flat , ref , polys ) : <EOL> """<STR_LIT>""" <EOL> self . flat = flat <EOL> self . ref = ref <EOL> self . polys = polys <EOL> @ property <EOL> def areal ( self ) : <EOL> """<STR_LIT>""" <EOL> def area ( pts , polys ) : <EOL> ppts = pts [ polys ] <EOL> cross = np . cross ( ppts [ : , <NUM_LIT:1> ] - ppts [ : , <NUM_LIT:0> ] , ppts [ : , <NUM_LIT:2> ] - ppts [ : , <NUM_LIT:0> ] ) <EOL> return np . sqrt ( ( cross ** <NUM_LIT:2> ) . sum ( - <NUM_LIT:1> ) ) <EOL> refarea = area ( self . ref , self . polys ) <EOL> flatarea = area ( self . flat , self . polys ) <EOL> tridists = np . log2 ( flatarea / refarea ) <EOL> vertratios = np . zeros ( ( len ( self . ref ) , ) ) <EOL> vertratios [ self . polys [ : , <NUM_LIT:0> ] ] += tridists <EOL> vertratios [ self . polys [ : , <NUM_LIT:1> ] ] += tridists <EOL> vertratios [ self . polys [ : , <NUM_LIT:2> ] ] += tridists <EOL> vertratios /= np . bincount ( self . polys . ravel ( ) ) <EOL> vertratios = np . nan_to_num ( vertratios ) <EOL> vertratios [ vertratios == <NUM_LIT:0> ] = <NUM_LIT:1> <EOL> return vertratios <EOL> @ property <EOL> def metric ( self ) : <EOL> """<STR_LIT>""" <EOL> import networkx as nx <EOL> def iter_surfedges ( tris ) : <EOL> for a , b , c in tris : <EOL> yield a , b <EOL> yield b , c <EOL> yield a , c <EOL> def make_surface_graph ( tris ) : <EOL> graph = nx . Graph ( ) <EOL> graph . add_edges_from ( iter_surfedges ( tris ) ) <EOL> return graph <EOL> G = make_surface_graph ( self . polys ) <EOL> selverts = np . unique ( self . polys . ravel ( ) ) <EOL> ref_dists = [ np . sqrt ( ( ( self . ref [ G . neighbors ( ii ) ] - self . ref [ ii ] ) ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) <EOL> for ii in selverts ] <EOL> flat_dists = [ np . sqrt ( ( ( self . flat [ G . neighbors ( ii ) ] - self . flat [ ii ] ) ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) <EOL> for ii in selverts ] <EOL> msdists = np . array ( [ ( f - r ) . mean ( ) for r , f in zip ( ref_dists , flat_dists ) ] ) <EOL> alldists = np . zeros ( ( len ( self . ref ) , ) ) <EOL> alldists [ selverts ] = msdists <EOL> return alldists <EOL> def tetra_vol ( pts ) : <EOL> '''<STR_LIT>''' <EOL> tetra = pts [ <NUM_LIT:1> : ] - pts [ <NUM_LIT:0> ] <EOL> return np . abs ( np . dot ( tetra [ <NUM_LIT:0> ] , np . cross ( tetra [ <NUM_LIT:1> ] , tetra [ <NUM_LIT:2> ] ) ) ) / <NUM_LIT:6> <EOL> def brick_vol ( pts ) : <EOL> '''<STR_LIT>''' <EOL> return tetra_vol ( pts [ [ <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:2> , <NUM_LIT:4> ] ] ) + tetra_vol ( pts [ [ <NUM_LIT:0> , <NUM_LIT:2> , <NUM_LIT:3> , <NUM_LIT:4> ] ] ) + tetra_vol ( pts [ [ <NUM_LIT:2> , <NUM_LIT:3> , <NUM_LIT:4> , <NUM_LIT:5> ] ] ) <EOL> def sort_polys ( polys ) : <EOL> amin = polys . argmin ( <NUM_LIT:1> ) <EOL> xind = np . arange ( len ( polys ) ) <EOL> return np . array ( [ polys [ xind , amin ] , polys [ xind , ( amin + <NUM_LIT:1> ) % <NUM_LIT:3> ] , polys [ xind , ( amin + <NUM_LIT:2> ) % <NUM_LIT:3> ] ] ) . T <EOL> def face_area ( pts ) : <EOL> '''<STR_LIT>''' <EOL> return <NUM_LIT:0.5> * np . sqrt ( ( np . cross ( pts [ : , <NUM_LIT:1> ] - pts [ : , <NUM_LIT:0> ] , pts [ : , <NUM_LIT:2> ] - pts [ : , <NUM_LIT:0> ] ) ** <NUM_LIT:2> ) . sum ( <NUM_LIT:1> ) ) <EOL> def face_volume ( pts1 , pts2 , polys ) : <EOL> '''<STR_LIT>''' <EOL> vols = np . zeros ( ( len ( polys ) , ) ) <EOL> for i , face in enumerate ( polys ) : <EOL> vols [ i ] = brick_vol ( np . append ( pts1 [ face ] , pts2 [ face ] , axis = <NUM_LIT:0> ) ) <EOL> if i % <NUM_LIT:1000> == <NUM_LIT:0> : <EOL> print ( i ) <EOL> return vols <EOL> def decimate ( pts , polys ) : <EOL> from tvtk . api import tvtk <EOL> pd = tvtk . PolyData ( points = pts , polys = polys ) <EOL> dec = tvtk . DecimatePro ( input = pd ) <EOL> dec . set ( preserve_topology = True , splitting = False , boundary_vertex_deletion = False , target_reduction = <NUM_LIT:1.0> ) <EOL> dec . update ( ) <EOL> dpts = dec . output . points . to_array ( ) <EOL> dpolys = dec . output . polys . to_array ( ) . reshape ( - <NUM_LIT:1> , <NUM_LIT:4> ) [ : , <NUM_LIT:1> : ] <EOL> return dpts , dpolys <EOL> def inside_convex_poly ( pts ) : <EOL> """<STR_LIT>""" <EOL> tri = Delaunay ( pts ) <EOL> return lambda x : tri . find_simplex ( x ) != - <NUM_LIT:1> <EOL> def make_cube ( center = ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) , size = <NUM_LIT:1> ) : <EOL> pts = np . array ( [ ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) , ( <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:0> ) , ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:0> ) , ( <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:0> ) , <EOL> ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:1> ) , ( <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:1> ) , ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:1> ) , ( <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:1> ) ] , dtype = float ) <EOL> pts -= ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) <EOL> polys = np . array ( [ ( <NUM_LIT:0> , <NUM_LIT:2> , <NUM_LIT:3> ) , ( <NUM_LIT:0> , <NUM_LIT:3> , <NUM_LIT:1> ) , ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:4> ) , ( <NUM_LIT:1> , <NUM_LIT:5> , <NUM_LIT:4> ) , <EOL> ( <NUM_LIT:1> , <NUM_LIT:3> , <NUM_LIT:5> ) , ( <NUM_LIT:3> , <NUM_LIT:7> , <NUM_LIT:5> ) , ( <NUM_LIT:2> , <NUM_LIT:7> , <NUM_LIT:3> ) , ( <NUM_LIT:2> , <NUM_LIT:6> , <NUM_LIT:7> ) , <EOL> ( <NUM_LIT:0> , <NUM_LIT:6> , <NUM_LIT:2> ) , ( <NUM_LIT:0> , <NUM_LIT:4> , <NUM_LIT:6> ) , ( <NUM_LIT:4> , <NUM_LIT:7> , <NUM_LIT:6> ) , ( <NUM_LIT:4> , <NUM_LIT:5> , <NUM_LIT:7> ) ] , dtype = np . uint32 ) <EOL> return pts * size + center , polys <EOL> def boundary_edges ( polys ) : <EOL> '''<STR_LIT>''' <EOL> edges = dict ( ) <EOL> for i , poly in enumerate ( np . sort ( polys ) ) : <EOL> for a , b in [ ( <NUM_LIT:0> , <NUM_LIT:1> ) , ( <NUM_LIT:1> , <NUM_LIT:2> ) , ( <NUM_LIT:0> , <NUM_LIT:2> ) ] : <EOL> key = poly [ a ] , poly [ b ] <EOL> if key not in edges : <EOL> edges [ key ] = [ ] <EOL> edges [ key ] . append ( i ) <EOL> epts = [ ] <EOL> for edge , faces in edges . items ( ) : <EOL> if len ( faces ) == <NUM_LIT:1> : <EOL> epts . append ( edge ) <EOL> return np . array ( epts ) <EOL> def trace_poly ( edges ) : <EOL> '''<STR_LIT>''' <EOL> idx = dict ( ( i , set ( [ ] ) ) for i in np . unique ( edges ) ) <EOL> for i , ( x , y ) in enumerate ( edges ) : <EOL> idx [ x ] . add ( i ) <EOL> idx [ y ] . add ( i ) <EOL> eset = set ( range ( len ( edges ) ) ) <EOL> while len ( eset ) > <NUM_LIT:0> : <EOL> eidx = eset . pop ( ) <EOL> poly = list ( edges [ eidx ] ) <EOL> stack = set ( [ eidx ] ) <EOL> while poly [ - <NUM_LIT:1> ] != poly [ <NUM_LIT:0> ] or len ( poly ) == <NUM_LIT:1> : <EOL> next = list ( idx [ poly [ - <NUM_LIT:1> ] ] - stack ) [ <NUM_LIT:0> ] <EOL> eset . remove ( next ) <EOL> stack . add ( next ) <EOL> if edges [ next ] [ <NUM_LIT:0> ] == poly [ - <NUM_LIT:1> ] : <EOL> poly . append ( edges [ next ] [ <NUM_LIT:1> ] ) <EOL> elif edges [ next ] [ <NUM_LIT:1> ] == poly [ - <NUM_LIT:1> ] : <EOL> poly . append ( edges [ next ] [ <NUM_LIT:0> ] ) <EOL> else : <EOL> raise Exception <EOL> yield poly <EOL> def rasterize ( poly , shape = ( <NUM_LIT> , <NUM_LIT> ) ) : <EOL> import subprocess as sp <EOL> import cStringIO <EOL> import shlex <EOL> from PIL import Image <EOL> polygon = "<STR_LIT:U+0020>" . join ( [ "<STR_LIT>" % tuple ( p [ : : - <NUM_LIT:1> ] ) for p in np . array ( poly ) - ( <NUM_LIT> , <NUM_LIT> ) ] ) <EOL> cmd = '<STR_LIT>' % ( shape [ <NUM_LIT:0> ] , shape [ <NUM_LIT:1> ] , polygon ) <EOL> proc = sp . Popen ( shlex . split ( cmd ) , stdout = sp . PIPE ) <EOL> png = cStringIO . StringIO ( proc . communicate ( ) [ <NUM_LIT:0> ] ) <EOL> im = Image . open ( png ) <EOL> return ( np . array ( im ) [ : , : , <NUM_LIT:0> ] > <NUM_LIT> ) . T <EOL> def voxelize ( pts , polys , shape = ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) , center = ( <NUM_LIT> , <NUM_LIT> , <NUM_LIT> ) , mp = True ) : <EOL> from tvtk . api import tvtk <EOL> pd = tvtk . PolyData ( points = pts + center + ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) , polys = polys ) <EOL> plane = tvtk . Planes ( normals = [ ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:1> ) ] , points = [ ( <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ) ] ) <EOL> clip = tvtk . ClipPolyData ( clip_function = plane , input = pd ) <EOL> feats = tvtk . FeatureEdges ( <EOL> manifold_edges = False , <EOL> non_manifold_edges = False , <EOL> feature_edges = False , <EOL> boundary_edges = True , <EOL> input = clip . output ) <EOL> def func ( i ) : <EOL> plane . points = [ ( <NUM_LIT:0> , <NUM_LIT:0> , i ) ] <EOL> feats . update ( ) <EOL> vox = np . zeros ( shape [ : <NUM_LIT:2> ] [ : : - <NUM_LIT:1> ] , np . uint8 ) <EOL> if feats . output . number_of_lines > <NUM_LIT:0> : <EOL> epts = feats . output . points . to_array ( ) <EOL> edges = feats . output . lines . to_array ( ) . reshape ( - <NUM_LIT:1> , <NUM_LIT:3> ) [ : , <NUM_LIT:1> : ] <EOL> for poly in trace_poly ( edges ) : <EOL> vox += rasterize ( epts [ poly ] [ : , : <NUM_LIT:2> ] + [ <NUM_LIT> , <NUM_LIT> ] , shape = shape [ : <NUM_LIT:2> ] [ : : - <NUM_LIT:1> ] ) <EOL> return vox % <NUM_LIT:2> <EOL> if mp : <EOL> from . import mp <EOL> layers = mp . map ( func , range ( shape [ <NUM_LIT:2> ] ) ) <EOL> else : <EOL> layers = map ( func , range ( shape [ <NUM_LIT:2> ] ) ) <EOL> return np . array ( layers ) . T <EOL> def measure_volume ( pts , polys ) : <EOL> from tvtk . api import tvtk <EOL> pd = tvtk . PolyData ( points = pts , polys = polys ) <EOL> mp = tvtk . MassProperties ( input = pd ) <EOL> return mp . volume <EOL> def marching_cubes ( volume , smooth = True , decimate = True , ** kwargs ) : <EOL> import tvtk <EOL> imgdata = tvtk . ImageData ( dimensions = volume . shape ) <EOL> imgdata . point_data . scalars = volume . flatten ( '<STR_LIT:F>' ) <EOL> contours = tvtk . ContourFilter ( input = imgdata , number_of_contours = <NUM_LIT:1> ) <EOL> contours . set_value ( <NUM_LIT:0> , <NUM_LIT:1> ) <EOL> if smooth : <EOL> smoothargs = dict ( number_of_iterations = <NUM_LIT> , feature_angle = <NUM_LIT> , pass_band = <NUM_LIT> ) <EOL> smoothargs . update ( kwargs ) <EOL> contours = tvtk . WindowedSincPolyDataFilter ( input = contours . output , ** smoothargs ) <EOL> if decimate : <EOL> contours = tvtk . QuadricDecimation ( input = contours . output , target_reduction = <NUM_LIT> ) <EOL> contours . update ( ) <EOL> pts = contours . output . points . to_array ( ) <EOL> polys = contours . output . polys . to_array ( ) . reshape ( - <NUM_LIT:1> , <NUM_LIT:4> ) [ : , <NUM_LIT:1> : ] <EOL> return pts , polys </s>
