src/methods/leo/dataset_creation.py

+import csv
+import requests
+import os
+
+# Chemin du fichier CSV
+csv_file_path = '..\\published_images.csv'
+
+# Dossier de destination pour enregistrer les images
+destination_folder = 'dataset/'
+
+# Nombre de lignes à télécharger
+nombre_lignes_a_telecharger = 500
+
+# Vérifiez si le dossier de destination existe, sinon créez-le
+if not os.path.exists(destination_folder):
+    os.makedirs(destination_folder)
+
+# Ouvrir le fichier CSV
+with open(csv_file_path, newline='', encoding='utf-8') as csvfile:
+    # Lire le fichier CSV
+    csv_reader = csv.reader(csvfile)
+    
+    # Ignorer l'en-tête si nécessaire
+    next(csv_reader, None)
+    
+    # Parcourir les 1000 premières lignes du CSV
+    for i, row in enumerate(csv_reader):
+        if i >= nombre_lignes_a_telecharger:
+            break  # Sortir de la boucle une fois que le nombre de lignes souhaité est atteint
+
+        print(f"Téléchargement de l'image {i + 1}...")
+
+        # Récupérer l'URL de l'image depuis la deuxième colonne (index 1)
+        image_url = row[2]
+        print("img url : ",image_url)
+        # Récupérer le nom de fichier depuis la première colonne (index 0)
+        image_filename = row[0] + '.jpg'
+        
+        # Chemin complet du fichier destination
+        destination_path = os.path.join(destination_folder, image_filename)
+
+        try:
+            # Télécharger l'image et l'enregistrer localement
+            response = requests.get(image_url)
+            response.raise_for_status()  # Vérifier si la requête a réussi
+
+            with open(destination_path, 'wb') as img_file:
+                img_file.write(response.content)
+                #print(f"Chemin complet du fichier destination : {destination_path}")
+
+            print(f"Image {i + 1} téléchargée avec succès.")
+        except Exception as e:
+            print(f"Erreur lors du téléchargement de l'image {i + 1}: {e}")
\ No newline at end of file

src/methods/leo/dataset_format.py

+import cv2
+import os
+import numpy as np
+from src.forward_model import CFA
+from src.utils import load_image, save_image, psnr, ssim
+image_path = 'images/img_1.png'
+
+img = load_image(image_path)
+op = CFA('bayer', img.shape)
+dossier_images = "dataset"
+
+
+images_redimensionnees = []
+
+
+for fichier in os.listdir(dossier_images):
+
+    if fichier.endswith(".jpg"):
+    
+        chemin_image = os.path.join(dossier_images, fichier)
+        image = cv2.imread(chemin_image)
+        
+        
+        image_redimensionnee = cv2.resize(image, (1024, 1024))
+        
+       
+        images_redimensionnees.append(image_redimensionnee)
+
+
+tableau_images_or = np.stack(images_redimensionnees, axis=-1)
+
+
+print("Dimensions du tableau d'images :", tableau_images_or.shape)
+with open('array.npy', 'wb') as f:
+    np.save(f, tableau_images_or)
+
+images_redimensionnees = []
+
+
+for k in range(500):
+    
+   
+    image = tableau_images_or[:,:,:,k]
+        
+    
+    y = op.direct(image)
+    z = op.adjoint(y)
+    
+    images_redimensionnees.append(z)
+tableau_images = np.empty((1024, 1024, 3, len(images_redimensionnees)), dtype=np.uint8)
+
+# Remplissez le tableau avec les images de la liste
+for i, image in enumerate(images_redimensionnees):
+    print(i)
+    tableau_images[:, :, :, i] = image
+with open('array_z.npy', 'wb') as f:
+    np.save(f, tableau_images)
\ No newline at end of file

src/methods/leo/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
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    # Performing the reconstruction.
+    # TODO
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+    z = op.adjoint(y)
+    interpolated_image = nearest_neighbor_interpolation(z)
+    return interpolated_image
+
+
+import numpy as np
+
+
+
+
+
+
+def nearest_neighbor_interpolation(image):
+    interpolated_image = np.zeros_like(image)
+    for i in range(image.shape[0]):
+        for j in range(image.shape[1]):
+            for k in range(3):
+                if (image[i, j,k] == 0):  
+                    closest_pixel = find_closest_nonzero_pixel(image[:,:,k], i, j)
+                    interpolated_image[i, j,k] = closest_pixel
+                else:
+                    interpolated_image[i, j,k] = image[i, j,k]  
+    return interpolated_image
+
+
+def find_closest_nonzero_pixel(image, i, j):
+    m, n = image.shape
+    nb_pixel = 1
+    while True:
+        for x in range(max(0, i - nb_pixel), min(m, i + nb_pixel + 1)):
+            for y in range(max(0, j - nb_pixel), min(n, j + nb_pixel + 1)):
+                if np.any(image[x, y] != 0):
+                    return image[x, y]
+        nb_pixel += 1
+
+
+
+
+# Afficher ou enregistrer l'image interpolée
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/leo/reconstruct_2.py

+"""The main file for the baseline reconstruction.
+This file should NOT be modified.
+"""
+
+
+import numpy as np
+
+from src.forward_model import CFA
+from src.methods.baseline.demo_reconstruction import naive_interpolation
+from tensorflow.keras.models import load_model
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+    z = op.adjoint(y)
+    model = load_model("model")
+    predictions = model.predict(z)
+
+    return predictions
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller
+
+
+import tensorflow as tf
+from tensorflow.keras import layers, models
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+with open('array.npy', 'rb') as f:
+    tableau_image_or = np.load(f)
+with open('array_z.npy', 'rb') as f:
+    tableau_images = np.load(f)
+    
+X = tableau_images.transpose((3,0,1,2))
+y = tableau_image_or.transpose((3,0,1,2))
+
+# Création du modèle
+model = models.Sequential()
+
+# Encoder
+model.add(layers.Conv2D(32, (3, 3), activation='relu', padding='same', input_shape=(1024, 1024, 3)))
+model.add(layers.MaxPooling2D((2, 2), padding='same'))
+
+# Decoder
+model.add(layers.Conv2D(32, (3, 3), activation='relu', padding='same'))
+model.add(layers.UpSampling2D((2, 2)))
+model.add(layers.Conv2D(3, (3, 3), activation='sigmoid', padding='same'))
+
+# Compiler le modèle
+model.compile(optimizer='adam', loss='mse')  # Mean Squared Error comme fonction de perte
+
+# Entraînement du modèle
+model.fit(X, y, epochs=3, batch_size=10, validation_split=0.2)
+
+# Évaluation du modèle
+loss = model.evaluate(X, y)
+print(f"Loss sur l'ensemble de données : {loss}")
+
+# Faire des prédictions
+predictions = model.predict(X)
+
+# Afficher quelques exemples
+for i in range(5):  # Afficher les 5 premières images par exemple
+    plt.figure(figsize=(10, 5))
+
+    # Image originale
+    plt.subplot(1, 3, 1)
+    plt.imshow(y[i])
+    plt.title('Image Originale')
+
+    # Image avec pixels masqués
+    plt.subplot(1, 3, 2)
+    plt.imshow(X[i])
+    plt.title('Image avec Pixels Masqués')
+
+    # Image prédite
+    plt.subplot(1, 3, 3)
+    plt.imshow(predictions[i])
+    plt.title('Image Prédite')
+
+    plt.show()
\ No newline at end of file

src/methods/leo/train_model.py

+import tensorflow as tf
+from tensorflow.keras import layers, models
+import numpy as np
+import matplotlib.pyplot as plt
+
+with open('array.npy', 'rb') as f:
+    tableau_image_or = np.load(f)
+    
+with open('array_z.npy', 'rb') as f:
+    tableau_images_z = np.load(f)
+
+X = tableau_images_z.transpose((3,0,1,2))
+y = tableau_image_or.transpose((3,0,1,2))
+
+# Création du modèle
+model = models.Sequential()
+
+# Encoder
+model.add(layers.Conv2D(32, (3, 3), activation='relu', padding='same', input_shape=(1024, 1024, 3)))
+model.add(layers.MaxPooling2D((2, 2), padding='same'))
+
+# Decoder
+model.add(layers.Conv2D(32, (3, 3), activation='relu', padding='same'))
+model.add(layers.UpSampling2D((2, 2)))
+model.add(layers.Conv2D(3, (3, 3), activation='sigmoid', padding='same'))
+
+# Compiler le modèle
+model.compile(optimizer='adam', loss='mse')  # Mean Squared Error comme fonction de perte
+
+# Entraînement du modèle
+model.fit(X, y, epochs=3, batch_size=10, validation_split=0.2)
+
+# Évaluation du modèle
+loss = model.evaluate(X, y)
+print(f"Loss sur l'ensemble de données : {loss}")
+
+# Faire des prédictions
+model.save("model")
+
+

src/methods/lioretn/README.txt

+Quite simple to use : reconstruct.py file calls demosaicking.py (where is made the demosaicking) and is called by the main.ipynb
\ No newline at end of file

