eic-opticks/ana/shape.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 TODO: get this to work with python3
0023 """
0024 import logging, copy
0025
0026 log = logging.getLogger(__name__)
0027
0028 import numpy as np, math
0029
0030 try:
0031     import matplotlib.pyplot as plt
0032     from matplotlib.patches import Rectangle, Circle, Ellipse, PathPatch
0033     import matplotlib.lines as mlines
0034     import matplotlib.path as mpath
0035 except ImportError:
0036     plt = None
0037 pass
0038
0039
0040 def ellipse_closest_approach_to_point( ex, ez, _c ):
0041     """
0042     Ellipse natural frame, semi axes ex, ez.  _c coordinates of point
0043
0044     :param ex: semi-major axis
0045     :param ez: semi-major axis
0046     :param c: xz coordinates of point
0047
0048     :return p: point on ellipse of closest approach to center of torus circle
0049
0050     Closest approach on the bulb ellipse to the center of torus "circle"
0051     is a good point to target for hype/cone/whatever neck,
0052     as are aiming to eliminate the cylinder neck anyhow
0053
0054     equation of RHS torus circle, in ellipse frame
0055
0056         (x - R)^2 + (z - z0)^2 - r^2 = 0
0057
0058     equation of ellipse
0059
0060         (x/ex)^2 + (z/ez)^2 - 1 = 0
0061
0062     """
0063     c = np.asarray( _c )   # center of RHS torus circle
0064     assert c.shape == (2,)
0065
0066     t = np.linspace( 0, 2*np.pi, 1000000 )
0067     e = np.zeros( [len(t), 2] )
0068     e[:,0] = ex*np.cos(t)
0069     e[:,1] = ez*np.sin(t)   # 1M parametric points on the ellipse
0070
0071     p = e[np.sum(np.square(e-c), 1).argmin()]   # point on ellipse closest to c
0072     return p
0073
0074
0075
0076 def ellipse_points( xy=[0,-5.], ex=254., ez=190., n=1000 ):
0077     """
0078     :param ec: center of ellipse
0079     :param ex: xy radius of ellipse
0080     :param ez: z radius of ellipse
0081     :param n: number of points
0082     :return e: array of shape (n,2) of points on the ellipse
0083     """
0084     t = np.linspace( 0, 2*np.pi, n )
0085     e = np.zeros([len(t), 2])
0086     e[:,0] = ex*np.cos(t) + xy[0]
0087     e[:,1] = ez*np.sin(t) + xy[1]
0088     return e
0089
0090 def circle_points( xy=[0,0], tr=80, n=1000 ):
0091     """
0092     :param tc: center of circle
0093     :param tr: radius of circle
0094     :param n: number of points
0095     :return c: array of shape (n,2) of points on the circle
0096     """
0097     t = np.linspace( 0, 2*np.pi, n )
0098     c = np.zeros([len(t), 2])
0099     c[:,0] = tr*np.cos(t) + xy[0]
0100     c[:,1] = tr*np.sin(t) + xy[1]
0101     return c
0102
0103
0104 def points_inside_circle(points, center, radius):
0105     """
0106     :param points: (n,2) array of points
0107     :param center: (2,) coordinates of circle center
0108     :param radius:
0109     :return mask: boolean array of dimension (n,2) indicating if points are within the circle
0110     """
0111     return np.sqrt(np.sum(np.square(points-center),1)) - radius < 0.
0112
0113
0114 def ellipse_points_inside_circle():
0115
0116     tc = np.array([torus_x,torus_z])
0117     tr = m4_torus_r
0118     e = ellipse_points( xy=[0,-5.], ex=254., ez=190., n=1000000 )
0119
0120
0121
0122
0123
0124
0125
0126 class X(object):
0127     def __init__(self, root):
0128         self.root = root
0129
0130     def __repr__(self):
0131         return "\n".join( map(repr, self.constituents()))
0132
0133     def find(self, shape):
0134         return self.root.find(shape)
0135
0136     def find_one(self, shape):
0137         ff = self.root.find(shape)
0138         assert len(ff) == 1
0139         return ff[0]
0140
0141     def constituents(self):
0142         return self.root.constituents()
0143
0144
0145     def replacement_cons(self):
0146         """
0147         """
0148         i = self.find_one("STorus")
0149         r = i.param[0]
0150         R = i.param[1]
0151
0152         d = self.find_one("SEllipsoid")
0153         ex = d.param[0]
0154         ez = d.param[1]
0155
0156         print("r %s R %s ex %s ez %s " % (r,R,ex,ez))
0157         print(" SEllipsoid d.xy %s " % repr(d.xy) )
0158         print(" STorus     i.xy %s " % repr(i.xy) )
0159
0160         z0 = i.xy[1]   # torus z-plane in ellipsoid frame
0161
0162         p = ellipse_closest_approach_to_point( ex, ez, [R,z0] )    # [R,z0] is center of torus circle
0163
0164         pr, pz = p    # at torus/ellipse closest point : no guarantee of intersection
0165         print(" ellipse closest approach to torus  %s " % repr(p) )
0166
0167         r2 = pr
0168         r1 = R - r
0169         mz = (z0 + pz)/2.   # mid-z cone coordinate (ellipsoid frame)
0170         hz = (pz - z0)/2.   # cons half height
0171
0172         f = SCons( "f", [r1,r2,hz] )
0173         B = np.array( [0, mz] )
0174
0175         print(" replacment SCons %s offset %s " % (repr(f),repr(B)))
0176
0177         return f, B
0178
0179
0180     def spawn_rationalized(self):
0181         """
0182
0183         ::
0184
0185                    UnionSolid
0186                    /         \
0187            Ellipsoid          Subtraction
0188                               /         \
0189                             Tubs        Torus
0190
0191
0192                    UnionSolid
0193                    /         \
0194            Ellipsoid           Cons
0195
0196
0197
0198         """
0199         name = self.__class__.__name__
0200
0201         x = copy.deepcopy(self)
0202
0203         # establish expectations for tree
0204         e = x.find_one("SEllipsoid")
0205         t = x.find_one("STorus")
0206         ss = t.parent
0207         assert ss is not None and ss.shape == "SSubtractionSolid"
0208         us = ss.parent
0209         assert us is not None and us.shape == "SUnionSolid"
0210         assert us.left is not None and us.left == e and us.right == ss and ss.right == t
0211         assert us.right is not None and us.right == ss
0212
0213
0214         if name == "x018":   # cathode vacuum cap
0215             assert x.root.shape == "SIntersectionSolid"
0216             x.root = e
0217             e.parent = None
0218         elif name == "x019":  # remainder vacuum
0219             assert x.root.shape == "SSubtractionSolid"
0220             left = x.root.left
0221             assert left.shape == "SUnionSolid"
0222             left.parent = None
0223             x.root = left
0224         else:
0225             pass
0226         pass
0227
0228         if name in ["x019","x020","x021"]:
0229             # calculate the parameters of the replacement cons
0230             cons, offset =  x.replacement_cons()
0231
0232             # tree surgery : replacing the right child of UnionSolid
0233             us.right = cons
0234             cons.parent = us
0235             cons.ltransform = offset
0236         pass
0237
0238         return x
0239
0240
0241
0242 class Shape(object):
0243     """
0244     matplotlib patches do not support deferred placement it seems,
0245     so do that here
0246     """
0247     KWA = dict(fill=False)
0248     dtype = np.float64
0249
0250     PRIMITIVE = ["SEllipsoid","STubs","STorus", "SCons", "SHype", "SBox", "SPolycone"]
0251     PRIMITIVE2 = ["Box", "Cylinder", "Tubs" ]
0252     ALL_PRIMITIVE = PRIMITIVE + PRIMITIVE2
0253
0254     COMPOSITE = ["SUnionSolid", "SSubtractionSolid", "SIntersectionSolid"]
0255     ALL = ALL_PRIMITIVE + COMPOSITE
0256
0257     def __repr__(self):
0258         return "%s : %20s : %s : %s " % (
0259                       self.name,
0260                       self.shape,
0261                       repr(self.ltransform),
0262                       repr(self.param)
0263                      )
0264
0265
0266     def __init__(self, name, param, **kwa ):
0267         shape = self.__class__.__name__
0268
0269         if not shape in self.ALL:
0270             log.error("shape class name %s is not in the list %s " % ( shape, str(self.ALL)))
0271         pass
0272         assert shape in self.ALL
0273         primitive = shape in self.ALL_PRIMITIVE
0274         composite = shape in self.COMPOSITE
0275
0276         d = self.KWA.copy()
0277         d.update(kwa)
0278         self.kwa = d
0279
0280         self.name = name
0281         self.shape = shape
0282         self.param = param
0283
0284         self.parent = None
0285         self.ltransform = None
0286         self.left = None
0287         self.right = None
0288
0289         if composite:
0290             left = self.param[0]
0291             right = self.param[1]
0292             right.ltransform = self.param[2]
0293
0294             left.parent = self
0295             right.parent = self
0296
0297             self.left = left
0298             self.right = right
0299         pass
0300
0301
0302     is_primitive = property(lambda self:self.left is None and self.right is None)
0303     is_composite = property(lambda self:self.left is not None and self.right is not None)
0304
0305     def _get_xy(self):
0306         """
0307         Assumes only translations, adds the node.ltransform obtained by following
0308         parent links up the tree of shapes.
0309
0310
0311
0312                            a                                            Intersection
0313                                                                        /            \
0314                       b             m(D)                          Union             m:Tubs
0315                                                                   /    \
0316                 c          k(C)                             Union       Tubs
0317                                                            /     \
0318             d     f(B)                            Ellipsoid   Subtraction
0319                                                               /          \
0320                 g(B)  i(B+A)                                Tubs         Torus
0321
0322
0323         """
0324         xy = np.array([0,0], dtype=self.dtype )
0325         node = self
0326         while node is not None:
0327             if node.ltransform is not None:
0328                 log.debug("adding ltransform %s " %  node.ltransform)
0329                 xy += node.ltransform
0330             pass
0331             node = node.parent
0332         pass
0333         return xy
0334
0335     xy = property(_get_xy)
0336
0337     def constituents(self):
0338         if self.is_primitive:
0339             return [self]
0340         else:
0341             assert self.is_composite
0342             cts = []
0343             cts.extend( self.left.constituents() )
0344             cts.extend( self.right.constituents() )
0345             return cts
0346         pass
0347
0348     def find(self, shape):
0349         cts = self.constituents()
0350         return filter( lambda ct:ct.shape == shape, cts )
0351
0352     def patches(self):
0353         """
0354         Positioning is relying on self.xy of the primitives
0355         with nothing being passed into composites.
0356
0357         For composites self.param[2] is the local right transform
0358         """
0359         if self.shape == "SEllipsoid":
0360             return self.make_ellipse( self.xy, self.param, **self.kwa )
0361         elif self.shape == "STubs" or self.shape == "Tubs":
0362             return self.make_rect( self.xy, self.param, **self.kwa)
0363         elif self.shape == "STorus":
0364             return self.make_torus( self.xy, self.param, **self.kwa)
0365         elif self.shape == "SCons":
0366             return self.make_cons( self.xy, self.param, **self.kwa)
0367         elif self.shape == "SHype":
0368             return self.make_hype( self.xy, self.param, **self.kwa)
0369         elif self.shape == "SBox" or self.shape == "Box":
0370             return self.make_rect( self.xy, self.param, **self.kwa)
0371         elif self.shape == "SPolycone":
0372             return self.make_polycone( self.xy, self.param, **self.kwa)
0373         else:
0374             if not self.is_composite:
0375                 log.error("shape :%s: not handled in patches()" % self.shape )
0376             pass
0377             assert self.is_composite
0378             pts = []
0379             pts.extend( self.left.patches() )
0380             pts.extend( self.right.patches() )
0381             return pts
0382         pass
0383
0384     @classmethod
0385     def create(cls, pt ):
0386         pass
0387
0388     @classmethod
0389     def make_rect(cls, xy , wh, **kwa ):
0390         """
0391         :param xy: center of rectangle
0392         :param wh: halfwidth, halfheight
0393         """
0394         ll = ( xy[0] - wh[0], xy[1] - wh[1] )
0395         return [Rectangle( ll,  2.*wh[0], 2.*wh[1], **kwa  )]
0396
0397     @classmethod
0398     def make_ellipse(cls, xy , param, **kwa ):
0399         return [Ellipse( xy,  width=2.*param[0], height=2.*param[1], **kwa  )]
0400
0401     @classmethod
0402     def make_circle(cls, xy , radius, **kwa ):
0403         return [Circle( xy,  radius=radius, **kwa  )]
0404
0405
0406     @classmethod
0407     def make_torus(cls, xy, param, **kwa ):
0408         r = param[0]
0409         R = param[1]
0410
0411         pts = []
0412         lhs = cls.make_circle( xy + [-R,0], r, **kwa)
0413         rhs = cls.make_circle( xy + [+R,0], r, **kwa)
0414         pts.extend(lhs)
0415         pts.extend(rhs)
0416         return pts
0417
0418     @classmethod
0419     def make_pathpatch(cls, xy, vtxs, **kwa ):
0420         """see analytic/pathpatch.py"""
0421         Path = mpath.Path
0422         path_data = []
0423
0424         for i, vtx in enumerate(vtxs):
0425             act = Path.MOVETO if i == 0 else Path.LINETO
0426             path_data.append( (act, (vtx[0]+xy[0], vtx[1]+xy[1])) )
0427         pass
0428         path_data.append( (Path.CLOSEPOLY, (vtxs[0,0]+xy[0], vtxs[0,1]+xy[1])) )
0429         pass
0430
0431         codes, verts = zip(*path_data)
0432         path = Path(verts, codes)
0433         patch = PathPatch(path, **kwa)
0434         return [patch]
0435
0436     @classmethod
0437     def make_cons(cls, xy , param, **kwa ):
0438         """
0439         (-r2,z2)      (r2,z2)
0440               1---------2
0441                \       /
0442                 0 ... 3
0443         (-r1,z1)     (r1,z1)
0444         """
0445         r1 = param[0]
0446         r2 = param[1]
0447         hz = param[2]
0448
0449         z2 =  hz + xy[1]
0450         z1 = -hz + xy[1]
0451
0452         vtxs = np.zeros( (4,2) )
0453         vtxs[0] = ( -r1, z1)
0454         vtxs[1] = ( -r2,  z2)
0455         vtxs[2] = (  r2,  z2)
0456         vtxs[3] = (  r1,  z1)
0457         return cls.make_pathpatch( xy, vtxs, **kwa )
0458
0459     @classmethod
0460     def make_polycone(cls, xy , param, **kwa ):
0461         """
0462         """
0463         zp = param
0464         nz = len(zp)
0465         assert zp.shape == (nz, 3), zp
0466         assert nz > 1 , zp
0467
0468         rmin = zp[:,0]
0469         rmax = zp[:,1]
0470         z = zp[:,2]
0471
0472         vtxs = np.zeros( (2*nz,2) )
0473         for i in range(nz):
0474             vtxs[i] = ( -rmax[i], z[i] )
0475             vtxs[2*nz-i-1] = ( rmax[i], z[i] )
0476         pass
0477         log.debug(" xy : %r " %  xy )
0478         return cls.make_pathpatch( xy, vtxs, **kwa )
0479
0480
0481     @classmethod
0482     def make_hype(cls, xy , param, **kwa ):
0483         """
0484              4----------- 5
0485               3          6
0486                2        7
0487               1          8
0488              0 ---------- 9
0489
0490            sqrt(x^2+y^2) =   r0 * np.sqrt(  (z/zf)^2  +  1 )
0491
0492         """
0493         r0 = param[0]
0494         stereo = param[1]
0495         hz = param[2]
0496         zf = r0/np.tan(stereo)
0497
0498         r_ = lambda z:r0*np.sqrt( np.square(z/zf) + 1. )
0499
0500         nz = 20
0501         zlhs = np.linspace( -hz, hz, nz )
0502         zrhs = np.linspace(  hz, -hz, nz )
0503
0504         vtxs = np.zeros( (nz*2,2) )
0505
0506         vtxs[:nz,0] = -r_(zlhs) + xy[0]
0507         vtxs[:nz,1] = zlhs + xy[1]
0508
0509         vtxs[nz:,0] = r_(zrhs) + xy[0]
0510         vtxs[nz:,1] = zrhs + xy[1]
0511
0512         return cls.make_pathpatch( xy, vtxs, **kwa )
0513
0514
0515 class SEllipsoid(Shape):pass
0516 class STubs(Shape):pass
0517 class STorus(Shape):pass
0518 class SCons(Shape):pass
0519 class SHype(Shape):pass
0520 class SPolycone(Shape):pass
0521 class SUnionSolid(Shape):pass
0522 class SSubtractionSolid(Shape):pass
0523 class SIntersectionSolid(Shape):pass
0524
0525
0526
0527 if __name__ == '__main__':
0528     pass
0529
0530
0531
0532
0533
0534
0535