diff --git a/doc_template/install.rst b/doc_template/install.rst
index d40ef71b..a289e6af 100644
--- a/doc_template/install.rst
+++ b/doc_template/install.rst
@@ -14,7 +14,7 @@ Two popular tools for managing Python environments are `anaconda <https://anacon
 
 .. code-block:: bash
 
-    conda create -y --name environment-name python=3.6
+    conda create -y --name environment-name python=3.7
     conda activate environment-name
 
 and using venv:
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/__main__.py b/neuron_morphology/snap_polygons/__main__.py
index f5d0c873..62f75f82 100644
--- a/neuron_morphology/snap_polygons/__main__.py
+++ b/neuron_morphology/snap_polygons/__main__.py
@@ -33,26 +33,29 @@ 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.")
     # 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 +72,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)        
 
@@ -83,20 +86,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"))
diff --git a/neuron_morphology/snap_polygons/_from_lims.py b/neuron_morphology/snap_polygons/_from_lims.py
index 4ad3e8a1..68c5291c 100644
--- a/neuron_morphology/snap_polygons/_from_lims.py
+++ b/neuron_morphology/snap_polygons/_from_lims.py
@@ -4,25 +4,26 @@
 Example Usage
 -------------
 python -m neuron_morphology.snap_polygons
-    --host <lims host> 
-    --port <lims port> 
-    --user <username> 
-    --password <password> 
-    --database <lims db> 
+    --host <lims host>
+    --port <lims port>
+    --user <username>
+    --password <password>
+    --database <lims db>
     --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 marshmallow as mm
-
+import numpy as np
+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 (
@@ -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
@@ -58,13 +59,11 @@ def query_for_layer_polygons(
     """
 
     query = f"""
-        select
+        select distinct
             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,50 +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}"
                     )
-            else:
-                warnings.warn(
-                    f"found multiple polygon records for {name} "
-                    "(identical paths)"
-                )
-
-        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:
@@ -127,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["Pia"], results["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')
     """
@@ -206,7 +218,7 @@ def query_for_image_dims(
     """
 
     query = f"""
-        select 
+        select
             im.height as height,
             im.width as width
         from specimens sp
@@ -229,30 +241,28 @@ 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)
-    image_width, image_height = query_for_image_dims(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,
         "pia_surface": pia_surface,
         "wm_surface": wm_surface,
-        "image_dimensions": {"width": image_width, "height": image_height}
     }
 
     if image_output_root is not None:
@@ -262,37 +272,40 @@ def get_inputs_from_lims(
     return results
 
 
-class PostgresInputConfigSchema(mm.Schema):
+class PostgresInputConfigSchema(DefaultSchema):
     """The parameters required to query a postgres database.
     """
 
     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="",
-        required=False,  # seems not to get hydrated from the default
+        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(
@@ -307,20 +320,20 @@ 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
         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
-            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")
         )
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=(
diff --git a/neuron_morphology/snap_polygons/geometries.py b/neuron_morphology/snap_polygons/geometries.py
index 6c6f4a1f..c43ab55c 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
 
@@ -512,14 +514,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 masks (of the same shape). 0 indicates
+        absence 
 
     Returns
     -------
@@ -576,8 +577,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 +604,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 +637,7 @@ def shared_faces(poly: Polygon, others: Iterable[Polygon]) -> LineString:
 
     Returns
     -------
-    LineString representing the shared face
+    LineString representing the shared face, or None if not found
     """
 
     faces_list = []
@@ -644,8 +651,15 @@ def shared_faces(poly: Polygon, others: Iterable[Polygon]) -> LineString:
         faces = shapely.ops.linemerge(backward)
 
         if not faces.is_empty:
+            if faces.type == 'MultiLineString':
+                try:
+                    faces = shapely.ops.linemerge([*faces, LineString(shapely.ops.nearest_points(*faces))])
+                except:
+                    raise ValueError("Can't find continuous shared face between polygons.")
             faces_list.append(faces)
 
+    if not faces_list:
+        return None
     merged_faces = shapely.ops.linemerge(faces_list)
     coordinates = list(merged_faces.coords)
     shared_line = ensure_linestring(coordinates)
diff --git a/neuron_morphology/snap_polygons/types.py b/neuron_morphology/snap_polygons/types.py
index ead97d4b..2a9fa289 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.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
 
diff --git a/tests/snap_polygons/test_geometries.py b/tests/snap_polygons/test_geometries.py
index 75392541..a1f393a3 100644
--- a/tests/snap_polygons/test_geometries.py
+++ b/tests/snap_polygons/test_geometries.py
@@ -208,7 +208,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)])