PSF Width Fitting Example

This example demonstrates fitting with variable PSF width using the GaussianXYNBS model, which includes the PSF width as a fitted parameter.

When to Use Variable PSF Fitting

Use GaussianXYNBS when:

• PSF width varies across your dataset (e.g., z-dependent defocus)
• PSF width is unknown and needs to be measured
• Quality control: filter localizations by fitted PSF width
• Detecting PSF changes during acquisition

Basic Usage

using GaussMLE
using Statistics

# Create variable-sigma model (no fixed parameters)
psf = GaussianXYNBS()

# Create fitter
fitter = GaussMLEConfig(psf_model = psf)

# Fit data
data = rand(Float32, 11, 11, 100)
smld, info = fit(data, fitter)

# Access fitted PSF width from emitters (Emitter2DFitSigma type)
sigmas = [e.σ for e in smld.emitters]
sigma_uncertainties = [e.σ_σ for e in smld.emitters]

println("Mean PSF width: $(mean(sigmas)) microns")
println("Mean sigma uncertainty: $(mean(sigma_uncertainties)) microns")

Complete Working Example

using GaussMLE
using Statistics

println("=== PSF Width Fitting Example ===\n")

# Generate synthetic data
n_rois = 100
data = rand(Float32, 11, 11, n_rois)

# Create variable-sigma fitter
fitter = GaussMLEConfig(
    psf_model = GaussianXYNBS(),
    iterations = 25  # More iterations for 5-parameter fit
)

# Fit the data
println("Fitting $n_rois ROIs with variable PSF width...")
smld, info = fit(data, fitter)

# The emitters are Emitter2DFitSigma type with sigma field
println("\n=== Results ===")
println("Fitted: $(length(smld.emitters)) localizations")

# Extract all fields
x_positions = [e.x for e in smld.emitters]
y_positions = [e.y for e in smld.emitters]
photons = [e.photons for e in smld.emitters]
backgrounds = [e.bg for e in smld.emitters]
sigmas = [e.σ for e in smld.emitters]

# Uncertainties
sigma_x = [e.σ_x for e in smld.emitters]
sigma_sigma = [e.σ_σ for e in smld.emitters]

println("\nPosition Statistics:")
println("  Mean x: $(round(mean(x_positions), digits=3)) microns")
println("  Mean y: $(round(mean(y_positions), digits=3)) microns")
println("  Mean precision: $(round(mean(sigma_x)*1000, digits=1)) nm")

println("\nPhotometry Statistics:")
println("  Mean photons: $(round(mean(photons), digits=1))")
println("  Mean background: $(round(mean(backgrounds), digits=1))")

println("\nPSF Width Statistics:")
println("  Mean sigma: $(round(mean(sigmas), digits=4)) microns")
println("  Sigma std: $(round(std(sigmas), digits=4)) microns")
println("  Sigma range: $(round.(extrema(sigmas), digits=4)) microns")
println("  Mean sigma uncertainty: $(round(mean(sigma_sigma), digits=4)) microns")

Understanding Emitter2DFitSigma

The GaussianXYNBS model returns Emitter2DFitSigma emitters with these fields:

Quality Control with PSF Width

One key use of variable PSF fitting is quality control:

using GaussMLE
using Statistics

# Fit with variable sigma
fitter = GaussMLEConfig(psf_model = GaussianXYNBS())
smld, info = fit(data, fitter)

# Extract PSF widths
sigmas = [e.σ for e in smld.emitters]

# Analyze distribution
println("PSF Width Analysis:")
println("  Mean: $(round(mean(sigmas), digits=4)) microns")
println("  Median: $(round(median(sigmas), digits=4)) microns")
println("  Std: $(round(std(sigmas), digits=4)) microns")
println("  IQR: $(round.(quantile(sigmas, [0.25, 0.75]), digits=4)) microns")

# Filter by PSF width using @filter (typical range: 80-180nm)
valid = @filter(smld, σ > 0.08 && σ < 0.18)
println("\nValid localizations (sigma in expected range):")
println("  $(length(valid.emitters)) / $(length(smld.emitters))")

Comparing Fixed vs Variable PSF Models

using GaussMLE
using Statistics

# Same data, two models
data = rand(Float32, 11, 11, 1000)

# Fixed PSF model
fitter_fixed = GaussMLEConfig(psf_model = GaussianXYNB(0.13f0))
smld_fixed, info_fixed = fit(data, fitter_fixed)

