How to improve a variable lens blur algorithm in Python OpenCV?

The only way to implement a filter where the kernel is different for every pixel is to create the kernel for each pixel and apply it in a loop, like OP's code does. The Fourier transform does not apply to this case. Python is a very slow language, the same algorithm implemented in a compiled language would be much faster. Unless there is some predefined structure in how the kernel is created at each pixel, there is no way to reduce the complexity of the algorithm.

For example, the uniform filter with a square kernel (commonly called the "box" filter) can be computed based on the integral image, using only 4 additions per pixel. This implementation should be able to choose a different kernel size at each pixel without any additional cost.

DIPlib has an implementation of an adaptive Gaussian filter [disclaimer: I'm an author of DIPlib, but I did not implement this functionality]. Here is the documentation. This filter applies a Gaussian filter, but the Gaussian kernel is scaled and rotated differently at every pixel.

Lens blur is not a Gaussian, but it's not easy to see the difference by eye in most cases; the difference matters only if there is a very small dot with high contrast.

OP's case would be implemented as follows:

import diplib as dip

img = dip.ImageRead('examples/trui.ics')

blur_map = dip.CreateRadiusSquareCoordinate(img.Sizes())
blur_map /= dip.Maximum(blur_map)

img_blur = dip.AdaptiveGauss(img, [0, blur_map], sigmas=[5])

(the blur_map here is defined differently, I chose a quadratic function of the distance to the center, because I think it looks really nice; use dip.CreateRadiusCoordinate() to reproduce OP's map).

I've chosen a maximum blur of 5 (this is the sigma, in pixels, of the Gaussian, not the footprint of the kernel), and blur_map here scales this sigma with a factor between 0 in the middle and 1 at the corners of the image.

Another interesting effect would be as follows, with increasing blur tangential to each circle centered in the middle of the image, with very little blur radially:

angle_map = dip.CreatePhiCoordinate(img.Sizes())
img_blur = dip.AdaptiveGauss(img, [angle_map, blur_map], sigmas=[8,1])

Here is one way to apply (uniform, non-varying) lens defocus blur in Python/OpenCV by transforming both the image and filter to the Fourier (frequency) domain.

• Read the input
• Take dft of input to transform to Fourier domain
• Draw a white filled circle on a black background the size of the input as a mask (filter kernel). This is the defocus kernel in the spatial domain, i.e. a circular rect function.
• Blur the circle slightly to anti-alias the edge
• Roll the mask so that the center is at the origin (top left corner) and normalize so that the sum of values = 1
• Take dft of mask to transform to Fourier domain where its amplitude profile is a jinx function.
• Multiply the two dft images to apply the blur
• Take the idft of the product to transform back to spatial domain
• Get the magnitude of the real and imaginary components of the product, clip and convert to uint8 as the result
• Save the result

Input:

import numpy as np
import cv2

# read input and convert to grayscale
img = cv2.imread('lena_512_gray.png', cv2.IMREAD_GRAYSCALE)

# do dft saving as complex output
dft_img = np.fft.fft2(img, axes=(0,1))

# create circle mask
radius = 32
mask = np.zeros_like(img)
cy = mask.shape[0] // 2
cx = mask.shape[1] // 2
cv2.circle(mask, (cx,cy), radius, 255, -1)[0]

# blur the mask slightly to antialias
mask = cv2.GaussianBlur(mask, (3,3), 0)

# roll the mask so that center is at origin and normalize to sum=1
mask_roll = np.roll(mask, (256,256), axis=(0,1))
mask_norm = mask_roll / mask_roll.sum() 

# take dft of mask
dft_mask_norm = np.fft.fft2(mask_norm, axes=(0,1))

# apply dft_mask to dft_img
dft_shift_product = np.multiply(dft_img, dft_mask_norm)

# do idft saving as complex output
img_filtered = np.fft.ifft2(dft_shift_product, axes=(0,1))

# combine complex real and imaginary components to form (the magnitude for) the original image again
img_filtered = np.abs(img_filtered).clip(0,255).astype(np.uint8)

cv2.imshow("ORIGINAL", img)
cv2.imshow("MASK", mask)
cv2.imshow("FILTERED DFT/IFT ROUND TRIP", img_filtered)
cv2.waitKey(0)
cv2.destroyAllWindows()

# write result to disk
cv2.imwrite("lena_512_gray_mask.png", mask)
cv2.imwrite("lena_dft_numpy_lowpass_filtered_rad32.jpg", img_filtered)

Mask - Filter Kernel In Spatial Domain:

Result for Circle Radius=4:

Result for Circle Radius=8:

Result for Circle Radius=16:

Result for Circle Radius=32

:

ADDITION

Using OpenCV for the dft, etc rather than Numpy, the above becomes:

import numpy as np
import cv2

# read input and convert to grayscale
img = cv2.imread('lena_512_gray.png', cv2.IMREAD_GRAYSCALE)

# do dft saving as complex output
dft_img = cv2.dft(np.float32(img), flags = cv2.DFT_COMPLEX_OUTPUT)

# create circle mask
radius = 32
mask = np.zeros_like(img)
cy = mask.shape[0] // 2
cx = mask.shape[1] // 2
cv2.circle(mask, (cx,cy), radius, 255, -1)[0]

# blur the mask slightly to antialias
mask = cv2.GaussianBlur(mask, (3,3), 0)

# roll the mask so that center is at origin and normalize to sum=1
mask_roll = np.roll(mask, (256,256), axis=(0,1))
mask_norm = mask_roll / mask_roll.sum() 

# take dft of mask
dft_mask_norm = cv2.dft(np.float32(mask_norm), flags = cv2.DFT_COMPLEX_OUTPUT)

# apply dft_mask to dft_img
dft_product = cv2.mulSpectrums(dft_img, dft_mask_norm, 0)

# do idft saving as complex output, then clip and convert to uint8
img_filtered = cv2.idft(dft_product, flags=cv2.DFT_SCALE+cv2.DFT_REAL_OUTPUT)
img_filtered = img_filtered.clip(0,255).astype(np.uint8)

cv2.imshow("ORIGINAL", img)
cv2.imshow("MASK", mask)
cv2.imshow("FILTERED DFT/IFT ROUND TRIP", img_filtered)
cv2.waitKey(0)
cv2.destroyAllWindows()

# write result to disk
cv2.imwrite("lena_512_gray_mask.png", mask)
cv2.imwrite("lena_dft_opencv_defocus_rad32.jpg", img_filtered)