src/methods/Yosra_Jelassi/new_reconstruction.py

+"""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 new_interpolation(op: CFA, y: np.ndarray) -> np.ndarray:
+    """Performs a simple interpolation of the lost pixels.
+
+    Args:
+        op (CFA): CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    z = op.adjoint(y)
+    if op.cfa == 'bayer':
+        demosaicked_red = np.empty_like(z[:,:,0])
+
+        weights_red = np.array([[0, 1],                      
+                                [0, 0]])   
+
+        for i in range(1, z.shape[0] - 1,2):
+            for j in range(1, z.shape[1] - 1,2):
+                demosaicked_red[i-1:i+1, j-1:j+1] = np.sum(z[i-1:i+1, j-1:j+1,0] * weights_red)
+   
+        demosaicked_green = np.empty_like(z[:,:,1])
+
+        weights_green = np.array([[1, 0],
+                                  [0, 1]]) 
+    
+        for i in range(1, z.shape[0] - 1,2):
+            for j in range(1, z.shape[1] - 1,2):
+                demosaicked_green[i-1:i+1, j-1:j+1] = np.sum(z[i-1:i+1, j-1:j+1,1] * weights_green)/2
+
+        weights_blue = np.array([[0, 0],                      
+                                 [1, 0]])  
+        demosaicked_blue = np.empty_like(z[:,:,2])
+
+        for i in range(1, z.shape[0] - 1,2):
+            for j in range(1, z.shape[1] - 1,2):
+                demosaicked_blue[i-1:i+1, j-1:j+1] = np.sum(z[i-1:i+1, j-1:j+1,2] * weights_blue)
+
+        demosaicked_image = np.dstack((demosaicked_red, demosaicked_green, demosaicked_blue))
+    else:
+        demosaicked_red = np.empty_like(z[:,:,0])
+
+        weights_red = np.array([[0, 0, 1, 1 ],                      
+                                [0, 0, 1, 1 ],
+                                [0, 0, 0, 0 ],
+                                [0, 0, 0, 0 ]])   
+
+        for i in range(2, z.shape[0] - 2, 4):
+            for j in range(2, z.shape[1] - 2, 4):
+                demosaicked_red[i-2:i+2, j-2:j+2] = np.sum(z[i-2:i+2, j-2:j+2,0] * weights_red)/4
+   
+        demosaicked_green = np.empty_like(z[:,:,1])
+
+        weights_green = np.array([[1, 1, 0, 0 ],                      
+                                  [1, 1, 0, 0 ],
+                                  [0, 0, 1, 1 ],
+                                  [0, 0, 1, 1 ]])   
+    
+        for i in range(2, z.shape[0] - 2, 4):
+            for j in range(2, z.shape[1] - 2, 4):
+                demosaicked_green[i-2:i+2, j-2:j+2] = np.sum(z[i-2:i+2, j-2:j+2,1] * weights_green)/8
+
+        weights_blue = np.array([[0, 0, 0, 0 ],                      
+                                 [0, 0, 0, 0 ],
+                                 [1, 1, 0, 0 ],
+                                 [1, 1, 0, 0 ]])  
+        demosaicked_blue = np.empty_like(z[:,:,2])
+
+        for i in range(2, z.shape[0] - 2, 4):
+            for j in range(2, z.shape[1] - 2, 4):
+                demosaicked_blue[i-2:i+2, j-2:j+2] = np.sum(z[i-2:i+2, j-2:j+2,2] * weights_blue)/4
+
+        demosaicked_image = np.dstack((demosaicked_red, demosaicked_green, demosaicked_blue))
+        
+    return demosaicked_image
+
+
+
+def extract_padded(M, size, i, j):
+    N_i, N_j = M.shape
+    res = np.zeros((size, size))
+    middle_size = int((size - 1) / 2)
+
+    for ii in range(- middle_size, middle_size + 1):
+        for jj in range(- middle_size, middle_size + 1):
+            if i + ii >= 0 and i + ii < N_i and j + jj >= 0 and j + jj < N_j:
+                res[middle_size + ii, middle_size + jj] = M[i + ii, j + jj]
+
+    return res
+
+
+def varying_kernel_convolution(M, K_list):
+    N_i, N_j = M.shape
+    res = np.zeros_like(M)
+
+    for i in range(N_i):
+        for j in range(N_j):
+            res[i, j] = np.sum(extract_padded(M, K_list[4 * (i % 4) + j % 4].shape[0], i, j) * K_list[4 * (i % 4) + j % 4])
+
+    np.clip(res, 0, 1, res)
+
+    return res
+
+
+K_identity = np.zeros((5, 5))
+K_identity[2, 2] = 1
+
+K_red_0 = np.zeros((5, 5))
+K_red_0[2, :] = np.array([-3, 13, 0, 0, 2]) / 12
+
+K_red_1 = np.zeros((5, 5))
+K_red_1[2, :] = np.array([2, 0, 0, 13, -3]) / 12
+
+K_red_8 = np.zeros((5, 5))
+K_red_8[:2, :2] = np.array([[-1, -1], [-1, 9]]) / 6
+
+K_red_9 = np.zeros((5, 5))
+K_red_9[:2, 3:] = np.array([[-1, -1], [9, -1]]) / 6
+
+K_red_10 = np.zeros((5, 5))
+K_red_10[:, 2] = np.array([-3, 13, 0, 0, 2]) / 12
+
+K_red_12 = np.zeros((5, 5))
+K_red_12[3:, :2] = np.array([[-1, 9], [-1, -1]]) / 6
+
+K_red_13 = np.zeros((5, 5))
+K_red_13[3:, 3:] = np.array([[9, -1], [-1, -1]]) / 6
+
+K_red_14 = np.zeros((5, 5))
+K_red_14[:, 2] = np.array([2, 0, 0, 13, -3]) / 12
+
+K_list_red = [K_red_0, K_red_1, K_identity, K_identity, K_red_0, K_red_1, K_identity, K_identity, K_red_8, K_red_9, K_red_10, K_red_10, K_red_12, K_red_13, K_red_14, K_red_14]
+
+
+K_green_2 = np.zeros((5, 5))
+K_green_2[2, :] = [-3, 13, 0, 0, 2]
+K_green_2[:, 2] = [-3, 13, 0, 0, 2]
+K_green_2 = K_green_2 / 24
+
+K_green_3 = np.zeros((5, 5))
+K_green_3[2, :] = [2, 0, 0, 13, -3]
+K_green_3[:, 2] = [-3, 13, 0, 0, 2]
+K_green_3 = K_green_3 / 24
+
+K_green_6 = np.zeros((5, 5))
+K_green_6[2, :] = [-3, 13, 0, 0, 2]
+K_green_6[:, 2] = [2, 0, 0, 13, -3]
+K_green_6 = K_green_6 / 24
+
+K_green_7 = np.zeros((5, 5))
+K_green_7[2, :] = [2, 0, 0, 13, -3]
+K_green_7[:, 2] = [2, 0, 0, 13, -3]
+K_green_7 = K_green_7 / 24
+
+K_list_green = [K_identity, K_identity, K_green_2, K_green_3, K_identity, K_identity, K_green_6, K_green_7, K_green_2, K_green_3, K_identity, K_identity, K_green_6, K_green_7, K_identity, K_identity]
+
+
+K_list_blue = [K_red_10, K_red_10, K_red_8, K_red_9, K_red_14, K_red_14, K_red_12, K_red_13, K_identity, K_identity, K_red_0, K_red_1, K_identity, K_identity, K_red_0, K_red_1]
+
+
+
+ker_bayer_red_blue = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 4
+
+ker_bayer_green = np.array([[0, 1, 0], [1, 4, 1], [0, 1, 0]]) / 4
+
+
+