src/methods/lioretn/demosaicking.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
+
+# High quality linear interpolation filters for bayer mosaic
+bayer_g_at_r = np.array([[ 0,  0, -1,  0,  0],
+                         [ 0,  0,  2,  0,  0],
+                         [-1,  2,  4,  2, -1],
+                         [ 0,  0,  2,  0,  0],
+                         [ 0,  0, -1,  0,  0]]) / 8
+bayer_g_at_b = bayer_g_at_r
+bayer_r_at_green_rrow_bcol = np.array([[ 0,  0, 1/2,  0,  0],
+                                       [ 0, -1,  0,  -1,  0],
+                                       [-1,  4,  5,   4, -1],
+                                       [ 0, -1,  0,  -1,  0],
+                                       [ 0,  0, 1/2,  0,  0]]) / 8
+bayer_r_at_green_brow_rcol = bayer_r_at_green_rrow_bcol.T
+bayer_b_at_green_rrow_bcol = bayer_r_at_green_brow_rcol
+bayer_b_at_green_brow_rcol = bayer_r_at_green_rrow_bcol
+bayer_r_at_b = np.array([[  0,  0, -3/2, 0,   0],
+                            [  0,  2,   0,  2,   0],
+                            [-3/2, 0,   6,  0, -3/2],
+                            [  0,  2,   0,  2,   0],
+                            [  0,  0, -3/2, 0,   0]]) / 8
+bayer_b_at_r = bayer_r_at_b
+
+# High quality linear interpolation filters for quad mosaic
+quad_g_at_r = np.array([[ 0,  0,  0, 0, -1, -1, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0, -1, -1, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0,  2,  2, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0,  2,  2, 0, 0,  0,  0],
+                        [-1, -1,  2, 2,  4,  4, 2, 2, -1, -1],
+                        [-1, -1,  2, 2,  4,  4, 2, 2, -1, -1],
+                        [ 0,  0,  0, 0,  2,  2, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0,  2,  2, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0, -1, -1, 0, 0,  0,  0],
+                        [ 0,  0,  0, 0, -1, -1, 0, 0,  0,  0]]) / 32
+quad_g_at_b = quad_g_at_r
+quad_r_at_green_rrow_bcol = np.array([[ 0,  0,  0,  0, 1/2, 1/2,  0,  0,  0,  0],
+                                      [ 0,  0,  0,  0, 1/2, 1/2,  0,  0,  0,  0],
+                                      [ 0,  0, -1, -1,  0,   0,  -1, -1,  0,  0],
+                                      [ 0,  0, -1, -1,  0,   0,  -1, -1,  0,  0],
+                                      [-1, -1,  4,  4,  5,   5,   4,  4, -1, -1],
+                                      [-1, -1,  4,  4,  5,   5,   4,  4, -1, -1],
+                                      [ 0,  0, -1, -1,  0,   0,  -1, -1,  0,  0],
+                                      [ 0,  0, -1, -1,  0,   0,  -1, -1,  0,  0],
+                                      [ 0,  0,  0,  0, 1/2, 1/2,  0,  0,  0,  0],
+                                      [ 0,  0,  0,  0, 1/2, 1/2,  0,  0,  0,  0]]) / 32
+quad_r_at_green_brow_rcol = quad_r_at_green_rrow_bcol.T
+quad_b_at_green_rrow_bcol = quad_r_at_green_brow_rcol
+quad_b_at_green_brow_rcol = quad_r_at_green_rrow_bcol
+quad_r_at_b = np.array([[  0,    0,  0,  0, -3/2, -3/2, 0,  0,   0,    0],
+                        [  0,    0,  0,  0, -3/2, -3/2, 0,  0,   0,    0],
+                        [  0,    0,  2,  2,   0,    0,  2,  2,   0,    0],
+                        [  0,    0,  2,  2,   0,    0,  2,  2,   0,    0],
+                        [-3/2, -3/2, 0,  0,   6,    6,  0,  0, -3/2, -3/2],
+                        [-3/2, -3/2, 0,  0,   6,    6,  0,  0, -3/2, -3/2],
+                        [  0,    0,  2,  2,   0,    0,  2,  2,   0,    0],
+                        [  0,    0,  2,  2,   0,    0,  2,  2,   0,    0],
+                        [  0,    0,  0,  0, -3/2, -3/2, 0,  0,   0,    0],
+                        [  0,    0,  0,  0, -3/2, -3/2, 0,  0,   0,    0]]) / 32
+quad_b_at_r = quad_r_at_b
+
+def high_quality_linear_interpolation(op : CFA, y : np.ndarray) -> np.ndarray:
+    z = op.adjoint(y)
+    res = np.empty(op.input_shape)
+    mask = op.mask
+    R = 0 ; G = 1 ; B = 2
+
+    if op.cfa == 'bayer':
+        y_pad = np.pad(y, 2, 'constant', constant_values=0)
+
+        for i in range(y.shape[0]):
+            for j in range(y.shape[1]):
+                i_pad = i+2 ; j_pad = j+2
+
+                if mask[i, j, G] == 1: # We must estimate R and B components
+                    res[i, j, G] = y[i, j]
+
+                    # Estimate R (B at left and right or B at top and bottom)
+                    if mask[i, max(0, j-1), R] == 1 or mask[i, min(y.shape[1]-1, j+1), R] == 1:
+                        res[i, j, R] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_r_at_green_rrow_bcol, mode='valid'))
+                    else:
+                        res[i, j, R] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_r_at_green_brow_rcol, mode='valid'))
+
+                    # Estimate B (R at left and right or R at top and bottom)
+                    if mask[i, max(0, j-1), B] == 1 or mask[i, min(y.shape[1]-1, j+1), B] == 1:
+                        res[i, j, B] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_b_at_green_brow_rcol, mode='valid'))
+                    else:
+                        res[i, j, B] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_b_at_green_rrow_bcol, mode='valid'))
+
+                elif mask[i, j, R] == 1: # We must estimate G and B components
+                    res[i, j, R] = y[i, j]
+                    res[i, j, G] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_g_at_r, mode='valid'))
+                    res[i, j, B] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_b_at_r, mode='valid'))
+
+                else: # We must estimate R and G components
+                    res[i, j, B] = y[i, j]
+                    res[i, j, G] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_g_at_b, mode='valid'))
+                    res[i, j, R] = float(convolve2d(y_pad[i_pad-2:i_pad+3, j_pad-2:j_pad+3], bayer_r_at_b, mode='valid'))
+
+    else:
+        y_pad = np.pad(y, 4, 'constant', constant_values=0)
+
+        for i in range(0, y.shape[0], 2):
+            for j in range(0, y.shape[1], 2):
+                i_pad = i+4 ; j_pad = j+4
+
+                if mask[i, j, G] == 1: # We must estimate R and B components
+                    res[i:i+2, j:j+2, G] = y[i:i+2, j:j+2]
+
+                    # Estimate R (B at left and right or B at top and bottom)
+                    if mask[i, max(0, j-1), R] == 1 or mask[i, min(y.shape[1]-1, j+1), R] == 1:
+                        res[i:i+2, j:j+2, R] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_r_at_green_rrow_bcol, mode='valid'))
+                    else:
+                        res[i:i+2, j:j+2, R] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_r_at_green_brow_rcol, mode='valid'))
+
+                    # Estimate B (R at left and right or R at top and bottom)
+                    if mask[i, max(0, j-1), B] == 1 or mask[i, min(y.shape[1]-1, j+1), B] == 1:
+                        res[i:i+2, j:j+2, B] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_b_at_green_brow_rcol, mode='valid'))
+                    else:
+                        res[i:i+2, j:j+2, B] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_b_at_green_rrow_bcol, mode='valid'))
+
+                elif mask[i, j, R] == 1: # We must estimate G and B components
+                    res[i:i+2, j:j+2, R] = y[i:i+2, j:j+2]
+                    res[i:i+2, j:j+2, G] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_g_at_r, mode='valid'))
+                    res[i:i+2, j:j+2, B] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_b_at_r, mode='valid'))
+
+                else: # We must estimate R and G components
+                    res[i:i+2, j:j+2, B] = y[i:i+2, j:j+2]
+                    res[i:i+2, j:j+2, G] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_g_at_b, mode='valid'))
+                    res[i:i+2, j:j+2, R] = float(convolve2d(y_pad[i_pad-4:i_pad+6, j_pad-4:j_pad+6], quad_r_at_b, mode='valid'))
+
+    return res
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

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

