Mappings from physical components (phase,retardance,orientation) to HSV /JCh

examples/map_to_hsv.py

@@ -0,0 +1,84 @@
+from mantis.analysis import visualization as vz
+from pathlib import Path
+from iohub import open_ome_zarr
+import numpy as np
+import glob
+from natsort import natsorted
+from mantis.cli import utils
+
+if __name__ == "__main__":
+    HSV_method = "JCh"  # "RO" or "PRO" or "JCh"
+    hsv_channels = ["Retardance - recon", "Orientation - recon"]
+    num_processes = 10
+    input_data_paths = "/hpc/projects/comp.micro/zebrafish/2023_02_02_zebrafish_casper/2-prototype-reconstruction/fov6-reconstruction.zarr/*/*/*"
+    output_data_path = f"./test_{HSV_method}_new_method2.zarr"
+
+    input_data_paths = [Path(path) for path in natsorted(glob.glob(input_data_paths))]
+    output_data_path = Path(output_data_path)
+
+    # Taking the input sample
+    with open_ome_zarr(input_data_paths[0], mode="r") as dataset:
+        T, C, Z, Y, X = dataset.data.shape
+        dataset_scale = dataset.scale
+        channel_names = dataset.channel_names
+
+        input_channel_idx = []
+        # FIXME: these should be part of a config
+        # FIXME:this is hardcoded to spit out 3 chans for RGB
+        output_channel_idx = [0, 1, 2]
+        time_indices = list(range(T))
+
+        if HSV_method == "PRO":
+            # hsv_channels = ["Orientation", "Retardance", "Phase"]
+            HSV_func = vz.HSV_PRO
+            hsv_func_kwargs = dict(
+                channel_order=output_channel_idx, max_val_V=0.5, max_val_S=1.0
+            )
+
+        elif HSV_method == "RO":
+            # hsv_channels = [
+            #     "Orientation",
+            #     "Retardance",
+            # ]
+            HSV_func = vz.HSV_RO
+            hsv_func_kwargs = dict(channel_order=output_channel_idx, max_val_V=0.5)
+
+        elif HSV_method == "JCh":
+            # hsv_channels = ["Orientation", "Retardance"]
+            HSV_func = vz.JCh_mapping
+            hsv_func_kwargs = dict(
+                channel_order=output_channel_idx, max_val_ret=150, noise_level=1
+            )
+        for channel in hsv_channels:
+            if channel in channel_names:
+                input_channel_idx.append(channel_names.index(channel))
+        rgb_channel_names = ["Red", "Green", "Blue"]
+
+    # Here the functions will output an RGB image
+    output_metadata = {
+        "shape": (len(time_indices), len(rgb_channel_names), Z, Y, X),
+        "chunks": None,
+        "scale": dataset_scale,
+        "channel_names": rgb_channel_names,
+        "dtype": np.float32,
+    }
+
+    utils.create_empty_hcs_zarr(
+        store_path=output_data_path,
+        position_keys=[p.parts[-3:] for p in input_data_paths],
+        **output_metadata,
+    )
+
+    for input_position_path in input_data_paths:
+        utils.process_single_position_v2(
+            HSV_func,
+            input_data_path=input_position_path,  # source store
+            output_path=output_data_path,  # target store
+            time_indices=time_indices,
+            input_channel_idx=input_channel_idx,
+            output_channel_idx=output_channel_idx,
+            num_processes=num_processes,  # parallel processing over time
+            **hsv_func_kwargs,
+        )
+
+# %%

mantis/analysis/register.py

@@ -267,7 +267,7 @@ def find_lir(registered_zyx: np.ndarray, plot: bool = False) -> Tuple:
     # Find the lir Z
     zyx_shape = registered_zyx.shape
     registered_zx_bool = registered_zyx.transpose((2, 0, 1)) > 0
