src/methods/brice_convers/main_template.py

+import src.methods.brice_convers.dataHandler as DataHandler
+import src.methods.brice_convers.dataEvaluation as DataEvaluation
+import time
+
+WORKING_DIRECOTRY_PATH = "SICOM_Image_Analysis/sicom_image_analysis_project/"
+
+DataHandler = DataHandler.DataHandler(WORKING_DIRECOTRY_PATH)
+DataEvaluation = DataEvaluation.DataEvaluation(DataHandler)
+
+def main(DataHandler):
+
+    IMAGE_PATH = WORKING_DIRECOTRY_PATH + "images/"
+
+    CFA_NAME = "quad_bayer"
+
+    METHOD = "menon"
+
+    startTime = time.time()
+
+    DataHandler.list_images(IMAGE_PATH)
+
+    DataHandler.print_list_images()
+
+    DataHandler.compute_CFA_images(CFA_NAME)
+
+    DataHandler.compute_reconstruction_images(METHOD, {"cfa": CFA_NAME})
+
+    DataHandler.plot_reconstructed_image(0, METHOD, {"cfa": CFA_NAME}, zoomSize="large")
+
+    DataEvaluation.print_metrics(0, METHOD)
+
+    endTime = time.time()
+
+    print("[INFO] Elapsed time: " + str(endTime - startTime) + "s")
+
+    print("[INFO] End")
+
+if __name__ == "__main__":
+    main(DataHandler)

src/methods/brice_convers/menon.py

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

src/methods/brice_convers/reconstruct.py

+"""The main file for the reconstruction.
+This file should NOT be modified except the body of the 'run_reconstruction' function.
+Students can call their functions (declared in others files of src/methods/your_name).
+"""
+
+
+import numpy as np
+import cv2
+from src.forward_model import CFA
+from src.methods.brice_convers.menon import demosaicing_CFA_Bayer_Menon2007
+import src.methods.brice_convers.configuration as configuration
+from src.methods.brice_convers.utilities import quad_bayer_to_bayer
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA("bayer", input_shape)
+
+    if cfa == "quad_bayer":
+        y = quad_bayer_to_bayer(y)
+        op = CFA("bayer", y.shape)
+
+    reconstructed_image = demosaicing_CFA_Bayer_Menon2007(y, op.mask, configuration.PIXEL_PATTERN, configuration.REFINING_STEP)
+
+    if cfa == "quad_bayer":
+        return cv2.resize(reconstructed_image, input_shape[:2], interpolation=cv2.INTER_CUBIC)
+    else:
+        return reconstructed_image
+
+####
+####
+####
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/brice_convers/utilities.py

