diff --git a/CLAUDE.md b/CLAUDE.md
new file mode 100644
index 0000000..c9d15ee
--- /dev/null
+++ b/CLAUDE.md
@@ -0,0 +1,86 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## Project Overview
+
+SMLMFrameConnection.jl performs frame-connection on 2D localization microscopy data, combining repeated localizations of blinking fluorophores into higher-precision localizations. Part of the JuliaSMLM ecosystem.
+
+**Reference:** Schodt & Lidke (2021), "Spatiotemporal Clustering of Repeated Super-Resolution Localizations via Linear Assignment Problem"
+
+## Commands
+
+```bash
+# Run tests
+julia --project -e 'using Pkg; Pkg.test()'
+
+# Run a specific test file
+julia --project test/test_frameconnect.jl
+
+# Interactive development
+julia --project
+using SMLMFrameConnection, SMLMData
+```
+
+## Architecture
+
+### Algorithm Pipeline
+
+The main entry point `frameconnect()` executes four stages:
+
+1. **Preclustering** (`precluster.jl`): Groups spatiotemporally nearby localizations using KD-tree nearest-neighbor search. Distance threshold: `nsigmadev × mean(σ_x, σ_y)`. Temporal threshold: `maxframegap` frames.
+
+2. **Parameter Estimation** (`estimateparams.jl`, `estimatedensities.jl`): Estimates photophysics (k_on, k_off, k_bleach, p_miss) and emitter densities from precluster statistics.
+
+3. **LAP Solving** (`create_costmatrix.jl`, `solveLAP.jl`, `linkclusters.jl`): Constructs cost matrices for each precluster and solves linear assignment problems via Hungarian algorithm to determine which localizations belong to the same physical emitter.
+
+4. **Combination** (`combinelocalizations.jl`): Merges connected localizations using MLE weighted mean: `x = Σ(x/σ²) / Σ(1/σ²)`, uncertainty: `σ = √(1/Σ(1/σ²))`.
+
+### Key Types
+
+- **Input**: `BasicSMLD{T, Emitter2DFit{T}}` from SMLMData.jl
+- **Output**: `(combined::BasicSMLD, info::ConnectInfo)` tuple
+- **ConnectInfo{T}** (`connectinfo.jl`): Secondary output with connected SMLD, statistics, and photophysics
+- **ParamStruct** (`structdefinitions.jl`): Internal algorithm parameters
+
+### File Organization
+
+```
+src/
+├── frameconnect.jl          # Main entry point
+├── connectinfo.jl           # ConnectInfo output struct
+├── structdefinitions.jl     # ParamStruct, helper conversions
+├── precluster.jl            # Stage 1: spatiotemporal clustering
+├── organizeclusters.jl      # Prepare cluster data structures
+├── estimateparams.jl        # Stage 2: photophysics estimation
+├── estimatedensities.jl     # Stage 2: density estimation
+├── create_costmatrix.jl     # Stage 3: LAP cost matrix construction
+├── solveLAP.jl              # Stage 3: Hungarian algorithm wrapper
+├── linkclusters.jl          # Stage 3: apply LAP solution
+├── combinelocalizations.jl  # Stage 4: MLE combination with covariance
+├── defineidealFC.jl         # Ground-truth validation utility
+└── compress_connectID.jl    # Normalize track_id to 1:N
+```
+
+### Dependencies
+
+- **SMLMData.jl** (v0.6+): Core types (BasicSMLD, Emitter2DFit with σ_xy)
+- **Hungarian.jl**: LAP solver
+- **NearestNeighbors.jl**: KD-tree for spatial queries
+- **Optim.jl**: Parameter optimization
+
+## Testing
+
+Tests use helper fixtures in `test/test_helpers.jl`:
+- `make_test_camera()`, `make_emitter()`, `make_test_smld()`
+- `make_blinking_molecule()`: Creates emitter sequence for one molecule
+- `make_two_molecules_smld()`, `make_single_molecule_smld()`
+
+## API
+
+Five exports:
+- `frameconnect(smld)` → `(combined::BasicSMLD, info::ConnectInfo)` tuple
+- `ConnectInfo{T}` → struct with connected SMLD, statistics (n_input, n_tracks, n_combined), photophysics (k_on, k_off, k_bleach, p_miss), timing
+- `combinelocalizations(smld)` → combines emitters by `track_id` using MLE with full covariance
+- `defineidealFC(smld)` → ground-truth frame connection for simulated data
+- `ParamStruct` → internal algorithm parameters
diff --git a/Project.toml b/Project.toml
index 28498b4..540df8e 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.2.0"
+version = "0.3.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.5"
+SMLMData = "0.6"
 StatsBase = "0.33, 0.34"
 julia = "1.6"
 
diff --git a/README.md b/README.md
index 2885d8d..bf5bb5a 100644
--- a/README.md
+++ b/README.md
@@ -24,9 +24,10 @@ Pkg.add("SMLMFrameConnection")
 using SMLMData, SMLMFrameConnection
 
 # Run frame connection on your data
-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 connected SMLD, statistics, and estimated photophysics
 ```
 
 ## Input Requirements
@@ -52,41 +53,53 @@ using SMLMData, SMLMFrameConnection
 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)
+# Constructor: Emitter2DFit{T}(x, y, photons, bg, σ_x, σ_y, σ_xy, σ_photons, σ_bg, frame, dataset, track_id, id)
 emitters = [
-    Emitter2DFit{Float64}(
-        5.0, 5.0,      # 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
-    ),
-    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.0, 5.0, 1000.0, 10.0, 0.02, 0.02, 0.0, 50.0, 2.0, 1, 1, 0, 1),
+    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)
+(combined, info) = frameconnect(smld)
 
 # Result: localizations connected based on spatial/temporal proximity
+println("Combined $(info.n_input) localizations into $(info.n_combined)")
+println("Formed $(info.n_tracks) tracks from $(info.n_preclusters) preclusters")
+
 # Combined uncertainty: σ_combined ≈ σ_individual / √n_connected
 ```
 
 ## Outputs Explained
 
 ```julia
-smld_connected, smld_preclustered, smld_combined, params = frameconnect(smld)
+(combined, info) = frameconnect(smld)
 ```
 
