diff options
| -rw-r--r-- | __init__.py | 6 | ||||
| -rw-r--r-- | color/chromaticity.py (renamed from chromaticity.py) | 1 | ||||
| -rw-r--r-- | color/ciexyz64_1.csv (renamed from ciexyz64_1.csv) | 0 | ||||
| -rw-r--r-- | color/sbrgb10w.csv (renamed from sbrgb10w.csv) | 0 | ||||
| -rw-r--r-- | detectors/lbne.py | 72 | ||||
| -rw-r--r-- | geometry.py | 158 | ||||
| -rw-r--r-- | make.py | 14 | ||||
| -rw-r--r-- | materials/OPTICS.ratdb (renamed from OPTICS.ratdb) | 0 | ||||
| -rw-r--r-- | materials/__init__.py (renamed from materials.py) | 43 | ||||
| -rw-r--r-- | materials/ratdb.py (renamed from ratdb.py) | 0 | ||||
| -rw-r--r-- | solid.py | 49 | ||||
| -rw-r--r-- | solids/__init__.py | 1 | ||||
| -rw-r--r-- | solids/r7081.py | 29 | ||||
| -rw-r--r-- | stl.py (renamed from mesh.py) | 31 | ||||
| -rwxr-xr-x | view.py | 19 | ||||
| -rw-r--r-- | viewpmt.py | 30 |
16 files changed, 203 insertions, 250 deletions
diff --git a/__init__.py b/__init__.py index f6d25eb..e69de29 100644 --- a/__init__.py +++ b/__init__.py @@ -1,6 +0,0 @@ -import geometry -import materials -import transform -import stl -import photon -import likelihood diff --git a/chromaticity.py b/color/chromaticity.py index 4080496..6229b51 100644 --- a/chromaticity.py +++ b/color/chromaticity.py @@ -1,6 +1,5 @@ import numpy as np -#f = open('ciexyz64_1.csv') f = open('sbrgb10w.csv') color_map = [] diff --git a/ciexyz64_1.csv b/color/ciexyz64_1.csv index f437db7..f437db7 100644 --- a/ciexyz64_1.csv +++ b/color/ciexyz64_1.csv diff --git a/sbrgb10w.csv b/color/sbrgb10w.csv index b654ad3..b654ad3 100644 --- a/sbrgb10w.csv +++ b/color/sbrgb10w.csv diff --git a/detectors/lbne.py b/detectors/lbne.py index a439682..7f1bb41 100644 --- a/detectors/lbne.py +++ b/detectors/lbne.py @@ -6,13 +6,12 @@ dir = os.path.split(os.path.realpath(__file__))[0] sys.path.append(dir + '/..') import models -from mesh import mesh_from_stl -from solid import Solid -from geometry import Geometry +from geometry import * +from materials import * +from solids import * from transform import rotate, make_rotation_matrix from itertools import product -from materials import * -from make import cylinder +import make endcap_spacing = .485 @@ -27,46 +26,19 @@ class LBNE(Geometry): super(LBNE, self).__init__() # outer cylinder - cylinder_solid = Solid(0, cylinder(radius, height+height/(pmts_per_string-1)), material1=lightwater_sno, material2=vacuum, surface=black_surface, color=0x0000ff) + cylinder_mesh = \ + make.cylinder(radius, height+height/(pmts_per_string-1)) + cylinder_solid = Solid(cylinder_mesh, lightwater_sno, vacuum, black_surface, 0x0000ff) self.add_solid(cylinder_solid) - pmt_inner_mesh = \ - mesh_from_stl(models.dir + '/hamamatsu_12inch_inner.stl') - pmt_outer_mesh = \ - mesh_from_stl(models.dir + '/hamamatsu_12inch_outer.stl') - - photocathode_triangles = np.mean(pmt_inner_mesh[:], axis=1)[:,1] > 0 - - inner_color = np.empty(len(pmt_inner_mesh.triangles), np.uint32) - inner_color[photocathode_triangles] = 0xff0000 - inner_color[~photocathode_triangles] = 0x00ff00 - - outer_color = np.empty(len(pmt_outer_mesh.triangles), np.uint32) - outer_color[:] = 0xffffff - - inner_surface = np.empty(len(pmt_inner_mesh.triangles), np.object) - inner_surface[photocathode_triangles] = black_surface - inner_surface[~photocathode_triangles] = shiny_surface - - self.pmtids = [1] + self.pmtids = [] # construct the barrel for i in range(pmts_per_string): for j in range(nstrings): rotation = make_rotation_matrix(j*2*np.pi/nstrings, (0,0,1)) displacement = rotate((0,-radius,-height/2+i*height/(pmts_per_string-1)), j*2*np.pi/nstrings, (0,0,1)) - self.add_solid(Solid(self.pmtids[-1], pmt_inner_mesh, - rotation, displacement, - material1=vacuum, - material2=glass, - surface=inner_surface, - color=inner_color)) - self.add_solid(Solid(self.pmtids[-1], pmt_outer_mesh, - rotation, displacement, - material1=glass, - material2=lightwater_sno, - color=outer_color)) - self.pmtids += [self.pmtids[-1]+1] + self.pmtids.append(self.add_solid(r7081, rotation, displacement)) # construct the top endcap for x, y in np.array(tuple(product(\ @@ -75,18 +47,7 @@ class LBNE(Geometry): if np.sqrt(x**2 + y**2) <= radius: rotation = make_rotation_matrix(+np.pi/2, (1,0,0)) displacement = (x,y,+height/2+height/(pmts_per_string-1)/2) - self.add_solid(Solid(self.pmtids[-1], pmt_inner_mesh, - rotation, displacement, - material1=vacuum, - material2=glass, - surface=inner_surface, - color=inner_color)) - self.add_solid(Solid(self.pmtids[-1], pmt_outer_mesh, - rotation, displacement, - material1=glass, - material2=lightwater_sno, - color=outer_color)) - self.pmtids += [self.pmtids[-1]+1] + self.pmtids.append(self.add_solid(r7081, rotation, displacement)) # construct the bottom endcap for x, y in np.array(tuple(product(\ @@ -95,15 +56,4 @@ class LBNE(Geometry): if np.sqrt(x**2 + y**2) <= radius: rotation = make_rotation_matrix(-np.pi/2, (1,0,0)) displacement = (x,y,-height/2-height/(pmts_per_string-1)/2) - self.add_solid(Solid(self.pmtids[-1], pmt_inner_mesh, - rotation, displacement, - material1=vacuum, - material2=glass, - surface=inner_surface, - color=inner_color)) - self.add_solid(Solid(self.pmtids[-1], pmt_outer_mesh, - rotation, displacement, - material1=glass, - material2=lightwater_sno, - color=outer_color)) - self.pmtids += [self.pmtids[-1]+1] + self.pmtids.append(self.add_solid(r7081, rotation, displacement)) diff --git a/geometry.py b/geometry.py index 5cd1bf2..be31e4b 100644 --- a/geometry.py +++ b/geometry.py @@ -1,9 +1,117 @@ import numpy as np -import numpy.ma as ma import pycuda.driver as cuda from pycuda import gpuarray -from mesh import Mesh -from materials import standard_wavelengths + +# all material/surface properties are interpolated at these +# wavelengths when they are sent to the gpu +standard_wavelengths = np.arange(200, 810, 20).astype(np.float32) + +class Mesh(object): + def __init__(self, vertices, triangles): + vertices = np.asarray(vertices, dtype=np.float32) + triangles = np.asarray(triangles, dtype=np.int32) + + if len(vertices.shape) != 2 or vertices.shape[1] != 3: + raise ValueError('shape mismatch') + + if len(triangles.shape) != 2 or triangles.shape[1] != 3: + raise ValueError('shape mismatch') + + if (triangles < 0).any(): + raise ValueError('indices in `triangles` must be positive.') + + if (triangles >= len(vertices)).any(): + raise ValueError('indices in `triangles` must be less than the ' + 'length of the vertex array.') + + self.vertices = vertices + self.triangles = triangles + + def build(self): + return self.vertices[self.triangles] + + def __getitem__(self, key): + return self.vertices[self.triangles[key]] + + def __len__(self): + return len(self.triangles) + + def __add__(self, other): + return Mesh(np.concatenate((self.vertices, other.vertices)), np.concatenate((self.triangles, other.triangles + len(self.vertices)))) + +class Solid(object): + def __init__(self, mesh, material1=None, material2=None, surface=None, color=0xffffffff): + self.mesh = mesh + + if np.iterable(material1): + if len(material1) != len(mesh): + raise ValueError('shape mismatch') + self.material1 = np.array(material1, dtype=np.object) + else: + self.material1 = np.tile(material1, len(self.mesh)) + + if np.iterable(material2): + if len(material2) != len(mesh): + raise ValueError('shape mismatch') + self.material2 = np.array(material2, dtype=np.object) + else: + self.material2 = np.tile(material2, len(self.mesh)) + + if np.iterable(surface): + if len(surface) != len(mesh): + raise ValueError('shape mismatch') + self.surface = np.array(surface, dtype=np.object) + else: + self.surface = np.tile(surface, len(self.mesh)) + + if np.iterable(color): + if len(color) != len(mesh): + raise ValueError('shape mismatch') + self.color = np.array(color, dtype=np.uint32) + else: + self.color = np.tile(color, len(self.mesh)).astype(np.uint32) + + def __len__(self): + return len(self.mesh) + + def __add__(self, other): + return Solid(self.mesh + other.mesh, np.concatenate((self.material1, other.material1)), np.concatenate((self.material2, other.material2)), np.concatenate((self.surface, other.surface)), np.concatenate((self.color, other.color))) + +class Material(object): + """Material optical properties.""" + def __init__(self, name='none'): + self.name = name + + self.refractive_index = None + self.absorption_length = None + self.scattering_length = None + + def set(self, name, value, wavelengths=standard_wavelengths): + if np.iterable(value): + if len(value) != len(wavelengths): + raise ValueError('shape mismatch') + else: + value = np.tile(value, len(wavelengths)) + + self.__dict__[name] = np.array(zip(wavelengths, value), dtype=np.float32) + +class Surface(object): + """Surface optical properties.""" + def __init__(self, name='none'): + self.name = name + + self.absorption = None + self.reflection_diffuse = None + self.reflection_specular = None + + def set(self, name, value, wavelengths=standard_wavelengths): + if np.iterable(value): + if len(value) != len(wavelengths): + raise ValueError('shape mismatch') + else: + value = np.tile(value, len(wavelengths)) + + self.__dict__[name] = np.array(zip(wavelengths, value), dtype=np.float32) def interleave(arr, bits): """ @@ -51,15 +159,30 @@ def morton_order(mesh, bits): return interleave(mean_positions, bits) class Geometry(object): - def __init__(self, solids=None): - if solids is not None: - self.solids = list(solids) - else: - self.solids = [] + def __init__(self): + self.solids = [] + self.solid_rotations = [] + self.solid_displacements = [] + + def add_solid(self, solid, rotation=np.identity(3), displacement=(0,0,0)): + rotation = np.asarray(rotation, dtype=np.float32) + + if rotation.shape != (3,3): + raise ValueError('shape mismatch') + + self.solid_rotations.append(rotation.astype(np.float32)) + + displacement = np.asarray(displacement, dtype=np.float32) + + if displacement.shape != (3,): + raise ValueError('shape mismatch') + + self.solid_displacements.append(displacement) - def add_solid(self, solid): self.solids.append(solid) + return len(self.solids)-1 + def build(self, bits=8): offsets = [ (0,0) ] for solid in self.solids: @@ -69,15 +192,14 @@ class Geometry(object): vertices = np.zeros(shape=(offsets[-1][0], 3), dtype=np.float32) triangles = np.zeros(shape=(offsets[-1][1],3), dtype=np.int32) - for solid, (vertex_offset, triangle_offset) in zip(self.solids, offsets[:-1]): + for i, (solid, (vertex_offset, triangle_offset)) in \ + enumerate(zip(self.solids, offsets[:-1])): triangles[triangle_offset:triangle_offset+len(solid.mesh.triangles),:] = \ solid.mesh.triangles + vertex_offset - vertices[vertex_offset:vertex_offset + len(solid.mesh.vertices),:] = \ - np.inner(solid.mesh.vertices, solid.rotation) + solid.displacement + vertices[vertex_offset:vertex_offset+len(solid.mesh.vertices),:] = \ + np.inner(solid.mesh.vertices, self.solid_rotations[i]) + self.solid_displacements[i] self.mesh = Mesh(vertices, triangles) - #del vertices - #del triangles zvalues_mesh = morton_order(self.mesh[:], bits) reorder = np.argsort(zvalues_mesh) @@ -121,10 +243,12 @@ class Geometry(object): else: self.surface_index[i] = -1 - self.colors = np.concatenate([solid.color for solid in self.solids])[reorder] + self.colors = \ + np.concatenate([solid.color for solid in self.solids])[reorder] - self.solid_id = np.concatenate([np.tile(solid.id, len(solid.mesh)) \ - for solid in self.solids])[reorder] + self.solid_id = \ + np.concatenate([np.tile(i, len(solid.mesh)) for i, solid in \ + enumerate(self.solids)])[reorder] unique_zvalues = np.unique(zvalues_mesh) zvalues = np.empty(unique_zvalues.size, dtype=np.uint64) @@ -1,5 +1,5 @@ import numpy as np -from mesh import Mesh +from geometry import Mesh def cylinder(radius=1.0, height=1.0, theta=np.pi/32): angles = np.arange(0, 2*np.pi, theta) @@ -41,15 +41,3 @@ def cylinder(radius=1.0, height=1.0, theta=np.pi/32): np.roll(np.arange(len(angles)),1) return Mesh(vertices, triangles) - -if __name__ == '__main__': - from solid import Solid - from geometry import Geometry - from view import view - - cylinder_mesh = cylinder() - cylinder_solid = Solid(0, cylinder_mesh) - geometry = Geometry([cylinder_solid]) - geometry.build(bits=8) - - view(geometry) diff --git a/OPTICS.ratdb b/materials/OPTICS.ratdb index c7551a5..c7551a5 100644 --- a/OPTICS.ratdb +++ b/materials/OPTICS.ratdb diff --git a/materials.py b/materials/__init__.py index 9feb866..d2644ad 100644 --- a/materials.py +++ b/materials/__init__.py @@ -1,43 +1,12 @@ +import os +import sys import numpy as np from ratdb import load -standard_wavelengths = np.arange(200, 810, 20).astype(np.float32) +dir = os.path.split(os.path.realpath(__file__))[0] +sys.path.append(dir + '/..') -class Material(object): - """Material optical properties.""" - def __init__(self, name='none'): - self.name = name - - self.refractive_index = None - self.absorption_length = None - self.scattering_length = None - - def set(self, name, value, wavelengths=standard_wavelengths): - if np.iterable(value): - if len(value) != len(wavelengths): - raise ValueError('shape mismatch') - else: - value = np.tile(value, len(wavelengths)) - - self.__dict__[name] = np.array(zip(wavelengths, value), dtype=np.float32) - -class Surface(object): - """Surface optical properties.""" - def __init__(self, name='none'): - self.name = name - - self.absorption = None - self.reflection_diffuse = None - self.reflection_specular = None - - def set(self, name, value, wavelengths=standard_wavelengths): - if np.iterable(value): - if len(value) != len(wavelengths): - raise ValueError('shape mismatch') - else: - value = np.tile(value, len(wavelengths)) - - self.__dict__[name] = np.array(zip(wavelengths, value), dtype=np.float32) +from geometry import Material, Surface vacuum = Material('vacuum') vacuum.set('refractive_index', 1) @@ -59,7 +28,7 @@ shiny_surface.set('absorption', 0) shiny_surface.set('reflection_specular', 1) shiny_surface.set('reflection_diffuse', 0) -db = load(open('OPTICS.ratdb'))['OPTICS'] +db = load(open(dir + '/OPTICS.ratdb'))['OPTICS'] glass = Material('pmt_glass') glass.set('refractive_index', db['glass']['RINDEX_value2'], db['glass']['RINDEX_value1']) diff --git a/ratdb.py b/materials/ratdb.py index 2b14761..2b14761 100644 --- a/ratdb.py +++ b/materials/ratdb.py diff --git a/solid.py b/solid.py deleted file mode 100644 index 707a786..0000000 --- a/solid.py +++ /dev/null @@ -1,49 +0,0 @@ -import numpy as np - -class Solid(object): - def __init__(self, id, mesh, rotation=np.identity(3), displacement=(0,0,0), material1=None, material2=None, surface=None, color=0xffffffff): - self.id = id - self.mesh = mesh - - if rotation.shape != (3,3): - raise ValueError('shape mismatch') - - self.rotation = rotation.astype(np.float32) - - displacement = np.asarray(displacement, dtype=np.float32) - - if displacement.shape != (3,): - raise ValueError('shape mismatch') - - self.displacement = displacement - - if np.iterable(material1): - if len(material1) != len(mesh): - raise ValueError('shape mismatch') - self.material1 = np.array(material1, dtype=np.object) - else: - self.material1 = np.tile(material1, len(self.mesh)) - - if np.iterable(material2): - if len(material2) != len(mesh): - raise ValueError('shape mismatch') - self.material2 = np.array(material2, dtype=np.object) - else: - self.material2 = np.tile(material2, len(self.mesh)) - - if np.iterable(surface): - if len(surface) != len(mesh): - raise ValueError('shape mismatch') - self.surface = np.array(surface, dtype=np.object) - else: - self.surface = np.tile(surface, len(self.mesh)) - - if np.iterable(color): - if len(color) != len(mesh): - raise ValueError('shape mismatch') - self.color = np.array(color, dtype=np.uint32) - else: - self.color = np.tile(color, len(self.mesh)).astype(np.uint32) - - def __len__(self): - return len(self.mesh) diff --git a/solids/__init__.py b/solids/__init__.py new file mode 100644 index 0000000..f01d97f --- /dev/null +++ b/solids/__init__.py @@ -0,0 +1 @@ +from r7081 import r7081 diff --git a/solids/r7081.py b/solids/r7081.py new file mode 100644 index 0000000..f91c2ab --- /dev/null +++ b/solids/r7081.py @@ -0,0 +1,29 @@ +import os +import sys +import numpy as np + +dir = os.path.split(os.path.realpath(__file__))[0] +sys.path.append(dir + '/..') + +import models +from mesh import mesh_from_stl +from geometry import * +from materials import * + +r7081_outer_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_outer.stl') +r7081_inner_mesh = mesh_from_stl(models.dir + '/hamamatsu_12inch_inner.stl') + +photocathode_triangles = np.mean(r7081_inner_mesh[:], axis=1)[:,1] > 0 + +inner_color = np.empty(len(r7081_inner_mesh.triangles), np.uint32) +inner_color[photocathode_triangles] = 0xff0000 +inner_color[~photocathode_triangles] = 0x00ff00 + +inner_surface = np.empty(len(r7081_inner_mesh.triangles), np.object) +inner_surface[photocathode_triangles] = black_surface +inner_surface[~photocathode_triangles] = shiny_surface + +r7081_inner_solid = Solid(r7081_inner_mesh, vacuum, glass, inner_surface, color=inner_color) +r7081_outer_solid = Solid(r7081_outer_mesh, glass, lightwater_sno) + +r7081 = r7081_inner_solid + r7081_outer_solid @@ -1,36 +1,7 @@ import numpy as np import string import struct - -class Mesh(object): - def __init__(self, vertices, triangles): - vertices = np.asarray(vertices, dtype=np.float32) - triangles = np.asarray(triangles, dtype=np.int32) - - if len(vertices.shape) != 2 or vertices.shape[1] != 3: - raise ValueError('shape mismatch') - - if len(triangles.shape) != 2 or triangles.shape[1] != 3: - raise ValueError('shape mismatch') - - if (triangles < 0).any(): - raise ValueError('indices in `triangles` must be positive.') - - if (triangles >= len(vertices)).any(): - raise ValueError('indices in `triangles` must be less than the ' - 'length of the vertex array.') - - self.vertices = vertices - self.triangles = triangles - - def build(self): - return self.vertices[self.triangles] - - def __getitem__(self, key): - return self.vertices[self.triangles[key]] - - def __len__(self): - return len(self.triangles) +from geometry import Mesh def mesh_from_stl(filename): f = open(filename) @@ -6,7 +6,6 @@ from pygame.locals import * import src from camera import * -from solid import * from geometry import * from transform import * @@ -182,6 +181,8 @@ if __name__ == '__main__': import optparse from mesh import * + + import solids import detectors parser = optparse.OptionParser('%prog filename.stl') @@ -197,20 +198,26 @@ if __name__ == '__main__': if ext.lower() == '.stl': geometry = Geometry() - geometry.add_solid(Solid(0, mesh_from_stl(args[0]))) + geometry.add_solid(Solid(mesh_from_stl(args[0]))) geometry.build(options.bits) view(geometry, tail) else: import inspect - members = dict(inspect.getmembers(detectors)) + members = dict(inspect.getmembers(detectors) + inspect.getmembers(solids)) if args[0] in members: - if issubclass(members[args[0]], Geometry): + if type(members[args[0]]) is type and \ + issubclass(members[args[0]], Geometry): geometry = members[args[0]]() geometry.build(options.bits) view(geometry, args[0]) + elif isinstance(members[args[0]], Solid): + geometry = Geometry() + geometry.add_solid(members[args[0]]) + geometry.build(options.bits) + view(geometry, args[0]) else: - sys.exit("%s doesn't inherit from a Geometry object" % args[0]) + sys.exit("can't display %s" % args[0]) else: - sys.exit("couldn't find a geometry with name %s" % args[0]) + sys.exit("couldn't find object %s" % args[0]) diff --git a/viewpmt.py b/viewpmt.py deleted file mode 100644 index 474b6d8..0000000 --- a/viewpmt.py +++ /dev/null @@ -1,30 +0,0 @@ -import numpy as np -from solid import Solid -from mesh import mesh_from_stl -from geometry import Geometry -from view import view -import models - -pmt_inner_mesh = \ - mesh_from_stl(models.dir + '/hamamatsu_12inch_inner.stl') -pmt_outer_mesh = \ - mesh_from_stl(models.dir + '/hamamatsu_12inch_outer.stl') - -pmt_outer_mesh.triangles = \ - pmt_outer_mesh.triangles[np.mean(pmt_outer_mesh[:], axis=1)[:,0] > 0] - -photocathode_triangles = np.mean(pmt_inner_mesh[:], axis=1)[:,1] > 0 - -inner_color = np.empty(len(pmt_inner_mesh.triangles), np.uint32) -inner_color[photocathode_triangles] = 0xff0000 -inner_color[~photocathode_triangles] = 0x00ff00 - -outer_color = np.empty(len(pmt_outer_mesh.triangles), np.uint32) -outer_color[:] = 0xffffff - -geometry = Geometry([Solid(0, pmt_inner_mesh, color=inner_color), - Solid(1, pmt_outer_mesh, color=outer_color)]) - -geometry.build(bits=8) - -view(geometry) |