84 lines
2.9 KiB
Python
84 lines
2.9 KiB
Python
# Based on https://web.archive.org/web/20231116015653/http://nrl.northumbria.ac.uk/id/eprint/29339/1/Paper_accepted.pdf IV. B..
|
|
|
|
from PIL import Image
|
|
from statistics import mean, median
|
|
from tqdm import tqdm
|
|
|
|
# What about other color channels? See #11.
|
|
MODE = 'L'
|
|
|
|
IImage = Image.open('9f04e2005fddb9d5512e2f42a3b826b019755717.jpg').convert(MODE)
|
|
I = IImage.load()
|
|
|
|
rImage = Image.new(MODE, (IImage.size[0] - 2, IImage.size[1] - 2))
|
|
r = rImage.load()
|
|
|
|
# This threshold is debatable. See #13.
|
|
THRESHOLD = 20
|
|
|
|
# Equation (10)
|
|
# Accelerate computation. See #15.
|
|
for m in range(1, IImage.size[0] - 1):
|
|
for n in range(1, IImage.size[1] - 1):
|
|
e = I[m, n + 1]
|
|
se = I[m + 1, n + 1]
|
|
s = I[m + 1, n]
|
|
sw = I[m + 1, n - 1]
|
|
w = I[m, n - 1]
|
|
nw = I[m - 1, n - 1]
|
|
no = I[m - 1, n]
|
|
ne = I[m - 1, n + 1]
|
|
A = [e, se, s, sw, w, nw, no, ne]
|
|
if max(A) - min(A) <= THRESHOLD:
|
|
newPixel = I[m, n] - mean(A)
|
|
elif abs(e - w) - abs(no - s) > THRESHOLD:
|
|
newPixel = I[m, n] - (s + no) / 2
|
|
elif abs(s - no) - abs(e - w) > THRESHOLD:
|
|
newPixel = I[m, n] - (e + w) / 2
|
|
elif abs(sw - ne) - abs(se - nw) > THRESHOLD:
|
|
newPixel = I[m, n] - (se + nw) / 2
|
|
elif abs(se - nw) - abs(sw - ne) > THRESHOLD:
|
|
newPixel = I[m, n] - (sw + ne) / 2
|
|
else:
|
|
newPixel = I[m, n] - median(A)
|
|
r[m - 1, n - 1] = round(newPixel)
|
|
|
|
# Why need to rotate the image? See #14.
|
|
rImage.rotate(-90).show()
|
|
|
|
Q = 3
|
|
# $\sigma_0^2$ is the noise variance.
|
|
sigma_0 = 9 ** 0.5
|
|
|
|
h_wImage = Image.new(MODE, (rImage.size[0], rImage.size[1]))
|
|
h_wImagePixels = h_wImage.load()
|
|
|
|
# Wiener filter.
|
|
def h_w(hImage, h, i, j):
|
|
# Equation (7)
|
|
return h[i, j] * sigma(hImage, h, i, j) / (sigma(hImage, h, i, j) + sigma_0 ** 2)
|
|
|
|
# Minimum of the considered variances.
|
|
def sigma(hImage, h, i, j):
|
|
# Equation (9)
|
|
return sigma_q(hImage, h, i, j, Q)
|
|
|
|
def getPixelIndexesAround(i, numberOfPixelsInEachDirection):
|
|
return range(i - numberOfPixelsInEachDirection, i + numberOfPixelsInEachDirection + 1)
|
|
|
|
# Expand image with border pixels.
|
|
def getPixelWithinImage(z, upperBound):
|
|
return max(min(z, upperBound - 1), 0)
|
|
|
|
# Local variance obtained by Maximum A Posteriori (MAP).
|
|
def sigma_q(hImage, h, i, j, q):
|
|
# Equation (8)
|
|
numberOfPixelsInEachDirection = (q - 1) // 2
|
|
B_q = [(x, z) for x in getPixelIndexesAround(i, numberOfPixelsInEachDirection) for z in getPixelIndexesAround(j, numberOfPixelsInEachDirection)]
|
|
return max(0, (1 / q ** 2) * sum([h[getPixelWithinImage(x, hImage.size[0]), getPixelWithinImage(z, hImage.size[1])] ** 2 - sigma_0 ** 2 for (x, z) in B_q]))
|
|
|
|
for i in tqdm(range(rImage.size[0])):
|
|
for j in range(rImage.size[1]):
|
|
h_wImagePixels[i, j] = round(h_w(rImage, r, i, j))
|
|
|
|
h_wImage.rotate(-90).show() |