Automated analysis of zebrafish vascular networks using the VISTA-Z pipeline.

Code ↔ Paper

6 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 6 matches

[1] § Materials and methods › Image pre-processing and segmentation ↔ Codes/VISTA-Z_ISV_loss.ipynb, lines 451–548 · score 0.89 · curated segmentation mask, manually curated, sigma range, segment labelling, Maximum intensity projections, binary masks
[2] § Materials and methods › Image pre-processing and segmentation ↔ Codes/VISTA-Z_lateral.ipynb, lines 328–423 · score 0.88 · curated segmentation mask, manually curated, sigma range, segment labelling, Maximum intensity projections, binary masks
[3] § Materials and methods › Automated image-based skeleton, vessel metric analysis and outlier filtering ↔ Codes/VISTA-Z_ROI.ipynb, lines 84–171 · score 0.84 · median absolute deviation, metric normalisation, vessel masks, vessel length, transform, windows
[4] § Materials and methods › Automated image-based skeleton, vessel metric analysis and outlier filtering ↔ Codes/VISTA-Z_lateral.ipynb, lines 76–139 · score 0.83 · median absolute deviation, metric normalisation, vessel masks, vessel length, transform, windows
[5] § Results › Development of VISTA-Z, an automated pipeline for vascular quantification ↔ Codes/VISTA-Z_ROI.ipynb, lines 84–171 · score 0.80 · Median Absolute Deviation, image pre processing, vessel masking, vessel length, noise, selection
[6] § Results › Development of VISTA-Z, an automated pipeline for vascular quantification ↔ Codes/VISTA-Z_lateral.ipynb, lines 76–139 · score 0.79 · Median Absolute Deviation, image pre processing, vessel masking, vessel length, noise, selection

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

Jupyter notebook · 1,054 lines · 61 KB · MIT · 3 matches

