diff --git a/src/methods/charpentier_laurine/charpentier.pdf b/src/methods/charpentier_laurine/charpentier.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..ee757501680296067200cb5d7d4c22a11dd6db99
Binary files /dev/null and b/src/methods/charpentier_laurine/charpentier.pdf differ
diff --git a/src/methods/charpentier_laurine/functions.py b/src/methods/charpentier_laurine/functions.py
new file mode 100644
index 0000000000000000000000000000000000000000..78890250d247f1bed19b9f986af9e5f3eb98961c
--- /dev/null
+++ b/src/methods/charpentier_laurine/functions.py
@@ -0,0 +1,79 @@
+import numpy as np
+import cv2
+from scipy.signal import convolve2d
+
+from src.forward_model import CFA
+
+
+def superpixel(op: CFA, y: np.ndarray) -> np.ndarray:
+    """Performs the method of variable number of gradients
+
+    Args:
+        op (CFA): CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    z = op.adjoint(y)
+
+    if op.cfa == 'bayer':
+        res = np.empty((op.input_shape[0]//2, op.input_shape[1]//2, op.input_shape[2]))
+
+        for i in range(1,op.input_shape[0],2):
+            for j in range(1, op.input_shape[1],2):
+                res[i//2,j//2,0] = z[ i-1,j ,0]
+                res[i//2,j//2,1] = (z[i,j,1] + z[i-1,j-1,1]) / 2
+                res[i//2,j//2,2] = z[ i,j-1 ,2]
+
+    else:
+        res = np.empty((op.input_shape[0]//2, op.input_shape[1]//2, op.input_shape[2]))
+
+        print()
+
+        for i in range(2,op.input_shape[0]-5,4):
+            for j in range(2, op.input_shape[1]-5,4):
+
+                res[i//2,j//2,0] = z[ i-2,j ,0]
+                res[i//2,j//2,1] = (z[i,j,1] + z[i-2,j-2,1]) / 2
+                res[i//2,j//2,2] = z[ i,j-2 ,2]
+
+                res[i//2+1,j//2,0] = z[ i-2+1,j ,0]
+                res[i//2+1,j//2,1] = (z[i+1,j,1] + z[i-2+1,j-2,1]) / 2
+                res[i//2+1,j//2,2] = z[ i+1,j-2 ,2]
+
+                res[i//2+1,j//2+1,0] = z[ i-2+1,j+1 ,0]
+                res[i//2+1,j//2+1,1] = (z[i+1,j+1,1] + z[i-2+1,j-2+1,1]) / 2
+                res[i//2+1,j//2+1,2] = z[ i+1,j-2+1 ,2]
+
+                res[i//2,j//2+1,0] = z[ i-2,j+1 ,0]
+                res[i//2,j//2+1,1] = (z[i,j+1,1] + z[i-2,j-2+1,1]) / 2
+                res[i//2,j//2+1,2] = z[ i,j-2+1 ,2]
+
+    return cv2.resize(res, (op.input_shape[0], op.input_shape[1]))
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller
diff --git a/src/methods/charpentier_laurine/reconstruct.py b/src/methods/charpentier_laurine/reconstruct.py
new file mode 100755
index 0000000000000000000000000000000000000000..062ada46b47d9d0cd006b645b423f2d9ae8a3dd2
--- /dev/null
+++ b/src/methods/charpentier_laurine/reconstruct.py
@@ -0,0 +1,54 @@
+"""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.charpentier_laurine.functions import superpixel
+
+
+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.
+    """
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+
+    res = superpixel(op, y)
+
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller