From d0825a136ff65b36069ff8b078b9fd97adeed0df Mon Sep 17 00:00:00 2001
From: Anthony LaTorre <telatorre@gmail.com>
Date: Wed, 1 Jun 2011 17:00:37 -0400
Subject: first step towards moving to a new mesh/solid/geometry structure

---
 detectors/lbne.py |   5 ++-
 geometry.py       | 126 +++++++++++++++++-------------------------------------
 mesh.py           | 113 ++++++++++++++++++++++++++++++++++++++++++++++++
 solid.py          |  58 +++++++++++++++++++++++++
 src/intersect.h   |   6 +--
 stl.py            |  43 -------------------
 transform.py      |  18 +++++---
 view.py           |  12 ++++--
 8 files changed, 239 insertions(+), 142 deletions(-)
 create mode 100644 mesh.py
 create mode 100644 solid.py
 delete mode 100644 stl.py

diff --git a/detectors/lbne.py b/detectors/lbne.py
index a9136dd..81f5dd7 100644
--- a/detectors/lbne.py
+++ b/detectors/lbne.py
@@ -10,9 +10,10 @@ dir = os.path.split(os.path.realpath(__file__))[0]
 sys.path.append(dir + '/..')
 
 import models
-from stl import read_stl
+from mesh import *
 from transform import rotate
-from geometry import Geometry, Solid
+from solid import Solid
+from geometry import Geometry
 from materials import glass, h2o
 
 endcap_spacing = .485
diff --git a/geometry.py b/geometry.py
index f50fe89..d1bec2a 100644
--- a/geometry.py
+++ b/geometry.py
@@ -48,103 +48,57 @@ def morton_order(mesh, bits):
 
     return interleave(mean_positions, bits)
 
-class Solid(object):
-    """
-    Object which stores a closed triangle mesh associated with a physically
-    distinct object.
-
-    Args:
-        - mesh, array
-            A closed triangle mesh.
-        - material1, Material
-            The material inside within the mesh.
-        - material2, Material
-            The material outside the mesh.
-        - surface1, Surface,
-            The surface on the inside of the mesh.
-        - surface2, Surface,
-            The surface on the outside of the mesh.
-
-    .. warning::
-        It is possible to define logically inconsistent geometries unless
-        you are careful. For example, solid A may define its inside material
-        as water but contain solid B which defines its outside material as air.
-        In this case, a photon traveling out of solid B will reflect/refract
-        assuming it's going into air, but upon reaching solid A's boundary will
-        calculate attenuation factors assuming it just traveled through water.
-    """
-    def __init__(self, mesh, material1, material2, \
-                     surface1=None, surface2=None):
-        if len(mesh.shape) != 3 or mesh.shape[1] != 3 or mesh.shape[2] != 3:
-            raise Exception('shape mismatch; mesh must be a triangle mesh')
+class Geometry(object):
+    def __init__(self, solids=[]):
+        self.solids = solids
 
-        self.mesh = mesh
-        self.material1 = material1
-        self.material2 = material2
-        self.surface1 = surface1
-        self.surface2 = surface2
+    def add_solid(self, solid):
+        self.solids.append(solid)
 
-    def __len__(self):
-        return self.mesh.shape[0]
+    def build(self, bits=8):
+        self.mesh = np.concatenate([solid.build() for solid in self.solids])
 
-class Geometry(object):
-    """Object which stores the global mesh for a geometry."""
+        zvalues_mesh = morton_order(self.mesh, bits)
+        reorder = np.argsort(zvalues_mesh)
+        zvalues_mesh = zvalues_mesh[reorder]
 
-    def __init__(self):
-        self.solids = []
-        self.materials = []
-        self.surfaces = []
+        if (np.diff(zvalues_mesh) < 0).any():
+            raise Exception('zvalues_mesh out of order.')
 
-    def add_solid(self, solid):
-        """Add a solid to the geometry."""
-        self.solids.append(solid)
+        self.mesh = self.mesh[reorder]
 
-        if solid.material1 not in self.materials:
-            self.materials.append(solid.material1)
+        material1 = np.concatenate([solid.material1 for solid in self.solids])
+        material2 = np.concatenate([solid.material2 for solid in self.solids])
+        
+        self.materials = \
+            list(np.unique(np.concatenate([material1, material2])))
 
