Generalized PSF Fitting Methods

• Generalized Point Spread Function Fitting (GPSF) is a framework that employs flexible PSF representations beyond single Gaussian or shift-invariant models.
• It integrates forward modeling with regularized inversion to estimate both object and PSF parameters across applications like microscopy, coherent imaging, and CMB analysis.
• Approaches range from sparse analog calibrations to neural-based phase fitting, offering enhanced accuracy and practical trade-offs in image reconstruction.

Generalized Point Spread Function Fitting (GPSF) denotes a family of inverse formulations in which the point spread function is fitted, calibrated, or embedded in a forward model under assumptions that are more general than a single Gaussian kernel, strict shift invariance, or a fixed low-dimensional basis. In the literature, the label has been used for sparse analog PSF calibration in fluorescence microscopy, spatially variant deconvolution in coherent imaging, semi-blind segmented PSF recovery from partially occulted images, basis-free neural phase-function fitting, PSF-inclusive beam-profile fitting for scintillating screens, and full-covariance multi-source template fitting in CMB polarization analysis (Samuylov et al., 2018, Lee et al., 2017, Hofmeister et al., 2022, Valouev, 2024, Novokshonov, 2024, Wang et al., 16 Aug 2025). This suggests that GPSF is best understood as a methodological category defined by generalized PSF representations and generalized inverse operators, rather than as a single canonical algorithm.

1. Scope of the term and major formulations

The literature uses GPSF for several non-identical but structurally related problems. In each case, the PSF enters a forward model whose parameters are estimated from measured or simulated data.

Domain	PSF representation	Inverse task
Fluorescence microscopy	Sparse Gaussian mixture	Variational calibration
Coherent imaging	Spatially variant kernel or transmission matrix	Generalized deconvolution
Partially occulted imaging	Segmented PSF weights	Semi-blind linear fitting
Neural PSF engineering	Implicit neural phase function	Basis-free PSF fitting
XFEL screen diagnostics	Simulated aberration and scintillator PSFs	Convolution-based beam fitting
CMB polarization	Beam-convolved source templates with full covariance	Multi-source subtraction

In fluorescence microscopy, the PSF is modeled as a finite mixture of multivariate Gaussians, with the number of active kernels controlling an explicit accuracy-efficiency trade-off (Samuylov et al., 2018). In coherent imaging, the generalization is spatial variance: the forward operator is a full kernel $H(x,y;x',y')$ rather than a shift-invariant convolution kernel (Lee et al., 2017). In partially occulted images, the PSF is partitioned into segments whose weights are directly estimated from a multi-linear system (Hofmeister et al., 2022). In basis-free neural PSF engineering, the fitted object is the pupil-plane phase function $\phi(u,v)=f_\theta(u,v)$, from which the PSF is generated by wave optics (Valouev, 2024). In XFEL diagnostics, the PSF is simulated from aberrations and scintillator geometry and then included in the fitting function for a Gaussian beam (Novokshonov, 2024). In CMB polarization, GPSF refers to pixel-domain fitting of overlapping polarized point sources using the full covariance, including off-diagonal terms (Wang et al., 16 Aug 2025).

A common misconception is that PSF fitting is synonymous with fitting a single Gaussian or a single shift-invariant blur kernel. The GPSF literature explicitly departs from that restriction.

2. Shared mathematical structure

Despite their different application domains, GPSF methods share a common pattern: a forward operator maps latent object or source parameters to observed data, and PSF parameters or PSF-induced templates are estimated by minimizing a likelihood or a regularized loss.

For sparse analog calibration in fluorescence microscopy, the fitted PSF is

$$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$

with $w_k\ge 0$ and $\sum_k w_k=1$. Calibration is formulated as

$$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$

and, under the convex relaxation $\phi(w_k)=|w_k|$, as

$\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$

Here $\lambda>0$ trades off fitting error versus sparsity of the mixture (Samuylov et al., 2018).

For coherent imaging without shift invariance, the forward model is

$$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$

or, after discretization,

$\phi(u,v)=f_\theta(u,v)$0