src/methods/marty/high_quality_interpolation.py

+import numpy as np
+
+
+def is_green(i, j):
+    return (not i%2 and not j%2) or (i%2 and j%2)
+
+def is_red(i, j):
+    return not i%2 and j%2
+
+def is_blue(i, j):
+    return i%2 and not j%2
+
+
+g_ker = 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
+
+rgrrow_ker = 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
+
+rgrcol_ker = 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]
+]) / 8
+
+rb_ker = 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
+
+
+
+def high_quality_interpolation(image):
+
+    res = np.empty((image.shape[0], image.shape[1], 3))
+    
+    estimated_green = np.zeros_like(image)
+    estimated_red = np.zeros_like(image)
+    estimated_blue = np.zeros_like(image)
+
+    for i in range(2, image.shape[0]-2):
+        for j in range(2, image.shape[1]-2):
+            if (is_red(i, j)):
+                estimated_red[i, j] = image[i, j]
+                estimated_green[i, j] = (image[i-2:i+3, j-2:j+3] * g_ker).sum()
+                estimated_blue[i, j] = (image[i-2:i+3, j-2:j+3] * rb_ker).sum()
+            elif (is_blue(i, j)):
+                estimated_red[i, j] = (image[i-2:i+3, j-2:j+3] * rb_ker).sum()
+                estimated_green[i, j] = (image[i-2:i+3, j-2:j+3] * g_ker).sum()
+                estimated_blue[i, j] = image[i, j]
+            else:
+                estimated_green[i, j] = image[i, j]
+                if (is_red(i-1, j)):
+                    estimated_red[i, j] = (image[i-2:i+3, j-2:j+3] * rgrcol_ker).sum()
+                    estimated_blue[i, j] = (image[i-2:i+3, j-2:j+3] * rgrrow_ker).sum()
+                else:
+                    estimated_red[i, j] = (image[i-2:i+3, j-2:j+3] * rgrrow_ker).sum()
+                    estimated_blue[i, j] = (image[i-2:i+3, j-2:j+3] * rgrcol_ker).sum()
+
+
+    res[:, :, 0] = estimated_red.clip(0, 1)
+    res[:, :, 1] = estimated_green.clip(0, 1)
+    res[:, :, 2] = estimated_blue.clip(0, 1)
+    
+    return res
+

