Repository: sanderslab/magellanmapper (language: python)

Issue:
Support ITK-Elastix as an alternative to SimpleITK with Elastix
We have used the SimpleITK with Elastix (formerly SimpleElastix) library for image registration (see #42, #379). Although this library has worked well, it is currently not available as a precompiled wheel with Elastix on PyPi. ITK-Elastix is a package available on PyPi that provides this access to Elastix, which would remove the requirement for the custom SimpleITK build that we currently provide. 

Since both SimpleITK and ITK-Elastix use ITK, they have similar APIs though with important caveats. Our goal will be to extend support for image registration using ITK-Elastix while maintaining support for SimpleITK and eventually making SimpleITK optional for MagellanMapper. Once we no longer require custom builds of SimpleITK, we can begin to make MagellanMapper available on PyPi as wel.

Goals of PRs in this series:
- [x] Support ITK-Elastix as an alternative to SimpleITK with Elastix (#495)
- [x] Maintain support for SimpleITK with Elastix
- [x] Make SimpleITK optional (#497)
- [x] Make ITK optional (#497)
- [x] Support groupwise registration in ITK-Elastix (see #438; #500)


Discussion / hints:

Fix patch:
diff --git a/docs/install.md b/docs/install.md
index b450a82c..7ec3b8d2 100644
--- a/docs/install.md
+++ b/docs/install.md
@@ -194,14 +194,16 @@ These package versions are used for automated testing (continuous integration).
 |-----|-----|-----|
 | `most` | Import and GUI tools | Has `import`, `gui` |
 | `all` | All groups plus `seaborn`, `scikit-learn` | Has all below |
-| `import` | Imports proprietary image formats | |
-| `gui` | Main graphical interface | |
-| `pandas_plus` | Exports styled and Excel formats | |
+| `3D` | 3D rendering | |
 | `aws` | Tools for accessing AWS | |
+| `classifer` | Tensorflow | |
 | `docs` | Tools for building docs | |
+| `gui` | Main graphical interface | |
+| `import` | Imports proprietary image formats | |
+| `itk` | ITK-Elastix | |
 | `jupyter` | Running Notebooks | |
-| `classifer` | Tensorflow | |
-| `3D` | 3D rendering | |
+| `pandas_plus` | Exports styled and Excel formats | |
+| `simpleitk` | Custom SimpleITK with Elastix | |
 
 ### Optional Dependency Build and Runtime Requirements
 
diff --git a/docs/release/release_v1.6.md b/docs/release/release_v1.6.md
index ad67f60c..55fcd6d0 100644
--- a/docs/release/release_v1.6.md
+++ b/docs/release/release_v1.6.md
@@ -81,6 +81,10 @@
 - Landmark distance measurements save the raw distances and no longer require spacing (#147)
 - Masks with angle planes can be constructed (#252)
 - `register.RegImgs` is a data class to track registered images (#335)
+- Supports image I/O and registration through ITK
+  - `ITK-Elastix` support has been added as an alternative to `SimpleITK` (#495, #500)
+  - Both libraries are now optional rather than required dependencies (#497)
+  - Better support for image direction metadata (#497)
 - Fixed changes to large label IDs during registration (#303)
 
 #### Cell detection
@@ -200,7 +204,8 @@
 - The `jupyter` install group installs packages for running the Jupyter sample commands notebook in a Bash kernel (#122)
 - Missing dependencies are starting to use more consistent error messages and instructions (#226)
 - Python 3.6 and 3.7 have been deprecated for removal in MM v1.7 and have separate pinned dependencies (`envs/requirements_py3<n>`) (#232, #379)
-- The SimpleElastix custom binaries have been switched to SimpleITK with Elastix (#379)
+- The SimpleElastix custom binaries are now built on SimpleITK with Elastix and is no longer a required dependency (#379, #501)
+- ITK-Elastix is now supported for image registration (#501)
 - Custom dependency binaries are now built for Python 3.8-3.11 (#379)
 - Fixed error on deprecated NumPy data type aliases (#364)
 - Fixed `qt4 backend` error by installing PyQt <= 5.15.7 (#431)
diff --git a/magmap/atlas/atlas_refiner.py b/magmap/atlas/atlas_refiner.py
index 77337759..cd7cff8b 100644
--- a/magmap/atlas/atlas_refiner.py
+++ b/magmap/atlas/atlas_refiner.py
@@ -1,5 +1,5 @@
 # Atlas refinement
-# Author: David Young, 2019
+# Author: David Young, 2019, 2023
 """Refine atlases in 3D.
 """
 from collections import OrderedDict
@@ -7,9 +7,16 @@
 import os
 from typing import Dict, List, Optional, Sequence, Tuple, Union, Any
 
-import SimpleITK as sitk
+try:
+    import itk
+except ImportError:
+    itk = None
 import numpy as np
 import pandas as pd
+try:
+    import SimpleITK as sitk
+except ImportError:
+    sitk = None
 from skimage import filters
 from skimage import measure
 from skimage import morphology
@@ -17,7 +24,7 @@
 
 from magmap.atlas import labels_meta
 from magmap.cv import chunking, cv_nd, segmenter
-from magmap.io import df_io, export_stack, importer, libmag, np_io, sitk_io
+from magmap.io import df_io, export_stack, importer, libmag, sitk_io
 from magmap.plot import colormaps, plot_3d, plot_support
 from magmap.settings import atlas_prof, config, profiles
 
@@ -279,7 +286,7 @@ def _curate_labels(img, img_ref, mirror=None, edge=None, expand=None,
     rotation = rotate["rotation"] if rotate else None
     
     # cast to int that takes the full range of the labels image
-    img_np = sitk.GetArrayFromImage(img)
+    img_np = sitk_io.convert_img(img)
     label_ids_orig = np.unique(img_np)
     try:
         signed = True if mirror_mult == -1 else None
@@ -292,7 +299,7 @@ def _curate_labels(img, img_ref, mirror=None, edge=None, expand=None,
         # fallback to large signed data type
         print(e)
         img_np.astype(np.int32)
-    img_ref_np = sitk.GetArrayFromImage(img_ref)
+    img_ref_np = sitk_io.convert_img(img_ref)
     tot_planes = len(img_np)
     
     # lateral edges of atlas labels (ie low z in sagittal orientation) are 
@@ -686,8 +693,8 @@ def _smoothing(img_np, img_np_orig, filter_size, spacing=None):
     
     print("\nMeasuring foreground overlap of labels after smoothing:")
     measure_overlap_labels(
-        make_labels_fg(sitk.GetImageFromArray(img_np)), 
-        make_labels_fg(sitk.GetImageFromArray(img_np_orig)))
+        make_labels_fg(sitk_io.convert_img(img_np)), 
+        make_labels_fg(sitk_io.convert_img(img_np_orig)))
     
     return filter_size, df_metrics, df_raw
 
@@ -908,8 +915,8 @@ def smooth_labels(
     # show DSC for labels
     print("\nMeasuring overlap of labels:")
     measure_overlap_labels(
-        sitk.GetImageFromArray(labels_img_np_orig), 
-        sitk.GetImageFromArray(labels_img_np))
+        sitk_io.convert_img(labels_img_np_orig), 
+        sitk_io.convert_img(labels_img_np))
     
     # weighted pixel ratio metric of volume change
     weighted_size_ratio = 0
@@ -1081,14 +1088,14 @@ def aggr_smoothing_metrics(df_pxs):
 
 
 def transpose_img(
-        img_sitk: sitk.Image, plane: Optional[str] = None,
+        img_sitk: "sitk.Image", plane: Optional[str] = None,
         rotate: Optional[int] = None,
         rotate_deg: Optional[Sequence[Dict[str, Any]]] = None,
         target_size: Optional[Union[float, Sequence[int]]] = None,
         target_size_res: Optional[Union[float, Sequence[int]]] = None,
         flip: Optional[Union[int, Sequence[int]]] = None,
         order: Optional[int] = None
-) -> sitk.Image:
+) -> "sitk.Image":
     """Transpose an image to a different plane or rotation.
     
     Supports the `ROTATE` and `FLIP` settings in
@@ -1154,9 +1161,9 @@ def transpose_img(
         return img_sitk
     
     # transform as a np array
-    img = sitk.GetArrayFromImage(img_sitk)
-    spacing = img_sitk.GetSpacing()[::-1]
-    origin = img_sitk.GetOrigin()[::-1]
+    img = sitk_io.convert_img(img_sitk)
+    spacing = list(img_sitk.GetSpacing())[::-1]
+    origin = list(img_sitk.GetOrigin())[::-1]
     transposed = img
     is_2d = transposed.ndim == 2
     
@@ -1228,30 +1235,37 @@ def transpose_img(
         spacing = np.multiply(spacing, target_size_res)
     
     # convert back to sitk image
-    transposed = sitk.GetImageFromArray(transposed)
+    transposed = sitk_io.convert_img(transposed)
     transposed.SetSpacing(spacing[::-1])
     transposed.SetOrigin(origin[::-1])
     return transposed
 
 
-def match_atlas_labels(img_atlas, img_labels, metrics=None):
-    """Apply register profile settings to labels and match atlas image 
+def match_atlas_labels(
+        img_atlas: Union["sitk.Image", "itk.Image"],
+        img_labels: Union["sitk.Image", "itk.Image"],
+        metrics: Optional[
+            Dict["config.AtlasMetrics", Union[float, Sequence[float]]]] = None
+) -> Tuple[Union["sitk.Image", "itk.Image"],
+           Union["sitk.Image", "itk.Image"],
+           Optional[pd.DataFrame], Optional[pd.DataFrame]]:
+    """Apply register profile settings to labels and match atlas image
     accordingly.
     
     Args:
-        img_atlas (:obj:`sitk.Image`): Reference image, such as histology.
-        img_labels (:obj:`sitk.Image`): Labels image.
-        metrics (:obj:`dict`): Dictionary to store metrics; defaults to 
-            None, in which case metrics will not be measured.
+        img_atlas: Reference image, such as histology.
+        img_labels: Labels image.
+        metrics: Dictionary of metrics. If None, metrics will not be measured.
     
     Returns:
-        Tuple: ``img_atlas``, the updated atlas; ``img_labels``, the 
-        updated labels; ``img_borders``, a new (:obj:`sitk.Image`) of the 
-        same shape as the prior images except an extra channels dimension 
-        as given by :func:``_curate_labels``; ``df_sm``, a 
-        data frame of smoothing stats, or None if smoothing was not performed; 
-        and ``df_sm_raw``, a data frame of raw smoothing stats, or 
-        None if smoothing was not performed.
+        Tuple of:
+        - ``img_atlas``: the updated atlas
+        - ``img_labels``: the updated labels
+        - ``df_sm``: a data frame of smoothing stats, or None if smoothing was
+          not performed
+        - ``df_sm_raw``: a data frame of raw smoothing stats, or None if
+          smoothing was not performed
+    
     """
     pre_plane = config.atlas_profile["pre_plane"]
     mirror = config.atlas_profile["labels_mirror"]
@@ -1267,11 +1281,11 @@ def match_atlas_labels(img_atlas, img_labels, metrics=None):
     far_hem_neg = config.atlas_profile["make_far_hem_neg"]
     
     if pre_plane:
-        # images in the correct desired orientation may need to be 
-        # transposed prior to label curation since mirroring assumes 
+        # images in the correct desired orientation may need to be
+        # transposed prior to label curation since mirroring assumes
         # a sagittal orientation
-        img_atlas_np = sitk.GetArrayFromImage(img_atlas)
-        img_labels_np = sitk.GetArrayFromImage(img_labels)
+        img_atlas_np = sitk_io.convert_img(img_atlas)
+        img_labels_np = sitk_io.convert_img(img_labels)
         arrs_3d, _ = plot_support.transpose_images(
             pre_plane, [img_atlas_np, img_labels_np])
         img_atlas_np = arrs_3d[0]
@@ -1286,12 +1300,12 @@ def match_atlas_labels(img_atlas, img_labels, metrics=None):
         # include any lateral extension and mirroring if ACTIVE flag is on
         img_labels_np, extis, df_sm, df_sm_raw = (
             _curate_labels(
-                img_labels, img_atlas, mirror, edge, expand, rotate, smooth, 
+                img_labels, img_atlas, mirror, edge, expand, rotate, smooth,
                 affine))
     else:
         # turn off lateral extension and/or mirroring
         img_labels_np, _, df_sm, df_sm_raw = _curate_labels(
-            img_labels, img_atlas, mirror if is_mirror else None, 
+            img_labels, img_atlas, mirror if is_mirror else None,
             edge if is_edge else None, expand, rotate, smooth,
             affine)
         if metrics or crop and (
@@ -1302,13 +1316,13 @@ def match_atlas_labels(img_atlas, img_labels, metrics=None):
                   "for measurements and any cropping:")
             resize = is_mirror and mirror["start"] is not None
             lbls_np_mir, extis, _, _ = _curate_labels(
-                img_labels, img_atlas, mirror, edge, expand, rotate, None, 
+                img_labels, img_atlas, mirror, edge, expand, rotate, None,
                 affine, resize)
             mask_lbls = lbls_np_mir != 0
             print()
     
     # adjust atlas with same settings
-    img_atlas_np = sitk.GetArrayFromImage(img_atlas)
+    img_atlas_np = sitk_io.convert_img(img_atlas)
     if rotation:
         for rot in rotation:
             img_atlas_np = cv_nd.rotate_nd(
@@ -1368,11 +1382,11 @@ def match_atlas_labels(img_atlas, img_labels, metrics=None):
     
     if metrics is not None:
         
-        # meas DSC of labeled hemisphere, using the sagittal midline 
+        # meas DSC of labeled hemisphere, using the sagittal midline
         # to define the hemispheric boundaries
         dsc = measure_overlap_combined_labels(
-            sitk.GetImageFromArray(img_atlas_np[:extis[1]]), 
-            sitk.GetImageFromArray(img_labels_np[:extis[1]]), 
+            sitk_io.convert_img(img_atlas_np[:extis[1]]),
+            sitk_io.convert_img(img_labels_np[:extis[1]]),
             config.atlas_profile["overlap_meas_add_lbls"])
         metrics[config.AtlasMetrics.DSC_ATLAS_LABELS_HEM] = [dsc]
         
@@ -1384,10 +1398,10 @@ def match_atlas_labels(img_atlas, img_labels, metrics=None):
             img_labels_np[:extis[0]] != 0, thresh_atlas[:extis[0]])
         # simply treat rest of hem as labeled to focus on unlabeled lat portion
         metrics[config.AtlasMetrics.LAT_UNLBL_VOL] = 1 - (
-            (np.sum(lbl_edge) + np.sum(thresh_atlas[extis[0]:extis[1]])) 
+            (np.sum(lbl_edge) + np.sum(thresh_atlas[extis[0]:extis[1]]))
             / np.sum(thresh_atlas[:extis[1]]))
 
-        # meas frac of planes that are at least partially labeled, using 
+        # meas frac of planes that are at least partially labeled, using
         # mask labels since from atlas portion that should be labeled
         frac = 0   # mask labels fully covered, so fully labeled by this def
         if mask_lbls is not None:
@@ -1440,20 +1454,24 @@ def match_atlas_labels(img_atlas, img_labels, metrics=None):
     return img_atlas, img_labels, df_sm, df_sm_raw
 
 
-def import_atlas(atlas_dir, show=True, prefix=None):
-    """Import atlas from the given directory, processing it according 
-    to the register settings specified at :attr:``config.register_settings``.
+def import_atlas(
+        atlas_dir: str, show: bool = True, prefix: Optional[str] = None):
+    """Import atlas from the given directory.
+     
+    Process the atlas according to the register settings specified at
+    :attr:``config.register_settings``.
     
-    The imported atlas will be saved to a directory of the same path as 
-    ``atlas_dir`` except with ``_import`` appended to the end. DSC 
+    The imported atlas will be saved to a directory of the same path as
+    ``atlas_dir`` except with ``_import`` appended to the end. DSC
     will be calculated and saved as a CSV file in this directory as well.
     
     Args:
-        atlas_dir (str): Path to atlas directory.
-        show (bool): True to show the imported atlas.
-        prefix (str): Output path; defaults to None to ignore. If ends with
+        atlas_dir: Path to atlas directory.
+        show: True to show the imported atlas. Only supported for SimpleITK.
+        prefix: Output path; defaults to None to ignore. If ends with
             a file separator,``atlas_dir`` will still be used for the output
             filename; otherwise, the basename will be used for this filename.
+    
     """
     # load atlas and corresponding labels
     img_atlas, path_atlas = sitk_io.read_sitk(
@@ -1510,17 +1528,17 @@ def import_atlas(atlas_dir, show=True, prefix=None):
     if truncate:
         # truncate labels
         img_labels_np = truncate_labels(
-            sitk.GetArrayFromImage(img_labels), *truncate)
+            sitk_io.convert_img(img_labels), *truncate)
         img_labels = sitk_io.replace_sitk_with_numpy(img_labels, img_labels_np)
     
     # show labels
     print("labels output data type:", img_labels.GetPixelIDTypeAsString())
-    img_labels_np = sitk.GetArrayFromImage(img_labels)
+    img_labels_np = sitk_io.convert_img(img_labels)
     label_ids = np.unique(img_labels_np)
     print("number of labels: {}".format(label_ids.size))
     print(label_ids)
     
-    # write images with atlas saved as MagellanMapper/Numpy format to 
+    # write images with atlas saved as MagellanMapper/Numpy format to
     # allow opening as an image within MagellanMapper alongside the labels image
     imgs_write = {
         config.RegNames.IMG_ATLAS.value: img_atlas, 
@@ -1529,7 +1547,7 @@ def import_atlas(atlas_dir, show=True, prefix=None):
         imgs_write, name_prefix, copy_to_suffix=True, 
         ext=os.path.splitext(path_atlas)[1])
     config.resolutions = [img_atlas.GetSpacing()[::-1]]
-    img_ref_np = sitk.GetArrayFromImage(img_atlas)
+    img_ref_np = sitk_io.convert_img(img_atlas)
     img_ref_np = img_ref_np[None]
     importer.save_np_image(img_ref_np, name_prefix, 0)
 
@@ -1560,7 +1578,7 @@ def import_atlas(atlas_dir, show=True, prefix=None):
     lbls_meta = labels_meta.LabelsMeta(name_prefix)
     lbls_meta.path_ref = config.atlas_labels[config.AtlasLabels.PATH_REF]
     lbls_meta.region_ids_orig = np.unique(
-        sitk.GetArrayFromImage(img_labels)).tolist()
+        sitk_io.convert_img(img_labels)).tolist()
     lbls_meta.save()
     # create relative symlink to only PATH_LABELS_META in the atlas dir for
     # registering by giving the dir alone
@@ -1568,14 +1586,15 @@ def import_atlas(atlas_dir, show=True, prefix=None):
         os.path.relpath(lbls_meta.save_path, target_dir),
         os.path.join(target_dir, labels_meta.LabelsMeta.PATH_LABELS_META))
     
-    if show:
+    if show and sitk:
+        # show the atlas using the default viewer in SimpleITK
         sitk.Show(img_atlas)
         sitk.Show(img_labels)
 
 
 def measure_atlas_refinement(
-        metrics: Dict[Enum, List], img_atlas: sitk.Image,
-        img_labels: sitk.Image, atlas_profile: atlas_prof.AtlasProfile,
+        metrics: Dict[Enum, List], img_atlas: "sitk.Image",
+        img_labels: "sitk.Image", atlas_profile: atlas_prof.AtlasProfile,
         path: str = None) -> pd.DataFrame:
     """
     
@@ -1600,9 +1619,9 @@ def measure_atlas_refinement(
     metrics[config.AtlasMetrics.VOL_ATLAS] = [np.sum(atlas_mask)]
     metrics[config.AtlasMetrics.VOL_LABELS] = [np.sum(labels_mask)]
     
-    # compactness of whole atlas (non-label) image; use lower threshold for 
+    # compactness of whole atlas (non-label) image; use lower threshold for
     # compactness measurement to minimize noisy surface artifacts
-    img_atlas_np = sitk.GetArrayFromImage(img_atlas)
+    img_atlas_np = sitk_io.convert_img(img_atlas)
     thresh = atlas_profile["atlas_threshold_all"]
     thresh_atlas = img_atlas_np > thresh
     compactness, _, _ = cv_nd.compactness_3d(
@@ -1611,21 +1630,26 @@ def measure_atlas_refinement(
 
     _logger.info("\nWhole atlas stats:")
     df = df_io.dict_to_data_frame(metrics, path, show=" ")
-    return df  
+    return df
 
 
 def measure_overlap(
-        img1: sitk.Image, img2: sitk.Image, thresh_img1: Optional[float] = None,
+        img1: Union["sitk.Image", "itk.Image", np.ndarray],
+        img2: Union["sitk.Image", "itk.Image", np.ndarray],
+        thresh_img1: Optional[float] = None,
         thresh_img2: Optional[float] = None,
         add_to_img1_mask: Optional[np.ndarray] = None,
         return_masks: bool = False
-) -> Union[float, Tuple[float, sitk.Image, np.ndarray]]:
-    """Measure the Dice Similarity Coefficient (DSC) between two foreground 
-    of two images.
+) -> Union[float, Tuple[float, Union[
+        "sitk.Image", "itk.Image", np.ndarray], np.ndarray]]:
+    """Measure the Dice Similarity Coefficient between two intensity images.
+     
+    Calculdate the DSC between the foreground of two images.
     
     Args:
-        img1: First image.
-        img2: Image to compare.
+        img1: Intensity image. The image type (:class:`itk.Image` or
+            :class:`sitk.Image`) will determine which library is used.
+        img2: Intensity mage to compare.
         thresh_img1: Threshold to determine the foreground of ``img1``;
             defaults to None to determine by a mean threshold.
         thresh_img2: Threshold to determine the foreground of ``img2``;
@@ -1641,12 +1665,16 @@ def measure_overlap(
         The DSC of the foreground of the two given images, or ``np.nan`` if
         the DSC cannot be computed. If`return_masks`` is True, also returns
         the masked images as NumPy arrays.
+    
     """
-    # upper threshold does not seem be set with max despite docs for 
+    # determine library based on first image
+    is_sitk = sitk and isinstance(img1, sitk.Image)
+    
+    # upper threshold does not seem be set with max despite docs for
     # sitk.BinaryThreshold, so need to set with max explicitly
-    img1_np = sitk.GetArrayFromImage(img1)
+    img1_np = sitk_io.convert_img(img1)
     thresh_img1_up = float(np.amax(img1_np))
-    img2_np = sitk.GetArrayFromImage(img2)
+    img2_np = sitk_io.convert_img(img2)
     thresh_img2_up = float(np.amax(img2_np))
     
     # use threshold mean if lower thresholds not given
@@ -1659,101 +1687,141 @@ def measure_overlap(
         thresh_img1, thresh_img2)
     
     # similar to simple binary thresholding via Numpy
-    binary_img1 = sitk.BinaryThreshold(
-        img1, thresh_img1, thresh_img1_up)
-    binary_img1_np = sitk.GetArrayFromImage(binary_img1)
+    if is_sitk:
+        binary_img1 = sitk.BinaryThreshold(
+            img1, thresh_img1, thresh_img1_up)
+    else:
+        # get ITK docs from:
+        # itk.itkBinaryThresholdImageFilterPython.binary_threshold_image_filter
+        binary_img1 = itk.binary_threshold_image_filter(
+            img1, lower_threshold=int(thresh_img1),
+            upper_threshold=int(thresh_img1_up))
+    
+    binary_img1_np = sitk_io.convert_img(binary_img1)
     if add_to_img1_mask is not None:
         # add mask to foreground of img1
         binary_img1_np[add_to_img1_mask] = True
         binary_img1 = sitk_io.replace_sitk_with_numpy(
             binary_img1, binary_img1_np)
-    binary_img2 = sitk.BinaryThreshold(
-        img2, thresh_img2, thresh_img2_up)
-
+    
+    # threshold 2nd image
+    if is_sitk:
+        binary_img2 = sitk.BinaryThreshold(
+            img2, thresh_img2, thresh_img2_up)
+    else:
+        binary_img2 = itk.binary_threshold_image_filter(
+            img2, lower_threshold=int(thresh_img2),
+            upper_threshold=int(thresh_img2_up))
+    
     # match world info in case of slight rounding that might prevent the
     # filter from executing
     sitk_io.match_world_info(binary_img1, binary_img2)
-    overlap_filter = sitk.LabelOverlapMeasuresImageFilter()
-    total_dsc = np.nan
-    try:
-        overlap_filter.Execute(binary_img1, binary_img2)
-        total_dsc = overlap_filter.GetDiceCoefficient()
-        _logger.info("Foreground DSC: %s\n", total_dsc)
-    except RuntimeError as e:
-        _logger.warn(e)
+    
+    # find overlap between thresholded images
+    total_dsc = measure_overlap_labels(binary_img1, binary_img2)
+    
     if return_masks:
-        return total_dsc, binary_img1_np, sitk.GetArrayFromImage(binary_img2)
+        return total_dsc, binary_img1_np, sitk_io.convert_img(binary_img2)
     return total_dsc
 
 
-def measure_overlap_labels(labels_img1, labels_img2):
-    """Measure the mean Dice Similarity Coefficient (DSC) between two 
-    labeled images.
+def measure_overlap_labels(
+        img1: Union["sitk.Image", "itk.Image"],
+        img2: Union["sitk.Image", "itk.Image"]) -> float:
+    """Measure the mean Dice Similarity Coefficient between two labeled images.
     
     Args:
-        labels_img1 (:obj:`sitk.Image`): Labels image 1.
-        labels_img2 (:obj:`sitk.Image`): Labels image 2.
+        img1: Labels image 1.
+        img2: Labels image 2.
     
     Returns:
-        float: The mean label-by-label DSC of the two given images.
+        The mean label-by-label DSC of the two given images.
+    
     """
-    overlap_filter = sitk.LabelOverlapMeasuresImageFilter()
-    overlap_filter.Execute(labels_img1, labels_img2)
-    mean_region_dsc = overlap_filter.GetDiceCoefficient()
-    _logger.info("Mean label overlap DSC: {}".format(mean_region_dsc))
+    is_sitk = sitk and isinstance(img1, sitk.Image)
+    mean_region_dsc = np.nan
+    try:
+        if is_sitk:
+            overlap_filter = sitk.LabelOverlapMeasuresImageFilter()
+            overlap_filter.Execute(img1, img2)
+        
+        else:
+            # ITK filter appears to only take unsigned char images, which
+            # should be set on instantiation to set the filter type
+            # get ITK docs from:
+            # itk.itkLabelOverlapMeasuresImageFilterPython.label_overlap_measures_image_filter
+            src_img = img1.astype(itk.US)
+            tgt_img = img2.astype(itk.US)
+            overlap_filter = itk.LabelOverlapMeasuresImageFilter.New(src_img)
+            
+            # set both images again since the image set in New is unclear
+            overlap_filter.SetSourceImage(src_img)
+            overlap_filter.SetTargetImage(tgt_img)
+        
+        mean_region_dsc = overlap_filter.GetDiceCoefficient()
+        _logger.info("Mean label overlap DSC: %s\n", mean_region_dsc)
+    except RuntimeError as e:
+        _logger.warn(e)
     return mean_region_dsc
 
 
-def make_labels_fg(labels_sitk):
+def make_labels_fg(
+        labels_sitk: Union["sitk.Image", "itk.Image"]
+) -> Union["sitk.Image", "itk.Image"]:
     """Make a labels foreground image.
     
     Args:
-        labels_sitk (:obj:`sitk.Image`): Labels image where
-            0 = background, and all other values are considered foreground.
+        labels_sitk: Labels image where 0 = background, and all other values
+            are considered foreground.
     
     Returns:
-        :obj:`sitk.Image`: Labels foreground image.
+        Labels foreground image.
+    
     """
-    fg_img = sitk.GetArrayFromImage(labels_sitk)
+    fg_img = sitk_io.convert_img(labels_sitk)
     fg_img[fg_img != 0] = 1
     fg_img_sitk = sitk_io.replace_sitk_with_numpy(labels_sitk, fg_img)
     return fg_img_sitk
 
 
-def measure_overlap_combined_labels(base_img, labels_img, add_lbls=None,
-                                    return_masks=False):
+def measure_overlap_combined_labels(
+        base_img: Union["sitk.Image", "itk.Image"],
+        labels_img: Union["sitk.Image", "itk.Image"],
+        add_lbls: Optional[Sequence[Tuple[int, int]]] = None,
+        return_masks: bool = False) -> float:
     """Measures the overlap of a combined labels images to another image
 
     Should be 1.0 by definition when using ``labels_img`` to curate another
     image.
 
     Args:
-        base_img (:obj:`sitk.Image`): Base image.
-        labels_img (:obj:`sitk.Image`): Labels image, where all non-zero
+        base_img: Base image.
+        labels_img: Labels image, where all non-zero
             values will be treated as foreground.
-        add_lbls (List[int]): Sequence of labels in ``labels_img`` to add
-            to the mask for ``base_img``; defaults to None.
-        return_masks (bool): True to return the thresholded, mask images
-            used for the overlap calculate; defaults to False.
+        add_lbls: Sequence of labels in ``labels_img`` to add
+            to the mask for ``base_img``.
+        return_masks: True to return the thresholded, mask images
+            used for the overlap calculate.
 
     Returns:
-        float: Dice Similarity Coefficient from :meth:`measure_overlap`.
+        Dice Similarity Coefficient from :meth:`measure_overlap`.
 
     """
     lbls_fg = make_labels_fg(labels_img)
     mask = None
     if add_lbls is not None:
-        # build mask from labels to add to fixed image's foreground, such 
-        # as labeled ventricles; TODO: get children of labels rather than 
-        # taking labels range, but would need to load labels reference; 
-        # TODO: consider using "atlas_threshold_all" profile setting 
-        # instead, but would need to ensure that fixed thresholds work 
-        # for both atlas and sample histology
-        labels_np_abs = np.absolute(sitk.GetArrayFromImage(labels_img))
+        # build mask from labels to add to fixed image's foreground, such
+        # as labeled ventricles
+        # TODO: get children of labels rather than taking labels range,
+        #   but would need to load labels reference;
+        # TODO: consider using "atlas_threshold_all" profile setting
+        #   instead, but would need to ensure that fixed thresholds work
+        #   for both atlas and sample histology
+        labels_np_abs = np.absolute(sitk_io.convert_img(labels_img))
         mask = np.zeros_like(labels_np_abs, dtype=bool)
         for lbl in add_lbls:
             print("adding abs labels within", lbl)
-            mask[np.all([labels_np_abs >= lbl[0], 
+            mask[np.all([labels_np_abs >= lbl[0],
                          labels_np_abs < lbl[1]], axis=0)] = True
     print("DSC of thresholded fixed image compared with combined labels:")
     return measure_overlap(
diff --git a/magmap/atlas/edge_seg.py b/magmap/atlas/edge_seg.py
index 2bb84475..78904dca 100644
--- a/magmap/atlas/edge_seg.py
+++ b/magmap/atlas/edge_seg.py
@@ -1,14 +1,17 @@
 # Segmentation based on edge detection
-# Author: David Young, 2019
+# Author: David Young, 2019, 2023
 """Re-segment atlases based on edge detections.
 """
 import os
 from time import time
 from typing import List, Optional
 
-import SimpleITK as sitk
 import numpy as np
 import pandas as pd
+try:
+    import SimpleITK as sitk
+except ImportError:
+    sitk = None
 from skimage import color
 
 from magmap.atlas import atlas_refiner
@@ -111,18 +114,18 @@ def make_edge_images(path_img, show=True, atlas=True, suffix=None,
         path_atlas = path_img
         path_labels = os.path.join(path_atlas_dir, labels_suffix)
         print("loading labels from", path_labels)
-        labels_sitk = sitk.ReadImage(path_labels)
+        labels_sitk = sitk_io.read_img(path_labels)
     elif labels_suffix:
         # load labels registered to sample image
         labels_sitk = sitk_io.load_registered_img(
             mod_path, labels_suffix, get_sitk=True)
-    labels_img_np = None if labels_sitk is None else sitk.GetArrayFromImage(
+    labels_img_np = None if labels_sitk is None else sitk_io.convert_img(
         labels_sitk)
     # load atlas image, set resolution from it
     atlas_sitk = sitk_io.load_registered_img(
         path_atlas, atlas_suffix, get_sitk=True)
     config.resolutions = np.array([atlas_sitk.GetSpacing()[::-1]])
-    atlas_np = sitk.GetArrayFromImage(atlas_sitk)
+    atlas_np = sitk_io.convert_img(atlas_sitk)
     
     if config.rgb:
         # convert RGB atlas image to grayscale for single channel
@@ -182,7 +185,7 @@ def make_edge_images(path_img, show=True, atlas=True, suffix=None,
         config.RegNames.IMG_LABELS_INTERIOR.value: labels_sitk_interior, 
         config.RegNames.IMG_LABELS_DIST.value: dist_sitk, 
     }
-    if show:
+    if show and sitk:
         for img in imgs_write.values():
             if img: sitk.Show(img)
     
@@ -296,10 +299,10 @@ def edge_aware_segmentation(
         load_path, config.RegNames.IMG_LABELS_MARKERS.value, get_sitk=True)
     
     # get Numpy arrays of images
-    atlas_img_np = sitk.GetArrayFromImage(atlas_sitk)
-    atlas_edge = sitk.GetArrayFromImage(atlas_sitk_edge)
-    labels_img_np = sitk.GetArrayFromImage(labels_sitk)
-    markers = sitk.GetArrayFromImage(labels_sitk_markers)
+    atlas_img_np = sitk_io.convert_img(atlas_sitk)
+    atlas_edge = sitk_io.convert_img(atlas_sitk_edge)
+    labels_img_np = sitk_io.convert_img(labels_sitk)
+    markers = sitk_io.convert_img(labels_sitk_markers)
     
     # segment image from markers
     sym_axis = atlas_refiner.find_symmetric_axis(atlas_img_np)
@@ -393,7 +396,7 @@ def edge_aware_segmentation(
     # show and write image to same directory as atlas with appropriate suffix
     sitk_io.write_reg_images(
         {config.RegNames.IMG_LABELS.value: labels_sitk_seg}, mod_path)
-    if show: sitk.Show(labels_sitk_seg)
+    if show and sitk: sitk.Show(labels_sitk_seg)
     return path_atlas
 
 
@@ -433,7 +436,7 @@ def merge_atlas_segmentations(img_paths, show=True, atlas=True, suffix=None):
         if erode["markers"]:
             # use default minimal post-erosion size (not setting erosion frac)
             markers, df = erode_labels(
-                sitk.GetArrayFromImage(labels_sitk), erosion,
+                sitk_io.convert_img(labels_sitk), erosion,
                 mirrored=mirrored, mirror_mult=mirror_mult)
             labels_sitk_markers = sitk_io.replace_sitk_with_numpy(
                 labels_sitk, markers)
@@ -474,7 +477,7 @@ def merge_atlas_segmentations(img_paths, show=True, atlas=True, suffix=None):
         if meas_edge_dist or erode_interior:
             labels_sitk = sitk_io.load_registered_img(
                 mod_path, config.RegNames.IMG_LABELS.value, get_sitk=True)
-            labels_np = sitk.GetArrayFromImage(labels_sitk)
+            labels_np = sitk_io.convert_img(labels_sitk)
             if meas_edge_dist:
                 # make edge distance images and stats
                 dist_to_orig, labels_edge = edge_distances(
@@ -501,7 +504,7 @@ def merge_atlas_segmentations(img_paths, show=True, atlas=True, suffix=None):
             config.RegNames.IMG_LABELS_INTERIOR.value: labels_sitk_interior, 
         }
         sitk_io.write_reg_images(imgs_write, mod_path)
-        if show:
+        if show and sitk:
             for img in imgs_write.values():
                 if img: sitk.Show(img)
         print("finished {}".format(path))
@@ -570,7 +573,7 @@ def make_sub_segmented_labels(img_path, suffix=None):
         img_path, config.RegNames.IMG_ATLAS_EDGE.value)
     
     # sub-divide the labels and save to file
-    labels_img_np = sitk.GetArrayFromImage(labels_sitk)
+    labels_img_np = sitk_io.convert_img(labels_sitk)
     labels_subseg = segmenter.sub_segment_labels(labels_img_np, atlas_edge)
     labels_subseg_sitk = sitk_io.replace_sitk_with_numpy(
         labels_sitk, labels_subseg)
diff --git a/magmap/atlas/register.py b/magmap/atlas/register.py
index c6281d10..6fe6e619 100644
--- a/magmap/atlas/register.py
+++ b/magmap/atlas/register.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 # Image registration
-# Author: David Young, 2017, 2019
+# Author: David Young, 2017, 2023
 """Register images to one another.
 
 The registration type can be given on the command-line (see :mod:``cli``) as 
@@ -43,11 +43,18 @@
 import os
 import shutil
 from time import time
-from typing import Any, Dict, List, Optional, Sequence, Tuple, Union
+from typing import Any, Callable, Dict, List, Optional, Sequence, Tuple, Union
 
+try:
+    import itk
+except ImportError:
+    itk = None
 import pandas as pd
 import numpy as np
-import SimpleITK as sitk
+try:
+    import SimpleITK as sitk
+except ImportError:
+    sitk = None
 from skimage import filters, measure, morphology, transform
 
 from magmap.atlas import atlas_refiner, edge_seg, labels_meta, ontology, \
@@ -79,12 +86,23 @@
 @dataclasses.dataclass
 class RegImgs:
     """Data class for tracking registered images."""
-    exp_orig: Optional[Union[np.ndarray, sitk.Image]] = None
-    exp: Optional[Union[np.ndarray, sitk.Image]] = None
-    atlas: Optional[Union[np.ndarray, sitk.Image]] = None
-    labels: Optional[Union[np.ndarray, sitk.Image]] = None
-    labels_markers: Optional[Union[np.ndarray, sitk.Image]] = None
-    borders: Optional[Union[np.ndarray, sitk.Image]] = None
+    #: Original experimental/fixed image.
+    exp_orig: Optional[Union[np.ndarray, "sitk.Image", "itk.Image"]] = None
+    #: Experimental/fixed image.
+    exp: Optional[Union[np.ndarray, "sitk.Image", "itk.Image"]] = None
+    #: Atlas/moving image.
+    atlas: Optional[Union[np.ndarray, "sitk.Image", "itk.Image"]] = None
+    #: Atlas/moving labels image.
+    labels: Optional[Union[np.ndarray, "sitk.Image", "itk.Image"]] = None
+    #: Labels markers image.
+    labels_markers: Optional[Union[
+        np.ndarray, "sitk.Image", "itk.Image"]] = None
+    #: Borders image.
+    borders: Optional[Union[np.ndarray, "sitk.Image", "itk.Image"]] = None
+    #: Experimental/fixed mask image.
+    exp_mask: Optional[Union[np.ndarray, "sitk.Image", "itk.Image"]] = None
+    #: Atlas/moving mask image.
+    atlas_mask: Optional[Union[np.ndarray, "sitk.Image", "itk.Image"]] = None
     
     @staticmethod
     def get_order(name: str) -> Optional[int]:
@@ -98,6 +116,22 @@ def get_order(name: str) -> Optional[int]:
 
         """
         return 0 if name in ("labels", "labels_markers", "borders") else None
+    
+    def update_fields(
+            self, fn: Callable[
+                [Optional[Union[np.ndarray, "sitk.Image", "itk.Image"]]],
+                Optional[Union[np.ndarray, "sitk.Image", "itk.Image"]]]):
+        """Update all fields.
+        
+        Args:
+            fn: Function to apply to each field if not None.
+
+        """
+        for field in dataclasses.fields(RegImgs):
+            reg_img = getattr(self, field.name)
+            if reg_img is not None:
+                reg_img = fn(reg_img)
+            setattr(self, field.name, reg_img)
 
 
 def _translation_adjust(orig, transformed, translation, flip=False):
@@ -157,13 +191,21 @@ def _show_overlays(imgs, translation, fixed_file, plane):
 
 
 def _handle_transform_file(fixed_file, transform_param_map=None):
+    # read or write transformation file
     base_name = sitk_io.reg_out_path(fixed_file, "")
     filename = base_name + "transform.txt"
-    param_map = None
     if transform_param_map is None:
-        param_map = sitk.ReadParameterFile(filename)
+        if sitk:
+            param_map = sitk.ReadParameterFile(filename)
+        else:
+            param_map = itk.ParameterObject.New()
+            param_map.ReadParameterFile(filename)
     else:
-        sitk.WriteParameterFile(transform_param_map[0], filename)
+        if sitk:
+            sitk.WriteParameterFile(transform_param_map[0], filename)
+        else:
+            param = itk.ParameterObject.New()
+            param.WriteParameterFile(transform_param_map[0], filename)
         param_map = transform_param_map[0]
     return param_map, None # TODO: not using translation parameters
     transform = np.array(param_map["TransformParameters"]).astype(float)
@@ -179,13 +221,14 @@ def _handle_transform_file(fixed_file, transform_param_map=None):
     else:
         print("Transform parameters do not match scaling dimensions")
     return param_map, translation
+    return param_map, None  # TODO: not using translation parameters
 
 
 def curate_img(
-        fixed_img: sitk.Image, labels_img: sitk.Image,
-        imgs: Optional[List[sitk.Image]] = None, inpaint: bool = True,
+        fixed_img: "sitk.Image", labels_img: "sitk.Image",
+        imgs: Optional[List["sitk.Image"]] = None, inpaint: bool = True,
         carve: bool = True, thresh: Optional[float] = None,
-        holes_area: Optional[int] = None) -> List[sitk.Image]:
+        holes_area: Optional[int] = None) -> List["sitk.Image"]:
     """Curate an image by the foreground of another.
     
     In-painting where corresponding pixels are present in fixed image but
@@ -210,8 +253,8 @@ def curate_img(
         with the curated ``labels_img``, followed by the images in ``imgs``.
     
     """
-    fixed_img_np = sitk.GetArrayFromImage(fixed_img)
-    labels_img_np = sitk.GetArrayFromImage(labels_img)
+    fixed_img_np = sitk_io.convert_img(fixed_img)
+    labels_img_np = sitk_io.convert_img(labels_img)
     
     # ensure that labels image is first
     if imgs:
@@ -234,7 +277,7 @@ def curate_img(
         # in-paint and remove pixels from result image where fixed image is
         # below threshold
         img = imgs[i]
-        result_img_np = sitk.GetArrayFromImage(img)
+        result_img_np = sitk_io.convert_img(img)
         if inpaint:
             result_img_np = cv_nd.in_paint(result_img_np, to_fill)
         if carve:
@@ -253,17 +296,19 @@ def curate_img(
 
 
 def register_repeat(
-        transformix_img_filter: sitk.TransformixImageFilter,
-        img: sitk.Image, preserve_idents: bool = False) -> sitk.Image:
+        transformix_img_filter: Union[
+            "sitk.TransformixImageFilter", "itk.TransformixFilter"],
+        img: Union["sitk.Image", "itk.Image"], preserve_idents: bool = False
+) -> Union["sitk.Image", "itk.Image"]:
     """Transform labels to match a prior registration.
     
     Uses an Elastix Transformix filter to reproduce a transformation on
-    a labels image. Ensures that the output pixel type remains the same as
-    the input.
+    a labels image. This filter can be from SimpleITK or ITK. Ensures that
+    the output pixel type remains the same as the input.
     
     Args:
         transformix_img_filter: Filter generated from a prior registration.
-        img: SimpleITK image.
+        img: SimpleITK or ITK image.
         preserve_idents: True to ensure that identities in ``img`` are
             preserved. Typically used for label images. Defaults to False.
 
@@ -271,30 +316,64 @@ def register_repeat(
         Transformed image.
 
     """
-    # apply atlas transformation to labels image
-    pixel_id = img.GetPixelID()  # now as signed int
+    # use SimpleITK if it has Elastix and filter is from it
+    is_sitk = sitk and "TransformixImageFilter" in dir(sitk) and isinstance(
+        transformix_img_filter, sitk.TransformixImageFilter)
+    
+    # check if input image is sitk type
+    is_sitk_img = sitk and isinstance(img, sitk.Image)
+    
+    if is_sitk:
+        # store data type
+        pixel_id = img.GetPixelID()
+    else:
+        if is_sitk_img:
+            # convert to ITK Image
+            img = sitk_io.sitk_to_itk_img(img)
+        
+        # cast to required data type for ITK transformation
+        pixel_id = type(img)
+        img = img.astype(itk.F)
     
     img_sitk = img
     img_unique = None
     if preserve_idents:
         # map values to indices in their unique array to minimize change of
         # rounding errors for large values
-        img_np = sitk.GetArrayFromImage(img_sitk)
+        img_np = sitk_io.convert_img(img_sitk)
         img_unique, img_inds = np.unique(img_np, return_inverse=True)
         
         # offset by 1 since background will be inserted as 0
         img_inds = img_inds.reshape(img_np.shape) + 1
-        img_sitk = sitk_io.replace_sitk_with_numpy(img, img_inds)
+        img_sitk = sitk_io.replace_sitk_with_numpy(
+            img, img_inds.astype(np.float32))
+    
+    if is_sitk:
+        # reapply transformation using sitk
+        transformix_img_filter.SetMovingImage(img_sitk)
+        transformix_img_filter.Execute()
+        transf_img = transformix_img_filter.GetResultImage()
+        transf_img = sitk.Cast(transf_img, pixel_id)
+    else:
+        # reapply transformation using ITK
+        
+        # set up new filter since it needs to be set up with the image
+        # to support 3D images
+        params = transformix_img_filter.GetTransformParameterObject()
+        transformix_img_filter = itk.TransformixFilter.New(img_sitk)
+        transformix_img_filter.SetTransformParameterObject(params)
+        
+        # perform the transformation
+        transformix_img_filter.UpdateLargestPossibleRegion()
+        transf_img = transformix_img_filter.GetOutput()
+        cast_filter = itk.CastImageFilter[type(transf_img), pixel_id].New()
+        cast_filter.SetInput(transf_img)
+        transf_img = cast_filter.GetOutput()
     
-    # apply transformation
-    transformix_img_filter.SetMovingImage(img_sitk)
-    transformix_img_filter.Execute()
-    transf_img = transformix_img_filter.GetResultImage()
-    transf_img = sitk.Cast(transf_img, pixel_id)
-
     if preserve_idents:
         # map indices back to original values
-        transf_inds = sitk.GetArrayFromImage(transf_img)
+        transf_inds = sitk_io.convert_img(transf_img)
+        transf_inds = transf_inds.astype(int)
         bkgdi = np.nonzero(img_unique == 0)
         if len(bkgdi) > 0:
             # convert inserted background (0) to background index
@@ -312,6 +391,9 @@ def register_repeat(
     variance = LabelStatistics.GetVariance(1)
     print("count: {}, mean: {}, variance: {}".format(count, mean, variance))
     '''
+    
+    if not is_sitk and is_sitk_img:
+        transf_img = sitk_io.itk_to_sitk_img(transf_img)
     return transf_img
 
 
@@ -331,25 +413,32 @@ def _config_reg_resolutions(grid_spacing_schedule, param_map, ndim):
 
 
 def register_duo(
-        fixed_img: sitk.Image, moving_img: sitk.Image,
-        path: Optional[str] = None, fixed_mask: Optional[sitk.Image] = None,
-        moving_mask: Optional[sitk.Image] = None,
+        fixed_img: Union["sitk.Image", "itk.Image"],
+        moving_img: Union["sitk.Image", "itk.Image"],
+        path: Optional[str] = None,
+        fixed_mask: Optional[Union["sitk.Image", "itk.Image"]] = None,
+        moving_mask: Optional[Union["sitk.Image", "itk.Image"]] = None,
         regs: Optional[Union[
             Sequence[str], Sequence["atlas_prof.RegParamMap"]]] = None
-) -> Tuple[sitk.Image, sitk.TransformixImageFilter]:
-    """Register two images to one another using ``SimpleElastix``.
+) -> Tuple[Union["sitk.Image", "itk.Image"],
+           Union["sitk.TransformixImageFilter", "itk.TransformixFilter"]]:
+    """Register two images to one another using ``Elastix``.
+    
+    Supports Elastix provided by SimpleITK or ITK. If SimpleITK images are
+    given with Elastix enabled in SimpleITK, registration will be performed
+    in SimpleITK. If Elastix is not enabled, however, the images will be
+    converted to ITK, and registration will use ITK-Elastix. If ITK images
+    are given, ITK-Elastix is assumed to be present; no conversion will take
+    place to SimpleITK.
     
     Args:
         fixed_img: The image to be registered to.
         moving_img: The image to register to ``fixed_img``.
         path: Path as string from whose parent directory the points-based
             registration files ``fix_pts.txt`` and ``mov_pts.txt`` will
-            be found; defaults to None, in which case points-based
-            reg will be ignored even if set.
-        fixed_mask: Mask for ``fixed_img``, typically a uint8 image.; defaults
-            to None.
-        moving_mask: Mask for ``moving_img``, typically a uint8 image.; defaults
-            to None.
+            be found. If None, points-based reg will be ignored even if set.
+        fixed_mask: Mask for ``fixed_img``, typically a uint8 image.
+        moving_mask: Mask for ``moving_img``, typically a uint8 image.
         regs: Sequence of atlas profile registration keys or registration
             parameter objects. The default of None gives all three major
             registration types, "reg_translation", "reg_affine", "reg_bspline".
@@ -363,29 +452,83 @@ def register_duo(
         # default to perform all the major registration types
         regs = ("reg_translation", "reg_affine", "reg_bspline")
     
-    # basic info from images just prior to SimpleElastix filtering for
-    # registration; to view raw images, show these images rather than merely
-    # turning all iterations to 0 since simply running through the filter
-    # will alter images
-    print("fixed image (type {}):\n{}".format(
-        fixed_img.GetPixelIDTypeAsString(), fixed_img))
-    print("moving image (type {}):\n{}".format(
-        moving_img.GetPixelIDTypeAsString(), moving_img))
-    
-    # set up SimpleElastix filter
-    elastix_img_filter = sitk.ElastixImageFilter()
-    elastix_img_filter.SetFixedImage(fixed_img)
-    elastix_img_filter.SetMovingImage(moving_img)
-
-    # add any masks
-    if fixed_mask is not None:
-        elastix_img_filter.SetFixedMask(fixed_mask)
-    if moving_mask is not None:
-        elastix_img_filter.SetMovingMask(moving_mask)
+    # collect images
+    reg_imgs = RegImgs(
+        exp=fixed_img, atlas=moving_img, exp_mask=fixed_mask,
+        atlas_mask=moving_mask)
+    
+    if not itk and not sitk:
+        raise ImportError(
+            config.format_import_err(
+                "itk-elastix", "ITK-Elastix or SimpleITK with Elastix",
+                "image registration"))
+    
+    is_sitk_img = sitk and isinstance(reg_imgs.exp, sitk.Image)
+    if itk and isinstance(reg_imgs.exp, itk.Image):
+        # use ITK-Elastix if images are of ITK type
+        is_sitk = False
+        
+    elif is_sitk_img:
+        # images are of SimpleITK type
+        if "ElastixImageFilter" in dir(sitk):
+            # use SimpleITK if Elastix is enabled
+            is_sitk = True
+        else:
+            # fall back to ITK-Elastix
+            is_sitk = False
+            _logger.debug(
+                "Converting SimpleITK to ITK images since Elastix was not "
+                "found in the SimpleITK library")
+            
+            # convert sitk to ITK images
+            reg_imgs.update_fields(
+                lambda x: sitk_io.sitk_to_itk_img(x))
+        
+    else:
+        raise TypeError(
+            f"Images must be ITK or SimpleITK Image types, but 'fixed_img' is "
+            f"of type {type(reg_imgs.exp)}")
+    
+    _logger.info(f"Fixed image:\n{reg_imgs.exp}")
+    _logger.info(f"Moving image:\n{reg_imgs.atlas}")
+    
+    param_map_vector = None  # for sitk
+    reg_params = None  # for ITK
+    if is_sitk:
+        # set up SimpleElastix filter
+        _logger.info("Registering images using SimpleITK with Elastix")
+        elastix_img_filter = sitk.ElastixImageFilter()
+        elastix_img_filter.SetFixedImage(reg_imgs.exp)
+        elastix_img_filter.SetMovingImage(reg_imgs.atlas)
+        param_map_vector = sitk.VectorOfParameterMap()
+        
+        # add any masks
+        if reg_imgs.exp_mask is not None:
+            elastix_img_filter.SetFixedMask(reg_imgs.exp_mask)
+        if reg_imgs.atlas_mask is not None:
+            elastix_img_filter.SetMovingMask(reg_imgs.atlas_mask)
+        
+    else:
+        # set up ITK-Elastix filter
+        _logger.info("Registering images using ITK-Elastix")
+        
+        # main images must be float
+        reg_imgs.exp = reg_imgs.exp.astype(itk.F)
+        reg_imgs.atlas = reg_imgs.atlas.astype(itk.F)
+        elastix_img_filter = itk.ElastixRegistrationMethod.New(
+            reg_imgs.exp, reg_imgs.atlas)
+        reg_params = itk.ParameterObject.New()
+        
+        # masks must be unsigned char
+        if reg_imgs.exp_mask:
+            reg_imgs.exp_mask = reg_imgs.exp_mask.astype(itk.UC)
+            elastix_img_filter.SetFixedMask(reg_imgs.exp_mask)
+        if reg_imgs.atlas_mask:
+            reg_imgs.atlas_mask = reg_imgs.atlas_mask.astype(itk.UC)
+            elastix_img_filter.SetMovingMask(reg_imgs.atlas_mask)
     
     # set up parameter maps for the included registration types
     settings = config.atlas_profile
-    param_map_vector = sitk.VectorOfParameterMap()
     for reg in regs:
         # get registration parameters from profile
         params = reg if isinstance(
@@ -393,10 +536,18 @@ def register_duo(
         if not params: continue
         
         max_iter = params["max_iter"]
-        # TODO: consider removing since does not skip if "0" and need at least
-        # one transformation for reg, even if 0 iterations
+        # TODO: Consider removing since does not skip if "0" and need at least
+        #   one transformation for reg, even if 0 iterations. Also, note that
+        #   turning all iterations to 0 since simply running through the filter
+        #   will still alter images.
         if not max_iter: continue
-        param_map = sitk.GetDefaultParameterMap(params["map_name"])
+        
+        reg_name = params["map_name"]
+        if reg_params:
+            param_map = reg_params.GetDefaultParameterMap(reg_name)
+        else:
+            param_map = sitk.GetDefaultParameterMap(reg_name)
+        
         similarity = params["metric_similarity"]
         if len(param_map["Metric"]) > 1:
             param_map["Metric"] = [similarity, *param_map["Metric"][1:]]
@@ -407,7 +558,7 @@ def register_duo(
         if grid_spacing_sched:
             # fine tune the spacing for multi-resolution registration
             _config_reg_resolutions(
-                grid_spacing_sched, param_map, fixed_img.GetDimension())
+                grid_spacing_sched, param_map, reg_imgs.exp.GetDimension())
         else:
             # num of resolutions is automatically set by spacing sched
             param_map["NumberOfResolutions"] = [params["num_resolutions"]]
@@ -429,27 +580,54 @@ def register_duo(
             fix_pts_path = os.path.join(os.path.dirname(path), "fix_pts.txt")
             move_pts_path = os.path.join(
                 os.path.dirname(path), "mov_pts.txt")
+            
             if os.path.isfile(fix_pts_path) and os.path.isfile(move_pts_path):
+                # set point files if both exist
                 metric = list(param_map["Metric"])
                 metric.append("CorrespondingPointsEuclideanDistanceMetric")
                 param_map["Metric"] = metric
                 #param_map["Metric2Weight"] = ["0.5"]
                 elastix_img_filter.SetFixedPointSetFileName(fix_pts_path)
                 elastix_img_filter.SetMovingPointSetFileName(move_pts_path)
-        param_map_vector.append(param_map)
+        
+        if param_map_vector is not None:
+            param_map_vector.append(param_map)
+        elif reg_params is not None:
+            reg_params.AddParameterMap(param_map)
+    
+    if is_sitk:
+        # perform registration in sitk
+        elastix_img_filter.SetParameterMap(param_map_vector)
+        elastix_img_filter.PrintParameterMap()
+        elastix_img_filter.Execute()
+        transformed_img = elastix_img_filter.GetResultImage()
+        
+        # prep filter to apply transformation to label files; turn off final
+        # bspline interpolation to avoid overshooting the interpolation for the
+        # labels image (see Elastix manual section 4.3)
+        transform_param_map = elastix_img_filter.GetTransformParameterMap()
+        transformix_img_filter = sitk.TransformixImageFilter()
+        transform_param_map[-1]["FinalBSplineInterpolationOrder"] = ["0"]
+        transformix_img_filter.SetTransformParameterMap(transform_param_map)
 
-    elastix_img_filter.SetParameterMap(param_map_vector)
-    elastix_img_filter.PrintParameterMap()
-    elastix_img_filter.Execute()
-    transformed_img = elastix_img_filter.GetResultImage()
-    
-    # prep filter to apply transformation to label files; turn off final
-    # bspline interpolation to avoid overshooting the interpolation for the
-    # labels image (see Elastix manual section 4.3)
-    transform_param_map = elastix_img_filter.GetTransformParameterMap()
-    transformix_img_filter = sitk.TransformixImageFilter()
-    transform_param_map[-1]["FinalBSplineInterpolationOrder"] = ["0"]
-    transformix_img_filter.SetTransformParameterMap(transform_param_map)
+    else:
+        # perform registration in ITK-Elastix
+        elastix_img_filter.SetParameterObject(reg_params)
+        elastix_img_filter.SetLogToConsole(config.verbose)
+        elastix_img_filter.UpdateLargestPossibleRegion()
+        transformed_img = elastix_img_filter.GetOutput()
+        transf_params = elastix_img_filter.GetTransformParameterObject()
+        
+        # set up transformix
+        transf_params.SetParameter(
+            transf_params.GetNumberOfParameterMaps() - 1,
+            "FinalBSplineInterpolationOrder", "0")
+        transformix_img_filter = itk.TransformixFilter.New()
+        transformix_img_filter.SetTransformParameterObject(transf_params)
+        
+        if is_sitk_img:
+            # convert back to sitk Image
+            transformed_img = sitk_io.itk_to_sitk_img(transformed_img)
     
     return transformed_img, transformix_img_filter
 
@@ -458,9 +636,10 @@ def register(
         fixed_file: str, moving_img_path: str, show_imgs: bool = True,
         write_imgs: bool = True, name_prefix: Optional[str] = None,
         new_atlas: bool = False,
-        transformix: Optional["sitk.TransformixImageFilter"] = None
-) -> "sitk.TransformixImageFilter":
-    """Register an atlas to a sample image using the SimpleElastix library.
+        transformix: Optional[Union[
+            "sitk.TransformixImageFilter", "itk.TransformixFilter"]] = None
+) -> Union["sitk.TransformixImageFilter", "itk.TransformixFilter"]:
+    """Register an atlas to a sample image using the Elastix library.
     
     Loads the images, applies any transformations to the moving image, and
     registers the moving to the sample images. Applies the identical
@@ -481,13 +660,13 @@ def register(
         * If the registration does not complete, the moving image origin and
           direction will be changed to that of the fixed image
         * If no reg, the spacing will also be matched
-        * If still no reg, the atlas will be output as-is 
+        * If still no reg, the atlas will be output as-is
         * If the atlas profile ``fallback`` parameter sets a DSC threshold
           and alternate registration settings, the image will be re-registered
           under these new settings if below the threshold
     
     Args:
-        fixed_file: The image to register, given as a Numpy archive file to 
+        fixed_file: The image to register, given as a Numpy archive file to
             be read by :importer:`read_file`.
         moving_img_path (str): Moving image base path from which an intensity
             and a labels image will be loaded using registered image suffixes.
@@ -495,18 +674,18 @@ def register(
             for the labels image. The intensity image will be used for
             registration. The atlas is currently used as the moving file
             since it is likely to be lower resolution than the Numpy file.
-        show_imgs: True if the output images should be displayed; defaults to 
+        show_imgs: True if the output images should be displayed; defaults to
             True.
-        write_imgs: True if the images should be written to file; defaults to 
+        write_imgs: True if the images should be written to file; defaults to
             False.
-        name_prefix: Path with base name where registered files will be output; 
+        name_prefix: Path with base name where registered files will be output;
             defaults to None, in which case the fixed_file path will be used.
-        new_atlas: True to generate registered images that will serve as a 
+        new_atlas: True to generate registered images that will serve as a
             new atlas; defaults to False.
-        transformix: SimpleElastix transformation filter from prior run;
+        transformix: Elastix transformation filter from prior run;
             defaults to None.
     Returns:
-        The SimpleElastix transformation filter generated by the registration.
+        The Elastix transformation filter generated by the registration.
     
     """
     def get_similarity_metric():
@@ -544,7 +723,7 @@ def get_similarity_metric():
     # TODO: assume fixed image is preprocessed before starting this reg?
     fixed_img_orig = fixed_img
     if settings["preprocess"]:
-        img_np = sitk.GetArrayFromImage(fixed_img)
+        img_np = sitk_io.convert_img(fixed_img)
         #img_np = plot_3d.saturate_roi(img_np)
         img_np = plot_3d.denoise_roi(img_np)
         fixed_img = sitk_io.replace_sitk_with_numpy(fixed_img, img_np)
@@ -591,7 +770,7 @@ def get_similarity_metric():
     truncate_labels = settings["truncate_labels"]
     if truncate_labels is not None:
         # generate a truncated/cropped version of the labels image
-        labels_trunc_np = sitk.GetArrayFromImage(labels_img)
+        labels_trunc_np = sitk_io.convert_img(labels_img)
         atlas_refiner.truncate_labels(labels_trunc_np, *truncate_labels)
         labels_trunc = sitk_io.replace_sitk_with_numpy(
             labels_img, labels_trunc_np)
@@ -613,11 +792,11 @@ def get_similarity_metric():
         labels_img, target_size_res=rescale, rotate_deg=rotate_deg, order=0)
 
     # get Numpy arrays of moving images for preprocessing
-    moving_img_np = sitk.GetArrayFromImage(moving_img)
-    labels_img_np = sitk.GetArrayFromImage(labels_img)
+    moving_img_np = sitk_io.convert_img(moving_img)
+    labels_img_np = sitk_io.convert_img(labels_img)
     moving_mask_np = None
     if moving_mask is not None:
-        moving_mask_np = sitk.GetArrayFromImage(moving_mask)
+        moving_mask_np = sitk_io.convert_img(moving_mask)
 
     crop_out_labels = config.atlas_profile["crop_out_labels"]
     if crop_out_labels is not None:
@@ -772,7 +951,7 @@ def make_labels(lbls_img):
             imgs_write[suffix] = make_labels(moving_imgs[suffix])[0]
             print("transforming label", suffix, imgs_write[suffix])
 
-    if show_imgs:
+    if show_imgs and sitk:
         # show individual SimpleITK images in default viewer
         for img in imgs_write.values():
             if img is not None:
@@ -797,19 +976,12 @@ def make_labels(lbls_img):
             lbls_meta.prefix = name_prefix
             lbls_meta.save()
 
-    # save transform parameters and attempt to find the original position 
-    # that corresponds to the final position that will be displayed
-    if transformix_filter is not None:
-        _, translation = _handle_transform_file(
-            name_prefix, transformix_filter.GetTransformParameterMap())
-        _translation_adjust(moving_img, img_moved, translation, flip=True)
-    
     # compare original atlas with registered labels taken as a whole
     dsc_labels = atlas_refiner.measure_overlap_combined_labels(
         fixed_img_orig, labels_moved)
     
     # measure compactness of fixed image
-    fixed_img_orig_np = sitk.GetArrayFromImage(fixed_img_orig)
+    fixed_img_orig_np = sitk_io.convert_img(fixed_img_orig)
     thresh_atlas = fixed_img_orig_np > filters.threshold_mean(fixed_img_orig_np)
     compactness, _, _ = cv_nd.compactness_3d(
         thresh_atlas, fixed_img_orig.GetSpacing()[::-1])
@@ -832,6 +1004,13 @@ def make_labels(lbls_img):
     df_io.dict_to_data_frame(metrics, df_path, show="\t")
 
     '''
+    # save transform parameters and attempt to find the original position
+    # that corresponds to the final position that will be displayed
+    if transformix_filter is not None:
+        _, translation = _handle_transform_file(
+            name_prefix, transformix_filter.GetTransformParameterMap())
+        _translation_adjust(moving_img, img_moved, translation, flip=True)
+    
     # show 2D overlays or registered image and atlas last since blocks until 
     # fig is closed
     imgs = [
@@ -886,14 +1065,14 @@ def register_rev(fixed_path, moving_path, reg_base=None, reg_names=None,
             defaults to None.
         show: True to show images after registration; defaults to True.
     """
-    fixed_img = sitk.ReadImage(fixed_path)
+    fixed_img = sitk_io.read_img(fixed_path)
     mod_path = moving_path
     if suffix is not None:
         # adjust image path to load with suffix
         mod_path = libmag.insert_before_ext(mod_path, suffix)
     if reg_base is None:
         # load the image directly from given path
-        moving_img = sitk.ReadImage(mod_path)
+        moving_img = sitk_io.read_img(mod_path)
     else:
         # treat the path as a base path to which a reg suffix will be combined
         moving_img = sitk_io.load_registered_img(
@@ -941,7 +1120,7 @@ def register_rev(fixed_path, moving_path, reg_base=None, reg_names=None,
         # distinct name to avoid overwriting previously registered images
         imgs_write[name] = img
     sitk_io.write_reg_images(imgs_write, output_base)
-    if show:
+    if show and sitk:
         for img in imgs_write.values(): sitk.Show(img)
 
 
@@ -1007,20 +1186,21 @@ def register_group(
     Registration parameters are assumed to be in a "b-spline"
     :class:`magmap.settings.atlas_prof.RegParamMap`.
     
+    Performs registration using SimpleITK if the package is present, otherwise
+    using ITK-Elastix.
+    
     Args:
         img_files: Paths to image files to register. A minimum of 4 images
             is required for groupwise registration.
         rotate: List of number of 90 degree rotations for images
             corresponding to ``img_files``; defaults to None, in which
             case the `config.transform` rotate attribute will be used.
-        show_imgs: True if the output images should be displayed; defaults to
-            True.
-        write_imgs: True if the images should be written to file; defaults to
-            True.
+        show_imgs: True if the output images should be displayed.
+        write_imgs: True if the images should be written to file.
         name_prefix: Path with base name where registered files will be output;
             defaults to None, in which case the fixed_file path will be used.
         scale: Rescaling factor as a scalar value, used to find the rescaled,
-            smaller images corresponding to ``img_files``. Defaults to None.
+            smaller images corresponding to ``img_files``.
     
     """
     start_time = time()
@@ -1030,23 +1210,26 @@ def register_group(
         rotate = config.transform[config.Transforms.ROTATE]
     target_size = config.atlas_profile["target_size"]
     
-    '''
-    # TESTING: assuming first file is a raw groupwise registered image,
-    # import it for post-processing
-    img = sitk.ReadImage(img_files[0])
-    img_np = sitk.GetArrayFromImage(img)
-    print("thresh mean: {}".format(filters.threshold_mean(img_np)))
-    carve_threshold = config.register_settings["carve_threshold"]
-    holes_area = config.register_settings["holes_area"]
-    img_np, img_np_unfilled = plot_3d.carve(
-        img_np, thresh=carve_threshold, holes_area=holes_area,
-        return_unfilled=True)
-    sitk.Show(sitk_io.replace_sitk_with_numpy(img, img_np_unfilled))
-    sitk.Show(sitk_io.replace_sitk_with_numpy(img, img_np))
-    return
-    '''
+    if not itk and not sitk:
+        raise ImportError(
+            config.format_import_err(
+                "itk-elastix", "ITK-Elastix or SimpleITK with Elastix",
+                "groupwise image registration"))
+    
+    imgs = None
+    img_vector = None  # for sitk
+    reg_params = None  # for ITK
+    if sitk:
+        # register with SimpleITK if present
+        _logger.info("Groupwise registering images using SimpleITK with Elastix")
+        img_vector = sitk.VectorOfImage()
+        
+    else:
+        # set up ITK-Elastix parameters
+        _logger.info("Groupwise registering images using ITK-Elastix")
+        reg_params = itk.ParameterObject.New()
+        imgs = []
     
-    img_vector = sitk.VectorOfImage()
     # image properties of 1st image, in SimpleITK format
     origin = None
     size_orig = None
@@ -1062,8 +1245,8 @@ def register_group(
         rot = rotate and libmag.get_if_within(rotate, i, 0) >= 2
         chl = config.channel[0] if config.channel else 0
         img = sitk_io.load_numpy_to_sitk(img_file, rot, chl)
-        size = img.GetSize()
-        img_np = sitk.GetArrayFromImage(img)
+        size = img.GetSize() if sitk else itk.size(img)
+        img_np = sitk_io.convert_img(img)
         if img_np_template is None:
             img_np_template = np.copy(img_np)
 
@@ -1097,7 +1280,7 @@ def register_group(
         if origin is None:
             origin = img.GetOrigin()
             size_orig = size
-            size_cropped = img.GetSize()
+            size_cropped = img.GetSize() if sitk else tuple(itk.size(img))
             spacing = img.GetSpacing()
             start_y = y_cropped
             print("size_cropped: ", size_cropped, ", size_orig", size_orig)
@@ -1108,52 +1291,91 @@ def register_group(
             # favor of comparing shapes of large original and registered images
             img.SetOrigin(origin)
             img.SetSpacing(spacing)
-        print("img_file: {}\n{}".format(img_file, img))
-        img_vector.push_back(img)
-        #sitk.Show(img)
+        
+        if img_vector is not None:
+            img_vector.push_back(img)
+        else:
+            imgs.append(img)
+    
     #sitk.ProcessObject.SetGlobalDefaultDirectionTolerance(1)
     #sitk.ProcessObject.SetGlobalDefaultCoordinateTolerance(100)
-    img_combined = sitk.JoinSeries(img_vector)
     
-    # add b-spline registration parameter map
+    if img_vector is not None:
+        # sitk provides an array to register images of different sizes
+        img_combined = sitk.JoinSeries(img_vector)
+        elastix_img_filter = sitk.ElastixImageFilter()
+        elastix_img_filter.SetFixedImage(img_combined)
+        elastix_img_filter.SetMovingImage(img_combined)
+    else:
+        # assume all images are of the same size for ITK
+        imgs_np = [sitk_io.convert_img(m) for m in imgs]
+        imgs_np = np.array(imgs_np)
+        img_combined = sitk_io.convert_img(imgs_np, False).astype(itk.F)
+        sitk_io.match_world_info(imgs[0], img_combined)
+        elastix_img_filter = itk.ElastixRegistrationMethod.New(
+            img_combined, img_combined)
+    
+    _logger.info("Images to groupwise register:\n%s", img_combined)
+    
+    # get custom parameters from b-spline registration profile
     settings = config.atlas_profile
     reg = settings["reg_bspline"]
-    elastix_img_filter = sitk.ElastixImageFilter()
-    elastix_img_filter.SetFixedImage(img_combined)
-    elastix_img_filter.SetMovingImage(img_combined)
-    param_map = sitk.GetDefaultParameterMap("groupwise")
-    param_map["FinalGridSpacingInVoxels"] = [
-        reg["grid_space_voxels"]]
+    
+    # get default parameters for groupwise registeration
+    reg_name = "groupwise"
+    if reg_params:
+        param_map = reg_params.GetDefaultParameterMap(reg_name)
+    else:
+        param_map = sitk.GetDefaultParameterMap(reg_name)
+    
+    # change spacing to voxels
+    param_map["FinalGridSpacingInVoxels"] = [reg["grid_space_voxels"]]
     del param_map["FinalGridSpacingInPhysicalUnits"]  # avoid conflict with vox
+    
+    # set iterations
     param_map["MaximumNumberOfIterations"] = [reg["max_iter"]]
     # TESTING:
-    #param_map["MaximumNumberOfIterations"] = ["0"]
+    # param_map["MaximumNumberOfIterations"] = ["0"]
+    
+    # add a set of axis resolutions for each resolution level
     _config_reg_resolutions(
         reg["grid_spacing_schedule"], param_map, img_np_template.ndim)
-    elastix_img_filter.SetParameterMap(param_map)
-    elastix_img_filter.PrintParameterMap()
-    transform_filter = elastix_img_filter.Execute()
-    transformed_img = elastix_img_filter.GetResultImage()
-    
-    # extract individual 3D images from 4D result image
-    extract_filter = sitk.ExtractImageFilter()
-    size = list(transformed_img.GetSize())
-    size[3] = 0 # set t to 0 to collapse this dimension
-    extract_filter.SetSize(size)
-    imgs = []
-    num_images = len(img_files)
-    for i in range(num_images):
-        extract_filter.SetIndex([0, 0, 0, i])  # x, y, z, t
-        img = extract_filter.Execute(transformed_img)
-        img_np = sitk.GetArrayFromImage(img)
-        # resize to original shape of first image, all aligned to position
-        # of subject within first image
-        img_large_np = np.zeros(size_orig[::-1])
-        img_large_np[:, start_y:start_y+img_np.shape[1]] = img_np
-        if show_imgs:
-            sitk.Show(sitk_io.replace_sitk_with_numpy(img, img_large_np))
-        imgs.append(img_large_np)
     
+    if sitk:
+        # groupwise register images in sitk
+        elastix_img_filter.SetParameterMap(param_map)
+        elastix_img_filter.PrintParameterMap()
+        elastix_img_filter.Execute()
+        transformed_img = elastix_img_filter.GetResultImage()
+        extract_filter = sitk.ExtractImageFilter()
+        
+        # extract individual 3D images from 4D result image
+        size = list(transformed_img.GetSize())
+        size[3] = 0  # set t to 0 to collapse this dimension
+        extract_filter.SetSize(size)
+        imgs = []
+        num_images = len(img_files)
+        for i in range(num_images):
+            extract_filter.SetIndex([0, 0, 0, i])  # x, y, z, t
+            img = extract_filter.Execute(transformed_img)
+            img_np = sitk_io.convert_img(img)
+            # resize to original shape of first image, all aligned to position
+            # of subject within first image
+            img_large_np = np.zeros(size_orig[::-1])
+            img_large_np[:, start_y:start_y + img_np.shape[1]] = img_np
+            if show_imgs and sitk:
+                sitk.Show(sitk_io.replace_sitk_with_numpy(img, img_large_np))
+            imgs.append(img_large_np)
+    
+    else:
+        # groupwise register images in ITK
+        reg_params.AddParameterMap(param_map)
+        elastix_img_filter.SetParameterObject(reg_params)
+        elastix_img_filter.SetLogToConsole(config.verbose)
+        elastix_img_filter.UpdateLargestPossibleRegion()
+        transformed_img = elastix_img_filter.GetOutput()
+        imgs = sitk_io.convert_img(transformed_img)
+
     # combine all images by taking their mean
     img_mean = np.mean(imgs, axis=0)
     extend_borders = settings["extend_borders"]
@@ -1172,8 +1394,7 @@ def register_group(
     img_raw = sitk_io.replace_sitk_with_numpy(transformed_img, img_mean)
     
     # carve groupwise registered image if given thresholds
-    imgs_to_show = []
-    imgs_to_show.append(img_raw)
+    imgs_to_show = [img_raw]
     holes_area = settings["holes_area"]
     if carve_threshold and holes_area:
         img_mean, _, img_mean_unfilled = cv_nd.carve(
@@ -1187,7 +1408,7 @@ def register_group(
         imgs_to_show.append(img_unfilled)
         imgs_to_show.append(transformed_img)
     
-    if show_imgs:
+    if show_imgs and sitk:
         for img in imgs_to_show: sitk.Show(img)
     
     #transformed_img = img_raw
@@ -1197,90 +1418,13 @@ def register_group(
         out_path = os.path.join(name_prefix, config.RegNames.IMG_GROUPED.value)
         if not os.path.exists(name_prefix):
             os.makedirs(name_prefix)
-        print("writing {}".format(out_path))
-        sitk.WriteImage(transformed_img, out_path, False)
-        img_np = sitk.GetArrayFromImage(transformed_img)
-        config.resolutions = [transformed_img.GetSpacing()[::-1]]
+        sitk_io.write_img(transformed_img, out_path)
+        img_np = sitk_io.convert_img(transformed_img)
+        config.resolutions = [tuple(transformed_img.GetSpacing())[::-1]]
         importer.save_np_image(img_np[None], out_path, config.series)
     
-    print("time elapsed for groupwise registration (s): {}"
-          .format(time() - start_time))
-
-
-def register_labels_to_atlas(path_fixed):
-    """Register annotation labels to its own atlas.
-    
-    Labels may not match their own underlying atlas image well, 
-    particularly in the orthogonal directions in which the labels 
-    were not constructed. To improve alignment between the labels 
-    and the atlas itself, register the edge-detected versions of the 
-    labels and its atlas.
-    
-    Edge files are assumed to have been generated by 
-    :func:``edge_seg.make_edge_images``.
-    
-    Args:
-        path_fixed: Path to the fixed file, typically the atlas file 
-            with stained sections. The corresponding edge and labels 
-            files will be loaded based on this path.
-    """
-    # load corresponding edge files
-    fixed_sitk = sitk_io.load_registered_img(
-        path_fixed, config.RegNames.IMG_ATLAS_EDGE.value, get_sitk=True)
-    moving_sitk = sitk_io.load_registered_img(
-        path_fixed, config.RegNames.IMG_LABELS_EDGE.value, get_sitk=True)
-    
-    # set up SimpleElastix filter
-    elastix_img_filter = sitk.ElastixImageFilter()
-    elastix_img_filter.SetFixedImage(fixed_sitk)
-    elastix_img_filter.SetMovingImage(moving_sitk)
-    
-    # set up registration parameters
-    settings = config.atlas_profile
-    param_map_vector = sitk.VectorOfParameterMap()
-    # bspline for non-rigid deformation
-    param_map = sitk.GetDefaultParameterMap("bspline")
-    param_map["FinalGridSpacingInVoxels"] = [
-        settings["bspline_grid_space_voxels"]]
-    del param_map["FinalGridSpacingInPhysicalUnits"] # avoid conflict with vox
-    param_map["MaximumNumberOfIterations"] = [settings["bspline_iter_max"]]
-    #param_map["MaximumNumberOfIterations"] = "0"
-    _config_reg_resolutions(
-        settings["grid_spacing_schedule"], param_map, fixed_sitk.GetDimension())
-    param_map_vector.append(param_map)
-    elastix_img_filter.SetParameterMap(param_map_vector)
-    elastix_img_filter.PrintParameterMap()
-    
-    # perform the registration
-    transform = elastix_img_filter.Execute()
-    transformed_img = elastix_img_filter.GetResultImage()
-    sitk.Show(transformed_img)
-    
-    # set up filter to apply the same transformation to label file, 
-    # ensuring that no new labels are interpolated
-    transform_param_map = elastix_img_filter.GetTransformParameterMap()
-    transform_param_map[-1]["FinalBSplineInterpolationOrder"] = ["0"]
-    transformix_img_filter = sitk.TransformixImageFilter()
-    transformix_img_filter.SetTransformParameterMap(transform_param_map)
-    
-    # apply transformation, manually resetting spacing in case of rounding
-    labels_sitk = sitk_io.load_registered_img(
-        path_fixed, config.RegNames.IMG_LABELS.value, get_sitk=True)
-    transformed_labels = register_repeat(
-        transformix_img_filter, labels_sitk, True)
-    transformed_labels.SetSpacing(transformed_img.GetSpacing())
-    #sitk.Show(transformed_labels)
-    
-    # write transformed labels file
-    out_path_base = os.path.splitext(path_fixed)[0] + "_edgereg.mhd"
-    imgs_write = {
-        config.RegNames.IMG_LABELS.value: transformed_labels}
-    sitk_io.write_reg_images(imgs_write, out_path_base)
-    # copy original atlas metadata file to allow opening this atlas 
-    # alongside new labels image for comparison
-    shutil.copy(
-        sitk_io.reg_out_path(path_fixed, config.RegNames.IMG_ATLAS.value), 
-        out_path_base)
+    _logger.info(
+        "Time elapsed for groupwise registration (s): %s", time() - start_time)
 
 
 def overlay_registered_imgs(fixed_file, moving_file_dir, plane=None, 
@@ -1307,18 +1451,17 @@ def overlay_registered_imgs(fixed_file, moving_file_dir, plane=None,
     image5d = img5d.img
     roi = image5d[0, ...]  # not using time dimension
     
-    # get the atlas file and transpose it to match the orientation of the 
+    # get the atlas file and transpose it to match the orientation of the
     # experiment image
     out_path = os.path.join(moving_file_dir, config.RegNames.IMG_ATLAS.value)
-    print("Reading in {}".format(out_path))
-    moving_sitk = sitk.ReadImage(out_path)
+    moving_sitk = sitk_io.read_img(out_path)
     moving_sitk = atlas_refiner.transpose_img(
         moving_sitk, plane, rotate)
-    moving_img = sitk.GetArrayFromImage(moving_sitk)
+    moving_img = sitk_io.convert_img(moving_sitk)
     
     # get the registered atlas file, which should already be transposed
     transformed_sitk = sitk_io.load_registered_img(name_prefix, get_sitk=True)
-    transformed_img = sitk.GetArrayFromImage(transformed_sitk)
+    transformed_img = sitk_io.convert_img(transformed_sitk)
     
     # get the registered labels file, which should also already be transposed
     labels_img = sitk_io.load_registered_img(
@@ -1590,7 +1733,7 @@ def volumes_by_id(
                 libmag.warn("will load atlas image instead")
                 img_sitk = sitk_io.load_registered_img(
                     mod_path, config.RegNames.IMG_ATLAS.value, get_sitk=True)
-            img_np = sitk.GetArrayFromImage(img_sitk)
+            img_np = sitk_io.convert_img(img_sitk)
             spacing = img_sitk.GetSpacing()[::-1]
             
             # load labels in order of priority: config > full labels
@@ -1781,7 +1924,7 @@ def volumes_by_id_compare(img_paths, labels_ref_paths, unit_factor=None,
                 img = sitk_io.load_registered_img(
                     img_path, config.RegNames.IMG_LABELS.value, get_sitk=True)
             labels_imgs_sitk.append(img)
-        labels_imgs = [sitk.GetArrayFromImage(img) for img in labels_imgs_sitk]
+        labels_imgs = [sitk_io.convert_img(img) for img in labels_imgs_sitk]
         spacing = labels_imgs_sitk[0].GetSpacing()[::-1]
         
         # load heat map of nuclei per voxel if available based on 1st path
@@ -1931,8 +2074,8 @@ def _test_region_from_id():
         # been given similarly to register fn
         path = os.path.join(
             config.filenames[1], config.RegNames.IMG_LABELS.value)
-        labels_img = sitk.ReadImage(path)
-        labels_img = sitk.GetArrayFromImage(labels_img)
+        labels_img = sitk_io.read_img(path)
+        labels_img = sitk_io.convert_img(labels_img)
         scaling = np.ones(3)
         print("loaded labels image from {}".format(path))
     else:
@@ -1940,8 +2083,8 @@ def _test_region_from_id():
         labels_img = sitk_io.load_registered_img(
             config.filename, config.RegNames.IMG_LABELS.value)
         if config.filename.endswith(".mhd"):
-            img = sitk.ReadImage(config.filename)
-            img = sitk.GetArrayFromImage(img)
+            img = sitk_io.read_img(config.filename)
+            img = sitk_io.convert_img(img)
             image5d = img[None]
         else:
             img5d = importer.read_file(config.filename, config.series)
@@ -1967,10 +2110,11 @@ def _test_curate_img(path, prefix):
     result_imgs = curate_img(
         fixed_img, labels_img, [atlas_img], inpaint=False, 
         holes_area=holes_area)
-    sitk.Show(fixed_img)
-    sitk.Show(labels_img)
-    sitk.Show(result_imgs[0])
-    sitk.Show(result_imgs[1])
+    if sitk:
+        sitk.Show(fixed_img)
+        sitk.Show(labels_img)
+        sitk.Show(result_imgs[0])
+        sitk.Show(result_imgs[1])
 
 
 def _test_smoothing_metric():
@@ -1983,16 +2127,18 @@ def _test_smoothing_metric():
 
 
 def main():
-    """Handle registration processing tasks as specified in 
-    :attr:`magmap.config.register_type`.
+    """Handle registration processing tasks
+    
+    Tasks are specified in :attr:`magmap.config.register_type`.
+    
     """
-    if hasattr(sitk.ProcessObject, "SetGlobalDefaultThreader"):
+    if sitk and hasattr(sitk.ProcessObject, "SetGlobalDefaultThreader"):
         # manually set threader for SimpleITK >= 2 to avoid potential hangs
         # during Python multiprocessing; set by default in SimpleITK 2.1
         sitk.ProcessObject.SetGlobalDefaultThreader("Platform")
     
-    # name prefix to use a different name from the input files, such as when 
-    # registering transposed/scaled images but outputting paths corresponding 
+    # name prefix to use a different name from the input files, such as when
+    # registering transposed/scaled images but outputting paths corresponding
     # to the original image
     if config.prefix is not None:
         print("Formatting registered filenames to match {}"
@@ -2209,10 +2355,6 @@ def main():
             edge_seg.make_edge_images(
                 img_path, config.show, atlas, config.suffix, config.load_labels)
     
-    elif reg is config.RegisterTypes.REG_LABELS_TO_ATLAS:
-        # register labels to its underlying atlas
-        register_labels_to_atlas(config.filename)
-    
     elif reg in (
             config.RegisterTypes.MERGE_ATLAS_SEGS,
             config.RegisterTypes.MERGE_ATLAS_SEGS_EXP):
diff --git a/magmap/io/export_regions.py b/magmap/io/export_regions.py
index 584d3f1a..0bedf41a 100644
--- a/magmap/io/export_regions.py
+++ b/magmap/io/export_regions.py
@@ -1,5 +1,5 @@
 # Region labels export to data frames and CSV files
-# Author: David Young, 2019
+# Author: David Young, 2019, 2023
 """Region labels export to data frames and CSV files.
 
 Convert regions from ontology files or atlases to data frames.
@@ -11,7 +11,10 @@
 from time import time
 from typing import Dict, Optional, Sequence, Tuple
 
-import SimpleITK as sitk
+try:
+    import SimpleITK as sitk
+except ImportError:
+    sitk = None
 import numpy as np
 import pandas as pd
 
@@ -186,7 +189,7 @@ def export_common_labels(img_paths, output_path):
 def make_density_image(
         img_path: str, scale: Optional[float] = None,
         shape: Optional[Sequence[int]] = None, suffix: Optional[str] = None, 
-        labels_img_sitk: Optional[sitk.Image] = None,
+        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,
@@ -257,7 +260,7 @@ def make_heat_map():
         scaling = np.divide(shape, blobs.roi_size)
         labels_spacing = np.divide(blobs.resolutions[0], scaling)
         labels_img = np.zeros(shape, dtype=np.uint8)
-        labels_img_sitk = sitk.GetImageFromArray(labels_img)
+        labels_img_sitk = sitk_io.convert_img(labels_img)
         labels_img_sitk.SetSpacing(labels_spacing[::-1])
     
     else:
@@ -265,7 +268,7 @@ def make_heat_map():
         if labels_img_sitk is None:
             labels_img_sitk = sitk_io.load_registered_img(
                 mod_path, config.RegNames.IMG_LABELS.value, get_sitk=True)
-        labels_img = sitk.GetArrayFromImage(labels_img_sitk)
+        labels_img = sitk_io.convert_img(labels_img_sitk)
         
         is_2d = labels_img.ndim == 2
         labels_res = list(labels_img_sitk.GetSpacing()[::-1])
@@ -447,7 +450,7 @@ def make_labels_diff_img(img_path, df_path, meas, fn_avg, prefix=None,
     reg_name = (config.RegNames.IMG_LABELS.value if level is None 
                 else config.RegNames.IMG_LABELS_LEVEL.value.format(level))
     labels_sitk = sitk_io.load_registered_img(img_path, reg_name, get_sitk=True)
-    labels_np = sitk.GetArrayFromImage(labels_sitk)
+    labels_np = sitk_io.convert_img(labels_sitk)
     df = pd.read_csv(df_path)
     labels_diff = vols.map_meas_to_labels(
         labels_np, df, meas, fn_avg, col_wt=col_wt)
@@ -466,7 +469,7 @@ def make_labels_diff_img(img_path, df_path, meas, fn_avg, prefix=None,
     imgs_write = {reg_diff: labels_diff_sitk}
     out_path = prefix if prefix else img_path
     sitk_io.write_reg_images(imgs_write, out_path)
-    if show:
+    if show and sitk:
         for img in imgs_write.values():
             if img: sitk.Show(img)
 
@@ -474,7 +477,7 @@ def make_labels_diff_img(img_path, df_path, meas, fn_avg, prefix=None,
 def make_labels_level_img(
         img_path: Optional[str], level: int, prefix: Optional[str] = None,
         show: bool = False
-) -> Dict[str, sitk.Image]:
+) -> Dict[str, "sitk.Image"]:
     """Replace labels in an image with their parents at the given level.
     
     Labels that do not fall within a parent at that level will remain in place.
@@ -507,7 +510,7 @@ def make_labels_level_img(
         labels_sitk = sitk_io.load_registered_img(
             img_path, config.RegNames.IMG_LABELS.value, get_sitk=True)
         ref = ontology.LabelsRef(config.load_labels).load()
-    labels_np = sitk.GetArrayFromImage(labels_sitk)
+    labels_np = sitk_io.convert_img(labels_sitk)
     
     # remap labels to given level
     labels_np = ontology.make_labels_level(labels_np, ref, level)
@@ -525,7 +528,7 @@ def make_labels_level_img(
     }
     out_path = prefix if prefix else img_path
     sitk_io.write_reg_images(imgs_write, out_path)
-    if show:
+    if show and sitk:
         for img in imgs_write.values():
             if img: sitk.Show(img)
     
diff --git a/magmap/io/export_rois.py b/magmap/io/export_rois.py
index 82b85701..caf6fd51 100644
--- a/magmap/io/export_rois.py
+++ b/magmap/io/export_rois.py
@@ -1,5 +1,5 @@
 # ROI exporter for MagellanMapper
-# Author: David Young, 2017, 2019
+# Author: David Young, 2017, 2023
 """ROI exporter for MagellanMapper.
 
 Convert images and corresponding database entries into formats for 
@@ -13,16 +13,12 @@
 import numpy as np
 import pandas as pd
 from matplotlib import pyplot as plt
-import SimpleITK as sitk
 
-from magmap.settings import config
-from magmap.cv import cv_nd
-from magmap.cv import detector
-from magmap.io import libmag
-from magmap.io import sqlite
-from magmap.plot import plot_3d
+from magmap.cv import cv_nd, detector
 from magmap.gui import roi_editor
-from magmap.io import df_io
+from magmap.io import df_io, libmag, sitk_io, sqlite
+from magmap.plot import plot_3d
+from magmap.settings import config
 from magmap.stats import vols
 
 if TYPE_CHECKING:
@@ -152,15 +148,15 @@ def export_rois(
             # WORKAROUND: for some reason SimpleITK gives a conversion error 
             # when converting from uint16 (>u2) Numpy array
             img3d = img3d.astype(np.float64)
-            img3d_sitk = sitk.GetImageFromArray(img3d)
+            img3d_sitk = sitk_io.convert_img(img3d)
             '''
             print(img3d_sitk)
             print("orig img:\n{}".format(img3d[0]))
-            img3d_back = sitk.GetArrayFromImage(img3d_sitk)
+            img3d_back = sitk_io.convert_img(img3d_sitk)
             print(img3d.shape, img3d.dtype, img3d_back.shape, img3d_back.dtype)
             print("sitk img:\n{}".format(img3d_back[0]))
             '''
-            sitk.WriteImage(img3d_sitk, path_img_nifti, False)
+            sitk_io.write_img(img3d_sitk, path_img_nifti)
             roi_ed = roi_editor.ROIEditor(img5d)
             print("shape", img3d.shape)
             roi_ed.plot_roi(
@@ -197,9 +193,8 @@ def export_rois(
                 img3d_truth[img3d_truth < 0.5] = 0
             _logger.debug("Exporting truth ROI of shape: %s", img3d_truth.shape)
             np.save(path_img_annot, img3d_truth)
-            sitk.WriteImage(
-                sitk.GetImageFromArray(img3d_truth), path_img_annot_nifti, 
-                False)
+            sitk_io.write_img(
+                sitk_io.convert_img(img3d_truth), path_img_annot_nifti)
             _logger.info(
                 "Wrote NIfTI formatted images: %s", path_img_annot_nifti)
             # avoid smoothing interpolation, using "nearest" instead
diff --git a/magmap/io/libmag.py b/magmap/io/libmag.py
index db51ac2f..c4c71c1a 100644
--- a/magmap/io/libmag.py
+++ b/magmap/io/libmag.py
@@ -117,30 +117,37 @@ def roll_elements(arr, shift, axis=None):
     return arr
 
 
-def pad_seq(seq, n, pad=None):
-    """Pad a sequence with a given value or truncate the sequence to fit 
-    a given length.
+def pad_seq(
+        seq: Sequence[Any], n: int, pad: Optional[Any] = None) -> Sequence[Any]:
+    """Pad a sequence with a given value or truncate to fit a given length.
     
     Args:
-        seq (List): Sequence to fill in-place.
-        n (int): Target length.
-        pad (float, List): Value with which to fill; defaults to None.
-            If a sequence, filling with start with the first corresponding
-            missing value in ``seq``.
+        seq: Sequence to fill in-place.
+        n: Target length.
+        pad: Value with which to fill; defaults to None.
+            If a sequence, missing values in ``seq`` will be filled with
+            corresponding values in ``pad``.
     
     Returns:
-        :List: A truncated view of ``seq`` if the sequence is longer than ``n``
-        or ``seq`` as a List with ``pad`` appended to reach a length of ``n``.
+        A truncated view of ``seq`` if the sequence is longer than ``n``
+        or ``seq`` with ``pad`` appended to reach a length of ``n``. ``seq``
+        is modified in-place if it is extended and not a NumPy array.
+    
     """
     len_seq = len(seq)
     if len_seq >= n:
         # truncate if seq is longer than n is
         seq = seq[:n]
+    
     else:
-        # pad with the given value
         if isinstance(seq, np.ndarray):
             # convert to list if ndarray to allow mixing with None values
             seq = seq.tolist()
+        elif isinstance(seq, tuple):
+            # convert to list to allow adding values
+            seq = list(seq)
+
+        # pad with the given value
         if is_seq(pad):
             if len(pad) > len_seq:
                 # fill starting with corresponding first missing value
@@ -151,9 +158,41 @@ def pad_seq(seq, n, pad=None):
         else:
             # fill with pad value repeats
             seq += [pad] * (n - len_seq)
+    
     return seq
 
 
+def replace_seq(
+        seq: Sequence[Any], replacement: Sequence[Any]) -> np.ndarray:
+    """Replace elements in one sequence with those from another.
+    
+    Sequences will be converted to NumPy arrays to support multiple dimensions.
+    
+    Args:
+        seq: Sequence to be replaced.
+        replacement: Sequence to replace in ``seq``. Must be of the same
+            shape and type as those of ``seq``.
+
+    Returns:
+        ``seq`` as a new NumPy array with as many values ass possible replaced
+        by ``replacement``.
+
+    """
+    # convert to np to support multi-dimensional arrays
+    seq_np = np.array(seq)
+    tgt_np = np.array(replacement)
+    
+    # truncate to smallest shape in all dimension
+    shape = np.minimum(seq_np.shape, tgt_np.shape)
+    slices = []
+    for s in shape:
+        slices.append(slice(s))
+    slices = tuple(slices)
+    seq_np[slices] = tgt_np[slices]
+    
+    return seq_np
+
+
 def combine_arrs(arrs, filter_none=True, fn=None, **kwargs):
     """Combine arrays with array filtering.
     
diff --git a/magmap/io/np_io.py b/magmap/io/np_io.py
index ad0ed516..6c95b3e2 100644
--- a/magmap/io/np_io.py
+++ b/magmap/io/np_io.py
@@ -458,7 +458,7 @@ def add_metadata():
         if bg_atlas:
             # extract labels image from BrainGlobe atlas
             config.labels_img = bg_atlas.annotation
-            config.labels_img_sitk = sitk_io.numpy_to_sitk(config.labels_img)
+            config.labels_img_sitk = sitk_io.convert_img(config.labels_img)
         else:
             try:
                 # load labels image
diff --git a/magmap/io/sitk_io.py b/magmap/io/sitk_io.py
index 9cee2c86..4706bae3 100644
--- a/magmap/io/sitk_io.py
+++ b/magmap/io/sitk_io.py
@@ -6,10 +6,17 @@
 """
 import os
 import shutil
-from typing import Any, Dict, List, Optional, Sequence, Tuple, Union
+from typing import Callable, Dict, List, Optional, Sequence, Tuple, Union
 
+try:
+    import itk
+except ImportError:
+    itk = None
 import numpy as np
-import SimpleITK as sitk
+try:
+    import SimpleITK as sitk
+except ImportError:
+    sitk = None
 from skimage import transform
 
 from magmap.settings import config
@@ -20,23 +27,6 @@
 EXTS_3D: Sequence[str] = (".mhd", ".mha", ".nii.gz", ".nii", ".nhdr", ".nrrd")
 # TODO: include all formats supported by SimpleITK
 
-# Mapping of NumPy dtypes to SimpleITK pixel IDs.
-_DTYPES_NP_SITK: Dict[str, Any] = {
-    np.dtype(np.character).name: sitk.sitkUInt8,
-    np.dtype(np.uint8).name: sitk.sitkUInt8,
-    np.dtype(np.uint16).name: sitk.sitkUInt16,
-    np.dtype(np.uint32).name: sitk.sitkUInt32,
-    np.dtype(np.uint64).name: sitk.sitkUInt64,
-    np.dtype(np.int8).name: sitk.sitkInt8,
-    np.dtype(np.int16).name: sitk.sitkInt16,
-    np.dtype(np.int32).name: sitk.sitkInt32,
-    np.dtype(np.int64).name: sitk.sitkInt64,
-    np.dtype(np.float32).name: sitk.sitkFloat32,
-    np.dtype(np.float64).name: sitk.sitkFloat64,
-    np.dtype(np.complex64).name: sitk.sitkComplexFloat32,
-    np.dtype(np.complex128).name: sitk.sitkComplexFloat64,
-}
-
 _logger = config.logger.getChild(__name__)
 
 
@@ -64,9 +54,54 @@ def reg_out_path(file_path, reg_name, match_ext=False):
         return file_path_base + "_" + reg_name
 
 
+def convert_img(
+        img: Union["sitk.Image", "itk.Image", np.ndarray],
+        multichannel: Optional[bool] = None, to_sitk: Optional[bool] = None
+) -> Union["sitk.Image", "itk.Image", np.ndarray]:
+    """Convert SimpleITK, ITK, and NumPy images.
+    
+    Args:
+        img: SimpleITK Image, ITK Image, or NumPy array. Non-NumPy arrays
+            will be converted to a NumPy array.
+        multichannel: True if the resulting SimpleITK or ITK Image should be
+            multichannel. Defaults to None, which will attempt to auto-detect
+            multichannel images. Only used if ``img`` is a NumPy array.
+        to_sitk: True to convert a NumPy ``img`` to a SimpleITK Image. If
+            False, an ITK Image is output instead. If None, the parameter will
+            be True if the SimpleITK library is found. Ignored if ``img`` is
+            not a NumPy array.
+
+    Returns:
+        The converted image.
+
+    """
+    if to_sitk is None:
+        # default to convert to SimpleITK if the library is present
+        to_sitk = sitk is not None
+    
+    conv = img
+    if sitk and isinstance(img, sitk.Image):
+        # convert an sitk Image to an np array
+        conv = sitk.GetArrayFromImage(img)
+    elif itk and isinstance(img, (itk.Image, itk.VectorImage)):
+        # convert an ITK image to an np array
+        conv = itk.GetArrayFromImage(img)
+    elif sitk and to_sitk:
+        # convert an np array to sitk Image
+        conv = sitk.GetImageFromArray(img, multichannel)
+    elif itk:
+        # convert an np array to ITK Image
+        if multichannel is None:
+            # default to treat 4D+ images as multichannel
+            multichannel = img.ndim > 3
+        conv = itk.GetImageFromArray(img, multichannel)
+    return conv
+
+
 def replace_sitk_with_numpy(
-        img_sitk: sitk.Image, img_np: np.ndarray, multichannel: bool = None
-) -> sitk.Image:
+        img_sitk: Union["sitk.Image", "itk.Image"], img_np: np.ndarray,
+        multichannel: bool = None
+) -> Union["sitk.Image", "itk.Image"]:
     """Generate a :class:``sitk.Image`` object based on another Image object,
     but replace its array with a new Numpy array.
     
@@ -77,54 +112,40 @@ def replace_sitk_with_numpy(
             the multichannel setting from ``img_sitk`` will be used.
     
     Returns:
-        New :class:``sitk.Image`` object with same spacing, origin, and 
-        direction as that of ``img_sitk`` and array replaced with ``img_np``.
+        New SimpleITK or ITK ``sitk.Image`` object with same spacing, origin,
+        and direction as that of ``img_sitk`` and array replaced by ``img_np``.
     
     """
-    # get original settings
-    spacing = img_sitk.GetSpacing()
-    origin = img_sitk.GetOrigin()
-    direction = img_sitk.GetDirection()
-    
     # treat as vector (multichannel) image if source is multichannel
     if multichannel is None:
         multichannel = (
             True if img_sitk.GetNumberOfComponentsPerPixel() > 1 else None)
-    img_sitk_back = sitk.GetImageFromArray(img_np, multichannel)
-    
-    # transfer original settings to new sitk Image
-    img_sitk_back.SetSpacing(spacing)
-    img_sitk_back.SetOrigin(origin)
+    img_sitk_back = convert_img(
+        img_np, multichannel, sitk and isinstance(img_sitk, sitk.Image))
     
-    try:
-        img_sitk_back.SetDirection(direction)
-    except RuntimeError:
-        # direction format and length may not be directly transferable, such
-        # as direction from groupwise reg image
-        _logger.warn(
-            "Could not replace image direction with: %s\n"
-            "Leaving default direction: %s",
-            direction, img_sitk_back.GetDirection())
+    # transfer original settings to new Image
+    match_world_info(img_sitk, img_sitk_back)
     
     return img_sitk_back
 
 
 def match_world_info(
-        source: sitk.Image, target: sitk.Image,
+        source: Union["sitk.Image", "itk.Image"],
+        target: Union["sitk.Image", "itk.Image"],
         spacing: Union[bool, Tuple[int], List[int]] = True,
         origin: Union[bool, Tuple[int], List[int]] = True,
         direction: Union[bool, Tuple[int], List[int]] = True):
-    """Copy world information (eg spacing, origin, direction) from one
-    image object to another.
+    """Transfer world information from one image object to another.
 
-    This matching is sometimes necessary for slight differences in
+    Supported world information is spacing, origin, direction. Matching
+    this information is sometimes necessary for slight differences in
     metadata perhaps from founding that may prevent ITK filters from executing.
 
     Args:
-        source (:obj:`sitk.Image`): Source object whose relevant metadata
-            will be copied into ``target``.
-        target (:obj:`sitk.Image`): Target object whose corresponding
-            metadata will be overwritten by that of ``source``.
+        source: Source object whose relevant metadata will be copied into
+            ``target``.
+        target: Target object whose corresponding metadata will be
+            overwritten by that of ``source``.
         spacing: True to copy the spacing from ``source`` to ``target``, or
             the spacing to set in ``target``; defaults to True.
         origin: True to copy the origin from ``source`` to ``target``, or
@@ -141,34 +162,86 @@ def match_world_info(
     if direction is True:
         direction = source.GetDirection()
     
-    # set the values in the target
-    _logger.debug(
-        "Adjusting spacing from %s to %s, origin from %s to %s, "
-        "direction from %s to %s", target.GetSpacing(), spacing,
-        target.GetOrigin(), origin, target.GetDirection(), direction)
     if spacing:
+        # transfer spacing, truncating if necessary
+        spacing = libmag.replace_seq(target.GetSpacing(), spacing)
         target.SetSpacing(spacing)
+    
     if origin:
+        # transfer origin, truncating if necessary
+        origin = libmag.replace_seq(target.GetOrigin(), origin)
         target.SetOrigin(origin)
+    
     if direction:
-        target.SetDirection(direction)
+        # convert directions to 2D arrays if necessary
+        dir_np = np.array(direction)
+        dir_back = np.array(target.GetDirection())
+        is_dir_1d = dir_np.ndim == 1
+        if is_dir_1d:
+            # sitk uses a 1D array; determine number of dimensions from spacing
+            dir_np = np.reshape(np.array(dir_np), [len(source.GetSpacing())] * 2)
+            dir_back = np.reshape(
+                np.array(dir_back), [len(target.GetSpacing())] * 2)
+        
+        # fill target direction with source, truncating it if necessary for ITK
+        dir_back = libmag.replace_seq(dir_back, dir_np)
+        if is_dir_1d:
+            dir_back = np.ravel(dir_back)
+        
+        try:
+            target.SetDirection(dir_back)
+        except RuntimeError:
+            # direction format and length may not be directly transferable, such
+            # as direction from groupwise reg image
+            _logger.warn(
+                "Could not replace image direction with: %s\n"
+                "Leaving default direction: %s",
+                direction, target.GetDirection())
 
 
-def read_sitk(path, dryrun=False):
-    """Read an image file into :class:``sitk.Image`` format, checking for 
-    alternative supported extensions if necessary.
+def read_img(path) -> Union["sitk.Image", "itk.Image"]:
+    """Read an image file into SimpleITK or ITK ``Image`` format.
     
     Args:
-        path (str): Path, including prioritized extension to check first.
-        dryrun (bool): True to find an existing path if available, without
+        path: Path, including prioritized extension to check first.
+    
+    Returns:
+        Image object located at ``path``. Loads into SimpleITK format if the
+        library is available, otherwise into ITK format.
+    
+    """
+    if sitk:
+        # load as sitk Image if available
+        _logger.debug("Loading image with SimpleITK: %s", path)
+        img_sitk = sitk.ReadImage(path)
+    else:
+        # load as ITK image
+        _logger.debug("Loading image with ITK: %s", path)
+        img_sitk = itk.imread(path)
+    return img_sitk
+
+
+def read_sitk(
+        path: str, dryrun: bool = False
+) -> Tuple[Optional[Union["sitk.Image", "itk.Image"]], Optional[str]]:
+    """Read an image file into SimpleITK or ITK ``Image`` format
+    
+    Checks for alternative supported extensions if necessary.
+    
+    Args:
+        path: Path, including prioritized extension to check first.
+        dryrun: True to find an existing path if available, without
             loading the image; defaults to False.
     
     Returns:
-        :obj:`sitk.Image`, str: Image object located at ``path`` with
-        the found extension, or None if unable to load; and the loaded path,
-        or None if no matching, existing path is found. If a file at
-        ``path`` cannot be found, its extension is replaced successively
-        with remaining extensions in :const:``EXTS_3D`` until a file is found.
+        Tuple of:
+        - ``img_sitk``: Image object located at ``path`` with the found
+          extension, or None if unable to load. Type is determined by
+          :meth:`read_img`
+        - ``path_loaded``: the loaded path, or None if no matching, existing
+          path is found. If a file at ``path`` cannot be found, its extension
+          is replaced successively with remaining extensions in
+          :const:``EXTS_3D`` until a file is found.
     
     """
     # prioritize given extension
@@ -184,14 +257,15 @@ def read_sitk(path, dryrun=False):
         img_path = path_split[0] + ext
         if os.path.exists(img_path):
             if not dryrun:
-                _logger.debug("Loading image with SimpleITK: %s", img_path)
-                img_sitk = sitk.ReadImage(img_path)
+                img_sitk = read_img(img_path)
             path_loaded = img_path
             break
+    
     if not dryrun and img_sitk is None:
         _logger.warn(
             "could not find image from %s and extensions %s", path_split[0],
             exts)
+    
     return img_sitk, path_loaded
 
 
@@ -204,7 +278,7 @@ def _load_reg_img_to_combine(path, reg_name, img_nps):
         img_np_base = img_nps[0]
     img_sitk, loaded_path = load_registered_img(
         path, reg_name, get_sitk=True, return_path=True)
-    img_np = sitk.GetArrayFromImage(img_sitk)
+    img_np = convert_img(img_sitk)
     if img_np_base is not None:
         if img_np_base.shape != img_np.shape:
             # resize to first image
@@ -218,7 +292,7 @@ def _load_reg_img_to_combine(path, reg_name, img_nps):
     return img_sitk, loaded_path
 
 
-def _make_3d(img_sitk: sitk.Image, spacing_z=1) -> sitk.Image:
+def _make_3d(img_sitk: "sitk.Image", spacing_z=1) -> "sitk.Image":
     """Make a 2D image into a 3D image with a single plane.
     
     Args:
@@ -232,9 +306,9 @@ def _make_3d(img_sitk: sitk.Image, spacing_z=1) -> sitk.Image:
     spacing = img_sitk.GetSpacing()
     if len(spacing) == 2:
         # prepend an additional axis for 2D images to make them 3D
-        img_np = sitk.GetArrayFromImage(img_sitk)[None]
+        img_np = convert_img(img_sitk)[None]
         spacing = list(spacing) + [spacing_z]
-        img_sitk = sitk.GetImageFromArray(img_np)
+        img_sitk = convert_img(img_np)
         img_sitk.SetSpacing(spacing)
     return img_sitk
 
@@ -243,7 +317,7 @@ def read_sitk_files(
         filename_sitk: str,
         reg_names: Optional[Union[str, Sequence[str]]] = None,
         return_sitk: bool = False, make_3d: bool = False
-) -> Union["np_io.Image5d", Tuple["np_io.Image5d", sitk.Image]]:
+) -> Union["np_io.Image5d", Tuple["np_io.Image5d", "sitk.Image"]]:
     """Read image files through SimpleITK.
     
     Supports identifying files based on registered suffixes and combining
@@ -308,14 +382,14 @@ def read_sitk_files(
             raise FileNotFoundError(
                 "could not find file {}".format(filename_sitk))
         img_sitk, _ = read_sitk(filename_sitk)
-        img_np = sitk.GetArrayFromImage(img_sitk)
+        img_np = convert_img(img_sitk)
     
     if make_3d:
         # convert 2D images to 3D
         # TODO: consider converting img_np to 3D regardless so array in
         #   Image5d is at least 4D
         img_sitk = _make_3d(img_sitk)
-        img_np = sitk.GetArrayFromImage(img_sitk)
+        img_np = convert_img(img_sitk)
     
     if config.resolutions is None:
         # fallback to determining metadata directly from sitk file
@@ -335,8 +409,8 @@ def read_sitk_files(
 
 def load_registered_imgs(
         img_path: str, reg_names: Sequence[str], *args, **kwargs
-) -> Dict[str, Union[Union[np.ndarray, sitk.Image],
-          Tuple[Union[np.ndarray, sitk.Image], str]]]:
+) -> Dict[str, Union[Union[np.ndarray, "sitk.Image"],
+          Tuple[Union[np.ndarray, "sitk.Image"], str]]]:
     """Load atlas-based images registered to another image.
     
     Args:
@@ -364,8 +438,8 @@ def load_registered_imgs(
 def load_registered_img(
         img_path: str, reg_name: str, get_sitk: bool = False,
         return_path: bool = False
-) -> Union[Union[np.ndarray, sitk.Image],
-           Tuple[Union[np.ndarray, sitk.Image], str]]:
+) -> Union[Union[np.ndarray, "sitk.Image"],
+           Tuple[Union[np.ndarray, "sitk.Image"], str]]:
     """Load atlas-based image that has been registered to another image.
     
     Args:
@@ -409,7 +483,7 @@ def load_registered_img(
     
     if not get_sitk:
         # convert to ndarray
-        reg_img = sitk.GetArrayFromImage(reg_img)
+        reg_img = convert_img(reg_img)
     
     # return image, including path if flagged
     return (reg_img, reg_img_path) if return_path else reg_img
@@ -440,39 +514,63 @@ def find_atlas_labels(
             os.path.dirname(labels_ref_path),
             config.RegNames.IMG_LABELS.value))
         if orig_labels_sitk is not None:
-            config.labels_img_orig = sitk.GetArrayFromImage(orig_labels_sitk)
+            config.labels_img_orig = convert_img(orig_labels_sitk)
             return np.unique(config.labels_img_orig).tolist()
     # fall back to using all labels in ontology reference, including any
     # hierarchical labels
     return list(labels_ref_lookup.keys())
 
 
-def write_registered_image(img_np, img_path, reg_name, img_sitk=None,
-                           load_reg_names=None, overwrite=False):
+def write_img(
+        img: Union["sitk.Image", "itk.Image"], path: str,
+        compression: bool = False, **kwargs):
+    """Write SimpleITK/ITK images.
+    
+    Args:
+        img: Image as SimpleITK or ITK type.
+        path: Path to save the image.
+        compression: True to use compression.
+        kwargs: Additional arguments to underlying :meth:`sitk.WriteImage`
+            or :meth:`itk.imwrite`.
+
+    """
+    # write using the available ITK-based library
+    if sitk and isinstance(img, sitk.Image):
+        sitk.WriteImage(img, path, compression, **kwargs)
+    else:
+        itk.imwrite(img, path, compression, **kwargs)
+    _logger.info("Wrote image to: %s", path)
+
+
+def write_registered_image(
+        img_np: np.ndarray, img_path: str, reg_name: str,
+        img_sitk: Optional["sitk.Image"] = None,
+        load_reg_names: Optional[Sequence[str]] = None,
+        overwrite: bool = False) -> "sitk.Image":
     """Write a Numpy array as a registered 3D image file through SimpleITK.
     
     To find metadata for the output image, another SimpleITK image must
     be given or discovered as a registered image.
     
     Args:
-        img_np (:class:`numpy.ndarray`): Image array to write.
-        img_path (str): Base path from which to construct output path.
-        reg_name (str): Registered image suffix, which will also specify the
+        img_np: Image array to write.
+        img_path: Base path from which to construct output path.
+        reg_name: Registered image suffix, which will also specify the
             output file format.
-        img_sitk (:class:`sitk.Image`): SimpleITK Image object to use as
+        img_sitk: SimpleITK Image object to use as
             a template for image metadata; defaults to None, in which case
             a registered image will be loaded through ``load_reg_names``.
-        load_reg_names (List[str]): Sequence of registered image suffixes
+        load_reg_names: Sequence of registered image suffixes
             from which to load a template image for metdata. Names are
             checked until the first image loads successfully. Defaults to
             None to use the atlas volume followed by the experimental image.
             If ``img_sitk`` is None and no registered image can be found,
             writing will be aborted.
-        overwrite (bool): True to overwrite any existing image at the output
+        overwrite: True to overwrite any existing image at the output
             path; defaults to False.
 
     Returns:
-        :class:`sitk.Image`: The saved image as a SimpleITK Image object.'
+        The saved image as a SimpleITK Image object.
     
     Raises:
         FileExistsError: if ``overwrite`` is false and existing file is at
@@ -505,7 +603,7 @@ def write_registered_image(img_np, img_path, reg_name, img_sitk=None,
     if img_sitk:
         # copy the template and replace its array with the given array to save
         reg_img = replace_sitk_with_numpy(img_sitk, img_np)
-        sitk.WriteImage(reg_img, reg_img_path, False)
+        write_img(reg_img, reg_img_path)
         print("wrote {} with current registered image".format(reg_img_path))
         return reg_img
     else:
@@ -514,24 +612,27 @@ def write_registered_image(img_np, img_path, reg_name, img_sitk=None,
             f"'{reg_img_path}'")
 
 
-def write_reg_images(imgs_write, prefix, copy_to_suffix=False, ext=None,
-                     prefix_is_dir=False):
+def write_reg_images(
+        imgs_write: Dict[str, Union["sitk.Image", "itk.Image"]], prefix: str,
+        copy_to_suffix: bool = False, ext: Optional[str] = None,
+        prefix_is_dir: bool = False):
     """Write registered images to file.
     
     Args:
-        imgs_write: Dictionary of ``{suffix: image}``, where ``suffix`` 
-            is a registered images suffix, such as :const:``IMAGE_LABELS``, 
-            and ``image`` is a SimpleITK image object. If the image does 
+        imgs_write: Dictionary of ``{suffix: image}``, where ``suffix``
+            is a registered images suffix, such as :const:``IMAGE_LABELS``,
+            and ``image`` is a SimpleITK or ITK image object. If the image does
             not exist, the file will not be written.
         prefix: Base path from which to construct registered file paths.
             Parent directories will be created if necessary.
-        copy_to_suffix: If True, copy the output path to a file in the 
-            same directory with ``suffix`` as the filename, which may 
-            be useful when setting the registered images as the 
-            main images in the directory. Defaults to False.
+        copy_to_suffix: If True, copy the output path to a file in the
+            same directory with ``suffix`` as the filename, which may
+            be useful when setting the registered images as the
+            main images in the directory.
         ext: Replace extension with this value if given; defaults to None.
-        prefix_is_dir (bool): True to treat ``prefix`` as a directory;
-            defaults to False to get the directory name of ``prefix``.
+        prefix_is_dir: True to treat ``prefix`` as a directory.
+            If False, will use the directory name of ``prefix``.
+    
     """
     # get parent directories and create them if necessary
     target_dir = prefix if prefix_is_dir else os.path.dirname(prefix)
@@ -544,14 +645,16 @@ def write_reg_images(imgs_write, prefix, copy_to_suffix=False, ext=None,
         if img is None: continue
         if ext: suffix = libmag.match_ext(ext, suffix)
         out_path = reg_out_path(prefix, suffix)
-        sitk.WriteImage(img, out_path, False)
-        print("wrote registered image to", out_path)
+        
+        # write using the available ITK-based library
+        write_img(img, out_path)
+        
         if copy_to_suffix:
-            # copy metadata file to allow opening images from bare suffix name, 
+            # copy metadata file to allow opening images from bare suffix name,
             # such as when this atlas becomes the new atlas for registration
             out_path_copy = os.path.join(target_dir, suffix)
             shutil.copy(out_path, out_path_copy)
-            print("also copied to", out_path_copy)
+            _logger.info("Also copied to: %s", out_path_copy)
 
 
 def write_pts(path: str, pts: Sequence[Sequence[int]], pt_type: str = "index"):
@@ -575,8 +678,11 @@ def write_pts(path: str, pts: Sequence[Sequence[int]], pt_type: str = "index"):
             f.write(f"{' '.join([str(p) for p in pt])}\n")
 
 
-def merge_images(img_paths, reg_name, prefix=None, suffix=None, 
-                 fn_combine=np.sum):
+def merge_images(
+        img_paths: Sequence[str], reg_name: str, prefix: Optional[str] = None,
+        suffix: Optional[str] = None,
+        fn_combine: Callable[[Sequence[np.ndarray]], np.ndarray] = np.sum
+) -> "sitk.Image":
     """Merge images from multiple paths.
     
     Assumes that the images are relatively similar in size, but will resize 
@@ -584,18 +690,18 @@ def merge_images(img_paths, reg_name, prefix=None, suffix=None,
     
     Args:
         img_paths: Paths from which registered paths will be found.
-        reg_name: Registration suffix to load for the given paths 
+        reg_name: Registration suffix to load for the given paths
             in ``img_paths``.
-        prefix: Start of output path; defaults to None to use the first 
-           path in ``img_paths`` instead.
-        suffix: Portion of path to be combined with each path 
-            in ``img_paths`` and output path; defaults to None.
-        fn_combine: Function to apply to combine images with ``axis=0``. 
-            Defaults to :func:``np.sum``. If None, each image will be 
-            inserted as a separate channel.
+        prefix: Start of output path. If None, the first path in ``img_paths``
+            will be used instead.
+        suffix: Portion of path to be combined with each path
+            in ``img_paths`` and output path.
+        fn_combine: Function to apply to combine images with ``axis=0``.
+            If None, each image will be inserted as a separate channel.
     
     Returns:
         The combined image in SimpleITK format.
+    
     """
     if len(img_paths) < 1: return None
     
@@ -631,7 +737,7 @@ def merge_images(img_paths, reg_name, prefix=None, suffix=None,
 
 def load_numpy_to_sitk(
         img5d: Union[str, "np_io.Image5d"], rotate: bool = False,
-        channel: Optional[Union[int, Sequence[int]]] = None) -> sitk.Image:
+        channel: Optional[Union[int, Sequence[int]]] = None) -> "sitk.Image":
     """Load Numpy image array to SimpleITK Image object.
 
     Use ``channel`` to extract a single channel before generating a
@@ -663,11 +769,11 @@ def load_numpy_to_sitk(
     
     if rotate:
         # rotate 180 deg
-        # TODO: consider deprecating, deferring to atlas_refier.transpose_img
+        # TODO: consider deprecating, deferring to atlas_refiner.transpose_img
         roi = np.rot90(roi, 2, (1, 2))
     
     # get sitk image and transfer metadata to it
-    sitk_img = sitk.GetImageFromArray(roi)
+    sitk_img = convert_img(roi)
     spacing = config.resolutions[0]
     sitk_img.SetSpacing(spacing[::-1])
     # TODO: consider setting z-origin to 0 since image generally as
@@ -678,21 +784,57 @@ def load_numpy_to_sitk(
     return sitk_img
 
 
-def numpy_to_sitk(img: np.ndarray):
-    """Wrapper to convert a NumPy array to a SimpleITK Image object.
+def sitk_to_itk_img(sitk_img: "sitk.Image") -> "itk.Image":
+    """Convert a SimpleITK image to an ITK Image.
     
     Args:
-        img: NumPy array.
+        sitk_img: SimpleITK Image. If None, will simply be returned.
 
     Returns:
-        ``img`` as a SimpleITK Image object with default settings.
+        ITK Image.
 
     """
-    try:
-        img_sitk = sitk.GetImageFromArray(img)
-    except TypeError:
-        # WORKAROUND: sitk may not match img's dtype to sitk pixel IDs; match
-        # based on dtype name instead of dtype object identity
-        img_sitk = sitk.Image(
-            img.shape[::-1], _DTYPES_NP_SITK[np.dtype(img.dtype).name])
-    return img_sitk
+    if sitk_img is None: return sitk_img
+    
+    # construct an ITK Image through the ndarray extracted from the sitk Image
+    itk_img = itk.GetImageFromArray(
+        convert_img(sitk_img),
+        is_vector=sitk_img.GetNumberOfComponentsPerPixel() > 1)
+    
+    # transfer metadata
+    itk_img.SetOrigin(sitk_img.GetOrigin())
+    itk_img.SetSpacing(sitk_img.GetSpacing())
+    
+    # sitk Image directions are flattened arrays, but ITK expects an
+    # (ndim x ndim) matrix
+    itk_img.SetDirection(itk.GetMatrixFromArray(np.reshape(
+        np.array(sitk_img.GetDirection()), [sitk_img.GetDimension()] * 2)))
+    
+    return itk_img
+
+
+def itk_to_sitk_img(itk_img: "itk.Image") -> "sitk.Image":
+    """Convert an ITK image to a SimpleITK Image.
+    
+    Args:
+        itk_img: ITK Image. If None, will simply be returned.
+
+    Returns:
+        SimpleITK Image.
+
+    """
+    if itk_img is None: return itk_img
+    
+    # construct a sitk Image through the ndarray extracted from the ITK Image
+    sitk_img = sitk.GetImageFromArray(
+        itk.GetArrayFromImage(itk_img),
+        isVector=itk_img.GetNumberOfComponentsPerPixel() > 1)
+    
+    # transfer metadata
+    sitk_img.SetOrigin(tuple(itk_img.GetOrigin()))
+    sitk_img.SetSpacing(tuple(itk_img.GetSpacing()))
+    
+    # ITK Image directions are 2D arrays, but ITK expects a 1D array
+    sitk_img.SetDirection(np.ravel(np.array(itk_img.GetDirection())))
+    
+    return sitk_img
diff --git a/magmap/settings/config.py b/magmap/settings/config.py
index 78becfc2..5010181c 100644
--- a/magmap/settings/config.py
+++ b/magmap/settings/config.py
@@ -32,7 +32,14 @@
 from magmap.settings import logs
 
 if TYPE_CHECKING:
-    import SimpleITK as sitk
+    try:
+        import itk
+    except ImportError:
+        itk = None
+    try:
+        import SimpleITK as sitk
+    except ImportError:
+        sitk = None
     from magmap.atlas import labels_meta, ontology
     from magmap.cv import detector
     from magmap.io import np_io
@@ -607,8 +614,6 @@ class RegNames(Enum):
     COMBINED = "combined.mhd"  # spliced into other registered names
 
 
-#: Path to labels image metadata file. 
-
 #: Loaded labels metadata.
 labels_metadata: Optional["labels_meta.LabelsMeta"] = None
 
@@ -620,10 +625,10 @@ class RegNames(Enum):
 #: Numpy array of a labels image file, typically corresponding to ``img5d``.
 labels_img: Optional = None
 #: Labels image as a SimpleITK Image instance.
-labels_img_sitk: Optional["sitk.Image"] = None
+labels_img_sitk: Optional[Union["sitk.Image", "itk.Image"]] = None
 #: Original labels image, before any processing.
 labels_img_orig: Optional[np.ndarray] = None
-#: Scaling factors from ``labels_img`` to ``img5d``. 
+#: Scaling factors from ``labels_img`` to ``img5d``.
 labels_scaling: Optional[Sequence[float]] = None
 #: Labels reference IDs corresponding to the labels image values.
 labels_ref: Optional["ontology.LabelsRef"] = None
@@ -632,8 +637,8 @@ class RegNames(Enum):
 
 VOL_KEY = "volume"
 BLOBS_KEY = "blobs"
-VARIATION_BLOBS_KEY = "var_blobs" # variation in blob density
-VARIATION_EXP_KEY = "var_exp" # variation in experiment intensity
+VARIATION_BLOBS_KEY = "var_blobs"  # variation in blob density
+VARIATION_EXP_KEY = "var_exp"  # variation in experiment intensity
 GENOTYPE_KEY = "Geno"
 SUB_SEG_MULT = 100  # labels multiplier for sub-segmentations
 REGION_ALL = "all"
@@ -690,7 +695,6 @@ class ABAKeys(Enum):
         "MAKE_EDGE_IMAGES",
         "MAKE_EDGE_IMAGES_EXP",
         "MERGE_ATLAS_SEGS",
-        "REG_LABELS_TO_ATLAS",
         "VOL_STATS",
         "VOL_COMPARE",
         "MAKE_DENSITY_IMAGES",
diff --git a/magmap/stats/atlas_stats.py b/magmap/stats/atlas_stats.py
index e69a5ecf..b22d9715 100644
--- a/magmap/stats/atlas_stats.py
+++ b/magmap/stats/atlas_stats.py
@@ -1,5 +1,5 @@
 # Atlas measurements and statistics
-# Author: David Young, 2019, 2021
+# Author: David Young, 2019, 2023
 """Low-level measurement of atlases and statistics generation.
 
 Typically applied to specific types of atlases and less generalizable
@@ -10,7 +10,6 @@
 
 import numpy as np
 import pandas as pd
-import SimpleITK as sitk
 
 from magmap.atlas import atlas_refiner, ontology
 from magmap.io import df_io, export_stack, libmag, np_io, sitk_io
@@ -562,7 +561,7 @@ def meas_landmark_dist(
     
     # transpose 2nd image with any config settings
     labels_imgs[1] = atlas_refiner.transpose_img(labels_imgs[1])
-    labels_imgs = [sitk.GetArrayFromImage(m) for m in labels_imgs]
+    labels_imgs = [sitk_io.convert_img(m) for m in labels_imgs]
     
     # set up output path and sample name based on 1st image's path
     out_path = libmag.make_out_path(libmag.combine_paths(
diff --git a/setup.py b/setup.py
index 4715e234..2eab8cd3 100644
--- a/setup.py
+++ b/setup.py
@@ -50,6 +50,19 @@
 #: Optional dependencies for the 3D viewer.
 _EXTRAS_3D = ["mayavi"]
 
+#: Optional pre-built SimpleITK with Elastix for image I/O and registration.
+_EXTRAS_SIMPLEITK = [
+    "simpleitk==2.0.2rc2.dev785+g8ac4f ; python_version < '3.8'",
+    "simpleitk==2.3.0.dev117+g0640d ; python_version >= '3.8'",
+]
+
+#: Optional ITK and Elastix for image I/O and registration.
+_EXTRAS_ITK = [
+    "itk",
+    "itk-elastix",
+]
+
+
 # installation configuration
 config = {
     "name": "magellanmapper",
@@ -73,11 +86,6 @@
         # PlotEditor performance regression with 3.3.0-3.3.1
         "matplotlib != 3.3.0, != 3.3.1",
         "pandas",
-        
-        # pre-built SimpleITK with Elastix
-        "simpleitk==2.0.2rc2.dev785+g8ac4f ; python_version < '3.8'",
-        "simpleitk==2.3.0.dev117+g0640d ; python_version >= '3.8'",
-        
         "PyYAML",
         "appdirs",
         # part of stdlib in Python >= 3.8
@@ -101,12 +109,15 @@
         "classifier": _EXTRAS_CLASSIFER,
         "gui": _EXTRAS_GUI,
         "3d": _EXTRAS_3D,
+        "itk": _EXTRAS_ITK,
+        "simplitk": _EXTRAS_SIMPLEITK,
         
         # dependencies for most common tasks
         "most": [
             "matplotlib_scalebar",
             "pyamg",  # for Random-Walker segmentation "cg_mg" mode
             *_EXTRAS_GUI,
+            *_EXTRAS_ITK,
             *_EXTRAS_IMPORT,
         ],
         
@@ -118,6 +129,7 @@
             "scikit-learn",
             *_EXTRAS_GUI,
             *_EXTRAS_PANDAS,
+            *_EXTRAS_ITK,
             *_EXTRAS_IMPORT,
             *_EXTRAS_AWS,
             *_EXTRAS_JUPYTER,


Test patch:
diff --git a/magmap/tests/test_libmag.py b/magmap/tests/test_libmag.py
index 2036bdd8..2108454e 100644
--- a/magmap/tests/test_libmag.py
+++ b/magmap/tests/test_libmag.py
@@ -4,12 +4,44 @@
 
 import unittest
 
+import numpy as np
 
 from magmap.io import libmag
 
 
 class TestLibmag(unittest.TestCase):
     
+    def test_pad_seq(self):
+        seq = ["a", "b", "c"]
+        self.assertSequenceEqual(libmag.pad_seq(list(seq), 2), ["a", "b"])
+        self.assertSequenceEqual(libmag.pad_seq(list(seq), 3), ["a", "b", "c"])
+        self.assertSequenceEqual(
+            libmag.pad_seq(list(seq), 4), ["a", "b", "c", None])
+        self.assertSequenceEqual(
+            libmag.pad_seq(list(seq), 4, 1), ["a", "b", "c", 1])
+        self.assertSequenceEqual(
+            libmag.pad_seq(list(seq), 5, [1, 2, 3, 4, 5]), ["a", "b", "c", 4, 5])
+        
+        # fills with last pad value since pad is shorter than seq is
+        self.assertSequenceEqual(
+            libmag.pad_seq(list(seq), 5, [1, 2]), ["a", "b", "c", 2, 2])
+        
+        # modifies the original sequence
+        self.assertSequenceEqual(libmag.pad_seq(tuple(seq), 4), ["a", "b", "c", None])
+        self.assertSequenceEqual(seq, ["a", "b", "c"])
+        
+        # modifies the original sequence
+        self.assertSequenceEqual(libmag.pad_seq(seq, 4), ["a", "b", "c", None])
+        self.assertSequenceEqual(seq, ["a", "b", "c", None])
+    
+    def test_replace_seq(self):
+        seq = [1, 2, 3]
+        np.testing.assert_array_equal(
+            libmag.replace_seq(seq, [10, 20]), [10, 20, 3])
+        np.testing.assert_array_equal(
+            libmag.replace_seq(seq, [10, 20, 30, 40]), [10, 20, 30])
+        np.testing.assert_array_equal(libmag.replace_seq(seq, []), [1, 2, 3])
+    
     def test_flatten(self):
         self.assertSequenceEqual(
             list(libmag.flatten(["a", "b"])), ["a", "b"])