eic-opticks/ana/simtrace_plot.py

0001 #!/usr/bin/env python
0002
0003 import os, logging, numpy as np
0004 log = logging.getLogger(__name__)
0005
0006 from opticks.ana.eget import efloatlist_, elookce_, elook_epsilon_, eint_
0007 from opticks.sysrap.sframe import sframe , X, Y, Z
0008
0009 SIZE = np.array([1280, 720])   ## SIZE*2 [2560, 1440]
0010 LEGEND =  not "NOLEGEND" in os.environ # when not MASK=pos legend often too many lines, so can switch it off
0011 GSPLOT = eint_("GSPLOT", "0")
0012
0013
0014
0015 MP =  not "NOMP" in os.environ
0016 if MP:
0017     try:
0018         import matplotlib.pyplot as mp
0019         from matplotlib.patches import Circle, Rectangle, Ellipse
0020     except ImportError:
0021         mp = None
0022     pass
0023 else:
0024     mp = None
0025 pass
0026
0027
0028 PV =  not "NOPV" in os.environ
0029 PVGRID = not "NOPVGRID" in os.environ
0030 if PV:
0031     try:
0032         import pyvista as pv
0033         themes = ["default", "dark", "paraview", "document" ]
0034         pv.set_plot_theme(themes[1])
0035     except ImportError:
0036         pv = None
0037     pass
0038 else:
0039     pv = None
0040 pass
0041
0042 from opticks.ana.pvplt import *
0043
0044
0045
0046 class SimtracePlot(object):
0047     def __init__(self, pl, feat, gs, frame, pos, outdir ):
0048         """
0049         :param pl: pyvista plotter instance, can be None
0050         :param feat: Feature instance
0051         :param gs: FrameGensteps instance
0052         :param frame: sframe instance
0053         :param pos: Positions instance
0054         :param outdir: str
0055
0056         ## hmm regarding annotation, what should come from remote and what local ?
0057
0058         XX,YY,ZZ
0059            lists of ordinates for drawing lines parallel to axes
0060
0061         """
0062         if not os.path.isdir(outdir):
0063             os.makedirs(outdir)
0064         pass
0065         self.pl = pl
0066         self.feat = feat
0067         self.gs = gs
0068         self.frame = frame
0069         self.pos = pos
0070         self.outdir = outdir
0071         self.pl = None
0072
0073         topline = os.environ.get("TOPLINE", "CSGOptiXSimtraceTest.py:PH")
0074         botline = os.environ.get("BOTLINE", "cxs")
0075         note = os.environ.get("NOTE", "")
0076         note1 = os.environ.get("NOTE1", "")
0077
0078         self.topline = topline
0079         self.botline = botline
0080         self.note = note
0081         self.note1 = note1
0082
0083         self.look = efloatlist_("LOOK", "0,0,0")
0084         self.look_ce = elookce_(extent=10)
0085
0086
0087         epsilon = self.frame.propagate_epsilon
0088         if epsilon == 0.: epsilon = 0.05
0089         self.look_epsilon = elook_epsilon_(epsilon)
0090
0091
0092         aa = {}
0093         aa[X] = efloatlist_("XX")
0094         aa[Y] = efloatlist_("YY")
0095         aa[Z] = efloatlist_("ZZ")
0096
0097         self.aa = aa
0098         self.sz = float(os.environ.get("SZ","1.0"))
0099
0100         log.info(" aa[X] %s " % str(self.aa[X]))
0101         log.info(" aa[Y] %s " % str(self.aa[Y]))
0102         log.info(" aa[Z] %s " % str(self.aa[Z]))
0103
0104     def outpath_(self, stem="positions", ptype="pvplt"):
0105         sisel = self.feat.sisel
0106         return os.path.join(self.outdir,"%s_%s_%s.png" % (stem, ptype, self.feat.name))
0107
0108     def positions_mpplt(self):
0109         axes = self.frame.axes
0110         if len(axes) == 2:
0111             self.positions_mpplt_2D(legend=LEGEND, gsplot=GSPLOT)
0112         else:
0113             log.info("mp skip 3D plotting as PV is so much better at that")
0114         pass
0115
0116     def positions_mpplt_2D(self, legend=True, gsplot=0):
0117         """
0118         (H,V) are the plotting axes
0119         (X,Y,Z) = (0,1,2) correspond to absolute axes which can be mapped to plotting axes in various ways
0120
0121         when Z is vertical lines of constant Z appear horizontal
0122         when Z is horizontal lines of constant Z appear vertical
0123         """
0124         upos = self.pos.upos  # may have mask applied
0125         ugsc = self.gs.ugsc
0126         lim = self.gs.lim
0127
0128         H,V = self.frame.axes       # traditionally H,V = X,Z  but now generalized
0129         _H,_V = self.frame.axlabels
0130         log.info(" frame.axes H:%s V:%s " % (_H, _V))
0131
0132         feat = self.feat
0133         sz = self.sz
0134         print("positions_mpplt feat.name %s " % feat.name )
0135
0136         igs = slice(None) if len(ugsc) > 1 else 0
0137
0138         #fig, ax = mp.subplots(figsize=SIZE/100.)  # mpl uses dpi 100
0139         #title = [self.topline, self.botline, self.frame.thirdline]
0140         #fig.suptitle("\n".join(title))    # now done in frame.mp_subplots
0141
0142         fig, ax = self.frame.mp_subplots(mp)   # this way is sensitive to FOCUS
0143
0144         self.ax = ax  # TODO: arrange to pass this in
0145         self.fig = fig
0146
0147         note = self.note
0148         note1 = self.note1
0149
0150         if len(note) > 0:
0151              mp.text(0.01, 0.99, note, horizontalalignment='left', verticalalignment='top', transform=ax.transAxes)
0152         pass
0153         if len(note1) > 0:
0154              mp.text(0.01, 0.95, note1, horizontalalignment='left', verticalalignment='top', transform=ax.transAxes)
0155         pass
0156
0157
0158         # loop over unique values of the feature
0159         for idesc in range(feat.unum):
0160             uval, selector, label, color, skip, msg = feat(idesc)
0161             if skip: continue
0162             pos = upos[selector]
0163             ## hmm any masking needs be applied to both upos and selector ?
0164             ## alternatively could apply the mask early and use that from the Feature machinery
0165             #log.info(" feat.unum %d idesc %d label %s color %s msg %s ." % (feat.unum, idesc, label, color, msg ))
0166             log.debug(" feat.unum %d msg %s " % (feat.unum, msg ))
0167             ax.scatter( pos[:,H], pos[:,V], label=label, color=color, s=sz )
0168         pass
0169
0170
0171         mpplt_parallel_lines(ax, self.gs.lim, self.aa, self.frame.axes, linestyle="dashed" )
0172
0173         if hasattr(self, 'x_lpos'):
0174             x_lpos = self.x_lpos
0175             ax.scatter( x_lpos[:,H], x_lpos[:,V], label="x_lpos", s=10 )
0176             mpplt_add_contiguous_line_segments(ax, x_lpos, axes=self.frame.axes, linewidths=2)
0177         pass
0178         if hasattr(self, 't_spos'):
0179             t_spos = self.t_spos
0180             ax.scatter( t_spos[:,H], t_spos[:,V], label="t_spos", s=5 )
0181         pass
0182         if hasattr(self, 'simtrace_selection'):
0183             sts = self.simtrace_selection
0184             mpplt_simtrace_selection_line(ax, sts, axes=self.frame.axes, linewidths=2)
0185         pass
0186         if not self.look_ce is None:
0187             mpplt_ce_multiple(ax, self.look_ce, axes=self.frame.axes)
0188         pass
0189         if not self.look_epsilon is None:
0190             mpplt_ce(ax, self.look_epsilon, axes=self.frame.axes, colors="yellow" )
0191         pass
0192         label = "gs_center XZ"
0193         if gsplot > 0:
0194             ax.scatter( ugsc[igs, H], ugsc[igs,V], label=None, s=sz )
0195         pass
0196
0197         # not setting xlim ylim here, as thats
0198         # now done in sframe.py:mp_subplots
0199         #
0200         # ax.set_aspect('equal')
0201         # ax.set_xlim( lim[H] )
0202         # ax.set_ylim( lim[V] )
0203         # ax.set_xlabel(_H)
0204         # ax.set_ylabel(_V)
0205
0206         if legend:
0207             ax.legend(loc="upper right", markerscale=4)
0208             ptype = "mpplt"
0209         else:
0210             ptype = "mpnky"
0211         pass
0212         if GUI:
0213             fig.show()
0214         pass
0215         #fig.savefig()
0216
0217
0218     def positions_pvplt(self):
0219         axes = self.frame.axes
0220         if len(axes) == 2:
0221             self.positions_pvplt_2D()
0222         else:
0223             self.positions_pvplt_3D()
0224         pass
0225
0226     @classmethod
0227     def MakePVPlotter(cls):
0228         log.info("MakePVPlotter")
0229         pl = pv.Plotter(window_size=SIZE*2 )  # retina 2x ?
0230         return pl
0231
0232     def get_pv_plotter(self):
0233         if self.pl is None:
0234             pl = self.MakePVPlotter()
0235             self.pl = pl
0236         else:
0237             pl = self.pl
0238             log.info("using preexisting plotter")
0239         pass
0240         return pl
0241
0242
0243     def positions_pvplt_3D(self):
0244         """
0245         Could try to reconstruct solid surface from the point cloud of intersects
0246         https://docs.pyvista.org/api/core/_autosummary/pyvista.PolyDataFilters.reconstruct_surface.html#pyvista.PolyDataFilters.reconstruct_surface
0247         """
0248         pass
0249         pl = self.get_pv_plotter()
0250
0251         feat = self.feat
0252         upos = self.pos.upos   ## typically local frame
0253
0254         log.info("feat.unum %d " % feat.unum)
0255
0256         for idesc in range(feat.unum):
0257             uval, selector, label, color, skip, msg = feat(idesc)
0258             if skip: continue
0259             pos = upos[selector]
0260             print(msg)
0261             pl.add_points( pos[:,:3], color=color, point_size=10  )
0262         pass
0263         pl.enable_eye_dome_lighting()
0264         ## improves depth peception for point cloud, especially from a distance
0265         ## https://www.kitware.com/eye-dome-lighting-a-non-photorealistic-shading-technique/
0266         pl.show_grid()
0267
0268
0269     def positions_pvplt_2D(self):
0270         """
0271         * actually better to use set_position reset=True after adding points to auto get into ballpark
0272
0273         * previously always starts really zoomed in, requiring two-finger upping to see the intersects
0274         * following hint from https://github.com/pyvista/pyvista/issues/863 now set an adhoc zoom factor
0275
0276         Positioning the eye with a simple global frame y-offset causes distortion
0277         and apparent untrue overlaps due to the tilt of the geometry.
0278         Need to apply the central transform to the gaze vector to get straight on view.
0279
0280         https://docs.pyvista.org/api/core/_autosummary/pyvista.Camera.zoom.html?highlight=zoom
0281
0282         In perspective mode, decrease the view angle by the specified factor.
0283
0284         In parallel mode, decrease the parallel scale by the specified factor.
0285         A value greater than 1 is a zoom-in, a value less than 1 is a zoom-out.
0286         """
0287
0288         lim = self.gs.lim
0289         ugsc = self.gs.ugsc
0290         upos = self.pos.upos
0291
0292         feat = self.feat
0293
0294         pl = self.get_pv_plotter()
0295
0296         pl.add_text(self.topline, position="upper_left")
0297         pl.add_text(self.botline, position="lower_left")
0298         pl.add_text(self.frame.thirdline, position="lower_right")
0299
0300         log.info("positions_pvplt_2D feat.name %s " % feat.name )
0301
0302         for idesc in range(feat.unum):
0303             uval, selector, label, color, skip, msg = feat(idesc)
0304             if skip: continue
0305             pos = upos[selector]
0306             log.info(msg)
0307             pl.add_points( pos[:,:3], color=color )
0308         pass
0309
0310         if hasattr(self, 'x_lpos'):
0311             pvplt_add_contiguous_line_segments(pl, self.x_lpos[:,:3], point_size=25 )
0312         pass
0313
0314         if hasattr(self, 'simtrace_selection'):
0315             sts = self.simtrace_selection
0316             pvplt_simtrace_selection_line(pl, sts)
0317         pass
0318
0319         if not self.look_ce is None:
0320             pvplt_ce_multiple(pl, self.look_ce, axes=self.frame.axes)
0321         pass
0322
0323         if not self.look_epsilon is None:
0324             pvplt_ce(pl, self.look_epsilon, axes=self.frame.axes, color="yellow" )
0325         pass
0326
0327         show_genstep_grid = len(self.frame.axes) == 2 # too obscuring with 3D
0328         if show_genstep_grid and PVGRID:
0329             pl.add_points( ugsc[:,:3], color="white" )
0330         pass
0331
0332         pvplt_parallel_lines(pl, self.gs.lim, self.aa, self.frame.axes, self.look )
0333
0334         self.frame.pv_compose(pl, local=True)
0335
0336         cp = pl.show()
0337
0338         return cp
0339
0340
0341     def mp_show(self):
0342         """
0343         Fatal Python error: Segmentation fault
0344         """
0345         if hasattr(self, 'img'):
0346             mp.imshow(self.img)
0347             mp.show()
0348         pass
0349
0350
0351
0352 if __name__ == '__main__':
0353     logging.basicConfig(level=logging.INFO)
0354
0355