src/methods/marty/mdwi.py

+import numpy as np
+
+
+def is_green(i, j):
+    return (not i%2 and not j%2) or (i%2 and j%2)
+
+def is_red(i, j):
+    return not i%2 and j%2
+
+def is_blue(i, j):
+    return i%2 and not j%2
+
+
+
+def MDWI(image):
+    
+    res = np.empty((image.shape[0], image.shape[1], 3))
+    estimated_green = np.zeros_like(image)
+    estimated_red = np.zeros_like(image)
+    estimated_blue = np.zeros_like(image)
+
+    EPSILON = 10e-4
+
+    # Estimation of green plan first
+
+    pre_estimation = {}
+    diag_gradient_factor = {}
+    weight = {}
+    H = [-8/256, 23/256, -48/256, 161/256, 161/256, -48/256, 23/256, -8/256]
+    coord_NE = [(-4, -3), (-3, -2), (-2, -1), (-1, 0), (0, 1), (1, 2), (2, 3), (3, 4)]
+    coord_SE = [(-3, -4), (-2, -3), (-1, -2), (0, -1), (1, 0), (2, 1), (3, 2), (4, 3)]
+    coord_NW = [(3, -4), (2, -3), (1, -2), (0, -1), (-1, 0), (-2, 1), (-3, 2), (-4, 3)]
+    coord_SW = [(4, -3), (3, -2), (2, -1), (1, 0), (0, 1), (-1, 2), (-2, 3), (-3, 4)]
+
+    for i in range(3, image.shape[0] - 4):
+        for j in range(3, image.shape[1] - 4):
+            
+
+            if (not is_green(i, j)):       
+                pre_estimation["N"] = image[i-1, j] + (image[i, j] - image[i-2, j]) / 2   
+                pre_estimation["S"] = image[i+1, j] + (image[i, j] - image[i+2, j]) / 2
+                pre_estimation["W"] = image[i, j-1] + (image[i, j] - image[i, j-2]) / 2
+                pre_estimation["E"] = image[i, j+1] + (image[i, j] - image[i, j+2]) / 2
+
+                diag_gradient_factor["N"] = abs(image[i-2, j-1] - image[i, j-1]) + abs(image[i-3, j] - image[i-1, j]) + abs(image[i-2, j+1] - image[i, j+1]) + abs(image[i-3, j-1] - image[i-1, j-1]) + abs(image[i-3, j+1] - image[i-1, j+1]) + abs(image[i-2, j] - image[i, j]) + EPSILON
+                diag_gradient_factor["S"] = abs(image[i+2, j-1] - image[i, j-1]) + abs(image[i+3, j] - image[i+1, j]) + abs(image[i+2, j+1] - image[i, j+1]) + abs(image[i+3, j-1] - image[i+1, j-1]) + abs(image[i+3, j+1] - image[i+1, j+1]) + abs(image[i+2, j] - image[i, j]) + EPSILON
+                diag_gradient_factor["W"] = abs(image[i-1, j-2] - image[i-1, j]) + abs(image[i, j-3] - image[i, j-1]) + abs(image[i+1, j-2] - image[i+1, j]) + abs(image[i-1, j-3] - image[i-1, j-1]) + abs(image[i+1, j-3] - image[i+1, j-1]) + abs(image[i, j-2] - image[i, j]) + EPSILON
+                diag_gradient_factor["E"] = abs(image[i-1, j+2] - image[i-1, j]) + abs(image[i, j+3] - image[i, j+1]) + abs(image[i+1, j+2] - image[i+1, j]) + abs(image[i-1, j+3] - image[i-1, j+1]) + abs(image[i+1, j+3] - image[i+1, j+1]) + abs(image[i, j+2] - image[i, j]) + EPSILON
+
+
+                pre_estimation["NW"] = 0
+                pre_estimation["SW"] = 0
+                pre_estimation["NE"] = 0
+                pre_estimation["SE"] = 0
+                for k in range(8):
+                    pre_estimation["NW"] += image[i + coord_NW[k][0], j + coord_NW[k][1]] * H[k]
+                    pre_estimation["SW"] += image[i + coord_SW[k][0], j + coord_SW[k][1]] * H[k]
+                    pre_estimation["NE"] += image[i + coord_NE[k][0], j + coord_NE[k][1]] * H[k]
+                    pre_estimation["SE"] += image[i + coord_SE[k][0], j + coord_SE[k][1]] * H[k]
+                
+                diag_gradient_factor["NW"] = abs(image[i-2, j-1] - image[i-1, j]) + abs(image[i-1, j] - image[i, j+1]) + abs(image[i-1, j-2] - image[i, j-1]) + abs(image[i, j-1] - image[i+1, j]) + abs(image[i-1, j-1] - image[i+1, j+1]) + abs(image[i-2, j-2] - image[i, j]) + EPSILON
+                diag_gradient_factor["SW"] = abs(image[i-1, j] - image[i, j-1]) + abs(image[i, j+1] - image[i+1, j]) + abs(image[i-1, j+1] - image[i+1, j-1]) + abs(image[i, j-1] - image[i+1, j-2]) + abs(image[i+1, j] - image[i+2, j-1]) + abs(image[i, j] - image[i+2, j-2]) + EPSILON
+                diag_gradient_factor["NE"] = abs(image[i-1, j] - image[i, j-1]) + abs(image[i, j+1] - image[i+1, j]) + abs(image[i-1, j+1] - image[i+1, j-1]) + abs(image[i-1, j] - image[i-2, j+1]) + abs(image[i, j+1] - image[i-1, j+2]) + abs(image[i, j] - image[i-2, j+2]) + EPSILON
+                diag_gradient_factor["SE"] = abs(image[i+2, j+1] - image[i+1, j]) + abs(image[i-1, j] - image[i, j+1]) + abs(image[i+1, j+2] - image[i, j+1]) + abs(image[i, j-1] - image[i+1, j]) + abs(image[i-1, j-1] - image[i+1, j+1]) + abs(image[i+2, j+2] - image[i, j]) + EPSILON
+
+                
+                for k, v in diag_gradient_factor.items():
+                    weight[k] =  1 / v
+
+                num = 0
+                den = 0
+                for k, v in pre_estimation.items():
+                    num += v * weight[k]
+                    den += weight[k]
+
+                estimated_green[i, j] = num / den
+
+            else:
+                estimated_green[i, j] = image[i, j]
+
+
+    # Then we work on blue and red plans
+    
+    diag = [(-1, -1), (-1, 1), (1, -1), (1, 1)]
+    
+    for i in range(3, image.shape[0] - 4):
+        for j in range(3, image.shape[1] - 4):
+            if (is_red(i, j) or is_blue(i, j)):
+                ksi_nw = image[i-1, j-1] - estimated_green[i-1, j-1]
+                ksi_ne = image[i-1, j+1] - estimated_green[i-1, j+1]
+                ksi_sw = image[i+1, j-1] - estimated_green[i+1, j-1]
+                ksi_se = image[i+1, j+1] - estimated_green[i+1, j+1]
+                
+                g_nw = abs(estimated_green[i-2, j-1] - estimated_green[i-1, j]) + abs(estimated_green[i-1, j-2] - estimated_green[i, j-1]) + abs(estimated_green[i-2, j-2] - estimated_green[i-1, j-1]) + abs(estimated_green[i-1, j-1] - estimated_green[i, j]) + abs(image[i-1, j-1] - image[i+1, j+1]) + EPSILON
+                g_ne = abs(estimated_green[i-2, j+1] - estimated_green[i-1, j]) + abs(estimated_green[i-1, j+2] - estimated_green[i, j+1]) + abs(estimated_green[i-2, j+2] - estimated_green[i-1, j+1]) + abs(estimated_green[i-1, j+1] - estimated_green[i, j]) + abs(image[i-1, j+1] - image[i+1, j-1]) + EPSILON
+                g_sw = abs(estimated_green[i+2, j-1] - estimated_green[i+1, j]) + abs(estimated_green[i+1, j-2] - estimated_green[i, j-1]) + abs(estimated_green[i+2, j-2] - estimated_green[i+1, j-1]) + abs(estimated_green[i+1, j-1] - estimated_green[i, j]) + abs(image[i+1, j-1] - image[i-1, j+1]) + EPSILON
+                g_se = abs(estimated_green[i+2, j+1] - estimated_green[i+1, j]) + abs(estimated_green[i+1, j+2] - estimated_green[i, j+1]) + abs(estimated_green[i+2, j+2] - estimated_green[i+1, j+1]) + abs(estimated_green[i+1, j+1] - estimated_green[i, j]) + abs(image[i+1, j+1] - image[i-1, j-1]) + EPSILON
+
+                pre_est = estimated_green[i, j] + (ksi_nw/g_nw + ksi_ne/g_ne + ksi_sw/g_sw + ksi_se/g_se) / (1/g_nw + 1/g_ne + 1/g_sw + 1/g_se)
+
+                sum_dif = 0
+                sum_sig = 0
+
+                for coord in diag:
+                    mu = image[i + coord[0], j + coord[1]] - pre_est
+                    sum_dif += mu / (1 + abs(mu))
+                    sum_sig += 1 / (1 + abs(mu))
+
+                if (is_red(i, j)):
+                    estimated_blue[i, j] = pre_est
+                    estimated_red[i, j] = image[i, j]
+                else:
+                    estimated_red[i, j] = pre_est
+                    estimated_blue[i, j] = image[i, j]
+
+            else:
+                if (is_red(i-1, j)):
+                    estimated_red[i, j] = (image[i-1, j] + image[i+1, j]) / 2
+                    estimated_blue[i, j] = (image[i, j-1] + image[i, j+1]) / 2
+                else:
+                    estimated_red[i, j] = (image[i, j-1] + image[i, j+1]) / 2
+                    estimated_blue[i, j] = (image[i-1, j] + image[i+1, j]) / 2
+
+
+    res[:, :, 0] = estimated_red
+    res[:, :, 1] = estimated_green
+    res[:, :, 2] = estimated_blue
+
+    return res
\ No newline at end of file