# %% [markdown]
# **Table of contents**<a id='toc0_'></a>
# - [Import packages](#toc1_)
# - [Set up working directory](#toc2_)
# - [Vessel analysis - Lateral images](#toc3_)
# - [Define parameters](#toc3_1_)
# - [Define functions](#toc3_2_)
# - [Run samples - Single channel](#toc3_3_)
# - [Run samples - Two or more channels](#toc3_4_)
#
# <!-- vscode-jupyter-toc-config
# numbering=false
# anchor=true
# flat=false
# minLevel=1
# maxLevel=6
# /vscode-jupyter-toc-config -->
#
# %% [markdown]
# # <a id='toc1_'></a>[Import packages](#toc0_)
# %%
import sys
import os
import numpy as np
import vessel_metrics as vm
import czifile
from aicspylibczi import CziFile
import matplotlib.pyplot as plt
import cv2
import imageio
import matplotlib.pyplot as plt
from matplotlib import rcParams
import matplotlib.colors
from matplotlib.pyplot import rc_context
import math
import seaborn as sns
from scipy.stats import median_abs_deviation
import glob
import time
import pandas as pd
import openpyxl
import seaborn as sns
import plotly.graph_objects as go
from skimage.measure import label, regionprops
from skimage import color
from skimage.color import label2rgb
from matplotlib.cm import get_cmap
from skimage.morphology import remove_small_objects
from aicsimageio import AICSImage, readers
from scipy.spatial import distance
from skimage.morphology import skeletonize
import traceback
import gc
# %% [markdown]
# # <a id='toc2_'></a>[Set up working directory](#toc0_)
# %%
"""
Define the path to the working directory and the output folder where the 3D volumes will be saved.
Define the name of the files to be processed.
"""
data_path = 'path_to_working_directory'
output_path = 'path_to_output_directory'
file_name = glob.glob(f'{data_path}*.*', recursive = True) # List of file names to process
# %%
"""
Ensure the output directory exists and that files are read correctly.
"""
total_files = len(file_name)
print(f"Wolking forlder contains {total_files} files")
# %% [markdown]
# # <a id='toc3_'></a>[Vessel analysis - Lateral images](#toc0_)
# %% [markdown]
# ## <a id='toc3_1_'></a>[Define parameters](#toc0_)
# %%
## Please README before running the following code ##
## Parameters that the user should define before running VISTA-Z ROI analysis ##
# clahe = (int, int) # CLAHE parameters (clipLimit, tileGridSize)
# channel = int [0,1,2] # Channel number to be analysed per sample (only if multi-channel images) [0]: mCherry, [1]: EGFP, [2]: DAPI, etc.
# name_start = str # Starting string of the file names to be processed
# pdf_dpi = int # DPI for saving PDF figures
# tiles = int # Number of tiles used to image the whole zebrafish trunk
# min_vessel_size = int # Minimum vessel size (in pixels) to be thresholded for analysis
# Vessel segmentation parameters:
# seg_method = str # Vessel segmentation method: 'meijering', 'frangi', 'sato' or 'jerman'
# sigma1 = range(int, int, int) # Sigma range for vessel enhancement filter
# hole_size = int # Maximum hole size to be filled in the vessel mask
# ditzle_size = int # Maximum size of small objects to be removed from the vessel mask
# thresh = int # Threshold value to binarize the vessel enhanced image
# skel_method = str # Skeletonisation method: 'lee' or 'zhang'
# Branch_dist = int # Maximum distance (in pixels) to consider two vessel branches as connected
# MAD filtering threshold: Threshold for median absolute deviation (MAD) filtering of vessel metrics
# threshold_low = float
# threshold_high = float
# Normalisation of vessel metrics:
# window_size = 10 # 10 is to make an 100µm² area of 10x10. Transform into the low int value to create the scale for the density values
# scale = float # Scale factor to normalise vessel metrics (pixels per micron)
# output_name = str # Name of the output files to be saved (excel file)
# %%
# Image pre-processing and ROI selection parameters
clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(4, 4))
channel = 0
name_start = 'kdrlmCh'
pdf_dpi = 300
tiles = 3
# Vessel segmentation parameters
min_vessel_size = 500 # Note: User are reccommended to set this parameter after initial analysis based on median(vessel_length) - std(vessel_length)
seg_method = 'meijering'
sigma1 = range(3, 8, 1)
hole_size = 200
ditzle_size = 50 # Note: Users are reccommended to set hole_size and ditzle_size based on the visualisation of the vessel mask after thresholding. Adjust these parameters to ensure that the vessel mask captures the vessels accurately while minimising noise and artefacts.
thresh = 10 # Note: Users are reccommended to set this parameter based on the visualisation of the vessel enhanced image and the vessel mask after thresholding. Adjust this parameter to ensure that the vessel mask captures the vessels accurately while minimising noise and artefacts.
skel_method = 'lee'
#Branchpoint refinement parameters
Branch_dist = 30 # Note: User are reccommended to set this parameter based median(vessel_length)
# MAD filtering threshold
threshold_low = 1
threshold_high = 3
# Diameter calculation and vessel metrics normalisation parameters
window_size = 10 # 10 is to make an 100µm² area of 10x10. Transform into the low int value to create the scale for the density values
scale = 1.2044
# Output file name
output_name = "VISTA-Z_lateral_Results"
output_tile_name = "VISTA-Z_lateral_tile_Results"
# %% [markdown]
# ## <a id='toc3_2_'></a>[Define functions](#toc0_)
# %%
def load_and_process_czi_single(file_path, file_name, output_path):
"""
Loads and processes a fluorescent image (.czi file) with a single channel,
Args:
file_path (str): Full path to the .czi file.
file_name (str): File name used to derive the sample ID.
output_path (str): Directory where output images will be saved.
Returns:
tuple:
image (ndarray): Raw image data after squeezing dimensions.
mip_image (ndarray): Normalized maximum-intensity projection (MIP) image in uint8.
gray_image (ndarray): Greyscale MIP image.
clahe_image (ndarray): CLAHE-enhanced greyscale image.
otsu_mask (ndarray): Binary mask from Otsu thresholding.
output_mip_path (str): Path for saving the MIP image.
output_mask_path (str): Path for saving the Otsu mask image.
Returns (None, None, None, None) if processing fails.
"""
try:
sampleID = file_name.split('_MaxInt')[0] # Name of the sample
with czifile.CziFile(file_path) as czi:
data = czi.asarray() # Read confocal image data from .czi file
data_squeezed = np.squeeze(data) # Squeeze data to remove unwanted dimensions/metadata and retrieve image data
image = data_squeezed
mip_image = cv2.normalize(image, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8) # Normalise image to uint8
mip_image_rgb = cv2.cvtColor(mip_image, cv2.COLOR_GRAY2RGB)
gray_image = cv2.cvtColor(mip_image_rgb, cv2.COLOR_RGB2GRAY)
clahe_image = clahe.apply(gray_image) # Apply CLAHE for contrast enhacement
_, otsu_mask = cv2.threshold(clahe_image, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU) # Otsu thresholding
output_mip_path = os.path.join(output_path, f"{sampleID}.tiff") # Define output MIP path
output_mask_path = os.path.join(output_path, f"{sampleID}_mask.tiff") # Define output mask path
return image, mip_image, gray_image, clahe_image, otsu_mask, output_mip_path, output_mask_path
except Exception as e:
print(f"Error processing {sampleID}: {e}") # Print error message if processing fails
return None, None, None, None
def load_and_process_czi_double(file_path, file_name, output_path, channel = channel):
"""
Loads and processes a fluorescent image (.czi file) with two or more channels,
Args:
file_path (str): Full path to the .czi file.
file_name (str): File name used to derive the sample ID.
output_path (str): Directory where output images will be saved.
channel (int): Channel index. For example mCherry = 0, EGFP = 1, DAPI = 2, etc.
Returns:
tuple:
image (ndarray): Raw image data after squeezing dimensions.
mip_image (ndarray): Normalized maximum-intensity projection (MIP) image in uint8.
gray_image (ndarray): Greyscale MIP image.
clahe_image (ndarray): CLAHE-enhanced greyscale image.
otsu_mask (ndarray): Binary mask from Otsu thresholding.
output_mip_path (str): Path for saving the MIP image.
output_mask_path (str): Path for saving the Otsu mask image.
Returns (None, None, None, None) if processing fails.
"""
try:
sampleID = file_name.split('_MaxInt')[0] # Name of the sample
with czifile.CziFile(file_path) as czi:
data = czi.asarray() # Read confocal image data from .czi file
data_squeezed = np.squeeze(data) # Squeeze data to remove unwanted dimensions/metadata
image = data_squeezed[channel] # Retrieve image data from the specified channel
mip_image = cv2.normalize(image, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8) # Normalise image to uint8
mip_image_rgb = cv2.cvtColor(mip_image, cv2.COLOR_GRAY2RGB)
gray_image = cv2.cvtColor(mip_image_rgb, cv2.COLOR_RGB2GRAY)
clahe_image = clahe.apply(gray_image) # Apply CLAHE for contrast enhacement
_, otsu_mask = cv2.threshold(clahe_image, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU) # Otsu thresholding
output_mip_path = os.path.join(output_path, f"{sampleID}.tiff") # Define output MIP path
output_mask_path = os.path.join(output_path, f"{sampleID}_mask.tiff") # Define output mask path
return image, mip_image, gray_image, clahe_image, otsu_mask, output_mip_path, output_mask_path
except Exception as e:
print(f"Error processing {sampleID}: {e}") # Print error message if processing fails
return None, None, None, None
def save_plot_as_pdf(fig, output_path, dpi=pdf_dpi):
"""
Save a Matplotlib figure to a PDF file.
Args:
fig (matplotlib.figure.Figure): Figure to save.
output_path (str): Output path to save the PDF file.
dpi (int): Resolution to use when saving.
Returns:
None
"""
try:
fig.savefig(output_path, format='pdf', dpi=dpi, bbox_inches='tight') # Save figure as PDF
print(f"Plot saved as PDF: {output_path}")
except Exception as e: # Print error message if saving fails
print(f"Error saving plot as PDF: {e}")
def display_and_save_mip_otsu(gray_image, clahe_image, otsu_mask, output_mip_path, output_mask_path):
"""
Display MIP, CLAHE, and Otsu mask panels and save images.
Args:
gray_image (ndarray): Greyscale MIP image.
clahe_image (ndarray): CLAHE-enhanced greyscale image.
otsu_mask (ndarray): Binary mask from Otsu thresholding.
output_mip_path (str): Path for saving the MIP image.
output_mask_path (str): Path for saving the Otsu mask image.
Returns:
None
"""
fig, axes = plt.subplots(1, 3, figsize=(8, 5)) # Create figure with 3 panels
axes[0].imshow(gray_image, cmap="gray") # Display MIP
axes[0].set_title("Maximum Intensity Projection (MIP)")
axes[0].axis("off")
axes[1].imshow(clahe_image, cmap="gray") # Display CLAHE processed image
axes[1].set_title("Clahe processed image")
axes[1].axis("off")
axes[2].imshow(otsu_mask, cmap="gray") # Display Otsu's threshold mask
axes[2].set_title("Otsu’s Threshold Mask")
axes[2].axis("off")
plt.show() # Show the figure
cv2.imwrite(output_mip_path, gray_image) # Save MIP image
cv2.imwrite(output_mask_path, otsu_mask) # Save Otsu mask image
print(f"Images saved: {output_mip_path} and {output_mask_path}")
def visualize_segmented_regions(gray_image, seg_im, cmap_name='hsv', min_vessel_size = min_vessel_size):
"""
Visualise labelled vessel segments overlaid on the greyscale image.
Args:
gray_image (ndarray): Greyscale MIP image for background.
seg_im (ndarray): Binary vessel segmentation mask.
cmap_name (str): Matplotlib colormap name for segment colouring.
min_vessel_size (int): Minimum segment size to keep (pixels).
Returns:
labeled_segments (ndarray): Labelled segmentation image with unique integer labels for each segment.
Returns None on failure.
"""
try:
gray_image_norm = (gray_image - gray_image.min()) / (gray_image.max() - gray_image.min()) # Normalise greyscale image for better visualisation
labeled_segments = label(seg_im) # Label each vessel segment
num_labels = np.max(labeled_segments) # Get number of unique segments
labeled_segments = remove_small_objects(labeled_segments, min_size=min_vessel_size) # Remove small vessel segments based on minimum vessel size
if num_labels == 0: # Check if any segments are detected
print("No segments detected.")
return None
cmap = matplotlib.colormaps.get_cmap(cmap_name) # Generate colours from the chosen colourmap
colors = cmap(np.linspace(0, 1, num_labels))[:, :3] # Get RGB values for each label
overlay = label2rgb(labeled_segments, image=gray_image_norm, bg_label=0, alpha=0.6, colors=colors) # Overlay coloured segments on greyscale image
# Create a figure
fig, ax = plt.subplots(figsize=(8, 6)) # Set figure size
ax.imshow(overlay, cmap='gray') # Display overlay image
ax.set_title("Vessel Segmentation Visualization")
for region in regionprops(labeled_segments): # Loop through each labeled region
centroid = region.centroid # Plot segment labels at centroids
ax.text(centroid[1], centroid[0], str(region.label), color='white', fontsize=8, ha='center', va='center', fontweight='bold')
ax.axis('off') # Remove axes
plt.show() # Show the figure
return labeled_segments # Return labeled segments for further processing
except Exception as e:
print(f"Error visualizing segmented regions: {e}") # Print error message if visualization fails
return None
def remove_selected_segments(labeled_segments, remove_list):
"""
Manual curation: Remove user-selected labelled segments from a segmentation mask.
Args:
labeled_segments (ndarray): Labelled segmentation image.
remove_list (list[int] or None): Segment labels to remove. Use the segment numbers displayed in the visualisation for reference. If None, no segments will be removed.
Returns:
cleaned_mask (ndarray): Binary mask with curated segments
Returns None on failure.
"""
try:
mask = np.isin(labeled_segments, remove_list, invert=True) # Create mask to exclude selected segments
cleaned_mask = labeled_segments * mask # Apply mask to labeled segments
return cleaned_mask > 0 # Return manually curated binary mask
except Exception as e:
print(f"Error removing selected segments: {e}") # Print error message if removal fails
return None
def segment_and_analyze_vessels(data_path, file_name, image, gray_image, clahe_image, otsu_mask, output_path,
tile_index, im_filter=seg_method, sigma1=sigma1, hole_size=hole_size,
ditzle_size=ditzle_size, thresh=thresh):
"""
Vessel segmentation, optional apply of manual curation and results visualisation
Args:
data_path (str): Directory containing the input image (used for labeling).
file_name (str): File name used to derive the sample ID.
image (ndarray): Raw image data used for segmentation.
gray_image (ndarray): Greyscale MIP image.
clahe_image (ndarray): CLAHE-enhanced greyscale image.
otsu_mask (ndarray): Binary mask from Otsu thresholding.
output_path (str): Directory for saving figures.
tile_index (int): Index of the tile being processed (used for labelling and individual tile vessel metrics analysis).
im_filter (str): Vessel enhancement filter name (meijering, frangi, sato or jerman).
sigma1 (range): Sigma range for vessel enhancement.
hole_size (int): Max hole size to fill in the mask.
ditzle_size (int): Max size of small objects to remove.
thresh (int): Threshold for binarisation of the enhanced image.
Returns:
tuple:
vessel_seg (ndarray): Final vessel segmentation mask (curated if manual curation is applied).
remove_list (list[int] or None): Labels removed by the user.
Returns None on failure.
"""
try:
# Segment vessels using specified filtering method and parameters
vessel_seg = vm.segment_image(image,
im_filter=im_filter, # Vessel segmentation method
sigma1=sigma1, # Sigma range for vessel enhancement filter
hole_size=hole_size, # Hole size to be filled in the vessel mask
ditzle_size=ditzle_size, # Size of small objects to be removed from the vessel mask
thresh=thresh) # Threshold value to binarise the vessel enhanced image
seg_im = vessel_seg.astype(np.uint8) # Convert binary mask to uint8 format
# Visualise initial segmentation results
fig, axes = plt.subplots(1, 4, figsize=(15, 5)) # Create figure with 4 panels
axes[0].imshow(gray_image, cmap="gray") # Display MIP
axes[0].set_title("Maximum Intensity Projection (MIP)")
axes[0].axis("off")
axes[1].imshow(clahe_image, cmap="gray") # Display CLAHE processed image
axes[1].set_title("Clahe processed image")
axes[1].axis("off")
axes[2].imshow(seg_im, cmap="gray") # Display vessel segmentation image
axes[2].set_title("Vessel segmentation image")
axes[2].axis("off")
axes[3].imshow(otsu_mask, cmap="gray") # Display Otsu's threshold mask
axes[3].set_title("Otsu’s Threshold Mask")
axes[3].axis("off")
plt.show() # Show the figure
sampleID = file_name.split('_MaxInt')[0]
fig_title = f"{sampleID}_tile{tile_index}_segmentation_mask_prefiltering" # Title for the figure
output_pdf_path = os.path.join(output_path, f"{fig_title}.pdf") # Define output PDF path
save_plot_as_pdf(fig, output_pdf_path) # Save figure as PDF
# Visualise segmented regions with labels for manual curation
labeled_segments = visualize_segmented_regions(gray_image, seg_im, cmap_name='hsv', min_vessel_size = min_vessel_size) # Visualise segmented regions
sampleID = file_name.split('_MaxInt')[0]
fig_title = f"{sampleID}_tile{tile_index}_segmentation_labelled" # Title for the figure
output_pdf_path = os.path.join(output_path, f"{fig_title}.pdf") # Define output PDF path
save_plot_as_pdf(fig, output_pdf_path) # Save figure as PDF
# Manual curation: Remove unwanted segments based on user input
remove_list = input("Enter segment numbers to remove (comma-separated) or press Enter to skip: ").strip() # Get user input for segments to remove
if remove_list:
remove_list = [int(x) for x in remove_list.split(",")] # Convert input string to list of integers
cleaned_mask = remove_selected_segments(labeled_segments, remove_list) # Remove selected segments
labeled_segments = label(cleaned_mask) # Relabel the manually curated mask
vessel_seg = labeled_segments # Update vessel segmentation mask
else:
remove_list = None # No segments removed
cleaned_mask = remove_selected_segments(labeled_segments, remove_list) # Keep original mask
labeled_segments = label(cleaned_mask)
vessel_seg = labeled_segments
# Visualise manually curated segmentation results
fig, axes = plt.subplots(1, 4, figsize=(15, 5)) # Create figure with 4 panels
axes[0].imshow(gray_image, cmap="gray")
axes[0].set_title("Maximum Intensity Projection (MIP)") # Display MIP
axes[0].axis("off")
axes[1].imshow(clahe_image, cmap="gray") # Display CLAHE processed image
axes[1].set_title("Clahe processed image")
axes[1].axis("off")
axes[2].imshow(cleaned_mask, cmap="gray") # Display manually filtered vessel segmentation image
axes[2].set_title("Filtered vessel segmentation image")
axes[2].axis("off")
axes[3].imshow(otsu_mask, cmap="gray") # Display Otsu's threshold mask
axes[3].set_title("Otsu’s Threshold Mask")
axes[3].axis("off")
plt.show() # Show the figure
sampleID = file_name.split('_MaxInt')[0]
fig_title = f"{sampleID}_tile{tile_index}_segmentation_mask" # Title for the figure
output_pdf_path = os.path.join(output_path, f"{fig_title}.pdf") # Define output PDF path
save_plot_as_pdf(fig, output_pdf_path) # Save figure as PDF
return vessel_seg, remove_list
except Exception as e:
print(f"Error segmenting vessels for {file_name}: {e}") # Print error message if segmentation fails
return None
def calculate_and_print_metrics(label_path, vessel_seg):
"""
Compute vessel metrics and overlap statistics.
Args:
label_path (str): Path to the labelled image.
vessel_seg (ndarray): Vessel segmentation mask.
Returns:
tuple:
length (float): Vessel length metric.
area (float): Vessel area metric.
conn (float): Connectivity metric.
Q (float): Combined metric of Q = connectivity * length * area.
jaccard (float): Jaccard index between label and segmentation.
label (ndarray): Raw labelled image.
label_im (ndarray): Labelled image cast to uint8.
Returns None on failure.
"""
try:
label = cv2.imread(label_path, 0) # Load label image
label_im = label.astype(np.uint8) # Convert label image to uint8
seg_im = vessel_seg.astype(np.uint8) # Convert segmentation mask to uint8, [i] is the index of each ROI
length, area, conn, Q = vm.cal(label_im, seg_im) # Calculate vessel metrics: length, area, connectivity and Q factor (combined metric of connectivity * length * area)
jaccard = vm.jaccard(label_im, seg_im) # Calculate Jaccard index between label and segmentation mask
return length, area, conn, Q, jaccard, label, label_im # Return vessel metrics and labelled images
except Exception as e:
print(f"Error calculating metrics: {e}") # Print error message if vessel metrics calculation fails
return None
def find_branchpoints(skel, min_distance=Branch_dist):
"""
Detect branchpoints from a skeleton image and removes nearby branchpoints based on a minimum distance threshold.
Args:
skel (ndarray): Skeletonised vessel image.
min_distance (int): Minimum distance between branchpoints (pixels). Users are reccommended to set this parameter based median(vessel_length)
Returns:
tuple:
edges (ndarray): Skeleton edges after removing branchpoints.
branchpoints (ndarray): Binary image of cleaned branchpoints.
"""
skel_binary = np.zeros_like(skel) # Create binary skeleton image
skel_binary[skel > 0] = 1 # Binarise skeleton image
skel_index = np.argwhere(skel_binary == True) # Get coordinates of skeleton pixels
tile_sum = [] # List to store sum of neighbourhood pixels
neighborhood_image = np.zeros(skel.shape) # Image to store neighbourhood sums
for i, j in skel_index: # Loop through each skeleton pixel
this_tile = skel_binary[i - 1 : i + 2, j - 1 : j + 2] # Extract 3x3 neighbourhood
tile_sum.append(np.sum(this_tile)) # Sum of neighbourhood pixels
neighborhood_image[i, j] = np.sum(this_tile) # Store sum in neighbourhood image
branch_points = np.zeros_like(neighborhood_image) # Create branchpoint image
branch_points[neighborhood_image > 3] = 1 # Identify branchpoints (neighbourhood sum > 3)
branch_points = branch_points.astype(np.uint8) # Convert branchpoint image to uint8
branch_coords = np.argwhere(branch_points == 1) # Get coordinates of branchpoints
to_keep = np.ones(len(branch_coords), dtype=bool) # Boolean array to track branchpoints to keep
dists = distance.squareform(distance.pdist(branch_coords)) # Calculate pairwise distances between branchpoints
for i in range(len(branch_coords)): # Loop through each branchpoint
if not to_keep[i]: # Skip if already marked for removal
continue
close = np.where((dists[i] < min_distance) & (dists[i] > 0))[0] # Find nearby branchpoints within min_distance
to_keep[close] = False # Keep only the first point and remove closeby branchpoints
cleaned_branch_points = np.zeros_like(skel, dtype=np.uint8) # Creates empty image to store filtered branchpoints
for idx in np.where(to_keep)[0]: # Loop through the index of kept branchpoints
y, x = branch_coords[idx] # Get coordinates for each branchpoint
cleaned_branch_points[y, x] = 1 # Mark the location of the branchpoint
edges = skel_binary - cleaned_branch_points # Edges are defined as skeleton pixels minus the cleaned branchpoints
edges[edges < 0] = 0 # Ensures no negative values are stored
branchpoints = cleaned_branch_points # Store filtered branchpoints
return edges.astype(np.uint8), branchpoints
def analyze_skeleton(vessel_seg, gray_image, file_name, tile_index):
"""
Analyse skeleton and vessel diameter. Visualise skeleton and branchpoints and save figures as PDF.
Args:
vessel_seg (ndarray): Vessel segmentation mask.
gray_image (ndarray): Greyscale MIP image.
file_name (str): File name used to derive the sample ID.
tile_index (int): Index of the tile being processed (used for labelling and individual tile vessel metrics analysis).
Returns:
tuple:
skel (ndarray): Skeletonised image.
edges (ndarray): Skeleton edges with branchpoints removed.
edge_labels (ndarray): Labelled edges image.
branchpoints (ndarray): Binary branchpoint image.
coords (list): Edge coordinates from endpoint analysis.
endpoints (list): Endpoint coordinates from endpoint analysis.
Returns None on failure.
"""
try:
skel, edges, _ = vm.skeletonize_vm(vessel_seg, method=skel_method) # Calculate skeleton metrics and edges of each image, [i] is the index of each ROI
_, edge_labels = cv2.connectedComponents(edges) # Label connected edges with a unique identifier
_, branchpoints = find_branchpoints(skel) # Define branhpoints (refined from find_branchpoints)
coords, endpoints = vm.find_endpoints(edges) # Find coordinates of each edge for further analysis and visualisation
# Visualise skeleton plot
cmap1 = matplotlib.colors.LinearSegmentedColormap.from_list("", [(1, 1, 1, 0), (0, 0, 0, 1)]) # Colourmap for skeleton (transparent to black)
cmap2 = matplotlib.colors.LinearSegmentedColormap.from_list("", [(1, 1, 1, 0), (1, 0, 0, 1)]) # Colourmap for overlay (transparent to red)
fig, axes = plt.subplots(1, 3, figsize=(15, 5)) # Create a figure with 3 panels
axes[0].imshow(gray_image, cmap="gray") # Display MIP
axes[0].set_title("Maximum Intensity Projection (MIP)")
axes[0].axis("off")
axes[1].imshow(skel, cmap=cmap1) # Display skeleton
axes[1].set_title("Skeleton Plot")
axes[1].axis("off")
axes[2].imshow(gray_image, cmap="gray") # Display overlay of skeleton over MIP
axes[2].imshow(skel, cmap=cmap2)
axes[2].set_title("Overlap")
axes[2].axis("off")
plt.show() # Show the figure
sampleID = file_name.split('_MaxInt')[0]
fig_title = f"{sampleID}_tile{tile_index}_skeleton_plot" # Title for the figure
output_pdf_path = os.path.join(output_path, f"{fig_title}.pdf") # Define output PDF path
save_plot_as_pdf(fig, output_pdf_path) # Save figure as PDF
# Visualise skeleton lines
kernel = np.ones((10, 10), np.uint8) # Adjust kernel size to make branchpoints bigger
branchpoints_dilated = cv2.dilate(branchpoints.astype(np.uint8), kernel, iterations=1) # Dilate branchpoints to make them more visible on the image (only for visualisation purposes)
fig, axes = plt.subplots(1, 3, figsize=(15, 5)) # Create a figure with 3 panels
axes[0].imshow(gray_image, cmap="gray") # Display MIP
axes[0].set_title("Maximum Intensity Projection (MIP)")
axes[0].axis("off")
axes[1].imshow(branchpoints_dilated, cmap='gray_r') # Display branchpoints
axes[1].set_title("Branchpoints")
axes[1].axis("off")
axes[2].imshow(gray_image, cmap="gray") # Display overlay of branchpoints over MIP
axes[2].imshow(branchpoints_dilated, cmap=cmap2)
axes[2].set_title("Overlap")
axes[2].axis("off")
plt.show() # Show the figure
sampleID = file_name.split('_MaxInt')[0]
fig_title = f"{sampleID}_tile{tile_index}_branchpoint_plot" # Title for the figure
output_pdf_path = os.path.join(output_path, f"{fig_title}.pdf") # Define output PDF path
save_plot_as_pdf(fig, output_pdf_path) # Save figure as PDF
return skel, edges, edge_labels, branchpoints, coords, endpoints
except Exception as e:
print(f"Error analyzing skeleton: {e}") # Print error message if skeletonisation fails
return None
def is_outlier(series, threshold_low=threshold_low, threshold_high=threshold_high):
"""
Identify outliers using Median Absolute Deviation (MAD).
Args:
series (pandas.Series): Series of numeric values to calculate central tendency.
threshold_low (float): Lower MAD threshold.
threshold_high (float): Upper MAD threshold.
Returns:
pandas.Series: Boolean mask where True indicates an outlier.
"""
median = series.median() # Calculates central tendency using median
mad = median_abs_deviation(series) # Calculates non-scaled MAD dispersion
# Compute assymetric bounds using provided thresholds
lower_bound = median - threshold_low * mad # Lower threshold
upper_bound = median + threshold_high * mad # Upper threshold
return (series < lower_bound) | (series > upper_bound) # Return a boolean with outliers
def DF_filtering(df):
"""
Filter outliers from each column of a DataFrame using MAD.
Args:
df (pandas.DataFrame): DataFrame to filter column-wise based on MAD outlier analysis.
Returns:
pandas.DataFrame: DataFrame with outliers removed per column.
"""
filtered_dict = {} # Empty dictionary to store filtered dataframes
for column in df.columns: # Iterate through each column
data = df[column].dropna() # Removes missing values to avoid propagating NaNs into MAD and comparisons
outliers = is_outlier(data) # Compute outlier mask for the current column
filtered_data = data[~outliers] # Keep only non-outlier values
filtered_dict[column] = filtered_data.reset_index(drop=True) # Reset index and store filtered data in the dictionary
filtered_values_df = pd.DataFrame(filtered_dict) # Build new dataframe with filtered results
return filtered_values_df
def analyze_vessel_diameters(image, vessel_seg, edge_labels, gray_image, file_name, tile_index):
"""
Compute vessel diameters and visualise diameter crosslines.
Args:
image (ndarray): Raw image data used for diameter extraction.
vessel_seg (ndarray): Vessel segmentation mask.
edge_labels (ndarray): Labelled edges image from skeletonisation.
gray_image (ndarray): Greyscale MIP image.
file_name (str): File name used to derive the sample ID.
tile_index (int): Index of the tile being processed (used for labelling and individual tile vessel metrics analysis).
Returns:
tuple:
diam_values (list[float]): Diameter values in pixels.
viz (ndarray): Visualisation image of diameter crosslines.
filtered_diam_values (pandas.DataFrame): MAD-filtered diameters.
Returns ([], None, None) on failure.
"""
try:
viz, diameters = vm.whole_anatomy_diameter(image, vessel_seg, edge_labels) # Extract diameter measurements from the entire ROI
cmap2 = matplotlib.colors.LinearSegmentedColormap.from_list("", [(1, 1, 1, 0), (1, 0, 0, 1)]) # Define colourmap for overlay (transparent to red)
fig, axes = plt.subplots(1, 3, figsize=(15, 5)) # Figure with 3 panels
axes[0].imshow(gray_image, cmap="gray") # Display MIP
axes[0].set_title("Maximum Intensity Projection (MIP)")
axes[0].axis("off")
axes[1].imshow(viz, cmap="gray_r") # Display diameter crosslines
axes[1].set_title("Diameters")
axes[1].axis("off")
axes[2].imshow(gray_image, cmap="gray") # Display overlay of diameter crosslines over MIP
axes[2].imshow(viz, cmap=cmap2)
axes[2].set_title("Overlap")
axes[2].axis("off")
plt.show() # Show the figure
sampleID = file_name.split('.tiff')[0]
fig_title = f"{sampleID}_tile{tile_index}_diameter_plot" # Title for the figure
output_pdf_path = os.path.join(output_path, f"{fig_title}.pdf") # Define output PDF path
save_plot_as_pdf(fig, output_pdf_path) # Save figure as PDF
num_diameters = diameters[-1][0] # Re-index diameters with unique identifiers after previous diameter removal
indexed_diameters = [] # Empty list to store new indexed diameters
for new_index in range(1, num_diameters + 1): # Loop through each diameter identifier
found = next((d for d in diameters if d[0] == new_index), None) # Find unique identifiers
if found:
indexed_diameters.append((new_index, found[1])) # If vessel ID was found, keep identifier
filtered_diameters = [(idx, diam) for idx, diam in indexed_diameters if diam > 0] # Remove crosslines that present diameters with zero values
diam_values = [diam for idx, diam in filtered_diameters] # Extract only positive diameter values
filtered_diam_values = DF_filtering(pd.DataFrame(diam_values, columns=['Diameter (px)'])) # MAD-based filtering of outliers
return diam_values, viz, filtered_diam_values # Return the diameter values
except Exception as e:
print(f"Error analyzing vessel diameters: {e}") # Print error message if diameter analysis fails
return [], None, None # Return an empty list in case of an error
def visualize_vessel_metrics(edge_labels, diam_values, scale = scale):
"""
Calculate vessel length, tortuosity, and diameter metrics and performs MAD filtering.
Args:
edge_labels (ndarray): Labelled edges image from skeletonisation.
diam_values (list[float]): Vessel diameter values in pixels.
scale (float): Pixels per micron for spatial conversion. Users need to set this parameter based on the imaging resolution of their dataset.
Returns:
tuple:
length (ndarray): Segment lengths in pixels.
vessel_len (float): Mean segment length in pixels.
tort (ndarray): Tortuosity values per segment.
vessel_tor (float): Mean tortuosity.
diam_values_um (ndarray): Diameters converted to microns.
length_um (ndarray): Segment lengths converted to microns.
filtered_length_df (pandas.DataFrame): MAD-filtered lengths (pixels).
filtered_length_um_df (pandas.DataFrame): MAD-filtered lengths (microns).
filtered_diam_values_df (pandas.DataFrame): MAD-filtered diameters (pixels).
filtered_diam_values_um_df (pandas.DataFrame): MAD-filtered diameters (microns).
Returns None on failure.
"""
try:
_, length = vm.vessel_length(edge_labels) # Calculate the length of each segment (pixels)
vessel_len = np.mean(length) # Compute the mean of segment lengths
tort, _ = vm.tortuosity(edge_labels) # Calculate tortuosity values
vessel_tor = np.mean(tort) # Compute the mean tortuosity
tort = np.array(tort) # Transform values into array
# Spatial unit convertion: pixels to microns
diam_values_um = np.array(diam_values) / scale # Convert diameter values
length_um = length / scale # Convert segment length values
filtered_length_df = DF_filtering(pd.DataFrame(length)) # Filter outliers for segment length in pixels
filtered_length_um_df = DF_filtering(pd.DataFrame(length_um)) # Filter outliers for segment length in microns
filtered_diam_values_df = DF_filtering(pd.DataFrame(diam_values)) # Filter outliers for vessel diameters in pixels
filtered_diam_values_um_df = DF_filtering(pd.DataFrame(diam_values_um)) # Filter outliers for vessel diameters in microns
return length, vessel_len, tort, vessel_tor, diam_values_um, length_um, filtered_length_df, filtered_length_um_df, filtered_diam_values_df, filtered_diam_values_um_df
except Exception as e:
print(f"Error visualizing vessel metrics: {e}") # Print error message if vessel measurement fails
return None
def analyze_vessel_network(image, vessel_seg, edges, label, scale = scale, window = window_size):
"""
Compute network-level metrics such as total length and density.
Args:
image (ndarray): Raw image data.
vessel_seg (ndarray): Vessel segmentation mask.
edges (ndarray): Skeleton edges image.
label (ndarray): Labelled image for branchpoint density.
scale (float): Pixels per micron for spatial conversion.
window (int): Window size (microns) for density binning.
Returns:
tuple:
net_length (float): Total network length in pixels.
net_length_um (float): Total network length in microns.
density (float): Mean vessel density.
density_array (ndarray): Bin density values.
overlay (ndarray): Density overlay image.
vessel_density (float): Fraction of vessel in micron².
bp_density (tuple): Branchpoint density outputs.
Returns None on failure.
"""
density_scale = np.floor(scale * window).astype(int) # Convert pixel scale to a 10×10 grid corresponding to 100µm² bins
try:
net_length = vm.network_length(edges) # Calculate total network length (sum of vessel segments)
net_length_um = net_length / scale # Convert total network length into microns
density, density_array, overlay = vm.vessel_density(image, vessel_seg, density_scale, density_scale) # Compute vessel density map
vessel_density = sum(density_array) # Number of vessel pixels versus total pixels
bp_density = vm.branchpoint_density(label) # Optional: Calculates branchpoint density
return net_length, net_length_um, density, density_array, overlay, vessel_density, bp_density
except Exception as e:
print(f"Error analyzing vessel network: {e}") # Print error message if network measurements fail
def save_results_to_excel(results, diam_results):
"""
Save vessel analysis metrics and raw values to an Excel file.
Args:
results (list[list]): Summary metrics per sample.
diam_results (dict): Raw and filtered diameters and length results (pixels and microns).
Returns:
None
"""
try:
# Sheet 1: Summary embryo metrics
embryo_metrics = pd.DataFrame(results) # Convert results into a dataframe
# Assign columns names to the summary results sheet
embryo_metrics.columns = ["File name", "Area", "Connectivity", "Q", "Jaccard",
"Mean raw vessel length (px)", "Mean raw vessel Length (um)",
"Mean vessel length (px)", "Mean vessel Length (um)",
"Mean raw vessel diameter (px)", "Mean raw vessel diameter (um)",
"Mean vessel diameter (px)", "Mean vessel diameter (um)",
"Vessel Tortuosity", "Network Length (px)", "Network Length (um)",
"Vessel_density", "Branchpoints"]
# Sheet 2: Raw and filtered diameter values (px)
diam_values_df = pd.DataFrame(diam_results["diam_values"]).transpose() # Convert entries into dataframe and transpose so each row corresponds to an image
diam_values_df.columns = [f"Diameter {i+1}" for i in range(diam_values_df.shape[1])] # Rename columns with diameter identifiers for each embryo
# Same for MAD-filtered diameter values
filtered_diam_values_df = pd.DataFrame(diam_results["filtered_diam_values"]).transpose()
filtered_diam_values_df.columns = [f"Diameter {i+1}" for i in range(filtered_diam_values_df.shape[1])]
# Sheet 3: Raw and filtered diameter values (µm)
diam_values_um_df = pd.DataFrame(diam_results["diam_values_um"]).transpose()
diam_values_um_df.columns = [f"Diameter {i+1} (um)" for i in range(diam_values_um_df.shape[1])]
filtered_diam_values_um_df = pd.DataFrame(diam_results["filtered_diam_values_um"]).transpose()
filtered_diam_values_um_df.columns = [f"Diameter {i+1}" for i in range(filtered_diam_values_um_df.shape[1])]
# Sheet 4: Vessel Lengths (px)
vessel_lengths_df = pd.DataFrame(diam_results["length"]).transpose()
vessel_lengths_df.columns = [f"Length {i+1}" for i in range(vessel_lengths_df.shape[1])]
filtered_vessel_lengths_df = pd.DataFrame(diam_results["filtered_length"]).transpose()
filtered_vessel_lengths_df.columns = [f"Length {i+1}" for i in range(filtered_vessel_lengths_df.shape[1])]
# Sheet 5: Vessel Lengths (µm)
vessel_lengths_um_df = pd.DataFrame(diam_results["length_um"]).transpose()
vessel_lengths_um_df.columns = [f"Length {i+1}" for i in range(vessel_lengths_um_df.shape[1])]
filtered_vessel_lengths_um_df = pd.DataFrame(diam_results["filtered_length_um"]).transpose()
filtered_vessel_lengths_um_df.columns = [f"Length {i+1}" for i in range(filtered_vessel_lengths_um_df.shape[1])]
# Save results to Excel
output_excel_path = os.path.join(output_path, f"{output_name}.xlsx") # Construct the output path under within the working directory
with pd.ExcelWriter(output_excel_path) as writer:
# Write each DataFrame to a separate sheet
embryo_metrics.to_excel(writer, sheet_name="Embryo Metrics", index=False)
diam_values_df.to_excel(writer, sheet_name="Diameter Values", index=False)
filtered_diam_values_df.to_excel(writer, sheet_name="MAD filtered diameter Values", index=False)
diam_values_um_df.to_excel(writer, sheet_name="Diameter Values (um)", index=False)
filtered_diam_values_um_df.to_excel(writer, sheet_name="MAD filtered diameter Values (um)", index=False)
vessel_lengths_df.to_excel(writer, sheet_name="Vessel Lengths", index=False)
filtered_vessel_lengths_df.to_excel(writer, sheet_name="MAD filtered vessel lengths", index=False)
vessel_lengths_um_df.to_excel(writer, sheet_name="Vessel Lengths (um)", index=False)
filtered_vessel_lengths_um_df.to_excel(writer, sheet_name="MAD filtered vessel lengths (um)", index=False)
print(f"Results saved to {output_excel_path}")
except Exception as e:
print(f"Error saving results to Excel: {e}") # Print error message if saving results fail
def process_all_files(data_path, output_path):
"""
Process all single-channel images from the input directory.
Args:
data_path (str): Directory containing input images.
output_path (str): Directory for saving outputs.
Returns:
None
"""
file_names = [f for f in os.listdir(data_path) if f.startswith(name_start) and f.endswith(".czi")] # Find all images matching naming pattern
all_results = [] # Empty list to store final metrics per embryo
# Empty dictionary to store vessel metrics across samples
all_diam_results = {
"diam_values": [],
"filtered_diam_values": [],
"diam_values_um": [],
"filtered_diam_values_um": [],
"length": [],
"filtered_length": [],
"length_um": [],
"filtered_length_um": [],
}
for file_name in file_names: # Loop through each file found in directory
file_path = os.path.join(data_path, file_name) # Find the full file path for each image
sampleID = file_name.split('_MaxInt')[0] # Retrieve sample identifier
print(f"\n--- Processing {sampleID} ---")
with czifile.CziFile(file_path) as czi:
data = czi.asarray() # Read confocal image data from .czi file
data_squeezed = np.squeeze(data) # Squeeze data to remove unwanted dimensions/metadata and retrieve image data
n_tiles = tiles # Split lateral image by the number of horizontal tiles (decrease computational requirements and processing time)
tile_width = data_squeezed.shape[1] // n_tiles # Calculate the width of each tile
tile_metrics = [] # Empty list to store vessel measurements for each tile
# Empty dictionary to store vessel metrics across tiles
tile_diam_results = {
"diam_values": [],
"filtered_diam_values": [],
"diam_values_um": [],
"filtered_diam_values_um": [],
"length": [],
"filtered_length": [],
"length_um": [],
"filtered_length_um": [],
}
for i in range(n_tiles): # Loop through each tile
# Find tile boundaries
start_x = i * tile_width
end_x = (i + 1) * tile_width if i < n_tiles - 1 else data_squeezed.shape[1]
tile = data_squeezed[:, start_x:end_x] # Extract tile from full image
image = tile # Rename tile image for convenience
mip_image = cv2.normalize(image, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8) # Normalise image to uint8
mip_image_rgb = cv2.cvtColor(mip_image, cv2.COLOR_GRAY2RGB)
gray_image = cv2.cvtColor(mip_image_rgb, cv2.COLOR_RGB2GRAY)
clahe_image = clahe.apply(gray_image) # Apply CLAHE for contrast enhacement
_, otsu_mask = cv2.threshold(clahe_image, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU) # Otsu thresholding
output_mip_path = os.path.join(output_path, f"{sampleID}_tile{i+1}.tiff") # Define output MIP path
output_mask_path = os.path.join(output_path, f"{sampleID}_tile{i+1}_mask.tiff") # Define output mask path
display_and_save_mip_otsu(gray_image, clahe_image, otsu_mask, output_mip_path, output_mask_path)
vessel_seg, _ = segment_and_analyze_vessels(data_path, file_name, image, gray_image, clahe_image, otsu_mask, output_path, i+1,
im_filter=seg_method, sigma1=sigma1, hole_size=hole_size,
ditzle_size=ditzle_size, thresh=thresh
) # Vessel segmentation and topology analysis
if vessel_seg is None: # Skip tile if segmentation failed
continue
label_path = output_mask_path # Use mask path for metrics function
length, area, conn, Q, jaccard, label, _ = calculate_and_print_metrics(label_path, vessel_seg) # Calculate vessel metrics
skel, edges, edge_labels, _, _, _ = analyze_skeleton(vessel_seg, gray_image, file_name, i+1) # Perform skeletonisation and branchpoint analysis
diam_values, viz, filtered_diam_values = analyze_vessel_diameters(image, vessel_seg, edge_labels, gray_image, file_name, i+1) # Measure vessel diameters
# Calculate diameter and vessel length metrics, including spatial conversion
(length, vessel_len, _, vessel_tor, diam_values_um, length_um, filtered_length_df,
filtered_length_um_df, filtered_diam_values_df, filtered_diam_values_um_df) = visualize_vessel_metrics(edge_labels, diam_values, scale = scale)
# Calculate overall network metrics
net_length, net_length_um, _, _, _, vessel_density, bp_density = analyze_vessel_network(image, vessel_seg, edges, label, scale = scale, window = window_size)
# Store metric results for selected tile
tile_metrics.append([
area, conn, Q, jaccard, vessel_len,
np.mean(length_um), filtered_length_df.values.mean(), filtered_length_um_df.values.mean(),
np.mean(diam_values), np.mean(diam_values_um), filtered_diam_values_df.values.mean(), filtered_diam_values_um_df.values.mean(),
vessel_tor, net_length, net_length_um, vessel_density, len(bp_density[0])
])
# Store raw metrics for further saving
tile_diam_results["diam_values"].append(diam_values)
tile_diam_results["filtered_diam_values"].append(filtered_diam_values.values.flatten().tolist() if filtered_diam_values is not None else [])
tile_diam_results["diam_values_um"].append(diam_values_um)
tile_diam_results["filtered_diam_values_um"].append(filtered_diam_values_um_df.values.flatten().tolist() if filtered_diam_values_um_df is not None else [])
tile_diam_results["length"].append(length)
tile_diam_results["filtered_length"].append(filtered_length_df.values.flatten().tolist() if filtered_length_df is not None else [])
tile_diam_results["length_um"].append(length_um)
tile_diam_results["filtered_length_um"].append(filtered_length_um_df.values.flatten().tolist() if filtered_length_um_df is not None else [])
if tile_metrics: # Ensure all tile metrics were processed successfully
tile_metrics_np = np.array(tile_metrics, dtype=np.float64) # Convert arrays for averaging
avg_metrics = np.nanmean(tile_metrics_np, axis=0) # Average metrics across tiles
sum_metrics = np.nansum(tile_metrics_np[:, 13:17], axis=0) # Sum network metrics
avg_metrics[13:17] = sum_metrics # Replace averaged entries with summed values
all_results.append([sampleID] + avg_metrics.tolist()) # Store final results per embryo
per_tile_results = [] # Save metrics per tile
for idx, metrics in enumerate(tile_metrics):
per_tile_results.append([sampleID, f"Tile {idx+1}"] + metrics)
# Store all tile-specific metrics
if 'all_tile_results' not in locals():
all_tile_results = []
all_tile_results.extend(per_tile_results)
# Flatten lists so each sample has 1 row per metric type
for key in all_diam_results.keys():
flat = [item for sublist in tile_diam_results[key] for item in (sublist if isinstance(sublist, list) else [sublist])]
all_diam_results[key].append(flat)
save_results_to_excel(all_results, all_diam_results) # Save summary results to Excel
try: # Save per-tile results to Excel
per_tile_df = pd.DataFrame(all_tile_results)
per_tile_df.columns = ["File name", "Tile", "Area", "Connectivity", "Q", "Jaccard",
"Mean raw vessel length (px)", "Mean raw vessel Length (um)",
"Mean vessel length (px)", "Mean vessel Length (um)",
"Mean raw vessel diameter (px)", "Mean raw vessel diameter (um)",
"Mean vessel diameter (px)", "Mean vessel diameter (um)",
"Vessel Tortuosity", "Network Length (px)", "Network Length (um)",
"Vessel_density", "Branchpoints"]
output_tile_excel_path = os.path.join(output_path, f"{output_tile_name}.xlsx")
per_tile_df.to_excel(output_tile_excel_path, index=False)
print(f"Per-tile results saved to {output_tile_excel_path}")
except Exception as e:
print(f"Error saving per-tile results to Excel: {e}") # Print error message if tile-specific results were not saved successfully
def process_all_files_double(data_path, output_path):
"""
Process all multi-channel images from the input directory.
Args:
data_path (str): Directory containing input images.
output_path (str): Directory for saving outputs.
Returns:
None
"""
file_names = [f for f in os.listdir(data_path) if f.startswith(name_start) and f.endswith(".czi")] # Find all images matching naming pattern
all_results = [] # Empty list to store final metrics per embryo
# Empty dictionary to store vessel metrics across samples
all_diam_results = {
"diam_values": [],
"filtered_diam_values": [],
"diam_values_um": [],
"filtered_diam_values_um": [],
"length": [],
"filtered_length": [],
"length_um": [],
"filtered_length_um": [],
}
for file_name in file_names[0:2]: # Loop through each file found in directory
file_path = os.path.join(data_path, file_name) # Find the full file path for each image
sampleID = file_name.split('_MaxInt')[0] # Retrieve sample identifier
print(f"\n--- Processing {sampleID} ---")
with czifile.CziFile(file_path) as czi:
data = czi.asarray() # Read confocal image data from .czi file
data_squeezed = np.squeeze(data)
data_squeezed = data_squeezed[channel]
n_tiles = tiles # Split lateral image by the number of horizontal tiles (decrease computational requirements and processing time)
tile_width = data_squeezed.shape[1] // n_tiles # Calculate the width of each tile
tile_metrics = [] # Empty list to store vessel measurements for each tile
# Empty dictionary to store vessel metrics across tiles
tile_diam_results = {
"diam_values": [],
"filtered_diam_values": [],
"diam_values_um": [],
"filtered_diam_values_um": [],
"length": [],
"filtered_length": [],
"length_um": [],
"filtered_length_um": [],
}
for i in range(n_tiles): # Loop through each tile
# Find tile boundaries
start_x = i * tile_width
end_x = (i + 1) * tile_width if i < n_tiles - 1 else data_squeezed.shape[1]
tile = data_squeezed[:, start_x:end_x] # Extract tile from full image
image = tile # Rename tile image for convenience
mip_image = cv2.normalize(image, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8) # Normalise image to uint8
mip_image_rgb = cv2.cvtColor(mip_image, cv2.COLOR_GRAY2RGB)
gray_image = cv2.cvtColor(mip_image_rgb, cv2.COLOR_RGB2GRAY)
clahe_image = clahe.apply(gray_image) # Apply CLAHE for contrast enhacement
_, otsu_mask = cv2.threshold(clahe_image, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU) # Otsu thresholding
output_mip_path = os.path.join(output_path, f"{sampleID}_tile{i+1}.tiff") # Define output MIP path
output_mask_path = os.path.join(output_path, f"{sampleID}_tile{i+1}_mask.tiff") # Define output mask path
display_and_save_mip_otsu(gray_image, clahe_image, otsu_mask, output_mip_path, output_mask_path)
vessel_seg, _ = segment_and_analyze_vessels(data_path, file_name, image, gray_image, clahe_image, otsu_mask, output_path, i+1,
im_filter=seg_method, sigma1=sigma1, hole_size=hole_size,
ditzle_size=ditzle_size, thresh=thresh
) # Vessel segmentation and topology analysis
if vessel_seg is None: # Skip tile if segmentation failed
continue
label_path = output_mask_path # Use mask path for metrics function
length, area, conn, Q, jaccard, label, _ = calculate_and_print_metrics(label_path, vessel_seg) # Calculate vessel metrics
skel, edges, edge_labels, _, _, _ = analyze_skeleton(vessel_seg, gray_image, file_name, i+1) # Perform skeletonisation and branchpoint analysis
diam_values, viz, filtered_diam_values = analyze_vessel_diameters(image, vessel_seg, edge_labels, gray_image, file_name, i+1) # Measure vessel diameters
# Calculate diameter and vessel length metrics, including spatial conversion
(length, vessel_len, _, vessel_tor, diam_values_um, length_um, filtered_length_df,
filtered_length_um_df, filtered_diam_values_df, filtered_diam_values_um_df) = visualize_vessel_metrics(edge_labels, diam_values, scale = scale)
# Calculate overall network metrics
net_length, net_length_um, _, _, _, vessel_density, bp_density = analyze_vessel_network(image, vessel_seg, edges, label, scale = scale, window = window_size)
# Store metric results for selected tile
tile_metrics.append([
area, conn, Q, jaccard, vessel_len,
np.mean(length_um), filtered_length_df.values.mean(), filtered_length_um_df.values.mean(),
np.mean(diam_values), np.mean(diam_values_um), filtered_diam_values_df.values.mean(), filtered_diam_values_um_df.values.mean(),
vessel_tor, net_length, net_length_um, vessel_density, len(bp_density[0])
])
# Store raw metrics for further saving
tile_diam_results["diam_values"].append(diam_values)
tile_diam_results["filtered_diam_values"].append(filtered_diam_values.values.flatten().tolist() if filtered_diam_values is not None else [])
tile_diam_results["diam_values_um"].append(diam_values_um)
tile_diam_results["filtered_diam_values_um"].append(filtered_diam_values_um_df.values.flatten().tolist() if filtered_diam_values_um_df is not None else [])
tile_diam_results["length"].append(length)
tile_diam_results["filtered_length"].append(filtered_length_df.values.flatten().tolist() if filtered_length_df is not None else [])
tile_diam_results["length_um"].append(length_um)
tile_diam_results["filtered_length_um"].append(filtered_length_um_df.values.flatten().tolist() if filtered_length_um_df is not None else [])
if tile_metrics: # Ensure all tile metrics were processed successfully
tile_metrics_np = np.array(tile_metrics, dtype=np.float64) # Convert arrays for averaging
avg_metrics = np.nanmean(tile_metrics_np, axis=0) # Average metrics across tiles
sum_metrics = np.nansum(tile_metrics_np[:, 13:17], axis=0) # Sum network metrics
avg_metrics[13:17] = sum_metrics # Replace averaged entries with summed values
all_results.append([sampleID] + avg_metrics.tolist()) # Store final results per embryo
per_tile_results = [] # Save metrics per tile
for idx, metrics in enumerate(tile_metrics):
per_tile_results.append([sampleID, f"Tile {idx+1}"] + metrics)
# Store all tile-specific metrics
if 'all_tile_results' not in locals():
all_tile_results = []
all_tile_results.extend(per_tile_results)
# Flatten lists so each sample has 1 row per metric type
for key in all_diam_results.keys():
flat = [item for sublist in tile_diam_results[key] for item in (sublist if isinstance(sublist, list) else [sublist])]
all_diam_results[key].append(flat)
save_results_to_excel(all_results, all_diam_results) # Save summary results to Excel
try: # Save per-tile results to Excel
per_tile_df = pd.DataFrame(all_tile_results)
per_tile_df.columns = ["File name", "Tile", "Area", "Connectivity", "Q", "Jaccard",
"Mean raw vessel length (px)", "Mean raw vessel Length (um)",
"Mean vessel length (px)", "Mean vessel Length (um)",
"Mean raw vessel diameter (px)", "Mean raw vessel diameter (um)",
"Mean vessel diameter (px)", "Mean vessel diameter (um)",
"Vessel Tortuosity", "Network Length (px)", "Network Length (um)",
"Vessel_density", "Branchpoints"]
output_tile_excel_path = os.path.join(output_path, f"{output_tile_name}.xlsx")
per_tile_df.to_excel(output_tile_excel_path, index=False)
print(f"Per-tile results saved to {output_tile_excel_path}")
except Exception as e:
print(f"Error saving per-tile results to Excel: {e}") # Print error message if tile-specific results were not saved successfully
# %% [markdown]
# ## <a id='toc3_3_'></a>[Run samples - Single channel](#toc0_)
# %%
if not os.path.exists(output_path):
os.makedirs(output_path)
process_all_files(data_path, output_path)
# %% [markdown]
# ## <a id='toc3_4_'></a>[Run samples - Two or more channels](#toc0_)
# %%
if not os.path.exists(output_path):
os.makedirs(output_path)
process_all_files_double(data_path, output_path)

