Layer depths update

neuron_morphology/layered_point_depths/__main__.py

@@ -139,12 +139,15 @@ def step_from_node(
 
     Returns
     -------
-    The depth of the intersection between the path walked and the given surface
+    depth: depth of the intersection between the path walked and the given surface
+    path_dist: the path distance from start to the intersection point
 
     """
 
     retry_step = False
     cur_pos = np.array(list(pos))
+    path_dist = 0
+    path = [cur_pos]
 
     for _ in range(max_iter):
         if not retry_step:
@@ -154,7 +157,7 @@ def step_from_node(
 
             base_step = np.squeeze([dx, dy])
             if np.any(np.isnan(base_step)):
-                return None
+                raise ValueError("Path is outside of interpolator domain.")
             base_step = base_step / np.linalg.norm(base_step)
         step = adaptive_scale * step_size * base_step
 
@@ -178,15 +181,19 @@ def step_from_node(
                         if test_dist < dist:
                             dist = test_dist
                             closest_pt = test_pt
-                    intersection_pt = list(closest_pt.coords)
+                    intersection_pt = list(closest_pt.coords)[0]
                 else:
-                    intersection_pt = list(intersection.coords)
-                return float(depth_interp(intersection_pt[0]))
+                    intersection_pt = list(intersection.coords)[0]
+                path_dist += np.linalg.norm(intersection_pt-cur_pos)
+                path.append(intersection_pt)
+                return float(depth_interp(intersection_pt)), path_dist
 
         cur_pos = next_pos
+        path_dist += np.linalg.norm(step)
+        path.append(next_pos)
         retry_step = False
-
-    return None # uneccessary, but satisfies linter :/
+    
+    raise ValueError("Path failed to intersect surface.")
 
 
 def get_node_intersections(
@@ -267,10 +274,10 @@ def get_node_intersections(
         "depth": depth,
         "local_layer_pia_side_depth": step_from_node(
             pos, depth_interp, dx_interp, dy_interp, pia, step_size, max_iter
-        ),
+        )[0],
         "local_layer_wm_side_depth": step_from_node(
             pos, depth_interp, dx_interp, dy_interp, wm, -step_size, max_iter
-        ),
+        )[0],
         "point_type": node["type"]
     }
 
@@ -362,7 +369,7 @@ def main():
     args = cp.deepcopy(parser.args)
     logging.getLogger().setLevel(args.pop("log_level"))
 
-    output = main(**args)
+    output = run_layered_point_depths(**args)
     output.update({"inputs": parser.args})
 
     parser.output(output)

neuron_morphology/transforms/pia_wm_streamlines/_schemas.py

@@ -9,13 +9,13 @@ class PiaWmStreamlineSchema(ArgSchema):
 
     pia_path_str = String(
         required=True,
-        description='string alternating x, y coordinates outlining the pia')
+        description='string alternating x, y coordinates outlining the pia (in pixels)')
     wm_path_str = String(
         required=True,
-        description='string alternating x, y coordinates outlining the wm')
+        description='string alternating x, y coordinates outlining the wm (in pixels)')
     soma_path_str = String(
         required=False,
-        description='string alternating x, y coordinates outlining the soma. '
+        description='string alternating x, y coordinates outlining the soma (in pixels). '
                     'If provided, streamlines will be translated so that '
                     'the origin is at the soma')
     resolution = Float(required=False,
@@ -29,7 +29,7 @@ class PiaWmStreamlineSchema(ArgSchema):
                            description='Fixed value wm boundary condition')
     mesh_res = Int(required=False,
                    default=20,
-                   description='Resolution for mesh for laplace solver')
+                   description='Resolution for mesh for laplace solver (microns)')
     output_dir = OutputDir(required=True,
                            description='Directory to write xarray results')
 

neuron_morphology/transforms/pia_wm_streamlines/calculate_pia_wm_streamlines.py

@@ -33,9 +33,35 @@ def run_streamlines(pia_path_str: str,
          wm_path_str: str,
          resolution: float,
          soma_path_str: Optional[str] = None,
-         mesh_res: int = 20,
+         mesh_res: float = 20,
          pia_fixed_value: float = 1.0,
-         wm_fixed_value: float = 0.0,):
+         wm_fixed_value: float = 0.0,
+         grid_res: float = 1.0,):
+    """
+        Continuously interpolate 'depth' coordinates and gradients
+        (directions of fastest ascent) via Laplace's equation,
+        for all points lying between two lines (pia and wm).
+        If pia/wm fixed_value defaults are used, depth_field
+        will be the depth from the top, normalized by region thickness.
+
+        Parameters
+        ----------
+        pia_path_str: path string for pia boundary
+        wm_path_str: path string for pia boundary
+        soma_path_str: path string for soma - if provided, 
+            result coordinates are centered about the soma
+        mesh_res: resolution of the mesh for equation solving
+        pia_fixed_value: 'depth' field value for pia/top boundary
+        wm_fixed_value: 'depth' field value for wm/bottom boundary
+        grid_res: resolution of coordinate grid for interpolated output
+
+        Returns
+        -------
+        depth_field: DataArray for depth values on grid
+        gradient_field: DataArray for gradient vectors on grid
+        translation: coordinate shift applied to center around soma
+
+    """
 
     pia_path = convert_path_str_to_list(pia_path_str, resolution)
     wm_path = convert_path_str_to_list(wm_path_str, resolution)
@@ -60,8 +86,8 @@ def run_streamlines(pia_path_str: str,
 
     x = [coord[0] for coord in mesh_coords]
     y = [coord[1] for coord in mesh_coords]
-    xx = np.arange(min(x)-1, max(x)+1, 1)
-    yy = np.arange(min(y)-1, max(y)+1, 1)
+    xx = np.arange(min(x)-grid_res, max(x)+grid_res, grid_res)
+    yy = np.arange(min(y)-grid_res, max(y)+grid_res, grid_res)
 
     grid_x, grid_y = np.meshgrid(xx, yy, indexing='ij')
     grid_u = griddata((x, y),

tests/layered_point_depths/test_main.py

@@ -85,13 +85,17 @@ def test_step_from_node(self):
 
         surface = LineString([(0, 2), (4, 2)])
 
-        depth = lpd.step_from_node(
+        depth, path = lpd.step_from_node(
             pos, depth_interp, dx_interp, dy_interp, surface, 1.0, 1000)
 
         self.assertAlmostEqual(
             depth, 
             depth_interp((1, 2))
         )
+        self.assertAlmostEqual(
+            path, 
+            np.sqrt(5)
+        )
 
     def test_get_node_intersections(self):