diff --git a/src/methods/antoine_vouillon/functions.py b/src/methods/antoine_vouillon/functions.py
new file mode 100644
index 0000000000000000000000000000000000000000..d94b192e597840c7dfec1022f2a493cd07faa838
--- /dev/null
+++ b/src/methods/antoine_vouillon/functions.py
@@ -0,0 +1,203 @@
+import numpy as np
+
+def find_neighbors (z,chan,i,j,N,M):
+    
+    """Finds all the neighbors of a pixel on a given channel
+    Args:
+        z: image.
+        chan: chosen chanel.
+        i: line of the pixel
+        j: column of the pixel
+        N: number of lines
+        M: number of columns
+
+    Returns:
+        np.ndarray: Neighbors of the pixel in a list (including the pixel itself)
+    """
+
+    P1 = z[(i-1)%N,(j-1)%M,chan]
+    P2 = z[(i-1)%N,j%M,chan]
+    P3 = z[(i-1)%N,(j+1)%M,chan]
+    P4 = z[i%N,(j-1)%M,chan]
+    P5 = z[i%N,j%M,chan]
+    P6 = z[i%N,(j+1)%M,chan]
+    P7 = z[(i+1)%N,(j-1)%M,chan]
+    P8 = z[(i+1)%N,j%M,chan]
+    P9 = z[(i+1)%N,(j+1)%M,chan]
+
+    return np.array([P1,P2,P3,P4,P5,P6,P7,P8,P9])
+
+def find_dir_deriv(neighbors):
+    """Calculates the directional derivative of a pixel.
+
+    Args:
+        neighbors: list of the neighbors of the pixel
+
+    Returns:
+        np.ndarray: directional derivatives in this order: Dx, Dy, Ddx, Ddy
+    """
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neighbors
+    Dx = (P4 - P6)/2
+    Dy = (P2 - P8)/2
+    Dxd = (P3 - P7)/(2*np.sqrt(2))
+    Dyd = (P1 - P9)/(2*np.sqrt(2))
+
+    return [Dx,Dy,Dxd,Dyd]
+
+def find_weights(z, neigh, dir_deriv,chan,i,j,N,M):
+    """Finds all the neighbors of a pixel on a given channel
+    Args:
+        z: image.
+        dir_deriv: directional derivatives
+        chan: chosen chanel.
+        i: line of the pixel
+        j: column of the pixel
+        N: number of lines
+        M: number of columns
+
+    Returns:
+        np.ndarray: Weights from E1 to E9
+    """
+    [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_neighbors(z,chan,i+k,j+k,N,M)
+            dd = find_dir_deriv(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(neigh,weights):
+    
+    """interpolates pixels from a grid where one of two pixels is missing regularly spaced
+    Args:
+        neigh: neighbors of the pixel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the interpolated pixel
+    """
+
+    [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_RB(neigh, neigh_G, weights):
+
+    """interpolates the central missing pixel from the red or blue channel from a bayer patern
+    Args:
+        neigh: neighbors of the pixel in the red or blue channel.
+        neigh_G: neighbors of the pixel in the green channel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the interpolated pixel in the red or blue channel
+    """
+
+    [P1,P2,P3,P4,P5,P6,P7,P8,P9] = neigh
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neigh_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
+
+def correction_G(neigh_G ,neigh_R ,neigh_B, weights):
+
+    """corrects the value of the green pixel in the third phase of the kimmel algorythme
+    Args:
+
+        neigh_G: neighbors of the pixel in the green channel.
+        neigh_R: neighbors of the pixel in the red channel.
+        neigh_B: neighbors of the pixel in the blue channel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the corrected pixel in the green channel
+    """
+
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neigh_G
+    [R1,R2,R3,R4,R5,R6,R7,R8,R9] = neigh_R
+    [B1,B2,B3,B4,B5,B6,B7,B8,B9] = neigh_B
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+
+    num_Gb5 = ((E2*G2)/B2) + ((E4*G4)/B4) + ((E6*G6)/B6) + ((E8*G8)/B8)
+    num_Gr5 = ((E2*G2)/R2) + ((E4*G4)/R4) + ((E6*G6)/R6) + ((E8*G8)/R8)
+    den5 = E2 + E4 + E6 + E8
+
+    Gb5 = B5 * num_Gb5/den5
+    Gr5 = R5 * num_Gr5/den5
+
+    G5 = (Gb5 + Gr5)/2
+
+    return Gr5
+
+def correction_R(neigh_G, neigh_R, weights):
+
+    """corrects the value of the red pixel in the third phase of the kimmel algorythme
+    Args:
+
+        neigh_G: neighbors of the pixel in the green channel.
+        neigh_R: neighbors of the pixel in the red channel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the corrected pixel in the red channel
+    """
+
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neigh_G
+    [R1,R2,R3,R4,R5,R6,R7,R8,R9] = neigh_R
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+
+    num_R5 = ((E1*R1)/G1) + ((E2*R2)/G2) + ((E3*R3)/G3) + ((E4*R4)/G4) + ((E6*R6)/G6) + ((E7*R7)/G7) + ((E8*R8)/G8) + ((E9*R9)/G9)
+
+    den5 = sum(weights)
+
+    R5 = G5 * num_R5/den5
+
+    return R5
+
+def correction_B(neigh_G ,neigh_B, weights):
+
+    """corrects the value of the blue pixel in the third phase of the kimmel algorythme
+    Args:
+
+        neigh_G: neighbors of the pixel in the green channel.
+        neigh_B: neighbors of the pixel in the blue channel.
+        weights: weight of the neighbors.
+
+
+    Returns:
+        np.ndarray: The value of the corrected pixel in the blue channel
+    """
+
+    [G1,G2,G3,G4,G5,G6,G7,G8,G9] = neigh_G
+    [B1,B2,B3,B4,B5,B6,B7,B8,B9] = neigh_B
+    [E1,E2,E3,E4,E6,E7,E8,E9] = weights
+
+    num_B5 = ((E1*B1)/G1) + ((E2*B2)/G2) + ((E3*B3)/G3) + ((E4*B4)/G4) + ((E6*B6)/G6) + ((E7*B7)/G7) + ((E8*B8)/G8) + ((E9*B9)/G9)
+
+    den5 = sum(weights)
+
+    B5 = G5 * num_B5/den5
+
+    return B5
diff --git a/src/methods/antoine_vouillon/image_analysis_report_antoine_vouillon.pdf b/src/methods/antoine_vouillon/image_analysis_report_antoine_vouillon.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..e29075fe0512fc1252bbc95dc443a3769d2a0c32
Binary files /dev/null and b/src/methods/antoine_vouillon/image_analysis_report_antoine_vouillon.pdf differ
diff --git a/src/methods/antoine_vouillon/reconstruct.py b/src/methods/antoine_vouillon/reconstruct.py
new file mode 100644
index 0000000000000000000000000000000000000000..7400260660831d17f023cc82d0210b43cf5e69d8
--- /dev/null
+++ b/src/methods/antoine_vouillon/reconstruct.py
@@ -0,0 +1,136 @@
+"""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.antoine_vouillon.functions import *
+
+#!!!!!!!! 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)
+
+    img_res = op.adjoint(y)
+    N = img_res[:,:,0].shape[0]
+    M = img_res[:,:,0].shape[1]
+
+    #interpolating green channel
+
+    for i in range (N):
+        for j in range (M):
+            if img_res[i,j,1] ==0:
+
+                neighbors = find_neighbors(img_res,1,i,j,N,M)
+                dir_deriv = find_dir_deriv(neighbors)
+                weights = find_weights(img_res, neighbors, dir_deriv,1,i,j,N,M)
+
+                img_res[i,j,1] = interpolate(neighbors,weights)
+                
+    
+    img_res[img_res>1] = 1
+    img_res[img_res<0] = 0
+
+    #first intepolation of red channel
+
+    for i in range (1,N,2):
+        for j in range (0,M,2):
+            
+            neighbors = find_neighbors(img_res,0,i,j,N,M)
+
+            neighbors_G = find_neighbors(img_res,1,i,j,N,M)
+            dir_deriv = find_dir_deriv(neighbors_G)
+            weights = find_weights(img_res,neighbors_G, dir_deriv,1,i,j,N,M) 
+
+
+
+            img_res[i,j,0] = interpolate_RB(neighbors, neighbors_G, weights)
+
+    # second interpolation of red channel
+
+    for i in range (N):
+        for j in range (M):
+            if img_res[i,j,0] ==0:
+
+                neighbors = find_neighbors(img_res,0,i,j,N,M)
+                dir_deriv = find_dir_deriv(neighbors)
+                weights = find_weights(img_res,neighbors, dir_deriv,0,i,j,N,M)
+
+                img_res[i,j,0] = interpolate(neighbors,weights)
+
+    img_res[img_res>1] = 1
+    img_res[img_res<0] = 0
+
+    #first interpolation of blue channel
+
+    for i in range (0,N,2):
+        for j in range (1,M,2):
+            
+            neighbors = find_neighbors(img_res,2,i,j,N,M)
+
+            neighbors_G = find_neighbors(img_res,1,i,j,N,M)
+            dir_deriv = find_dir_deriv(neighbors_G)
+            weights = find_weights(img_res,neighbors_G, dir_deriv,1,i,j,N,M) 
+
+
+
+            img_res[i,j,2] = interpolate_RB(neighbors, neighbors_G, weights)
+
+    #second interpolation of blue channel
+
+    for i in range (N):
+        for j in range (M):
+            if img_res[i,j,2] ==0:
+
+                neighbors = find_neighbors(img_res,2,i,j,N,M)
+                dir_deriv = find_dir_deriv(neighbors)
+                weights = find_weights(img_res, neighbors, dir_deriv,2,i,j,N,M)
+
+                img_res[i,j,2] = interpolate(neighbors,weights)
+
+    img_res[img_res>1] = 1
+    img_res[img_res<0] = 0
+
+    return img_res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller
diff --git a/src/methods/template/reconstruct.py b/src/methods/template/reconstruct.py
deleted file mode 100755
index a97bd3f6e3c68df763b36c46b2727461af078bd2..0000000000000000000000000000000000000000
--- a/src/methods/template/reconstruct.py
+++ /dev/null
@@ -1,53 +0,0 @@
-"""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
-
-
-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
-    input_shape = (y.shape[0], y.shape[1], 3)
-    op = CFA(cfa, input_shape)
-
-    return np.zeros(op.input_shape)
-
-
-####
-####
-####
-
-####      ####                ####        #############
-####      ######              ####      ##################
-####      ########            ####      ####################
-####      ##########          ####      ####        ########
-####      ############        ####      ####            ####
-####      ####  ########      ####      ####            ####
-####      ####    ########    ####      ####            ####
-####      ####      ########  ####      ####            ####
-####      ####  ##    ######  ####      ####          ######
-####      ####  ####      ##  ####      ####    ############
-####      ####  ######        ####      ####    ##########
-####      ####  ##########    ####      ####    ########
-####      ####      ########  ####      ####
-####      ####        ############      ####
-####      ####          ##########      ####
-####      ####            ########      ####
-####      ####              ######      ####
-
-# 2023
-# Authors: Mauro Dalla Mura and Matthieu Muller