where $\phi(u,v)=f_\theta(u,v)$1 is generally dense. A standard regularized solution is

$\phi(u,v)=f_\theta(u,v)$2

This is the generalized image-deconvolution formulation based on spatially variant PSFs or an equivalent transmission matrix (Lee et al., 2017).

For partially occulted images, the PSF is partitioned into segments:

$\phi(u,v)=f_\theta(u,v)$3

which yields a linear system

$\phi(u,v)=f_\theta(u,v)$4

over occulted pixels. The weights are obtained by a non-negative least-squares or regularized fit, embedded in an iterative fit-deconvolution loop (Hofmeister et al., 2022).

For basis-free neural fitting, the wave-optics forward model is

$\phi(u,v)=f_\theta(u,v)$5

with $\phi(u,v)=f_\theta(u,v)$6 represented by a small multi-layer perceptron. Training minimizes point-wise MSE on the PSF (Valouev, 2024).

For CMB polarization, the fitted object is a local source-template amplitude vector in a Gaussian likelihood,

$\phi(u,v)=f_\theta(u,v)$7

where $\phi(u,v)=f_\theta(u,v)$8 is the full pixel-pixel covariance for $\phi(u,v)=f_\theta(u,v)$9 data, including off-diagonal terms (Wang et al., 16 Aug 2025).

This recurrence of forward modeling plus regularized or covariance-weighted inversion is the main mathematical unifier of GPSF.

3. Sparse analog GPSF in fluorescence microscopy

In fluorescence microscopy, GPSF was introduced to improve reconstruction accuracy beyond the Gaussian model while preserving computational efficiency. The core representation is a sparse combination of Gaussian kernels that retains analog, continuous modeling needed for sub-pixel localization accuracy. A single Gaussian typically models only the central lobe, whereas real PSFs have “fat” tails, including side-lobes and asymmetries; mixing several Gaussians with different covariances $$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$0 and shifts $$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$1 provides a flexible model with closed-form convolution and fast evaluation, for example via FFT or IFGT (Samuylov et al., 2018).

The inverse problem is posed variationally. The ideal penalty on the number of active kernels is combinatorial and nonconvex, so it is replaced by an $$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$2 surrogate on the mixture weights. After solving the convex problem, small $$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$3 are thresholded to zero and the remaining weights may be refit on the pure $$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$4 term under positivity. To solve the large-scale convex program efficiently, the method uses a fully-split formulation with an Alternating Split Bregman scheme. The split variables separately handle data fidelity, $$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$5 shrinkage, and nonnegativity. For fixed $$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$6, the split-Bregman or ADMM theory gives convergence; in practice, $$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$7–$$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$8 ASB iterations suffice (Samuylov et al., 2018).

The key design variable is the support size

$$h(x;\Theta)=\sum_{k=1}^M w_k\,\mathcal N(x;\mu_k,\Sigma_k),$$9

the number of active kernels. Accuracy may be measured by

$w_k\ge 0$0

or normalized deviance. Varying $w_k\ge 0$1 and the dictionary size $w_k\ge 0$2 traces out a curve $w_k\ge 0$3: typically $w_k\ge 0$4 decreases as $w_k\ge 0$5 grows, but computation time $w_k\ge 0$6 grows roughly linearly with $w_k\ge 0$7. The paper recommends inspecting the “elbow” in $w_k\ge 0$8 when choosing the model size (Samuylov et al., 2018).

Experimental validation is reported on synthetic Born–Wolf stacks and on real widefield/confocal data. The single-Gaussian model gives $w_k\ge 0$9–$\sum_k w_k=1$0 but under-/over-estimates intensity and fails in tails. GPSF with $\sum_k w_k=1$1–$\sum_k w_k=1$2 kernels halves the photometry error and slightly improves $\sum_k w_k=1$3. A full Born–Wolf model is best in accuracy but approximately $\sum_k w_k=1$4 slower. A representative timing-and-deviance comparison lists SG with $\sum_k w_k=1$5 and $\sum_k w_k=1$6, GPSF with small $\sum_k w_k=1$7 at $\sum_k w_k=1$8, $\sum_k w_k=1$9, $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$0, GPSF with medium $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$1 at $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$2, $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$3, $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$4, and BW at $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$5 with $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$6 (Samuylov et al., 2018).

