diff --git a/algorithms/context_adaptative_interpolator/context_adaptive_interpolator.py b/algorithms/context_adaptative_interpolator/context_adaptive_interpolator.py
new file mode 100644
index 0000000..da680b1
--- /dev/null
+++ b/algorithms/context_adaptative_interpolator/context_adaptive_interpolator.py
@@ -0,0 +1,53 @@
+# 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 wiener_filter import wienerFilter
+
+# Assume greyscale PIL image passed.
+# What about other color channels? See #11.
+def contextAdaptiveInterpolator(I, IImage):
+    rImage = Image.new('L', (IImage.size[0] - 2, IImage.size[1] - 2))
+    r = rImage.load()
+
+    # This threshold is debatable. See #13.
+    THRESHOLD = 20
+
+    print('before for loops')
+    # 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)
+    print('after for loops')
+
+    # 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
+
+    print('before wiener filter')
+    return wienerFilter(r, rImage, Q, sigma_0)
\ No newline at end of file
diff --git a/algorithms/context_adaptative_interpolator/test.py b/algorithms/context_adaptative_interpolator/test.py
new file mode 100644
index 0000000..ddf1c0d
--- /dev/null
+++ b/algorithms/context_adaptative_interpolator/test.py
@@ -0,0 +1,9 @@
+from PIL import Image
+from context_adaptive_interpolator import contextAdaptiveInterpolator
+
+MODE = 'L'
+
+IImage = Image.open('9f04e2005fddb9d5512e2f42a3b826b019755717.jpg').convert(MODE)
+I = IImage.load()
+
+contextAdaptiveInterpolator(I, IImage).show()
\ No newline at end of file
diff --git a/algorithms/context_adaptative_interpolator/wiener_filter.py b/algorithms/context_adaptative_interpolator/wiener_filter.py
new file mode 100644
index 0000000..33d3577
--- /dev/null
+++ b/algorithms/context_adaptative_interpolator/wiener_filter.py
@@ -0,0 +1,38 @@
+from PIL import Image
+from tqdm import tqdm
+
+def wienerFilter(r, rImage, Q, sigma_0):
+    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]))
+
+    print('wiener filter start for loops')
+    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))
+    print('wiener filter end for loops')
+
+    return h_wImage.rotate(-90)
\ No newline at end of file