eic-opticks/ana/pvprim1.py

0001 #!/usr/bin/env python
0002 """
0003 pvprim1.py
0004 =============
0005
0006 Places a cylinder on a sphere and then calulates (eye,look,up)
0007 viewpoint using tangential frame coordinates with origin
0008 at the target cylinder.
0009
0010 https://docs.pyvista.org/examples/00-load/create-geometric-objects.html
0011 """
0012 from collections import OrderedDict as odict
0013 import numpy as np
0014 import pyvista as pv
0015
0016 _white = "ffffff"
0017 _red = "ff0000"
0018 _green = "00ff00"
0019 _blue = "0000ff"
0020
0021 DTYPE = np.float64
0022 SIZE = np.array([1280, 720])
0023
0024
0025 def make_transform(r, t, p ):
0026     """
0027
0028            -Y    Z            -T   +P
0029             .  /               .  /
0030             . /                . /
0031             ./                 ./
0032             +------ X          +------ +R
0033             |                 .|
0034             |                . |
0035             |               .  |
0036             Y             -P   +T
0037
0038
0039
0040     With the below transform get units vector mapping of (X, Y, Z)  =>  (R, T, P)::
0041
0042        c2s = np.array([
0043                      [ np.sin(theta)*np.cos(phi) , np.sin(theta)*np.sin(phi) , np.cos(theta)  ,  0. ],
0044                      [ np.cos(theta)*np.cos(phi) , np.cos(theta)*np.sin(phi) , -np.sin(theta) ,  0. ],
0045                      [ -np.sin(phi)              , np.cos(phi)               , 0.             ,  0. ],
0046                      [ 0.                        , 0.                        , 0.             ,  1. ]
0047                    ], dtype=DTYPE )
0048
0049
0050     (red)   +X arrows  =>  +R radial-outwards
0051
0052     (green) +Y arrows  =>  +T theta-tangent (North to South direction)
0053
0054     (blue)  +Z arrows  =>  +P phi-tangent (West to East direction)
0055
0056
0057     It is best to get used to that (R P T) tangential frame rather
0058     than suffering the confusion of trying to rotate to something else like (P T R)
0059     which might initially seem more intuitive.
0060
0061
0062        Z                      X
0063        |  Y                   |  Y
0064        | /                    | /
0065        |/                     |/
0066        +----- X        Z -----+
0067
0068     """
0069     radius = r
0070     theta = t*np.pi
0071     phi = p*np.pi
0072
0073     tx = radius*np.sin(theta)*np.cos(phi)
0074     ty = radius*np.sin(theta)*np.sin(phi)
0075     tz = radius*np.cos(theta)
0076
0077     tra = np.array([[1,   0,  0,  0],
0078                    [0,   1,  0,  0],
0079                    [0,   0,  1,  0],
0080                    [tx, ty, tz,  1]])
0081
0082     rot = np.array([
0083                      [ np.sin(theta)*np.cos(phi) , np.sin(theta)*np.sin(phi) , np.cos(theta)  ,  0. ],
0084                      [ np.cos(theta)*np.cos(phi) , np.cos(theta)*np.sin(phi) , -np.sin(theta) ,  0. ],
0085                      [ -np.sin(phi)              , np.cos(phi)               , 0.             ,  0. ],
0086                      [ 0.                        , 0.                        , 0.             ,  1. ]
0087                    ], dtype=DTYPE )
0088
0089
0090     tra_rot = np.array([
0091                      [ np.sin(theta)*np.cos(phi) , np.sin(theta)*np.sin(phi) , np.cos(theta)  ,  0. ],
0092                      [ np.cos(theta)*np.cos(phi) , np.cos(theta)*np.sin(phi) , -np.sin(theta) ,  0. ],
0093                      [ -np.sin(phi)              , np.cos(phi)               , 0.             ,  0. ],
0094                      [ tx                        , ty                        , tz             ,  1. ]
0095                    ], dtype=DTYPE )
0096
0097
0098
0099     #return np.dot(rot, tra )
0100     return tra_rot
0101
0102
0103 class TangentialFrame(object):
0104     """
0105     Tangential frame origin is a fixed point on the surface of a sphere.
0106     For example a point on the X axis (like null_island on equator)::
0107
0108         xyz :(radius,0,0)
0109         theta:0.5 phi:0
0110         rtpw:(0,0,0,1)
0111
0112
0113     TODO: no extent scaling yet
0114     TODO: pyvista view does not appear until interact with the view somehow
0115
0116     """
0117     def __init__(self, rtp, extent=1):
0118         r,t,p = rtp
0119         radius = r
0120         theta = t*np.pi
0121         phi = p*np.pi
0122
0123         x = radius*np.sin(theta)*np.cos(phi)
0124         y = radius*np.sin(theta)*np.sin(phi)
0125         z = radius*np.cos(theta)
0126         w = 1.
0127
0128         rot_tra = np.array([
0129                      [ np.sin(theta)*np.cos(phi)        , np.sin(theta)*np.sin(phi)        , np.cos(theta)         ,  0. ],
0130                      [ np.cos(theta)*np.cos(phi)        , np.cos(theta)*np.sin(phi)        , -np.sin(theta)        ,  0. ],
0131                      [ -np.sin(phi)                     , np.cos(phi)                      , 0.                    ,  0. ],
0132                      [ radius*np.sin(theta)*np.cos(phi) , radius*np.sin(theta)*np.sin(phi) , radius*np.cos(theta)  ,  1. ]
0133                      ], dtype=DTYPE )
0134
0135         self.rot_tra = rot_tra
0136         self.xyzw = np.array([x,y,z,w], dtype=DTYPE)
0137         self.ce = np.array([x,y,z,extent], dtype=DTYPE)
0138
0139
0140     def tangential_to_conventional(self, rtpw):
0141         """
0142         aka tangential_to_conventional
0143
0144                    z:p
0145                    |
0146                    |  y:t
0147                    | /
0148                    |/
0149                    +-----x:r
0150
0151
0152         :param rtpw: tangential cartesian coordinate in frame with origin at point on surface of sphere
0153                       and with tangent-to-sphere and normal-to-sphere as "theta" "phi" unit vectors
0154         :return  xyzw: conventional cartesian coordinate in frame with  origin at center of sphere
0155         """
0156         return np.dot( rtpw, self.rot_tra )
0157
0158
0159
0160
0161 if __name__ == '__main__':
0162
0163
0164     pl = pv.Plotter(window_size=SIZE*2 )
0165
0166     edges = False
0167     radius = 10.
0168
0169     sphere = pv.Sphere()   # curious choice of radius 0.5 for Sphere
0170     tr_sphere = np.array([[2*radius,          0,        0, 0],
0171                           [0,          2*radius,        0, 0],
0172                           [0,                 0, 2*radius, 0],
0173                           [0,                 0,        0, 1]])
0174
0175
0176
0177
0178     rtp = odict()
0179     rtp["null_island"] = np.array( [radius, 0.5, 0.0] )
0180     rtp["midlat"] = np.array(      [radius, 0.25, 0.0] )
0181
0182     #k = "null_island"
0183     k = "midlat"
0184
0185     tf = TangentialFrame(rtp[k], extent=1.)
0186     print("tf.xyzw %s " % str(tf.xyzw))
0187
0188     tr_cyl = tf.rot_tra
0189
0190     #view = "from_above"
0191     view = "from_side"
0192
0193     if view == "from_above":
0194         look_rtpw = np.array( [0.,  0.,  0., 1.] )
0195         eye_rtpw  = np.array( [1.,  0.,  0., 1.] )
0196         up_rtpw   = np.array( [0., -1.,  0., 0.] )
0197     elif view == "from_side":
0198         look_rtpw = np.array( [0.,   0.,  0., 1.] )
0199         eye_rtpw  = np.array( [0.,   5.,  0., 1.] )
0200         up_rtpw   = np.array( [0.,   0.,  1., 0.] )
0201     else:
0202         assert 0
0203     pass
0204
0205     look_xyzw = tf.tangential_to_conventional( look_rtpw )
0206     eye_xyzw  = tf.tangential_to_conventional( eye_rtpw )
0207     up_xyzw   = tf.tangential_to_conventional( up_rtpw )
0208     zoom = 1.
0209
0210     print("look_xyzw : %s " % (str(look_xyzw)))
0211     print("eye_xyzw : %s " % (str(eye_xyzw)))
0212     print("up_xyzw : %s " % (str(up_xyzw)))
0213     print("zoom : %s " % (str(zoom)))
0214
0215     sphere.transform(tr_sphere)
0216     pl.add_mesh(sphere, color=_blue, show_edges=edges, style="wireframe")
0217     pl.show_grid()
0218
0219     cyl = pv.Cylinder(direction=(1,0,0))
0220     cyl.transform(tr_cyl.T)
0221     pl.add_mesh(cyl, color=_white, show_edges=True )
0222
0223     for t in np.linspace(0,1,10):
0224         for p in np.linspace(0,2,10):
0225             tr = make_transform( radius, t, p )
0226
0227             x_arrow = pv.Arrow(direction=(1,0,0))  # X->R
0228             y_arrow = pv.Arrow(direction=(0,1,0))  # Y->T
0229             z_arrow = pv.Arrow(direction=(0,0,1))  # Z->P
0230
0231             x_arrow.transform(tr.T)
0232             y_arrow.transform(tr.T)
0233             z_arrow.transform(tr.T)
0234
0235             pl.add_mesh(x_arrow, color=_red , show_edges=False)
0236             pl.add_mesh(y_arrow, color=_green, show_edges=False)
0237             pl.add_mesh(z_arrow, color=_blue, show_edges=False)
0238
0239             plane = pv.Plane(direction=(1,0,0))   # ends up tangential to sphere
0240             plane.transform(tr.T)
0241             pl.add_mesh(plane, color=_white, show_edges=edges)
0242         pass
0243     pass
0244
0245     pl.set_focus(    look_xyzw )
0246     pl.set_viewup(   up_xyzw[:3] )
0247     pl.set_position( eye_xyzw )    # reset=True changes to get everything into frame, so its at odds with trying to carefully control the viewpoint
0248     pl.camera.Zoom(zoom)
0249     #pl.update()
0250
0251     cpos = pl.show()
0252     #pl.update()
0253     print(cpos)
0254