src/methods/CASA_NOVA_Vincent/functions.py

+from scipy import ndimage
+import numpy as np
+
+############################################################
+
+    
+def color_pixel(i,j,cfa = "bayer"):
+    if (cfa == "quad_bayer"):
+        i = i//2
+        j = j//2
+    
+    if ((i+j)%2==0):
+        return 'green'
+    else:
+        if (i%2==0):
+            return 'red'
+        else:
+            return 'blue'
+
+def rmse_pixel(pixel_raw,pixel_extrapolate):
+    return np.sqrt(np.mean((pixel_raw-pixel_extrapolate)**2))
+
+######### Method extrapolation with edge detection #########
+def compute_orientation_matrix(img_raw):
+    vertical = ndimage.sobel(img_raw, 0)
+    horizontal = ndimage.sobel(img_raw, 1)
+    orientation_matrix = np.zeros(img_raw.shape)
+    orientation_matrix[vertical < horizontal] = 1
+    return orientation_matrix
+
+## Green Channel ##
+##Formulas for etrapolation of pixels:
+def extrapolate_green_top(img_raw,i,j):
+    return img_raw[i-1,j] + 3/4*(img_raw[i,j]-img_raw[i-2,j])-1/4*(img_raw[i-1,j]-img_raw[i-3,j])
+
+def extrapolate_green_bottom(img_raw,i,j):
+    return img_raw[i+1,j] + 3/4*(img_raw[i,j]-img_raw[i+2,j])-1/4*(img_raw[i+1,j]-img_raw[i+3,j])
+
+def extrapolate_green_left(img_raw,i,j):
+    return img_raw[i,j-1] + 3/4*(img_raw[i,j]-img_raw[i,j-2])-1/4*(img_raw[i,j-1]-img_raw[i,j-3])
+
+def extrapolate_green_right(img_raw,i,j):
+    return img_raw[i,j+1] + 3/4*(img_raw[i,j]-img_raw[i,j+2])-1/4*(img_raw[i,j+1]-img_raw[i,j+3])
+
+## Extrapolation method:
+def median_extrapolate_green_pixel(img_raw,i,j,orientations_to_drop):
+    list_extrapolate_pixel = []
+
+    if ("top" not in orientations_to_drop):
+        list_extrapolate_pixel.append(extrapolate_green_top(img_raw,i,j))
+    if ("bottom" not in orientations_to_drop):
+        list_extrapolate_pixel.append(extrapolate_green_bottom(img_raw,i,j))  
+    if("left" not in orientations_to_drop):
+        list_extrapolate_pixel.append(extrapolate_green_left(img_raw,i,j))
+    if("right" not in orientations_to_drop):
+        list_extrapolate_pixel.append(extrapolate_green_right(img_raw,i,j))
+
+    return np.median(list_extrapolate_pixel)
+
+def extrapolate_green_pixel(img_raw,i,j,orientation):
+    # First the borders:
+    orientations_to_drop = []
+    if (i<2):
+        orientations_to_drop.append('top')
+    if (i>img_raw.shape[0]-4):
+        orientations_to_drop.append('bottom')
+    if (j<2):
+        orientations_to_drop.append('left')
+    if (j>img_raw.shape[1]-4):
+        orientations_to_drop.append('right')
+    
+    
+    # Then the rest of the image:
+    else:
+        if (orientation == 1): # V < H so we gonna eliminate one horizontal pixel. 
+            if ("right" not in orientations_to_drop and "left" not in orientations_to_drop):
+                rmse_pixel_left = rmse_pixel(img_raw[i,j],extrapolate_green_left(img_raw,i,j))
+                rmse_pixel_right = rmse_pixel(img_raw[i,j],extrapolate_green_right(img_raw,i,j))
+                if (rmse_pixel_left > rmse_pixel_right):
+                    orientations_to_drop.append('left')
+                else:
+                    orientations_to_drop.append('right')
+    
+        else: # V > H so we gonna eliminate one vertical pixel.
+            if ("top" not in orientations_to_drop and "bottom" not in orientations_to_drop):
+                rmse_pixel_top = rmse_pixel(img_raw[i,j],extrapolate_green_top(img_raw,i,j))
+                rmse_pixel_bottom = rmse_pixel(img_raw[i,j],extrapolate_green_bottom(img_raw,i,j))
+                if (rmse_pixel_top > rmse_pixel_bottom):
+                    orientations_to_drop.append('top')
+                else:
+                    orientations_to_drop.append('bottom')
+    return median_extrapolate_green_pixel(img_raw,i,j,orientations_to_drop)
+
+def extrapolate_green(img_raw,extrapolate_img):
+    orientation_matrix = compute_orientation_matrix(img_raw)
+    for i in range(img_raw.shape[0]):
+        for j in range(img_raw.shape[1]):
+            if (color_pixel(i,j)!= "green"):
+                extrapolate_img[i,j,1] = extrapolate_green_pixel(img_raw,i,j,orientation_matrix[i,j])
+            else:
+                extrapolate_img[i,j,1] = img_raw[i,j]
+    return extrapolate_img
+
+## Red and Blue Channels ##
+def extrapolate_top(img_raw,img_extrapolate,i,j):
+    return (img_raw[i-1,j] + img_raw[i,j]-img_extrapolate[i-1,j,1])
+def extrapolate_left(img_raw,img_extrapolate,i,j):
+    return (img_raw[i,j-1] + img_raw[i,j]-img_extrapolate[i,j-1,1])
+def extrapolate_right(img_raw,img_extrapolate,i,j):
+    return (img_raw[i,j+1] + img_raw[i,j]-img_extrapolate[i,j+1,1])
+def extrapolate_bottom(img_raw,img_extrapolate,i,j):
+    return (img_raw[i+1,j] + img_raw[i,j]-img_extrapolate[i+1,j,1])
+
+def extrapolate_top_left(img_raw,img_extrapolate,i,j):
+    return (img_raw[i-1,j-1] + img_extrapolate[i,j,1]-img_extrapolate[i-1,j-1,1])
+def extrapolate_top_right(img_raw,img_extrapolate,i,j):
+    return (img_raw[i-1,j+1] + img_extrapolate[i,j,1]-img_extrapolate[i-1,j+1,1])
+def extrapolate_bottom_left(img_raw,img_extrapolate,i,j):
+    return (img_raw[i+1,j-1] + img_extrapolate[i,j,1]-img_extrapolate[i+1,j-1,1])
+def extrapolate_bottom_right(img_raw,img_extrapolate,i,j):
+    return (img_raw[i+1,j+1] + img_extrapolate[i,j,1]-img_extrapolate[i+1,j+1,1])
+
+def median_pixel(img_raw,img_extrapolate,i,j,orientations_to_drop):
+    list_extrapolate = []
+    if (color_pixel(i,j) != "green"):
+        if("top_left" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_top_left(img_raw,img_extrapolate,i,j))
+        if("top_right" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_top_right(img_raw,img_extrapolate,i,j))
+        if("bottom_left" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_bottom_left(img_raw,img_extrapolate,i,j))
+        if("bottom_right" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_bottom_right(img_raw,img_extrapolate,i,j))
+
+    elif (color_pixel(i,j) == "green"):
+        if("top" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_top(img_raw,img_extrapolate,i,j))
+        if("left" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_left(img_raw,img_extrapolate,i,j))
+        if("right" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_right(img_raw,img_extrapolate,i,j))
+        if("bottom" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_bottom(img_raw,img_extrapolate,i,j))
+    return np.median(list_extrapolate)
+
+def extrapolate_pixel(img_raw,img_extrapolate,i,j,color):
+    orientations_to_drop = []
+
+    if (color_pixel(i,j)!='green'):
+        if (i<1):
+            orientations_to_drop.append("top_left")
+            orientations_to_drop.append("top_right")
+        if (i>img_raw.shape[0]-2):
+            orientations_to_drop.append("bottom_left")
+            orientations_to_drop.append("bottom_right")
+        if (j<1):
+            orientations_to_drop.append("top_left")
+            orientations_to_drop.append("bottom_left")
+        if (j>img_raw.shape[1]-2):
+            orientations_to_drop.append("top_right")
+            orientations_to_drop.append("bottom_right")
+        
+        
+        if ("top_left" not in orientations_to_drop and "top_right" not in orientations_to_drop and "bottom_left" not in orientations_to_drop and "bottom_right" not in orientations_to_drop):
+            rmse_top_left = rmse_pixel(img_raw[i,j],extrapolate_top_left(img_raw,img_extrapolate,i,j))
+            rmse_top_right = rmse_pixel(img_raw[i,j],extrapolate_top_right(img_raw,img_extrapolate,i,j))
+            rmse_bottom_left = rmse_pixel(img_raw[i,j],extrapolate_bottom_left(img_raw,img_extrapolate,i,j))
+            rmse_bottom_right = rmse_pixel(img_raw[i,j],extrapolate_bottom_right(img_raw,img_extrapolate,i,j))
+            if (rmse_bottom_left> rmse_bottom_right and rmse_bottom_left> rmse_top_left and rmse_bottom_left> rmse_top_right):
+                orientations_to_drop.append("bottom_left")
+            elif (rmse_bottom_right> rmse_bottom_left and rmse_bottom_right> rmse_top_left and rmse_bottom_right> rmse_top_right):
+                orientations_to_drop.append("bottom_right")
+            elif (rmse_top_left> rmse_bottom_left and rmse_top_left> rmse_bottom_right and rmse_top_left> rmse_top_right):
+                orientations_to_drop.append("top_left")
+            else:
+                orientations_to_drop.append("top_right")
+    elif(color_pixel(i,j)=="green"):
+        if (i<1):
+            orientations_to_drop.append("top")
+        if (i>img_raw.shape[0]-2):
+            orientations_to_drop.append("bottom")
+        if (j<1):
+            orientations_to_drop.append("left")
+        if (j>img_raw.shape[1]-2):
+            orientations_to_drop.append("right")
+        
+        if ((i%2!=0 and color == "red") or (i%2==0 and color == "blue")):
+            if ("right" not in orientations_to_drop and "left" not in orientations_to_drop):
+                rmse_pixel_left = rmse_pixel(img_raw[i,j],extrapolate_left(img_raw,img_extrapolate,i,j))
+                rmse_pixel_right = rmse_pixel(img_raw[i,j],extrapolate_right(img_raw,img_extrapolate,i,j))
+                if (rmse_pixel_left > rmse_pixel_right):
+                    orientations_to_drop.append('left')
+                else:
+                    orientations_to_drop.append('right')
+        else:
+            if ("top" not in orientations_to_drop and "bottom" not in orientations_to_drop):
+                rmse_pixel_top = rmse_pixel(img_raw[i,j],extrapolate_top(img_raw,img_extrapolate,i,j))
+                rmse_pixel_bottom = rmse_pixel(img_raw[i,j],extrapolate_bottom(img_raw,img_extrapolate,i,j))
+                if (rmse_pixel_top > rmse_pixel_bottom):
+                    orientations_to_drop.append('top')
+                else:
+                    orientations_to_drop.append('bottom')
+    return median_pixel(img_raw,img_extrapolate,i,j,orientations_to_drop)
+
+def extrapolate_red(img_raw,img_extrapolate):
+    for i in range(img_raw.shape[0]):
+        for j in range(img_raw.shape[1]):
+            if (color_pixel(i,j)!="red"):
+                img_extrapolate[i,j,0] = extrapolate_pixel(img_raw,img_extrapolate,i,j,"red")
+            else:
+                img_extrapolate[i,j,0] = img_raw[i,j]
+
+def extrapolate_blue(img_raw,img_extrapolate):
+    for i in range(img_raw.shape[0]):
+        for j in range(img_raw.shape[1]):
+            if (color_pixel(i,j)!="blue"):
+                img_extrapolate[i,j,2] = extrapolate_pixel(img_raw,img_extrapolate,i,j,"blue")
+            else:
+                img_extrapolate[i,j,2] = img_raw[i,j]
+
+def extrapolate_img(img_cfa):
+    extapolate_img = np.zeros(img_cfa.shape + (3,))
+    extrapolate_green(img_cfa,extapolate_img)
+    extrapolate_red(img_cfa,extapolate_img)
+    extrapolate_blue(img_cfa,extapolate_img)
+    return extapolate_img
+
+
+#################################################
+##QUAD BAYER
+#################################################
+
+
+## Green Channel ##
+
+### Formulas extrapolation of pixels:
+def extrapolate_green_top_quad(img_raw,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] = img_raw[i+m-1*2,j+n] + 3/4*(img_raw[i+m,j+n]-img_raw[i+m-2*2,j+n])-1/4*(img_raw[i+m-1*2,j+n]-img_raw[i+m-3*2,j+n])
+    return extrapolate_quad
+
+def extrapolate_green_bottom_quad(img_raw,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] = img_raw[i+m+1*2,j+n] + 3/4*(img_raw[i+m,j+n]-img_raw[i+m+2*2,j+n])-1/4*(img_raw[i+m+1*2,j+n]-img_raw[i+m+3*2,j+n])
+    return extrapolate_quad
+
+def extrapolate_green_left_quad(img_raw,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] = img_raw[i+m,j+n-1*2] + 3/4*(img_raw[i+m,j+n]-img_raw[i+m,j+n-2*2])-1/4*(img_raw[i+m,j+n-1*2]-img_raw[i+m,j+n-3*2])
+    return extrapolate_quad
+
+def extrapolate_green_right_quad(img_raw,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] = img_raw[i+m,j+n+1*2] + 3/4*(img_raw[i+m,j+n]-img_raw[i+m,j+n+2*2])-1/4*(img_raw[i+m,j+n+1*2]-img_raw[i+m,j+n+3*2])
+    return extrapolate_quad
+
+### Extrapolation method:
+def median_extrapolate_green_pixel_quad(img_raw,i,j,orientations_to_drop):
+    list_extrapolate_pixel = []
+
+    if ("top" not in orientations_to_drop):
+        list_extrapolate_pixel.append(extrapolate_green_top_quad(img_raw,i,j))
+    if ("bottom" not in orientations_to_drop):
+        list_extrapolate_pixel.append(extrapolate_green_bottom_quad(img_raw,i,j))  
+    if("left" not in orientations_to_drop):
+        list_extrapolate_pixel.append(extrapolate_green_left_quad(img_raw,i,j))
+    if("right" not in orientations_to_drop):
+        list_extrapolate_pixel.append(extrapolate_green_right_quad(img_raw,i,j))
+    median_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            median_quad[m,n] = np.median([list_extrapolate_pixel[k] for k in range(len(list_extrapolate_pixel))])
+    return median_quad
+
+def extrapolate_green_pixel_quad(img_raw,i,j,orientation):
+    # First the borders:
+    orientations_to_drop = []
+    if (i<2):
+        orientations_to_drop.append('top')
+    if (i>img_raw.shape[0]-4*2):
+        orientations_to_drop.append('bottom')
+    if (j<2):
+        orientations_to_drop.append('left')
+    if (j>img_raw.shape[1]-4*2):
+        orientations_to_drop.append('right')
+    
+    
+    # Then the rest of the image:
+    else:
+        if (orientation >0.5): # V < H so we gonna eliminate one horizontal pixel. 
+            if ("right" not in orientations_to_drop and "left" not in orientations_to_drop):
+                rmse_pixel_left = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_green_left_quad(img_raw,i,j))
+                rmse_pixel_right = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_green_right_quad(img_raw,i,j))
+                if (np.sum(rmse_pixel_left) > np.sum(rmse_pixel_right)):
+                    orientations_to_drop.append('left')
+                else:
+                    orientations_to_drop.append('right')
+    
+        else: # V > H so we gonna eliminate one vertical pixel.
+            if ("top" not in orientations_to_drop and "bottom" not in orientations_to_drop):
+                rmse_pixel_top = rmse_pixel(img_raw[i+2,j+2],extrapolate_green_top_quad(img_raw,i,j))
+                rmse_pixel_bottom = rmse_pixel(img_raw[i+2,j+2],extrapolate_green_bottom_quad(img_raw,i,j))
+                if (np.sum(rmse_pixel_top) > np.sum(rmse_pixel_bottom)):
+                    orientations_to_drop.append('top')
+                else:
+                    orientations_to_drop.append('bottom')
+    return median_extrapolate_green_pixel_quad(img_raw,i,j,orientations_to_drop)
+
+def extrapolate_green_quad(img_raw,extrapolate_img):
+    orientation_matrix = compute_orientation_matrix(img_raw)
+    for i in range(0,img_raw.shape[0],2):
+        for j in range(0,img_raw.shape[1],2):
+            if (color_pixel(i,j,'quad_bayer')!= "green"):
+                extrapolate_img[i:i+2,j:j+2,1] = extrapolate_green_pixel_quad(img_raw,i,j,(1/4) *np.sum(orientation_matrix[i:i+2,j:j+2]))
+            else:
+                extrapolate_img[i:i+2,j:j+2,1] = img_raw[i:i+2,j:j+2]
+    return extrapolate_img
+
+## Red and Blue Channels ##
+
+def extrapolate_top_quad(img_raw,img_extrapolate,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] = (img_raw[i+m-1*2,j+n] + img_raw[i+m,j+n]-img_extrapolate[i+m-1*2,j+n,1])
+    return extrapolate_quad
+def extrapolate_left_quad(img_raw,img_extrapolate,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] = (img_raw[i+m,j+n-1*2] + img_raw[i+m,j+n]-img_extrapolate[i+m,j+n-1*2,1])
+    return extrapolate_quad
+def extrapolate_right_quad(img_raw,img_extrapolate,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] = (img_raw[i+m,j+n+1*2] + img_raw[i+m,j+n]-img_extrapolate[i+m,j+n+1*2,1])
+    return extrapolate_quad
+def extrapolate_bottom_quad(img_raw,img_extrapolate,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] = (img_raw[i+m+1*2,j+n] + img_raw[i+m,j+n]-img_extrapolate[i+m+1*2,j+n,1])
+    return extrapolate_quad
+
+def extrapolate_top_left_quad(img_raw,img_extrapolate,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] =(img_raw[i+m-1*2,j+n-1*2] + img_extrapolate[i+m,j+n,1]-img_extrapolate[i+m-1*2,j+n-1*2,1])
+    return extrapolate_quad
+def extrapolate_top_right_quad(img_raw,img_extrapolate,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] =(img_raw[i+m-1*2,j+n+1*2] + img_extrapolate[i+m,j+n,1]-img_extrapolate[i+m-1*2,j+n+1*2,1])
+    return extrapolate_quad
+def extrapolate_bottom_left_quad(img_raw,img_extrapolate,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] =(img_raw[i+m+1*2,j+n-1*2] + img_extrapolate[i+m,j+n,1]-img_extrapolate[i+m+1*2,j+n-1*2,1])
+    return extrapolate_quad
+def extrapolate_bottom_right_quad(img_raw,img_extrapolate,i,j):
+    extrapolate_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            extrapolate_quad[m,n] =(img_raw[i+m+1*2,j+n+1*2] + img_extrapolate[i+m,j+n,1]-img_extrapolate[i+m+1*2,j+n+1*2,1])
+    return extrapolate_quad
+
+def median_pixel_quad(img_raw,img_extrapolate,i,j,orientations_to_drop):
+    list_extrapolate = []
+    if (color_pixel(i,j,"quad_bayer") != "green"):
+        if("top_left" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_top_left_quad(img_raw,img_extrapolate,i,j))
+        if("top_right" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_top_right_quad(img_raw,img_extrapolate,i,j))
+        if("bottom_left" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_bottom_left_quad(img_raw,img_extrapolate,i,j))
+        if("bottom_right" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_bottom_right_quad(img_raw,img_extrapolate,i,j))
+
+    elif (color_pixel(i,j,"quad_bayer") == "green"):
+        if("top" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_top_quad(img_raw,img_extrapolate,i,j))
+        if("left" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_left_quad(img_raw,img_extrapolate,i,j))
+        if("right" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_right_quad(img_raw,img_extrapolate,i,j))
+        if("bottom" not in orientations_to_drop):
+            list_extrapolate.append(extrapolate_bottom_quad(img_raw,img_extrapolate,i,j))
+    
+    median_quad = np.zeros((2,2))
+    for m in range(2):
+        for n in range(2):
+            median_quad[m,n] = np.median([list_extrapolate[k][m,n] for k in range(len(list_extrapolate))])
+    return median_quad
+
+def extrapolate_pixel_quad(img_raw,img_extrapolate,i,j,color):
+    orientations_to_drop = []
+
+    if (color_pixel(i,j,"quad_bayer")!='green'):
+        if (i<1):
+            orientations_to_drop.append("top_left")
+            orientations_to_drop.append("top_right")
+        if (i>img_raw.shape[0]-2*2):
+            orientations_to_drop.append("bottom_left")
+            orientations_to_drop.append("bottom_right")
+        if (j<1):
+            orientations_to_drop.append("top_left")
+            orientations_to_drop.append("bottom_left")
+        if (j>img_raw.shape[1]-2*2):
+            orientations_to_drop.append("top_right")
+            orientations_to_drop.append("bottom_right")
+        
+        
+        if ("top_left" not in orientations_to_drop and "top_right" not in orientations_to_drop and "bottom_left" not in orientations_to_drop and "bottom_right" not in orientations_to_drop):
+            rmse_top_left = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_top_left_quad(img_raw,img_extrapolate,i,j))
+            rmse_top_right = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_top_right_quad(img_raw,img_extrapolate,i,j))
+            rmse_bottom_left = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_bottom_left_quad(img_raw,img_extrapolate,i,j))
+            rmse_bottom_right = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_bottom_right_quad(img_raw,img_extrapolate,i,j))
+            if (rmse_bottom_left> rmse_bottom_right and rmse_bottom_left> rmse_top_left and rmse_bottom_left> rmse_top_right):
+                orientations_to_drop.append("bottom_left")
+            elif (rmse_bottom_right> rmse_bottom_left and rmse_bottom_right> rmse_top_left and rmse_bottom_right> rmse_top_right):
+                orientations_to_drop.append("bottom_right")
+            elif (rmse_top_left> rmse_bottom_left and rmse_top_left> rmse_bottom_right and rmse_top_left> rmse_top_right):
+                orientations_to_drop.append("top_left")
+            else:
+                orientations_to_drop.append("top_right")
+    elif(color_pixel(i,j,"quad_bayer")=="green"):
+        if (i<1):
+            orientations_to_drop.append("top")
+        if (i>img_raw.shape[0]-2*2):
+            orientations_to_drop.append("bottom")
+        if (j<1):
+            orientations_to_drop.append("left")
+        if (j>img_raw.shape[1]-2*2):
+            orientations_to_drop.append("right")
+        
+        if (((i/2)%2!=0 and color == "red") or ((i/2)%2==0 and color == "blue")):
+            if ("right" not in orientations_to_drop and "left" not in orientations_to_drop):
+                rmse_pixel_left = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_left_quad(img_raw,img_extrapolate,i,j))
+                rmse_pixel_right = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_right_quad(img_raw,img_extrapolate,i,j))
+                if (rmse_pixel_left > rmse_pixel_right):
+                    orientations_to_drop.append('left')
+                else:
+                    orientations_to_drop.append('right')
+        else:
+            if ("top" not in orientations_to_drop and "bottom" not in orientations_to_drop):
+                rmse_pixel_top = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_top_quad(img_raw,img_extrapolate,i,j))
+                rmse_pixel_bottom = rmse_pixel(img_raw[i:i+2,j:j+2],extrapolate_bottom_quad(img_raw,img_extrapolate,i,j))
+                if (rmse_pixel_top > rmse_pixel_bottom):
+                    orientations_to_drop.append('top')
+                else:
+                    orientations_to_drop.append('bottom')
+    return median_pixel_quad(img_raw,img_extrapolate,i,j,orientations_to_drop)
+
+def extrapolate_red_quad(img_raw,img_extrapolate):
+    for i in range(0,img_raw.shape[0],2):
+        for j in range(0,img_raw.shape[1],2):
+            if (color_pixel(i,j,"quad_bayer")!="red"):
+                img_extrapolate[i:i+2,j:j+2,0] = extrapolate_pixel_quad(img_raw,img_extrapolate,i,j,"red")
+            else:
+                img_extrapolate[i:i+2,j:j+2,0] = img_raw[i:i+2,j:j+2]
+
+def extrapolate_blue_quad(img_raw,img_extrapolate):
+    for i in range(0,img_raw.shape[0],2):
+        for j in range(0,img_raw.shape[1],2):
+            if (color_pixel(i,j,"quad_bayer")!="blue"):
+                img_extrapolate[i:i+2,j:j+2,2] = extrapolate_pixel_quad(img_raw,img_extrapolate,i,j,"blue")
+            else:
+                img_extrapolate[i:i+2,j:j+2,2] = img_raw[i:i+2,j:j+2]
+
+def extrapolate_img_quad(img_cfa):
+    extapolate_img = np.zeros(img_cfa.shape + (3,))
+    extrapolate_green_quad(img_cfa,extapolate_img)
+    extrapolate_red_quad(img_cfa,extapolate_img)
+    extrapolate_blue_quad(img_cfa,extapolate_img)
+    return extapolate_img
+
+def extrapolate_cfa(img_cfa,cfa):
+    if (cfa=="bayer"):
+        return extrapolate_img(img_cfa)
+    elif(cfa=="quad_bayer"):
+        return extrapolate_img_quad(img_cfa)
+    else:
+        print("Error: cfa not recognized")
+        return None
\ No newline at end of file

