Added a sequential ray-tracing function to raypier.core.tracer module

raypier/core/ctracer.pxd

@@ -229,6 +229,7 @@ cdef class FaceList(object):
 
     cpdef void sync_transforms(self)
     cdef int intersect_c(self, ray_t *ray, vector_t end_point)
+    cdef int intersect_one_face_c(self, ray_t *ray, vector_t end_point, int face_idx)
     cdef int intersect_para_c(self, para_t *ray, vector_t ray_end, Face face)
     cdef orientation_t compute_orientation_c(self, Face face, vector_t point)
 
@@ -243,10 +244,24 @@ cdef RayCollection trace_segment_c(RayCollection rays,
                                     list decomp_faces,
                                     float max_length)
 
+cdef RayCollection trace_one_face_segment_c(RayCollection rays,
+                                    FaceList face_set,
+                                    int face_idx,
+                                    list all_faces,
+                                    list decomp_faces,
+                                    float max_length)
+
 cdef GaussletCollection trace_gausslet_c(GaussletCollection gausslets, 
                                     list face_sets, 
                                     list all_faces,
                                     list decomp_faces,
                                     double max_length)
 
+cdef GaussletCollection trace_one_face_gausslet_c(GaussletCollection gausslets, 
+                                    FaceList face_set,
+                                    int face_idx, 
+                                    list all_faces,
+                                    list decomp_faces,
+                                    double max_length)
+
 cdef double ray_power_(ray_t ray)

raypier/core/ctracer.pyx

@@ -1797,6 +1797,27 @@ cdef class FaceList(object):
 
     def __getitem__(self, intidx):
         return self.faces[intidx]
+
+
+    cdef int intersect_one_face_c(self, ray_t *ray, vector_t ray_end, int face_idx):
+        """Finds the intersection with the face specified by face_idx.
+        """
+        cdef:
+            vector_t p1 = transform_c(self.inv_trans, ray.origin)
+            vector_t p2 = transform_c(self.inv_trans, ray_end)
+            list faces=self.faces
+            unsigned int i
+            int all_idx=-1
+            double dist
+            Face face
+
+        face = faces[face_idx]
+        dist = face.intersect_c(p1, p2, 1)
+        if face.tolerance < dist < ray.length:
+            ray.length = dist
+            all_idx = face.idx
+            ray.end_face_idx = all_idx
+        return all_idx
 
 
     cdef int intersect_c(self, ray_t *ray, vector_t ray_end):
