Merge pull request #27 from jobh/trace

Simplify trace_ray, remove heap mallocs

pykonal/fields.pyx

@@ -33,7 +33,6 @@ from . import transformations
 # Cython built-in imports.
 from libc.math cimport sqrt, sin
 from libcpp.vector cimport vector as cpp_vector
-from libc.stdlib   cimport malloc, free
 
 # Third-party Cython imports
 cimport numpy as np
@@ -418,26 +417,33 @@ cdef class ScalarField3D(Field3D):
         :rtype:  numpy.ndarray(shape=(N,3), dtype=numpy.float)
         """
 
-        cdef cpp_vector[constants.REAL_t *]       ray
+        cdef cpp_vector[constants.REAL_t]         ray
         cdef constants.REAL_t                     norm, step_size, value, value_1back
-        cdef constants.REAL_t                     *point_new
-        cdef constants.REAL_t[3]                  gg, point, point_1back
-        cdef Py_ssize_t                           idx, jdx
-        cdef np.ndarray[constants.REAL_t, ndim=2] ray_np
+        cdef constants.REAL_t[3]                  gg
+        cdef constants.REAL_t[:]                  point
+        cdef Py_ssize_t                           idx
         cdef str                                  coord_sys
+        cdef np.ndarray[constants.REAL_t, ndim=1] ray_np
         cdef VectorField3D                        grad
 
         coord_sys = self.coord_sys
         step_size = self.step_size
         grad      = self.gradient
 
-        point_new = <constants.REAL_t *> malloc(3 * sizeof(constants.REAL_t))
-        point_new[0], point_new[1], point_new[2] = end
-        ray.push_back(point_new)
-        point = ray.back()
-        value = self.value(point)
+        for idx in range(3):
+            ray.push_back(end[idx])
+        value = np.infty
 
         while True:
+            point = <constants.REAL_t[:3]>&ray[ray.size()-3]
+            try:
+                value_1back = value
+                value = self.value(point)
+                if value >= value_1back or np.isnan(value):
+                    ray.resize(ray.size()-3) # drop bad pt
+                    break
+            except ValueError:
+                return (None)
             gg   = grad.value(point)
             norm = sqrt(gg[0]**2 + gg[1]**2 + gg[2]**2)
             if np.isnan(norm):
@@ -447,27 +453,13 @@ cdef class ScalarField3D(Field3D):
             if coord_sys == "spherical":
                 gg[1] /= point[0]
                 gg[2] /= point[0] * sin(point[1])
-            point_new = <constants.REAL_t *> malloc(3 * sizeof(constants.REAL_t))
             for idx in range(3):
-                point_new[idx] = point[idx] - step_size * gg[idx]
-            point_1back = ray.back()
-            ray.push_back(point_new)
-            point  = ray.back()
-            try:
-                value_1back = value
-                value = self.value(point)
-                if value_1back <= value or np.isnan(value):
-                    break
-            except ValueError:
-                for idx in range(ray.size()-1):
-                    free(ray[idx])
-                return (None)
-        ray_np = np.zeros((ray.size()-1, 3), dtype=constants.DTYPE_REAL)
-        for idx in range(ray.size()-1):
-            for jdx in range(3):
-                ray_np[idx, jdx] = ray[idx][jdx]
-            free(ray[idx])
-        return (np.flipud(ray_np))
+                ray.push_back(point[idx] - step_size * gg[idx])
+
+        ray_np = np.empty(ray.size(), dtype=constants.DTYPE_REAL)
+        for idx in range(ray_np.size):
+            ray_np[idx] = ray[idx]
+        return (np.flipud(ray_np.reshape(-1,3)))
 
 
     cpdef constants.REAL_t value(