eic-opticks/ana/makeflight.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 makeflight.py : creates eye-look-up input NumPy arrays used by okc/FlightPath.cc
0023 ===================================================================================
0024
0025 See docs/misc/making_flightpath_raytrace_movies.rst
0026
0027 This creates input eye-look-up-ctrl .npy files read by optickscore/FlightPath
0028
0029 ::
0030
0031     ipython -i -- makeflight.py
0032
0033     ana/makeflight.sh
0034
0035
0036 """
0037 import os, inspect, logging, argparse, numpy as np
0038 log = logging.getLogger(__name__)
0039
0040 dtype = np.float32
0041
0042 class Flight(object):
0043     DEFAULT_BASE = os.path.expanduser("~/.opticks/flight")
0044
0045     @classmethod
0046     def EyeLine(cls, axis='X', scale=1, start=500, stop=1, steps=32):
0047         """
0048         Move eye in straight line along an axis
0049
0050
0051                 Z
0052                 |   Y
0053                 |  /
0054                 | /
0055                 |/
0056                 +-----<---X
0057
0058         """
0059         method = inspect.currentframe().f_code.co_name
0060         name = "{method}{axis}".format(**locals())
0061         log.info("axis %s scale %s steps %s name %s " % (axis, scale, steps, name))
0062
0063         f = cls.Make( name, steps)
0064         v = scale*np.linspace(start, stop,steps )
0065
0066         f.l[:] = [0,0,0,1]
0067
0068         f.e[:,0] = v if axis == 'X' else 0
0069         f.e[:,1] = v if axis == 'Y' else 0
0070         f.e[:,2] = v if axis == 'Z' else 0
0071         f.e[:,3] = 1
0072
0073         f.u[:] = [0,0,1,0] if axis in "XY" else [1,0,0,0]
0074
0075         return f
0076
0077     @classmethod
0078     def Roundabout(cls, plane='XY', scale=1, steps=32):
0079         """
0080
0081         Move eye in circle in XY/XZ plane whilst looking towards the center, up +Z/+Y
0082
0083         **plane XY, +Z up**                          **plane XZ, +Y up**
0084
0085                                                       Y
0086                  Y                                    |
0087                 3|                                    |
0088                  |   2                                |
0089                  |                                    |
0090                  |      1                             |
0091                  |                                    O--------0----X
0092                  O--------0----X                     /
0093                 /                                   /       1
0094                /                                       2
0095               /                                   Z
0096              Z
0097
0098         Note that the flightpath can be scaled after loading
0099         into Opticks executables using the --flightpathscale option.
0100
0101         Its generally more flexible to change scale that way, avoiding
0102         the need to recreate flightpath.npy files for simple scale changes.
0103         """
0104
0105         method = inspect.currentframe().f_code.co_name
0106         name = "{method}{plane}".format(**locals())
0107
0108         log.info("plane %s scale %s steps %s name %s " % (plane, scale, steps, name))
0109
0110         f = cls.Make( name, steps)
0111         ta = np.linspace( 0, 2*np.pi, steps )
0112         st = np.sin(ta)
0113         ct = np.cos(ta)
0114         n = len(ta)
0115
0116         if plane == 'XY':        # starts at  (scale, 0, 0) -> (0, scale, 0)
0117             f.e[:,0] = ct*scale  # X
0118             f.e[:,1] = st*scale  # Y
0119             f.e[:,2] = 0
0120             f.e[:,3] = 1
0121
0122             f.l[:] = [0,0,0,1]   # always looking at center
0123             f.u[:] = [0,0,1,0]   # always up +Z
0124
0125         elif plane == 'XZ':       # starts at  (scale, 0, 0) -> (0, 0, scale)
0126
0127             f.e[:,0] = ct*scale  # X
0128             f.e[:,1] = 0
0129             f.e[:,2] = st*scale  # Z
0130             f.e[:,3] = 1
0131
0132             f.l[:] = [0,0,0,1]   # always looking at center
0133
0134             f.u[:,0] = -st*scale
0135             f.u[:,1] = 0
0136             f.u[:,2] = ct*scale
0137             f.u[:,3] = 0
0138
0139         else:
0140             pass
0141         pass
0142         return f
0143
0144     @classmethod
0145     def Path(cls, name):
0146         return os.path.join(cls.DEFAULT_BASE, "%s.npy" % name)
0147
0148     @classmethod
0149     def Make(cls, name, n):
0150         eluc = np.zeros( (n,4,4), dtype=np.float32)
0151         return cls(name, eluc)
0152
0153     @classmethod
0154     def Load(cls, name):
0155         path = cls.Path(name)
0156         eluc = np.load(path)
0157         return cls(name, eluc)
0158
0159     @classmethod
0160     def Combine(cls, names, combined_name):
0161         arrs = []
0162         for name in names:
0163             path = cls.Path(name)
0164             arr = np.load(path)
0165             print( " %15s : %s " % (str(arr.shape), path ))
0166             arrs.append(arr)
0167         pass
0168         return cls.CombineArrays(arrs, combined_name)
0169
0170     @classmethod
0171     def CombineArrays(cls, arrs, combined_name):
0172         eluc = np.concatenate(tuple(arrs))
0173         return cls(combined_name, eluc)
0174
0175     def __init__(self, name, eluc ):
0176         self.name = name
0177         self.eluc = eluc
0178
0179     e = property(lambda self:self.eluc[:,0,:4] )
0180     l = property(lambda self:self.eluc[:,1,:4] )
0181     g = property(lambda self:self.l - self.e)
0182     u = property(lambda self:self.eluc[:,2,:4] )
0183     c = property(lambda self:self.eluc[:,3,:4] )
0184
0185     e3 = property(lambda self:self.eluc[:,0,:3] )
0186     g3 = property(lambda self:self.eluc[:,1,:3] - self.eluc[:,0,:3] )
0187     u3 = property(lambda self:self.eluc[:,2,:3] )
0188     r3 = property(lambda self:np.cross( self.g3, self.u3 ))
0189
0190     def save(self):
0191         path = self.Path(self.name)
0192         fold = os.path.dirname(path)
0193         if not os.path.isdir(fold):
0194             log.info("creating directory %s " % fold)
0195             os.makedirs(fold)
0196         pass
0197         log.info("saving to %s " % path )
0198         np.save(path, self.eluc )
0199
0200     def print_cmds(self):
0201         print(self.c.copy().view("|S2"))
0202
0203     def __len__(self):
0204         return len(self.eluc)
0205
0206     def __repr__(self):
0207         return "Flight %s eluc.shape %s " % (self.name, str(self.eluc.shape))
0208
0209
0210     def quiver_plot(self, ax, sc):
0211         e = self.e
0212         l = self.l
0213         u = self.u
0214         g = l - e
0215
0216         x = sc*e[:,0]
0217         y = sc*e[:,1]
0218         z = sc*e[:,2]
0219
0220         u0 = g[:, 0]
0221         v0 = g[:, 1]
0222         w0 = g[:, 2]
0223
0224         u1 = u[:, 0]
0225         v1 = u[:, 1]
0226         w1 = u[:, 2]
0227
0228         #ax.plot( x,z )
0229         ax.quiver( x, y, z, u0, v0, w0  )
0230         ax.quiver( x, y, z, u1, v1, w1  )
0231
0232         labels = False
0233         if labels:
0234             for i in range(len(e)):
0235                 ax.text( x[i], y[i], z[i], i , "z" )
0236             pass
0237         pass
0238
0239
0240 def parse_args(doc, **kwa):
0241     np.set_printoptions(suppress=True, precision=3, linewidth=200)
0242     parser = argparse.ArgumentParser(doc)
0243     parser.add_argument( "--level", default="info", help="logging level" )
0244     parser.add_argument( "--steps", default=32, type=int, help="Number of steps in flightpath that are interpolated between in InterpolatedView " )
0245     parser.add_argument( "--scale", default=1, type=float, help="scale of the flightpath, for example the radius for circles" )
0246     parser.add_argument( "--start", default=100, type=float, help="start input for flightpath, for example EyeLine start along an axis" )
0247     parser.add_argument( "--stop",  default=1,  type=float, help="stop input for flightpath, for example EyeLine stop along an axis" )
0248
0249     args = parser.parse_args()
0250     fmt = '[%(asctime)s] p%(process)s {%(pathname)s:%(lineno)d} %(levelname)s - %(message)s'
0251     logging.basicConfig(level=getattr(logging,args.level.upper()), format=fmt)
0252     return args
0253
0254 if __name__ == '__main__':
0255     pass
0256     np.set_printoptions(suppress=True)
0257     args = parse_args(__doc__)
0258
0259     ff = {}
0260     for p in ['XY', 'XZ' ]:
0261         ff[p] = Flight.Roundabout(plane=p, scale=args.scale, steps=args.steps)
0262     pass
0263
0264     p = 'XY_XZ'
0265     ff[p] = Flight.CombineArrays( [ff['XY'].eluc, ff['XZ'].eluc], 'Roundabout%s' % p )
0266
0267     for axis in "XYZ":
0268         ff[axis] = Flight.EyeLine(axis, scale=args.scale, steps=args.steps)
0269     pass
0270
0271     for p in ff.keys():
0272         f = ff[p]
0273         print(f)
0274         f.save()
0275     pass
0276