eic-opticks/CSG/phicut.py

0001 #!/usr/bin/env python
0002 """
0003 phicut.py
0004 ============
0005
0006
0007 ~/o/CSG/phicut.sh
0008
0009
0010 """
0011 import numpy as np
0012 import logging
0013 from opticks.CSG.scan import Scan
0014
0015 log = logging.getLogger(__name__)
0016
0017 try:
0018     import pyvista as pv
0019 except ImportError:
0020     pv = None
0021 pass
0022
0023
0024 class HalfPlane(object):
0025
0026     @classmethod
0027     def phi_quadrant(cls, phi):
0028         """
0029         This way of defining the quadrant avoids precision of pi issues.
0030         Yes, but its better to use the same definition that can be applied to the
0031         intersect positions without having to calculate phi from x,y.
0032         Problems likely at boundaries between quadrants
0033
0034                          0.5
0035                  0.75     |     0.25
0036                           |
0037                           |
0038                1.0  ------+------- 0.0 , 2.0
0039                           |
0040                           |
0041                 1.25      |     1.75
0042                          1.5
0043
0044
0045            cosPhi < 0  |  cosPhi > 0
0046                        |
0047                -+      |     ++
0048                01      |     11
0049                        |              sinPhi > 0
0050                - - - - +----------------------------
0051                        |              sinPhi < 0
0052                --      |     +-
0053                00      |     10
0054
0055         """
0056         if phi >= 0 and phi < 0.5:
0057             quadrant = 3
0058         elif phi >= 0.5 and phi < 1.0:
0059             quadrant = 1
0060         elif phi >= 1.0 and phi < 1.5:
0061             quadrant = 0
0062         elif phi >= 1.5 and phi <= 2.0:
0063             quadrant = 2
0064         else:
0065             quadrant = -1
0066         pass
0067         return quadrant
0068
0069     @classmethod
0070     def xy_quadrant(cls, x, y ):
0071         xpos = int(x >= 0.)
0072         ypos = int(y >= 0.)
0073         return 2*xpos + ypos
0074
0075     @classmethod
0076     def test_quadrant(cls, n=17):
0077         """
0078         In [3]: np.linspace(0,2,17)
0079         Out[3]: array([0.   , 0.125, 0.25 , 0.375, 0.5  , 0.625, 0.75 , 0.875, 1.   , 1.125, 1.25 , 1.375, 1.5  , 1.625, 1.75 , 1.875, 2.   ])
0080         """
0081         phi = np.linspace(0, 2, n)
0082         for p in phi:
0083             x = np.cos(p*np.pi)
0084             y = np.sin(p*np.pi)
0085             _phi_quadrant = cls.phi_quadrant(p)
0086             _xy_quadrant = cls.xy_quadrant(x,y)
0087             print( " p %10.4f  phi_quadrant %d   cosPhi %10.4f sinPhi %10.4f  xy_quadrant %d " % (p, _phi_quadrant, x, y, _xy_quadrant ))
0088         pass
0089
0090     def __init__(self, phi, debug=False ):
0091         cosPhi, sinPhi = np.cos(phi*np.pi), np.sin(phi*np.pi)
0092
0093         #n_quadrant = self.phi_quadrant(phi)
0094         n_quadrant = self.xy_quadrant(cosPhi, sinPhi)
0095         n = np.array( [  sinPhi, -cosPhi, 0. ]  )
0096
0097         if debug:
0098             log.info( " phi %10.3f cosPhi %10.3f sinPhi %10.3f  n_quadrant %d " % (phi, cosPhi, sinPhi, n_quadrant ))
0099         pass
0100
0101         self.n_quadrant = n_quadrant
0102         self.n = n
0103         self.sinPhi = sinPhi
0104         self.cosPhi = cosPhi
0105         self.debug = debug
0106
0107     def intersect(self, ray_origin, ray_direction, t_min ):
0108         """
0109
0110
0111         phi=1  -(*)- - - - - O-----*--+-----*------ X   phi = 0.
0112                 /                 /   |      \
0113                /                 /    |       \
0114               /                 /     n        \
0115                                /                \
0116                               /
0117
0118         How to generalize disqualification ? Possible techniques
0119
0120         1. quadrant match
0121
0122            * check the signs of (x,y) at the intersect
0123              corresponds to signs of (cosPhi, sinPhi) of the plane
0124
0125            * issues with this at phi 1.5, 2.0 : presumably from depending on the sign of small values
0126            * using "side > 0" avoids the problems and is simpler
0127
0128         2. rotate intersect (x,y) by -phi, so sign of xrot can be required +ve
0129
0130            * xrot = x*cosPhi + y*sinPhi
0131            * yrot = -x*sinPhi + y*cosPhi
0132
0133         3. intersect must fulfil p.n = 0 yes but disqualified does too.. so no use
0134
0135         4. dot product between the vector from origin to intersect and the vector
0136            (cosPhi, sinPhi, 0 ) that is within the plane and identifies the phi
0137            should be +ve
0138
0139            * side = x*cosPhi + y*sinPhi
0140            * actually the maths of this is the same as backwards rotation
0141              its just thinking of the same thing in two different frames
0142
0143
0144         """
0145         n = self.n
0146         sinPhi = self.sinPhi
0147         cosPhi = self.cosPhi
0148         debug = self.debug
0149
0150         dn = np.dot(ray_direction, n )
0151         on = np.dot(ray_origin,    n )
0152         t_cand = t_min if dn == 0. else -on/dn
0153
0154         xyz = ray_origin + t_cand*ray_direction
0155
0156         x = ray_origin[0] + t_cand*ray_direction[0]
0157         y = ray_origin[1] + t_cand*ray_direction[1]
0158         side0 = x*cosPhi + y*sinPhi
0159
0160         side =     ray_origin[0]*cosPhi + ray_origin[1]*sinPhi + ( ray_direction[0]*cosPhi + ray_direction[1]*sinPhi )*t_cand
0161
0162         sidecheck = True
0163         if sidecheck:
0164             valid_intersect = t_cand > t_min and side > 0.
0165         else:
0166             valid_intersect = t_cand > t_min
0167         pass
0168
0169         if debug:
0170             afmt_ = lambda a:"%7.3f %7.3f %7.3f" % (a[0], a[1], a[2])
0171             print(" ray_origin %s ray_direction %s xyz %s t_cand %7.3f  side0 %10.4f side %10.4f valid_intersect %d" %
0172                (afmt_(ray_origin), afmt_(ray_direction), afmt_(xyz), t_cand, side0, side, valid_intersect ))
0173         pass
0174         if valid_intersect:
0175             isect = np.zeros(4)
0176             isect[:3] = n
0177             isect[3] = t_cand
0178         else:
0179             isect = None
0180         pass
0181         return isect
0182
0183
0184 class PhiCut(object):
0185     @classmethod
0186     def plane(cls, phi):
0187         return np.array( [np.sin(np.pi*phi), -np.cos(np.pi*phi), 0, 0 ], dtype=np.float32 )
0188
0189     def __init__(self, phiStart, phiDelta ):
0190         phi0 = phiStart
0191         phi1 = phiStart + phiDelta
0192
0193         cosPhi0, sinPhi0 = np.cos(phi0*np.pi), np.sin(phi0*np.pi)
0194         cosPhi1, sinPhi1 = np.cos(phi1*np.pi), np.sin(phi1*np.pi)
0195         n0 = np.array( [  sinPhi0, -cosPhi0, 0. ]  )
0196         n1 = np.array( [ -sinPhi1,  cosPhi1, 0. ]  )
0197         self.n0 = n0
0198         self.n1 = n1
0199         self.cosPhi0 = cosPhi0
0200         self.sinPhi0 = sinPhi0
0201         self.cosPhi1 = cosPhi1
0202         self.sinPhi1 = sinPhi1
0203
0204
0205     def intersect(self, ray_origin, ray_direction, t_min ):
0206         """
0207         TODO: add handling of pure-z rays and then mixed 3D rays, add unbounded_exit signally
0208
0209
0210         These are currently full planes need to disqualify half of them.
0211         What exactly is the convention for the allowed intersects ?
0212
0213         The solid lines below denote the valid half planes.
0214         For the infinite plane intersection points [0] [1] (0) (1)
0215
0216
0217                                p1 [cosPhi1, sinPhi1]
0218                                 /
0219                                /
0220                              (1)
0221                              /:
0222                             / :
0223                            /  :
0224            . . .{0}. . .  +--(0)--------- p0   [ cosPhi0, sinPhi0 ]
0225                  :       .    :
0226                  :      .     :
0227                  :     .      :
0228                  :    .       :
0229                  :   .        :
0230                  :  .         :
0231                  : .          :
0232                  :.           :
0233                 {1}           :
0234                 .:            :
0235                . :            :
0236               .  :            :
0237                  o            o
0238
0239
0240         """
0241         cosPhi0 = self.cosPhi0
0242         sinPhi0 = self.sinPhi0
0243         cosPhi1 = self.cosPhi1
0244         sinPhi1 = self.sinPhi1
0245
0246         # dot products for direction and origin with::
0247         # normal0 [  sinPhi0, -cosPhi0, 0. ]
0248         # normal1 [ -sinPhi1,  cosPhi1, 0. ]
0249
0250         d_n0 = ray_direction[0]*sinPhi0    + ray_direction[1]*(-cosPhi0)
0251         d_n1 = ray_direction[0]*(-sinPhi1) + ray_direction[1]*(cosPhi1)
0252         o_n0 = ray_origin[0]*sinPhi0       + ray_origin[1]*(-cosPhi0)
0253         o_n1 = ray_origin[0]*(-sinPhi1)    + ray_origin[1]*(cosPhi1)
0254         t0 = t_min if d_n0 == 0. else -o_n0/d_n0    # intersect distance to halfPlane0
0255         t1 = t_min if d_n1 == 0. else -o_n1/d_n1    # intersect distance to halfPlane1
0256
0257         # Disqualify intersects with the wrong halves of the planes.
0258         # When the dot product of the hit point with the unit vector
0259         # within the plane (in valid direction) is positive it means the
0260         # intersect is with the valid half of the plane.  When negative
0261         # are intersecting with the non-allowed half of the plane.
0262         #
0263         # sideI = (ray_origin + ray_direction * tI) · (cosPhiI, sinPhiI)    I={0,1}
0264
0265         side0 = ray_origin[0]*cosPhi0 + ray_origin[1]*sinPhi0 + ( ray_direction[0]*cosPhi0 + ray_direction[1]*sinPhi0 )*t0
0266         side1 = ray_origin[0]*cosPhi1 + ray_origin[1]*sinPhi1 + ( ray_direction[0]*cosPhi1 + ray_direction[1]*sinPhi1 )*t1
0267         if side0 < 0.: t0 = t_min
0268         if side1 < 0.: t1 = t_min
0269
0270         t_near = min(t0,t1)
0271         t_far  = max(t0,t1)
0272
0273         #t_cand = t_near > t_min  ?  t_near : ( t_far > t_min ? t_far : t_min ) ;
0274         t_cand = t_min
0275         if t_near > t_min:
0276             t_cand = t_near
0277         else:
0278             if t_far > t_min:
0279                 t_cand = t_far
0280             else:
0281                 t_cand = t_min
0282             pass
0283         pass
0284
0285         valid_intersect = t_cand > t_min
0286
0287         if valid_intersect:
0288             isect = np.zeros(4)
0289             isect[0]= -sinPhi1 if t_cand == t1 else sinPhi0
0290             isect[1]=  cosPhi1 if t_cand == t1 else -cosPhi0
0291             isect[2]=  0.
0292             isect[3] = t_cand
0293         else:
0294             isect = None
0295         pass
0296         return isect
0297
0298
0299 def check_normals():
0300     #phi = np.linspace(0,2,100)
0301     phi = np.linspace(0,2,10)
0302
0303     planes = np.zeros( (len(phi), 4), dtype=np.float32 )
0304     pos = np.zeros( (len(phi), 3), dtype=np.float32 )
0305
0306     for i in range(len(phi)):
0307         planes[i] = PhiCut.plane(phi[i])
0308         pos[i] = (10*np.cos(phi[i]*np.pi), 10*np.sin(phi[i]*np.pi), 0 )
0309     pass
0310
0311     size = np.array([1280, 720])
0312     pl = pv.Plotter(window_size=size*2 )
0313     #pl.camera.ParallelProjectionOn()
0314     pl.view_xy()
0315     pl.add_arrows( pos,  planes[:,:3],  color="red" )   # suitable for phi0
0316     pl.add_arrows( pos, -planes[:,:3],  color="blue" )  # suitable for phi1
0317
0318     look = (0,0,0)
0319     up = (0,1,0)
0320     eye = (0,0, 10)
0321
0322     pl.set_focus(    look )
0323     pl.set_viewup(   up )
0324     pl.set_position( eye, reset=True )   ## for reset=True to succeed to auto-set the view, must do this after add_points etc..
0325     pl.camera.Zoom(1)
0326     pl.show_grid()
0327     cp = pl.show()
0328
0329
0330
0331 if __name__ == '__main__':
0332     logging.basicConfig(level=logging.INFO)
0333
0334     #HalfPlane.test_quadrant()
0335
0336     #geom = HalfPlane( 0.0 )
0337     #geom = HalfPlane( 0.25 )
0338     #geom = HalfPlane( 0.50 )
0339     #geom = HalfPlane( 0.75 )
0340     #geom = HalfPlane( 1.00 )
0341     #geom = HalfPlane( 1.25 )
0342     #geom = HalfPlane( 1.5, debug=False)
0343     #geom = HalfPlane( 1.75, debug=False)
0344     #geom = HalfPlane( 2.0, debug=False)
0345
0346     #geom = PhiCut( 0.25, 0.1 )
0347     geom = PhiCut( 0.0, 0.5 )
0348     #geom = PhiCut( 0.0, 0.75 )
0349     #geom = PhiCut( 0.0, 1.5 )
0350
0351     modes = [0,1,2,3]
0352     #modes = [2,]
0353
0354     scan = Scan.XY(geom, 100, modes=modes)
0355     Scan.Plot(scan)
0356
0357