Implementing MATLAB’s psf2otf Function in Python: Why and How

The Point Spread Function (PSF) and the Optical Transfer Function (OTF) are fundamental tools in image processing and the design and analysis of imaging systems. They allow us to understand and manipulate the way that light from the object space is transformed into an image, and to correct for distortions and blurring caused by the imaging process.

The PSF describes the response of an imaging system to a point source or point object. It represents the distribution of light from a single point in the object space, as it appears in the image. The PSF often represents the blurring caused by the imaging system, such as the blur caused by out-of-focus optics or motion blur in a camera, or the diffraction pattern caused by the aperture in a microscope.

The OTF, on the other hand, is the Fourier transform of the PSF. It describes how different spatial frequencies are handled by the imaging system. The magnitude of the OTF (the Modulation Transfer Function, or MTF) describes how much each frequency is attenuated, while the phase of the OTF describes how each frequency is phase-shifted.

In MATLAB, the psf2otf function is used to convert a PSF to an OTF. This function is often used in image processing tasks, particularly in deconvolution operations. In this article, we will discuss how to implement an equivalent function in Python.

The psf2otf function in MATLAB works by padding, circularly shifting (to ensure the center of the PSF is at (0,0)), and taking the Fourier Transform of the PSF to generate the OTF.

Let’s break down these steps and translate them into Python code using the NumPy and SciPy libraries.

1. Padding: The PSF needs to be padded to the same size as the image that it will be applied to. This can be done using the numpy.pad function.

import numpy as np

def pad_psf_to_shape(psf, shape):
    pad_sizes = [(s - p) // 2 for s, p in zip(shape, psf.shape)]
    pad_sizes = [(p, s - p) for p, s in zip(pad_sizes, shape)]
    return np.pad(psf, pad_sizes, mode='constant')

2. Circular Shift: The PSF needs to be circularly shifted so that its center is at (0,0). This can be done using the numpy.roll function.

def circular_shift_psf(psf, shift):
    return np.roll(psf, shift, axis=(0, 1))

3. Fourier Transform: Finally, the Fourier Transform of the PSF is taken to generate the OTF. This can be done using the numpy.fft.fft2 function.

def compute_otf(psf):
    return np.fft.fft2(psf)

Now, let’s put these steps together to create a Python function equivalent to MATLAB’s psf2otf.

def psf2otf(psf, shape):
    # Pad the PSF to the desired shape
    psf = pad_psf_to_shape(psf, shape)

    # Compute the shift needed to move the center of the PSF to (0,0)
    shift = [dim // 2 for dim in psf.shape]

    # Circularly shift the PSF
    psf = circular_shift_psf(psf, shift)

    # Compute the OTF
    otf = compute_otf(psf)

    return otf

This Python function should provide equivalent functionality to MATLAB’s psf2otf. It’s important to note that the PSF should be normalized (i.e., its sum should be 1) before being passed to this function, as is the case with the MATLAB function.

This function can be used in image processing tasks in Python in the same way that the psf2otf function is used in MATLAB. For example, it can be used in deconvolution operations to recover the original image from a blurred image and the PSF used to blur it.

Applications of the PSF and OTF are vast. In astronomy, the PSF can be used to correct for atmospheric distortion in images taken by ground-based telescopes. The PSF representing the distortion is measured (often by observing a known star), and then deconvolution is used to correct the images. In medical imaging, such as MRI or CT scans, the PSF can be used to correct for the blurring caused by the imaging process. This can help to produce sharper images, which can aid in diagnosis.

In conclusion, the PSF and OTF are fundamental tools in image processing and the design and analysis of imaging systems. By implementing MATLAB’s psf2otf function in Python, we can leverage these tools in a Python environment, opening up a wide range of possibilities for image analysis and manipulation.