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/nada_kouddane/functions.py

+"""A file containing the main function for the chosen interpolation (Adams-Hamilton Algorithm).
+"""
+
+### BE CAREFUL TO ADD THE MODULE TORCH TO THE FILE requirements.txt
+
+import numpy as np
+from scipy.signal import convolve2d
+import torch
+import torch.nn.functional as F
+
+from src.forward_model import CFA
+
+
+# Initial kernels for the bilinear interpolation
+ker_bayer_red_blue = np.array([[1, 2, 1], 
+                               [2, 4, 2], 
+                               [1, 2, 1]]) / 4
+ker_bayer_green = np.array([[0, 1, 0], 
+                            [1, 4, 1], 
+                            [0, 1, 0]]) / 4
+
+# Kernels for the quad_bayer filter
+ker_quad_red_blue = np.array([[1, 1, 2, 2, 0, 1], 
+                              [1, 1, 2, 2, 1, 1], 
+                              [2, 2, 4, 4, 2, 2], 
+                              [2, 2, 4, 4, 2, 2], 
+                              [1, 1, 2, 2, 1, 1], 
+                              [1, 1, 2, 2, 1, 1]])/16
+ker_quad_green = np.array([[0, 0, 1, 1, 0, 0], 
+                           [0, 0, 1, 1, 0, 0], 
+                           [1, 1, 4, 4, 1, 1], 
+                           [1, 1, 4, 4, 1, 1], 
+                           [0, 0, 1, 1, 0, 0], 
+                           [0, 0, 1, 1, 0, 0]])/16
+
+
+def second_naive_interpolation(op: CFA, y: np.ndarray) -> np.ndarray:
+    """Performs a naive interpolation of the lost pixels. for the quad_bayer, performs a convolution 
+    with kernels (above: ker_quad_red_blue and ker_quad_green) inspired from the classic bilnear interpolation
+
+    Args:
+        op (CFA): CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    z = op.adjoint(y)
+
+    if op.cfa == 'bayer':
+        res = np.empty(op.input_shape)
+
+        res[:, :, 0] = convolve2d(z[:, :, 0], ker_bayer_red_blue, mode='same')
+        res[:, :, 1] = convolve2d(z[:, :, 1], ker_bayer_green, mode='same')
+        res[:, :, 2] = convolve2d(z[:, :, 2], ker_bayer_red_blue, mode='same')
+
+    else:
+        res = np.empty(op.input_shape)
+
+        res[:, :, 0] = convolve2d(z[:, :, 0], ker_quad_red_blue, mode='same')
+        res[:, :, 1] = convolve2d(z[:, :, 1], ker_quad_green, mode='same')
+        res[:, :, 2] = convolve2d(z[:, :, 2], ker_quad_red_blue, mode='same')
+        #res[:, :, 0] = convolution_pad_stride(z[:, :, 0], ker_quad_red_blue)
+        #res[:, :, 1] = convolution_pad_stride(z[:, :, 1], ker_quad_green)
+        #res[:, :, 2] = convolution_pad_stride(z[:, :, 2], ker_quad_red_blue)
+
+
+    return res
+
+
+import torch
+import torch.nn.functional as F
+
+def convolution_pad_stride(input, kernel):
+    """An attempt of convolution with stride and padding.
+
+    Args:
+        input(np.ndarray): input (mosaicked) image.
+        kernel (np.ndarray): convolution kernel.
+
+    Returns:
+        np.ndarray: Convolution of the image with the kernel with stride of 2 and padding of 512.
+    """
+    input_tensor = torch.tensor(np.expand_dims(input, axis=(0,1)))
+    kernel_tensor = torch.tensor(np.expand_dims(kernel, axis=(0,1)))
+    output = F.conv2d(input_tensor, kernel_tensor, stride=2, padding=514)
+    output = output.numpy().squeeze()
+
+    return output
+
+
+import numpy as np
+
+
+
+
+
+def take(array2d, i, j):
+    """
+    Helper function that returns the indices of an array, prevents index outbound
+
+    Args:
+        array2d (np.ndarray): input array
+        i: row index
+        j: column index
+
+    Returns:
+        np.float64(0): array value at position[i,j]
+    """
+    if 0 <= i < array2d.shape[0] and 0 <= j < array2d.shape[1]:
+        return array2d[i, j]
+    return np.float64(0)
+
+
+def red_blue_positions(img):
+    """
+    Helper function that yields the red and blue positions in the mosaicked image
+
+    Args:
+        img (np.ndarray): mosaicked image 
+
+    Returns:
+        None
+    """
+    first_non_green = 1
+    for i in range(img.shape[0]):
+        for j in range(first_non_green, img.shape[1], 2):
+            yield i, j
+        first_non_green = 1 - first_non_green
+
+
+def directional_green_interpolation(img):
+    """
+    Function that performs green interpolation in horizontal and vertical directions.
+
+    Args: 
+        img (np.ndarray): input (mosaicked) image
+
+    Returns:
+        green_h, green_v (np.ndarray): horizontally and vertically interpolated green components
+    """
+
+    green_h = img.copy()  # green positions are copied
+    green_v = img.copy()  # other values will be replaced
+
+    for i, j in red_blue_positions(img):
+        r = lambda k: take(img, i, j + k)  # r - relative indexing
+        green_h[i, j] = (r(1) + r(-1) + r(0)) / 2 - (r(2) + r(-2)) / 4
+        r = lambda k: take(img, i + k, j)
+        green_v[i, j] = (r(1) + r(-1) + r(0)) / 2 - (r(2) + r(-2)) / 4
+
+    return green_h, green_v
+
+
+def green_decision(img, green_h, green_v, cardinal_directions_improvement = True):
+    
+    """
+    Function that performs the green decision between the chrominance components based on the color difference uniformity
+    by calculating the horizontal and the vertical gradients.
+
+    Args: 
+        img (np.ndarray): input (mosaicked) image
+        green_h (np.ndarray): horizontally interpolated green component
+        green_v (np.ndarray): vertically interpolated green component
+        cardinal_directions_improvement (bool) = True (default) : parameter that allows to adjust a given window weight in order
+        to improve the gradients calculation
+
+    Returns:
+        green (np.ndarray): interpolated green image
+        delta_h (np.ndarray): horizontal gradient image
+        delta_v (np.ndarray): vertical gradient image
+    """
+
+    height, width = img.shape
+
+    # "chrominance" is R - G in red locations, B - G in blue locations
+    # and 0 in green locations
+    chrominance_h = img - green_h
+    chrominance_v = img - green_v
+
+    # also 0 in green locations, this will be useful
+    gradient_h = chrominance_h.copy()
+    gradient_v = chrominance_v.copy()
+
+    for i, j in red_blue_positions(img):
+        gradient_h[i, j] -= take(chrominance_h, i, j + 2)
+        gradient_v[i, j] -= take(chrominance_v, i + 2, j)
+    gradient_h = np.abs(gradient_h)
+    gradient_v = np.abs(gradient_v)
+    # could be easily rewritten without loops
+
+    window = np.ones(shape=(5, 5), dtype=np.float64)
+    if cardinal_directions_improvement:
+        window[2, :] = 3
+        window[:, 2] = 3
+
+    delta_h = np.zeros(shape=(img.shape), dtype=np.float64)
+    delta_v = delta_h.copy()
+    padded_grad_h = np.zeros(shape=(img.shape[0] + 4, img.shape[1] + 4), dtype=np.float64)
+    padded_grad_v = padded_grad_h.copy()
+    padded_grad_h[2 : img.shape[0] + 2, 2 : img.shape[1] + 2] = gradient_h
+    padded_grad_v[2 : img.shape[0] + 2, 2 : img.shape[1] + 2] = gradient_v
+    green = green_h.copy()
+    for i, j in red_blue_positions(img):
+        delta_h[i, j] = np.sum(window * padded_grad_h[i : i + 5, j : j + 5])
+        delta_v[i, j] = np.sum(window * padded_grad_v[i : i + 5, j : j + 5])
+        if delta_v[i, j] < delta_h[i, j]:
+            green[i, j] = green_v[i, j]
+
+    return green, delta_h, delta_v
+
+
+def red_blue_interpolation(img, green, delta_h, delta_v):
+
+    """
+    Function that performs the red and blue components interpolation.
+
+    Args: 
+        img (np.ndarray): input (mosaicked) image
+        green (np.ndarray): interpolated green image
+        delta_h (np.ndarray): horizontal gradient image
+        delta_v (np.ndarray): vertical gradient image
+
+    Returns:
+        red, blue (np.ndarray): interpolated red and blue image components
+    """
+
+    height, width = img.shape
+    red = img.copy()
+    blue = img.copy()
+
+    # green positions first
+    for i in range(0, height, 2):  # green-red rows
+        for j in range(0, width, 2):
+            red[i, j] = (take(img, i, j - 1) +
+                         take(img, i, j + 1)) / 2
+            blue[i, j] = (take(img, i - 1, j) +
+                          take(img, i + 1, j)) / 2
+    for i in range(1, height, 2):  # green-blue rows
+        for j in range(1, width, 2):
+            blue[i, j] = (take(img, i, j - 1) +
+                          take(img, i, j + 1)) / 2
+            red[i, j] = (take(img, i - 1, j) +
+                         take(img, i + 1, j)) / 2
+    
+    # now red in blue positions, blue in red positions
+    red_minus_blue = red - blue
+    for i in range(1, height, 2):
+        for j in range(0, width, 2):
+            if delta_v[i, j] < delta_h[i, j]:
+                red[i, j] = blue[i, j] + (take(red_minus_blue, i - 1, j) +
+                                          take(red_minus_blue, i + 1, j)) / 2
+            else:
+                red[i, j] = blue[i, j] + (take(red_minus_blue, i, j - 1) +
+                                          take(red_minus_blue, i, j + 1)) / 2
+    for i in range(0, height, 2):
+        for j in range(1, width, 2):
+            if delta_v[i, j] < delta_h[i, j]:
+                blue[i, j] = red[i, j] - (take(red_minus_blue, i - 1, j) +
+                                          take(red_minus_blue, i + 1, j)) / 2
+            else:
+                blue[i, j] = red[i, j] - (take(red_minus_blue, i, j - 1) +
+                                          take(red_minus_blue, i, j + 1)) / 2
+
+    return red, blue
+
+
+def demosaicking_algorithm(img):
+    """
+    Main function of the Daniele Menon demosaicking algorithm.
+
+    Args:
+        img (np.ndarray): input (mosaicked image)
+
+    Returns:
+        np.ndarray: reconstructed image
+
+    """
+    green_h, green_v = directional_green_interpolation(img)
+    green, delta_h, delta_v = green_decision(img, green_h, green_v)
+    red, blue = red_blue_interpolation(img, green, delta_h, delta_v)
+    return np.clip(np.dstack((red, green, blue)), 0, 1)
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller
\ No newline at end of file

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

src/methods/quelletl/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.quelletl.some_function import interpolation
+import cv2
+from sklearn.ensemble import RandomForestRegressor
+from sklearn.model_selection import train_test_split
+from sklearn.metrics import mean_squared_error
+from scipy.signal import convolve2d
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+    res = interpolation(y, op)
+    
+    return res
+
+def quad_to_bayer(y):
+    for i in range(1, y.shape[0], 4):
+        save = np.copy(y[:,i])
+        y[:,i] = y[:,i+1]
+        y[:,i+1] = save
+    for j in range(1, y.shape[0], 4):
+        save = np.copy(y[j,:])
+        y[j,:] = y[j+1,:]
+        y[j+1,:] = save
+    for i in range(1, y.shape[0], 4):
+        for j in range(1, y.shape[0], 4):
+            save = np.copy(y[i,j])
+            y[i,j] = y[i+1,j+1]
+            y[i+1,j+1] = save
+    return y
+
+ker_bayer_green = np.array([[0, 1, 0], [1, 4, 1], [0, 1, 0]]) / 4
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/quelletl/some_function.py

+import numpy as np
+
+def bayer_blue_red(res, y, i, j) : 
+    """
+    Compute estimated blue/red pixel in red/blue bayer pixel
+
+    Args : 
+        res     : estimated image
+        y       : image to reconstruct
+        i, j    : indices
+    Return : 
+        value   : value of the estimated pixels
+    """
+    K2 = res[i-1,j-1,1] - y[i-1,j-1]
+    K4 = res[i-1,j+1,1] - y[i-1,j+1]
+    K10 = res[i+1,j-1,1] - y[i+1,j-1]
+    K12 = res[i+1,j+1,1] - y[i+1,j+1]
+    value = res[i,j,1]-1/4*(K2+K4+K10+K12)
+    return value
+
+def bayer_green_vert(res, y, i, j) : 
+    """
+    Compute estimated blue/red pixel in green bayer pixel in vertical direction
+
+    Args : 
+        res     : estimated image
+        y       : image to reconstruct
+        i, j    : indices
+    Return : 
+        value   : value of the estimated pixels
+    """
+    k1 = res[i-1,j,1] - y[i-1,j]
+    k2 = res[i+1,j,1] - y[i+1,j]
+    value = y[i,j] - 1/2*(k1+k2)
+    return value
+
+def bayer_green_hor(res, y, i, j): 
+    """
+    Compute estimated blue/red pixel in green bayer pixel in horizontal direction
+
+    Args : 
+        res     : estimated image
+        y       : image to reconstruct
+        i, j    : indices
+    Return : 
+        value   : value of the estimated pixels
+    """
+    k1 = res[i,j-1,1] - y[i,j-1]
+    k2 = res[i,j+1,1] - y[i,j+1]
+    value = y[i,j] - 1/2*(k1+k2)
+    return value
+
+def interpolate_green(res, y, z):
+    """
+    Directional interpolation of the green channel
+
+    Args : 
+        res     : estimated image
+        y       : image to reconstruct
+        z       : bayer pattern
+    Return : 
+        res     : reconstructed image
+    """
+    for i in range(2,y.shape[0]-1):
+        for j in range(2,y.shape[1]-1):
+        # Vertical and horizontal gradient
+            if z[i,j,1]==0:
+                d_h = np.abs(y[i,j-1]-y[i,j+1])
+                d_v = np.abs(y[i-1,j]-y[i+1,j])
+                if d_h > d_v:
+                    green = (y[i-1,j]+y[i+1,j])/2
+                elif d_v > d_h:
+                    green = (y[i,j-1]+y[i,j+1])/2
+                else :
+                    green = (y[i,j-1]+y[i,j+1]+y[i-1,j]+y[i+1,j])/4
+            else : 
+                green = y[i,j]
+            res[i,j,1] = green
+    return res
+
+def quad_to_bayer(y):
+    """
+    Convert Quad Bayer to Bayer
+
+    Args : 
+        res     : estimated image
+        y       : image to reconstruct
+        i, j    : indices
+    Return : 
+        value   : value of the estimated pixels
+    """
+    for i in range(1, y.shape[0], 4):
+        save = np.copy(y[:,i])
+        y[:,i] = y[:,i+1]
+        y[:,i+1] = save
+    for j in range(1, y.shape[0], 4):
+        save = np.copy(y[j,:])
+        y[j,:] = y[j+1,:]
+        y[j+1,:] = save
+    for i in range(1, y.shape[0], 4):
+        for j in range(1, y.shape[0], 4):
+            save = np.copy(y[i,j])
+            y[i,j] = y[i+1,j+1]
+            y[i+1,j+1] = save
+    return y
+
+def interpolation(y, op):
+    """
+    Reconstruct image
+
+    Args : 
+        y       : image to reconstruct
+        op      : CFA operator
+    Return : 
+        np.ndarray: Demosaicked image.
+    """
+    if op.cfa == 'quad_bayer':
+        y = quad_to_bayer(y)
+        op.mask = quad_to_bayer(op.mask)
+
+    z = op.adjoint(y)
+    res = np.empty(op.input_shape)
+
+    # Interpolation of green channel
+    res = interpolate_green(res, y, z)
+    # Interpolation of R and B channels using channel correlation
+    for i in range(2,y.shape[0]-2):
+        for j in range(2, y.shape[1]-2):
+            # Bayer is Green
+            if z[i,j,1] != 0 : 
+                # Green is between 2 vertical bleu
+                if z[i+1,j,0] == 0:
+                    red = bayer_green_hor(res, y, i, j) # Compute Red channel
+                    blue = bayer_green_vert(res, y, i, j) # Compute Blue channel 
+                # Green is between 2 vertical red
+                else:
+                    blue = bayer_green_hor(res, y, i, j) # Compute Blue channel 
+                    red = bayer_green_vert(res, y, i, j) # Compute Red channel
+            # Bayer is red
+            elif z[i,j,0] != 0 :
+                red = y[i,j]    # Red channel
+                blue = bayer_blue_red(res, y, i, j) # Blue channel
+            # Bayer is bleue
+            elif z[i,j,2] != 0 :
+                blue = y[i,j] # Bleu channel
+                red = bayer_blue_red(res, y, i, j) # Res channel
+
+            res[i,j,0] = np.clip(red, 0, 255)
+            res[i,j,2] = np.clip(blue,0,255)
+    return res
+

src/methods/ramanantsitonta_harizo/fonctions.py

+import numpy as np 
+from src.forward_model import CFA
+import cv2 as cv2
+
+def hamilton_adams(y, input_shape):
+    n,p = input_shape[0], input_shape[1]
+    z = np.copy(y)
+    for i in range(2 , n-2): #green interpolation by gradient comparison for every red and blue pixels
+        for j in range(2, p-2):
+            if z[i,j,1] == 0 and z[i,j,0] != 0:   #red pixel
+                #Vertical and horizontal gradient
+                grad_y = np.abs(z[i-1,j,1] - z[i+1,j,1]) + np.abs(2*z[i,j,0] - z[i-2,j,0] - z[i+2,j,0])
+                grad_x = np.abs(z[i,j-1,1] - z[i,j+1,1]) + np.abs(2*z[i,j,0] - z[i,j-2,0] - z[i,j+2,0])
+
+                if grad_x < grad_y:
+                    z[i,j,1] = (z[i,j-1,1] + z[i,j+1,1])/2 + (2*z[i,j,0] - z[i,j-2,0] - z[i,j+2,0])/4
+                elif  grad_x > grad_y:
+                    z[i,j,1] = (z[i-1,j,1] + z[i+1,j,1])/2 + (2*z[i,j,0] - z[i-2,j,0] - z[i+2,j,0])/4
+                else:
+                    z[i,j,1] = (z[i-1,j,1] + z[i+1,j,1] + z[i,j-1,1] + z[i,j+1,1])/4 + (2*z[i,j,0] - z[i,j-2,0] - z[i,j+2,0] + 2*z[i,j,0] - z[i-2,j,0] - z[i+2,j,0])/8
+                
+            elif z[i,j,1] == 0 and z[i,j,2] != 0:    #blue pixel
+                #Vertical and horizontal gradient
+                grad_y = np.abs(z[i-1,j,1] - z[i+1,j,1]) + np.abs(2*z[i,j,2] - z[i-2,j,2] - z[i+2,j,2])
+                grad_x = np.abs(z[i,j-1,1] - z[i,j+1,1]) + np.abs(2*z[i,j,2] - z[i,j-2,2] - z[i,j+2,2])
+
+                if grad_x < grad_y:
+                    z[i,j,1] = (z[i,j-1,1] + z[i,j+1,1])/2 + (2*z[i,j,2] - z[i,j-2,2] - z[i,j+2,2])/4
+                elif  grad_x > grad_y:
+                    z[i,j,1] = (z[i-1,j,1] + z[i+1,j,1])/2 + (2*z[i,j,2] - z[i-2,j,2] - z[i+2,j,2])/4
+                else:
+                    z[i,j,1] = (z[i-1,j,1] + z[i+1,j,1] + z[i,j-1,1] + z[i,j+1,1])/4 + (2*z[i,j,2] - z[i,j-2,2] - z[i,j+2,2] + 2*z[i,j,2] - z[i-2,j,2] - z[i+2,j,2])/8
+        
+    for i in range(1 , n-1): #red/blue interpolation by bilinear interpolation on blue/red pixels
+        for j in range(1, p-1):
+            if z[i,j,2] != 0 :
+                 z[i,j,0] = (z[i-1,j-1,0] + z[i-1,j+1,0] + z[i+1,j-1,0] + z[i+1,j+1,0]) / 4
+
+            elif z[i,j,0] != 0:
+                 z[i,j,2] = (z[i-1,j-1,2] + z[i-1,j+1,2] + z[i+1,j-1,2] + z[i+1,j+1,2]) / 4
+            
+            else:
+                z[i,j] = z[i,j]
+
+    for i in range(1 , n-1): #blue and red interpolation by bilinear interpolation on green pixels
+        for j in range(1, p-1):    
+            if z[i,j,0] == z[i,j,2]:
+                z[i,j,0] = (z[i-1,j,0] + z[i,j-1,0] + z[i+1,j,0] + z[i,j+1,0]) / 4             
+                z[i,j,2] = (z[i-1,j,2] + z[i,j-1,2] + z[i+1,j,2] + z[i,j+1,2]) / 4  
+            
+    return z
+
+# def SSD(y, cfa_img, input_shape):  #SSD algo
+#     n,p = input_shape[0], input_shape[1]
+#     hlist = [16,4,1]
+#     res = np.copy(y)
+#     for h in hlist:
+#         res = NLh(res, cfa_img ,n,p,h)
+#         res = CR(res,n,p)
+    
+#     return res
+
+# def NLh(y, cfa_img, n, p, h):  #NLh part
+#     res = np.copy(cfa_img)
+#     for i in range(4,n-3):
+#         for j in range(4,p-3):
+#             if cfa_img[i,j,0] != 0:
+#                 res[i,j,1] = NLh_calc(y, cfa_img, i, j, 1, h)
+#                 res[i,j,2] = NLh_calc(y, cfa_img, i, j, 2, h)
+#                 print((i,j))
+
+#             elif cfa_img[i,j,1] != 0:
+#                 res[i,j,0] = NLh_calc(y, cfa_img, i, j, 0, h)
+#                 res[i,j,2] = NLh_calc(y, cfa_img, i, j, 2, h)
+#                 print((i,j))
+
+#             else:
+#                 res[i,j,0] = NLh_calc(y, cfa_img, i, j, 0, h)
+#                 res[i,j,1] = NLh_calc(y, cfa_img, i, j, 1, h)
+#                 print((i,j))
+
+#     return res
+
+# def NLh_calc(y, cfa_img, i, j, channel, h):   #aux function to calculate the main sum for each pixel
+#     sum = 0
+#     norm = 0
+#     for k in range(-2,3):
+#         for l in range(-2,3):
+#             if k!=0 and j!=0:
+#                 a = poids(y, channel, i,j,k,l,h)
+#                 sum += a * cfa_img[i+k,j+l,channel]
+#                 norm += a
+
+#     return sum / norm
+
+# def poids(y, channel, i,j,k,l,h):  #aux function to calcultate the weight for a given p=(i,j) , q=(k,l)
+#     som = 0
+#     # for tx in range(-1,2):
+#     #     for ty in range(-1,2):
+
+#     som += (np.abs(y[i-1 , j-1, channel] - y[k+1 , l-1, channel]))**2
+#     som += (np.abs(y[i-1 , j, channel] - y[k+1 , l, channel]))**2
+#     som += (np.abs(y[i-1 , j+1, channel] - y[k+1 , l+1, channel]))**2
+
+#     som += (np.abs(y[i, j-1, channel] - y[k, l-1, channel]))**2
+#     som += (np.abs(y[i, j, channel] - y[k, l, channel]))**2
+#     som += (np.abs(y[i, j+1, channel] - y[k, l+1, channel]))**2   
+
+#     som += (np.abs(y[i+1 , j-1, channel] - y[k+1 , l-1, channel]))**2
+#     som += (np.abs(y[i+1 , j, channel] - y[k+1 , l, channel]))**2
+#     som += (np.abs(y[i+1 , j+1, channel] - y[k+1 , l+1, channel]))**2
+
+
+#     res = np.exp((-1/h**2) * som)
+#     return res
+
+# def CR(img_rgb):  #Chrominance median 
+#     res = cv2.cvtColor(img_rgb, cv2.COLOR_RGB2YUV)
+#     y, u, v = cv2.split(res)
+#     U_med = cv2.medianBlur(u, 3)
+#     V_med = cv2.medianBlur(v, 3)
+
+#     y = cv2.cvtColor(y, cv2.COLOR_GRAY2RGB)
+#     u = cv2.cvtColor(u, cv2.COLOR_GRAY2RGB)
+#     v = cv2.cvtColor(v, cv2.COLOR_GRAY2RGB)
+
+#     res = np.vstack([y, u, v])
+
+#     return res
+
+