-    registered_zx_bool = registered_zx_bool[zyx_shape[0] // 2].copy()
+    registered_zx_bool = registered_zx_bool[zyx_shape[2] // 2].copy()
     rectangle_coords_zx = lir.lir(registered_zx_bool)
     x = rectangle_coords_zx[0]
     z = rectangle_coords_zx[1]

mantis/analysis/visualization.py

@@ -0,0 +1,147 @@
+import numpy as np
+
+from colorspacious import cspace_convert
+from skimage.color import hsv2rgb
+from skimage.exposure import rescale_intensity
+
+
+def HSV_PRO(czyx, channel_order, max_val_V: float = 1.0, max_val_S: float = 1.0):
+    """
+    HSV encoding of retardance + orientation + phase image with hsv colormap (orientation in h, retardance in s, phase in v)
+    Parameters
+    ----------
+        czyx : numpy.ndarray
+        channel_order: list
+                    the order in which the channels should be stacked i.e([orientation_c_idx, retardance_c_idx,phase_c_idx])
+                    the 0 index corresponds to the orientation image (range from 0 to pi)
+                    the 1 index corresponds to the retardance image
+                    the 2 index corresponds to the Phase image
+        max_val_V   : float
+                      raise the brightness of the phase channel by 1/max_val_V
+        max_val_S   : float
+                      raise the brightness of the retardance channel by 1/max_val_S
+    Returns:
+        RGB with HSV (Orientation, Retardance, Phase)
+    """
+
+    C, Z, Y, X = czyx.shape
+    assert C == 3, "The input array must have 3 channels"
+    print(f"channel_order: {channel_order}")
+
+    czyx_out = np.zeros((3, Z, Y, X), dtype=np.float32)
+    # Normalize the stack
+    ordered_stack = np.stack(
+        (
+            # Normalize the first channel by dividing by pi
+            czyx[channel_order[0]] / np.pi,
+            # Normalize the second channel and rescale intensity
+            rescale_intensity(
+                czyx[channel_order[1]],
+                in_range=(
+                    np.min(czyx[channel_order[1]]),
+                    np.max(czyx[channel_order[1]]),
+                ),
+                out_range=(0, 1),
+            )
+            / max_val_S,
+            # Normalize the third channel and rescale intensity
+            rescale_intensity(
+                czyx[channel_order[2]],
+                in_range=(
+                    np.min(czyx[channel_order[2]]),
+                    np.max(czyx[channel_order[2]]),
+                ),
+                out_range=(0, 1),
+            )
+            / max_val_V,
+        ),
+        axis=0,
+    )
+    czyx_out = hsv2rgb(ordered_stack, channel_axis=0)
+    return czyx_out
+
+
+def HSV_RO(czyx, channel_order: list[int], max_val_V: int = 1):
+    """
+    visualize retardance + orientation with hsv colormap (orientation in h, saturation=1 s, retardance in v)
+    Parameters
+    ----------
+        czyx : numpy.ndarray
+        channel_order: list
+                    the order in which the channels should be stacked i.e([orientation_c_idx, retardance_c_idx])
+                    the 0 index corresponds to the orientation image (range from 0 to pi)
+                    the 1 index corresponds to the retardance image
+        max_val_V   : float
+                    raise the brightness of the phase channel by 1/max_val_V
+    Returns:
+
+        RGB with HSV (Orientation, _____ , Retardance)
+    """
+    C, Z, Y, X = czyx.shape
+    assert C == 2, "The input array must have 2 channels"
+    czyx_out = np.zeros((3, Z, Y, X), dtype=np.float32)
+    ordered_stack = np.stack(
+        (
+            # Normalize the first channel by dividing by pi and then rescale intensity
+            czyx[channel_order[0]] / np.pi,
+            # Set the second channel to ones = Saturation 1
+            np.ones_like(czyx[channel_order[0]]),
+            # Normalize the third channel and rescale intensity
+            np.minimum(
+                1,
+                rescale_intensity(
+                    czyx[channel_order[1]],
+                    in_range=(
+                        np.min(czyx[channel_order[1]]),
+                        np.max(czyx[channel_order[1]]),
+                    ),
+                    out_range=(0, max_val_V),
+                ),
+            ),
+        ),
+        axis=0,
+    )
+    # HSV-RO encoding
+    czyx_out = hsv2rgb(ordered_stack, channel_axis=0)
+    return czyx_out
+
+
+def JCh_mapping(czyx, channel_order: list[int], max_val_ret: int = None, noise_level: int = 1):
+    """
+    JCh retardance + orientation + phase image with hsv colormap (orientation in h, retardance in s, phase in v)
+    Parameters
+    ----------
+        czyx : numpy.ndarray
+        channel_order: list
+                     the order in which the channels should be stacked i.e([retardance_c_idx, orientation_c_idx])
+                      the 0 index corresponds to the retardance image
+                      the 1 index corresponds to the orientation image (range from 0 to pi)
+
+        max_val_V   : float
+                      raise the brightness of the phase channel by 1/max_val_ret
+    Returns:
+        RGB with JCh (Retardance, Orientation)
+    """
+    # retardance, orientation
+    C, Z, Y, X = czyx.shape
+    assert C == 2, "The input array must have 2 channels"
+
+    # Retardance,chroma,Hue
+    czyx_out = np.zeros((3, Z, Y, X), dtype=np.float32)
+    for z_idx in range(Z):
+        # Retardance
+        if max_val_ret is None:
+            max_val_ret = np.max(czyx[channel_order[0], z_idx])
+        retardance = np.clip(czyx[channel_order[0], z_idx], 0, max_val_ret)
+        # Chroma of each pixel, set to 60 by default, with noise handling
+        chroma = np.where(czyx[channel_order[0], z_idx] < noise_level, 0, 60)
+        # Orientation 180 to 360 to match periodic hue
+        hue = czyx[channel_order[1], z_idx] * 360 / np.pi
+        # Stack arrays in the correct order (Y, X, 3)
+        I_JCh = np.stack((retardance, chroma, hue), axis=-1)
+        # Transpose to shape for the skimage or colorspace functions
+        JCh_rgb = cspace_convert(I_JCh, "JCh", "sRGB1")
+        JCh_rgb = np.clip(JCh_rgb, 0, 1)
+        czyx_out[:, z_idx] = np.transpose(JCh_rgb, (2, 0, 1))
+
+    return czyx_out

mantis/cli/apply_affine.py

@@ -108,7 +108,8 @@ def apply_affine(
         )
         # Overwrite the previous target shape
         Z_target, Y_target, X_target = cropped_target_shape_zyx[-3:]
-        click.echo(f'Shape of cropped output dataset: {target_shape_zyx}\n')
+        cropped_target_shape_zyx = Z_target, Y_target, X_target
+        click.echo(f'Shape of cropped output dataset: {cropped_target_shape_zyx}\n')
     else:
         Z_slice, Y_slice, X_slice = (
             slice(0, Z_target),