From d824208f9c6d5845fade6d91cbe2c937de164894 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 23 Feb 2021 11:30:51 -0800 Subject: [PATCH 1/3] save streamline intersection path dist, add errors --- .../layered_point_depths/__main__.py | 27 ++++++++++++------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/neuron_morphology/layered_point_depths/__main__.py b/neuron_morphology/layered_point_depths/__main__.py index 63e1dd90..717a11ea 100644 --- a/neuron_morphology/layered_point_depths/__main__.py +++ b/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) From ea58c25a74b09b0aee70d172970ee6848fb8621f Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 23 Feb 2021 14:19:35 -0800 Subject: [PATCH 2/3] test for step_from_node path length --- tests/layered_point_depths/test_main.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/layered_point_depths/test_main.py b/tests/layered_point_depths/test_main.py index cd5b7ae3..0332528c 100644 --- a/tests/layered_point_depths/test_main.py +++ b/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): From 56fa7d469227a86a6515affa8c30bb7d3054e083 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Thu, 25 Feb 2021 11:01:01 -0800 Subject: [PATCH 3/3] run_streamlines docstring and schema clarification --- .../transforms/pia_wm_streamlines/_schemas.py | 8 ++--- .../calculate_pia_wm_streamlines.py | 34 ++++++++++++++++--- 2 files changed, 34 insertions(+), 8 deletions(-) diff --git a/neuron_morphology/transforms/pia_wm_streamlines/_schemas.py b/neuron_morphology/transforms/pia_wm_streamlines/_schemas.py index 2b9c044d..990e3bde 100644 --- a/neuron_morphology/transforms/pia_wm_streamlines/_schemas.py +++ b/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') diff --git a/neuron_morphology/transforms/pia_wm_streamlines/calculate_pia_wm_streamlines.py b/neuron_morphology/transforms/pia_wm_streamlines/calculate_pia_wm_streamlines.py index 3bf38cae..a94b369b 100644 --- a/neuron_morphology/transforms/pia_wm_streamlines/calculate_pia_wm_streamlines.py +++ b/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),