src/reconstruct.py → src/methods/Yosra_Jelassi/reconstruct.py

 """The main file for the reconstruction.
 This file should NOT be modified except the body of the 'run_reconstruction' function.
-Students should only call their functions (declared in others files of src/) in it.
+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.Yosra_Jelassi.new_reconstruction import new_interpolation
 
-# Import the files needed for the reconstruction like:
-# from src.file import function_1
-from src.demo_reconstructions import naive_interpolation
 
 
 def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
@@ -18,20 +16,16 @@ def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
 
     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
-
-    ######################################## Demo code ########################################
     input_shape = (y.shape[0], y.shape[1], 3)
     op = CFA(cfa, input_shape)
-
-
-    res = naive_interpolation(op, y)
-    ###########################################################################################
+    res = new_interpolation(op, y)
 
     return res
 

src/methods/anouk_enz/reconstruct.py

+"""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
+import cv2 as cv
+from skimage.filters import rank
+from src.forward_model import CFA
+from skimage.morphology import disk
+import matplotlib.pyplot as plt
+
+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)
+
+    z = op.adjoint(y)
+
+
+    ################# MY CODE ############################
+    if op.cfa == 'quad_bayer' : 
+        # Transform into a Bayer-patterned image 
+        y = QuadBayer_to_Bayer(y, op, vizualize=False)
+
+        # Demosaïck the image 
+        res = demosaicking_bayer(y)
+       
+    if op.cfa == 'bayer':
+        res = demosaicking_bayer(y)
+    ####################################################
+
+    return res
+
+def get_color(mask, i, j):
+    """
+        Returns the index of the chanel that contains already a color at a pixel localization. 
+
+        Args :
+            mask (np.ndarray) : mask of the used cfa.
+            i (int): index of the line of the pixel.
+            j (int): index of the column of the pixel.
+
+        Returns 
+            int: the index of the color avaible at the indicated pixel.
+    """
+    if mask[i,j,0] == 1:
+        return 0
+    if mask[i,j,1] == 1:
+        return 1
+    if mask[i,j,2] == 1:
+        return 2
+    return None
+
+def demosaicking_bayer(y):
+    """
+        Performs High-Quality Linear Interpolation on a Bayer-patterned image.
+
+        Args : 
+            y (np.ndarray): Mosaicked image to be reconstructed.
+       
+            Returns 
+                np.ndarray: Demosaicked image. 
+    """
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA("bayer", input_shape)
+    res = np.empty(op.input_shape)
+
+    for i in range (2, input_shape[0]-2):
+        for j in range (2, input_shape[0]-2):
+            patch = y[i-2:i+3, j-2:j+3]
+            channel = get_color(op.mask, i, j)
+            # If the Red value is avaible 
+            if channel == 0:
+                res[i, j, 0] = y[i,j]
+                # Interpolate the Green value
+                res[i, j, 1] = np.sum(patch*G_at_R)/8
+                #Interpolate the Blue value 
+                res[i, j, 2] = np.sum(patch* B_at_R)/8
+
+            # If the Green value is avaible 
+            if channel == 1 : 
+                res[i, j, 1] = y[i,j]
+                # Interpolation of the Red value
+                if get_color(op.mask, i, j+1) == 0:
+                    res[i,j,0] = np.sum(patch*R_at_G_row)/8
+                else : 
+                    res[i,j,0] = np.sum(patch*R_at_G_column)/8
+                
+                # Interpolation of the Blue value
+                if get_color(op.mask, i, j+1) == 2:
+                    res[i,j,2] = np.sum(patch* B_at_G_row)/8
+                else : 
+                    res[i,j,2] = np.sum(patch * B_at_G_column)/8
+
+            # If the Blue value is avaible 
+            if channel == 2:
+                res[i, j, 2] = y[i,j]
+                # Interpolate the Red value
+                res[i, j, 1] = np.sum(patch* G_at_B)/8
+                #Interpolate the Red value 
+                res[i, j, 0] = np.sum(patch* R_at_B)/8
+    return res
+
+def QuadBayer_to_Bayer(y, op, vizualize=False):
+    """
+    Applies the swapping method to transform a QuadBayer-patterned image into a Bayer-patterned image
+
+    Args : 
+        y (np.ndarray): Mosaicked image to be tranformed 
+        op (CFA) : the CFA object 
+        vizualize (Boolean) : show the evolution of the mask if True
+
+    Returns :
+        (np.ndarray): a bayer-patterned image 
+    """
+    input_shape = (y.shape[0], y.shape[1], 3)
+
+    if vizualize : 
+        fig, axs = plt.subplots(1, 4, figsize=(15,15))
+        axs[0].imshow(op.mask[0:10, 0:10, :])
+        axs[0].set_title('Original mask')
+
+    # Step 1 : Swap 2 columns every 2 columns
+    temp = np.zeros((input_shape[0], 1))
+    temp_mask = np.zeros((input_shape[0], 1, 3))
+    for col in range (1, input_shape[0]-1, 4):
+        temp = np.copy(y[:, col])
+        y[:, col] =  np.copy(y[:, col+1])
+        y[:, col+1] = np.copy(temp)
+
+        if vizualize: 
+            temp_mask = np.copy(op.mask[:, col,:])
+            op.mask[:, col,:] =  np.copy(op.mask[:, col+1,:])
+            op.mask[:, col+1,:] = np.copy(temp_mask)
+    if vizualize:         
+        axs[1].imshow(op.mask[0:10, 0:10, :])
+        axs[1].set_title('Mask after first step')
+        
+    #Step 2 : Swap 2 lines every 2 lines 
+    temp = np.zeros((1, input_shape[1]))
+    temp_mask = np.zeros((1, input_shape[1],3))
+    for line in range (1, input_shape[1], 4):
+        temp = np.copy(y[line, :])
+        y[line, :] =  np.copy(y[line+1, :])
+        y[line+1, :] = np.copy(temp)
+
+        if vizualize: 
+            temp_mask = np.copy(op.mask[line, :, :])
+            op.mask[line, :, :] =  np.copy(op.mask[line+1, :, :])
+            op.mask[line+1, :, :] = np.copy(temp_mask)
+    if vizualize: 
+        axs[2].imshow(op.mask[0:10, 0:10, :])
+        axs[2].set_title('Mask after second step')
+
+    # 3: Swap back some diagonal greens
+    temp_mask = np.zeros((1,1,3))
+    for i in range(0, input_shape[0], 4):
+        for j in range(2, input_shape[1], 4):
+            temp = y[i, j]
+            y[i, j] = y[i+1, j-1]
+            y[i+1, j-1] = temp
+            if vizualize: 
+                temp_mask = op.mask[i, j, :]
+                op.mask[i, j, :] = np.copy(op.mask[i+1, j-1, :])
+                op.mask[i+1, j-1, :] = np.copy(temp_mask)
+    if vizualize: 
+        axs[3].imshow(op.mask[0:10, 0:10, :])
+        axs[3].set_title('Mask after third step')
+        plt.tight_layout()
+        plt.show()
+
+    return y
+
+
+# Defining the mask
+
+# Interpolate the Green 
+# Green at Red 
+G_at_R = [[ 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]]
+# Green at Blue
+G_at_B = [[ 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]]
+
+#Interpolate the Red 
+# R at G in R row, B column
+R_at_G_row = [[ 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]] 
+# R at G in B row, R column
+R_at_G_column = [[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]]
+# Red at Blue
+R_at_B = [[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]]
+
+#Intelropation for Blue 
+# Blue at Green in B row, R column
+B_at_G_row = [[ 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]] 
+
+B_at_G_column = [[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]]
+
+B_at_R = [[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]]
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/antoine_vouillon/functions.py

