First PRNU distance test

algorithms/context_adaptive_interpolator/wiener_filter.py

@@ -0,0 +1,41 @@
+from PIL import Image
+from tqdm import tqdm
+
+def wienerFilter(r, rImage, Q, sigma_0, showProgress):
+    h_wImage = Image.new('L', (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]))
+
+    if showProgress:
+        print('wiener filter start for loops')
+    rImageSize0Range = range(rImage.size[0])
+    for i in tqdm(rImageSize0Range) if showProgress else rImageSize0Range:
+        for j in range(rImage.size[1]):
+            h_wImagePixels[i, j] = round(h_w(rImage, r, i, j))
+    if showProgress:
+        print('wiener filter end for loops')
+
+    return h_wImage.rotate(-90)