From 95796b53fc5093e8db9035dd5c56b0b5a1acb874 Mon Sep 17 00:00:00 2001 From: Benjamin Loison <12752145+Benjamin-Loison@users.noreply.github.com> Date: Wed, 15 May 2024 17:41:19 +0200 Subject: [PATCH] #69: Remove artifacts from low pass denoiser See [Robust_image_source_identification_on_modern_smartphones/issues/69#issuecomment-1892](https://gitea.lemnoslife.com/Benjamin_Loison/Robust_image_source_identification_on_modern_smartphones/issues/69#issuecomment-1892). --- datasets/raise/utils.py | 70 ++++++++++++++++++++++++++++++++++------- 1 file changed, 58 insertions(+), 12 deletions(-) diff --git a/datasets/raise/utils.py b/datasets/raise/utils.py index 1250784..a8f8110 100644 --- a/datasets/raise/utils.py +++ b/datasets/raise/utils.py @@ -3,12 +3,15 @@ import skimage.restoration import numpy as np import rawpy from tqdm import tqdm -from PIL import Image +from PIL import Image, ImageDraw from skimage import img_as_float from datetime import datetime import builtins as __builtin__ from scipy.ndimage import gaussian_filter import os +from skimage.restoration import estimate_sigma +from scipy import fftpack +import imageio class Color(Enum): RED = auto() @@ -24,16 +27,23 @@ class Denoiser(Enum): BILATERAL = auto() TV_CHAMBOLLE = auto() MEAN = auto() + NL_MEANS = auto() + LOW_PASS = auto() def __str__(self): return self.name.lower() -# Among: -# - `wavelet` -# - `bilateral` -# - `tv_chambolle` +# Source: https://stackoverflow.com/a/43346070 +def getGaussianKernel(size, sigma): + axis = np.linspace(-(size - 1) / 2., (size - 1) / 2., size) + gauss = np.exp(-0.5 * np.square(axis) / np.square(sigma)) + kernel = np.outer(gauss, gauss) + return kernel / np.sum(kernel) + +# `denoiser` can be all values of `Denoiser` except `MEAN`. def denoise(imageNpArray, denoiser): - skImageRestorationDenoise = getattr(skimage.restoration, f'denoise_{denoiser}') + if denoiser != Denoiser.LOW_PASS: + skImageRestorationDenoise = getattr(skimage.restoration, f'denoise_{denoiser}') match denoiser: case Denoiser.WAVELET: @@ -42,6 +52,42 @@ def denoise(imageNpArray, denoiser): imageDenoisedNpArray = skImageRestorationDenoise(imageNpArray, sigma_color=0.05, sigma_spatial=15) case Denoiser.TV_CHAMBOLLE: imageDenoisedNpArray = skImageRestorationDenoise(imageNpArray, weight=0.2) + case Denoiser.NL_MEANS: + sigmaEst = np.mean(estimate_sigma(imageNpArray, channel_axis = -1)) + imageDenoisedNpArray = skImageRestorationDenoise(imageNpArray, patch_size = 7, patch_distance = 4, h = 0.8 * sigmaEst, fast_mode = True) + case Denoiser.LOW_PASS: + # Based on [the Stack Overflow answer 54662336](https://stackoverflow.com/a/54662336). + # fft of image + fft1 = fftpack.fftshift(fftpack.fft2(imageNpArray)) + + GAUSSIAN_SIGMA = 50 + #low_pass_np = np.zeros(image1_np.shape, dtype = np.uint8) + + low_pass_np_size = min(image1_np.shape) + low_pass_np = getGaussianKernel(low_pass_np_size, GAUSSIAN_SIGMA) + + #rectangle_low_pass_np = np.zeros(fft1.shape) + #rectangle_low_pass_np[] + + requiredH, requiredW = image1_np.shape + vert_pad_size = (requiredH - low_pass_np_size) // 2 + hor_pad_size = (requiredW - low_pass_np_size) // 2 + rectangle_low_pass_np = np.pad(low_pass_np, [(vert_pad_size, vert_pad_size), (hor_pad_size, hor_pad_size)]) + + # Unsure that doing such *normalization* is appropriate. + rectangle_low_pass_np /= rectangle_low_pass_np.max() + + #plt.imshow(rectangle_low_pass_np) + #plt.show() + + # multiply both the images + filtered = np.multiply(fft1, lowPassNp) + + # inverse fft + ifft2 = np.real(fftpack.ifft2(fftpack.ifftshift(filtered))) + ifft2 = np.maximum(0, np.minimum(ifft2, 255)) + + imageDenoisedNpArray = ifft2.astype(np.uint8) return imageDenoisedNpArray class iterativeMean: @@ -91,7 +137,7 @@ def isARawImage(imageFilePath): def getColorChannel(imageFilePath, color): if isARawImage(imageFilePath): numpyFilePath = f'{imageFilePath}.{color}.npy' - if os.path.isfile(numpyFilePath): + if False and os.path.isfile(numpyFilePath): imageNpArray = np.load(numpyFilePath) else: with rawpy.imread(imageFilePath) as raw: @@ -144,7 +190,8 @@ def updateExtremes(imageNpArray, minColor, maxColor): return minColor, maxColor def print(*toPrint): - __builtin__.print(datetime.now(), *toPrint) + #__builtin__.print(datetime.now(), *toPrint) + tqdm.write(f'{datetime.now()} {str(toPrint)}') def getColorMeans(imagesFileNames, colors, singleColorChannelCropResolution = None): colorMeans = {} @@ -159,19 +206,18 @@ def getColorMeans(imagesFileNames, colors, singleColorChannelCropResolution = No colorMeans[color] = colorIterativeMean.mean return colorMeans -def getImageNpArray(imageFilePath, computeExtremes, color, denoiser): - global minColor, maxColor +def getImageNpArray(imageFilePath, computeExtremes, color, denoiser, minColor, maxColor): imageNpArray = getColorChannel(imageFilePath, color) if computeExtremes: minColor, maxColor = updateExtremes(imageNpArray, minColor, maxColor) - return + return None, minColor, maxColor if isARawImage(imageFilePath) and denoiser != Denoiser.MEAN: imageNpArray = rescaleRawImageForDenoiser(imageNpArray, minColor, maxColor) # Pay attention to range of values expected by the denoiser. # Indeed if provide the thousands valued raw image, then the denoiser only returns values between 0 and 1 and making the difference between both look pointless. - return imageNpArray + return imageNpArray, minColor, maxColor def getImageCrop(image, cropResolution): imageCrop = image[:cropResolution[0], :cropResolution[1]]