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/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 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
+
+
+
+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 = new_interpolation(op, y)
+
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

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/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/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

src/methods/brice_convers/dataEvaluation.py

+from src.methods.brice_convers.dataHandler import DataHandler
+from src.utils import psnr, ssim
+from sklearn.metrics import f1_score, mean_squared_error
+import numpy as np
+
+class DataEvaluation:
+    def __init__(self, DataHandler: DataHandler):
+        DataEvaluation.DataHandler = DataHandler
+
+    def print_metrics(self, indexImage, method):
+        DataEvaluation.DataHandler.indexImageExists(indexImage)
+
+        img = DataEvaluation.DataHandler.load_image(indexImage)
+        res = DataEvaluation.DataHandler.get_reconstructed_image(indexImage, method)
+
+        ssimMetric = ssim(img, res)
+        psnrMetrc = psnr(img, res)
+        mse = mean_squared_error(img.flatten(), res.flatten())
+        mseRedPixels = mean_squared_error(img[:,:,0], res[:,:,0])
+        mseGreenPixels = mean_squared_error(img[:,:,1], res[:,:,1])
+        mseBluePixels = mean_squared_error(img[:,:,2], res[:,:,2])
+        miMetric = DataEvaluation.MI(img, res)
+        ccMetric = DataEvaluation.CC(img, res)
+        sadMetric = DataEvaluation.SAD(img, res)
+        lsMetric = DataEvaluation.LS(img, res)
+
+        print("[INFO] Metrics for image {}".format(indexImage))
+        print("#" * 30)
+        print(" SSIM: {:.6}  ".format(ssimMetric))
+        print(" PSNR: {:.6}   ".format(psnrMetrc))
+        print(" MSE : {:.3e}    ".format(mse))
+        print(" MSE (R): {:.3e}    ".format(mseRedPixels))
+        print(" MSE (G): {:.3e}    ".format(mseGreenPixels))
+        print(" MSE (B): {:.3e}    ".format(mseBluePixels))
+        print(" MI: {:.6}     ".format(miMetric))
+        print(" CC: {:.6}    ".format(ccMetric))
+        print(" SAD: {:.6} ".format(sadMetric))
+        print(" LS: {:.3e} ".format(lsMetric))
+        print("#" * 30)
+
+    #Mutual Information
+    def MI(img_mov, img_ref):
+        hgram, x_edges, y_edges = np.histogram2d(img_mov.ravel(), img_ref.ravel(), bins=20)
+        pxy = hgram / float(np.sum(hgram))
+        px = np.sum(pxy, axis=1) # marginal for x over y
+        py = np.sum(pxy, axis=0) # marginal for y over x
+        px_py = px[:, None] * py[None, :] # Broadcast to multiply marginals
+        # Now we can do the calculation using the pxy, px_py 2D arrays
+        nzs = pxy > 0 # Only non-zero pxy values contribute to the sum
+        return np.sum(pxy[nzs] * np.log(pxy[nzs] / px_py[nzs]))
+
+    # Cross Correlation
+    def CC(img_mov, img_ref):
+        # Vectorized versions of c,d,e
+        a = img_mov.astype('float64')
+        b = img_ref.astype('float64')
+
+        # Calculating mean values
+        AM = np.mean(a)
+        BM = np.mean(b)
+
+        c_vect = (a - AM) * (b - BM)
+        d_vect = (a - AM) ** 2
+        e_vect = (b - BM) ** 2
+
+        # Finally get r using those vectorized versions
+        r_out = np.sum(c_vect) / float(np.sqrt(np.sum(d_vect) * np.sum(e_vect)))
+        return r_out
+
+    #Sum of Absolute Differences
+    def SAD(img_mov, img_ref):
+        img1 = img_mov.astype('float64')
+        img2 = img_ref.astype('float64')
+        ab = np.abs(img1 - img2)
+        sav = np.sum(ab.ravel())
+        sav /= ab.ravel().shape[0]
+        return sav
+
+    #Sum of Least Squared Errors
+    def LS(img_mov, img_ref):
+        img1 = img_mov.astype('float64')
+        img2 = img_ref.astype('float64')
+        r = (img1 - img2)**2
+        sse = np.sum(r.ravel())
+        sse /= r.ravel().shape[0]
+        return sse
+

src/methods/brice_convers/dataHandler.py