-| 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 |
+| Output | Type | Description |
+|--------|------|-------------|
+| `combined` | `BasicSMLD` | **Main output.** Combined high-precision localizations |
+| `info` | `ConnectInfo` | Metadata: connected SMLD, statistics, estimated photophysics |
+
+### ConnectInfo Fields
+
+| Field | Type | Description |
+|-------|------|-------------|
+| `info.connected` | `BasicSMLD` | Original localizations with `track_id` populated |
+| `info.n_input` | `Int` | Number of input localizations |
+| `info.n_tracks` | `Int` | Number of tracks formed |
+| `info.n_combined` | `Int` | Number of output localizations |
+| `info.n_preclusters` | `Int` | Number of preclusters |
+| `info.k_on` | `Float64` | Estimated on rate (1/frame) |
+| `info.k_off` | `Float64` | Estimated off rate (1/frame) |
+| `info.k_bleach` | `Float64` | Estimated bleach rate (1/frame) |
+| `info.p_miss` | `Float64` | Estimated miss probability |
+| `info.initialdensity` | `Vector{Float64}` | Density estimate per cluster (emitters/μm²) |
+| `info.elapsed_ns` | `UInt64` | Processing time in nanoseconds |
+| `info.algorithm` | `Symbol` | Algorithm used (`:lap`) |
 
 ## Parameters
 
