src/methods/dzik_erwan/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.dzik_erwan.functions import malvar2004, quad_bayer_to_bayer
+
+
+
+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)
+
+    if cfa == 'bayer':
+        res = malvar2004(y,op)
+    else:
+        bayer_equivalent = quad_bayer_to_bayer(y)
+        res = malvar2004(bayer_equivalent, CFA("bayer", input_shape))
+    
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/fortin_come/DemosaickingInterChannel.py

+######################################################################################
+# Demosaicking algorithm based on "A Demosaicking Algorithme with 
+# Adaptive Inter-Channel Correlation" paper from Joan Duran and Antoni Buades
+######################################################################################
+
+import numpy as np
+
+
+
+
+
+
+def cvt_YUV_space(R: np.ndarray, G: np.ndarray, B: np.ndarray):
+    """Convert RGB into YUV
+
+    Args:
+        R (np.ndarray): Red channel of image
+        G (np.ndarray): Green channel of image
+        B (np.ndarray): Blue channel of image
+
+    Returns:
+        np.ndarray: Image in the YUV space
+    """
+    Y = 0.299*R + 0.587*G + 0.114*B
+    U = R-Y
+    V = B-Y
+    return Y, U, V
+
+
+
+
+
+
+def compute_gradient(L: int, U: np.ndarray, V: np.ndarray, i: int, j: int, dir: str):
+    """Compute the gradient
+
+    Args:
+        L (int): Size of the local neighborhood
+        U (np.ndarray): U channel of the image
+        V (np.ndarray): V channel of the image
+        i (int): position in line
+        j (int): position in column
+        dir (string) : direction of interpolation
+
+    Returns:
+        float: Gradient
+    """
+    height, width = U.shape
+    gradient_U = 0
+    gradient_V = 0
+    
+    if dir == "north":
+        if j < L+1: max = j+1
+        else : max = L+1
+        for l in range(1, max):
+            gradient_U += (U[j-l, i] - U[j, i])**2
+            gradient_V += (V[j-l, i] - V[j, i])**2
+    elif dir == "south":
+        if height-j < L+1: max = height-j
+        else: max = L+1
+        for l in range(1, max):
+            gradient_U += (U[j+l, i] - U[j, i])**2
+            gradient_V += (V[j+l, i] - V[j, i])**2
+    elif dir == "east":
+        if width-i < L+1: max = width-i
+        else: max = L+1
+        for l in range(1, max):
+            gradient_U += (U[j, i+l] - U[j, i])**2
+            gradient_V += (V[j, i+l] - V[j, i])**2
+    elif dir == "west":
+        if i < L+1: max = i+1
+        else: max = L+1
+        for l in range(1, max):
+            gradient_U += (U[j, i-l] - U[j, i])**2
+            gradient_V += (V[j, i-l] - V[j, i])**2
+    
+    if max == L+1: max = L
+    
+    gradient = (np.sqrt(gradient_U) + np.sqrt(gradient_V)) / max
+    
+    return gradient
+
+
+
+
+
+
+
+
+
+
+
+def local_int(R0: np.ndarray, G0: np.ndarray, B0: np.ndarray, beta: float, L: int , eps: float):
+    """Local direction interpolation algorithm
+
+    Args:
+        R0 (np.ndarray): Red mosaiked image
+        G0 (np.ndarray): Green mosaiked image
+        B0 (np.ndarray): Blue mosaiked image
+        beta (float): Channel correlation parameter [0; 1]
+        L (int): Size of the local neighborhood
+        eps (float): tresholding parameter > 0
+
+    Returns:
+        np.ndarray: The interpolated image (R, G, B)
+    """
+    
+    assert R0.shape == G0.shape and R0.shape == B0.shape
+    assert eps > 0
+    
+    lin, col = G0.shape
+    
+    # Initalisation
+    
+    Rn, Rs, Re, Rw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
+    Gn, Gs, Ge, Gw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
+    Bn, Bs, Be, Bw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
+    
+    RGn, RGs, RGe, RGw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
+    BGn, BGs, BGe, BGw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
+    
+    Gradn, Grads, Grade, Gradw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
+    
+    low_wn, low_ws, low_we, low_ww = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
+    W = np.zeros((lin, col))
+    
+    R, G, B = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
+    
+    
+    # Computation
+    
+    for i in range(col):
+        for j in range(lin):
+            
+            # Directionnal interpolation of green channel
+            if G0[j, i] != 0:
+                Gn[j, i], Gs[j, i], Ge[j, i], Gw[j, i] = G0[j, i], G0[j, i], G0[j, i], G0[j, i]
+            else:
+                if j-2 >= 0: # North
+                    # If we are on a Red square: B0[j, i] = B0[j-2, i] = 0
+                    # The same if we are on a Blue square
+                    Gn[j, i] = G0[j-1, i] + 0.5*beta*(R0[j, i] - R0[j-2, i] + B0[j, i] - B0[j-2, i])
+                else: # R0[j-2, i] = 0 or B0[j-2, i] = 0
+                    if j-1 >= 0:
+                        Gn[j, i] = G0[j-1, i] + beta*(R0[j, i] + B0[j, i])
+                    else: # G0[j-1, i] = 0
+                        Gn[j, i] = beta*(R0[j, i] + B0[j, i])
+                
+                if j+2 < lin: # South
+                    Gs[j, i] = G0[j+1, i] + 0.5*beta*(R0[j, i] - R0[j+2, i] + B0[j, i] - B0[j+2, i])
+                else:
+                    if j+1 < lin:
+                        Gs[j, i] = G0[j+1, i] + beta*(R0[j, i] + B0[j, i])
+                    else:
+                        Gs[j, i] = beta*(R0[j, i] + B0[j, i])
+                
+                if i+2 < col: # East
+                    Ge[j, i] = G0[j, i+1] + 0.5*beta*(R0[j, i] - R0[j, i+2] + B0[j, i] - B0[j, i+2])
+                else:
+                    if i+1 < col:
+                        Ge[j, i] = G0[j, i+1] + beta*(R0[j, i] + B0[j, i])
+                    else:
+                        Ge[j, i] = beta*(R0[j, i] + B0[j, i])
+                
+                if i-2 >= 0: # West
+                    Gw[j, i] = G0[j, i-1] + 0.5*beta*(R0[j, i] - R0[j, i-2] + B0[j, i] - B0[j, i-2])
+                else:
+                    if i-1 >= 0:
+                        Gw[j, i] = G0[j, i-1] + beta*(R0[j, i] + B0[j, i])
+                    else:
+                        Gw[j, i] = beta*(R0[j, i] + B0[j, i])
+            
+            
+            
+            # Bilinear interpolation of red and blue channels
+            if R0[j, i] != 0:
+                Rn[j, i], Rs[j, i], Re[j, i], Rw[j, i] = R0[j, i], R0[j, i], R0[j, i], R0[j, i]
+                
+                RGn[j, i] = R0[j, i] - beta*Gn[j, i]
+                RGs[j, i] = R0[j, i] - beta*Gs[j, i]
+                RGe[j, i] = R0[j, i] - beta*Ge[j, i]
+                RGw[j, i] = R0[j, i] - beta*Gw[j, i]
+            if B0[j, i] != 0:
+                Bn[j, i], Bs[j, i], Be[j, i], Bw[j, i] = B0[j, i], B0[j, i], B0[j, i], B0[j, i]
+                
+                BGn[j, i] = B0[j, i] - beta*Gn[j, i]
+                BGs[j, i] = B0[j, i] - beta*Gs[j, i]
+                BGe[j, i] = B0[j, i] - beta*Ge[j, i]
+                BGw[j, i] = B0[j, i] - beta*Gw[j, i]
+            
+            if i+1 < col:
+                if G0[j, i] != 0 and R0[j, i+1] != 0: # Red square next to us
+                    if i-1 >=0 :
+                        Rn[j, i] = 0.5*(RGn[j, i-1] + RGn[j, i+1]) + beta*Gn[j, i]
+                        Rs[j, i] = 0.5*(RGs[j, i-1] + RGs[j, i+1]) + beta*Gs[j, i]
+                        Re[j, i] = 0.5*(RGe[j, i-1] + RGe[j, i+1]) + beta*Ge[j, i]
+                        Rw[j, i] = 0.5*(RGw[j, i-1] + RGw[j, i+1]) + beta*Gw[j, i]
+
+                    
+                    # The case nest to us is red, so it can't be blue
+                    if j-1 >=0:
+                        if j+1 < lin:
+                            Bn[j, i] = 0.5*(BGn[j-1, i] + BGn[j+1, i]) + beta*Gn[j, i]
+                            Bs[j, i] = 0.5*(BGs[j-1, i] + BGs[j+1, i]) + beta*Gs[j, i]
+                            Be[j, i] = 0.5*(BGe[j-1, i] + BGe[j+1, i]) + beta*Ge[j, i]
+                            Bw[j, i] = 0.5*(BGw[j-1, i] + BGw[j+1, i]) + beta*Gw[j, i]
+                        else:
+                            Bn[j, i] = BGn[j-1, i] + beta*Gn[j, i]
+                            Bs[j, i] = BGs[j-1, i] + beta*Gs[j, i]
+                            Be[j, i] = BGe[j-1, i] + beta*Ge[j, i]
+                            Bw[j, i] = BGw[j-1, i] + beta*Gw[j, i]
+                    else:
+                        Bn[j, i] = BGn[j+1, i] + beta*Gn[j, i]
+                        Bs[j, i] = BGs[j+1, i] + beta*Gs[j, i]
+                        Be[j, i] = BGe[j+1, i] + beta*Ge[j, i]
+                        Bw[j, i] = BGw[j+1, i] + beta*Gw[j, i]
+                
+                if G0[j, i] != 0 and B0[j, i+1] != 0: # Blue square next to us
+                    if i-1 >=0 :
+                        Bn[j, i] = 0.5*(BGn[j, i-1] + BGn[j, i+1]) + beta*Gn[j, i]
+                        Bs[j, i] = 0.5*(BGs[j, i-1] + BGs[j, i+1]) + beta*Gs[j, i]
+                        Be[j, i] = 0.5*(BGe[j, i-1] + BGe[j, i+1]) + beta*Ge[j, i]
+                        Bw[j, i] = 0.5*(BGw[j, i-1] + BGw[j, i+1]) + beta*Gw[j, i]
+                        
+                        # The case nest to us is blue, so it can't be red
+                        if j-1 >=0:
+                            if j+1 < lin:
+                                Rn[j, i] = 0.5*(RGn[j-1, i] + RGn[j+1, i]) + beta*Gn[j, i]
+                                Rs[j, i] = 0.5*(RGs[j-1, i] + RGs[j+1, i]) + beta*Gs[j, i]
+                                Re[j, i] = 0.5*(RGe[j-1, i] + RGe[j+1, i]) + beta*Ge[j, i]
+                                Rw[j, i] = 0.5*(RGw[j-1, i] + RGw[j+1, i]) + beta*Gw[j, i]
+                            else:
+                                Rn[j, i] = RGn[j-1, i] + beta*Gn[j, i]
+                                Rs[j, i] = RGs[j-1, i] + beta*Gs[j, i]
+                                Re[j, i] = RGe[j-1, i] + beta*Ge[j, i]
+                                Rw[j, i] = RGw[j-1, i] + beta*Gw[j, i]
+            
+            if R0[j, i] == 0 and G0[j, i] == 0:
+                if i-1 >= 0:
+                    if i+1 < col:
+                        if j-1 >= 0:
+                            if j+1 < lin:
+                                Rn[j, i] = 0.25*(RGn[j-1, i-1] + RGn[j-1, i+1] + RGn[j+1, i-1] + RGn[j+1, i+1]) + beta*Gn[j, i]
+                                Rs[j, i] = 0.25*(RGs[j-1, i-1] + RGs[j-1, i+1] + RGs[j+1, i-1] + RGs[j+1, i+1]) + beta*Gs[j, i]
+                                Re[j, i] = 0.25*(RGe[j-1, i-1] + RGe[j-1, i+1] + RGe[j+1, i-1] + RGe[j+1, i+1]) + beta*Ge[j, i]
+                                Rw[j, i] = 0.25*(RGw[j-1, i-1] + RGw[j-1, i+1] + RGw[j+1, i-1] + RGw[j+1, i+1]) + beta*Gw[j, i]
+                            else:
+                                Rn[j, i] = 0.5*(RGn[j-1, i-1] + RGn[j-1, i+1]) + beta*Gn[j, i]
+                                Rs[j, i] = 0.5*(RGs[j-1, i-1] + RGs[j-1, i+1]) + beta*Gs[j, i]
+                                Re[j, i] = 0.5*(RGe[j-1, i-1] + RGe[j-1, i+1]) + beta*Ge[j, i]
+                                Rw[j, i] = 0.5*(RGw[j-1, i-1] + RGw[j-1, i+1]) + beta*Gw[j, i]
+                        else:
+                            Rn[j, i] = 0.5*(RGn[j+1, i-1] + RGn[j+1, i+1]) + beta*Gn[j, i]
+                            Rs[j, i] = 0.5*(RGs[j+1, i-1] + RGs[j+1, i+1]) + beta*Gs[j, i]
+                            Re[j, i] = 0.5*(RGe[j+1, i-1] + RGe[j+1, i+1]) + beta*Ge[j, i]
+                            Rw[j, i] = 0.5*(RGw[j+1, i-1] + RGw[j+1, i+1]) + beta*Gw[j, i]
+                    else:
+                        if j-1 >= 0:
+                            if j+1 < lin:
+                                Rn[j, i] = 0.5*(RGn[j-1, i-1] + RGn[j+1, i-1]) + beta*Gn[j, i]
+                                Rs[j, i] = 0.5*(RGs[j-1, i-1] + RGs[j+1, i-1]) + beta*Gs[j, i]
+                                Re[j, i] = 0.5*(RGe[j-1, i-1] + RGe[j+1, i-1]) + beta*Ge[j, i]
+                                Rw[j, i] = 0.5*(RGw[j-1, i-1] + RGw[j+1, i-1]) + beta*Gw[j, i]
+                            else:
+                                Rn[j, i] = RGn[j-1, i-1] + beta*Gn[j, i]
+                                Rs[j, i] = RGs[j-1, i-1] + beta*Gs[j, i]
+                                Re[j, i] = RGe[j-1, i-1] + beta*Ge[j, i]
+                                Rw[j, i] = RGw[j-1, i-1] + beta*Gw[j, i]
+                        else:
+                            Rn[j, i] = beta*Gn[j, i]
+                            Rs[j, i] = beta*Gs[j, i]
+                            Re[j, i] = beta*Ge[j, i]
+                            Rw[j, i] = beta*Gw[j, i]
+                else:
+                    if j-1 >= 0:
+                        if j+1 < lin:
+                            Rn[j, i] = 0.5*(RGn[j-1, i+1] + RGn[j+1, i+1]) + beta*Gn[j, i]
+                            Rs[j, i] = 0.5*(RGs[j-1, i+1] + RGs[j+1, i+1]) + beta*Gs[j, i]
+                            Re[j, i] = 0.5*(RGe[j-1, i+1] + RGe[j+1, i+1]) + beta*Ge[j, i]
+                            Rw[j, i] = 0.5*(RGw[j-1, i+1] + RGw[j+1, i+1]) + beta*Gw[j, i]
+                        else:
+                            Rn[j, i] = RGn[j-1, i+1] + beta*Gn[j, i]
+                            Rs[j, i] = RGs[j-1, i+1] + beta*Gs[j, i]
+                            Re[j, i] = RGe[j-1, i+1] + beta*Ge[j, i]
+                            Rw[j, i] = RGw[j-1, i+1] + beta*Gw[j, i]
+                    else:
+                        Rn[j, i] = beta*Gn[j, i]
+                        Rs[j, i] = beta*Gs[j, i]
+                        Re[j, i] = beta*Ge[j, i]
+                        Rw[j, i] = beta*Gw[j, i]
+                
+            if B0[j, i] == 0 and G0[j, i] == 0:
+                if i-1 >= 0:
+                    if i+1 < col:
+                        if j-1 >= 0:
+                            if j+1 < lin:
+                                Bn[j, i] = 0.25*(BGn[j-1, i-1] + BGn[j-1, i+1] + BGn[j+1, i-1] + BGn[j+1, i+1]) + beta*Gn[j, i]
+                                Bs[j, i] = 0.25*(BGs[j-1, i-1] + BGs[j-1, i+1] + BGs[j+1, i-1] + BGs[j+1, i+1]) + beta*Gs[j, i]
+                                Be[j, i] = 0.25*(BGe[j-1, i-1] + BGe[j-1, i+1] + BGe[j+1, i-1] + BGe[j+1, i+1]) + beta*Ge[j, i]
+                                Bw[j, i] = 0.25*(BGw[j-1, i-1] + BGw[j-1, i+1] + BGw[j+1, i-1] + BGw[j+1, i+1]) + beta*Gw[j, i]
+                            else:
+                                Bn[j, i] = 0.5*(BGn[j-1, i-1] + BGn[j-1, i+1]) + beta*Gn[j, i]
+                                Bs[j, i] = 0.5*(BGs[j-1, i-1] + BGs[j-1, i+1]) + beta*Gs[j, i]
+                                Be[j, i] = 0.5*(BGe[j-1, i-1] + BGe[j-1, i+1]) + beta*Ge[j, i]
+                                Bw[j, i] = 0.5*(BGw[j-1, i-1] + BGw[j-1, i+1]) + beta*Gw[j, i]
+                        else:
+                            Bn[j, i] = 0.5*(BGn[j+1, i-1] + BGn[j+1, i+1]) + beta*Gn[j, i]
+                            Bs[j, i] = 0.5*(BGs[j+1, i-1] + BGs[j+1, i+1]) + beta*Gs[j, i]
+                            Be[j, i] = 0.5*(BGe[j+1, i-1] + BGe[j+1, i+1]) + beta*Ge[j, i]
+                            Bw[j, i] = 0.5*(BGw[j+1, i-1] + BGw[j+1, i+1]) + beta*Gw[j, i]
+                    else:
+                        if j-1 >= 0:
+                            if j+1 < lin:
+                                Bn[j, i] = 0.5*(BGn[j-1, i-1] + BGn[j+1, i-1]) + beta*Gn[j, i]
+                                Bs[j, i] = 0.5*(BGs[j-1, i-1] + BGs[j+1, i-1]) + beta*Gs[j, i]
+                                Be[j, i] = 0.5*(BGe[j-1, i-1] + BGe[j+1, i-1]) + beta*Ge[j, i]
+                                Bw[j, i] = 0.5*(BGw[j-1, i-1] + BGw[j+1, i-1]) + beta*Gw[j, i]
+                            else:
+                                Bn[j, i] = BGn[j-1, i-1] + beta*Gn[j, i]
+                                Bs[j, i] = BGs[j-1, i-1] + beta*Gs[j, i]
+                                Be[j, i] = BGe[j-1, i-1] + beta*Ge[j, i]
+                                Bw[j, i] = BGw[j-1, i-1] + beta*Gw[j, i]
+                        else:
+                            Bn[j, i] = beta*Gn[j, i]
+                            Bs[j, i] = beta*Gs[j, i]
+                            Be[j, i] = beta*Ge[j, i]
+                            Bw[j, i] = beta*Gw[j, i]
+                else:
+                    if j-1 >= 0:
+                        if j+1 < lin:
+                            Bn[j, i] = 0.5*(BGn[j-1, i+1] + BGn[j+1, i+1]) + beta*Gn[j, i]
+                            Bs[j, i] = 0.5*(BGs[j-1, i+1] + BGs[j+1, i+1]) + beta*Gs[j, i]
+                            Be[j, i] = 0.5*(BGe[j-1, i+1] + BGe[j+1, i+1]) + beta*Ge[j, i]
+                            Bw[j, i] = 0.5*(BGw[j-1, i+1] + BGw[j+1, i+1]) + beta*Gw[j, i]
+                        else:
+                            Bn[j, i] = BGn[j-1, i+1] + beta*Gn[j, i]
+                            Bs[j, i] = BGs[j-1, i+1] + beta*Gs[j, i]
+                            Be[j, i] = BGe[j-1, i+1] + beta*Ge[j, i]
+                            Bw[j, i] = BGw[j-1, i+1] + beta*Gw[j, i]
+                    else:
+                        Bn[j, i] = beta*Gn[j, i]
+                        Bs[j, i] = beta*Gs[j, i]
+                        Be[j, i] = beta*Ge[j, i]
+                        Bw[j, i] = beta*Gw[j, i]
+
+
+
+            # Pixel-level fusion of full color interpolated images
+            
+            # Computation of the YUV space
+            
+            _, Un, Vn = cvt_YUV_space(Rn, Gn, Bn)
+            _, Us, Vs = cvt_YUV_space(Rs, Gs, Bs)
+            _, Ue, Ve = cvt_YUV_space(Re, Ge, Be)
+            _, Uw, Vw = cvt_YUV_space(Rw, Gw, Bw)
+            
+            # Computation of the gradients
+            
+            Gradn[j, i] = compute_gradient(L, Un, Vn, i, j, "north")
+            Grads[j, i] = compute_gradient(L, Us, Vs, i, j, "south")
+            Grade[j, i] = compute_gradient(L, Ue, Ve, i, j, "east")
+            Gradw[j, i] = compute_gradient(L, Uw, Vw, i, j, "west")
+            
+            low_wn[j, i] = 1 / (Gradn[j, i] + eps)
+            low_ws[j, i] = 1 / (Grads[j, i] + eps)
+            low_we[j, i] = 1 / (Grade[j, i] + eps)
+            low_ww[j, i] = 1 / (Gradw[j, i] + eps)
+            
+            W[j, i] = low_wn[j, i] + low_ws[j, i] + low_we[j, i] + low_ww[j, i]
+            
+            low_wn[j, i] = low_wn[j, i] / W[j, i]
+            low_ws[j, i] = low_ws[j, i] / W[j, i]
+            low_we[j, i] = low_we[j, i] / W[j, i]
+            low_ww[j, i] = low_ww[j, i] / W[j, i]
+            
+            R[j, i] = low_wn[j, i]*Rn[j, i] + low_ws[j, i]*Rs[j, i] + low_we[j, i]*Re[j, i] + low_ww[j, i]*Rw[j, i]
+            G[j, i] = low_wn[j, i]*Gn[j, i] + low_ws[j, i]*Gs[j, i] + low_we[j, i]*Ge[j, i] + low_ww[j, i]*Gw[j, i]
+            B[j, i] = low_wn[j, i]*Bn[j, i] + low_ws[j, i]*Bs[j, i] + low_we[j, i]*Be[j, i] + low_ww[j, i]*Bw[j, i]
+    
+    return R, G, B
\ No newline at end of file

src/methods/fortin_come/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.fortin_come.DemosaickingInterChannel import local_int
+
+
+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)
+    z = op.adjoint(y)
+    
+    R0 = z[:, :, 0]
+    G0 = z[:, :, 1]
+    B0 = z[:, :, 2]
+    
+    beta = 0.9 # Inter-channel parameter
+    L = 5 # Size of the local neighborhood
+    eps = 0.00000001 # Tresholding parameter > 0
+    
+    R1, G1, B1 = local_int(R0, G0, B0, beta, L, eps)
+    
+    img_demosaicked = np.zeros((R1.shape[0], R1.shape[1], 3))
+    img_demosaicked[:, :, 0] = R1
+    img_demosaicked[:, :, 1] = G1
+    img_demosaicked[:, :, 2] = B1
+    
+
+    return img_demosaicked
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: FORTIN Côme

src/methods/girardey/methods.py

+import numpy as np
+from src.forward_model import CFA
+
+def HQ_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] = varying_kernel_convolution_HQ(y[:, :], K_red_list)
+        res[:, :, 1] = varying_kernel_convolution_HQ(y[:, :], K_green_list)
+        res[:, :, 2] = varying_kernel_convolution_HQ(y[:, :], K_blue_list)
+ 
+    else:
+        res = np.empty(op.input_shape)
+
+        bayer_y = quad_bayer_to_bayer(y)
+        res[:, :, 0] = varying_kernel_convolution_HQ(bayer_y[:, :], K_red_list)
+        res[:, :, 1] = varying_kernel_convolution_HQ(bayer_y[:, :], K_green_list)
+        res[:, :, 2] = varying_kernel_convolution_HQ(bayer_y[:, :], K_blue_list)
+
+    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
+
+def varying_kernel_convolution_HQ(M, K_list):
+
+    res = np.zeros_like(M)
+    M_padded = np.pad(M,2,"constant",constant_values=0)
+    N_i, N_j = M_padded.shape
+    for i in range(2,N_i-2):
+        for j in range(2,N_j-2):
+            res[i-2, j-2] = np.sum(extract_padded(M_padded, 5, i ,j) * K_list[2* (i % 2) + j % 2])
+
+    np.clip(res, 0, 1, res)
+
+    return res
+
+def quad_bayer_to_bayer(y):
+    res = y
+    N_i, N_j = y.shape
+    for j in range(int(N_j/4)):
+        buff = res[:,j*4+1]
+        res[:,j*4+1] = res[:,j*4+2]
+        res[:,j*4+2] = buff
+    for i in range(int(N_i/4)):
+        buff = res[i*4+1,:]
+        res[i*4+1,:] = res[i*4+2,:]
+        res[i*4+2,:] = buff
+    for i in range(int(N_i/4)):
+        for j in range(int(N_j/4)):
+            buff = res[i*4+2,j*4+1]
+            res[i*4+2,j*4+1] = res[i*4+1,j*4+2]
+            res[i*4+1,j*4+2] = buff   
+    return res
+
+
+
+K_identity = np.zeros((5,5))
+K_identity[2,2] = 1
+
+K_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
+
+K_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
+K_r_and_b_g_brow_rcol = K_r_and_b_g_rrow_bcol.T
+K_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
+
+K_green_list = [K_identity, 
+                K_g_r_and_b,
+                K_g_r_and_b,
+                K_identity]
+
+K_red_list = [K_r_and_b_g_rrow_bcol,
+              K_identity,
+              K_r_b,
+              K_r_and_b_g_brow_rcol]
+
+K_blue_list = [K_r_and_b_g_brow_rcol,
+               K_r_b,
+               K_identity,
+               K_r_and_b_g_rrow_bcol]
\ No newline at end of file

src/methods/girardey/reconstruct.py

+
+import numpy as np
+
+from src.forward_model import CFA
+from src.methods.girardey.methods import HQ_interpolation
+
+
+def run_reconstruction_HQ(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 = HQ_interpolation(op, y)
+
+    return res

src/methods/inssaf_cherki/other_functions.py

+import numpy as np
+from scipy import signal
+from src.forward_model import CFA
+
+def compute_gradient(channel, points):
+    """ Compute the gradient for the channel based on the specified points. """
+    gradient = np.zeros_like(channel)
+    for point in points:
+        gradient += np.roll(channel, shift=point, axis=(0, 1))
+    gradient /= len(points)
+    gradient -= channel
+    return gradient
+
+def gradient_correction_interpolation(op: CFA, z: np.ndarray,alpha=0.5, beta=0.5, gamma=0.5) -> np.ndarray:
+    # defining bilinear interpolation filters for the different channels (bayer pattern case)
+    bilinear_filter_red = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 4
+    bilinear_filter_green = np.array([[0, 1, 0], [1, 4, 1], [0, 1, 0]]) / 4
+    bilinear_filter_blue = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 4
+
+    interpolated_img = np.empty(op.input_shape)
+    
+    # applying the convolution product with the bilinear filters
+    interpolated_img[:, :, 0] = signal.convolve2d(z[:, :, 0], bilinear_filter_red,  boundary='symm',mode='same')
+    interpolated_img[:, :, 1] = signal.convolve2d(z[:, :, 1], bilinear_filter_green,  boundary='symm',mode='same')
+    interpolated_img[:, :, 2] = signal.convolve2d(z[:, :, 2], bilinear_filter_blue, boundary='symm', mode='same')
+
+    """ Applying gradient correction to the bilinearly interpolated image. """
+
+
+    R = interpolated_img[:, :, 0]
+    G = interpolated_img[:, :, 1]
+    B = interpolated_img[:, :, 2]
+
+    # computing gradients for R and B channels
+    points_R = [(0, -2), (0, 2), (-2, 0), (2, 0)] 
+    points_B = [(-1, -1), (-1, 1), (1, -1), (1, 1), (0, 0)]  
+    grad_R = compute_gradient(R, points_R)
+    grad_B = compute_gradient(B, points_B)
+
+    G_corrected = np.copy(G)
+    G_corrected[0::2, 0::2] = G[0::2, 0::2] + alpha * grad_R[0::2, 0::2]
+    
+    G_corrected[1::2, 1::2] = G[1::2, 1::2] + gamma * grad_B[1::2, 1::2]
+    
+    # correction of R at G locations and B at G locations
+    R_corrected = np.copy(R)
+    B_corrected = np.copy(B)
+    R_corrected[1::2, 0::2] = R[1::2, 0::2] + beta * grad_R[1::2, 0::2]  
+    R_corrected[0::2, 1::2] = R[0::2, 1::2] + beta * grad_R[0::2, 1::2] 
+    B_corrected[1::2, 0::2] = B[1::2, 0::2] + beta * grad_B[1::2, 0::2] 
+    B_corrected[0::2, 1::2] = B[0::2, 1::2] + beta * grad_B[0::2, 1::2]  
+
+   
+    corrected_image = np.stack((R_corrected, G_corrected, B_corrected), axis=-1)
+
+    return corrected_image
+
+def swapping(op, y):
+    """
+    Convert both the CFA operator's mask and an image from Quad Bayer to Bayer by swapping method.
+
+    Args:
+        op: CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        tuple: A tuple containing the updated CFA operator and the converted image.
+    """
+
+    if op.cfa == 'quad_bayer':
+        #  swapping 2 columns every 2 columns
+        op.mask[:, 1::4], op.mask[:, 2::4] = op.mask[:, 2::4].copy(), op.mask[:, 1::4].copy()
+        #  swapping 2 lines every 2 lines
+        op.mask[1::4, :], op.mask[2::4, :] = op.mask[2::4, :].copy(), op.mask[1::4, :].copy()
+        # swapping the diagonal pairs
+        op.mask[1::4, 1::4], op.mask[2::4, 2::4] = op.mask[2::4, 2::4].copy(), op.mask[1::4, 1::4].copy()
+
+    converted_image = y.copy()
+    converted_image[:, 1::4], converted_image[:, 2::4] = converted_image[:, 2::4].copy(), converted_image[:, 1::4].copy()
+    converted_image[1::4, :], converted_image[2::4, :] = converted_image[2::4, :].copy(), converted_image[1::4, :].copy()
+    converted_image[1::4, 1::4], converted_image[2::4, 2::4] = converted_image[2::4, 2::4].copy(), converted_image[1::4, 1::4].copy()
+
+    # updating to bayer pattern
+    op.cfa = 'bayer'
+
+    return op, converted_image    
\ No newline at end of file

src/methods/inssaf_cherki/reconstruct.py

+import numpy as np
+from src.methods.inssaf_cherki.other_functions import gradient_correction_interpolation, swapping, compute_gradient
+from src.forward_model import CFA
+
+
+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)
+
+    if op.cfa == 'quad_bayer':
+        op, y = swapping(op, y)   
+
+    z = op.adjoint(y)
+    reconstructed_image = gradient_correction_interpolation(op,z,alpha=0.5, beta=0.5, gamma=0.5)
+
+    return reconstructed_image
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/khattabi/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.khattabi.utils import HA
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    # Performing the reconstruction.
+    # TODO
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+    y_adjoint=op.adjoint(y)
+    res = HA(y_adjoint)
+
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/khattabi/utils.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 HA(y_adjoint):
+    height, width = y_adjoint.shape[:2]
+    # Green channel interpolation
+
+
+    for i in range(2 , height-2): 
+        for j in range(2, width-2):
+            # Red pixels
+            if y_adjoint[i,j,1] == 0 and y_adjoint[i,j,0] != 0:  
+                grad_y = np.abs(y_adjoint[i-1,j,1] - y_adjoint[i+1,j,1]) + np.abs(2*y_adjoint[i,j,0] - y_adjoint[i-2,j,0] - y_adjoint[i+2,j,0])
+                grad_x = np.abs(y_adjoint[i,j-1,1] - y_adjoint[i,j+1,1]) + np.abs(2*y_adjoint[i,j,0] - y_adjoint[i,j-2,0] - y_adjoint[i,j+2,0])
+
+                if grad_x < grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/2 + (2*y_adjoint[i,j,0] - y_adjoint[i,j-2,0] - y_adjoint[i,j+2,0])/4
+                elif  grad_x > grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1])/2 + (2*y_adjoint[i,j,0] - y_adjoint[i-2,j,0] - y_adjoint[i+2,j,0])/4
+                else:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1] + y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/4 + \
+                                       (2*y_adjoint[i,j,0] - y_adjoint[i,j-2,0] - y_adjoint[i,j+2,0] + 2*y_adjoint[i,j,0] - y_adjoint[i-2,j,0] - y_adjoint[i+2,j,0])/8
+            # Blue Pixels
+            elif y_adjoint[i,j,1] == 0 and y_adjoint[i,j,2] != 0:   
+                grad_y = np.abs(y_adjoint[i-1,j,1] - y_adjoint[i+1,j,1]) + np.abs(2*y_adjoint[i,j,2] - y_adjoint[i-2,j,2] - y_adjoint[i+2,j,2])
+                grad_x = np.abs(y_adjoint[i,j-1,1] - y_adjoint[i,j+1,1]) + np.abs(2*y_adjoint[i,j,2] - y_adjoint[i,j-2,2] - y_adjoint[i,j+2,2])
+
+                if grad_x < grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/2 + (2*y_adjoint[i,j,2] - y_adjoint[i,j-2,2] - y_adjoint[i,j+2,2])/4
+                elif  grad_x > grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1])/2 + (2*y_adjoint[i,j,2] - y_adjoint[i-2,j,2] - y_adjoint[i+2,j,2])/4
+                else:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1] + y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/4 + \
+                                       (2*y_adjoint[i,j,2] - y_adjoint[i,j-2,2] - y_adjoint[i,j+2,2] + 2*y_adjoint[i,j,2] - y_adjoint[i-2,j,2] - y_adjoint[i+2,j,2])/8
+                    
+    # Red and blue channel interpolation
+    for i in range(1 , height-1):  
+        for j in range(1, width-1):
+            if y_adjoint[i,j,2] != 0 :
+                 y_adjoint[i,j,0] = (y_adjoint[i-1,j-1,0] + y_adjoint[i-1,j+1,0] + y_adjoint[i+1,j-1,0] + y_adjoint[i+1,j+1,0]) / 4
+
+            elif y_adjoint[i,j,0] != 0:
+                 y_adjoint[i,j,2] = (y_adjoint[i-1,j-1,2] + y_adjoint[i-1,j+1,2] + y_adjoint[i+1,j-1,2] + y_adjoint[i+1,j+1,2]) / 4
+            
+            else:
+                y_adjoint[i,j] = y_adjoint[i,j]
+
+    for i in range(1 , height-1):
+        for j in range(1, width-1):    
+            if y_adjoint[i,j,0] == y_adjoint[i,j,2]:
+                y_adjoint[i,j,0] = (y_adjoint[i-1,j,0] + y_adjoint[i,j-1,0] + y_adjoint[i+1,j,0] + y_adjoint[i,j+1,0]) / 4             
+                y_adjoint[i,j,2] = (y_adjoint[i-1,j,2] + y_adjoint[i,j-1,2] + y_adjoint[i+1,j,2] + y_adjoint[i,j+1,2]) / 4  
+            
+    return y_adjoint
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/laporte_auriane/fct_to_demosaick.py