@@ -2016,24 +2037,94 @@ cdef RayCollection trace_segment_c(RayCollection rays,
     return new_rays
 
 
+cdef RayCollection trace_one_face_segment_c(RayCollection rays, 
+                                    FaceList face_set,
+                                    int face_idx, 
+                                    list all_faces,
+                                    list decomp_faces,
+                                    float max_length):
+    cdef:
+        Face face
+        size_t i, j
+        vector_t point
+        orientation_t orient
+        int idx, nearest_idx=-1
+        ray_t *ray
+        RayCollection new_rays
+
+    #need to allocate the output rays here 
+    new_rays = RayCollection(rays.n_rays)
+    new_rays.parent = rays
+
+    for i in range(rays.n_rays):
+        ray = rays.rays + i
+        ray.length = max_length
+        ray.end_face_idx = -1
+        nearest_idx=-1
+        point = addvv_(ray.origin, 
+                            multvs_(ray.direction, 
+                                    max_length))
+        #print "points", P1, P2
+        #intersect_c returns the face idx of the intersection, or -1 otherwise
+        idx = (<FaceList>face_set).intersect_one_face_c(ray, point, face_idx)
+        if idx >= 0:
+            nearest_idx = idx
+        if nearest_idx >= 0:
+            #print "GET FACE", nearest.face_idx, len(all_faces)
+            face = all_faces[nearest_idx]
+            face.count += 1
+            #print "ray length", ray.length
+            point = addvv_(ray.origin, multvs_(ray.direction, ray.length))
+            orient = (<FaceList>face_set).compute_orientation_c(face, point)
+            #print "s normal", normal
+            (<InterfaceMaterial>(face.material)).eval_child_ray_c(ray, i, 
+                                                    point,
+                                                    orient,
+                                                    new_rays
+                                                    )
+    return new_rays
+
+
 def trace_segment(RayCollection rays, 
                     list face_sets, 
                     list all_faces,
                     max_length=100,
                     decomp_faces=[]):
+    cdef:
+        FaceList fs
     for fs in face_sets:
         fs.sync_transforms()
     return trace_segment_c(rays, face_sets, all_faces, decomp_faces, max_length)
 
+def trace_one_face_segment(RayCollection rays, 
+                    FaceList face_set,
+                    int face_idx, 
+                    list all_faces,
+                    max_length=100,
+                    decomp_faces=[]):
+    face_set.sync_transforms()
+    return trace_one_face_segment_c(rays, face_set, face_idx, all_faces, decomp_faces, max_length)
+
 def trace_gausslet(GaussletCollection rays, 
                     list face_sets, 
                     list all_faces,
                     max_length=100,
                     decomp_faces=[]):
+    cdef:
+        FaceList fs
     for fs in face_sets:
         (<FaceList>fs).sync_transforms()
     return trace_gausslet_c(rays, face_sets, all_faces, decomp_faces, max_length)
 
+def trace_one_face_gausslet(GaussletCollection rays, 
+                    FaceList face_set,
+                    int face_idx, 
+                    list all_faces,
+                    max_length=100,
+                    decomp_faces=[]):
+    face_set.sync_transforms()
+    return trace_one_face_gausslet_c(rays, face_set, face_idx, all_faces, decomp_faces, max_length)
+
 
 cdef GaussletCollection trace_gausslet_c(GaussletCollection gausslets, 
                                     list face_sets, 
@@ -2102,6 +2193,69 @@ cdef GaussletCollection trace_gausslet_c(GaussletCollection gausslets,
     return new_gausslets
 
 
+cdef GaussletCollection trace_one_face_gausslet_c(GaussletCollection gausslets, 
+                                    FaceList face_set,
+                                    int face_idx, 
+                                    list all_faces,
+                                    list decomp_faces,
+                                    double max_length):
+    cdef:
+        Face face
+        size_t i, j, n_decomp = len(decomp_faces)
+        vector_t point
+        orientation_t orient
+        int idx, nearest_idx=-1
+        gausslet_t *gausslet
+        ray_t *ray
+        GaussletCollection new_gausslets
+        RayCollection child_rays
+
+    #need to allocate the output rays here 
+    new_gausslets = GaussletCollection(gausslets.n_rays)
+    new_gausslets.parent = gausslets
+
+    child_rays = RayCollection(2)
+
+    for i in range(gausslets.n_rays):
+        gausslet = gausslets.rays + i
+        ray = &gausslet.base_ray
+        ray.end_face_idx = -1
+        nearest_idx=-1
+        point = addvv_(ray.origin, 
+                            multvs_(ray.direction, 
+                                    max_length))
+        #print "points", P1, P2
+        #intersect_c returns the face idx of the intersection, or -1 otherwise
+        idx = (<FaceList>face_set).intersect_one_face_c(ray, point, face_idx)
+        if idx >= 0:
+            nearest_idx = idx
+        if nearest_idx >= 0:
+            #print "GET FACE", nearest.face_idx, len(all_faces)
+            face = all_faces[nearest_idx]
+            face.count += 1
+            #print "ray length", ray.length
+            point = addvv_(ray.origin, multvs_(ray.direction, ray.length))
+            orient = face_set.compute_orientation_c(face, point)
+            #print "s normal", normal
+            ### Clear the child_rays structure
+            child_rays.n_rays = 0
+            (<InterfaceMaterial>(face.material)).eval_child_ray_c(ray, i, 
+                                                    point,
+                                                    orient,
+                                                    child_rays
+                                                    )
+            trace_parabasal_rays(gausslet, child_rays, face, face_set, new_gausslets, max_length)
+
+    for j in range(n_decomp):
+        face = decomp_faces[j]
+        if face.count:
+            (<InterfaceMaterial>(face.material)).eval_decomposed_rays_c(new_gausslets)
+            face.count = 0
+
+    new_gausslets.reset_length_c(max_length)
+    return new_gausslets
+
+
 cdef void trace_parabasal_rays(gausslet_t *g_in, RayCollection base_rays, Face face, FaceList face_set, 
                           GaussletCollection new_gausslets, double max_length):
     cdef:

raypier/core/tracer.py

@@ -1,5 +1,6 @@
 
-from .ctracer import trace_segment, trace_gausslet, RayCollection
+from .ctracer import trace_segment, trace_gausslet, RayCollection, trace_one_face_segment,
+        trace_one_face_gausslet
 
 from itertools import chain
 import numpy
@@ -45,3 +46,53 @@ def trace_rays(input_rays, face_lists, recursion_limit=100, max_length=100.0):
 
     return traced_rays, all_faces 
 
+
+def trace_ray_sequence(input_rays, face_lists, face_idx_list, recursion_limit=100, max_length=100.0):
+    """
+    Tracing routine for a sequential ray-trace. Takes either a RayCollection or GaussletCollection
+    and traces the rays sequentially through the given list of FaceList objects.
+    
+    If you want to trace multiple faces within one FaceList object, then the FaceList should have 
+    multiple entries in the face_lists list with corresponding face_idx entries in the face_idx_list list. 
+    
+     - face_idx_list : list of indices specifying which face in each FaceList should be traced.
+    
+    The input_rays should already have a consistent wavelengths property set.
+    
+    returns - (traced_rays, all_faces)
+            where traced_rays is a list of RayCollection or GaussletCollections 
+            representing the sequence of ray generations. The 'all_faces' list
+            is a list of cfaces.Face objects. The end_face_idx member of each
+            traced ray indexes into this list to give the face where it terminates.
+    """
+    input_rays.reset_length(max_length)
+    traced_rays = []
+    trace_func = trace_one_face_segment if isinstance(input_rays, RayCollection) else trace_one_face_gausslet
+    count = 0
+    wavelengths = numpy.asarray(input_rays.wavelengths)
+
+    all_faces = list(chain(*(fs.faces for fs in face_lists)))
+    for i, f in enumerate(all_faces):
+        f.idx = i
+        f.count = 0 #reset intersection count
+        f.update() #This was a bad idea.
+        f.material.wavelengths = wavelengths
+        f.max_length = max_length
+
+    face_sets = list(face_lists)
+    decomp_faces = [f for f in all_faces if f.material.is_decomp_material()]
+
+    rays = input_rays
+
+    face_iter = zip(face_lists, face_idx_list)
+
+    while rays.n_rays>0 and count<recursion_limit:
+        traced_rays.append(rays)
+        next_face_list, next_face_idx = next(face_iter)
+        rays = trace_func(rays, next_face_list, next_face_idx, 
+                          all_faces, 
+                         max_length=max_length,
+                         decomp_faces=decomp_faces)
+        count += 1
+
+    return traced_rays, all_faces