-        if solid.material2 not in self.materials:
-            self.materials.append(solid.material2)
+        self.material1_index = np.empty(len(self.mesh), dtype=np.uint32)
+        self.material2_index = np.empty(len(self.mesh), dtype=np.uint32)
 
-        if solid.surface1 not in self.surfaces:
-            self.surfaces.append(solid.surface1)
+        for i, material in enumerate(material1[reorder]):
+            self.material1_index[i] = self.materials.index(material)
 
-        if solid.surface2 not in self.surfaces:
-            self.surfaces.append(solid.surface1)
+        for i, material in enumerate(material2[reorder]):
+            self.material2_index[i] = self.materials.index(material)
 
-    def build(self, bits=8):
-        """Build the bounding volume hierarchy of the geometry."""
-        mesh = np.concatenate([solid.mesh for solid in self.solids])
-
-        # lookup solid/material/surface index from triangle index
-        solid_index = \
-            np.concatenate([np.tile(self.solids.index(solid), \
-                                        len(solid)) for solid in self.solids])
-        material1_index = \
-            np.concatenate([np.tile(self.materials.index(solid.material1), \
-                                        len(solid)) for solid in self.solids])
-        material2_index = \
-            np.concatenate([np.tile(self.materials.index(solid.material2), \
-                                        len(solid)) for solid in self.solids])
-        surface1_index = \
-            np.concatenate([np.tile(self.surfaces.index(solid.surface1), \
-                                        len(solid)) for solid in self.solids])
-        surface2_index = \
-            np.concatenate([np.tile(self.surfaces.index(solid.surface2), \
-                                        len(solid)) for solid in self.solids])
-
-        zvalues_mesh = morton_order(mesh, bits)
-        reorder = np.argsort(zvalues_mesh)
-        zvalues_mesh = zvalues_mesh[reorder]
+        surface1 = np.concatenate([solid.surface1 for solid in self.solids])
+        surface2 = np.concatenate([solid.surface2 for solid in self.solids])
 
-        if (np.diff(zvalues_mesh) < 0).any():
-            raise Exception('zvalues_mesh out of order')
-
-        self.mesh = mesh[reorder]
-        self.solid_index = solid_index[reorder]
-        self.material1_index = material1_index[reorder]
-        self.material2_index = material2_index[reorder]
-        self.surface1_index = surface1_index[reorder]
-        self.surface2_index = surface2_index[reorder]
+        self.surfaces = list(np.unique(np.concatenate([surface1, surface2])))
+
+        self.surface1_index = np.empty(len(self.mesh), dtype=np.uint32)
+        self.surface2_index = np.empty(len(self.mesh), dtype=np.uint32)
+
+        for i, surface in enumerate(surface1[reorder]):
+            self.surface1_index[i] = self.surfaces.index(surface)
+
+        for i, surface in enumerate(surface2[reorder]):
+            self.surface2_index[i] = self.surfaces.index(surface)
+
+        self.solid_index = \
+            np.concatenate([np.tile(i, len(solid.mesh)) \
+                                for i in range(len(self.solids))])[reorder]
 
         unique_zvalues = np.unique(zvalues_mesh)
         zvalues = np.empty(unique_zvalues.size, dtype=np.uint64)
