diff --git a/src/methods/lioretn/demosaicking.py b/src/methods/lioretn/demosaicking.py
new file mode 100644
index 0000000000000000000000000000000000000000..7b2ed87c4059aa36a5366c366a51fba9fc8e83f5
--- /dev/null
+++ b/src/methods/lioretn/demosaicking.py
@@ -0,0 +1,162 @@
+"""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
+
+# High quality linear interpolation filters for bayer mosaic
+bayer_g_at_r = 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
+bayer_g_at_b = bayer_g_at_r
+bayer_r_at_green_rrow_bcol = np.array([[ 0,  0, 1/2,  0,  0],
+                                       [ 0, -1,  0,  -1,  0],
+                                       [-1,  4,  5,   4, -1],
+                                       [ 0, -1,  0,  -1,  0],
+                                       [ 0,  0, 1/2,  0,  0]]) / 8
+bayer_r_at_green_brow_rcol = bayer_r_at_green_rrow_bcol.T
+bayer_b_at_green_rrow_bcol = bayer_r_at_green_brow_rcol
+bayer_b_at_green_brow_rcol = bayer_r_at_green_rrow_bcol
+bayer_r_at_b = np.array([[  0,  0, -3/2, 0,   0],
+                            [  0,  2,   0,  2,   0],
+                            [-3/2, 0,   6,  0, -3/2],
+                            [  0,  2,   0,  2,   0],
+                            [  0,  0, -3/2, 0,   0]]) / 8
+bayer_b_at_r = bayer_r_at_b
+
+# High quality linear interpolation filters for quad mosaic
+quad_g_at_r = np.array([[ 0,  0,  0, 0, -1, -1, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0, -1, -1, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0,  2,  2, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0,  2,  2, 0, 0,  0,  0],
+                        [-1, -1,  2, 2,  4,  4, 2, 2, -1, -1],
+                        [-1, -1,  2, 2,  4,  4, 2, 2, -1, -1],
+                        [ 0,  0,  0, 0,  2,  2, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0,  2,  2, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0, -1, -1, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0, -1, -1, 0, 0,  0,  0]]) / 32
+quad_g_at_b = quad_g_at_r
+quad_r_at_green_rrow_bcol = np.array([[ 0,  0,  0,  0, 1/2, 1/2,  0,  0,  0,  0],
+                                      [ 0,  0,  0,  0, 1/2, 1/2,  0,  0,  0,  0],
+                                      [ 0,  0, -1, -1,  0,   0,  -1, -1,  0,  0],
+                                      [ 0,  0, -1, -1,  0,   0,  -1, -1,  0,  0],
+                                      [-1, -1,  4,  4,  5,   5,   4,  4, -1, -1],
+                                      [-1, -1,  4,  4,  5,   5,   4,  4, -1, -1],
+                                      [ 0,  0, -1, -1,  0,   0,  -1, -1,  0,  0],
+                                      [ 0,  0, -1, -1,  0,   0,  -1, -1,  0,  0],
+                                      [ 0,  0,  0,  0, 1/2, 1/2,  0,  0,  0,  0],
+                                      [ 0,  0,  0,  0, 1/2, 1/2,  0,  0,  0,  0]]) / 32
+quad_r_at_green_brow_rcol = quad_r_at_green_rrow_bcol.T
+quad_b_at_green_rrow_bcol = quad_r_at_green_brow_rcol
+quad_b_at_green_brow_rcol = quad_r_at_green_rrow_bcol
+quad_r_at_b = np.array([[  0,    0,  0,  0, -3/2, -3/2, 0,  0,   0,    0],
+                        [  0,    0,  0,  0, -3/2, -3/2, 0,  0,   0,    0],
+                        [  0,    0,  2,  2,   0,    0,  2,  2,   0,    0],
+                        [  0,    0,  2,  2,   0,    0,  2,  2,   0,    0],
+                        [-3/2, -3/2, 0,  0,   6,    6,  0,  0, -3/2, -3/2],
+                        [-3/2, -3/2, 0,  0,   6,    6,  0,  0, -3/2, -3/2],
+                        [  0,    0,  2,  2,   0,    0,  2,  2,   0,    0],
+                        [  0,    0,  2,  2,   0,    0,  2,  2,   0,    0],
+                        [  0,    0,  0,  0, -3/2, -3/2, 0,  0,   0,    0],
+                        [  0,    0,  0,  0, -3/2, -3/2, 0,  0,   0,    0]]) / 32
+quad_b_at_r = quad_r_at_b
+
+def high_quality_linear_interpolation(op : CFA, y : np.ndarray) -> np.ndarray:
+    z = op.adjoint(y)
+    res = np.empty(op.input_shape)
+    mask = op.mask
+    R = 0 ; G = 1 ; B = 2
+
+    if op.cfa == 'bayer':
+        y_pad = np.pad(y, 2, 'constant', constant_values=0)
+
+        for i in range(y.shape[0]):
+            for j in range(y.shape[1]):
+                i_pad = i+2 ; j_pad = j+2
+
+                if mask[i, j, G] == 1: # We must estimate R and B components
+                    res[i, j, G] = y[i, j]
+
+                    if mask[i, max(0, j-1), R] == 1 or mask[i, min(y.shape[1]-1, j+1), R] == 1:
+                        res[i, j, R] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_r_at_green_rrow_bcol, mode='valid'))
+                    else:
+                        res[i, j, R] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_r_at_green_brow_rcol, mode='valid'))
+
+                    if mask[i, max(0, j-1), B] == 1 or mask[i, min(y.shape[1]-1, j+1), B] == 1:
+                        res[i, j, B] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_b_at_green_brow_rcol, mode='valid'))
+                    else:
+                        res[i, j, B] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_b_at_green_rrow_bcol, mode='valid'))
+
+                elif mask[i, j, R] == 1: # We must estimate G and B components
+                    res[i, j, R] = y[i, j]
+                    res[i, j, G] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_g_at_r, mode='valid'))
+                    res[i, j, B] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_b_at_r, mode='valid'))
+
+                else: # We must estimate R and G components
+                    res[i, j, B] = y[i, j]
+                    res[i, j, G] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_g_at_b, mode='valid'))
+                    res[i, j, R] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_r_at_b, mode='valid'))
+
+    else:
+        y_pad = np.pad(y, 4, 'constant', constant_values=0)
+
+        for i in range(0, y.shape[0], 2):
+            for j in range(0, y.shape[1], 2):
+                i_pad = i+4 ; j_pad = j+4
+
+                if mask[i, j, G] == 1: # We must estimate R and B components
+                    res[i:i+2, j:j+2, G] = y[i:i+2, j:j+2]
+
+                    if mask[i, max(0, j-1), R] == 1 or mask[i, min(y.shape[1]-1, j+1), R] == 1:
+                        res[i:i+2, j:j+2, R] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_r_at_green_rrow_bcol, mode='valid'))
+                    else:
+                        res[i:i+2, j:j+2, R] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_r_at_green_brow_rcol, mode='valid'))
+
+                    if mask[i, max(0, j-1), B] == 1 or mask[i, min(y.shape[1]-1, j+1), B] == 1:
+                        res[i:i+2, j:j+2, B] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_b_at_green_brow_rcol, mode='valid'))
+                    else:
+                        res[i:i+2, j:j+2, B] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_b_at_green_rrow_bcol, mode='valid'))
+
+                elif mask[i, j, R] == 1: # We must estimate G and B components
+                    res[i:i+2, j:j+2, R] = y[i:i+2, j:j+2]
+                    res[i:i+2, j:j+2, G] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_g_at_r, mode='valid'))
+                    res[i:i+2, j:j+2, B] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_b_at_r, mode='valid'))
+
+                else: # We must estimate R and G components
+                    res[i:i+2, j:j+2, B] = y[i:i+2, j:j+2]
+                    res[i:i+2, j:j+2, G] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_g_at_b, mode='valid'))
+                    res[i:i+2, j:j+2, R] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_r_at_b, mode='valid'))
+
+    return res
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller
diff --git a/src/methods/lioretn/lioretn.pdf b/src/methods/lioretn/lioretn.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..6b152311853467d297d10a3b641f500b411bb6d3
Binary files /dev/null and b/src/methods/lioretn/lioretn.pdf differ
diff --git a/src/methods/lioretn/reconstruct.py b/src/methods/lioretn/reconstruct.py
new file mode 100644
index 0000000000000000000000000000000000000000..80f771adc5b319cc3b90196af398a5469dc3e743
--- /dev/null
+++ b/src/methods/lioretn/reconstruct.py
@@ -0,0 +1,57 @@
+"""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.lioretn.demoisaicing_fct import High_Quality_Linear_Interpolation
+from src.methods.lioretn.demosaicking import high_quality_linear_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.
+    # TODO
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+
+    res = high_quality_linear_interpolation(op, y)
+    # res = High_Quality_Linear_Interpolation(op, y)
+
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller