Merge branch 'withoutMain' into 'master'

Image Analysis Project

See merge request !35

src/methods/brice_convers/README.md

+# Image Analysis Project
+
+Numerous RGB cameras in the commercial sector employ Color Filter Array (CFA) technology.
+This technology involves an array of red, green, or blue filters placed atop the sensors, usually
+organized in periodic patterns. The incident light is consequently filtered by each filter, before
+being captured by the sensor. The typical process of acquiring images utilizes a predefined CFA
+pattern to allocate a color to each pixel of the sensor. The goal of this code project is to demosaicing these CFA Image.
+
+You can check the report called **Image_Analysis_Project_Report_Brice_Convers** to find more information about it.
+
+## How to use
+
+To use the code project you can move the **main_template.py** outside the **src** file and execute it.
+Or you can also call the **run_reconstruction** method in you programme.
+
+```Python
+
+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})
+
+    #The first agurment (ex: 3) is the image index in the list print by "DataHandler.print_list_images()"
+    DataHandler.plot_reconstructed_image(3, METHOD, {"cfa": CFA_NAME}, zoomSize="large")
+
+    DataEvaluation.print_metrics(3, METHOD)
+
+    endTime = time.time()
+
+    print("[INFO] Elapsed time: " + str(endTime - startTime) + "s")
+
+    print("[INFO] End")
+
+if __name__ == "__main__":
+    main(DataHandler)
+
+```
+
+## TODO List:
+
+- Fix menon method pour quad bayer pattern with landscape picture
+
+
+
+
+
+
+
+## References:
+
+[1] [*Research Paper:*](https://ieeexplore.ieee.org/document/4032820) Used for Menon Method.
+
+## Authors
+- [Brice Convers](https://briceconvers.com)

src/methods/brice_convers/configuration.py

+
+PIXEL_PATTERN = "RGB"
+REFINING_STEP = True

src/methods/brice_convers/dataEvaluation.py

+from src.methods.brice_convers.dataHandler import DataHandler
+from src.utils import psnr, ssim
+from sklearn.metrics import f1_score, mean_squared_error
+import numpy as np
+
+class DataEvaluation:
+    def __init__(self, DataHandler: DataHandler):
+        DataEvaluation.DataHandler = DataHandler
+
+    def print_metrics(self, indexImage, method):
+        DataEvaluation.DataHandler.indexImageExists(indexImage)
+
+        img = DataEvaluation.DataHandler.load_image(indexImage)
+        res = DataEvaluation.DataHandler.get_reconstructed_image(indexImage, method)
+
+        ssimMetric = ssim(img, res)
+        psnrMetrc = psnr(img, res)
+        mse = mean_squared_error(img.flatten(), res.flatten())
+        mseRedPixels = mean_squared_error(img[:,:,0], res[:,:,0])
+        mseGreenPixels = mean_squared_error(img[:,:,1], res[:,:,1])
+        mseBluePixels = mean_squared_error(img[:,:,2], res[:,:,2])
+        miMetric = DataEvaluation.MI(img, res)
+        ccMetric = DataEvaluation.CC(img, res)
+        sadMetric = DataEvaluation.SAD(img, res)
+        lsMetric = DataEvaluation.LS(img, res)
+
+        print("[INFO] Metrics for image {}".format(indexImage))
+        print("#" * 30)
+        print(" SSIM: {:.6}  ".format(ssimMetric))
+        print(" PSNR: {:.6}   ".format(psnrMetrc))
+        print(" MSE : {:.3e}    ".format(mse))
+        print(" MSE (R): {:.3e}    ".format(mseRedPixels))
+        print(" MSE (G): {:.3e}    ".format(mseGreenPixels))
+        print(" MSE (B): {:.3e}    ".format(mseBluePixels))
+        print(" MI: {:.6}     ".format(miMetric))
+        print(" CC: {:.6}    ".format(ccMetric))
+        print(" SAD: {:.6} ".format(sadMetric))
+        print(" LS: {:.3e} ".format(lsMetric))
+        print("#" * 30)
+
+    #Mutual Information
+    def MI(img_mov, img_ref):
+        hgram, x_edges, y_edges = np.histogram2d(img_mov.ravel(), img_ref.ravel(), bins=20)
+        pxy = hgram / float(np.sum(hgram))
+        px = np.sum(pxy, axis=1) # marginal for x over y
+        py = np.sum(pxy, axis=0) # marginal for y over x
+        px_py = px[:, None] * py[None, :] # Broadcast to multiply marginals
+        # Now we can do the calculation using the pxy, px_py 2D arrays
+        nzs = pxy > 0 # Only non-zero pxy values contribute to the sum
+        return np.sum(pxy[nzs] * np.log(pxy[nzs] / px_py[nzs]))
+
+    # Cross Correlation
+    def CC(img_mov, img_ref):
+        # Vectorized versions of c,d,e
+        a = img_mov.astype('float64')
+        b = img_ref.astype('float64')
+
+        # Calculating mean values
+        AM = np.mean(a)
+        BM = np.mean(b)
+
+        c_vect = (a - AM) * (b - BM)
+        d_vect = (a - AM) ** 2
+        e_vect = (b - BM) ** 2
+
+        # Finally get r using those vectorized versions
+        r_out = np.sum(c_vect) / float(np.sqrt(np.sum(d_vect) * np.sum(e_vect)))
+        return r_out
+
+    #Sum of Absolute Differences
+    def SAD(img_mov, img_ref):
+        img1 = img_mov.astype('float64')
+        img2 = img_ref.astype('float64')
+        ab = np.abs(img1 - img2)
+        sav = np.sum(ab.ravel())
+        sav /= ab.ravel().shape[0]
+        return sav
+
+    #Sum of Least Squared Errors
+    def LS(img_mov, img_ref):
+        img1 = img_mov.astype('float64')
+        img2 = img_ref.astype('float64')
+        r = (img1 - img2)**2
+        sse = np.sum(r.ravel())
+        sse /= r.ravel().shape[0]
+        return sse
+

src/methods/brice_convers/dataHandler.py

+import os
+import matplotlib.pyplot as plt
+import matplotlib.image as mpimg
+from src.forward_model import CFA
+from pathlib import Path
+from src.methods.baseline.reconstruct import run_reconstruction
+from src.methods.brice_convers.reconstruct import run_reconstruction as run_reconstruction_brice_convers
+from src.methods.brice_convers.utilities import folderExists
+from src.utils import normalise_image, save_image, psnr, ssim
+from skimage.io import imread
+import numpy as np
+
+class DataHandler:
+    def __init__(self, workingDirectoryPath = ""):
+        DataHandler.imagePaths =[]
+        DataHandler.imagePathsLabels = []
+        DataHandler.CFA_images = {}
+        DataHandler.reconstruction_images = {"interpolation": {}, "menon": {}}
+
+        DataHandler.IMAGE_TYPES = (".jpg", ".jpeg", ".png", ".bmp", ".tif", ".tiff")
+        DataHandler.WD_PATH = workingDirectoryPath
+
+    def list_files(self, basePath, interval, validExts=None, contains=None):
+        sum = 0
+
+        if interval is not None:
+            counter = interval[0]
+            lastFile = interval[1]
+        else:
+            counter = 0
+            lastFile = -1
+
+        # loop over the directory structure
+        for (rootDir, dirNames, filenames) in os.walk(basePath):
+            # loop over the filenames in the current directory
+            for filename in filenames:
+                # if the contains string is not none and the filename does not contain
+                # the supplied string, then ignore the file
+                if contains is not None and filename.find(contains) == -1:
+                    continue
+
+                # determine the file extension of the current file
+                ext = filename[filename.rfind("."):].lower()
+
+                # check to see if the file is an image and should be processed
+                if validExts is None or ext.endswith(validExts):
+                    if sum >= counter and (sum <= lastFile) or lastFile == -1:
+                        # construct the path to the image and yield it
+                        imagePath = os.path.join(rootDir, filename)
+                        imageName = Path(imagePath).stem
+                        DataHandler.imagePaths.append(imagePath)
+                        DataHandler.imagePathsLabels.append(imageName)
+
+                    sum += 1
+
+    def list_images(self, basePath, interval = None, contains=None):
+        # return the set of files that are valid
+        DataHandler.list_files(self, basePath, interval, validExts=DataHandler.IMAGE_TYPES, contains=contains)
+        print("[INFO] There are {} images.".format(len(DataHandler.imagePaths)))
+
+    def print_list_images(self):
+        print("[INFO] This is the order or your image in the list. IndexImage is the index in this table:")
+        print(DataHandler.imagePathsLabels)
+
+    def plot_input_images(self, rows=3, cols=3, figsize=(15, 5), max_images = 6, title="Input Images"):
+        fig = plt.figure(figsize=figsize)
+        fig.suptitle(title, fontsize=16)
+
+        for i, imagePath in enumerate(DataHandler.imagePaths):
+            if i >= max_images:
+                break
+
+            img = mpimg.imread(imagePath)
+            fig.add_subplot(rows, cols, i+1)
+            # image title
+            plt.title(imagePath.split(os.path.sep)[-1])
+            plt.text(0, -0.1, f"Shape of the image: {img.shape}.", fontsize=12, transform=plt.gca().transAxes)
+            plt.imshow(img)
+
+        plt.show()
+
+    def plot_raw_transformation(self, zoom = True, specificIndex = None, rows=8, cols=3, figsize=(15, 5), max_images = 24, title="Raw Transformation"):
+        fig, axs = plt.subplots(rows, cols, figsize=figsize)
+        fig.suptitle(title, fontsize=16)
+
+        for i, imagePath in enumerate(DataHandler.imagePaths):
+            if i >= max_images and max_images != -1:
+                break
+
+            if specificIndex is not None:
+                if i != specificIndex:
+                    continue
+
+            nameImage = Path(imagePath).stem
+
+            if DataHandler.CFA_images.get(nameImage) is None:
+                print("[ERROR] There is no CFA image for the image {}.".format(nameImage))
+                continue
+
+            img = mpimg.imread(imagePath)
+
+            y = DataHandler.CFA_images[nameImage].direct(img)
+            z = DataHandler.CFA_images[nameImage].adjoint(y)
+
+            if specificIndex is None:
+                line = 2*i
+                subLine = 2*i+1
+            else:
+                line = 0
+                subLine = 1
+
+            axs[line, 0].imshow(img)
+            axs[line, 0].set_title('Input image')
+            axs[line, 1].imshow(y, cmap='gray')
+            axs[line, 1].set_title('Output image')
+            axs[line, 2].imshow(z)
+            axs[line, 2].set_title('Adjoint image')
+
+            if zoom:
+                axs[subLine, 0].imshow(img[800:864, 450:514])
+                axs[subLine, 0].set_title('Zoomed input image')
+                axs[subLine, 1].imshow(y[800:864, 450:514], cmap='gray')
+                axs[subLine, 1].set_title('Zoomed output image')
+                axs[subLine, 2].imshow(z[800:864, 450:514])
+                axs[subLine, 2].set_title('Zoomed adjoint image')
+
+        print("[INFO] There are {} ploted images.".format(max_images))
+
+    def plot_specific_raw_transformation(self, indexImage, zoom = True, rows=2, cols=3, figsize=(15, 5), title="Raw Transformation"):
+        DataHandler.plot_raw_transformation(self, zoom, indexImage, rows, cols, figsize, -1, title)
+
+    def compute_CFA_images(self, CFA_NAME):
+        for i, imagePath in enumerate(DataHandler.imagePaths):
+            nameImage = Path(imagePath).stem
+            DataHandler.compute_CFA_image(self, CFA_NAME, nameImage, i)
+
+        print("[INFO] There are {} CFA images.".format(len(DataHandler.CFA_images)))
+
+    def compute_reconstruction_image(self, method, indexImage, options = None):
+        if len(DataHandler.imagePaths) == 0:
+            print("[ERROR] There is no image in imagePaths")
+            return
+
+        # Test key method
+        if method not in DataHandler.reconstruction_images.keys():
+            print("[ERROR] The method {} is not valid.".format(method))
+            exit(1)
+
+        nameImage = Path(DataHandler.imagePaths[indexImage]).stem
+
+        if DataHandler.CFA_images.get(nameImage) is None:
+            print("[ERROR] There is no CFA image for the image {}.".format(nameImage))
+            exit(1)
+
+        img = DataHandler.load_image(self, indexImage)
+
+        img_CFA = DataHandler.CFA_images[nameImage].direct(img)
+
+        cfa = options.get("cfa")
+
+        if cfa is None:
+            print("[ERROR] You must specify the cfa.")
+            exit(1)
+
+        if method == "interpolation":
+
+            DataHandler.reconstruction_images[method].setdefault(nameImage, run_reconstruction(img_CFA, cfa))
+
+        if method == "menon":
+
+            DataHandler.reconstruction_images[method].setdefault(nameImage, run_reconstruction_brice_convers(img_CFA, cfa))
+
+    def compute_reconstruction_images(self, method, options = None):
+
+        for i in range(len(DataHandler.imagePaths)):
+            DataHandler.compute_reconstruction_image(self, method, i, options)
+
+        print("[INFO] There are {} images which have been reconstructed.".format(len(DataHandler.reconstruction_images[method])))
+
+    def plot_reconstructed_image(self, indexImage, method, cfa = {"cfa": "bayer"}, zoomSize = "small", rows=1, cols=4, figsize=(15, 5)):
+        # Test key method
+        if method not in DataHandler.reconstruction_images.keys():
+            print("[ERROR] The method {} is not valid.".format(method))
+            exit(1)
+
+        res = DataHandler.get_reconstructed_image(self, indexImage, method)
+
+        fig, axs = plt.subplots(rows, cols, figsize=figsize)
+        fig.suptitle("Reconstructed Image with method: {} and pattern type: {}".format(method, cfa["cfa"]), fontsize=16)
+
+        axs[0].imshow(DataHandler.load_image(self, indexImage))
+        axs[0].set_title('Original Image')
+        axs[1].imshow(res)
+        axs[1].set_title('Reconstructed Image')
+
+        if zoomSize == "small":
+            axs[2].imshow(DataHandler.load_image(self, indexImage)[800:864, 450:514])
+            axs[2].set_title('Zoomed Input Image')
+            axs[3].imshow(res[800:864, 450:514])
+            axs[3].set_title('Zoomed Reconstructed Image')
+        else:
+            axs[2].imshow(DataHandler.load_image(self, indexImage)[2000:2064, 2000:2064])
+            axs[2].set_title('Zoomed Input Image')
+            axs[3].imshow(res[2000:2064, 2000:2064])
+            axs[3].set_title('Zoomed Reconstructed Image')
+
+        outputPath = os.path.join(DataHandler.WD_PATH, "output")
+        folderExists(outputPath)
+
+        imageName = Path(DataHandler.imagePaths[indexImage]).stem
+        plotCompImagePath = os.path.join(outputPath, "Compararison_" + imageName + "_" + method + "_" + cfa["cfa"] + ".png")
+
+        #save_image( plotReconstructedImagePath, res)
+        fig.savefig(plotCompImagePath)
+
+
+    def get_reconstructed_image(self, indexImage, method):
+        # Test key method
+        if method not in DataHandler.reconstruction_images.keys():
+            print("[ERROR] The method {} is not valid.".format(method))
+            exit(1)
+
+        DataHandler.indexImageExists(self, indexImage)
+
+        nameImage = Path(DataHandler.imagePaths[indexImage]).stem
+
+        if DataHandler.reconstruction_images[method].get(nameImage) is None:
+            print("[ERROR] There is no reconstruction image for the image {} and for the method {}.".format(nameImage, method))
+            exit(1)
+
+        return DataHandler.reconstruction_images[method][nameImage]
+
+    def shape_of_recontructed_image(self, indexImage, method):
+        if DataHandler.reconstruction_images[method].get(Path(DataHandler.imagePaths[indexImage]).stem) is None:
+            print("[ERROR] There is no reconstruction image for the image {} and for the method {}.".format(Path(DataHandler.imagePaths[indexImage]).stem, method))
+            exit(1)
+
+        return DataHandler.reconstruction_images[method][Path(DataHandler.imagePaths[indexImage]).stem].shape
+
+    def compute_CFA_image(self, CFA_NAME, nameImage, indexImage):
+        DataHandler.CFA_images.setdefault(nameImage, CFA(CFA_NAME, DataHandler.shape_of_image(self, indexImage)))
+
+    def shape_of_image(self, indexImage):
+        img = mpimg.imread(DataHandler.imagePaths[indexImage])
+        return img.shape
+
+    def load_image(self, indexImage):
+        img = imread(DataHandler.imagePaths[indexImage])
+        img = normalise_image(img)
+        return img
+
+    def indexImageExists(self, indexImage):
+        if indexImage >= len(DataHandler.imagePaths):
+            print("[ERROR] The index {} is not valid.".format(indexImage))

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/template/reconstruct.py → src/methods/brice_convers/reconstruct.py

@@ -5,8 +5,11 @@ Students can call their functions (declared in others files of src/methods/your_
 
 
 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:
@@ -19,18 +22,24 @@ def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
     Returns:
         np.ndarray: Demosaicked image.
     """
-    # Performing the reconstruction.
-    # TODO
+
     input_shape = (y.shape[0], y.shape[1], 3)
-    op = CFA(cfa, input_shape)
+    op = CFA("bayer", input_shape)
 
-    return np.zeros(op.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
 
 ####
 ####
 ####
-
 ####      ####                ####        #############
 ####      ######              ####      ##################
 ####      ########            ####      ####################

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