diff --git a/mesh.py b/mesh.py
new file mode 100644
index 0000000..a4b9048
--- /dev/null
+++ b/mesh.py
@@ -0,0 +1,113 @@
+import numpy as np
+import string
+import struct
+
+class Mesh(object):
+    def __init__(self, vertices, triangles):
+        vertices = np.asarray(vertices, dtype=float)
+        triangles = np.asarray(triangles, dtype=int)
+
+        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 __add__(self, other):
+        vertices = np.concatenate((self.vertices, other.vertices))
+        triangles = np.concatenate((self.triangles, other.triangles+len(self.vertices)))
+
+        return Mesh(vertices, triangles)
+
+    def __len__(self):
+        return len(self.triangles)
+
+def mesh_from_stl(filename):
+    f = open(filename)
+    buf = f.read(200)
+    f.close()
+
+    for char in buf:
+        if char not in string.printable:
+            return mesh_from_binary_stl(filename)
+
+    return mesh_from_ascii_stl(filename)
+
+def mesh_from_ascii_stl(filename):
+    f = open(filename)
+
+    vertices = []
+    triangles = []
+    vertex_map = {}
+
+    while True:
+        line = f.readline()
+
+        if line == '':
+            break
+
+        if not line.strip().startswith('vertex'):
+            continue
+
+        triangle = [None]*3
+        for i in range(3):
+            vertex = tuple([float(s) for s in line.strip().split()[1:]])
+
+            if vertex not in vertex_map:
+                vertices.append(vertex)
+                vertex_map[vertex] = len(vertices) - 1
+
+            triangle[i] = vertex_map[vertex]
+
+            if i < 3:
+                line = f.readline()
+
+        triangles.append(triangle)
+
+    f.close()
+
+    return Mesh(np.array(vertices), np.array(triangles, dtype=np.uint32))
+
+def mesh_from_binary_stl(filename):
+    f = open(filename)
+
+    vertices = []
+    triangles = []
+    vertex_map = {}
+
+    f.read(80)
+    ntriangles = struct.unpack('<I', f.read(4))[0]
+
+    for i in range(ntriangles):
+        f.read(12)
+
+        triangle = [None]*3
+        for j in range(3):
+            vertex = tuple(struct.unpack('<fff', f.read(12)))
+
+            if vertex not in vertex_map:
+                vertices.append(vertex)
+                vertex_map[vertex] = len(vertices) - 1
+
+            triangle[j] = vertex_map[vertex]
+
+        triangles.append(triangle)
+
+        f.read(2)
+
+    f.close()
+
+    return Mesh(np.array(vertices), np.array(triangles, dtype=np.uint32))
diff --git a/solid.py b/solid.py
new file mode 100644
index 0000000..fe22463
--- /dev/null
+++ b/solid.py
@@ -0,0 +1,58 @@
+import numpy as np
+
+class Solid(object):
+    def __init__(self, id, mesh, rotation=np.identity(3), displacement=(0,0,0),
+                 material1=None, material2=None, surface1=None, surface2=None,
+                 color=None):
+        self.id = id
+        self.mesh = mesh
+
+        if rotation.shape != (3,3):
+            raise ValueError('shape mismatch')
+
+        self.rotation = rotation
+
+        displacement = np.asarray(displacement)
+
+        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(surface1):
+            if len(surface1) != len(mesh):
+                raise ValueError('shape mismatch')
+            self.surface1 = np.array(surface1, dtype=np.object)
+        else:
+            self.surface1 = np.tile(surface1, len(self.mesh))
+
+        if np.iterable(surface2):
+            if len(surface2) != len(mesh):
+                raise ValueError('shape mismatch')
+            self.surface2 = np.array(surface2, dtype=np.object)
+        else:
+            self.surface2 = np.tile(surface2, 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))
+
+    def build(self):
+        return np.inner(self.mesh.build(), self.rotation) + self.displacement
diff --git a/src/intersect.h b/src/intersect.h
index b984612..a968f09 100644
--- a/src/intersect.h
+++ b/src/intersect.h
@@ -26,21 +26,21 @@ __device__ bool intersect_triangle(const float3 &origin, const float3 &direction
 		    (m.a02*m.a21 - m.a01*m.a22)*b.y +
 		    (m.a01*m.a12 - m.a02*m.a11)*b.z)/determinant;
 
-	if (u1 < 0.0f)
+	if (u1 < 0.0f || u1 > 1.0f)
 		return false;
 
 	float u2 = ((m.a12*m.a20 - m.a10*m.a22)*b.x +
 		    (m.a00*m.a22 - m.a02*m.a20)*b.y +
 		    (m.a02*m.a10 - m.a00*m.a12)*b.z)/determinant;
 
-	if (u2 < 0.0f)
+	if (u2 < 0.0f || u2 > 1.0f)
 		return false;
 
 	float u3 = ((m.a10*m.a21 - m.a11*m.a20)*b.x +
 		    (m.a01*m.a20 - m.a00*m.a21)*b.y +
 		    (m.a00*m.a11 - m.a01*m.a10)*b.z)/determinant;
 
-	if (u3 < 0.0f || (1-u1-u2) < 0.0f)
+	if (u3 < 0.0f || (1.0f-u1-u2) < 0.0f)
 		return false;
 
 	distance = u3;
