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
penaltyis 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
penaltyarray 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
penaltyoperator. 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
regularizationparameter, and \(W\) is the optionalwindowfunction (defaulting to 1 if not provided). \(|S|^{2}\) is the power spectrum of the source OTF.When the
penaltyis set to'laplacian', the regularization reproduces the approach used by thepypherpackage (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_psfandtarget_psfmust 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_psfandtarget_psfmust 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
penaltyisNone, this is expressed as a fraction of the peak power in the source PSF’s Fourier transform. Whenpenaltyis 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.- penalty
None,'laplacian','biharmonic', or 2Dndarray, 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 thepypherpackage (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
shapeparameter (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.
- source_psf2D
- Returns:
- kernel2D
ndarray The matching kernel to go from
source_psftotarget_psf. The output matching kernel is normalized such that it sums to 1.
- kernel2D
- 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
regularizationis not positive, or ifpenaltyis not a valid value.- TypeError
If the input
windowis not callable.
See also
make_kernelMake 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