   1# -*- coding: utf-8 -*-
      2"""
      3Created on Thu Apr 13 14:02:27 2023
      4
      5@author: Bart Steemans. Govers Lab.
      6"""
      7
      8import matplotlib.pyplot as plt
      9import numpy as np
     10import pandas as pd
     11import random
     12from pathlib import Path
     13from matplotlib.colors import LinearSegmentedColormap, ListedColormap
     14from matplotlib.lines import Line2D
     15from skimage.color import label2rgb
     16from skimage.measure import label
     17from . import utilities as u
     18import tifffile
     19from skimage.segmentation import find_boundaries
     20from scipy.spatial import cKDTree
     21dpi = 300
     22
     23import numpy as np
     24import matplotlib.pyplot as plt
     25
     26
     27def _tutorial_contrast(image):
     28    """Return a display range without changing the measured image."""
     29    low, high = np.nanpercentile(image, (1, 99.7))
     30    return (low, high if high > low else low + 1)
     31
     32
     33def _tutorial_crop(shape, center, window_px):
     34    half = window_px // 2
     35    y = int(np.clip(center[0] - half, 0, max(shape[0] - window_px, 0)))
     36    x = int(np.clip(center[1] - half, 0, max(shape[1] - window_px, 0)))
     37    return slice(y, min(y + window_px, shape[0])), slice(x, min(x + window_px, shape[1]))
     38
     39
     40def pick_representative_cell(image_objects, margin_px=20):
     41    """Pick an interior, isolated, median-sized cell with valid geometry.
     42
     43    The rule is deterministic and uses geometry only, never fluorescence intensity.
     44    Returns the image object and cell ID so the same cell can be revisited.
     45    """
     46    candidates = []
     47    for image in image_objects:
     48        height, width = image.image.shape
     49        for cell in getattr(image, "cells", []):
     50            contour = np.asarray(cell.contour)
     51            mesh = np.asarray(cell.mesh)
     52            if (contour.ndim != 2 or contour.shape[1] != 2 or len(contour) < 6
     53                    or mesh.ndim != 2 or mesh.shape[1] != 4 or len(mesh) < 4
     54                    or not np.isfinite(contour).all() or not np.isfinite(mesh).all()):
     55                continue
     56            low = contour.min(axis=0)
     57            high = contour.max(axis=0)
     58            if (low[0] < margin_px or low[1] < margin_px
     59                    or high[0] > height - margin_px or high[1] > width - margin_px):
     60                continue
     61            box = high - low
     62            long_side, short_side = max(box), min(box)
     63            if short_side < 5 or long_side / short_side < 1.5:
     64                continue
     65            center = contour.mean(axis=0)
     66            center_distance = np.linalg.norm(center / [height, width] - 0.5)
     67            candidates.append((image, cell.cell_id, long_side, center_distance, center))
     68    if not candidates:
     69        raise ValueError("No interior cells with valid elongated contours and meshes")
     70    median_length = np.median([item[2] for item in candidates])
     71    typical = [item for item in candidates if abs(item[2] - median_length) <= 0.25 * median_length]
     72    neighbors = {}
     73    for image in image_objects:
     74        group = [item for item in candidates if item[0] is image]
     75        if not group:
     76            continue
     77        all_centers = [np.asarray(cell.contour).mean(axis=0)
     78                       for cell in image.cells
     79                       if np.asarray(cell.contour).ndim == 2
     80                       and np.asarray(cell.contour).shape[1] == 2
     81                       and len(cell.contour) > 0]
     82        if len(all_centers) < 2:
     83            for item in group:
     84                neighbors[(image.image_name, item[1])] = np.inf
     85            continue
     86        distances = cKDTree(all_centers).query([item[4] for item in group], k=2)[0]
     87        for item, distance in zip(group, distances):
     88            neighbors[(image.image_name, item[1])] = distance[1] / item[2]
     89    winner = min(
     90        typical,
     91        key=lambda item: (-neighbors[(item[0].image_name, item[1])],
     92                          item[3], abs(item[2] - median_length),
     93                          item[0].image_name, item[1]),
     94    )
     95    return winner[0], winner[1]
     96
     97
     98def select_example_cell(image_objects, field_id=None, cell_id=None):
     99    """Choose an existing field/cell, or a representative cell in the chosen field."""
    100    objects = list(image_objects)
    101    if field_id is not None:
    102        objects = [image for image in objects if image.image_name.rsplit("_", 1)[0] == field_id]
    103        if not objects:
    104            raise ValueError(f"Field {field_id!r} was not found")
    105    if cell_id is None:
    106        return pick_representative_cell(objects)
    107    if field_id is None:
    108        raise ValueError("Set field_id when selecting a cell_id")
    109    for image in objects:
    110        if any(cell.cell_id == cell_id for cell in image.cells):
    111            return image, cell_id
    112    raise ValueError(f"Cell {cell_id} was not found in field {field_id}")
    113
    114
    115def _tutorial_channel_color(label):
    116    label = label.casefold()
    117    if "nucleoid" in label or "dapi" in label:
    118        return "#1DBCD3"
    119    if "rna" in label or "gfp" in label:
    120        return "#2BAE66"
    121    if "fm" in label or "membrane" in label:
    122        return "#D85755"
    123    return "#197E89"
    124
    125
    126def _tutorial_channel_cmap(label):
    127    color = _tutorial_channel_color(label)
    128    return LinearSegmentedColormap.from_list(
    129        f"bactoscoop_{label.replace(' ', '_')}", ["#000000", color]
    130    )
    131
    132
    133def plot_field_channels(dataset_dir, field_id, channel_labels, crop_size=650):
    134    """Preview the same pixel window across the named TIFF channels."""
    135    dataset_dir = Path(dataset_dir)
    136    labels = dict(channel_labels)
    137    if not labels:
    138        raise ValueError("channel_labels is empty")
    139    first = tifffile.imread(dataset_dir / f"{field_id}_{next(iter(labels))}.tiff")
    140    crop = _tutorial_crop(first.shape, (first.shape[0] // 2, first.shape[1] // 2), crop_size)
    141    fig, axes = plt.subplots(1, len(labels), figsize=(3.8 * len(labels), 4))
    142    for ax, (channel, label) in zip(np.atleast_1d(axes), labels.items()):
    143        image = tifffile.imread(dataset_dir / f"{field_id}_{channel}.tiff")
    144        if image.shape != first.shape:
    145            raise ValueError(f"{channel} dimensions do not match the first channel")
    146        data = image[crop]
    147        low, high = _tutorial_contrast(data)
    148        cmap = "gray" if channel == next(iter(labels)) else _tutorial_channel_cmap(label)
    149        ax.imshow(data, cmap=cmap, vmin=low, vmax=high)
    150        ax.set_title(f"{channel}  {label.upper()}", loc="left", fontsize=10, fontweight="bold", color="#193348")
    151        ax.set_xticks([])
    152        ax.set_yticks([])
    153        for spine in ax.spines.values():
    154            spine.set_visible(False)
    155    fig.suptitle(f"CHANNELS  /  {field_id}", x=0.02, y=0.98, ha="left", fontsize=14, fontweight="bold", color="#193348")
    156    fig.tight_layout(rect=(0, 0, 1, 0.98))
    157    return fig
    158
    159
    160def plot_segmentation_review(image_objects, crop_size=500, mask_label=None, center=None):
    161    """Compare phase and mask outlines in a selectable cell-centered crop."""
    162    objects = list(image_objects)
    163    if not objects:
    164        raise ValueError("image_objects is empty")
    165    fig, axes = plt.subplots(len(objects), 2, figsize=(10, 4.5 * len(objects)), squeeze=False)
    166    for row, image in enumerate(objects):
    167        if image.mask is None:
    168            raise ValueError(f"No mask is loaded for {image.image_name}")
    169        labels = np.unique(image.mask)
    170        labels = labels[labels > 0]
    171        if mask_label is not None and mask_label not in labels:
    172            raise ValueError(f"Mask label {mask_label} was not found in {image.image_name}")
    173        chosen = mask_label
    174        if chosen is None and center is None and len(labels):
    175            midpoint = np.asarray(image.image.shape) / 2
    176            nonzero = np.argwhere(image.mask > 0)
    177            closest = np.argmin(np.sum((nonzero - midpoint) ** 2, axis=1))
    178            chosen = image.mask[tuple(nonzero[closest])]
    179        if chosen is not None:
    180            focus = np.argwhere(image.mask == chosen).mean(axis=0)
    181        else:
    182            focus = center if center is not None else np.asarray(image.image.shape) / 2
    183        crop = _tutorial_crop(image.image.shape, focus, crop_size)
    184        phase, mask = image.image[crop], image.mask[crop]
    185        low, high = _tutorial_contrast(phase)
    186        for ax in axes[row]:
    187            ax.imshow(phase, cmap="gray", vmin=low, vmax=high)
    188            ax.set_xticks([])
    189            ax.set_yticks([])
    190            for spine in ax.spines.values():
    191                spine.set_visible(False)
    192        outlines = np.ma.masked_where(~find_boundaries(mask, mode="outer"), np.ones(mask.shape))
    193        axes[row, 1].imshow(outlines, cmap="autumn", vmin=0, vmax=1, alpha=0.9, interpolation="nearest")
    194        if chosen is not None:
    195            selected_outline = find_boundaries(mask == chosen, mode="outer")
    196            axes[row, 1].imshow(np.ma.masked_where(~selected_outline, np.ones(mask.shape)),
    197                                cmap="autumn", vmin=0, vmax=1,
    198                                interpolation="nearest")
    199        field = image.image_name.rsplit("_", 1)[0]
    200        axes[row, 0].set_title(f"{field}  |  phase", loc="left", fontsize=11)
    201        axes[row, 1].set_title(f"Mask outlines  |  selected label {chosen}", loc="left", fontsize=11)
    202    fig.suptitle("SEGMENTATION REVIEW", x=0.04, y=0.98, ha="left", fontsize=15, fontweight="bold", color="#193348")
    203    fig.tight_layout(rect=(0, 0, 1, 0.91))
    204    return fig
    205
    206
    207def plot_cell_gallery(image_object, cell_id, channel="C2", window_px=110,
    208                      channel_label=None, align_signal=False, show_objects=False):
    209    """Show phase and mesh, optionally with a registered signal crop and objects.
    210
    211    The registered crop uses the same local crop/phase registration called by
    212    ``object_detection(align=True)``. Plot it after detection to compare its
    213    detected object contours with the phase-cell contour.
    214    """
    215    cell = next((item for item in image_object.cells if item.cell_id == cell_id), None)
    216    if cell is None:
    217        raise ValueError(f"Cell {cell_id} was not found")
    218    if channel is not None and (not image_object.channels or channel not in image_object.channels):
    219        raise ValueError(f"Channel {channel} is not loaded")
    220    contour = np.asarray(cell.contour)
    221    mesh = np.asarray(cell.mesh)
    222    midline = np.asarray(cell.midline)
    223    extent = int(np.ceil(np.max(np.ptp(contour, axis=0)))) + 24
    224    crop = _tutorial_crop(image_object.image.shape, contour.mean(axis=0), max(window_px, extent))
    225    phase = image_object.image[crop]
    226    label = channel_label or channel
    227    panels = [(phase, "gray", "01  PHASE"), (phase, "gray", "02  CONTOUR + MESH")]
    228    signal_origin = (crop[1].start, crop[0].start)
    229    if channel is not None:
    230        signal = image_object.channels[channel][crop]
    231        if align_signal:
    232            signal, _, _, row_offset, column_offset = u.crop_image(
    233                image_object.channels[channel], cell.contour,
    234                phase_img=image_object.image,
    235            )
    236            signal_origin = (column_offset, row_offset)
    237        heading = "ALIGNED " if align_signal else ""
    238        panels.append((signal, _tutorial_channel_cmap(label),
    239                       f"03  {heading}{channel} {label.upper()}"))
    240    fig, axes = plt.subplots(1, len(panels), figsize=(3.7 * len(panels), 3.9))
    241    for ax, (data, cmap, title) in zip(axes, panels):
    242        low, high = _tutorial_contrast(data)
    243        ax.imshow(data, cmap=cmap, vmin=low, vmax=high)
    244        ax.set_title(title, loc="left", fontsize=10, fontweight="bold", color="#193348")
    245        ax.set_xticks([])
    246        ax.set_yticks([])
    247        for spine in ax.spines.values():
    248            spine.set_visible(False)
    249    x0, y0 = crop[1].start, crop[0].start
    250    axes[1].plot(contour[:, 1] - x0, contour[:, 0] - y0, color="#FFCC66", lw=1.6)
    251    if channel is not None:
    252        axes[2].plot(contour[:, 1] - signal_origin[0], contour[:, 0] - signal_origin[1], color="#FFCC66", lw=1.6)
    253        if show_objects:
    254            for detected in cell.object_meshdata.get(channel, {}).get("object_contour", []):
    255                detected = np.asarray(detected)
    256                if detected.ndim == 2 and detected.shape[1] == 2:
    257                    axes[2].plot(detected[:, 1] - signal_origin[0],
    258                                 detected[:, 0] - signal_origin[1],
    259                                 color="#4DE6F2", lw=1.2)
    260    if mesh.ndim == 2 and mesh.shape[1] == 4:
    261        axes[1].plot(mesh[:, 1] - x0, mesh[:, 0] - y0, color="#36D6D0", lw=0.7)
    262        axes[1].plot(mesh[:, 3] - x0, mesh[:, 2] - y0, color="#36D6D0", lw=0.7)
    263        for segment in mesh[::3]:
    264            axes[1].plot([segment[1] - x0, segment[3] - x0], [segment[0] - y0, segment[2] - y0], color="#36D6D0", lw=0.55)
    265    if midline.ndim == 2 and midline.shape[1] == 2:
    266        axes[1].plot(midline[:, 1] - x0, midline[:, 0] - y0, color="#F26463", lw=1.3)
    267    max_scale_px = phase.shape[1] * 0.45
    268    scale_um = next(
    269        (value for value in (5, 2, 1, 0.5, 0.2, 0.1)
    270         if value / image_object.px <= max_scale_px),
    271        0.1,
    272    )
    273    scale_px = max(1, round(scale_um / image_object.px))
    274    axes[0].plot([7, 7 + scale_px], [phase.shape[0] - 9] * 2, color="white", lw=3, solid_capstyle="butt")
    275    axes[0].text(7, phase.shape[0] - 14, f"{scale_um:g} µm", color="white", fontsize=8)
    276    fig.suptitle(f"CELL {cell_id}  /  {image_object.image_name}", x=0.02, y=0.98, ha="left", fontsize=13, fontweight="bold", color="#193348")
    277    fig.tight_layout(rect=(0, 0, 1, 0.93))
    278    return fig
    279
    280def plot_C2_alignment(image_object, cell_id, window_px=110,
    281                      channel="C2", channel_label=None,
    282                      show_objects=False):
    283    """
    284    Compare the raw C2 signal with the phase-aligned C2 signal.
    285
    286    BEFORE: original fluorescence image, cropped around the cell.
    287    AFTER:  fluorescence crop returned by u.crop_image(..., phase_img=...).
    288
    289    Cell contours are overlaid using the correct origin for each crop.
    290    Detected object contours can optionally be shown on the aligned panel.
    291    """
    292
    293    cell = next(
    294        (item for item in image_object.cells if item.cell_id == cell_id),
    295        None,
    296    )
    297    if cell is None:
    298        raise ValueError(f"Cell {cell_id} was not found")
    299
    300    if not image_object.channels or channel not in image_object.channels:
    301        raise ValueError(f"Channel {channel} is not loaded")
    302
    303    contour = np.asarray(cell.contour)
    304    label = channel_label or channel
    305
    306    # BEFORE: use the original image and the same cropping method
    307    # as plot_cell_gallery().
    308    extent = int(np.ceil(np.max(np.ptp(contour, axis=0)))) + 24
    309
    310    crop = _tutorial_crop(
    311        image_object.image.shape,
    312        contour.mean(axis=0),
    313        max(window_px, extent),
    314    )
    315
    316    signal_before = image_object.channels[channel][crop]
    317
    318    # Crop origin in image coordinates: (column, row).
    319    origin_before = (
    320        crop[1].start,
    321        crop[0].start,
    322    )
    323
    324    # AFTER: apply the same local registration as in plot_cell_gallery().
    325    signal_after, _, _, row_offset, column_offset = u.crop_image(
    326        image_object.channels[channel],
    327        cell.contour,
    328        phase_img=image_object.image,
    329    )
    330
    331    # The function returns row offset first, column offset second.
    332    # For plotting, x = column and y = row.
    333    origin_after = (
    334        column_offset,
    335        row_offset,
    336    )
    337
    338    panels = [
    339        (signal_before, origin_before, f"01  UNALIGNED {channel}"),
    340        (signal_after, origin_after, f"02  ALIGNED {channel}"),
    341    ]
    342
    343    fig, axes = plt.subplots(
    344        1, 2,
    345        figsize=(10, 4.8),
    346    )
    347
    348    # Use identical contrast limits in both panels.
    349    combined = np.concatenate([
    350        signal_before.ravel(),
    351        signal_after.ravel(),
    352    ])
    353    combined = combined[np.isfinite(combined)]
    354
    355    if combined.size == 0:
    356        raise ValueError("Both signal crops contain no finite values")
    357
    358    low, high = np.percentile(combined, [1, 99.5])
    359
    360    if low == high:
    361        high = low + 1
    362
    363    for ax, (signal, origin, title) in zip(axes, panels):
    364
    365        ax.imshow(
    366            signal,
    367            cmap=_tutorial_channel_cmap(label),
    368            vmin=low,
    369            vmax=high,
    370        )
    371
    372        ax.set_title(
    373            title,
    374            loc="left",
    375            fontsize=11,
    376            fontweight="bold",
    377            color="#193348",
    378        )
    379
    380        ax.set_xticks([])
    381        ax.set_yticks([])
    382
    383        for spine in ax.spines.values():
    384            spine.set_visible(False)
    385
    386        # Cell contours use (row, column).
    387        # Matplotlib uses (x=column, y=row).
    388        x0, y0 = origin
    389
    390        ax.plot(
    391            contour[:, 1] - x0,
    392            contour[:, 0] - y0,
    393            color="#FFCC66",
    394            lw=1.6,
    395        )
    396
    397    # Show detected objects only on the aligned panel.
    398    # These contours are taken directly from the cell object.
    399    if show_objects:
    400        x0, y0 = origin_after
    401
    402        for detected in cell.object_meshdata.get(
    403            channel, {}
    404        ).get("object_contour", []):
    405
    406            detected = np.asarray(detected)
    407
    408            if detected.ndim == 2 and detected.shape[1] == 2:
    409                axes[1].plot(
    410                    detected[:, 1] - x0,
    411                    detected[:, 0] - y0,
    412                    color="#4DE6F2",
    413                    lw=1.2,
    414                )
    415
    416    fig.suptitle(
    417        f"{channel} ALIGNMENT  /  CELL {cell_id}  /  "
    418        f"{image_object.image_name}",
    419        x=0.02,
    420        y=0.98,
    421        ha="left",
    422        fontsize=13,
    423        fontweight="bold",
    424        color="#193348",
    425    )
    426
    427    fig.tight_layout(rect=(0, 0, 1, 0.93))
    428
    429    return fig
    430def snapshot_cell_contours(image_objects):
    431    """Copy cell IDs and contours before curation changes the collection."""
    432    return {
    433        image.image_name: {
    434            cell.cell_id: np.asarray(cell.contour).copy() for cell in image.cells
    435        }
    436        for image in image_objects
    437    }
    438
    439
    440def plot_curation_review(before_contours, image_objects, crop_size=700):
    441    """Show retained cells in green and removed cells in red on phase images.
    442
    443    ``before_contours`` is made by :func:`snapshot_cell_contours` before
    444    curation. Cells are matched by image name and cell ID.
    445    """
    446    objects = list(image_objects)
    447    if not objects:
    448        raise ValueError("image_objects is empty")
    449    fig, axes = plt.subplots(1, len(objects), figsize=(5.2 * len(objects), 5.2), squeeze=False)
    450    for ax, image in zip(axes[0], objects):
    451        original = before_contours.get(image.image_name)
    452        if original is None:
    453            raise ValueError(f"No pre-curation contours for {image.image_name}")
    454        kept = {cell.cell_id for cell in image.cells}
    455        if not kept.issubset(original):
    456            raise ValueError(f"Curated cell IDs do not match the original meshes for {image.image_name}")
    457        for cell in image.cells:
    458            old_center = np.asarray(original[cell.cell_id]).mean(axis=0)
    459            new_center = np.asarray(cell.contour).mean(axis=0)
    460            if np.linalg.norm(old_center - new_center) > 1:
    461                raise ValueError(f"Curated cell geometry does not match the original mesh for {image.image_name}")
    462        crop = _tutorial_crop(
    463            image.image.shape,
    464            (image.image.shape[0] // 2, image.image.shape[1] // 2),
    465            crop_size,
    466        )
    467        phase = image.image[crop]
    468        low, high = _tutorial_contrast(phase)
    469        ax.imshow(phase, cmap="gray", vmin=low, vmax=high)
    470        x0, y0 = crop[1].start, crop[0].start
    471        for cell_id, contour in original.items():
    472            contour = np.asarray(contour)
    473            margin = 3
    474            fully_inside = (
    475                (contour[:, 0] >= crop[0].start + margin)
    476                & (contour[:, 0] < crop[0].stop - margin)
    477                & (contour[:, 1] >= crop[1].start + margin)
    478                & (contour[:, 1] < crop[1].stop - margin)
    479            )
    480            if not fully_inside.all():
    481                continue
    482            keep = cell_id in kept
    483            ax.plot(
    484                contour[:, 1] - x0, contour[:, 0] - y0,
    485                color="#34D986" if keep else "#FF686B",
    486                lw=1.15 if keep else 0.85,
    487                alpha=0.95 if keep else 0.82,
    488            )
    489        field = image.image_name.rsplit("_", 1)[0]
    490        ax.set_title(
    491            f"{field}\n{len(kept):,} kept  ·  {len(original) - len(kept):,} removed",
    492            loc="left", fontsize=11, color="#193348",
    493        )
    494        ax.set_xticks([])
    495        ax.set_yticks([])
    496        for spine in ax.spines.values():
    497            spine.set_visible(False)
    498    fig.legend(
    499        handles=[
    500            Line2D([], [], color="#34D986", lw=2, label="Kept"),
    501            Line2D([], [], color="#FF686B", lw=2, label="Removed"),
    502        ],
    503        loc="lower center", ncol=2, frameon=False, bbox_to_anchor=(0.5, -0.01),
    504    )
    505    fig.suptitle("CURATION REVIEW", x=0.02, y=1.04, ha="left", fontsize=15, fontweight="bold", color="#193348")
    506    fig.tight_layout(rect=(0, 0.035, 1, 1))
    507    return fig
    508
    509
    510def plot_signal_gallery(image_object, cell_id, channel_labels=None,
    511                        window_px=120, show_objects=False):
    512    """Show the phase image beside raw channel crops, with optional detections.
    513
    514    Cell contours are displayed on all panels. Detected object contours
    515    are optionally overlaid on the corresponding fluorescence channels.
    516
    517    Crops share pixel coordinates. This function does not register channels.
    518    """
    519    labels = dict(channel_labels or {"C4": "RNA", "C5": "DAPI / nucleoid"})
    520    if not labels:
    521        raise ValueError("channel_labels is empty")
    522
    523    cell = next(
    524        (item for item in image_object.cells if item.cell_id == cell_id),
    525        None,
    526    )
    527    if cell is None:
    528        raise ValueError(f"Cell {cell_id} was not found")
    529
    530    missing = [
    531        channel for channel in labels
    532        if not image_object.channels or channel not in image_object.channels
    533    ]
    534    if missing:
    535        raise ValueError(f"Channels not loaded: {missing}")
    536
    537    contour = np.asarray(cell.contour)
    538    crop = _tutorial_crop(
    539        image_object.image.shape, contour.mean(axis=0), window_px
    540    )
    541    x0, y0 = crop[1].start, crop[0].start
    542
    543    # Prepare panels
    544    panels = [("C1  PHASE", image_object.image, "gray", None)]
    545
    546    for channel, label in labels.items():
    547        cmap = _tutorial_channel_cmap(label)
    548        panels.append((
    549            f"{channel}  {label.upper()}",
    550            image_object.channels[channel],
    551            cmap,
    552            channel,
    553        ))
    554
    555    fig, axes = plt.subplots(
    556        1, len(panels),
    557        figsize=(3.7 * len(panels), 3.9),
    558    )
    559
    560    for ax, (title, full_image, cmap, channel) in zip(axes, panels):
    561        data = full_image[crop]
    562        low, high = _tutorial_contrast(data)
    563
    564        ax.imshow(data, cmap=cmap, vmin=low, vmax=high)
    565        ax.set_title(
    566            title, loc="left", fontsize=10,
    567            fontweight="bold", color="#193348",
    568        )
    569        ax.set_xticks([])
    570        ax.set_yticks([])
    571
    572        for spine in ax.spines.values():
    573            spine.set_visible(False)
    574
    575        # Plot cell contour on every panel
    576        ax.plot(
    577            contour[:, 1] - x0,
    578            contour[:, 0] - y0,
    579            color="#FFCF70",
    580            lw=1.5,
    581        )
    582
    583        # Optionally plot detected objects on fluorescence panels
    584        if channel is not None and show_objects:
    585            contours = cell.object_meshdata.get(
    586                channel, {}
    587            ).get("object_contour", [])
    588
    589            if contours is not None:
    590                for item in contours:
    591                    outline = np.asarray(item)
    592
    593                    if outline.ndim == 2 and outline.shape[1] == 2:
    594                        ax.plot(
    595                            outline[:, 1] - x0,
    596                            outline[:, 0] - y0,
    597                            color="#4DE6F2",
    598                            lw=1.2,
    599                        )
    600
    601    fig.suptitle(
    602        f"CHANNEL EXAMPLE  /  CELL {cell_id}",
    603        x=0.02, y=0.98, ha="left",
    604        fontsize=14, fontweight="bold", color="#193348",
    605    )
    606
    607    fig.tight_layout(rect=(0, 0, 1, 0.93))
    608    return fig
    609
    610
    611def plot_object_area_vs_cell_area(features, channel_labels):
    612    """Compare detected object area with cell area for each named channel."""
    613    labels = dict(channel_labels)
    614    if not labels:
    615        raise ValueError("channel_labels is empty")
    616    if "cell_area" not in features:
    617        raise ValueError("Missing cell_area")
    618    fig, axes = plt.subplots(1, len(labels), figsize=(4.3 * len(labels), 4), squeeze=False)
    619    x = np.asarray(features["cell_area"], dtype=float)
    620    for ax, (channel, label) in zip(axes[0], labels.items()):
    621        column = f"{channel}_cell_total_obj_area"
    622        if column not in features:
    623            raise ValueError(f"Missing {column}; calculate object features first")
    624        y = np.asarray(features[column], dtype=float)
    625        valid = np.isfinite(x) & np.isfinite(y)
    626        ax.scatter(x[valid], y[valid], s=12, color=_tutorial_channel_color(label), alpha=0.42, linewidths=0)
    627        ax.set_title(f"{channel}  {label}", loc="left", fontsize=11, fontweight="bold", color="#193348")
    628        ax.set_xlabel("Cell area (µm²)")
    629        ax.set_ylabel("Detected object area (µm²)")
    630        ax.text(0.97, 0.96, f"n = {valid.sum():,}", transform=ax.transAxes, ha="right", va="top", fontsize=9, color="#516676")
    631        ax.spines[["top", "right"]].set_visible(False)
    632        ax.grid(color="#E7EDF1", lw=0.8)
    633        ax.set_axisbelow(True)
    634        ax.tick_params(colors="#516676", labelsize=9)
    635    fig.suptitle("SUBCELLULAR AREA  /  CELL AREA", x=0.02, y=1.04, ha="left", fontsize=14, fontweight="bold", color="#193348")
    636    fig.tight_layout()
    637    return fig
    638
    639
    640def plot_feature_overview(features, signal_column="C4_NC_ratio"):
    641    """
    642    Exploratory overview of cell morphology, constriction,
    643    detected objects, and relative signal localization.
    644
    645    C2 = outer membrane
    646    C3 = inner membrane
    647    C4 = RNA
    648    C5 = DAPI
    649
    650    Missing features are indicated in their panels.
    651    Nothing is saved.
    652    """
    653
    654    def numeric(column):
    655        if column not in features.columns:
    656            return None
    657        return pd.to_numeric(features[column], errors="coerce").to_numpy()
    658
    659    def unavailable(ax, *columns):
    660        ax.text(
    661            0.5, 0.5,
    662            "Feature unavailable:\n" + "\n".join(columns),
    663            ha="center", va="center",
    664            transform=ax.transAxes,
    665            fontsize=9, color="gray",
    666        )
    667        ax.set_axis_off()
    668
    669    length = numeric("cell_length")
    670    width = numeric("cell_width")
    671    constriction = numeric("cell_constriction_degree_width_axial")
    672
    673    fig, axes = plt.subplots(
    674        2, 3,
    675        figsize=(15, 9),
    676        constrained_layout=True,
    677    )
    678
    679    axes = axes.ravel()
    680
    681    for ax in axes:
    682        ax.spines[["top", "right"]].set_visible(False)
    683        ax.grid(axis="y", color="#E7EDF1", lw=0.8)
    684        ax.set_axisbelow(True)
    685        ax.tick_params(labelsize=9)
    686
    687    # 1. Cell morphology, coloured by constriction
    688    ax = axes[0]
    689
    690    if length is not None and width is not None:
    691        valid = np.isfinite(length) & np.isfinite(width)
    692
    693        if constriction is not None:
    694            valid &= np.isfinite(constriction)
    695
    696            scatter = ax.scatter(
    697                length[valid],
    698                width[valid],
    699                c=constriction[valid],
    700                cmap="viridis",
    701                s=10,
    702                alpha=0.55,
    703                linewidths=0,
    704            )
    705
    706            fig.colorbar(
    707                scatter,
    708                ax=ax,
    709                label="Constriction degree",
    710            )
    711        else:
    712            ax.scatter(
    713                length[valid],
    714                width[valid],
    715                color="#197E89",
    716                s=10,
    717                alpha=0.4,
    718                linewidths=0,
    719            )
    720
    721        ax.set(
    722            xlabel="Cell length (µm)",
    723            ylabel="Cell width (µm)",
    724            title="Cell shape and constriction",
    725        )
    726    else:
    727        unavailable(ax, "cell_length", "cell_width")
    728
    729    # 2. Constriction-degree distribution
    730    ax = axes[1]
    731
    732    if constriction is not None:
    733        valid = constriction[np.isfinite(constriction)]
    734
    735        ax.hist(
    736            valid,
    737            bins=35,
    738            color="#197E89",
    739            edgecolor="white",
    740        )
    741
    742        ax.set(
    743            xlabel="Constriction degree",
    744            ylabel="Cells",
    745            title="Cell constriction",
    746        )
    747    else:
    748        unavailable(ax, "cell_constriction_degree_width_axial")
    749
    750    # 3. Number of detected DAPI objects per cell
    751    ax = axes[2]
    752    object_number = numeric("C5_object_number")
    753
    754    if object_number is not None:
    755        counts = (
    756            pd.Series(object_number)
    757            .dropna()
    758            .value_counts()
    759            .sort_index()
    760        )
    761
    762        ax.bar(
    763            counts.index,
    764            counts.values,
    765            color="#7F3F98",
    766            width=0.75,
    767        )
    768
    769        ax.set(
    770            xlabel="Detected DAPI objects per cell",
    771            ylabel="Cells",
    772            title="DAPI object distribution",
    773        )
    774    else:
    775        unavailable(ax, "C5_object_number")
    776
    777    # 4. Nucleoid-to-cell ratio
    778    ax = axes[3]
    779    nc_ratio = numeric("C5_NC_ratio")
    780
    781    if nc_ratio is not None:
    782        valid = nc_ratio[np.isfinite(nc_ratio)]
    783
    784        ax.hist(
    785            valid,
    786            bins=35,
    787            color="#7F3F98",
    788            edgecolor="white",
    789        )
    790
    791        ax.set(
    792            xlabel="DAPI NC ratio",
    793            ylabel="Cells",
    794            title="Nucleoid occupancy",
    795        )
    796    else:
    797        unavailable(ax, "C5_NC_ratio")
    798
    799    # 5. Relative object occupancy versus cell length
    800    ax = axes[4]
    801    signal = numeric(signal_column)
    802
    803    if length is not None and signal is not None:
    804        valid = np.isfinite(length) & np.isfinite(signal)
    805
    806        ax.scatter(
    807            length[valid],
    808            signal[valid],
    809            color="#E17953",
    810            s=10,
    811            alpha=0.4,
    812            linewidths=0,
    813        )
    814
    815        ax.set(
    816            xlabel="Cell length (µm)",
    817            ylabel=signal_column.replace("_", " "),
    818            title="Object occupancy versus cell length",
    819        )
    820    else:
    821        unavailable(ax, "cell_length", signal_column)
    822
    823    # 6. Outer- versus inner-membrane signal constriction
    824    ax = axes[5]
    825
    826    outer = numeric("C2_signal_constriction_degree_mesh")
    827    inner = numeric("C3_signal_constriction_degree_mesh")
    828
    829    if outer is not None and inner is not None:
    830        valid = np.isfinite(outer) & np.isfinite(inner)
    831
    832        ax.scatter(
    833            outer[valid],
    834            inner[valid],
    835            color="#36A5A0",
    836            s=10,
    837            alpha=0.4,
    838            linewidths=0,
    839        )
    840
    841        ax.set(
    842            xlabel="Outer membrane constriction degree (C2)",
    843            ylabel="Inner membrane constriction degree (C3)",
    844            title="Membrane signal localization",
    845        )
    846    else:
    847        unavailable(
    848            ax,
    849            "C2_signal_constriction_degree_mesh",
    850            "C3_signal_constriction_degree_mesh",
    851        )
    852
    853    fig.suptitle(
    854        f"CELL FEATURE OVERVIEW  /  {len(features):,} cells",
    855        fontsize=15,
    856        fontweight="bold",
    857        color="#193348",
    858    )
    859
    860    return fig
    861import numpy as np
    862import pandas as pd
    863import matplotlib.pyplot as plt
    864
    865
    866def plot_feature_histogram(
    867    df,
    868    feature,
    869    bins=40,
    870    xlabel=None,
    871    color="#197E89",
    872    figsize=(7, 4.5),
    873):
    874    """Plot a feature distribution with mean, median and interquartile range."""
    875
    876    if feature not in df.columns:
    877        raise ValueError(f"Feature not found: {feature}")
    878
    879    values = pd.to_numeric(df[feature], errors="coerce").to_numpy()
    880    values = values[np.isfinite(values)]
    881
    882    if len(values) == 0:
    883        raise ValueError(f"No valid values for {feature}")
    884
    885    mean = np.mean(values)
    886    median = np.median(values)
    887    q1, q3 = np.percentile(values, [25, 75])
    888
    889    fig, ax = plt.subplots(figsize=figsize)
    890
    891    # Histogram
    892    ax.hist(
    893        values,
    894        bins=bins,
    895        color=color,
    896        alpha=0.85,
    897        edgecolor="white",
    898        linewidth=0.6,
    899    )
    900
    901    # Interquartile range
    902    ax.axvspan(
    903        q1, q3,
    904        color=color,
    905        alpha=0.12,
    906        label=f"IQR: {q1:.2f}–{q3:.2f}",
    907    )
    908
    909    # Mean and median
    910    ax.axvline(
    911        mean,
    912        color="#E17953",
    913        linestyle="--",
    914        linewidth=1.8,
    915        label=f"Mean: {mean:.2f}",
    916    )
    917
    918    ax.axvline(
    919        median,
    920        color="#193348",
    921        linestyle="-",
    922        linewidth=1.8,
    923        label=f"Median: {median:.2f}",
    924    )
    925
    926    ax.set_xlabel(xlabel or feature.replace("_", " "), fontsize=11)
    927    ax.set_ylabel("Number of cells", fontsize=11)
    928    ax.set_title(
    929        f"{feature.replace('_', ' ')}\nN = {len(values):,} cells",
    930        fontsize=13,
    931        fontweight="bold",
    932        loc="left",
    933    )
    934
    935    ax.legend(frameon=False, fontsize=9)
    936
    937    ax.spines[["top", "right"]].set_visible(False)
    938    ax.grid(axis="y", alpha=0.15)
    939    ax.set_axisbelow(True)
    940
    941    fig.tight_layout()
    942    return fig
    943
    944
    945def plot_feature_scatter(
    946    df,
    947    x_feature,
    948    y_feature,
    949    color_by=None,
    950    xlabel=None,
    951    ylabel=None,
    952    trendline=True,
    953    figsize=(7, 5),
    954):
    955    """
    956    Plot two features against each other.
    957
    958    Optionally colour points by a third numeric feature.
    959    A linear regression line and Pearson correlation can be shown.
    960    """
    961
    962    columns = [x_feature, y_feature]
    963
    964    if color_by is not None:
    965        columns.append(color_by)
    966
    967    missing = [col for col in columns if col not in df.columns]
    968    if missing:
    969        raise ValueError(f"Features not found: {missing}")
    970
    971    # Convert selected features to numeric values
    972    data = df[columns].apply(pd.to_numeric, errors="coerce")
    973    data = data.replace([np.inf, -np.inf], np.nan).dropna()
    974
    975    if data.empty:
    976        raise ValueError("No valid cells remain for the selected features")
    977
    978    x = data[x_feature].to_numpy()
    979    y = data[y_feature].to_numpy()
    980
    981    fig, ax = plt.subplots(figsize=figsize)
    982
    983    # Scatter plot
    984    if color_by is None:
    985        ax.scatter(
    986            x, y,
    987            s=13,
    988            alpha=0.3,
    989            color="#197E89",
    990            edgecolors="none",
    991            rasterized=True,
    992        )
    993    else:
    994        scatter = ax.scatter(
    995            x, y,
    996            c=data[color_by],
    997            cmap="viridis",
    998            s=13,
    999            alpha=0.5,
   1000            edgecolors="none",
   1001            rasterized=True,
   1002        )
   1003
   1004        fig.colorbar(
   1005            scatter,
   1006            ax=ax,
   1007            label=color_by.replace("_", " "),
   1008        )
   1009
   1010    # Linear trend and Pearson correlation
   1011    if len(x) >= 2 and np.std(x) > 0 and np.std(y) > 0:
   1012
   1013        r = np.corrcoef(x, y)[0, 1]
   1014
   1015        if trendline:
   1016            slope, intercept = np.polyfit(x, y, 1)
   1017
   1018            x_line = np.linspace(x.min(), x.max(), 200)
   1019
   1020            ax.plot(
   1021                x_line,
   1022                slope * x_line + intercept,
   1023                color="#E17953",
   1024                linewidth=2,
   1025                linestyle="--",
   1026                label=f"Linear fit (r = {r:.2f})",
   1027            )
   1028
   1029            ax.legend(frameon=False)
   1030
   1031        else:
   1032            ax.text(
   1033                0.03, 0.97,
   1034                f"Pearson r = {r:.2f}",
   1035                transform=ax.transAxes,
   1036                ha="left",
   1037                va="top",
   1038                fontsize=10,
   1039            )
   1040
   1041    ax.set_xlabel(xlabel or x_feature.replace("_", " "), fontsize=11)
   1042    ax.set_ylabel(ylabel or y_feature.replace("_", " "), fontsize=11)
   1043
   1044    ax.set_title(
   1045        f"{x_feature.replace('_', ' ')} vs. {y_feature.replace('_', ' ')}"
   1046        f"\nN = {len(data):,} cells",
   1047        fontsize=13,
   1048        fontweight="bold",
   1049        loc="left",
   1050    )
   1051
   1052    ax.spines[["top", "right"]].set_visible(False)
   1053    ax.grid(alpha=0.15)
   1054    ax.set_axisbelow(True)
   1055
   1056    fig.tight_layout()
   1057    return fig
   1058
   1059def plot_contour(
   1060    cell_id,
   1061    image,
   1062    base_channel=None,
   1063    verbose=False,
   1064    add_scalebar=True,
   1065    scalebar_um=1,
   1066    um_per_px=0.065,
   1067    scalebar_color="white",
   1068    scalebar_lw=6,
   1069    window_px=None,        # <-- NEW: fixed window size (square)
   1070    save_svg_path=None,
   1071):
   1072    """
   1073    Plot the contour of a cell with optional scalebar.
   1074
   1075    Parameters
   1076    ----------
   1077    window_px : int or None
   1078        If provided, creates a fixed square window (in pixels) centered on the cell.
   1079        If None, tight crop is used (default behavior).
   1080    """
   1081
   1082    cellobj = next((cell for cell in image.cells if cell.cell_id == cell_id), None)
   1083
   1084    if cellobj is None:
   1085        if verbose:
   1086            print(f"Cell with ID {cell_id} not found.")
   1087        return
   1088
   1089    contour = cellobj.contour
   1090    full_img = image.image
   1091
   1092    # ---------------------------
   1093    # Fixed window crop
   1094    # ---------------------------
   1095    if window_px is not None:
   1096
   1097        # center of cell
   1098        cy = int(np.mean(contour[:, 0]))
   1099        cx = int(np.mean(contour[:, 1]))
   1100
   1101        half = window_px // 2
   1102
   1103        y_min = max(cy - half, 0)
   1104        y_max = min(cy + half, full_img.shape[0])
   1105        x_min = max(cx - half, 0)
   1106        x_max = min(cx + half, full_img.shape[1])
   1107
   1108        cropped_img = full_img[y_min:y_max, x_min:x_max]
   1109
   1110        adjusted_contour_x = contour[:, 1] - x_min
   1111        adjusted_contour_y = contour[:, 0] - y_min
   1112
   1113    # ---------------------------
   1114    # Default tight crop
   1115    # ---------------------------
   1116    else:
   1117        cropped_img, _, _, x_offset, y_offset = u.crop_image(
   1118            image=full_img, contour=contour, mask_to_crop=None, phase_img=None
   1119        )
   1120
   1121        adjusted_contour_x = contour[:, 1] - y_offset
   1122        adjusted_contour_y = contour[:, 0] - x_offset
   1123
   1124    # ---------------------------
   1125    # SVG settings (editable)
   1126    # ---------------------------
   1127    plt.rcParams.update({
   1128        "font.family": "Arial",
   1129        "svg.fonttype": "none",
   1130    })
   1131
   1132    fig, ax = plt.subplots(figsize=(6, 6))
   1133    ax.imshow(cropped_img, cmap="gist_gray")
   1134
   1135    ax.plot(adjusted_contour_x, adjusted_contour_y, "-", c="w", lw=3)
   1136
   1137    # ---------------------------
   1138    # Scalebar (vector)
   1139    # ---------------------------
   1140    if add_scalebar:
   1141        bar_px = scalebar_um / um_per_px
   1142        h, w = cropped_img.shape[:2]
   1143
   1144        x_start = int(w * 0.05)
   1145        y_start = int(h * 0.95)
   1146
   1147        ax.plot(
   1148            [x_start, x_start + bar_px],
   1149            [y_start, y_start],
   1150            color=scalebar_color,
   1151            lw=scalebar_lw,
   1152            solid_capstyle="butt",
   1153        )
   1154
   1155        ax.text(
   1156            x_start + bar_px / 2,
   1157            y_start - h * 0.05,
   1158            f"{scalebar_um:.1f} µm",
   1159            color=scalebar_color,
   1160            ha="center",
   1161            va="bottom",
   1162            fontsize=12,
   1163            fontfamily="Arial",
   1164        )
   1165
   1166    ax.set_title(
   1167        f"Cell {cellobj.cell_id} Frame {image.frame}",
   1168        fontfamily="Arial"
   1169    )
   1170
   1171    ax.axis("off")
   1172    fig.tight_layout(pad=0)
   1173
   1174    if save_svg_path:
   1175        fig.savefig(
   1176            save_svg_path,
   1177            format="svg",
   1178            bbox_inches="tight",
   1179            pad_inches=0
   1180        )
   1181
   1182    plt.show()
   1183    plt.close(fig)
   1184
   1185
   1186
   1187    # plt.plot(midline.T[1], midline.T[0], '-', c='blue', lw=5)
   1188    # plt.plot([msh[:,1][::2], msh[:,3][::2]], [msh[:,0][::2], msh[:,2][::2]], '-', c='cyan', lw=2)
   1189import os
   1190import numpy as np
   1191import matplotlib.pyplot as plt
   1192
   1193def plot_cell_contour_gallery(
   1194    cell_id,
   1195    image,
   1196    channels,
   1197    window_px=100,
   1198    phase_first=True,
   1199    phase_cmap="gist_grey",
   1200    channel_cmap="gist_grey",
   1201    contour_color="yellow",
   1202    contour_lw=1.5,
   1203    add_scalebar=True,
   1204    scalebar_um=1,
   1205    um_per_px=0.065,
   1206    scalebar_color="white",
   1207    scalebar_lw=4,
   1208    title=True,
   1209    figsize_per_panel=3.0,
   1210    save_path=None,
   1211    save_individual_dir=None,
   1212    dpi=500,
   1213    verbose=False,
   1214):
   1215    """
   1216    Plot a cropped gallery for one cell using:
   1217      - image.image for phase contrast
   1218      - image.channels[ch] for channel images
   1219      - only the cell contour overlay
   1220
   1221    Parameters
   1222    ----------
   1223    cell_id : int or str
   1224    image : object
   1225        Typically ic.image_objects[frame_id]
   1226    channels : list of str
   1227        Example: ["C1", "C2", "C3", "C4", "C5"]
   1228    """
   1229
   1230    cellobj = next((cell for cell in image.cells if cell.cell_id == cell_id), None)
   1231    if cellobj is None:
   1232        if verbose:
   1233            print(f"Cell with ID {cell_id} not found.")
   1234        return None
   1235
   1236    contour = np.asarray(cellobj.contour)   # expected shape: [N, 2] with [y, x]
   1237    full_phase = image.image                # phase contrast image
   1238
   1239    if full_phase is None:
   1240        if verbose:
   1241            print("Phase image not found at image.image")
   1242        return None
   1243
   1244    # fixed square crop centered on the cell
   1245    cy = int(np.mean(contour[:, 0]))
   1246    cx = int(np.mean(contour[:, 1]))
   1247    half = int(window_px // 2)
   1248
   1249    y_min = max(cy - half, 0)
   1250    y_max = min(cy + half, full_phase.shape[0])
   1251    x_min = max(cx - half, 0)
   1252    x_max = min(cx + half, full_phase.shape[1])
   1253
   1254    adjusted_contour_x = contour[:, 1] - x_min
   1255    adjusted_contour_y = contour[:, 0] - y_min
   1256
   1257    panels = []
   1258
   1259    if phase_first:
   1260        panels.append(("Phase", full_phase[y_min:y_max, x_min:x_max], phase_cmap))
   1261
   1262    for ch in channels:
   1263        if ch not in image.channels:
   1264            if verbose:
   1265                print(f"Channel {ch} not found in image.channels, skipping.")
   1266            continue
   1267        ch_img = image.channels[ch]
   1268        cropped = ch_img[y_min:y_max, x_min:x_max]
   1269        panels.append((str(ch), cropped, channel_cmap))
   1270
   1271    if len(panels) == 0:
   1272        if verbose:
   1273            print("No panels available to plot.")
   1274        return None
   1275
   1276    plt.rcParams.update({
   1277        "font.family": "Arial",
   1278        "svg.fonttype": "none",
   1279    })
   1280
   1281    n_panels = len(panels)
   1282    fig, axes = plt.subplots(
   1283        1, n_panels,
   1284        figsize=(figsize_per_panel * n_panels, figsize_per_panel),
   1285        squeeze=False
   1286    )
   1287    axes = axes.ravel()
   1288
   1289    def _add_scalebar(ax, img_shape):
   1290        if not add_scalebar:
   1291            return
   1292        h, w = img_shape[:2]
   1293        bar_px = scalebar_um / um_per_px
   1294        x_start = int(w * 0.05)
   1295        y_start = int(h * 0.93)
   1296
   1297        ax.plot(
   1298            [x_start, x_start + bar_px],
   1299            [y_start, y_start],
   1300            color=scalebar_color,
   1301            lw=scalebar_lw,
   1302            solid_capstyle="butt",
   1303        )
   1304        ax.text(
   1305            x_start + bar_px / 2,
   1306            y_start - h * 0.04,
   1307            f"{scalebar_um:.1f} µm",
   1308            color=scalebar_color,
   1309            ha="center",
   1310            va="bottom",
   1311            fontsize=10,
   1312            fontfamily="Arial",
   1313        )
   1314
   1315    for ax, (panel_name, cropped_img, cmap) in zip(axes, panels):
   1316        ax.imshow(cropped_img, cmap=cmap)
   1317        ax.plot(adjusted_contour_x, adjusted_contour_y, color=contour_color, lw=contour_lw)
   1318        _add_scalebar(ax, cropped_img.shape)
   1319
   1320        if title:
   1321            ax.set_title(panel_name)
   1322
   1323        ax.axis("off")
   1324
   1325    fig.tight_layout(pad=0.05)
   1326
   1327    if save_path is not None:
   1328        os.makedirs(os.path.dirname(save_path), exist_ok=True)
   1329        fig.savefig(save_path, bbox_inches="tight", pad_inches=0.02, dpi=dpi)
   1330
   1331    if save_individual_dir is not None:
   1332        os.makedirs(save_individual_dir, exist_ok=True)
   1333
   1334        for panel_name, cropped_img, cmap in panels:
   1335            fig_i, ax_i = plt.subplots(figsize=(figsize_per_panel, figsize_per_panel))
   1336            ax_i.imshow(cropped_img, cmap=cmap)
   1337            ax_i.plot(adjusted_contour_x, adjusted_contour_y, color=contour_color, lw=contour_lw)
   1338            _add_scalebar(ax_i, cropped_img.shape)
   1339
   1340            if title:
   1341                ax_i.set_title(panel_name)
   1342
   1343            ax_i.axis("off")
   1344            fig_i.tight_layout(pad=0.05)
   1345
   1346            safe_name = "".join(c if c.isalnum() or c in ("-", "_") else "_" for c in str(panel_name))
   1347            out_path = os.path.join(save_individual_dir, f"cell_{cell_id}_{safe_name}.png")
   1348            fig_i.savefig(out_path, bbox_inches="tight", pad_inches=0.02, dpi=dpi)
   1349            plt.close(fig_i)
   1350
   1351    plt.show()
   1352
   1353    return {
   1354        "cell_id": cell_id,
   1355        "crop_bounds": (y_min, y_max, x_min, x_max),
   1356        "panel_names": [p[0] for p in panels],
   1357        "figure": fig,
   1358    }
   1359
   1360def plot_aligned_contours(cell_id, image, channel, verbose=False):
   1361    """
   1362    Plot the contour of a cell on an image.
   1363
   1364    This function takes a cell object and an image object and plots the contour of the specified cell on the image.
   1365    It also displays the cell ID in the title of the plot.
   1366
   1367    Parameters
   1368    ----------
   1369    cell : Cell
   1370        The cell object to be plotted.
   1371    image : Image
   1372        The image object on which to plot the cell contour.
   1373    """
   1374    cellobj = next((cell for cell in image.cells if cell.cell_id == cell_id), None)
   1375    if cellobj is None:
   1376        if verbose:
   1377            print(f"Cell with ID {cell_id} not found.")
   1378        return
   1379    fig = plt.figure(figsize=(12, 12))
   1380
   1381    contour = cellobj.contour
   1382    midline = cellobj.midline
   1383    aligned_contour = cellobj.shifted_contour
   1384    xmin = np.min(contour.T[1] - 10)
   1385    xmax = np.max(contour.T[1] + 10)
   1386    ymin = np.min(contour.T[0] - 10)
   1387    ymax = np.max(contour.T[0] + 10)
   1388    plt.xlim(xmin, xmax)
   1389    plt.ylim(ymin, ymax)
   1390    plt.imshow(image.channels[channel], cmap="gist_gray")
   1391    plt.plot(contour.T[1], contour.T[0], "-", c="w", lw=5)
   1392    if aligned_contour is not None:
   1393        plt.plot(aligned_contour.T[1], aligned_contour.T[0], "-", c="orange", lw=5)
   1394    # plt.plot(midline.T[1], midline.T[0], '-', c='blue', lw=5)
   1395
   1396    plt.title(f"Cell {cellobj.cell_id} Frame {image.frame}")
   1397    plt.show()
   1398
   1399
   1400def plot_svm_controls(frame_cell_pairs, image_objects, message=None):
   1401    """
   1402    Plot cells based on frame and cell ID pairs from SVM control analysis.
   1403
   1404    This function takes frame and cell ID pairs, along with a list of image objects, and plots the cells from SVM control analysis.
   1405    It allows adding an optional message to the title of the plot.
   1406
   1407    Parameters
   1408    ----------
   1409    frame_cell_pairs : list of tuples
   1410        A list of (cell_id, frame) pairs for cells to be plotted.
   1411    image_objects : list
   1412        A list of image objects corresponding to the frames in frame_cell_pairs.
   1413    message : str, optional
   1414        An optional message to be added to the title of the plot.
   1415    """
   1416    images_by_frame = {}
   1417    for image_obj in image_objects:
   1418        if image_obj.frame in images_by_frame:
   1419            raise ValueError(f"Duplicate image frame {image_obj.frame}.")
   1420        images_by_frame[image_obj.frame] = image_obj
   1421
   1422    for cell_id, frame in frame_cell_pairs:
   1423        image_obj = images_by_frame.get(frame)
   1424        if image_obj is None:
   1425            raise ValueError(f"No image object has frame {frame}.")
   1426        cellobj = next(
   1427            (cell for cell in image_obj.cells if cell.cell_id == cell_id), None
   1428        )
   1429        if cellobj is None:
   1430            raise ValueError(f"Cell {cell_id} is absent from frame {frame}.")
   1431        contour = cellobj.contour
   1432        # ori_contour = image.mesh_dataframe['contour'][cell]
   1433        xmin = np.min(contour.T[1] - 10)
   1434        xmax = np.max(contour.T[1] + 10)
   1435        ymin = np.min(contour.T[0] - 10)
   1436        ymax = np.max(contour.T[0] + 10)
   1437
   1438        plt.figure(figsize=(12, 12))
   1439        plt.xlim(xmin, xmax)
   1440        plt.ylim(ymin, ymax)
   1441        plt.imshow(image_obj.image, cmap="gist_gray")
   1442        plt.plot(contour.T[1], contour.T[0], "-", c="w", lw=5)
   1443        if message is not None:
   1444            plt.title(f"Cell: {cellobj.cell_id} ---Frame: {frame} ---Type: {message}")
   1445            plt.show()
   1446
   1447
   1448def plot_mask(image_object, joined=False, dpi=400, alpha=0.4):
   1449    """
   1450    Plot the mask overlay on an image.
   1451
   1452    This function takes an image object and optionally a joined mask, and plots the mask overlay on the image.
   1453    The alpha parameter controls the transparency of the mask overlay.
   1454
   1455    Parameters
   1456    ----------
   1457    image_object : Image object
   1458        An image object containing the image and mask data.
   1459    joined : bool, optional
   1460        Whether to use the joined mask for overlay (default is False).
   1461    dpi : int, optional
   1462        Dots per inch for the plot (default is 400).
   1463    alpha : float, optional
   1464        Alpha (transparency) value for the mask overlay (default is 0.4).
   1465    """
   1466
   1467    plt.figure(
   1468        figsize=(image_object.image.shape[1] / dpi, image_object.image.shape[0] / dpi),
   1469        dpi=dpi,
   1470    )
   1471    plt.imshow(image_object.image, cmap="gray")
   1472
   1473    if joined is True:
   1474        plt.imshow(
   1475            label2rgb(label(image_object.joined_mask, connectivity=1), bg_label=0),
   1476            alpha=alpha,
   1477        )
   1478    else:
   1479
   1480        plt.imshow(
   1481            label2rgb(label(image_object.mask, connectivity=1), bg_label=0), alpha=alpha
   1482        )
   1483
   1484    plt.axis("off")  # Turn off the axis
   1485    plt.show()
   1486
   1487
   1488def plot_random_contours(image, num_cells_to_plot=4, crop_size=200, scalebar=False):
   1489    """
   1490    Plot random cell contours within an image.
   1491
   1492    This function selects and plots random cell contours within an image, showing a specified number of cells.
   1493
   1494    Parameters
   1495    ----------
   1496    image : Image object
   1497        An image object containing the image and cell contour data.
   1498    num_cells_to_plot : int, optional
   1499        Number of random cell contours to plot (default is 4).
   1500    crop_size : int, optional
   1501        Size of the cropped area around each cell (default is 200).
   1502    scalebar : bool, optional
   1503        Whether to add a scale bar to the plots (default is False).
   1504
   1505    Returns
   1506    -------
   1507    None
   1508    """
   1509
   1510    fig, axes = plt.subplots(2, 2, figsize=(10, 10))
   1511
   1512    plotted_cells = set()  # To keep track of which cells have been plotted
   1513
   1514    for i in range(num_cells_to_plot):
   1515        row = i // 2
   1516        col = i % 2
   1517
   1518        # Find a random cell that hasn't been plotted yet
   1519        while True:
   1520            cell_idx = random.randint(0, len(image.cells) - 1)
   1521            if cell_idx not in plotted_cells:
   1522                break
   1523
   1524        cellobj = image.cells[cell_idx]
   1525        contour = cellobj.contour
   1526        center_x = np.mean(contour.T[1])
   1527        center_y = np.mean(contour.T[0])
   1528        half_crop_size = crop_size / 2
   1529        xmin = center_x - half_crop_size
   1530        xmax = center_x + half_crop_size
   1531        ymin = center_y - half_crop_size
   1532        ymax = center_y + half_crop_size
   1533
   1534        axes[row, col].set_xlim(xmin, xmax)
   1535        axes[row, col].set_ylim(ymin, ymax)
   1536        axes[row, col].imshow(image.image, cmap="gist_gray")
   1537        axes[row, col].plot(contour.T[1], contour.T[0], "-", c="w", lw=2)
   1538        axes[row, col].set_title(f"Cell {cellobj.cell_id}")
   1539
   1540        plotted_cells.add(cell_idx)  # Mark this cell as plotted
   1541        if scalebar:
   1542            scale_length_um = 1
   1543            scale_length_px = scale_length_um / 0.065  # Convert µm to pixels
   1544            scale_x = (
   1545                xmin + (xmax - xmin) * 0.02
   1546            )  # Adjust the X position of the scale bar
   1547            scale_y = (
   1548                ymin + (ymax - ymin) * 0.02
   1549            )  # Adjust the Y position of the scale bar
   1550            axes[row, col].plot(
   1551                [scale_x, scale_x + scale_length_px],
   1552                [scale_y, scale_y],
   1553                color="white",
   1554                lw=5,
   1555            )
   1556
   1557    # Remove axes for all subplots
   1558    for ax in axes.flat:
   1559        ax.axis("off")
   1560
   1561    plt.tight_layout()
   1562    plt.show()
   1563
   1564
   1565def plot_random_objects(
   1566    image, chann, num_objects_to_plot=4, crop_size=300, scalebar=True
   1567):
   1568    num_rows = 2  # Number of rows in the grid
   1569    num_cols = 2  # Number of columns in the grid
   1570
   1571    fig, axes = plt.subplots(num_rows, num_cols, figsize=(10, 10))
   1572
   1573    plotted_objects = set()  # To keep track of which objects have been plotted
   1574
   1575    for i in range(num_objects_to_plot):
   1576        # Find a random cell that hasn't been plotted yet
   1577        while True:
   1578            object_idx = random.randint(0, len(image.cells) - 1)
   1579            if object_idx not in plotted_objects:
   1580                break
   1581
   1582        cellobj = image.cells[object_idx]
   1583        contour = cellobj.contour
   1584
   1585        # Calculate cropping boundaries based on the fixed crop size
   1586        center_x = np.mean(contour.T[1])
   1587        center_y = np.mean(contour.T[0])
   1588        half_crop_size = crop_size / 2
   1589        xmin = center_x - half_crop_size
   1590        xmax = center_x + half_crop_size
   1591        ymin = center_y - half_crop_size
   1592        ymax = center_y + half_crop_size
   1593
   1594        row = i // num_cols
   1595        col = i % num_cols
   1596
   1597        axes[row, col].set_xlim(xmin, xmax)
   1598        axes[row, col].set_ylim(ymin, ymax)
   1599        axes[row, col].imshow(image.channels[chann], cmap="gist_gray")
   1600        axes[row, col].plot(contour.T[1], contour.T[0], "-", c="yellow", lw=2)
   1601        if cellobj.object_meshdata[chann]["object_contour"] is not None:
   1602            for nuc_contour in cellobj.object_meshdata[chann]["object_contour"]:
   1603                axes[row, col].plot(
   1604                    nuc_contour.T[1], nuc_contour.T[0], "-", c="cyan", lw=2
   1605                )
   1606
   1607        axes[row, col].set_title(f"Cell {cellobj.cell_id}")
   1608        axes[row, col].axis("off")
   1609
   1610        # Add scale bar
   1611        if scalebar:
   1612            scale_length_um = 1
   1613            scale_length_px = scale_length_um / 0.065  # Convert µm to pixels
   1614            scale_x = (
   1615                xmin + 10
   1616            )  # Adjust the position of the scale bar (e.g., 10 pixels from the left)
   1617            scale_y = (
   1618                ymin + 10
   1619            )  # Adjust the position of the scale bar (e.g., 10 pixels from the bottom)
   1620            axes[row, col].plot(
   1621                [scale_x, scale_x + scale_length_px],
   1622                [scale_y, scale_y],
   1623                color="white",
   1624                lw=5,
   1625            )
   1626
   1627        plotted_objects.add(object_idx)  # Mark this object as plotted
   1628
   1629    plt.tight_layout()
   1630    plt.show()
   1631
   1632
   1633def plot_random_axial(feature_dataframe, channels, method, num_objects_to_plot=2):
   1634    num_plots = 2  # Number of vertical plots
   1635    num_objects_per_plot = num_objects_to_plot // num_plots  # Objects per plot
   1636    fig, axes = plt.subplots(num_plots, 1, figsize=(12, 10))  # Make the plot taller
   1637
   1638    # Get a list of unique cell IDs
   1639    unique_cell_ids = feature_dataframe["cell_id"].unique()
   1640
   1641    # Randomly select cell IDs to plot
   1642    random_cell_ids = random.sample(list(unique_cell_ids), num_objects_to_plot)
   1643
   1644    for plot_index in range(num_plots):
   1645        for i in range(num_objects_per_plot):
   1646            cell_id = random_cell_ids[plot_index * num_objects_per_plot + i]
   1647
   1648            for channel in channels:
   1649                axial_intensity_column = f"{channel}_axial_intensity"
   1650                cell_frame_condition = feature_dataframe["cell_id"] == cell_id
   1651
   1652                if feature_dataframe[cell_frame_condition].empty:
   1653                    continue
   1654
   1655                axial_intensity = feature_dataframe[axial_intensity_column][
   1656                    cell_frame_condition
   1657                ].iloc[0]
   1658
   1659                # Normalize the axial intensity curve by dividing by the maximum value
   1660                normalized_midline_intensity_data = (
   1661                    axial_intensity - np.min(axial_intensity)
   1662                ) / (np.max(axial_intensity) - np.min(axial_intensity))
   1663
   1664                color = f"C{channels.index(channel)}"  # Assign a unique color to each channel
   1665                ax = axes[plot_index]
   1666
   1667                ax.plot(
   1668                    normalized_midline_intensity_data,
   1669                    linestyle="-",
   1670                    color=color,
   1671                    label=f"Channel {channel}",
   1672                )
   1673
   1674            # Retrieve the corresponding frame for the cell
   1675            frame = feature_dataframe["frame"][cell_frame_condition].iloc[0]
   1676
   1677            ax = axes[plot_index]
   1678            ax.set_title(f"Cell {cell_id}, Frame {frame}")
   1679            ax.set_xlabel("Position along Axial")
   1680            ax.set_ylabel("Normalized Axial Intensity")
   1681            ax.set_ylim(0, 1.2)
   1682            ax.legend(frameon=False)  # Remove the frame around the legend
   1683
   1684    plt.tight_layout()
   1685    plt.show()
   1686
   1687
   1688def plot_normalized_axial_intensity(
   1689    feature_dataframe, channels, selected_frame=None, selected_cell_id=None,
   1690    show=True,
   1691):
   1692    """Plot normalized axial profiles for selected cells."""
   1693    profile_columns = [f"{channel}_axial_intensity" for channel in channels]
   1694    missing = set(profile_columns) - set(feature_dataframe.columns)
   1695    if missing:
   1696        raise ValueError(f"Missing axial profile columns: {sorted(missing)}")
   1697    mask = np.ones(len(feature_dataframe), dtype=bool)
   1698    if selected_frame is not None:
   1699        mask &= feature_dataframe["frame"].to_numpy() == selected_frame
   1700    if selected_cell_id is not None:
   1701        mask &= feature_dataframe["cell_id"].to_numpy() == selected_cell_id
   1702    selected = feature_dataframe.loc[mask, ["frame", "cell_id", *profile_columns]]
   1703    figures = []
   1704    colors = {"C2": "#E3484D", "C3": "#CF42BA", "C4": "#27AA59", "C5": "#3186DE"}
   1705    for _, row in selected.drop_duplicates(["frame", "cell_id"]).iterrows():
   1706        fig, ax = plt.subplots(figsize=(8, 4))
   1707        for channel, column in zip(channels, profile_columns):
   1708            values = np.asarray(row[column], dtype=float)
   1709            finite = values[np.isfinite(values)]
   1710            if not len(finite):
   1711                continue
   1712            low, high = finite.min(), finite.max()
   1713            scaled = (values - low) / (high - low) if high > low else np.zeros_like(values)
   1714            ax.plot(np.linspace(0, 1, len(scaled)), scaled, label=channel,
   1715                    color=colors.get(channel), lw=1.8)
   1716        ax.set(xlabel="Relative position along cell", ylabel="Normalized intensity",
   1717               ylim=(-0.03, 1.08), title=f"Frame {row['frame']}  |  Cell {row['cell_id']}")
   1718        ax.spines[["top", "right"]].set_visible(False)
   1719        ax.grid(color="#E7EDF1", lw=0.8)
   1720        ax.legend(frameon=False, ncol=len(channels))
   1721        fig.tight_layout()
   1722        figures.append(fig)
   1723        if show:
   1724            plt.show()
   1725    return figures
   1726
   1727
   1728def plot_object_mesh(cell_id, image, chann, verbose=False):
   1729
   1730    # Look up the cell by its ID
   1731    cellobj = next((cell for cell in image.cells if cell.cell_id == cell_id), None)
   1732
   1733    if cellobj is None:
   1734        if verbose:
   1735            print(f"Cell with ID {cell_id} not found.")
   1736        return
   1737    fig = plt.figure(figsize=(12, 12))
   1738    contour = cellobj.contour
   1739    xmin = np.min(contour.T[1] - 10)
   1740    xmax = np.max(contour.T[1] + 10)
   1741    ymin = np.min(contour.T[0] - 10)
   1742    ymax = np.max(contour.T[0] + 10)
   1743
   1744    plt.xlim(xmin, xmax)
   1745    plt.ylim(ymin, ymax)
   1746    plt.imshow(image.channels[chann], cmap="gist_gray")
   1747    plt.plot(contour.T[1], contour.T[0], "-", c="yellow", lw=2)
   1748
   1749    if cellobj.object_meshdata[chann]["object_mesh"] is not None:
   1750        for cnt, msh, mdl in zip(
   1751            cellobj.object_meshdata[chann]["object_contour"],
   1752            cellobj.object_meshdata[chann]["object_mesh"],
   1753            cellobj.object_meshdata[chann]["object_midline"],
   1754        ):
   1755            if np.any(msh):
   1756                plt.plot(
   1757                    [msh[:, 1][::2], msh[:, 3][::2]],
   1758                    [msh[:, 0][::2], msh[:, 2][::2]],
   1759                    "-",
   1760                    c="cyan",
   1761                    lw=2,
   1762                )
   1763            plt.plot(cnt.T[1], cnt.T[0], "-", c="white", lw=2)
   1764
   1765    plt.title(f"Cell {cellobj.cell_id} ----- Channel {chann} ----- Frame {image.frame}")
   1766    plt.axis("off")
   1767    plt.show()
   1768
   1769
   1770import numpy as np
   1771import matplotlib.pyplot as plt
   1772
   1773
   1774import numpy as np
   1775import matplotlib.pyplot as plt
   1776
   1777
   1778import numpy as np
   1779import matplotlib.pyplot as plt
   1780
   1781
   1782def plot_object_contours(
   1783    cell_id,
   1784    image,
   1785    contour_channels,
   1786    shift=False,
   1787    base_channel=None,
   1788    verbose=False,
   1789    add_scalebar=True,
   1790    scalebar_um=1,
   1791    um_per_px=0.065,
   1792    scalebar_color="white",
   1793    scalebar_lw=8,
   1794    window_px=None,          # fixed square window size (pixels); if None -> original crop_image behavior
   1795    save_svg_path=None,
   1796    save_tiff_path=None    # SVG only (editable in Illustrator)
   1797):
   1798    """
   1799    Plot the contours of a cell and its objects for multiple channels on a single base channel image.
   1800
   1801    Behaves like your original version:
   1802      - base_channel None -> uses image.image
   1803      - shift True -> crop_image(... phase_img=image.image)
   1804      - else -> crop_image(... phase_img=None)
   1805      - plots cell contour + object contours from cropped_object_contour
   1806      - adds scalebar the same way
   1807
   1808    Added:
   1809      - window_px: optional fixed square window around cell center (bypasses crop_image).
   1810      - save_svg_path: save editable SVG (scalebar/text/contours editable in Illustrator).
   1811    """
   1812
   1813    cellobj = next((cell for cell in image.cells if cell.cell_id == cell_id), None)
   1814    if cellobj is None:
   1815        if verbose:
   1816            print(f"Cell with ID {cell_id} not found.")
   1817        return
   1818
   1819    contour = cellobj.contour
   1820
   1821    # Keep SVG text editable
   1822    plt.rcParams.update({
   1823        "font.family": "Arial",
   1824        "svg.fonttype": "none",
   1825    })
   1826
   1827    # ---------------------------------------------------------
   1828    # 1) Get cropped_img and cropped_contour exactly like before
   1829    #    unless window_px is requested.
   1830    # ---------------------------------------------------------
   1831    if window_px is None:
   1832        # Original logic (respects base_channel and shift exactly)
   1833        if base_channel is None:
   1834            cropped_img, _, cropped_contour, x_offset, y_offset = u.crop_image(
   1835                image=image.image,
   1836                contour=contour,
   1837                mask_to_crop=None,
   1838                phase_img=None,
   1839            )
   1840        elif shift:
   1841            cropped_img, _, cropped_contour, x_offset, y_offset = u.crop_image(
   1842                image=image.channels[base_channel],
   1843                contour=contour,
   1844                mask_to_crop=None,
   1845                phase_img=image.image,
   1846            )
   1847        else:
   1848            cropped_img, _, cropped_contour, x_offset, y_offset = u.crop_image(
   1849                image=image.channels[base_channel],
   1850                contour=contour,
   1851                mask_to_crop=None,
   1852                phase_img=None,
   1853            )
   1854
   1855        # Object contours: use exactly what your original code used
   1856        # (already in cropped coords)
   1857        cropped_object_contours = {}
   1858        for chann in contour_channels:
   1859            # keep it robust if key missing
   1860            cropped_object_contours[chann] = (
   1861                cellobj.object_meshdata.get(chann, {}).get("cropped_object_contour", None)
   1862            )
   1863
   1864    else:
   1865        # ---------------------------------------------------------
   1866        # Fixed window mode (consistent output size across cells)
   1867        # This mode does NOT call crop_image, so we must:
   1868        #  - choose the same base image selection logic
   1869        #  - shift the cell contour + object contours into window coords
   1870        # ---------------------------------------------------------
   1871        if base_channel is None:
   1872            base_img = image.image
   1873        else:
   1874            base_img = image.channels[base_channel]
   1875
   1876        # center of cell in full image coords
   1877        cy = int(np.mean(contour[:, 0]))
   1878        cx = int(np.mean(contour[:, 1]))
   1879        half = int(window_px // 2)
   1880
   1881        y_min = max(cy - half, 0)
   1882        y_max = min(cy + half, base_img.shape[0])
   1883        x_min = max(cx - half, 0)
   1884        x_max = min(cx + half, base_img.shape[1])
   1885
   1886        cropped_img = base_img[y_min:y_max, x_min:x_max]
   1887
   1888        # Cell contour in window coords
   1889        cropped_contour = np.column_stack([contour[:, 0] - y_min, contour[:, 1] - x_min])
   1890
   1891        # Object contours:
   1892        # In many pipelines, cropped_object_contour is in crop_image coords,
   1893        # so for window mode we instead try to use the full-res object contour
   1894        # if available, otherwise we fall back to cropped_object_contour (unchanged).
   1895        cropped_object_contours = {}
   1896        for chann in contour_channels:
   1897            meshdata = cellobj.object_meshdata.get(chann, {})
   1898
   1899            # Try common keys for full-image object contours
   1900            full_key_candidates = ["object_contour", "object_contours", "raw_object_contour", "raw_object_contours"]
   1901            full_contours = None
   1902            for k in full_key_candidates:
   1903                if k in meshdata and meshdata[k] is not None:
   1904                    full_contours = meshdata[k]
   1905                    break
   1906
   1907            if full_contours is not None:
   1908                # Shift full contours into window coords
   1909                shifted = []
   1910                for obj in full_contours:
   1911                    obj = np.asarray(obj)
   1912                    shifted.append(np.column_stack([obj[:, 0] - y_min, obj[:, 1] - x_min]))
   1913                cropped_object_contours[chann] = shifted
   1914            else:
   1915                # Fallback: use cropped_object_contour as-is (may not align in window mode,
   1916                # but prevents hard failure if only that exists)
   1917                cropped_object_contours[chann] = meshdata.get("cropped_object_contour", None)
   1918
   1919    # ---------------------------------------------------------
   1920    # 2) Plot (same operations as your original version)
   1921    # ---------------------------------------------------------
   1922    fig, ax = plt.subplots(figsize=(6, 6))
   1923    ax.imshow(cropped_img, cmap="gist_grey")
   1924    ax.plot(cropped_contour.T[1], cropped_contour.T[0], color="yellow")
   1925
   1926    for chann in contour_channels:
   1927        ccont = cropped_object_contours.get(chann, None)
   1928        if ccont is not None:
   1929            for obj_contour in ccont:
   1930                ax.plot(
   1931                    obj_contour.T[1],
   1932                    obj_contour.T[0],
   1933                    "-",
   1934                    lw=2,
   1935                    label=f"Channel {chann}",
   1936                )
   1937
   1938    if contour_channels:
   1939        ax.legend(frameon=False)
   1940
   1941    # ---------------------------------------------------------
   1942    # 3) Scalebar (same logic; editable in SVG)
   1943    # ---------------------------------------------------------
   1944    if add_scalebar:
   1945        bar_px = scalebar_um / um_per_px
   1946        h, w = cropped_img.shape[:2]
   1947
   1948        x_start = int(w * 0.05)
   1949        y_start = int(h * 0.95)
   1950
   1951        ax.plot(
   1952            [x_start, x_start + bar_px],
   1953            [y_start, y_start],
   1954            color=scalebar_color,
   1955            lw=scalebar_lw,
   1956            solid_capstyle="butt",
   1957        )
   1958        ax.text(
   1959            x_start + bar_px / 2,
   1960            y_start - h * 0.03,
   1961            f"{scalebar_um:.1f} µm",
   1962            color=scalebar_color,
   1963            ha="center",
   1964            va="bottom",
   1965            fontsize=20,
   1966            fontfamily="Arial",
   1967        )
   1968
   1969    ax.set_title(f"Cell {cellobj.cell_id} - Frame {image.frame}— Base Channel {base_channel}", fontfamily="Arial")
   1970    ax.axis("off")
   1971    fig.tight_layout(pad=0)
   1972    if save_tiff_path:
   1973        plt.imsave(
   1974            save_tiff_path,
   1975            cropped_img,
   1976            cmap="gist_grey",
   1977            format="tiff"
   1978        )
   1979    # SVG only (editable)
   1980    if save_svg_path:
   1981        fig.savefig(save_svg_path, format="svg", bbox_inches="tight", pad_inches=0)
   1982
   1983    plt.show()
   1984
   1985
   1986
   1987def plot_mesh_scientific(cell_id, image, crop_size=300, verbose=False):
   1988
   1989    cellobj = next((cell for cell in image.cells if cell.cell_id == cell_id), None)
   1990
   1991    if cellobj is None:
   1992        if verbose:
   1993            print(f"Cell with ID {cell_id} not found.")
   1994        return
   1995
   1996    fig = plt.figure(figsize=(12, 12))
   1997
   1998    profile_mesh = cellobj.profile_mesh
   1999    # Calculate the center of the contour
   2000    center_x = np.mean(cellobj.contour.T[1])
   2001    center_y = np.mean(cellobj.contour.T[0])
   2002
   2003    # Calculate the new cropping boundaries based on the fixed crop size
   2004    half_crop_size = crop_size / 2
   2005    xmin = center_x - half_crop_size
   2006    xmax = center_x + half_crop_size
   2007    ymin = center_y - half_crop_size
   2008    ymax = center_y + half_crop_size
   2009
   2010    plt.xlim(xmin, xmax)
   2011    plt.ylim(ymin, ymax)
   2012    plt.imshow(image.image, cmap="gist_gray")
   2013    plt.plot(cellobj.contour.T[1], cellobj.contour.T[0], "-", c="w", lw=3)
   2014
   2015    plt.plot(
   2016        [cellobj.y1[::2], cellobj.y2[::2]],
   2017        [cellobj.x1[::2], cellobj.x2[::2]],
   2018        c="cyan",
   2019        lw=3,
   2020    )
   2021    plt.plot(cellobj.midline.T[1], cellobj.midline.T[0], "-", c="orange", lw=3)
   2022    # Add scale bar
   2023    scale_length_um = 1
   2024    scale_length_px = scale_length_um / 0.065  # Convert µm to pixels
   2025    scale_x = (xmin + xmax) / 2 - scale_length_px / 2  # Center the scale bar
   2026    scale_y = ymin + (ymax - ymin) * 0.02  # Adjust the position of the scale bar
   2027    plt.plot(
   2028        [scale_x, scale_x + scale_length_px], [scale_y, scale_y], color="white", lw=5
   2029    )
   2030    # plt.scatter(profile_mesh[1][::2].T, profile_mesh[0][::2].T, s = 3, c = 'cyan')
   2031    plt.title(
   2032        f"Cell number: {cellobj.cell_id} with contour and profiling mesh", fontsize=16
   2033    )
   2034    plt.show()
   2035
   2036
   2037def plot_width(cell, image, crop_size=70):
   2038    cellobj = image.cells[cell]
   2039    other_cellobj = image.cells[cell + 44]
   2040
   2041    fig1, axs1 = plt.subplots(figsize=(12, 12))
   2042    profile_mesh = cellobj.profile_mesh
   2043
   2044    center_x = np.mean(cellobj.contour.T[1])
   2045    center_y = np.mean(cellobj.contour.T[0])
   2046
   2047    # Calculate the new cropping boundaries based on the fixed crop size
   2048    half_crop_size = crop_size / 2
   2049    xmin = center_x - half_crop_size
   2050    xmax = center_x + half_crop_size
   2051    ymin = center_y - half_crop_size
   2052    ymax = center_y + half_crop_size
   2053
   2054    axs1.set_position([0.1, 0.1, 0.5, 0.8])
   2055    axs1.imshow(image.image, cmap="gist_gray")
   2056    axs1.plot(cellobj.contour.T[1], cellobj.contour.T[0], lw=3, c="orange")
   2057
   2058    axs1.set_xlim(xmin, xmax)
   2059    axs1.set_ylim(ymin, ymax)
   2060    # Add scale bar
   2061    scale_length_um = 1
   2062    scale_length_px = scale_length_um / 0.065  # Convert µm to pixels
   2063    scale_x = (xmin + xmax) / 2 - scale_length_px / 2  # Center the scale bar
   2064    scale_y = ymin + (ymax - ymin) * 0.02  # Adjust the position of the scale bar
   2065    axs1.plot(
   2066        [scale_x, scale_x + scale_length_px], [scale_y, scale_y], color="white", lw=5
   2067    )
   2068    axs1.set_xticks([])
   2069    axs1.set_yticks([])
   2070    axs1.set_title("Cell " + str(cell), fontname="Arial", fontsize=14)
   2071
   2072    fig2, axs2 = plt.subplots(figsize=(6, 4))
   2073
   2074    # Plot the cell width for the first cell (red curve)
   2075    axs2.plot(cellobj.profiling_data["cell_widthno"], c="orange")
   2076
   2077    # Plot the cell width for the other cell (green curve)
   2078    axs2.plot(other_cellobj.profiling_data["cell_widthno"], c="b")
   2079
   2080    # Set x-axis ticks and labels
   2081    newxticks = [
   2082        np.round(
   2083            cellobj.contour_features["cell_length"] * (tick / profile_mesh.shape[1]), 1
   2084        )
   2085        for tick in axs2.get_xticks()
   2086    ]
   2087    axs2.set_xticks(axs2.get_xticks())
   2088    axs2.set_xticklabels(newxticks, fontname="Arial", fontsize=12)
   2089
   2090    # Set y-axis ticks and label
   2091    newyticks = [np.round(tick, 1) for tick in axs2.get_yticks()]
   2092    axs2.set_yticks(axs2.get_yticks())
   2093    axs2.set_yticklabels(newyticks, fontname="Arial", fontsize=12)
   2094
   2095    # Set x-axis to start at 0
   2096    axs2.set_xlim(left=0)
   2097    axs2.set_ylim(bottom=0)
   2098    # Set x and y axis labels
   2099    axs2.set_xlabel("cell length [µm]", fontname="Arial", fontsize=12)
   2100    axs2.set_ylabel("cell width [µm]", fontname="Arial", fontsize=12)
   2101
   2102    axs2.legend(["overlapping cell", "single cell"], frameon=False)
   2103
   2104    fig3, axs3 = plt.subplots(figsize=(12, 12))
   2105    other_profile_mesh = other_cellobj.profile_mesh
   2106
   2107    center_x = np.mean(other_cellobj.contour.T[1])
   2108    center_y = np.mean(other_cellobj.contour.T[0])
   2109
   2110    # Calculate the new cropping boundaries based on the fixed crop size
   2111    half_crop_size = crop_size / 2
   2112    xmin = center_x - half_crop_size
   2113    xmax = center_x + half_crop_size
   2114    ymin = center_y - half_crop_size
   2115    ymax = center_y + half_crop_size
   2116
   2117    axs3.set_position([0.1, 0.1, 0.5, 0.8])
   2118    axs3.imshow(image.image, cmap="gist_gray")
   2119    axs3.plot(other_cellobj.contour.T[1], other_cellobj.contour.T[0], lw=3, c="b")
   2120
   2121    axs3.set_xlim(xmin, xmax)
   2122    axs3.set_ylim(ymin, ymax)
   2123
   2124    axs3.set_xticks([])
   2125    axs3.set_yticks([])
   2126    axs3.set_title("Other Cell", fontname="Arial", fontsize=14)
   2127    scale_length_um = 1
   2128    scale_length_px = scale_length_um / 0.065  # Convert µm to pixels
   2129    scale_x = (xmin + xmax) / 2 - scale_length_px / 2  # Center the scale bar
   2130    scale_y = ymin + (ymax - ymin) * 0.02  # Adjust the position of the scale bar
   2131    axs3.plot(
   2132        [scale_x, scale_x + scale_length_px], [scale_y, scale_y], color="white", lw=5
   2133    )
   2134
   2135    plt.show()
   2136
   2137
   2138def plot_signal_profile(cell, image):
   2139    cellobj = image.cells[cell]
   2140    profile_mesh = cellobj.profile_mesh
   2141    profiling_mesh = cellobj.profiling_data["phaco_mesh_intensity"]
   2142    fig, ax = plt.subplots(2, 1, figsize=(10, 10))
   2143    xmin = np.min(profile_mesh[1] - 10)
   2144    xmax = np.max(profile_mesh[1] + 10)
   2145    ymin = np.min(profile_mesh[0] - 10)
   2146    ymax = np.max(profile_mesh[0] + 10)
   2147    ax[0].imshow(profiling_mesh, aspect="auto")
   2148    ax[0].get_xticks()
   2149    xticks = ax[0].get_xticks()
   2150    newxticks = [
   2151        np.round(
   2152            cellobj.contour_features["cell_length"] * (tick / profiling_mesh.shape[1]),
   2153            1,
   2154        )
   2155        for tick in xticks
   2156    ]
   2157    ax[0].set_xticklabels(newxticks, fontname="Arial", fontsize=12)
   2158    ax[0].set_yticks([])
   2159    ax[0].set_ylabel("signal\nstraighten image\n", fontname="Arial", fontsize=12)
   2160    ax[0].set_xlabel("cell length [µm]", fontname="Arial", fontsize=12)
   2161
   2162    ax[1].imshow(image, cmap="gist_gray", aspect="auto")
   2163    ax[1].plot(cellobj.contour.T[1], cellobj.contour.T[0], "-", c="r")
   2164    ax[1].plot(cellobj.midline.T[1], cellobj.midline.T[0], "-", c="y")
   2165    ax[1].set_xlim(xmin, xmax)
   2166    ax[1].set_ylim(ymin, ymax)
   2167    plt.show()
   2168
   2169from scipy.interpolate import interp1d
   2170
   2171def plot_demograph_rotated(cell_lengths, normalized_average_mesh_intensity, title='Demograph Plot',
   2172                           y_label='Normalized Distance From Midcell (µm)', x_label='Cell Length Percentile',
   2173                           cmap='rainbow', cbar_title='DnaN-msfGFP Normalized Intensity'):
   2174    # Filter out NaN values from cell_lengths and corresponding intensities
   2175    valid_indices = [i for i, arr in enumerate(normalized_average_mesh_intensity) if arr is not None and not np.isnan(np.nanmean(arr))]
   2176    cell_lengths = cell_lengths[valid_indices].reset_index(drop=True)
   2177    normalized_average_mesh_intensity = normalized_average_mesh_intensity[valid_indices].reset_index(drop=True)
   2178
   2179
   2180    # Sort the arrays based on cell_lengths
   2181    sorted_indices = np.argsort(cell_lengths)
   2182    sorted_lengths = cell_lengths[sorted_indices]
   2183    sorted_intensities = normalized_average_mesh_intensity[sorted_indices]
   2184
   2185    # Reset the indices to ensure correct alignment
   2186    sorted_lengths = np.array(sorted_lengths)
   2187    sorted_intensities = np.array(sorted_intensities)
   2188
   2189    # Find the maximum length of the arrays
   2190    max_length = max(len(arr) for arr in sorted_intensities)
   2191
   2192    # Interpolate missing values for shorter arrays
   2193    interpolated_arrays = []
   2194    for array in sorted_intensities:
   2195        x = np.arange(len(array))  # x-coordinates for the existing intensity values
   2196        f = interp1d(x, array, kind='linear', fill_value='extrapolate')  # Interpolation function
   2197        interpolated_array = f(np.linspace(0, len(array), max_length))  # Interpolated array
   2198        interpolated_arrays.append(interpolated_array)
   2199    stacked_demograph = np.vstack(interpolated_arrays)
   2200
   2201    fig1 = plt.figure(figsize=(10, 8))  # Adjust figure size as needed
   2202
   2203    ax = plt.subplot(111)
   2204    image = ax.imshow(stacked_demograph.T, aspect='auto', cmap=cmap)  # Transpose the array
   2205    cbar = plt.colorbar(image)  # Use the image as the mappable object for the colorbar
   2206    cbar.set_label(cbar_title, rotation=90, labelpad=20, fontsize=14)  # Set colorbar label
   2207    cbar.ax.tick_params(labelsize=14)  # Adjust colorbar tick label size
   2208
   2209    # Calculate y-axis values and middle index
   2210    y_axis_values = np.linspace(-1, 1, max_length)  # Ensure values are between -1 and 1
   2211    middle_index = len(y_axis_values) // 2
   2212
   2213    # Set y-axis ticks and labels
   2214    plt.yticks([-0.5, middle_index, len(y_axis_values)-0.5], [-1, 0, 1])
   2215    plt.xticks([0, len(sorted_lengths) * 0.25, len(sorted_lengths) * 0.5, len(sorted_lengths) * 0.75, len(sorted_lengths)],
   2216               ['0', '25', '50', '75', '100'])
   2217    plt.xlabel(x_label, fontsize=16)
   2218    plt.ylabel(y_label, fontsize=16)
   2219
   2220    plt.title(title, fontsize=18)  # Set the plot title
   2221
   2222    plt.show()
   2223
   2224
   2225def list_field_ids(dataset_dir, phase_channel="C1"):
   2226    dataset_dir = Path(dataset_dir)
   2227    return sorted(
   2228        {
   2229            path.stem.rsplit("_", 1)[0]
   2230            for path in dataset_dir.glob(f"*_{phase_channel}.tiff")
   2231        }
   2232    )
   2233
   2234
   2235def plot_multichannel_overview(dataset_dir, field_id, channel_labels, figsize=(14, 3.8)):
   2236    """Return aligned display crops for each channel in a field of view."""
   2237    fig = plot_field_channels(dataset_dir, field_id, channel_labels)
   2238    if figsize is not None:
   2239        fig.set_size_inches(*figsize)
   2240        fig.tight_layout(rect=(0, 0, 1, 0.98))
   2241    return fig
   2242
   2243
   2244def plot_mesh_overlay(image_object, cell_id, crop_size=150, verbose=False):
   2245    cellobj = next((cell for cell in image_object.cells if cell.cell_id == cell_id), None)
   2246    if cellobj is None:
   2247        if verbose:
   2248            print(f"Cell with ID {cell_id} not found.")
   2249        return
   2250
   2251    contour = np.asarray(cellobj.contour)
   2252    mesh = np.asarray(cellobj.mesh)
   2253    midline = np.asarray(cellobj.midline)
   2254
   2255    center_x = contour[:, 1].mean()
   2256    center_y = contour[:, 0].mean()
   2257    half_crop = crop_size / 2
   2258    xmin = center_x - half_crop
   2259    xmax = center_x + half_crop
   2260    ymin = center_y - half_crop
   2261    ymax = center_y + half_crop
   2262
   2263    fig, ax = plt.subplots(figsize=(5.5, 5.5))
   2264    ax.imshow(image_object.image, cmap="gray")
   2265    ax.plot(contour[:, 1], contour[:, 0], color="white", lw=2, label="Cell contour")
   2266    if mesh.size:
   2267        ax.plot(
   2268            [mesh[:, 1][::2], mesh[:, 3][::2]],
   2269            [mesh[:, 0][::2], mesh[:, 2][::2]],
   2270            color="cyan",
   2271            lw=1,
   2272        )
   2273    ax.plot(midline[:, 1], midline[:, 0], color="magenta", lw=2, label="Midline")
   2274    ax.set_xlim(xmin, xmax)
   2275    ax.set_ylim(ymin, ymax)
   2276    ax.set_title(f"Mesh overlay for cell {cell_id}")
   2277    ax.axis("off")
   2278    ax.legend(frameon=False, loc="upper right")
   2279    plt.tight_layout()
   2280    plt.show()
   2281
   2282
   2283def pick_detected_cell(image_object, preferred_channels=("C4", "C5")):
   2284    """Choose an interior, isolated cell with a detected object when possible."""
   2285    candidates = []
   2286    all_centers = []
   2287    for cell in image_object.cells:
   2288        contour = np.asarray(cell.contour)
   2289        if contour.ndim != 2 or contour.shape[1] != 2 or not len(contour):
   2290            continue
   2291        center = contour.mean(axis=0)
   2292        all_centers.append(center)
   2293        detections = sum(
   2294            len(contours) for channel in preferred_channels
   2295            for contours in [cell.object_meshdata.get(channel, {}).get("object_contour")]
   2296            if contours is not None
   2297        )
   2298        if detections:
   2299            margin = min(center[0], center[1],
   2300                         image_object.image.shape[0] - center[0],
   2301                         image_object.image.shape[1] - center[1])
   2302            length = np.ptp(contour, axis=0).max()
   2303            candidates.append((cell.cell_id, center, margin, length))
   2304    if not candidates:
   2305        return pick_representative_cell([image_object])[1]
   2306    median_length = np.median([item[3] for item in candidates])
   2307    tree = cKDTree(all_centers) if len(all_centers) > 1 else None
   2308    def score(item):
   2309        distance = tree.query(item[1], k=2)[0][1] if tree is not None else np.inf
   2310        return (item[2] >= 60, distance / max(item[3], 1),
   2311                -abs(item[3] - median_length), item[2])
   2312    return max(candidates, key=score)[0]
   2313
   2314
   2315def plot_morphology_summary(df):
   2316    fig, axes = plt.subplots(2, 2, figsize=(13, 10))
   2317
   2318    axes[0, 0].hist(df["cell_length"].dropna(), bins=30, color="#4c72b0", alpha=0.9)
   2319    axes[0, 0].set_title("Cell length distribution")
   2320    axes[0, 0].set_xlabel("Length (um)")
   2321    axes[0, 0].set_ylabel("Cell count")
   2322
   2323    axes[0, 1].hist(df["cell_width"].dropna(), bins=30, color="#55a868", alpha=0.9)
   2324    axes[0, 1].set_title("Cell width distribution")
   2325    axes[0, 1].set_xlabel("Width (um)")
   2326    axes[0, 1].set_ylabel("Cell count")
   2327
   2328    scatter = axes[1, 0].scatter(
   2329        df["cell_length"],
   2330        df["cell_area"],
   2331        c=df["C5_cell_total_obj_area"].fillna(0),
   2332        cmap="magma",
   2333        s=20,
   2334        alpha=0.75,
   2335    )
   2336    axes[1, 0].set_title("Cell area vs. cell length")
   2337    axes[1, 0].set_xlabel("Length (um)")
   2338    axes[1, 0].set_ylabel("Area (um^2)")
   2339    cbar = plt.colorbar(scatter, ax=axes[1, 0])
   2340    cbar.set_label("Nucleoid object area")
   2341
   2342    image_names = sorted(df["image_name"].dropna().unique())
   2343    width_groups = [
   2344        df.loc[df["image_name"] == image_name, "cell_width"].dropna()
   2345        for image_name in image_names
   2346    ]
   2347    axes[1, 1].boxplot(width_groups, labels=image_names, patch_artist=True)
   2348    axes[1, 1].set_title("Cell width by field of view")
   2349    axes[1, 1].set_ylabel("Width (um)")
   2350    axes[1, 1].tick_params(axis="x", rotation=30)
   2351
   2352    plt.tight_layout()
   2353    plt.show()
   2354
   2355
   2356def plot_signal_summary(df):
   2357    fig, axes = plt.subplots(1, 3, figsize=(16, 4.5))
   2358
   2359    axes[0].hist(
   2360        df["C2_signal_constriction_degree_mesh"].dropna(),
   2361        bins=30,
   2362        color="#dd8452",
   2363        alpha=0.9,
   2364    )
   2365    axes[0].set_title("Outer-membrane signal localization")
   2366    axes[0].set_xlabel("C2 signal constriction degree")
   2367    axes[0].set_ylabel("Cell count")
   2368
   2369    max_objects = int(
   2370        np.nanmax([df["C4_object_number"].max(), df["C5_object_number"].max()])
   2371    )
   2372    bins = np.arange(0.5, max_objects + 1.6, 1)
   2373    axes[1].hist(
   2374        df["C4_object_number"].dropna(),
   2375        bins=bins,
   2376        alpha=0.7,
   2377        label="RNA objects",
   2378        color="#55a868",
   2379    )
   2380    axes[1].hist(
   2381        df["C5_object_number"].dropna(),
   2382        bins=bins,
   2383        alpha=0.7,
   2384        label="Nucleoid objects",
   2385        color="#c44e52",
   2386    )
   2387    axes[1].set_title("Detected object counts per cell")
   2388    axes[1].set_xlabel("Objects per cell")
   2389    axes[1].set_ylabel("Cell count")
   2390    axes[1].legend(frameon=False)
   2391
   2392    scatter = axes[2].scatter(
   2393        df["C4_NC_ratio"],
   2394        df["C5_NC_ratio"],
   2395        c=df["cell_length"],
   2396        cmap="viridis",
   2397        s=20,
   2398        alpha=0.75,
   2399    )
   2400    axes[2].set_title("RNA versus nucleoid occupancy")
   2401    axes[2].set_xlabel("C4 object area / cell area")
   2402    axes[2].set_ylabel("C5 object area / cell area")
   2403    cbar = plt.colorbar(scatter, ax=axes[2])
   2404    cbar.set_label("Cell length (um)")
   2405
   2406    plt.tight_layout()
   2407    plt.show()
   2408
   2409
   2410def plot_pairwise_metric_heatmap(
   2411    df,
   2412    metric_suffix,
   2413    channel_order=("C2", "C3", "C4", "C5"),
   2414):
   2415    matching_columns = [column for column in df.columns if column.endswith(metric_suffix)]
   2416    if not matching_columns:
   2417        print(f"No columns found for metric suffix: {metric_suffix}")
   2418        return
   2419
   2420    matrix = pd.DataFrame(np.nan, index=channel_order, columns=channel_order, dtype=float)
   2421    for column in matching_columns:
   2422        ch1, ch2, _ = column.split("_", 2)
   2423        value = df[column].median(skipna=True)
   2424        matrix.loc[ch1, ch2] = value
   2425        matrix.loc[ch2, ch1] = value
   2426
   2427    if "correlation" in metric_suffix:
   2428        np.fill_diagonal(matrix.values, 1.0)
   2429        vmin, vmax, cmap = -1, 1, "coolwarm"
   2430    else:
   2431        vmin = np.nanmin(matrix.values)
   2432        vmax = np.nanmax(matrix.values)
   2433        cmap = "viridis"
   2434
   2435    fig, ax = plt.subplots(figsize=(6, 5))
   2436    image = ax.imshow(matrix.values.astype(float), cmap=cmap, vmin=vmin, vmax=vmax)
   2437    ax.set_xticks(range(len(channel_order)))
   2438    ax.set_yticks(range(len(channel_order)))
   2439    ax.set_xticklabels(channel_order)
   2440    ax.set_yticklabels(channel_order)
   2441    ax.set_title(metric_suffix.replace("_", " "))
   2442
   2443    for row in range(len(channel_order)):
   2444        for col in range(len(channel_order)):
   2445            value = matrix.iloc[row, col]
   2446            if not np.isnan(value):
   2447                ax.text(col, row, f"{value:.2f}", ha="center", va="center", color="black")
   2448
   2449    plt.colorbar(image, ax=ax, fraction=0.046, pad=0.04)
   2450    plt.tight_layout()
   2451    plt.show()