moved class definitions for Solid, Mesh, Material, and Surface into geometry.py and moved instances of these classes into separate folders. the Solid object no longer contains a rotation, displacement, or id variable; instead, they are passed to a geometry object when calling add_solid().

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)
diff --git a/make.py b/make.py
index 72b8406..e69b28e 100644
--- a/make.py
+++ b/make.py
@@ -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
diff --git a/mesh.py b/stl.py
index 0f1e8b6..9bd3a38 100644
--- a/mesh.py
+++ b/stl.py
@@ -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)
diff --git a/view.py b/view.py
index 8a61513..da33b75 100755
--- a/view.py
+++ b/view.py
@@ -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)