4. Spatially variant and semi-blind GPSF

A second major line of work generalizes PSF fitting by abandoning shift invariance. In coherent optical imaging, the imaging system is represented by a spatially variant kernel $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$7 or, equivalently, by a transmission matrix. Full characterization is achieved by successively recording optical responses under various laser illumination angles controlled by a digital micro-mirror device. Input spectra vectors and output spectra are stacked into matrices $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$8 and $$\min_{\Theta} \|h(\Theta)*u-f\|_2^2+\lambda\sum_{k=1}^M \phi(w_k),$$9, and the transmission matrix is computed as

$\phi(w_k)=|w_k|$0

The inversion may be performed directly in the transmission-matrix formulation by SVD, or in object space through regularized least squares (Lee et al., 2017).

This generalized deconvolution is experimentally contrasted with conventional Fourier-domain deconvolution. For optical defocus, conventional and GPSF both achieve field-field correlation greater than $\phi(w_k)=|w_k|$1. For a tilted lens, conventional deconvolution gives approximately $\phi(w_k)=|w_k|$2, whereas GPSF gives approximately $\phi(w_k)=|w_k|$3. For a scattering layer, conventional deconvolution falls below $\phi(w_k)=|w_k|$4 and can be worse than the raw image, whereas GPSF achieves approximately $\phi(w_k)=|w_k|$5–$\phi(w_k)=|w_k|$6. The restored effective resolution reaches $\phi(w_k)=|w_k|$7 line pairs, approaching the diffraction limit. Practical limitations include transmission-matrix measurement time, reported as approximately $\phi(w_k)=|w_k|$8 for $\phi(w_k)=|w_k|$9 angle steps at $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$0, the memory and compute cost of a dense $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$1, and noise amplification from small singular values, handled by $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$2-regularization (Lee et al., 2017).

A different semi-blind generalization derives the instrumental PSF from partially occulted images. Here the PSF is partitioned into nonoverlapping regions $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$3, with unknown weights $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$4. The occulted pixels provide linear equations in those weights, and the algorithm alternates between fitting $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$5, reassembling the PSF, deconvolving the full image, and updating the current approximation of the true image under the mask constraint (Hofmeister et al., 2022).

The reported guarantees are unusually explicit. The algorithm is described as guaranteed to converge towards the correct instrumental PSF for a large class of occultations, and, under mild conditions on the mask and segmentation, the solution is unique in the noise-free setting (Hofmeister et al., 2022). Quantitatively, the central weight of the PSF, interpreted as the fraction of photons not scattered by the instrument, is accurate to better than $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$6. The mean absolute percentage error between reconstructed and true PSF is usually between $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$7 and $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$8 for the entire PSF, between $\min_{w\ge 0}\|h(w)*u-f\|_2^2+\lambda\|w\|_1.$9 and $\lambda>0$0 for the PSF core, and between $\lambda>0$1 and $\lambda>0$2 for the PSF tail (Hofmeister et al., 2022).

Together, these two formulations show that “generalized” may refer either to a full space-variant forward operator or to a semi-blind segmented PSF estimate without a predefined functional form.

5. Basis-free neural and simulation-constrained fitting

A more recent extension generalizes PSF fitting at the level of the pupil function. In basis-free point spread function engineering, the phase is represented implicitly as

$\lambda>0$3

where $\lambda>0$4 is a small multi-layer perceptron. The predicted PSF is

$\lambda>0$5

with uniform pupil amplitude. No explicit Tikhonov or TV regularizer is added; the network’s smoothness serves as an implicit regularization. The input uses “radial encoding” with $\lambda>0$6 and encoding length $\lambda>0$7, and the network uses $\lambda>0$8–$\lambda>0$9 hidden fully connected layers of width $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$0 with Leaky ReLU. Training is performed with Adam at initial learning rate $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$1, cosine-annealed to $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$2 over $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$3 epochs, with batch size $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$4 PSF per parameter update (Valouev, 2024).

