src/methods/domer/image-colorization/network.ipynb

+%% Cell type:code id: tags:
+
+``` 
+from network import ColorizeNet
+
+model = ColorizeNet()
+model
+```
+
+%% Output
+
+    ColorizeNet(
+      (encoder): Sequential(
+        (0): Conv2d(3, 64, kernel_size=(7, 7), stride=(2, 2), padding=(3, 3), bias=False)
+        (1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+        (2): ReLU(inplace=True)
+        (3): MaxPool2d(kernel_size=3, stride=2, padding=1, dilation=1, ceil_mode=False)
+        (4): Sequential(
+          (0): BasicBlock(
+            (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+          )
+          (1): BasicBlock(
+            (conv1): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+          )
+        )
+        (5): Sequential(
+          (0): BasicBlock(
+            (conv1): Conv2d(64, 128, kernel_size=(3, 3), stride=(2, 2), padding=(1, 1), bias=False)
+            (bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (downsample): Sequential(
+              (0): Conv2d(64, 128, kernel_size=(1, 1), stride=(2, 2), bias=False)
+              (1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            )
+          )
+          (1): BasicBlock(
+            (conv1): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn1): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(128, 128, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(128, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+          )
+        )
+      )
+      (decoder): Sequential(
+        (0): Sequential(
+          (0): BasicBlock(
+            (conv1): Conv2d(128, 64, kernel_size=(5, 5), stride=(1, 1), padding=(2, 2), bias=False)
+            (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (res_conv): Conv2d(128, 64, kernel_size=(1, 1), stride=(1, 1), bias=False)
+            (res_bn): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (upsample): Upsample(scale_factor=2.0, mode=nearest)
+          )
+          (1): BasicBlock(
+            (conv1): Conv2d(64, 64, kernel_size=(5, 5), stride=(1, 1), padding=(2, 2), bias=False)
+            (bn1): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(64, 64, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(64, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (activation): ReLU(inplace=True)
+          )
+        )
+        (1): Sequential(
+          (0): BasicBlock(
+            (conv1): Conv2d(64, 32, kernel_size=(5, 5), stride=(1, 1), padding=(2, 2), bias=False)
+            (bn1): BatchNorm2d(32, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(32, 32, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(32, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (res_conv): Conv2d(64, 32, kernel_size=(1, 1), stride=(1, 1), bias=False)
+            (res_bn): BatchNorm2d(32, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (upsample): Upsample(scale_factor=2.0, mode=nearest)
+          )
+          (1): BasicBlock(
+            (conv1): Conv2d(32, 32, kernel_size=(5, 5), stride=(1, 1), padding=(2, 2), bias=False)
+            (bn1): BatchNorm2d(32, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(32, 32, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(32, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (activation): ReLU(inplace=True)
+          )
+        )
+        (2): Sequential(
+          (0): BasicBlock(
+            (conv1): Conv2d(32, 2, kernel_size=(5, 5), stride=(1, 1), padding=(2, 2), bias=False)
+            (bn1): BatchNorm2d(2, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(2, 2, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(2, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (res_conv): Conv2d(32, 2, kernel_size=(1, 1), stride=(1, 1), bias=False)
+            (res_bn): BatchNorm2d(2, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (upsample): Upsample(scale_factor=2.0, mode=nearest)
+          )
+          (1): BasicBlock(
+            (conv1): Conv2d(2, 2, kernel_size=(5, 5), stride=(1, 1), padding=(2, 2), bias=False)
+            (bn1): BatchNorm2d(2, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (relu): ReLU(inplace=True)
+            (conv2): Conv2d(2, 2, kernel_size=(3, 3), stride=(1, 1), padding=(1, 1), bias=False)
+            (bn2): BatchNorm2d(2, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)
+            (activation): Sigmoid()
+          )
+        )
+      )
+    )
+
+%% Cell type:code id: tags:
+
+``` 
+import torch
+
+# to check if our model is working as expected
+model(torch.rand((2, 1, 224, 224))).shape
+```
+
+%% Output
+
+    torch.Size([2, 2, 224, 224])
+
+%% Cell type:code id: tags:
+
+``` 
+from utils import count_params
+
+count_params(model)  # no of trainable parameters
+```
+
+%% Output
+
+    1166016

src/methods/domer/image-colorization/network.py

+import torch.nn as nn
+import torchvision.models as models
+
+
+class BasicBlock(nn.Module):
+    def __init__(self, in_channels, out_channels,
+                 activation=None, upsample=None):
+        super().__init__()
+
+        self.conv1 = nn.Conv2d(
+            in_channels=in_channels, out_channels=out_channels,
+            kernel_size=5, stride=1, padding=2, bias=False
+        )
+        self.bn1 = nn.BatchNorm2d(num_features=out_channels)
+        self.relu = nn.ReLU(inplace=True)
+        self.conv2 = nn.Conv2d(
+            in_channels=out_channels, out_channels=out_channels,
+            kernel_size=3, stride=1, padding=1, bias=False
+        )
+        self.bn2 = nn.BatchNorm2d(num_features=out_channels)
+
+        if activation is not None:
+            self.activation = activation
+        else:
+            self.res_conv = nn.Conv2d(
+                in_channels=in_channels, out_channels=out_channels,
+                kernel_size=1, stride=1, bias=False
+            )
+            self.res_bn = nn.BatchNorm2d(num_features=out_channels)
+
+        self.upsample = upsample
+
+    def forward(self, t):
+
+        res = t
+        t = self.conv1(t)
+        t = self.bn1(t)
+        t = self.relu(t)
+        t = self.conv2(t)
+        t = self.bn2(t)
+
+        if self.upsample is not None:
+            res = self.res_conv(res)
+            res = self.res_bn(res)
+            t += res
+            t = self.relu(t)
+            t = self.upsample(t)
+        else:
+            t += res
+            t = self.activation(t)
+
+        return t
+
+
+class ColorizeNet(nn.Module):
+    def __init__(self):
+        super().__init__()
+
+        # make pretrained=True before starting the training process
+        resnet18 = models.resnet18(pretrained=False)
+        # change first conv layer to accept single channel (grayscale)
+        resnet18.conv1.weight = nn.Parameter(
+            resnet18.conv1.weight.mean(dim=1).unsqueeze(dim=1))
+
+        # use first 3 layers of ResNet-18 as encoder
+        self.encoder = nn.Sequential(
+            *list(resnet18.children())[:6]
+        )
+        self.decoder = nn.Sequential(
+            self._make_layer(BasicBlock, 128, 64, nn.ReLU(inplace=True)),
+            self._make_layer(BasicBlock, 64, 32, nn.ReLU(inplace=True)),
+            self._make_layer(BasicBlock, 32, 2, nn.Sigmoid())
+        )
+
+    def _make_layer(self, block, in_channels, out_channels, activation):
+        upsample = nn.Upsample(scale_factor=2, mode='nearest')
+        layers = []
+        layers.append(block(in_channels, out_channels, upsample=upsample))
+        layers.append(block(out_channels, out_channels, activation=activation))
+
+        return nn.Sequential(*layers)
+
+    def forward(self, t):
+        t = self.encoder(t)
+        t = self.decoder(t)
+
+        return t

src/methods/domer/image-colorization/requirements.txt

+matplotlib==3.3.2
+numpy==1.19.4
+Pillow==8.0.1
+scikit-image==0.17.2
+torch==1.7.0+cpu
+torchvision==0.8.1+cpu
+tqdm==4.51.0
\ No newline at end of file

src/methods/domer/image-colorization/train.ipynb

+%% Cell type:code id: tags:
+
+``` 
+# from google.colab import drive
+# drive.mount('/content/drive')
+```
+
+%% Output
+
+    Mounted at /content/drive
+
+%% Cell type:code id: tags:
+
+``` 
+# Download and unzip (2.2GB)
+# !wget http://data.csail.mit.edu/places/places205/testSetPlaces205_resize.tar.gz
+# !mkdir images/
+# !tar -xzf drive/MyDrive/testSetPlaces205_resize.tar.gz -C 'images/'
+```
+
+%% Cell type:code id: tags:
+
+``` 
+import torch
+import torch.nn as nn
+import torch.optim as optim
+import torchvision.transforms as transforms
+
+from tqdm import tqdm
+from network import ColorizeNet
+from utils import GrayscaleImageFolder
+```
+
+%% Cell type:code id: tags:
+
+``` 
+device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')
+```
+
+%% Cell type:code id: tags:
+
+``` 
+transform = transforms.Compose([
+    transforms.RandomResizedCrop(224),
+    transforms.RandomHorizontalFlip()
+])
+
+train_set = GrayscaleImageFolder('images/', transform=transform)
+train_loader = torch.utils.data.DataLoader(train_set, batch_size=128, shuffle=True, num_workers=4)
+```
+
+%% Cell type:code id: tags:
+
+``` 
+criterion = nn.MSELoss()
+
+model = ColorizeNet().to(device)
+optimizer = optim.Adam(model.parameters(), lr=0.001)
+# scheduler = optim.lr_scheduler.MultiStepLR(optimizer, milestones=[8, 24], gamma=1/3)
+scheduler = optim.lr_scheduler.ReduceLROnPlateau(optimizer, factor=1/3, patience=5, verbose=True)
+```
+
+%% Cell type:code id: tags:
+
+``` 
+try:
+    checkpoint = torch.load('drive/MyDrive/checkpoint_0.001.pth', map_location=device)
+    start_epoch = checkpoint['next_epoch']
+    model.load_state_dict(checkpoint['model'])
+    optimizer.load_state_dict(checkpoint['optimizer'])
+    scheduler.load_state_dict(checkpoint['scheduler'])
+
+except:
+    start_epoch = 0
+```
+
+%% Cell type:code id: tags:
+
+``` 
+num_epochs = 64
+for epoch in range(start_epoch, num_epochs):
+    train_loss = 0
+    loop = tqdm(train_loader)
+    for batch in loop:
+        in_gray, in_ab = batch[0].to(device), batch[1].to(device)
+        out_ab = model(in_gray)
+        loss = criterion(out_ab, in_ab)
+        optimizer.zero_grad()
+        loss.backward()
+        optimizer.step()
+
+        train_loss += loss.item()*in_gray.size(0)
+        loop.set_description(f'Epoch [{epoch+1:2d}/{num_epochs}]')
+        loop.set_postfix(loss=train_loss)
+
+    scheduler.step(train_loss)
+    checkpoint = {
+        'next_epoch': epoch + 1,
+        'model': model.state_dict(),
+        'optimizer': optimizer.state_dict(),
+        'scheduler': scheduler.state_dict()
+    }
+    torch.save(checkpoint, f'drive/MyDrive/checkpoint_{scheduler._last_lr[0]}.pth')
+```
+
+%% Output
+
+    Epoch [39/64]: 100%|██████████| 321/321 [12:00<00:00,  2.24s/it, loss=97.4]
+    Epoch [40/64]: 100%|██████████| 321/321 [11:30<00:00,  2.15s/it, loss=98]
+    Epoch [41/64]: 100%|██████████| 321/321 [11:18<00:00,  2.12s/it, loss=97.4]
+    Epoch [42/64]: 100%|██████████| 321/321 [11:28<00:00,  2.15s/it, loss=103]
+    Epoch [43/64]: 100%|██████████| 321/321 [11:23<00:00,  2.13s/it, loss=98.4]
+    Epoch [44/64]: 100%|██████████| 321/321 [11:15<00:00,  2.10s/it, loss=97.3]
+    Epoch [45/64]: 100%|██████████| 321/321 [11:29<00:00,  2.15s/it, loss=97.6]
+    Epoch [46/64]: 100%|██████████| 321/321 [11:25<00:00,  2.14s/it, loss=96.5]
+    Epoch [47/64]: 100%|██████████| 321/321 [11:40<00:00,  2.18s/it, loss=96.7]
+    Epoch [48/64]: 100%|██████████| 321/321 [11:44<00:00,  2.19s/it, loss=96.3]
+    Epoch [49/64]: 100%|██████████| 321/321 [11:13<00:00,  2.10s/it, loss=96.5]
+    Epoch [50/64]: 100%|██████████| 321/321 [11:09<00:00,  2.09s/it, loss=96.4]
+    Epoch [51/64]: 100%|██████████| 321/321 [11:09<00:00,  2.09s/it, loss=95.7]
+    Epoch [52/64]: 100%|██████████| 321/321 [11:13<00:00,  2.10s/it, loss=95.9]
+    Epoch [53/64]: 100%|██████████| 321/321 [11:08<00:00,  2.08s/it, loss=95.8]
+    Epoch [54/64]: 100%|██████████| 321/321 [11:31<00:00,  2.15s/it, loss=95.4]
+    Epoch [55/64]: 100%|██████████| 321/321 [11:54<00:00,  2.23s/it, loss=95.5]
+    Epoch [56/64]: 100%|██████████| 321/321 [12:11<00:00,  2.28s/it, loss=95.5]
+    Epoch [57/64]: 100%|██████████| 321/321 [12:00<00:00,  2.24s/it, loss=94.9]
+    Epoch [58/64]: 100%|██████████| 321/321 [11:54<00:00,  2.23s/it, loss=94.9]
+    Epoch [59/64]: 100%|██████████| 321/321 [11:56<00:00,  2.23s/it, loss=94.9]
+    Epoch [60/64]: 100%|██████████| 321/321 [12:06<00:00,  2.26s/it, loss=94.3]
+    Epoch [61/64]: 100%|██████████| 321/321 [12:03<00:00,  2.26s/it, loss=95]
+    Epoch [62/64]: 100%|██████████| 321/321 [12:01<00:00,  2.25s/it, loss=95]
+    Epoch [63/64]: 100%|██████████| 321/321 [12:05<00:00,  2.26s/it, loss=96.7]
+    Epoch [64/64]: 100%|██████████| 321/321 [12:10<00:00,  2.27s/it, loss=95]
+
+%% Cell type:code id: tags:
+
+``` 
+torch.save(model.state_dict(), f'./models/model.pth')  # save the model separately from the checkpoint file
+```

src/methods/domer/image-colorization/utils.py

+import torch
+from torchvision import transforms, datasets
+import numpy as np
+from PIL import Image
+from skimage.color import rgb2lab, rgb2gray, lab2rgb
+
+
+def count_params(model):
+    '''
+    returns the number of trainable parameters in some model
+    '''
+    return sum(p.numel() for p in model.parameters() if p.requires_grad)
+
+
+class GrayscaleImageFolder(datasets.ImageFolder):
+    '''
+    Custom dataloader for various operations on images before loading them.
+    '''
+
+    def __getitem__(self, index):
+        path, target = self.imgs[index]
+        img = self.loader(path)
+        if self.transform is not None:
+            img_orig = self.transform(img)  # apply transforms
+            img_orig = np.asarray(img_orig)  # convert to numpy array
+            img_lab = rgb2lab(img_orig)  # convert RGB image to LAB
+            img_ab = img_lab[:, :, 1:3]  # separate AB channels from LAB
+            img_ab = (img_ab + 128) / 255  # normalize the pixel values
+            # transpose image from HxWxC to CxHxW and turn it into a tensor
+            img_ab = torch.from_numpy(img_ab.transpose((2, 0, 1))).float()
+            img_orig = rgb2gray(img_orig)  # convert RGB to grayscale
+            # add a channel axis to grascale image and turn it into a tensor
+            img_orig = torch.from_numpy(img_orig).unsqueeze(0).float()
+
+        if self.target_transform is not None:
+            target = self.target_transform(target)
+
+        return img_orig, img_ab, target
+
+
+def load_gray(path, max_size=360, shape=None):
+    '''
+    load an image as grayscale, change the shape as per input,
+    perform transformations and convert it to model compatable shape.
+    '''
+    img_gray = Image.open(path).convert('L')
+
+    if max(img_gray.size) > max_size:
+        size = max_size
+    else:
+        size = max(img_gray.size)
+
+    if shape is not None:
+        size = shape
+
+    img_transform = transforms.Compose([
+        transforms.Resize(size),
+        transforms.ToTensor()
+    ])
+
+    img_gray = img_transform(img_gray).unsqueeze(0)
+    return img_gray
+
+
+def to_rgb(img_l, img_ab):
+    '''
+    concatinates Lightness (grayscale) and AB channels,
+    and converts the resulting LAB image to RGB
+    '''
+    if img_l.shape == img_ab.shape:
+        img_lab = torch.cat((img_l, img_ab), 1).numpy().squeeze()
+    else:
+        img_lab = torch.cat(
+            (img_l, img_ab[:, :, :img_l.size(2), :img_l.size(3)]),
+            dim=1
+        ).numpy().squeeze()
+
+    img_lab = img_lab.transpose(1, 2, 0)  # transpose image to HxWxC
+    img_lab[:, :, 0] = img_lab[:, :, 0] * 100  # range pixel values from 0-100
+    img_lab[:, :, 1:] = img_lab[:, :, 1:] * 255 - 128  # un-normalize
+    img_rgb = lab2rgb(img_lab.astype(np.float64))  # convert LAB image to RGB
+
+    return img_rgb

src/methods/domer/reconstruct.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Jan 30 19:35:12 2024
+
+@author: stani1
+"""
+
+"""The main file for the baseline reconstruction.
+This file should NOT be modified.
+"""
+
+
+import numpy as np
+from PIL import Image
+import subprocess
+import os
+
+
+from src.forward_model import CFA
+from src.methods.baseline.demo_reconstruction import naive_interpolation
+
+
+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.
+    """
+    current_directory = os.getcwd() #we take the current directory to be sure we will have the correct one
+    input_path='src/methods/domer/image-colorization/input/0.jpg'
+    input_path = os.path.join(current_directory, input_path)#location to store the image
+    output_path='src/methods/domer/image-colorization/output/0.jpg'
+    output_path = os.path.join(current_directory, output_path)
+    # venv_path='src/methods/domer/image-colorization/image-colorization/bin/activate'
+    # venv_path = os.path.join(current_directory, venv_path) #activation of the venv, but no need for a lot of versions of libraries
+    y2=Image.fromarray(y)
+    y2=y2.convert('L')
+    y2.save(input_path)
+    resolution=np.shape(y)[0]
+    # activation_venv = ['source ', venv_path]
+    # subprocess.run(activation_venv, shell=True)
+    n=y2.size[0]
+    # for i in range (n) :
+    #     for j in range (n) :
+    #         y2.putpixel((i,j),y.getpixel((i,j))*255)
+    command = [
+        'python', 'src/methods/domer/image-colorization/colorize.py',
+        '-i', input_path,
+        '-o', output_path,
+        '-r', str(resolution)
+    ]
+    subprocess.run(command) #calling of the NN with the location of the input and output and the resolution
+    # deactivation_venv = f"deactivate"
+    # subprocess.run(deactivation_venv, shell=True)
+    res=Image.open(output_path) #We take the result of the NN
+    res=np.asarray(res)#reseting the result in a numpy array
+
+    return res
+
+

src/methods/khattabi/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.khattabi.utils import HA
+
+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)
+    y_adjoint=op.adjoint(y)
+    res = HA(y_adjoint)
+
+    return res
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/khattabi/utils.py

+"""A file containing a (pretty useless) reconstruction.
+It serves as example of how the project works.
+This file should NOT be modified.
+"""
+
+
+import numpy as np
+from scipy.signal import convolve2d
+
+from src.forward_model import CFA
+
+def HA(y_adjoint):
+    height, width = y_adjoint.shape[:2]
+    # Green channel interpolation
+
+
+    for i in range(2 , height-2): 
+        for j in range(2, width-2):
+            # Red pixels
+            if y_adjoint[i,j,1] == 0 and y_adjoint[i,j,0] != 0:  
+                grad_y = np.abs(y_adjoint[i-1,j,1] - y_adjoint[i+1,j,1]) + np.abs(2*y_adjoint[i,j,0] - y_adjoint[i-2,j,0] - y_adjoint[i+2,j,0])
+                grad_x = np.abs(y_adjoint[i,j-1,1] - y_adjoint[i,j+1,1]) + np.abs(2*y_adjoint[i,j,0] - y_adjoint[i,j-2,0] - y_adjoint[i,j+2,0])
+
+                if grad_x < grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/2 + (2*y_adjoint[i,j,0] - y_adjoint[i,j-2,0] - y_adjoint[i,j+2,0])/4
+                elif  grad_x > grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1])/2 + (2*y_adjoint[i,j,0] - y_adjoint[i-2,j,0] - y_adjoint[i+2,j,0])/4
+                else:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1] + y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/4 + \
+                                       (2*y_adjoint[i,j,0] - y_adjoint[i,j-2,0] - y_adjoint[i,j+2,0] + 2*y_adjoint[i,j,0] - y_adjoint[i-2,j,0] - y_adjoint[i+2,j,0])/8
+            # Blue Pixels
+            elif y_adjoint[i,j,1] == 0 and y_adjoint[i,j,2] != 0:   
+                grad_y = np.abs(y_adjoint[i-1,j,1] - y_adjoint[i+1,j,1]) + np.abs(2*y_adjoint[i,j,2] - y_adjoint[i-2,j,2] - y_adjoint[i+2,j,2])
+                grad_x = np.abs(y_adjoint[i,j-1,1] - y_adjoint[i,j+1,1]) + np.abs(2*y_adjoint[i,j,2] - y_adjoint[i,j-2,2] - y_adjoint[i,j+2,2])
+
+                if grad_x < grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/2 + (2*y_adjoint[i,j,2] - y_adjoint[i,j-2,2] - y_adjoint[i,j+2,2])/4
+                elif  grad_x > grad_y:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1])/2 + (2*y_adjoint[i,j,2] - y_adjoint[i-2,j,2] - y_adjoint[i+2,j,2])/4
+                else:
+                    y_adjoint[i,j,1] = (y_adjoint[i-1,j,1] + y_adjoint[i+1,j,1] + y_adjoint[i,j-1,1] + y_adjoint[i,j+1,1])/4 + \
+                                       (2*y_adjoint[i,j,2] - y_adjoint[i,j-2,2] - y_adjoint[i,j+2,2] + 2*y_adjoint[i,j,2] - y_adjoint[i-2,j,2] - y_adjoint[i+2,j,2])/8
+                    
+    # Red and blue channel interpolation
+    for i in range(1 , height-1):  
+        for j in range(1, width-1):
+            if y_adjoint[i,j,2] != 0 :
+                 y_adjoint[i,j,0] = (y_adjoint[i-1,j-1,0] + y_adjoint[i-1,j+1,0] + y_adjoint[i+1,j-1,0] + y_adjoint[i+1,j+1,0]) / 4
+
+            elif y_adjoint[i,j,0] != 0:
+                 y_adjoint[i,j,2] = (y_adjoint[i-1,j-1,2] + y_adjoint[i-1,j+1,2] + y_adjoint[i+1,j-1,2] + y_adjoint[i+1,j+1,2]) / 4
+            
+            else:
+                y_adjoint[i,j] = y_adjoint[i,j]
+
+    for i in range(1 , height-1):
+        for j in range(1, width-1):    
+            if y_adjoint[i,j,0] == y_adjoint[i,j,2]:
+                y_adjoint[i,j,0] = (y_adjoint[i-1,j,0] + y_adjoint[i,j-1,0] + y_adjoint[i+1,j,0] + y_adjoint[i,j+1,0]) / 4             
+                y_adjoint[i,j,2] = (y_adjoint[i-1,j,2] + y_adjoint[i,j-1,2] + y_adjoint[i+1,j,2] + y_adjoint[i,j+1,2]) / 4  
+            
+    return y_adjoint
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller

src/methods/laporte_auriane/fct_to_demosaick.py

+"""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

src/methods/laporte_auriane/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.laporte_auriane.fct_to_demosaick import *
+
+
+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.
+    """
+    
+    #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 z
+
+
+####
+####
+####
+
+####      ####                ####        #############
+####      ######              ####      ##################
+####      ########            ####      ####################
+####      ##########          ####      ####        ########
+####      ############        ####      ####            ####
+####      ####  ########      ####      ####            ####
+####      ####    ########    ####      ####            ####
+####      ####      ########  ####      ####            ####
+####      ####  ##    ######  ####      ####          ######
+####      ####  ####      ##  ####      ####    ############
+####      ####  ######        ####      ####    ##########
+####      ####  ##########    ####      ####    ########
+####      ####      ########  ####      ####
+####      ####        ############      ####
+####      ####          ##########      ####
+####      ####            ########      ####
+####      ####              ######      ####
+
+# 2023
+# Authors: Mauro Dalla Mura and Matthieu Muller