+import os
+import matplotlib.pyplot as plt
+import matplotlib.image as mpimg
+from src.forward_model import CFA
+from pathlib import Path
+from src.methods.baseline.reconstruct import run_reconstruction
+from src.methods.brice_convers.reconstruct import run_reconstruction as run_reconstruction_brice_convers
+from src.methods.brice_convers.utilities import folderExists
+from src.utils import normalise_image, save_image, psnr, ssim
+from skimage.io import imread
+import numpy as np
+
+class DataHandler:
+    def __init__(self, workingDirectoryPath = ""):
+        DataHandler.imagePaths =[]
+        DataHandler.imagePathsLabels = []
+        DataHandler.CFA_images = {}
+        DataHandler.reconstruction_images = {"interpolation": {}, "menon": {}}
+
+        DataHandler.IMAGE_TYPES = (".jpg", ".jpeg", ".png", ".bmp", ".tif", ".tiff")
+        DataHandler.WD_PATH = workingDirectoryPath
+
+    def list_files(self, basePath, interval, validExts=None, contains=None):
+        sum = 0
+
+        if interval is not None:
+            counter = interval[0]
+            lastFile = interval[1]
+        else:
+            counter = 0
+            lastFile = -1
+
+        # loop over the directory structure
+        for (rootDir, dirNames, filenames) in os.walk(basePath):
+            # loop over the filenames in the current directory
+            for filename in filenames:
+                # if the contains string is not none and the filename does not contain
+                # the supplied string, then ignore the file
+                if contains is not None and filename.find(contains) == -1:
+                    continue
+
+                # determine the file extension of the current file
+                ext = filename[filename.rfind("."):].lower()
+
+                # check to see if the file is an image and should be processed
+                if validExts is None or ext.endswith(validExts):
+                    if sum >= counter and (sum <= lastFile) or lastFile == -1:
+                        # construct the path to the image and yield it
+                        imagePath = os.path.join(rootDir, filename)
+                        imageName = Path(imagePath).stem
+                        DataHandler.imagePaths.append(imagePath)
+                        DataHandler.imagePathsLabels.append(imageName)
+
+                    sum += 1
+
+    def list_images(self, basePath, interval = None, contains=None):
+        # return the set of files that are valid
+        DataHandler.list_files(self, basePath, interval, validExts=DataHandler.IMAGE_TYPES, contains=contains)
+        print("[INFO] There are {} images.".format(len(DataHandler.imagePaths)))
+
+    def print_list_images(self):
+        print("[INFO] This is the order or your image in the list. IndexImage is the index in this table:")
+        print(DataHandler.imagePathsLabels)
+
+    def plot_input_images(self, rows=3, cols=3, figsize=(15, 5), max_images = 6, title="Input Images"):
+        fig = plt.figure(figsize=figsize)
+        fig.suptitle(title, fontsize=16)
+
+        for i, imagePath in enumerate(DataHandler.imagePaths):
+            if i >= max_images:
+                break
+
+            img = mpimg.imread(imagePath)
+            fig.add_subplot(rows, cols, i+1)
+            # image title
+            plt.title(imagePath.split(os.path.sep)[-1])
+            plt.text(0, -0.1, f"Shape of the image: {img.shape}.", fontsize=12, transform=plt.gca().transAxes)
+            plt.imshow(img)
+
+        plt.show()
+
+    def plot_raw_transformation(self, zoom = True, specificIndex = None, rows=8, cols=3, figsize=(15, 5), max_images = 24, title="Raw Transformation"):
+        fig, axs = plt.subplots(rows, cols, figsize=figsize)
+        fig.suptitle(title, fontsize=16)
+
+        for i, imagePath in enumerate(DataHandler.imagePaths):
+            if i >= max_images and max_images != -1:
+                break
+
+            if specificIndex is not None:
+                if i != specificIndex:
+                    continue
+
+            nameImage = Path(imagePath).stem
+
+            if DataHandler.CFA_images.get(nameImage) is None:
+                print("[ERROR] There is no CFA image for the image {}.".format(nameImage))
+                continue
+
+            img = mpimg.imread(imagePath)
+
+            y = DataHandler.CFA_images[nameImage].direct(img)
+            z = DataHandler.CFA_images[nameImage].adjoint(y)
+
+            if specificIndex is None:
+                line = 2*i
+                subLine = 2*i+1
+            else:
+                line = 0
+                subLine = 1
+
+            axs[line, 0].imshow(img)
+            axs[line, 0].set_title('Input image')
+            axs[line, 1].imshow(y, cmap='gray')
+            axs[line, 1].set_title('Output image')
+            axs[line, 2].imshow(z)
+            axs[line, 2].set_title('Adjoint image')
+
+            if zoom:
+                axs[subLine, 0].imshow(img[800:864, 450:514])
+                axs[subLine, 0].set_title('Zoomed input image')
+                axs[subLine, 1].imshow(y[800:864, 450:514], cmap='gray')
+                axs[subLine, 1].set_title('Zoomed output image')
+                axs[subLine, 2].imshow(z[800:864, 450:514])
+                axs[subLine, 2].set_title('Zoomed adjoint image')
+
+        print("[INFO] There are {} ploted images.".format(max_images))
+
+    def plot_specific_raw_transformation(self, indexImage, zoom = True, rows=2, cols=3, figsize=(15, 5), title="Raw Transformation"):
+        DataHandler.plot_raw_transformation(self, zoom, indexImage, rows, cols, figsize, -1, title)
+
+    def compute_CFA_images(self, CFA_NAME):
+        for i, imagePath in enumerate(DataHandler.imagePaths):
+            nameImage = Path(imagePath).stem
+            DataHandler.compute_CFA_image(self, CFA_NAME, nameImage, i)
+
+        print("[INFO] There are {} CFA images.".format(len(DataHandler.CFA_images)))
+
+    def compute_reconstruction_image(self, method, indexImage, options = None):
+        if len(DataHandler.imagePaths) == 0:
+            print("[ERROR] There is no image in imagePaths")
+            return
+
+        # Test key method
+        if method not in DataHandler.reconstruction_images.keys():
+            print("[ERROR] The method {} is not valid.".format(method))
+            exit(1)
+
+        nameImage = Path(DataHandler.imagePaths[indexImage]).stem
+
+        if DataHandler.CFA_images.get(nameImage) is None:
+            print("[ERROR] There is no CFA image for the image {}.".format(nameImage))
+            exit(1)
+
+        img = DataHandler.load_image(self, indexImage)
+
+        img_CFA = DataHandler.CFA_images[nameImage].direct(img)
+
+        cfa = options.get("cfa")
+
+        if cfa is None:
+            print("[ERROR] You must specify the cfa.")
+            exit(1)
+
+        if method == "interpolation":
+
+            DataHandler.reconstruction_images[method].setdefault(nameImage, run_reconstruction(img_CFA, cfa))
+
+        if method == "menon":
+
+            DataHandler.reconstruction_images[method].setdefault(nameImage, run_reconstruction_brice_convers(img_CFA, cfa))
+
+    def compute_reconstruction_images(self, method, options = None):
+
+        for i in range(len(DataHandler.imagePaths)):
+            DataHandler.compute_reconstruction_image(self, method, i, options)
+
+        print("[INFO] There are {} images which have been reconstructed.".format(len(DataHandler.reconstruction_images[method])))
+
+    def plot_reconstructed_image(self, indexImage, method, cfa = {"cfa": "bayer"}, zoomSize = "small", rows=1, cols=4, figsize=(15, 5)):
+        # Test key method
+        if method not in DataHandler.reconstruction_images.keys():
+            print("[ERROR] The method {} is not valid.".format(method))
+            exit(1)
+
+        res = DataHandler.get_reconstructed_image(self, indexImage, method)
+
+        fig, axs = plt.subplots(rows, cols, figsize=figsize)
+        fig.suptitle("Reconstructed Image with method: {} and pattern type: {}".format(method, cfa["cfa"]), fontsize=16)
+
+        axs[0].imshow(DataHandler.load_image(self, indexImage))
+        axs[0].set_title('Original Image')
+        axs[1].imshow(res)
+        axs[1].set_title('Reconstructed Image')
+
+        if zoomSize == "small":
+            axs[2].imshow(DataHandler.load_image(self, indexImage)[800:864, 450:514])
+            axs[2].set_title('Zoomed Input Image')
+            axs[3].imshow(res[800:864, 450:514])
+            axs[3].set_title('Zoomed Reconstructed Image')
+        else:
+            axs[2].imshow(DataHandler.load_image(self, indexImage)[2000:2064, 2000:2064])
+            axs[2].set_title('Zoomed Input Image')
+            axs[3].imshow(res[2000:2064, 2000:2064])
+            axs[3].set_title('Zoomed Reconstructed Image')
+
+        outputPath = os.path.join(DataHandler.WD_PATH, "output")
+        folderExists(outputPath)
+
+        imageName = Path(DataHandler.imagePaths[indexImage]).stem
+        plotCompImagePath = os.path.join(outputPath, "Compararison_" + imageName + "_" + method + "_" + cfa["cfa"] + ".png")
+
+        #save_image( plotReconstructedImagePath, res)
+        fig.savefig(plotCompImagePath)
+
+
+    def get_reconstructed_image(self, indexImage, method):
+        # Test key method
+        if method not in DataHandler.reconstruction_images.keys():
+            print("[ERROR] The method {} is not valid.".format(method))
+            exit(1)
+
+        DataHandler.indexImageExists(self, indexImage)
+
+        nameImage = Path(DataHandler.imagePaths[indexImage]).stem
+
+        if DataHandler.reconstruction_images[method].get(nameImage) is None:
+            print("[ERROR] There is no reconstruction image for the image {} and for the method {}.".format(nameImage, method))
+            exit(1)
+
+        return DataHandler.reconstruction_images[method][nameImage]
+
+    def shape_of_recontructed_image(self, indexImage, method):
+        if DataHandler.reconstruction_images[method].get(Path(DataHandler.imagePaths[indexImage]).stem) is None:
+            print("[ERROR] There is no reconstruction image for the image {} and for the method {}.".format(Path(DataHandler.imagePaths[indexImage]).stem, method))
+            exit(1)
+
+        return DataHandler.reconstruction_images[method][Path(DataHandler.imagePaths[indexImage]).stem].shape
+
+    def compute_CFA_image(self, CFA_NAME, nameImage, indexImage):
+        DataHandler.CFA_images.setdefault(nameImage, CFA(CFA_NAME, DataHandler.shape_of_image(self, indexImage)))
+
+    def shape_of_image(self, indexImage):
+        img = mpimg.imread(DataHandler.imagePaths[indexImage])
+        return img.shape
+
+    def load_image(self, indexImage):
+        img = imread(DataHandler.imagePaths[indexImage])
+        img = normalise_image(img)
+        return img
+
+    def indexImageExists(self, indexImage):
+        if indexImage >= len(DataHandler.imagePaths):
+            print("[ERROR] The index {} is not valid.".format(indexImage))

