ana/pmt/gcsg.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 GCSG
0023 ======
0024
0025 * THIS IS SCHEDULED FOR ERADICATION,
0026   ONCE NCSG BINARY-TREE APPROACH HAS PERCOLATED THRU FULL CHAIN
0027
0028 * GCSG uses a first child, last child approach which is
0029   not convenient for CSG evaluation on GPU
0030
0031 * Renamed from "CSG" to "GCSG" to avoid confusion with "NCSG" and dev.csg.csg:CSG
0032   and associate this with corresponding C++ side class ggeo/GCSG.
0033
0034
0035 Where/how GCSG is used
0036 -------------------------
0037
0038 ddpart.py
0039     Parts lists retain association to the basis shapes
0040     they came from via gcsg attributes created
0041     within the Elem.partition methods
0042
0043     ::
0044
0045             350     def partition_intersection(self, material=None):
0046             ...
0047             372         if len(spheres) == 3:
0048             373             pts = self.partition_intersection_3spheres(spheres, material=material)
0049             374         elif len(spheres) == 2:
0050             375             pts = self.partition_intersection_2spheres(spheres, material=material)
0051             376         else:
0052             ...
0053             380         pts.gcsg = GCSG(self, spheres)  # so pts keep hold of reference to the basis Elem shapes list they came from
0054             382         return pts
0055
0056             385     def partition_union_intersection_tubs(self, comps, material=None, verbose=False):
0057             ...
0058             390         ipts = comps[0].partition_intersection()
0059             391         sparts = Part.ascending_bbox_zmin(ipts)
0060             ...
0061             428         rparts = Parts()
0062             429         rparts.extend(sparts)
0063             430         rparts.extend([tpart])
0064             ...
0065             436         rparts.gcsg = GCSG(self, [ipts.gcsg, comps[1]] )
0066             437
0067             438         return rparts
0068
0069
0070     As the parts lists get formed and combined the gcsg attributes
0071     are passed along for the ride.
0072
0073     This complicated way of doing things was needed
0074     in order to retain connection between the parts and the source
0075     volumes (for correct material/boundary/surface assignment for example)
0076     but without having a tree structure to place the info upon. As the
0077     partlist was the primary structure, with no node tree to work with.
0078
0079
0080
0081 analytic.py
0082     at the tail of analytic.py main GPmt.save
0083     invokes GCSG.serialize_list if the Buf instance
0084     has a .gcsg attribute
0085
0086
0087 Buffer Layout
0088 ---------------
0089
0090 ::
0091
0092     In [12]: p.csgbuf.view(np.int32)[:,2:]
0093     Out[12]:
0094
0095                                       typecode, nodeindex, parentindex, base
0096                                       csgindex, nchild,  , firstchild,  lastchild_inc_recursive
0097
0098
0099     array([[[10,  1,  0,  0],        10:Union
0100             [ 0,  2,  1,  5]],       0:idx   2:nc  1:fc 5:lcir         lcir - fc + 1 = 4 != 2 nc
0101
0102                [[20,  0,  0,  0],        20:Intersection
0103                 [ 1,  3,  2,  4]],       1:idx   3:nc  2:fc  4:lcir         lcir - fc + 1 = 3 = nc  (means no recursion, ie immediate primitives)
0104
0105                    [[ 3,  0,  0,  0],        3:sph
0106                     [ 2,  0,  0,  0]],       2:idx
0107
0108                    [[ 3,  0,  0,  0],        3:sph
0109                     [ 3,  0,  0,  0]],       3:idx
0110
0111                    [[ 3,  0,  0,  0],        3:sph
0112                     [ 4,  0,  0,  0]],       4:idx
0113
0114            [[ 4,  0,  0,  0],        4:tubs
0115             [ 5,  0,  0,  0]],       5:idx   (lc of 0:idx)
0116
0117
0118                                       typecode, nodeindex, parentindex, base
0119                                       csgindex, nchild,  , firstchild,  lastchild_inc_recursive
0120
0121            [[10,  2,  1,  0],        10:Union
0122             [ 6,  2,  7, 11]],        6:idx  2:nc 7:fc 11:lcir
0123
0124                [[20,  0,  0,  0],        20:Intersection
0125                 [ 7,  3,  8, 10]],       7:idx fc of 6:idx
0126
0127                    [[ 3,  0,  0,  0],        3:sph
0128                     [ 8,  0,  0,  0]],       8: fc of 7:idx
0129
0130                    [[ 3,  0,  0,  0],        3:sph
0131                     [ 9,  0,  0,  0]],
0132
0133                    [[ 3,  0,  0,  0],        3:sph
0134                     [10,  0,  0,  0]],       10: lc of 7:idx
0135
0136            [[ 4,  0,  0,  0],        4:tubs
0137             [11,  0,  0,  0]],       11: lc of 6:idx
0138
0139
0140
0141                [[10,  3,  2,  0],      10:Union         2:pidx
0142                 [12,  2, 13, 14]],
0143
0144                    [[ 3,  0,  0,  0],      3:sph
0145                     [13,  0,  0,  0]],
0146
0147                    [[ 3,  0,  0,  0],      3:sph
0148                     [14,  0,  0,  0]],
0149
0150
0151                [[ 3,  4,  2,  0],      3:sph           2:pidx
0152                 [15,  0,  0,  0]],
0153
0154                [[ 4,  5,  2,  0],      4:tubs          2:pidx
0155                 [16,  0,  0,  0]]],
0156
0157             dtype=int32)
0158
0159
0160
0161
0162
0163
0164 """
0165 import logging, hashlib, sys, os
0166 import numpy as np
0167 np.set_printoptions(precision=2)
0168 from opticks.ana.nbase import Buf
0169
0170 log = logging.getLogger(__name__)
0171
0172
0173 class Progeny(list):
0174     pass
0175
0176
0177 TYPCODE = {'Union':10, 'Intersection':20, 'Sphere':3, 'Tubs':4 }
0178
0179
0180 class GCSG(object):
0181     """
0182     """
0183     lvnodes = []
0184     def __init__(self, ele, children=[]):
0185         """
0186         :param ele: dd.py Elem instance wrapping an lxml one
0187         :param children: comprised of either base elements or other CSG nodes
0188
0189         Hmm need to retain original node association for material assignment
0190         """
0191         self.ele = ele
0192         self.children = children
0193         self.typ = self.ele.__class__.__name__
0194         self.lv = None
0195         self.node = None
0196
0197
0198     def progeny_count(self):
0199         n = 0
0200         for c in self.children:
0201             if type(c) is GCSG:
0202                 n += c.progeny_count()
0203             else:
0204                 n += 1
0205             pass
0206         return n
0207
0208     def as_gcsg(self):
0209         return [[0,0,0,0],[0,0,0,0],[0,0,0,0],[0,0,0,0]]
0210
0211
0212     @classmethod
0213     def serialize_list(cls, gcsg):
0214         """
0215         :param gcsg: list of GCSG nodes
0216         :return gcsgbuf: Buf instance allowing attributes
0217
0218         Hmm this was an early attempt to serialize a GCSG tree
0219         prior to having the ability to render on GPU.
0220
0221         Is it actually used from C++ side ?
0222         YES, used for Geant4 test geometry construction in cfg4/CMaker etc..
0223
0224         (2,0) typecode
0225         (2,1) nodeindex
0226
0227         (3,0) current offset, ie index of the serialized record
0228         (3,1) number of children
0229         (3,2) index corresponding to first child
0230         (3,3) index corresponding to last child (including recursion)
0231
0232         Formerly Tree.csg_serialize
0233
0234         * array csg nodes and their progeny into flat array, with n elements
0235         * allocate numpy buffer (n,4,4)
0236         * huh, looks like flat is just to find the length
0237
0238         """
0239         flat = []
0240         for cn in gcsg:
0241             flat.extend([cn])
0242             pr = cn.progeny()
0243             flat.extend(pr)
0244         pass
0245         for k,p in enumerate(flat):
0246             log.debug(" %s:%s " % (k, repr(p)))
0247
0248         data = np.zeros([len(flat),4,4],dtype=np.float32)
0249
0250         log.info(" dump gcsg repr, count  %d " % len(gcsg) )
0251         for i,cn in enumerate(gcsg):
0252             print i, repr(cn)
0253         log.info(" dump gcsg type, count  %d " % len(gcsg) )
0254         for i,cn in enumerate(gcsg):
0255             is_gcsg = type(cn).__name__ == "GCSG"
0256             print i, type(cn), type(cn).__name__, "YES" if is_gcsg else "NO"
0257
0258         # curious "type(cn) is GCSG" is failing for : <class 'opticks.ana.pmt.gcsg.GCSG'>
0259         # so check on the name instead
0260
0261         offset = 0
0262         for cn in gcsg:
0263             assert type(cn).__name__ == "GCSG", (cn, type(cn).__name__)
0264             offset = cls.serialize_r(data, offset, cn)
0265         pass
0266
0267         log.info("GCSG.serialize tot len(flat) %s final offset %s " % (len(flat), offset))
0268         assert offset == len(flat)
0269         buf = data.view(Buf)
0270         return buf
0271
0272
0273     @classmethod
0274     def serialize_r(cls, data, offset, obj):
0275         """
0276         :param data: npy buffer to fill
0277         :param offset: primary index of buffer at which to start writing
0278         :param obj: to serialize
0279         :return: offset updated to next place to write
0280         """
0281         pass
0282         nchild = 0
0283         nodeindex = 0
0284         parentindex = 0
0285
0286         if type(obj).__name__ == "GCSG":
0287             log.debug("**serialize offset %s typ %s [%s] (%s)" % (offset,obj.typ,repr(obj),repr(obj.node)))
0288             nchild = len(obj.children)
0289             payload = obj.ele if nchild == 0 else obj    # CSG nodes wrapping single elem, kinda different
0290             if obj.lv is not None:
0291                cls.lvnodes.append(obj.lv)
0292                nodeindex = cls.lvnodes.index(obj.lv) + 1
0293                if obj.node.parent is not None:
0294                    parentindex = cls.lvnodes.index(obj.node.parent.lv) + 1
0295                pass
0296                log.info("serialize nodeindex %s parentindex %s " % (nodeindex, parentindex))
0297                log.debug("serialize nodeindex %s lv %s no %s " % (nodeindex, repr(obj.lv), repr(obj.node)))
0298                log.debug("serialize parentindex %s parent %s " % (parentindex, repr(obj.node.parent)))
0299         else:
0300             payload = obj
0301         pass
0302         base = offset
0303
0304         data[base] = payload.as_gcsg()
0305         data[base].view(np.int32)[2,0] = TYPCODE.get(payload.typ, -1 )   # name -> code using above dict
0306         data[base].view(np.int32)[2,1] = nodeindex
0307         data[base].view(np.int32)[2,2] = parentindex
0308         data[base].view(np.int32)[3,0] = base
0309         offset += 1
0310
0311         if nchild > 0:
0312             data[base].view(np.int32)[3,1] = nchild
0313             data[base].view(np.int32)[3,2] = base + 1
0314             for c in obj.children:
0315                 offset = cls.serialize_r(data, offset, c)
0316             pass
0317             data[base].view(np.int32)[3,3] = offset - 1    # after the recursion
0318         pass
0319         return offset
0320
0321
0322     def progeny(self):
0323         pgs = Progeny()
0324         nchild = len(self.children)
0325
0326         if nchild == 0:
0327             pgs.extend([self.ele])
0328         else:
0329             for c in self.children:
0330                 if type(c) is GCSG:
0331                     pg = c.progeny()
0332                     pgs.extend(pg)
0333                 else:
0334                     pgs.extend([c])
0335                 pass
0336             pass
0337         pass
0338         return pgs
0339
0340
0341
0342     def __repr__(self):
0343         pc = self.progeny_count()
0344         return "GCSG node, lv %s, %s children, %s progeny [%s] " % ( repr(self.lv),len(self.children), pc, repr(self.ele)) + "\n" + "\n".join(map(lambda _:"    " + repr(_), self.children))
0345
0346
0347
0348 if __name__ == '__main__':
0349     pass
0350