CSG/tests/CSGIntersectSolidTest.py

0001 #!/usr/bin/env python
0002 """
0003 CSGIntersectSolidTest.py
0004 ==========================
0005
0006
0007 AnnulusFourBoxUnion
0008     notice that spurious intersects all from 2nd circle roots
0009
0010
0011 IXYZ=-6,0,0 SPUR=1 ./csg_geochain.sh ana
0012
0013
0014
0015 """
0016 import os, logging, numpy as np
0017 from opticks.ana.fold import Fold
0018 from opticks.sysrap.SCenterExtentGenstep import SCenterExtentGenstep
0019
0020 from opticks.ana.gridspec import GridSpec, X, Y, Z
0021 from opticks.ana.npmeta import NPMeta
0022 from opticks.ana.eget import efloat_, efloatlist_, eint_, eintlist_
0023
0024
0025 SIZE = np.array([1280, 720])
0026
0027 log = logging.getLogger(__name__)
0028 np.set_printoptions(suppress=True, edgeitems=5, linewidth=200,precision=3)
0029
0030 try:
0031     import matplotlib.pyplot as mp
0032 except ImportError:
0033     mp = None
0034 pass
0035 #mp = None
0036
0037 try:
0038     import vtk
0039     import pyvista as pv
0040     from pyvista.plotting.colors import hexcolors
0041     themes = ["default", "dark", "paraview", "document" ]
0042     pv.set_plot_theme(themes[1])
0043 except ImportError:
0044     pv = None
0045     hexcolors = None
0046 pass
0047 #pv = None
0048
0049
0050
0051 def make_arrows(cent, direction, mag=1):
0052     direction = direction.copy()
0053     if cent.ndim != 2:
0054         cent = cent.reshape((-1, 3))
0055
0056     if direction.ndim != 2:
0057         direction = direction.reshape((-1, 3))
0058
0059     direction[:,0] *= mag
0060     direction[:,1] *= mag
0061     direction[:,2] *= mag
0062
0063     pdata = pv.vector_poly_data(cent, direction)
0064     # Create arrow object
0065     arrow = vtk.vtkArrowSource()
0066     arrow.Update()
0067     glyph3D = vtk.vtkGlyph3D()
0068     glyph3D.SetSourceData(arrow.GetOutput())
0069     glyph3D.SetInputData(pdata)
0070     glyph3D.SetVectorModeToUseVector()
0071     glyph3D.Update()
0072
0073     arrows = pv.wrap(glyph3D.GetOutput())
0074     return arrows
0075
0076
0077
0078 def lines_rectangle_YX(center, halfside):
0079     """
0080
0081           1                0
0082            +------+------+
0083            |      |      |
0084            |      |      |
0085            |- - - + - - -|
0086            |      |      |
0087            |      |      |
0088            +------+------+
0089          2                3
0090
0091
0092          X
0093          |
0094          +-- Y
0095
0096     """
0097     p0 = np.array( [center[0]+halfside[0], center[1]+halfside[1], center[2] ])
0098     p1 = np.array( [center[0]+halfside[0], center[1]-halfside[1], center[2] ])
0099     p2 = np.array( [center[0]-halfside[0], center[1]-halfside[1], center[2] ])
0100     p3 = np.array( [center[0]-halfside[0], center[1]+halfside[1], center[2] ])
0101
0102     ll = np.zeros( (4,2,3),  dtype=np.float32 )
0103
0104     ll[0] = p0, p1
0105     ll[1] = p1, p2
0106     ll[2] = p2, p3
0107     ll[3] = p3, p0
0108
0109     return ll.reshape(-1,3)
0110
0111
0112 def AnnulusFourBoxUnion_YX(pl, opt="+X -X +Y -Y circ", color='cyan', width=1, radius=45.):
0113     if opt.find("circ") > -1:
0114         arc = pv.CircularArc([0,  radius, 0], [0, radius, 0], [0, 0, 0], negative=True)
0115         pl.add_mesh(arc, color='cyan', line_width=1)
0116     pass
0117     lls = []
0118     if opt.find("+X") > -1:
0119         llpx = lines_rectangle_YX([ 52., 0.,0.],[10.,11.5, 6.5])
0120         lls.append(llpx)
0121     pass
0122     if opt.find("-X") > -1:
0123         llnx = lines_rectangle_YX([-52., 0.,0.],[10.,11.5, 6.5])
0124         lls.append(llnx)
0125     pass
0126     if opt.find("+Y") > -1:
0127         llpy = lines_rectangle_YX([0., 50.,0.],[15.,15.,  6.5])
0128         lls.append(llpy)
0129     pass
0130     if opt.find("-Y") > -1:
0131         llny = lines_rectangle_YX([0.,-50.,0.],[15.,15.,  6.5])
0132         lls.append(llny)
0133     pass
0134     for ll in lls:
0135         pl.add_lines(ll, color=color, width=width )
0136     pass
0137
0138 def ReferenceGeometry(pl):
0139     geom = os.environ.get("GEOM",None)
0140     if geom is None: return
0141
0142     print("ReferenceGeometry for GEOM %s " % geom)
0143     if geom == "AnnulusFourBoxUnion_YX":
0144         AnnulusFourBoxUnion_YX(pl, opt="+X -X +Y -Y circ")
0145     elif geom == "BoxCrossTwoBoxUnion_YX":
0146         AnnulusFourBoxUnion_YX(pl, opt="+Y +X")
0147
0148
0149     else:
0150         print("ReferenceGeometry not implemented for GEOM %s " % geom)
0151     pass
0152
0153
0154
0155 class Rays(object):
0156     """
0157
0158     ## different presentation of selected intersects
0159     ## red position points (usually invisible under arrow) and normal arrows
0160
0161     #pl.add_arrows( s_ray_origin, s_ray_direction, color="blue", mag=s_t )
0162     # drawing arrow from ray origins to intersects works, but the arrows are too big
0163     # TODO: find way to do this with different presentation
0164     # as it avoids having to python loop over the intersects as done
0165     # below which restricts to small numbers
0166     """
0167     @classmethod
0168     def Draw(cls, pl, p0, p1, n1, n0=None, ori_color="magenta", nrm_color="red", ray_color="blue", ray_width=1, pos_shift=None ):
0169         """
0170         :param pl: pyvista plotter
0171         :param p0: start position
0172         :param p1: end position
0173         :param n1: direction vec at end
0174         :param n0: direction vec at start or None
0175         :param pos_shift: shift vector or None
0176         """
0177         assert len(p0) == len(p1)
0178
0179         if not pos_shift is None:
0180             q0 = p0 + pos_shift
0181             q1 = p1 + pos_shift
0182         else:
0183             q0 = p0
0184             q1 = p1
0185         pass
0186
0187         pl.add_points( q0, color=nrm_color )
0188
0189         if not n1 is None:
0190             assert len(q0) == len(n1)
0191             pl.add_arrows( q1, n1, color=nrm_color, mag=10 )
0192         pass
0193
0194         if not n0 is None:
0195             assert len(q0) == len(n0)
0196             pl.add_arrows( q0, n0, color=nrm_color, mag=10 )
0197         pass
0198
0199         ll = np.zeros( (len(q0), 2, 3), dtype=np.float32 )
0200         ll[:,0] = q0
0201         ll[:,1] = q1
0202
0203         pl.add_points( q0, color=ori_color, point_size=16.0 )
0204         for i in range(len(ll)):
0205             pl.add_lines( ll[i].reshape(-1,3), color=ray_color, width=ray_width )
0206         pass
0207     pass
0208
0209
0210 class CSGRecord(object):
0211     @classmethod
0212     def Draw0(cls, pl, recs):
0213         """
0214         # works but then difficult for labelling
0215         """
0216         if len(recs) > 0:
0217             arrows = make_arrows( recs[:,2,:3], recs[:,0,:3], mag=10 )
0218             pl.add_mesh( arrows, color="pink" )
0219             print(arrows)
0220         pass
0221     @classmethod
0222     def Draw1(cls, pl, recs):
0223         """
0224         Old record layout
0225         """
0226         rec_tloop_nodeIdx_ctrl = rec[3,:3].view(np.int32)
0227         arrows = make_arrows( rec_pos, rec_dir, mag=10 )
0228         pl.add_mesh( arrows, color="pink" )
0229         points = arrows.points
0230         mask = slice(0,1)
0231         pl.add_point_labels( points[mask], [ str(rec_tloop_nodeIdx_ctrl) ] , point_size=20, font_size=36 )
0232
0233     @classmethod
0234     def Draw(cls, pl, recs, ray_origin, ray_direction, off):
0235         """
0236         :param pl: pyvista plotter
0237         :param recs: CSGRecord array for a single ray CSG tree intersection
0238         :param ray_origin:
0239         :param ray_direction:
0240         :param off: normalized vector out of the genstep plane, pointing to viewpoint
0241
0242         # less efficient but gives easy access to each arrow position
0243         """
0244
0245         pos_shifts = []
0246         for i in range(len(recs)):
0247             left  = recs[i,0]
0248             n_left = left[:3]
0249             t_left  = np.abs(left[3])
0250
0251             right = recs[i,1]
0252             n_right = right[:3]
0253             t_right = np.abs(right[3])
0254
0255             tmin = recs[i,3,0]
0256             t_min = recs[i,3,1]
0257             tminAdvanced = recs[i-1,3,2] if i > 0 else 0.
0258
0259             ## tminAdvanced is bit confusing as collected
0260             ## into record prior to the decision  (for leaf looping anyhow)
0261
0262             print("CSGRecord.Draw rec %d : tmin %10.4f t_min %10.4f tminAdvanced %10.4f (from prior rec) " % (i, tmin, t_min, tminAdvanced ))
0263
0264             result = recs[i,4]
0265             t_result = np.abs(result[3])
0266             n_result = result[:3]
0267
0268             pos_t_min = ray_origin + t_min*ray_direction                # base
0269             pos_tmin  = ray_origin + tmin*ray_direction
0270             pos_tminAdvanced = ray_origin + tminAdvanced*ray_direction
0271             pos_left  = ray_origin + t_left*ray_direction
0272             pos_right = ray_origin + t_right*ray_direction
0273             pos_result = ray_origin + t_result*ray_direction
0274
0275             pos_shift_0 = float(i*10+0)*off
0276             pos_shift_1 = float(i*10+1)*off
0277             pos_shift_2 = float(i*10+2)*off
0278             pos_shift_3 = float(i*10+3)*off
0279
0280             pos_shifts.append(pos_shift_0)
0281
0282             if tminAdvanced > 0.:
0283                 Rays.Draw(pl, pos_t_min, pos_tminAdvanced ,  ray_direction, ray_width=4, pos_shift=pos_shift_0, ray_color="green"  )
0284             pass
0285
0286             Rays.Draw(pl, pos_tmin,  pos_left ,  n_left,        ray_width=4, pos_shift=pos_shift_1  )
0287             Rays.Draw(pl, pos_tmin,  pos_right,  n_right,       ray_width=4, pos_shift=pos_shift_2  )
0288
0289             if t_result > 0.:
0290                 Rays.Draw(pl, pos_tmin,  pos_result, n_result,      ray_width=4, pos_shift=pos_shift_3, ray_color="cyan", nrm_color="cyan"  )
0291             pass
0292         pass
0293         return pos_shifts
0294     pass
0295 pass
0296
0297
0298
0299
0300 def apply_phicut(select):
0301     print(" SPHI : selecting intersects based on phi angle of the position ")
0302     phi0,phi1 = efloatlist_("SPHI", "0.25,1.75")
0303     cosPhi0,sinPhi0 = np.cos(phi0*np.pi),np.sin(phi0*np.pi)
0304     cosPhi1,sinPhi1 = np.cos(phi1*np.pi),np.sin(phi1*np.pi)
0305     PQ = cosPhi0*sinPhi1 - cosPhi1*sinPhi0
0306     PR = cosPhi0*pos[:,1] - pos[:,0]*sinPhi0
0307     QR = cosPhi1*pos[:,1] - pos[:,0]*sinPhi1   #Q ^ R = cosPhi1*d.y - d.x*sinPhi1
0308     select_phi = np.logical_and( PR > 0., QR < 0. ) if PQ > 0. else np.logical_or( PR > 0., QR < 0. )
0309     select_phi = np.logical_not( select_phi )
0310     count_phi = np.count_nonzero(select_phi)
0311     print( "%40s : %d " % ("count_phi", count_phi) )
0312     select = np.logical_and( select, select_phi )
0313     count_select = np.count_nonzero(select)
0314     print("%40s : %d   SPHI %s  phi0 %s phi1 %s PQ %s \n" % ("count_select", count_select, os.environ["SPHI"], phi0, phi1, PQ ))
0315     return select
0316 pass
0317
0318 def CSGIntersectFold():
0319     """
0320     eg /tmp/blyth/opticks/GeoChain_Darwin/nmskSolidMaskTail/CSGIntersectSolidTest/nmskSolidMaskTail_XZ
0321     """
0322     identdir = "$GEOM" if "GEOM" in os.environ else "$SOPR"
0323     isectFold = os.path.expandvars("$CFBASE/CSGIntersectSolidTest/%s" % identdir)
0324     print("CSGIntersectFolder : %s " % isectFold )
0325     return isectFold
0326
0327 def X4IntersectFold():
0328     """
0329     eg /tmp/blyth/opticks/extg4/X4IntersectSolidTest/nmskSolidMaskTail_XZ/X4Intersect/
0330     """
0331     identdir = "$GEOM" if "GEOM" in os.environ else "$SOPR"
0332     isectFold = os.path.expandvars("/tmp/$USER/opticks/extg4/X4IntersectSolidTest/%s/X4Intersect" % identdir)
0333     print("X4IntersectFolder : %s " % isectFold )
0334     return isectFold
0335
0336 def select_ixyz_genstep(select, ixyz):
0337     """
0338     when envvar IXYZ is defined restrict to a single genstep source of photons identified by grid coordinates
0339     """
0340     ix,iy,iz = ixyz
0341     print(isect_gsid)
0342     select_isect_gsid = np.logical_and( np.logical_and( isect_gsid[:,0] == ix , isect_gsid[:,1] == iy ), isect_gsid[:,2] == iz )
0343     count_isect_gsid = np.count_nonzero(select_isect_gsid)
0344
0345     select = np.logical_and( select, select_isect_gsid )
0346     count_select = np.count_nonzero(select)
0347
0348     print("%40s : %d   IXYZ %s  ix:%d iy:%d iz:%d" %   ("count_isect_gsid", count_isect_gsid, os.environ.get("IXYZ",None), ix,iy,iz ))
0349     print("%40s : %d  \n" % ("count_select", count_select))
0350
0351     return select
0352
0353 def select_iw_genstep(select, iw ):
0354     """
0355     when envvar IW is defined restrict to a single photon index, usually used together with IXYZ to pick one intersect
0356     by selecting a single intersect as well as single genstep
0357     """
0358     select_iw = isect_gsid[:,3] == iw
0359     count_iw = np.count_nonzero(select_iw)
0360     select = np.logical_and( select, select_iw )
0361     count_select = np.count_nonzero(select)
0362
0363     print("%40s : %d   IW %s " %   ("count_iw", count_iw, os.environ["IW"] ))
0364     print("%40s : %d  \n" % ("count_select", count_select))
0365
0366     return select
0367
0368
0369
0370
0371
0372
0373 if __name__ == '__main__':
0374     logging.basicConfig(level=logging.INFO)
0375
0376     pl = pv.Plotter(window_size=SIZE*2 )
0377
0378     if 'EDL' in os.environ:
0379         print("EDL : enable_eye_dome_lighting :  improves depth perception with point clouds")
0380         pl.enable_eye_dome_lighting()
0381     pass
0382
0383     if 'X4' in os.environ:
0384         x4_fold_ = X4IntersectFold()
0385         x4_cegs = SCenterExtentGenstep(x4_fold_)
0386         x4_fold = x4_cegs.fold
0387     else:
0388         x4_fold = None
0389     pass
0390
0391
0392     csg_fold_ = CSGIntersectFold()
0393     cegs = SCenterExtentGenstep(csg_fold_)
0394     fold = cegs.fold
0395     rel_name = "saveCenterExtentGenstepIntersect"
0396     recs_path = os.path.join(csg_fold_, rel_name )
0397
0398
0399     gspos = fold.gs[:,5,:3]
0400     gsid_ =  fold.gs[:,0,2].view(np.uint32).copy()
0401     gsid =  gsid_.view(np.int8).reshape(-1,4)            # q0.u.z (4 packed signed char)
0402
0403     no_gs = 'NO_GS' in os.environ
0404     if no_gs:
0405         print(" skip gspos points due to NO_GS " )
0406     else:
0407         pl.add_points( gspos, color="yellow" )
0408     pass
0409
0410     gsmeta = NPMeta(fold.gs_meta)  # TODO cxs_geochain uses genstep_meta : standardize
0411     gridspec = GridSpec(fold.peta, gsmeta)
0412     axes = gridspec.axes
0413     off = gridspec.off
0414     print(" gridspec.axes : %s   gridspec.off %s  " % (str(axes), str(off)) )
0415
0416     gridspec.pv_compose(pl)  ## composition is sensitive to envvars RESET PARA ZOOM
0417
0418     no_isect = not hasattr(fold, 'isect')
0419     gs_only = 'GS_ONLY' in os.environ
0420     quiet = True
0421
0422     if no_isect or gs_only:
0423        log.fatal("early exit just showing gensteps as no_isect:%s in fold OR gs_only:%d selected by GS_ONLY envvar" % (no_isect, gs_only))
0424        pl.show_grid()
0425        cp = pl.show()
0426        sys.exit(0)
0427     pass
0428
0429     dir_, t = fold.isect[:,0,:3], fold.isect[:,0,3]   # intersect normal and ray distance
0430     pos, sd = fold.isect[:,1,:3], fold.isect[:,1,3]   # intersect position and surface distance
0431     ray_origin    = fold.isect[:, 2, :3]
0432     ray_direction = fold.isect[:, 3, :3]
0433
0434     if not x4_fold is None:
0435         x4_pos = x4_fold.isect[:,0,:3]   # TODO: standardize isect layout
0436     else:
0437         x4_pos = None
0438     pass
0439
0440     isect_gsid_ = fold.isect[:,3,3].copy().view(np.uint32)
0441     isect_gsid = isect_gsid_.view(np.int8).reshape(-1,4)
0442     # genstep (ix,iy,iz,iw) from which each intersect originated from,
0443     # iw is photon index within genstep
0444
0445     ## bitwise_and to remove bytes holding the photon index, for easier connection to gs level
0446     isect_gsoid_ = np.bitwise_and( fold.isect[:,3,3].view(np.uint32), 0x00ffffff ).copy()
0447     isect_gsoid = isect_gsoid_.view(np.int8).reshape(-1,4)
0448     assert np.all( isect_gsid[:,:3]  == isect_gsoid[:,:3] )
0449
0450
0451     asd = np.abs(sd)   # absolute distance to surface : should be close to zero for all intersects
0452     asd_cut = efloat_("ASD", 1e-2)
0453     select_spurious = asd > asd_cut
0454     select_all = t > 0.
0455
0456
0457     gss_ = np.unique(isect_gsoid_[select_spurious])
0458     gss = gss_.view(np.int8).reshape(-1,4)         # gensteps which have spurious intersects
0459     gsid_idx = np.where(np.in1d(gsid_, gss_))[0]   # gsid_ indices of
0460     gss_select = np.isin(gsid_, gss_)              # genstep mask selecting gensteps the rays from which have spurious intersects
0461
0462     if no_gs:
0463         print(" skip gspos points due to NO_GS " )
0464     else:
0465         gss_select_count = np.count_nonzero(gss_select)
0466         if gss_select_count > 0:
0467             pl.add_points( gspos[gss_select], color="green" )
0468         pass
0469         ## gs with spurious intersects presented differently
0470     pass
0471
0472     select = select_all
0473
0474     count_all = np.count_nonzero(select_all)
0475     count_spurious = np.count_nonzero(select_spurious)
0476     cmdline = os.environ.get("CMDLINE", "no-CMDLINE")
0477
0478     print("\n\n")
0479     print(cmdline)
0480     print("\n")
0481     print( "%40s : %d " % ("count_all", count_all) )
0482     print( "%40s : %d " % ("count_spurious", count_spurious) )
0483
0484     if 'SPHI' in os.environ:
0485         select = apply_phicut(select)
0486     pass
0487
0488     if 'SPUR' in os.environ:  # when envvar SPUR is defined restrict to intersects that have sd less than sd_cut
0489         select = np.logical_and( select, select_spurious )
0490         count_select = np.count_nonzero(select)
0491         print("%40s : %d    SPUR %s " % ("count_spurious", count_spurious, os.environ["SPUR"] ))
0492         print("%40s : %d  \n" % ("count_select", count_select))
0493     pass
0494
0495     sxyzw = eintlist_("SXYZW", None)
0496     _sxyzw = np.array( sxyzw, dtype=np.int8 ).view(np.uint32)[0] if not sxyzw is None else None
0497     _sxyzw_idx = np.where( isect_gsid_ == _sxyzw )  if not _sxyzw is None else None
0498
0499     ixyz = eintlist_("IXYZ", None)
0500     if not ixyz is None:
0501         log.info("select_ixyz_genstep")
0502         select = select_ixyz_genstep(select, ixyz)
0503     pass
0504
0505     iw = eint_("IW","-1")
0506     if iw > -1:
0507         log.info("select_iw_genstep")
0508         select = select_iw_genstep(select, iw)
0509     pass
0510
0511     gsid_label = None
0512     if not ixyz is None and not iw is None:
0513         ix,iy,iz = ixyz
0514         gsid_label = "gsid_%d_%d_%d_%d" % (ix,iy,iz,iw)
0515         print(" IXYZ/IW defined =>  gsid_label : %s " % gsid_label )
0516     elif not sxyzw is None:
0517         sx,sy,sz,sw = sxyzw
0518         gsid_label = "gsid_%d_%d_%d_%d" % (sx,sy,sz,sw)
0519         print(" SXYZW defined  =>  gsid_label : %s " % gsid_label )
0520     pass
0521
0522     if not gsid_label is None:
0523         recs_fold = Fold.Load(recs_path, gsid_label, quiet=True)
0524         recs = [] if recs_fold is None else recs_fold.CSGRecord
0525         if len(recs) > 0:
0526             print("%40s : %s   recs_fold.base %s " % ("recs", str(recs.shape), recs_fold.base ))
0527         else:
0528             if quiet == False:
0529                 print(" no recs to load at recs_path/gsid_label %s/%s " % (recs_path, gsid_label))
0530             pass
0531         pass
0532     else:
0533         recs = []
0534     pass
0535
0536     pos_shifts = []
0537     if len(recs) > 0:
0538         assert not _sxyzw_idx is None
0539         _ray_origin = ray_origin[_sxyzw_idx][0]
0540         _ray_direction = ray_direction[_sxyzw_idx][0]
0541         pos_shifts = CSGRecord.Draw(pl, recs, _ray_origin, _ray_direction, off )
0542     pass
0543
0544     count_select = np.count_nonzero(select)
0545     print("%40s : %d  \n" % ("count_select", count_select))
0546
0547     s_t = t[select]
0548     s_pos = pos[select]
0549     s_sd = sd[select]
0550     s_isect_gsid = isect_gsid[select]
0551     s_dir = dir_[select]
0552     s_isect = fold.isect[select]
0553     s_ray_origin = ray_origin[select]
0554     s_ray_direction = ray_direction[select]
0555     s_pos_r = np.sqrt(np.sum(s_pos*s_pos, axis=1))
0556
0557     s_count = len(s_pos)
0558     s_limited = min( s_count, 100 )
0559     # default number of photons per genstep is 100, so this gets all of those when selecting single genstep
0560     selected_isect = s_isect[:s_limited]
0561
0562     print("%40s : %d  \n" % ("s_count", s_count))
0563     print("%40s : %d  \n" % ("s_limited", s_limited))
0564     print("%40s : %d  \n" % ("selected_isect", len(selected_isect)))
0565
0566     ## sub selection of the selected that are also spurious as judged by signed distance
0567
0568     select_and_spurious = np.logical_and( select, select_spurious )
0569     ss_pos = pos[select_and_spurious]
0570     ss_dir = dir_[select_and_spurious]
0571     ss_ray_origin = ray_origin[select_and_spurious]
0572
0573     ss_count = len(ss_pos)
0574     ss_limited = min( ss_count, 100 )
0575
0576     print("%40s : %d  \n" % ("ss_count", ss_count))
0577     print("%40s : %d  \n" % ("ss_limited", ss_limited))
0578
0579     def fmt(i):
0580         _s_isect_gsid = "SXYZW=%d,%d,%d,%d  " % tuple(s_isect_gsid[i])
0581         _s_t = " s_t ( %10.4f ) " % s_t[i]
0582         _s_pos = " s_pos ( %10.4f %10.4f %10.4f ) " % tuple(s_pos[i])
0583         _s_pos_r = " s_pos_r ( %10.4f ) " % s_pos_r[i]
0584         _s_sd = " s_sd ( %10.4f ) " % s_sd[i]
0585         return " %-20s %s %s %s %s " % (_s_isect_gsid, _s_t, _s_pos, _s_pos_r, _s_sd )
0586     pass
0587     desc = "\n".join(map(fmt,np.arange(s_limited)))
0588     print(desc)
0589     print(" Use IXYZ to select gensteps, IW to select photons within the genstep ")
0590
0591     log.info( "asd_cut %10.4g sd.min %10.4g sd.max %10.4g num select %d " % (asd_cut, sd.min(), sd.max(), len(s_pos)))
0592
0593
0594     selected_isect_path = "/tmp/selected_isect.npy"
0595     key = 'SAVE_SELECTED_ISECT'
0596     save_selected_isect = key in os.environ
0597     if save_selected_isect and len(selected_isect) > 0 :
0598         print(" %s : saving selected_isect to %s " % (key, selected_isect_path)  )
0599         np.save( selected_isect_path, selected_isect )
0600     else:
0601         print(" define key %s to save selected isect to %s " % (key, selected_isect_path ))
0602     pass
0603
0604     if 'PLOT_SELECTED_ONLY' in os.environ:   # normally all positions are plotted not just selected
0605         pl.add_points( s_pos, color="white" )
0606     else:
0607         if len(pos_shifts) == 0:
0608             pl.add_points( pos, color="white" )
0609         else:
0610             for pos_shift in pos_shifts:
0611                 pl.add_points( pos+pos_shift, color="white" )
0612             pass
0613         pass
0614     pass
0615
0616     if not x4_pos is None:
0617         print("x4_pos %s" % str(x4_pos.shape))
0618         pl.add_points( x4_pos+10*off, color="cyan" )
0619     else:
0620         print("x4_pos None")
0621     pass
0622
0623
0624     ## different presentation of selected intersects and normals
0625     if s_count > 0 and sxyzw is None:
0626         Rays.Draw( pl, s_ray_origin[:s_limited], s_pos[:s_limited], s_dir[:s_limited] )
0627     pass
0628     if ss_count > 0 and sxyzw is None:
0629         Rays.Draw( pl, ss_ray_origin[:ss_limited], ss_pos[:ss_limited], ss_dir[:ss_limited], ray_color="yellow", nrm_color="yellow" )
0630     pass
0631
0632     pl.show_grid()
0633     ReferenceGeometry(pl)
0634     topline = os.environ.get("TOPLINE", "topline")
0635     botline = os.environ.get("BOTLINE", "botline")
0636     pl.add_text(topline, position="upper_left")
0637     pl.add_text(botline, position="lower_left")
0638
0639     figsdir = os.path.join(csg_fold_, "figs")
0640     if not os.path.isdir(figsdir):
0641         os.makedirs(figsdir)
0642     pass
0643
0644     outname = "out.png"
0645     outpath = os.path.join(figsdir, outname)
0646     print(" outpath : %s " % outpath )
0647
0648     cp = pl.show(screenshot=outpath)
0649
0650