eic-opticks/ana/cylinder.py

0001 #!/usr/bin/env python
0002 """
0003 cylinder.py
0004 =============
0005
0006
0007 """
0008
0009 import numpy as np
0010 #import pyvista as pv
0011 #pv.set_plot_theme("dark")
0012 from opticks.ana.pvplt import *
0013 from opticks.ana.make_rotation_matrix import make_rotation_matrix
0014
0015 DTYPE = np.float64
0016
0017 X = np.array([1,0,0,0], dtype=DTYPE)
0018 Y = np.array([0,1,0,0], dtype=DTYPE)
0019 Z = np.array([0,0,1,0], dtype=DTYPE)
0020 O = np.array([0,0,0,1], dtype=DTYPE)
0021
0022
0023 def make_transforms_cube_corners():
0024     trs = np.zeros( (8,4,4) )
0025     for i in range(len(trs)):
0026         m = np.eye(4)
0027         tx = sc if i & 0x1 else -sc
0028         ty = sc if i & 0x2 else -sc
0029         tz = sc if i & 0x4 else -sc
0030         m[3] = (tx,ty,tz,1)
0031         trs[i] = m.T
0032     pass
0033     return trs
0034
0035 def make_transforms_around_cylinder(radius, halfheight, num_ring=10, num_in_ring=16):
0036     zz = np.linspace( -halfheight, halfheight, num_ring )
0037
0038     phi = np.linspace( 0, 2*np.pi, num_in_ring+1 )[:-1]
0039     xy_outwards = np.zeros( (len(phi), 4) )
0040     xy_outwards[:,0] = np.cos(phi)
0041     xy_outwards[:,1] = np.sin(phi)
0042     xy_outwards[:,2] = 0.
0043     xy_outwards[:,3] = 0.
0044
0045     assert len(zz) == num_ring
0046     assert len(xy_outwards) == num_in_ring
0047
0048     num_tr = num_ring*num_in_ring
0049     trs = np.zeros( (num_tr, 4, 4) )
0050     for i in range(num_ring):
0051         z = zz[i]
0052         for j in range(num_in_ring):
0053             outwards = xy_outwards[j]
0054             m4 = make_rotation_matrix(Z, outwards)
0055             # HUH: expected -outwards to point the arrows inwards
0056             # presumably this is due to the transpose
0057             # done for pyvista
0058             pos = outwards*radius
0059             pos[2] = z
0060             m4[3,:3] = pos[:3]
0061             ## HMM: somewhat dirty stuffing the translation in like that
0062             ## as its a bit unclear in what multiplication order it corresponds to
0063             ## should really multiply matrices
0064             idx = i*num_in_ring + j
0065             trs[idx] = m4.T       ## seems pyvista needs transposed
0066         pass
0067     pass
0068     return trs
0069
0070
0071 if __name__ == '__main__':
0072
0073     sc = 5
0074     radius = sc
0075     halfheight = sc
0076     trs0 = make_transforms_around_cylinder(radius, halfheight)
0077     trs1 = np.load("/tmp/SPlace_test/AroundCylinder1.npy")   ## with the flip
0078     dt = np.abs( trs0 - trs1 ).max()
0079     assert dt < 1e-9
0080
0081     trs2 = np.load("/tmp/SPlace_test/AroundSphere1.npy")
0082     trs = trs2
0083
0084     arrows = []
0085     for i in range(len(trs)):
0086         arr = pv.Arrow(direction=Z[:3])
0087         arr.transform(trs[i])
0088         arrows.append(arr)
0089     pass
0090     cyl = pv.Cylinder(center=O[:3], direction=Z[:3], radius=radius, height=2*halfheight)
0091
0092     pl = pvplt_plotter()
0093     pl.add_mesh(cyl, style="wireframe")
0094     for arr in arrows:
0095         pl.add_mesh(arr)
0096     pass
0097     pl.show()
0098
0099