Lucy Richardson output mostly white lines

[BLOO123]

Description

Way to reproduce

# Place the full code we need to recreate your issue here
# upload all necessary images to github too!

Code

    # Lucy Richardson
    
    import matplotlib.pyplot as graph
    import numpy as np
    from numpy import fft
    import math
    import cv2
    import skimage.restoration
    
    #Simulation motion blur
    def motion_process(image_size, motion_angle):
        psf = np.zeros(image_size)
        center_position = (image_size [0] -1)/2
        slope_tan =math.tan(motion_angle * math.pi/180)
        slope_cot = 1/slope_tan
        if slope_tan<= 1:
            for i in range(15):
                offset= round(i * slope_tan) #((center_position-i) * slope_tan)
                psf[int(center_position + offset), int(center_position-offset)]=1
            return psf/psf.sum() #normalize the brightness of the point spread function
        else:
            for i in range (15):
                offset=round(i*slope_cot)
                psf[int(center_position-offset), int(center_position + offset)]=1
            return psf/psf.sum()
    
    # Motion blur the picture
    def make_blurred(input, psf, eps):
        input_fft=fft.fft2(input) #Fourier transform a two-dimensional array
        psf_fft=fft.fft2(psf) + eps # Added to remove zeroes
        blurred=fft.ifft2(input_fft*psf_fft)
        blurred=np.abs(fft.fftshift(blurred))
        return blurred
    
    image=cv2.imread("face.jpg",0)
    img_h=image.shape[0]
    img_w=image.shape[1]
    
    graph.figure(1)
    graph.xlabel("original image")
    graph.gray()
    graph.imshow(image) #Show the original image
    
    graph.figure(2)
    graph.gray()
    
    #Perform motion blur
    psf= motion_process((img_h, img_w), 60) # Get PSF
    blurred = np.abs(make_blurred(image, psf, 1e-3)) # Get power spectrum of image after convolving with psf
    graph.subplot(241)
    graph.xlabel("motion blurred")
    graph.imshow(blurred)
    
    # Get the inverse filter of motion blurred image
    result = skimage.restoration.richardson_lucy(blurred, psf, iterations=80,clip=True,filter_epsilon=1e-3)
    graph.subplot(242)
    graph.xlabel("LR deblurred (NUMIT = 50)")
    graph.imshow(result)
    graph.show()
    
    # Add noise to blurred image
    blurred_noisy=blurred + 0.1 * blurred.std() * \
       np.random.standard_normal(blurred.shape) #Add noise,standard_normal generates a random function
    graph.subplot(243)
    graph.xlabel("motion&noisy blurred")
    graph.imshow(blurred_noisy) #Show images with added noise and motion blur
    
    # Get inverse filter of final image
    result = skimage.restoration.richardson_lucy(blurred_noisy, psf, iterations=80,clip=True,filter_epsilon=1e-3)
    graph.subplot(244)
    graph.xlabel("LR deblurred (NUMIT = 50)")
    graph.imshow(result)
    graph.show()

Output

The image below shows the output from the code above.

If the filter_epsilon was set to 0, the output image would be completely white.

The image below shows the output when iteration is set to 30. The image has more white lines than at lower iteration.

Is there a mistake in my implementation of the lucy richardson restoration? Thank you!

Version information

# Paste the output of the following python commands
from __future__ import print_function
import sys; print(sys.version)
import platform; print(platform.platform())
import skimage; print("scikit-image version: {}".format(skimage.__version__))
import numpy; print("numpy version: {}".format(numpy.__version__))

# your output here
3.8.8 (default, Apr 13 2021, 15:08:03) [MSC v.1916 64 bit (AMD64)]
Windows-10-10.0.19041-SP0
scikit-image version: 0.18.1
numpy version: 1.20.2

[mkcor]

Hi @BLOO123,

Is there a mistake in my implementation of the lucy richardson restoration? Thank you!

I believe your question belongs better in the user forum: https://forum.image.sc/ 🙏