For random pupils, the method reports mean MSSIM $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$5 versus $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$6 with pixel-wise optimization, median PSNR $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$7 versus $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$8, and mean PSNR $$g(x,y)=\iint H(x,y;x',y')\,f(x',y')\,dx'\,dy'+n(x,y),$$9 versus $\phi(u,v)=f_\theta(u,v)$00. For Zernike pupils, the reported median MSSIM is $\phi(u,v)=f_\theta(u,v)$01 versus $\phi(u,v)=f_\theta(u,v)$02, while mean PSNR is $\phi(u,v)=f_\theta(u,v)$03 versus $\phi(u,v)=f_\theta(u,v)$04. An ablation on EDoF reports that pixel-wise optimization wins on very high-frequency phase, but with much more noise and artifacts; the MLP produces smooth, manufacturable $\phi(u,v)=f_\theta(u,v)$05 (Valouev, 2024). The same summary presents this as a GPSF pipeline: collect target PSFs, build the differentiable forward operator, train with Adam plus cosine annealing for $\phi(u,v)=f_\theta(u,v)$06–$\phi(u,v)=f_\theta(u,v)$07 epochs, monitor MSSIM or PSNR, and export the learned phase at the target SLM resolution.

Another formulation is simulation-constrained rather than basis-free. For European XFEL scintillating screens, the PSF is decomposed into optical aberrations and scintillator-induced geometric spreading. The aberration term $\phi(u,v)=f_\theta(u,v)$08 is approximated by a 2D Gaussian with width $\phi(u,v)=f_\theta(u,v)$09 for the Schneider-Kreuznach $\phi(u,v)=f_\theta(u,v)$10 objective and $\phi(u,v)=f_\theta(u,v)$11 for the $\phi(u,v)=f_\theta(u,v)$12 objective. The scintillator term $\phi(u,v)=f_\theta(u,v)$13 is fitted by a super-Gaussian comb of four lobes with power $\phi(u,v)=f_\theta(u,v)$14. Assuming a 2D Gaussian beam profile

$\phi(u,v)=f_\theta(u,v)$15

the observed image is modeled as

$\phi(u,v)=f_\theta(u,v)$16

or, with background,

$\phi(u,v)=f_\theta(u,v)$17

The free parameters are the beam centroid, beam widths, and background, while the PSF parameters are held fixed from Zemax fits and estimated by non-linear least squares, for example Levenberg–Marquardt (Novokshonov, 2024).

These two cases enlarge GPSF in different directions: one removes the basis restriction in the pupil plane, and the other makes PSF simulation itself part of the fitting function.

6. Covariance-aware multi-source GPSF in CMB polarization

In CMB B-mode analysis, GPSF refers to a map-domain method for removing polarized point-source contamination. The data vector collects local $\phi(u,v)=f_\theta(u,v)$18 and $\phi(u,v)=f_\theta(u,v)$19 pixels within a fitting region, and the background is modeled as Gaussian with full pixel-pixel covariance

$\phi(u,v)=f_\theta(u,v)$20

including off-diagonal terms between pixels and between Stokes components. A source at known position is described by a beam-convolved template, and overlapping sources are fitted jointly by minimizing the generalized $\phi(u,v)=f_\theta(u,v)$21 over the source amplitudes (Wang et al., 16 Aug 2025).

The implementation proceeds in two stages. First, candidate sources are identified from an input catalog, sorted by polarized flux $\phi(u,v)=f_\theta(u,v)$22, and tested in a local model including neighbors within $\phi(u,v)=f_\theta(u,v)$23. Sources above a significance threshold $\phi(u,v)=f_\theta(u,v)$24–$\phi(u,v)=f_\theta(u,v)$25 are retained. Second, nearby sources are partitioned into disjoint groups and fit simultaneously with the same covariance $\phi(u,v)=f_\theta(u,v)$26; the best-fit templates are then subtracted from the observed maps. The method is designed for realistic conditions in which source blending is important, especially for small-aperture telescopes with large beams (Wang et al., 16 Aug 2025).