+import numpy as np
+
+def find_neighbors (z,chan,i,j,N,M):
+    
+    """Finds all the neighbors of a pixel on a given channel
+    Args:
+        z: image.
+        chan: chosen chanel.
+        i: line of the pixel
+        j: column of the pixel
+        N: number of lines
+        M: number of columns
+
+    Returns:
+        np.ndarray: Neighbors of the pixel in a list (including the pixel itself)
+    """
+
+    P1 = z[(i-1)%N,(j-1)%M,chan]
+    P2 = z[(i-1)%N,j%M,chan]
+    P3 = z[(i-1)%N,(j+1)%M,chan]
+    P4 = z[i%N,(j-1)%M,chan]
+    P5 = z[i%N,j%M,chan]
+    P6 = z[i%N,(j+1)%M,chan]
+    P7 = z[(i+1)%N,(j-1)%M,chan]
+    P8 = z[(i+1)%N,j%M,chan]
+    P9 = z[(i+1)%N,(j+1)%M,chan]
+
+    return np.array([P1,P2,P3,P4,P5,P6,P7,P8,P9])
+
+def find_dir_deriv(neighbors):
+    """Calculates the directional derivative of a pixel.
+
+    Args:
+        neighbors: list of the neighbors of the pixel
+
+    Returns:
+        np.ndarray: directional derivatives in this order: Dx, Dy, Ddx, Ddy
+    """
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neighbors
+    Dx = (P4 - P6)/2
+    Dy = (P2 - P8)/2
+    Dxd = (P3 - P7)/(2*np.sqrt(2))
+    Dyd = (P1 - P9)/(2*np.sqrt(2))
+
+    return [Dx,Dy,Dxd,Dyd]
+
+def find_weights(z, neigh, dir_deriv,chan,i,j,N,M):
+    """Finds all the neighbors of a pixel on a given channel
+    Args:
+        z: image.
+        dir_deriv: directional derivatives
+        chan: chosen chanel.
+        i: line of the pixel
+        j: column of the pixel
+        N: number of lines
+        M: number of columns
+
+    Returns:
+        np.ndarray: Weights from E1 to E9
+    """
+    [Dx,Dy,Dxd,Dyd] = dir_deriv
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neigh
+    E = []
+    c = 1
+    for k in range (-1,2):
+        for k in range (-1,2):
+
+            n = find_neighbors(z,chan,i+k,j+k,N,M)
+            dd = find_dir_deriv(n)
+            if c == 1 or c == 9:
+                E.append(1/np.sqrt(1 + Dyd**2 + dd[3]**2))
+            elif c == 2 or c == 8:
+                E.append(1/np.sqrt(1 + Dy**2 + dd[1]**2))
+            elif c == 3 or c == 7:
+                E.append(1/np.sqrt(1 + Dxd**2 + dd[2]**2))
+            elif c == 4 or c == 6:
+                E.append(1/np.sqrt(1 + Dx**2 + dd[0]**2))
+            c += 1
+    return E       
+
+def interpolate(neigh,weights):
+    
+    """interpolates pixels from a grid where one of two pixels is missing regularly spaced
+    Args:
+        neigh: neighbors of the pixel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the interpolated pixel
+    """
+
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neigh
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+    num5 = E2*P2 + E4*P4 + E6*P6 + E8*P8
+    den5 = E2 + E4 + E6 + E8
+    I5 = num5/den5
+    return I5
+
+def interpolate_RB(neigh, neigh_G, weights):
+
+    """interpolates the central missing pixel from the red or blue channel from a bayer patern
+    Args:
+        neigh: neighbors of the pixel in the red or blue channel.
+        neigh_G: neighbors of the pixel in the green channel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the interpolated pixel in the red or blue channel
+    """
+
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neigh
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neigh_G
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+    num5 = ((E1*P1)/G1) + ((E3*P3)/G3) + ((E7*P7)/G7) + ((E9*P9)/G9)
+    den5 = E1 + E3 + E7 + E9
+    I5 = G5 * num5/den5
+
+    return I5
+
+def correction_G(neigh_G ,neigh_R ,neigh_B, weights):
+
+    """corrects the value of the green pixel in the third phase of the kimmel algorythme
+    Args:
+
+        neigh_G: neighbors of the pixel in the green channel.
+        neigh_R: neighbors of the pixel in the red channel.
+        neigh_B: neighbors of the pixel in the blue channel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the corrected pixel in the green channel
+    """
+
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neigh_G
+    [R1,R2,R3,R4,R5,R6,R7,R8,R9] = neigh_R
+    [B1,B2,B3,B4,B5,B6,B7,B8,B9] = neigh_B
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+
+    num_Gb5 = ((E2*G2)/B2) + ((E4*G4)/B4) + ((E6*G6)/B6) + ((E8*G8)/B8)
+    num_Gr5 = ((E2*G2)/R2) + ((E4*G4)/R4) + ((E6*G6)/R6) + ((E8*G8)/R8)
+    den5 = E2 + E4 + E6 + E8
+
+    Gb5 = B5 * num_Gb5/den5
+    Gr5 = R5 * num_Gr5/den5
+
+    G5 = (Gb5 + Gr5)/2
+
+    return Gr5
+
+def correction_R(neigh_G, neigh_R, weights):
+
+    """corrects the value of the red pixel in the third phase of the kimmel algorythme
+    Args:
+
+        neigh_G: neighbors of the pixel in the green channel.
+        neigh_R: neighbors of the pixel in the red channel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the corrected pixel in the red channel
+    """
+
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neigh_G
+    [R1,R2,R3,R4,R5,R6,R7,R8,R9] = neigh_R
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+
+    num_R5 = ((E1*R1)/G1) + ((E2*R2)/G2) + ((E3*R3)/G3) + ((E4*R4)/G4) + ((E6*R6)/G6) + ((E7*R7)/G7) + ((E8*R8)/G8) + ((E9*R9)/G9)
+
+    den5 = sum(weights)
+
+    R5 = G5 * num_R5/den5
+
+    return R5
+
+def correction_B(neigh_G ,neigh_B, weights):
+
+    """corrects the value of the blue pixel in the third phase of the kimmel algorythme
+    Args:
+
+        neigh_G: neighbors of the pixel in the green channel.
+        neigh_B: neighbors of the pixel in the blue channel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the corrected pixel in the blue channel
+    """
+
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neigh_G
+    [B1,B2,B3,B4,B5,B6,B7,B8,B9] = neigh_B
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+
+    num_B5 = ((E1*B1)/G1) + ((E2*B2)/G2) + ((E3*B3)/G3) + ((E4*B4)/G4) + ((E6*B6)/G6) + ((E7*B7)/G7) + ((E8*B8)/G8) + ((E9*B9)/G9)
+
+    den5 = sum(weights)
+
+    B5 = G5 * num_B5/den5
+
+    return B5

src/methods/antoine_vouillon/reconstruct.py

+"""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.antoine_vouillon.functions import *
+
+#!!!!!!!! It is normal that the reconstructions lasts several minutes (3min on my computer)
+
+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.
+    """
+
+    # Define constants and operators
+    cfa_name = 'bayer' # bayer or quad_bayer
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa_name, input_shape)
+
+    img_res = op.adjoint(y)
+    N = img_res[:,:,0].shape[0]
+    M = img_res[:,:,0].shape[1]
+
+    #interpolating green channel
+
+    for i in range (N):
+        for j in range (M):
+            if img_res[i,j,1] ==0:
+
+                neighbors = find_neighbors(img_res,1,i,j,N,M)
+                dir_deriv = find_dir_deriv(neighbors)
+                weights = find_weights(img_res, neighbors, dir_deriv,1,i,j,N,M)
+
+                img_res[i,j,1] = interpolate(neighbors,weights)
+                
+    
+    img_res[img_res>1] = 1
+    img_res[img_res<0] = 0
+
+    #first intepolation of red channel
+
+    for i in range (1,N,2):
+        for j in range (0,M,2):
+            
+            neighbors = find_neighbors(img_res,0,i,j,N,M)
+
+            neighbors_G = find_neighbors(img_res,1,i,j,N,M)
+            dir_deriv = find_dir_deriv(neighbors_G)
+            weights = find_weights(img_res,neighbors_G, dir_deriv,1,i,j,N,M) 
+
+
+
+            img_res[i,j,0] = interpolate_RB(neighbors, neighbors_G, weights)
+
+    # second interpolation of red channel
+
+    for i in range (N):
+        for j in range (M):
+            if img_res[i,j,0] ==0:
+
+                neighbors = find_neighbors(img_res,0,i,j,N,M)
+                dir_deriv = find_dir_deriv(neighbors)
+                weights = find_weights(img_res,neighbors, dir_deriv,0,i,j,N,M)
+
+                img_res[i,j,0] = interpolate(neighbors,weights)
+
+    img_res[img_res>1] = 1
+    img_res[img_res<0] = 0
+
+    #first interpolation of blue channel
+
+    for i in range (0,N,2):
+        for j in range (1,M,2):
+            
+            neighbors = find_neighbors(img_res,2,i,j,N,M)
+
+            neighbors_G = find_neighbors(img_res,1,i,j,N,M)
+            dir_deriv = find_dir_deriv(neighbors_G)
+            weights = find_weights(img_res,neighbors_G, dir_deriv,1,i,j,N,M) 
+
+
+
+            img_res[i,j,2] = interpolate_RB(neighbors, neighbors_G, weights)
+
+    #second interpolation of blue channel
+
+    for i in range (N):
+        for j in range (M):
+            if img_res[i,j,2] ==0:
+
+                neighbors = find_neighbors(img_res,2,i,j,N,M)
+                dir_deriv = find_dir_deriv(neighbors)
+                weights = find_weights(img_res, neighbors, dir_deriv,2,i,j,N,M)
+
+                img_res[i,j,2] = interpolate(neighbors,weights)
+
+    img_res[img_res>1] = 1
+    img_res[img_res<0] = 0
+
+    return img_res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/babacar/demoisaicing_fct.py