# Variable PSF model
fitter_var = GaussMLEConfig(psf_model = GaussianXYNBS())
smld_var, info_var = fit(data, fitter_var)

# Compare position estimates
x_fixed = [e.x for e in smld_fixed.emitters]
x_var = [e.x for e in smld_var.emitters]

rms_diff = sqrt(mean((x_fixed .- x_var).^2))
println("Position difference (RMS): $(round(rms_diff*1000, digits=2)) nm")

# Compare uncertainties
sigma_x_fixed = mean([e.σ_x for e in smld_fixed.emitters])
sigma_x_var = mean([e.σ_x for e in smld_var.emitters])

println("Mean x uncertainty:")
println("  Fixed PSF:    $(round(sigma_x_fixed*1000, digits=2)) nm")
println("  Variable PSF: $(round(sigma_x_var*1000, digits=2)) nm")
println("  Ratio: $(round(sigma_x_var/sigma_x_fixed, digits=2))x")

# Performance comparison
t_fixed = @elapsed fit(data, fitter_fixed)
t_var = @elapsed fit(data, fitter_var)

println("\nPerformance:")
println("  Fixed PSF:    $(round(1000/t_fixed)) fits/second")
println("  Variable PSF: $(round(1000/t_var)) fits/second")
println("  Speed ratio: $(round(t_fixed/t_var, digits=2))x")

Using with ROIBatch

In a typical SMLM pipeline, ROIBatch comes from SMLMBoxer.jl. For testing:

using GaussMLE

# Create camera
camera = IdealCamera(0:1023, 0:1023, 0.1)  # 100nm pixels

# Generate test data
batch = generate_roi_batch(
    camera,
    GaussianXYNBS(),  # Variable sigma model
    n_rois = 100,
    roi_size = 11
)

# Fit
fitter = GaussMLEConfig(psf_model = GaussianXYNBS())
smld, info = fit(batch, fitter)

# Extract results - positions in camera coordinates
sigmas = [e.σ for e in smld.emitters]
println("Fitted PSF widths: $(round.(extrema(sigmas), digits=4)) microns")

Anisotropic PSF (GaussianXYNBSXSY)

For elliptical PSFs, use GaussianXYNBSXSY:

using GaussMLE
using Statistics

# Anisotropic model - fits sigma_x and sigma_y independently
fitter = GaussMLEConfig(psf_model = GaussianXYNBSXSY())
smld, info = fit(data, fitter)

# Returns Emitter2DFitSigmaXY emitters
# Fitted PSF widths (distinct from position uncertainties):
sigma_x_psf = [e.σx for e in smld.emitters]   # Fitted PSF width in x (microns)
sigma_y_psf = [e.σy for e in smld.emitters]   # Fitted PSF width in y (microns)

# Position uncertainties (CRLB):
precision_x = [e.σ_x for e in smld.emitters]  # x position uncertainty (microns)
precision_y = [e.σ_y for e in smld.emitters]   # y position uncertainty (microns)

# PSF width uncertainties:
sigma_sigma_x = [e.σ_σx for e in smld.emitters]  # Uncertainty in fitted σx
sigma_sigma_y = [e.σ_σy for e in smld.emitters]  # Uncertainty in fitted σy

Troubleshooting

PSF Width Convergence Issues

If fitted PSF widths are unreasonable:

# Check for extreme values
sigmas = [e.σ for e in smld.emitters]

extreme_low = count(s -> s < 0.05, sigmas)   # < 50nm
extreme_high = count(s -> s > 0.5, sigmas)   # > 500nm

println("Extreme low sigma: $extreme_low")
println("Extreme high sigma: $extreme_high")

if extreme_low + extreme_high > length(sigmas) * 0.1
    println("Warning: >10% of fits have extreme PSF widths")
    println("Consider:")
    println("  - Using fixed PSF model (GaussianXYNB)")
    println("  - Checking data quality")
    println("  - Increasing iterations")
end

Slow Convergence

# Use more iterations for variable PSF
fitter = GaussMLEConfig(
    psf_model = GaussianXYNBS(),
    iterations = 30  # Default is 20
)

Next Steps

• Learn about 3D astigmatic fitting with AstigmaticXYZNB
• Explore GPU Support for large datasets
• See the API Reference for all options