5.6. PSF Convolution

5.6.1. Convolution of Colored Images

5.6.1.1. Polychromatic Convolution

In general, both the image and the PSF are wavelength-dependent and the image convolution must be performed individually for each wavelength:

(5.154)\[\text{im}_2(x, y, \lambda) = \text{im}(x, y, \lambda) \otimes \text{psf}(x, y, \lambda)\]

5.6.1.2. Special Case Spectral Homogenity

In the case of so-called spectral homogeneity, the spectrum \(s\) is consistent across the entire image but is scaled by a location-dependent intensity factor \(\text{im}_s(x, y)\). The concept of spectral homogeneity goes back to [1], who also describes this mathematical simplification. [2] also uses this to generate color images.

For an example \(\text{im}_s(\lambda)\), it holds that:

(5.155)\[\begin{split}\text{im}_2(x, y, \lambda) &= \left(\text{im}_r(\lambda) \text{im}_s(x, y)\right) \otimes \text{psf}(x, y, \lambda)\\ &= \text{im}_s(x, y) \otimes \left(\text{im}_s(\lambda) \text{psf}(x, y, \lambda)\right)\\ &= \text{im}_s(x, y) \otimes \text{psf}_s(x, y, \lambda)\\\end{split}\]

The spectrum is a constant with respect to the convolution expression and can therefore also be multiplied with the PSF. We define the new spectrally weighted PSF \(\text{psf}_s(x, y, \lambda)\).

5.6.1.3. Convolution of RGB Images

Typically, the spectral profile of the image result is not of interest. Instead, the color coordinates for display on a monitor are more relevant. For instance, to calculate the red stimulus produced by the convolved image, it can be multiplied by the sRGB red color matching function \(r(\lambda)\) and integrated over all wavelengths:

(5.156)\[\begin{split}\text{im}_{2,r\rightarrow r}(x, y) &= \int \Big[\text{im}_r(x, y) \otimes \text{psf}_r(x, y, \lambda)\Big] r(\lambda)\;\text{d}\lambda\\ &= \text{im}_r(x, y) \otimes \int\text{psf}_r(x, y, \lambda) r(\lambda)\;\text{d}\lambda\\ &= \text{im}_r(x, y) \otimes \text{psf}_{r\rightarrow r}(x, y)\\\end{split}\]

Since only the PSF has a wavelength dependency, the image can be factored out of the integral expression. Because the convolution is independent of wavelength, the integration can be performed prior to convolution. This results in the red-to-red PSF \(\text{psf}_{r \rightarrow r}(x, y)\).

If the light spectra \(\text{im}_r(\lambda), \text{im}_g(\lambda), \text{im}_b(\lambda)\) are selected such that they correspond to the primary spectra of sRGB, they form linearly independent color channels. These channels can be used to compose all colors within the sRGB color space via linear combination. The PSF can similarly be decomposed into three individual channel PSFs that are linearly independent.

The RGB color image from the red image \(\text{im}_r(x, y)\) with spatially homogeneous \(\text{im}_r(\lambda)\) for convolution with the PSF can then be simplified and summarized as:

(5.157)\[\begin{split}\text{im}_{2,r\rightarrow rgb}(x, y) = \left[\begin{array}\, \text{im}_r(x, y) \otimes \text{psf}_{r\rightarrow r}(x, y)\\ \text{im}_r(x, y) \otimes \text{psf}_{r\rightarrow g}(x, y)\\ \text{im}_r(x, y) \otimes \text{psf}_{r\rightarrow b}(x, y)\\ \end{array}\right]\end{split}\]

Even if the original red spectrum generates a pure red in the sRGB color space, this might not be the case after convolution with the PSF. Chromatic aberration could lead to more transverse chromatic aberration at smaller wavelengths, resulting in a yellow fringe in the PSF and introducing components in \(\text{psf}_{r \rightarrow g}(x, y)\).

Similarly, this applies to the G-channel with the sRGB color matching function \(g(\lambda)\):

(5.158)\[\begin{split}\text{im}_{2,g\rightarrow rgb}(x, y) = \left[\begin{array}\, \text{im}_g(x, y) \otimes \text{psf}_{g\rightarrow r}(x, y)\\ \text{im}_g(x, y) \otimes \text{psf}_{g\rightarrow g}(x, y)\\ \text{im}_g(x, y) \otimes \text{psf}_{g\rightarrow b}(x, y)\\ \end{array}\right]\end{split}\]

And the blue channel with matching function \(b(\lambda)\):

(5.159)\[\begin{split}\text{im}_{2,b\rightarrow rgb}(x, y) = \left[\begin{array}\, \text{im}_b(x, y) \otimes \text{psf}_{b\rightarrow r}(x, y)\\ \text{im}_b(x, y) \otimes \text{psf}_{b\rightarrow g}(x, y)\\ \text{im}_b(x, y) \otimes \text{psf}_{b\rightarrow b}(x, y)\\ \end{array}\right]\end{split}\]

The overall image \(\text{im}_{2,rgb \rightarrow rgb}\) is obtained from the sum of all convolved R, G, B color components in the image. However, the mixing ratio of all channels must be considered: If the color PSFs were all simulated with a power of one watt, this does not correspond to the correct mixing ratio for white in the sRGB color space. This must be adjusted so that equal parts in \(\text{im}_r, \text{im}_g, \text{im}_b\) produce white in the color space.

