diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index bf89f33..6c1d294 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -10,19 +10,19 @@ jobs: fail-fast: false matrix: version: - - '1.6' + - '1' - 'nightly' os: - ubuntu-latest arch: - x64 steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: actions/cache@v1 + - uses: actions/cache@v4 env: cache-name: cache-artifacts with: @@ -35,15 +35,15 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v4 with: file: lcov.info docs: name: Documentation runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 with: version: '1' - run: | diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..4797afb --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,72 @@ +# CLAUDE.md + +This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository. + +## What This Package Does + +SMLMFrameConnection performs frame-connection on 2D SMLM localization data: combining repeated localizations of a single blinking fluorophore across frames into a single higher-precision localization. Implements the spatiotemporal LAP algorithm from Schodt & Lidke (2021). + +## Development Commands + +```bash +# Run tests +julia --project=. -e 'using Pkg; Pkg.test()' + +# Run a single test file +julia --project=. -e 'using Test, SMLMFrameConnection, SMLMData; include("test/test_helpers.jl"); include("test/test_frameconnect.jl")' + +# Build docs locally +julia --project=docs -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' +julia --project=docs docs/make.jl + +# Quick REPL usage +julia --project=. +``` + +Test files must `include("test/test_helpers.jl")` first since `runtests.jl` loads shared fixtures from it. + +## Architecture + +### Algorithm Pipeline + +`frameconnect()` in `frameconnect.jl` is the main entry point. It orchestrates a 4-stage pipeline: + +1. **Precluster** (`precluster.jl`): Spatiotemporal clustering using KDTree nearest-neighbor search. Groups localizations within `max_sigma_dist * mean(σ)` distance and `max_frame_gap` frames. Stores cluster assignments in emitter `track_id` fields. + +2. **Estimate parameters** (`estimateparams.jl`, `estimatedensities.jl`): Estimates photophysics rates (k_on, k_off, k_bleach, p_miss) from precluster statistics. Uses Optim.jl NelderMead to fit cumulative localization counts. `estimatedensities.jl` estimates local emitter density per cluster using KDTree neighbor distances. + +3. **Connect via LAP** (`connectlocalizations.jl` -> `create_costmatrix.jl` -> `solveLAP.jl` -> `linkclusters.jl`): For each multi-emitter precluster, builds a 2Nx2N cost matrix with connection/birth/death blocks using negative log-likelihoods from spatial separation, observation probability, and photophysics. Solves with Hungarian.jl. `linkclusters.jl` updates `connectID` from LAP assignments. + +4. **Combine** (`combinelocalizations.jl`): MLE weighted mean using full 2x2 covariance (precision-weighted). Produces higher-precision output localizations. + +### Data Flow + +- Input/output: `BasicSMLD` from SMLMData.jl containing `Emitter2DFit` emitters +- Internal representation: `organizeclusters()` converts to `Vector{Matrix{Float32}}` where each matrix is a cluster with columns: `[x, y, σ_x, σ_y, σ_xy, frame, dataset, connectID, sortindex]` +- `ParamStruct`: mutable struct accumulating estimated parameters through the pipeline +- Return: tuple `(combined::BasicSMLD, info::FrameConnectInfo)` + +### Key Types + +- `FrameConnectConfig <: AbstractSMLMConfig`: User-facing config (keyword-constructible) +- `FrameConnectInfo{T} <: AbstractSMLMInfo`: Algorithm output metadata (track assignments, rates, timing) +- `ParamStruct`: Internal mutable state for pipeline parameters (not exported for direct use) + +### Dual API Pattern + +`frameconnect()` accepts both kwargs and a `FrameConnectConfig` struct. The kwargs form constructs the config internally. Both `FrameConnectConfig` and `FrameConnectInfo` inherit from SMLMData abstract types. + +## Dependencies + +- **SMLMData.jl** (v0.7): Provides `BasicSMLD`, `Emitter2DFit`, `IdealCamera`, abstract config/info types +- **Hungarian.jl**: LAP solver in `solveLAP.jl` +- **NearestNeighbors.jl**: KDTree for preclustering and density estimation +- **Optim.jl**: Parameter estimation (Fminbox + NelderMead) + +## Conventions + +- Positions and uncertainties in microns +- Frame numbers are 1-based integers +- `track_id=0` means unconnected; populated values are compressed to `1:n_tracks` +- Mutating functions use `!` suffix with non-mutating wrappers (e.g., `connectlocalizations!`/`connectlocalizations`) +- Emitter precision type `ET` is derived from emitter field values, not SMLD type parameter diff --git a/Project.toml b/Project.toml index 540df8e..9e9e07f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SMLMFrameConnection" uuid = "1517b1a9-81e8-461f-b994-92eb29599690" authors = ["klidke@unm.edu"] -version = "0.3.0" +version = "0.5.0" [deps] Hungarian = "e91730f6-4275-51fb-a7a0-7064cfbd3b39" @@ -15,7 +15,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Hungarian = "0.6, 0.7" NearestNeighbors = "0.4" Optim = "1, 2" -SMLMData = "0.6" +SMLMData = "0.7" StatsBase = "0.33, 0.34" julia = "1.6" diff --git a/README.md b/README.md index 2885d8d..91933da 100644 --- a/README.md +++ b/README.md @@ -5,11 +5,7 @@ [![Build Status](https://github.com/JuliaSMLM/SMLMFrameConnection.jl/workflows/CI/badge.svg)](https://github.com/JuliaSMLM/SMLMFrameConnection.jl/actions) [![Coverage](https://codecov.io/gh/JuliaSMLM/SMLMFrameConnection.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaSMLM/SMLMFrameConnection.jl) -## Overview - -SMLMFrameConnection performs **frame-connection** on 2D localization microscopy data: combining repeated localizations of a single blinking fluorophore into a single higher-precision localization. - -Uses the spatiotemporal clustering algorithm from [Schodt & Lidke 2021](https://doi.org/10.3389/fbinf.2021.724325). +Frame-connection for single-molecule localization microscopy: linking localizations from the same fluorophore blinking event across consecutive frames into single, higher-precision localizations. Uses spatiotemporal LAP assignment to optimally connect temporally adjacent detections based on spatial proximity and estimated blinking kinetics. ## Installation @@ -21,124 +17,128 @@ Pkg.add("SMLMFrameConnection") ## Quick Start ```julia -using SMLMData, SMLMFrameConnection +using SMLMFrameConnection + +# Frame connection on localization data +(combined, info) = frameconnect(smld) + +# Output: combined high-precision localizations +println("$(info.n_input) → $(info.n_combined) localizations") +``` + +For complete SMLM workflows (detection + fitting + frame-connection + rendering), see [SMLMAnalysis.jl](https://github.com/JuliaSMLM/SMLMAnalysis.jl). + +## Configuration + +`frameconnect()` accepts keyword arguments or a `FrameConnectConfig` struct: + +| Parameter | Default | Description | +|-----------|---------|-------------| +| `n_density_neighbors` | 2 | Nearest preclusters for local density estimation | +| `max_sigma_dist` | 5.0 | Sigma multiplier for preclustering distance threshold | +| `max_frame_gap` | 5 | Maximum frame gap for temporal adjacency | +| `max_neighbors` | 2 | Maximum nearest-neighbors for precluster membership | -# Run frame connection on your data -smld_connected, smld_preclustered, smld_combined, params = frameconnect(smld) +```julia +# Keyword form (most common) +(combined, info) = frameconnect(smld; max_frame_gap=10, max_sigma_dist=3.0) -# smld_combined is the main output - higher precision localizations +# Config struct form (reusable settings) +config = FrameConnectConfig(max_frame_gap=10, max_sigma_dist=3.0) +(combined, info) = frameconnect(smld, config) ``` -## Input Requirements +**Parameter guidance:** Default values work well for standard dSTORM/PALM data. For dense samples, reduce `max_sigma_dist` to 3.0. For long dark states (dSTORM), increase `max_frame_gap` to 10-20. + +## Output Format + +`frameconnect()` returns `(combined::BasicSMLD, info::FrameConnectInfo)`. -Input must be a `BasicSMLD{T, Emitter2DFit{T}}` from [SMLMData.jl](https://github.com/JuliaSMLM/SMLMData.jl). +| Output | Description | +|--------|-------------| +| `combined` | High-precision combined localizations (main output) | +| `info.connected` | Input with `track_id` assigned (per-frame data with labels) | +| `info.n_input` | Number of input localizations | +| `info.n_tracks` | Number of tracks formed | +| `info.n_combined` | Number of output localizations | +| `info.k_on` | Estimated on rate (1/frame) | +| `info.k_off` | Estimated off rate (1/frame) | +| `info.k_bleach` | Estimated bleach rate (1/frame) | +| `info.p_miss` | Probability of missed detection | +| `info.initial_density` | Density estimate per cluster (emitters/μm²) | +| `info.elapsed_s` | Wall time (seconds) | +| `info.algorithm` | Algorithm used (`:lap`) | +| `info.n_preclusters` | Number of preclusters formed | -**Required fields** (algorithm fails without these): -- `x`, `y`: Position coordinates in microns -- `σ_x`, `σ_y`: Position uncertainties in microns (used for MLE weighting; must be > 0) -- `frame`: Frame number (1-based integer) +### Combination Method -**Optional fields** (combined in output if present): -- `photons`, `σ_photons`: Photon count and uncertainty (summed across connected localizations) -- `bg`, `σ_bg`: Background and uncertainty -- `dataset`: Dataset identifier (defaults to 1; for multi-ROI or multi-acquisition data) +Connected localizations are combined using maximum likelihood estimation (MLE) weighted mean: +- Position: `x_combined = Σ(x/σ²) / Σ(1/σ²)` (inverse-variance weighted) +- Uncertainty: `σ_combined = √(1/Σ(1/σ²))` ≈ `σ_individual / √n` +- Photons: summed across connected localizations -## Complete Example +## Algorithm Pipeline + +1. **Precluster**: Spatiotemporal clustering using KDTree nearest-neighbor search within `max_sigma_dist * σ` distance and `max_frame_gap` frames +2. **Estimate parameters**: Fit photophysics rates (k_on, k_off, k_bleach, p_miss) from precluster statistics using Optim.jl NelderMead +3. **Connect via LAP**: Build cost matrix from spatial separation and photophysics likelihoods, solve with Hungarian.jl +4. **Combine**: MLE weighted mean using full 2×2 covariance (precision-weighted) + +## Example ```julia -using SMLMData, SMLMFrameConnection +using SMLMFrameConnection # Create camera (pixel_ranges, pixel_size in μm) cam = IdealCamera(1:512, 1:512, 0.1) # Create emitters representing the same molecule blinking across 3 frames -# Constructor: Emitter2DFit{T}(x, y, photons, bg, σ_x, σ_y, σ_photons, σ_bg; frame, dataset, track_id) emitters = [ Emitter2DFit{Float64}( - 5.0, 5.0, # x, y position (μm) + 5.00, 5.00, # x, y position (μm) 1000.0, 10.0, # photons, background - 0.02, 0.02, # σ_x, σ_y uncertainties (μm) - 50.0, 2.0; # σ_photons, σ_bg - frame=1 + 0.02, 0.02, 0.0, # σ_x, σ_y, σ_xy + 50.0, 2.0, # σ_photons, σ_bg + 1, 1, 0, 1 # frame, dataset, track_id, id ), - Emitter2DFit{Float64}(5.01, 5.01, 1200.0, 12.0, 0.02, 0.02, 60.0, 2.0; frame=2), - Emitter2DFit{Float64}(5.02, 4.99, 1100.0, 11.0, 0.02, 0.02, 55.0, 2.0; frame=3), + Emitter2DFit{Float64}(5.01, 5.01, 1200.0, 12.0, 0.02, 0.02, 0.0, 60.0, 2.0, 2, 1, 0, 2), + Emitter2DFit{Float64}(5.02, 4.99, 1100.0, 11.0, 0.02, 0.02, 0.0, 55.0, 2.0, 3, 1, 0, 3), ] -# Create SMLD: BasicSMLD(emitters, camera, n_frames, n_datasets) smld = BasicSMLD(emitters, cam, 3, 1) # Run frame connection -_, _, smld_combined, _ = frameconnect(smld) - -# Result: localizations connected based on spatial/temporal proximity -# Combined uncertainty: σ_combined ≈ σ_individual / √n_connected -``` +(combined, info) = frameconnect(smld) -## Outputs Explained - -```julia -smld_connected, smld_preclustered, smld_combined, params = frameconnect(smld) +# Result: Combined uncertainty ≈ σ_individual / √n_connected +println("$(info.n_input) → $(info.n_combined) localizations in $(info.elapsed_s)s") ``` -| Output | Description | When to use | -|--------|-------------|-------------| -| `smld_combined` | **Main output.** Combined high-precision localizations | Standard analysis | -| `smld_connected` | Original localizations with `track_id` populated | When you need per-frame data with connection labels | -| `smld_preclustered` | Intermediate preclustering result | Debugging, algorithm tuning | -| `params` | Estimated photophysics + input parameters | Inspecting estimated k_on, k_off, density | - -## Parameters - -```julia -frameconnect(smld; - nnearestclusters = 2, # Clusters used for local density estimation - nsigmadev = 5.0, # Distance threshold = nsigmadev × localization uncertainty - maxframegap = 5, # Max frames between connected localizations - nmaxnn = 2 # Nearest neighbors checked during preclustering -) -``` - -**Parameter guidance:** -- `nsigmadev`: Higher values allow connections over larger distances. Default (5.0) works for typical SMLM data. Reduce for dense samples. -- `maxframegap`: Set based on expected blinking duration. For dSTORM with long dark states, increase to 10-20. -- Defaults work well for standard dSTORM/PALM data with typical blinking kinetics. - -## Estimated Parameters (ParamStruct) - -The algorithm estimates fluorophore photophysics from your data: - -| Field | Description | -|-------|-------------| -| `k_on` | Rate of transitioning from dark to visible state (1/frame) | -| `k_off` | Rate of transitioning from visible to dark state (1/frame) | -| `k_bleach` | Photobleaching rate (1/frame) | -| `p_miss` | Probability of missing a localization when fluorophore is on | -| `initialdensity` | Estimated emitter density per cluster (emitters/μm²) | - -## Combination Method - -Connected localizations are combined using **maximum likelihood estimation (MLE) weighted mean**: -- Position: inverse-variance weighted average → `x_combined = Σ(x/σ²) / Σ(1/σ²)` -- Uncertainty: `σ_combined = √(1/Σ(1/σ²))` ≈ `σ_individual / √n` -- Photons: summed across connected localizations - ## Utility Functions ### combinelocalizations ```julia smld_combined = combinelocalizations(smld) ``` -Combines emitters that share the same `track_id`. Use when you have pre-labeled data. +Combines emitters with the same `track_id`. Use when you have pre-labeled data. ### defineidealFC ```julia -smld_connected, smld_combined = defineidealFC(smld; maxframegap=5) +smld_connected, smld_combined = defineidealFC(smld; max_frame_gap=5) ``` -For **simulated data** where `track_id` already contains ground-truth emitter IDs. Useful for validating frame-connection performance against known truth. +For simulated data where `track_id` contains ground-truth emitter IDs. Validates frame-connection performance against known truth. + +## Algorithm Reference + +> Schodt, D.J. and Lidke, K.A. "Spatiotemporal Clustering of Repeated Super-Resolution Localizations via Linear Assignment Problem." *Frontiers in Bioinformatics*, 2021. [DOI: 10.3389/fbinf.2021.724325](https://doi.org/10.3389/fbinf.2021.724325) + +## Related Packages -## Citation +- **[SMLMAnalysis.jl](https://github.com/JuliaSMLM/SMLMAnalysis.jl)** - Complete SMLM workflow (detection + fitting + frame-connection + rendering) +- **[SMLMData.jl](https://github.com/JuliaSMLM/SMLMData.jl)** - Core data types for SMLM +- **[GaussMLE.jl](https://github.com/JuliaSMLM/GaussMLE.jl)** - GPU-accelerated Gaussian PSF fitting +- **[SMLMSim.jl](https://github.com/JuliaSMLM/SMLMSim.jl)** - SMLM data simulation -David J. Schodt and Keith A. Lidke, "Spatiotemporal Clustering of Repeated Super-Resolution Localizations via Linear Assignment Problem", Frontiers in Bioinformatics, 2021 +## License -https://doi.org/10.3389/fbinf.2021.724325 +MIT License - see [LICENSE](LICENSE) file for details. diff --git a/api_overview.md b/api_overview.md index aa34ece..48ceb65 100644 --- a/api_overview.md +++ b/api_overview.md @@ -1,60 +1,94 @@ # SMLMFrameConnection API Overview -Frame-connection for 2D SMLM data: combines repeated localizations of blinking fluorophores into higher-precision localizations. +Frame-connection for 2D SMLM data: linking localizations from the same fluorophore blinking event across consecutive frames into single, higher-precision localizations. Uses spatiotemporal LAP assignment to optimally connect temporally adjacent detections based on spatial proximity and estimated blinking kinetics. ## Main Functions ### frameconnect ```julia -frameconnect(smld::BasicSMLD{T,Emitter2DFit{T}}; - nnearestclusters::Int=2, - nsigmadev::Float64=5.0, - maxframegap::Int=5, - nmaxnn::Int=2 -) where T -> (smld_connected, smld_preclustered, smld_combined, params) +# Kwargs form (most common) +(combined, info) = frameconnect(smld::BasicSMLD; kwargs...) + +# Config form (reusable settings) +(combined, info) = frameconnect(smld::BasicSMLD, config::FrameConnectConfig) ``` Main entry point. Connects repeated localizations and combines them. -**Parameters:** -- `nnearestclusters`: Nearest preclusters for local density estimation -- `nsigmadev`: Sigma multiplier defining preclustering distance threshold (higher = larger connection radius) -- `maxframegap`: Maximum frame gap for temporal adjacency in preclusters -- `nmaxnn`: Maximum nearest-neighbors inspected for precluster membership - -**Returns:** -- `smld_connected`: Input with `track_id` populated (localizations uncombined) -- `smld_preclustered`: Intermediate preclustering result (for debugging) -- `smld_combined`: **Main output** - combined high-precision localizations -- `params::ParamStruct`: Algorithm parameters (input + estimated photophysics) +**Returns tuple `(combined, info)`:** +- `combined::BasicSMLD`: **Main output** - combined high-precision localizations +- `info::FrameConnectInfo`: Track assignments and algorithm metadata ### combinelocalizations ```julia -combinelocalizations(smld::BasicSMLD{T,Emitter2DFit{T}}) where T -> BasicSMLD +combinelocalizations(smld::BasicSMLD) -> BasicSMLD ``` Combines emitters sharing the same `track_id` using MLE weighted mean. Use when `track_id` is already populated. ### defineidealFC ```julia -defineidealFC(smld::BasicSMLD{T,Emitter2DFit{T}}; - maxframegap::Int=5 -) where T -> (smld_connected, smld_combined) +(smld_connected, smld_combined) = defineidealFC(smld::BasicSMLD; max_frame_gap::Int=5) ``` For simulated data where `track_id` indicates ground-truth emitter ID. Useful for validation/benchmarking. ## Types +### FrameConnectConfig +```julia +@kwdef struct FrameConnectConfig <: AbstractSMLMConfig + n_density_neighbors::Int = 2 # Nearest preclusters for local density estimation + max_sigma_dist::Float64 = 5.0 # Sigma multiplier for preclustering distance threshold + max_frame_gap::Int = 5 # Maximum frame gap for temporal adjacency + max_neighbors::Int = 2 # Maximum nearest-neighbors for precluster membership +end +``` +Configuration parameters for frame connection. Use with `frameconnect(smld, config)` or pass fields as kwargs to `frameconnect(smld; kwargs...)`. + +**Example:** +```julia +# Create config with custom settings +config = FrameConnectConfig(max_frame_gap=10, max_sigma_dist=3.0) +(combined, info) = frameconnect(smld, config) + +# Equivalent kwargs form +(combined, info) = frameconnect(smld; max_frame_gap=10, max_sigma_dist=3.0) +``` + +### FrameConnectInfo{T} +```julia +struct FrameConnectInfo{T} <: AbstractSMLMInfo + connected::BasicSMLD{T} # Input with track_id assigned (uncombined) + n_input::Int # Number of input localizations + n_tracks::Int # Number of tracks formed + n_combined::Int # Number of output localizations + k_on::Float64 # Rate: dark → visible (1/frame) + k_off::Float64 # Rate: visible → dark (1/frame) + k_bleach::Float64 # Rate: photobleaching (1/frame) + p_miss::Float64 # Probability of missed detection + initial_density::Vector{Float64} # Emitter density per cluster (emitters/μm²) + elapsed_s::Float64 # Wall time in seconds + algorithm::Symbol # Algorithm used (:lap) + n_preclusters::Int # Number of preclusters formed +end +``` + +**Access track assignments:** +```julia +(combined, info) = frameconnect(smld) +track_ids = [e.track_id for e in info.connected.emitters] +``` + ### ParamStruct ```julia mutable struct ParamStruct - initialdensity::Vector{Float64} # Emitter density per cluster (emitters/μm²) - nnearestclusters::Int # Clusters for density estimation + initial_density::Vector{Float64} # Emitter density per cluster (emitters/μm²) + n_density_neighbors::Int # Clusters for density estimation k_on::Float64 # Rate: dark → visible (1/frame) k_off::Float64 # Rate: visible → dark (1/frame) k_bleach::Float64 # Rate: photobleaching (1/frame) p_miss::Float64 # Probability of missing localization when on - nsigmadev::Float64 # Preclustering threshold multiplier - maxframegap::Int # Max frame gap in preclusters - nmaxnn::Int # Max nearest-neighbors for preclustering + max_sigma_dist::Float64 # Preclustering threshold multiplier + max_frame_gap::Int # Max frame gap in preclusters + max_neighbors::Int # Max nearest-neighbors for preclustering end ``` @@ -66,7 +100,7 @@ end ## Input Requirements -Input `BasicSMLD` must contain `Emitter2DFit` emitters. +Input `BasicSMLD` must contain emitters with position uncertainties. **Required fields:** - `x`, `y`: Position (microns) @@ -81,7 +115,7 @@ Input `BasicSMLD` must contain `Emitter2DFit` emitters. ## Output -Output `smld_combined` contains `Emitter2DFit` emitters with: +Output `combined` contains emitters with: - Combined position via MLE weighted mean: `x = Σ(x/σ²) / Σ(1/σ²)` - Reduced uncertainties: `σ = √(1/Σ(1/σ²))` - Summed photons @@ -89,8 +123,8 @@ Output `smld_combined` contains `Emitter2DFit` emitters with: ## Dependencies -- SMLMData.jl (v0.5+): BasicSMLD, Emitter2DFit types -- Hungarian.jl (v0.6-0.7): Linear assignment problem solver +- SMLMData.jl (v0.7): BasicSMLD, Emitter types, AbstractSMLMConfig, AbstractSMLMInfo +- Hungarian.jl: Linear assignment problem solver - NearestNeighbors.jl: Spatial clustering - Optim.jl: Parameter estimation diff --git a/docs/src/index.md b/docs/src/index.md index b3165e2..368acf5 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,7 +4,7 @@ CurrentModule = SMLMFrameConnection # SMLMFrameConnection -Frame-connection for 2D single molecule localization microscopy (SMLM) data. Combines repeated localizations of blinking fluorophores into higher-precision localizations using the algorithm from [Schodt & Lidke 2021](https://doi.org/10.3389/fbinf.2021.724325). +Frame-connection for 2D single molecule localization microscopy (SMLM) data: linking localizations from the same fluorophore blinking event across consecutive frames into single, higher-precision localizations. Uses spatiotemporal LAP assignment to optimally connect temporally adjacent detections based on spatial proximity and estimated blinking kinetics. See [Schodt & Lidke 2021](https://doi.org/10.3389/fbinf.2021.724325). ## Installation @@ -16,12 +16,13 @@ Pkg.add("SMLMFrameConnection") ## Quick Start ```julia -using SMLMData, SMLMFrameConnection +using SMLMFrameConnection # Run frame connection on your BasicSMLD with Emitter2DFit emitters -smld_connected, smld_preclustered, smld_combined, params = frameconnect(smld) +(combined, info) = frameconnect(smld) -# smld_combined is the main output - higher precision localizations +# combined is the main output - higher precision localizations +# info contains track assignments and algorithm metadata ``` ## Input Requirements @@ -48,24 +49,24 @@ Input `BasicSMLD` must contain `Emitter2DFit` emitters with: | Output | Description | |--------|-------------| -| `smld_combined` | **Main output** - combined high-precision localizations | -| `smld_connected` | Original localizations with `track_id` labels | -| `smld_preclustered` | Intermediate result (debugging) | -| `params` | Estimated photophysics parameters | +| `combined` | **Main output** - combined high-precision localizations | +| `info.connected` | Original localizations with `track_id` assigned | +| `info.n_tracks` | Number of tracks formed | +| `info.elapsed_s` | Wall time in seconds | ## Parameters ```julia frameconnect(smld; - nnearestclusters = 2, # Clusters for density estimation - nsigmadev = 5.0, # Distance threshold multiplier - maxframegap = 5, # Max frame gap for connections - nmaxnn = 2 # Nearest neighbors for preclustering + n_density_neighbors = 2, # Clusters for density estimation + max_sigma_dist = 5.0, # Distance threshold multiplier + max_frame_gap = 5, # Max frame gap for connections + max_neighbors = 2 # Nearest neighbors for preclustering ) ``` -- `nsigmadev`: Higher values allow connections over larger distances -- `maxframegap`: Increase for dyes with long dark states (dSTORM: 10-20) +- `max_sigma_dist`: Higher values allow connections over larger distances +- `max_frame_gap`: Increase for dyes with long dark states (dSTORM: 10-20) ## API Reference diff --git a/examples/benchmark.jl b/examples/benchmark.jl index ae61d70..960671b 100644 --- a/examples/benchmark.jl +++ b/examples/benchmark.jl @@ -62,20 +62,21 @@ Benchmark a single run with timing and allocation tracking. function benchmark_single(smld; warmup::Bool = false) # Warmup run (if needed) if warmup - frameconnect(smld; nnearestclusters=1) + frameconnect(smld; n_density_neighbors=1) end # Timed run GC.gc() # Clean up before measurement - stats = @timed frameconnect(smld; nnearestclusters=2, nsigmadev=5.0, maxframegap=5, nmaxnn=2) + stats = @timed frameconnect(smld; n_density_neighbors=2, max_sigma_dist=5.0, max_frame_gap=5, max_neighbors=2) + (combined, info) = stats.value return ( time = stats.time, bytes = stats.bytes, gctime = stats.gctime, n_input = length(smld.emitters), - n_output = length(stats.value[3].emitters) + n_output = length(combined.emitters) ) end @@ -156,30 +157,31 @@ function profile_detailed(; n_molecules::Int = 50, n_frames::Int = 100) println("Input localizations: $(length(smld.emitters))") # Warmup - frameconnect(smld; nnearestclusters=1) + frameconnect(smld; n_density_neighbors=1) println("\nRunning with @time macro:") println("-" ^ 70) GC.gc() - @time result = frameconnect( + @time (combined, info) = frameconnect( smld; - nnearestclusters = 2, - nsigmadev = 5.0, - maxframegap = 5, - nmaxnn = 2 + n_density_neighbors = 2, + max_sigma_dist = 5.0, + max_frame_gap = 5, + max_neighbors = 2 ) println("\nOutput:") - println(" Final tracks: $(length(unique(e.track_id for e in result.connected.emitters)))") - println(" Combined emitters: $(length(result.combined.emitters))") + println(" Final tracks: $(info.n_tracks)") + println(" Combined emitters: $(info.n_combined)") + println(" Algorithm time: $(info.elapsed_s)s") # Precision improvement input_σ = mean([e.σ_x for e in smld.emitters]) - combined_σ = mean([e.σ_x for e in result.combined.emitters]) + combined_σ = mean([e.σ_x for e in combined.emitters]) println("\nPrecision improvement: $(round(input_σ / combined_σ, digits=2))x") - return smld, result + return smld, (combined, info) end #= @@ -194,14 +196,14 @@ function allocation_breakdown(; n_molecules::Int = 30, n_frames::Int = 50) println("\nDataset: $(length(smld.emitters)) localizations") # Warmup - frameconnect(smld; nnearestclusters=1) + frameconnect(smld; n_density_neighbors=1) # Prepare params params = SMLMFrameConnection.ParamStruct() - params.nnearestclusters = 2 - params.nsigmadev = 5.0 - params.maxframegap = 5 - params.nmaxnn = 2 + params.n_density_neighbors = 2 + params.max_sigma_dist = 5.0 + params.max_frame_gap = 5 + params.max_neighbors = 2 println("\n@allocated for each step:") println("-" ^ 50) @@ -237,9 +239,9 @@ function allocation_breakdown(; n_molecules::Int = 30, n_frames::Int = 50) # Step 4: Estimate densities GC.gc() alloc_densities = @allocated begin - params.initialdensity = SMLMFrameConnection.estimatedensities(smld_pre, clusterdata, params) + params.initial_density = SMLMFrameConnection.estimatedensities(smld_pre, clusterdata, params) end - params.initialdensity = SMLMFrameConnection.estimatedensities(smld_pre, clusterdata, params) + params.initial_density = SMLMFrameConnection.estimatedensities(smld_pre, clusterdata, params) @printf(" estimatedensities: %10.2f KB\n", alloc_densities / 1024) # Step 5: Connect localizations (LAP) diff --git a/examples/frameconnection_example.jl b/examples/frameconnection_example.jl index 8c3e009..8f0d718 100644 --- a/examples/frameconnection_example.jl +++ b/examples/frameconnection_example.jl @@ -54,7 +54,7 @@ function generate_synthetic_smld(; e = Emitter2DFit{Float64}( obs_x, obs_y, photons, 10.0, # photons, bg - σ_loc, σ_loc, # σ_x, σ_y + σ_loc, σ_loc, 0.0, # σ_x, σ_y, σ_xy σ_photons, 3.0, # σ_photons, σ_bg frame, 1, mol_id, 0 # frame, dataset, track_id, id ) @@ -68,7 +68,7 @@ function generate_synthetic_smld(; # Reset track_id to 0 for input (algorithm will populate it) input_emitters = [Emitter2DFit{Float64}( - e.x, e.y, e.photons, e.bg, e.σ_x, e.σ_y, e.σ_photons, e.σ_bg, + e.x, e.y, e.photons, e.bg, e.σ_x, e.σ_y, e.σ_xy, e.σ_photons, e.σ_bg, e.frame, e.dataset, 0, e.id # track_id = 0 ) for e in emitters] @@ -95,39 +95,40 @@ println(" True molecules: $(length(unique(e.track_id for e in smld_truth.emitte # Run frame connection println("\nRunning frame connection...") -result = frameconnect( +(combined, info) = frameconnect( smld_input; - nnearestclusters = 2, - nsigmadev = 5.0, - maxframegap = 5, - nmaxnn = 2 + n_density_neighbors = 2, + max_sigma_dist = 5.0, + max_frame_gap = 5, + max_neighbors = 2 ) -println(" Combined localizations: $(length(result.combined.emitters))") +println(" Combined localizations: $(length(combined.emitters))") +println(" Time: $(info.elapsed_s) seconds") # Compare with ideal result (using ground truth track_id) println("\nComputing ideal frame connection (ground truth)...") -smld_ideal_connected, smld_ideal_combined = defineidealFC(smld_truth; maxframegap = 5) +smld_ideal_connected, smld_ideal_combined = defineidealFC(smld_truth; max_frame_gap = 5) println(" Ideal combined: $(length(smld_ideal_combined.emitters))") # Print estimated parameters println("\nEstimated photophysics parameters:") -println(" k_on (dark→visible rate): $(round(result.params.k_on, digits=4)) /frame") -println(" k_off (visible→dark rate): $(round(result.params.k_off, digits=4)) /frame") -println(" k_bleach (bleaching rate): $(round(result.params.k_bleach, digits=4)) /frame") -println(" p_miss (miss probability): $(round(result.params.p_miss, digits=4))") +println(" k_on (dark→visible rate): $(round(info.k_on, digits=4)) /frame") +println(" k_off (visible→dark rate): $(round(info.k_off, digits=4)) /frame") +println(" k_bleach (bleaching rate): $(round(info.k_bleach, digits=4)) /frame") +println(" p_miss (miss probability): $(round(info.p_miss, digits=4))") # Compute precision improvement input_σ = mean([e.σ_x for e in smld_input.emitters]) -combined_σ = mean([e.σ_x for e in result.combined.emitters]) +combined_σ = mean([e.σ_x for e in combined.emitters]) println("\nPrecision improvement:") println(" Mean input σ: $(round(input_σ * 1000, digits=1)) nm") println(" Mean combined σ: $(round(combined_σ * 1000, digits=1)) nm") println(" Improvement factor: $(round(input_σ / combined_σ, digits=2))x") # Summary statistics -n_unique_tracks = length(unique(e.track_id for e in smld_connected.emitters)) -avg_locs_per_track = length(smld_connected.emitters) / n_unique_tracks +n_unique_tracks = length(unique(e.track_id for e in info.connected.emitters)) +avg_locs_per_track = length(info.connected.emitters) / n_unique_tracks println("\nConnection statistics:") println(" Unique tracks identified: $n_unique_tracks") println(" Avg localizations per track: $(round(avg_locs_per_track, digits=1))") diff --git a/examples/fullscale_benchmark.jl b/examples/fullscale_benchmark.jl index 55974c7..d17108f 100644 --- a/examples/fullscale_benchmark.jl +++ b/examples/fullscale_benchmark.jl @@ -83,21 +83,19 @@ function benchmark_single_dataset(smld::BasicSMLD) stats = @timed begin frameconnect( smld; - nnearestclusters = 2, - nsigmadev = 5.0, - maxframegap = 10, - nmaxnn = 2 + n_density_neighbors = 2, + max_sigma_dist = 5.0, + max_frame_gap = 10, + max_neighbors = 2 ) end - result = stats.value - n_output = length(result.combined.emitters) - n_tracks = length(unique(e.track_id for e in result.connected.emitters)) + (combined, info) = stats.value return ( n_input = n_input, - n_output = n_output, - n_tracks = n_tracks, + n_output = info.n_combined, + n_tracks = info.n_tracks, time = stats.time, bytes = stats.bytes, gctime = stats.gctime diff --git a/examples/realistic_frameconnection.jl b/examples/realistic_frameconnection.jl index 4f2fb23..eab8847 100644 --- a/examples/realistic_frameconnection.jl +++ b/examples/realistic_frameconnection.jl @@ -101,7 +101,7 @@ function main() # Clear track_id for input (algorithm should estimate it) input_emitters = [Emitter2DFit{Float64}( - e.x, e.y, e.photons, e.bg, e.σ_x, e.σ_y, e.σ_photons, e.σ_bg, + e.x, e.y, e.photons, e.bg, e.σ_x, e.σ_y, e.σ_xy, e.σ_photons, e.σ_bg, e.frame, e.dataset, 0, e.id # track_id = 0 ) for e in smld_noisy.emitters] @@ -119,15 +119,14 @@ function main() stats = @timed begin frameconnect( smld_input; - nnearestclusters = 2, - nsigmadev = 5.0, - maxframegap = 10, # Allow gaps since k_on=0.03 means reblinking - nmaxnn = 2 + n_density_neighbors = 2, + max_sigma_dist = 5.0, + max_frame_gap = 10, # Allow gaps since k_on=0.03 means reblinking + max_neighbors = 2 ) end - result = stats.value - fc_params = result.params + (combined, info) = stats.value @printf(" Time: %.2f seconds\n", stats.time) @printf(" Allocations: %.1f MB\n", stats.bytes / 1e6) @@ -135,15 +134,13 @@ function main() # ========================================================================= # 5. Results summary # ========================================================================= - n_tracks = length(unique(e.track_id for e in result.connected.emitters)) - n_combined = length(result.combined.emitters) - println("\n5. Frame Connection Results") println("-" ^ 70) - @printf(" Input localizations: %d\n", length(smld_input.emitters)) - @printf(" Connected tracks: %d\n", n_tracks) - @printf(" Combined emitters: %d\n", n_combined) - @printf(" Reduction ratio: %.1fx\n", length(smld_input.emitters) / n_combined) + @printf(" Input localizations: %d\n", info.n_input) + @printf(" Connected tracks: %d\n", info.n_tracks) + @printf(" Combined emitters: %d\n", info.n_combined) + @printf(" Preclusters formed: %d\n", info.n_preclusters) + @printf(" Reduction ratio: %.1fx\n", info.n_input / info.n_combined) # ========================================================================= # 6. Compare estimated vs true parameters @@ -159,12 +156,12 @@ function main() println(" Parameter True (Hz) Estimated (/frame) True (/frame)") println(" " * "-" ^ 60) @printf(" k_on %.2f %.4f %.6f\n", - 0.03, fc_params.k_on, true_k_on_pf) + 0.03, info.k_on, true_k_on_pf) @printf(" k_off %.2f %.4f %.4f\n", - 3.0, fc_params.k_off, true_k_off_pf) + 3.0, info.k_off, true_k_off_pf) @printf(" k_bleach N/A %.5f N/A (2-state model)\n", - fc_params.k_bleach) - @printf(" p_miss - %.4f -\n", fc_params.p_miss) + info.k_bleach) + @printf(" p_miss - %.4f -\n", info.p_miss) # ========================================================================= # 7. Precision improvement @@ -173,8 +170,8 @@ function main() println("-" ^ 70) input_σ = mean([e.σ_x for e in smld_input.emitters]) - combined_σ = mean([e.σ_x for e in result.combined.emitters]) - avg_locs_per_track = length(result.connected.emitters) / n_tracks + combined_σ = mean([e.σ_x for e in combined.emitters]) + avg_locs_per_track = info.n_input / info.n_tracks theoretical_improvement = sqrt(avg_locs_per_track) @printf(" Mean input σ: %.1f nm\n", input_σ * 1000) @@ -190,26 +187,26 @@ function main() println("-" ^ 70) # Use original smld_noisy which has true track_id from simulation - smld_ideal_connected, smld_ideal_combined = defineidealFC(smld_noisy; maxframegap = 10) + smld_ideal_connected, smld_ideal_combined = defineidealFC(smld_noisy; max_frame_gap = 10) n_ideal = length(smld_ideal_combined.emitters) @printf(" True emitters (SMLMSim): %d\n", n_true_emitters) @printf(" Ideal combined: %d\n", n_ideal) - @printf(" Algorithm combined: %d\n", n_combined) + @printf(" Algorithm combined: %d\n", info.n_combined) - if n_combined > n_ideal + if info.n_combined > n_ideal @printf(" Over-segmentation: %.1f%% (algorithm found more)\n", - 100.0 * (n_combined - n_ideal) / n_ideal) + 100.0 * (info.n_combined - n_ideal) / n_ideal) else @printf(" Under-segmentation: %.1f%% (algorithm merged too many)\n", - 100.0 * (n_ideal - n_combined) / n_ideal) + 100.0 * (n_ideal - info.n_combined) / n_ideal) end println("\n" * "=" ^ 70) println("Example complete!") println("=" ^ 70) - return smld_input, result, smld_noisy + return smld_input, (combined, info), smld_noisy end # Run diff --git a/src/SMLMFrameConnection.jl b/src/SMLMFrameConnection.jl index b2d34fd..b282f8d 100644 --- a/src/SMLMFrameConnection.jl +++ b/src/SMLMFrameConnection.jl @@ -1,15 +1,18 @@ module SMLMFrameConnection using SMLMData +import SMLMData: AbstractSMLMConfig, AbstractSMLMInfo using Hungarian using NearestNeighbors using Optim using StatsBase using Statistics -export frameconnect, defineidealFC, combinelocalizations, ParamStruct +export frameconnect, defineidealFC, combinelocalizations +export FrameConnectConfig, FrameConnectInfo, ParamStruct include("structdefinitions.jl") +include("connectinfo.jl") include("precluster.jl") include("defineidealFC.jl") include("organizeclusters.jl") diff --git a/src/connectinfo.jl b/src/connectinfo.jl new file mode 100644 index 0000000..b98474e --- /dev/null +++ b/src/connectinfo.jl @@ -0,0 +1,52 @@ +# FrameConnectInfo struct for tuple-pattern return + +""" + FrameConnectInfo{T} + +Secondary output from `frameconnect()` containing track assignments and algorithm metadata. + +# Fields +- `connected::BasicSMLD{T}`: Input SMLD with track_id assigned (localizations uncombined) +- `n_input::Int`: Number of input localizations +- `n_tracks::Int`: Number of tracks formed +- `n_combined::Int`: Number of output localizations +- `k_on::Float64`: Estimated on rate (1/frame) +- `k_off::Float64`: Estimated off rate (1/frame) +- `k_bleach::Float64`: Estimated bleach rate (1/frame) +- `p_miss::Float64`: Probability of missed detection +- `initial_density::Vector{Float64}`: Initial density estimate per cluster (emitters/μm²) +- `elapsed_s::Float64`: Wall time in seconds +- `algorithm::Symbol`: Algorithm used (`:lap`) +- `n_preclusters::Int`: Number of preclusters formed + +# Rate Parameter Interpretation + +The rate parameters describe the photophysics of blinking fluorophores: +- `k_on`: Rate at which dark emitters convert to visible state +- `k_off`: Rate at which visible emitters convert to reversible dark state +- `k_bleach`: Rate at which visible emitters are irreversibly photobleached +- Duty cycle = k_on / (k_on + k_off) +- For typical dSTORM: k_on << k_off (low duty cycle, mostly dark with brief blinks) + +# Example +```julia +(combined, info) = frameconnect(smld) +println("Connected \$(info.n_input) → \$(info.n_combined) localizations") +println("Formed \$(info.n_tracks) tracks in \$(info.elapsed_s)s") +# Access track assignments via info.connected +``` +""" +struct FrameConnectInfo{T} <: AbstractSMLMInfo + connected::BasicSMLD{T} + n_input::Int + n_tracks::Int + n_combined::Int + k_on::Float64 + k_off::Float64 + k_bleach::Float64 + p_miss::Float64 + initial_density::Vector{Float64} + elapsed_s::Float64 + algorithm::Symbol + n_preclusters::Int +end diff --git a/src/connectlocalizations.jl b/src/connectlocalizations.jl index df9bc8b..62f957c 100644 --- a/src/connectlocalizations.jl +++ b/src/connectlocalizations.jl @@ -1,8 +1,8 @@ using StatsBase """ - connectlocalizations!(connectID::Vector{Int64}, - clusterdata::Vector{Matrix{Float32}}, + connectlocalizations!(connectID::Vector{Int64}, + clusterdata::Vector{Matrix{Float32}}, params::ParamStruct, nframes::Int64) Connect localizations in `clusterdata` by solving a linear assignment problem. diff --git a/src/create_costmatrix.jl b/src/create_costmatrix.jl index 656aa38..187f912 100644 --- a/src/create_costmatrix.jl +++ b/src/create_costmatrix.jl @@ -31,8 +31,8 @@ function create_costmatrix(clusterdata::Vector{Matrix{Float32}}, k_off = params.k_off k_bleach = params.k_bleach p_miss = params.p_miss - rho_0 = params.initialdensity[clusterind] - maxframegap = params.maxframegap + rho_0 = params.initial_density[clusterind] + max_frame_gap = params.max_frame_gap # Populate the upper-left "connection" block. nlocalizations = length(framenum) @@ -89,8 +89,8 @@ function create_costmatrix(clusterdata::Vector{Matrix{Float32}}, indices = 1:nlocalizations framesint = Int32.(framenum) for nn in indices - deltaframe_past = minimum([maxframegap; framesint[nn]-startframe]) - deltaframe_future = minimum([maxframegap; nframes-framesint[nn]]) + deltaframe_past = minimum([max_frame_gap; framesint[nn]-startframe]) + deltaframe_future = minimum([max_frame_gap; nframes-framesint[nn]]) # Localization uncertainty ellipse area (μm²) to make density dimensionless # Area = π√det(Σ) where det(Σ) = σ_x²σ_y² - σ_xy² det_loc = x_se[nn]^2 * y_se[nn]^2 - xy_cov[nn]^2 diff --git a/src/defineidealFC.jl b/src/defineidealFC.jl index c41271d..abc73ba 100644 --- a/src/defineidealFC.jl +++ b/src/defineidealFC.jl @@ -3,7 +3,7 @@ using SMLMData """ connect1DS(smld::BasicSMLD{T,E}, dataset::Int, connectID::Vector{Int}, maxID::Int; - maxframegap::Int = 5) where {T, E<:SMLMData.AbstractEmitter} + max_frame_gap::Int = 5) where {T, E<:SMLMData.AbstractEmitter} Define the "ideal" frame-connection result for simulated `smld` with one dataset. @@ -17,7 +17,7 @@ call defineidealFC(). - `dataset`: Dataset number to be connected. - `connectID`: Current connection IDs for all emitters. - `maxID`: Current maximum ID value. -- `maxframegap`: Maximum frame gap allowed between localizations connected in +- `max_frame_gap`: Maximum frame gap allowed between localizations connected in the "ideal" result. # Returns @@ -25,7 +25,7 @@ call defineidealFC(). """ function connect1DS(smld::BasicSMLD{T,E}, dataset::Int, connectID::Vector{Int}, maxID::Int; - maxframegap::Int = 5) where {T, E<:SMLMData.AbstractEmitter} + max_frame_gap::Int = 5) where {T, E<:SMLMData.AbstractEmitter} emitters = smld.emitters # Find indices for current dataset @@ -58,13 +58,13 @@ function connect1DS(smld::BasicSMLD{T,E}, dataset::Int, # If all of these localizations are within the framegap, no action is # needed (they already share the same connectID). framediff = diff(sortedframes) - if all(framediff .<= maxframegap) + if all(framediff .<= max_frame_gap) continue end # Determine which localizations we can combine. for ff = 1:length(framediff) - if framediff[ff] <= maxframegap + if framediff[ff] <= max_frame_gap # Connect these localizations. connectID_DS[sorted_local_inds[ff+1]] = connectID_DS[sorted_local_inds[ff]] @@ -85,7 +85,7 @@ end """ smld_connected, smld_combined = defineidealFC( smld::BasicSMLD{T,E}; - maxframegap::Int = 5) where {T, E<:SMLMData.AbstractEmitter} + max_frame_gap::Int = 5) where {T, E<:SMLMData.AbstractEmitter} Define the "ideal" frame-connection result for a simulated `smld`. @@ -94,15 +94,15 @@ This function defines the "ideal" frame connection result from a simulation. That is to say, for a simulated BasicSMLD structure `smld` with `track_id` field populated to indicate emitter membership of localizations, this function will generate an "ideal" FC result which combines all blinking events that appeared -with frame gaps less than `maxframegap` of one another. Note that for very +with frame gaps less than `max_frame_gap` of one another. Note that for very high duty cycles, multiple blinking events might be mistakingly combined by -this method (i.e., if the emitter blinks back on within `maxframegap` frames +this method (i.e., if the emitter blinks back on within `max_frame_gap` frames of its previous blink). Note that localizations are not allowed to be connected across datasets. # Inputs - `smld`: BasicSMLD with track_id populated to indicate emitter ID. -- `maxframegap`: Maximum frame gap allowed between localizations connected in +- `max_frame_gap`: Maximum frame gap allowed between localizations connected in the "ideal" result. # Outputs @@ -111,7 +111,7 @@ connected across datasets. - `smld_combined`: Ideal frame-connection result with localizations combined. """ function defineidealFC(smld::BasicSMLD{T,E}; - maxframegap::Int = 5) where {T, E<:SMLMData.AbstractEmitter} + max_frame_gap::Int = 5) where {T, E<:SMLMData.AbstractEmitter} emitters = smld.emitters # Initialize connectID from track_id @@ -123,7 +123,7 @@ function defineidealFC(smld::BasicSMLD{T,E}; # Loop through datasets and combine localizations as appropriate. for ds in datasets - connectID, maxID = connect1DS(smld, ds, connectID, maxID; maxframegap = maxframegap) + connectID, maxID = connect1DS(smld, ds, connectID, maxID; max_frame_gap = max_frame_gap) end # Compress connectID diff --git a/src/estimatedensities.jl b/src/estimatedensities.jl index 72250a7..23d2cb4 100644 --- a/src/estimatedensities.jl +++ b/src/estimatedensities.jl @@ -2,13 +2,13 @@ using SMLMData using NearestNeighbors """ - initialdensity = estimatedensities(smld::BasicSMLD{T,E}, + initial_density = estimatedensities(smld::BasicSMLD{T,E}, clusterdata::Vector{Matrix{Float32}}, params::ParamStruct) where {T, E<:SMLMData.AbstractEmitter} Estimate local emitter densities for clusters in `smld` and `clusterdata`. # Description -The initial local densities `initialdensity` around each pre-cluster present in +The initial local densities `initial_density` around each pre-cluster present in `smld`/`clusterdata` are estimated based on the local density of pre-clusters throughout the entire set of data as well as some of the rate parameters provided in `params`. @@ -31,15 +31,15 @@ function estimatedensities(smld::BasicSMLD{T,E}, # If only one cluster is present, we should return an answer right away and stop. if nclusters == 1 if size(clusterdata[1], 1) == 1 - initialdensity = 1.0 + initial_density = 1.0 else clusterarea = (maximum(clusterdata[1][:, 1]) - minimum(clusterdata[1][:, 1])) * (maximum(clusterdata[1][:, 2]) - minimum(clusterdata[1][:, 2])) - initialdensity = (1 / clusterarea) * + initial_density = (1 / clusterarea) * ((params.k_bleach / params.k_off) / (1 - params.p_miss)) / (1 - exp(-params.k_bleach * dutycycle * (maxframe - 1))) end - return [initialdensity] # Return as vector to match ParamStruct.initialdensity type + return [initial_density] # Return as vector to match ParamStruct.initial_density type end # Determine the center of all clusters assuming each arose from the same @@ -79,7 +79,7 @@ function estimatedensities(smld::BasicSMLD{T,E}, # Estimate the local cluster density based on the distance to # nearest-neighbors. - kneighbors = minimum([params.nnearestclusters; nclusters - 1]) + kneighbors = minimum([params.n_density_neighbors; nclusters - 1]) kneighbors = max(kneighbors, 1) # Ensure at least 1 neighbor # Ensure we don't request more neighbors than available (including self) kneighbors = min(kneighbors, nclusters - 1) @@ -109,10 +109,10 @@ function estimatedensities(smld::BasicSMLD{T,E}, # Estimate the density of underlying emitters based on cluster density. lambda1 = params.k_bleach * dutycycle lambda2 = (params.k_on + params.k_off + params.k_bleach) - lambda1 - initialdensity = clusterdensity * + initial_density = clusterdensity * (1.0 / dutycycle) * (1.0 / params.k_off) * (1.0 / (1.0 - params.p_miss)) ./ ((1.0 / lambda1) * (1.0 - exp(-lambda1 * (maxframe - 1.0))) - (1.0 / lambda2) * (1.0 - exp(-lambda2 * (maxframe - 1.0)))) - return initialdensity + return initial_density end diff --git a/src/frameconnect.jl b/src/frameconnect.jl index b33cfea..b697926 100644 --- a/src/frameconnect.jl +++ b/src/frameconnect.jl @@ -1,65 +1,83 @@ using SMLMData """ - result = frameconnect(smld::BasicSMLD{T,E}; - nnearestclusters::Int=2, nsigmadev::Float64=5.0, - maxframegap::Int=5, nmaxnn::Int=2) where {T, E<:SMLMData.AbstractEmitter} + (combined, info) = frameconnect(smld::BasicSMLD, config::FrameConnectConfig) + (combined, info) = frameconnect(smld::BasicSMLD; kwargs...) Connect repeated localizations of the same emitter in `smld`. # Description Repeated localizations of the same emitter present in `smld` are connected and -combined into higher precision localizations of that emitter. This is done by +combined into higher precision localizations of that emitter. This is done by 1) forming pre-clusters of localizations, 2) estimating rate parameters from the pre-clusters, 3) solving a linear assignment problem for connecting localizations in each pre-cluster, and 4) combining the connected localizations using their MLE position estimate assuming Gaussian noise. -# Inputs -- `smld`: BasicSMLD containing the localizations that should be connected. - Must contain emitters with valid position uncertainties (σ_x, σ_y). -- `nnearestclusters`: Number of nearest preclusters used for local density - estimates. (default = 2)(see estimatedensities()) -- `nsigmadev`: Multiplier of localization errors that defines a pre-clustering - distance threshold. (default = 5)(see precluster())(microns) -- `maxframegap`: Maximum frame gap between temporally adjacent localizations in - a precluster. (default = 5)(see precluster())(frames) -- `nmaxnn`: Maximum number of nearest-neighbors inspected for precluster - membership. Ideally, this would be set to inf, but that's not - feasible for most data. (default = 2)(see precluster()) - -# Outputs -Returns a NamedTuple with fields: -- `combined`: Final frame-connection result (i.e., `smld` with localizations that - seem to be from the same blinking event combined into higher - precision localizations). -- `connected`: Input `smld` with track_id field updated to reflect connected - localizations (localizations remain uncombined). -- `params`: Structure of parameters used in the algorithm, with some copied - directly from the option kwargs to this function, and others - calculated internally (see SMLMFrameConnection.ParamStruct). +# Arguments +- `smld::BasicSMLD`: Localizations to connect. Must contain emitters with valid + position uncertainties (σ_x, σ_y). +- `config::FrameConnectConfig`: Configuration parameters (optional, can use kwargs instead) + +# Keyword Arguments (equivalent to FrameConnectConfig fields) +- `n_density_neighbors::Int=2`: Number of nearest preclusters used for local density + estimates (see `estimatedensities`) +- `max_sigma_dist::Float64=5.0`: Multiplier of localization errors that defines a + pre-clustering distance threshold (see `precluster`) +- `max_frame_gap::Int=5`: Maximum frame gap between temporally adjacent localizations + in a precluster (see `precluster`) +- `max_neighbors::Int=2`: Maximum number of nearest-neighbors inspected for precluster + membership (see `precluster`) + +# Returns +A tuple `(combined, info)`: +- `combined::BasicSMLD`: Connected localizations combined into higher precision results +- `info::FrameConnectInfo`: Track assignments and algorithm metadata (see [`FrameConnectInfo`](@ref)) + +# Example +```julia +# Using kwargs (most common) +(combined, info) = frameconnect(smld) +(combined, info) = frameconnect(smld; max_frame_gap=10) + +# Using config struct +config = FrameConnectConfig(max_frame_gap=10, max_sigma_dist=3.0) +(combined, info) = frameconnect(smld, config) + +println("Connected \$(info.n_input) → \$(info.n_combined) localizations") +println("Formed \$(info.n_tracks) tracks from \$(info.n_preclusters) preclusters") + +# Access track assignments for downstream analysis +track_ids = [e.track_id for e in info.connected.emitters] +``` """ -function frameconnect(smld::BasicSMLD{T,E}; - nnearestclusters::Int = 2, nsigmadev::Float64 = 5.0, - maxframegap::Int = 5, nmaxnn::Int = 2) where {T, E<:SMLMData.AbstractEmitter} +function frameconnect(smld::BasicSMLD{T,E}; kwargs...) where {T, E<:SMLMData.AbstractEmitter} + # kwargs form forwards to config form + config = FrameConnectConfig(; kwargs...) + return frameconnect(smld, config) +end + +function frameconnect(smld::BasicSMLD{T,E}, config::FrameConnectConfig) where {T, E<:SMLMData.AbstractEmitter} + t_start = time() # Prepare a ParamStruct to keep track of parameters used. params = ParamStruct() - params.nnearestclusters = nnearestclusters - params.nsigmadev = nsigmadev - params.maxframegap = maxframegap - params.nmaxnn = nmaxnn + params.n_density_neighbors = config.n_density_neighbors + params.max_sigma_dist = config.max_sigma_dist + params.max_frame_gap = config.max_frame_gap + params.max_neighbors = config.max_neighbors # Generate pre-clusters of localizations in `smld`. smld_preclustered = precluster(smld, params) clusterdata = organizeclusters(smld_preclustered) + n_preclusters = length(clusterdata) # Estimate rate parameters. params.k_on, params.k_off, params.k_bleach, params.p_miss = estimateparams(smld_preclustered, clusterdata) # Estimate the underlying density of emitters. - params.initialdensity = + params.initial_density = estimatedensities(smld_preclustered, clusterdata, params) # Get nframes @@ -89,5 +107,24 @@ function frameconnect(smld::BasicSMLD{T,E}; # Combine the connected localizations into higher precision localizations. smld_combined = combinelocalizations(smld_connected) - return (combined=smld_combined, connected=smld_connected, params=params) + elapsed_s = time() - t_start + + # Build FrameConnectInfo + n_tracks = length(unique(connectID_final)) + info = FrameConnectInfo{T}( + smld_connected, + length(smld.emitters), + n_tracks, + length(smld_combined.emitters), + params.k_on, + params.k_off, + params.k_bleach, + params.p_miss, + params.initial_density, + elapsed_s, + :lap, + n_preclusters + ) + + return (smld_combined, info) end diff --git a/src/precluster.jl b/src/precluster.jl index 4f2e5b7..8938d59 100644 --- a/src/precluster.jl +++ b/src/precluster.jl @@ -11,8 +11,8 @@ Cluster localizations in `smld` based on distance and time thresholds in `params # Description Localizations in the input structure `smld` are clustered together based on their spatiotemporal separations. All localizations within a spatial -threshold of `params.nsigmadev*mean([σ_x σ_y])` and a temporal -threshold of `params.maxframegap` of one another will be clustered together, +threshold of `params.max_sigma_dist*mean([σ_x σ_y])` and a temporal +threshold of `params.max_frame_gap` of one another will be clustered together, meaning that these localizations now share the same unique integer value for their track_id field. @@ -46,9 +46,9 @@ function precluster(smld::BasicSMLD{T,E}, mean_se = Statistics.mean([σ_x σ_y]; dims = 2) # Isolate some parameters from params. - maxframegap = params.maxframegap - nsigmadev = params.nsigmadev - nmaxnn = params.nmaxnn + max_frame_gap = params.max_frame_gap + max_sigma_dist = params.max_sigma_dist + max_neighbors = params.max_neighbors # Initialize a connectID array, with each localization being considered # a unique cluster. @@ -81,7 +81,7 @@ function precluster(smld::BasicSMLD{T,E}, # Determine which localizations should be considered for # clustering. currentindFN = (1:nperframe[ff]) .+ ncumulativeFN[ff] - candidateind = findall((framenumCDS .>= (frames[ff] .- maxframegap)) .& + candidateind = findall((framenumCDS .>= (frames[ff] .- max_frame_gap)) .& (framenumCDS .<= frames[ff])) if length(candidateind) < 2 maxID += 1 @@ -94,14 +94,14 @@ function precluster(smld::BasicSMLD{T,E}, # which is a candidate for clustering. kdtree = NearestNeighbors.KDTree(xyCDS[:, candidateind]) nnindices, nndist = NearestNeighbors.knn(kdtree, xyCDS[:, currentindFN], - min(nmaxnn + 1, length(candidateind)), true) + min(max_neighbors + 1, length(candidateind)), true) # Assign localizations to clusters based on `nndist`. for ii in 1:nperframe[ff] # Determine which candidates meet our distance cutoff. se_sum = mean_seCDS[currentindFN[ii]] .+ mean_seCDS[candidateind[nnindices[ii]]] - validnninds = nnindices[ii][nndist[ii].<=(nsigmadev*se_sum)] + validnninds = nnindices[ii][nndist[ii].<=(max_sigma_dist*se_sum)] # Update connectIDCDS to reflect the new clusters. updateinds = unique([currentindFN[ii] diff --git a/src/structdefinitions.jl b/src/structdefinitions.jl index a652b3d..caf3bd0 100644 --- a/src/structdefinitions.jl +++ b/src/structdefinitions.jl @@ -1,26 +1,62 @@ # This file defines some struct types used in the FrameConnection package. """ - ParamStruct(initialdensity::Vector{Float64}, - k_on::Float64, k_off::Float64, k_bleach::Float64, p_miss::Float64, - nsigmadev::Float64, maxframegap::Int, nnearestclusters::Int) + FrameConnectConfig + +Configuration parameters for frame connection algorithm. + +# Fields +- `n_density_neighbors::Int=2`: Number of nearest preclusters used for local density + estimates (see `estimatedensities`) +- `max_sigma_dist::Float64=5.0`: Multiplier of localization errors that defines a + pre-clustering distance threshold (see `precluster`) +- `max_frame_gap::Int=5`: Maximum frame gap between temporally adjacent localizations + in a precluster (see `precluster`) +- `max_neighbors::Int=2`: Maximum number of nearest-neighbors inspected for precluster + membership (see `precluster`) + +# Example +```julia +# Using default config +config = FrameConnectConfig() +(combined, info) = frameconnect(smld, config) + +# Custom config +config = FrameConnectConfig(max_frame_gap=10, max_sigma_dist=3.0) +(combined, info) = frameconnect(smld, config) + +# Kwargs form (equivalent to Config form) +(combined, info) = frameconnect(smld; max_frame_gap=10, max_sigma_dist=3.0) +``` +""" +Base.@kwdef struct FrameConnectConfig <: AbstractSMLMConfig + n_density_neighbors::Int = 2 + max_sigma_dist::Float64 = 5.0 + max_frame_gap::Int = 5 + max_neighbors::Int = 2 +end + +""" + ParamStruct(initial_density::Vector{Float64}, + k_on::Float64, k_off::Float64, k_bleach::Float64, p_miss::Float64, + max_sigma_dist::Float64, max_frame_gap::Int, n_density_neighbors::Int) Structure of parameters needed for frame-connection. # Fields -- `initialdensity`: Density of emitters at the start of the experiment. +- `initial_density`: Density of emitters at the start of the experiment. (see estimatedensities()) (emitters/μm²) -- `nnearestclusters`: Number of nearest preclusters used for local density +- `n_density_neighbors`: Number of nearest preclusters used for local density estimates. (default = 2)(see estimatedensities()) - `k_on`: Rate at which dark emitters convert to the visible state (1/frame). - `k_off`: Rate at which visible emitters convert to the reversible dark state (1/frame). - `k_bleach`: Rate at which visible emitters are irreversibly photobleached (1/frame). - `p_miss`: Probability of missing a localization of a visible emitter. -- `nsigmadev`: Multiplier of localization errors that defines a pre-clustering +- `max_sigma_dist`: Multiplier of localization errors that defines a pre-clustering distance threshold. (default = 5)(see precluster())(unitless) -- `maxframegap`: Maximum frame gap between temporally adjacent localizations in +- `max_frame_gap`: Maximum frame gap between temporally adjacent localizations in a precluster. (default = 5)(see precluster())(frames) -- `nmaxnn`: Maximum number of nearest-neighbors inspected for precluster +- `max_neighbors`: Maximum number of nearest-neighbors inspected for precluster membership. Ideally, this would be set to inf, but that's not feasible for most data. (default = 2)(see precluster()) @@ -29,15 +65,15 @@ The duty cycle (fraction of time in ON state) is k_on/(k_on + k_off). For typica dSTORM, k_on << k_off gives low duty cycle (mostly dark, brief blinks). """ mutable struct ParamStruct - initialdensity::Vector{Float64} - nnearestclusters::Int + initial_density::Vector{Float64} + n_density_neighbors::Int k_on::Float64 k_off::Float64 k_bleach::Float64 p_miss::Float64 - nsigmadev::Float64 - maxframegap::Int - nmaxnn::Int + max_sigma_dist::Float64 + max_frame_gap::Int + max_neighbors::Int end ParamStruct() = ParamStruct([], 2, 0.0, 0.0, 0.0, 0.0, 5.0, 5, 2) diff --git a/test/test_defineidealFC.jl b/test/test_defineidealFC.jl index 6c73775..4f4cfac 100644 --- a/test/test_defineidealFC.jl +++ b/test/test_defineidealFC.jl @@ -18,14 +18,14 @@ ] smld = make_test_smld(emitters; n_frames=3) - smld_connected, smld_combined = defineidealFC(smld; maxframegap=5) + smld_connected, smld_combined = defineidealFC(smld; max_frame_gap=5) # Should have 2 combined localizations (one per unique emitter) @test length(smld_combined.emitters) == 2 end - @testset "respects maxframegap" begin - # Same emitter (track_id=1) but with frame gap > maxframegap + @testset "respects max_frame_gap" begin + # Same emitter (track_id=1) but with frame gap > max_frame_gap emitters = [ make_emitter(5.0, 5.0, 1; track_id=1), make_emitter(5.0, 5.0, 2; track_id=1), @@ -34,8 +34,8 @@ ] smld = make_test_smld(emitters; n_frames=11) - # With maxframegap=5, should split into 2 blinking events - smld_connected, smld_combined = defineidealFC(smld; maxframegap=5) + # With max_frame_gap=5, should split into 2 blinking events + smld_connected, smld_combined = defineidealFC(smld; max_frame_gap=5) @test length(smld_combined.emitters) == 2 end @@ -45,7 +45,7 @@ emitters = [make_emitter(5.0, 5.0, i; track_id=1) for i in 1:5] smld = make_test_smld(emitters; n_frames=5) - _, smld_combined = defineidealFC(smld; maxframegap=5) + _, smld_combined = defineidealFC(smld; max_frame_gap=5) @test length(smld_combined.emitters) == 1 end @@ -60,7 +60,7 @@ ] smld = make_test_smld(emitters; n_frames=4) - _, smld_combined = defineidealFC(smld; maxframegap=5) + _, smld_combined = defineidealFC(smld; max_frame_gap=5) # Should have 2 combined localizations (one per emitter identity) @test length(smld_combined.emitters) == 2 diff --git a/test/test_frameconnect.jl b/test/test_frameconnect.jl index 7ac4794..f9ce90e 100644 --- a/test/test_frameconnect.jl +++ b/test/test_frameconnect.jl @@ -1,19 +1,23 @@ @testset "frameconnect" begin - @testset "return types" begin + @testset "return types - tuple pattern" begin # Use multiple molecules for proper density estimation emitters = vcat( [make_emitter(Float64(i), Float64(i), j) for i in 1:5 for j in 1:3]... ) smld = make_test_smld(emitters; n_frames=3) - result = frameconnect(smld) - - @test result isa NamedTuple - @test haskey(result, :combined) - @test haskey(result, :connected) - @test haskey(result, :params) - @test result.combined isa BasicSMLD - @test result.connected isa BasicSMLD - @test result.params isa ParamStruct + (combined, info) = frameconnect(smld) + + @test combined isa BasicSMLD + @test info isa FrameConnectInfo + + # FrameConnectInfo fields + @test info.connected isa BasicSMLD + @test info.n_input == length(smld.emitters) + @test info.n_tracks > 0 + @test info.n_combined == length(combined.emitters) + @test info.elapsed_s > 0 + @test info.algorithm == :lap + @test info.n_preclusters > 0 end @testset "single molecule connection" begin @@ -25,11 +29,13 @@ ) smld = make_test_smld(emitters; n_frames=5) - result = frameconnect(smld; maxframegap=5, nnearestclusters=1) + (combined, info) = frameconnect(smld; max_frame_gap=5, n_density_neighbors=1) # Target molecule should be combined (may have 5 background singles) # At minimum, we should have fewer emitters than input - @test length(result.combined.emitters) < length(smld.emitters) + @test length(combined.emitters) < length(smld.emitters) + @test info.n_input == length(smld.emitters) + @test info.n_combined == length(combined.emitters) end @testset "multiple separated molecules" begin @@ -42,11 +48,11 @@ ) smld = make_test_smld(emitters; n_frames=3) - result = frameconnect(smld; maxframegap=5) + (combined, info) = frameconnect(smld; max_frame_gap=5) # Should combine into ~4 molecules (fewer than 12 input emitters) - @test length(result.combined.emitters) <= 8 - @test length(result.combined.emitters) < length(smld.emitters) + @test length(combined.emitters) <= 8 + @test length(combined.emitters) < length(smld.emitters) end @testset "track_id population" begin @@ -57,15 +63,15 @@ ) smld = make_test_smld(emitters; n_frames=3) - result = frameconnect(smld) + (combined, info) = frameconnect(smld) - # All emitters should have non-zero track_id - for e in result.connected.emitters + # All emitters in info.connected should have non-zero track_id + for e in info.connected.emitters @test e.track_id > 0 end end - @testset "parameter estimation" begin + @testset "rate parameter estimation" begin emitters = vcat( make_blinking_molecule(5.0, 5.0, [1, 2, 3]), make_blinking_molecule(10.0, 10.0, [1, 2, 3]), @@ -74,17 +80,17 @@ ) smld = make_test_smld(emitters; n_frames=3) - result = frameconnect(smld) + (combined, info) = frameconnect(smld) # Parameters should be estimated (non-default values) - @test result.params.k_on >= 0 - @test result.params.k_off >= 0 - @test result.params.k_bleach >= 0 - @test 0 <= result.params.p_miss <= 1 - @test !isempty(result.params.initialdensity) + @test info.k_on >= 0 + @test info.k_off >= 0 + @test info.k_bleach >= 0 + @test 0 <= info.p_miss <= 1 + @test !isempty(info.initial_density) end - @testset "respects maxframegap" begin + @testset "respects max_frame_gap" begin # Create molecules with gaps, plus background for density estimation emitters = vcat( make_blinking_molecule(5.0, 5.0, [1, 2, 3]), # Frames 1-3 @@ -94,11 +100,11 @@ ) smld = make_test_smld(emitters; n_frames=12) - # With maxframegap=5, the two groups at (5,5) should NOT connect - result = frameconnect(smld; maxframegap=5) + # With max_frame_gap=5, the two groups at (5,5) should NOT connect + (combined, info) = frameconnect(smld; max_frame_gap=5) # Should have at least 3 combined localizations (2 for split molecule + 2 background) - @test length(result.combined.emitters) >= 3 + @test length(combined.emitters) >= 3 end @testset "single localization" begin @@ -106,20 +112,22 @@ emitters = [make_emitter(5.0, 5.0, 1)] smld = make_test_smld(emitters; n_frames=1) - result = frameconnect(smld) + (combined, info) = frameconnect(smld) - @test length(result.combined.emitters) == 1 - @test result.combined.emitters[1].x ≈ 5.0 atol=1e-10 + @test length(combined.emitters) == 1 + @test combined.emitters[1].x ≈ 5.0 atol=1e-10 + @test info.n_input == 1 + @test info.n_combined == 1 end @testset "preserves camera and metadata" begin smld = make_single_molecule_smld() - result = frameconnect(smld) + (combined, info) = frameconnect(smld) - @test result.connected.camera == smld.camera - @test result.combined.camera == smld.camera - @test result.connected.n_datasets == smld.n_datasets + @test info.connected.camera == smld.camera + @test combined.camera == smld.camera + @test info.connected.n_datasets == smld.n_datasets end @testset "precision improvement" begin @@ -136,10 +144,10 @@ emitters = vcat(target_emitters, background) smld = make_test_smld(emitters; n_frames=n_locs) - result = frameconnect(smld; maxframegap=n_locs, nnearestclusters=1) + (combined, info) = frameconnect(smld; max_frame_gap=n_locs, n_density_neighbors=1) # Find the combined target molecule (around position 5,5) - target_combined = filter(e -> e.x < 10.0 && e.y < 10.0, result.combined.emitters) + target_combined = filter(e -> e.x < 10.0 && e.y < 10.0, combined.emitters) @test length(target_combined) == 1 σ_output = target_combined[1].σ_x @@ -149,19 +157,31 @@ @test σ_output ≈ expected_σ rtol=0.1 end - @testset "custom parameters" begin + @testset "timing capture" begin smld = make_single_molecule_smld() - result = frameconnect(smld; - nnearestclusters=3, - nsigmadev=4.0, - maxframegap=10, - nmaxnn=3 - ) + (combined, info) = frameconnect(smld) + + # elapsed_s should be positive and reasonable (< 60 seconds) + @test info.elapsed_s > 0 + @test info.elapsed_s < 60.0 + end - @test result.params.nnearestclusters == 3 - @test result.params.nsigmadev == 4.0 - @test result.params.maxframegap == 10 - @test result.params.nmaxnn == 3 + @testset "FrameConnectInfo type parameter" begin + # Test with Float64 + smld64 = make_single_molecule_smld() + (_, info64) = frameconnect(smld64) + @test info64 isa FrameConnectInfo{Float64} + + # Test with Float32 + emitters32 = [SMLMData.Emitter2DFit{Float32}( + 5.0f0, 5.0f0, 1000.0f0, 10.0f0, + 0.02f0, 0.02f0, 0.0f0, 10.0f0, 1.0f0, + 1, 1, 0, 1 + )] + camera32 = SMLMData.IdealCamera(1:64, 1:64, 0.1f0) + smld32 = BasicSMLD(emitters32, camera32, 1, 1, Dict{String,Any}()) + (_, info32) = frameconnect(smld32) + @test info32 isa FrameConnectInfo{Float32} end end diff --git a/test/test_types.jl b/test/test_types.jl index 178fded..2d383c0 100644 --- a/test/test_types.jl +++ b/test/test_types.jl @@ -1,17 +1,44 @@ @testset "Types" begin + @testset "FrameConnectInfo" begin + # Create minimal test data + emitters = [SMLMData.Emitter2DFit{Float64}( + 5.0, 5.0, 1000.0, 10.0, 0.02, 0.02, 0.0, 10.0, 1.0, 1, 1, 1, 1 + )] + camera = SMLMData.IdealCamera(1:64, 1:64, 0.1) + smld = BasicSMLD(emitters, camera, 1, 1, Dict{String,Any}()) + + info = FrameConnectInfo{Float64}( + smld, 10, 5, 5, 0.1, 0.5, 0.01, 0.05, [1.0, 2.0], 1.5, :lap, 3 + ) + + @test info isa FrameConnectInfo{Float64} + @test info.connected === smld + @test info.n_input == 10 + @test info.n_tracks == 5 + @test info.n_combined == 5 + @test info.k_on == 0.1 + @test info.k_off == 0.5 + @test info.k_bleach == 0.01 + @test info.p_miss == 0.05 + @test info.initial_density == [1.0, 2.0] + @test info.elapsed_s == 1.5 + @test info.algorithm == :lap + @test info.n_preclusters == 3 + end + @testset "ParamStruct" begin @testset "default constructor" begin params = ParamStruct() @test params isa ParamStruct - @test params.initialdensity == [] - @test params.nnearestclusters == 2 + @test params.initial_density == [] + @test params.n_density_neighbors == 2 @test params.k_on == 0.0 @test params.k_off == 0.0 @test params.k_bleach == 0.0 @test params.p_miss == 0.0 - @test params.nsigmadev == 5.0 - @test params.maxframegap == 5 - @test params.nmaxnn == 2 + @test params.max_sigma_dist == 5.0 + @test params.max_frame_gap == 5 + @test params.max_neighbors == 2 end @testset "field mutation" begin @@ -24,21 +51,21 @@ params.k_off = 0.5 @test params.k_off == 0.5 - params.initialdensity = [1.0, 2.0] - @test params.initialdensity == [1.0, 2.0] + params.initial_density = [1.0, 2.0] + @test params.initial_density == [1.0, 2.0] end @testset "full constructor" begin params = ParamStruct([1.0], 3, 0.1, 0.5, 0.01, 0.05, 4.0, 10, 3) - @test params.initialdensity == [1.0] - @test params.nnearestclusters == 3 + @test params.initial_density == [1.0] + @test params.n_density_neighbors == 3 @test params.k_on == 0.1 @test params.k_off == 0.5 @test params.k_bleach == 0.01 @test params.p_miss == 0.05 - @test params.nsigmadev == 4.0 - @test params.maxframegap == 10 - @test params.nmaxnn == 3 + @test params.max_sigma_dist == 4.0 + @test params.max_frame_gap == 10 + @test params.max_neighbors == 3 end end end