eic-opticks/ana/geom2d.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 os, json, logging, numpy as np
0024 log = logging.getLogger(__name__)
0025
0026 from opticks.ana.key import keydir
0027 from opticks.ana.prim import Solid
0028
0029 class Geom2d(object):
0030     """
0031     Scratch geometry, for designing flight paths
0032
0033     Note that the GParts are not standardly saved, need to run
0034     an Opticks executable such as OpSnapTest with option --savegparts::
0035
0036        OpSnapTest --savegparts
0037
0038     """
0039     def __init__(self, kd, ridx=0):
0040         fp = open(os.path.join(kd, "cachemeta.json"), "r")
0041         meta = json.load(fp)
0042         gcv = meta["GEOCACHE_CODE_VERSION"]
0043         log.info("GEOCACHE_CODE_VERSION:%s" % gcv)
0044
0045         self.ridx = str(ridx)
0046         self.kd = kd
0047         self.meta = meta
0048
0049         self._pvn = None
0050         self._lvn = None
0051         self.ce = np.load(os.path.join(self.kd, "GNodeLib", "all_volume_center_extent.npy"))
0052
0053         dir_ = os.path.expandvars(os.path.join("$TMP/GParts",self.ridx))
0054         self.d = Solid(dir_, kd)     ## mm0 analytic
0055         self.select_gprim()
0056
0057     def _get_pvn(self):
0058         if self._pvn is None:
0059             self._pvn =  np.loadtxt(os.path.join(self.kd, "GNodeLib/all_volume_PVNames.txt" ), dtype="|S100" )
0060         return self._pvn
0061     pvn = property(_get_pvn)
0062
0063     def _get_lvn(self):
0064         if self._lvn is None:
0065             self._lvn =  np.loadtxt(os.path.join(self.kd, "GNodeLib/all_volume_LVNames.txt" ), dtype="|S100" )
0066         return self._lvn
0067     lvn = property(_get_lvn)
0068
0069
0070     def pvfind(self, pvname_start, encoding='utf-8'):
0071         """
0072         :param pvname_start: string start of PV name
0073         :return indices: array of matching indices in pvname array
0074
0075         Examples::
0076
0077             In [1]: mm0.pvfind("pTarget")
0078             Out[1]: array([67843])
0079
0080         """
0081         return np.flatnonzero(np.char.startswith(self.pvn, pvname_start.encode(encoding)))
0082
0083
0084     def select_gprim(self, names=False):
0085         pp = self.d.prims
0086         sli = slice(0,None)
0087         gprim = []
0088         for p in pp[sli]:
0089             if p.lvName.startswith("sPlane"): continue
0090             if p.lvName.startswith("sStrut"): continue
0091             if p.lvName.startswith("sWall"): continue
0092             if p.numParts > 1: continue     # skip compounds
0093             gprim.append(p)
0094             log.info(repr(p))
0095             vol = p.idx[0]     # global volume index
0096             log.info(self.ce[vol])
0097             #print(str(p))
0098             if names:
0099                 pvn = self.pvn[vol]
0100                 lvn = self.lvn[vol]
0101                 log.info(pvn)
0102                 log.info(lvm)
0103             pass
0104         pass
0105         self.gprim = gprim
0106
0107
0108     def dump(self):
0109         for i,p in enumerate(self.gprim):
0110             assert len(p.parts) == 1
0111             pt = p.parts[0]
0112             print(repr(p))
0113             #print(str(p))
0114             #print(pt)
0115             #print(pt.tran)
0116
0117     def render(self, ax, art3d=None):
0118         sc = 1000
0119         for i,p in enumerate(self.gprim):
0120             assert len(p.parts) == 1
0121             pt = p.parts[0]
0122             sh = pt.as_shape("prim%s" % i, sc=sc )
0123             if sh is None:
0124                print(str(p))
0125                continue
0126             #print(sh)
0127             for pa in sh.patches():
0128                 ax.add_patch(pa)
0129                 if not art3d is None:
0130                     art3d.pathpatch_2d_to_3d(pa, z=0, zdir="y")
0131                 pass
0132             pass
0133         pass
0134
0135
0136
0137
0138
0139 if __name__ == '__main__':
0140
0141     logging.basicConfig(level=logging.INFO)
0142     kd = keydir()
0143     log.info(kd)
0144     assert os.path.exists(kd), kd
0145
0146     os.environ["IDPATH"] = kd    ## TODO: avoid having to do this, due to prim internals
0147
0148     mm0 = Geom2d(kd, ridx=0)
0149
0150     import matplotlib.pyplot as plt
0151     from mpl_toolkits.mplot3d import Axes3D
0152     import mpl_toolkits.mplot3d.art3d as art3d
0153
0154     plt.ion()
0155     fig = plt.figure(figsize=(6,5.5))
0156     ax = fig.add_subplot(111,projection='3d')
0157     plt.title("mm0 geom2d")
0158     sz = 25
0159
0160     ax.set_xlim([-sz,sz])
0161     ax.set_ylim([-sz,sz])
0162     ax.set_zlim([-sz,sz])
0163
0164     ax.set_xlabel("x")
0165     ax.set_ylabel("y")
0166     ax.set_zlabel("z")
0167
0168
0169     mm0.render(ax, art3d=art3d)
0170
0171
0172
0173
0174
0175