Let \(a_r, a_g, a_b\) be the relative mixing factors. The final result can be expressed as:

(5.160)\[\text{im}_{2,rgb\rightarrow rgb}(x, y) = a_r \text{im}_{2,r\rightarrow rgb}(x, y) + a_g \text{im}_{2,g\rightarrow rgb}(x, y) + a_b \text{im}_{2,b\rightarrow rgb}(x, y)\]

Equation (5.188) illustrates the RGB color spectra and respective weighting factors. Rather than applying the rescaling subsequently, the power ratios might be integrated into the source rendering process of the PSFs. By doing this, the power ratios are inherently incorporated into the relative R, G, B PSFs.

5.6.1.4. Convolution of a spectral homogeneous image and an RGB PSF

In the special case where the image is spectrally homogeneous, let the spectrum be denoted as \(s\). The transformation can then be represented by:

(5.161)\[\begin{split}\text{im}_{2,s\rightarrow rgb}(x, y) = \left[\begin{array}\, \text{im}_s(x, y) \otimes \text{psf}_{s\rightarrow r}(x, y)\\ \text{im}_s(x, y) \otimes \text{psf}_{s\rightarrow g}(x, y)\\ \text{im}_s(x, y) \otimes \text{psf}_{s\rightarrow b}(x, y)\\ \end{array}\right]\end{split}\]

In this context, \(\text{im}_s\) describes the spatial distribution of the intensity of the source, which emits the spectrum \(s\).

5.6.1.5. Convolution of spectral homogeneous image and PSF

In the special scenario where both the image and the PSF are spectrally homogeneous, the relationship can be defined as:

(5.162)\[\text{im}_{2,s_1\rightarrow s_2}(x, y) = \text{im}_{s_1}(x, y) \otimes \text{psf}_{s_1\rightarrow s_2}(x, y)\]

In this equation, \(s_1\) represents the source spectrum, and \(s_2\) represents the detector spectrum. Due to potential absorption effects, these spectra do not need to be identical. As explained earlier, it is crucial that both the source and the PSF exhibit spectral homogeneity. A typical example involves a black and white image in conjunction with a wavelength-independent PSF.

5.6.2. Limitations

The limitations are detailed in Section 4.11.3.

5.6.3. Processing Steps

1. Convert the image and PSF to linear sRGB values, while including negative values.

2. For a grayscale PSF, normalize the PSF so that its sum (total power) equals one.

3. Downscale/interpolate the PSF so that the physical pixel sizes of the PSF and the image (after scaling with the magnification factor) are identical.

4. Pad the PSF with zeros to ensure a defined fall-off

5. Flip the image if the magnification factor is negative

6. Pad the image according to the chosen padding method.

7. Convolve image and PSF according to the methods in Section Section 5.6.1.

8. Convert the image back to sRGB, applying a selected gamut mapping.

9. Slice the image to remove padding, or trim it back to its original size

The convolution is conducted in sRGB coordinates because the channels are orthogonal, and this color space is the target for monitors. Nevertheless, the convolution must be executed as a linear operation using linear sRGB values. It is also crucial to include colors outside the color space (negative coordinates) to maintain linearity. If negative values persist in the image post-convolution, gamut mapping needs to be applied.

In the case of a grayscale PSF, it is automatically normalized. This ensures that, when used in conjunction with the normalize=False parameter of the convolve function, the brightness and color values remain unchanged. For colored PSFs, normalization poses a greater challenge because it requires knowing the amount of light from the source that actually reached the detector for the PSF. This could potentially be achieved using metadata from RenderImage. Nonetheless, the relevance of this option is debatable, as normalized images are typically preferred.

Downscaling the PSF needs to be executed in a manner that conserves energy. Moreover, it is essential to choose a method that prevents aliasing. We utilize the scaling with the INTER_AREA option from OpenCV within the resize function. The PSF must be rescaled so that the physical pixel dimensions of the image and the PSF align in both dimensions. Consequently, it becomes sufficient to convolve the image as a pixel matrix, even if the pixels are not square.

The convolution is performed in the Fourier space employing the convolution theorem. The Fourier transformation is calculated using the scipy.signal.fftconvolve() function. Due to the nature of this method, areas outside the image are assumed to be black. Consequently, a fall-off region appears in the resulting image, where the PSF increasingly convolves with the black portions at the edges. This transition area is as wide as the PSF region where its intensities exceed zero. For simplicity, we assume the entire PSF width for this. If the user desires a different padding method, additional padding must be applied to the image. First, the image should be padded according to the user-chosen method, and then again to prevent dark fall-off edges.

---

References

[1]

R. Barnden. Calculation of axial polychromatic optical transfer function. Optica Acta: International Journal of Optics, 21(12):981–1003, dec 1974. doi:10.1080/713818864.

[2]

Sowmya Ravikumar, Larry N. Thibos, and Arthur Bradley. Calculation of retinal image quality for polychromatic light. Journal of the Optical Society of America A, 25(10):2395, sep 2008. URL: https://doi.org/10.1364%2Fjosaa.25.002395, doi:10.1364/josaa.25.002395.