src/methods/CASA_NOVA_Vincent/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 functions as fu
+from src.forward_model import CFA
+import importlib 
+importlib.reload(fu)
+
+
+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
+
+    return fu.extrapolate_cfa(y,cfa)
+

src/methods/Chardon_tom/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.Chardon_tom.utils import *
+import pywt
+
+#!!!!!!!! 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)
+
+    res = op.adjoint(y)
+
+    N,M = input_shape[0], input_shape[1]
+
+    
+
+    #interpolating green channel
+
+    for i in range (N):
+        for j in range (M):
+            if res[i,j,1] ==0:
+
+                neighbors = get_neighbors(res,1,i,j,N,M)
+                weights = get_weights(res,i,j,1,N,M)
+                res[i,j,1] = interpolate_green(weights, neighbors)
+
+
+
+    #first intepolation of red channel
+
+    for i in range (1,N,2):
+        for j in range (0,M,2):
+            neighbors = get_neighbors(res,0,i,j,N,M)
+            neighbors_G = get_neighbors(res,1,i,j,N,M)
+            weights = get_weights(res,i,j,0,N,M) 
+            res[i,j,0] = interpolate_red_blue(weights,neighbors, neighbors_G)
+
+    # second interpolation of red channel
+
+    for i in range (N):
+        for j in range (M):
+            if res[i,j,0] ==0:
+                neighbors = get_neighbors(res,0,i,j,N,M)
+                weights = get_weights(res,i,j,0,N,M)
+                res[i,j,0] = interpolate_green(weights, neighbors)
+
+
+    #first interpolation of blue channel
+
+    for i in range (0,N,2):
+        for j in range (1,M,2):
+            neighbors = get_neighbors(res,2,i,j,N,M)
+            neighbors_G = get_neighbors(res,1,i,j,N,M)
+            weights = get_weights(res,i,j,2,N,M) 
+            res[i,j,2] = interpolate_red_blue(weights, neighbors, neighbors_G)
+
+    #second interpolation of blue channel
+
+    for i in range (N):
+        for j in range (M):
+            if res[i,j,2] ==0:
+                neighbors = get_neighbors(res,2,i,j,N,M)
+                weights = get_weights(res,i,j,2,N,M)
+                res[i,j,2] = interpolate_green(weights,neighbors)
+
+
+
+    # k=0
+    # while k<2 :
+    #     for i in range(input_shape[0]):
+    #         for j in range(input_shape[1]):
+    #             res[i][j][1] = correction_green(res,i,j,N,M) 
+    #     for i in range(input_shape[0]):
+    #         for j in range(input_shape[1]):
+    #             res[i][j][0] = correction_red(res,i,j,N,M) 
+    #     for i in range(input_shape[0]):
+    #         for j in range(input_shape[1]):
+    #             res[i][j][2] = correction_blue(res,i,j,N,M) 
+    #     k+=1
+
+    res[res>1] = 1
+    res[res<0] = 0
+
+
+    return res
+

