eic-opticks/analytic/sc.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 import numpy as np
0022 import os, logging, collections, sys
0023
0024 log = logging.getLogger(__name__)
0025
0026
0027 from opticks.ana.base import opticks_main, expand_, json_load_, json_save_, json_save_pretty_, splitlines_
0028 from opticks.analytic.treebase import Tree
0029 from opticks.analytic.treebuilder import TreeBuilder
0030 from opticks.analytic.gdml import GDML, Boolean
0031 from opticks.analytic.csg import CSG
0032 from opticks.analytic.glm import mdot_, mdotr_,  to_pyline
0033 from opticks.analytic.polyconfig import PolyConfig
0034
0035 def log_info(msg):
0036     sys.stdout.write(msg)
0037
0038
0039 class Mh(object):
0040     def __init__(self, lvIdx, lvName, soName, uri=""):
0041         self.lvIdx = lvIdx
0042         self.lvName = lvName
0043         self.soName = soName
0044         self.extras = dict(lvIdx=self.lvIdx, soName=self.soName, uri=uri )
0045
0046     def __repr__(self):
0047         return "Mh %4d : %30s %s " % (self.lvIdx, self.soName, self.lvName )
0048
0049     def _get_gltf(self):
0050         d = {}
0051         d["name"] = self.lvName
0052         d["extras"] = self.extras
0053         d["primitives"] = [dict(attributes=[])]
0054         return d
0055     gltf = property(_get_gltf)
0056
0057
0058 class Nd(object):
0059     identity = np.eye(4, dtype=np.float32)
0060     suppress_identity = True
0061     def __init__(self, ndIdx, soIdx, transform, boundary, pvname, depth, scene, selected):
0062         """
0063         :param ndIdx: local within subtree nd index, used for child/parent Nd referencing
0064         :param soIdx: local within substree so index, used for referencing to distinct solids/meshes
0065
0066         :param boundary: string with shortnames omat/osur/isur/imat eg "OwsWater///UnstStainlessSteel"
0067
0068         """
0069         self.ndIdx = ndIdx
0070         self.soIdx = soIdx
0071         self.transform = transform
0072
0073         #self.extras = dict(boundary=boundary, pvname=pvname, selected=selected)
0074         #self.extras = dict(boundary=boundary, pvname=pvname)
0075         # got rid of duplicated and unused, to slim the gltf
0076         self.extras = dict(boundary=boundary)
0077
0078         self.name = pvname
0079         self.depth = depth
0080         self.scene = scene
0081
0082         self.children = []
0083         self.parent = -1
0084
0085     def _get_matrix(self):
0086         m = list(map(float,self.transform.ravel()))
0087         if self.suppress_identity:
0088             is_identity = np.all( self.transform == self.identity )
0089             if is_identity:
0090                 m = None
0091             pass
0092         pass
0093         return m
0094     matrix = property(_get_matrix)
0095
0096
0097     brief = property(lambda self:"Nd ndIdx:%3d soIdx:%d nch:%d par:%d matrix:%s " % (self.ndIdx, self.soIdx,  len(self.children), self.parent, self.matrix))
0098
0099
0100     def _get_boundary(self):
0101         return self.extras['boundary']
0102     def _set_boundary(self, bnd):
0103         self.extras['boundary'] = bnd
0104     boundary = property(_get_boundary, _set_boundary)
0105
0106
0107     def _get_parent_node(self):
0108         return self.scene.get_node(self.parent) if self.parent > -1 else None
0109     parent_node = property(_get_parent_node)
0110     p = property(_get_parent_node)
0111
0112     def _get_hierarchy(self):
0113         """ starting from top, going down tree to self"""
0114         ancestors_ = []
0115         p = self.parent_node
0116         while p:
0117             ancestors_.append(p)
0118             p = p.parent_node
0119         pass
0120         return ancestors_[::-1] + [self]
0121     hierarchy = property(_get_hierarchy)
0122
0123     def _get_trs(self):
0124         hier = self.hierarchy
0125         trs_ = map(lambda n:n.transform, hier )
0126         return trs_
0127     trs = property(_get_trs)
0128
0129     gtr_mdot = property(lambda self:mdot_(self.trs))
0130     gtr_mdotr = property(lambda self:mdotr_(self.trs))
0131
0132     gtr_mdot_r = property(lambda self:mdot_(self.trs[::-1]))
0133     gtr_mdotr_r = property(lambda self:mdotr_(self.trs[::-1]))
0134
0135
0136
0137     def __repr__(self):
0138         #indent = ".. " * self.depth
0139         return "%30s %s " % (self.name, self.brief)
0140
0141     def __str__(self):
0142         indent = ".. " * self.depth
0143         return "\n".join([indent + self.brief] + map(repr,map(self.scene.get_node,self.children)))
0144
0145     def find_nodes(self, lvIdx):
0146         nodes = []
0147         def find_nodes_r(node):
0148             if node.mesh.lvIdx == lvIdx:
0149                 nodes.append(node)
0150             pass
0151             for child in node.children:
0152                 find_nodes_r(self.scene.get_node(child))
0153             pass
0154         pass
0155         find_nodes_r(self)
0156         return nodes
0157
0158
0159     def _get_gltf(self):
0160         d = {}
0161         d["mesh"] = self.soIdx
0162         d["name"] = self.name
0163         d["extras"] = self.extras
0164         if len(self.children) > 0:
0165             d["children"] = self.children
0166         pass
0167         m = self.matrix
0168         if m is not None:
0169             d["matrix"] = m
0170         pass
0171         return d
0172     gltf = property(_get_gltf)
0173
0174
0175 class Sc(object):
0176     def __init__(self, maxcsgheight=4):
0177         self.ulv = set()
0178         self.uso = set()
0179         self.nodes = collections.OrderedDict()
0180         self.meshes = collections.OrderedDict()
0181         self.extras = {}
0182         self.maxcsgheight = maxcsgheight
0183         self.translate_node_count = 0
0184         self.add_node_count = 0
0185         self.selected = []
0186
0187
0188     def _get_gltf(self):
0189         root = 0
0190         d = {}
0191         d["scene"] = 0
0192         d["scenes"] = [{ "nodes":[root] }]
0193         d["asset"] = { "version":"2.0", "extras":self.extras }
0194         d["nodes"] = [node.gltf for node in self.nodes.values()]
0195         d["meshes"] = [mesh.gltf for mesh in self.meshes.values()]
0196         return d
0197     gltf = property(_get_gltf)
0198
0199     brief = property(lambda self:"Sc nodes:%d meshes:%d len(ulv):%d len(uso):%d " % (len(self.nodes), len(self.meshes), len(self.ulv), len(self.uso)))
0200
0201     def __repr__(self):
0202          return "\n".join([self.brief])
0203
0204     def __str__(self):
0205          return "\n".join([self.brief] +  map(repr, self.meshes.items()))
0206
0207
0208     def lv2so(self, lvIdx):
0209         """
0210         Convert from an external "mesh" index lvIdx into
0211         local mesh index, using lvIdx identity
0212         """
0213         soIdx = list(self.meshes.iterkeys()).index(lvIdx)
0214         return soIdx
0215
0216     def add_mesh(self, lvIdx, lvName, soName):
0217         if not lvIdx in self.meshes:
0218             self.meshes[lvIdx] = Mh(lvIdx, lvName, soName)
0219             self.meshes[lvIdx].soIdx = self.lv2so(lvIdx)
0220         pass
0221         return self.meshes[lvIdx]
0222
0223     def get_mesh(self, lvIdx):
0224         return self.meshes[lvIdx]
0225
0226     def find_meshes_so(self, pfx):
0227         return filter(lambda mesh:mesh.soName.startswith(pfx),self.meshes.values())
0228
0229     def find_meshes_lv(self, pfx):
0230         return filter(lambda mesh:mesh.lvName.startswith(pfx),self.meshes.values())
0231
0232
0233     def add_node(self, lvIdx, lvName, pvName, soName, transform, boundary, depth, selected):
0234
0235         mesh = self.add_mesh(lvIdx, lvName, soName)
0236         soIdx = mesh.soIdx
0237
0238         ndIdx = len(self.nodes)
0239         name = "ndIdx:%3d,soIdx:%3d,lvName:%s" % (ndIdx, soIdx, lvName)
0240
0241         #log.info("add_node %s " % name)
0242         assert transform is not None
0243
0244         nd = Nd(ndIdx, soIdx, transform, boundary, pvName, depth, self, selected )
0245         nd.mesh = mesh
0246
0247
0248         assert not ndIdx in self.nodes
0249         self.nodes[ndIdx] = nd
0250         return nd
0251
0252     def get_node(self, ndIdx):
0253         return self.nodes[ndIdx]
0254
0255     def get_transform(self, ndIdx):
0256         nd = self.get_node(ndIdx)
0257         return nd.gtr_mdot_r
0258
0259     def add_node_gdml(self, node, depth, debug=False):
0260         """
0261         :param node: treeified ie (PV,LV) collapsed GDML node
0262         :param depth: integer
0263         :return nd: GLTF translation of the input node
0264         """
0265         lvIdx = node.lv.idx
0266         lvName = node.lv.name
0267         pvName = node.pv.name
0268         soName = node.lv.solid.name
0269         transform = node.pv.transform
0270         boundary = node.boundary
0271         nodeIdx = node.index
0272         selected = node.selected
0273
0274         msg = "sc.py:add_node_gdml nodeIdx:%4d lvIdx:%2d soName:%30s lvName:%s " % (nodeIdx, lvIdx, soName, lvName )
0275         #print msg
0276
0277         if debug:
0278             solidIdx = node.lv.solid.idx
0279             self.ulv.add(lvIdx)
0280             self.uso.add(solidIdx)
0281             assert len(self.ulv) == len(self.uso)
0282             sys.stderr.write(msg+"\n" + repr(transform)+"\n")
0283         pass
0284
0285         nd = self.add_node( lvIdx, lvName, pvName, soName, transform, boundary, depth, selected )
0286
0287         ## hmm: why handle csg translation at node level, its more logical to do at mesh level ?
0288         ##      Presumably done here as it is then easy to access the lv ?
0289         ##
0290
0291         if getattr(nd.mesh,'csg',None) is None:
0292             #print msg
0293             csg = self.translate_lv( node.lv, self.maxcsgheight )
0294             nd.mesh.csg = csg
0295             self.translate_node_count += 1
0296
0297             if csg.meta.get('skip',0) == 1:
0298                 log.warning("tlv(%3d): csg.skip as height %2d > %d lvn %s lvidx %s " % (self.translate_node_count, csg.height, self.maxcsgheight, node.lv.name, node.lv.idx ))
0299             pass
0300         pass
0301
0302
0303         if selected:
0304             #log.info("\n\nselected nd %s \n\n%s\n\n" % (nd,  str(nd.mesh.csg.txt) ))
0305             self.selected.append(nd)
0306         pass
0307
0308         return nd
0309
0310
0311     @classmethod
0312     def translate_lv(cls, lv, maxcsgheight, maxcsgheight2=0 ):
0313         """
0314         NB dont be tempted to convert to node here as CSG is a mesh level thing, not node level
0315
0316         :param lv:
0317         :param maxcsgheight:  CSG trees greater than this are balanced
0318         :param maxcsgheight2:  required post-balanced height to avoid skipping
0319
0320         There are many `solid.as_ncsg` implementations, one for each the supported GDML solids,
0321         some of them return single primitives others return boolean composites, some
0322         such as the Polycone invokes treebuilder to provide uniontree composites.
0323
0324         """
0325
0326         if maxcsgheight2 == 0 and maxcsgheight != 0:
0327             maxcsgheight2 = maxcsgheight + 1
0328         pass
0329
0330         solid = lv.solid
0331         log.debug("translate_lv START %-15s %s  " % (solid.__class__.__name__, lv.name ))
0332
0333         rawcsg = solid.as_ncsg()
0334
0335         if rawcsg is None:
0336             err = "translate_lv solid.as_ncsg failed for solid %r lv %r " % ( solid, lv )
0337             log.fatal(err)
0338             rawcsg = CSG.MakeUndefined(err=err,lv=lv)
0339         pass
0340         rawcsg.analyse()
0341
0342         log.debug("translate_lv DONE %-15s height %3d csg:%s " % (solid.__class__.__name__, rawcsg.height, rawcsg.name))
0343
0344         csg = cls.optimize_csg(rawcsg, maxcsgheight, maxcsgheight2 )
0345
0346         polyconfig = PolyConfig(lv.shortname)
0347         csg.meta.update(polyconfig.meta )
0348
0349         # lv.solid.idx gives an incorrect (too high) soIdx ??
0350         csg.meta.update(lvname=lv.name, soname=lv.solid.name, lvIdx=lv.idx, height=csg.height)
0351
0352         ### Nope pvname is not appropriate in the CSG, CSG is a mesh level tink not a node/volume level thing
0353
0354         return csg
0355
0356
0357     @classmethod
0358     def optimize_csg(self, rawcsg, maxcsgheight, maxcsgheight2):
0359         """
0360         :param rawcsg:
0361         :param maxcsgheight:  tree balancing is for height > maxcsgheight
0362         :param maxcsgheight2: error is raised if balanced tree height reamains > maxcsgheight2
0363         :return csg:  balanced csg tree
0364         """
0365         overheight_ = lambda csg,maxheight:csg.height > maxheight and maxheight != 0
0366
0367         is_balance_disabled = rawcsg.is_balance_disabled()
0368
0369         #log.info(" %s %s " % ( is_balance_disabled, rawcsg.name ))
0370
0371         is_overheight = overheight_(rawcsg, maxcsgheight)
0372         if is_overheight:
0373             if is_balance_disabled:
0374                 log.warning("tree is_overheight but marked balance_disabled leaving raw : %s " % rawcsg.name )
0375                 return rawcsg
0376             else:
0377                 log.debug("proceed to balance")
0378         else:
0379             return rawcsg
0380         pass
0381         log.debug("optimize_csg OVERHEIGHT h:%2d maxcsgheight:%d maxcsgheight2:%d %s " % (rawcsg.height,maxcsgheight, maxcsgheight2, rawcsg.name))
0382
0383         rawcsg.positivize()
0384
0385         csg = TreeBuilder.balance(rawcsg)
0386
0387         log.debug("optimize_csg compressed tree from height %3d to %3d " % (rawcsg.height, csg.height ))
0388
0389         #assert not overheight_(csg, maxcsgheight2)
0390         if overheight_(csg, maxcsgheight2):
0391             csg.meta.update(err="optimize_csg.overheight csg.height %s maxcsgheight:%s maxcsgheight2:%s " % (csg.height,maxcsgheight,maxcsgheight2) )
0392         pass
0393
0394         return csg
0395
0396
0397
0398     def add_tree_gdml(self, target, maxdepth=0):
0399         """
0400         :param target: treebase.Node instance, typically the root node
0401
0402         invoked from gdml2gltf_main, notice the two different types of node:
0403
0404         node
0405             input treebase.Node instances, derived from the GDML parse and treeification
0406             to de-stripe from PV-LV-PV-LV-.. to (PV,LV)-(PV,LV)-..
0407             node.children is used to traverse the tree
0408
0409         nd
0410             output sc.Nd instances, which correspond to the GLTF output
0411
0412         """
0413         self.add_node_count = 0
0414         def build_r(node, depth=0):
0415             self.add_node_count += 1
0416             if self.add_node_count % 1000 == 0:
0417                 log.info("add_tree_gdml count %s depth %s maxdepth %s " % (self.add_node_count,depth,maxdepth ))
0418             pass
0419             if maxdepth == 0 or depth < maxdepth:
0420                 nd = self.add_node_gdml(node, depth)
0421                 assert nd is not None
0422                 for child in node.children:
0423                     ch = build_r(child, depth+1)
0424                     if ch is not None:
0425                         ch.parent = nd.ndIdx
0426                         nd.children.append(ch.ndIdx)
0427                     pass
0428                 pass
0429             else:
0430                 nd = None
0431             pass
0432             return nd
0433         pass
0434         log.info("add_tree_gdml START maxdepth:%d maxcsgheight:%d nodesCount:%5d" % (maxdepth, self.maxcsgheight, len(self.nodes)))
0435         #log.info("add_tree_gdml targetNode: %r " % (target))
0436         tg = build_r(target)
0437         log.info("add_tree_gdml DONE maxdepth:%d maxcsgheight:%d nodesCount:%5d tlvCount:%d addNodeCount:%d tgNd:%r " %
0438              (maxdepth, self.maxcsgheight, len(self.nodes),self.translate_node_count, self.add_node_count, tg))
0439         return tg
0440
0441     def save_extras(self, gdir):
0442         gdir = expand_(gdir)
0443         self.dump_extras()
0444         extras_dir = os.path.join( gdir, "extras" )
0445         log.debug("save_extras %s " % extras_dir )
0446         if not os.path.exists(extras_dir):
0447             os.makedirs(extras_dir)
0448         pass
0449         btxt = []
0450         count = 0
0451         for lvIdx, mesh in self.meshes.items():
0452             soIdx = mesh.soIdx
0453             lvdir = os.path.join( extras_dir, "%d" % lvIdx )
0454             uri = os.path.relpath(lvdir, gdir)
0455             mesh.extras["uri"] = uri
0456             mesh.csg.save(lvdir, soIdx, lvIdx)
0457             btxt.append(uri)
0458             count += 1
0459         pass
0460
0461         log.info("save_extras %s  : saved %d " % (extras_dir, count) )
0462
0463         csgtxt_path = os.path.join(extras_dir, "csg.txt")
0464         log.info("write %d lines to %s " % (len(btxt), csgtxt_path))
0465         file(csgtxt_path,"w").write("\n".join(btxt))
0466
0467     def dump_extras(self):
0468         log.info("dump_extras %d " %  len(self.meshes))
0469         for lvIdx, mesh in self.meshes.items():
0470             soIdx = mesh.soIdx
0471             log.debug( "lv %5d so %5d " % (lvIdx, soIdx)  )
0472         pass
0473
0474     def save(self, path, load_check=True, pretty_also=False):
0475         gdir = os.path.dirname(path)
0476         log.info("saving extras in %s " % gdir )
0477         self.save_extras(gdir)    # sets uri for extra external files, so must come before the json gltf save
0478
0479         log.info("saving gltf into %s " % path )
0480         gltf = self.gltf
0481         json_save_(path, gltf)
0482
0483         if pretty_also:
0484             pretty_path = path.replace(".gltf",".pretty.gltf")
0485             log.info("also saving to %s " % pretty_path )
0486             json_save_pretty_(pretty_path, gltf)
0487         pass
0488
0489         if load_check:
0490             gltf2 = json_load_(path)
0491         pass
0492         return gltf
0493
0494
0495
0496
0497
0498     def dump_all(self, lvns):
0499         log.info("dump_all lvns %d " %  len(lvns))
0500         for lvn in lvns:
0501             self.dump(lvn)
0502         pass
0503
0504     def dump(self, lvn):
0505         mss = self.find_meshes_lv(lvn)
0506         assert len(mss) == 1
0507         ms = mss[0]
0508
0509         tree = ms.csg
0510         #tree.analyse()
0511
0512         tree.dump(lvn)
0513
0514         if not tree.is_positive_form():
0515             tree.positivize()
0516             tree.dump(lvn + " (converted to positive form)")
0517         pass
0518
0519         if TreeBuilder.can_balance(tree):
0520             balanced = TreeBuilder.balance(tree)
0521             balanced.dump(lvn + " (TreeBuilder balanced form)")
0522         else:
0523             log.warning("cannot balance")
0524         pass
0525
0526
0527
0528
0529
0530 def gdml2gltf_main( args ):
0531     """
0532     main used by bin/gdml2gltf.py
0533     """
0534     # envvars are set within opticks_main
0535     gdmlpath = os.environ['OPTICKS_GDMLPATH']
0536     gltfpath = os.environ['OPTICKS_GLTFPATH']
0537
0538     #assert gltfpath.startswith("/tmp")
0539
0540     if gltfpath.startswith("/tmp"):
0541         pass
0542     else:
0543         assert gdmlpath.replace('.gdml','.gltf') == gltfpath
0544         assert gltfpath.replace('.gltf','.gdml') == gdmlpath
0545     pass
0546
0547     log.info(" gdmlpath %s " % gdmlpath )
0548     log.info(" gltfpath %s " % gltfpath )
0549
0550     #return 1
0551
0552     log.info("start GDML parse")
0553     gdml = GDML.parse(gdmlpath)
0554
0555     log.info("start treeify")
0556     tree = Tree(gdml.world)
0557
0558     log.info("start apply_selection")
0559     tree.apply_selection(args.query)   # sets node.selected "volume mask"
0560
0561     log.info("start Sc.ctor")
0562     sc = Sc(maxcsgheight=3)
0563
0564     sc.extras["verbosity"] = 1
0565     sc.extras["targetnode"] = 0   # args.query.query_range[0]   # hmm get rid of this ?
0566
0567     log.info("start Sc.add_tree_gdml")
0568
0569     tg = sc.add_tree_gdml( tree.root, maxdepth=0)
0570
0571     log.info("start Sc.add_tree_gdml DONE")
0572
0573     #path = args.gltfpath
0574     gltf = sc.save(gltfpath)
0575
0576     sc.gdml = gdml
0577     sc.tree = tree
0578
0579     Boolean.SaveBuffer()
0580
0581     return sc
0582
0583
0584 def test_range():
0585     #q = "index:3159,depth:1"
0586     q = "range:3158:3160"
0587     if q is not None:
0588         os.environ['OPTICKS_QUERY']=q
0589     pass
0590
0591     args = opticks_main()
0592     sc = gdml2gltf_main( args )
0593
0594     nd = sc.get_node(3159)
0595     tx = sc.get_transform(3159)
0596
0597     print nd.mesh.csg.txt
0598     print to_pyline(nd.gtr_mdot_r, "gtr")
0599
0600     print tx
0601
0602
0603 if __name__ == '__main__':
0604     from opticks.ana.base import opticks_main
0605     args = opticks_main()
0606     sc = gdml2gltf_main( args )
0607
0608
0609