eic-opticks/ana/gridspec.py

0001 #!/usr/bin/env python
0002
0003 import logging, os
0004 import numpy as np
0005 log = logging.getLogger(__name__)
0006
0007 eary_ = lambda ekey, edef:np.array( list(map(float, os.environ.get(ekey,edef).split(","))) )
0008 #efloat_ = lambda ekey, edef: float( os.environ.get(ekey,edef) )
0009
0010 from opticks.ana.eget import efloat_
0011
0012 X,Y,Z = 0,1,2
0013
0014 from opticks.ana.axes import *
0015
0016
0017 class GridSpec(object):
0018     @classmethod
0019     def DemoPeta(cls, ix0=-2, ix1=2, iy0=-2, iy1=2, iz0=0, iz1=0, photons_per_genstep=100):
0020         peta = np.zeros( (1,4,4), dtype=np.float32 )
0021         peta.view(np.int32)[0,0] = [ix0,ix1,iy0,iy1]
0022         peta.view(np.int32)[0,1] = [iz0,iz1,photons_per_genstep, 0]
0023         ce = [ -492.566,  -797.087, 19285.   ,   264.   ]
0024         peta[0,2] = ce
0025         return peta
0026
0027
0028     def __init__(self, peta, gsmeta ):
0029         """
0030         :param peta:
0031         :param gsmeta:
0032
0033
0034
0035         """
0036         moi = gsmeta.find("moi:", None)
0037         midx = gsmeta.find("midx:", None)
0038         mord = gsmeta.find("mord:", None)
0039         iidx = gsmeta.find("iidx:", None)
0040
0041         coords = "RTP" if not iidx is None and int(iidx) == -3 else "XYZ"   ## NB RTP IS CORRECT ORDERING radiusUnitVec:thetaUnitVec:phiUnitVec
0042         log.info(" moi %s midx %s mord %s iidx %s coords %s " % (moi, midx, mord, iidx, coords))
0043
0044
0045         ## CSGOptiX::setCEGS tests/CSGOptiXSimtraceTest.cc
0046
0047         ix0,ix1,iy0,iy1 = peta[0,0].view(np.int32)
0048         iz0,iz1,photons_per_genstep,_ = peta[0,1].view(np.int32)
0049         gridscale = peta[0,1,3]
0050
0051         ce = tuple(peta[0,2])
0052         sce = (" %7.2f" * 4 ) % ce
0053
0054         assert photons_per_genstep != 0
0055         nx = (ix1 - ix0)//2
0056         ny = (iy1 - iy0)//2
0057         nz = (iz1 - iz0)//2
0058
0059         log.info(" ix0 %d ix1 %d nx %d  " % (ix0, ix1, nx))
0060         log.info(" iy0 %d iy1 %d ny %d  " % (iy0, iy1, ny))
0061         log.info(" iz0 %d iz1 %d nz %d  " % (iz0, iz1, nz))
0062         log.info(" gridscale %10.4f " % gridscale )
0063         log.info(" sce %s " % sce )
0064
0065         # below default from envvars are overridden for planar data
0066         eye = eary_("EYE","1.,1.,1.")
0067
0068         #look = eary_("LOOK","0.,0.,0.")
0069         look = ce[:3]
0070
0071         up  = eary_("UP","0.,0.,1.")
0072         off  = eary_("OFF","0.,0.,1.")
0073
0074
0075         EYES = efloat_("EYES", "6.")   # TODO: why 6 ? how to control FOV to gain more control of this
0076
0077
0078         axes = self.determine_axes(nx, ny, nz)
0079         planar = len(axes) == 2
0080
0081         if planar:
0082             H, V = axes
0083             axlabels =  coords[H], coords[V]
0084             HV = "%s%s" % (coords[H],coords[V])
0085             up  = Axes.Up(H,V)
0086             off = Axes.Off(H,V)
0087             eye = look + ce[3]*off*EYES
0088
0089         else:
0090             H, V, D = axes
0091             HV = None
0092             axlabels =  coords[H], coords[V], coords[D]
0093             up = XYZ.up
0094             off = XYZ.off
0095             ## hmm in 3D case makes less sense : better to just use the input EYE
0096
0097             eye = ce[3]*eye*EYES
0098
0099             pass
0100         pass
0101         log.info(" planar %d  eye %s  EYES %s " % (planar, str(eye), EYES))
0102
0103
0104
0105         self.coords = coords
0106         self.eye = eye
0107         self.look = look
0108         self.up  = up
0109         self.off = off
0110         self.HV = HV
0111         self.peta = peta
0112         self.gsmeta = gsmeta
0113
0114         self.axes = axes
0115         self.planar = planar
0116         self.axlabels = axlabels
0117
0118         self.nx = nx
0119         self.ny = ny
0120         self.nz = nz
0121         self.ce = ce
0122         self.sce = sce
0123         self.thirdline = " ce: " + sce
0124         self.photons_per_genstep = photons_per_genstep
0125
0126
0127     def pv_compose(self, pl ):
0128         """
0129         :param pl: pyvista plotter instance
0130         :param reset: for reset=True to succeed to auto-set the view, must do this after add_points etc..
0131
0132         Note for greater control of the view it is better to use reset=False
0133         """
0134
0135         PARA = "PARA" in os.environ
0136         RESET = "RESET" in os.environ
0137         ZOOM = efloat_("ZOOM", "1.")
0138
0139         #eye = look + self.off
0140
0141         look = self.look
0142         eye = self.eye
0143         up = self.up
0144
0145         print("pv_arrange_viewpoint look:%s eye: %s up:%s  PARA:%s RESET:%d ZOOM:%s  " % (str(look), str(eye), str(up), RESET, PARA, ZOOM ))
0146
0147         if PARA:
0148             pl.camera.ParallelProjectionOn()
0149         pass
0150
0151         pl.set_focus(    look )
0152         pl.set_viewup(   up )
0153         pl.set_position( eye, reset=RESET )   ## for reset=True to succeed to auto-set the view, must do this after add_points etc..
0154         pl.camera.Zoom(ZOOM)
0155
0156
0157
0158
0159     def __str__(self):
0160         return "GridSpec (nx ny nz) (%d %d %d) axes %s axlabels %s " % (self.nx, self.ny, self.nz, str(self.axes), str(self.axlabels) )
0161
0162     def determine_axes(self, nx, ny, nz):
0163         """
0164         :param nx:
0165         :param nx:
0166         :param nx:
0167
0168         With planar axes the order is arranged to make the longer axis the first horizontal one
0169         followed by the shorter axis as the vertical one.
0170
0171             +------+
0172             |      |  nz     ->  ( Y, Z )    ny_over_nz > 1
0173             +------+
0174                ny
0175
0176             +------+
0177             |      |  ny     ->  ( Z, Y )    ny_over_nz < 1
0178             +------+
0179                nz
0180
0181         """
0182         if nx == 0 and ny > 0 and nz > 0:
0183             ny_over_nz = float(ny)/float(nz)
0184             axes = (Y,Z) if ny_over_nz > 1 else (Z,Y)
0185         elif nx > 0 and ny == 0 and nz > 0:
0186             nx_over_nz = float(nx)/float(nz)
0187             axes = (X,Z) if nx_over_nz > 1 else (Z,X)
0188         elif nx > 0 and ny > 0 and nz == 0:
0189             nx_over_ny = float(nx)/float(ny)
0190             axes = (X,Y) if nx_over_ny > 1 else (Y,X)
0191         else:
0192             axes = (X,Y,Z)
0193         pass
0194         return axes
0195
0196
0197 class CrossHairs(object):
0198     @classmethod
0199     def draw(cls, pl, sz=100):
0200
0201         nll = 3 ;
0202         ll = np.zeros( (nll, 2, 3), dtype=np.float32 )
0203         ll[0,0] = (-sz,0,0)
0204         ll[0,1] = (sz,0,0)
0205
0206         ll[1,0] = (0,-sz,0)
0207         ll[1,1] = (0,sz,0)
0208
0209         ll[2,0] = (0,0,-sz)
0210         ll[2,1] = (0,0,sz)
0211
0212         pl.add_points( ll.reshape(-1,3), color="magenta", point_size=16.0 )
0213
0214         for i in range(len(ll)):
0215              pl.add_lines( ll[i].reshape(-1,3), color="blue" )
0216         pass
0217
0218
0219
0220
0221 if __name__ == '__main__':
0222      logging.basicConfig(level=logging.INFO)
0223      peta = GridSpec.DemoPeta()
0224      grid = GridSpec(peta)
0225      print(grid)
0226
0227