Merge pull request #3 from DigitalHolography/port_pipeline_waveform_shape_metrics

Rough implementation of old doppler view pipeline

Author: AdrienGordon

.gitignore

@@ -1,3 +1,9 @@
+# Temporary workspace files
+system_design/
+AGENTS.md
+uv.lock
+test/
+
 # Byte-compiled / cache
 __pycache__/
 *.py[cod]
@@ -21,6 +27,7 @@ dist/
 wheelhouse/
 
 # Test / coverage
+.tmp/
 .pytest_cache/
 .tox/
 .nox/
@@ -63,4 +70,4 @@ Desktop.ini
 process_result.csv
 pipelines.txt
 requirements-optional.txt
-eyeflow_matlab/
+eyeflow_matlab/

pyproject.toml

@@ -4,15 +4,24 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "EyeFlow"
-version = "0.1.0"
+version = "0.2.0"
 description = "Cohort-analysis engine for retinal Doppler holography"
 readme = "README.md"
 requires-python = ">=3.10"
 license = { text = "GPL-3.0-only" }
 
 # =============== [ DEPENDENCIES ] ===============
 
-dependencies = ["numpy>=1.24", "h5py>=3.9", "sv-ttk>=2.6", "tkinterdnd2"]
+dependencies = [
+    "numpy>=1.24",
+    "h5py>=3.9",
+    "pydantic>=2.6",
+    "scipy>=1.10",
+    "scikit-image>=0.20",
+    "joblib>=1.3",
+    "sv-ttk>=2.6",
+    "tkinterdnd2",
+]
 
 [project.optional-dependencies]
 # For specific pipelines

src/app_settings.py

@@ -12,7 +12,7 @@
 APP_NAME = "EyeFlow"
 SETTINGS_FILENAME = "settings.json"
 DEFAULT_SETTINGS_FILENAME = "default_settings.json"
-LAST_BATCH_LOG_FILENAME = "last_batch_log.txt"
+LAST_BATCH_LOG_FILENAME = "last_EF_log.txt"
 VERSION_PATTERN = re.compile(r'^version\s*=\s*"([^"]+)"\s*$')
 INVALID_PATH_CHARS_PATTERN = re.compile(r'[<>:"/\\|?*]+')
 

src/cli.py

@@ -23,6 +23,10 @@
 from collections.abc import Callable, Sequence
 from pathlib import Path
 
+from runtime_limits import configure_numeric_threads
+
+configure_numeric_threads()
+
 import h5py
 
 from pipelines import (
@@ -31,7 +35,7 @@
     load_pipeline_catalog,
 )
 from pipelines.core.errors import format_pipeline_exception
-from pipelines.core.utils import write_combined_results_h5
+from input_output import write_combined_results_h5
 
 
 def _build_pipeline_registry() -> dict[str, PipelineDescriptor]:

src/domain/__init__.py

@@ -0,0 +1,13 @@
+"""Pure EyeFlow domain calculations."""
+
+from .steps import (
+    ArterialWaveformAnalysisStep,
+    DomainStep,
+    VesselVelocityEstimatorStep,
+)
+
+__all__ = [
+    "ArterialWaveformAnalysisStep",
+    "DomainStep",
+    "VesselVelocityEstimatorStep",
+]

src/domain/blood_flow_velocity/__init__.py

@@ -0,0 +1,13 @@
+"""Pure blood-flow velocity calculations ported from Matlab EyeFlow."""
+
+from .per_beat import PerBeatAnalysisInput, PerBeatAnalysisResult, run_per_beat_analysis
+from .per_beat_signal import PerBeatSignalAnalysisResult, per_beat_signal_analysis
+
+__all__ = [
+    "PerBeatAnalysisInput",
+    "PerBeatAnalysisResult",
+    "PerBeatSignalAnalysisResult",
+    "per_beat_signal_analysis",
+    "run_per_beat_analysis",
+]
+

src/domain/blood_flow_velocity/_signal_utils.py

