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

Issue

I want to emulate the blur of a cheap camera lens (like Holga).
Blur is very weak close to the photo center.
And it's getting more decisive close to corners.

I wrote the code and it works in general.

Input image:

Result image:

.

But I feel that it could be done better and faster.
I've found a similar question but it still has no answer.

How to improve an algorithm speed and avoid iteration over pixels?

UPDATE:
It's not the same as standard Gaussian or 2D filter blur with constant kernel size.

import cv2
import numpy as np
import requests
from tqdm import tqdm
import warnings
warnings.filterwarnings("ignore")

def blur(img=None, blur_radius=None, test=False):
    # test image loading
    if img is None:
        test=True
        print('test mode ON')
        print('loading image...')
        url = r'http://www.lenna.org/lena_std.tif'
        resp = requests.get(url, stream=True).raw
        img = np.asarray(bytearray(resp.read()), dtype="uint8")
        img = cv2.imdecode(img, cv2.IMREAD_COLOR)
        cv2.imwrite('img_input.png', img)
        print('image loaded')

    # channels splitting
    img_lab = cv2.cvtColor(img, cv2.COLOR_BGR2LAB)
    l, a, b = cv2.split(img_lab)
    if test:
        cv2.imwrite('l_channel.png', l)
        print('l channel saved')

    # make blur map 
    height, width = l.shape[:2]
    center = np.array([height/2, width/2])
    diag = ((height / 2) ** 2 + (width / 2) ** 2) ** 0.5
    blur_map = np.linalg.norm(
        np.indices(img.shape[:2]) - center[:,None,None] + 0.5,
        axis = 0
    )

    if blur_radius is None:
        blur_radius = int(max(height, width) * 0.03)

    blur_map = blur_map / diag 
    blur_map = blur_map * blur_radius
    if test:
        blur_map_norm = cv2.normalize(blur_map, None, 0, 255, cv2.NORM_MINMAX, cv2.CV_32F)
        cv2.imwrite('blur_map.png', blur_map_norm)
        print('blur map saved')

    # very inefficient blur algorithm!!!
    l_blur = np.copy(l)
    for x in tqdm(range(width)):
        for y in range(height):
            kernel_size = int(blur_map[y, x])
       
            if kernel_size == 0:
                l_blur[y, x] = l[y, x]
                continue
            
            kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (kernel_size, kernel_size))
            cut = l[
                max(0, y - kernel_size):min(height, y + kernel_size),
                max(0, x - kernel_size):min(width, x + kernel_size)
            ]

            if cut.shape == kernel.shape:
                cut = (cut * kernel).mean()
            else:
                cut = cut.mean()

            l_blur[y, x] = cut
    if test: cv2.imwrite('l_blur.png', l_blur); print('l_blur saved')
    if test: print('done')
    return l_blur

blur() 

Solution

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])

Answered By - Cris Luengo