diff --git a/src/methods/ELAMRANI_Mouna/main.ipynb b/src/methods/ELAMRANI_Mouna/main.ipynb
old mode 100755
new mode 100644
diff --git a/src/methods/digeronimo/reconstruct.py b/src/methods/digeronimo/reconstruct.py
old mode 100755
new mode 100644
diff --git a/src/methods/marty/Report.pdf b/src/methods/marty/Report.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..24c0a7f93fc289f9dd00312886c774fcbb5c5bde
Binary files /dev/null and b/src/methods/marty/Report.pdf differ
diff --git a/src/methods/marty/high_quality_interpolation.py b/src/methods/marty/high_quality_interpolation.py
new file mode 100644
index 0000000000000000000000000000000000000000..64f52087ba29ec525f3f63e859b761987bf42572
--- /dev/null
+++ b/src/methods/marty/high_quality_interpolation.py
@@ -0,0 +1,81 @@
+import numpy as np
+
+
+def is_green(i, j):
+    return (not i%2 and not j%2) or (i%2 and j%2)
+
+def is_red(i, j):
+    return not i%2 and j%2
+
+def is_blue(i, j):
+    return i%2 and not j%2
+
+
+g_ker = np.array([
+    [0, 0, -1, 0, 0],
+    [0, 0, 2, 0, 0],
+    [-1, 2, 4, 2, -1],
+    [0, 0, 2, 0, 0],
+    [0, 0, -1, 0, 0]
+]) / 8
+
+rgrrow_ker = np.array([
+    [0, 0, 0.5, 0, 0],
+    [0, -1, 0, -1, 0],
+    [-1, 4, 5, 4, -1],
+    [0, -1, 0, -1, 0],
+    [0, 0, 0.5, 0, 0]
+]) / 8
+
+rgrcol_ker = np.array([
+    [0, 0, -1, 0, 0],
+    [0, -1, 4, -1, 0],
+    [0.5, 0, 5, 0, 0.5],
+    [0, -1, 4, -1, 0],
+    [0, 0, -1, 0, 0]
+]) / 8
+
+rb_ker = np.array([
+    [0, 0, -1.5, 0, 0],
+    [0, 2, 0, 2, 0],
+    [-1.5, 0, 6, 0, -1.5],
+    [0, 2, 0, 2, 0],
+    [0, 0, -1.5, 0, 0]
+]) / 8
+
+
+
+def high_quality_interpolation(image):
+
+    res = np.empty((image.shape[0], image.shape[1], 3))
+    
+    estimated_green = np.zeros_like(image)
+    estimated_red = np.zeros_like(image)
+    estimated_blue = np.zeros_like(image)
+
+    for i in range(2, image.shape[0]-2):
+        for j in range(2, image.shape[1]-2):
+            if (is_red(i, j)):
+                estimated_red[i, j] = image[i, j]
+                estimated_green[i, j] = (image[i-2:i+3, j-2:j+3] * g_ker).sum()
+                estimated_blue[i, j] = (image[i-2:i+3, j-2:j+3] * rb_ker).sum()
+            elif (is_blue(i, j)):
+                estimated_red[i, j] = (image[i-2:i+3, j-2:j+3] * rb_ker).sum()
+                estimated_green[i, j] = (image[i-2:i+3, j-2:j+3] * g_ker).sum()
+                estimated_blue[i, j] = image[i, j]
+            else:
+                estimated_green[i, j] = image[i, j]
+                if (is_red(i-1, j)):
+                    estimated_red[i, j] = (image[i-2:i+3, j-2:j+3] * rgrcol_ker).sum()
+                    estimated_blue[i, j] = (image[i-2:i+3, j-2:j+3] * rgrrow_ker).sum()
+                else:
+                    estimated_red[i, j] = (image[i-2:i+3, j-2:j+3] * rgrrow_ker).sum()
+                    estimated_blue[i, j] = (image[i-2:i+3, j-2:j+3] * rgrcol_ker).sum()
+
+
+    res[:, :, 0] = estimated_red.clip(0, 1)
+    res[:, :, 1] = estimated_green.clip(0, 1)
+    res[:, :, 2] = estimated_blue.clip(0, 1)
+    
+    return res
+
diff --git a/src/methods/marty/mdwi.py b/src/methods/marty/mdwi.py
new file mode 100644
index 0000000000000000000000000000000000000000..800b8fa967c203493d2fb76cfd24ee23a9bfcb63
--- /dev/null
+++ b/src/methods/marty/mdwi.py
@@ -0,0 +1,129 @@
+import numpy as np
+
+
+def is_green(i, j):
+    return (not i%2 and not j%2) or (i%2 and j%2)
+
+def is_red(i, j):
+    return not i%2 and j%2
+
+def is_blue(i, j):
+    return i%2 and not j%2
+
+
+
+def MDWI(image):
+    
+    res = np.empty((image.shape[0], image.shape[1], 3))
+    estimated_green = np.zeros_like(image)
+    estimated_red = np.zeros_like(image)
+    estimated_blue = np.zeros_like(image)
+
+    EPSILON = 10e-4
+
+    # Estimation of green plan first
+
+    pre_estimation = {}
+    diag_gradient_factor = {}
+    weight = {}
+    H = [-8/256, 23/256, -48/256, 161/256, 161/256, -48/256, 23/256, -8/256]
+    coord_NE = [(-4, -3), (-3, -2), (-2, -1), (-1, 0), (0, 1), (1, 2), (2, 3), (3, 4)]
+    coord_SE = [(-3, -4), (-2, -3), (-1, -2), (0, -1), (1, 0), (2, 1), (3, 2), (4, 3)]
+    coord_NW = [(3, -4), (2, -3), (1, -2), (0, -1), (-1, 0), (-2, 1), (-3, 2), (-4, 3)]
+    coord_SW = [(4, -3), (3, -2), (2, -1), (1, 0), (0, 1), (-1, 2), (-2, 3), (-3, 4)]
+
+    for i in range(3, image.shape[0] - 4):
+        for j in range(3, image.shape[1] - 4):
+            
+
+            if (not is_green(i, j)):       
+                pre_estimation["N"] = image[i-1, j] + (image[i, j] - image[i-2, j]) / 2   
+                pre_estimation["S"] = image[i+1, j] + (image[i, j] - image[i+2, j]) / 2
+                pre_estimation["W"] = image[i, j-1] + (image[i, j] - image[i, j-2]) / 2
+                pre_estimation["E"] = image[i, j+1] + (image[i, j] - image[i, j+2]) / 2
+
+                diag_gradient_factor["N"] = abs(image[i-2, j-1] - image[i, j-1]) + abs(image[i-3, j] - image[i-1, j]) + abs(image[i-2, j+1] - image[i, j+1]) + abs(image[i-3, j-1] - image[i-1, j-1]) + abs(image[i-3, j+1] - image[i-1, j+1]) + abs(image[i-2, j] - image[i, j]) + EPSILON
+                diag_gradient_factor["S"] = abs(image[i+2, j-1] - image[i, j-1]) + abs(image[i+3, j] - image[i+1, j]) + abs(image[i+2, j+1] - image[i, j+1]) + abs(image[i+3, j-1] - image[i+1, j-1]) + abs(image[i+3, j+1] - image[i+1, j+1]) + abs(image[i+2, j] - image[i, j]) + EPSILON
+                diag_gradient_factor["W"] = abs(image[i-1, j-2] - image[i-1, j]) + abs(image[i, j-3] - image[i, j-1]) + abs(image[i+1, j-2] - image[i+1, j]) + abs(image[i-1, j-3] - image[i-1, j-1]) + abs(image[i+1, j-3] - image[i+1, j-1]) + abs(image[i, j-2] - image[i, j]) + EPSILON
+                diag_gradient_factor["E"] = abs(image[i-1, j+2] - image[i-1, j]) + abs(image[i, j+3] - image[i, j+1]) + abs(image[i+1, j+2] - image[i+1, j]) + abs(image[i-1, j+3] - image[i-1, j+1]) + abs(image[i+1, j+3] - image[i+1, j+1]) + abs(image[i, j+2] - image[i, j]) + EPSILON
+
+
+                pre_estimation["NW"] = 0
+                pre_estimation["SW"] = 0
+                pre_estimation["NE"] = 0
+                pre_estimation["SE"] = 0
+                for k in range(8):
+                    pre_estimation["NW"] += image[i + coord_NW[k][0], j + coord_NW[k][1]] * H[k]
+                    pre_estimation["SW"] += image[i + coord_SW[k][0], j + coord_SW[k][1]] * H[k]
+                    pre_estimation["NE"] += image[i + coord_NE[k][0], j + coord_NE[k][1]] * H[k]
+                    pre_estimation["SE"] += image[i + coord_SE[k][0], j + coord_SE[k][1]] * H[k]
+                
+                diag_gradient_factor["NW"] = abs(image[i-2, j-1] - image[i-1, j]) + abs(image[i-1, j] - image[i, j+1]) + abs(image[i-1, j-2] - image[i, j-1]) + abs(image[i, j-1] - image[i+1, j]) + abs(image[i-1, j-1] - image[i+1, j+1]) + abs(image[i-2, j-2] - image[i, j]) + EPSILON
+                diag_gradient_factor["SW"] = abs(image[i-1, j] - image[i, j-1]) + abs(image[i, j+1] - image[i+1, j]) + abs(image[i-1, j+1] - image[i+1, j-1]) + abs(image[i, j-1] - image[i+1, j-2]) + abs(image[i+1, j] - image[i+2, j-1]) + abs(image[i, j] - image[i+2, j-2]) + EPSILON
+                diag_gradient_factor["NE"] = abs(image[i-1, j] - image[i, j-1]) + abs(image[i, j+1] - image[i+1, j]) + abs(image[i-1, j+1] - image[i+1, j-1]) + abs(image[i-1, j] - image[i-2, j+1]) + abs(image[i, j+1] - image[i-1, j+2]) + abs(image[i, j] - image[i-2, j+2]) + EPSILON
+                diag_gradient_factor["SE"] = abs(image[i+2, j+1] - image[i+1, j]) + abs(image[i-1, j] - image[i, j+1]) + abs(image[i+1, j+2] - image[i, j+1]) + abs(image[i, j-1] - image[i+1, j]) + abs(image[i-1, j-1] - image[i+1, j+1]) + abs(image[i+2, j+2] - image[i, j]) + EPSILON
+
+                
+                for k, v in diag_gradient_factor.items():
+                    weight[k] =  1 / v
+
+                num = 0
+                den = 0
+                for k, v in pre_estimation.items():
+                    num += v * weight[k]
+                    den += weight[k]
+
+                estimated_green[i, j] = num / den
+
+            else:
+                estimated_green[i, j] = image[i, j]
+
+
+    # Then we work on blue and red plans
+    
+    diag = [(-1, -1), (-1, 1), (1, -1), (1, 1)]
+    
+    for i in range(3, image.shape[0] - 4):
+        for j in range(3, image.shape[1] - 4):
+            if (is_red(i, j) or is_blue(i, j)):
+                ksi_nw = image[i-1, j-1] - estimated_green[i-1, j-1]
+                ksi_ne = image[i-1, j+1] - estimated_green[i-1, j+1]
+                ksi_sw = image[i+1, j-1] - estimated_green[i+1, j-1]
+                ksi_se = image[i+1, j+1] - estimated_green[i+1, j+1]
+                
+                g_nw = abs(estimated_green[i-2, j-1] - estimated_green[i-1, j]) + abs(estimated_green[i-1, j-2] - estimated_green[i, j-1]) + abs(estimated_green[i-2, j-2] - estimated_green[i-1, j-1]) + abs(estimated_green[i-1, j-1] - estimated_green[i, j]) + abs(image[i-1, j-1] - image[i+1, j+1]) + EPSILON
+                g_ne = abs(estimated_green[i-2, j+1] - estimated_green[i-1, j]) + abs(estimated_green[i-1, j+2] - estimated_green[i, j+1]) + abs(estimated_green[i-2, j+2] - estimated_green[i-1, j+1]) + abs(estimated_green[i-1, j+1] - estimated_green[i, j]) + abs(image[i-1, j+1] - image[i+1, j-1]) + EPSILON
+                g_sw = abs(estimated_green[i+2, j-1] - estimated_green[i+1, j]) + abs(estimated_green[i+1, j-2] - estimated_green[i, j-1]) + abs(estimated_green[i+2, j-2] - estimated_green[i+1, j-1]) + abs(estimated_green[i+1, j-1] - estimated_green[i, j]) + abs(image[i+1, j-1] - image[i-1, j+1]) + EPSILON
+                g_se = abs(estimated_green[i+2, j+1] - estimated_green[i+1, j]) + abs(estimated_green[i+1, j+2] - estimated_green[i, j+1]) + abs(estimated_green[i+2, j+2] - estimated_green[i+1, j+1]) + abs(estimated_green[i+1, j+1] - estimated_green[i, j]) + abs(image[i+1, j+1] - image[i-1, j-1]) + EPSILON
+
+                pre_est = estimated_green[i, j] + (ksi_nw/g_nw + ksi_ne/g_ne + ksi_sw/g_sw + ksi_se/g_se) / (1/g_nw + 1/g_ne + 1/g_sw + 1/g_se)
+
+                sum_dif = 0
+                sum_sig = 0
+
+                for coord in diag:
+                    mu = image[i + coord[0], j + coord[1]] - pre_est
+                    sum_dif += mu / (1 + abs(mu))
+                    sum_sig += 1 / (1 + abs(mu))
+
+                if (is_red(i, j)):
+                    estimated_blue[i, j] = pre_est
+                    estimated_red[i, j] = image[i, j]
+                else:
+                    estimated_red[i, j] = pre_est
+                    estimated_blue[i, j] = image[i, j]
+
+            else:
+                if (is_red(i-1, j)):
+                    estimated_red[i, j] = (image[i-1, j] + image[i+1, j]) / 2
+                    estimated_blue[i, j] = (image[i, j-1] + image[i, j+1]) / 2
+                else:
+                    estimated_red[i, j] = (image[i, j-1] + image[i, j+1]) / 2
+                    estimated_blue[i, j] = (image[i-1, j] + image[i+1, j]) / 2
+
+
+    res[:, :, 0] = estimated_red
+    res[:, :, 1] = estimated_green
+    res[:, :, 2] = estimated_blue
+
+    return res
\ No newline at end of file
diff --git a/src/methods/marty/quad_bayer.py b/src/methods/marty/quad_bayer.py
new file mode 100644
index 0000000000000000000000000000000000000000..c25953a132b9eecdc6094acfd615b1c409f2814b
--- /dev/null
+++ b/src/methods/marty/quad_bayer.py
@@ -0,0 +1,36 @@
+import numpy as np
+
+from scipy.signal import convolve2d
+
+
+def down_sample(image):    
+
+    down_sampled = np.empty((int(image.shape[0] / 2), int(image.shape[1] / 2)))
+
+    down_sampled[:, :] = (image[::2, ::2] + image[1::2, ::2] + image[::2, 1::2] + image[1::2, 1::2]) / 4
+
+    return down_sampled
+
+
+def up_sample(image):
+
+    up_sampled = np.empty((int(image.shape[0] * 2), int(image.shape[1] * 2), image.shape[2]))
+
+    up_sampled[::2, ::2, :] = image[::1, ::1, :]
+    up_sampled[1::2, ::2, :] = image[::1, ::1, :]
+    up_sampled[::2, 1::2, :] = image[::1, ::1, :]
+    up_sampled[1::2, 1::2, :] = image[::1, ::1, :]
+
+    return up_sampled
+
+
+
+def refine(image):
+    
+    res = np.empty((image.shape[0], image.shape[1], 3))
+    
+    ker = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 16
+    for i in range(3):
+        res[:, :, i] = convolve2d(image[:, :, i], ker, mode='same')
+
+    return res
diff --git a/src/methods/marty/reconstruct.py b/src/methods/marty/reconstruct.py
new file mode 100644
index 0000000000000000000000000000000000000000..729446e128b82d1556f73c1eddae35f939f0d153
--- /dev/null
+++ b/src/methods/marty/reconstruct.py
@@ -0,0 +1,73 @@
+"""The main file for the reconstruction.
+This file should NOT be modified except the body of the 'run_reconstruction' function.
+Students can call their functions (declared in others files of src/methods/your_name).
+"""
+
+
+import numpy as np
+
+from src.forward_model import CFA
+
+from src.methods.marty.mdwi import MDWI
+from src.methods.marty.high_quality_interpolation import high_quality_interpolation
+from src.methods.marty.quad_bayer import down_sample, up_sample, refine
+
+
+# Either MDWI or high_quality_interpolation
+method = high_quality_interpolation
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    
+    # Performing the reconstruction.
+    if (cfa == "bayer") :
+        
+        res = method(y)
+        
+
+    else :
+        down_sample_image = down_sample(y)
+
+        bayer_image = method(down_sample_image)
+        
+        res = up_sample(bayer_image)
+
+        res = refine(res)
+    return res
+
+
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller