diff --git a/src/methods/laporte_auriane/Project_report_LAPORTE.pdf b/src/methods/laporte_auriane/Project_report_LAPORTE.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..a60362af40c78f667e3dbb1b1cf15a849c8ef585
Binary files /dev/null and b/src/methods/laporte_auriane/Project_report_LAPORTE.pdf differ
diff --git a/src/methods/laporte_auriane/fct_to_demosaick.py b/src/methods/laporte_auriane/fct_to_demosaick.py
new file mode 100644
index 0000000000000000000000000000000000000000..22c383c97dd108a10b58c5be6cbc5b991eafd422
--- /dev/null
+++ b/src/methods/laporte_auriane/fct_to_demosaick.py
@@ -0,0 +1,332 @@
+"""This file contains all the functions that are used in order to perform a Kimmel algorithm demosaicking"""
+
+##### Importations #####
+import numpy as np
+import matplotlib.pyplot as plt
+from src.forward_model import CFA
+
+##### Some functions used #####
+def neigh(i, j, colour_chan, img):
+
+ """Finds the neighbours of the pixel (i, j) of a colour channel
+ Args:
+ i: column of the pixel
+ j: line of the pixel
+ colour_chan: colour channel chosen
+ img: image
+ Returns:
+ List of the neighbours of the pixel (i, j) and the pixel (i, j)
+ """
+
+ X = img[:, :, colour_chan].shape[0]
+ Y = img[:, :, colour_chan].shape[1]
+
+ # Place each Pi compared to P5 (i horizontal and j vertical)
+ # See P1 to P9 drawing, page 2 of http://ultra.sdk.free.fr/docs/Image-Processing/Demosaic/An%20Improved%20Demosaicing%20Algorithm.pdf
+ # % (modulo) to make sure it is a multiple (avoid error index 1024 out of bounds)
+ P1 = img[(i-1)%Y, (j-1)%X, colour_chan]
+ P2 = img[i%Y, (j-1)%X, colour_chan]
+ P3 = img[(i+1)%Y, (j-1)%X, colour_chan]
+ P4 = img[(i-1)%Y, j%X, colour_chan]
+ P5 = img[i%Y, j%X, colour_chan]
+ P6 = img[(i+1)%Y, j%X, colour_chan]
+ P7 = img[(i-1)%Y, (j+1)%X, colour_chan]
+ P8 = img[i%Y, (j+1)%X, colour_chan]
+ P9 = img[(i+1)%Y, (j+1)%X, colour_chan]
+ #print(f'Size of P1: {P1.size}') # 1
+
+ return [P1, P2, P3, P4, P5, P6, P7, P8, P9]
+
+def derivat(neighbour_list):
+
+ """Compute directional derivatives of pixel
+ Args:
+ neighbour_list: list of neighbours
+ Returns:
+ List of the directional derivatives
+ """
+
+ # Must be 9 (9P neigbours)
+ #print(f'neighbour_list len to derivate: {len(neighbour_list)}') # 9
+
+ # Assign P neighbours to the neighbour_list (to then compute the derivatives)
+ [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbour_list
+
+ # Compute the derivatives (independant of the colour component)
+ Dx = (P4 - P6) / 2
+ Dy = (P2 - P8) / 2
+ Dxd = (P3 - P7) / (2 * np.sqrt(2))
+ Dyd = (P1 - P9) / (2 * np.sqrt(2))
+ # Formula given at the top of page 3 (but very long computing time, and does not provide better results)
+ #Dxd = np.max( ( np.abs((P3 - P5) / (np.sqrt(2))), np.abs((P7 - P5) / (np.sqrt(2))) ) )
+ #Dyd = np.max( ( np.abs((P1 - P5) / (np.sqrt(2))), np.abs((P9 - P5) / (np.sqrt(2))) ) )
+
+ # Store the computed values in deriv_list
+ deriv_list = [Dx, Dy, Dxd, Dyd]
+ #print(f'Len de deriv_list: {len(deriv_list)}') # 4
+
+ return deriv_list
+
+def weight_fct(i, j, neighbour_list, deriv_list, colour_chan, img):
+
+ """Compute the weights Ei
+ Args:
+ i: column of the pixel
+ j: line of the pixel
+ neighbour_list: neighbour_list: list of neighbours
+ deriv_list: list of derivatives
+ colour_chan:colour channel chosen
+ img: image
+ Returns:
+ List of the weights
+ """
+
+ # Assignment
+ [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbour_list
+ [Dx, Dy, Dxd, Dyd] = deriv_list
+ X = img[:,:,colour_chan].shape[0]
+ Y = img[:,:,colour_chan].shape[1]
+
+ # E list to complete (weighting function)
+ E = []
+
+ # Fix the start
+ cur = 1
+
+ # Neighbourhood
+ for step in range(-1, 2): # from -1 to 1 (2-1)
+ for step in range(-1, 2): # otherwise only 3 values in E
+ # Find the neigbours
+ n = neigh(i+step, j+step, colour_chan, img)
+ # Derivatives
+ d = derivat(n)
+
+ # See P1 to P9 drawing, page 2 of http://ultra.sdk.free.fr/docs/Image-Processing/Demosaic/An%20Improved%20Demosaicing%20Algorithm.pdf
+ if cur==4 or cur==6:
+ E.append( 1 / np.sqrt(1 + Dx**2 + d[0]**2) )
+ elif cur==2 or cur==8:
+ E.append( 1 / np.sqrt(1 + Dy**2 + d[1]**2) )
+ elif cur==7 or cur==3:
+ E.append( 1 / np.sqrt(1 + Dxd**2 + d[2]**2) )
+ else: # 1 or 9
+ E.append( 1 / np.sqrt(1 + Dyd**2 + d[3]**2) )
+ cur = cur+1
+
+ #print(f'Len of E: {len(E)}') # 9
+ # other way to code (Page 3/6 of http://elynxsdk.free.fr/ext-docs/Demosaicing/more/news0/Optimal%20Recovery%20Demosaicing.pdf, method 1)), but I did not understand what is or how T was chosen?
+
+ return E
+
+def clip_extreme(img):
+
+ """Clip the values of the image, so they are between 0-1
+ Args:
+ img: image
+ Returns:
+ The clipped image
+ """
+
+ #max_val = np.amax(img)
+ #min_val = np.amin(img)
+ #print(f'Max value before clipping: {max_val}')
+ #print(f'Min value before clipping: {min_val}')
+
+ # Put extreme values between 0 and 1
+ img = np.clip(img, 0, 1)
+
+ return img
+
+def quad_to_bayer(img, op):
+
+ """Performs a conversion from Quad_bayer to Bayer pattern
+ Args:
+ img: image
+ op: CFA op
+ Returns:
+ The image Bayer pattern
+ """
+
+ # Based on p28/32 https://pyxalis.com/wp-content/uploads/2021/12/PYX-ImageViewer-User_Guide.pdf
+ ### 1. Swap 2 col every 2 columns
+ # Start: 2nd col --> 1 (bc start counting at 0) / Step: 4 (look at drawing of 1st swap p28 Pyxalis) / End: input_shape[0]-2 --> input_shape[0]-1
+ for j in range (1, img.shape[0]-1, 4):
+ store = np.copy(img[:, j]) # Store col j of the img
+ img[:, j] = np.copy(img[:, j+1]) # Col j becomes like col j+1
+ img[:, j+1] = store # Col j+1 become like col j (previously saved otherwise value lost)
+ store_op = np.copy(op.mask[:, j])
+ op.mask[:, j] = np.copy(op.mask[:, j+1])
+ op.mask[:, j+1] = store_op
+
+ ### 2. Swap 2 lines every 2 lines (same process for the lines)
+ for i in range (1, img.shape[1]-1, 4):
+ store = np.copy(img[i, :])
+ img[i, :] = np.copy(img[i+1, :])
+ img[i+1, :] = store
+ store_op = np.copy(op.mask[i, :])
+ op.mask[i, :] = np.copy(op.mask[i+1, :])
+ op.mask[i+1, :] = store_op
+
+ ### 3. Swap back some diagonal greens
+ # Starts: 0 and 2 / Steps: 4 / End: img.shape[nb]-1 --> img.shape[nb]
+ for k in range(0, img.shape[0], 4):
+ for l in range(2, img.shape[1], 4):
+ store = np.copy(img[k, l])
+ img[k, l] = np.copy(img[k+1, l-1])
+ img[k+1, l-1] = store
+ store_op = op.mask[k, l]
+ op.mask[k, l] = np.copy(op.mask[k+1, l-1])
+ op.mask[k+1, l-1] = store_op
+
+ return img
+
+
+##### 1. Interpolation of green colour #####
+def interpol_green(neighbour_list, weights_list):
+
+ """Interpolates missing green pixel
+ Args:
+ neighbour_list: list of neighbours
+ weights_list: list of weights
+ Returns:
+ The interpolated pixel
+ """
+
+ # Assignment
+ [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbour_list
+ [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+
+ # Compute missing green pixel (combination of 4 nearest neighbours)
+ G5 = (E2*P2 + E4*P4 + E6*P6 + E8*P8) / (E2+E4+E6+E8)
+ #print(f"G5 size: {G5.size}") # 1
+
+ return G5
+
+
+##### 2. Interpolation of red and blue colours #####
+def interpol_red_blue(neighbour_list, weights_list, green_neighbour):
+
+ """Interpolates missing blue/red pixel
+ Args:
+ neighbour_list: list of neighbours in blue/red channel
+ weights_list: list of weights
+ green_neighbour: list of the neighbours in green channel
+ Returns:
+ The interpolated pixel (in blue/red channel)
+ """
+
+ # Assignment
+ [P1, P2, P3, P4, P5, P6, P7, P8, P9] = neighbour_list
+ [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+ [G1, G2, G3, G4, G5, G6, G7, G8, G9] = green_neighbour
+
+ # Compute missing red/blue pixel (combination of 4 nearest neighbours)
+ R5_num = E1*(P1/G1) + E3*(P3/G3) + E7*(P7/G7) + E9*(P9/G9)
+ R5_denom = E1+E3+E7+E9
+ R5 = G5 * R5_num / R5_denom
+ #print(f"R5 size: {R5.size}") # 1
+
+ return R5
+
+
+##### 3. Correction stage #####
+def green_correction(red_neighbour, green_neighbour, blue_neighbour, weights_list):
+
+ """Correction of green pixels
+ Args:
+ red_neighbour: neighbours in red channel
+ green_neighbour: neighbours in green channel
+ blue_neighbour: neighbours in blue channel
+ weights_list: list of weights
+ Returns:
+ Corrected green pixel
+ """
+
+ # Assignment
+ [G1, G2, G3, G4, G5, G6, G7, G8, G9] = green_neighbour
+ [R1, R2, R3, R4, R5, R6, R7, R8, R9] = red_neighbour
+ [B1, B2, B3, B4, B5, B6, B7, B8, B9] = blue_neighbour
+ [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+
+ # Compute correction
+ GB5_num = E2*(G2/B2) + E4*(G4/B4) + E6*(G6/B6) + E8*(G8/B8)
+ GR5_num = E2*(G2/R2) + E4*(G4/R4) + E6*(G6/R6) + E8*(G8/R8)
+ G5_denom = E2+E4+E6+E8
+ GB5 = B5 * GB5_num / G5_denom
+ GR5 = R5 * GR5_num / G5_denom
+
+ #print(f"GB5 size: {GB5.size}") # 1
+ print(f"GB5: {GB5}")
+ #print(f"GR5 size: {GR5.size}") # 1
+ print(f"GR5: {GR5}")
+
+ G5 = (GR5 + GB5) / 2
+ print(f"G5 size: {G5.size}") # 1
+ print(f"G5: {G5}")
+
+ return G5
+
+def blue_correction(green_neighbour, blue_neighbour, weights_list):
+
+ """Correction of blue pixels
+ Args:
+ green_neighbour: neighbours in green channel
+ blue_neighbour: neighbours in blue channel
+ weights_list: list of weights
+ Returns:
+ Corrected blue pixel
+ """
+
+ # Assignment
+ [G1, G2, G3, G4, G5, G6, G7, G8, G9] = green_neighbour
+ [B1, B2, B3, B4, B5, B6, B7, B8, B9] = blue_neighbour
+ [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+
+ # Compute correction
+ B5_num = 0
+ B5_denom = 0
+ print(f'Len de weights_list: {len(weights_list)}')
+
+ for i in range(len(weights_list)):
+ if i != 5:
+ B5_num += weights_list[i] * blue_neighbour[i] / green_neighbour[i]
+ B5_denom += weights_list[i]
+
+ B5 = G5 * B5_num / B5_denom
+
+ #print(f"B5_denom size: {B5_denom.size}") # 1
+ print(f"B5_denom: {B5_denom}") # number
+ #print(f"B5_num size: {B5_num.size}") # 1
+ print(f"B5_num: {B5_num}") # nan
+ #print(f"B5 size: {B5.size}") # 1
+ print(f"B5: {B5}")
+
+ return B5
+
+def red_correction(red_neighbour, green_neighbour, weights_list):
+
+ """Correction of red pixels
+ Args:
+ red_neighbour: neighbours in red channel
+ green_neighbour: neighbours in green channel
+ weights_list: list of weights
+ Returns:
+ Corrected red pixel
+ """
+
+ # Assignment
+ [G1, G2, G3, G4, G5, G6, G7, G8, G9] = green_neighbour
+ [R1, R2, R3, R4, R5, R6, R7, R8, R9] = red_neighbour
+ [E1, E2, E3, E4, E5, E6, E7, E8, E9] = weights_list
+
+ # Compute correction
+ R5_num = 0
+ R5_denom = 0
+ for i in range(len(weights_list)):
+ if i != 5:
+
+ R5_num += weights_list[i] * red_neighbour[i] / green_neighbour[i]
+ R5_denom += weights_list[i]
+
+ R5 = G5 * R5_num / R5_denom
+
+ return R5
\ No newline at end of file
diff --git a/src/methods/laporte_auriane/reconstruct.py b/src/methods/laporte_auriane/reconstruct.py
index a97bd3f6e3c68df763b36c46b2727461af078bd2..bac96a08738de8d3540462cc6617090ea910a210 100644
--- a/src/methods/laporte_auriane/reconstruct.py
+++ b/src/methods/laporte_auriane/reconstruct.py
@@ -3,10 +3,9 @@ This file should NOT be modified except the body of the 'run_reconstruction' fun
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.laporte_auriane.fct_to_demosaick import *
def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
@@ -19,12 +18,121 @@ def run_reconstruction(y: np.ndarray, cfa: str) -> np.ndarray:
Returns:
np.ndarray: Demosaicked image.
"""
- # Performing the reconstruction.
- # TODO
+
+ #print(f'y shape: {y.shape}') # (1024, 1024)
input_shape = (y.shape[0], y.shape[1], 3)
op = CFA(cfa, input_shape)
+
+ ########## TODO ##########
+
+ # Colour channels
+ red = 0
+ green = 1
+ blue = 2
+ correction = 0
+
+ # In case of Quad_bayer op
+ if op.cfa == 'quad_bayer':
+ print("Transformation into Bayer pattern")
+ y = quad_to_bayer(y, op)
+ print('Quad_bayer to Bayer done\n\nDemoisaicking processing...')
+
+ z = op.adjoint(y)
+ # 1st and 2nd dim of a colour channel
+ X = z[:,:,0].shape[0]
+ Y = z[:,:,0].shape[1]
+ #print(f"X: {X}") # 1024
+ #print(f"Y: {Y}") # 1024
+
+
+ ### Kimmel algorithm ###
+ # About 3 to 4 minutes
+
+
+ # Interpolation of green colour
+ for i in range (X):
+ for j in range (Y):
+ # Control if value 0 in green channel
+ if z[i,j,1] == 0:
+
+ # Find neigbours
+ n = neigh(i, j, green, z)
+ # Derivatives
+ d = derivat(n)
+ # Weighting fct
+ w = weight_fct(i, j, n, d, green, z)
+ # Green interpolation (cross)
+ z[i,j,1] = interpol_green(n, w)
+ # Clip (avoid values our of range 0-1)
+ z = clip_extreme(z)
+
+
+ # 2 steps for the blues
+ # Interpolate missing blues at the red locations
+ for i in range (0, X, 2):
+ for j in range (1, Y, 2):
+
+ # Find neigbours + green
+ n = neigh(i, j, blue, z)
+ n_G = neigh(i, j, green, z)
+ d = derivat(n_G)
+ w = weight_fct(i, j, n, d, green, z)
+ # Blue interpolation (4 corners)
+ z[i,j,2] = interpol_red_blue(n, w, n_G)
+ # Interpolate missing blues at the green locations
+ for i in range (X):
+ for j in range (Y):
+ # Control if value 0 in blue channel
+ if z[i,j,2] == 0:
+
+ # Find neigbours
+ n_B = neigh(i, j, blue, z)
+ d = derivat(n_B)
+ w = weight_fct(i, j, n_B, d, blue, z)
+ # Blue interpolation (cross)
+ z[i,j,2] = interpol_green(n_B,w)
+ z = clip_extreme(z)
+
+
+ # 2 steps for the reds
+ # Interpolate missing reds at the blue locations
+ for i in range (1, X, 2):
+ for j in range (0, Y, 2):
+
+ # Find neigbours + green
+ n = neigh(i, j, red, z)
+ n_G = neigh(i, j, green, z)
+ d = derivat(n_G)
+ w = weight_fct(i, j, n_G, d, green, z)
+ # Red interpolation (4 corners)
+ z[i,j,0] = interpol_red_blue(n, w, n_G)
+ # Interpolate missing reds at the green locations
+ for i in range (X):
+ for j in range (Y):
+ # Control if value 0 in red channel
+ if z[i,j,0] == 0:
+
+ # Find neigbours
+ n_R = neigh(i, j, red, z)
+ d = derivat(n_R)
+ w = weight_fct(i, j, n_R, d, red, z)
+ # Red interpolation (cross)
+ z[i,j,0] = interpol_green(n_R,w)
+ z = clip_extreme(z)
+
+
+ # Correction step (repeated 3 times)
+ if correction == 1:
+ for i in range(3):
+
+ n_G[4] = green_correction(n_R, n_G, n_B, w)
+ n_B[4] = blue_correction(n_G, n_B, w)
+ n_R[4] = red_correction(n_R, n_G, w)
+ # nan
+
+ #print(f'Image reconstructed shape: {z.shape}') # 1024, 1024, 3
- return np.zeros(op.input_shape)
+ return z
####