make_wiener_kernel#

photutils.psf_matching.make_wiener_kernel(source_psf, target_psf, *, regularization=0.0001, penalty=None, window=None)[source]#

Make a convolution kernel that matches an input PSF to a target PSF using Wiener regularization in Fourier space.

This function computes a Wiener-regularized PSF-matching kernel in the Fourier domain. The denominator includes a regularization term that stabilizes inversion of the source OTF (Optical Transfer Function, the Fourier transform of the PSF) by preventing division by small values, thereby suppressing noise amplification at spatial frequencies where the source response is weak.

When no penalty is provided, the regularization is a frequency-independent (zero-order scalar) Tikhonov term expressed as a fraction of the peak power in the source OTF. In this case, the kernel is computed as:

\[K = \mathcal{F}^{-1} \left[ W \cdot \frac{T \cdot S^{*}} {|S|^{2} + \lambda \cdot \max(|S|^{2})} \right]\]

When a penalty array is provided (e.g., a Laplacian operator), the regularization becomes frequency-dependent:

\[K = \mathcal{F}^{-1} \left[ W \cdot \frac{T \cdot S^{*}} {|S|^{2} + \lambda \cdot |P|^{2}} \right]\]

where \(P\) is the OTF of the penalty operator. This penalizes high spatial frequencies more heavily, which is particularly effective at suppressing noise amplification.

In both equations, \(\mathcal{F}^{-1}\) is the inverse Fourier transform, \(S\) and \(T\) are the Fourier transforms of the source and target PSFs (the OTFs), \(S^{*}\) is the complex conjugate of \(S\), \(\lambda\) is the regularization parameter, and \(W\) is the optional window function (defaulting to 1 if not provided). \(|S|^{2}\) is the power spectrum of the source OTF.

When the penalty is set to 'laplacian', the regularization reproduces the approach used by the pypher package (Boucaud et al. 2016), which applies a discrete Laplacian operator as the penalty. This provides stronger suppression of high spatial frequencies, which can be beneficial when working with noisy or undersampled PSFs.

Compared to make_kernel, which uses a hard threshold on Fourier amplitude, this approach provides continuous, smooth regularization that is better suited for PSFs that have near-zero power at high spatial frequencies. The hard-cutoff approach zeros out frequencies where the source amplitude is below a threshold, which can introduce discontinuities in Fourier space. Wiener regularization instead smoothly down-weights those frequencies, producing matching kernels with less ringing.

Parameters:
source_psf2D ndarray

The source PSF. The source PSF should have higher resolution (i.e., narrower) than the target PSF. source_psf and target_psf must have the same shape and pixel scale. It is assumed to be centered on the central pixel.

target_psf2D ndarray

The target PSF. The target PSF should have lower resolution (i.e., broader) than the source PSF. source_psf and target_psf must have the same shape and pixel scale. It is assumed to be centered on the central pixel.

regularizationfloat, optional

The regularization parameter that controls the strength of the Wiener (Tikhonov) regularization. When penalty is None, this is expressed as a fraction of the peak power in the source PSF’s Fourier transform. When penalty is provided, this scales the penalty operator’s power spectrum directly. Larger values produce smoother but less accurate matching kernels; smaller values preserve more detail but may amplify noise. Must be a positive number.

penaltyNone, 'laplacian', 'biharmonic', or 2D ndarray, optional

The regularization penalty operator. This controls the structure of the regularization term in the denominator:

  • None (default): Scalar Tikhonov regularization. The denominator is \(|S|^2 + \lambda \cdot \max(|S|^2)\), providing uniform regularization across all spatial frequencies. Use this for well-behaved PSFs or when you want simple, frequency-independent smoothing.

  • 'laplacian': Uses a discrete Laplacian operator (second derivative) as the penalty, producing frequency-dependent regularization that penalizes high spatial frequencies more heavily. The denominator becomes \(|S|^2 + \lambda \cdot |L|^2\) where \(L\) is the OTF of the Laplacian kernel [[0, -1, 0], [-1, 4, -1], [0, -1, 0]]. This reproduces the regularization used by the pypher package (Boucaud et al. 2016). This is the most commonly used penalty for PSF matching and works well for most applications. Requires PSFs to be at least 3x3.

  • 'biharmonic': Uses a biharmonic operator (fourth derivative, Laplacian of the Laplacian) as the penalty, producing very strong suppression of high spatial frequencies. Uses the kernel [[0, 0, 1, 0, 0], [0, 2, -8, 2, 0], [1, -8, 20, -8, 1], [0, 2, -8, 2, 0], [0, 0, 1, 0, 0]]. This produces the smoothest matching kernels and is useful when working with very noisy or poorly sampled PSFs, at the cost of reduced accuracy in matching. Requires PSFs to be at least 5x5.

  • 2D ndarray: A custom penalty operator array. Its OTF will be computed and used in the denominator as \(|S|^2 + \lambda \cdot |P|^2\). The PSFs must be at least as large as the penalty array in both dimensions.

windowcallable, optional

The window (taper) function or callable class instance used to remove high frequency noise from the PSF matching kernel. The window function should be a callable that accepts a single shape parameter (a tuple defining the 2D array shape) and returns a 2D array of the same shape. The returned window values must be in the range [0, 1], where 1.0 indicates full preservation of that spatial frequency and 0.0 indicates complete suppression. The window must be centered on the central pixel. Built-in window classes include:

A window function is generally not needed when using Wiener regularization because the regularization itself suppresses high-frequency noise. However, a window may still be useful when working with noisy or undersampled PSFs.

For more information on window functions, custom windows, and example usage, see PSF Matching.

Returns:
kernel2D ndarray

The matching kernel to go from source_psf to target_psf. The output matching kernel is normalized such that it sums to 1.

Raises:
ValueError

If the PSFs are not 2D arrays, have even dimensions, do not have the same shape, are too small for the specified penalty, if regularization is not positive, or if penalty is not a valid value.

TypeError

If the input window is not callable.

See also

make_kernel

Make a matching kernel using a hard frequency cutoff instead of Wiener regularization.

Examples

Make a matching kernel between two Gaussian PSFs:

>>> import numpy as np
>>> from astropy.modeling.models import Gaussian2D
>>> from photutils.psf_matching import make_wiener_kernel
>>> y, x = np.mgrid[0:51, 0:51]
>>> psf1 = Gaussian2D(100, 25, 25, 3, 3)(x, y)
>>> psf2 = Gaussian2D(100, 25, 25, 5, 5)(x, y)
>>> psf1 /= psf1.sum()
>>> psf2 /= psf2.sum()
>>> kernel = make_wiener_kernel(psf1, psf2)
>>> print(f'{kernel.sum():.1f}')
1.0

Use the Laplacian penalty for frequency-dependent regularization:

>>> kernel = make_wiener_kernel(psf1, psf2, penalty='laplacian')
>>> print(f'{kernel.sum():.1f}')
1.0

Use the biharmonic penalty for maximum smoothness:

>>> kernel = make_wiener_kernel(psf1, psf2, penalty='biharmonic')
>>> print(f'{kernel.sum():.1f}')
1.0