@@ -0,0 +1,58 @@
+"""Shared helpers for blood-flow velocity signal analysis."""
+
+from __future__ import annotations
+
+import numpy as np
+
+
+def next_power_of_two(value: int) -> int:
+    if value < 1:
+        raise ValueError("next_power_of_two expects a strictly positive integer.")
+    return 1 << (value - 1).bit_length()
+
+
+def normalize_cycle_boundaries(
+    cycle_boundaries,
+    signal_length: int,
+    *,
+    index_base: int | None = None,
+) -> np.ndarray:
+    boundaries = np.asarray(cycle_boundaries, dtype=np.int64).reshape(-1)
+    if boundaries.size < 2:
+        raise ValueError("At least two cycle boundaries are required.")
+    if signal_length <= 0:
+        raise ValueError("signal_length must be positive.")
+
+    inferred_index_base = _infer_index_base(boundaries, signal_length, index_base)
+    normalized = boundaries - int(inferred_index_base)
+    _validate_cycle_boundaries(normalized, signal_length, inferred_index_base)
+    return normalized.astype(np.int64, copy=False)
+
+
+def _infer_index_base(
+    boundaries: np.ndarray,
+    signal_length: int,
+    index_base: int | None,
+) -> int:
+    if index_base is not None:
+        return int(index_base)
+    if np.any(boundaries == 0):
+        return 0
+    if np.any(boundaries == signal_length):
+        return 1
+    return 1
+
+
+def _validate_cycle_boundaries(
+    normalized: np.ndarray,
+    signal_length: int,
+    index_base: int,
+) -> None:
+    if np.any(np.diff(normalized) <= 0):
+        raise ValueError("Cycle boundaries must be strictly increasing.")
+    if normalized[0] >= 0 and normalized[-1] < signal_length:
+        return
+    raise ValueError(
+        "Cycle boundaries fall outside the available signal length after "
+        f"normalization (index_base={index_base})."
+    )

src/domain/blood_flow_velocity/per_beat.py

