Repository: sanderslab/magellanmapper (language: python)

Issue:
Back up blob archives and fix saving blob columns
We're splitting out the blob changes from the classifier branch (#199) here for modularity and to integrate them sooner.

- Blob archives are backed up before overwriting
- Fixed saving additional blob columns
- Fixed blob accessors/setters
- Fixed blob detections in the GUI with no selected channels or loaded blobs but none in the ROI
- Fixed exporting ROIs
- Reduced debug messages during blob detection


Discussion / hints:

Fix patch:
diff --git a/bin/sample_cmds_bash.ipynb b/bin/sample_cmds_bash.ipynb
index e42b1adc2..a3dd50561 100644
--- a/bin/sample_cmds_bash.ipynb
+++ b/bin/sample_cmds_bash.ipynb
@@ -315,8 +315,97 @@
   },
   {
    "cell_type": "markdown",
-   "id": "26f5d7fc-e42d-41ed-8b7a-f7cac9ae0426",
+   "id": "195fb19a-fe6c-4ea6-afe1-0070788da70d",
+   "metadata": {
+    "tags": []
+   },
+   "source": [
+    "### Register Atlas to Image\n",
+    "\n",
+    "Now that we have downsampled our sample image, we can register an atlas to it. Image registration shifts and morphs a \"moving\" image to align it to a \"fixed\" image. To preserve the morphology of the sample image, we make the atlas the \"moving\" image that moves to fit the sample.\n",
+    "\n",
+    "Let's set up the atlas. The registration toolkit supports 2D to 2D or 3D to 3D image registration. Following the conventions for the Allen Insitute atlases, MagellanMapper will look for an atlas intensity image named `atlasVolume.<ext>`, where `ext` can be any extension supported by ITK, such as MHD, NRRD, and NIFTI. The annotated image should be named `annotation.<ext>`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "59c949bd-2cdb-451f-a842-de9b7afd8052",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# path to atlas directory, assumed to have an atlasVolume and annotation file\n",
+    "atlas_dir=\"/path/to/your/atlas\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "51fc570f-935c-4004-a9d8-40f3e743ac67",
+   "metadata": {},
+   "source": [
+    "We now register the images by listing the moving image followed by the fixed image.\n",
+    "- `single` registration registers single images as opposed to a group of images\n",
+    "- `--prefix \"$img\"` outputs image filenames formatted according to the original rather than the resized image"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "74a61265-63c5-4ea8-bb5b-d645e9bf9bd0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "./run.py --img \"$img_resized\" \"$atlas_dir\" --register single --prefix \"$img\" -v"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "33f0954e-b4b1-436b-96c3-90c2f1d7d47d",
+   "metadata": {},
+   "source": [
+    "Atlases are often need reorientation to fit a given sample. Here, we add several example options:\n",
+    "- `--transform rotate=2` rotates the atas by 90 deg x2 (= 180 deg) counter-clockwise\n",
+    "- `--atlas_profile abaccfv3` uses settings in an atlas profile designed for the Allen CCFv3 atlas\n",
+    "- `--channel 1` uses channel 1 (ie the 2nd channel) of the sample for registration"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bbbfbbff-5065-4904-93cc-b6ebddd047a4",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "./run.py --img \"$img_resized\" \"$atlas_dir\" --register single -v \\\n",
+    "  --transform rotate=2 --atlas_profile abaccfv3 --prefix \"$img\" --channel 1"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7191b3e8-5a55-4cef-bd0d-c4c3813a5d22",
+   "metadata": {},
+   "source": [
+    "To view the registered atlas on the sample, we use the base name of the image and add \"registered image suffixes\" to specify which registered images to load. Registered images use the base image name and add a suffix depending on the image type, such as `_atlasVolume.mhd` and `_annotation.mhd` for the registered atlas intensity and annotation images, respectively. `_exp.mhd` is the same as the sample image.\n",
+    "- Registered images suffixes are given as: `--reg_suffixes atlas=<suffix> annotation=<suffix> borders=<suffix> fixed_mask=<suffix> moving_mask=<suffix>`, where the type can be omitted if in the given order\n",
+    "- `--roi_profile atlas` can be added to use a grayscale color profile for the intensity image"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8cc8cd6d-eb09-4580-acb7-14e375d08f41",
    "metadata": {},
+   "outputs": [],
+   "source": [
+    "./run.py --img \"$img\" --reg_suffixes exp.mhd annotation.mhd"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "26f5d7fc-e42d-41ed-8b7a-f7cac9ae0426",
+   "metadata": {
+    "tags": []
+   },
    "source": [
     "## Cell Detection\n",
     "\n",
@@ -403,93 +492,110 @@
   },
   {
    "cell_type": "markdown",
-   "id": "195fb19a-fe6c-4ea6-afe1-0070788da70d",
+   "id": "92bedce5-8fec-4310-8d18-f56671eb6f3b",
    "metadata": {},
    "source": [
-    "### Register Atlas to Image\n",
+    "### Classify cells\n",
     "\n",
-    "Now that we have downsampled our sample image, we can register an atlas to it. Image registration shifts and morphs a \"moving\" image to align it to a \"fixed\" image. To preserve the morphology of the sample image, we make the atlas the \"moving\" image that moves to fit the sample.\n",
+    "The blob detector identifies objects based on contrast to surrounding areas but cannot distinguish between types of objects. A classifier can help to filter out unwanted detections, such as background noise, or identify different cell types.\n",
     "\n",
-    "Let's set up the atlas. The registration toolkit supports 2D to 2D or 3D to 3D image registration. Following the conventions for the Allen Insitute atlases, MagellanMapper will look for an atlas intensity image named `atlasVolume.<ext>`, where `ext` can be any extension supported by ITK, such as MHD, NRRD, and NIFTI. The annotated image should be named `annotation.<ext>`."
+    "#### Setup\n",
+    "\n",
+    "MagellanMapper uses additional dependencies for the classifier, which can be installed with the `classifier` group:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "59c949bd-2cdb-451f-a842-de9b7afd8052",
+   "id": "da6e3ee9-e941-4367-a636-46e8615bdbbe",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# path to atlas directory, assumed to have an atlasVolume and annotation file\n",
-    "atlas_dir=\"/path/to/your/atlas\""
+    "pip install -e .[classifier]"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "51fc570f-935c-4004-a9d8-40f3e743ac67",
+   "id": "dfa83fc8-792c-4825-aeaf-c7db5c3cfb47",
    "metadata": {},
    "source": [
-    "We now register the images by listing the moving image followed by the fixed image.\n",
-    "- `single` registration registers single images as opposed to a group of images\n",
-    "- `--prefix \"$img\"` outputs image filenames formatted according to the original rather than the resized image"
+    "#### Train classifier\n",
+    "\n",
+    "The classifier trains on labeled images as ground truth. The blob detector can identify objects in your image, which you can label using the annotation tools in the ROI Editor. The classifier learns from these annotations and the images patches around them using a classifier network in Pytorch.\n",
+    "\n",
+    "TODO\n",
+    "\n",
+    "#### Apply classifier\n",
+    "\n",
+    "The trained classifier can be applied to new images through the GUI or CLI.\n",
+    "\n",
+    "##### GUI\n",
+    "\n",
+    "1. In the \"Detect\" panel > \"Classifier model\" field, browse to the trained `.h5` classifier model\n",
+    "1. Check the desired channels in the \"Chl\" area.\n",
+    "1. Presse the \"Detect\" button to detect blobs, which will automatically classify them afterward.\n",
+    "\n",
+    "Blobs are colored by their class in the ROI Editor. These classes also appear in the \"confirmed\" column of the blobs table.\n",
+    "\n",
+    "##### CLI\n",
+    "\n",
+    "To detect and classify blobs, specify the path to your classifier model and run the `detect` task:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "74a61265-63c5-4ea8-bb5b-d645e9bf9bd0",
+   "id": "11156681-0095-47d3-a7e9-07731c998329",
    "metadata": {},
    "outputs": [],
    "source": [
-    "./run.py --img \"$img_resized\" \"$atlas_dir\" --register single --prefix \"$img\" -v"
+    "roi_profile=\"lightsheet\"\n",
+    "model_path=\"/path/to/your/model.h5\"\n",
+    "mm \"$img\" --roi_profile \"$roi_profile\" --classifier model=\"$model_path\" --proc detect"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "33f0954e-b4b1-436b-96c3-90c2f1d7d47d",
+   "id": "d5414d4b-8739-48d3-becf-c501a8b8ed5b",
    "metadata": {},
    "source": [
-    "Atlases are often need reorientation to fit a given sample. Here, we add several example options:\n",
-    "- `--transform rotate=2` rotates the atas by 90 deg x2 (= 180 deg) counter-clockwise\n",
-    "- `--atlas_profile abaccfv3` uses settings in an atlas profile designed for the Allen CCFv3 atlas\n",
-    "- `--channel 1` uses channel 1 (ie the 2nd channel) of the sample for registration"
+    "You can also classify previously detections, which allows repeatedly re-classifying blobs while tweaking the model:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "bbbfbbff-5065-4904-93cc-b6ebddd047a4",
+   "id": "07ad6d67-5a91-4b19-a936-728fa8c8e072",
    "metadata": {},
    "outputs": [],
    "source": [
-    "./run.py --img \"$img_resized\" \"$atlas_dir\" --register single -v \\\n",
-    "  --transform rotate=2 --atlas_profile abaccfv3 --prefix \"$img\" --channel 1"
+    "mm \"$img\" --classifier model=\"$model_path\" --proc classify --load blobs"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "7191b3e8-5a55-4cef-bd0d-c4c3813a5d22",
+   "id": "5362966b-234a-4380-8b73-b57ae646bc64",
    "metadata": {},
    "source": [
-    "To view the registered atlas on the sample, we use the base name of the image and add \"registered image suffixes\" to specify which registered images to load. Registered images use the base image name and add a suffix depending on the image type, such as `_atlasVolume.mhd` and `_annotation.mhd` for the registered atlas intensity and annotation images, respectively. `_exp.mhd` is the same as the sample image.\n",
-    "- Registered images suffixes are given as: `--reg_suffixes atlas=<suffix> annotation=<suffix> borders=<suffix> fixed_mask=<suffix> moving_mask=<suffix>`, where the type can be omitted if in the given order\n",
-    "- `--roi_profile atlas` can be added to use a grayscale color profile for the intensity image"
+    "Load the blobs to view their classifications in the GUI:"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "8cc8cd6d-eb09-4580-acb7-14e375d08f41",
+   "id": "c85448f0-040f-4800-9438-a2a4badd4695",
    "metadata": {},
    "outputs": [],
    "source": [
-    "./run.py --img \"$img\" --reg_suffixes exp.mhd annotation.mhd"
+    "mm \"$img\" --load blobs"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "d762cc46-a5b8-4e46-b925-c159762258e6",
-   "metadata": {},
+   "metadata": {
+    "tags": []
+   },
    "source": [
     "## 3D Visualization\n",
     "\n",
diff --git a/magmap/cv/classifier.py b/magmap/cv/classifier.py
new file mode 100644
index 000000000..111e28229
--- /dev/null
+++ b/magmap/cv/classifier.py
@@ -0,0 +1,350 @@
+"""Blob classifier."""
+from enum import Enum, auto
+from time import time
+from typing import Tuple, Sequence, Optional
+
+import numpy as np
+
+from magmap.cv import chunking, detector
+from magmap.io import cli
+from magmap.plot import plot_3d
+from magmap.settings import config
+
+_logger = config.logger.getChild(__name__)
+
+
+def extract_patches(roi: np.ndarray, blobs: np.ndarray, size: int = 16):
+    """Extract image patches for blobs.
+    
+    Patches are 2D, centered on each blob but offset by one pixel in width
+    and height for even-numbered patch dimensions.
+    
+    Args:
+        roi: Image region of interst as a 3/4D array (``z, y, x, [c]``).
+        blobs: 2D blobs array with each blob as a row in ``z, y, x, ...``.
+        size: Patch size as an int for both width and height; defaults to 16.
+
+    Returns:
+
+    """
+    # px backward from blob center
+    size_back = size // 2
+    # px forward
+    size_fwd = -(size // -2)
+    
+    patches = []
+    for blob in blobs[:, :3].astype(int):
+        # extract 2D patch around blob
+        blob = blob
+        z = blob[0]
+        y = blob[1]
+        x = blob[2]
+        patch = roi[z, y - size_back:y + size_fwd,
+                    x - size_back:x + size_back, ...]
+        
+        # normalize patch to itself
+        patch = patch / np.max(patch)
+        patches.append(patch)
+    
+    # combine patches and add a channel axis
+    x = np.stack(patches)
+    shape = list(x.shape)
+    shape.append(1)
+    x = x.reshape(shape)
+    
+    return x
+
+
+def classify_patches(model, x: np.ndarray, thresh: float = 0.5
+) -> Tuple[np.ndarray, np.ndarray]:
+    """Classify patches with a model.
+    
+    Args:
+        model: Keras model.
+        x: 2D array of image patches, each in ``y, x, c`` format.
+        thresh: Score threshold to classify as 1, otherwise 0. Defaults to 0.5.
+
+    Returns:
+        Tuple of:
+        - ``y_pred``: Integer array of class predictions.
+        - ``y_score``: Float array of raw prediction scores.
+
+    """
+    # calculate prediction scores and assign predictions based on threshold
+    y_score = model.predict(x).squeeze()
+    y_pred = (y_score > thresh).astype(int).squeeze()
+    
+    return y_pred, y_score
+
+
+def setup_classification_roi(
+        image5d: np.ndarray, subimg_offset: Sequence[int],
+        subimg_size: Sequence[int],
+        blobs: "detector.Blobs", patch_size: int, blobs_relative: bool
+) -> Tuple[np.ndarray, np.ndarray, Sequence[int]]:
+    """Set up ROI for blob classification.
+    
+    Sets up an image ROI with a border to provide uniformly sized image
+    patches around blobs, including blobs on the ROI edge. If the border
+    would exceed an edge of the image, a padding area is created where
+    blobs are excluded so that edge blobs can still have the full-sized
+    patches.
+    
+    Args:
+        subimg_offset: Subimage offset in ``z, y, x``.
+        subimg_size: Subimage size in ``z, y, x``.
+        blobs: Blobs instance.
+        patch_size: Patch size as an int for both width and height.
+    
+    Returns:
+        Tuple of:
+        - ``roi``: region of interest as ``z, y, x, [c]``
+        - ``blobs_roi_mask``: mask for ``blobs`` in the ROI.
+        - ``blobs_shift``: Offset of blobs relative to ``subimg_offset``.
+
+    """
+    
+    # reduce subimage size if it would exceed image boundaries
+    img_shape = image5d.shape[1:4]
+    border_far_roi = np.add(subimg_offset, subimg_size)
+    border_far_roi = np.where(
+        np.greater_equal(border_far_roi, img_shape), img_shape, border_far_roi)
+    subimg_size = border_far_roi - subimg_offset
+    border_far_roi = np.add(subimg_offset, subimg_size)
+    
+    # initialize size of ROI border along each axis and offset
+    border = (0, patch_size // 2, patch_size // 2)
+    border_offset = np.subtract(subimg_offset, border)
+    
+    # set ROI bounding box by defining opposite corners, preventing them from
+    # exceeding image boundaries
+    border_near = np.where(border_offset < 0, 0, border_offset)
+    border_far_full = border_far_roi + border
+    border_far = np.where(
+        border_far_full > img_shape, img_shape, border_far_full)
+    roi = plot_3d.prepare_subimg(
+        image5d, border_near, np.subtract(border_far, border_near))
+    
+    # blobs ROI defaults to ROI without border, but any ROI border truncation
+    # is converted to blob ROI padding to make up for the border loss
+    blobs_near = np.where(border_offset < 0, -border_offset, subimg_offset)
+    blobs_far = np.where(
+        border_far_full > img_shape,
+        np.multiply(img_shape, 2) - border_far_full, border_far_roi)
+    
+    # convert blob offsets to positions relative to ROI without border
+    blobs_rel_offset = np.subtract(
+        blobs_near, subimg_offset) if blobs_relative else blobs_near
+    blobs_size = np.subtract(blobs_far, blobs_near)
+    blobs_shift = np.subtract(subimg_offset, border_near)
+    blobs_roi, blobs_roi_mask = detector.get_blobs_in_roi(
+        blobs.blobs, blobs_rel_offset, blobs_size, reverse=False)
+    
+    return roi, blobs_roi_mask, blobs_shift
+
+
+def classify_blobs(
+        path: str, image5d: np.ndarray, subimg_offset: Sequence[int],
+        subimg_size: Sequence[int], channels: Sequence[int],
+        blobs: "detector.Blobs", patch_size: int = 16, blobs_relative=False
+) -> Tuple[np.ndarray, Sequence[int]]:
+    """Classify blobs based on surrounding image patches.
+    
+    Args:
+        path: Path to Keras model.
+        image5d: Image in ``t, z, y, x, [c]`` order.
+        subimg_offset: Subimage offset in ``z, y, x``.
+        subimg_size: Subimage size in ``z, y, x``.
+        channels: Sequence of channels in ``image5d``.
+        blobs: Blobs instance.
+        patch_size: Patch size as an int for both width and height; defaults
+            to 16.
+        blobs_relative: True to treat blob coordinates as relative to
+            ``subimg_offset``; defaults to False.
+        
+    Returns:
+        Tuple of:
+        - ``blobs_mask``: row mask for blobs in sub-image.
+        - ``classifications``: corresponding blob classifications.
+    
+    Raises:
+        `ModuleNotFoundError` if Tensorflow is not installed.
+
+    """
+    
+    # set up image and blobs ROIs
+    _logger.info("Setting up classification ROI")
+    roi_class, blobs_roi_mask, blobs_shift = setup_classification_roi(
+        image5d, subimg_offset, subimg_size, blobs, patch_size, blobs_relative)
+    
+    # load model with Keras
+    try:
+        from tensorflow.keras.models import load_model
+    except ModuleNotFoundError as e:
+        raise ModuleNotFoundError(
+            "Tensorflow is required for classification. Please install, eg "
+            "with 'pip install tensorflow'.")
+    model = load_model(path)
+
+    for chl in channels:
+        _logger.info("Classifying blobs in channel: %s", chl)
+        
+        # get blobs in channel for mask and to ensure blobs are present
+        blobs_chl, blobs_chl_mask = blobs.blobs_in_channel(
+            blobs.blobs, chl, True)
+        
+        # combine masks to set blobs array view without intermediate copy
+        blobs_mask = np.logical_and(blobs_roi_mask, blobs_chl_mask)
+        if np.sum(blobs_mask) < 1: continue
+        
+        # shift blob coordinates by border or padding
+        blobs_chl = np.add(blobs.blobs[blobs_mask, :3], blobs_shift)
+        if not blobs_relative:
+            # convert absolute to relative coordinates
+            blobs_chl -= subimg_offset
+        
+        # classify blobs based on surrounding image patches
+        roi_chl = (roi_class if roi_class.ndim < 4
+                   else roi_class[..., chl])
+        patches = extract_patches(roi_chl, blobs_chl, patch_size)
+        y_pred, y_score = classify_patches(model, patches)
+        blobs.set_blob_confirmed(blobs.blobs, y_pred, mask=blobs_mask)
+        # blobs.set_blob_truth(blobs.blobs, y_score)
+    
+    classifications = blobs.get_blob_confirmed(blobs.blobs)[blobs_roi_mask]
+    return blobs_roi_mask, classifications
+
+
+class ClassifyImage:
+    """Convert a label to an edge with class methods as an encapsulated 
+    way to use in multiprocessing without requirement for global variables.
+
+    """
+    image5d = None
+    blobs = None
+    
+    @classmethod
+    def classify_whole_image(
+            cls, model_path: Optional[str] = None,
+            image5d: Optional[np.ndarray] = None,
+            channels: Optional[Sequence[int]] = None,
+            blobs: Optional["detector.Blobs"] = None, **kwargs):
+        """Classify blobs in the whole image through multiprocessing
+        
+        Args:
+            model_path: Path to Keras model. Defaults to None, in which case
+                :attr:`magmap.settings.config.classifier` is accessed.
+            image5d: Image in ``t, z, y, x, [c]`` order. Defaults to None, in which
+                case :attr:`magmap.settings.config.img5d.img` is accessed.
+            channels: Sequence of channels in ``image5d``. Defaults to None, in
+                which case :attr:`magmap.settings.config.chanel` is accessed.
+            blobs: Blobs instance. Defaults to None, in which case
+                :attr:`magmap.settings.config.blobs` is accessed.
+            kwargs: Additional arguments to :meth:`classify_blobs`.
+        
+        Raises:
+            `FileNotFoundError` if a classifier model, image, or blobs are not
+            found.
+    
+        """
+        
+        if not model_path:
+            # get model path from config
+            model_path = config.classifier.model
+            if not model_path:
+                raise FileNotFoundError("No classifier model found")
+        
+        if image5d is None and config.img5d is not None:
+            # get main image from config
+            image5d = config.img5d.img
+            if image5d is None:
+                raise FileNotFoundError("No image found")
+        
+        if channels is None:
+            # set up channels from config
+            channels = plot_3d.setup_channels(
+                config.image5d, config.channel, 4)[1]
+        
+        if blobs is None:
+            # get blobs from config
+            blobs = config.blobs
+            if blobs is None:
+                raise FileNotFoundError("No blobs found")
+
+        is_fork = chunking.is_fork()
+        if is_fork:
+            # share large objects as class variables in forked processes
+            cls.image5d = image5d
+            cls.blobs = blobs
+        
+        # set up multiprocessing
+        start_time = time()
+        pool = chunking.get_mp_pool()
+        pool_results = []
+        
+        # chunk image simply by planes
+        step = 100
+        img_shape = image5d.shape[1:4]
+        size = list(img_shape)
+        size[0] = step
+        nsteps = -(img_shape[0] // -step)
+        
+        _logger.info("Classifying whole image in %s steps", nsteps)
+        for i in range(nsteps):
+            # increment offset z-plane by step size
+            offset = (i * step, 0, 0)
+            pool_results.append(
+                pool.apply_async(
+                    cls.classify_chunk,
+                    args=(model_path, offset, size, channels), kwds=kwargs))
+        
+        for result in pool_results:
+            offset, blobs_mask, predictions = result.get()
+            _logger.info(
+                "Classified blobs at offset %s of %s", offset, img_shape)
+            
+            # update blobs since modifications in separate processes are lost
+            blobs.set_blob_confirmed(blobs.blobs, predictions, blobs_mask)
+            print(blobs.blobs[blobs_mask][:5])
+        pool.close()
+        pool.join()
+
+        _logger.debug(
+            "Time elapsed to classify whole image: %s", time() - start_time)
+    
+    @classmethod
+    def classify_chunk(
+            cls, model_path: str, subimg_offset: Sequence[int],
+            subimg_size: Sequence[int], channels: Sequence[int], **kwargs):
+        """Classify blobs in an image chunk.
+        
+        Args:
+            model_path: Path to Keras model
+            subimg_offset: Subimage offset in ``z, y, x``.
+            subimg_size: Subimage size in ``z, y, x``.
+            channels: Sequence of channels in :attr:`image5d`.
+            kwargs: Additional arguments to :meth:`classify_blobs`.
+        
+        Returns:
+            Tuple of:
+            - ``subimg_offset``: the ``subimg_offset`` argument to track
+              during multiprocessing.
+            - ``blobs_mask``: row mask for blobs in sub-image.
+            - ``classifications``: corresponding blob classifications.
+    
+        """
+        
+        if cls.image5d is None:
+            # reopen main image in each spawned process
+            cli.process_cli_args()
+            cli.setup_image(config.filename)
+            cls.image5d = config.img5d.img
+            cls.blobs = config.blobs
+        
+        # classify blobs using model
+        blobs_mask, classifications = classify_blobs(
+            model_path, cls.image5d, subimg_offset, subimg_size, channels,
+            cls.blobs, **kwargs)
+        
+        return subimg_offset, blobs_mask, classifications
+        
diff --git a/magmap/cv/stack_detect.py b/magmap/cv/stack_detect.py
index 6743f34bf..dec77519b 100644
--- a/magmap/cv/stack_detect.py
+++ b/magmap/cv/stack_detect.py
@@ -8,12 +8,12 @@
 from enum import Enum
 import os
 from time import time
-from typing import NamedTuple, Sequence, TYPE_CHECKING
+from typing import NamedTuple, Sequence, TYPE_CHECKING, Optional, Tuple
 
 import numpy as np
 import pandas as pd
 
-from magmap.cv import chunking, colocalizer, detector, verifier
+from magmap.cv import chunking, classifier, colocalizer, detector, verifier
 from magmap.io import cli, df_io, importer, libmag, naming
 from magmap.plot import plot_3d
 from magmap.settings import config, roi_prof
@@ -320,9 +320,13 @@ def setup_blocks(
         exclude_border, tol, overlap_base, overlap, overlap_padding, max_pixels)
 
 
-def detect_blobs_blocks(filename_base, image5d, offset, size, channels,
-                        verify=False, save_dfs=True, full_roi=False,
-                        coloc=False):
+def detect_blobs_blocks(
+        filename_base: str, image5d: np.ndarray,
+        offset: Optional[Sequence[int]] = None,
+        size: Optional[Sequence[int]] = None,
+        channels: Optional[Sequence[int]] = None, verify: bool = False,
+        save_dfs: bool = True, full_roi: bool = False, coloc: bool = False
+) -> Tuple[Tuple[int, int, int], str, "detector.Blobs"]:
     """Detect blobs by block processing of a large image.
     
     All channels are processed in the same blocks.
@@ -332,19 +336,20 @@ def detect_blobs_blocks(filename_base, image5d, offset, size, channels,
         image5d: Large image to process as a Numpy array of t,z,y,x,[c]
         offset: Sub-image offset given as coordinates in z,y,x.
         size: Sub-image shape given in z,y,x.
-        channels (Sequence[int]): Sequence of channels, where None detects
+        channels: Sequence of channels, where None detects
             in all channels.
         verify: True to verify detections against truth database; defaults 
             to False.
         save_dfs: True to save data frames to file; defaults to True.
-        full_roi (bool): True to treat ``image5d`` as the full ROI; defaults
-            to False.
-        coloc (bool): True to perform blob co-localizations; defaults to False.
+        full_roi: True to treat ``image5d`` as the full ROI; defaults to False.
+        coloc: True to perform blob co-localizations; defaults to False.
     
     Returns:
-        tuple[int, int, int], str, :class:`magmap.cv.detector.Blobs`:
-        Accuracy metrics from :class:`magmap.cv.detector.verify_rois`,
-        feedback message from this same function, and detected blobs.
+        Tuple of:
+        - ``stats_detection``: accuracy metrics from
+          :class:`magmap.cv.detector.verify_rois`
+        - ``fdbk``: feedback message from this same function
+        - ``blobs``: detected blobs
     
     """
     time_start = time()
@@ -489,14 +494,17 @@ def detect_blobs_blocks(filename_base, image5d, offset, size, channels,
     return stats_detection, fdbk, blobs
 
 
-def detect_blobs_stack(filename_base, subimg_offset, subimg_size, coloc=False):
+def detect_blobs_stack(
+        filename_base: str, subimg_offset: Optional[Sequence[int]] = None,
+        subimg_size: Optional[Sequence[int]] = None, coloc: bool = False
+) -> Tuple[Tuple[int, int, int], str, "detector.Blobs"]:
     """Detect blobs in a full stack, such as a whole large image.
     
     Process channels in separate sets of blocks if their profiles specify
     different block sizes.
     
     Args:
-        filename_base (str): 
+        filename_base (str): Base image filename, for saving files.
         subimg_offset (Sequence[int]): Sub-image offset as ``z,y,x`` to load
             from :attr:`config.image5d`; defaults to None.
         subimg_size (Sequence[int]): Sub-image size as ``z,y,x`` to load
@@ -506,23 +514,26 @@ def detect_blobs_stack(filename_base, subimg_offset, subimg_size, coloc=False):
             use the ``coloc_match`` task
             (:meth:`magmap.colocalizer.StackColocalizer.colocalize_stack`)
             instead.
-
+    
     Returns:
-        tuple[int, int, int], str, :class:`magmap.cv.detector.Blobs`:
-        Combined ccuracy metrics from :class:`magmap.cv.detector.verify_rois`,
-        feedback message from this same function, and detected blobs across
-        all channels in :attr:`magmap.settings.config.channel`.
+        Tuple of:
+        - ``stats_detection``: combined accuracy metrics from
+          :class:`magmap.cv.detector.verify_rois`
+        - ``fdbk``: feedback message from this same function
+        - ``blobs``: detected blobs across all channels in
+          :attr:`magmap.settings.config.channel`
 
     """
     channels = plot_3d.setup_channels(config.image5d, config.channel, 4)[1]
     if roi_prof.ROIProfile.is_identical_settings(
             [config.get_roi_profile(c) for c in channels],
             roi_prof.ROIProfile.BLOCK_SIZES):
-        print("Will process channels together in the same blocks")
+        _logger.info("Will process channels together in the same blocks")
+        # combine channels into nested list
         channels = [channels]
     else:
-        print("Will process channels in separate blocks defined by their "
-              "profiles")
+        _logger.info(
+            "Will process channels in separate blocks defined by profiles")
     
     cols = ("stats", "fdbk", "blobs")
     detection_out = {}
@@ -556,6 +567,15 @@ def detect_blobs_stack(filename_base, subimg_offset, subimg_size, coloc=False):
         blobs_all.colocalizations = libmag.combine_arrs(
             [b.colocalizations for b in detection_out["blobs"]
              if b.colocalizations is not None])
+
+        if channels:
+            try:
+                # classify blobs if model is set in config
+                classifier.ClassifyImage.classify_whole_image(
+                    channels=channels[0], blobs=blobs_all)
+            except FileNotFoundError as e:
+                _logger.debug(e)
+        
         blobs_all.save_archive()
         print()
         
diff --git a/magmap/gui/visualizer.py b/magmap/gui/visualizer.py
index 9b1a3f571..eca9dd3cb 100644
--- a/magmap/gui/visualizer.py
+++ b/magmap/gui/visualizer.py
@@ -60,7 +60,8 @@
 
 from magmap.atlas import ontology
 from magmap.brain_globe import bg_controller
-from magmap.cv import colocalizer, cv_nd, detector, segmenter, verifier
+from magmap.cv import classifier, colocalizer, cv_nd, detector, segmenter, \
+    verifier
 from magmap.gui import atlas_editor, atlas_threads, import_threads, \
     roi_editor, verifier_editor, vis_3d, vis_handler
 from magmap.io import cli, importer, libmag, load_env, naming, np_io, sitk_io, \
@@ -345,6 +346,8 @@ class Visualization(HasTraits):
         ["-1"],
         tooltip="All blob labels will be set to this value when running Detect",
     )
+    _segs_model_path = Str  # blobs classifier model path
+    _segs_model_btn = Button("Browse")  # button to select model file
     _segs_visible = List  # blob visibility options
     _colocalize = List  # blob co-localization options
     _blob_color_style = List  # blob coloring
@@ -733,6 +736,10 @@ class Visualization(HasTraits):
                      format_func=lambda x: x)),
             Item("_btn_verifier", show_label=False),
         ),
+        HGroup(
+            Item("_segs_model_path", label="Classifier model", style="simple"),
+            Item("_segs_model_btn", show_label=False),
+        ),
         HGroup(
             Item("_segs_visible", style="custom", show_label=False,
                  editor=CheckListEditor(
@@ -1038,6 +1045,8 @@ def __init__(self):
         # set up blobs
         self.blobs = detector.Blobs()
         self._blob_color_style = [BlobColorStyles.ATLAS_LABELS.value]
+        model = config.classifier.model
+        self._segs_model_path = model if model else ""
 
         # set up profiles selectors
         self._profiles_cats = [ProfileCats.ROI.value]
@@ -2635,11 +2644,11 @@ def get_atlas_cmap(coords):
             
             else:
                 # default to color by channel
-                chls = detector.Blobs.get_blobs_channel(
+                blob_chls = detector.Blobs.get_blobs_channel(
                     segs_all[self.segs_in_mask]).astype(np.int)
                 cmap = colormaps.discrete_colormap(
-                    max(chls) + 1, alpha, True, config.seed)
-                self.segs_cmap = cmap[chls]
+                    max(blob_chls) + 1, alpha, True, config.seed)
+                self.segs_cmap = cmap[blob_chls]
         
         if self._DEFAULTS_2D[2] in self._check_list_2d:
             blobs = self.segments[self.segs_in_mask]
@@ -2666,6 +2675,16 @@ def get_atlas_cmap(coords):
             '''
         #detector.show_blob_surroundings(self.segments, self.roi)
         
+        if self._segs_model_path and self.blobs.blobs is not None:
+            try:
+                # classify blobs if model is set
+                classifier.classify_blobs(
+                    self._segs_model_path, config.image5d, offset[::-1],
+                    roi_size[::-1], chls, self.blobs, blobs_relative=True)
+            except ModuleNotFoundError as e:
+                _logger.exception(e)
+                self.segs_feedback += f"\n{e.msg}"
+        
         if (self.selected_viewer_tab is vis_handler.ViewerTabs.ROI_ED or
                 self.selected_viewer_tab is vis_handler.ViewerTabs.MAYAVI and
                 self.stale_viewers[vis_handler.ViewerTabs.MAYAVI]):
@@ -2721,7 +2740,17 @@ def _colocalize_blobs(self):
         if self.roi_ed:
             self.roi_ed.show_colocalized_blobs(
                 ColocalizeOptions.INTENSITY.value in self._colocalize)
-        
+    
+    @on_trait_change("_segs_model_btn")
+    def _segs_model_path_updated(self):
+        """Open a Pyface file dialog with path set to the classifer model path.
+        """
+        open_dialog = FileDialog(
+            action="open", default_path=os.path.dirname(self._filename))
+        if open_dialog.open() == OK:
+            # get user selected path
+            self._segs_model_path = open_dialog.path
+    
     @on_trait_change("_segs_visible")
     def _update_blob_visibility(self):
         """Change blob visibilty based on toggle check box."""
diff --git a/magmap/io/cli.py b/magmap/io/cli.py
index 98d96b305..8a6e45c5e 100644
--- a/magmap/io/cli.py
+++ b/magmap/io/cli.py
@@ -24,7 +24,8 @@
 """
 
 import argparse
-from enum import Enum 
+import dataclasses
+from enum import Enum
 import logging
 import os
 import sys
@@ -35,7 +36,7 @@
 
 from magmap.atlas import register, transformer
 from magmap.cloud import notify
-from magmap.cv import chunking, colocalizer, stack_detect
+from magmap.cv import chunking, classifier, colocalizer, stack_detect
 from magmap.io import df_io, export_stack, importer, libmag, naming, np_io, \
     sqlite
 from magmap.plot import colormaps, plot_2d
@@ -119,24 +120,27 @@ def args_with_dict(args: List[str]) -> List[Union[Dict[str, Any], str, int]]:
     return parsed
 
 
-def args_to_dict(args: List[str], keys_enum: Type[Enum],
-                 args_dict: Optional[Dict[Enum, Any]] = None,
-                 sep_args: str = "=", sep_vals: str = ",",
-                 default: Optional[Any] = None) -> Dict[Enum, Any]:
+def args_to_dict(
+        args: List[str], keys: Union[Type[Enum], Type["config.DataClass"]],
+        args_dict: Optional[Dict[Enum, Any]] = None, sep_args: str = "=",
+        sep_vals: str = ",", default: Optional[Any] = None
+) -> Union[Dict[Enum, Any], Type["config.DataClass"]]:
     """Parse positional and keyword-based arguments to an enum-keyed dictionary.
     
     Args:
         args: List of arguments with positional values followed by
             ``sep_args``-delimited values. Positional values will be entered
-            in the existing order of ``keys_enum`` based on member values, 
-            while keyword-based values will be entered if an enum 
-            member corresponding to the keyword exists.
-            Entries can also be ``sep_vals``-delimited to specify lists.
-        keys_enum: Enum class to use as keys for dictionary. Values are
-            assumed to range from 1 to number of members as output 
-            by the default Enum functional API.
+            in the existing order of ``keys`` based on member values,
+            while keyword-based values will be entered if a member
+            corresponding to the keyword exists. Entries can also be
+            ``sep_vals``-delimited to specify lists.
+        keys: Enum class or data class instance whose fields will be used as
+            keys for dictionary. Enum values are assumed to range from 1 to
+            number of members as output by the default Enum functional API.
         args_dict: Dictionary to be filled or updated with keys from
-            ``keys_enum``; defaults to None, which will assign an empty dict.
+            ``keys_enum``. Defaults to None, which will assign an empty dict.
+            Ignored if ``keys_class`` is a data class instance, which will
+            be updated instead.
         sep_args: Separator between arguments and values; defaults to "=".
         sep_vals: Separator within values; defaults to ",".
         default: Default value for each argument. Effectively turns off
@@ -145,48 +149,61 @@ def args_to_dict(args: List[str], keys_enum: Type[Enum],
             undergo splitting by ``sep_vals``.
     
     Returns:
-        Dictionary filled with arguments. Values that contain commas
-        will be split into comma-delimited lists. All values will be 
-        converted to ints if possible.
+        Dictionary or data class corresponding to ``keys``, filled with
+        arguments. Values that contain commas will be split into
+        comma-delimited lists. All values will be converted to ints if possible.
     
     """
-    if args_dict is None:
-        args_dict = {}
+    is_data = dataclasses.is_dataclass(keys)
+    if is_data:
+        # use fields as keys
+        keys_data = [f.name for f in dataclasses.fields(keys)]
+        nkeys = len(keys_data)
+        out = keys
+    else:
+        # use Enum members as keys
+        keys_data = None
+        nkeys = len(keys)
+        out = {} if args_dict is None else args_dict
+    
     by_position = True
-    num_enums = len(keys_enum)
     for i, arg in enumerate(args):
         arg_split = arg.split(sep_args)
         if default and len(arg_split) < 2:
             # add default value unless another value is given
             arg_split.append(default)
         len_arg_split = len(arg_split)
+        
         # assume by position until any keyword given
         by_position = by_position and len_arg_split < 2
         key = None
         vals = arg
         if by_position:
-            # positions are based on enum vals, assumed to range from 
+            # positions are based on enum vals, assumed to range from
             # 1 to num of members
             n = i + 1
-            if n > num_enums:
+            if n > nkeys:
                 _logger.warn(
-                    "No further parameters in {} to assign \"{}\" by "
-                    "position, skipping".format(keys_enum, arg))
+                    "No further parameters in '%s' to assign '%s' by "
+                    "position, skipping", keys, arg)
                 continue
-            key = keys_enum(n)
+            key = keys_data[n] if is_data else keys(n)
+        
         elif len_arg_split < 2:
-            _logger.warn("parameter {} does not contain a keyword, skipping"
-                         .format(arg))
+            _logger.warn(
+                "Parameter '%s' does not contain a keyword, skipping", arg)
+        
         else:
             # assign based on keyword if its equivalent enum exists
+            key_str = arg_split[0]
             vals = arg_split[1]
-            key_str = arg_split[0].upper()
             try:
-                key = keys_enum[key_str]
+                key = key_str.lower() if is_data else keys[key_str.upper()]
             except KeyError:
-                _logger.warn("Unable to find '{}' in {}, skipping".format(
-                    key_str, keys_enum))
+                _logger.warn(
+                    "Unable to find '%s' in %s, skipping", key_str, keys)
                 continue
+        
         if key:
             if isinstance(vals, str):
                 # split delimited strings
@@ -195,12 +212,17 @@ def args_to_dict(args: List[str], keys_enum: Type[Enum],
                     vals = vals_split
                 # cast to numeric type if possible
                 vals = libmag.get_int(vals)
-            # assign to found enum
-            args_dict[key] = vals
-    return args_dict
+            # assign to found enum to data class
+            if is_data:
+                setattr(out, key, vals)
+            else:
+                out[key] = vals
+    
+    return out
 
 
-def _get_args_dict_help(msg: str, keys: Type[Enum]) -> str:
+def _get_args_dict_help(
+        msg: str, keys: Union[Type[Enum], Type["config.DataClass"]]) -> str:
     """Get the help message for command-line arguments that are converted
     to a dictionary.
     
@@ -212,9 +234,15 @@ def _get_args_dict_help(msg: str, keys: Type[Enum]) -> str:
         Help message with available keys.
 
     """
-    return ("{} Available keys follow this order until the first "
-            "key=value pair is given: {}".format(
-                msg, libmag.enum_names_aslist(keys)))
+    if dataclasses.is_dataclass(keys):
+        # get all fields from data class
+        names = [f.name for f in dataclasses.fields(keys)]
+    else:
+        # use enum names
+        names = libmag.enum_names_aslist(keys)
+    
+    return (f"{msg} Available keys, which follow this order positionally "
+            f"until the first key=value pair is given: {names}")
 
 
 def process_cli_args():
@@ -353,6 +381,11 @@ def process_cli_args():
     parser.add_argument(
         "--set_meta", nargs="*",
         help="Set metadata values; see config.MetaKeys for settings")
+    parser.add_argument(
+        "--classifier", nargs="*",
+        help=_get_args_dict_help(
+            "Classifier values; see config.ClassifierKeys for settings.",
+            config.ClassifierData))
 
     # image and figure display arguments
     parser.add_argument(
@@ -727,6 +760,11 @@ def process_cli_args():
     # set up Matplotlib styles/themes
     plot_2d.setup_style()
     
+    if args.classifier is not None:
+        # classifier settings
+        args_to_dict(args.classifier, config.classifier)
+        print("Set classifier to {}".format(config.classifier))
+
     if args.db:
         # set main database path to user arg
         config.db_path = args.db
@@ -754,18 +792,58 @@ def process_cli_args():
             f"ignored: {args_unknown}")
 
 
+def setup_image(
+        path: str, series: Optional[int] = None,
+        proc_tasks: Optional[Dict["config.ProcessTypes", Any]] = None):
+    """Set up the main image from CLI args and process any tasks.
+    
+    Args:
+        path: Image path.
+        series: Image series, such as a tile; defaults to None.
+        proc_tasks: Dictionary of processing tasks; defaults to None.
+
+    """
+    # deconstruct user-supplied image filename
+    filename, offset, size, reg_suffixes = importer.deconstruct_img_name(path)
+    set_subimg, _ = importer.parse_deconstructed_name(
+        filename, offset, size, reg_suffixes)
+    
+    if not set_subimg:
+        # sub-image parameters set in filename takes precedence for
+        # the loaded image, but fall back to user-supplied args
+        offset = config.subimg_offsets[0] if config.subimg_offsets else None
+        size = config.subimg_sizes[0] if config.subimg_sizes else None
+    
+    if proc_tasks:
+        for proc_task, proc_val in proc_tasks.items():
+            # set up image for the given task
+            np_io.setup_images(
+                filename, series, offset, size, proc_task,
+                fallback_main_img=False)
+            process_file(
+                filename, proc_task, proc_val, series, offset, size,
+                config.roi_offsets[0] if config.roi_offsets else None,
+                config.roi_sizes[0] if config.roi_sizes else None)
+        
+    else:
+        # set up image without a task specified, eg for display
+        np_io.setup_images(filename, series, offset, size)
+    
+
 def process_proc_tasks(
         path: Optional[str] = None,
         series_list: Optional[Sequence[int]] = None
-) -> Optional[Dict[config.ProcessTypes, Any]]:
+) -> Optional[Dict["config.ProcessTypes", Any]]:
     """Apply processing tasks.
     
     Args:
         path: Base path to main image file; defaults to None, in which case
             :attr:`config.filename` will be used.
-        series_list: 
+        series_list: Sequence of images series, such as tiles; defaults
+            to None.
 
     Returns:
+        Dictionary of set processing types.
 
     """
     if path is None:
@@ -781,30 +859,8 @@ def process_proc_tasks(
     for series in series_list:
         # process files for each series, typically a tile within a
         # microscopy image set or a single whole image
-        filename, offset, size, reg_suffixes = \
-            importer.deconstruct_img_name(path)
-        set_subimg, _ = importer.parse_deconstructed_name(
-            filename, offset, size, reg_suffixes)
-        if not set_subimg:
-            # sub-image parameters set in filename takes precedence for
-            # the loaded image, but fall back to user-supplied args
-            offset = (config.subimg_offsets[0] if config.subimg_offsets
-                      else None)
-            size = (config.subimg_sizes[0] if config.subimg_sizes
-                    else None)
-        if proc_tasks:
-            for proc_task, proc_val in proc_tasks.items():
-                # set up image for the given task
-                np_io.setup_images(
-                    filename, series, offset, size, proc_task,
-                    fallback_main_img=False)
-                process_file(
-                    filename, proc_task, proc_val, series, offset, size,
-                    config.roi_offsets[0] if config.roi_offsets else None,
-                    config.roi_sizes[0] if config.roi_sizes else None)
-        else:
-            # set up image without a task specified, eg for display
-            np_io.setup_images(filename, series, offset, size)
+        setup_image(path, series, proc_tasks)
+    
     return proc_tasks
     
 
@@ -1116,7 +1172,7 @@ def process_file(
     Assumes that the image has already been set up.
     
     Args:
-        path: Path to image from which MagellanMapper-style paths will 
+        path: Path to image from which MagellanMapper-style paths will
             be generated.
         proc_type: Processing type, which should be a one of
             :class:`config.ProcessTypes`.
@@ -1131,8 +1187,9 @@ def process_file(
             ``(x, y, z)``; defaults to None.
     
     Returns:
-        Tuple of stats from processing, or None if no stats, and 
+        Tuple of stats from processing, or None if no stats, and
         text feedback from the processing, or None if no feedback.
+    
     """
     # PROCESS BY TYPE
     stats = None
@@ -1149,9 +1206,8 @@ def process_file(
         print("imported {}, will exit".format(path))
     
     elif proc_type is config.ProcessTypes.EXPORT_ROIS:
-        # export ROIs; assumes that info_proc was already loaded to 
-        # give smaller region from which smaller ROIs from the truth DB 
-        # will be extracted
+        # export ROIs; assumes that metadata was already loaded to give smaller
+        # region from which smaller ROIs from the truth DB will be extracted
         from magmap.io import export_rois
         db = config.db if config.truth_db is None else config.truth_db
         export_path = naming.make_subimage_name(
@@ -1165,7 +1221,7 @@ def process_file(
     elif proc_type is config.ProcessTypes.TRANSFORM:
         # transpose, rescale, and/or resize whole large image
         transformer.transpose_img(
-            path, series, plane=config.plane, 
+            path, series, plane=config.plane,
             rescale=config.transform[config.Transforms.RESCALE],
             target_size=config.roi_size)
         
@@ -1205,6 +1261,14 @@ def process_file(
             colocalizer.insert_matches(config.db, matches)
         else:
             print("No blobs loaded to colocalize, skipping")
+    
+    elif proc_type is config.ProcessTypes.CLASSIFY:
+        # classify blobs
+        try:
+            classifier.ClassifyImage.classify_whole_image()
+            config.blobs.save_archive()
+        except FileNotFoundError as e:
+            _logger.debug(e)
 
     elif proc_type in (config.ProcessTypes.EXPORT_PLANES,
                        config.ProcessTypes.EXPORT_PLANES_CHANNELS):
diff --git a/magmap/io/export_regions.py b/magmap/io/export_regions.py
index de395cc83..5948abf7d 100644
--- a/magmap/io/export_regions.py
+++ b/magmap/io/export_regions.py
@@ -185,7 +185,8 @@ def make_density_image(
         labels_img_sitk: Optional[sitk.Image] = None,
         channel: Optional[Sequence[int]] = None,
         matches: Dict[Tuple[int, int], "colocalizer.BlobMatch"] = None,
-        atlas_profile: Optional["atlas_prof.AtlasProfile"] = None
+        atlas_profile: Optional["atlas_prof.AtlasProfile"] = None,
+        include: Sequence[int] = None
 ) -> Tuple[np.ndarray, str]:
     """Make a density image based on associated blobs.
     
@@ -218,6 +219,8 @@ def make_density_image(
         matches: Dictionary of channel combinations to blob matches; defaults
             to None.
         atlas_profile: Atlas profile, used for scaling; defaults to None.
+        include: Sequence of blob ``confirmed`` flags to include; defaults
+            to None, in which case all flags will be included.
     
     Returns:
         Tuple of the density image as a Numpy array in the
@@ -288,14 +291,25 @@ def make_heat_map():
             labels_img_sitk.SetSpacing(labels_spacing[::-1])
     _logger.debug("Using image scaling: {}".format(scaling))
     
-    # annotate blobs based on position
     blobs_chl = blobs.blobs
+    _logger.info("Initial number of blobs: %s", len(blobs_chl))
     if channel is not None:
+        # filter blobs by channel
         _logger.info(
             "Using blobs from channel(s), combining if multiple channels: %s",
             channel)
         blobs_chl = blobs_chl[np.isin(detector.Blobs.get_blobs_channel(
             blobs_chl), channel)]
+        _logger.info("Number of remaining blobs: %s", len(blobs_chl))
+
+    if include is not None:
+        # filter blobs by confirmation flag
+        _logger.info("Using blobs with confirmed flag(s): %s", include)
+        blobs_chl = blobs_chl[np.isin(detector.Blobs.get_blob_confirmed(
+            blobs_chl), include)]
+        _logger.info("Number of remaining blobs: %s", len(blobs_chl))
+    
+    # annotate blobs based on position
     heat_map = make_heat_map()
     if is_2d:
         # convert back to 3D
@@ -386,7 +400,8 @@ def make_density_images_mp(img_paths, scale=None, shape=None, suffix=None,
         pool_results.append(pool.apply_async(
             make_density_image,
             args=(img_path, scale, shape, suffix, None, channel, matches,
-                  config.atlas_profile)))
+                  config.atlas_profile,
+                  config.classifier.include)))
     for result in pool_results:
         _, path = result.get()
         print("finished {}".format(path))
diff --git a/magmap/io/np_io.py b/magmap/io/np_io.py
index abbd7caca..086e92c2b 100644
--- a/magmap/io/np_io.py
+++ b/magmap/io/np_io.py
@@ -513,6 +513,7 @@ def add_metadata():
         # crop full image to bounds of sub-image
         config.image5d = plot_3d.prepare_subimg(
             config.image5d, offset, size)[None]
+        config.img5d.img = config.image5d
         config.image5d_is_roi = True
 
     # add any additional image5d thresholds for multichannel images, such 
diff --git a/magmap/plot/plot_3d.py b/magmap/plot/plot_3d.py
index 9ccb57c60..7ecdc4008 100644
--- a/magmap/plot/plot_3d.py
+++ b/magmap/plot/plot_3d.py
@@ -4,7 +4,7 @@
 
 Prepare volumetric image stacks for plotting.
 """
-from typing import Optional, Sequence, Tuple
+from typing import Optional, Sequence, Tuple, Union
 
 import numpy as np
 from skimage import draw
@@ -322,61 +322,69 @@ def get_isotropic_vis(settings):
     return isotropic
 
 
-def prepare_subimg(image5d, offset, size, ndim_base=5):
+def prepare_subimg(
+        image5d: np.ndarray, offset: Sequence[int], size: Sequence[int],
+        ndim_base: int = 5, border: Optional[Union[int, Sequence[int]]] = None
+) -> np.ndarray:
     """Extracts a subimage from a larger image.
     
     Args:
-        image5d (:obj:`np.ndarray`): 5D image array in the order,
-            ``t,z,y,x[,c]``, where the final dimension is optional as with
-            many one channel images.
-        offset (List[int]): Tuple of offset given as ``z,y,x`` for the region
-            of interest. Defaults to ``(0, 0, 0)``.
-        size (List[int]): Size of the region of interest as ``z,y,x``.
-        ndim_base (int): Number of dimensions on which ``image5d`` is based.
-            Typically 3 or 5, defaulting to 5 as ``t,z,y,x[,c]``.
-            If 3, the ``t`` dimension is removed.
+        image5d: 5D image array in the order, ``t, z, y, x, [c]`` unless
+            ``ndim_base`` is given.
+        offset: Tuple of offset given as ``z, y, x`` for the region of interest.
+        size: Size of the region of interest as ``z, y, x``.
+        ndim_base: Number of dimensions in ``image5d``. Defaults to 5 for
+            `t, z, y, x, [c]``; if fewer, the ``t`` dimension is removed.
+        border: Border size as either a scalar for all dimensions of values
+            for ``z, y, x``. The resulting subimage will be expanded to
+            include this border outside ``offset``, with the same extra size
+            on opposite sides along each axis.
     
     Returns:
-        :obj:`np.ndarray`: The sub-imge without separate time dimension as
-        a 3D (or 4-D array if channel dimension exists) array.
+        The sub-imge without separate time dimension as a 3D (or 4-D array
+        if channel dimension exists) array.
     
     """
-    cube_slices = [slice(o, o + s) for o, s in zip(offset, size)]
-    libmag.printv("preparing sub-image at offset: {}, size: {}, slices: {}"
-                  .format(offset, size, cube_slices))
+    if border is not None:
+        # add border outside offset with same extra space on opposite sides
+        offset = np.subtract(offset, border)
+        size = np.add(size, np.multiply(border, 2))
+
+    # ROI with corner at offset and shape given by size
+    roi_slices = [slice(o, o + s) for o, s in zip(offset, size)]
     
-    # cube with corner at offset and shape given by size
     img = image5d
     if ndim_base >= 5:
         # remove time axis
         img = image5d[0]
-    return img[cube_slices[0], cube_slices[1], cube_slices[2]]
+    return img[roi_slices[0], roi_slices[1], roi_slices[2]]
 
 
-def prepare_roi(image5d, roi_offset, roi_size, ndim_base=5):
+def prepare_roi(
+        image5d: np.ndarray, roi_offset: Sequence[int],
+        roi_size: Sequence[int], *args, **kwargs) -> np.ndarray:
     """Extracts a region of interest (ROI).
 
     Calls :meth:`prepare_subimage` but expects size and offset variables to
     be in x,y,z order following this software's legacy convention.
 
     Args:
-        image5d (:obj:`np.ndarray`): 5D image array in the order,
+        image5d: 5D image array in the order,
             ``t,z,y,x[,c]``, where the final dimension is optional as with
             many one channel images.
-        roi_offset (List[int]): Tuple of offset given as ``x,y,z`` for the region
+        roi_offset: Tuple of offset given as ``x,y,z`` for the region
             of interest. Defaults to ``(0, 0, 0)``.
-        roi_size (List[int]): Size of the region of interest as ``x,y,z``.
-        ndim_base (int): Number of dimensions on which ``image5d`` is based.
-            Typically 3 or 5, defaulting to 5 as ``t,z,y,x[,c]``.
-            If 3, the ``t`` dimension is removed.
+        roi_size: Size of the region of interest as ``x,y,z``.
+        args: Positional arguments to :meth:`prepare_subroi`.
+        kwargs: Named arguments to :meth:`prepare_subroi`.
 
     Returns:
-        :obj:`np.ndarray`: The region of interest without separate time
+        The region of interest without separate time
         dimension as a 3D (or 4-D array if channel dimension exists) array.
     
     """
-    libmag.printv("preparing ROI at x,y,z:")
-    return prepare_subimg(image5d, roi_offset[::-1], roi_size[::-1], ndim_base)
+    return prepare_subimg(
+        image5d, roi_offset[::-1], roi_size[::-1], *args, **kwargs)
 
 
 def roi_center_to_offset(
diff --git a/magmap/settings/config.py b/magmap/settings/config.py
index 551742dad..11b154ab6 100644
--- a/magmap/settings/config.py
+++ b/magmap/settings/config.py
@@ -9,9 +9,13 @@
 
 """
 
+import dataclasses
 from enum import Enum, auto
 import pathlib
-from typing import Any, Dict, Optional, Sequence, TYPE_CHECKING, Union
+from typing import Any, Dict, Optional, Sequence, TypeVar, TYPE_CHECKING, Union
+
+# imports from typing directly if available
+from typing_extensions import Protocol
 
 try:
     from appdirs import AppDirs
@@ -42,6 +46,15 @@
 DNS_REVERSE = f"io.github.sanderslab.{APP_NAME}"
 
 
+class DataClassProtocol(Protocol):
+    """Typing protocol for data classes."""
+    __dataclass_fields__: Dict[str, Any]
+
+
+#: Type hint for data classes.
+DataClass = TypeVar("DataClass", bound=DataClassProtocol)
+
+
 class DocsURLs(Enum):
     """URLs to online documentation."""
     #: Docs base URL.
@@ -257,6 +270,8 @@ class ProcessTypes(Enum):
     DETECT_COLOC = auto()
     #: Detect blobs along with match-based colocalization.
     COLOC_MATCH = auto()
+    #: Classify blobs.
+    CLASSIFY = auto()
     #: Load previously processed images and blobs. DEPRECATED: use ``--load``
     #: CLI parameter instead.
     LOAD = auto()
@@ -382,6 +397,19 @@ class Transforms(Enum):
 transform = dict.fromkeys(Transforms, None)
 
 
+@dataclasses.dataclass
+class ClassifierData:
+    """Classifier data class for CLI arguments."""
+    #: Path to pre-trained classifier model.
+    model: Optional[str] = None
+    #: Classification flags (confirmed column) to include.
+    include: Optional[Union[int, Sequence[int]]] = None
+
+
+#: Classifier settings.
+classifier: ClassifierData = ClassifierData()
+
+
 # extensions for saving figures.
 
 #: tuple[str, ...]: Extension for 3D renderings.
diff --git a/setup.py b/setup.py
index fa0387359..b42745ce5 100644
--- a/setup.py
+++ b/setup.py
@@ -31,6 +31,9 @@
 # optional dependencies for Jupyter notebooks
 _EXTRAS_JUPYTER = ["jupyterlab", "bash_kernel"]
 
+# optional dependencies for classification
+_EXTRAS_CLASSIFER = ["tensorflow"]
+
 # installation configuration
 config = {
     "name": "magellanmapper",
@@ -71,6 +74,7 @@
         "dataclasses ; python_version < '3.7'",
         # BrainGlobe dependencies for access to cloud-hosted atlases
         "bg-atlasapi @ https://github.com/brainglobe/bg-atlasapi/archive/refs/heads/master.zip",
+        "typing_extensions",
     ],
     "extras_require": {
         "import": _EXTRAS_IMPORT,
@@ -78,6 +82,7 @@
         "pandas_plus": _EXTRAS_PANDAS,
         "docs": _EXTRAS_DOCS,
         "jupyter": _EXTRAS_JUPYTER,
+        "classifier": _EXTRAS_CLASSIFER,
         
         # dependencies for most common tasks
         "most": [


Test patch:
diff --git a/magmap/tests/test_classifier.py b/magmap/tests/test_classifier.py
new file mode 100644
index 000000000..19f03d1b2
--- /dev/null
+++ b/magmap/tests/test_classifier.py
@@ -0,0 +1,48 @@
+# MagellanMapper unit testing for classification
+"""Unit testing for the MagellanMapper classifier module."""
+
+import unittest
+
+import numpy as np
+
+from magmap.cv import classifier
+from magmap.tests import test_detector
+
+
+class TestClassifier(unittest.TestCase):
+    
+    def test_setup_classification_roi(self):
+        
+        # set up blobs
+        blobs = test_detector.TestDetector.make_random_blobs()
+        
+        # set up image
+        image5d = np.zeros((1, 10, 15, 20, 2))
+        
+        # test ROI with complete borders
+        roi, blobs_mask, blobs_shift = classifier.setup_classification_roi(
+            image5d, (2, 3, 4), (5, 6, 7), blobs, 2)
+        np.testing.assert_array_equal(roi.shape, (5, 8, 9, 2))
+        np.testing.assert_array_equal(blobs_shift, (0, 1, 1))
+        
+        # test ROI exceeding the image's far edge
+        roi, blobs_mask, blobs_shift = classifier.setup_classification_roi(
+            image5d, (2, 3, 4), (5, 14, 7), blobs, 2)
+        np.testing.assert_array_equal(roi.shape, (5, 13, 9, 2))
+        np.testing.assert_array_equal(blobs_shift, (0, 1, 1))
+        
+        # test ROI aligned with the image's near edge
+        roi, blobs_mask, blobs_shift = classifier.setup_classification_roi(
+            image5d, (0, 0, 0), (5, 6, 7), blobs, 2)
+        np.testing.assert_array_equal(roi.shape, (5, 7, 8, 2))
+        np.testing.assert_array_equal(blobs_shift, (0, 0, 0))
+        
+        # test ROI whose border partially exceeds image's near edge
+        roi, blobs_mask, blobs_shift = classifier.setup_classification_roi(
+            image5d, (0, 1, 1), (5, 6, 7), blobs, 4)
+        np.testing.assert_array_equal(roi.shape, (5, 9, 10, 2))
+        np.testing.assert_array_equal(blobs_shift, (0, 1, 1))
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/magmap/tests/test_detector.py b/magmap/tests/test_detector.py
index b63e38262..b92f827e6 100644
--- a/magmap/tests/test_detector.py
+++ b/magmap/tests/test_detector.py
@@ -10,14 +10,25 @@
 
 class TestDetector(unittest.TestCase):
     
-    def test_blob_accessors(self):
-        # set up random blobs
+    @staticmethod
+    def make_random_blobs() -> "detector.Blobs":
+        """Make blobs with random coordinates.
+        
+        Returns:
+            Blobs instance.
+
+        """
         blobs = np.random.rand(20).reshape((5, 4))
-        blobs[:, 0:3] = np.multiply(blobs[:, 0:3], 100).astype(int)
+        blobs[:, :3] = np.multiply(blobs[:, :3], 100).astype(int)
         blobs[:, 3] = blobs[:, 3] * 10
         bl = detector.Blobs(blobs)
         bl.blobs = bl.format_blobs(bl.blobs)
         print(f"Blobs:\n{bl.blobs}")
+        return bl
+    
+    def test_blob_accessors(self):
+        # set up random blobs
+        bl = self.make_random_blobs()
 
         # test blob confirmation flag
         np.testing.assert_array_equal(