diff --git a/stl.py b/stl.py
deleted file mode 100644
index a458708..0000000
--- a/stl.py
+++ /dev/null
@@ -1,43 +0,0 @@
-import numpy as np
-import string
-import struct
-
-def read_stl(filename):
-    """Return a triangle mesh from `filename`."""
-    f = open(filename)
-    buf = f.read(200)
-    f.close()
-
-    for char in buf:
-        if char not in string.printable:
-            return read_binary_stl(filename)
-
-    return read_ascii_stl(filename)
-
-def read_ascii_stl(filename):
-    f = open(filename)
-
-    vertex = []
-    for line in f:
-        if not line.strip().startswith('vertex'):
-            continue
-        vertex.append([float(s) for s in line.strip().split()[1:]])
-
-    f.close()
-    return np.array(vertex).reshape(len(vertex)//3,3,3)
-
-def read_binary_stl(filename):
-    f = open(filename)
-
-    f.read(80)
-    triangles = struct.unpack('<I', f.read(4))[0]
-
-    vertex = []
-    for i in range(triangles):
-        f.read(12)
-        for j in range(3):
-            vertex.append(struct.unpack('<fff', f.read(12)))
-        f.read(2)
-
-    f.close()
-    return np.array(vertex).reshape(len(vertex)//3,3,3)
diff --git a/transform.py b/transform.py
index 1321c8c..ef8c96c 100644
--- a/transform.py
+++ b/transform.py
@@ -1,16 +1,24 @@
 import numpy as np
 
+def make_rotation_matrix(phi, n):
+    """
+    Make the rotation matrix to rotate points through an angle `phi`
+    counter-clockwise around the axis `n` (when looking towards +infinity).
+
+    Source: Weissten, Eric W. "Rotation Formula." Mathworld.
+    """
+    n = np.asarray(n)/np.linalg.norm(n)
+    
+    return np.cos(phi)*np.identity(3) + (1-np.cos(phi))*np.outer(n,n) + \
+        np.sin(phi)*np.array([[0,n[2],-n[1]],[-n[2],0,n[0]],[n[1],-n[0],0]])
+
 def rotate(x, phi, n):
     """
     Rotate an array of points `x` through an angle phi counter-clockwise
     around the axis `n` (when looking towards +infinity).
-
-    Source: Weisstein, Eric W. "Rotation Formula." Mathworld.
     """
     x = np.asarray(x)
-    n = np.asarray(n)/np.linalg.norm(n)
 
-    r = np.cos(phi)*np.identity(3) + (1-np.cos(phi))*np.outer(n,n) + \
-        np.sin(phi)*np.array([[0,n[2],-n[1]],[-n[2],0,n[0]],[n[1],-n[0],0]])
+    r = make_rotation_matrix(phi, n)
 
     return np.inner(x,r)
diff --git a/view.py b/view.py
index cc28355..1859195 100755
--- a/view.py
+++ b/view.py
@@ -6,9 +6,12 @@ from pygame.locals import *
 
 import src
 from camera import *
+from solid import *
 from geometry import *
 from transform import *
 
+import time
+
 from pycuda import autoinit
 from pycuda.compiler import SourceModule
 from pycuda import gpuarray
@@ -81,8 +84,12 @@ def view(geometry, name=''):
 
     def render():
         """Render the mesh and display to screen."""
+        t0 = time.time()
         cuda_raytrace(np.int32(pixels.size), origins_gpu, directions_gpu, np.int32(geometry.node_map.size-1), np.int32(geometry.first_node), pixels_gpu, texrefs=texrefs, **gpu_kwargs)
         cuda.Context.synchronize()
+        elapsed = time.time() - t0
+
+        print 'elapsed %f sec' % elapsed
 
         cuda.memcpy_dtoh(pixels, pixels_gpu)
         pygame.surfarray.blit_array(screen, pixels.reshape(size))
@@ -174,8 +181,7 @@ if __name__ == '__main__':
     import sys
     import optparse
 
-    from stl import *
-    from materials import *
+    from mesh import *
     import detectors
 
     parser = optparse.OptionParser('%prog filename.stl')
@@ -191,7 +197,7 @@ if __name__ == '__main__':
 
     if ext.lower() == '.stl':
         geometry = Geometry()
-        geometry.add_solid(Solid(read_stl(args[0]), vacuum, vacuum))
+        geometry.add_solid(Solid(0, mesh_from_stl(args[0])))
         geometry.build(options.bits)
         view(geometry, tail)
     else:
-- 
cgit