eic-opticks/analytic/prism.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022 """
0023 import numpy as np
0024
0025
0026 def make_normal( a, b, c, dtype=np.float32):
0027     """
0028               C
0029         C-A  / \
0030          /  /   \
0031            /     \
0032           /       \
0033          A_________B
0034               ->
0035              B-A
0036
0037          n = (B - A) x (C - A)
0038
0039       * normal outwards
0040       * order A, B, C as anti-clockwise as viewed from outside
0041
0042     """
0043     a = np.asarray(a, dtype=dtype)
0044     b = np.asarray(b, dtype=dtype)
0045     c = np.asarray(c, dtype=dtype)
0046
0047     assert a.shape == (3,)
0048     assert b.shape == (3,)
0049     assert c.shape == (3,)
0050
0051     ba = b - a
0052     ca = c - a
0053
0054     n = np.cross( ba, ca )
0055     nn = np.sqrt(np.dot(n,n))
0056
0057     #print "a %40s b %40s c %40s " % (a,b,c)
0058     n /= nn
0059     #print "ba %40s ca %40s n %40s nn %s " % (ba, ca, n, nn )
0060
0061     return n
0062
0063 def make_plane( normal, point, dtype=np.float32 ):
0064     """
0065     Form plane from its normal and any point that lies on the plane
0066     """
0067     normal = np.asarray(normal, dtype=dtype)
0068     point = np.asarray(point, dtype=dtype)
0069
0070     assert normal.shape == (3,)
0071     assert point.shape == (3,)
0072     normal /= np.sqrt(np.dot(normal,normal))
0073
0074     d = np.dot(normal, point)  # distance of point from origin
0075     plane = np.array( [normal[0], normal[1], normal[2], d], dtype=np.float32 )
0076     return plane
0077
0078 def make_plane3( a, b, c, dtype=np.float32):
0079     """
0080     Form plane from three points that lie in the plane, NB the winding
0081     order is anti-clockwise, ie right hand rule around a->b->c
0082     should point outwards
0083     """
0084     n = make_normal( a, b, c, dtype=dtype )
0085     p = make_plane( n, a )
0086     return p
0087
0088
0089
0090 class ConvexPolyhedronSrc(object):
0091     """
0092     Starting from ctor argument with the vertices
0093     collects faces and planes via calls with vertex indices.
0094     """
0095     def __init__(self, **kwa):
0096         self.srcmeta = kwa
0097
0098     def _get_verts(self):
0099         return self._verts
0100
0101     def _set_verts(self, v ):
0102         """
0103         :param v: vertex array of shape (nv, 3)
0104         """
0105         self.dtype = v.dtype
0106
0107         assert len(v.shape) == 2
0108         assert v.shape[0] > 0
0109         assert v.shape[1] == 3
0110         assert len(v) == v.shape[0]
0111
0112         self.nv = len(v)
0113         self._verts = v
0114         self.planes_ = []
0115         self.faces_ = []
0116
0117     verts = property(_get_verts, _set_verts)
0118
0119     def add_quad_ccw(self, i0, i1, i2, i3):
0120         """
0121         *ccw-order, NOT z-order*
0122
0123         ::
0124
0125             105 template <typename T>
0126             106 void NOpenMeshBuild<T>::add_face_(typename T::VertexHandle v0,typename T::VertexHandle v1, typename T::VertexHandle v2, typename T::VertexHandle v3 )
0127             107 {
0128             108    /*
0129             109               3-------2
0130             110               |     . |
0131             111               |   .   |
0132             112               | .     |
0133             113               0-------1
0134             114    */
0135             115
0136             116     add_face_(v0,v1,v2, -1);
0137             117     add_face_(v2,v3,v0, -1);
0138             118 }
0139
0140         """
0141         self(i0,i1,i2)
0142         self(i2,i3,i0)
0143
0144
0145     def __call__(self, i0, i1, i2):
0146         """
0147         :param i0,i1,i2: index of tri vertices
0148
0149         1. vertices are looked up
0150         2. plane through them is constructed and collected into planes_ list
0151         3. indices are collected into faces_ list
0152
0153         """
0154         nv = self.nv
0155         assert i0 < nv
0156         assert i1 < nv
0157         assert i2 < nv
0158
0159         v = self._verts
0160         a = v[i0]
0161         b = v[i1]
0162         c = v[i2]
0163
0164         p = make_plane3( a, b, c, dtype=self.dtype )
0165         self.planes_.append( p )
0166         self.faces_.append( [i0,i1,i2,-1] )
0167
0168
0169     def _get_planes(self):
0170         """
0171         Converts the list of planes into (n,4) array
0172         """
0173         p = np.zeros( (len(self.planes_),4), dtype=self.dtype )
0174         for i in range(len(p)):
0175             p[i] = self.planes_[i]
0176         pass
0177         return p
0178     planes = property(_get_planes)
0179
0180     def _get_faces(self):
0181         """
0182         Converts the list of faces into an np.int32 (n,4) array
0183         """
0184         f = np.zeros( (len(self.faces_),4), dtype=np.int32 )
0185         for i in range(len(f)):
0186             f[i] = self.faces_[i]
0187         pass
0188         return f
0189     faces = property(_get_faces)
0190
0191     def _get_bbox(self):
0192         """
0193         Calulate the minimum and maximum vertex coordinates to provide bounding box
0194         """
0195         v = self._verts
0196         b = np.zeros( (2,3), dtype=self.dtype )
0197         for i in range(3):
0198             b[0,i] = np.min(v[:,i])
0199             b[1,i] = np.max(v[:,i])
0200         pass
0201         return b
0202     bbox = property(_get_bbox)
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212 def make_prism( angle, height, depth, dtype=np.float32, layout=0, crosscheck=True):
0213     #angle = 30
0214     #height = 200
0215     #depth = 200
0216     #dtype = np.float32
0217     #layout = 0
0218     #crosscheck = True
0219     """
0220     Mid line of the symmetric prism spanning along z from -depth/2 to depth/2
0221
0222                            v0
0223                             A  (0,height,0)     Y
0224                            /|\                  |
0225                           / | \                 |
0226                          /  |  \                +---- X
0227                         /   h   \              Z
0228                        /    |    \ (x,y)
0229                       M     |     N
0230                      /      |      \
0231                     L-------O-------R
0232
0233
0234          (-hwidth,0, 0)           (hwidth, 0, 0)
0235
0236
0237     For apex angle 90 degrees, hwidth = height
0238     """
0239
0240
0241     def a_(xyz, dtype=np.float32):
0242         a = np.asarray(xyz, dtype=dtype)
0243         a /= np.sqrt(np.dot(a,a))
0244         return a
0245
0246
0247     src = ConvexPolyhedronSrc(
0248                 src_type="prism",
0249                 src_angle=angle,
0250                 src_height=height,
0251                  src_depth=depth)
0252
0253
0254     hwidth = height*np.tan((np.pi/180.)*angle/2.)
0255     v = np.zeros( (6,3), dtype=dtype)
0256     p = np.zeros( (5,4), dtype=dtype)
0257
0258     if crosscheck:
0259         p2 = np.zeros( (5,4), dtype=dtype)
0260         n1 = np.zeros( (5,3), dtype=dtype)
0261         n2 = np.zeros( (5,3), dtype=dtype)
0262
0263
0264     if layout == 0:
0265         """
0266         Front and back faces
0267
0268                                     v3
0269                                      +
0270                                     / \
0271                            v0      /   \                Y
0272                             +     /     \               |
0273                            / \   /       \              |
0274                           /   \ /         \             |
0275                          /     \           \            +---- X
0276                         /     / \           \          /
0277                        /     /   \           \        /
0278                       /     +-----\-----------+      Z
0279                      /    v4       \          v5
0280                     /               \
0281                    +-----------------+
0282                   v1                 v2
0283
0284                 (x,y)
0285
0286         """
0287
0288         xmin = -hwidth
0289         xmax =  hwidth
0290         ymin = -height/2.
0291         ymax =  height/2.
0292         zmin = -depth/2.
0293         zmax =  depth/2.
0294
0295         v[0] = [      0, ymax,  zmax ]   # front apex
0296         v[1] = [   xmin, ymin,  zmax ]   # front left
0297         v[2] = [   xmax, ymin,  zmax ]   # front right
0298
0299         v[3] = [      0, ymax,  zmin ]   # back apex
0300         v[4] = [   xmin, ymin,  zmin ]   # back left
0301         v[5] = [   xmax, ymin,  zmin ]   # back right
0302
0303
0304         p[0] = make_plane3( v[5], v[3], v[0] )
0305         p[1] = make_plane3( v[1], v[0], v[3] )
0306         p[2] = make_plane3( v[5], v[2], v[1] )
0307         p[3] = make_plane3( v[1], v[2], v[0] )
0308         p[4] = make_plane3( v[4], v[3], v[5] )
0309
0310         src.verts = v
0311
0312         src(5,3,0)
0313         src(1,0,3)
0314         src(5,2,1)
0315         src(1,2,0)
0316         src(4,3,5)
0317
0318
0319         if crosscheck:
0320
0321             n2[0] = a_([ height, hwidth, 0])
0322             n2[1] = a_([-height, hwidth, 0])
0323             n2[2] = a_([ 0, -1,  0])
0324             n2[3] = a_([ 0,  0,  1])
0325             n2[4] = a_([ 0,  0, -1])
0326
0327             # anti-clockwise from outside winding
0328             n1[0] = make_normal( v[5], v[3], v[0] )
0329             n1[1] = make_normal( v[1], v[0], v[3] )
0330             n1[2] = make_normal( v[5], v[2], v[1] )
0331             n1[3] = make_normal( v[1], v[2], v[0] )
0332             n1[4] = make_normal( v[4], v[3], v[5] )
0333
0334             assert np.allclose(n1, n2)
0335
0336             p2[0] = make_plane( n1[0], v[0] )  # +X+Y
0337             p2[1] = make_plane( n1[1], v[0] )  # -X+Y
0338             p2[2] = make_plane( n1[2], v[1] )  # -Y
0339             p2[3] = make_plane( n1[3], v[0] )  # +Z
0340             p2[4] = make_plane( n1[4], v[3] )  # -Z
0341
0342             assert np.allclose(p, p2)
0343         pass
0344
0345
0346
0347     elif layout == 1:
0348         """
0349         Front and back faces
0350
0351                                     v3
0352                                      +
0353                                     / \
0354                            v0      /   \                Z
0355                             +     /     \               |  Y
0356                            / \   /       \              | /
0357                           /   \ /         \             |/
0358                          /     \           \            +---- X
0359                         /     / \           \
0360                        /     /   \           \
0361                       /     +-----\-----------+
0362                      /    v4       \          v5
0363                     /               \
0364                    +-----------------+
0365                   v1                 v2
0366
0367                 (x,z)
0368
0369         """
0370
0371         xmin = -hwidth
0372         xmax =  hwidth
0373         ymin = -depth/2.
0374         ymax =  depth/2.
0375         zmin = -height/2.
0376         zmax =  height/2.
0377
0378         v[0] = [       0,  ymin, zmax ]   # front apex
0379         v[1] = [    xmin,  ymin, zmin ]   # front left
0380         v[2] = [    xmax,  ymin, zmin ]   # front right
0381
0382         v[3] = [       0,  ymax, zmax ]   # back apex
0383         v[4] = [    xmin,  ymax, zmin ]   # back left
0384         v[5] = [    xmax,  ymax, zmin ]   # back right
0385
0386         p[0] = make_plane3( v[5], v[3], v[0] )
0387         p[1] = make_plane3( v[1], v[0], v[3] )
0388         p[2] = make_plane3( v[5], v[2], v[1] )
0389         p[3] = make_plane3( v[1], v[2], v[0] )
0390         p[4] = make_plane3( v[4], v[3], v[5] )
0391
0392         src.verts = v
0393
0394         src(5,3,0)
0395         src(1,0,3)
0396         src(5,2,1)
0397         src(1,2,0)
0398         src(4,3,5)
0399
0400
0401         if crosscheck:
0402
0403             n2[0] = a_([ height, 0, hwidth])  # +X+Z
0404             n2[1] = a_([-height, 0, hwidth])  # -X+Z
0405             n2[2] = a_([ 0, 0,-1])            # -Z
0406             n2[3] = a_([ 0,-1, 0])            # -Y
0407             n2[4] = a_([ 0, 1, 0])            # +Y
0408
0409             # anti-clockwise from outside winding
0410             n1[0] = make_normal( v[5], v[3], v[0] )
0411             n1[1] = make_normal( v[1], v[0], v[3] )
0412             n1[2] = make_normal( v[5], v[2], v[1] )
0413             n1[3] = make_normal( v[1], v[2], v[0] )
0414             n1[4] = make_normal( v[4], v[3], v[5] )
0415
0416             assert np.allclose(n1, n2)
0417
0418             p2[0] = make_plane( n1[0], v[0] )  # +X+Z
0419             p2[1] = make_plane( n1[1], v[0] )  # -X+Z
0420             p2[2] = make_plane( n1[2], v[1] )  # -Z
0421             p2[3] = make_plane( n1[3], v[0] )  # -Y
0422             p2[4] = make_plane( n1[4], v[3] )  # +Y
0423
0424             assert np.allclose(p, n2)
0425         pass
0426     else:
0427         assert 0
0428
0429
0430
0431     bb = np.zeros( (2,3), dtype=dtype )
0432     for i in range(3):
0433         bb[0,i] = np.min(v[:,i])
0434         bb[1,i] = np.max(v[:,i])
0435     pass
0436
0437
0438     planes = src.planes
0439     verts = src.verts
0440     faces = src.faces
0441     bbox = src.bbox
0442
0443     assert np.all( planes == p )
0444     assert np.all( verts == v )
0445     assert np.all( bbox == bb )
0446
0447
0448
0449     return p, v, bb, src
0450
0451
0452
0453
0454
0455 def make_segment( phi0, phi1, sz, sr, dtype=np.float32 ):
0456     """
0457
0458     Prism intended for deltaphi intersecting with
0459     vertices 0 and 3 on the z-axis at -sz/2 and sz/2
0460
0461     ::
0462
0463                5
0464               / \
0465              /   \
0466          sr /     \
0467            /       \
0468           /         \
0469          3-----------4       top plane at z = sz/2
0470               sr
0471
0472
0473                2
0474               / \
0475              /   \
0476          sr /     \
0477            /       \
0478           /         \
0479          0-----------1        base plane at z = -sz/2
0480               sr   (x1,y1)
0481
0482
0483
0484                                        Z
0485                                        |  Y
0486                                        | /
0487                                        |/
0488                                        +------ X
0489
0490     """
0491
0492     src = ConvexPolyhedronSrc(
0493                     src_type="segment",
0494                     src_phi0=phi0,
0495                     src_phi1=phi1,
0496                     src_sz=sz,
0497                     src_sr=sr)
0498
0499
0500
0501     xy_ = lambda phi:np.array([np.cos(phi*np.pi/180.), np.sin(phi*np.pi/180.)], dtype=dtype)
0502
0503     v = np.zeros( (6,3), dtype=dtype)   # verts
0504
0505     x0,y0 = 0.,0.
0506     x1,y1 = sr*xy_(phi0)
0507     x2,y2 = sr*xy_(phi1)
0508
0509     v[0] = [    x0,     y0 , -sz/2. ]
0510     v[1] = [    x1,     y1 , -sz/2. ]
0511     v[2] = [    x2,     y2 , -sz/2. ]
0512
0513     v[3] = [    x0,     y0 ,  sz/2. ]
0514     v[4] = [    x1,     y1 ,  sz/2. ]
0515     v[5] = [    x2,     y2 ,  sz/2. ]
0516
0517
0518
0519     p = np.zeros( (5,4), dtype=dtype)   # planes
0520     p[0] = make_plane3( v[0], v[2], v[1] ) # -Z
0521     p[1] = make_plane3( v[3], v[4], v[5] ) # +Z
0522     p[2] = make_plane3( v[0], v[1], v[3] ) #
0523     p[3] = make_plane3( v[0], v[3], v[5] ) #
0524     p[4] = make_plane3( v[2], v[5], v[4] ) #
0525
0526     src.verts = v
0527
0528     src(0,2,1)   # add planes by referencing vertex indices
0529     src(3,4,5)
0530     src(0,1,3)
0531     src(0,3,5)
0532     src(2,5,4)
0533
0534     b = np.zeros( (2,3), dtype=dtype )  # bbox
0535     for i in range(3):
0536         b[0,i] = np.min(v[:,i])
0537         b[1,i] = np.max(v[:,i])
0538     pass
0539
0540
0541     planes = src.planes
0542     verts = src.verts
0543     bbox = src.bbox
0544
0545     assert np.all( planes == p )
0546     assert np.all( verts == v )
0547     assert np.all( bbox == b )
0548
0549     return p, v, b, src
0550
0551
0552
0553
0554
0555
0556
0557 def make_cubeplanes(hx=0.5, hy=0.5, hz=0.5, dtype=np.float32):
0558     """
0559     z-order verts
0560
0561
0562                   6----------7
0563                  /|         /|
0564                 / |        / |
0565                4----------5  |
0566                |  |       |  |
0567                |  |       |  |         Z
0568                |  2-------|--3         |  Y
0569                | /        | /          | /
0570                |/         |/           |/
0571                0----------1            +------ X
0572
0573     """
0574
0575     src = ConvexPolyhedronSrc(
0576               src_type="cubeplanes",
0577               src_hx=hx,
0578               src_hy=hy,
0579               src_hz=hz
0580            )
0581
0582     v = np.zeros( (8,3), dtype=dtype)   # verts
0583
0584                                 # ZYX
0585     v[0] = [ -hx , -hy , -hz ]  # 000
0586     v[1] = [  hx , -hy , -hz ]  # 001
0587     v[2] = [ -hx ,  hy , -hz ]  # 010
0588     v[3] = [  hx ,  hy , -hz ]  # 011
0589
0590     v[4] = [ -hx , -hy ,  hz ]  # 100
0591     v[5] = [  hx , -hy ,  hz ]  # 101
0592     v[6] = [ -hx ,  hy ,  hz ]  # 110
0593     v[7] = [  hx ,  hy ,  hz ]  # 111
0594
0595
0596     src.verts = v
0597
0598     #src(3,7,5)
0599     #src(0,4,6)
0600     #src(2,6,7)
0601     #src(1,5,4)
0602     #src(5,7,6)
0603     #src(3,1,0)
0604
0605     # visualize the cube from different faces, and read off the ccw cycles
0606
0607     src.add_quad_ccw(0,1,5,4)   # from -Y
0608     src.add_quad_ccw(0,4,6,2)   # from -X
0609     src.add_quad_ccw(0,2,3,1)   # from -Z
0610     src.add_quad_ccw(5,7,6,4)   # from +Z
0611     src.add_quad_ccw(3,7,5,1)   # from +X
0612     src.add_quad_ccw(6,7,3,2)   # from +Y
0613
0614
0615     planes = src.planes
0616     faces = src.faces
0617     verts = src.verts
0618     bbox = src.bbox
0619
0620     return src
0621
0622
0623
0624
0625 def make_trapezoid( z, x1, y1, x2, y2, dtype=np.float32 ):
0626     """
0627     z-order verts
0628
0629
0630                   6----------7
0631                  /|         /|
0632                 / |        / |
0633                4----------5  |
0634                |  |       |  |
0635                |  |       |  |         Z
0636                |  2-------|--3         |  Y
0637                | /        | /          | /
0638                |/         |/           |/
0639                0----------1            +------ X
0640
0641
0642     x1: x length at -z
0643     y1: y length at -z
0644
0645     x2: x length at +z
0646     y2: y length at +z
0647
0648     z:  z length
0649
0650     """
0651     src = ConvexPolyhedronSrc(
0652               src_type="trapezoid",
0653               src_z=z,
0654              src_x1=x1,
0655              src_y1=y1,
0656              src_x2=x2,
0657              src_y2=y2
0658            )         # prefix as BParameters/BList did not supporting multilevel ? Now are moving to BMeta which does
0659
0660     v = np.zeros( (8,3), dtype=dtype)   # verts
0661                                     # ZYX
0662     v[0] = [ -x1/2., -y1/2. , -z/2. ]  # 000
0663     v[1] = [  x1/2., -y1/2. , -z/2. ]  # 001
0664     v[2] = [ -x1/2.,  y1/2. , -z/2. ]  # 010
0665     v[3] = [  x1/2.,  y1/2. , -z/2. ]  # 011
0666
0667     v[4] = [ -x2/2., -y2/2. ,  z/2. ]  # 100
0668     v[5] = [  x2/2., -y2/2. ,  z/2. ]  # 101
0669     v[6] = [ -x2/2.,  y2/2. ,  z/2. ]  # 110
0670     v[7] = [  x2/2.,  y2/2. ,  z/2. ]  # 111
0671
0672     p = np.zeros( (6,4), dtype=dtype)   # planes
0673     p[0] = make_plane3( v[3], v[7], v[5] ) # +X
0674     p[1] = make_plane3( v[0], v[4], v[6] ) # -X
0675     p[2] = make_plane3( v[2], v[6], v[7] ) # +Y
0676     p[3] = make_plane3( v[1], v[5], v[4] ) # -Y
0677     p[4] = make_plane3( v[5], v[7], v[6] ) # +Z
0678     p[5] = make_plane3( v[3], v[1], v[0] ) # -Z
0679
0680     b = np.zeros( (2,3), dtype=dtype )  # bbox
0681     for i in range(3):
0682         b[0,i] = np.min(v[:,i])
0683         b[1,i] = np.max(v[:,i])
0684     pass
0685
0686
0687     src.verts = v
0688
0689     src(3,7,5)
0690     src(0,4,6)
0691     src(2,6,7)
0692     src(1,5,4)
0693     src(5,7,6)
0694     src(3,1,0)
0695
0696     planes = src.planes
0697     faces = src.faces
0698     verts = src.verts
0699     bbox = src.bbox
0700
0701     assert np.all( planes == p )
0702     assert np.all( verts == v )
0703     assert np.all( bbox == b )
0704     assert len(faces) == len(planes)
0705
0706     return p, v, b, src
0707
0708
0709
0710
0711 def make_icosahedron(scale=500., dtype=np.float32):
0712     """
0713     dtype = np.float32
0714
0715     * 20 faces (triangular)
0716     * 12 verts
0717     * 30 edges
0718
0719     Translation of npy-/NTrianglesNPY::icosahedron()
0720
0721     """
0722     CZ = 2./np.sqrt(5)
0723     SZ = 1./np.sqrt(5)
0724
0725     src = ConvexPolyhedronSrc(
0726                     src_type="icosahedron",
0727                    src_scale=scale,
0728                     src_cz=CZ,
0729                     src_sz=SZ)
0730
0731
0732     C1 = np.cos( np.pi*18./180. )
0733     S1 = np.sin( np.pi*18./180. )
0734     C2 = np.cos( np.pi*54./180. )
0735     S2 = np.sin( np.pi*54./180. )
0736
0737     SC = 1.
0738     X1 = C1*CZ
0739     Y1 = S1*CZ
0740     X2 = C2*CZ
0741     Y2 = S2*CZ
0742
0743     SC *= scale
0744     X2 *= scale
0745     Y2 *= scale
0746     SZ *= scale
0747     X1 *= scale
0748     Y1 *= scale
0749     CZ *= scale
0750
0751     vec3_ = lambda x,y,z:np.array([x,y,z], dtype=dtype)
0752
0753     v = np.zeros( (12,3), dtype=dtype)   # verts
0754     Ip0 = v[0] = vec3_(0,0,SC)
0755     Ip1 = v[1] = vec3_(-X2,-Y2,SZ)
0756     Ip2 = v[2] = vec3_( X2,-Y2,SZ)
0757     Ip3 = v[3] = vec3_( X1, Y1,SZ)
0758     Ip4 = v[4] = vec3_(  0, CZ,SZ)
0759     Ip5 = v[5] = vec3_(-X1, Y1,SZ)
0760
0761     p0,p1,p2,p3,p4,p5 = range(0,6)
0762
0763     Im0 = v[6] = vec3_(-X1, -Y1,-SZ)
0764     Im1 = v[7] = vec3_(  0, -CZ,-SZ)
0765     Im2 = v[8] = vec3_( X1, -Y1,-SZ)
0766     Im3 = v[9] = vec3_( X2,  Y2,-SZ)
0767     Im4 = v[10] = vec3_(-X2,  Y2,-SZ)
0768     Im5 = v[11] = vec3_(0,0,-SC)
0769
0770     m0,m1,m2,m3,m4,m5 = range(6,12)
0771
0772
0773     vv = np.hstack( [Ip0,Ip1,Ip2,Ip3,Ip4,Ip5, Im0,Im1,Im2,Im3,Im4,Im5] ).reshape(-1,3)
0774
0775     assert np.all( v == vv )
0776
0777
0778     p = np.zeros( (20,4), dtype=dtype)   # planes
0779
0780     # front pole
0781     p[ 0] = make_plane3(Ip0, Ip1, Ip2)
0782     p[ 1] = make_plane3(Ip0, Ip5, Ip1)
0783     p[ 2] = make_plane3(Ip0, Ip4, Ip5)
0784     p[ 3] = make_plane3(Ip0, Ip3, Ip4)
0785     p[ 4] = make_plane3(Ip0, Ip2, Ip3)
0786
0787     # mid
0788     p[ 5] = make_plane3(Ip1, Im0, Im1)
0789     p[ 6] = make_plane3(Im0, Ip1, Ip5)
0790     p[ 7] = make_plane3(Ip5, Im4, Im0)
0791     p[ 8] = make_plane3(Im4, Ip5, Ip4)
0792     p[ 9] = make_plane3(Ip4, Im3, Im4)
0793     p[10] = make_plane3(Im3, Ip4, Ip3)
0794     p[11] = make_plane3(Ip3, Im2, Im3)
0795     p[12] = make_plane3(Im2, Ip3, Ip2)
0796     p[13] = make_plane3(Ip2, Im1, Im2)
0797     p[14] = make_plane3(Im1, Ip2, Ip1)
0798
0799     # back pole
0800     p[15] = make_plane3(Im3, Im2, Im5)
0801     p[16] = make_plane3(Im4, Im3, Im5)
0802     p[17] = make_plane3(Im0, Im4, Im5)
0803     p[18] = make_plane3(Im1, Im0, Im5)
0804     p[19] = make_plane3(Im2, Im1, Im5)
0805
0806
0807     src.verts = v
0808
0809     # front pole
0810     src(p0,p1,p2)
0811     src(p0,p5,p1)
0812     src(p0,p4,p5)
0813     src(p0,p3,p4)
0814     src(p0,p2,p3)
0815
0816     # mid
0817     src(p1, m0, m1)
0818     src(m0, p1, p5)
0819     src(p5, m4, m0)
0820     src(m4, p5, p4)
0821     src(p4, m3, m4)
0822     src(m3, p4, p3)
0823     src(p3, m2, m3)
0824     src(m2, p3, p2)
0825     src(p2, m1, m2)
0826     src(m1, p2, p1)
0827
0828     # back pole
0829     src(m3, m2, m5)
0830     src(m4, m3, m5)
0831     src(m0, m4, m5)
0832     src(m1, m0, m5)
0833     src(m2, m1, m5)
0834
0835
0836     b = np.zeros( (2,3), dtype=dtype )  # bbox
0837     for i in range(3):
0838         b[0,i] = np.min(v[:,i])
0839         b[1,i] = np.max(v[:,i])
0840     pass
0841
0842     planes = src.planes
0843     faces = src.faces
0844     bbox = src.bbox
0845
0846     assert np.all( planes == p )
0847     assert np.all( bbox == b )
0848     assert len(faces) == len(planes)
0849
0850
0851     return p, v, b, src
0852
0853
0854
0855
0856 def test_prism():
0857     p, v, b, src = make_prism( 45, 400,  400 )
0858
0859 def test_trapezoid():
0860     p, v, b, src = make_trapezoid(z=50.02, x1=100, y1=27, x2=237.2, y2=27 )
0861
0862 def test_icosahedron():
0863     p, v, b, src = make_icosahedron()
0864
0865 def test_segment():
0866     p, v, b, src = make_segment(0,45,100,200)
0867
0868
0869 if __name__ == '__main__':
0870     pass
0871
0872     test_prism()
0873     test_trapezoid()
0874     test_icosahedron()
0875     test_segment()
0876
0877
0878
0879
0880
0881