eic-opticks/ana/analytic_boundary_dbg.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 tboolean.py: CSG boolean geometry tracetest
0023 =======================================================
0024
0025 For analysis of photon_buffer written by oxrap/cu/generate.cu:tracetest
0026 which is intended for testing geometrical intersection only, with
0027 no propagation.
0028
0029 Unsigned boundaries 1 and 124 all from CSG intersection, the 215 was a
0030 abitrary marker written by evaluative_csg::
0031
0032     In [43]: flg[:100]
0033     Out[43]:
0034     A()sliced
0035     A([[  0,   0,   1, 215],
0036            [  0,   0, 123,   0],
0037            [  0,   0,   1, 215],
0038            [  0,   0, 123,   0],
0039            [  0,   0, 123,   0],
0040            [  0,   0,   1, 215],
0041            [  0,   0, 123,   0],
0042            [  0,   0,   1, 215],
0043            [  0,   0, 123,   0],
0044            [  0,   0, 123,   0],
0045            [  0,   0, 123,   0],
0046            [  0,   0,   1, 215],
0047            [  0,   0, 123,   0],
0048            [  0,   0, 124, 215],
0049            [  0,   0,   1, 215],
0050            [  0,   0,   1, 215],
0051            [  0,   0,   1, 215],
0052            [  0,   0,   1, 215],
0053            [  0,   0,   1, 215],
0054            [  0,   0,   1, 215],
0055
0056
0057     In [46]: np.unique(flg[flg[:,3] == 215][:,2])
0058     Out[46]:
0059     A()sliced
0060     A([  1, 124], dtype=uint32)
0061
0062
0063     In [48]: count_unique_sorted(flg[:,2])   # unsigned 0-based boundaries
0064     Out[48]:
0065     array([[    1, 58822],
0066            [  123, 35537],
0067            [  124,  5641]], dtype=uint64)
0068
0069
0070     In [57]: t0[np.where(flg[:,2] == 1)]
0071     A([ 149.8889,  149.8889,  349.8889, ...,  348.1284,  349.8889,  349.8889], dtype=float32)
0072
0073     In [58]: t0[np.where(flg[:,2] == 124)]
0074     A([ 51.6622,  47.394 ,  61.4086, ...,  61.0178,  66.7235,  47.1538], dtype=float32)
0075
0076     In [59]: t0[np.where(flg[:,2] == 123)]
0077     A([ 1400.    ,  1400.    ,  1400.    , ...,  2021.8895,  1400.    ,  1400.    ], dtype=float32)
0078
0079
0080 ::
0081
0082     In [13]: count_unique(ib)   # signed 1-based boundaries encode inner/outer
0083     Out[13]:
0084     array([[  -125.,   4696.],
0085            [    -2.,  58316.],
0086            [     2.,    506.],
0087            [   124.,  35537.],
0088            [   125.,    945.]])
0089
0090     In [35]: count_unique(ub)   # unsigned 0-based boundaries, for lookup against the blib,   58316 + 506 = 58822,  4696 + 945 = 5641
0091     Out[35]:
0092     array([[    1, 58822],
0093            [  123, 35537],
0094            [  124,  5641]], dtype=uint64)
0095
0096
0097
0098     In [33]: print "\n".join(["%3d : %s " % (i, n) for i,n in enumerate(blib.names)])
0099       0 : Vacuum///Vacuum
0100       1 : Vacuum///Rock
0101       2 : Rock///Air
0102       3 : Air/NearPoolCoverSurface//PPE
0103       4 : Air///Aluminium
0104       5 : Aluminium///Foam
0105       6 : Foam///Bakelite
0106     ...
0107     119 : OwsWater/NearOutInPiperSurface//PVC
0108     120 : OwsWater/NearOutOutPiperSurface//PVC
0109     121 : DeadWater/LegInDeadTubSurface//ADTableStainlessSteel
0110     122 : Rock///RadRock
0111
0112
0113 """
0114
0115 import os, sys, logging, numpy as np
0116 log = logging.getLogger(__name__)
0117
0118 import matplotlib.pyplot as plt
0119 #from opticks.ana.nbase import count_unique_sorted  # doenst work with signed
0120 from opticks.ana.nbase import count_unique
0121 from opticks.ana.evt import Evt
0122 from opticks.ana.proplib import PropLib
0123
0124 X,Y,Z,W = 0,1,2,3
0125
0126
0127 if __name__ == '__main__':
0128     from opticks.ana.main import opticks_main
0129     #args = opticks_main(tag="1", det="boolean", src="torch")
0130     ok = opticks_main()
0131
0132
0133     blib = PropLib("GBndLib")
0134
0135     evt = Evt(tag=ok.tag, det=ok.det, src=ok.src, pfx=ok.pfx, args=ok)
0136
0137     if not evt.valid:
0138        log.fatal("failed to load evt %s " % repr(args))
0139        sys.exit(1)
0140
0141     ox = evt.ox
0142
0143     ## assume layout written by oxrap/cu/generate.cu:tracetest ##
0144
0145     p0 = ox[:,0,:W]
0146     d0 = ox[:,1,:W]
0147     t0 = ox[:,1,W]
0148     flg = ox[:,3].view(np.uint32)  # u-flags
0149
0150     p1 = p0 + np.repeat(t0,3).reshape(-1,3)*d0   # intersect position
0151
0152
0153     ub = flg[:,2]                   # unsigned boundary
0154     ib = ox[:,2,W].view(np.int32)   # signed boundary
0155
0156     b = 1
0157     #b = 123
0158     #b = 124
0159
0160
0161     thin = slice(0,None,100)
0162     s = np.where(ub == b)[0]
0163
0164     plt.ion()
0165     plt.close()
0166
0167
0168     #plt.scatter( p0[s,X], p0[s,Y] )
0169     plt.scatter( p1[s,X], p1[s,Y] )
0170
0171
0172     plt.show()
0173