src/methods/Chardon_tom/reconstruct.py

+"""The main file for the reconstruction.
+This file should NOT be modified except the body of the 'run_reconstruction' function.
+Students can call their functions (declared in others files of src/methods/your_name).
+"""
+
+
+import numpy as np
+
+from src.forward_model import CFA
+from src.methods.Chardon_tom.utils import *
+import pywt
+
+#!!!!!!!! It is normal that the reconstructions lasts several minutes (3min on my computer)
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+
+    # Define constants and operators
+    cfa_name = 'bayer' # bayer or quad_bayer
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa_name, input_shape)
+
+    res = op.adjoint(y)
+
+    N,M = input_shape[0], input_shape[1]
+
+    
+
+    #interpolating green channel
+
+    for i in range (N):
+        for j in range (M):
+            if res[i,j,1] ==0:
+
+                neighbors = get_neighbors(res,1,i,j,N,M)
+                weights = get_weights(res,i,j,1,N,M)
+                res[i,j,1] = interpolate_green(weights, neighbors)
+
+
+
+    #first intepolation of red channel
+
+    for i in range (1,N,2):
+        for j in range (0,M,2):
+            neighbors = get_neighbors(res,0,i,j,N,M)
+            neighbors_G = get_neighbors(res,1,i,j,N,M)
+            weights = get_weights(res,i,j,0,N,M) 
+            res[i,j,0] = interpolate_red_blue(weights,neighbors, neighbors_G)
+
+    # second interpolation of red channel
+
+    for i in range (N):
+        for j in range (M):
+            if res[i,j,0] ==0:
+                neighbors = get_neighbors(res,0,i,j,N,M)
+                weights = get_weights(res,i,j,0,N,M)
+                res[i,j,0] = interpolate_green(weights, neighbors)
+
+
+    #first interpolation of blue channel
+
+    for i in range (0,N,2):
+        for j in range (1,M,2):
+            neighbors = get_neighbors(res,2,i,j,N,M)
+            neighbors_G = get_neighbors(res,1,i,j,N,M)
+            weights = get_weights(res,i,j,2,N,M) 
+            res[i,j,2] = interpolate_red_blue(weights, neighbors, neighbors_G)
+
+    #second interpolation of blue channel
+
+    for i in range (N):
+        for j in range (M):
+            if res[i,j,2] ==0:
+                neighbors = get_neighbors(res,2,i,j,N,M)
+                weights = get_weights(res,i,j,2,N,M)
+                res[i,j,2] = interpolate_green(weights,neighbors)
+
+
+
+    # k=0
+    # while k<2 :
+    #     for i in range(input_shape[0]):
+    #         for j in range(input_shape[1]):
+    #             res[i][j][1] = correction_green(res,i,j,N,M) 
+    #     for i in range(input_shape[0]):
+    #         for j in range(input_shape[1]):
+    #             res[i][j][0] = correction_red(res,i,j,N,M) 
+    #     for i in range(input_shape[0]):
+    #         for j in range(input_shape[1]):
+    #             res[i][j][2] = correction_blue(res,i,j,N,M) 
+    #     k+=1
+
+    res[res>1] = 1
+    res[res<0] = 0
+
+
+    return res
+

src/methods/Chardon_tom/utils.py

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

src/methods/EL-MURR-Theresa/reconstruct.py → src/methods/EL_MURR_Theresa/reconstruct.py

@@ -7,7 +7,7 @@ Students can call their functions (declared in others files of src/methods/your_
 import numpy as np
 
 from src.forward_model import CFA
-from src.methods.el_murr.malvar import malvar_he_cutler
+from src.methods.EL_MURR_Theresa.malvar import malvar_he_cutler
 
 def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
     """Performs demosaicking on y.

src/methods/GERAYELI_Kourosh/GERAYELI_Kourosh.py

+
+# Importing libraries
+import os
+import colour
+from colour_demosaicing import (
+    demosaicing_CFA_Bayer_bilinear,
+    demosaicing_CFA_Bayer_Malvar2004,
+    demosaicing_CFA_Bayer_Menon2007,
+    mosaicing_CFA_Bayer)
+
+from src.utils import psnr,ssim
+
+# Image path
+image_pathes = ['images/img_1.png','images/img_2.png','images/img_3.png','images/img_4.png']
+
+for i in image_pathes:
+    LIGHTHOUSE_IMAGE = colour.io.read_image(i)
+    img = LIGHTHOUSE_IMAGE
+    colour.plotting.plot_image(
+        colour.cctf_encoding(LIGHTHOUSE_IMAGE))
+
+
+    # Mosaicing
+    CFA = mosaicing_CFA_Bayer(LIGHTHOUSE_IMAGE)
+
+    colour.plotting.plot_image(
+        colour.cctf_encoding(CFA),
+        text_kwargs={'text': 'Lighthouse - CFA - RGGB'})
+
+    colour.plotting.plot_image(
+        colour.cctf_encoding(mosaicing_CFA_Bayer(LIGHTHOUSE_IMAGE, 'BGGR')), 
+        text_kwargs={'text': 'Lighthouse - CFA - BGGR'});
+
+
+    # Demosaicing bilinear
+    colour.plotting.plot_image(
+        colour.cctf_encoding(demosaicing_CFA_Bayer_bilinear(CFA)), 
+        text_kwargs={'text': 'Demosaicing - Bilinear'});
+
+    recons_bilinear = colour.cctf_encoding(demosaicing_CFA_Bayer_bilinear(CFA))
+
+
+    # demosaicing Malvar
+    recons_malvar = colour.cctf_encoding(demosaicing_CFA_Bayer_Malvar2004(CFA))
+    colour.plotting.plot_image(
+        colour.cctf_encoding(demosaicing_CFA_Bayer_Malvar2004(CFA)), 
+        text_kwargs={'text': 'Demosaicing - Malvar (2004)'});
+
+    # demosaicing Menon
+    recons_menon = colour.cctf_encoding(demosaicing_CFA_Bayer_Menon2007(CFA))
+    colour.plotting.plot_image(
+        colour.cctf_encoding(demosaicing_CFA_Bayer_Menon2007(CFA)), 
+        text_kwargs={'text': 'Demosaicing - Menon (2007)'});
+
+
+    print('bilinear : ')
+    print(f'PSNR: {psnr(img, recons_bilinear):.2f}')
+    print(f'SSIM: {ssim(img, recons_bilinear):.4f}')
+
+    print('Malvar : ')
+    print(f'PSNR: {psnr(img, recons_malvar):.2f}')
+    print(f'SSIM: {ssim(img, recons_malvar):.4f}')
+
+    print('Menon')
+    print(f'PSNR: {psnr(img, recons_menon):.2f}')
+    print(f'SSIM: {ssim(img, recons_menon):.4f}')
\ No newline at end of file

src/methods/RHOUCH_Oussama/clip_values.py

+import numpy as np
+
+def clip_values(img: np.ndarray) -> np.ndarray:
+    """Clip values of an image between 0 and 1
+    
+    Args:
+        img (np.ndarray): image to clip
+        
+    Returns:
+        np.ndarray: clipped image
+    """
+    # Clip values between 0 and 1
+    img[img < 0] = 0
+    img[img > 1] = 1 
+    return img

src/methods/RHOUCH_Oussama/find_direct_neighbors.py

+import numpy as np  
+
+def find_direct_neighbors(neighbors: list) -> list:
+    """
+    Calculate the gradients in horizontal, vertical, and diagonal directions for a pixel
+    based on its neighboring pixels.
+
+    Args:
+        neighbors (list): The neighbors of a pixel in a 3x3 grid. The order of neighbors
+                          is assumed to be [top-left, top, top-right, right, center, left,
+                          bottom-left, bottom, bottom-right].
+
+    Returns:
+        list: The gradients [horizontal, vertical, diagonal-x, diagonal-y].
+    """
+    
+    # Horizontal gradient: Difference between the right and left neighbors
+    horizontal_gradient = (neighbors[3] - neighbors[5]) / 2
+
+    # Vertical gradient: Difference between the top and bottom neighbors
+    vertical_gradient = (neighbors[1] - neighbors[7]) / 2
+
+    # Diagonal gradient along x-axis: Difference between top-right and bottom-left neighbors
+    diagonal_gradient_x = (neighbors[2] - neighbors[6]) / (2 * np.sqrt(2))
+
+    # Diagonal gradient along y-axis: Difference between top-left and bottom-right neighbors
+    diagonal_gradient_y = (neighbors[0] - neighbors[8]) / (2 * np.sqrt(2))
+    
+    return [horizontal_gradient, vertical_gradient, diagonal_gradient_x, diagonal_gradient_y]

src/methods/RHOUCH_Oussama/find_neighbors.py

+import numpy as np
+
+def find_neighbors(img: np.ndarray, channel: int, i: int, j: int, N: int, M: int) -> list:
+    """
+    Retrieve the 3x3 neighborhood of a pixel in a specified channel of an image.
+    The function wraps around the image edges to handle border pixels.
+
+    Args:
+        img (np.ndarray): The image to process.
+        channel (int): The index of the channel to process.
+        i (int): The row index of the pixel.
+        j (int): The column index of the pixel.
+        N (int): Height of the image.
+        M (int): Width of the image.
+        
+    Returns:
+        np.ndarray: An array of neighbor pixel values in the specified channel.
+    """
+    
+    neighbors = []
+
+    for di in [-1, 0, 1]:
+        for dj in [-1, 0, 1]:
+            neighbor_i = (i + di) % N
+            neighbor_j = (j + dj) % M
+
+            neighbors.append(img[neighbor_i, neighbor_j, channel])
+
+    return np.array(neighbors)

src/methods/RHOUCH_Oussama/find_weights.py

+import numpy as np
+from .find_neighbors import find_neighbors
+from .find_direct_neighbors import find_direct_neighbors
+
+def find_weights(img: np.ndarray, direct_neighbors: list, channel: int, i: int, j: int, N: int, M: int) -> list:
+    """
+    Find the weights of the neighbors of a pixel in the image.
+    
+    Args:
+        img (np.ndarray): The image to process.
+        direct_neighbors (list): The list of direct neighbors of the pixel.
+        channel (int): The index of the channel to process.
+        i (int): The row index of the pixel.
+        j (int): The column index of the pixel.
+        N (int): Height of the image.
+        M (int): Width of the image.
+        
+    Returns:
+        list: The list of weights of the neighbors of the pixel.
+    """
+    
+    [Dx, Dy, Dxx, Dyy] = direct_neighbors
+    E = []
+    c = 1
+    
+    for k in [-1, 0, 1]:
+        for l in [-1, 0, 1]:
+            n = find_neighbors(img, channel, i + k, j + l, N, M)
+            dd = find_direct_neighbors(n)
+
+            sqrt_arguments = {
+                1: 1 + Dyy * 2 + dd[3] * 2,
+                3: 1 + Dxx * 2 + dd[2] * 2,
+                2: 1 + Dy * 2 + dd[1] * 2,
+                4: 1 + Dx * 2 + dd[0] * 2
+            }
+
+            value = sqrt_arguments.get(c, 1)
+            if value < 0:
+                E.append(0)
+            else:
+                E.append(1 / np.sqrt(value))
+
+            c += 1
+            
+    return E

src/methods/RHOUCH_Oussama/interpolate.py

+def interpolate(neighbors: list, weights: list) -> float:
+    """
+    Interpolate the pixel value.
+    
+    Args:
+        neighbors (list): The neighbors of the pixel.
+        weights (list): The weights of the direct neighbors.
+        
+    Returns:
+        float: The interpolated value.
+    """
+    
+    return (weights[1] * neighbors[1] + weights[3] * neighbors[3] + weights[4] * neighbors[5] + weights[6] * neighbors[7]) / (weights[1] + weights[1] + weights[4] + weights[6])
+  
\ No newline at end of file

src/methods/RHOUCH_Oussama/interpolate_neighboring.py

+def interpolate_neighboring(neighbors_0: list, neighbors_1: list, weights: list) -> float:
+    """
+    Interpolate the neighboring pixels.
+    
+    Args:
+        neighbors_0 (list): The neighbors of the pixel in the channel to process.
+        neighbors_1 (list): The neighbors of the pixel in the green channel.
+        weights (list): The weights of the direct neighbors.
+        
+    Returns:
+        float: The interpolated value.
+    """
+    
+    return neighbors_1[4] * (((weights[0] * neighbors_0[0]) / neighbors_1[0]) + ((weights[2] * neighbors_0[2]) / neighbors_1[2]) + ((weights[5] * neighbors_0[6]) / neighbors_1[6]) + ((weights[7] * neighbors_0[8]) / neighbors_1[8])) / (weights[0] + weights[2] + weights[5] + weights[7])
+    
\ No newline at end of file

src/methods/RHOUCH_Oussama/main.ipynb

+%% Cell type:markdown id: tags:
+
+## Main notebook for experimenting with the algorithms
+
+Here is a simple notebook to test your code. You can modify it as you please.
+
+Remember to restart the jupyter kernel each time you modify a file.
+
+%% Cell type:code id: tags:
+
+``` python
+%%capture
+
+%pip install -r requirements.txt
+```
+
+%% Cell type:code id: tags:
+
+``` python
+import matplotlib.pyplot as plt
+
+from src.utils import load_image, save_image, psnr, ssim
+from src.forward_model import CFA
+from src.methods.RHOUCH_Oussama.reconstruct import run_reconstruction
+```
+
+%% Cell type:markdown id: tags:
+
+### Load the input image
+
+%% Cell type:code id: tags:
+
+``` python
+image_path = 'images/img_4.png'
+
+img = load_image(image_path)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# Shows some information on the image
+plt.imshow(img)
+plt.show()
+print(f'Shape of the image: {img.shape}.')
+```
+
+%% Cell type:markdown id: tags:
+
+### Definition of the forward model
+
+To setup the forward operator we just need to instanciate the `CFA` class. This class needs two arguments: `cfa_name` being the kind of pattern (bayer or kodak), and `input_shape` the shape of the inputs of the operator.
+
+This operation is linear and can be represented by a matrix $A$ but no direct access to this matrix is given (one can create it if needed). However the method `direct` allows to perform $A$'s operation. Likewise the method `adjoint` will perform the operation of $A^T$.
+
+For example let $X \in \mathbb R^{M \times N \times 3}$ the input RGB image in natural shape. Then we got $x \in \mathbb R^{3MN}$ (vectorized version of $X$) and $A \in \mathbb R^{MN \times 3MN}$, leading to:
+
+\begin{equation*}
+    y = Ax \in \mathbb R^{MN} \quad \text{and} \quad z = A^Ty  \in \mathbb R^{3MN}
+\end{equation*}
+
+However thanks to `direct` and `adjoint` there is no need to work with vectorized images, except if it is interesting to create the matrix $A$ explicitly.
+
+%% Cell type:code id: tags:
+
+``` python
+cfa_name = 'bayer' # bayer or quad_bayer
+op = CFA(cfa_name, img.shape)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# Shows the mask
+plt.imshow(op.mask[:10, :10])
+plt.show()
+print(f'Shape of the mask: {op.mask.shape}.')
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# Applies the mask to the image
+y = op.direct(img)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# Applies the adjoint operation to y
+z = op.adjoint(y)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+fig, axs = plt.subplots(2, 3, figsize=(15, 10))
+axs[0, 0].imshow(img)
+axs[0, 0].set_title('Input image')
+axs[0, 1].imshow(y, cmap='gray')
+axs[0, 1].set_title('Output image')
+axs[0, 2].imshow(z)
+axs[0, 2].set_title('Adjoint image')
+axs[1, 0].imshow(img[800:864, 450:514])
+axs[1, 0].set_title('Zoomed input image')
+axs[1, 1].imshow(y[800:864, 450:514], cmap='gray')
+axs[1, 1].set_title('Zoomed output image')
+axs[1, 2].imshow(z[800:864, 450:514])
+axs[1, 2].set_title('Zoomed adjoint image')
+plt.show()
+```
+
+%% Cell type:markdown id: tags:
+
+### Run the reconstruction
+
+Here the goal is to reconstruct the image `img` using only `y` and `op` (using `img` is forbidden).
+
+To run the reconstruction we simply call the function `run_reconstruction`. This function takes in argument the image to reconstruct and the kind of CFA used (bayer or kodak). All the parameters related to the reconstruction itself must be written inside `run_reconstruction`.
+
+%% Cell type:code id: tags:
+
+``` python
+res = run_reconstruction(y, cfa_name)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# Prints some information on the reconstruction
+print(f'Size of the reconstruction: {res.shape}.')
+```
+
+%% Cell type:markdown id: tags:
+
+### Quantitative and qualitative results
+
+%% Cell type:code id: tags:
+
+``` python
+fig, axs = plt.subplots(1, 2, figsize=(15, 15))
+axs[0].imshow(img)
+axs[0].set_title('Original image')
+axs[1].imshow(res)
+axs[1].set_title('Reconstructed image')
+plt.show()
+```
+
+%% Cell type:code id: tags:
+
+``` python
+# Computes some metrics
+print(f'PSNR: {psnr(img, res):.2f}')
+print(f'SSIM: {ssim(img, res):.4f}')
+```
+
+%% Cell type:code id: tags:
+
+``` python
+reconstructed_path = 'output/reconstructed_image.png'
+
+save_image(reconstructed_path, res)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+```