+"""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 High_Quality_Linear_Interpolation(op: CFA, y: np.ndarray) -> np.ndarray:
+    """Performs a High-Quality Linear Interpolation of the lost pixels.
+
+    Args:
+        op (CFA): CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    z = op.adjoint(y)
+    #operation for the bayer pattern
+    if op.cfa == 'bayer':
+        res = np.empty(op.input_shape)
+        mask = op.mask
+        R  = 0
+        G = 1
+        B = 2
+
+        y_padded = np.pad(y,2, 'constant', constant_values=0)
+
+        for i in range(2,y_padded.shape[0]-2):
+            for j in range(2,y_padded.shape[1]-2):
+                for c in range(op.input_shape[2]):
+                    i_real = i-2
+                    j_real = j-2
+                    if (c==R and mask[i_real,j_real,R] == 1) or (c==G and mask[i_real,j_real,G] == 1)or (c==B and mask[i_real,j_real,B] == 1): #no need to find the pixel value
+                        res[i_real,j_real,c] = y_padded[i,j]
+                    elif (c == G) and (mask[i_real,j_real,B] == 1 or mask[i_real,j_real,R] == 1): #green channel and pixel value = B or R
+                        res[i_real,j_real,c] = float(convolve2d(y_padded[i-2:i+3, j-2:j+3], hq_g_r_and_b, mode='valid'))
+
+                    elif (c == R or c== B )and (mask[i_real,j_real,G] == 1) and (mask[i_real,j_real-1,c] == 0): #red or blue channel and pixel value = G and left pixel isn't in the recovered layer
+                        res[i_real,j_real,c] = float(convolve2d(y_padded[i-2:i+3, j-2:j+3], hq_r_and_b_g_brow_rcol, mode='valid'))
+
+                    elif (c == R or c == B) and (mask[i_real,j_real,G] == 1) and (mask[i_real-1,j_real,c] == 0): #red or blue channel and pixel value = G and up pixel isn't in the recovered layer
+                        res[i_real,j_real,c] = float(convolve2d(y_padded[i-2:i+3, j-2:j+3], hq_r_and_b_g_rrow_bcol, mode='valid'))
+
+                    elif (c == R and mask[i_real,j_real,B] == 1) or (c == B and mask[i_real,j_real,R] == 1) : #red channel and pixel value = B or inverse
+                        res[i_real,j_real,c] = float(convolve2d(y_padded[i-2:i+3, j-2:j+3], hq_r_b, mode='valid'))
+                    else : print("error")
+
+    #operation for the quad bayer pattern
+    else:
+        res = np.empty(op.input_shape)
+        mask = op.mask
+        R  = 0
+        G = 1
+        B = 2
+
+        y_padded = np.pad(y,4, 'constant', constant_values=0)
+
+        for i in range(4,y_padded.shape[0]-4, 2):
+            for j in range(4,y_padded.shape[1]-4, 2):
+                for c in range(op.input_shape[2]):
+                    i_real = i-4
+                    j_real = j-4
+                    if (c==R and mask[i_real,j_real,R] == 1) or (c==G and mask[i_real,j_real,G] == 1)or (c==B and mask[i_real,j_real,B] == 1): #no need to find the pixel value
+                        res[i_real:i_real+2,j_real:j_real+2,c] = y_padded[i:i+2,j:j+2]
+                        
+                    elif (c == G) and (mask[i_real,j_real,B] == 1 or mask[i_real,j_real,R] == 1): #green channel and pixel value = B or R
+                        # res[i_real:i_real+2,j_real:j_real+2,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], hq_quad_g_r_and_b, mode='valid'))
+
+                        res[i_real,j_real,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], A_hq_quad_g_r_and_b, mode='valid'))
+                        res[i_real,j_real+1,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], B_hq_quad_g_r_and_b, mode='valid'))
+                        res[i_real+1,j_real,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], C_hq_quad_g_r_and_b, mode='valid'))
+                        res[i_real+1,j_real+1,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], D_hq_quad_g_r_and_b, mode='valid'))
+
+                    elif (c == R or c== B )and (mask[i_real,j_real,G] == 1) and (mask[i_real,j_real-2,c] == 0): #red or blue channel and pixel value = G and left pixel isn't in the recovered layer
+                        # res[i_real:i_real+2,j_real:j_real+2,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], hq_quad_r_and_b_g_brow_rcol, mode='valid'))
+
+
+                        res[i_real,j_real,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], A_hq_quad_r_and_b_g_brow_rcol, mode='valid'))
+                        res[i_real,j_real+1,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], B_hq_quad_r_and_b_g_brow_rcol, mode='valid'))
+                        res[i_real+1,j_real,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], C_hq_quad_r_and_b_g_brow_rcol, mode='valid'))
+                        res[i_real+1,j_real+1,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], D_hq_quad_r_and_b_g_brow_rcol, mode='valid'))
+
+                    elif (c == R or c == B) and (mask[i_real,j_real,G] == 1) and (mask[i_real-2,j_real,c] == 0): #red or blue channel and pixel value = G and up pixel isn't in the recovered layer
+                        # res[i_real:i_real+2,j_real:j_real+2,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], hq_quad_r_and_b_g_rrow_bcol, mode='valid'))
+
+                        res[i_real,j_real,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], A_hq_quad_r_and_b_g_rrow_bcol, mode='valid'))
+                        res[i_real,j_real+1,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], B_hq_quad_r_and_b_g_rrow_bcol, mode='valid'))
+                        res[i_real+1,j_real,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], C_hq_quad_r_and_b_g_rrow_bcol, mode='valid'))
+                        res[i_real+1,j_real+1,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], D_hq_quad_r_and_b_g_rrow_bcol, mode='valid'))
+                        
+
+                    elif (c == R and mask[i_real,j_real,B] == 1) or (c == B and mask[i_real,j_real,R] == 1) : #red channel and pixel value = B or inverse
+                        # res[i_real:i_real+2,j_real:j_real+2,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], hq_quad_r_b, mode='valid'))
+
+                        res[i_real,j_real,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], A_hq_quad_r_b, mode='valid'))
+                        res[i_real,j_real+1,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], B_hq_quad_r_b, mode='valid'))
+                        res[i_real+1,j_real,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], C_hq_quad_r_b, mode='valid'))
+                        res[i_real+1,j_real+1,c] = float(convolve2d(y_padded[i-4:i+6, j-4:j+6], D_hq_quad_r_b, mode='valid'))
+                    else : print("error")
+    np.clip(res, 0, 1, res)
+    return res
+
+
+# Bayer masks
+
+hq_g_r_and_b = 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
+
+hq_r_and_b_g_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
+
+hq_r_and_b_g_brow_rcol = hq_r_and_b_g_rrow_bcol.T
+
+hq_r_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
+
+
+# Quad bayer masks
+
+## Calcul the value of 4 pixels with one mask, then the value is apply for the 4 pixels
+
+
+hq_quad_g_r_and_b = 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]]) / 16
+
+hq_quad_r_and_b_g_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]]) / 16
+
+hq_quad_r_and_b_g_brow_rcol = hq_quad_r_and_b_g_rrow_bcol.T
+
+hq_quad_r_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]]) / 16
+
+
+## Calcul the value of one pixel of the 2x2 matrix independentely
+# A = up left pixel / B = up right pixel / C = down left pixel / D = down right pixel
+
+A_hq_quad_g_r_and_b = np.array([[0,0, 0,0, -1,0, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                              [0,0, 0,0, 2,0, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                              [-1,0, 2,0, 4,0, 2,0, -1,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                              [0,0, 0,0, 2,0, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                              [0,0, 0,0, -1,0, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0]]) / 8
+
+B_hq_quad_g_r_and_b = np.array([[0,0, 0,0, 0,-1, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                              [0,0, 0,0, 0,2, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                              [0,-1, 0,2, 0,4, 0,2, 0,-1],[0,0, 0,0, 0,0, 0,0, 0,0],
+                              [0,0, 0,0, 0,2, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                              [0,0, 0,0, 0,-1, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0]]) / 8
+
+C_hq_quad_g_r_and_b = np.array([[0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, -1,0, 0,0, 0,0],
+                              [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 2,0, 0,0, 0,0],
+                              [0,0, 0,0, 0,0, 0,0, 0,0],[-1,0, 2,0, 4,0, 2,0, -1,0],
+                              [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 2,0, 0,0, 0,0],
+                              [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, -1,0, 0,0, 0,0]]) / 8
+
+D_hq_quad_g_r_and_b = np.array([[0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 0,-1, 0,0, 0,0],
+                              [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 0,2, 0,0, 0,0],
+                              [0,0, 0,0, 0,0, 0,0, 0,0],[0,-1, 0,2, 0,4, 0,2, 0,-1],
+                              [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 0,2, 0,0, 0,0],
+                              [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 0,-1, 0,0, 0,0]]) / 8
+
+
+A_hq_quad_r_and_b_g_rrow_bcol = np.array([[0,0, 0,0, 1/2,0, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                                        [0,0, -1,0, 0,0, -1,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                                        [-1,0, 4,0, 5,0, 4,0, -1,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                                        [0,0, -1,0, 0,0, -1,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                                        [0,0, 0,0, 1/2,0, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0]]) / 8
+
+B_hq_quad_r_and_b_g_rrow_bcol = np.array([[0,0, 0,0, 0,1/2, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                                        [0,0, 0,-1, 0,0, 0,-1, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                                        [0,-1, 0,4, 0,5, 0,4, 0,-1],[0,0, 0,0, 0,0, 0,0, 0,0],
+                                        [0,0, 0,-1, 0,0, 0,-1, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                                        [0,0, 0,0, 0,-1/2, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0]]) / 8
+
+C_hq_quad_r_and_b_g_rrow_bcol = np.array([[0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 1/2,0, 0,0, 0,0],
+                                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, -1,0, 0,0, -1,0, 0,0],
+                                        [0,0, 0,0, 0,0, 0,0, 0,0],[-1,0, 4,0, 5,0, 4,0, -1,0],
+                                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, -1,0, 0,0, -1,0, 0,0],
+                                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 1/2,0, 0,0, 0,0]]) / 8
+
+D_hq_quad_r_and_b_g_rrow_bcol = np.array([[0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 0,1/2, 0,0, 0,0],
+                                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,-1, 0,0, 0,-1, 0,0],
+                                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,-1, 0,4, 0,5, 0,4, 0,-1],
+                                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,-1, 0,0, 0,-1, 0,0],
+                                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 0,-1/2, 0,0, 0,0]]) / 8
+
+
+A_hq_quad_r_and_b_g_brow_rcol = A_hq_quad_r_and_b_g_rrow_bcol.T
+B_hq_quad_r_and_b_g_brow_rcol = B_hq_quad_r_and_b_g_rrow_bcol.T
+C_hq_quad_r_and_b_g_brow_rcol = C_hq_quad_r_and_b_g_rrow_bcol.T
+D_hq_quad_r_and_b_g_brow_rcol = D_hq_quad_r_and_b_g_rrow_bcol.T
+
+
+
+
+A_hq_quad_r_b = np.array([[0,0, 0,0, -3/2,0, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                        [0,0, 2,0, 0,0, 2,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                        [-3/2,0, 0,0, 6,0, 0,0, -3/2,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                        [0,0, 2,0, 0,0, 2,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                        [0,0, 0,0, -3/2,0, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0]]) / 8
+
+B_hq_quad_r_b = np.array([[0,0, 0,0, 0,-3/2, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                        [0,0, 0,2, 0,0, 0,2, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                        [0,-3/2, 0,0, 0,6, 0,0, 0,-3/2],[0,0, 0,0, 0,0, 0,0, 0,0],
+                        [0,0, 0,2, 0,0, 0,2, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0],
+                        [0,0, 0,0, 0,-3/2, 0,0, 0,0],[0,0, 0,0, 0,0, 0,0, 0,0]]) / 8
+
+
+C_hq_quad_r_b = np.array([[0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, -3/2,0, 0,0, 0,0],
+                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 2,0, 0,0, 2,0, 0,0],
+                        [0,0, 0,0, 0,0, 0,0, 0,0],[-3/2,0, 0,0, 6,0, 0,0, -3/2,0],
+                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 2,0, 0,0, 2,0, 0,0],
+                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, -3/2,0, 0,0, 0,0]]) / 8
+
+D_hq_quad_r_b = np.array([[0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 0,-3/2, 0,0, 0,0],
+                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,2, 0,0, 0,2, 0,0],
+                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,-3/2, 0,0, 0,6, 0,0, 0,-3/2],
+                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,2, 0,0, 0,2, 0,0],
+                        [0,0, 0,0, 0,0, 0,0, 0,0],[0,0, 0,0, 0,-3/2, 0,0, 0,0]]) / 8
+
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/babacar/reconstruct.py

+"""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
+import src.methods.babacar.demoisaicing_fct as b
+
+
+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.
+
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+    res  = b.High_Quality_Linear_Interpolation(op, y)
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/bach_antoine/MLRI_reconstruction.py

+"""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 src.forward_model import CFA
+import cv2
+
+def MLRI_interpolation(op: CFA, y: np.ndarray) -> np.ndarray:
+    z = op.adjoint(y)
+    res = np.empty(op.input_shape)
+    if op.cfa == 'bayer':
+        green_estimation=green_estimate(y, z)        
+        red_estimation=red_estimate(z,green_estimation)
+        blue_estimation=blue_estimate(z,green_estimation)
+        res[:, :, 0] = red_estimation
+        res[:, :, 1] = green_estimation
+        res[:, :, 2] = blue_estimation
+    else :
+        new_z = cv2.resize(z, (z.shape[1] // 2, z.shape[0] // 2), interpolation=cv2.INTER_AREA)
+        new_y=np.sum(new_z, axis=2)
+        new_z_upsampled = np.repeat(np.repeat(new_z, 2, axis=0), 2, axis=1)
+        new_y_upsampled=np.sum(new_z_upsampled, axis=2)
+        mask_out=cv2.merge([y-new_y_upsampled]*3)
+
+        res_downsampled = np.zeros_like(new_z)
+        green_estimation=green_estimate(new_y,new_z)        
+        red_estimation=red_estimate(new_z,green_estimation)
+        blue_estimation=blue_estimate(new_z,green_estimation)
+        res_downsampled[:, :, 0] = red_estimation
+        res_downsampled[:, :, 1] = green_estimation
+        res_downsampled[:, :, 2] = blue_estimation
+        res_downsampled_upsampled= np.repeat(np.repeat(res_downsampled, 2, axis=0), 2, axis=1)
+        res=res_downsampled_upsampled+mask_out
+    return res
+
+
+    
+def green_estimate(y : np.ndarray, z : np.ndarray)-> np.ndarray:
+    # estimation of green by gradient based threshold free color filter array interpolation(GBTF)
+    (grh_est, grv_est, 
+     gbh_est, gbv_est, 
+     rgh_est, rgv_est, 
+     bgh_est, bgv_est)=hamilton_and_adam_interpolation(y)
+    
+    (HCDE,VCDE)=color_difference(grh_est, grv_est , 
+                     gbh_est , gbv_est , 
+                     rgh_est , rgv_est ,
+                     bgh_est , bgv_est, z)
+    
+    (gradient_H,gradient_V)=gradient_compute(HCDE, VCDE, y)
+    Weight_W, Weight_E, Weight_S, Weight_N=directionnal_weight(gradient_H, gradient_V)
+    green_estimation=final_estimation(Weight_W, Weight_E, Weight_S, Weight_N,
+                     HCDE, VCDE,y, z)
+    return green_estimation  
+
+def red_estimate(z : np.ndarray, green_estimation : np.ndarray)-> np.ndarray:
+    # estimation of red by minimizing laplacian residuals interpolation (MLRI)
+    return red_blue_estimate(z, green_estimation, 0)
+
+def blue_estimate(z : np.ndarray, green_estimation : np.ndarray)-> np.ndarray:
+    # estimation of blue by minimizing laplacian residuals interpolation (MLRI)
+    return red_blue_estimate(z, green_estimation, 2)
+
+def red_blue_estimate(z : np.ndarray, green_estimation : np.ndarray, blue_red)-> np.ndarray:
+    
+    # estimation of blue or red by minimizing laplacian residuals interpolation
+    
+    # laplacian interpolation
+    F = np.array([
+        [0, 0, -1, 0, 0],
+        [0, 0, 0, 0, 0],
+        [-1, 0, 4, 0, -1],
+        [0, 0, 0, 0, 0],
+        [0, 0, -1, 0, 0]
+    ], dtype=np.float32)
+
+    lap_rb = cv2.filter2D(z[:,:,blue_red], -1, F, borderType=cv2.BORDER_REPLICATE)
+    lap_green = cv2.filter2D(green_estimation * (z[:,:,blue_red]!=0), -1, F, borderType=cv2.BORDER_REPLICATE)
+    
+    #estimate the residuals
+    mean_a,mean_b = guided_filter(green_estimation*(z[:,:,blue_red]!=0), z[:,:,blue_red], 
+                                  lap_green*(z[:,:,blue_red]!=0) , lap_rb*(z[:,:,blue_red]!=0))
+
+    tentativeRB = mean_a * green_estimation + mean_b
+
+    residualRB = (z[:,:,blue_red]!=0) * (z[:,:,blue_red] - tentativeRB)
+
+    # Bilinear interpolation of the residuals
+    H = np.array([
+        [1/4, 1/2, 1/4],
+        [1/2, 1, 1/2],
+        [1/4, 1/2, 1/4]
+    ], dtype=np.float32)
+
+    residualRB2 = cv2.filter2D(residualRB, -1, H, borderType=cv2.BORDER_REPLICATE)
+    rb_estimation = residualRB2 + tentativeRB
+    return rb_estimation
+
+def hamilton_and_adam_interpolation(y : np.ndarray)-> (np.ndarray,np.ndarray,
+                                                       np.ndarray,np.ndarray,
+                                                       np.ndarray,np.ndarray,
+                                                       np.ndarray,np.ndarray):
+    
+    # interpolation method to all pixels in both vertical and horizontal directions.
+    grh_est=np.zeros_like(y)
+    grv_est=np.zeros_like(y)
+    gbh_est=np.zeros_like(y)
+    gbv_est=np.zeros_like(y)
+
+    rgh_est=np.zeros_like(y)
+    rgv_est=np.zeros_like(y)
+    bgh_est=np.zeros_like(y)
+    bgv_est=np.zeros_like(y)
+    rows,cols=y.shape
+    for i in range(rows):
+        for j in range(cols):
+            if (i%2==0 and j%2!=0):
+                # estimate horizontal green with red
+                if (j==1):
+                    grh_est[i,j] =0
+                elif (j==cols-2) or (j==cols-1):
+                    grh_est[i,j] = 2*y[i,j-2]-y[i,j-4]
+                else:
+                    grh_est[i,j] = (y[i, j-1] + y[i, j+1]) / 2 + (2 * y[i, j] - y[i, j-2] - y[i, j+2]) / 4
+                # estimate vertical green with red
+                if (i==0):
+                    grv_est[i,j] = 0
+                elif (i==rows-2)or (i==rows-1):
+                    grv_est[i,j] = 2*y[i-2,j]-y[i-4,j]
+                else:
+                    grv_est[i,j] = (y[i-1, j] + y[i+1, j]) / 2 + (2 * y[i, j] - y[i-2, j] - y[i+2, j]) / 4
+        
+            if (i%2!=0 and j%2==0):
+                # estimate horizontal green with blue
+                if (j==0):
+                    gbh_est[i,j] =0
+                elif (j==cols-2) or (j==cols-1):
+                    gbh_est[i,j] = 2*y[i,j-2]-y[i,j-4]
+                else:
+                    gbh_est[i,j] = (y[i, j-1] + y[i, j+1]) / 2 + (2 * y[i, j] - y[i, j-2] - y[i, j+2]) / 4
+                # estimate vertical green with blue
+                if (i==1):
+                    gbv_est[i,j] = 0
+                elif (i==rows-2)or (i==rows-1):
+                    gbv_est[i,j] = 2*y[i-2,j]-y[i-4,j]
+                else:
+                    gbv_est[i,j] = (y[i-1, j] + y[i+1, j]) / 2 + (2 * y[i, j] - y[i-2, j] - y[i+2, j]) / 4
+            
+            if (i%2==0 and j%2==0):    
+                # estimate horizontal red with green
+                if (j==0):
+                    rgh_est[i,j] = 0
+                elif (j==cols-2) or (j==cols-1):
+                    rgh_est[i,j] = 2*y[i,j-2]-y[i,j-4]
+                else:
+                    rgh_est[i,j] = (y[i, j-1] + y[i, j+1]) / 2 + (2 * y[i, j] - y[i, j-2] - y[i, j+2]) / 4
+                # estimate vertical blue with green
+                if (i==0):
+                    bgv_est[i,j] = 0
+                elif (i==rows-2)or (i==rows-1):
+                    bgv_est[i,j] = 2*y[i-2,j]-y[i-4,j]
+                else:
+                    bgv_est[i,j] = (y[i-1, j] + y[i+1, j]) / 2 + (2 * y[i, j] - y[i-2, j] - y[i+2, j]) / 4
+        
+            if (i%2!=0 and j%2!=0):  
+                # estimate vertical red with green
+                if (i==1):
+                    rgv_est[i,j] = 0
+                elif (i==rows-2)or (i==rows-1):
+                    rgv_est[i,j] = 2*y[i-2,j]-y[i-4,j]
+                else:
+                    rgv_est[i,j] = (y[i-1, j] + y[i+1, j]) / 2 + (2 * y[i, j] - y[i-2, j] - y[i+2, j]) / 4
+                # estimate horizontal blue with green
+                if (j==1):
+                    bgh_est[i,j] = 0
+                elif (j==cols-2) or (j==cols-1):
+                    bgh_est[i,j] = 2*y[i,j-2]-y[i,j-4]
+                else:
+                    bgh_est[i,j] = (y[i, j-1] + y[i, j+1]) / 2 + (2 * y[i, j] - y[i, j-2] - y[i, j+2]) / 4 
+                    
+    grh_est[:,1] = 2*grh_est[:,3]-grh_est[:,5]
+    grv_est[0,:] = 2*grv_est[2,:]-grv_est[4,:]
+    gbh_est[:,0] = 2*gbh_est[:,2]-gbh_est[:,4]
+    gbv_est[1,:] = 2*gbv_est[3,:]-gbv_est[5,:]
+
+    rgh_est[:,0] = 2*rgh_est[:,2]-rgh_est[:,4]
+    rgv_est[1,:] = 2*rgv_est[3,:]-rgv_est[5,:]
+    bgh_est[:,1] = 2*bgh_est[:,3]-bgh_est[:,5]
+    bgv_est[0,:] = 2*bgv_est[2,:]-bgv_est[4,:]
+    
+    return grh_est, grv_est, gbh_est, gbv_est, rgh_est, rgv_est, bgh_est, bgv_est
+
+def color_difference(grh_est : np.ndarray, grv_est : np.ndarray, 
+                     gbh_est : np.ndarray, gbv_est : np.ndarray, 
+                     rgh_est : np.ndarray, rgv_est : np.ndarray,
+                     bgh_est : np.ndarray, bgv_est : np.ndarray,
+                     z : np.ndarray)-> (np.ndarray, np.ndarray) :
+    
+    # estimate the color difference array in horizontal and vertical directions
+    HCDE=z[:,:,1]-z[:,:,0]-z[:,:,2]
+    VCDE=HCDE.copy()
+    difference_rg_h=grh_est-rgh_est
+    difference_rg_v=grv_est-rgv_est
+    difference_bg_h=gbh_est-bgh_est
+    difference_bg_v=gbv_est-bgv_est
+    HCDE=HCDE+difference_rg_h+difference_bg_h
+    VCDE=VCDE+difference_rg_v+difference_bg_v
+    
+    return HCDE,VCDE
+def gradient_compute(HCDE : np.ndarray, VCDE : np.ndarray, y: np.ndarray)-> (np.ndarray, np.ndarray) :
+    
+    # estimate the color difference gradient array in horizontal and vertical directions
+    new_gradient_vcde= cv2.copyMakeBorder(VCDE, 1,1,1,1, cv2.BORDER_CONSTANT,None,0)
+    new_gradient_hcde= cv2.copyMakeBorder(HCDE, 1,1,1,1, cv2.BORDER_CONSTANT,None,0)
+    gradient_V=np.zeros_like(y)
+    gradient_H=np.zeros_like(y)
+    rows,cols=y.shape
+    for i in range(rows):
+        for j in range(cols):
+            gradient_V[i,j]=np.abs(new_gradient_vcde[i,j+1]-new_gradient_vcde[i+2,j+1])
+            gradient_H[i,j]=np.abs(new_gradient_hcde[i+1,j]-new_gradient_hcde[i+1,j+2])
+    return gradient_H,gradient_V
+
+def directionnal_weight(gradient_H : np.ndarray, gradient_V : np.ndarray)-> (np.ndarray, np.ndarray, 
+                                                                             np.ndarray, np.ndarray):
+
+    # Calculate the weights for each direction by computing the reciprocal of power gradients 
+    # in this direction within a local window 3*3
+    Kernel_Weight = np.ones((3, 3), dtype=np.float32)
+    Weight_H = cv2.filter2D(gradient_H, -1, Kernel_Weight, borderType=cv2.BORDER_REPLICATE)
+    Weight_V = cv2.filter2D(gradient_V, -1, Kernel_Weight, borderType=cv2.BORDER_REPLICATE)
+
+    Weight_W=cv2.copyMakeBorder(Weight_H,0,0,1,0, cv2.BORDER_REPLICATE,None,0)[:,:-1]
+    Weight_E=cv2.copyMakeBorder(Weight_H,0,0,0,1, cv2.BORDER_REPLICATE,None,0)[:,1:]
+    Weight_S=cv2.copyMakeBorder(Weight_V,0,1,0,0, cv2.BORDER_REPLICATE,None,0)[1:,:]
+    Weight_N=cv2.copyMakeBorder(Weight_V,1,0,0,0, cv2.BORDER_REPLICATE,None,0)[:-1,:]
+
+
+    Weight_W = 1.0 / (np.power(Weight_W, 2))
+    Weight_E = 1.0 / (np.power(Weight_E, 2))
+    Weight_S = 1.0 / (np.power(Weight_S, 2))
+    Weight_N = 1.0 / (np.power(Weight_N, 2))
+    return Weight_W, Weight_E, Weight_S, Weight_N
+
+def final_estimation(Weight_W : np.ndarray, Weight_E : np.ndarray, 
+                     Weight_S : np.ndarray, Weight_N : np.ndarray,
+                     HCDE : np.ndarray, VCDE : np.ndarray,
+                     y : np.ndarray, z : np.ndarray, 
+                     size=9, sigma=1)-> np.ndarray:
+    
+    # directional color differences and weight is computing together to estimate green
+    h = cv2.getGaussianKernel(size, sigma)
+    Ke = np.array([0, 0, 0, 0, 1, 1, 1, 1, 1], dtype=np.float32) * h.T
+    Kw = np.array([1, 1, 1, 1, 1, 0, 0, 0, 0], dtype=np.float32) * h.T
+    Ke /= np.sum(Ke, axis=1)
+    Kw /= np.sum(Kw, axis=1)
+    Ks = np.transpose(Ke)
+    Kn = np.transpose(Kw)
+
+    difn = cv2.filter2D(VCDE, -1, Kn, borderType=cv2.BORDER_REPLICATE)
+    difs = cv2.filter2D(VCDE, -1, Ks, borderType=cv2.BORDER_REPLICATE)
+    difw = cv2.filter2D(HCDE, -1, Kw, borderType=cv2.BORDER_REPLICATE)
+    dife = cv2.filter2D(HCDE, -1, Ke, borderType=cv2.BORDER_REPLICATE)
+
+    Wt = Weight_W + Weight_E + Weight_N + Weight_S
+    final_estimation = (Weight_N * difn + Weight_S * difs + Weight_W * difw + Weight_E * dife) / Wt
+    green = final_estimation + y
+    green_estimate = green * (z[:,:,1]==0) + z[:,:,1]
+    return green_estimate
+
+def guided_filter(guide : np.ndarray, to_interpolate : np.ndarray, 
+                  guide_laplacian : np.ndarray, to_interpolate_laplacian : np.ndarray, 
+                  r=11, eps=0)-> (np.ndarray, np.ndarray) :
+    
+    # interpolate the red or blue value with a guided filter create with the green estimation
+    kernel = np.ones((r, r), dtype=np.float32)
+    mask_normalize = cv2.filter2D(np.ones_like(guide)*(to_interpolate!=0), -1, kernel, borderType=cv2.BORDER_CONSTANT)
+    mask_normalize_ab = cv2.filter2D(np.ones_like(guide), -1, kernel, borderType=cv2.BORDER_CONSTANT)
+
+    mean_p = cv2.filter2D(to_interpolate_laplacian, -1, kernel, borderType=cv2.BORDER_CONSTANT)/mask_normalize
+    mean_I = cv2.filter2D(guide_laplacian, -1, kernel, borderType=cv2.BORDER_CONSTANT)/mask_normalize
+    corr_I = cv2.filter2D(guide_laplacian*guide_laplacian, -1, kernel, borderType=cv2.BORDER_CONSTANT)/mask_normalize
+    corr_Ip = cv2.filter2D(guide_laplacian*to_interpolate_laplacian, -1, kernel, borderType=cv2.BORDER_CONSTANT)/mask_normalize
+
+    var_I = corr_I - mean_I*mean_I
+    cov_Ip = corr_Ip - mean_I*mean_p
+
+    a = cov_Ip / (var_I + eps)
+
+    mean_guide = cv2.filter2D(guide, -1, kernel, borderType=cv2.BORDER_CONSTANT)/mask_normalize
+    mean_to_interpolate = cv2.filter2D(to_interpolate, -1, kernel, borderType=cv2.BORDER_CONSTANT)/mask_normalize
+
+    b = mean_to_interpolate - a * mean_guide
+    mean_a = cv2.filter2D(a, -1, kernel, borderType=cv2.BORDER_CONSTANT)/mask_normalize_ab
+    mean_b = cv2.filter2D(b, -1, kernel, borderType=cv2.BORDER_CONSTANT)/mask_normalize_ab
+    return mean_a,mean_b
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 28/01/2024
+# Authors: BACH Antoine

src/methods/bach_antoine/reconstruct.py

+"""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.bach_antoine.MLRI_reconstruction import MLRI_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 = MLRI_interpolation(op, y)
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 28/01/2024
+# Authors: BACH Antoine
+

src/methods/baseline/demo_reconstruction.py

+"""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 naive_interpolation(op: CFA, y: np.ndarray) -> np.ndarray:
+    """Performs a simple interpolation of the lost pixels.
+
+    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)
+
+        res[:, :, 0] = convolve2d(z[:, :, 0], ker_bayer_red_blue, mode='same')
+        res[:, :, 1] = convolve2d(z[:, :, 1], ker_bayer_green, mode='same')
+        res[:, :, 2] = convolve2d(z[:, :, 2], ker_bayer_red_blue, mode='same')
+
+    else:
+        res = np.empty(op.input_shape)
+
+        res[:, :, 0] = varying_kernel_convolution(z[:, :, 0], K_list_red)
+        res[:, :, 1] = varying_kernel_convolution(z[:, :, 1], K_list_green)
+        res[:, :, 2] = varying_kernel_convolution(z[:, :, 2], K_list_blue)
+
+    return res
+
+
+def extract_padded(M, size, i, j):
+    N_i, N_j = M.shape
+    res = np.zeros((size, size))
+    middle_size = int((size - 1) / 2)
+
+    for ii in range(- middle_size, middle_size + 1):
+        for jj in range(- middle_size, middle_size + 1):
+            if i + ii >= 0 and i + ii < N_i and j + jj >= 0 and j + jj < N_j:
+                res[middle_size + ii, middle_size + jj] = M[i + ii, j + jj]
+
+    return res
+
+
+def varying_kernel_convolution(M, K_list):
+    N_i, N_j = M.shape
+    res = np.zeros_like(M)
+
+    for i in range(N_i):
+        for j in range(N_j):
+            res[i, j] = np.sum(extract_padded(M, K_list[4 * (i % 4) + j % 4].shape[0], i, j) * K_list[4 * (i % 4) + j % 4])
+
+    np.clip(res, 0, 1, res)
+
+    return res
+
+
+K_identity = np.zeros((5, 5))
+K_identity[2, 2] = 1
+
+K_red_0 = np.zeros((5, 5))
+K_red_0[2, :] = np.array([-3, 13, 0, 0, 2]) / 12
+
+K_red_1 = np.zeros((5, 5))
+K_red_1[2, :] = np.array([2, 0, 0, 13, -3]) / 12
+
+K_red_8 = np.zeros((5, 5))
+K_red_8[:2, :2] = np.array([[-1, -1], [-1, 9]]) / 6
+
+K_red_9 = np.zeros((5, 5))
+K_red_9[:2, 3:] = np.array([[-1, -1], [9, -1]]) / 6
+
+K_red_10 = np.zeros((5, 5))
+K_red_10[:, 2] = np.array([-3, 13, 0, 0, 2]) / 12
+
+K_red_12 = np.zeros((5, 5))
+K_red_12[3:, :2] = np.array([[-1, 9], [-1, -1]]) / 6
+
+K_red_13 = np.zeros((5, 5))
+K_red_13[3:, 3:] = np.array([[9, -1], [-1, -1]]) / 6
+
+K_red_14 = np.zeros((5, 5))
+K_red_14[:, 2] = np.array([2, 0, 0, 13, -3]) / 12
+
+K_list_red = [K_red_0, K_red_1, K_identity, K_identity, K_red_0, K_red_1, K_identity, K_identity, K_red_8, K_red_9, K_red_10, K_red_10, K_red_12, K_red_13, K_red_14, K_red_14]
+
+
+K_green_2 = np.zeros((5, 5))
+K_green_2[2, :] = [-3, 13, 0, 0, 2]
+K_green_2[:, 2] = [-3, 13, 0, 0, 2]
+K_green_2 = K_green_2 / 24
+
+K_green_3 = np.zeros((5, 5))
+K_green_3[2, :] = [2, 0, 0, 13, -3]
+K_green_3[:, 2] = [-3, 13, 0, 0, 2]
+K_green_3 = K_green_3 / 24
+
+K_green_6 = np.zeros((5, 5))
+K_green_6[2, :] = [-3, 13, 0, 0, 2]
+K_green_6[:, 2] = [2, 0, 0, 13, -3]
+K_green_6 = K_green_6 / 24
+
+K_green_7 = np.zeros((5, 5))
+K_green_7[2, :] = [2, 0, 0, 13, -3]
+K_green_7[:, 2] = [2, 0, 0, 13, -3]
+K_green_7 = K_green_7 / 24
+
+K_list_green = [K_identity, K_identity, K_green_2, K_green_3, K_identity, K_identity, K_green_6, K_green_7, K_green_2, K_green_3, K_identity, K_identity, K_green_6, K_green_7, K_identity, K_identity]
+
+
+K_list_blue = [K_red_10, K_red_10, K_red_8, K_red_9, K_red_14, K_red_14, K_red_12, K_red_13, K_identity, K_identity, K_red_0, K_red_1, K_identity, K_identity, K_red_0, K_red_1]
+
+
+
+ker_bayer_red_blue = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 4
+
+ker_bayer_green = np.array([[0, 1, 0], [1, 4, 1], [0, 1, 0]]) / 4
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/baseline/reconstruct.py

+"""The main file for the baseline reconstruction.
+This file should NOT be modified.
+"""
+
+
+import numpy as np
+
+from src.forward_model import CFA
+from src.methods.baseline.demo_reconstruction import naive_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.
+    """
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+
+    res = naive_interpolation(op, y)
+
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/brice_convers/README.md

+# Image Analysis Project
+
+Numerous RGB cameras in the commercial sector employ Color Filter Array (CFA) technology.
+This technology involves an array of red, green, or blue filters placed atop the sensors, usually
+organized in periodic patterns. The incident light is consequently filtered by each filter, before
+being captured by the sensor. The typical process of acquiring images utilizes a predefined CFA
+pattern to allocate a color to each pixel of the sensor. The goal of this code project is to demosaicing these CFA Image.
+
+You can check the report called **Image_Analysis_Project_Report_Brice_Convers** to find more information about it.
+
+## How to use
+
+To use the code project you can move the **main_template.py** outside the **src** file and execute it.
+Or you can also call the **run_reconstruction** method in you programme.
+
+```Python
+
+import src.methods.brice_convers.dataHandler as DataHandler
+import src.methods.brice_convers.dataEvaluation as DataEvaluation
+import time
+
+WORKING_DIRECOTRY_PATH = "SICOM_Image_Analysis/sicom_image_analysis_project/"
+
+DataHandler = DataHandler.DataHandler(WORKING_DIRECOTRY_PATH)
+DataEvaluation = DataEvaluation.DataEvaluation(DataHandler)
+
+def main(DataHandler):
+
+    IMAGE_PATH = WORKING_DIRECOTRY_PATH + "images/"
+
+    CFA_NAME = "quad_bayer"
+
+    METHOD = "menon"
+
+    startTime = time.time()
+
+    DataHandler.list_images(IMAGE_PATH)
+
+    DataHandler.print_list_images()
+
+    DataHandler.compute_CFA_images(CFA_NAME)
+
+    DataHandler.compute_reconstruction_images(METHOD, {"cfa": CFA_NAME})
+
+    #The first agurment (ex: 3) is the image index in the list print by "DataHandler.print_list_images()"
+    DataHandler.plot_reconstructed_image(3, METHOD, {"cfa": CFA_NAME}, zoomSize="large")
+
+    DataEvaluation.print_metrics(3, METHOD)
+
+    endTime = time.time()
+
+    print("[INFO] Elapsed time: " + str(endTime - startTime) + "s")
+
+    print("[INFO] End")
+
+if __name__ == "__main__":
+    main(DataHandler)
+
+```
+
+## TODO List:
+
+- Fix menon method pour quad bayer pattern with landscape picture
+
+
+
+
+
+
+
+## References:
+
+[1] [*Research Paper:*](https://ieeexplore.ieee.org/document/4032820) Used for Menon Method.
+
+## Authors
+- [Brice Convers](https://briceconvers.com)

src/methods/brice_convers/configuration.py

+
+PIXEL_PATTERN = "RGB"
+REFINING_STEP = True