eic-opticks/ana/geometry.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 geometry.py : Shape, Ray, Plane, Intersect, IntersectFrame
0023 ============================================================
0024
0025 Simple geometry intersection calculations.
0026
0027
0028
0029 """
0030 import os, logging
0031 import numpy as np
0032 from opticks.ana.boundary import Boundary
0033 from StringIO import StringIO
0034
0035 log = logging.getLogger(__name__)
0036
0037
0038 def mat4_tostring(m):
0039     from StringIO import StringIO
0040     sio = StringIO()
0041     np.savetxt(sio,m,delimiter=',',fmt="%.3f")
0042     s = sio.getvalue().replace("\n",",")
0043     return s
0044
0045 def mat4_fromstring(s):
0046     m = np.fromstring(s,sep=",").reshape(4,4)
0047     return m
0048
0049 class Shape(object):
0050     def __init__(self, parameters, boundary):
0051          self.parameters = np.fromstring(parameters, sep=",")
0052          if type(boundary) is str:
0053              self.boundary = Boundary(boundary)
0054          else:
0055              self.boundary = boundary
0056          pass
0057
0058     def refractive_index(self, wavelength):
0059          return self.boundary.imat.refractive_index(wavelength)
0060
0061     def __repr__(self):
0062          return "%s %s %s " % ( self.__class__.__name__ , self.parameters, self.boundary )
0063
0064
0065 class Ray(object):
0066     def __init__(self, o, d):
0067         self.origin = np.array(o)
0068         self.direction = np.array(d)
0069
0070     def position(self, t):
0071         return self.origin + t*self.direction
0072
0073
0074 class Plane(object):
0075     def __init__(self, n, p):
0076         """
0077         :param n: normal
0078         :param p: point in the plane
0079         """
0080         n = np.array(n)
0081         p = np.array(p)
0082         n = n/np.linalg.norm(n,2)
0083         d = -np.dot(n,p)
0084         pass
0085         self.n = n
0086         self.p = p
0087         self.d = d
0088
0089     def ntile(self, num):
0090         return np.tile( self.n, num).reshape(-1,3)
0091
0092     def intersect(self, ray):
0093         denom_ = np.dot( self.n, ray.direction )
0094         return denom_, -(self.d + np.dot(self.n, ray.origin))/denom_
0095
0096
0097 class Intersect(object):
0098     def __init__(self, i, t, n, p):
0099         self.i = i
0100         self.t = t
0101         self.n = n
0102         self.p = p
0103
0104     def __repr__(self):
0105         return "i %2d t %10.4f n %25s p %25s  " % (self.i, self.t, self.n,self.p)
0106
0107
0108 class IntersectFrame(object):
0109     """
0110     world_to_intersect 4x4 homogenous matrix
0111     intersect frame basis
0112
0113     Arrange frame with origin at
0114     intersection point and with basis vectors
0115
0116     U: in plane of face
0117     V: normal to the intersected plane
0118     W: perpendicular to both the above
0119
0120                  V
0121                  |
0122                  |
0123           -------*---U------A-
0124                 W           another
0125     """
0126     def __init__(self, isect, a):
0127
0128         p = isect.p
0129         n = isect.n
0130
0131         norm_ = lambda v:v/np.linalg.norm(v)
0132
0133         U = norm_( a - p )
0134         V = n
0135         W = np.cross(U,V)
0136
0137         ro = np.identity(4)
0138         ro[:3,0] = U
0139         ro[:3,1] = V
0140         ro[:3,2] = W
0141
0142         tr = np.identity(4)
0143         tr[:3,3] = p
0144
0145         itr = np.identity(4)
0146         itr[:3,3] = -p
0147
0148         self.w2i = np.dot(ro.T,itr)
0149         self.i2w = np.dot(tr, ro)
0150
0151         self.p = p
0152         self.n = n
0153         self.a = a
0154
0155
0156
0157     def homogenize(self, v, w=1):
0158         assert len(v) == 3
0159         vh = np.zeros(4)
0160         vh[:3] = v
0161         vh[3] = w
0162         return vh
0163
0164     def world_to_intersect(self, v, w=1):
0165         wfp = self.homogenize(v, w)
0166         ifp = np.dot( self.w2i, wfp )
0167         log.info("world_to_intersect  wfp %s -> ifp %s " % (wfp, ifp))
0168         return ifp
0169
0170     def intersect_to_world(self, v, w=1):
0171         ifp = self.homogenize(v, w)
0172         wfp = np.dot( self.i2w, ifp )
0173         log.info("intersect_to_world  ifp %s -> wfp %s " % (ifp, wfp))
0174         return wfp
0175
0176     def i2w_string(self):
0177         return mat4_tostring(self.i2w.T)
0178
0179     def id_string(self):
0180         return mat4_tostring(np.identity(4))
0181
0182
0183
0184
0185 if __name__ == '__main__':
0186     pass
0187     logging.basicConfig(level=logging.INFO)
0188
0189
0190
0191