From 94b57fcbd0f952fca9f07e13523d42827cea03cd Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 23 Feb 2021 11:30:51 -0800 Subject: [PATCH 01/13] 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 354675574cd809392c7d818ee70a89cd075e9579 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 23 Feb 2021 14:19:35 -0800 Subject: [PATCH 02/13] 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 796b668201f57b60b618c1133abb387afae3273e Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 6 Oct 2020 10:49:35 -0700 Subject: [PATCH 03/13] update FromLimsSource for current argschema --- neuron_morphology/snap_polygons/_from_lims.py | 25 ++++++++++--------- 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/neuron_morphology/snap_polygons/_from_lims.py b/neuron_morphology/snap_polygons/_from_lims.py index 4ad3e8a1..7e5aa569 100644 --- a/neuron_morphology/snap_polygons/_from_lims.py +++ b/neuron_morphology/snap_polygons/_from_lims.py @@ -19,10 +19,10 @@ import warnings import logging -import marshmallow as mm +from argschema.schemas import DefaultSchema from argschema.fields import Int, OutputDir, String -from argschema.sources import ArgSource +from argschema.sources import ConfigurableSource from allensdk.internal.core import lims_utilities as lu from neuron_morphology.snap_polygons.types import ( @@ -262,7 +262,7 @@ def get_inputs_from_lims( return results -class PostgresInputConfigSchema(mm.Schema): +class PostgresInputConfigSchema(DefaultSchema): """The parameters required to query a postgres database. """ @@ -285,7 +285,7 @@ class PostgresInputConfigSchema(mm.Schema): ) port = Int( description="", - required=False, # seems not to get hydrated from the default + required=False, default=5432 ) @@ -307,20 +307,21 @@ class FromLimsSchema(PostgresInputConfigSchema): ) -class FromLimsSource(ArgSource): +class FromLimsSource(ConfigurableSource): """ An alternate argschema source which gets its inputs from lims directly """ ConfigSchema = FromLimsSchema def get_dict(self): - image_output = getattr(self, "image_output_root", None) + config = self.config + image_output = getattr(config, "image_output_root", None) return get_inputs_from_lims( - self.host, # pylint: disable=no-member - self.port, # pylint: disable=no-member - self.database, # pylint: disable=no-member - self.user, # pylint: disable=no-member - self.password, # pylint: disable=no-member, - self.focal_plane_image_series_id, # pylint: disable=no-member + config["host"], # pylint: disable=no-member + config["port"], # pylint: disable=no-member + config["database"], # pylint: disable=no-member + config["user"], # pylint: disable=no-member + config["password"], # pylint: disable=no-member, + config["focal_plane_image_series_id"], # pylint: disable=no-member image_output # pylint: disable=no-member ) From 134946b286185f11702cd5dcb706d7518753e189 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 6 Oct 2020 10:50:14 -0700 Subject: [PATCH 04/13] use lims env variables for defaults --- neuron_morphology/lims_apical_queries.py | 8 ++++---- neuron_morphology/snap_polygons/_from_lims.py | 11 +++++++---- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/neuron_morphology/lims_apical_queries.py b/neuron_morphology/lims_apical_queries.py index a58e2826..6320a795 100644 --- a/neuron_morphology/lims_apical_queries.py +++ b/neuron_morphology/lims_apical_queries.py @@ -20,10 +20,10 @@ def convert_coords_str(coords_str: str, resolution=None): def get_data(query): - host = os.getenv("DBHOST") - dbname = os.getenv("DBNAME") - user = os.getenv("DBREADER") - password = os.getenv("DBPASSWORD") + host = os.getenv("LIMS_HOST") + dbname = os.getenv("LIMS_DBNAME") + user = os.getenv("LIMS_USER") + password = os.getenv("LIMS_PASSWORD") conn_str = f'host={host} dbname={dbname} user={user} password={password}' data = {} diff --git a/neuron_morphology/snap_polygons/_from_lims.py b/neuron_morphology/snap_polygons/_from_lims.py index 7e5aa569..b5f4dc99 100644 --- a/neuron_morphology/snap_polygons/_from_lims.py +++ b/neuron_morphology/snap_polygons/_from_lims.py @@ -268,20 +268,23 @@ class PostgresInputConfigSchema(DefaultSchema): host = String( description="", - required=True + required=False, + default=os.environ.get("LIMS_HOST") ) database = String( description="", - required=True + required=False, + default=os.environ.get("LIMS_DBNAME") ) user = String( description="", - required=True + required=False, + default=os.environ.get("LIMS_USER") ) password = String( description="", required=False, - default=os.environ.get("POSTGRES_SOURCE_PASSWORD") + default=os.environ.get("LIMS_PASSWORD") ) port = Int( description="", From aeb2b98b344482f5474b5ba431403339958961c8 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 6 Oct 2020 11:48:27 -0700 Subject: [PATCH 05/13] cleanup FromLimsSource edits --- neuron_morphology/snap_polygons/_from_lims.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/neuron_morphology/snap_polygons/_from_lims.py b/neuron_morphology/snap_polygons/_from_lims.py index b5f4dc99..0d434a83 100644 --- a/neuron_morphology/snap_polygons/_from_lims.py +++ b/neuron_morphology/snap_polygons/_from_lims.py @@ -246,13 +246,11 @@ def get_inputs_from_lims( layer_polygons = query_for_layer_polygons(engine, imser_id) pia_surface, wm_surface = query_for_cortical_surfaces(engine, imser_id) - image_width, image_height = query_for_image_dims(engine, imser_id) results = { "layer_polygons": layer_polygons, "pia_surface": pia_surface, "wm_surface": wm_surface, - "image_dimensions": {"width": image_width, "height": image_height} } if image_output_root is not None: @@ -318,13 +316,12 @@ class FromLimsSource(ConfigurableSource): def get_dict(self): config = self.config - image_output = getattr(config, "image_output_root", None) return get_inputs_from_lims( - config["host"], # pylint: disable=no-member - config["port"], # pylint: disable=no-member - config["database"], # pylint: disable=no-member - config["user"], # pylint: disable=no-member - config["password"], # pylint: disable=no-member, - config["focal_plane_image_series_id"], # pylint: disable=no-member - image_output # pylint: disable=no-member + config["host"], + config["port"], + config["database"], + config["user"], + config["password"], + config["focal_plane_image_series_id"], + config.get("image_output_root") ) From 3e1bcffa58ed42e595b830d3df90289bd6d5ed54 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 6 Oct 2020 11:55:43 -0700 Subject: [PATCH 06/13] parser updates for lims source --- neuron_morphology/snap_polygons/__main__.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/neuron_morphology/snap_polygons/__main__.py b/neuron_morphology/snap_polygons/__main__.py index f5d0c873..3649a844 100644 --- a/neuron_morphology/snap_polygons/__main__.py +++ b/neuron_morphology/snap_polygons/__main__.py @@ -83,20 +83,19 @@ def run_snap_polygons( return results +class Parser(ArgSchemaParser): + """An ArgschemaParser that can pull data from LIMS + """ + default_sources = \ + ArgSchemaParser.default_sources + (FromLimsSource,) + default_schema=InputParameters + default_output_schema=OutputParameters + def main(): """CLI entrypoint for snapping polygons """ - class Parser(ArgSchemaParser): - """An ArgschemaParser that can pull data from LIMS - """ - default_configurable_sources = \ - ArgSchemaParser.default_configurable_sources + [FromLimsSource] - - parser = Parser( - schema_type=InputParameters, - output_schema_type=OutputParameters - ) + parser = Parser() args = cp.deepcopy(parser.args) logging.getLogger().setLevel(args.pop("log_level")) From 6c7118432ab0ad5be5542dedefabaf8ecd4e644a Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 13 Oct 2020 13:10:52 -0700 Subject: [PATCH 07/13] query_for_layer_polygons select distinct --- neuron_morphology/snap_polygons/_from_lims.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/neuron_morphology/snap_polygons/_from_lims.py b/neuron_morphology/snap_polygons/_from_lims.py index 0d434a83..402a742b 100644 --- a/neuron_morphology/snap_polygons/_from_lims.py +++ b/neuron_morphology/snap_polygons/_from_lims.py @@ -58,7 +58,7 @@ def query_for_layer_polygons( """ query = f""" - select + select distinct st.acronym as name, polygon.path as path, polygon.id as polygon_id @@ -105,11 +105,6 @@ def query_for_layer_polygons( raise ValueError( f"found multiple distinct layer drawings for {name}" ) - else: - warnings.warn( - f"found multiple polygon records for {name} " - "(identical paths)" - ) polygons.append({ "name": name, From 97827810e5235e58891632ad51358d279917a953 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 23 Feb 2021 15:24:06 -0800 Subject: [PATCH 08/13] clarify find_vertical_surfaces behavior for pia/wm --- neuron_morphology/snap_polygons/geometries.py | 33 +++++++++++-------- tests/snap_polygons/test_geometries.py | 2 +- 2 files changed, 21 insertions(+), 14 deletions(-) diff --git a/neuron_morphology/snap_polygons/geometries.py b/neuron_morphology/snap_polygons/geometries.py index 333bf95b..3b7d52cb 100644 --- a/neuron_morphology/snap_polygons/geometries.py +++ b/neuron_morphology/snap_polygons/geometries.py @@ -512,14 +512,13 @@ def clear_overlaps(stack: Dict[str, np.ndarray]): def closest_from_stack(stack: Dict[str, np.ndarray]): - """ Given a stack of images describing distance from several objects, find - the closest object to each pixel. + """ Given a stack of image masks representing several objects, find + the closest object to each pixel in the image space. Parameters ---------- - stack : Keys are names, values are ndarrays (of the same shape). Each pixel - in the values describes the distance from that pixel to the named - object + stack : Keys are names, values are ndarrays (of the same shape) containing + binary masks of objects. Returns ------- @@ -576,8 +575,10 @@ def find_vertical_surfaces( pia: Optional[LineString] = None, white_matter: Optional[LineString] = None ): - """ Given a set of polygons describing cortical layer boundaries, find the + """ Given a set of polygons describing cortical layers, find the boundaries between each layer. + Pia and white matter surfaces are optional, used for pia-side of top layer + and wm-side of bottom layer. Otherwise, these are not assigned. Parameters ---------- @@ -601,18 +602,22 @@ def find_vertical_surfaces( for index, name in enumerate(names): current = polygons[name] # up side - if index == 0 and pia is not None: - results[f"{name}_pia"] = pia + if index == 0: + top = pia else: above_layers = [polygons[name] for name in names[:index]] - results[f"{name}_pia"] = shared_faces(current, above_layers) + top = shared_faces(current, above_layers) + if top is not None: + results[f"{name}_pia"] = top # down side - if index == len(names) - 1 and white_matter is not None: - results[f"{name}_wm"] = white_matter + if index == len(names) - 1: + bottom = white_matter else: below_layers = [polygons[name] for name in names[index + 1:]] - results[f"{name}_wm"] = shared_faces(current, below_layers) + bottom = shared_faces(current, below_layers) + if bottom is not None: + results[f"{name}_wm"] = bottom return results @@ -630,7 +635,7 @@ def shared_faces(poly: Polygon, others: Iterable[Polygon], snap_tolerance=10) -> Returns ------- - LineString representing the shared face + LineString representing the shared face, or None if not found """ merged_others = shapely.ops.unary_union(others) @@ -658,6 +663,8 @@ def shared_faces(poly: Polygon, others: Iterable[Polygon], snap_tolerance=10) -> flat_faces_list[i + 1] = shapely.geometry.LineString(list(f.coords) + list(f_prev.boundary.geoms[prev_ind].coords)) faces = shapely.ops.linemerge(flat_faces_list) + if not faces: + return None coordinates = list(faces.coords) shared_line = ensure_linestring(coordinates) return shared_line diff --git a/tests/snap_polygons/test_geometries.py b/tests/snap_polygons/test_geometries.py index e1f3319d..2decce73 100644 --- a/tests/snap_polygons/test_geometries.py +++ b/tests/snap_polygons/test_geometries.py @@ -213,7 +213,7 @@ def test_get_snapped_polys(self): [(1.0, 1.0), (1.0, 3.0), (3.0, 3.0), (3.0, 1.0), (1.0, 1.0)] ) - def find_vertical_surfaces(self): + def test_find_vertical_surfaces(self): polys = { "layer1": Polygon([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]), "layer2": Polygon([(0, 1), (1, 1), (1, 2), (0, 2), (0, 1)]) From c9bca9cac9dfed037b081b3fabcae984eeed8f70 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Tue, 23 Feb 2021 17:04:09 -0800 Subject: [PATCH 09/13] allow missing pia/wm paths --- neuron_morphology/snap_polygons/__main__.py | 40 +++++++++++-------- neuron_morphology/snap_polygons/_from_lims.py | 2 +- neuron_morphology/snap_polygons/_schemas.py | 8 +++- 3 files changed, 30 insertions(+), 20 deletions(-) diff --git a/neuron_morphology/snap_polygons/__main__.py b/neuron_morphology/snap_polygons/__main__.py index 3649a844..118466fe 100644 --- a/neuron_morphology/snap_polygons/__main__.py +++ b/neuron_morphology/snap_polygons/__main__.py @@ -33,26 +33,32 @@ def run_snap_polygons( """Finds and returns close fit boundaries. May write diagnostic images as a side effect. """ - + if len(layer_polygons)==0: + raise ValueError("No polygons provided.") + layer_names = [layer['name'] for layer in layer_polygons] + if len(layer_names) != len(set(layer_names)): + raise ValueError("Duplicate layer names.") # setup input geometries geometries = Geometries() geometries.register_polygons(layer_polygons) # setup cortex boundaries - hull = geometries.convex_hull() - pia = trim_to_close(hull, surface_distance_threshold, pia_surface["path"]) - white_matter = trim_to_close( - hull, surface_distance_threshold, wm_surface["path"] - ) - - geometries.register_surface("pia", pia) - geometries.register_surface("wm", white_matter) - - pia_wm_vertices = get_vertices_from_two_lines( - pia.coords[:], white_matter.coords[:] - ) - bounds = shapely.geometry.polygon.Polygon(pia_wm_vertices) - + if pia_surface is not None and wm_surface is not None: + hull = geometries.convex_hull() + pia = trim_to_close(hull, surface_distance_threshold, pia_surface["path"]) + white_matter = trim_to_close( + hull, surface_distance_threshold, wm_surface["path"] + ) + geometries.register_surface("pia", pia) + geometries.register_surface("wm", white_matter) + pia_wm_vertices = get_vertices_from_two_lines( + pia.coords[:], white_matter.coords[:] + ) + # is this ever a good idea? seems like it may cut too much + bounds = shapely.geometry.polygon.Polygon(pia_wm_vertices) + else: + bounds = geometries.convex_hull() + multipolygon_resolver = partial( select_largest_subpolygon, error_threshold=multipolygon_error_threshold @@ -69,8 +75,8 @@ def run_snap_polygons( boundaries = find_vertical_surfaces( result_geos.polygons, layer_order, - pia=geometries.surfaces["pia"], - white_matter=geometries.surfaces["wm"] + pia=geometries.surfaces.get("pia"), + white_matter=geometries.surfaces.get("wm") ) result_geos.register_surfaces(boundaries) diff --git a/neuron_morphology/snap_polygons/_from_lims.py b/neuron_morphology/snap_polygons/_from_lims.py index 402a742b..77be8f56 100644 --- a/neuron_morphology/snap_polygons/_from_lims.py +++ b/neuron_morphology/snap_polygons/_from_lims.py @@ -154,7 +154,7 @@ def query_for_cortical_surfaces( "name": item["name"], "path": ensure_path(item["path"]) } - return results["Pia"], results["White Matter"] + return results.get("Pia"), results.get("White Matter") def query_for_images( diff --git a/neuron_morphology/snap_polygons/_schemas.py b/neuron_morphology/snap_polygons/_schemas.py index 86cb64e7..ff50a2b6 100644 --- a/neuron_morphology/snap_polygons/_schemas.py +++ b/neuron_morphology/snap_polygons/_schemas.py @@ -68,14 +68,18 @@ class InputParameters(ArgSchema): pia_surface = Nested( SimpleGeometry, description="A path defining the pia-side surface of the cortex", - required=True + required=True, + default=None, + allow_none=True ) wm_surface = Nested( SimpleGeometry, description=( "A path defining the white matter-side surface of the cortex" ), - required=True + required=False, + default=None, + allow_none=True ) working_scale = Float( description=( From 03f1ad18471151071e91fe6342d7169fc51b0573 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Fri, 26 Feb 2021 16:42:28 -0800 Subject: [PATCH 10/13] remove duplicate layer error to pass tests --- neuron_morphology/snap_polygons/__main__.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/neuron_morphology/snap_polygons/__main__.py b/neuron_morphology/snap_polygons/__main__.py index 118466fe..62f75f82 100644 --- a/neuron_morphology/snap_polygons/__main__.py +++ b/neuron_morphology/snap_polygons/__main__.py @@ -35,9 +35,6 @@ def run_snap_polygons( """ if len(layer_polygons)==0: raise ValueError("No polygons provided.") - layer_names = [layer['name'] for layer in layer_polygons] - if len(layer_names) != len(set(layer_names)): - raise ValueError("Duplicate layer names.") # setup input geometries geometries = Geometries() geometries.register_polygons(layer_polygons) From a4c5cc7a81aa824c23162258cf1830f428c06423 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Mon, 18 Oct 2021 19:06:23 -0700 Subject: [PATCH 11/13] improve handling of invalid shapes, logging --- neuron_morphology/snap_polygons/geometries.py | 6 ++++-- neuron_morphology/snap_polygons/types.py | 9 ++++++--- neuron_morphology/transforms/geometry.py | 2 +- 3 files changed, 11 insertions(+), 6 deletions(-) diff --git a/neuron_morphology/snap_polygons/geometries.py b/neuron_morphology/snap_polygons/geometries.py index 3b7d52cb..ef42fd9c 100644 --- a/neuron_morphology/snap_polygons/geometries.py +++ b/neuron_morphology/snap_polygons/geometries.py @@ -372,12 +372,14 @@ def cut( for key, polygon in self.polygons.items(): polygon = polygon.intersection(template) polygon = multipolygon_resolver(polygon) - result.register_polygon(key, polygon) + if not polygon.is_empty: + result.register_polygon(key, polygon) for key, surface in self.surfaces.items(): surface = surface.intersection(template) surface = multisurface_resolver(surface) - result.register_surface(key, surface) + if not surface.is_empty: + result.register_surface(key, surface) return result diff --git a/neuron_morphology/snap_polygons/types.py b/neuron_morphology/snap_polygons/types.py index a17c6258..8941a3f6 100644 --- a/neuron_morphology/snap_polygons/types.py +++ b/neuron_morphology/snap_polygons/types.py @@ -40,13 +40,16 @@ def ensure_polygon(candidate: PolyType) -> Polygon: candidate = ensure_path(candidate) if isinstance(candidate, Polygon): - return candidate + poly = candidate elif isinstance(candidate, LinearRing): - return Polygon(candidate) + poly = Polygon(candidate) elif isinstance(candidate, collections.abc.Sequence): - return Polygon([item for item in map(tuple, candidate)]) + poly = Polygon([item for item in map(tuple, candidate)]) else: raise TypeError(f"did not understand type: {type(candidate)}") + if not poly.is_valid: + poly = poly.buffer(0) + return poly def ensure_linestring(candidate: LineType) -> LineString: diff --git a/neuron_morphology/transforms/geometry.py b/neuron_morphology/transforms/geometry.py index 90b419dd..7cb1266c 100644 --- a/neuron_morphology/transforms/geometry.py +++ b/neuron_morphology/transforms/geometry.py @@ -55,7 +55,7 @@ def prune_two_lines(line1: List[Tuple], line2: List[Tuple]): prune = True if prune: - warnings.warn(f"lines are modified \nline1: {line1}\nline2: {line2}", UserWarning) + warnings.warn(f"Trimming pia/wm to simplify region, lines are likely too curved.", UserWarning) return line1, line2 From 07fc8a7d16e94d3df718a50208e0d268eba99516 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Mon, 18 Oct 2021 19:09:27 -0700 Subject: [PATCH 12/13] lims queries skip invalid shapes, allow multiples Will now return multiple distinct shapes with same label if validate_polys=False. --- neuron_morphology/snap_polygons/_from_lims.py | 143 ++++++++++-------- 1 file changed, 80 insertions(+), 63 deletions(-) diff --git a/neuron_morphology/snap_polygons/_from_lims.py b/neuron_morphology/snap_polygons/_from_lims.py index 77be8f56..68c5291c 100644 --- a/neuron_morphology/snap_polygons/_from_lims.py +++ b/neuron_morphology/snap_polygons/_from_lims.py @@ -4,22 +4,23 @@ Example Usage ------------- python -m neuron_morphology.snap_polygons - --host - --port - --user - --password - --database + --host + --port + --user + --password + --database --focal_plane_image_series_id 522408212 # for instance --image_output_root /some_directory """ from typing import Callable, List, Dict, Tuple, Union, Optional +from collections import defaultdict from functools import partial import os import warnings import logging - +import numpy as np from argschema.schemas import DefaultSchema from argschema.fields import Int, OutputDir, String from argschema.sources import ConfigurableSource @@ -34,7 +35,7 @@ def query_for_layer_polygons( - query_engine: QueryEngineType, + query_engine: QueryEngineType, focal_plane_image_series_id: int, validate_polys: bool = True, treatment: str = "Biocytin" @@ -43,13 +44,13 @@ def query_for_layer_polygons( Parameters ---------- - query_engine : executes a query, passed in as a string. Must not require + query_engine : executes a query, passed in as a string. Must not require any additional database information. focal_plane_image_series_id : used to determine which polygons to fetch validate_polys : if True, fail when - a label is associated with multiple distinct valid geometries - a label is associated with one or more geometries, but none are valid - treatment: The layer polygons are associated with Biocytin and DAPI + treatment: The layer polygons are associated with Biocytin and DAPI treatments. We only need one. Returns @@ -62,9 +63,7 @@ def query_for_layer_polygons( st.acronym as name, polygon.path as path, polygon.id as polygon_id - from specimens sp - join specimens spp on spp.id = sp.parent_id - join image_series imser on imser.specimen_id = spp.id + from image_series imser join sub_images si on si.image_series_id = imser.id join images im on im.id = si.image_id join treatments tm on tm.id = im.treatment_id @@ -72,45 +71,57 @@ def query_for_layer_polygons( join avg_group_labels label on label.id = layer.group_label_id join avg_graphic_objects polygon on polygon.parent_id = layer.id join structures st on st.id = polygon.cortex_layer_id - where + where imser.id = {focal_plane_image_series_id} and label.name in ('Cortical Layers') and tm.name = '{treatment}' """ - + results = query_engine(query) + validate_fcn = ensure_polygon + polygons = sort_lims_graphics_objects( + results, validate_fcn=validate_fcn, validate_polys=validate_polys + ) + return polygons + +def sort_lims_graphics_objects(results, validate_fcn=None, validate_polys=True): polygons = [] candidate_names = set() - found_names: Dict[str, NicePathType] = {} + found_names: Dict[str, list(NicePathType)] = defaultdict(list) - for result in query_engine(query): + for result in results: name = result["name"] path = ensure_path(result["path"]) poly_id = result["polygon_id"] candidate_names.add(name) - if validate_polys: + if validate_fcn: try: - ensure_polygon(path) + if not validate_fcn(path): + continue except (ValueError, TypeError, IndexError): warnings.warn( - "unable to build shapely object from avg graphic " - f"object {poly_id} (label: {name})" + f"invalid graphic object {poly_id} (label: {name})" ) continue if name in found_names: - if path != found_names[name]: + if not any([np.array_equal(path, x) for x in found_names[name]]): if validate_polys: raise ValueError( - f"found multiple distinct layer drawings for {name}" + f"found multiple distinct polygons by name: {name}" ) - - polygons.append({ - "name": name, - "path": path - }) - found_names[name] = path + polygons.append({ + "name": name, + "path": path + }) + found_names[name].append(path) + else: + polygons.append({ + "name": name, + "path": path + }) + found_names[name].append(path) invalid = candidate_names - set(found_names.keys()) if validate_polys and invalid: @@ -122,60 +133,66 @@ def query_for_layer_polygons( def query_for_cortical_surfaces( query_engine: QueryEngineType, - focal_plane_image_series_id: int + focal_plane_image_series_id: int, + cell_specimen_id: int = None, + validate_polys: bool = True, ) -> Tuple[ - Dict[str, Union[NicePathType, str]], + Dict[str, Union[NicePathType, str]], Dict[str, Union[NicePathType, str]] ]: """ Return the pia and white matter surface drawings for this image series """ query = f""" - select + select distinct polygon.path as path, - label.name as name - from specimens sp - join specimens spp on spp.id = sp.parent_id - join image_series imser on imser.specimen_id = spp.id + label.name as name, + polygon.id as polygon_id + from image_series imser join sub_images si on si.image_series_id = imser.id join images im on im.id = si.image_id join treatments tm on tm.id = im.treatment_id join avg_graphic_objects layer on layer.sub_image_id = si.id join avg_graphic_objects polygon on polygon.parent_id = layer.id join avg_group_labels label on label.id = layer.group_label_id + JOIN biospecimen_polygons bsp ON bsp.polygon_id = polygon.id where imser.id = {focal_plane_image_series_id} and label.name in ('Pia', 'White Matter') and tm.name = 'Biocytin' """ - results = {} - for item in query_engine(query): - results[item["name"]] = { - "name": item["name"], - "path": ensure_path(item["path"]) - } - return results.get("Pia"), results.get("White Matter") - + if cell_specimen_id is not None: + query += f"and bsp.biospecimen_id = {cell_specimen_id}" + results = query_engine(query) + validate_fcn = lambda path: (len([x for x in path if not x==[0,0]]) > 1) + results_valid = sort_lims_graphics_objects( + results, validate_fcn=validate_fcn, validate_polys=validate_polys + ) + results_dict = {} + for x in results_valid: + results_dict[x['name']] = x + return results_dict.get("Pia"), results_dict.get("White Matter") + def query_for_images( - query_engine: QueryEngineType, + query_engine: QueryEngineType, focal_plane_image_series_id: int, output_dir: str ) -> List[Dict[str, str]]: - """ Return Biocytin and DAPI images associated with a focal plane image + """ Return Biocytin and DAPI images associated with a focal plane image series """ query = f""" - select - im.jp2, + select + im.jp2, sl.storage_directory, tm.name - from sub_images si - join images im on im.id = si.image_id - join slides sl on sl.id = im.slide_id + from sub_images si + join images im on im.id = si.image_id + join slides sl on sl.id = im.slide_id join treatments tm on tm.id = im.treatment_id - where + where image_series_id = {focal_plane_image_series_id} and tm.name in ('Biocytin', 'DAPI') """ @@ -201,7 +218,7 @@ def query_for_image_dims( """ query = f""" - select + select im.height as height, im.width as width from specimens sp @@ -224,23 +241,23 @@ def get_inputs_from_lims( database: str, user: str, password: str, - imser_id: int, + imser_id: int, image_output_root: Optional[str] ): """ Utility for building module inputs from a direct LIMS query """ engine = partial( - lu.query, - host=host, - port=port, - database=database, - user=user, + lu.query, + host=host, + port=port, + database=database, + user=user, password=password ) - layer_polygons = query_for_layer_polygons(engine, imser_id) - pia_surface, wm_surface = query_for_cortical_surfaces(engine, imser_id) + layer_polygons = query_for_layer_polygons(engine, imser_id, validate_polys=True) + pia_surface, wm_surface = query_for_cortical_surfaces(engine, imser_id, validate_polys=False) results = { "layer_polygons": layer_polygons, @@ -281,14 +298,14 @@ class PostgresInputConfigSchema(DefaultSchema): ) port = Int( description="", - required=False, + required=False, default=5432 ) class FromLimsSchema(PostgresInputConfigSchema): - """The parameters required to query LIMS for a set of cortical layer - polygons and cortical surface boundaries. + """The parameters required to query LIMS for a set of cortical layer + polygons and cortical surface boundaries. """ focal_plane_image_series_id = Int( From 1111cb562156ec617ae10a4756878eca0e715d20 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Thu, 13 Apr 2023 16:51:52 -0700 Subject: [PATCH 13/13] soma layer depths module --- neuron_morphology/soma_layer_depths.py | 412 +++++++++++++++++++++++++ 1 file changed, 412 insertions(+) create mode 100644 neuron_morphology/soma_layer_depths.py diff --git a/neuron_morphology/soma_layer_depths.py b/neuron_morphology/soma_layer_depths.py new file mode 100644 index 00000000..7d022295 --- /dev/null +++ b/neuron_morphology/soma_layer_depths.py @@ -0,0 +1,412 @@ +import logging +from functools import partial +from itertools import chain +from multiprocessing import Pool + +import numpy as np +import pandas as pd +from shapely.geometry import Polygon, Point, LineString +from scipy.interpolate import CloughTocher2DInterpolator + +import argschema as ags +from allensdk.internal.core import lims_utilities as lu +from neuron_morphology.lims_apical_queries import get_data +from neuron_morphology.transforms.pia_wm_streamlines.calculate_pia_wm_streamlines import generate_laplace_field +from neuron_morphology.transforms.geometry import get_vertices_from_two_lines +from neuron_morphology.snap_polygons.__main__ import run_snap_polygons, InputParameters +from neuron_morphology.snap_polygons._from_lims import query_for_layer_polygons, query_for_cortical_surfaces +from neuron_morphology.snap_polygons.types import ensure_linestring, ensure_path +from neuron_morphology.snap_polygons.geometries import Geometries, select_largest_subpolygon +from neuron_morphology.features.layer.reference_layer_depths import DEFAULT_HUMAN_MTG_REFERENCE_LAYER_DEPTHS, DEFAULT_MOUSE_REFERENCE_LAYER_DEPTHS +from neuron_morphology.layered_point_depths.__main__ import step_from_node + +logger = logging.getLogger(__name__) + +WELL_KNOWN_REFERENCE_LAYER_DEPTHS = { + "human": DEFAULT_HUMAN_MTG_REFERENCE_LAYER_DEPTHS, + "mouse": DEFAULT_MOUSE_REFERENCE_LAYER_DEPTHS, +} +class LayerDepthError(Exception): + pass + +def get_cell_soma_data(cell_specimen_ids): + # based on query in lims_apical_queries but remove requirement of reconstruction + # there is probably a better way to do this, and should be in neuron_morphology + ids_str = ', '.join([str(sid) for sid in cell_specimen_ids]) + query_for_soma = f""" + SELECT DISTINCT sp.id as specimen_id, 'null', layert.name as path_type, poly.path, sc.resolution, 'null', 'null' + FROM specimens sp + JOIN biospecimen_polygons AS bsp ON bsp.biospecimen_id=sp.id + JOIN avg_graphic_objects poly ON poly.id=bsp.polygon_id + JOIN avg_graphic_objects layer ON layer.id=poly.parent_id + JOIN avg_group_labels layert ON layert.id=layer.group_label_id + AND layert.prevent_missing_polygon_structure=false + JOIN sub_images AS si ON si.id=layer.sub_image_id + AND si.failed=false + JOIN images AS im ON im.id=si.image_id + JOIN slides AS s ON s.id=im.slide_id + JOIN scans AS sc ON sc.slide_id=s.id + AND sc.superseded=false + JOIN treatments t ON t.id = im.treatment_id AND t.id = 300080909 --Why? + WHERE sp.id IN ({ids_str}) + ORDER BY sp.id + """ + # all results returned as 'invalid_data' with only soma coords and resolution + _, cell_data = get_data(query_for_soma) + soma_centers = {k: cell_data[k]["soma_center"] for k in cell_specimen_ids + if cell_data[k]["soma_center"] is not None} + resolution = cell_data[int(cell_specimen_ids[1])]["resolution"] + return soma_centers, resolution + +def lims_slice_cells_info(query_engine, specimen_ids): + ids_str = ", ".join([str(s) for s in specimen_ids]) + query = f""" + select + sp.id as specimen_id, + slice.id as slice_id, + image_series.id as image_series_id + from specimens sp + join specimens slice on slice.id = sp.parent_id + join image_series on image_series.specimen_id=slice.id + where image_series.type='FocalPlaneImageSeries' + and sp.id in ({ids_str}) + """ + results = query_engine(query) + return results + +def query_for_soma_centers(query_engine, cell_specimen_ids): + """Return soma center coordinates for a list of specimens in a single image, + converted to microns (along with the resolution). + """ + ids_str = ', '.join([str(sid) for sid in cell_specimen_ids]) + + query = f""" + select distinct + polygon.path as path, + label.name as name, + sc.resolution as resolution, + bsp.biospecimen_id as specimen_id + from specimens sp + join specimens spp on spp.id = sp.parent_id + join image_series imser on imser.specimen_id = spp.id + join sub_images si on si.image_series_id = imser.id + join images im on im.id = si.image_id + join slides s on s.id=im.slide_id + join scans sc on sc.slide_id=s.id + join treatments tm on tm.id = im.treatment_id + join avg_graphic_objects layer on layer.sub_image_id = si.id + join avg_graphic_objects polygon on polygon.parent_id = layer.id + join avg_group_labels label on label.id = layer.group_label_id + join biospecimen_polygons bsp on bsp.polygon_id = polygon.id + where + label.name = 'Soma' + and tm.name = 'Biocytin' + and bsp.biospecimen_id in ({ids_str}) + """ + results = query_engine(query) + soma_centers = {} + if len(results)==0: + raise ValueError('No soma centers found.') + for item in results: + resolution = item["resolution"] + path = ensure_path(item["path"]) + soma_centers[item["specimen_id"]] = np.asarray(path).mean(axis=0)*resolution + return soma_centers, resolution + +def trim_layers(layer_polygons, pia_surface, wm_surface, + multipolygon_error_threshold=100, **kwargs): + """Follow the snap_polygons approach to trim a list of layer polygons to lie + within the bounds of a pia/wm-defined region.""" + if pia_surface is None or wm_surface is None: + raise ValueError("Missing pia/wm.") + geometries = Geometries() + geometries.register_polygons(layer_polygons) + geometries.register_surface("pia", pia_surface["path"]) + geometries.register_surface("wm", wm_surface["path"]) + pia_wm_vertices = get_vertices_from_two_lines( + pia_surface['path'], wm_surface['path']) + bounds = Polygon(pia_wm_vertices) + + multipolygon_resolver = partial( + select_largest_subpolygon, + error_threshold=multipolygon_error_threshold + ) + + result_geos = geometries.cut(bounds, multipolygon_resolver=multipolygon_resolver) + return result_geos.to_json() + +def layer_info_from_snap_polygons_output(output, resolution=1): + layers = {} + pia_path = None + wm_path = None + for polygon in output["polygons"]: + layers[polygon['name']] = {'bounds': Polygon(resolution*np.array(polygon['path']))} + for surface in output["surfaces"]: + name = surface['name'] + path = list(resolution*np.array(surface['path'])) + if name=='pia': + pia_path = path + elif name=='wm': + wm_path = path + else: + path = LineString(path) + layer, side = name.split('_') + layers[layer][f"{side}_surface"] = path + return layers, pia_path, wm_path + +def get_missing_layer_info(layers, species, infer_missing_layers=True): + ref_layer_depths = WELL_KNOWN_REFERENCE_LAYER_DEPTHS[species].copy() + # don't want to include wm as a layer! + ref_layer_depths.pop('wm') + all_layers_ordered = sorted(ref_layer_depths.keys()) + complete_layers = sorted(( + layer.replace("Layer", '') for layer, layer_poly in layers.items() + if 'pia_surface' in layer_poly and 'wm_surface' in layer_poly + )) + if not complete_layers: + raise LayerDepthError("No layer boundaries found.") + first = complete_layers[0] + last = complete_layers[-1] + top_path = layers[f"Layer{first}"]['pia_surface'] + top_path = list(top_path.coords) + bottom_path = layers[f"Layer{last}"]['wm_surface'] + bottom_path = list(bottom_path.coords) + missing_above = all_layers_ordered[:all_layers_ordered.index(first)] + missing_below = all_layers_ordered[all_layers_ordered.index(last)+1:] + pia_extra_dist = wm_extra_dist = 0 + if len(missing_above) > 0: + if infer_missing_layers: + pia_extra_dist = sum(ref_layer_depths[layer].thickness for layer in missing_above) + else: + pia_extra_dist = np.nan + wm_extra_dist = sum(ref_layer_depths[layer].thickness for layer in missing_below) + if len(missing_below) > 0: + if infer_missing_layers: + wm_extra_dist = sum(ref_layer_depths[layer].thickness for layer in missing_below) + else: + wm_extra_dist = np.nan + return top_path, bottom_path, pia_extra_dist, wm_extra_dist + +def get_layer_for_point(point, layer_polys): + in_layer = [ + layer for layer in layer_polys if + layer_polys[layer]['bounds'].intersects(Point(*point)) + # checks for common boundary or interior + ] + + if len(in_layer) == 0: + raise LayerDepthError("Point not found in any layer") + elif len(in_layer) == 1: + layer_name = in_layer[0] + else: + # overlap means point is likely on a boundary + # choose upper layer, avoiding L1 + for layer in sorted(in_layer): + if not layer=="Layer1": + layer_name = layer + logger.warning(f"Overlapping layers: {in_layer}. Choosing {layer_name}") + layer_poly = layer_polys[layer_name] + return layer_name, layer_poly + +def get_layer_depths(point, layer_polys, pia_path, wm_path, depth_interp, dx_interp, dy_interp, + step_size=1.0, max_iter=1000, + pia_extra_dist=0, wm_extra_dist=0): + + def dist_to_boundary(boundary_path, direction): + try: + _, dist = step_from_node( + point, depth_interp, dx_interp, dy_interp, + boundary_path, direction*step_size, max_iter, adaptive_scale=1 + ) + except ValueError as e: + logger.warning(e) + dist = np.nan + return dist + + layer_name, layer_poly = get_layer_for_point(point, layer_polys) + if pia_path is not None: + pia_path = ensure_linestring(pia_path) + pia_distance = dist_to_boundary(pia_path, 1) + else: + pia_distance = np.nan + if wm_path is not None: + wm_path = ensure_linestring(wm_path) + wm_distance = dist_to_boundary(wm_path, -1) + else: + wm_distance = np.nan + if 'pia_surface' in layer_poly: + pia_side_dist = dist_to_boundary(layer_poly['pia_surface'], 1) + else: + pia_side_dist = dist_to_boundary(layer_poly['bounds'].boundary, 1) + if 'wm_surface' in layer_poly: + wm_side_dist = dist_to_boundary(layer_poly['wm_surface'], -1) + else: + wm_side_dist = dist_to_boundary(layer_poly['bounds'].boundary, -1) + pia_distance += pia_extra_dist + wm_distance += wm_extra_dist + + layer_thickness = wm_side_dist + pia_side_dist + cortex_thickness = pia_distance + wm_distance + out = { + 'layer_depth': pia_side_dist, + 'layer_thickness': layer_thickness, + 'normalized_layer_depth': pia_side_dist/layer_thickness, + 'normalized_depth': pia_distance/cortex_thickness, #vs depth_interp(point)? + 'absolute_depth': pia_distance, + 'cortex_thickness': cortex_thickness, + 'wm_distance': wm_distance, + 'layer':layer_name, + } + return out + +def resample_line(coords, distance_delta=50): + line = LineString(coords) + distances = np.arange(0, line.length, distance_delta) + points = [line.interpolate(distance) for distance in distances] + [line.boundary[1]] + line_coords = [point.coords[0] for point in points] + return line_coords + +def get_depths_from_processed_layers(layers, pia_path, wm_path, soma_centers, + step_size=2.0, max_iter=1000, + species=None, orient_by_layer_bounds=False, + infer_missing_layers=True + ): + """Get depths (within-layer and absolute) for a list of soma positions, + given processed layer and pia/wm records. + + The key difference from the layered_point_depths approach is to not require + annotated pia_surface/wm_surface for each layer. Instead you just need the + pia/wm for up/down orientation and use the whole polygon as a boundary.""" + errors = [] + try: + top_path = pia_path + bottom_path = wm_path + pia_extra_dist = wm_extra_dist = 0 + if (pia_path is None) or (wm_path is None): + if not orient_by_layer_bounds: + raise ValueError("Pia or WM is missing, can't find depths.") + top_path, bottom_path, pia_extra_dist, wm_extra_dist = get_missing_layer_info( + layers, species, infer_missing_layers=infer_missing_layers) + + top_path = resample_line(top_path) + bottom_path = resample_line(bottom_path) + + (_, _, _, mesh_coords, mesh_values, mesh_gradients) = generate_laplace_field( + top_path, + bottom_path, + ) + interp = CloughTocher2DInterpolator + depth_interp = interp(mesh_coords, mesh_values) + dx_interp = interp(mesh_coords, mesh_gradients[:,0]) + dy_interp = interp(mesh_coords, mesh_gradients[:,1]) + except LayerDepthError as exc: + top_path = bottom_path = None + pia_extra_dist = wm_extra_dist = np.nan + depth_interp = dx_interp = dy_interp = lambda x: np.nan + logger.error(exc) + errors.append(str(exc)) + + outputs = {} + cell_errors = {} + for name, point in soma_centers.items(): + try: + outputs[name] = get_layer_depths( + point, layers, top_path, bottom_path, depth_interp, dx_interp, dy_interp, + step_size=step_size, max_iter=max_iter, + pia_extra_dist=pia_extra_dist, wm_extra_dist=wm_extra_dist + ) + except (LayerDepthError,) as exc: + error = f"Failure getting depth info for cell {name}: {exc}" + logger.error(error) + cell_errors[name] = str(exc) + if ((len(layers) >= 3) | ((pia_path is not None) & (wm_path is not None))): + if len(soma_centers) == len(cell_errors): + raise ValueError(f"All cells in slice failed unexpectedly: {cell_errors}") + return outputs, errors, cell_errors + +def run_layer_depths(image_series_id, cell_specimen_ids, species=None, grouped_cells=True, + use_snap_polygons_boundaries=False, infer_missing_layers=True): + if grouped_cells: + try: + layer_polygons = query_for_layer_polygons(lu.query, image_series_id) + pia_surface, wm_surface = query_for_cortical_surfaces(lu.query, image_series_id) + except ValueError: + # error may be from multiple pia/wm drawings, running cells individually + # lets them each find the correct drawing + slice_records = chain.from_iterable( + [run_layer_depths(image_series_id, [specimen_id], grouped_cells=False) + for specimen_id in cell_specimen_ids] + ) + return slice_records + else: + layer_polygons = query_for_layer_polygons(lu.query, image_series_id, validate_polys=False) + pia_surface, wm_surface = query_for_cortical_surfaces( + lu.query, image_series_id, cell_specimen_id=cell_specimen_ids[0] + ) + + try: + soma_centers, resolution = query_for_soma_centers(lu.query, cell_specimen_ids) + if use_snap_polygons_boundaries: + input_data=dict(layer_polygons=layer_polygons, pia_surface=pia_surface, wm_surface=wm_surface) + args = InputParameters().load(input_data) + args.pop('log_level') + output = run_snap_polygons(**args) + else: + # trim layers to pia/wm bounds (subset of snap_polygons process) + # helps run streamlines in L1, in case boundary extends slightly above pia + output = trim_layers(layer_polygons, pia_surface, wm_surface) + layers, pia_path, wm_path = layer_info_from_snap_polygons_output(output, resolution) + + slice_results, slice_errors, cell_errors = get_depths_from_processed_layers( + layers, pia_path, wm_path, soma_centers, species=species, max_iter=4000, + orient_by_layer_bounds=use_snap_polygons_boundaries, + infer_missing_layers=infer_missing_layers) + + slice_records = [dict(specimen_id=specimen_id, image_series_id=image_series_id, + errors=slice_errors, **results) + for specimen_id, results in slice_results.items()] + missing_records = [dict(specimen_id=specimen_id, image_series_id=image_series_id, + errors="Missing soma center") + for specimen_id in cell_specimen_ids if specimen_id not in soma_centers] + error_records = [dict(specimen_id=specimen_id, image_series_id=image_series_id, errors=err) + for specimen_id, err in cell_errors.items()] + slice_records += missing_records + error_records + except ValueError as exc: + logging.warning(f'Failed to get depths: image series {image_series_id}. {exc}') + slice_records = [dict(specimen_id=specimen_id, image_series_id=image_series_id, + errors=str(exc)) + for specimen_id in cell_specimen_ids] + except Exception as exc: + logging.exception(f'Failed to get depths: image series {image_series_id}') + slice_records = [dict(specimen_id=specimen_id, image_series_id=image_series_id, + errors=repr(exc)) + for specimen_id in cell_specimen_ids] + return slice_records + +def main(input_file, output_file, **kwargs): + cells = pd.read_csv(input_file, index_col=0).index + + df = pd.DataFrame.from_records(lims_slice_cells_info(lu.query, cells)) + cell_groups = df.groupby('image_series_id')['specimen_id'].apply(list) + image_series_ids = cell_groups.index.values + cell_specimen_ids = cell_groups.values + + pool = Pool() + fcn = partial(run_layer_depths, **kwargs) + results = pool.starmap(fcn, zip(image_series_ids, cell_specimen_ids)) + + depth_df = pd.DataFrame.from_records(chain(*results)) + depth_df.to_csv(output_file, index=False) + +class SomaLayerDepthSchema(ags.ArgSchema): + input_file = ags.fields.InputFile( + description="input cell list as a csv with the specimen ids in the first column" + ) + output_file = ags.fields.OutputFile(default="soma_depths.csv") + use_snap_polygons_boundaries = ags.fields.Boolean(default=False) + +if __name__ == "__main__": + module = ags.ArgSchemaParser(schema_type=SomaLayerDepthSchema) + module.args.pop("log_level") + # logger.setLevel(module.args.pop("log_level")) + main(**module.args)