diff --git a/src/methods/khattabi/KHATTABI.pdf b/src/methods/khattabi/KHATTABI.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..2f942c78a40e0ed68fd40fe4e721f9e2869b1f47
Binary files /dev/null and b/src/methods/khattabi/KHATTABI.pdf differ
diff --git a/src/methods/khattabi/reconstruct.py b/src/methods/khattabi/reconstruct.py
new file mode 100755
index 0000000000000000000000000000000000000000..c907d3703e9a265790638f24d46342fe1e5f2ad2
--- /dev/null
+++ b/src/methods/khattabi/reconstruct.py
@@ -0,0 +1,55 @@
+"""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.khattabi.utils import HA
+
+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.
+    # TODO
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+    y_adjoint=op.adjoint(y)
+    res = HA(y_adjoint)
+
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller
diff --git a/src/methods/khattabi/utils.py b/src/methods/khattabi/utils.py
new file mode 100755
index 0000000000000000000000000000000000000000..a87dfc0b223601a0653ca204bcfe5808359f8efd
--- /dev/null
+++ b/src/methods/khattabi/utils.py
@@ -0,0 +1,88 @@
+"""A file containing a (pretty useless) reconstruction.
+It serves as example of how the project works.
+This file should NOT be modified.
+"""
+
+
+import numpy as np
+from scipy.signal import convolve2d
+
+from src.forward_model import CFA
+
+def HA(y_adjoint):
+    height, width = y_adjoint.shape[:2]
+    # Green channel interpolation
+
+
+    for i in range(2 , height-2): 
+        for j in range(2, width-2):
+            # Red pixels
+            if y_adjoint[i,j,1] == 0 and y_adjoint[i,j,0] != 0:  
+                grad_y = np.abs(y_adjoint[i-1,j,1] - y_adjoint[i+1,j,1]) + np.abs(2*y_adjoint[i,j,0] - y_adjoint[i-2,j,0] - y_adjoint[i+2,j,0])
+                grad_x = np.abs(y_adjoint[i,j-1,1] - y_adjoint[i,j+1,1]) + np.abs(2*y_adjoint[i,j,0] - y_adjoint[i,j-2,0] - y_adjoint[i,j+2,0])
+
+                if grad_x < grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/2 + (2*y_adjoint[i,j,0] - y_adjoint[i,j-2,0] - y_adjoint[i,j+2,0])/4
+                elif  grad_x > grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1])/2 + (2*y_adjoint[i,j,0] - y_adjoint[i-2,j,0] - y_adjoint[i+2,j,0])/4
+                else:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1] + y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/4 + \
+                                       (2*y_adjoint[i,j,0] - y_adjoint[i,j-2,0] - y_adjoint[i,j+2,0] + 2*y_adjoint[i,j,0] - y_adjoint[i-2,j,0] - y_adjoint[i+2,j,0])/8
+            # Blue Pixels
+            elif y_adjoint[i,j,1] == 0 and y_adjoint[i,j,2] != 0:   
+                grad_y = np.abs(y_adjoint[i-1,j,1] - y_adjoint[i+1,j,1]) + np.abs(2*y_adjoint[i,j,2] - y_adjoint[i-2,j,2] - y_adjoint[i+2,j,2])
+                grad_x = np.abs(y_adjoint[i,j-1,1] - y_adjoint[i,j+1,1]) + np.abs(2*y_adjoint[i,j,2] - y_adjoint[i,j-2,2] - y_adjoint[i,j+2,2])
+
+                if grad_x < grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/2 + (2*y_adjoint[i,j,2] - y_adjoint[i,j-2,2] - y_adjoint[i,j+2,2])/4
+                elif  grad_x > grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1])/2 + (2*y_adjoint[i,j,2] - y_adjoint[i-2,j,2] - y_adjoint[i+2,j,2])/4
+                else:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1] + y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/4 + \
+                                       (2*y_adjoint[i,j,2] - y_adjoint[i,j-2,2] - y_adjoint[i,j+2,2] + 2*y_adjoint[i,j,2] - y_adjoint[i-2,j,2] - y_adjoint[i+2,j,2])/8
+                    
+    # Red and blue channel interpolation
+    for i in range(1 , height-1):  
+        for j in range(1, width-1):
+            if y_adjoint[i,j,2] != 0 :
+                 y_adjoint[i,j,0] = (y_adjoint[i-1,j-1,0] + y_adjoint[i-1,j+1,0] + y_adjoint[i+1,j-1,0] + y_adjoint[i+1,j+1,0]) / 4
+
+            elif y_adjoint[i,j,0] != 0:
+                 y_adjoint[i,j,2] = (y_adjoint[i-1,j-1,2] + y_adjoint[i-1,j+1,2] + y_adjoint[i+1,j-1,2] + y_adjoint[i+1,j+1,2]) / 4
+            
+            else:
+                y_adjoint[i,j] = y_adjoint[i,j]
+
+    for i in range(1 , height-1):
+        for j in range(1, width-1):    
+            if y_adjoint[i,j,0] == y_adjoint[i,j,2]:
+                y_adjoint[i,j,0] = (y_adjoint[i-1,j,0] + y_adjoint[i,j-1,0] + y_adjoint[i+1,j,0] + y_adjoint[i,j+1,0]) / 4             
+                y_adjoint[i,j,2] = (y_adjoint[i-1,j,2] + y_adjoint[i,j-1,2] + y_adjoint[i+1,j,2] + y_adjoint[i,j+1,2]) / 4  
+            
+    return y_adjoint
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller