CSG/tests/cross2D_angle_range_without_trig.py

0001 #!/usr/bin/env python
0002 """
0003 cross2D_angle_range_without_trig.py
0004 =====================================
0005
0006 Working out how to use 2D cross product to select
0007 a range of phi directions without using arc-trig funcs
0008 like atan2.
0009
0010 """
0011
0012
0013 def cross2D( a , b ):
0014     """
0015     https://www.nagwa.com/en/explainers/175169159270/
0016
0017
0018     2D cross product definition::
0019
0020          A = Ax i + Ay j
0021
0022          B = Bx i + By j
0023
0024
0025          A ^ B =  (Ax By - Bx Ay ) k
0026
0027                 = |A||B| sin(AB_angle) k
0028
0029                 = sin(AB_angle) k         when A and B are normalized
0030
0031
0032     Anti-commutative::
0033
0034           A ^ B = - B ^ A
0035
0036     """
0037     assert a.shape[-1] == 3
0038     assert b.shape[-1] == 3
0039
0040     if a.ndim == 1 and b.ndim == 1:
0041         c = a[0]*b[1] - b[0]*a[1]
0042     elif a.ndim == 1 and b.ndim == 2:
0043         c = a[0]*b[:,1] - b[:,0]*a[1]
0044     elif a.ndim == 2 and b.ndim == 1:
0045         c = a[:,0]*b[1] - b[0]*a[:,1]
0046     else:
0047         c = None
0048     pass
0049     assert not c is None
0050     print("a.ndim : %d a.shape %s " % (a.ndim, str(a.shape)))
0051     print("b.ndim : %d b.shape %s " % (b.ndim, str(b.shape)))
0052     print("c.ndim : %d c.shape %s " % (c.ndim, str(c.shape)))
0053     return c
0054
0055
0056
0057
0058
0059 import numpy as np
0060 import pyvista as pv
0061
0062 def compose_XY( pl ):
0063      """
0064
0065
0066
0067                Y
0068                |
0069                U
0070                |
0071                |
0072                |
0073                L---------- X
0074               /
0075              /
0076             E
0077            /
0078           Z
0079
0080      """
0081      look = np.array( [0,0,0], dtype=np.float32 )
0082      up = np.array( [0,1,0], dtype=np.float32 )
0083      eye = np.array( [0,0,10], dtype=np.float32 )
0084      zoom = 1
0085
0086      pl.set_focus(    look )
0087      pl.set_viewup(   up )
0088      pl.set_position( eye, reset=False )
0089      pl.camera.Zoom(zoom)
0090      pl.show_grid()
0091
0092 if __name__ == '__main__':
0093
0094
0095      SIZE = np.array([1280, 720])
0096      pl = pv.Plotter(window_size=SIZE*2 )
0097
0098      #phiMode = "pacman"
0099      #phiMode = "pacman_small"
0100      #phiMode = "ppquad"
0101      phiMode = "nnquad"
0102      #phiMode = "uphemi"
0103
0104      if phiMode == "pacman":
0105          phiStart = 0.25
0106          phiDelta = 1.50
0107      elif phiMode == "pacman_small":
0108          phiStart = 0.10
0109          phiDelta = 1.80
0110      elif phiMode == "ppquad":
0111          phiStart = 0.00
0112          phiDelta = 0.25
0113      elif phiMode == "uphemi":
0114          phiStart = 0.00
0115          phiDelta = 1.00
0116      elif phiMode == "nnquad":
0117          phiStart = 1.00
0118          phiDelta = 0.5
0119      else:
0120          assert 0
0121      pass
0122
0123      phi0_pi = phiStart
0124      phi1_pi = phiStart+phiDelta
0125      phi_pi = np.linspace( phi0_pi,  phi1_pi, 100 )
0126      phi = np.pi*phi_pi
0127      cosPhi = np.cos( phi )
0128      sinPhi = np.sin( phi )
0129
0130
0131      ## edge vectors : as standin for geometry
0132      evec0 = np.array( [cosPhi[0], sinPhi[0], 0], dtype=np.float32 )
0133      evec1 = np.array( [cosPhi[-1], sinPhi[-1], 0], dtype=np.float32 )
0134
0135      ## edge lines from origin
0136      ll = np.zeros( (2, 2, 3), dtype=np.float32 )
0137      ll[:,0] = (0,0,0)
0138      ll[0,1] = evec0
0139      ll[1,1] = evec1
0140
0141      for i in range(len(ll)):
0142          pl.add_lines( ll[i].reshape(-1,3), color="blue" )
0143      pass
0144
0145      ## around the circular arc points
0146      pos = np.zeros( (len(phi), 3 ), dtype=np.float32 )
0147      pos[:,0] = cosPhi
0148      pos[:,1] = sinPhi
0149      pos[:,2] = 0.
0150      pl.add_points( pos, color="white" )
0151
0152
0153
0154      ## random phi in full range
0155      rnum = 100
0156      #rphi = 2.*np.random.random_sample(rnum)
0157      rphi = np.linspace( 0., 2., rnum)
0158
0159      rpos = np.zeros( (rnum, 3), dtype=np.float32 )
0160      rpos[:,0] = np.cos(np.pi*rphi)
0161      rpos[:,1] = np.sin(np.pi*rphi)
0162      rpos[:,2] = 0
0163
0164      s01 = cross2D(evec0, evec1)  # sine of angle evec0 -> evec1
0165      s0 = cross2D(evec0, rpos )   # sine of angle evec0 -> rpos
0166      s1 = cross2D(evec1, rpos )   # sine of angle evec1 -> rpos
0167
0168      ## hmm could flip the s1 definition
0169      ## to make the rpos inbetween case have both signs the same
0170
0171
0172      # check anti-commutativity
0173      s0_check = -cross2D(rpos, evec0) ; assert np.all( s0 == s0_check )
0174      s1_check = -cross2D(rpos, evec1) ; assert np.all( s1 == s1_check )
0175
0176
0177      """
0178      cross2D(evec0, evec1) > 0   means angle less than 1. (pi)
0179
0180      The sign-of-the-cross2D-sine tells you which "side of the line",
0181      so the *logical_and* of two such signs gives the region between the lines.
0182      This works fine when the angle between evec0 and evec1 is less than pi,
0183      selecting less than half the cake.
0184
0185
0186
0187                 Y
0188                 .        /  s0 > 0
0189                 .       / + + + + + + +
0190                 .      / + + + + + + +
0191                 .    evec0 + + + + + +
0192                 .    / + + + + + + +
0193                 .   / + + + + + + +
0194                 .  / + + + + + + +
0195                 . / + + + + + + +
0196                 ./ + + + + + + + +
0197                 O . . . . . . . . . . . . . . . X
0198                .| + + + + + + + +
0199               . | + + + + + + + +
0200              .  | + + + + + + + +
0201             .   | + + + + + + + +
0202            .    | + + + + + + + +
0203           .     | + + + + + + + +
0204          .     evec1  + + + + + +
0205                 | + + + + + + + +
0206                 |   s1 < 0
0207
0208
0209       When the angle between evec0 and evec1 is more than pi
0210       need to use *logical_or* of between the two signs
0211       to be more inclusive and thus greedily select more than half of the cake.
0212
0213        evec0
0214         \ + +   .
0215          \ + +  .
0216           \ + + .
0217            \ + +.
0218             \ + .
0219              \ +.+
0220               \ . +
0221                \. +
0222                 O + +
0223                 |. + +
0224                 | . + +
0225                 |+ . + +
0226                 |+ +. + +
0227                 |+ + . + +
0228                 |+ + +. + +
0229                 |+ + + . + +
0230              evec1
0231
0232      """
0233
0234      rr = np.zeros( (rnum,2,3), dtype=np.float32 )
0235      rr[:,0] = (0,0,0)
0236      rr[:,1] = rpos
0237
0238      #mask = np.logical_and( s0 > 0. , s1 > 0. )    # +Y quad
0239      #mask = np.logical_or(  s0 > 0. , s1 > 0. )     # not -Y quad
0240      #mask = np.logical_xor(  s0 > 0. , s1 > 0. )     # +X and -X
0241
0242      #mask = np.logical_and( s0 < 0. , s1 < 0. )    # -Y quad
0243      #mask = np.logical_or( s0 < 0. , s1 < 0. )      # not +Y quad
0244      #mask = np.logical_xor( s0 < 0. , s1 < 0. )      # +X and -X
0245
0246      #mask = np.logical_and( s0 < 0. , s1 > 0. )    # +X quad
0247      #mask = np.logical_or( s0 < 0. , s1 > 0. )      # not -X quad
0248      #mask = np.logical_xor( s0 < 0. , s1 > 0. )      # +Y and -Y
0249
0250      #mask = np.logical_and( s0 > 0. , s1 < 0. )     # -X quad           ## one to use for ppquad
0251      #mask = np.logical_or(  s0 > 0. , s1 < 0. )      # not +X quad    ## one to use for pacman, pacman_small
0252      #mask = np.logical_xor(  s0 > 0. , s1 < 0. )      # +Y and -Y
0253
0254      if s01 > 0.:   # angle less than pi
0255          mask = np.logical_and( s0 > 0. , s1 < 0. )
0256      else:          # angle greater than pi
0257          mask = np.logical_or( s0 > 0. , s1 < 0. )
0258      pass
0259
0260      urr = rr[mask]
0261
0262      for i in range(len(urr)):
0263          pl.add_lines( urr[i].reshape(-1,3), color="green" )
0264      pass
0265
0266
0267      compose_XY(pl)
0268      cp = pl.show()
0269
0270