Skip to content
Snippets Groups Projects
Commit dee7c15c authored by Matthieu Muller's avatar Matthieu Muller
Browse files

Merge branch 'master' into 'master'

Rendu projet

See merge request !14
parents 0ded0e79 eb661f59
No related branches found
No related tags found
1 merge request!14Rendu projet
This diff is collapsed.
######################################################################################
# Demosaicking algorithm based on "A Demosaicking Algorithme with
# Adaptive Inter-Channel Correlation" paper from Joan Duran and Antoni Buades
######################################################################################
import numpy as np
def cvt_YUV_space(R: np.ndarray, G: np.ndarray, B: np.ndarray):
"""Convert RGB into YUV
Args:
R (np.ndarray): Red channel of image
G (np.ndarray): Green channel of image
B (np.ndarray): Blue channel of image
Returns:
np.ndarray: Image in the YUV space
"""
Y = 0.299*R + 0.587*G + 0.114*B
U = R-Y
V = B-Y
return Y, U, V
def compute_gradient(L: int, U: np.ndarray, V: np.ndarray, i: int, j: int, dir: str):
"""Compute the gradient
Args:
L (int): Size of the local neighborhood
U (np.ndarray): U channel of the image
V (np.ndarray): V channel of the image
i (int): position in line
j (int): position in column
dir (string) : direction of interpolation
Returns:
float: Gradient
"""
height, width = U.shape
gradient_U = 0
gradient_V = 0
if dir == "north":
if j < L+1: max = j+1
else : max = L+1
for l in range(1, max):
gradient_U += (U[j-l, i] - U[j, i])**2
gradient_V += (V[j-l, i] - V[j, i])**2
elif dir == "south":
if height-j < L+1: max = height-j
else: max = L+1
for l in range(1, max):
gradient_U += (U[j+l, i] - U[j, i])**2
gradient_V += (V[j+l, i] - V[j, i])**2
elif dir == "east":
if width-i < L+1: max = width-i
else: max = L+1
for l in range(1, max):
gradient_U += (U[j, i+l] - U[j, i])**2
gradient_V += (V[j, i+l] - V[j, i])**2
elif dir == "west":
if i < L+1: max = i+1
else: max = L+1
for l in range(1, max):
gradient_U += (U[j, i-l] - U[j, i])**2
gradient_V += (V[j, i-l] - V[j, i])**2
if max == L+1: max = L
gradient = (np.sqrt(gradient_U) + np.sqrt(gradient_V)) / max
return gradient
def local_int(R0: np.ndarray, G0: np.ndarray, B0: np.ndarray, beta: float, L: int , eps: float):
"""Local direction interpolation algorithm
Args:
R0 (np.ndarray): Red mosaiked image
G0 (np.ndarray): Green mosaiked image
B0 (np.ndarray): Blue mosaiked image
beta (float): Channel correlation parameter [0; 1]
L (int): Size of the local neighborhood
eps (float): tresholding parameter > 0
Returns:
np.ndarray: The interpolated image (R, G, B)
"""
assert R0.shape == G0.shape and R0.shape == B0.shape
assert eps > 0
lin, col = G0.shape
# Initalisation
Rn, Rs, Re, Rw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
Gn, Gs, Ge, Gw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
Bn, Bs, Be, Bw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
RGn, RGs, RGe, RGw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
BGn, BGs, BGe, BGw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
Gradn, Grads, Grade, Gradw = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
low_wn, low_ws, low_we, low_ww = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
W = np.zeros((lin, col))
R, G, B = np.zeros((lin, col)), np.zeros((lin, col)), np.zeros((lin, col))
# Computation
for i in range(col):
for j in range(lin):
# Directionnal interpolation of green channel
if G0[j, i] != 0:
Gn[j, i], Gs[j, i], Ge[j, i], Gw[j, i] = G0[j, i], G0[j, i], G0[j, i], G0[j, i]
else:
if j-2 >= 0: # North
# If we are on a Red square: B0[j, i] = B0[j-2, i] = 0
# The same if we are on a Blue square
Gn[j, i] = G0[j-1, i] + 0.5*beta*(R0[j, i] - R0[j-2, i] + B0[j, i] - B0[j-2, i])
else: # R0[j-2, i] = 0 or B0[j-2, i] = 0
if j-1 >= 0:
Gn[j, i] = G0[j-1, i] + beta*(R0[j, i] + B0[j, i])
else: # G0[j-1, i] = 0
Gn[j, i] = beta*(R0[j, i] + B0[j, i])
if j+2 < lin: # South
Gs[j, i] = G0[j+1, i] + 0.5*beta*(R0[j, i] - R0[j+2, i] + B0[j, i] - B0[j+2, i])
else:
if j+1 < lin:
Gs[j, i] = G0[j+1, i] + beta*(R0[j, i] + B0[j, i])
else:
Gs[j, i] = beta*(R0[j, i] + B0[j, i])
if i+2 < col: # East
Ge[j, i] = G0[j, i+1] + 0.5*beta*(R0[j, i] - R0[j, i+2] + B0[j, i] - B0[j, i+2])
else:
if i+1 < col:
Ge[j, i] = G0[j, i+1] + beta*(R0[j, i] + B0[j, i])
else:
Ge[j, i] = beta*(R0[j, i] + B0[j, i])
if i-2 >= 0: # West
Gw[j, i] = G0[j, i-1] + 0.5*beta*(R0[j, i] - R0[j, i-2] + B0[j, i] - B0[j, i-2])
else:
if i-1 >= 0:
Gw[j, i] = G0[j, i-1] + beta*(R0[j, i] + B0[j, i])
else:
Gw[j, i] = beta*(R0[j, i] + B0[j, i])
# Bilinear interpolation of red and blue channels
if R0[j, i] != 0:
Rn[j, i], Rs[j, i], Re[j, i], Rw[j, i] = R0[j, i], R0[j, i], R0[j, i], R0[j, i]
RGn[j, i] = R0[j, i] - beta*Gn[j, i]
RGs[j, i] = R0[j, i] - beta*Gs[j, i]
RGe[j, i] = R0[j, i] - beta*Ge[j, i]
RGw[j, i] = R0[j, i] - beta*Gw[j, i]
if B0[j, i] != 0:
Bn[j, i], Bs[j, i], Be[j, i], Bw[j, i] = B0[j, i], B0[j, i], B0[j, i], B0[j, i]
BGn[j, i] = B0[j, i] - beta*Gn[j, i]
BGs[j, i] = B0[j, i] - beta*Gs[j, i]
BGe[j, i] = B0[j, i] - beta*Ge[j, i]
BGw[j, i] = B0[j, i] - beta*Gw[j, i]
if i+1 < col:
if G0[j, i] != 0 and R0[j, i+1] != 0: # Red square next to us
if i-1 >=0 :
Rn[j, i] = 0.5*(RGn[j, i-1] + RGn[j, i+1]) + beta*Gn[j, i]
Rs[j, i] = 0.5*(RGs[j, i-1] + RGs[j, i+1]) + beta*Gs[j, i]
Re[j, i] = 0.5*(RGe[j, i-1] + RGe[j, i+1]) + beta*Ge[j, i]
Rw[j, i] = 0.5*(RGw[j, i-1] + RGw[j, i+1]) + beta*Gw[j, i]
# The case nest to us is red, so it can't be blue
if j-1 >=0:
if j+1 < lin:
Bn[j, i] = 0.5*(BGn[j-1, i] + BGn[j+1, i]) + beta*Gn[j, i]
Bs[j, i] = 0.5*(BGs[j-1, i] + BGs[j+1, i]) + beta*Gs[j, i]
Be[j, i] = 0.5*(BGe[j-1, i] + BGe[j+1, i]) + beta*Ge[j, i]
Bw[j, i] = 0.5*(BGw[j-1, i] + BGw[j+1, i]) + beta*Gw[j, i]
else:
Bn[j, i] = BGn[j-1, i] + beta*Gn[j, i]
Bs[j, i] = BGs[j-1, i] + beta*Gs[j, i]
Be[j, i] = BGe[j-1, i] + beta*Ge[j, i]
Bw[j, i] = BGw[j-1, i] + beta*Gw[j, i]
else:
Bn[j, i] = BGn[j+1, i] + beta*Gn[j, i]
Bs[j, i] = BGs[j+1, i] + beta*Gs[j, i]
Be[j, i] = BGe[j+1, i] + beta*Ge[j, i]
Bw[j, i] = BGw[j+1, i] + beta*Gw[j, i]
if G0[j, i] != 0 and B0[j, i+1] != 0: # Blue square next to us
if i-1 >=0 :
Bn[j, i] = 0.5*(BGn[j, i-1] + BGn[j, i+1]) + beta*Gn[j, i]
Bs[j, i] = 0.5*(BGs[j, i-1] + BGs[j, i+1]) + beta*Gs[j, i]
Be[j, i] = 0.5*(BGe[j, i-1] + BGe[j, i+1]) + beta*Ge[j, i]
Bw[j, i] = 0.5*(BGw[j, i-1] + BGw[j, i+1]) + beta*Gw[j, i]
# The case nest to us is blue, so it can't be red
if j-1 >=0:
if j+1 < lin:
Rn[j, i] = 0.5*(RGn[j-1, i] + RGn[j+1, i]) + beta*Gn[j, i]
Rs[j, i] = 0.5*(RGs[j-1, i] + RGs[j+1, i]) + beta*Gs[j, i]
Re[j, i] = 0.5*(RGe[j-1, i] + RGe[j+1, i]) + beta*Ge[j, i]
Rw[j, i] = 0.5*(RGw[j-1, i] + RGw[j+1, i]) + beta*Gw[j, i]
else:
Rn[j, i] = RGn[j-1, i] + beta*Gn[j, i]
Rs[j, i] = RGs[j-1, i] + beta*Gs[j, i]
Re[j, i] = RGe[j-1, i] + beta*Ge[j, i]
Rw[j, i] = RGw[j-1, i] + beta*Gw[j, i]
if R0[j, i] == 0 and G0[j, i] == 0:
if i-1 >= 0:
if i+1 < col:
if j-1 >= 0:
if j+1 < lin:
Rn[j, i] = 0.25*(RGn[j-1, i-1] + RGn[j-1, i+1] + RGn[j+1, i-1] + RGn[j+1, i+1]) + beta*Gn[j, i]
Rs[j, i] = 0.25*(RGs[j-1, i-1] + RGs[j-1, i+1] + RGs[j+1, i-1] + RGs[j+1, i+1]) + beta*Gs[j, i]
Re[j, i] = 0.25*(RGe[j-1, i-1] + RGe[j-1, i+1] + RGe[j+1, i-1] + RGe[j+1, i+1]) + beta*Ge[j, i]
Rw[j, i] = 0.25*(RGw[j-1, i-1] + RGw[j-1, i+1] + RGw[j+1, i-1] + RGw[j+1, i+1]) + beta*Gw[j, i]
else:
Rn[j, i] = 0.5*(RGn[j-1, i-1] + RGn[j-1, i+1]) + beta*Gn[j, i]
Rs[j, i] = 0.5*(RGs[j-1, i-1] + RGs[j-1, i+1]) + beta*Gs[j, i]
Re[j, i] = 0.5*(RGe[j-1, i-1] + RGe[j-1, i+1]) + beta*Ge[j, i]
Rw[j, i] = 0.5*(RGw[j-1, i-1] + RGw[j-1, i+1]) + beta*Gw[j, i]
else:
Rn[j, i] = 0.5*(RGn[j+1, i-1] + RGn[j+1, i+1]) + beta*Gn[j, i]
Rs[j, i] = 0.5*(RGs[j+1, i-1] + RGs[j+1, i+1]) + beta*Gs[j, i]
Re[j, i] = 0.5*(RGe[j+1, i-1] + RGe[j+1, i+1]) + beta*Ge[j, i]
Rw[j, i] = 0.5*(RGw[j+1, i-1] + RGw[j+1, i+1]) + beta*Gw[j, i]
else:
if j-1 >= 0:
if j+1 < lin:
Rn[j, i] = 0.5*(RGn[j-1, i-1] + RGn[j+1, i-1]) + beta*Gn[j, i]
Rs[j, i] = 0.5*(RGs[j-1, i-1] + RGs[j+1, i-1]) + beta*Gs[j, i]
Re[j, i] = 0.5*(RGe[j-1, i-1] + RGe[j+1, i-1]) + beta*Ge[j, i]
Rw[j, i] = 0.5*(RGw[j-1, i-1] + RGw[j+1, i-1]) + beta*Gw[j, i]
else:
Rn[j, i] = RGn[j-1, i-1] + beta*Gn[j, i]
Rs[j, i] = RGs[j-1, i-1] + beta*Gs[j, i]
Re[j, i] = RGe[j-1, i-1] + beta*Ge[j, i]
Rw[j, i] = RGw[j-1, i-1] + beta*Gw[j, i]
else:
Rn[j, i] = beta*Gn[j, i]
Rs[j, i] = beta*Gs[j, i]
Re[j, i] = beta*Ge[j, i]
Rw[j, i] = beta*Gw[j, i]
else:
if j-1 >= 0:
if j+1 < lin:
Rn[j, i] = 0.5*(RGn[j-1, i+1] + RGn[j+1, i+1]) + beta*Gn[j, i]
Rs[j, i] = 0.5*(RGs[j-1, i+1] + RGs[j+1, i+1]) + beta*Gs[j, i]
Re[j, i] = 0.5*(RGe[j-1, i+1] + RGe[j+1, i+1]) + beta*Ge[j, i]
Rw[j, i] = 0.5*(RGw[j-1, i+1] + RGw[j+1, i+1]) + beta*Gw[j, i]
else:
Rn[j, i] = RGn[j-1, i+1] + beta*Gn[j, i]
Rs[j, i] = RGs[j-1, i+1] + beta*Gs[j, i]
Re[j, i] = RGe[j-1, i+1] + beta*Ge[j, i]
Rw[j, i] = RGw[j-1, i+1] + beta*Gw[j, i]
else:
Rn[j, i] = beta*Gn[j, i]
Rs[j, i] = beta*Gs[j, i]
Re[j, i] = beta*Ge[j, i]
Rw[j, i] = beta*Gw[j, i]
if B0[j, i] == 0 and G0[j, i] == 0:
if i-1 >= 0:
if i+1 < col:
if j-1 >= 0:
if j+1 < lin:
Bn[j, i] = 0.25*(BGn[j-1, i-1] + BGn[j-1, i+1] + BGn[j+1, i-1] + BGn[j+1, i+1]) + beta*Gn[j, i]
Bs[j, i] = 0.25*(BGs[j-1, i-1] + BGs[j-1, i+1] + BGs[j+1, i-1] + BGs[j+1, i+1]) + beta*Gs[j, i]
Be[j, i] = 0.25*(BGe[j-1, i-1] + BGe[j-1, i+1] + BGe[j+1, i-1] + BGe[j+1, i+1]) + beta*Ge[j, i]
Bw[j, i] = 0.25*(BGw[j-1, i-1] + BGw[j-1, i+1] + BGw[j+1, i-1] + BGw[j+1, i+1]) + beta*Gw[j, i]
else:
Bn[j, i] = 0.5*(BGn[j-1, i-1] + BGn[j-1, i+1]) + beta*Gn[j, i]
Bs[j, i] = 0.5*(BGs[j-1, i-1] + BGs[j-1, i+1]) + beta*Gs[j, i]
Be[j, i] = 0.5*(BGe[j-1, i-1] + BGe[j-1, i+1]) + beta*Ge[j, i]
Bw[j, i] = 0.5*(BGw[j-1, i-1] + BGw[j-1, i+1]) + beta*Gw[j, i]
else:
Bn[j, i] = 0.5*(BGn[j+1, i-1] + BGn[j+1, i+1]) + beta*Gn[j, i]
Bs[j, i] = 0.5*(BGs[j+1, i-1] + BGs[j+1, i+1]) + beta*Gs[j, i]
Be[j, i] = 0.5*(BGe[j+1, i-1] + BGe[j+1, i+1]) + beta*Ge[j, i]
Bw[j, i] = 0.5*(BGw[j+1, i-1] + BGw[j+1, i+1]) + beta*Gw[j, i]
else:
if j-1 >= 0:
if j+1 < lin:
Bn[j, i] = 0.5*(BGn[j-1, i-1] + BGn[j+1, i-1]) + beta*Gn[j, i]
Bs[j, i] = 0.5*(BGs[j-1, i-1] + BGs[j+1, i-1]) + beta*Gs[j, i]
Be[j, i] = 0.5*(BGe[j-1, i-1] + BGe[j+1, i-1]) + beta*Ge[j, i]
Bw[j, i] = 0.5*(BGw[j-1, i-1] + BGw[j+1, i-1]) + beta*Gw[j, i]
else:
Bn[j, i] = BGn[j-1, i-1] + beta*Gn[j, i]
Bs[j, i] = BGs[j-1, i-1] + beta*Gs[j, i]
Be[j, i] = BGe[j-1, i-1] + beta*Ge[j, i]
Bw[j, i] = BGw[j-1, i-1] + beta*Gw[j, i]
else:
Bn[j, i] = beta*Gn[j, i]
Bs[j, i] = beta*Gs[j, i]
Be[j, i] = beta*Ge[j, i]
Bw[j, i] = beta*Gw[j, i]
else:
if j-1 >= 0:
if j+1 < lin:
Bn[j, i] = 0.5*(BGn[j-1, i+1] + BGn[j+1, i+1]) + beta*Gn[j, i]
Bs[j, i] = 0.5*(BGs[j-1, i+1] + BGs[j+1, i+1]) + beta*Gs[j, i]
Be[j, i] = 0.5*(BGe[j-1, i+1] + BGe[j+1, i+1]) + beta*Ge[j, i]
Bw[j, i] = 0.5*(BGw[j-1, i+1] + BGw[j+1, i+1]) + beta*Gw[j, i]
else:
Bn[j, i] = BGn[j-1, i+1] + beta*Gn[j, i]
Bs[j, i] = BGs[j-1, i+1] + beta*Gs[j, i]
Be[j, i] = BGe[j-1, i+1] + beta*Ge[j, i]
Bw[j, i] = BGw[j-1, i+1] + beta*Gw[j, i]
else:
Bn[j, i] = beta*Gn[j, i]
Bs[j, i] = beta*Gs[j, i]
Be[j, i] = beta*Ge[j, i]
Bw[j, i] = beta*Gw[j, i]
# Pixel-level fusion of full color interpolated images
# Computation of the YUV space
_, Un, Vn = cvt_YUV_space(Rn, Gn, Bn)
_, Us, Vs = cvt_YUV_space(Rs, Gs, Bs)
_, Ue, Ve = cvt_YUV_space(Re, Ge, Be)
_, Uw, Vw = cvt_YUV_space(Rw, Gw, Bw)
# Computation of the gradients
Gradn[j, i] = compute_gradient(L, Un, Vn, i, j, "north")
Grads[j, i] = compute_gradient(L, Us, Vs, i, j, "south")
Grade[j, i] = compute_gradient(L, Ue, Ve, i, j, "east")
Gradw[j, i] = compute_gradient(L, Uw, Vw, i, j, "west")
low_wn[j, i] = 1 / (Gradn[j, i] + eps)
low_ws[j, i] = 1 / (Grads[j, i] + eps)
low_we[j, i] = 1 / (Grade[j, i] + eps)
low_ww[j, i] = 1 / (Gradw[j, i] + eps)
W[j, i] = low_wn[j, i] + low_ws[j, i] + low_we[j, i] + low_ww[j, i]
low_wn[j, i] = low_wn[j, i] / W[j, i]
low_ws[j, i] = low_ws[j, i] / W[j, i]
low_we[j, i] = low_we[j, i] / W[j, i]
low_ww[j, i] = low_ww[j, i] / W[j, i]
R[j, i] = low_wn[j, i]*Rn[j, i] + low_ws[j, i]*Rs[j, i] + low_we[j, i]*Re[j, i] + low_ww[j, i]*Rw[j, i]
G[j, i] = low_wn[j, i]*Gn[j, i] + low_ws[j, i]*Gs[j, i] + low_we[j, i]*Ge[j, i] + low_ww[j, i]*Gw[j, i]
B[j, i] = low_wn[j, i]*Bn[j, i] + low_ws[j, i]*Bs[j, i] + low_we[j, i]*Be[j, i] + low_ww[j, i]*Bw[j, i]
return R, G, B
\ No newline at end of file
File added
"""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.fortin_come.DemosaickingInterChannel import local_int
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)
z = op.adjoint(y)
R0 = z[:, :, 0]
G0 = z[:, :, 1]
B0 = z[:, :, 2]
beta = 0.9 # Inter-channel parameter
L = 5 # Size of the local neighborhood
eps = 0.00000001 # Tresholding parameter > 0
R1, G1, B1 = local_int(R0, G0, B0, beta, L, eps)
img_demosaicked = np.zeros((R1.shape[0], R1.shape[1], 3))
img_demosaicked[:, :, 0] = R1
img_demosaicked[:, :, 1] = G1
img_demosaicked[:, :, 2] = B1
return img_demosaicked
####
####
####
#### #### #### #############
#### ###### #### ##################
#### ######## #### ####################
#### ########## #### #### ########
#### ############ #### #### ####
#### #### ######## #### #### ####
#### #### ######## #### #### ####
#### #### ######## #### #### ####
#### #### ## ###### #### #### ######
#### #### #### ## #### #### ############
#### #### ###### #### #### ##########
#### #### ########## #### #### ########
#### #### ######## #### ####
#### #### ############ ####
#### #### ########## ####
#### #### ######## ####
#### #### ###### ####
# 2023
# Authors: FORTIN Côme
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment