|
4 | 4 | import numpy as np |
5 | 5 | import pandas as pd |
6 | 6 | import matplotlib.pyplot as plt |
| 7 | +import tifffile |
| 8 | +from matplotlib import colors |
| 9 | +from skimage.segmentation import find_boundaries |
7 | 10 |
|
8 | 11 | from util import literature_reference_values |
9 | 12 |
|
10 | 13 | png_dpi = 300 |
11 | 14 |
|
12 | 15 |
|
| 16 | +def scramble_instance_labels(arr): |
| 17 | + """Scramble indexes of instance segmentation to avoid neighboring colors. |
| 18 | + """ |
| 19 | + unique = list(np.unique(arr)[1:]) |
| 20 | + rng = np.random.default_rng(seed=42) |
| 21 | + new_list = rng.uniform(1, len(unique) + 1, size=(len(unique))) |
| 22 | + new_arr = np.zeros(arr.shape) |
| 23 | + for old_id, new_id in zip(unique, new_list): |
| 24 | + new_arr[arr == old_id] = new_id |
| 25 | + return new_arr |
| 26 | + |
| 27 | + |
| 28 | +def plot_seg_crop(img_path, seg_path, save_path, xlim1, xlim2, ylim1, ylim2, boundary_rgba=[0, 0, 0, 0.5], plot=False): |
| 29 | + seg = tifffile.imread(seg_path) |
| 30 | + if len(seg.shape) == 3: |
| 31 | + seg = seg[10, xlim1:xlim2, ylim1:ylim2] |
| 32 | + else: |
| 33 | + seg = seg[xlim1:xlim2, ylim1:ylim2] |
| 34 | + |
| 35 | + img = tifffile.imread(img_path) |
| 36 | + img = img[10, xlim1:xlim2, ylim1:ylim2] |
| 37 | + |
| 38 | + # create color map with random distribution for coloring instance segmentation |
| 39 | + unique = list(np.unique(seg)[1:]) |
| 40 | + n_instances = len(unique) |
| 41 | + |
| 42 | + seg = scramble_instance_labels(seg) |
| 43 | + |
| 44 | + rng = np.random.default_rng(seed=42) # fixed seed for reproducibility |
| 45 | + colors_array = rng.uniform(0, 1, size=(n_instances, 4)) # RGBA values in [0,1] |
| 46 | + colors_array[:, 3] = 1.0 # full alpha |
| 47 | + colors_array[0, 3] = 0.0 # make label 0 transparent (background) |
| 48 | + cmap = colors.ListedColormap(colors_array) |
| 49 | + |
| 50 | + boundaries = find_boundaries(seg, mode="inner") |
| 51 | + boundary_overlay = np.zeros((*boundaries.shape, 4)) |
| 52 | + |
| 53 | + boundary_overlay[boundaries] = boundary_rgba # RGBA = black |
| 54 | + |
| 55 | + fig, ax = plt.subplots(figsize=(6, 6)) |
| 56 | + ax.imshow(img, cmap="gray") |
| 57 | + ax.imshow(seg, cmap=cmap, alpha=0.5, interpolation="nearest") |
| 58 | + ax.imshow(boundary_overlay) |
| 59 | + ax.axis("off") |
| 60 | + plt.tight_layout() |
| 61 | + plt.savefig(save_path, bbox_inches="tight", pad_inches=0.1, dpi=png_dpi) |
| 62 | + |
| 63 | + if plot: |
| 64 | + plt.show() |
| 65 | + else: |
| 66 | + plt.close() |
| 67 | + |
| 68 | + |
| 69 | +def fig_02b_sgn(save_dir, plot=False): |
| 70 | + """Plot crops of SGN segmentation of CochleaNet, Cellpose and micro-sam. |
| 71 | + """ |
| 72 | + cochlea_dir = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet" |
| 73 | + val_sgn_dir = f"{cochlea_dir}/predictions/val_sgn" |
| 74 | + image_dir = f"{cochlea_dir}/AnnotatedImageCrops/F1ValidationSGNs/for_consensus_annotation" |
| 75 | + |
| 76 | + crop_name = "MLR169R_PV_z3420_allturns_full" |
| 77 | + img_path = os.path.join(image_dir, f"{crop_name}.tif") |
| 78 | + |
| 79 | + xlim1 = 2000 |
| 80 | + xlim2 = 2500 |
| 81 | + ylim1 = 3100 |
| 82 | + ylim2 = 3600 |
| 83 | + boundary_rgba = [1, 1, 1, 0.5] |
| 84 | + |
| 85 | + for seg_net in ["distance_unet", "cellpose-sam", "micro-sam"]: |
| 86 | + save_path = os.path.join(save_dir, f"fig_02b_sgn_{seg_net}.png") |
| 87 | + seg_dir = os.path.join(val_sgn_dir, seg_net) |
| 88 | + seg_path = os.path.join(seg_dir, f"{crop_name}_seg.tif") |
| 89 | + |
| 90 | + plot_seg_crop(img_path, seg_path, save_path, xlim1, xlim2, ylim1, ylim2, boundary_rgba, plot=plot) |
| 91 | + |
| 92 | + |
| 93 | +def fig_02b_ihc(save_dir, plot=False): |
| 94 | + """Plot crops of IHC segmentation of CochleaNet, Cellpose and micro-sam. |
| 95 | + """ |
| 96 | + cochlea_dir = "/mnt/vast-nhr/projects/nim00007/data/moser/cochlea-lightsheet" |
| 97 | + val_sgn_dir = f"{cochlea_dir}/predictions/val_ihc" |
| 98 | + image_dir = f"{cochlea_dir}/AnnotatedImageCrops/F1ValidationIHCs" |
| 99 | + |
| 100 | + crop_name = "MLR226L_VGlut3_z1200_3turns_full" |
| 101 | + img_path = os.path.join(image_dir, f"{crop_name}.tif") |
| 102 | + |
| 103 | + xlim1 = 1900 |
| 104 | + xlim2 = 2400 |
| 105 | + ylim1 = 2000 |
| 106 | + ylim2 = 2500 |
| 107 | + boundary_rgba = [1, 1, 1, 0.5] |
| 108 | + |
| 109 | + for seg_net in ["distance_unet_v4b", "cellpose-sam", "micro-sam"]: |
| 110 | + save_path = os.path.join(save_dir, f"fig_02b_ihc_{seg_net}.png") |
| 111 | + seg_dir = os.path.join(val_sgn_dir, seg_net) |
| 112 | + seg_path = os.path.join(seg_dir, f"{crop_name}_seg.tif") |
| 113 | + |
| 114 | + plot_seg_crop(img_path, seg_path, save_path, xlim1, xlim2, ylim1, ylim2, boundary_rgba, plot=plot) |
| 115 | + |
| 116 | + |
13 | 117 | def fig_02c(save_path, plot=False, all_versions=False): |
14 | 118 | """Scatter plot showing the precision, recall, and F1-score of SGN (distance U-Net, manual), |
15 | 119 | IHC (distance U-Net, manual), and synapse detection (U-Net). |
@@ -299,6 +403,8 @@ def main(): |
299 | 403 | os.makedirs(args.figure_dir, exist_ok=True) |
300 | 404 |
|
301 | 405 | # Panel C: Evaluation of the segmentation results: |
| 406 | + fig_02b_sgn(save_dir=args.figure_dir, plot=args.plot) |
| 407 | + fig_02b_ihc(save_dir=args.figure_dir, plot=args.plot) |
302 | 408 | fig_02c(save_path=os.path.join(args.figure_dir, "fig_02c"), plot=args.plot, all_versions=False) |
303 | 409 |
|
304 | 410 | # Panel D: The number of SGNs, IHCs and average number of ribbon synapses per IHC |
|
0 commit comments