eic-opticks/ana/framegensteps.py

0001 #!/usr/bin/env python
0002 """
0003 framegensteps.py
0004 =================
0005
0006 Uses per-genstep transforms to give world_frame_centers and
0007 then uses the frame.w2m to get the grid centers in the target frame.
0008
0009 Related:
0010
0011 * sysrap/sframe.h
0012 * sysrap/SFrameGenstep.cc:SFrameGenstep::MakeCenterExtentGensteps
0013
0014
0015 """
0016 import os, logging, numpy as np
0017 from opticks.ana.axes import X,Y,Z
0018
0019 log = logging.getLogger(__name__)
0020
0021 class FrameGensteps(object):
0022     """
0023     Transform enabled gensteps:
0024
0025     * gs[igs,0,3] photons to generate for genstep *igs*
0026     * gs[igs,1] local frame center position
0027     * gs[igs,2:] 4x4 transform
0028
0029     From SEvent::ConfigureGenstep::
0030
0031     * gsid was MOVED from (1,3) to (0,2) when changing genstep to carry transform
0032
0033     Notice that every genstep has its own transform with slightly
0034     different translations according to the different grid points
0035     Conversely there is only one overall frame transform
0036     which corresponds to the targetted piece of geometry.
0037     """
0038     def __init__(self, genstep, frame, local=True, symbol="gs"):
0039         """
0040         :param genstep: (num_gs,6,4) array with grid transforms in 2: and position in 1
0041         :param frame: sframe instance (replacing former metatran array of 3 transforms and grid GridSpec instance)
0042         :param local: bool
0043         :param local_extent_scale: SUSPECT THIS SHOULD NO LONGER EVER BE TRUE
0044
0045         The way *local* is used implies that the simtrace genstep in q1 contains global centers,
0046         that are transformed here by frame.w2m to give the genstep centers *ugsc* in the local frame
0047
0048         BUT t.genstep[:,1] they are all origins [0., 0., 0., 1.]
0049
0050         YES, but the genstep transfrom is applied to that origin to give the world_frame_centers
0051         hence need the frame.w2m transform to give the local frame grid centers.
0052
0053         """
0054         gs = genstep
0055         local_extent_scale = frame.coords == "RTP"  ## KINDA KLUDGE DUE TO EXTENT HANDLING BEING DONE BY THE RTP TRANSFORM
0056
0057         numpho = gs.view(np.int32)[:,0,3]        # q0.i.w  top right values from all gensteps
0058         gsid = gs.view(np.int32)[:,0,2].copy()   # q0.i.z  SEvent::ConfigureGenstep
0059         all_one = np.all( gs[:,1,3] == 1. )      # q1.f.w
0060         assert all_one   # from SEvent::MakeCenterExtentGensteps that q1.f.w should always to 1.f
0061
0062
0063         ## apply the 4x4 transform in rows 2: to the position in row 1
0064         world_frame_centers = np.zeros( (len(gs), 4 ), dtype=np.float32 )
0065         for igs in range(len(gs)):
0066             gs_pos = gs[igs,1]          ## normally origin (0,0,0,1)
0067             gs_tran = gs[igs,2:]        ## m2w with grid translation
0068             gs_tran[:,3] = [0,0,0,1]   ## fixup 4th column, as may contain identity info
0069             world_frame_centers[igs] = np.dot( gs_pos, gs_tran )
0070             #   world_frame_centers = m2w * grid_translation * model_frame_positon
0071         pass
0072
0073         centers_local = np.dot( world_frame_centers, frame.w2m )
0074         # use w2m to transform global frame centers back to local frame
0075
0076         if local and local_extent_scale:
0077             extent = frame.ce[3]
0078             centers_local[:,:3] *= extent
0079             assert 0
0080             ## HUH:confusing, surely the horse has left the stable ?
0081             ## what was done on device is what matters, so why scale here ?
0082             ## vague recollection this is due to some issue with RTP tangential transforms
0083             ## TODO: ELIMINATE once get RTP tangential operational again
0084         pass
0085
0086
0087         ugsc = centers_local if local else world_frame_centers
0088
0089         lim = {}
0090         lim[X] = np.array([ugsc[:,X].min(), ugsc[:,X].max()])
0091         lim[Y] = np.array([ugsc[:,Y].min(), ugsc[:,Y].max()])
0092         lim[Z] = np.array([ugsc[:,Z].min(), ugsc[:,Z].max()])
0093
0094         self.gs = gs
0095         self.gsid = gsid
0096         self.numpho = numpho
0097         self.totpho = np.sum(numpho)
0098         self.centers = world_frame_centers
0099         self.centers_local = centers_local
0100         self.ugsc = ugsc
0101         self.lim = lim
0102         self.symbol = symbol
0103         log.info(repr(self))
0104
0105     def __repr__(self):
0106         symbol = self.symbol
0107         return "\n".join([
0108                    "FrameGensteps",
0109                    "%s.gs %s " % (symbol, str(self.gs.shape)),
0110                    "%s.centers %s " % (symbol, str(self.centers.shape)),
0111                    "%s.centers_local %s " % (symbol, str(self.centers_local.shape)),
0112                    "%s.numpho[0] %d " % (symbol, self.numpho[0]) ,
0113                    "%s.totpho    %d " % (symbol, self.totpho) ,
0114                    "%s.lim[X] %s " % (symbol,str(self.lim[X])),
0115                    "%s.lim[Y] %s " % (symbol,str(self.lim[Y])),
0116                    "%s.lim[Z] %s " % (symbol,str(self.lim[Z])),
0117               ])
0118
0119     @classmethod
0120     def CombineLim(cls, stuv_ ):
0121         """
0122         """
0123         stuv = list(filter(None, stuv_))
0124         lim = {}
0125         if len(stuv) == 1:
0126             s,t,u = stuv[0], None,None
0127             lim = s.lim
0128         elif len(stuv) == 2:
0129             s,t,u = stuv[0], stuv[1],None
0130             for d in [X,Y,Z]:
0131                 sl = s.lim[d]
0132                 tl = t.lim[d]
0133                 assert tl.shape == sl.shape == (2,)
0134                 lim[d] = np.array( [min(tl[0],sl[0]), max(tl[1],sl[1])] )
0135             pass
0136         elif len(stuv) == 3:
0137             s,t,u = stuv[0], stuv[1], stuv[2]
0138             for d in [X,Y,Z]:
0139                 sl = s.lim[d]
0140                 tl = t.lim[d]
0141                 ul = u.lim[d]
0142                 assert sl.shape == tl.shape == ul.shape == (2,)
0143                 lim[d] = np.array( [min(sl[0],tl[0],ul[0]), max(sl[1],tl[1],ul[0])] )
0144             pass
0145         elif len(stuv) == 4:
0146             s,t,u,v = stuv[0], stuv[1], stuv[2], stuv[3]
0147             for d in [X,Y,Z]:
0148                 sl = s.lim[d]
0149                 tl = t.lim[d]
0150                 ul = u.lim[d]
0151                 vl = v.lim[d]
0152                 assert sl.shape == tl.shape == ul.shape == vl.shape == (2,)
0153                 lim[d] = np.array( [min(sl[0],tl[0],ul[0],vl[0]), max(sl[1],tl[1],ul[1],vl[1])] )
0154             pass
0155         else:
0156             lim = None
0157         pass
0158         return lim
0159
0160