summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--detectors/lbne.py5
-rw-r--r--geometry.py126
-rw-r--r--mesh.py113
-rw-r--r--solid.py58
-rw-r--r--src/intersect.h6
-rw-r--r--stl.py43
-rw-r--r--transform.py18
-rwxr-xr-xview.py12
8 files changed, 239 insertions, 142 deletions
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:
generated by cgit (git 2.34.1) at 2025-05-08 05:34:34 +0000