96 lines
3.6 KiB
Python
96 lines
3.6 KiB
Python
import numpy as np
|
|
from astropy.convolution import Gaussian2DKernel, interpolate_replace_nans
|
|
from PIL import Image
|
|
import os
|
|
import sys
|
|
from scipy import fftpack
|
|
|
|
sys.path.insert(0, '../')
|
|
|
|
from utils import Color, mergeSingleColorChannelImagesAccordingToBayerFilter
|
|
import matplotlib.pyplot as plt
|
|
|
|
PREFIX = 'means/mean_rafael_230424_mean_'
|
|
X_STDDEV = 2
|
|
|
|
def getImageByColor(color):
|
|
filePath = PREFIX + f'{color}.npy'
|
|
image = np.load(filePath)
|
|
return image
|
|
|
|
singleColorChannelImages = {color: getImageByColor(color) for color in Color}
|
|
multipleColorsImage = mergeSingleColorChannelImagesAccordingToBayerFilter(singleColorChannelImages)
|
|
image = multipleColorsImage
|
|
|
|
fft1 = fftpack.fftshift(fftpack.fft2(image))
|
|
|
|
def removePeriodicPatterns(fft1Part):
|
|
# This example is intended to demonstrate how astropy.convolve and
|
|
# scipy.convolve handle missing data, so we start by setting the brightest
|
|
# pixels to NaN to simulate a "saturated" data set
|
|
# Unclear why the Gaussian kernel usage does not work as expected, even when just consider the *real* part.
|
|
# See [Benjamin_Loison/astropy/issues/2](https://codeberg.org/Benjamin_Loison/astropy/issues/2).
|
|
fixedImagePart = fft1Part.copy()
|
|
height, width = fft1Part.shape
|
|
for x in range(width):
|
|
#fft1Part[height // 2, x] = np.nan
|
|
fixedImagePart[height // 2, x] = np.mean([fft1Part[height // 2 - 1, x], fft1Part[height // 2 + 1, x]])
|
|
|
|
middleX = width // 2
|
|
RANGE = 1
|
|
for verticalLineX in [width // 4, width // 2, round(3 * width / 4)]:
|
|
for y in range(height):
|
|
for x in range(verticalLineX - RANGE, verticalLineX + RANGE + 1):
|
|
#fft1Part[y, x] = np.nan
|
|
fixedImagePart[y, x] = np.mean([fft1Part[y, x - RANGE - 1], fft1Part[y, x + RANGE + 1]])
|
|
|
|
# We smooth with a Gaussian kernel
|
|
#kernel = Gaussian2DKernel(x_stddev = X_STDDEV)
|
|
|
|
# create a "fixed" image with NaNs replaced by interpolated values
|
|
# See [Benjamin_Loison/astropy/issues/1](https://codeberg.org/Benjamin_Loison/astropy/issues/1).
|
|
#fixedImagePart = interpolate_replace_nans(fft1Part, kernel)
|
|
return fixedImagePart
|
|
|
|
realFixedImage = removePeriodicPatterns(np.real(fft1))
|
|
imaginaryFixedImage = removePeriodicPatterns(np.imag(fft1))
|
|
fixedImage = realFixedImage + 1j * imaginaryFixedImage
|
|
|
|
figure, axes = plt.subplots(1, 2, sharex = True, sharey = True)
|
|
|
|
plt.suptitle('Attenuating FFT significant lines')
|
|
axes[0].set_title('Original FFT')
|
|
firstImage = np.log10(1 + abs(fft1))
|
|
secondImage = np.log10(1 + abs(fixedImage))
|
|
images = [firstImage, secondImage]
|
|
vMin = np.min(images)
|
|
vMax = np.max(images)
|
|
axes[0].imshow(firstImage, vmin = vMin, vmax = vMax)
|
|
axes[1].set_title('FFT with significant lines attenuated')
|
|
axes[1].imshow(secondImage, vmin = vMin, vmax = vMax)
|
|
plt.tight_layout()
|
|
plt.show()
|
|
|
|
figure, axes = plt.subplots(1, 3, sharex = True, sharey = True)
|
|
|
|
def inverseFft(fft):
|
|
ifft2 = np.real(fftpack.ifft2(fftpack.ifftshift(fft)))
|
|
return ifft2
|
|
|
|
invFft1 = inverseFft(fft1)
|
|
invFixedImage = inverseFft(fixedImage)
|
|
images = [invFft1, invFixedImage]
|
|
vMin = np.min(images)
|
|
vMax = np.max(images)
|
|
plt.suptitle('Rafael 23/04/24 PRNU mean denoiser with periodic patterns attenuated')
|
|
axes[0].set_title('Original PRNU')
|
|
axes[0].imshow(invFft1, vmin = vMin, vmax = vMax)
|
|
axes[1].set_title('PRNU with periodic patterns attenuated')
|
|
axes[1].imshow(invFixedImage, vmin = vMin, vmax = vMax)
|
|
axes[2].set_title('Difference between both left images')
|
|
differenceBetweenBothImages = invFft1 - invFixedImage
|
|
# `np.log10(1 + abs(differenceBetweenBothImages))` does not seem more interesting.
|
|
axes[2].imshow(differenceBetweenBothImages)
|
|
plt.tight_layout()
|
|
plt.show()
|