src/methods/brice_convers/main_template.py

+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})
+
+    DataHandler.plot_reconstructed_image(0, METHOD, {"cfa": CFA_NAME}, zoomSize="large")
+
+    DataEvaluation.print_metrics(0, METHOD)
+
+    endTime = time.time()
+
+    print("[INFO] Elapsed time: " + str(endTime - startTime) + "s")
+
+    print("[INFO] End")
+
+if __name__ == "__main__":
+    main(DataHandler)

src/methods/brice_convers/menon.py

+"""
+DDFAPD - Menon (2007) Bayer CFA Demosaicing
+===========================================
+
+*Bayer* CFA (Colour Filter Array) DDFAPD - *Menon (2007)* demosaicing.
+
+References
+----------
+-   :cite:`Menon2007c` : Menon, D., Andriani, S., & Calvagno, G. (2007).
+    Demosaicing With Directional Filtering and a posteriori Decision. IEEE
+    Transactions on Image Processing, 16(1), 132-141.
+    doi:10.1109/TIP.2006.884928
+"""
+
+import numpy as np
+from colour.hints import ArrayLike, Literal, NDArrayFloat
+from colour.utilities import as_float_array, ones, tsplit, tstack
+from scipy.ndimage.filters import convolve, convolve1d
+from src.forward_model import CFA
+
+def tensor_mask_to_RGB_mask(mask: ArrayLike, pixelPattern: str = "RGB"):
+        # We extract image chanels from mask
+    for i, letter in enumerate(pixelPattern):
+        if letter == "R":
+            R_m = mask[:, :, i]
+        elif letter == "G":
+            G_m = mask[:, :, i]
+        elif letter == "B":
+            B_m = mask[:, :, i]
+
+    return R_m, G_m, B_m
+
+def _cnv_h(x: ArrayLike, y: ArrayLike) -> NDArrayFloat:
+    """Perform horizontal convolution."""
+    # we go through the rows because axis = -1
+    return convolve1d(x, y, mode="mirror")
+
+
+def _cnv_v(x: ArrayLike, y: ArrayLike) -> NDArrayFloat:
+    """Perform vertical convolution."""
+
+    return convolve1d(x, y, mode="mirror", axis=0)
+
+
+def demosaicing_CFA_Bayer_Menon2007(
+    rawImage: ArrayLike,
+    mask: ArrayLike,
+    pixelPattern: str = "RGB",
+    refining_step: bool = True,
+):
+    """
+    Return the demosaiced *RGB* colourspace array from given *Bayer* CFA using
+    DDFAPD - *Menon (2007)* demosaicing algorithm.
+
+    Parameters
+    ----------
+    CFA
+        *Bayer* CFA.
+    pattern
+        Arrangement of the colour filters on the pixel array.
+    refining_step
+        Perform refining step.
+
+    Returns
+    -------
+    :class:`numpy.ndarray`
+        *RGB* colourspace array.
+
+    Notes
+    -----
+    -   The definition output is not clipped in range [0, 1] : this allows for
+        direct HDRI image generation on *Bayer* CFA data and post
+        demosaicing of the high dynamic range data as showcased in this
+        `Jupyter Notebook <https://github.com/colour-science/colour-hdri/\
+blob/develop/colour_hdri/examples/\
+examples_merge_from_raw_files_with_post_demosaicing.ipynb>`__.
+
+    References
+    ----------
+    :cite:`Menon2007c`
+    """
+
+    # We extract image chanels from mask
+    R_m, G_m, B_m = tensor_mask_to_RGB_mask(mask, pixelPattern)
+
+    # We extract known pixel intensities: when we have a zero in the mask, we have an unknown pixel intensity for the color
+    R = rawImage * R_m
+    G = rawImage * G_m
+    B = rawImage * B_m
+
+    # We define the horizontal and vertical filters
+    h_0 = as_float_array([0.0, 0.5, 0.0, 0.5, 0.0])
+    h_1 = as_float_array([-0.25, 0.0, 0.5, 0.0, -0.25])
+
+    # Green components interpolation along both horizontal and veritcal directions:
+    # For each unkown green pixel, we compute the gradient along both horizontal and vertical directions
+    G_H = np.where(G_m == 0, _cnv_h(rawImage, h_0) + _cnv_h(rawImage, h_1), G)
+    G_V = np.where(G_m == 0, _cnv_v(rawImage, h_0) + _cnv_v(rawImage, h_1), G)
+
+    # We calculate the chrominance differences along both horizontal and vertical directions
+    # For each unknown red and blue pixel, we compute the difference between the pixel intensity and the horizontal green component
+    C_H = np.where(R_m == 1, R - G_H, 0)
+    C_H = np.where(B_m == 1, B - G_H, C_H)
+
+    # Sale method with vertical green component
+    C_V = np.where(R_m == 1, R - G_V, 0)
+    C_V = np.where(B_m == 1, B - G_V, C_V)
+
+    # We compute the directional gradients along both horizontal and vertical directions
+    # First we pad our arrayes with zeros to avoid boundary effects. Acxtually, we pad with the last value of the array
+    # We add two columns to the right of the horizontal array and two rows at the bottom of the  vertical array, with the reflect mode.
+    # Then we remove the first two columns of the horizontal array and the first two rows of the vertical array.
+    paded_D_H = np.pad(C_H, ((0, 0), (0, 2)), mode="reflect")[:, 2:]
+    paded_D_V = np.pad(C_V, ((0, 2), (0, 0)), mode="reflect")[2:, :]
+
+    # We compute the difference between the original array and the padded array.
+    # With the paded array, we have a difference between each pixel and the right neigborhood. We do not have issue with boundaries.
+    # It gives a measure of pixel intensity variation along the horizontal and vertical directions.
+    D_H = np.abs(C_H - paded_D_H)
+    D_V = np.abs(C_V - paded_D_V)
+
+    del h_0, h_1, C_V, C_H, paded_D_V, paded_D_H
+
+    # We define a sufficiently large neighborhood with a size of (5, 5).
+    k = as_float_array(
+        [
+            [0.0, 0.0, 1.0, 0.0, 1.0],
+            [0.0, 0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, 3.0, 0.0, 3.0],
+            [0.0, 0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, 1.0, 0.0, 1.0],
+        ]
+    )
+
+    # We convolve the difference component with the neighborhood. This method is used to highlight directional variations in the image, in two direction.
+    d_H = convolve(D_H, k, mode="constant")
+    d_V = convolve(D_V, np.transpose(k), mode="constant")
+
+    del D_H, D_V
+
+    # We estimate the green channel with our classifier
+    mask = d_V >= d_H
+    G = np.where(mask, G_H, G_V)
+    # We estimate the mask which represents the best directional reconstruction
+    M = np.where(mask, 1, 0)
+
+    del d_H, d_V, G_H, G_V
+
+    ## The, we estimate the red and blue channels
+
+    # We arrays with ones at the line where there is at least one red (blue) pixel in the red (blue) mask
+    R_r = np.transpose(np.any(R_m == 1, axis=1)[None]) * ones(R.shape)
+    B_r = np.transpose(np.any(B_m == 1, axis=1)[None]) * ones(B.shape)
+
+    # We define a new filter
+    k_b = as_float_array([0.5, 0, 0.5])
+
+    # We fill R array with the condition: if we are in a line where there is at least one red pixel in the red mask and we are on a green pixel in the green mask, we apply the filter horizontaly to the red channel.
+    # If not it means we are on a red pixel (only two possiblity) in the red mask, so we do not apply the filter because we know the red pixel
+    R = np.where(
+        np.logical_and(G_m == 1, R_r == 1),
+        G + _cnv_h(R, k_b) - _cnv_h(G, k_b),
+        R,
+    )
+
+    # Same but we test only the line where there is at least one blue pixel in the blue mask.
+    # When the condition is true, we apply the filter vertically because this time red pixel are aline vertically.
+    R = np.where(
+        np.logical_and(G_m == 1, B_r == 1) == 1,
+        G + _cnv_v(R, k_b) - _cnv_v(G, k_b),
+        R,
+    )
+
+    # It is the same logic for the blue image
+    B = np.where(
+        np.logical_and(G_m == 1, B_r == 1),
+        G + _cnv_h(B, k_b) - _cnv_h(G, k_b),
+        B,
+    )
+
+    B = np.where(
+        np.logical_and(G_m == 1, R_r == 1) == 1,
+        G + _cnv_v(B, k_b) - _cnv_v(G, k_b),
+        B,
+    )
+
+    # To finish R image we need to interpolate blue pixel. We use M to know wich direction is the best and then we interpolate the blue pixel with the filter.
+    R_b = np.where(
+            M == 1,
+            B + _cnv_h(R, k_b) - _cnv_h(B, k_b),
+            B + _cnv_v(R, k_b) - _cnv_v(B, k_b),
+        )
+
+    # Then we put the condition: if we are on a line where there is at least one blue pixel and we are on a blue pixel we take the previous interpolated value.
+    # If not we know the red pixel value and we keep it.
+    R = np.where(
+        np.logical_and(B_r == 1, B_m == 1),
+        R_b,
+        R,
+    )
+
+    # Same idea for the blue image.
+    B = np.where(
+        np.logical_and(R_r == 1, R_m == 1),
+        np.where(
+            M == 1,
+            R + _cnv_h(B, k_b) - _cnv_h(R, k_b),
+            R + _cnv_v(B, k_b) - _cnv_v(R, k_b),
+        ),
+        B,
+    )
+
+    # We stack the channels in the last dimension to get the final image
+    RGB = tstack([R, G, B])
+
+    del R, G, B, k_b, R_r, B_r
+
+    # We optionally perform the refining step
+    if refining_step:
+        RGB = refining_step_Menon2007(RGB, tstack([R_m, G_m, B_m]), M)
+
+    del M, R_m, G_m, B_m
+
+    return RGB
+
+
+def refining_step_Menon2007(
+    RGB: ArrayLike, RGB_m: ArrayLike, M: ArrayLike
+) -> NDArrayFloat:
+    """
+    Perform the refining step on given *RGB* colourspace array.
+
+    Parameters
+    ----------
+    RGB
+        *RGB* colourspace array.
+    RGB_m
+        *Bayer* CFA red, green and blue masks.
+    M
+        Estimation for the best directional reconstruction.
+
+    Returns
+    -------
+    :class:`numpy.ndarray`
+        Refined *RGB* colourspace array.
+    """
+    # Unpacking the RGB and RGB_m arrays.
+    R, G, B = tsplit(RGB)
+    R_m, G_m, B_m = tsplit(RGB_m)
+    M = as_float_array(M)
+
+    del RGB, RGB_m
+
+    # Updating of the green component.
+    R_G = R - G
+    B_G = B - G
+
+    # Definition of the low-pass filter.
+    FIR = ones(3) / 3
+
+    # When we are on a blue pixel, we convolve the pixel with the filter in function of the best direction.
+    B_G_m = np.where(
+        B_m == 1,
+        np.where(M == 1, _cnv_h(B_G, FIR), _cnv_v(B_G, FIR)),
+        0,
+    )
+
+    # Same for the red pixel.
+    R_G_m = np.where(
+        R_m == 1,
+        np.where(M == 1, _cnv_h(R_G, FIR), _cnv_v(R_G, FIR)),
+        0,
+    )
+
+    del B_G, R_G
+
+    # We update the green component for known red and blue pixels with the difference between the red or blue pixel intensity and the filtered pixel intensity.
+    G = np.where(R_m == 1, R - R_G_m, G)
+    G = np.where(B_m == 1, B - B_G_m, G)
+
+    # Updating of the red and blue components in the green locations.
+
+    # R_r is an array with ones at the line where there is at least one red pixel in the red mask.
+    R_r = np.transpose(np.any(R_m == 1, axis=1)[None]) * ones(R.shape)
+    # R_c is an array with ones at the column where there is at least one red pixel in the red mask.
+    R_c = np.any(R_m == 1, axis=0)[None] * ones(R.shape)
+    # B_r is an array with ones at the line where there is at least one blue pixel in the blue mask.
+    B_r = np.transpose(np.any(B_m == 1, axis=1)[None]) * ones(B.shape)
+    # B_c is an array with ones at the column where there is at least one blue pixel in the blue mask.
+    B_c = np.any(B_m == 1, axis=0)[None] * ones(B.shape)
+
+    R_G = R - G
+    B_G = B - G
+
+    k_b = as_float_array([0.5, 0.0, 0.5])
+
+    R_G_m = np.where(
+        np.logical_and(G_m == 1, B_r == 1),
+        _cnv_v(R_G, k_b),
+        R_G_m,
+    )
+    R = np.where(np.logical_and(G_m == 1, B_r == 1), G + R_G_m, R)
+    R_G_m = np.where(
+        np.logical_and(G_m == 1, B_c == 1),
+        _cnv_h(R_G, k_b),
+        R_G_m,
+    )
+    R = np.where(np.logical_and(G_m == 1, B_c == 1), G + R_G_m, R)
+
+    del B_r, R_G_m, B_c, R_G
+
+    B_G_m = np.where(
+        np.logical_and(G_m == 1, R_r == 1),
+        _cnv_v(B_G, k_b),
+        B_G_m,
+    )
+    B = np.where(np.logical_and(G_m == 1, R_r == 1), G + B_G_m, B)
+    B_G_m = np.where(
+        np.logical_and(G_m == 1, R_c == 1),
+        _cnv_h(B_G, k_b),
+        B_G_m,
+    )
+    B = np.where(np.logical_and(G_m == 1, R_c == 1), G + B_G_m, B)
+
+    del B_G_m, R_r, R_c, G_m, B_G
+
+    # Updating of the red (blue) component in the blue (red) locations.
+    R_B = R - B
+    R_B_m = np.where(
+        B_m == 1,
+        np.where(M == 1, _cnv_h(R_B, FIR), _cnv_v(R_B, FIR)),
+        0,
+    )
+    R = np.where(B_m == 1, B + R_B_m, R)
+
+    R_B_m = np.where(
+        R_m == 1,
+        np.where(M == 1, _cnv_h(R_B, FIR), _cnv_v(R_B, FIR)),
+        0,
+    )
+    B = np.where(R_m == 1, R - R_B_m, B)
+
+    del R_B, R_B_m, R_m
+
+    return tstack([R, G, B])