+import os
+import cv2
+
+def folderExists(path):
+	CHECK_FOLDER = os.path.isdir(path)
+
+	# If folder doesn't exist, it creates it.
+	if not CHECK_FOLDER:
+		os.makedirs(path)
+		print("[DATA] You created a new folder : " +  str(path))
+
+import numpy as np
+
+def quad_bayer_to_bayer(quad_bayer_pattern):
+    # We test that thee quad bayer size fit with a multiple of 2 for width and height). If not, we pad it.
+	if quad_bayer_pattern.shape[0] % 2 != 0 or quad_bayer_pattern.shape[1] % 2 != 0:
+		print("[INFO] The quad bayer pattern size is not valid. We need to pad it.")
+
+		pad_schema = []
+
+		if quad_bayer_pattern.shape[0] % 2 != 0:
+				pad_schema.append([0, 1])
+
+		if quad_bayer_pattern.shape[1] % 2 != 0:
+				pad_schema.append([0, 1])
+		else:
+				pad_schema.append([0, 0])
+
+		quad_bayer_pattern = np.pad(quad_bayer_pattern, pad_schema, mode="reflect")[:, 2:]
+
+	# we create a new bayer pattern with the good size
+	bayer_pattern = np.zeros((quad_bayer_pattern.shape[0] // 2, quad_bayer_pattern.shape[1] // 2))
+
+	# We combine adjacent pixels to create the Bayer pattern
+	for i in range(0, quad_bayer_pattern.shape[0], 2):
+		for j in range(0, quad_bayer_pattern.shape[1], 2):
+			bayer_pattern[i // 2, j // 2] = (
+			quad_bayer_pattern[i, j] +
+			quad_bayer_pattern[i, j + 1] +
+			quad_bayer_pattern[i + 1, j] +
+			quad_bayer_pattern[i + 1, j + 1]
+		) / 4
+
+	# We resize bayer iamge to the original image size
+
+	#return cv2.resize(bayer_pattern, quad_bayer_pattern.shape, interpolation=cv2.INTER_CUBIC)
+	return bayer_pattern

src/methods/chaari_mohamed/fonctions.py

+import numpy as np
+from scipy.signal import convolve2d
+from src.forward_model import CFA
+
+def bilinear_demosaicing(op: CFA, y: np.ndarray) -> np.ndarray:
+    """
+    Bilinear demosaicing method.
+
+    Args:
+        op (CFA): CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    # Copie des valeurs directement connues pour chaque canal
+    red = y[:, :, 0]
+    green = y[:, :, 1]
+    blue = y[:, :, 2]
+
+    # Création des masques pour chaque couleur selon le motif CFA
+    mask_red = (op.mask == 0)  # Supposons que 0 correspond au rouge dans le masque
+    mask_green = (op.mask == 1)  # Supposons que 1 correspond au vert
+    mask_blue = (op.mask == 2)  # Supposons que 2 correspond au bleu
+
+    # Interpolation bilinéaire pour le rouge et le bleu
+    # Note: np.multiply multiplie les éléments correspondants des tableaux, c'est pourquoi nous utilisons np.multiply au lieu de *
+    red_interp = convolve2d(np.multiply(red, mask_red), [[1/4, 1/2, 1/4], [1/2, 1, 1/2], [1/4, 1/2, 1/4]], mode='same')
+    blue_interp = convolve2d(np.multiply(blue, mask_blue), [[1/4, 1/2, 1/4], [1/2, 1, 1/2], [1/4, 1/2, 1/4]], mode='same')
+
+    # Interpolation bilinéaire pour le vert
+    # Pour le vert, nous utilisons un autre noyau car il y a plus de pixels verts
+    green_interp = convolve2d(np.multiply(green, mask_green), [[0, 1/4, 0], [1/4, 1, 1/4], [0, 1/4, 0]], mode='same')
+
+    # Création de l'image interpolée
+    demosaicked_image = np.stack((red_interp, green_interp, blue_interp), axis=-1)
+
+    # Correction des valeurs interpolées: on réapplique les valeurs connues pour éviter le flou
+    demosaicked_image[:, :, 0][mask_red] = red[mask_red]
+    demosaicked_image[:, :, 1][mask_green] = green[mask_green]
+    demosaicked_image[:, :, 2][mask_blue] = blue[mask_blue]
+
+    # Clip pour s'assurer que toutes les valeurs sont dans la plage [0, 1]
+    demosaicked_image = np.clip(demosaicked_image, 0, 1)
+
+    return demosaicked_image
+
+
+def quad_bayer_demosaicing(op: CFA, y: np.ndarray) -> np.ndarray:
+    """
+    Demosaicing method for Quad Bayer CFA pattern.
+
+    Args:
+        op (CFA): CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    
+    # Interpolation bilinéaire pour chaque canal
+    red_interp = convolve2d(np.multiply(y[:, :, 0], op.mask == 0), [[1/4, 1/2, 1/4], [1/2, 1, 1/2], [1/4, 1/2, 1/4]], mode='same')
+    green_interp = convolve2d(np.multiply(y[:, :, 1], op.mask == 1), [[0, 1/4, 0], [1/4, 1, 1/4], [0, 1/4, 0]], mode='same')
+    blue_interp = convolve2d(np.multiply(y[:, :, 2], op.mask == 2), [[1/4, 1/2, 1/4], [1/2, 1, 1/2], [1/4, 1/2, 1/4]], mode='same')
+
+    # Assemblage de l'image interpolée
+    demosaicked_image = np.stack((red_interp, green_interp, blue_interp), axis=-1)
+
+    # Réapplication des valeurs connues
+    demosaicked_image[:, :, 0][op.mask == 0] = y[:, :, 0][op.mask == 0]
+    demosaicked_image[:, :, 1][op.mask == 1] = y[:, :, 1][op.mask == 1]
+    demosaicked_image[:, :, 2][op.mask == 2] = y[:, :, 2][op.mask == 2]
+
+    # Clip des valeurs pour les maintenir dans la plage appropriée
+    demosaicked_image = np.clip(demosaicked_image, 0, 1)
+
+    return demosaicked_image

src/methods/chaari_mohamed/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.chaari_mohamed.fonctions import bilinear_demosaicing,quad_bayer_demosaicing
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """
+    Performs demosaicking on y based on the CFA pattern.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA pattern. 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)
+    
+    if cfa == 'bayer':
+        res = bilinear_demosaicing(op, y)
+    elif cfa == 'quad_bayer':
+        res = quad_bayer_demosaicing(op, y)
+    else:
+        raise ValueError("Unsupported CFA pattern. Supported patterns are 'bayer' and 'quad_bayer'.")
+    
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/charpentier_laurine/functions.py

+import numpy as np
+import cv2
+from scipy.signal import convolve2d
+
+from src.forward_model import CFA
+
+
+def superpixel(op: CFA, y: np.ndarray) -> np.ndarray:
+    """Performs the method of variable number of gradients
+
+    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[0]//2, op.input_shape[1]//2, op.input_shape[2]))
+
+        for i in range(1,op.input_shape[0],2):
+            for j in range(1, op.input_shape[1],2):
+                res[i//2,j//2,0] = z[ i-1,j ,0]
+                res[i//2,j//2,1] = (z[i,j,1] + z[i-1,j-1,1]) / 2
+                res[i//2,j//2,2] = z[ i,j-1 ,2]
+
+    else:
+        res = np.empty((op.input_shape[0]//2, op.input_shape[1]//2, op.input_shape[2]))
+
+        print()
+
+        for i in range(2,op.input_shape[0]-5,4):
+            for j in range(2, op.input_shape[1]-5,4):
+
+                res[i//2,j//2,0] = z[ i-2,j ,0]
+                res[i//2,j//2,1] = (z[i,j,1] + z[i-2,j-2,1]) / 2
+                res[i//2,j//2,2] = z[ i,j-2 ,2]
+
+                res[i//2+1,j//2,0] = z[ i-2+1,j ,0]
+                res[i//2+1,j//2,1] = (z[i+1,j,1] + z[i-2+1,j-2,1]) / 2
+                res[i//2+1,j//2,2] = z[ i+1,j-2 ,2]
+
+                res[i//2+1,j//2+1,0] = z[ i-2+1,j+1 ,0]
+                res[i//2+1,j//2+1,1] = (z[i+1,j+1,1] + z[i-2+1,j-2+1,1]) / 2
+                res[i//2+1,j//2+1,2] = z[ i+1,j-2+1 ,2]
+
+                res[i//2,j//2+1,0] = z[ i-2,j+1 ,0]
+                res[i//2,j//2+1,1] = (z[i,j+1,1] + z[i-2,j-2+1,1]) / 2
+                res[i//2,j//2+1,2] = z[ i,j-2+1 ,2]
+
+    return cv2.resize(res, (op.input_shape[0], op.input_shape[1]))
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/charpentier_laurine/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.charpentier_laurine.functions import superpixel
+
+
+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 = superpixel(op, y)
+
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/dai/function.py

+import numpy as np
+from scipy.signal import convolve2d
+
+from src.forward_model import CFA
+
+
+def naive_interpolation(op: CFA, y: np.ndarray) -> np.ndarray:
+    """Performs a simple interpolation of the lost pixels.
+
+    Args:
+        op (CFA): CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    
+
+    res = np.empty(op.input_shape)
+    z = op.adjoint(y)
+
+    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')
+    
+    return res
+
+
+
+
+def Spectral_difference(op,y):
+
+    z = op.adjoint(y)
+    y_hat  = naive_interpolation(op,y)
+    res = np.empty(op.input_shape)
+    for l in range(3):
+        res[:,:,0]+=z[:,:,l]+convolve2d(z[:,:,0]-y_hat[:,:,l]*op.mask[:,:,0], ker_bayer_red_blue, mode='same')*op.mask[:,:,l]
+        res[:,:,1]+=z[:,:,l]+convolve2d(z[:,:,1]-y_hat[:,:,l]*op.mask[:,:,1], ker_bayer_green, mode='same')*op.mask[:,:,l]
+        res[:,:,2]+=z[:,:,l]+convolve2d(z[:,:,2]-y_hat[:,:,l]*op.mask[:,:,2], ker_bayer_red_blue, mode='same')*op.mask[:,:,l]
+
+    return res
+
+def quad_bayer_to_bayer(y_quad):
+    y_bayer=np.copy(y_quad)
+
+    for i in range(1,y_quad.shape[0],4):
+        temp = np.copy(y_bayer[i,:])
+        y_bayer[i,:]=y_bayer[i+1,:]
+        y_bayer[i+1,:]=temp
+    
+    for j in range(1,y_quad.shape[1],4):
+        temp = np.copy(y_bayer[:,j])
+        y_bayer[:,j]=y_bayer[:,j+1]
+        y_bayer[:,j+1]=temp
+
+
+
+    return y_bayer
+
+
+
+
+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
\ No newline at end of file

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

src/methods/david_alexis/functions.py

+
+import numpy as np
+
+
+# interpolation kernels
+w1 = 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 # G at R location
+w2 = w1 #G at B location
+w3 = 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 # R at G location (Brow, Rcol)
+w4 = 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 # R at G location (Rrow, Bcol)
+w5 = np.array([[0,0,-3/2,0,0], [0,2,0,2,0], [-3/2,0,6,0,-1], [0,2,0,2,0], [0,0,-3/2,0,0]])/8 # R at B location (Brow, Bcol)
+w6 = w3 #R at G location (Brow, Rcol)
+w7 = w4 #B at G location (Rrow, Bcol)
+w8 = w5 #B at R location (Rrow, Rcol)
+w_s = [w1, w2, w3, w4, w5, w6, w7, w8]
+
+def conv(A,B):
+    return np.sum(A*B)
+
+def bayer_gradient_interpolation(y, op, w=w_s):
+    """
+    Second interpolation method: convolution with a 2D kernels (5x5) using multi-channel interpolations
+    y: input image
+    w: interpolation kernels
+    """
+    y_padded = np.pad(y, ((2,2),(2,2)), 'constant', constant_values=0)
+    size_padded = y_padded.shape
+
+    #separate channels
+    y_0 = y*op.mask[:, :, 0]
+    y_1 = y*op.mask[:, :, 1]
+    y_2 = y*op.mask[:, :, 2]
+
+    for i in range(2,size_padded[0]-2):
+        for j in range(2,size_padded[1]-2):
+
+            # determine the location of the pixel
+            if op.mask[i-2,j-2,1] == 1: # at G
+                if op.mask[i-2,j-3,0] == 1: #R row
+                    y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w3)
+                    y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w7)
+                else: #B row
+                    y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w4)
+                    y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w6)
+
+            elif op.mask[i-2,j-2,0] == 1: # at R
+                y_1[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w1)
+                y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w8)
+
+            else: # at B
+                y_1[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w2)
+                y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w5)
+
+    y_res = np.zeros((y.shape[0], y.shape[1], 3))
+    y_res[:, :, 0] = y_0
+    y_res[:, :, 1] = y_1
+    y_res[:, :, 2] = y_2
+
+    return y_res
+
+
+
+import numpy as np
+def conv(A,B):
+    return np.sum(A*B)
+def quad_gradient_interpolation(y, op):
+    """
+    Second interpolation method: convolution with a 2D kernels (5x5) using multi-channel interpolations
+    y: input image
+    w: interpolation kernels
+    """
+    y_padded = np.pad(y, ((2,2),(2,2)), 'constant', constant_values=0)
+    size_padded = y_padded.shape
+
+    mask_padded = np.pad(op.mask, ((2,2),(2,2),(0,0)), 'constant', constant_values=-1)
+
+    #separate channels
+    y_0 = y*op.mask[:, :, 0]
+    y_1 = y*op.mask[:, :, 1]
+    y_2 = y*op.mask[:, :, 2]
+
+    #interpolation kernels
+    # green quad kernels
+    wg11 = np.zeros((5,5))
+    wg11[0,2] = -3
+    wg11[1,2] = 13
+    wg11[2,4] = 2
+    wg11 += wg11.T
+    wg11 = wg11/24
+
+    # rotate the kernel
+    wg12 = np.rot90(wg11,3)
+    wg22 = np.rot90(wg11,2)
+    wg21 = np.rot90(wg11,1)
+
+    #red/blue quad kernels
+        # center part 
+    w_c11 = np.zeros((5,5))
+    w_c11[0:2,0:2] = np.array([[-1,-1],[-1,9]])/6
+    w_c21 = np.rot90(w_c11,1)
+    w_c22 = np.rot90(w_c11,2)
+    w_c12 = np.rot90(w_c11,3)
+        # left part
+    w_l11 = np.zeros((5,5))
+    w_l11[0:2, 2] = np.array([-3,13])
+    w_l11[4,2] = 2
+    w_l11 = w_l11/12
+    w_l22 = np.rot90(w_l11,2)
+        # top part
+    w_t11 = np.rot90(w_l11,1)
+    w_t22 = np.rot90(w_t11,2)
+
+
+    for i in range(2,size_padded[0]-2):
+        for j in range(2,size_padded[1]-2):
+
+            # determine the location of the pixel
+            if mask_padded[i,j,1] == 1: # at G
+                if mask_padded[i-1,j,0] + mask_padded[i+1,j,0] == 1: # at R column
+                    if mask_padded[i,j-1,2] == 1:
+                        y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_t11)
+                    elif mask_padded[i,j+1,2] == 1:
+                        y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_t22)
+                    if mask_padded[i-1,j,0] == 1:
+                        y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_l11)
+                    elif mask_padded[i+1,j,0] == 1:
+                        y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_l22)
+
+                else:
+                    if mask_padded[i,j-1,0] == 1:
+                        y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_t11)
+                    elif mask_padded[i,j+1,0] == 1:
+                        y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_t22)
+                    if mask_padded[i-1,j,2] == 1:
+                        y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_l11)
+                    elif mask_padded[i+1,j,2] == 1:
+                        y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_l22)
+
+
+            elif mask_padded[i,j,1] == 0: # at R or B
+                if mask_padded[i-1,j,1] ==1:
+                    if mask_padded[i,j-1,1] ==1:
+                        y_1[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],wg11)
+                    else:
+                        y_1[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],wg12)
+                else:
+                    if mask_padded[i,j-1,1] ==1:
+                        y_1[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],wg21)
+                    else:
+                        y_1[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],wg22)
+                
+                if mask_padded[i,j,0] == 1: #at R
+                    if mask_padded[i-1,j-1,2] ==1: #center (1,1)
+                        y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_c11)
+                    elif mask_padded[i-1,j+1,2] ==1: #center (1,2)
+                        y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_c12)
+                    elif mask_padded[i+1,j-1,2] ==1: #center (2,1)
+                        y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_c21)
+                    elif mask_padded[i+1,j+1,2] ==1: #center (2,2)
+                        y_2[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_c22)
+                    
+                elif mask_padded[i,j,2] == 1: #at B
+                    if mask_padded[i-1,j-1,0] ==1: #center (1,1)
+                        y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_c11)   
+                    elif mask_padded[i-1,j+1,0] ==1: #center (1,2)
+                        y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_c12)
+                    elif mask_padded[i+1,j-1,0] ==1: #center (2,1)
+                        y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_c21)
+                    elif mask_padded[i+1,j+1,0] ==1: #center (2,2)
+                        y_0[i-2,j-2] = conv(y_padded[i-2:i+3, j-2:j+3],w_c22)
+
+    y_res = np.zeros((y.shape[0], y.shape[1], 3))
+    y_res[:, :, 0] = y_0
+    y_res[:, :, 1] = y_1
+    y_res[:, :, 2] = y_2
+
+    return y_res
+
+from scipy import signal
+def interpolate(y, op, mode = 'bayer'):
+    """
+    First interpolation method: convolution with a 2D kernel
+    y: input image
+    op: object of class CFA
+    """
+    #separate channels
+    y_0 = y*op.mask[:, :, 0]
+    y_1 = y*op.mask[:, :, 1]
+    y_2 = y*op.mask[:, :, 2]
+
+    #convolution 2D
+    
+    if mode == 'bayer':
+        a = np.array([1,2,1])
+    elif mode == 'quad_bayer':
+        a = np.array([1,2,3,2,1])
+    else:
+        raise Exception('Mode not supported')
+    
+    # create the 2D kernel
+    w = a[:, None] * a[None, :]
+    w = w / np.sum(w)
+
+    #interpolate green
+    y_1_interpolated = signal.convolve2d(y_1, w, mode='same', boundary='fill', fillvalue=0)
+
+    #interpolate red
+    y_0_interpolated = signal.convolve2d(y_0, w, mode='same', boundary='fill', fillvalue=0)
+
+    # interpolate blue
+    y_2_interpolated = signal.convolve2d(y_2, w, mode='same', boundary='fill', fillvalue=0)
+
+    y_res = np.zeros((y.shape[0], y.shape[1], 3))   
+    y_res[:,:,0] = y_0_interpolated*4
+    y_res[:,:,1] = y_1_interpolated*2
+    y_res[:,:,2] = y_2_interpolated*4
+
+    return y_res

src/methods/david_alexis/reconstruct.py

+"""The main file for the reconstruction.
+"""
+
+import numpy as np
+
+from src.forward_model import CFA
+from src.methods.david_alexis.functions import bayer_gradient_interpolation, quad_gradient_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.
+    """
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+
+    if cfa == "bayer": 
+        res = bayer_gradient_interpolation(y, op)
+    else:
+        res = quad_gradient_interpolation(y, op)
+
+    return res
+
+
+
+# Author: Alexis DAVID

src/methods/digeronimo/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.digeronimo.somefunc import spectral_difference, normalization, quad_bayer_to_bayer_pattern, quad_bayer_to_bayer_image, bilinear_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.
+    
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa, input_shape)
+
+    if op.cfa != 'bayer' and op.cfa != 'quad_bayer':
+        raise ValueError("CFA must be bayer or quad_bayer")
+    
+    if op.cfa == 'quad_bayer':
+        quad_bayer_to_bayer_pattern(op)
+        quad_bayer_to_bayer_image(y)
+
+    # Apply adjoint operation on raw aquisition
+    z = op.adjoint(y)
+
+    # Demosaicing process
+    res_bi = bilinear_interpolation(op, z)
+    res_sd = spectral_difference(op, z, res_bi)
+    res = normalization(res_sd) # min-max normalization
+
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/digeronimo/somefunc.py

+import numpy as np
+from scipy.signal import convolve2d
+
+from src.forward_model import CFA
+
+def bilinear_interpolation(op: CFA, z: np.ndarray) -> np.ndarray:
+    """Perform bilinear interpolation for demosaicing
+
+    Args:
+        op (CFA): CFA operator.
+        z (np.ndarray): Adjoint image.
+
+    Returns:
+        np.ndarray: Interpolated image.
+    """
+    # Bi-linear 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
+
+    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')
+
+    return res
+
+def spectral_difference(op: CFA, z: np.ndarray, res: np.ndarray) -> np.ndarray:
+    """Perform spectral difference method for demosaicing
+
+    Args:
+        op (CFA): CFA operator.
+        z (np.ndarray): Adjoint image.
+        res (np.ndarray): Interpolated image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+
+    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
+
+    # Computation of spectral differences
+
+    delta_red_green_quad =  z[:, :, 0] - np.multiply(res[:, :, 1],op.mask[:,:,0])
+    delta_red_blue_quad =  z[:, :, 0] - np.multiply(res[:, :, 2],op.mask[:,:,0])
+
+    delta_blue_green_quad =  z[:, :, 2] - np.multiply(res[:, :, 1],op.mask[:,:,2])
+    delta_blue_red_quad =  z[:, :, 2] - np.multiply(res[:, :, 0],op.mask[:,:,2])
+
+    delta_green_red_quad =  z[:, :, 1] - np.multiply(res[:, :, 0],op.mask[:,:,1])
+    delta_green_blue_quad =  z[:, :, 1] - np.multiply(res[:, :, 2],op.mask[:,:,1])
+
+    # Estimation
+    res_sd = np.empty(op.input_shape)
+
+    res_sd[:,:,0] = res[:, :, 1] + convolve2d(delta_red_green_quad,ker_bayer_red_blue,mode='same') + res[:, :, 2] + convolve2d(delta_red_blue_quad,ker_bayer_red_blue,mode='same') 
+    res_sd[:,:,2] = res[:, :, 1] + convolve2d(delta_blue_green_quad,ker_bayer_red_blue,mode='same') + res[:, :, 0] + convolve2d(delta_blue_red_quad,ker_bayer_red_blue,mode='same') 
+    res_sd[:,:,1] = res[:, :, 0] + convolve2d(delta_green_red_quad,ker_bayer_green,mode='same') + res[:, :, 2] + convolve2d(delta_green_blue_quad,ker_bayer_green,mode='same') 
+
+    return res_sd
+
+def normalization(res_sd:np.ndarray)-> np.ndarray:
+    """Perform a min-max normalization
+
+    Args:
+        res_sd (np.ndarray): Demosaicked image.
+
+    Returns:
+        np.ndarray: Normalized image.
+    """
+
+    res = (res_sd - np.min(res_sd)) / (np.max(res_sd) - np.min(res_sd))
+    
+    return res
+
+def quad_bayer_to_bayer_pattern(op: CFA):
+    """Quad to Bayer pattern conversion by swapping method
+
+    Args:
+        op (CFA): CFA operator.
+
+    Returns:
+        CFA: Bayer pettern.
+    """
+    if op.cfa == 'quad_bayer':
+        for j in range(1, op.mask.shape[1], 4):
+            op.mask[:,j], op.mask[:,j+1] = op.mask[:,j+1].copy(), op.mask[:,j].copy()
+
+        for i in range(1, op.mask.shape[0], 4):
+            op.mask[i, :], op.mask[i+1,:] = op.mask[i+1,:].copy(), op.mask[i,:].copy()
+            
+        for i in range(1, op.mask.shape[0], 4):
+            for j in range(1, op.mask.shape[1], 4):
+                op.mask[i,j], op.mask[i+1,j+1] = op.mask[i+1,j+1].copy(), op.mask[i,j].copy() 
+    return 0
+
+def quad_bayer_to_bayer_image(y:np.array):
+    """Quad to Bayer conversion of a mosaicked image by swapping method
+
+    Args:
+        y (np.array): Mosaicked image.
+    """
+        
+    for j in range(1, y.shape[1], 4):
+        y[:,j], y[:,j+1] = y[:,j+1].copy(), y[:,j].copy()
+
+    for i in range(1, y.shape[0], 4):
+        y[i, :], y[i+1,:] = y[i+1,:].copy(), y[i,:].copy()
+        
+    for i in range(1, y.shape[0], 4):
+        for j in range(1, y.shape[1], 4):
+            y[i,j], y[i+1,j+1] = y[i+1,j+1].copy(), y[i,j].copy() 
+            
+    return 0
\ No newline at end of file