eic-opticks/ana/proplib.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 PropLib : Geocache Access
0023 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0024
0025 *PropLib* provides access to optical properties of materials and surfaces
0026 as a function of wavelength.
0027
0028 Material properties::
0029
0030     REFRACTIVE_INDEX
0031     ABSORPTION_LENGTH
0032     SCATTERING_LENGTH
0033     REEMISSION_PROB
0034
0035 Surface properties::
0036
0037     DETECT
0038     ABSORB
0039     REFLECT_SPECULAR
0040     REFLECT_DIFFUSE
0041
0042
0043 Example data shapes::
0044
0045     In [5]: mlib.data.shape
0046     Out[5]: (38, 2, 39, 4)
0047
0048     In [6]: slib.data.shape
0049     Out[6]: (48, 2, 39, 4)
0050
0051     In [7]: blib.data.shape
0052     Out[7]: (123, 4, 2, 39, 4)
0053
0054
0055 """
0056 import os, logging, numpy as np
0057 import sys, codecs
0058 b_ = lambda _:codecs.latin_1_encode(_)[0] if sys.version_info.major > 2 else _
0059 log = logging.getLogger(__name__)
0060
0061
0062 from opticks.ana.base import stamp_
0063 from opticks.ana.dat import Dat
0064 from opticks.ana.nload import np_load
0065 from opticks.ana.key import keydir
0066
0067 KEYDIR = keydir()
0068
0069
0070 class Bnd(object):
0071     def __init__(self, names):
0072         nam = np.zeros( (len(names),4 ), dtype="|S64")
0073         for i,name in enumerate(names):
0074             nam[i] = name.split("/")
0075         pass
0076         self.nam = nam
0077
0078     mats = property(lambda self:np.unique(np.hstack([self.nam[:,0],self.nam[:,3]])))
0079     surs = property(lambda self:np.unique(np.hstack([self.nam[:,1],self.nam[:,2]])))
0080
0081     def oms(self, mat):
0082         """bnd indices of omat"""
0083         return np.where(self.nam[:,0] == mat)[0]
0084
0085     def ims(self, mat):
0086         """bnd indices of imat"""
0087         return np.where(self.nam[:,3] == mat)[0]
0088
0089
0090
0091 class PropLib(object):
0092
0093     # 1st set of 4 [0]
0094     M_REFRACTIVE_INDEX = 0,0
0095     M_ABSORPTION_LENGTH = 0,1
0096     M_SCATTERING_LENGTH = 0,2
0097     M_REEMISSION_PROB = 0,3
0098
0099     L_GROUP_VELOCITY = 1,0
0100
0101     M_LABELS = "RINDEX ABSLEN RAYLEIGH REEMPROB GROUPVEL".split()
0102
0103
0104     # 2nd set of 4 [0] startswith GROUPVEL, currently not used
0105
0106     S_DETECT = 0,0
0107     S_ABSORB = 0,1
0108     S_REFLECT_SPECULAR = 0,2
0109     S_REFLECT_DIFFUSE = 0,3
0110
0111     B_OMAT = 0
0112     B_OSUR = 1
0113     B_ISUR = 2
0114     B_IMAT = 3
0115
0116     COARSE_DOMAIN = np.linspace(60.,820., 39)
0117     FINE_DOMAIN = np.linspace(60., 820., 761)
0118
0119     @classmethod
0120     def load_GBndLib(cls, base):
0121         t, t_paths = np_load(base,"GBndLib/GBndLib.npy")
0122         o, o_paths = np_load(base,"GBndLib/GBndLibOptical.npy")
0123         if t is None or o is None:
0124             log.warning("missing GBndLib data : cannot create blib Proplib")
0125             blib = None
0126         else:
0127             blib = cls("GBndLib", data=t, names=os.path.join(base,"GItemList/GBndLib.txt"), optical=o )
0128             ## hmm GBndLib.txt no longer exists ? presumably due top more dynamic nature of boundaries?
0129         pass
0130         return blib
0131
0132     def __init__(self, kls="GMaterialLib", data=None, names=None, optical=None):
0133         """
0134         :param kls:
0135         :param data:
0136         :param names:
0137
0138         *data* and *names* arguments can take several types:
0139
0140         None
0141             load standard paths
0142         str
0143             load the given path
0144         ndarray
0145              for data, adopt the array
0146         list
0147              for names, adopt the list
0148
0149         """
0150         self.kls = kls
0151         self.paths = []
0152         log.info("names : %s " % names)
0153
0154         if names is None:
0155             npath=os.path.join(KEYDIR, "GItemList/%(kls)s.txt" % locals())
0156         elif type(names) is str:
0157             npath = os.path.expandvars(names)
0158         elif type(names) is np.ndarray or type(names) is list:
0159             npath = None
0160         pass
0161
0162         log.info("npath : %s " % npath)
0163         if npath is None:
0164             log.warning("direct names override")
0165         else:
0166             self.paths.append(npath)
0167             names = list(map(lambda _:_[:-1],open(npath,"r").readlines()))
0168         pass
0169         log.info("names : %s " % names)
0170         self.names = names
0171
0172         if data is None:
0173             dpath=os.path.join(KEYDIR, "%(kls)s/%(kls)s.npy" % locals())
0174         elif type(data) is str:
0175             dpath = data
0176         elif type(data) is np.ndarray:
0177             dpath = None
0178         pass
0179
0180         if dpath is None:
0181             log.warning("direct data override")
0182             data_ = data
0183         else:
0184             dpath = os.path.expandvars(dpath)
0185             self.paths.append(dpath)
0186             data_ = np.load(dpath)
0187         pass
0188
0189         if len(self.names) != data_.shape[0]:
0190             data = data_.reshape(-1,4,2,data_.shape[0],4)
0191             log.warning("reshaped %s from  %r -> %r  " % (dpath, data_.shape, data.shape) )
0192         else:
0193             data = data_
0194         pass
0195         self.data = data
0196         #self.domain = self.COARSE_DOMAIN
0197         self.domain = self.FINE_DOMAIN
0198
0199         assert len(self.names) == self.data.shape[0]
0200         pass
0201
0202         if kls == "GBndLib":
0203             if optical is None:
0204                 opticalpath = idp_("%(kls)s/%(kls)sOptical.npy" % locals())
0205                 self.optical = np.load(opticalpath)
0206                 self.paths.append(opticalpath)
0207             else:
0208                 self.optical = optical
0209             pass
0210             self.dat = Dat(data, names, "omat osur isur imat".split(), "g0 g1".split() )
0211             self.bnd = Bnd(names)
0212         pass
0213
0214     brief = property(lambda self:"%s" % ( self.kls ))
0215
0216
0217     def _set_names(self, names):
0218         """
0219         ndarray has advantage of easy grabbing masks and index lists
0220         """
0221         self._names = np.zeros( (len(names)), dtype="|S128")   # 64 truncates
0222         self._names[:] = names
0223     def _get_names(self):
0224         return self._names
0225     names = property(_get_names, _set_names)
0226
0227
0228     def __repr__(self):
0229         return "\n".join([self.brief] + self.paths + map(stamp_, self.paths))
0230
0231     def index(self, name):
0232         names = list(map(lambda _:_.decode("utf-8"), self.names ))
0233         return names.index(name)
0234         #return np.where( self._names == b_(name) )[0][0]
0235
0236     def interp(self, name, wavelengths, jl):
0237         """
0238         :param name: eg material name "Water"
0239         :param wavelengths: values or array
0240         :param jl: (j,l) tuple
0241
0242                            i  j   k  l
0243         data shape is eg (17, 2, 39, 4)
0244
0245         i: material index 0:16
0246         j: 0:1 property group
0247         k: 0:38 wavelenth samples
0248         l: 0:3 property from the "float4"
0249
0250         """
0251         j,l = jl
0252
0253         #bname = b_(name)
0254         #if not bname in self.names:
0255         #    log.fatal("bname %s is not in list %s " % (bname, self.names))
0256         #pass
0257         #idx = list(self.names).index(bname)
0258         idx = self.index(name)
0259         return np.interp( wavelengths, self.domain, self.data[idx,j,:,l] )
0260
0261     def __call__(self, name):
0262         idx = self.index(name)
0263         return self.data[idx]
0264
0265     def as_mlib(self, order=None):
0266         return self.asmateriallib(self, order=order)
0267
0268     @classmethod
0269     def asmateriallib(cls, blib, order=None):
0270         """
0271         :param blib: boundary lib
0272         :return mlib: material lib
0273
0274         Convert a boundary lib with shape like (123, 4, 2, 39, 4)
0275         into material lib with shape like (38, 2, 39, 4)
0276
0277         Duplicated material property values from
0278         boundary buffer are collapsed into material buffer layout,
0279         checking that all occurences of the material in the boundary buffer
0280         are exact matches.
0281         """
0282         assert blib.kls == "GBndLib"
0283         bsh = blib.data.shape
0284         assert len(bsh) == 5
0285         names = blib.bnd.mats
0286
0287         if not order is None:
0288             order = map(str, order)  # convert from ndarray to simple list
0289             names = sorted(names, key=lambda name:order.index(name) if name in order else len(order))
0290         pass
0291
0292         mdat = np.zeros( ( len(names), bsh[2], bsh[3], bsh[4] ), dtype=np.float32 )
0293
0294         for imat,mat in enumerate(names):
0295             oms = blib.bnd.oms(mat)   # indices of omat
0296             ims = blib.bnd.ims(mat)   # indices of imat
0297
0298             for om in oms:
0299                 bom = blib.data[om,0]
0300
0301                 if np.all(mdat[imat] == 0.):
0302                     mdat[imat] = bom
0303                 else:
0304                     assert np.all(mdat[imat] == bom)
0305                 pass
0306             pass
0307
0308             for im in ims:
0309                 bim = blib.data[im,3]
0310
0311                 if np.all(mdat[imat] == 0.):
0312                     mdat[imat] = bim
0313                 else:
0314                     assert np.all(mdat[imat] == bim)
0315                 pass
0316             pass
0317         pass
0318         return cls("GMaterialLib", data=mdat, names=names)
0319
0320
0321
0322
0323 if __name__ == '__main__':
0324     from opticks.ana.main import opticks_main
0325     ok = opticks_main()
0326
0327     log.info("---")
0328
0329     mlib = PropLib("GMaterialLib")
0330     slib = PropLib("GSurfaceLib")
0331
0332     idx = mlib.index("Water")
0333     print("idx:%s " % idx)
0334
0335     m1 = "Water"
0336     wavelength = 442
0337     ri = mlib.interp(m1, wavelength, mlib.M_REFRACTIVE_INDEX)
0338
0339     print("m1 %s wl %s ri %s " % (m1, wavelength, ri))
0340
0341
0342     # not working , perhaps due to the move to dynamic long ago
0343     #blib = PropLib("GBndLib")
0344     #op = blib.optical
0345     #print("op:%s" % op)
0346
0347