eic-opticks/ana/prim.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 prim.py
0023 =========
0024
0025 Loads a list of primitives from a GParts persisted directory.
0026 Used for debugging things such as missing transforms.
0027
0028 See: notes/issues/OKX4Test_partBuffer_difference.rst
0029
0030 TODO: consolidate ana/prim.py with dev/csg/GParts.py
0031
0032 """
0033 import sys, os, logging, numpy as np
0034 log = logging.getLogger(__name__)
0035 from opticks.sysrap.OpticksCSG import CSG_
0036 from opticks.ana.blib import BLib
0037 from opticks.ana.mesh import Mesh
0038
0039 from opticks.ana.shape import Shape
0040
0041 class Ellipsoid(Shape):pass
0042 class Tubs(Shape):pass
0043 class Torus(Shape):pass
0044 class Cons(Shape):pass
0045 class Hype(Shape):pass
0046 class Box(Shape):pass
0047
0048 class UnionSolid(Shape):pass
0049 class SubtractionSolid(Shape):pass
0050 class IntersectionSolid(Shape):pass
0051
0052
0053 class Part(object):
0054     """
0055     Parts are CSG constituents, aka nodes of the CSG trees that make up each solid
0056     """
0057     part_idx = 0
0058
0059     def __init__(self, part, trans, d, prim):
0060         """
0061         Part instances are created within the parent Prim instance
0062         by mapping this Part ctor over elements of the parts (nodes)
0063         and global transforms arrays.
0064
0065         :param part: single csg node of shape (4,4)
0066         :param trans: 1 or more transforms, shape (ntran,4,4)
0067         :param d: Dir instance
0068         """
0069
0070         sidx = d.sidx
0071         self.d = d
0072         self.prim = prim
0073         self.sidx = sidx
0074
0075         assert part.shape == (4,4), part
0076         assert trans.shape[1:] == (4,4), trans
0077         assert trans.shape[0] < 16, trans
0078         assert d.__class__.__name__ == 'Solid'   # formerly 'Dir'
0079         ntran = trans.shape[0]
0080         assert ntran > 0
0081
0082         log.debug( "Part : trans.shape %s " % repr(trans.shape))
0083         #print("trans", trans)
0084
0085         f = part.view(np.float32)
0086         u = part.view(np.uint32)
0087
0088         fc = part.copy().view(np.float32)
0089         fc[1][2] = 0  # scrub boundary in copy, as it is known discrepant : due to lack of surfaces
0090
0091
0092         tc = u[2][3]          ## typecode eg CSG_.UNION
0093         tcn = CSG_.desc(tc)   ## typename
0094         bnd = u[1][2]         ## boundary index
0095
0096         # check the complements, viewing as float otherwise lost "-0. ok but not -0"
0097         comp = np.signbit(f[3,3])
0098          # shift everything away, leaving just the signbit
0099         comp2 = u[3,3] >> 31
0100         assert comp == comp2
0101
0102         # recover the gtransform index by getting rid of the complement signbit
0103         gt = u[3][3] & 0x7fffffff
0104
0105         if tc in [CSG_.DIFFERENCE, CSG_.INTERSECTION, CSG_.UNION]:
0106             assert gt == 0, "operators are expected to not have a gtransform"
0107         elif tc == CSG_.ZERO:
0108             assert gt == 0, "zeros are not expected to have a gtransform"
0109         else:
0110             assert gt > 0, "primitives are expected to have a gtransform "
0111
0112             #log.info(" sidx:%2d part_idx:%5d ntran:%2d gt:%2d tcn:%s " % ( sidx, self.__class__.part_idx, ntran, gt, tcn ))
0113
0114             if not gt <= ntran:
0115                 log.warn("1-based gt expected to be <= ntran (local index, not global)  gt:%d ntran:%d " % (gt,ntran))
0116             pass
0117             #assert gt <= ntran
0118             if ntran > 5:
0119                 pass
0120                 #log.info(" part_idx:%5d ntran:%2d gt:%2d tcn:%s " % ( self.__class__.part_idx, ntran, gt, tcn ))
0121             pass
0122         pass
0123
0124
0125         self.f = f
0126         self.fc = fc
0127
0128         self.tc = tc
0129         self.tcn = tcn
0130         self.comp = comp
0131         self.gt = gt        # 1-based gtransform pointer into trans
0132
0133         self.bnd = bnd
0134         self.bname = d.blib.bname(bnd)
0135
0136         try:
0137             tran = trans[gt-1] if gt > 0 else np.eye(4)
0138         except IndexError:
0139             log.error("trans issue gt %s trans.shape %s " % ( gt, repr(trans.shape)))
0140             tran = np.eye(4)
0141         pass
0142
0143         self.tran = tran
0144
0145         self.idx = self.__class__.part_idx
0146         self.__class__.part_idx += 1
0147
0148
0149
0150     def __repr__(self):
0151         return "    Part %1s%2s %2s  %15s   %3d %25s   tz:%10.3f    %s  " % ( "!" if self.comp else " ", self.tc, self.gt, self.tcn, self.bnd, self.bname, self.tz, self.detail() )
0152
0153     def maxdiff(self, other):
0154         """
0155         :param other: Part instance
0156         :return float: max absolute difference between float param values of the CSG constituent part
0157         (actually the boundary index is excluded by comparing a copy and scrubbing that)
0158         """
0159         return np.max( np.abs(self.fc - other.fc) )
0160
0161     r = property(lambda self:self.f[0][3])
0162     tz = property(lambda self:self.tran[3][2])
0163
0164     xbox = property(lambda self:self.f[0][0])
0165     ybox = property(lambda self:self.f[0][1])
0166     zbox = property(lambda self:self.f[0][2])
0167
0168     r1co = property(lambda self:self.f[0][0])
0169     z1co = property(lambda self:self.f[0][1])
0170     r2co = property(lambda self:self.f[0][2])
0171     z2co = property(lambda self:self.f[0][3])
0172
0173     z1 = property(lambda self:self.f[1][0])  # cy or zsphere
0174     z2 = property(lambda self:self.f[1][1])
0175     r1 = property(lambda self:self.f[0][2])
0176     r2 = property(lambda self:self.f[0][3])
0177     dz = property(lambda self:self.z2 - self.z1)
0178     dr = property(lambda self:self.r2 - self.r1)
0179
0180
0181     def as_shape(self, name, sc):
0182         if self.tc == CSG_.BOX3:
0183             sh = Box(name, [self.xbox/sc, self.zbox/sc ] )
0184         elif self.tc == CSG_.CYLINDER:
0185             sh = Tubs(name, [self.r/sc, abs(self.z1/sc) ] )
0186         elif self.tc == CSG_.TORUS:
0187             sh = Ellipsoid(name, [self.r/sc, self.r/sc ] )
0188         else:
0189             sh = None
0190         pass
0191         if not sh is None:
0192             sh.ltransform = [0, self.tz/sc ]
0193         pass
0194         return sh
0195
0196     def detail(self):
0197         tz = self.tz
0198         if self.tc == CSG_.ZSPHERE:
0199             msg = " r: %10.3f z1:%10.3f z2:%10.3f " % ( self.r, self.z1, self.z2 )
0200         elif self.tc == CSG_.SPHERE:
0201             msg = " r: %10.3f " % ( self.f[0][3]  )
0202         elif self.tc == CSG_.CYLINDER:
0203             msg = "   z1:%10.3f z2:%10.3f r :%10.3f               z1+tz:%10.3f z2+tz:%10.3f" % ( self.z1, self.z2, self.r, self.z1 + tz, self.z2 + tz)
0204         elif self.tc == CSG_.CONE:
0205             msg = "   z1:%10.3f z2:%10.3f r1:%10.3f r2:%10.3f z1+tz:%10.3f z2+tz:%10.3f" % ( self.z1co, self.z2co, self.r1co, self.r2co, self.z1co+tz, self.z2co+tz )
0206         elif self.tc == CSG_.BOX3:
0207             msg = "   x:%10.3f y:%10.3f z:%10.3f  " % ( self.xbox, self.ybox, self.zbox  )
0208         else:
0209             msg = ""
0210         pass
0211         return msg
0212
0213
0214
0215 class Prim(object):
0216     """
0217
0218     In [21]: gp[2].primbuf
0219     Out[21]:
0220     array([[ 0,  1,  0,  0],
0221            [ 1,  7,  1,  0],
0222            [ 8,  7,  5,  0],
0223            [15,  7,  8,  0],
0224            [22,  1, 11,  0],
0225            [23,  7, 12,  0]], dtype=int32)
0226              |   |   |   |
0227      partOffset  |   |   |
0228                  |   |   |
0229            numParts  |   |
0230                      |   |
0231              tranOffset  |
0232                          |
0233                  planOffset
0234     """
0235     def __init__(self, primIdx, prim, d):
0236         """
0237         """
0238         assert primIdx > -1 and primIdx < 10000, primIdx
0239         assert prim.shape == (4,), "unexpected prim.shape %s " % repr(prim.shape)
0240
0241         self.sidx = d.sidx
0242         self.primIdx = primIdx
0243         self.prim = prim
0244
0245         idx = d.idxbuf[primIdx] if d.idxbuf is not None else -np.ones(4,dtype=np.uint32)
0246         lvIdx = idx[2]
0247
0248         self.idx = idx
0249         self.lvIdx = lvIdx
0250
0251         self.lvName = d.ma.idx2name.get(lvIdx, "-") if d.ma is not None else "--"
0252
0253         partOffset, numParts, tranOffset, planOffset = prim
0254         numTran = d.ntran[primIdx]
0255
0256         parts_ = d.partbuf[partOffset:partOffset+numParts]
0257         trans_ = d.tranbuf[tranOffset:tranOffset+numTran,0]  # eg shape (2, 4, 4)  plucking the first from the t,v,q triplet of transforms
0258
0259         self.parts_ = parts_
0260         self.trans_ = trans_    ## without the python class wrapping
0261
0262         self.parts = list(map(lambda _:Part(_,trans_,d, self), parts_))   ## note that every part gets passed all the trans_ need to use the gt to determine which one to use
0263
0264         self.partOffset = partOffset
0265         self.numParts= numParts
0266         self.numTran = numTran
0267         self.tranOffset = tranOffset
0268         self.planOffset = planOffset
0269         self.d = d
0270
0271     def maxdiff(self, other):
0272         """
0273         :return float: max difference over the constituent parts, from Part.maxdiff
0274         """
0275         assert len(self.parts) == len(other.parts)
0276         return max(map( lambda ab:ab[0].maxdiff(ab[1]), zip( self.parts, other.parts)))
0277
0278     def tr_maxdiff(self, other):
0279         """
0280         :param other: Prim instance to compare self with
0281         :return value: max absolute difference between the (numtran,4,4) elements
0282         """
0283         return np.max(np.abs(self.trans_ - other.trans_))
0284
0285     def __repr__(self):
0286         fmt = "sidx %2d primIdx %3s idx %30s lvIdx %3d lvName %30s partOffset %3s numParts %3s tranOffset %3s numTran %3s planOffset %3s  "
0287         return fmt % (self.sidx, self.primIdx, str(self.idx), self.lvIdx, self.lvName, self.partOffset, self.numParts, self.tranOffset, self.numTran, self.planOffset )
0288
0289     def __str__(self):
0290         return "\n".join(["",repr(self)] + list(map(str,list(filter(lambda pt:pt.tc > 0, self.parts)))) + [repr(self.trans_)])
0291
0292
0293 class Solid(object):
0294     """
0295     This is Solid in the Opticks composite sense, formerly poorly named Dir
0296
0297     Handles the concatenated prims from $TMP/GParts/{0,1,2,...}
0298
0299     TODO: rename to "Prims"
0300     """
0301     def __init__(self, base, kd):
0302         """
0303         :param base: directory containing primBuffer.npy etc..  eg $TMP/GParts/0
0304
0305
0306         ::
0307
0308             In [5]: p.primBuf
0309             Out[5]:
0310             array([[0, 1, 0, 3],
0311                    [1, 7, 1, 0]], dtype=uint32)
0312
0313             In [6]: p.partBuf.shape
0314             Out[6]: (8, 4, 4)
0315
0316             In [7]: p.partBuf
0317             Out[7]:
0318             array([[[    0.    ,     0.    ,     0.    ,  1000.    ],
0319                     [    0.    ,     0.    ,     0.    ,     0.    ],
0320                     [-1000.01  , -1000.01  , -1000.01  ,     0.    ],
0321                     [ 1000.01  ,  1000.01  ,  1000.01  ,     0.    ]],
0322
0323                    [[    0.    ,     0.    ,     0.    ,   500.    ],
0324                     [    0.    ,     0.    ,     0.    ,     0.    ],
0325                     [ -500.01  ,  -500.01  ,  -500.01  ,     0.    ],
0326                     [  500.01  ,   500.01  ,   500.01  ,     0.    ]],
0327
0328         """
0329         sidx = int(os.path.basename(base))
0330         self.sidx = sidx
0331         self.base = base
0332         self.blib = BLib(kd)
0333
0334         primbuf = np.load(os.path.join(base, "primBuffer.npy"))  # "solid" tree level index into part and tran buffers
0335         idxbuf  = np.load(os.path.join(base, "idxBuffer.npy"))
0336         partbuf = np.load(os.path.join(base, "partBuffer.npy"))
0337         tranbuf = np.load(os.path.join(base, "tranBuffer.npy"))
0338
0339         assert len(primbuf) == len(idxbuf)
0340         assert len(tranbuf) <= len(partbuf)
0341
0342         tot_prim = len(primbuf)
0343         tot_tran = len(tranbuf)
0344         tot_part = len(partbuf)   # aka node
0345
0346         #assert tot_part in [1,3,7,15,31]   nope these are totals for all the prim
0347         """
0348                             1                        ->  1
0349                  10                  11              ->  3     +2
0350            100       101        110       111        ->  7     +4
0351         1000 1001 1010 1011  1100 1101 1110  1111    -> 15     +8
0352                                                      -> 31     +16
0353         """
0354         self.primbuf = primbuf
0355         self.partbuf = partbuf
0356         self.tranbuf = tranbuf
0357         self.idxbuf = idxbuf
0358
0359         log.info(str(self))
0360
0361
0362         ma = Mesh(kd)   # GItemList/GMeshLib.txt LV name2idx idx2name
0363
0364         lvidxs = idxbuf[:,1]
0365         lvns = list(map(lambda lvidx:ma.idx2name[lvidx],lvidxs))
0366         self.lvns = lvns
0367
0368         ntran = np.zeros( tot_prim, dtype=np.uint32)
0369         ntran[0:tot_prim-1] = primbuf[1:,2] - primbuf[:-1,2]    ## differencing the tranOffsets to give numtran
0370
0371         ntran_tot_excluding_last = ntran.sum()
0372         ntran_last = tot_tran - ntran_tot_excluding_last
0373         ntran[tot_prim-1] = ntran_last
0374
0375         self.ntran = ntran
0376
0377         self.ma = ma
0378         self.prims = self.get_prims()
0379
0380     def get_prims(self):
0381         """
0382         :return pp: python array of Prim instances deserialized from self.prim array
0383         """
0384         pp = []
0385
0386         tot_prim = len(self.primbuf)
0387         log.debug("get_prims sidx %d tot_prim %d " % (self.sidx, tot_prim))
0388
0389         for primIdx in range(tot_prim):
0390             prim_item = self.primbuf[primIdx]
0391             p = Prim(primIdx, prim_item, self)
0392             pp.append(p)
0393         pass
0394         return pp
0395
0396     def enumerate_prim_zip(self, other):
0397         """
0398         :return i_ab: (idx,(this_prim,other_prim))  i_ab[0]
0399         """
0400         assert len(self.prims) == len(other.prims)
0401         return enumerate(zip(self.prims, other.prims))
0402
0403     def where_discrepant_tr(self, other, cut=0.1):
0404         assert len(self.prims) == len(other.prims)
0405         return map(lambda i_ab:i_ab[0], filter(lambda i_ab:i_ab[1][0].tr_maxdiff(i_ab[1][1]) > cut , enumerate(zip(self.prims, other.prims)) ))
0406
0407     def where_discrepant_prims(self, other, cut=0.1):
0408         assert len(self.prims) == len(other.prims)
0409         return map(lambda i_ab:i_ab[0], filter(lambda i_ab:i_ab[1][0].maxdiff(i_ab[1][1]) > cut , enumerate(zip(self.prims, other.prims)) ))
0410
0411     def get_discrepant_prims(self, other, cut=0.1):
0412         assert len(self.prims) == len(other.prims)
0413         return filter(lambda i_ab:i_ab[1][0].maxdiff(i_ab[1][1]) > cut , enumerate(zip(self.prims, other.prims)) )
0414
0415     def __repr__(self):
0416         fmt = "Solid %d : %s : primbuf %s partbuf %s tranbuf %s idxbuf %s "
0417         return fmt % (self.sidx, self.base, repr(self.primbuf.shape), repr(self.partbuf.shape), repr(self.tranbuf.shape), repr(self.idxbuf.shape))
0418
0419
0420 if __name__ == '__main__':
0421     logging.basicConfig(level=logging.INFO)
0422
0423     ridx = 0
0424     ddir = os.path.expandvars(os.path.join("$TMP/GParts",str(ridx)))
0425     dir_ = sys.argv[1] if len(sys.argv) > 1 else ddir
0426     sli_ = sys.argv[2] if len(sys.argv) > 2 else "0:10"
0427     sli = slice(*map(int, sli_.split(":")))
0428
0429     if dir_ == ddir:
0430         log.warning("using hardcoded dir" ) ;
0431     pass
0432
0433     from opticks.ana.key import keydir
0434     kd = keydir(os.environ["OPTICKS_KEY"])
0435
0436     d = Dir(dir_, kd)
0437     print("Dir(dir_)", d)
0438
0439     pp = d.prims
0440     print("dump sliced prims from the dir slice %s " % repr(sli))
0441     for p in pp[sli]:
0442         print(p)
0443     pass
0444     #print(d.tran)
0445