src/methods/Chardon_tom/utils.py

+import numpy as np
+import pywt
+
+def get_neighbors (img,channel,i,j,N,M):
+
+    P1 = img[(i-1)%N,(j-1)%M,channel]
+    P2 = img[(i-1)%N,j%M,channel]
+    P3 = img[(i-1)%N,(j+1)%M,channel]
+    P4 = img[i%N,(j-1)%M,channel]
+    P5 = img[i%N,j%M,channel]
+    P6 = img[i%N,(j+1)%M,channel]
+    P7 = img[(i+1)%N,(j-1)%M,channel]
+    P8 = img[(i+1)%N,j%M,channel]
+    P9 = img[(i+1)%N,(j+1)%M,channel]
+
+    return np.array([P1,P2,P3,P4,P5,P6,P7,P8,P9])
+
+
+def get_derivatives(neighbors):
+
+    [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbors
+
+    D_x = (P4 - P6)/2
+    D_y = (P2 - P8)/2
+    D_xd = (P3 - P7)/(2*np.sqrt(2))
+    D_yd = (P1 - P9)/(2*np.sqrt(2))
+
+    return ([D_x, D_y, D_xd, D_yd])
+
+
+def get_weights(mosaic_image, i, j, channel, N, M):
+
+    derivatives_neigbors = []
+    for l in range(-1, 2):
+        for L in range(-1, 2):
+            derivatives_neigbors.append(get_derivatives(
+                get_neighbors(mosaic_image, channel, i+l, j+L, N, M)))
+
+    [Dx, Dy, Dxd, Dyd] = derivatives_neigbors[4]
+    E1 = 1/np.sqrt(1 + Dyd**2 + derivatives_neigbors[0][3]**2)
+    E2 = 1/np.sqrt(1 + Dy**2 + derivatives_neigbors[1][1]**2)
+    E3 = 1/np.sqrt(1 + Dxd**2 + derivatives_neigbors[2][2]**2)
+    E4 = 1/np.sqrt(1 + Dx**2 + derivatives_neigbors[3][0]**2)
+    E6 = 1/np.sqrt(1 + Dxd**2 + derivatives_neigbors[5][2]**2)
+    E7 = 1/np.sqrt(1 + Dy**2 + derivatives_neigbors[6][1]**2)
+    E8 = 1/np.sqrt(1 + Dyd**2 + derivatives_neigbors[7][3]**2)
+    E9 = 1/np.sqrt(1 + Dx**2 + derivatives_neigbors[8][0]**2)
+    E = [E1, E2, E3, E4, E6, E7, E8, E9]
+
+    return E    
+
+
+def interpolate_green(weights, neighbors):
+
+    [E1, E2, E3, E4, E6, E7, E8, E9] = weights
+    [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbors
+
+    I5 = (E2*P2 + E4*P4 + E6*P6 + E8*P8)/(E2 + E4 + E6 + E8)
+
+    return (I5)
+
+
+def interpolate_red_blue(weights, neighbors, green_neighbors):
+
+    [E1, E2, E3, E4, E6, E7, E8, E9] = weights
+    [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbors
+    [G1, G2, G3, G4, G5, G6, G7, G8, G9] = green_neighbors
+
+    I5 = G5*(E1*P1/G1 + E3*P3/G3 + E7*P7/G7 + E9*P9/G9)/(E1 + E3 + E7 + E9)
+
+    return (I5)
+
+
+def correction_green(res,i,j,N,M):
+
+
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = get_neighbors(res,1,i,j,N,M)
+    [R1,R2,R3,R4,R5,R6,R7,R8,R9] = get_neighbors(res,0,i,j,N,M)
+    [B1,B2,B3,B4,B5,B6,B7,B8,B9] = get_neighbors(res,2,i,j,N,M)
+    [E1,E2,E3,E4,E6,E7,E8,E9] = get_weights(res,i,j,1,N,M)
+
+    Gb5 = R5*((E2*G2)/B2 + (E4*G4)/B4 + (E6*G6)/B6 + (E8*G8)/B8)/(E2 + E4 + E6 + E8)
+    Gr5 = B5*((E2*G2)/R2 + (E4*G4)/R4 + (E6*G6)/R6 + (E8*G8)/R8)/(E2 + E4 + E6 + E8)
+
+    G5 = (Gb5 + Gr5)/2
+
+    return G5
+
+def correction_red(res,i,j,N,M) :
+
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = get_neighbors(res,1,i,j,N,M)
+    [R1,R2,R3,R4,R5,R6,R7,R8,R9] = get_neighbors(res,0,i,j,N,M)
+    [E1,E2,E3,E4,E6,E7,E8,E9] = get_weights(res,i,j,0,N,M)
+
+    R5 = G5*((E1*R1)/G1 + (E2*R2)/G2 + (E3*R3)/G3 + (E4*R4)/G4 + (E6*R6)/G6 + (E7*R7)/G7 + (E8*R8)/G8 + (E9*R9)/G9)/(E1 + E2 + E3 + E4 + E6 + E7 + E8 + E9)
+
+    return R5
+
+def correction_blue(res,i,j,N,M) :
+
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = get_neighbors(res,1,i,j,N,M)
+    [B1,B2,B3,B4,B5,B6,B7,B8,B9] = get_neighbors(res,2,i,j,N,M)
+    [E1,E2,E3,E4,E6,E7,E8,E9] = get_weights(res,i,j,2,N,M)
+
+    B5 = G5*((E1*B1)/G1 + (E2*B2)/G2 + (E3*B3)/G3 + (E4*B4)/G4 + (E6*B6)/G6 + (E7*B7)/G7 + (E8*B8)/G8 + (E9*B9)/G9)/(E1 + E2 + E3 + E4 + E6 + E7 + E8 + E9)
+
+    return B5
+
+
+
+

src/methods/ELAMRANI_Mouna/functions.py

+import numpy as np
+
+def find_Knearest_neighbors(z, chan, i, j, N, M):
+    """Finds a pixel's neighbors on a channel"""
+    return np.array([z[(i+di)%N, (j+dj)%M, chan] for di in range(-1, 2) for dj in range(-1, 2)])
+
+def calculate_directional_gradients(neighbors):
+    """Gives the directional derivative"""
+    P1, P2, P3, P4, P5, P6, P7, P8, P9 = neighbors
+    Dx, Dy = (P4 - P6)/2, (P2 - P8)/2
+    Dxd, Dyd = (P3 - P7)/(2*np.sqrt(2)), (P1 - P9)/(2*np.sqrt(2))
+    return [Dx, Dy, Dxd, Dyd]
+
+def calculate_adaptive_weights(z, neigh, dir_deriv,chan,i,j,N,M):
+
+    [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_Knearest_neighbors(z,chan,i+k,j+k,N,M)
+            dd = calculate_directional_gradients(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_pixel(neigh,weights):
+    """This function performs interpolation for a single pixel by calculating a weighted average of its neighboring pixels"""
+    [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_RedBlue(neighbors, neighbors_G, weights):
+    """This function specifically interpolates a pixel in the red or blue channels"""
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neighbors
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neighbors_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

src/methods/ELAMRANI_Mouna/reconstruct.py

+import numpy as np
+from src.forward_model import CFA
+from src.methods.ELAMRANI_Mouna.functions import *
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+
+    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]
+
+    def interpolate_channel(img_res, channel, first_pass, N, M):
+        for i in range(N):
+            for j in range(M):
+                if first_pass and ((channel == 0 and i % 2 == 1 and j % 2 == 0) or
+                               (channel == 2 and i % 2 == 0 and j % 2 == 1)):
+                    neighbors = find_Knearest_neighbors(img_res, channel, i, j, N, M)
+                    neighbors_G = find_Knearest_neighbors(img_res, 1, i, j, N, M)
+                    dir_deriv = calculate_directional_gradients(neighbors_G)
+                    weights = calculate_adaptive_weights(img_res, neighbors_G, dir_deriv, 1, i, j, N, M)
+                    img_res[i, j, channel] = interpolate_RedBlue(neighbors, neighbors_G, weights)
+                elif not first_pass and img_res[i, j, channel] == 0:
+                    neighbors = find_Knearest_neighbors(img_res, channel, i, j, N, M)
+                    dir_deriv = calculate_directional_gradients(neighbors)
+                    weights = calculate_adaptive_weights(img_res, neighbors, dir_deriv, channel, i, j, N, M)
+                    img_res[i, j, channel] = interpolate_pixel(neighbors, weights)
+        return img_res
+
+# Interpolation pour chaque canal
+    img_res = interpolate_channel(img_res, 1, False, N, M)  # Interpolation du canal vert
+    img_res = interpolate_channel(img_res, 0, True, N, M)   # Première interpolation du canal rouge
+    img_res = interpolate_channel(img_res, 0, False, N, M)  # Seconde interpolation du canal rouge
+    img_res = interpolate_channel(img_res, 2, True, N, M)   # Première interpolation du canal bleu
+    img_res = interpolate_channel(img_res, 2, False, N, M)  # Seconde interpolation du canal bleu
+
+    img_res[img_res > 1] = 1
+    img_res[img_res < 0] = 0
+
+    return img_res
+

src/methods/EL_MURR_Theresa/malvar.py

+import numpy as np
+from scipy.signal import correlate2d
+from src.forward_model import CFA
+
+def malvar_he_cutler(y: np.ndarray, op: CFA ) -> np.ndarray:
+    """Performs demosaicing using the malvar-he-cutler algorithm
+
+    Args:
+        op (CFA): CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+
+    red_mask, green_mask, blue_mask =  [op.mask[:, :, 0], op.mask[:, :, 1], op.mask[:, :, 2]]
+    mosaicked_image = np.float32(y)
+    demosaicked_image = np.empty(op.input_shape)
+
+    if op.cfa == 'quad_bayer':
+        filters = get_quad_bayer_filters()
+    else:
+        filters = get_default_filters()
+
+    demosaicked_image = apply_demosaicking_filters(
+        mosaicked_image,demosaicked_image, red_mask, green_mask, blue_mask, filters
+    )
+
+    return demosaicked_image
+
+def get_quad_bayer_filters():
+    coefficient_scale = 0.03125
+    return {
+        "G_at_R_and_B": np.array([
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, 2, 2, 0, 0, 0, 0],
+            [0, 0, 0, 0, 2, 2, 0, 0, 0, 0],
+            [-1, -1, 2, 2, 4, 4, 2, 2, -1, -1],
+            [-1, -1, 2, 2, 4, 4, 2, 2, -1, -1],
+            [0, 0, 0, 0, 2, 2, 0, 0, 0, 0],
+            [0, 0, 0, 0, 2, 2, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0]
+        ]) * coefficient_scale,
+        "R_at_GR_and_B_at_GB": np.array([
+            [0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0],
+            [0, 0, -1, -1, 0, 0, -1, -1, 0, 0],
+            [0, 0, -1, -1, 0, 0, -1, -1, 0, 0],
+            [-1, -1, 4, 4, 5, 5, 4, 4, -1, -1],
+            [-1, -1, 4, 4, 5, 5, 4, 4, -1, -1],
+            [0, 0, -1, -1, 0, 0, -1, -1, 0, 0],
+            [0, 0, -1, -1, 0, 0, -1, -1, 0, 0],
+            [0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0]
+        ]) * coefficient_scale,
+        "R_at_GB_and_B_at_GR": np.array([
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, -1, -1, 4, 4, -1, -1, 0, 0],
+            [0, 0, -1, -1, 4, 4, -1, -1, 0, 0],
+            [0.5, 0.5, 0, 0, 5, 5, 0, 0, 0.5, 0.5],
+            [0.5, 0.5, 0, 0, 5, 5, 0, 0, 0.5, 0.5],
+            [0, 0, -1, -1, 4, 4, -1, -1, 0, 0],
+            [0, 0, -1, -1, 4, 4, -1, -1, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0]
+        ]) * coefficient_scale,
+        "R_at_B_and_B_at_R": np.array([
+            [0, 0, 0, 0, -1.5, -1.5, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1.5, -1.5, 0, 0, 0, 0],
+            [0, 0, 2, 2, 0, 0, 2, 2, 0, 0],
+            [0, 0, 2, 2, 0, 0, 2, 2, 0, 0],
+            [-1.5, -1.5, 0, 0, 6, 6, 0, 0, -1.5, -1.5],
+            [-1.5, -1.5, 0, 0, 6, 6, 0, 0, -1.5, -1.5],
+            [0, 0, 2, 2, 0, 0, 2, 2, 0, 0],
+            [0, 0, 2, 2, 0, 0, 2, 2, 0, 0],
+            [0, 0, 0, 0, -1.5, -1.5, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1.5, -1.5, 0, 0, 0, 0]
+        ]) * coefficient_scale,
+    }
+
+def get_default_filters():
+    coefficient_scale = 0.125
+    return {
+        "G_at_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]
+        ]) * coefficient_scale,
+        "R_at_GR_and_B_at_GB": np.array([
+            [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]
+        ]) * coefficient_scale,
+        "R_at_GB_and_B_at_GR": np.array([
+            [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]
+        ]) * coefficient_scale,
+        "R_at_B_and_B_at_R": np.array([
+            [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]
+        ]) * coefficient_scale,
+    }
+
+def apply_demosaicking_filters(image, res, red_mask, green_mask, blue_mask, filters):
+    red_channel = image * red_mask
+    green_channel = image * green_mask
+    blue_channel = image * blue_mask
+
+    # Create the green channel after applying a filter
+    green_channel = np.where(
+        np.logical_or(red_mask == 1, blue_mask == 1),
+        correlate2d(image, filters['G_at_R_and_B'], mode="same", boundary="symm"),
+        green_channel
+    )
+
+
+    # Define masks for extracting pixel values
+    red_row_mask = np.any(red_mask == 1, axis=1)[:, np.newaxis].astype(np.float32)
+    red_col_mask = np.any(red_mask == 1, axis=0)[np.newaxis].astype(np.float32)
+
+    blue_row_mask = np.any(blue_mask == 1, axis=1)[:, np.newaxis].astype(np.float32)
+    blue_col_mask = np.any(blue_mask == 1, axis=0)[np.newaxis].astype(np.float32)
+
+    def update_channel(channel, row_mask, col_mask, filter_key):
+        return np.where(
+            np.logical_and(row_mask == 1, col_mask == 1),
+            correlate2d(image, filters[filter_key], mode="same", boundary="symm"),
+            channel
+        )
+
+# Update the red channel and blue channel
+    red_channel = update_channel(red_channel, red_row_mask, blue_col_mask, 'R_at_GR_and_B_at_GB')
+    red_channel = update_channel(red_channel, blue_row_mask, red_col_mask, 'R_at_GB_and_B_at_GR')
+
+    blue_channel = update_channel(blue_channel, blue_row_mask, red_col_mask, 'R_at_GR_and_B_at_GB')
+    blue_channel = update_channel(blue_channel, red_row_mask, blue_col_mask, 'R_at_GB_and_B_at_GR')
+
+    # Update R channel and B channel again
+    red_channel = update_channel(red_channel, blue_row_mask, blue_col_mask, 'R_at_B_and_B_at_R')
+    blue_channel = update_channel(blue_channel, red_row_mask, red_col_mask, 'R_at_B_and_B_at_R')
+    res[:, :, 0] = red_channel
+    res[:, :, 1] = green_channel
+    res[:, :, 2] = blue_channel
+    return res
\ No newline at end of file

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

src/methods/Elmehdi_lahmar/malvar.py

+import numpy as np
+from scipy.signal import convolve2d
+from src.forward_model import CFA
+
+def malvar(y: np.ndarray, op: CFA) -> np.ndarray:
+    """
+    Malvar-He-Cutler demosaicing algorithm for Bayer pattern.
+    """
+    # Convert the mosaicked image to the initial estimated channels.
+    z = op.adjoint(y)
+
+    # Define convolution kernels
+    kernel_G_at_RB = 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
+    kernel_RB_at_G = np.array([[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]]) / 8
+    kernel_RB_at_RB = np.array([[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]]) / 8
+
+    # Interpolate each channel
+    R = z[:, :, 0]
+    G = z[:, :, 1]
+    B = z[:, :, 2]
+
+    mask = op.mask
+    R_m, G_m, B_m = mask[:, :, 0], mask[:, :, 1], mask[:, :, 2]
+
+    # Interpolate G at R and B locations
+    G = np.where(np.logical_or(R_m == 1, B_m == 1), convolve2d(y, kernel_G_at_RB, mode='same'), G)
+
+    # Interpolate R at G and B locations, B at R and G locations
+    R = np.where(np.logical_or(G_m == 1, B_m == 1), convolve2d(y, kernel_RB_at_G, mode='same'), R)
+    B = np.where(np.logical_or(R_m == 1, G_m == 1), convolve2d(y, kernel_RB_at_G, mode='same'), B)
+
+    # Interpolate R at B locations and B at R locations
+    R = np.where(B_m == 1, convolve2d(y, kernel_RB_at_RB, mode='same'), R)
+    B = np.where(R_m == 1, convolve2d(y, kernel_RB_at_RB, mode='same'), B)
+
+    # Combine channels
+    return np.clip(np.stack((R, G, B), axis=-1), 0, 1)
+
+
+
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller
\ No newline at end of file