src/methods/marty/quad_bayer.py

+import numpy as np
+
+from scipy.signal import convolve2d
+
+
+def down_sample(image):    
+
+    down_sampled = np.empty((int(image.shape[0] / 2), int(image.shape[1] / 2)))
+
+    down_sampled[:, :] = (image[::2, ::2] + image[1::2, ::2] + image[::2, 1::2] + image[1::2, 1::2]) / 4
+
+    return down_sampled
+
+
+def up_sample(image):
+
+    up_sampled = np.empty((int(image.shape[0] * 2), int(image.shape[1] * 2), image.shape[2]))
+
+    up_sampled[::2, ::2, :] = image[::1, ::1, :]
+    up_sampled[1::2, ::2, :] = image[::1, ::1, :]
+    up_sampled[::2, 1::2, :] = image[::1, ::1, :]
+    up_sampled[1::2, 1::2, :] = image[::1, ::1, :]
+
+    return up_sampled
+
+
+
+def refine(image):
+    
+    res = np.empty((image.shape[0], image.shape[1], 3))
+    
+    ker = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 16
+    for i in range(3):
+        res[:, :, i] = convolve2d(image[:, :, i], ker, mode='same')
+
+    return res

src/methods/marty/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.marty.mdwi import MDWI
+from src.methods.marty.high_quality_interpolation import high_quality_interpolation
+from src.methods.marty.quad_bayer import down_sample, up_sample, refine
+
+
+# Either MDWI or high_quality_interpolation
+method = high_quality_interpolation
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    
+    # Performing the reconstruction.
+    if (cfa == "bayer") :
+        
+        res = method(y)
+        
+
+    else :
+        down_sample_image = down_sample(y)
+
+        bayer_image = method(down_sample_image)
+        
+        res = up_sample(bayer_image)
+
+        res = refine(res)
+    return res
+
+
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/maxime_torre/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).
+"""
+
+from src.forward_model import CFA
+import numpy as np
+from scipy.signal import convolve2d
+import cv2 
+
+def is_green(z,i, j):
+
+    return z[i, j, 1] != 0
+
+def hamilton_adams_interpolation(y, op, z):
+    height, width = y.shape
+    green_channel = np.copy(z[:, :, 1])  
+
+    for i in range(1, height-1):
+        for j in range(1, width-1):
+            if not is_green(z,i, j) :
+                delta_H = abs(z[i, j-1, 1] - z[i, j+1, 1]) + abs(z[i, j-1, 0] - z[i, j+1, 0] + z[i, j-1, 2] - z[i, j+1, 2]) / 2
+                # print(f"delta_H : {delta_H}")
+                delta_V = abs(z[i-1, j, 1] - z[i+1, j, 1]) + abs(z[i-1, j, 0] - z[i+1, j, 0] + z[i-1, j, 2] - z[i+1, j, 2]) / 2
+
+                if delta_H > delta_V:
+                    green_channel[i, j] = (z[i-1, j, 1] + z[i+1, j, 1]) / 2 + (z[i, j-1, 0] - z[i, j+1, 0] + z[i, j-1, 2] - z[i, j+1, 2]) / 4
+                elif delta_H < delta_V:
+                    green_channel[i, j] = (z[i, j-1, 1] + z[i, j+1, 1]) / 2 + (z[i-1, j, 0] - z[i+1, j, 0] + z[i-1, j, 2] - z[i+1, j, 2]) / 4
+                else:
+                    green_channel[i, j] = (z[i-1, j, 1] + z[i+1, j, 1] + z[i, j-1, 1] + z[i, j+1, 1]) / 4 + \
+                                            (z[i, j-1, 0] - z[i, j+1, 0] + z[i, j-1, 2] - z[i, j+1, 2] + \
+                                            z[i-1, j, 0] - z[i+1, j, 0] + z[i-1, j, 2] - z[i+1, j, 2]) / 8
+
+    return green_channel
+
+def interpolate_channel_difference(mosaicked_channel, green_channel_interpolated):
+    ker_bayer_red_blue = np.array([[1, 2, 1], [2, 4, 2], [1, 2, 1]]) / 4
+
+    print(mosaicked_channel.shape, green_channel_interpolated.shape)
+    difference = mosaicked_channel - green_channel_interpolated
+    difference_interpolated = convolve2d(difference, np.ones((3, 3)) / 9, mode='same', boundary='wrap')
+    channel_interpolated = green_channel_interpolated + difference_interpolated
+
+    channel_interpolated = convolve2d(channel_interpolated, ker_bayer_red_blue, mode='same')
+
+    return channel_interpolated
+
+def Constant_difference_based_interpolation_reconstruction(op, y, z):
+    if op.cfa == 'bayer':
+        print("bayer")
+
+        red_channel = z[:, :, 0]
+        green_channel = z[:, :, 1]
+        blue_channel = z[:, :, 2]
+
+        green_channel_reconstruct = hamilton_adams_interpolation(y, op, z)
+            
+        red_channel_interpolated = interpolate_channel_difference(red_channel, green_channel_reconstruct)
+        blue_channel_interpolated = interpolate_channel_difference(blue_channel, green_channel_reconstruct)
+
+        reconstructed_image = np.stack((red_channel_interpolated, green_channel_reconstruct, blue_channel_interpolated), axis=-1)
+
+        return reconstructed_image
+    
+    elif op.cfa == "quad_bayer":
+        print(f"quad_bayer")
+        new_z = cv2.resize(z, (z.shape[1] // 2, z.shape[0] // 2), interpolation=cv2.INTER_AREA)
+        new_y=np.sum(new_z, axis=2)
+        op.mask = op.mask[::2, ::2]
+        green_channel_reconstruct_new = hamilton_adams_interpolation(new_y, op, new_z)
+        red_channel_new = new_z[:, :, 0]
+        blue_channel_new = new_z[:, :, 2]
+        red_channel_interpolated_new = interpolate_channel_difference(red_channel_new, green_channel_reconstruct_new)
+        blue_channel_interpolated_new = interpolate_channel_difference(blue_channel_new, green_channel_reconstruct_new)
+        reconstructed_image_new = np.stack((red_channel_interpolated_new, green_channel_reconstruct_new, blue_channel_interpolated_new), axis=-1)
+        reconstructed_image_upsampled = cv2.resize(reconstructed_image_new, (z.shape[1], z.shape[0]), interpolation=cv2.INTER_LINEAR)
+        
+        return reconstructed_image_upsampled
+
+    else :
+        raise ValueError("CFA pattern not recognized")
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+    z = op.adjoint(y)
+    reconstructed_image = Constant_difference_based_interpolation_reconstruction(op, y, z)
+    
+    return reconstructed_image
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller