ana/pmt/plot.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 Plotting from the serialized PMT analytic geometry data
0023
0024 The very thin CATHODE and BOTTOM are composed of sphere parts (sparts)
0025 with close inner and outer radii.
0026 Inner sparts have parent attribute that points to outer sparts
0027
0028 To plot that, need to clip away the inside.
0029
0030 http://matplotlib.1069221.n5.nabble.com/removing-paths-inside-polygon-td40632.html
0031
0032 Could also just set a width, but thats cheating and the point of the
0033 plotting is to check the parts...
0034
0035
0036 """
0037 import numpy as np
0038 import logging, os
0039 log = logging.getLogger(__name__)
0040
0041 import matplotlib.pyplot as plt
0042 import matplotlib.lines as mlines
0043 import matplotlib.patches as mpatches
0044
0045 from opticks.ana.base import opticks_main
0046 from opticks.sysrap.OpticksCSG import CSG_
0047
0048 #TYPECODE = {'Sphere':1, 'Tubs':2, 'Box':3 }  ## equivalent to pre-unified hardcoded and duplicitous approach
0049 TYPECODE = {'Sphere':CSG_.SPHERE, 'Tubs':CSG_.TUBS, 'Box':CSG_.BOX }
0050
0051
0052 path_ = lambda _:os.path.expandvars("$IDPATH/GMergedMesh/1/%s.npy" % _)
0053
0054
0055 X = 0
0056 Y = 1
0057 Z = 2
0058
0059 ZX = [Z,X]
0060 ZY = [Z,Y]
0061 XY = [X,Y]
0062
0063
0064 class Mesh(object):
0065     def __init__(self):
0066         self.v = np.load(path_("vertices"))
0067         self.f = np.load(path_("indices"))
0068         self.i = np.load(path_("nodeinfo"))
0069         self.vc = np.zeros( self.i.shape[0]+1 )
0070         np.cumsum(self.i[:,1], out=self.vc[1:])
0071     def verts(self, solid):
0072         return self.v[self.vc[solid]:self.vc[solid+1]]
0073
0074 class Sphere(object):
0075     def __init__(self, center, radius):
0076         self.center = center
0077         self.radius = radius
0078     def __repr__(self):
0079         return "Sphere %s %s  " % (repr(self.center), self.radius)
0080
0081     def as_patch(self, axes):
0082         circle = mpatches.Circle(self.center[axes],self.radius)
0083         return circle
0084
0085 class ZTubs(object):
0086     def __init__(self, position, radius, sizeZ):
0087         self.position = position
0088         self.radius = radius
0089         self.sizeZ = sizeZ
0090
0091     def __repr__(self):
0092         return "ZTubs pos %s rad %s sizeZ %s " % (repr(self.position), self.radius, self.sizeZ)
0093
0094     def as_patch(self, axes):
0095         if Z == axes[0]:
0096             width = self.sizeZ
0097             height = 2.*self.radius
0098             botleft = self.position[axes] - np.array([self.sizeZ/2., self.radius])
0099             patch = mpatches.Rectangle(botleft, width, height)
0100         elif Z == axes[1]:
0101             assert 0
0102         else:
0103             patch = mpatches.Circle(self.position[axes],self.radius)
0104         return patch
0105
0106
0107 class Bbox(object):
0108     def __init__(self, min_, max_ ):
0109         self.min_ = np.array(min_)
0110         self.max_ = np.array(max_)
0111         self.dim  = max_ - min_
0112
0113     def as_patch(self, axes):
0114          width = self.dim[axes[0]]
0115          height = self.dim[axes[1]]
0116          botleft = self.min_[axes]
0117          rect = mpatches.Rectangle( botleft, width, height)
0118          return rect
0119
0120     def __repr__(self):
0121         return "Bbox %s %s %s " % (self.min_, self.max_, self.dim )
0122
0123
0124 class Pmt(object):
0125     def __init__(self, path):
0126         path = os.path.expandvars(path)
0127         log.info("loading Pmt from %s " % path)
0128         self.data = np.load(path).reshape(-1,4,4)
0129         self.num_parts = len(self.data)
0130         self.all_parts = range(self.num_parts)
0131         self.partcode = self.data[:,2,3].view(np.int32)
0132         self.partnode = self.data[:,3,3].view(np.int32)
0133         self.index    = self.data[:,1,1].view(np.int32)
0134         self.parent   = self.data[:,1,2].view(np.int32)
0135         self.flags    = self.data[:,1,3].view(np.int32)
0136
0137     def parts(self, solid):
0138         """
0139         :param solid: index of node/solid
0140         :return parts array:
0141         """
0142         pts = np.arange(len(self.partnode))
0143         if solid is not None:
0144             pts = pts[self.partnode == solid]
0145         pass
0146         return pts
0147
0148     def bbox(self, p):
0149         part = self.data[p]
0150         return Bbox(part[2,:3], part[3,:3])
0151
0152     def shape(self, p):
0153         """
0154         :param p: part index
0155         :return shape instance: Sphere or ZTubs
0156         """
0157         code = self.partcode[p]
0158         if code == TYPECODE['Sphere']:
0159             return self.sphere(p)
0160         elif code == TYPECODE['Tubs']:
0161             return self.ztubs(p)
0162         else:
0163             log.warning("Pmt.shape typecode %d not recognized, perhaps an old pre-enum-unification .npy ?" % code )
0164             return None
0165
0166     def sphere(self, p):
0167         """
0168         Creates *Shape* instance from Part data identified by index
0169
0170         :param p: part index
0171         :return Sphere:
0172         """
0173         part = self.data[p]
0174         sp = Sphere( part[0][:3], part[0][3])
0175
0176         log.debug("p %2d sp %s " % (p, repr(sp)))
0177
0178         #if p == 10:
0179         #    sp.center = np.array([ 0.,  0.,  -69.], dtype=np.float32)
0180
0181         return sp
0182
0183     def ztubs(self, p):
0184         """
0185         Creates *ZTubs* instance from Part data index p
0186         """
0187         q0,q1,q2,q3 = self.data[p]
0188         return ZTubs( q0[:3], q0[3], q1[0])
0189
0190
0191 class PmtPlot(object):
0192     def __init__(self, ax, pmt, axes):
0193         self.ax = ax
0194         self.axes = axes
0195         self.pmt = pmt
0196         self.patches = []
0197         self.ec = 'none'
0198         self.edgecolor = ['r','g','b','c','m','y','k']
0199         self.highlight = {}
0200
0201     def color(self, i, other=False):
0202         n = len(self.edgecolor)
0203         idx = (n-i-1)%n if other else i%n
0204         return self.edgecolor[idx]
0205
0206     def plot_bbox(self, parts=[]):
0207         for i,p in enumerate(parts):
0208             bb = self.pmt.bbox(p)
0209             _bb = bb.as_patch(self.axes)
0210             self.add_patch(_bb, self.color(i))
0211
0212     def plot_shape_simple(self, parts=[]):
0213         for i,p in enumerate(parts):
0214             sh = self.pmt.shape(p)
0215             _sh = sh.as_patch(self.axes)
0216             self.add_patch(_sh, self.color(i))
0217
0218     def plot_shape(self, parts=[], clip=True):
0219         log.info("plot_shape parts %r " % parts)
0220         for i,p in enumerate(parts):
0221             is_inner = self.pmt.parent[p] > 0
0222             bb = self.pmt.bbox(p)
0223             _bb = bb.as_patch(self.axes)
0224
0225             ec = self.color(i)
0226             #ec = 'none'
0227             self.add_patch(_bb, ec)
0228
0229             sh = self.pmt.shape(p)
0230             _sh = sh.as_patch(self.axes)
0231             ec = self.color(i,other=True)
0232             fc = self.highlight.get(p, 'none')
0233
0234             if is_inner:
0235                 fc = 'w'
0236
0237             self.add_patch(_sh, ec, fc)
0238             if clip:
0239                 _sh.set_clip_path(_bb)
0240
0241     def add_patch(self, patch, ec, fc='none'):
0242         patch.set_fc(fc)
0243         patch.set_ec(ec)
0244         self.patches.append(patch)
0245         self.ax.add_artist(patch)
0246
0247     def limits(self, sx=200, sy=150):
0248         self.ax.set_xlim(-sx,sx)
0249         self.ax.set_ylim(-sy,sy)
0250
0251
0252
0253
0254
0255 def mug_plot(fig, pmt, pts):
0256     for i, axes in enumerate([ZX,XY]):
0257         ax = fig.add_subplot(1,2,i+1, aspect='equal')
0258         pp = PmtPlot(ax, pmt, axes=axes)
0259         pp.plot_shape(pts, clip=True)
0260         pp.limits()
0261
0262 def clipped_unclipped_plot(fig, pmt, pts):
0263     for i, clip in enumerate([False, True]):
0264         ax = fig.add_subplot(1,2,i+1, aspect='equal')
0265         pp = PmtPlot(ax, pmt, axes=ZX)
0266         pp.plot_shape(pts, clip=clip)
0267         pp.limits()
0268
0269 def solids_plot(fig, pmt, solids=range(5)):
0270
0271     if len(solids)>4:
0272         ny,nx = 3,2
0273     else:
0274         ny,nx = 2,2
0275
0276     for i,solid in enumerate(solids):
0277         pts = pmt.parts(solid)
0278         ax = fig.add_subplot(nx,ny,i+1, aspect='equal')
0279         pp = PmtPlot(ax, pmt, axes=ZX)
0280         pp.plot_shape(pts, clip=True)
0281         pp.limits()
0282     pass
0283
0284
0285 def one_plot(fig, pmt, pts, clip=True, axes=ZX, highlight={}):
0286     ax = fig.add_subplot(1,1,1, aspect='equal')
0287     pp = PmtPlot(ax, pmt, axes=axes)
0288     pp.highlight = highlight
0289     pp.plot_shape(pts, clip=clip)
0290     pp.limits()
0291
0292
0293
0294 if __name__ == '__main__':
0295     logging.basicConfig(level=logging.INFO)
0296
0297     args = opticks_main(apmtidx=2)
0298     apmtpath = args.apmtpath
0299
0300
0301     # 0:4  PYREX
0302     # 4:8  VACUUM
0303     # 8:12 CATHODE
0304     # 12:14 BOTTOM
0305     # 14:15 DYNODE
0306
0307     highlight = {}
0308     highlight[8] = 'r'
0309     highlight[9] = 'r'
0310     highlight[10] = 'r'
0311     highlight[11] = 'b'
0312
0313     ALL, PYREX, VACUUM, CATHODE, BOTTOM, DYNODE = None,0,1,2,3,4
0314
0315     mesh = Mesh()
0316
0317     pmt = Pmt(apmtpath)
0318
0319     fig = plt.figure()
0320
0321     axes = ZX
0322
0323     #solid = CATHODE
0324     #solid = BOTTOM
0325     #solid = DYNODE
0326     solid = ALL
0327
0328     pts = pmt.parts(solid)
0329
0330
0331     #pts = np.arange(8)
0332
0333     #mug_plot(fig, pmt, pts)
0334     #clipped_unclipped_plot(fig, pmt, pts)
0335
0336     #one_plot(fig, pmt, pts, highlight=highlight)
0337     one_plot(fig, pmt, pts, axes=axes, clip=True)
0338
0339     # hmm not possible to split at part level, as those are sub solid
0340     #if mesh:
0341     #    vv = mesh.verts(solid)
0342     #    plt.scatter(vv[:,axes[0]],vv[:,axes[1]],c=vv[:,Y])
0343
0344
0345     #solids_plot(fig, pmt, solids=range(5))
0346
0347
0348     fig.show()
0349     fig.savefig(os.path.expandvars("$TMP/pmtplot.png"))
0350