@@ -104,36 +117,24 @@ frameconnect(smld;
 - `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`
+Connected localizations are combined using **maximum likelihood estimation (MLE) weighted mean** with full covariance propagation:
+- Position: inverse-variance weighted average using 2x2 covariance matrix
+- Uncertainty: properly propagated including σ_xy correlation
 - Photons: summed across connected localizations
 
 ## Utility Functions
 
 ### combinelocalizations
 ```julia
-smld_combined = combinelocalizations(smld)
+combined = combinelocalizations(smld)
 ```
 Combines emitters that share the same `track_id`. Use when you have pre-labeled data.
 
 ### defineidealFC
 ```julia
-smld_connected, smld_combined = defineidealFC(smld; maxframegap=5)
+(connected, combined) = defineidealFC(smld; maxframegap=5)
 ```
 For **simulated data** where `track_id` already contains ground-truth emitter IDs. Useful for validating frame-connection performance against known truth.
 
diff --git a/api_overview.md b/api_overview.md
index aa34ece..ce3e7ef 100644
--- a/api_overview.md
+++ b/api_overview.md
@@ -2,16 +2,20 @@
 
 Frame-connection for 2D SMLM data: combines repeated localizations of blinking fluorophores into higher-precision localizations.
 
+## Exports
+
+4 exports: `frameconnect`, `combinelocalizations`, `defineidealFC`, `ConnectInfo`, `ParamStruct`
+
 ## Main Functions
 
 ### frameconnect
 ```julia
-frameconnect(smld::BasicSMLD{T,Emitter2DFit{T}};
+(combined, info) = 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)
+) where T
 ```
 Main entry point. Connects repeated localizations and combines them.
 
@@ -21,28 +25,44 @@ Main entry point. Connects repeated localizations and combines them.
 - `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::BasicSMLD`: **Main output** - combined high-precision localizations
+- `info::ConnectInfo`: Metadata including connected SMLD, statistics, and estimated photophysics
 
 ### combinelocalizations
 ```julia
-combinelocalizations(smld::BasicSMLD{T,Emitter2DFit{T}}) where T -> BasicSMLD
+combined = combinelocalizations(smld::BasicSMLD{T,Emitter2DFit{T}}) where T -> BasicSMLD
 ```
-Combines emitters sharing the same `track_id` using MLE weighted mean. Use when `track_id` is already populated.
+Combines emitters sharing the same `track_id` using MLE weighted mean with full covariance propagation. Use when `track_id` is already populated.
 
 ### defineidealFC
 ```julia
-defineidealFC(smld::BasicSMLD{T,Emitter2DFit{T}};
+(connected, combined) = defineidealFC(smld::BasicSMLD{T,Emitter2DFit{T}};
     maxframegap::Int=5
-) where T -> (smld_connected, smld_combined)
+) where T
 ```
 For simulated data where `track_id` indicates ground-truth emitter ID. Useful for validation/benchmarking.
 
 ## Types
 
+### ConnectInfo{T}
+```julia
+struct ConnectInfo{T}
+    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                 # 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
+    initialdensity::Vector{Float64}  # Density estimate per cluster (emitters/μm²)
+    elapsed_ns::UInt64            # Processing time in nanoseconds
+    algorithm::Symbol             # Algorithm used (:lap)
+    n_preclusters::Int            # Number of preclusters formed
+end
+```
+
 ### ParamStruct
 ```julia
 mutable struct ParamStruct
@@ -76,20 +96,41 @@ Input `BasicSMLD` must contain `Emitter2DFit` emitters.
 **Optional fields:**
 - `photons`, `σ_photons`: Photon count (summed in output)
 - `bg`, `σ_bg`: Background
+- `σ_xy`: Position covariance (microns², propagated in output)
 - `dataset`: Dataset identifier (default: 1)
 - `track_id`: Set to 0 for input (populated by algorithm)
 
 ## Output
 
-Output `smld_combined` contains `Emitter2DFit` emitters with:
-- Combined position via MLE weighted mean: `x = Σ(x/σ²) / Σ(1/σ²)`
-- Reduced uncertainties: `σ = √(1/Σ(1/σ²))`
+Output `combined` contains `Emitter2DFit` emitters with:
+- Combined position via MLE weighted mean with full 2x2 covariance
+- Properly propagated uncertainties including σ_xy
 - Summed photons
 - `track_id` indicating connection group
 
+## Example
+
+```julia
+using SMLMData, SMLMFrameConnection
+
+# Load or create SMLD
+smld = ...
+
+# Run frame connection
+(combined, info) = frameconnect(smld; maxframegap=10)
+
+# Access results
+println("Combined $(info.n_input) → $(info.n_combined) localizations")
+println("Estimated k_on=$(info.k_on), k_off=$(info.k_off)")
+println("Processing time: $(info.elapsed_ns / 1e9) seconds")
+
+# Access connected (uncombined) localizations
+connected_smld = info.connected
+```
+
 ## Dependencies
 
-- SMLMData.jl (v0.5+): BasicSMLD, Emitter2DFit types
+- SMLMData.jl (v0.6+): BasicSMLD, Emitter2DFit types
 - Hungarian.jl (v0.6-0.7): 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..0a810c0 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -19,9 +19,10 @@ Pkg.add("SMLMFrameConnection")
 using SMLMData, 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 connected SMLD, statistics, and estimated photophysics
 ```
 
 ## Input Requirements
@@ -35,6 +36,7 @@ Input `BasicSMLD` must contain `Emitter2DFit` emitters with:
 
 **Optional:**
 - `photons`, `bg`: Photometry (summed in output)
+- `σ_xy`: Position covariance (propagated in output)
 - `dataset`: Dataset identifier (default: 1)
 
 ## Algorithm Overview
@@ -42,16 +44,32 @@ Input `BasicSMLD` must contain `Emitter2DFit` emitters with:
 1. **Preclustering**: Groups spatially and temporally adjacent localizations into candidate clusters
 2. **Parameter Estimation**: Estimates fluorophore blinking kinetics (`k_on`, `k_off`, `k_bleach`) and emitter density from the data
 3. **Frame Connection**: Uses Linear Assignment Problem (LAP) to optimally assign localizations to emitters based on spatial proximity and estimated photophysics
-4. **Combination**: Combines connected localizations using MLE weighted mean for improved precision
+4. **Combination**: Combines connected localizations using MLE weighted mean with full covariance propagation
 
 ## Outputs
 
-| 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 |
+```julia
+(combined, info) = frameconnect(smld)
+```
+
+| Output | Type | Description |
+|--------|------|-------------|
+| `combined` | `BasicSMLD` | **Main output** - combined high-precision localizations |
+| `info` | `ConnectInfo` | Metadata: connected SMLD, statistics, photophysics |
+
+### ConnectInfo Fields
+
+| Field | Description |
+|-------|-------------|
+| `info.connected` | Original localizations with `track_id` labels |
+| `info.n_input` | Number of input localizations |
+| `info.n_tracks` | Number of tracks formed |
+| `info.n_combined` | Number of output localizations |
+| `info.n_preclusters` | Number of preclusters |
+| `info.k_on`, `k_off`, `k_bleach` | Estimated photophysics |
+| `info.p_miss` | Estimated miss probability |
+| `info.elapsed_ns` | Processing time (nanoseconds) |
+| `info.algorithm` | Algorithm used (`:lap`) |
 
 ## Parameters
 
@@ -67,6 +85,23 @@ frameconnect(smld;
 - `nsigmadev`: Higher values allow connections over larger distances
 - `maxframegap`: Increase for dyes with long dark states (dSTORM: 10-20)
 
+## Example
+
+```julia
+using SMLMData, SMLMFrameConnection
+
+# Run frame connection
+(combined, info) = frameconnect(smld; maxframegap=10)
+
+# Access results
+println("Combined $(info.n_input) → $(info.n_combined) localizations")
+println("Formed $(info.n_tracks) tracks")
+println("Estimated k_on=$(info.k_on), k_off=$(info.k_off)")
+
+# Access connected (uncombined) localizations with track_id
+connected_smld = info.connected
+```
+
 ## API Reference
 
 ```@index
diff --git a/src/SMLMFrameConnection.jl b/src/SMLMFrameConnection.jl
index b2d34fd..1dc1af4 100644
--- a/src/SMLMFrameConnection.jl
+++ b/src/SMLMFrameConnection.jl
@@ -7,9 +7,10 @@ using Optim
 using StatsBase
 using Statistics
 
-export frameconnect, defineidealFC, combinelocalizations, ParamStruct
+export frameconnect, defineidealFC, combinelocalizations, ParamStruct, ConnectInfo
 
 include("structdefinitions.jl")
+include("connectinfo.jl")
 include("precluster.jl")
 include("defineidealFC.jl")
 include("organizeclusters.jl")
diff --git a/src/combinelocalizations.jl b/src/combinelocalizations.jl
index 2af1bd4..eff4258 100644
--- a/src/combinelocalizations.jl
+++ b/src/combinelocalizations.jl
@@ -21,6 +21,7 @@ function combinelocalizations(smld::BasicSMLD{T,E}) where {T, E<:SMLMData.Abstra
     y = [e.y for e in emitters][sortindices]
     σ_x = [e.σ_x for e in emitters][sortindices]
     σ_y = [e.σ_y for e in emitters][sortindices]
+    σ_xy = [e.σ_xy for e in emitters][sortindices]
     photons = [e.photons for e in emitters][sortindices]
     σ_photons = [e.σ_photons for e in emitters][sortindices]
     bg = [e.bg for e in emitters][sortindices]
@@ -42,6 +43,7 @@ function combinelocalizations(smld::BasicSMLD{T,E}) where {T, E<:SMLMData.Abstra
     y_combined = Vector{ET}(undef, nclusters)
     σ_x_combined = Vector{ET}(undef, nclusters)
     σ_y_combined = Vector{ET}(undef, nclusters)
+    σ_xy_combined = Vector{ET}(undef, nclusters)
     photons_combined = Vector{ET}(undef, nclusters)
     σ_photons_combined = Vector{ET}(undef, nclusters)
     bg_combined = Vector{ET}(undef, nclusters)
@@ -53,14 +55,46 @@ function combinelocalizations(smld::BasicSMLD{T,E}) where {T, E<:SMLMData.Abstra
 
     for nn in 1:nclusters
         indices = (1:nperID[nn]) .+ ncumulative[nn]
-        x_combined[nn] = StatsBase.mean(x[indices], weights(1.0 ./ σ_x[indices] .^ 2))
-        y_combined[nn] = StatsBase.mean(y[indices], weights(1.0 ./ σ_y[indices] .^ 2))
-        σ_x_combined[nn] = sqrt(1.0 / sum(1.0 ./ σ_x[indices] .^ 2))
-        σ_y_combined[nn] = sqrt(1.0 / sum(1.0 ./ σ_y[indices] .^ 2))
+
+        # Combine positions using precision-weighted mean with full 2x2 covariance.
+        # For each measurement i: Σᵢ = [σ_x² σ_xy; σ_xy σ_y²], Pᵢ = Σᵢ⁻¹
+        # Combined: Σ_comb = (Σ Pᵢ)⁻¹, pos_comb = Σ_comb * Σ(Pᵢ * posᵢ)
+        P_sum = zeros(ET, 2, 2)
+        μ_weighted = zeros(ET, 2)
+        for i in indices
+            det_i = σ_x[i]^2 * σ_y[i]^2 - σ_xy[i]^2
+            det_i = max(det_i, eps(ET))  # Ensure positive definite
+            # Precision matrix: P = [σ_y² -σ_xy; -σ_xy σ_x²] / det
+            P11 = σ_y[i]^2 / det_i
+            P22 = σ_x[i]^2 / det_i
+            P12 = -σ_xy[i] / det_i
+            P_sum[1,1] += P11
+            P_sum[2,2] += P22
+            P_sum[1,2] += P12
+            P_sum[2,1] += P12
+            μ_weighted[1] += P11 * x[i] + P12 * y[i]
+            μ_weighted[2] += P12 * x[i] + P22 * y[i]
+        end
+
+        # Combined covariance = inverse of summed precision
+        det_P = P_sum[1,1] * P_sum[2,2] - P_sum[1,2]^2
+        Σ_comb_11 = P_sum[2,2] / det_P  # σ_x²
+        Σ_comb_22 = P_sum[1,1] / det_P  # σ_y²
+        Σ_comb_12 = -P_sum[1,2] / det_P  # σ_xy
+
+        # Combined position
+        x_combined[nn] = Σ_comb_11 * μ_weighted[1] + Σ_comb_12 * μ_weighted[2]
+        y_combined[nn] = Σ_comb_12 * μ_weighted[1] + Σ_comb_22 * μ_weighted[2]
+        σ_x_combined[nn] = sqrt(Σ_comb_11)
+        σ_y_combined[nn] = sqrt(Σ_comb_22)
+        σ_xy_combined[nn] = Σ_comb_12
+
+        # Photons and background: sum (independent measurements)
         photons_combined[nn] = sum(photons[indices])
         σ_photons_combined[nn] = sqrt(sum(σ_photons[indices] .^ 2))
         bg_combined[nn] = sum(bg[indices])
         σ_bg_combined[nn] = sqrt(sum(σ_bg[indices] .^ 2))
+
         connectID_combined[nn] = connectID[indices[1]]
         framenum_combined[nn] = framenum[indices[1]]
         datasetnum_combined[nn] = datasetnum[indices[1]]
@@ -74,6 +108,7 @@ function combinelocalizations(smld::BasicSMLD{T,E}) where {T, E<:SMLMData.Abstra
             x_combined[nn], y_combined[nn],
             photons_combined[nn], bg_combined[nn],
             σ_x_combined[nn], σ_y_combined[nn],
+            σ_xy_combined[nn],
             σ_photons_combined[nn], σ_bg_combined[nn],
             framenum_combined[nn], datasetnum_combined[nn],
             connectID_combined[nn], id_combined[nn]
diff --git a/src/computeclusterinfo.jl b/src/computeclusterinfo.jl
index d35ab4a..ce9b789 100644
--- a/src/computeclusterinfo.jl
+++ b/src/computeclusterinfo.jl
@@ -15,7 +15,8 @@ function computeclusterinfo(clusterdata::Vector{Matrix{Float32}})
     nobservations = ones(nclusters)
     for nn = 1:nclusters
         # Compute the duration of this cluster.
-        currentframes = clusterdata[nn][:, 5]
+        # Column layout: 1:x 2:y 3:σ_x 4:σ_y 5:σ_xy 6:frame 7:dataset 8:connectID 9:sortindex
+        currentframes = clusterdata[nn][:, 6]
         clusterdurations[nn] = maximum(currentframes) - minimum(currentframes) + 1
 
         # Compute the number of observations of this cluster excluding the 
diff --git a/src/connectinfo.jl b/src/connectinfo.jl
new file mode 100644
index 0000000..9201033
--- /dev/null
+++ b/src/connectinfo.jl
@@ -0,0 +1,33 @@
+"""
+    ConnectInfo{T}
+
+Secondary output from `frameconnect()` containing connected localizations and 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 (combined) 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
+- `initialdensity::Vector{Float64}`: Initial density estimate per precluster (emitters/μm²)
+- `elapsed_ns::UInt64`: Wall time in nanoseconds
+- `algorithm::Symbol`: Algorithm used (`:lap`)
+- `n_preclusters::Int`: Number of preclusters formed
+"""
+struct ConnectInfo{T}
+    connected::BasicSMLD{T}
+    n_input::Int
+    n_tracks::Int
+    n_combined::Int
+    k_on::Float64
+    k_off::Float64
+    k_bleach::Float64
+    p_miss::Float64
+    initialdensity::Vector{Float64}
+    elapsed_ns::UInt64
+    algorithm::Symbol
+    n_preclusters::Int
+end
diff --git a/src/connectlocalizations.jl b/src/connectlocalizations.jl
index 176aadb..df9bc8b 100644
--- a/src/connectlocalizations.jl
+++ b/src/connectlocalizations.jl
@@ -28,7 +28,7 @@ function connectlocalizations!(connectID::Vector{Int64},
 
         # Update `connectID` to reflect the LAP solution.
         _, maxconnectID = linkclusters!(connectID, maxconnectID,
-            Int64.(clusterdata[nn][:, 8]), assignment[1])
+            Int64.(clusterdata[nn][:, 9]), assignment[1])  # sortindex is column 9
     end
 
     # Ensure `connectID` is contains the set of integers 1:nclusters.
diff --git a/src/create_costmatrix.jl b/src/create_costmatrix.jl
index f9215d2..656aa38 100644
--- a/src/create_costmatrix.jl
+++ b/src/create_costmatrix.jl
@@ -20,11 +20,13 @@ function create_costmatrix(clusterdata::Vector{Matrix{Float32}},
     params::ParamStruct, clusterind::Int64, nframes::Int64)
     # Extract some entries from clusterdata and params to make the code more
     # readable.
+    # Column layout: 1:x 2:y 3:σ_x 4:σ_y 5:σ_xy 6:frame 7:dataset 8:connectID 9:sortindex
     x = clusterdata[clusterind][:, 1]
     y = clusterdata[clusterind][:, 2]
     x_se = clusterdata[clusterind][:, 3]
     y_se = clusterdata[clusterind][:, 4]
-    framenum = clusterdata[clusterind][:, 5]
+    xy_cov = clusterdata[clusterind][:, 5]
+    framenum = clusterdata[clusterind][:, 6]
     k_on = params.k_on
     k_off = params.k_off
     k_bleach = params.k_bleach
@@ -54,11 +56,19 @@ function create_costmatrix(clusterdata::Vector{Matrix{Float32}},
             #       needed.
             continue
         end
-        sigma_x = sqrt(x_se[mm]^2 + x_se[nn]^2)
-        sigma_y = sqrt(y_se[mm]^2 + y_se[nn]^2)
-        separationcost = (log(2*pi*sigma_x*sigma_y)
-            + (x[mm]-x[nn])^2 / (2*sigma_x^2)
-            + (y[mm]-y[nn])^2 / (2*sigma_y^2));
+        # Combined covariance matrix for two localizations:
+        # Σ_comb = Σ_mm + Σ_nn (variances add for independent measurements)
+        σ_x² = x_se[mm]^2 + x_se[nn]^2
+        σ_y² = y_se[mm]^2 + y_se[nn]^2
+        σ_xy = xy_cov[mm] + xy_cov[nn]
+        det_Σ = σ_x² * σ_y² - σ_xy^2
+
+        # Mahalanobis distance: Δᵀ Σ⁻¹ Δ where Σ⁻¹ = [σ_y² -σ_xy; -σ_xy σ_x²] / det
+        Δx, Δy = x[mm] - x[nn], y[mm] - y[nn]
+        mahal = (σ_y² * Δx^2 - 2*σ_xy*Δx*Δy + σ_x² * Δy^2) / det_Σ
+
+        # Negative log-likelihood: -log(p) = 0.5*log(2π) + 0.5*log(det) + 0.5*mahal
+        separationcost = log(2*pi) + 0.5*log(det_Σ) + 0.5*mahal
         observationcost = -log((p_miss^(deltaframe-1)) * (1-p_miss))
         stilloncost = (k_off+k_bleach) * deltaframe
         costmatrix[mm, nn] = (separationcost+observationcost+stilloncost) / 2.0
@@ -81,8 +91,10 @@ function create_costmatrix(clusterdata::Vector{Matrix{Float32}},
     for nn in indices
         deltaframe_past = minimum([maxframegap; framesint[nn]-startframe])
         deltaframe_future = minimum([maxframegap; nframes-framesint[nn]])
-        # Localization uncertainty area (μm²) to make density dimensionless
-        loc_area = pi * x_se[nn] * y_se[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
+        loc_area = pi * sqrt(max(det_loc, eps()))
         birthcost = -log(1-p_miss) -
             log(rho_off[framesint[nn]] * loc_area * (1-exp(-k_on)) * exp(-deltaframe_past*k_on) +
                 rho_on[framesint[nn]-deltaframe_past] * loc_area * (p_miss^deltaframe_past))
diff --git a/src/defineidealFC.jl b/src/defineidealFC.jl
index 6a6f8b7..c41271d 100644
--- a/src/defineidealFC.jl
+++ b/src/defineidealFC.jl
@@ -136,7 +136,7 @@ function defineidealFC(smld::BasicSMLD{T,E};
     for i in 1:length(emitters)
         e = emitters[i]
         new_emitters[i] = SMLMData.Emitter2DFit{ET}(
-            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, connectID[i], e.id
         )
     end
diff --git a/src/estimatedensities.jl b/src/estimatedensities.jl
index ba4cb29..72250a7 100644
--- a/src/estimatedensities.jl
+++ b/src/estimatedensities.jl
@@ -43,18 +43,38 @@ function estimatedensities(smld::BasicSMLD{T,E},
     end
 
     # Determine the center of all clusters assuming each arose from the same
-    # emitter.
+    # emitter using precision-weighted mean with full covariance.
     clustercenters = zeros(2, nclusters)
     for nn = 1:nclusters
         # Isolate some arrays to improve readability.
+        # Column layout: 1:x 2:y 3:σ_x 4:σ_y 5:σ_xy 6:frame ...
         x = clusterdata[nn][:, 1]
         y = clusterdata[nn][:, 2]
         σ_x = clusterdata[nn][:, 3]
         σ_y = clusterdata[nn][:, 4]
+        σ_xy = clusterdata[nn][:, 5]
 
-        # Compute the cluster centers based on MLE of position.
-        clustercenters[1:2, nn] = [sum(x ./ σ_x .^ 2) / sum(1 ./ σ_x .^ 2)
-            sum(y ./ σ_y .^ 2) / sum(1 ./ σ_y .^ 2)]
+        # Compute the cluster centers using precision-weighted mean with full covariance.
+        # P = Σ⁻¹ = [σ_y² -σ_xy; -σ_xy σ_x²] / det where det = σ_x²σ_y² - σ_xy²
+        # Combined: center = (Σ Pᵢ)⁻¹ * Σ(Pᵢ * posᵢ)
+        P_sum = zeros(2, 2)
+        μ_weighted = zeros(2)
+        for i in eachindex(x)
+            det_i = σ_x[i]^2 * σ_y[i]^2 - σ_xy[i]^2
+            det_i = max(det_i, eps())  # Ensure positive definite
+            # Precision matrix elements
+            P11 = σ_y[i]^2 / det_i
+            P22 = σ_x[i]^2 / det_i
+            P12 = -σ_xy[i] / det_i
+            P_sum[1,1] += P11
+            P_sum[2,2] += P22
+            P_sum[1,2] += P12
+            P_sum[2,1] += P12
+            μ_weighted[1] += P11 * x[i] + P12 * y[i]
+            μ_weighted[2] += P12 * x[i] + P22 * y[i]
+        end
+        # Solve P_sum * center = μ_weighted
+        clustercenters[1:2, nn] = P_sum \ μ_weighted
     end
 
     # Estimate the local cluster density based on the distance to
diff --git a/src/frameconnect.jl b/src/frameconnect.jl
index 7f60e45..346595f 100644
--- a/src/frameconnect.jl
+++ b/src/frameconnect.jl
@@ -1,7 +1,7 @@
 using SMLMData
 
 """
-    result = frameconnect(smld::BasicSMLD{T,E};
+    (combined, info) = frameconnect(smld::BasicSMLD{T,E};
         nnearestclusters::Int=2, nsigmadev::Float64=5.0,
         maxframegap::Int=5, nmaxnn::Int=2) where {T, E<:SMLMData.AbstractEmitter}
 
@@ -29,20 +29,36 @@ using their MLE position estimate assuming Gaussian noise.
             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).
+Returns a tuple `(combined, info)`:
+- `combined::BasicSMLD`: Final frame-connection result with localizations combined
+                         into higher precision localizations.
+- `info::ConnectInfo`: Metadata including connected SMLD (with track_id assigned),
+                       estimated photophysics parameters, timing, and statistics.
+
+# Example
+```julia
+using SMLMFrameConnection, SMLMData
+
+# Load or create SMLD with localizations
+smld = ...
+
+# Run frame connection
+(combined, info) = frameconnect(smld)
+
+# Access results
+println("Combined \$(info.n_input) localizations into \$(info.n_combined)")
+println("Formed \$(info.n_tracks) tracks from \$(info.n_preclusters) preclusters")
+
+# Access connected (uncombined) localizations with track_id
+connected_smld = info.connected
+```
 """
 function frameconnect(smld::BasicSMLD{T,E};
     nnearestclusters::Int = 2, nsigmadev::Float64 = 5.0,
     maxframegap::Int = 5, nmaxnn::Int = 2) where {T, E<:SMLMData.AbstractEmitter}
 
+    t_start = time_ns()
+
     # Prepare a ParamStruct to keep track of parameters used.
     params = ParamStruct()
     params.nnearestclusters = nnearestclusters
@@ -54,6 +70,9 @@ function frameconnect(smld::BasicSMLD{T,E};
     smld_preclustered = precluster(smld, params)
     clusterdata = organizeclusters(smld_preclustered)
 
+    # Count preclusters
+    n_preclusters = length(unique(e.track_id for e in smld_preclustered.emitters))
+
     # Estimate rate parameters.
     params.k_on, params.k_off, params.k_bleach, params.p_miss =
         estimateparams(smld_preclustered, clusterdata)
@@ -79,7 +98,7 @@ function frameconnect(smld::BasicSMLD{T,E};
     for i in 1:length(emitters)
         e = emitters[i]
         new_emitters[i] = SMLMData.Emitter2DFit{ET}(
-            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, connectID_final[i], e.id
         )
     end
@@ -89,5 +108,25 @@ 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_ns = time_ns() - t_start
+
+    # Count tracks and build info
+    n_tracks = length(unique(e.track_id for e in smld_connected.emitters))
+
+    info = ConnectInfo{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.initialdensity,
+        elapsed_ns,
+        :lap,
+        n_preclusters
+    )
+
+    return (smld_combined, info)
 end
diff --git a/src/organizeclusters.jl b/src/organizeclusters.jl
index dab4426..5e3c2cc 100644
--- a/src/organizeclusters.jl
+++ b/src/organizeclusters.jl
@@ -18,12 +18,14 @@ function organizeclusters(smld::BasicSMLD{T,E}) where {T, E<:SMLMData.AbstractEm
     y = [e.y for e in emitters]
     σ_x = [e.σ_x for e in emitters]
     σ_y = [e.σ_y for e in emitters]
+    σ_xy = [e.σ_xy for e in emitters]
     framenum = [e.frame for e in emitters]
     datasetnum = [e.dataset for e in emitters]
 
     # Sort by connectID
+    # Column layout: 1:x 2:y 3:σ_x 4:σ_y 5:σ_xy 6:frame 7:dataset 8:connectID
     sortindices = sortperm(connectID)
-    combineddata = [x y σ_x σ_y framenum datasetnum connectID]
+    combineddata = [x y σ_x σ_y σ_xy framenum datasetnum connectID]
     combineddata = combineddata[sortindices, :]
 
     # Loop over preclusters and organize member localizations into a more
diff --git a/src/precluster.jl b/src/precluster.jl
index b3bf6dd..4f2e5b7 100644
--- a/src/precluster.jl
+++ b/src/precluster.jl
@@ -130,7 +130,7 @@ function precluster(smld::BasicSMLD{T,E},
     for i in 1:n_emitters
         e = emitters[i]
         new_emitters[i] = SMLMData.Emitter2DFit{ET}(
-            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, connectID_compressed[i], e.id
         )
     end
diff --git a/src/structdefinitions.jl b/src/structdefinitions.jl
index 771503b..a652b3d 100644
--- a/src/structdefinitions.jl
+++ b/src/structdefinitions.jl
@@ -51,7 +51,7 @@ consistent Emitter2DFit output.
 function to_emitter2dfit(e::SMLMData.AbstractEmitter, track_id::Int)
     T = typeof(e.x)
     SMLMData.Emitter2DFit{T}(
-        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, track_id, e.id
     )
 end
\ No newline at end of file
diff --git a/test/test_combinelocalizations.jl b/test/test_combinelocalizations.jl
index 67f8cf3..85466f5 100644
--- a/test/test_combinelocalizations.jl
+++ b/test/test_combinelocalizations.jl
@@ -20,8 +20,8 @@
         # Two localizations at same position with same uncertainty
         σ = 0.02
         emitters = [
-            make_emitter(5.0, 5.0, 1; σ_xy=σ, track_id=1, photons=1000.0),
-            make_emitter(5.0, 5.0, 2; σ_xy=σ, track_id=1, photons=1200.0),
+            make_emitter(5.0, 5.0, 1; σ_pos=σ, track_id=1, photons=1000.0),
+            make_emitter(5.0, 5.0, 2; σ_pos=σ, track_id=1, photons=1200.0),
         ]
         smld = make_test_smld(emitters)
         result = combinelocalizations(smld)
@@ -49,8 +49,8 @@
         # Two localizations with different uncertainties
         # Weighted mean should favor lower uncertainty
         emitters = [
-            make_emitter(5.0, 5.0, 1; σ_xy=0.01, track_id=1),  # High precision
-            make_emitter(5.1, 5.1, 2; σ_xy=0.10, track_id=1),  # Low precision
+            make_emitter(5.0, 5.0, 1; σ_pos=0.01, track_id=1),  # High precision
+            make_emitter(5.1, 5.1, 2; σ_pos=0.10, track_id=1),  # Low precision
         ]
         smld = make_test_smld(emitters)
         result = combinelocalizations(smld)
@@ -68,9 +68,9 @@
     @testset "three localizations" begin
         σ = 0.02
         emitters = [
-            make_emitter(5.00, 5.00, 1; σ_xy=σ, track_id=1),
-            make_emitter(5.02, 5.02, 2; σ_xy=σ, track_id=1),
-            make_emitter(4.98, 4.98, 3; σ_xy=σ, track_id=1),
+            make_emitter(5.00, 5.00, 1; σ_pos=σ, track_id=1),
+            make_emitter(5.02, 5.02, 2; σ_pos=σ, track_id=1),
+            make_emitter(4.98, 4.98, 3; σ_pos=σ, track_id=1),
         ]
         smld = make_test_smld(emitters)
         result = combinelocalizations(smld)
diff --git a/test/test_frameconnect.jl b/test/test_frameconnect.jl
index 4f91ff3..1c05e45 100644
--- a/test/test_frameconnect.jl
+++ b/test/test_frameconnect.jl
@@ -5,15 +5,15 @@
             [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 ConnectInfo
+        @test info.connected isa BasicSMLD
+        @test info.algorithm == :lap
+        @test info.elapsed_ns > 0
+        @test info.n_input == length(smld.emitters)
+        @test info.n_combined == length(combined.emitters)
     end
 
     @testset "single molecule connection" begin
@@ -25,11 +25,11 @@
         )
         smld = make_test_smld(emitters; n_frames=5)
 
-        result = frameconnect(smld; maxframegap=5, nnearestclusters=1)
+        (combined, info) = frameconnect(smld; maxframegap=5, nnearestclusters=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)
     end
 
     @testset "multiple separated molecules" begin
@@ -42,11 +42,11 @@
         )
         smld = make_test_smld(emitters; n_frames=3)
 
-        result = frameconnect(smld; maxframegap=5)
+        (combined, info) = frameconnect(smld; maxframegap=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,10 +57,10 @@
         )
         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
+        for e in info.connected.emitters
             @test e.track_id > 0
         end
     end
@@ -74,14 +74,16 @@
         )
         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.initialdensity)
+        @test info.n_preclusters > 0
+        @test info.n_tracks > 0
     end
 
     @testset "respects maxframegap" begin
@@ -95,10 +97,10 @@
         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)
+        (combined, info) = frameconnect(smld; maxframegap=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 +108,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
@@ -128,7 +132,7 @@
         σ_input = 0.02
 
         # Target molecule to test precision improvement
-        target_emitters = [make_emitter(5.0, 5.0, i; σ_xy=σ_input) for i in 1:n_locs]
+        target_emitters = [make_emitter(5.0, 5.0, i; σ_pos=σ_input) for i in 1:n_locs]
         # Background molecules for proper density estimation
         background = vcat(
             [make_emitter(Float64(20+i), Float64(20+i), 1) for i in 1:4]...
@@ -136,10 +140,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; maxframegap=n_locs, nnearestclusters=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 +153,20 @@
         @test σ_output ≈ expected_σ rtol=0.1
     end
 
-    @testset "custom parameters" begin
-        smld = make_single_molecule_smld()
-
-        result = frameconnect(smld;
-            nnearestclusters=3,
-            nsigmadev=4.0,
-            maxframegap=10,
-            nmaxnn=3
+    @testset "ConnectInfo statistics" begin
+        emitters = vcat(
+            make_blinking_molecule(5.0, 5.0, [1, 2, 3]),
+            make_blinking_molecule(10.0, 10.0, [1, 2]),
         )
+        smld = make_test_smld(emitters; n_frames=3)
+
+        (combined, info) = frameconnect(smld)
 
-        @test result.params.nnearestclusters == 3
-        @test result.params.nsigmadev == 4.0
-        @test result.params.maxframegap == 10
-        @test result.params.nmaxnn == 3
+        @test info.n_input == 5
+        @test info.n_combined <= info.n_input
+        @test info.n_tracks <= info.n_input
+        @test info.n_preclusters <= info.n_input
+        @test info.elapsed_ns isa UInt64
+        @test info.algorithm == :lap
     end
 end
diff --git a/test/test_helpers.jl b/test/test_helpers.jl
index ce96a9a..6548969 100644
--- a/test/test_helpers.jl
+++ b/test/test_helpers.jl
@@ -12,12 +12,14 @@ function make_test_camera(; xpixels=512, ypixels=512, pixelsize=0.1)
 end
 
 """
-    make_emitter(x, y, frame; σ_xy=0.02, photons=1000.0, track_id=0)
+    make_emitter(x, y, frame; σ_pos=0.02, σ_xy_cov=0.0, photons=1000.0, track_id=0)
 
 Create a single Emitter2DFit for testing with sensible defaults.
+σ_pos is used for both σ_x and σ_y, σ_xy_cov is the x-y covariance.
 """
 function make_emitter(x::Real, y::Real, frame::Int;
-                      σ_xy::Real=0.02,
+                      σ_pos::Real=0.02,  # position uncertainty (used for both σ_x and σ_y)
+                      σ_xy_cov::Real=0.0,  # x-y covariance
                       photons::Real=1000.0,
                       bg::Real=10.0,
                       track_id::Int=0,
@@ -25,7 +27,8 @@ function make_emitter(x::Real, y::Real, frame::Int;
     return Emitter2DFit{Float64}(
         Float64(x), Float64(y),
         Float64(photons), Float64(bg),
-        Float64(σ_xy), Float64(σ_xy),
+        Float64(σ_pos), Float64(σ_pos),
+        Float64(σ_xy_cov),  # σ_xy covariance
         Float64(sqrt(photons)), Float64(sqrt(bg)),
         frame, dataset, track_id, 0
     )
@@ -45,20 +48,20 @@ function make_test_smld(emitters::Vector{<:Emitter2DFit};
 end
 
 """
-    make_blinking_molecule(x, y, frames; σ_xy=0.02, position_jitter=0.002)
+    make_blinking_molecule(x, y, frames; σ_pos=0.02, position_jitter=0.002)
 
 Create emitters representing a single molecule blinking across multiple frames.
 Uses deterministic small offsets to simulate localization uncertainty.
 """
 function make_blinking_molecule(x::Real, y::Real, frames::AbstractVector{Int};
-                                σ_xy::Real=0.02,
+                                σ_pos::Real=0.02,
                                 position_jitter::Real=0.002,
                                 track_id::Int=0)
     emitters = Emitter2DFit{Float64}[]
     for (i, f) in enumerate(frames)
         # Deterministic small offset pattern (not random for reproducibility)
         offset = position_jitter * ((i % 3) - 1)  # -jitter, 0, +jitter pattern
-        push!(emitters, make_emitter(x + offset, y + offset, f; σ_xy=σ_xy, track_id=track_id))
+        push!(emitters, make_emitter(x + offset, y + offset, f; σ_pos=σ_pos, track_id=track_id))
     end
     return emitters
 end
diff --git a/test/test_types.jl b/test/test_types.jl
index 178fded..e398cc2 100644
--- a/test/test_types.jl
+++ b/test/test_types.jl
@@ -1,4 +1,42 @@
 @testset "Types" begin
+    @testset "ConnectInfo" begin
+        # Create minimal test data
+        camera = IdealCamera(512, 512, 0.1)
+        emitter = SMLMData.Emitter2DFit{Float64}(
+            1.0, 2.0, 1000.0, 50.0, 0.02, 0.02, 0.0, 30.0, 5.0, 1, 1, 1, 1
+        )
+        smld = BasicSMLD([emitter], camera, 1, 1)
+
+        info = ConnectInfo{Float64}(
+            smld,
+            10,      # n_input
+            5,       # n_tracks
+            5,       # n_combined
+            0.1,     # k_on
+            0.5,     # k_off
+            0.01,    # k_bleach
+            0.05,    # p_miss
+            [1.0],   # initialdensity
+            UInt64(1000000),  # elapsed_ns
+            :lap,    # algorithm
+            3        # n_preclusters
+        )
+
+        @test info isa ConnectInfo{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.initialdensity == [1.0]
+        @test info.elapsed_ns == UInt64(1000000)
+        @test info.algorithm == :lap
+        @test info.n_preclusters == 3
+    end
+
     @testset "ParamStruct" begin
         @testset "default constructor" begin
             params = ParamStruct()