The method is integrated with Needlet Internal Linear Combination (NILC): GPSF is applied first on the multi-frequency $\phi(u,v)=f_\theta(u,v)$27 maps, after which the maps are smoothed to a common beam of $\phi(u,v)=f_\theta(u,v)$28 FWHM, corrected for EB leakage via the recycling method, decomposed into needlets, and combined using the standard ILC weight

$\phi(u,v)=f_\theta(u,v)$29

The purpose is to remove resolved polarized sources before diffuse-component cleaning (Wang et al., 16 Aug 2025).

On simulations, the reported performance is explicit. For the tensor-to-scalar ratio $\phi(u,v)=f_\theta(u,v)$30, “PS unmitigated” yields mean bias $\phi(u,v)=f_\theta(u,v)$31 with $\phi(u,v)=f_\theta(u,v)$32. GPSF reduces the bias to $\phi(u,v)=f_\theta(u,v)$33, an $\phi(u,v)=f_\theta(u,v)$34 reduction, with $\phi(u,v)=f_\theta(u,v)$35, described as only a $\phi(u,v)=f_\theta(u,v)$36 increase. Masking gives $\phi(u,v)=f_\theta(u,v)$37 with $\phi(u,v)=f_\theta(u,v)$38, while inpainting yields $\phi(u,v)=f_\theta(u,v)$39 with $\phi(u,v)=f_\theta(u,v)$40. At the spectrum level, GPSF reduces residual $\phi(u,v)=f_\theta(u,v)$41 power by approximately $\phi(u,v)=f_\theta(u,v)$42 at $\phi(u,v)=f_\theta(u,v)$43 with nearly unchanged variance, and, before NILC, yields the lowest residual across $\phi(u,v)=f_\theta(u,v)$44 and $\phi(u,v)=f_\theta(u,v)$45 with standard deviations at the ideal level (Wang et al., 16 Aug 2025).

This usage of GPSF is notable because the generalized object is neither the PSF itself nor the pupil function, but the statistical treatment of PSF-convolved source templates through a full covariance model.

7. Recurring design choices, limitations, and directions

Across the literature, GPSF methods recurrently exchange restrictive optics assumptions for richer representations and more expensive inversions. In sparse Gaussian mixtures, the principal control knob is the support size $\phi(u,v)=f_\theta(u,v)$46, with an elbow search on $\phi(u,v)=f_\theta(u,v)$47 recommended for choosing the accuracy-efficiency operating point (Samuylov et al., 2018). In spatially variant deconvolution, the main burden is the dense operator or transmission matrix, together with measurement time and the conditioning of the inverse (Lee et al., 2017). In segmented semi-blind recovery, the practical issues are segmentation granularity, mask design, regularization choice, clustering of distant occulted pixels into super-pixels, bootstrap averaging, and hierarchical refinement (Hofmeister et al., 2022). In neural phase fitting, the stated limitation is the decay of the standard MLP plus positional-encoding spectrum at high frequencies; suggested remedies include periodic SIREN activations, higher-frequency Fourier features, transfer-learned initializations, and end-to-end optimization of upstream imaging metrics (Valouev, 2024). In CMB applications, performance depends on accurate beam knowledge, reliable input catalogs, and stable covariance estimation, and may degrade if source density is extremely high or beam non-Gaussianities are large beyond the symmetric Gaussian PSF assumption (Wang et al., 16 Aug 2025).

The most important interpretive point is terminological. “Generalized Point Spread Function Fitting” does not identify a single universally standardized formalism. Instead, it names a set of strategies that generalize PSF fitting along one or more axes: representation, spatial variance, blindness, physical forward modeling, basis choice, or statistical weighting. This suggests that the term functions as an umbrella for methods that keep the PSF—or PSF-convolved templates—inside the inverse problem at a higher fidelity than classical single-kernel fitting.