README.md

@@ -4,7 +4,7 @@
 
 As a final exam, this project aims to evaluate 3A SICOM students. 
 
-Students will work individually, and each student must submit its work (report and code) as a marge request on Gitlab by **Tuesday, June 6**.
+Students will work individually, and each student must submit its work (report and code) as a merge request on Gitlab by **Wednesday, January 31st**.
 
 ## 2. Presentation
 
@@ -18,15 +18,15 @@ Because of the filters the pixels on the sensor will only acquire one color (bei
 
 The goal of this project is to perform *demosaicking*: recover all the missing colors for each pixel. This way we will recover the full RGB image.
 
-We propose you to reconstruct 4 images (you can find them in the `images` folder of this project). These images are from the open dataset of the **National Gallery of Art**, USA, which can be found [here](https://github.com/NationalGalleryOfArt/opendata). To reconstruct these images we provide you the forward operator, modeling the effect of a CFA camera with 2 different CFA patterns:either Bayer of Quad Bayer pattern. This operation is described in `src/forward_operator.py`.
+We propose you to reconstruct 4 images (you can find them in the `images` folder of this project). These images are from the open dataset of the **National Gallery of Art**, USA, which can be found [here](https://github.com/NationalGalleryOfArt/opendata). To reconstruct these images we provide you the forward operator, modeling the effect of a CFA camera with 2 different CFA patterns: either Bayer of Quad Bayer pattern. This operation is described in `src/forward_operator.py`.
 
 ## 3. How to proceed?
 
-All of the images to reconstruct have the same size (1024x1024) and are RGB. You have to use methods seen during this semester to recover a fully colored RGB image using **only** the raw acquisition (gray-scale image) and the forward operator.
+All the images to reconstruct have the same size (1024x1024) and are RGB. You have to use methods seen during this semester to recover a fully colored RGB image using **only** the raw acquisition (gray-scale image) and the forward operator (but you are not forced to use it). Of course using directly the ground truth (original image) is forbiden, it can only be used to compute metrics (SSIM, PSNR).
 
 ## 4. Some hints to start
 
-We provide in `src/methods/baseline` a basic interpolation image which can help you to start. These methods are described in the PhD Thesis _Model Based Signal Processing Techniques for Nonconventional Optical Imaging Systems_ and the user's manual _Pyxalis Image Viewer_ (see references). You can also find in the academic litterature techniques (interpolation, inverse problem, machine learning, etc) that solves demosaicking.
+We provide in `src/methods/baseline` a basic interpolation image which can help you to start. These methods are described in the PhD Thesis _Model Based Signal Processing Techniques for Nonconventional Optical Imaging Systems_ and the user's manual _Pyxalis Image Viewer_ (see references). You can also find in the academic litterature techniques (interpolation, inverse problem, machine learning, etc) that solves problems close to this one.
 
 *The time of computations might be very long depending of your machine. To verify that your method works, try to work with smaller test images.*
 
@@ -44,48 +44,46 @@ sicom_image_analysis_project/
 ├─readme_imgs/                 # Images for the Readme
 ├─requirements.txt             # Requirement file for packages installation
 └─src/                         # All of the source files for the project
-  ├─reconstruct.py             # File containing the main reconstruction fonction
   ├─checks.py                  # File containing some sanity checks
   ├─forward_model.py           # File containing the CFA operator
   ├─utils.py                   # Some utilities
   └─methods/
-    ├─baseline/                # Example of demosaicking
+    ├─baseline/                # Example of reconstruction
     └─template/                # Template of your project (to be copied)
 ~~~
 
-Copy `src/methods/template` and rename it with your name. You are expected to write your own code in new files inside `src/methods/your_name`, and call them in the `src/methods/your_name/reconstruct.py` file. You are **not allowed to modify any file outside of `src/methods/your_name`** apart from `main.ipynb`, see instructions bellow.
-
-Please have a look in `src/methods/baseline`, your projet should work in the same way.
-
 ## 6. Instructions:
 
 Each student will fork this Git repository on its own account. You will then work in your own version of this repository.
 
-To submit your work you just need to create a merge request in Gitlab (see [here](https://docs.gitlab.com/ee/user/project/merge_requests/creating_merge_requests.html#when-you-work-in-a-fork) for a detailed explanation). This merge request will **only encompass the changes made to `src/methods/your_name`, without `main.ipynb`**. Because the merge request will only create `src/methods/your_name` there will not be any conflicts between each student's project.
-
 Along with your code you must submit a report **as a pdf file** with the name **name.pdf** inside the folder `src/methods/your_name`.
 
 The code and the report must be **written in English**.
 
-### 6.1. The report:
-
-Your report **must be a pdf file** written in English and should not be longer than 5 pages. In this file you are asked to explain:
-
-1. The problem statement, show us that you've understood the project.
-2. The solution that you've chosen, explaining the theory.
-3. With tools that **you've built** show us that your solution is relevant and working
-4. Results, give us your results of the images `img_1`, `img_2`, `img_3`, `img_4`.
-5. Conclusion, take a step back about your results. What can be improved?
-
-### 6.2. The code:
+### 6.1. The code:
 
 - You should first **copy** the folder `src/methods/template` and rename it with your name. It is in this folder that you will work, **nothing else should be modified** apart from `main.ipynb` for experimenting with the code.
 - You will add as many files you want in this folder but the only interface to run your code is `run_reconstruction` in `src/methods/your_name`. You can modify this function as you please, as long as it takes in argument the image to reconstruct (`y`), the name of the CFA (`cfa`) and returns a demosaicked image.
+- Have a look in `src/methods/baseline`, your projet should work in the same way.
 - You can use all the functions defined in the project, even the ones that you should not modify (`utils.py`, `forward_model.py`, etc).
 - Your code must be operational through `run_reconstruction`, as we will test it on new and private images.
 - The notebook provides a working bench. It should **not be included in the merge request**, it is just a work document for you.
 - Comment your code when needed.
 
+### 6.2. The report:
+
+Your report **must be a pdf file** written in English and should not be longer than 5 pages. In this file you are asked to explain:
+
+- The problem statement, show us that you've understood the project.
+- The solution that you've chosen, explaining the theory.
+- With tools that **you've built** show us that your solution is relevant and working
+- Results, give us your results of the images `img_1`, `img_2`, `img_3`, `img_4`.
+- Conclusion, take a step back about your results. What can be improved?
+
+### 6.3. Submission:
+
+To submit your work you just need to create a merge request in Gitlab (see [here](https://docs.gitlab.com/ee/user/project/merge_requests/creating_merge_requests.html#when-you-work-in-a-fork) for a detailed explanation). This merge request will **only encompass the changes made to `src/methods/your_name`, without `main.ipynb`**. Because the merge request will only create `src/methods/your_name` there will not be any conflicts between each student's project.
+
 ## 7. Supervisors
 
 - Mauro Dalla Mura: mauro.dalla-mura@gipsa-lab.grenoble-inp.fr

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 ..utils import load_image, save_image, psnr, ssim
+from src.utils import load_image, save_image, psnr, ssim
 from src.forward_model import CFA
 from src.methods.baseline.reconstruct import run_reconstruction
 ```
-%% Output
-
-    ---------------------------------------------------------------------------
-    ImportError                               Traceback (most recent call last)
-    /home/mullemat/code/sicom_image_analysis_project/src/methods/baseline/main.ipynb Cell 3 line 3
-          <a href='vscode-notebook-cell:/home/mullemat/code/sicom_image_analysis_project/src/methods/baseline/main.ipynb#W2sZmlsZQ%3D%3D?line=0'>1</a> import matplotlib.pyplot as plt
-    ----> <a href='vscode-notebook-cell:/home/mullemat/code/sicom_image_analysis_project/src/methods/baseline/main.ipynb#W2sZmlsZQ%3D%3D?line=2'>3</a> from ..utils import load_image, save_image, psnr, ssim
-          <a href='vscode-notebook-cell:/home/mullemat/code/sicom_image_analysis_project/src/methods/baseline/main.ipynb#W2sZmlsZQ%3D%3D?line=3'>4</a> from src.forward_model import CFA
-          <a href='vscode-notebook-cell:/home/mullemat/code/sicom_image_analysis_project/src/methods/baseline/main.ipynb#W2sZmlsZQ%3D%3D?line=4'>5</a> from src.methods.baseline.reconstruct import run_reconstruction
-    ImportError: attempted relative import with no known parent package
-
 %% Cell type:markdown id: tags:
 ### Load the input image
 %% Cell type:code id: tags:
 ``` python
 image_path = 'images/img_1.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}.')
 ```
 %% Output
     Shape of the image: (1024, 1024, 3).
 %% 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}.')
 ```
 %% Output
     Shape of the mask: (1024, 1024, 3).
 %% 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()
 ```
 %% Output
 %% 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}.')
 ```
 %% Output
     Size of the reconstruction: (1024, 1024, 3).
 %% 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()
 ```
 %% Output
 %% Cell type:code id: tags:
 ``` python
 # Computes some metrics
 print(f'PSNR: {psnr(img, res):.2f}')
 print(f'SSIM: {ssim(img, res):.4f}')
 ```
 %% Output
     PSNR: 34.63
     SSIM: 0.9502
 %% Cell type:markdown id: tags:
 ### Save the reconstructed image
 %% Cell type:code id: tags:
 ``` python
 reconstructed_path = 'output/reconstructed_image.png'
 save_image(reconstructed_path, res)
 ```

requirements.txt

+matplotlib
 numpy
-scipy
 scikit-image
-matplotlib
+scipy

src/checks.py

@@ -71,25 +71,14 @@ def check_data_range(img: np.ndarray) -> None:
         raise Exception(f'Pixel\'s values must be in range [0, 1]. Got range [{np.min(img)}, {np.max(img)}].')
 
 
-def check_len(img1_list: list | np.ndarray, img2_list: list | np.ndarray) -> None:
-    """Checks if two lists have the same length.
-
-    Args:
-        img1_list (list | np.ndarray): First list.
-        img2_list (list | np.ndarray): Secind list.
-
-    Raises:
-        Exception: Exception if the two list do not have the same length.
-    """
-    if len(img1_list) != len(img2_list):
-        raise Exception(f'The two lists must have the same length. Got {len(img1_list)} and {len(img2_list)}.')
-
-
 def check_cfa(cfa: str) -> None:
     """Checks if the CFA's name is correct.
 
     Args:
         cfa (str): CFA name.
+
+    Raises:
+        Exception: Exception if the name of the CFA is not correct.
     """
     if cfa not in ['bayer', 'quad_bayer']:
         raise Exception(f'Unknown CFA name. Got {cfa} but expected either bayer or quad_bayer.')

src/methods/CASA_NOVA_Vincent/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 functions as fu
+from src.forward_model import CFA
+import importlib 
+importlib.reload(fu)
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+    """Performs demosaicking on y.
+
+    Args:
+        y (np.ndarray): Mosaicked image to be reconstructed.
+        cfa (str): Name of the CFA. Can be bayer or quad_bayer.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+    # Performing the reconstruction.
+    # TODO
+
+    return fu.extrapolate_cfa(y,cfa)
+

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

+import numpy as np
+
+def find_Knearest_neighbors(z, chan, i, j, N, M):
+    """Finds a pixel's neighbors on a channel"""
+    return np.array([z[(i+di)%N, (j+dj)%M, chan] for di in range(-1, 2) for dj in range(-1, 2)])
+
+def calculate_directional_gradients(neighbors):
+    """Gives the directional derivative"""
+    P1, P2, P3, P4, P5, P6, P7, P8, P9 = neighbors
+    Dx, Dy = (P4 - P6)/2, (P2 - P8)/2
+    Dxd, Dyd = (P3 - P7)/(2*np.sqrt(2)), (P1 - P9)/(2*np.sqrt(2))
+    return [Dx, Dy, Dxd, Dyd]
+
+def calculate_adaptive_weights(z, neigh, dir_deriv,chan,i,j,N,M):
+
+    [Dx,Dy,Dxd,Dyd] = dir_deriv
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neigh
+    E = []
+    c = 1
+    for k in range (-1,2):
+        for k in range (-1,2):
+
+            n = find_Knearest_neighbors(z,chan,i+k,j+k,N,M)
+            dd = calculate_directional_gradients(n)
+            if c == 1 or c == 9:
+                E.append(1/np.sqrt(1 + Dyd**2 + dd[3]**2))
+            elif c == 2 or c == 8:
+                E.append(1/np.sqrt(1 + Dy**2 + dd[1]**2))
+            elif c == 3 or c == 7:
+                E.append(1/np.sqrt(1 + Dxd**2 + dd[2]**2))
+            elif c == 4 or c == 6:
+                E.append(1/np.sqrt(1 + Dx**2 + dd[0]**2))
+            c += 1
+    return E       
+
+def interpolate_pixel(neigh,weights):
+    """This function performs interpolation for a single pixel by calculating a weighted average of its neighboring pixels"""
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neigh
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+    num5 = E2*P2 + E4*P4 + E6*P6 + E8*P8
+    den5 = E2 + E4 + E6 + E8
+    I5 = num5/den5
+    return I5
+
+def interpolate_RedBlue(neighbors, neighbors_G, weights):
+    """This function specifically interpolates a pixel in the red or blue channels"""
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neighbors
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neighbors_G
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+    num5 = ((E1*P1)/G1) + ((E3*P3)/G3) + ((E7*P7)/G7) + ((E9*P9)/G9)
+    den5 = E1 + E3 + E7 + E9
+    I5 = G5 * num5/den5
+
+    return I5

src/methods/ELAMRANI_Mouna/reconstruct.py

+import numpy as np
+from src.forward_model import CFA
+from src.methods.ELAMRANI_Mouna.functions import *
+
+
+def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
+
+    cfa_name = 'bayer' # bayer or quad_bayer
+    input_shape = (y.shape[0], y.shape[1], 3)
+    op = CFA(cfa_name, input_shape)
+
+    img_res = op.adjoint(y)
+    N = img_res[:,:,0].shape[0]
+    M = img_res[:,:,0].shape[1]
+
+    def interpolate_channel(img_res, channel, first_pass, N, M):
+        for i in range(N):
+            for j in range(M):
+                if first_pass and ((channel == 0 and i % 2 == 1 and j % 2 == 0) or
+                               (channel == 2 and i % 2 == 0 and j % 2 == 1)):
+                    neighbors = find_Knearest_neighbors(img_res, channel, i, j, N, M)
+                    neighbors_G = find_Knearest_neighbors(img_res, 1, i, j, N, M)
+                    dir_deriv = calculate_directional_gradients(neighbors_G)
+                    weights = calculate_adaptive_weights(img_res, neighbors_G, dir_deriv, 1, i, j, N, M)
+                    img_res[i, j, channel] = interpolate_RedBlue(neighbors, neighbors_G, weights)
+                elif not first_pass and img_res[i, j, channel] == 0:
+                    neighbors = find_Knearest_neighbors(img_res, channel, i, j, N, M)
+                    dir_deriv = calculate_directional_gradients(neighbors)
+                    weights = calculate_adaptive_weights(img_res, neighbors, dir_deriv, channel, i, j, N, M)
+                    img_res[i, j, channel] = interpolate_pixel(neighbors, weights)
+        return img_res
+
+# Interpolation pour chaque canal
+    img_res = interpolate_channel(img_res, 1, False, N, M)  # Interpolation du canal vert
+    img_res = interpolate_channel(img_res, 0, True, N, M)   # Première interpolation du canal rouge
+    img_res = interpolate_channel(img_res, 0, False, N, M)  # Seconde interpolation du canal rouge
+    img_res = interpolate_channel(img_res, 2, True, N, M)   # Première interpolation du canal bleu
+    img_res = interpolate_channel(img_res, 2, False, N, M)  # Seconde interpolation du canal bleu
+
+    img_res[img_res > 1] = 1
+    img_res[img_res < 0] = 0
+
+    return img_res
+

src/methods/EL_MURR_Theresa/malvar.py

+import numpy as np
+from scipy.signal import correlate2d
+from src.forward_model import CFA
+
+def malvar_he_cutler(y: np.ndarray, op: CFA ) -> np.ndarray:
+    """Performs demosaicing using the malvar-he-cutler algorithm
+
+    Args:
+        op (CFA): CFA operator.
+        y (np.ndarray): Mosaicked image.
+
+    Returns:
+        np.ndarray: Demosaicked image.
+    """
+
+    red_mask, green_mask, blue_mask =  [op.mask[:, :, 0], op.mask[:, :, 1], op.mask[:, :, 2]]
+    mosaicked_image = np.float32(y)
+    demosaicked_image = np.empty(op.input_shape)
+
+    if op.cfa == 'quad_bayer':
+        filters = get_quad_bayer_filters()
+    else:
+        filters = get_default_filters()
+
+    demosaicked_image = apply_demosaicking_filters(
+        mosaicked_image,demosaicked_image, red_mask, green_mask, blue_mask, filters
+    )
+
+    return demosaicked_image
+
+def get_quad_bayer_filters():
+    coefficient_scale = 0.03125
+    return {
+        "G_at_R_and_B": np.array([
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, 2, 2, 0, 0, 0, 0],
+            [0, 0, 0, 0, 2, 2, 0, 0, 0, 0],
+            [-1, -1, 2, 2, 4, 4, 2, 2, -1, -1],
+            [-1, -1, 2, 2, 4, 4, 2, 2, -1, -1],
+            [0, 0, 0, 0, 2, 2, 0, 0, 0, 0],
+            [0, 0, 0, 0, 2, 2, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0]
+        ]) * coefficient_scale,
+        "R_at_GR_and_B_at_GB": np.array([
+            [0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0],
+            [0, 0, -1, -1, 0, 0, -1, -1, 0, 0],
+            [0, 0, -1, -1, 0, 0, -1, -1, 0, 0],
+            [-1, -1, 4, 4, 5, 5, 4, 4, -1, -1],
+            [-1, -1, 4, 4, 5, 5, 4, 4, -1, -1],
+            [0, 0, -1, -1, 0, 0, -1, -1, 0, 0],
+            [0, 0, -1, -1, 0, 0, -1, -1, 0, 0],
+            [0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0],
+            [0, 0, 0, 0, 0.5, 0.5, 0, 0, 0, 0]
+        ]) * coefficient_scale,
+        "R_at_GB_and_B_at_GR": np.array([
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, -1, -1, 4, 4, -1, -1, 0, 0],
+            [0, 0, -1, -1, 4, 4, -1, -1, 0, 0],
+            [0.5, 0.5, 0, 0, 5, 5, 0, 0, 0.5, 0.5],
+            [0.5, 0.5, 0, 0, 5, 5, 0, 0, 0.5, 0.5],
+            [0, 0, -1, -1, 4, 4, -1, -1, 0, 0],
+            [0, 0, -1, -1, 4, 4, -1, -1, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1, -1, 0, 0, 0, 0]
+        ]) * coefficient_scale,
+        "R_at_B_and_B_at_R": np.array([
+            [0, 0, 0, 0, -1.5, -1.5, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1.5, -1.5, 0, 0, 0, 0],
+            [0, 0, 2, 2, 0, 0, 2, 2, 0, 0],
+            [0, 0, 2, 2, 0, 0, 2, 2, 0, 0],
+            [-1.5, -1.5, 0, 0, 6, 6, 0, 0, -1.5, -1.5],
+            [-1.5, -1.5, 0, 0, 6, 6, 0, 0, -1.5, -1.5],
+            [0, 0, 2, 2, 0, 0, 2, 2, 0, 0],
+            [0, 0, 2, 2, 0, 0, 2, 2, 0, 0],
+            [0, 0, 0, 0, -1.5, -1.5, 0, 0, 0, 0],
+            [0, 0, 0, 0, -1.5, -1.5, 0, 0, 0, 0]
+        ]) * coefficient_scale,
+    }
+
+def get_default_filters():
+    coefficient_scale = 0.125
+    return {
+        "G_at_R_and_B": 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]
+        ]) * coefficient_scale,
+        "R_at_GR_and_B_at_GB": 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]
+        ]) * coefficient_scale,
+        "R_at_GB_and_B_at_GR": 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]
+        ]) * coefficient_scale,
+        "R_at_B_and_B_at_R": np.array([
+            [0, 0, -1.5, 0, 0],
+            [0, 2, 0, 2, 0],
+            [-1.5, 0, 6, 0, -1.5],
+            [0, 2, 0, 2, 0],
+            [0, 0, -1.5, 0, 0]
+        ]) * coefficient_scale,
+    }
+
+def apply_demosaicking_filters(image, res, red_mask, green_mask, blue_mask, filters):
+    red_channel = image * red_mask
+    green_channel = image * green_mask
+    blue_channel = image * blue_mask
+
+    # Create the green channel after applying a filter
+    green_channel = np.where(
+        np.logical_or(red_mask == 1, blue_mask == 1),
+        correlate2d(image, filters['G_at_R_and_B'], mode="same", boundary="symm"),
+        green_channel
+    )
+
+
+    # Define masks for extracting pixel values
+    red_row_mask = np.any(red_mask == 1, axis=1)[:, np.newaxis].astype(np.float32)
+    red_col_mask = np.any(red_mask == 1, axis=0)[np.newaxis].astype(np.float32)
+
+    blue_row_mask = np.any(blue_mask == 1, axis=1)[:, np.newaxis].astype(np.float32)
+    blue_col_mask = np.any(blue_mask == 1, axis=0)[np.newaxis].astype(np.float32)
+
+    def update_channel(channel, row_mask, col_mask, filter_key):
+        return np.where(
+            np.logical_and(row_mask == 1, col_mask == 1),
+            correlate2d(image, filters[filter_key], mode="same", boundary="symm"),
+            channel
+        )
+
+# Update the red channel and blue channel
+    red_channel = update_channel(red_channel, red_row_mask, blue_col_mask, 'R_at_GR_and_B_at_GB')
+    red_channel = update_channel(red_channel, blue_row_mask, red_col_mask, 'R_at_GB_and_B_at_GR')
+
+    blue_channel = update_channel(blue_channel, blue_row_mask, red_col_mask, 'R_at_GR_and_B_at_GB')
+    blue_channel = update_channel(blue_channel, red_row_mask, blue_col_mask, 'R_at_GB_and_B_at_GR')
+
+    # Update R channel and B channel again
+    red_channel = update_channel(red_channel, blue_row_mask, blue_col_mask, 'R_at_B_and_B_at_R')
+    blue_channel = update_channel(blue_channel, red_row_mask, red_col_mask, 'R_at_B_and_B_at_R')
+    res[:, :, 0] = red_channel
+    res[:, :, 1] = green_channel
+    res[:, :, 2] = blue_channel
+    return res
\ No newline at end of file

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