@@ -0,0 +1,83 @@
+"""Port of BloodFlowVelocity/perBeatAnalysis.m."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from .per_beat_signal import PerBeatSignalAnalysisResult, per_beat_signal_analysis
+
+
+@dataclass(frozen=True)
+class VesselPerBeatAnalysisResult:
+    signal: PerBeatSignalAnalysisResult
+    vmax_band_limited: np.ndarray
+    vmin_band_limited: np.ndarray
+    vti_per_beat: np.ndarray
+
+
+@dataclass(frozen=True)
+class PerBeatAnalysisInput:
+    arterial_velocity_signal: np.ndarray
+    venous_velocity_signal: np.ndarray
+    systolic_acceleration_peak_indexes: np.ndarray
+    band_limited_signal_harmonic_count: int
+    dt_seconds: float
+    beat_period_seconds: np.ndarray | None = None
+    index_base: int | None = None
+
+
+@dataclass(frozen=True)
+class PerBeatAnalysisResult:
+    beat_period_idx: np.ndarray
+    beat_period_seconds: np.ndarray
+    artery: VesselPerBeatAnalysisResult
+    vein: VesselPerBeatAnalysisResult
+
+
+def run_per_beat_analysis(inputs: PerBeatAnalysisInput) -> PerBeatAnalysisResult:
+    sys_idx_list = np.asarray(
+        inputs.systolic_acceleration_peak_indexes,
+        dtype=np.int64,
+    ).reshape(-1)
+    beat_period_idx = np.diff(sys_idx_list).astype(np.int32, copy=False)
+    return PerBeatAnalysisResult(
+        beat_period_idx=beat_period_idx,
+        beat_period_seconds=_beat_period_seconds(inputs, beat_period_idx),
+        vein=_run_vessel(inputs.venous_velocity_signal, sys_idx_list, inputs),
+        artery=_run_vessel(inputs.arterial_velocity_signal, sys_idx_list, inputs),
+    )
+
+
+def _beat_period_seconds(
+    inputs: PerBeatAnalysisInput,
+    beat_period_idx: np.ndarray,
+) -> np.ndarray:
+    if inputs.beat_period_seconds is not None:
+        periods = np.asarray(inputs.beat_period_seconds, dtype=np.float64).reshape(-1)
+        if periods.size == beat_period_idx.size:
+            return periods
+    return beat_period_idx.astype(np.float64, copy=False) * float(inputs.dt_seconds)
+
+
+def _run_vessel(
+    velocity_signal,
+    sys_idx_list: np.ndarray,
+    inputs: PerBeatAnalysisInput,
+) -> VesselPerBeatAnalysisResult:
+    signal = per_beat_signal_analysis(
+        velocity_signal,
+        sys_idx_list,
+        inputs.band_limited_signal_harmonic_count,
+        index_base=inputs.index_base,
+    )
+    return VesselPerBeatAnalysisResult(
+        signal=signal,
+        vmax_band_limited=np.max(signal.velocity_signal_per_beat_band_limited, axis=1),
+        vmin_band_limited=np.min(signal.velocity_signal_per_beat_band_limited, axis=1),
+        vti_per_beat=(
+            np.sum(signal.velocity_signal_per_beat, axis=1) * float(inputs.dt_seconds)
+        ),
+    )
+

src/domain/blood_flow_velocity/per_beat_signal.py

@@ -0,0 +1,122 @@
+"""Port of BloodFlowVelocity/perBeatSignalAnalysis.m."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from ._signal_utils import next_power_of_two, normalize_cycle_boundaries
+
+
+@dataclass(frozen=True)
+class PerBeatSignalAnalysisResult:
+    velocity_signal_per_beat: np.ndarray
+    velocity_signal_per_beat_fft: np.ndarray
+    velocity_signal_per_beat_band_limited: np.ndarray
+
+
+def _interpft_real(signal: np.ndarray, target_length: int) -> np.ndarray:
+    source = np.asarray(signal, dtype=np.float64).reshape(-1)
+    source_length = int(source.size)
+    if source_length == 0:
+        raise ValueError("interpft requires a non-empty signal.")
+    if target_length <= 0:
+        raise ValueError("interpft target_length must be positive.")
+    if target_length == source_length:
+        return source.copy()
+
+    spectrum = np.fft.fft(source)
+    resized = np.zeros(int(target_length), dtype=np.complex128)
+    _copy_resized_spectrum(spectrum, resized, source_length)
+    interpolated = np.fft.ifft(resized) * (float(target_length) / float(source_length))
+    return interpolated.real
+
+
+def _copy_resized_spectrum(
+    spectrum: np.ndarray,
+    resized: np.ndarray,
+    source_length: int,
+) -> None:
+    if source_length % 2 == 0:
+        half = source_length // 2
+        resized[:half] = spectrum[:half]
+        resized[-(source_length - half - 1) :] = spectrum[half + 1 :]
+        resized[half] = spectrum[half] / 2.0
+        resized[resized.size - half] = spectrum[half] / 2.0
+        return
+
+    pivot = source_length // 2 + 1
+    resized[:pivot] = spectrum[:pivot]
+    resized[-(source_length // 2) :] = spectrum[pivot:]
+
+
+def per_beat_signal_analysis(
+    signal,
+    sys_idx_list,
+    band_limited_signal_harmonic_count: int,
+    *,
+    index_base: int | None = None,
+) -> PerBeatSignalAnalysisResult:
+    signal_array = np.asarray(signal, dtype=np.float64).reshape(-1)
+    if signal_array.size == 0:
+        raise ValueError("signal must contain at least one sample.")
+    if band_limited_signal_harmonic_count < 1:
+        raise ValueError("band_limited_signal_harmonic_count must be positive.")
+
+    cycle_boundaries = normalize_cycle_boundaries(
+        sys_idx_list,
+        signal_array.size,
+        index_base=index_base,
+    )
+    return _analyze_cycles(
+        signal_array,
+        cycle_boundaries,
+        int(band_limited_signal_harmonic_count),
+    )
+
+
+def _analyze_cycles(
+    signal_array: np.ndarray,
+    cycle_boundaries: np.ndarray,
+    harmonic_count: int,
+) -> PerBeatSignalAnalysisResult:
+    number_of_beats = int(cycle_boundaries.size - 1)
+    n_fft = next_power_of_two(int(np.max(np.diff(cycle_boundaries))))
+    per_beat, per_beat_fft, band_limited = _empty_outputs(number_of_beats, n_fft)
+
+    for beat_index in range(number_of_beats):
+        start = int(cycle_boundaries[beat_index])
+        stop = int(cycle_boundaries[beat_index + 1]) + 1
+        beat_interp = _interpft_real(signal_array[start:stop], n_fft + 1)[:-1]
+        beat_fft = np.fft.fft(beat_interp, n=n_fft)
+        per_beat[beat_index, :] = beat_interp
+        per_beat_fft[beat_index, :] = beat_fft
+        band_limited[beat_index, :] = _band_limited_signal(
+            beat_fft,
+            n_fft,
+            harmonic_count,
+        )
+
+    return PerBeatSignalAnalysisResult(per_beat, per_beat_fft, band_limited)
+
+
+def _empty_outputs(
+    number_of_beats: int,
+    n_fft: int,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    per_beat = np.full((number_of_beats, n_fft), np.nan, dtype=np.float64)
+    per_beat_fft = np.full((number_of_beats, n_fft), np.nan + 0j, dtype=np.complex128)
+    band_limited = np.full((number_of_beats, n_fft), np.nan, dtype=np.float64)
+    return per_beat, per_beat_fft, band_limited
+
+
+def _band_limited_signal(
+    beat_fft: np.ndarray,
+    n_fft: int,
+    harmonic_count: int,
+) -> np.ndarray:
+    count = min(int(harmonic_count), n_fft)
+    band_limited_spectrum = beat_fft[:count].copy() * 2.0
+    band_limited_spectrum[0] = beat_fft[0]
+    return np.abs(np.fft.ifft(band_limited_spectrum, n=n_fft))

src/domain/steps/__init__.py

@@ -0,0 +1,15 @@
+"""Pure domain steps migrated from DopplerView calculation steps."""
+
+from .arterial_waveform_analysis import (
+    ArterialWaveformAnalysisStep,
+)
+from .base import DomainStep
+from .vessel_velocity_estimator import (
+    VesselVelocityEstimatorStep,
+)
+
+__all__ = [
+    "ArterialWaveformAnalysisStep",
+    "DomainStep",
+    "VesselVelocityEstimatorStep",
+]