+"""This file contains all the functions that are used in order to perform a Kimmel algorithm demosaicking"""
+
+##### Importations #####
+import numpy as np
+import matplotlib.pyplot as plt
+from src.forward_model import CFA
+
+##### Some functions used #####
+def neigh(i, j, colour_chan, img):
+    
+    """Finds the neighbours of the pixel (i, j) of a colour channel
+    Args:
+        i: column of the pixel
+        j: line of the pixel
+        colour_chan: colour channel chosen
+        img: image
+    Returns:
+        List of the neighbours of the pixel (i, j) and the pixel (i, j)
+    """
+    
+    X = img[:, :, colour_chan].shape[0]
+    Y = img[:, :, colour_chan].shape[1]
+    
+    # Place each Pi compared to P5 (i horizontal and j vertical)
+    # See P1 to P9 drawing, page 2 of http://ultra.sdk.free.fr/docs/Image-Processing/Demosaic/An%20Improved%20Demosaicing%20Algorithm.pdf
+    # % (modulo) to make sure it is a multiple (avoid error index 1024 out of bounds)
+    P1 = img[(i-1)%Y, (j-1)%X, colour_chan]
+    P2 = img[i%Y, (j-1)%X, colour_chan]
+    P3 = img[(i+1)%Y, (j-1)%X, colour_chan]
+    P4 = img[(i-1)%Y, j%X, colour_chan]
+    P5 = img[i%Y, j%X, colour_chan]
+    P6 = img[(i+1)%Y, j%X, colour_chan]
+    P7 = img[(i-1)%Y, (j+1)%X, colour_chan]
+    P8 = img[i%Y, (j+1)%X, colour_chan]
+    P9 = img[(i+1)%Y, (j+1)%X, colour_chan]
+    #print(f'Size of P1: {P1.size}') # 1
+    
+    return [P1, P2, P3, P4, P5, P6, P7, P8, P9]
+
+def derivat(neighbour_list):
+    
+    """Compute directional derivatives of pixel
+    Args:
+        neighbour_list: list of neighbours
+    Returns:
+        List of the directional derivatives
+    """
+    
+    # Must be 9 (9P neigbours)
+    #print(f'neighbour_list len to derivate: {len(neighbour_list)}') # 9
+    
+    # Assign P neighbours to the neighbour_list (to then compute the derivatives)
+    [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbour_list
+    
+    # Compute the derivatives (independant of the colour component)
+    Dx = (P4 - P6) / 2
+    Dy = (P2 - P8) / 2
+    Dxd = (P3 - P7) / (2 * np.sqrt(2))
+    Dyd = (P1 - P9) / (2 * np.sqrt(2))
+    # Formula given at the top of page 3 (but very long computing time, and does not provide better results)
+    #Dxd = np.max( ( np.abs((P3 - P5) / (np.sqrt(2))), np.abs((P7 - P5) / (np.sqrt(2))) ) )
+    #Dyd = np.max( ( np.abs((P1 - P5) / (np.sqrt(2))), np.abs((P9 - P5) / (np.sqrt(2))) ) )
+    
+    # Store the computed values in deriv_list
+    deriv_list = [Dx, Dy, Dxd, Dyd]
+    #print(f'Len de deriv_list: {len(deriv_list)}') # 4
+    
+    return deriv_list
+
+def weight_fct(i, j, neighbour_list, deriv_list, colour_chan, img):
+    
+    """Compute the weights Ei
+    Args:
+        i: column of the pixel
+        j: line of the pixel
+        neighbour_list: neighbour_list: list of neighbours
+        deriv_list: list of derivatives
+        colour_chan:colour channel chosen
+        img: image
+    Returns:
+        List of the weights
+    """
+    
+    # Assignment
+    [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbour_list
+    [Dx, Dy, Dxd, Dyd] = deriv_list
+    X = img[:,:,colour_chan].shape[0]
+    Y = img[:,:,colour_chan].shape[1]
+    
+    # E list to complete (weighting function)
+    E = []
+    
+    # Fix the start
+    cur = 1
+    
+    # Neighbourhood
+    for step in range(-1, 2): # from -1 to 1 (2-1)
+        for step in range(-1, 2): # otherwise only 3 values in E
+            # Find the neigbours
+            n = neigh(i+step, j+step, colour_chan, img)
+            # Derivatives
+            d = derivat(n)       
+            
+            # See P1 to P9 drawing, page 2 of http://ultra.sdk.free.fr/docs/Image-Processing/Demosaic/An%20Improved%20Demosaicing%20Algorithm.pdf
+            if cur==4 or cur==6:
+                E.append( 1 / np.sqrt(1 + Dx**2 + d[0]**2) )
+            elif cur==2 or cur==8:
+                E.append( 1 / np.sqrt(1 + Dy**2 + d[1]**2) )
+            elif cur==7 or cur==3:
+                E.append( 1 / np.sqrt(1 + Dxd**2 + d[2]**2) )
+            else: # 1 or 9
+                E.append( 1 / np.sqrt(1 + Dyd**2 + d[3]**2) )
+            cur = cur+1
+
+    #print(f'Len of E: {len(E)}') # 9
+    # other way to code (Page 3/6 of http://elynxsdk.free.fr/ext-docs/Demosaicing/more/news0/Optimal%20Recovery%20Demosaicing.pdf, method 1)), but I did not understand what is or how T was chosen?
+    
+    return E
+
+def clip_extreme(img):
+    
+    """Clip the values of the image, so they are between 0-1
+    Args:
+        img: image
+    Returns:
+        The clipped image
+    """
+
+    #max_val = np.amax(img)
+    #min_val = np.amin(img)
+    #print(f'Max value before clipping: {max_val}')
+    #print(f'Min value before clipping: {min_val}')
+    
+    # Put extreme values between 0 and 1
+    img = np.clip(img, 0, 1)
+    
+    return img
+ 
+def quad_to_bayer(img, op):
+    
+    """Performs a conversion from Quad_bayer to Bayer pattern
+    Args:
+        img: image
+        op: CFA op
+    Returns:
+        The image Bayer pattern
+    """
+    
+    # Based on p28/32 https://pyxalis.com/wp-content/uploads/2021/12/PYX-ImageViewer-User_Guide.pdf
+    ### 1. Swap 2 col every 2 columns
+    # Start: 2nd col --> 1 (bc start counting at 0) / Step: 4 (look at drawing of 1st swap p28 Pyxalis) / End: input_shape[0]-2 --> input_shape[0]-1
+    for j in range (1, img.shape[0]-1, 4):
+        store = np.copy(img[:, j]) # Store col j of the img
+        img[:, j] =  np.copy(img[:, j+1])  # Col j becomes like col j+1
+        img[:, j+1] = store # Col j+1 become like col j (previously saved otherwise value lost)
+        store_op = np.copy(op.mask[:, j])
+        op.mask[:, j] =  np.copy(op.mask[:, j+1])
+        op.mask[:, j+1] = store_op
+        
+    ### 2. Swap 2 lines every 2 lines (same process for the lines)
+    for i in range (1, img.shape[1]-1, 4):
+        store = np.copy(img[i, :])
+        img[i, :] =  np.copy(img[i+1, :])
+        img[i+1, :] = store
+        store_op = np.copy(op.mask[i, :])
+        op.mask[i, :] =  np.copy(op.mask[i+1, :])
+        op.mask[i+1, :] = store_op
+
+    ### 3. Swap back some diagonal greens
+    # Starts: 0 and 2 / Steps: 4 / End: img.shape[nb]-1 --> img.shape[nb]
+    for k in range(0, img.shape[0], 4):
+        for l in range(2, img.shape[1], 4):
+            store = np.copy(img[k, l])
+            img[k, l] = np.copy(img[k+1, l-1])
+            img[k+1, l-1] = store
+            store_op = op.mask[k, l]
+            op.mask[k, l] = np.copy(op.mask[k+1, l-1])
+            op.mask[k+1, l-1] = store_op
+            
+    return img
+
+
+##### 1. Interpolation of green colour #####
+def interpol_green(neighbour_list, weights_list):
+    
+    """Interpolates missing green pixel 
+    Args:
+        neighbour_list: list of neighbours
+        weights_list: list of weights
+    Returns:
+        The interpolated pixel
+    """
+    
+    # Assignment
+    [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbour_list
+    [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+    
+    # Compute missing green pixel (combination of 4 nearest neighbours)
+    G5 = (E2*P2 + E4*P4 + E6*P6 + E8*P8) / (E2+E4+E6+E8)
+    #print(f"G5 size: {G5.size}") # 1
+    
+    return G5
+
+
+##### 2. Interpolation of red and blue colours #####
+def interpol_red_blue(neighbour_list, weights_list, green_neighbour):
+    
+    """Interpolates missing blue/red pixel
+    Args:
+        neighbour_list: list of neighbours in blue/red channel
+        weights_list: list of weights
+        green_neighbour: list of the neighbours in green channel
+    Returns:
+        The interpolated pixel (in blue/red channel)
+    """
+    
+    # Assignment
+    [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbour_list
+    [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+    [G1, G2, G3, G4, G5, G6, G7, G8, G9] = green_neighbour
+    
+    # Compute missing red/blue pixel (combination of 4 nearest neighbours)
+    R5_num = E1*(P1/G1) + E3*(P3/G3) + E7*(P7/G7) + E9*(P9/G9)
+    R5_denom = E1+E3+E7+E9
+    R5 = G5 * R5_num / R5_denom
+    #print(f"R5 size: {R5.size}") # 1
+    
+    return R5
+
+
+##### 3. Correction stage #####
+def green_correction(red_neighbour, green_neighbour, blue_neighbour, weights_list):
+    
+    """Correction of green pixels
+    Args:
+        red_neighbour: neighbours in red channel
+        green_neighbour: neighbours in green channel
+        blue_neighbour: neighbours in blue channel
+        weights_list: list of weights
+    Returns:
+        Corrected green pixel
+    """
+    
+    # Assignment
+    [G1, G2, G3, G4, G5, G6, G7, G8, G9] = green_neighbour
+    [R1, R2, R3, R4, R5, R6, R7, R8, R9] = red_neighbour
+    [B1, B2, B3, B4, B5, B6, B7, B8, B9] = blue_neighbour
+    [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+
+    # Compute correction
+    GB5_num = E2*(G2/B2) + E4*(G4/B4) + E6*(G6/B6) + E8*(G8/B8)
+    GR5_num = E2*(G2/R2) + E4*(G4/R4) + E6*(G6/R6) + E8*(G8/R8)
+    G5_denom = E2+E4+E6+E8
+    GB5 = B5 * GB5_num / G5_denom
+    GR5 = R5 * GR5_num / G5_denom
+    
+    #print(f"GB5 size: {GB5.size}") # 1
+    print(f"GB5: {GB5}")
+    #print(f"GR5 size: {GR5.size}") # 1
+    print(f"GR5: {GR5}")
+    
+    G5 = (GR5 + GB5) / 2
+    print(f"G5 size: {G5.size}") # 1
+    print(f"G5: {G5}")
+    
+    return G5
+
+def blue_correction(green_neighbour, blue_neighbour, weights_list):
+    
+    """Correction of blue pixels
+    Args:
+        green_neighbour: neighbours in green channel
+        blue_neighbour: neighbours in blue channel
+        weights_list: list of weights
+    Returns:
+        Corrected blue pixel
+    """
+    
+    # Assignment
+    [G1, G2, G3, G4, G5, G6, G7, G8, G9] = green_neighbour
+    [B1, B2, B3, B4, B5, B6, B7, B8, B9] = blue_neighbour
+    [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+
+    # Compute correction
+    B5_num = 0
+    B5_denom = 0
+    print(f'Len de weights_list: {len(weights_list)}')
+    
+    for i in range(len(weights_list)):
+        if i != 5:
+            B5_num += weights_list[i] * blue_neighbour[i] / green_neighbour[i]
+            B5_denom += weights_list[i]
+    
+    B5 = G5 * B5_num / B5_denom
+    
+    #print(f"B5_denom size: {B5_denom.size}") # 1
+    print(f"B5_denom: {B5_denom}") # number
+    #print(f"B5_num size: {B5_num.size}") # 1
+    print(f"B5_num: {B5_num}") # nan
+    #print(f"B5 size: {B5.size}") # 1
+    print(f"B5: {B5}")
+    
+    return B5
+
+def red_correction(red_neighbour, green_neighbour, weights_list):
+    
+    """Correction of red pixels
+    Args:
+        red_neighbour: neighbours in red channel
+        green_neighbour: neighbours in green channel
+        weights_list: list of weights
+    Returns:
+        Corrected red pixel
+    """
+    
+    # Assignment
+    [G1, G2, G3, G4, G5, G6, G7, G8, G9] = green_neighbour
+    [R1, R2, R3, R4, R5, R6, R7, R8, R9] = red_neighbour
+    [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+
+    # Compute correction
+    R5_num = 0
+    R5_denom = 0
+    for i in range(len(weights_list)):
+        if i != 5:
+        
+            R5_num += weights_list[i] * red_neighbour[i] / green_neighbour[i]
+            R5_denom += weights_list[i]
+    
+    R5 = G5 * R5_num / R5_denom
+
+    return R5
\ No newline at end of file

src/methods/laporte_auriane/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.laporte_auriane.fct_to_demosaick import *
+
+
+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.
+    """
+    
+    #print(f'y shape: {y.shape}') # (1024, 1024)
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+    
+    ########## TODO ##########
+    
+     # Colour channels
+    red = 0
+    green = 1
+    blue = 2
+    correction = 0
+    
+    # In case of Quad_bayer op
+    if op.cfa == 'quad_bayer':
+        print("Transformation into Bayer pattern")
+        y = quad_to_bayer(y, op)
+        print('Quad_bayer to Bayer done\n\nDemoisaicking processing...')
+ 
+    z = op.adjoint(y)
+    # 1st and 2nd dim of a colour channel
+    X = z[:,:,0].shape[0]
+    Y = z[:,:,0].shape[1]
+    #print(f"X: {X}") # 1024
+    #print(f"Y: {Y}") # 1024
+    
+    
+    ### Kimmel algorithm ###
+    # About 3 to 4 minutes
+    
+    
+    # Interpolation of green colour
+    for i in range (X):
+        for j in range (Y):
+            # Control if value 0 in green channel
+            if z[i,j,1] == 0:
+
+                # Find neigbours
+                n = neigh(i, j, green, z)
+                # Derivatives
+                d = derivat(n)
+                # Weighting fct
+                w = weight_fct(i, j, n, d, green, z)
+                # Green interpolation (cross)
+                z[i,j,1] = interpol_green(n, w)
+    # Clip (avoid values our of range 0-1)
+    z = clip_extreme(z)
+    
+    
+    # 2 steps for the blues
+    # Interpolate missing blues at the red locations
+    for i in range (0, X, 2):
+        for j in range (1, Y, 2):
+            
+            # Find neigbours + green
+            n = neigh(i, j, blue, z)
+            n_G = neigh(i, j, green, z)
+            d = derivat(n_G)
+            w = weight_fct(i, j, n, d, green, z)
+            # Blue interpolation (4 corners)
+            z[i,j,2] = interpol_red_blue(n, w, n_G)
+    # Interpolate missing blues at the green locations
+    for i in range (X):
+        for j in range (Y):
+            # Control if value 0 in blue channel
+            if z[i,j,2] == 0:
+
+                # Find neigbours
+                n_B = neigh(i, j, blue, z)
+                d = derivat(n_B)
+                w = weight_fct(i, j, n_B, d, blue, z)
+                # Blue interpolation (cross)
+                z[i,j,2] = interpol_green(n_B,w)
+    z = clip_extreme(z)
+
+
+    # 2 steps for the reds
+    # Interpolate missing reds at the blue locations
+    for i in range (1, X, 2):
+        for j in range (0, Y, 2):
+            
+            # Find neigbours + green
+            n = neigh(i, j, red, z)
+            n_G = neigh(i, j, green, z)
+            d = derivat(n_G)
+            w = weight_fct(i, j, n_G, d, green, z)
+            # Red interpolation (4 corners)
+            z[i,j,0] = interpol_red_blue(n, w, n_G)
+    # Interpolate missing reds at the green locations
+    for i in range (X):
+        for j in range (Y):
+            # Control if value 0 in red channel
+            if z[i,j,0] == 0:
+
+                # Find neigbours
+                n_R = neigh(i, j, red, z)
+                d = derivat(n_R)
+                w = weight_fct(i, j, n_R, d, red, z)
+                # Red interpolation (cross)
+                z[i,j,0] = interpol_green(n_R,w)
+    z = clip_extreme(z)
+    
+    
+    # Correction step (repeated 3 times)
+    if correction == 1:
+        for i in range(3):
+            
+            n_G[4] = green_correction(n_R, n_G, n_B, w)
+            n_B[4] = blue_correction(n_G, n_B, w)
+            n_R[4] = red_correction(n_R, n_G, w)
+            # nan
+    
+    #print(f'Image reconstructed shape: {z.shape}') # 1024, 1024, 3
+
+    return z
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller