# -*- coding: utf-8 -*-
# Copyright (c) 2025 University Medical Center Göttingen, Germany.
# All rights reserved.
#
# Patent Pending: DE 10 2024 112 939.5
# SPDX-License-Identifier: LicenseRef-Proprietary-See-LICENSE
#
# This software is licensed under a custom license. See the LICENSE file
# in the root directory for full details.
#
# **Commercial use is prohibited without a separate license.**
# Contact MBM ScienceBridge GmbH (https://sciencebridge.de/en/) for licensing.
"""Domain motion analysis module.
This module provides functions for analyzing sarcomere contraction dynamics
within sarcomere domains over time, computing per-domain summary statistics
and contraction parameters.
"""
import logging
import os
from typing import Dict, List, Optional, Tuple, Union
import numpy as np
from scipy.ndimage import binary_closing, binary_opening, label
from scipy.signal import savgol_filter
from skimage.segmentation import clear_border
from contraction_net.prediction import predict_contractions
from sarcasm.structure_modules.domain_clustering import (
analyze_domains,
assign_vectors_to_domains,
)
logger = logging.getLogger(__name__)
[docs]
def compute_domain_timeseries(
pos_vectors_all: List[np.ndarray],
sarcomere_length_vectors_all: List[np.ndarray],
domain_mask: np.ndarray,
pixelsize: float,
n_domains: int,
) -> Dict[str, np.ndarray]:
"""
Compute per-domain sarcomere length statistics over time.
For each frame, assigns sarcomere vectors to domains based on their position
and computes summary statistics (mean, std, quartiles) of sarcomere lengths
within each domain.
Parameters
----------
pos_vectors_all : List[np.ndarray]
List of position vectors for each frame. Each element is shape (n_vectors, 2) in µm.
sarcomere_length_vectors_all : List[np.ndarray]
List of sarcomere length vectors for each frame. Each element is shape (n_vectors,) in µm.
domain_mask : np.ndarray
Integer-labeled domain mask from the reference frame. Domain IDs are 1, 2, 3, etc.
Background pixels have value 0.
pixelsize : float
Pixel size in µm.
n_domains : int
Number of domains in the mask (excluding background).
Returns
-------
dict
Dictionary containing per-domain time-series:
- 'domain_slen_timeseries': np.ndarray, shape (n_domains, n_frames), mean sarcomere length
- 'domain_slen_median_timeseries': np.ndarray, shape (n_domains, n_frames), median sarcomere length
- 'domain_slen_std_timeseries': np.ndarray, shape (n_domains, n_frames), std of sarcomere length
- 'domain_slen_q25_timeseries': np.ndarray, shape (n_domains, n_frames), 25th percentile
- 'domain_slen_q75_timeseries': np.ndarray, shape (n_domains, n_frames), 75th percentile
- 'domain_n_vectors_timeseries': np.ndarray, shape (n_domains, n_frames), number of vectors
"""
n_frames = len(pos_vectors_all)
# Initialize output arrays
domain_slen_mean = np.full((n_domains, n_frames), np.nan)
domain_slen_median = np.full((n_domains, n_frames), np.nan)
domain_slen_std = np.full((n_domains, n_frames), np.nan)
domain_slen_q25 = np.full((n_domains, n_frames), np.nan)
domain_slen_q75 = np.full((n_domains, n_frames), np.nan)
domain_n_vectors = np.zeros((n_domains, n_frames), dtype=np.int32)
# Process each frame
for frame_idx, (pos_vectors, sarcomere_lengths) in enumerate(
zip(pos_vectors_all, sarcomere_length_vectors_all)
):
if pos_vectors is None or len(pos_vectors) == 0:
continue
# Assign vectors to domains
domain_ids = assign_vectors_to_domains(pos_vectors, domain_mask, pixelsize)
# Compute statistics for each domain
for domain_id in range(1, n_domains + 1):
mask = domain_ids == domain_id
n_vec = np.sum(mask)
domain_n_vectors[domain_id - 1, frame_idx] = n_vec
if n_vec > 0:
lengths = sarcomere_lengths[mask]
domain_slen_mean[domain_id - 1, frame_idx] = np.nanmean(lengths)
domain_slen_median[domain_id - 1, frame_idx] = np.nanmedian(lengths)
domain_slen_std[domain_id - 1, frame_idx] = np.nanstd(lengths)
domain_slen_q25[domain_id - 1, frame_idx] = np.nanpercentile(lengths, 25)
domain_slen_q75[domain_id - 1, frame_idx] = np.nanpercentile(lengths, 75)
return {
'domain_slen_timeseries': domain_slen_mean,
'domain_slen_median_timeseries': domain_slen_median,
'domain_slen_std_timeseries': domain_slen_std,
'domain_slen_q25_timeseries': domain_slen_q25,
'domain_slen_q75_timeseries': domain_slen_q75,
'domain_n_vectors_timeseries': domain_n_vectors,
}
[docs]
def detect_domain_contractions(
domain_slen_timeseries: np.ndarray,
frametime: float,
model_path: str,
threshold: float = 0.3,
contr_time_min: float = 0.2,
merge_time_max: float = 0.05,
buffer_frames: int = 3,
min_valid_frames: float = 0.5,
) -> Dict[str, np.ndarray]:
"""
Detect contraction cycles from domain-averaged sarcomere length time-series using ContractionNet.
Uses the ContractionNet neural network to predict contraction states from the
mean sarcomere length signal of each domain, then applies morphological operations
to clean up the predictions.
Parameters
----------
domain_slen_timeseries : np.ndarray
Per-domain mean sarcomere length time-series. Shape (n_domains, n_frames).
frametime : float
Time between frames in seconds.
model_path : str
Path to the ContractionNet model weights (.pt file).
threshold : float, optional
Binary threshold for contraction state prediction. Default 0.3.
contr_time_min : float, optional
Minimal time of contraction in seconds. Shorter contractions are removed. Default 0.2.
merge_time_max : float, optional
Maximal time between two contractions. Closer contractions are merged. Default 0.05.
buffer_frames : int, optional
Remove contraction cycles within this many frames of start/end of time-series. Default 3.
min_valid_frames : float, optional
Minimum fraction of valid (non-NaN) frames required for a domain to be analyzed. Default 0.5.
Returns
-------
dict
Dictionary containing per-domain contraction detection results:
- 'domain_contr': np.ndarray, shape (n_domains, n_frames), binary contraction state
- 'domain_n_contr': np.ndarray, shape (n_domains,), number of contractions per domain
- 'domain_labels_contr': np.ndarray, shape (n_domains, n_frames), contraction cycle labels
- 'domain_beating_rate': np.ndarray, shape (n_domains,), beating rate in Hz
- 'domain_beating_rate_variability': np.ndarray, shape (n_domains,), std of inter-beat interval
"""
n_domains, n_frames = domain_slen_timeseries.shape
# Initialize output arrays
domain_contr = np.zeros((n_domains, n_frames), dtype=bool)
domain_n_contr = np.zeros(n_domains, dtype=np.int32)
domain_labels_contr = np.zeros((n_domains, n_frames), dtype=np.int32)
domain_beating_rate = np.full(n_domains, np.nan)
domain_beating_rate_var = np.full(n_domains, np.nan)
# Morphological structuring elements
structure_closing = np.ones(max(1, int(merge_time_max / frametime)))
structure_opening = np.ones(max(1, int(contr_time_min / frametime)))
# Process each domain
for domain_idx in range(n_domains):
slen_timeseries = domain_slen_timeseries[domain_idx]
# Check if domain has enough valid data
valid_fraction = np.sum(~np.isnan(slen_timeseries)) / n_frames
if valid_fraction < min_valid_frames:
logger.debug(f"Domain {domain_idx + 1} has insufficient valid data ({valid_fraction:.1%}), skipping.")
continue
# Interpolate NaN values for prediction
slen_interp = _interpolate_nans(slen_timeseries)
if np.all(np.isnan(slen_interp)):
continue
# Predict contractions using ContractionNet
try:
contr_pred = predict_contractions(slen_interp, model_path)
contr = contr_pred[0] > threshold
except Exception as e:
logger.warning(f"ContractionNet prediction failed for domain {domain_idx + 1}: {e}")
continue
# Apply morphological operations to clean up predictions
contr = binary_opening(binary_closing(contr, structure=structure_closing), structure=structure_opening)
# Remove incomplete contractions at beginning/end
contr = clear_border(contr, buffer_size=buffer_frames)
# Store binary contraction state
domain_contr[domain_idx] = contr
# Label contraction cycles
labels, n_contr = label(contr)
domain_labels_contr[domain_idx] = labels
domain_n_contr[domain_idx] = n_contr
# Calculate beating rate
if n_contr > 1:
start_frames = np.where(np.diff(contr.astype('float32')) > 0.5)[0]
if len(start_frames) > 1:
inter_beat_intervals = np.diff(start_frames) * frametime
domain_beating_rate[domain_idx] = 1 / np.mean(inter_beat_intervals)
domain_beating_rate_var[domain_idx] = np.std(inter_beat_intervals)
else:
logger.warning(f"Domain {domain_idx + 1}: Only {n_contr} contraction cycle(s) detected. "
f"Cannot compute beating rate (requires >= 2 cycles).")
return {
'domain_contr': domain_contr,
'domain_n_contr': domain_n_contr,
'domain_labels_contr': domain_labels_contr,
'domain_beating_rate': domain_beating_rate,
'domain_beating_rate_variability': domain_beating_rate_var,
}
[docs]
def analyze_domain_contraction_parameters(
domain_slen_timeseries: np.ndarray,
domain_labels_contr: np.ndarray,
domain_n_contr: np.ndarray,
frametime: float,
filter_params: Tuple[int, int] = (13, 5),
) -> Dict[str, np.ndarray]:
"""
Analyze contraction parameters for domain-averaged sarcomere length trajectories.
Computes per-domain, per-contraction-cycle parameters including maximum contraction,
maximum elongation, velocities, and timing parameters.
Parameters
----------
domain_slen_timeseries : np.ndarray
Per-domain mean sarcomere length time-series. Shape (n_domains, n_frames).
domain_labels_contr : np.ndarray
Per-domain contraction cycle labels. Shape (n_domains, n_frames).
domain_n_contr : np.ndarray
Number of contractions per domain. Shape (n_domains,).
frametime : float
Time between frames in seconds.
filter_params : Tuple[int, int], optional
Savitzky-Golay filter parameters (window_length, polyorder) for velocity smoothing.
Default (13, 5).
Returns
-------
dict
Dictionary containing per-domain contraction parameters:
- 'domain_equ': np.ndarray, shape (n_domains,), equilibrium/resting sarcomere length
- 'domain_contr_max': np.ndarray, shape (n_domains, max_n_contr), max contraction per cycle
- 'domain_elong_max': np.ndarray, shape (n_domains, max_n_contr), max elongation per cycle
- 'domain_vel_contr_max': np.ndarray, shape (n_domains, max_n_contr), max shortening velocity
- 'domain_vel_elong_max': np.ndarray, shape (n_domains, max_n_contr), max elongation velocity
- 'domain_time_to_peak': np.ndarray, shape (n_domains, max_n_contr), time to maximal contraction
- 'domain_time_to_relax': np.ndarray, shape (n_domains, max_n_contr), time from peak to relaxation
- 'domain_time_contr': np.ndarray, shape (n_domains, max_n_contr), contraction duration
"""
n_domains = domain_slen_timeseries.shape[0]
max_n_contr = int(np.max(domain_n_contr)) if np.max(domain_n_contr) > 0 else 1
# Initialize output arrays
domain_equ = np.full(n_domains, np.nan)
domain_contr_max = np.full((n_domains, max_n_contr), np.nan)
domain_elong_max = np.full((n_domains, max_n_contr), np.nan)
domain_vel_contr_max = np.full((n_domains, max_n_contr), np.nan)
domain_vel_elong_max = np.full((n_domains, max_n_contr), np.nan)
domain_time_to_peak = np.full((n_domains, max_n_contr), np.nan)
domain_time_to_relax = np.full((n_domains, max_n_contr), np.nan)
domain_time_contr = np.full((n_domains, max_n_contr), np.nan)
window_length, polyorder = filter_params
for domain_idx in range(n_domains):
slen = domain_slen_timeseries[domain_idx]
labels = domain_labels_contr[domain_idx]
n_contr = domain_n_contr[domain_idx]
if n_contr == 0 or np.all(np.isnan(slen)):
continue
# Calculate equilibrium length (median of non-NaN values)
valid_slen = slen[~np.isnan(slen)]
if len(valid_slen) > 0:
domain_equ[domain_idx] = np.median(valid_slen)
# Calculate velocity using Savitzky-Golay filter
slen_interp = _interpolate_nans(slen)
if len(slen_interp) >= window_length:
vel = savgol_filter(slen_interp, window_length, polyorder, deriv=1, delta=frametime)
else:
vel = np.gradient(slen_interp, frametime)
# Calculate delta (change from equilibrium)
delta_slen = slen_interp - domain_equ[domain_idx]
# Analyze each contraction cycle
for contr_idx in range(n_contr):
cycle_mask = labels == (contr_idx + 1)
if not np.any(cycle_mask):
continue
delta_cycle = delta_slen[cycle_mask]
vel_cycle = vel[cycle_mask]
# Contraction duration
domain_time_contr[domain_idx, contr_idx] = np.sum(cycle_mask) * frametime
# Max contraction (most negative delta) and elongation (most positive delta)
if len(delta_cycle) > 0:
domain_contr_max[domain_idx, contr_idx] = np.nanmin(delta_cycle)
domain_elong_max[domain_idx, contr_idx] = np.nanmax(delta_cycle)
# Max velocities
domain_vel_contr_max[domain_idx, contr_idx] = np.nanmin(vel_cycle)
domain_vel_elong_max[domain_idx, contr_idx] = np.nanmax(vel_cycle)
# Time to peak (time from start to minimum)
if not np.all(np.isnan(delta_cycle)):
peak_idx = np.nanargmin(delta_cycle)
domain_time_to_peak[domain_idx, contr_idx] = peak_idx * frametime
domain_time_to_relax[domain_idx, contr_idx] = (len(delta_cycle) - peak_idx) * frametime
return {
'domain_equ': domain_equ,
'domain_contr_max': domain_contr_max,
'domain_elong_max': domain_elong_max,
'domain_vel_contr_max': domain_vel_contr_max,
'domain_vel_elong_max': domain_vel_elong_max,
'domain_time_to_peak': domain_time_to_peak,
'domain_time_to_relax': domain_time_to_relax,
'domain_time_contr': domain_time_contr,
}
[docs]
def _interpolate_nans(arr: np.ndarray) -> np.ndarray:
"""
Interpolate NaN values in a 1D array using linear interpolation.
Parameters
----------
arr : np.ndarray
1D array potentially containing NaN values.
Returns
-------
np.ndarray
Array with NaN values interpolated.
"""
arr = arr.copy()
nans = np.isnan(arr)
if np.all(nans):
return arr
if np.any(nans):
indices = np.arange(len(arr))
arr[nans] = np.interp(indices[nans], indices[~nans], arr[~nans])
return arr