VISTA-Z_lateral.ipynb at commit 6dd30c4, under MIT · at the source

Overview

Authors: Ignacio Rodriguez-Pastrana¹, Joanna Richens¹, Robert N Wilkinson¹

ORCID iDs: Ignacio Rodriguez-Pastrana

School of Life Sciences, University of Nottingham, Nottingham, NG7 2UH UK

Institutions: University of Nottingham (United Kingdom)

Journal: Scientific reports, volume 16, issue 1, article 15611

Dates: received 19 December 2025; accepted 3 March 2026; published online 1 April 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1038/s41598-026-43301-5 · PMID 41922411 · PMCID PMC13187487 · OpenAlex W7147264012

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: zebrafish (organism), methods / tools (subfield)

Methods: Statistics

Keywords: Zebrafish, Vascular imaging, Automated image analysis, Vessel segmentation, Angiogenesis, High-throughput phenotyping, Biological techniques, Biotechnology, Computational biology and bioinformatics

MeSH: Angiogenesis*, Blood Vessels*, Image Processing, Computer-Assisted*, Neovascularization, Physiologic*, Zebrafish*, Animals, Animals, Genetically Modified, Embryo, Nonmammalian, Forkhead Transcription Factors, Zebrafish Proteins (* major topic)

Topic: Zebrafish Biomedical Research Applications (Cell Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex

Funding: Medical Research Council (MR/X008215/1)

Citations: cited by 1 paper (Europe PMC); 57 references in the paper

Abstract

High-resolution vascular imaging in zebrafish embryos offers unparalleled insight into angiogenesis, yet quantitative analysis remains limited by manual workflows and inconsistent segmentation. To overcome this, we developed VISTA-Z (Vascular Imaging and Segmentation for Topology Analysis in Zebrafish), an automated Python-based pipeline designed to standardise vascular quantification across diverse fluorescent datasets. The workflow combines adaptive contrast enhancement, vessel segmentation using Meijering filtering, artefact removal through segment labelling, and skeleton-based topology analysis with refined branchpoint detection. Metrics are normalised to physical units and imaging depth, enabling reproducible comparisons across experiments. We validated VISTA-Z using multiple endothelial transgenic lines and developmental stages. The pipeline detected subtle and severe phenotypes, including brain vessel loss in foxc1a mutants, brain and trunk vessel loss in kdrl mutants and widespread hyper-angiogenesis in plxnd1 crispants. VISTA-Z is an open-source, scalable platform for reproducible high-throughput quantification of zebrafish vascular architecture, providing a standardised framework for developmental research and preclinical screening

Supplementary Information: The online version contains supplementary material available at 10.1038/s41598-026-43301-5.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 6 matches between paragraphs and lines of code.

fishyvessels/VISTA-Z

License: MIT

State: the link answers, verified on 28 September 2026

Evidence: files inventoried

Commit: 6dd30c47f0092c3becbecf6e8fca0d28e6a16bd5, 3 April 2026

Languages: Jupyter (6)

Size: 12 files, 6 scripts

Software Heritage: not archived

Found in: “Data availability”

Holds: README, license file, environment (Codes/requirements.txt), 6 notebooks

Not found: CITATION.cff, tests, continuous integration, documentation

Tools: NumPy (6 files), SciPy (6 files), Matplotlib (5 files), OpenCV (5 files), Plotly (5 files), scikit-image (5 files), imageio (4 files), pandas (4 files), seaborn (4 files), napari (1 file)

Availability: 1 check, the latest on 28 September 2026: the link answers

28 September 2026: the link answers

8 files

Codes/3D_volume_rendering.ipynb, Jupyter, 260 lines
Codes/Colocalisation_analysis.ipynb, Jupyter, 277 lines
Codes/VISTA-Z_ISV_loss.ipynb, Jupyter, 769 lines, 1 match
Codes/VISTA-Z_ROI.ipynb, Jupyter, 1,479 lines, 2 matches
Codes/VISTA-Z_dorsal.ipynb, Jupyter, 929 lines
Codes/VISTA-Z_lateral.ipynb, Jupyter, 1,054 lines, 3 matches
LICENSE, License, 21 lines
README.md, Text, 127 lines

At the source: github.com/fishyvessels/VISTA-Z

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
6 scripts, each with its path and the digest of its content;
6 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

biostudies:S-BIAD2914, at BioStudies; found in “Data availability”

Data availability

Code is available via GitHub (https://github.com/fishyvessels/VISTA-Z). Test image datasets are available on BioImages Archive (S-BIAD2914): https://www.ebi.ac.uk/biostudies/bioimages/studies/S-BIAD2914. Other data available upon reasonable request from the corresponding author.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 3 authors, 9 keywords, 10 MeSH terms, 1 funder, 55 references.

Cite

This paper

Rodriguez-Pastrana, I., Richens, J., & Wilkinson, R. N. (2026). Automated analysis of zebrafish vascular networks using the VISTA-Z pipeline. Scientific reports, 16(1), 15611. https://doi.org/10.1038/s41598-026-43301-5

BibTeX

@article{rodriguezpastrana2026automated, author = {Rodriguez-Pastrana, Ignacio and Richens, Joanna and Wilkinson, Robert N}, title = {{Automated analysis of zebrafish vascular networks using the VISTA-Z pipeline}}, journal = {Scientific reports}, year = {2026}, month = apr, volume = {16}, number = {1}, pages = {15611}, publisher = {Nature Publishing Group}, issn = {2045-2322}, doi = {10.1038/s41598-026-43301-5}, url = {https://doi.org/10.1038/s41598-026-43301-5}, pmid = {41922411}, pmcid = {PMC13187487} }

RIS

TY - JOUR AU - Rodriguez-Pastrana, Ignacio AU - Richens, Joanna AU - Wilkinson, Robert N TI - Automated analysis of zebrafish vascular networks using the VISTA-Z pipeline T2 - Scientific reports J2 - Sci Rep PY - 2026 DA - 2026/04/01 VL - 16 IS - 1 SP - 15611 SN - 2045-2322 PB - Nature Publishing Group DO - 10.1038/s41598-026-43301-5 UR - https://doi.org/10.1038/s41598-026-43301-5 LA - en ER -

CSL-JSON

{ "id": "10.1038/s41598-026-43301-5", "type": "article-journal", "title": "Automated analysis of zebrafish vascular networks using the VISTA-Z pipeline", "container-title": "Scientific reports", "author": [ { "family": "Rodriguez-Pastrana", "given": "Ignacio" }, { "family": "Richens", "given": "Joanna" }, { "family": "Wilkinson", "given": "Robert N" } ], "container-title-short": "Sci Rep", "volume": "16", "issue": "1", "page": "15611", "DOI": "10.1038/s41598-026-43301-5", "PMID": "41922411", "PMCID": "PMC13187487", "ISSN": "2045-2322", "publisher": "Nature Publishing Group", "URL": "https://doi.org/10.1038/s41598-026-43301-5", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 1 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.