Make SimpleITK dependency optional

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
 

docs/release/release_v1.6.md

@@ -81,9 +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 (#495)
-  - `ITK-Elastix` support has been added as an alternative to `SimpleITK`
-  - Both libraries are now optional rather than required dependencies
+- Supports image I/O and registration through ITK
+  - `ITK-Elastix` support has been added as an alternative to `SimpleITK` (#495)
+  - 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

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)

magmap/atlas/register.py

@@ -216,10 +216,10 @@ def _handle_transform_file(fixed_file, transform_param_map=None):
 
 
 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
@@ -244,8 +244,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:
@@ -268,7 +268,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:
@@ -714,7 +714,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)
@@ -761,7 +761,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)
@@ -783,11 +783,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:
@@ -942,7 +942,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:
@@ -972,7 +972,7 @@ def make_labels(lbls_img):
         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])
@@ -1056,14 +1056,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(
@@ -1111,7 +1111,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)
 
 
@@ -1233,7 +1233,7 @@ def register_group(
         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)
+        img_np = sitk_io.convert_img(img)
         if img_np_template is None:
             img_np_template = np.copy(img_np)
 
@@ -1315,12 +1315,12 @@ def register_group(
     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)
+        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:
+        if show_imgs and sitk:
             sitk.Show(sitk_io.replace_sitk_with_numpy(img, img_large_np))
         imgs.append(img_large_np)
 
@@ -1357,7 +1357,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
@@ -1367,9 +1367,8 @@ 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)
+        sitk_io.write_img(transformed_img, out_path)
+        img_np = sitk_io.convert_img(transformed_img)
         config.resolutions = [transformed_img.GetSpacing()[::-1]]
         importer.save_np_image(img_np[None], out_path, config.series)
 
@@ -1424,7 +1423,8 @@ def register_labels_to_atlas(path_fixed):
     # perform the registration
     transform = elastix_img_filter.Execute()
     transformed_img = elastix_img_filter.GetResultImage()
-    sitk.Show(transformed_img)
+    if sitk:
+        sitk.Show(transformed_img)
 
     # set up filter to apply the same transformation to label file, 
     # ensuring that no new labels are interpolated
@@ -1477,18 +1477,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(
@@ -1760,7 +1759,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
@@ -1951,7 +1950,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
@@ -2101,17 +2100,17 @@ 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:
         # registered labels image and associated experiment file
         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)
@@ -2137,10 +2136,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():