Trying to blur highest variance point of an image but some conversion problem exist in code

Question

I am trying to blur of highest variance point from the image. I wrote code below. 1st part finds the variance of the image. I checked the resultant variance of an image and it is correct. (I used Lena's image) In 2nd part, I find the highest variance coordinates and send to this Function which finds gaussian blur. When I execute this code, it throws "C:\Tmp\blur_highest_variance.py", line 66, in sigma=15) numpy.core._exceptions.UFuncTypeError: Cannot cast ufunc 'subtract' output from dtype('float64') to dtype('uint8') with casting rule 'same_kind' I tried a few conversions between types but no avail. Can you show me some direction?

import numpy as np
import matplotlib.pylab as plt    
from skimage import measure
from PIL import Image, ImageChops
import math
import cv2
from skimage.morphology import rectangle
import skimage.filters as filters

######################Calculate Variance#######################
# Variance = mean of square of image - square of mean of image
# See # see https://en.wikipedia.org/wiki/Variance

# read the image
# convert to 16-bits grayscale since mean filter below is limited 
# to single channel 8 or 16-bits, not float
# and variance will be larger than 8-bit range
img = cv2.imread(r".\lena_std512_512.jpg", cv2.IMREAD_GRAYSCALE).astype(np.uint16)

# compute square of image
img_sq = cv2.multiply(img, img)

# compute local mean in 5x5 rectangular region of each image
# note: python will give warning about slower performance when processing 16-bit images
region = rectangle(10,10)
mean_img = filters.rank.mean(img, selem=region)
mean_img_sq = filters.rank.mean(img_sq, selem=region)

# compute square of local mean of img
sq_mean_img = cv2.multiply(mean_img, mean_img)

# compute variance using float versions of images
var = cv2.add(mean_img_sq.astype(np.float32), -sq_mean_img.astype(np.float32))

# compute standard deviation and convert to 8-bit format
std = cv2.sqrt(var).clip(0,255).astype(np.uint8)

# multiply by 2 to make brighter as an example
cv2.imwrite('lena_std_local_variance.jpg',std)    

#################Gaussian Blur Function###############
def gaussian_mask(x, y, shape, amp=1, sigma=15):
    """
    Returns an array of shape, with values based on

    amp * exp(-((i-x)**2 +(j-y)**2) / (2 * sigma ** 2))

    :param x: float
    :param y: float
    :param shape: tuple
    :param amp: float
    :param sigma: float
    :return: array
    """
    xv, yv = np.meshgrid(np.arange(shape[1]), np.arange(shape[0]))
    g = amp * np.exp(-((xv - x) ** 2 + (yv - y) ** 2) / (2 * sigma ** 2))
    return g

#################Find Gaussian Blur and Subtract###############

y, x = np.unravel_index(np.argmax(std), std.shape)       

std -= gaussian_mask(x, y,
                    shape=std.shape[:2],
                    amp=1,
                    sigma=15)

cv2.imwrite(r'.\gaussian\lena_std_local_variance.jpg',std)

Answer 1

The error message tells us the line and the reason for the error:

Traceback (most recent call last):
File "C:\Tmp\blur_highest_variance.py", line 66, in sigma=15)
numpy.core._exceptions.UFuncTypeError: Cannot cast ufunc 'subtract' output from dtype('float64') to dtype('uint8') with casting rule 'same_kind'

It is more simple to debug the code using intermediate variables:
For example, use an intermediate named gmask:

gmask = gaussian_mask(x, y,
                      shape=std.shape[:2],
                      amp=1,
                      sigma=15)

print('gmask.dtype = ' + str(gmask.dtype))  # gmask.dtype = float64
print('std.dtype = ' + str(std.dtype))  # std.dtype = uint8

We don't really need to print the dtype, we may use the debugger, but printing demonstrates the reason for the error.
We can't subtract float64 array from uint8 using -= operator!

Using std = std - gmask, is not an error, but the type of the result is float64.

---

Suggested solution:

Cast gmask to uint8 and use cv2.subtract:

gmask = gaussian_mask(x, y,
                      shape=std.shape[:2],
                      amp=1,
                      sigma=15)

std = cv2.subtract(std, gmask.clip(0, 255).astype(np.uint8))

Using cv2.subtract is a safe way for subtracting two uint8 matrices, because it clips the result to [0, 255] (includes overflow protection).

---

Result (lena_std_local_variance.jpg):