Keywords: Image processing, Signal processing, Optimization, Sparse Methods, Superresolution, Inverse problems, Plug-and-play, Proximal operator, Image deconvolution, Image denoising, Herschel observation maps, Variable splitting, Potential function Image deconvolution using Approximate Proximal Operators

In the framework of the acquisition chain and devices built by i2S, the objective of this PhD is to provide efficient algorithms able to merge different acquired images corresponding to slight spatial translational displacements to get a wider super-resolved image. Consequently, this PhD takes place within the general subject of super-resolution.

In this thesis, super-resolution is performed using different acquisitions. We propose the use of sensor displacement within the cameras of I2S. We propose a scheme to achieve an image with up to two times higher resolution using this technique. Furthermore, we also propose an additional image deconvolution algorithm that helps to improve the image quality further and to address any degradation problem that may occur through the super-resolution scheme. Our image deconvolution algorithm is based on variable splitting and takes advantage of the proximal operator and Fourier transform. We also proposed the use of new potential functions that can be used as prior information in image inverse problems for the first time used in image processing. Experimental results show promising capabilities of the proposed algorithm. The algorithm is successfully implemented within various cameras and devices of I2S. The practical experiments on real-world data prove the effectiveness and flexibility of our image-deconvolution. Experiments were conducted on Herschel observation maps, and promising results were obtained on such imaging data. In the last part of the thesis, the idea of plug-and-play priors for image denoising and deconvolution is presented. This thesis proposes the implementation of plug-and play-priors in an alternating minimisation scheme. The early result has shown potential to be adequate for image denoising/deconvolution application.

Résumé

Dans le cadre de la chaîne d'acquisition et des dispositifs construits par I2S, l'objectif de cette thèse est de fournir des algorithmes efficaces capables de fusionner différentes images acquises correspondant à de légers déplacements de translation spatiale pour obtenir une image super-résolue plus large. Par conséquent, cette thèse s'inscrit dans le cadre général de la super-résolution. Dans cette thèse, la super-résolution est réalisée à l'aide de différentes acquisitions. Nous proposons l'utilisation du capteur de déplacement au sein des caméras d'I2S. Nous proposons un schéma pour obtenir une image avec une résolution jusqu'à deux fois plus élevée en utilisant cette technique. En outre, nous proposons également un algorithme de déconvolution d'image supplémentaire qui aide à améliorer davantage la qualité de l'image et à résoudre les problémes de dégradation pouvant survenir via le schéma de super-résolution. Notre algorithme de déconvolution d'image est basé sur le fractionnement variable et tire parti de l'opérateur proximal et de la transformée de Fourier. Nous introduisons également de nouvelles fonctions potentielles qui peuvent être utilisées comme information préalable dans les problémes inverses qui sont utilisés pour la premiére fois en traitement d'image. Les résultats expérimentaux montrent des capacités prometteuses de l'algorithme proposé.

L'algorithme est mis en oeuvre avec succès dans diverses caméras et appareils d'I2S. Les expériences pratiques sur des données du monde réel prouvent l'efficacité et la flexibilité de notre méthode de déconvolution. Des expériences ont été menées sur des cartes d'observation de Herschel, et des résultats prometteurs ont été obtenus sur de telles données d'imagerie. Dans la dernière partie de la thèse, l'idée de priors plug-and-play pour le débruitage et la déconvolution des images est présentée. Cette thése propose l'implémentation de priors plug-and-play dans un schéma de minimisation alternée. Le premier résultat a montré un potentiel adéquat pour l'application de débruitage/déconvolution d'image.

I am extremely grateful to my supervisor, Dr Hussein Yahia, for his invaluable advice, continuous support, and patience during my PhD study. His immense knowledge and great experience have encouraged me in my academic research and daily life, and I am indebted to his efforts and kindness.

I would also like to thank Dr Alexis Cailly for his technical support at I2S. He has always been a great friend and role model and guided me throughout this work. Moreover, I would like to thank Dr Anne Cherif for her encouraging words, and thoughtful and detailed feedback has been very important to me. You have always cared for my personal and professional life.

I would like to express my gratitude to all the jury members who honoured me with their presence and feedback: Dr Jérôme Idier, Chercheur au CNRS (Laboratoire des Sciences du Numérique de Nantes); Professeur Stéphane Binczak (Université de Bourgogne); Dr Pierrick Legrand, Maître de Conférences, Université de Bordeaux; Professeur Sabir Jacquir (Université Paris-Saclay).

Thank you to my parents, Ali and Shiva, for your endless support and love. You have always stood behind me, and this was no exception. Thank you for all of your love and for constantly reminding me of the end goal.

I would like to thank Minoo for her love and constant support, for all the late nights and early mornings, and for keeping me sane over the past few months. Who has stood by me through all my absences, my fits of pique and impatience. She gave me support and help, discussed ideas and prevented several wrong turns.

I would like to express my gratitude to Xavier Datin, CEO of I2S and Christophe Lacroix, Director of innovation at I2S. Thank you for placing your trust and confidence in my abilities and allowing me to pursue this PhD at I2S and INRIA.

From the INRIA side, I am incredibly thankful for Sabrina and Chris for their kindness in supporting me. They have become part of my family and great friends who have always kindly offered their support. I thank my great colleagues and friends, Michael and Delphine. They have made my life in Bordeaux enjoyable and memorable. Michael, thank you for being my rock i Climbing partner and teaching me many professional life skills. I am also very grateful to all my friends at the lab who have supported me during the last few years: Anass, Chiheb, Biswajit, Anis, Guillaume and Brendan. I also thank our previous Post-Doctoral fellow, Augusto, for his help. fixed beta and dynamic beta. β i = κ * β (i-1) , i: iteration value. We note that larger β values converge slower than small ones. Fixed β can provide comparable performance to dynamic if µ is large. . . . . . . . . LIST OF FIGURES 3.5 Visual demonstration of the structural similarity (SSIM) measurement system [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]). . . . . . . . . . . . . . . . . . . . . . . . 3.8 PSNR values after deconvolution using the regularizer ϕ p q (x) as a function of p and q. The best ϕ p q (x) function for Lena lies around the region where p = 1.1 and q = 2.9. The best function for Cameraman lies around the region where p = 0.9 and q = 3 . . . . . . . . . . . . . . . 3.9 (a) Visual illustration of 19×19 kernel used inexperiment of Figure (3.10) of [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF]) The kernel used. . . . . . . . . . . . . 3.10 (a) Visual illustration of deconvolution for Lena and a piecewise constant image. (top-left) Original, (top-right) blurred with 19 × 19 kernel, (bottom-left) 3/4 and the proposed regularizer (bottom-right) ϕ p q (x) . For both images, an improvement can be assessed visually. The PSNR values for all images have been provided in Table 3.3 Top: Display of the singularity exponents of the raw, unfiltered map.

List of Figures

Below: after noise reduction, p = 1.5, q = -1, λ = 0.1. . . . . . . . . . 97 5.3 Singularity exponents of the Musca observation Herschel at 250 µm, zoom on subregion, after noise reduction, p = 1.5, q = -1, λ = 0.1.

Compare with (Yahia et al. (2021)), fig. 7. . . . . . . . . . . . . . . . . 98 iii LIST OF FIGURES 5.4 Visualization of the log of gradient's norms, Musca observation Herschel map, after noise reduction, p = 1.5, q = -1, λ = 0.1. Although the log of gradient's norms enhances the striations, the resulting spectrum, shown in Figure 5.5 is very different than the one obtained with the correlation measure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.5 Left: comparison of obtained singularity spectra for Musca observation map: raw data, log of gradient's norms, correlation measure. The red and green curves are generated from the filtered observation map, the blue spectrum is generated from the raw map.The green spectrum has been computed from the values of h(x) given by Equation (5.6). The two other spectra (blue and red) make use of the values h(x) from formula (B.7) in [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]. Right: comparison between the singularity spectrum of Musca unfiltered (blue) computed with h(x)

given by Equation (5.6) and the singularity spectrum of Musca with noise reduction p = 1.5, q = -1, λ = 0.1 computed with h(x) given by Equation (5.6)

3.3

Comparison of PSNR and SSIM of reconstructed images using regularizers: l 1 ,l p and ϕ p q (x). The optimal value of µ is also reported. There exists one ϕ p q (x)) able to provide better performance in both scales.The values in bold correspond to the best performance. . . . . . . . . . . .

Chapter 1: General Introduction

This chapter gives a brief introduction and information on the super-resolution technique. In this chapter, we present our super-resolution framework as the basis of our main work and contribution.

Motivation

Images and image processing play an important role in this modern era. Cameras can be found everywhere; it has become the most common tool to create data in the form of an image. With the development of smartphones and tablets, it has become both fun and practical to capture images. I2S is a world-leading company in the field of image processing. The products consist of book scanners, microscopes, 3D image construction, image processing algorithms and cameras for various applications. Our main goals and objectives are to design and implement a super-resolution technique that will allow the scanners of I2S to capture images with higher resolution and image quality. This work is a CIFRE thesis, and the main goal is to provide software solutions by scientific research considering industry demands and needs. The CIFRE (Conventions Industrielles de Formation par la Recherche) system allows French companies, local authorities or associations to entrust a doctoral candidate with an assignment in the framework of a research collaboration with an academic research laboratory affiliated to a doctoral school .

Super-resolution is a new area that has been considered to be implemented within the new camera developed by I2S with the intention of integrating it into various book scanners. This technique will allow I2S to produce images with up to 2 times higher resolution with higher quality and accuracy. In the beginning, the idea and structure of super-resolution will be presented.

Super-resolution

The digital images acquired have a restricted spatial resolution due to the limitations of solid-state sensors such as CCD or CMOS in the camera, which reduces their practical value. Furthermore, since the information each pixel records are easily polluted during the acquisition and transmission process, digital images can be degraded by blurring or noise due to lens limitations, electronics and photonic noise. Super-resolution is a technique that considers both resolution limitations and normal image degradation. The super-resolution technique can create a higher-resolution image or sequence by using redundant information between multi-frame images with relative sub-pixel motion. It has been widely used with medical images (Liu et al. (2013)), remote sensing images (Zhang et al. (2012a,b)) and video surveillance [START_REF] Cheng | Classification-based video superresolution using artificial neural networks[END_REF]; [START_REF] Ng | A total variation regularization based superresolution reconstruction algorithm for digital video[END_REF]).

Imaging techniques have been rapidly developed and expanded in the last decades, and image resolution has reached a new level. The question is, therefore: are image resolution enhancement techniques still meaningful?

Although the high-definition displays in recent years have reached a new level of resolution (e.g., 1920×10804 pixels for HDTV, 3840×2160 pixels for some ultra-HDTV, and 2048 × 1536 pixels for some mobile devices), the need for resolution enhancement is still a key component in many applications [START_REF] Takeda | Super-resolution without explicit subpixel motion estimation[END_REF]).

For instance, to guarantee the long-term stable operation of the image capturing devices and the appropriate frame rate for dynamic scenes, digital image systems tend to sacrifice resolution to some extent. A similar situation exists in the remote sensing field: there is always a trade-off between a captured image's spatial, spectral, and temporal resolutions; an image can not have all the resolutions at the highest level at once. In the case of medical imaging, within each imaging modality, specific physical laws are in control, defining the meaning of noise and the sensitivity of the imaging process. How to extract 3D models of the human structure with high-resolution images while reducing the level of radiation remains a challenge (Greenspan (2009a); [START_REF] Kennedy | Super-resolution in pet imaging[END_REF]).

Based on these facts, the current techniques cannot satisfy the demands. Resolution enhancement is therefore still necessary, especially in fields such as video surveillance, medical diagnosis, and remote sensing applications. Considering the high cost and the limitations of resolution enhancement through "hardware" techniques, especially for large-scale imaging devices, signal processing methods which are known as superresolution (SR), have become a potential way to obtain high-resolution (HR) images.

With SR methods, we can go beyond the limit of the low-resolution (LR) observations, rather than improving the hardware devices with difficulties such as limited space of sensor for an embedded solution and obstacles of reducing the cell size for a sensor.

SR is a technique that reconstructs a higher-resolution image or sequence from the observed LR images. Technically, SR can be considered as multi-frame or single-frame based on the input LR images [START_REF] Borman | Spatial resolution enhancement of low-resolution image sequencesa comprehensive review with directions for future research[END_REF]; [START_REF] Elad | Super-resolution reconstruction of image sequences[END_REF]; Farsiu et al. (2004a); [START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF]). Suppose multiple images of the same scene with sub-pixel misalignment and misplacement can be obtained. In that case, we can reconstruct a higher-resolution image or image sequence by using the complementary information between them, as Figure (1.1) shows. By developing SR using software solutions, we have the advantage of not needing to produce a new, more prominent and more expensive sensor. Furthermore, we can use an existing sensor (generally low cost).

Another constrain is that we can not reduce the size of sensor cells indefinitely (due to physical limitation). SR with software solution needs only a few images from existing acquisition hardware.

A Brief Review of Super-resolution Methods

In this section, we tend to review the strategies and current issues for SR with multiple observations. The key problem is how to use the supplementary data among the non-inheritable low-resolution pictures. [START_REF] Harris | Diffraction and resolving power[END_REF][START_REF] Harris | Diffraction and resolving power[END_REF]) established the theoretical foundation for the SR problem by introducing the theorems to solve the diffraction problem within an optical system. Twenty years later, Tsai and Huang (Tsai (1984)) first used the idea of SR to improve the spatial resolution of Landsat TM images.

Since then, several researchers have started to concentrate on SR, either in theoretical analysis or practical applications [START_REF] Altunbasak | Super-resolution still and video reconstruction from mpeg-coded video[END_REF]; [START_REF] Elad | A fast super-resolution reconstruction algorithm for pure translational motion and common space-invariant blur[END_REF]; [START_REF] Farsiu | Robust shift and add approach to superresolution[END_REF][START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF]; Greenspan (2009a); [START_REF] Keren | Image sequence enhancement using sub-pixel displacements[END_REF]; [START_REF] Liu | On bayesian adaptive video super resolution[END_REF]; [START_REF] Zhang | Super-resolution reconstruction for multi-angle remote sensing images considering resolution differences[END_REF][START_REF] Zhang | A blind super-resolution reconstruction method considering image registration errors[END_REF]Zhang et al. ( , 2010a[START_REF] Zhang | Application of tikhonov regularization to super-resolution reconstruction of brain mri images[END_REF]; [START_REF] Zhuang | Hallucinating faces: Lph super-resolution and neighbor reconstruction for residue compensation[END_REF]).

The number of input and output photos involved in the SR pipeline can characterize SR methods. Single-image single-output (SISO) super-resolution occurs when a single high-resolution (HR) image is created from a single low-resolution (LR) image. The SISO is not really regarded as a super-resolution technique. It is more like an upsampling with a blur prevention algorithm (which is better than simple interpolation).

The capacity to achieve resolution and quality increase is a potential use of SISO superresolution. Other SR methods use numerous LR frames to estimate a distinct HR image:

Multiple-input single-output (MISO) super-resolution is the term for this situation. An example application area is license plate recognition from a video stream to increase the alphanumeric recognition rates. Recent research has been conducted at reconstructing a set of HR images from an equivalent set of LR frames. This approach takes the name of

Super-resolution class Number of inputs Number of outputs

Applications

SISO

One One

Object recognition Image quality restoration

MISO

Multi One

Medical imaging Industrial imaging Astronomy imaging

MIMO

Multi Multi

Video enhancement Video surveillance As discussed, SR techniques can be classified by the number of the LR images involved, the actual reconstruction method, and the employed domain (spatial or frequency). Many survey papers on SR algorithms have mostly considered these parameters in their classification. In a primary taxonomy, SR algorithms can be classified based on their domain, i.e., the spatial domain or the frequency domain. Despite having the very first SR algorithms emerging from signal processing techniques in the frequency domain, it can be seen in many survey papers that the majority of these algorithms have been developed in the spatial domain.

In the beginning, most of the methods concentrated on the frequency domain [START_REF] Aizawa | Acquisition of very high resolution images using stereo cameras[END_REF]; [START_REF] Kim | Recursive reconstruction of high resolution image from noisy undersampled multiframes[END_REF]; [START_REF] Nguyen | A wavelet-based interpolation-restoration method for superresolution (wavelet superresolution)[END_REF]; [START_REF] Rhee | Discrete cosine transform based regularized high-resolution image reconstruction algorithm[END_REF]; [START_REF] Tsai | Multiframe image restoration and registration[END_REF]). Frequency domain algorithms can use the relationship between the HR image and the LR observations based on a simple theoretical basis with high computational efficiency. However, the methods have notable limitations, such as being sensitive to model errors and facing difficulty to handle more complicated motion models. As a result, this has prevented them from further development and progress.

Due to the drawbacks of the frequency domain algorithms, spatial domain methods have become the main trend [START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF].) The famous spatial domain methods can be categorized as: non-uniform interpolation [START_REF] Nguyen | A computationally efficient superresolution image reconstruction algorithm[END_REF]), iterative back-projection (IBP) [START_REF] Irani | Improving resolution by image registration[END_REF]), projection onto convex sets (POCS) [START_REF] Patti | High-resolution image reconstruction from a low-resolution image sequence in the presence of time-varying motion blur[END_REF]Stark and Oskoui (1989); [START_REF] Tekalp | High-resolution image reconstruction from lower-resolution image sequences and space-varying image restoration[END_REF]), the regularized methods [START_REF] Farsiu | Robust shift and add approach to superresolution[END_REF][START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF]; [START_REF] Liu | On bayesian adaptive video super resolution[END_REF]; [START_REF] Ng | A total variation regularization based superresolution reconstruction algorithm for digital video[END_REF]; [START_REF] Ng | A fast map algorithm for high-resolution image reconstruction with multisensors[END_REF]; [START_REF] Schultz | A bayesian approach to image expansion for improved definition[END_REF]; [START_REF] Ng | A total variation regularization based superresolution reconstruction algorithm for digital video[END_REF]; [START_REF] Takeda | Kernel regression for image processing and reconstruction[END_REF]), and a number of hybrid algorithms [START_REF] Elad | Restoration of a single superresolution image from several blurred, noisy, and undersampled measured images[END_REF]). Early review papers have provided specific descriptions, and explanations of those methods [START_REF] Nasrollahi | Super-resolution: a comprehensive survey[END_REF][START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF]; [START_REF] Tian | A survey on super-resolution imaging[END_REF]). Among them, the regularized methods are the most popular due to their effectiveness and flexibility. Therefore, most of the recent representative articles about SR have focused on regularized frameworks [START_REF] Babacan | Variational bayesian super resolution[END_REF]; [START_REF] Liu | On bayesian adaptive video super resolution[END_REF]; [START_REF] Ng | A total variation regularization based superresolution reconstruction algorithm for digital video[END_REF]Su et al. (2011);[START_REF] Takeda | Kernel regression for image processing and reconstruction[END_REF][START_REF] Takeda | Super-resolution without explicit subpixel motion estimation[END_REF]; [START_REF] Yuan | Adaptive multiple-frame image super-resolution based on u-curve[END_REF]; [START_REF] Zhang | A blind super-resolution reconstruction method considering image registration errors[END_REF]).

Observation Model

The imaging model, which corresponds to the observation model, is essential to SR, especially when using a regularized scheme. Before reconstruction, we need to observe the process by which the observed images have been gained. The image acquisition process is confronted with a set of degrading factors, such as optical diffraction, undersampling, relative motion, and system noise. Consequently, the following imaging model can be introduced:

y k = O k D k B k M k x + n k (1.1)
where there are k LR images participating in the reconstruction. As N 1k × N 2k is defined as the size of the kth input LR image, L 1k N 1k × L 2k N 2k is set as the size of the reconstructed HR data, which is determined by the horizontal and vertical magnification factors L 1k and L 2k . In Equation (1.1), x is the vector form of the reconstructed image with a constant size of L 1k N 1k L 2k N 2k × 1, y k is the vector form of the kth input dataset.

D k is the down-sampling matrix of size N 1k N 2k × L 1k N 1k L 2k N 2k . B k is the blur matrix of size L 1k N 1k L 2k N 2k × L 1k N 1k L 2k N 2k .
M k is the warp matrix describing the motion information (e.g. translation, rotation, etc.), n k (N 1k N 2k × 1) represents additive noise.

O k is the operator excluding the unobservable pixels from the kth image (Shen et al. (2015a)). We can obtain the observation model for single-frame SR when k = 1 in Equation (1.1). If D k and M k are excluded, it is a model for image restoration, only dealing with the problems of noise, blurring, or missing information:

y k = O k B k x + n k (1.2)

Frequency Domain Approach

A group of multi-frame SR methods use a frequency domain formulation for the SR problem. The main idea is that indications about high frequencies are spread across the numerous LR images in the form of aliased spectral frequencies. The first frequencydomain SR method can be traced back to Tsai and Huang (Tsai (1984)), who studied SR reconstruction from noise-free LR images. Their method consisted of first transforming the LR image data into the Discrete Fourier Transform (DFT) domain and then combining them following the relationship between the aliased DFT coefficients of the observed LR images. The approach is based on the following factors:

1. The shifting property of the Fourier transform.

2. The aliasing relationship between the continuous Fourier transform (CFT) and the DFT of observed LR images.

3.

The assumption that an original HR image is band-limited.

These factors make it possible to formulate an equation to express the relationship between the aliased DFT coefficients of the observed LR images to a sample of the CFT of an unknown HR image.

To reduce memory requirements and to bring down the computational costs, Rhee and Kang (Rhee and Kang (1999a)) exploited the Discrete Cosine Transform (DCT) instead of DFT. [START_REF] Woods | Stochastic methods for joint registration, restoration, and interpolation of multiple undersampled images[END_REF]), on the other hand, offered an iterative expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]) for simultaneously performing the registration, blind deconvolution, and interpolation operations.

The frequency-domain-based SR approaches have two main advantages. The first advantage is their theoretical simplicity: the relationship between the LR input images and the HR image can be clearly demonstrated. Consequently, frequency-domain-based SR approaches represent an intuitive way to enhance the details of the images by generalizing the high-frequency information presented in the LR images. Secondly, these approaches have low computational complexity and enable parallel implementations.

However, frequency-domain-based SR methods are insufficient to handle realworld practical applications since these methods assume that there only exists a global displacement between the observed images and the linear space invariant blur during the image acquisition process (i.e. B k is supposed to be the same for all the LR images).

Lately, researchers have begun to investigate the use of the wavelet transform to target the SR problem to recover the detailed information (generally the high-frequency information) lost or degraded during the image capturing process. The rationale behind this is the fact that the wavelet transform provides a powerful and efficient multi-scale representation of the image to recover the high-frequency information (Nguyen and Milanfar (2000a)). These approaches usually consider the observed LR images the lowpass filtered sub-bands of the unknown wavelet-transformed HR image. The aim is to estimate the finer scale sub-band coefficients, then apply the inverse wavelet transform to produce the HR image.

In Ji and Fermüller [START_REF] Ji | Wavelet-based super-resolution reconstruction: theory and algorithm[END_REF]; [START_REF] Ji | Robust wavelet-based super-resolution reconstruction: theory and algorithm[END_REF]) proposed a robust wavelet SR approach to handle the error encountered in both the registration computation and the blur identification computation process. Chappelli and Bose [START_REF] Chappalli | Simultaneous noise filtering and super-resolution with secondgeneration wavelets[END_REF]) incorporated a denoising stage into the conventional wavelet-domain SR framework to develop simultaneous denoising and SR reconstruction methodology.

Regularization-based approach

Based on the observation model explained above in Equation (1.1), the target is to reconstruct the HR image from a set of warped, blurred, noisy, and under-sampled measured images. As the model in (1.1) is ill-conditioned, SR turns out to be an ill-posed inverse problem due to the insufficient number of LR images and ill-conditioned blur operators. For undefined image x given its ill conditioned observed version y, Bayesian estimation uses the posterior conditional probability, p (x|y), to infer x. Maximum a posteriori probability (MAP) is the famous estimator for finding x. Using Bayes's rule, and (MAP) theory, the problem we need to solve can be transformed to the minimization problem as in [START_REF] Bouman | A generalized gaussian image model for edge-preserving map estimation[END_REF]; [START_REF] Ng | A total variation regularization based superresolution reconstruction algorithm for digital video[END_REF]):

x = argmax x p (x|y) = argmax x p (y|x) p(x) p(y) = argmax x p (y|x) p(x) = argmin x -log {p (y|x)} -log {p (x)} = argmin x k k=1 ρ(y k -H k x) + λU (x) = F (x) + λU (x) (1.3)
In the above formulation we used the fact that p(y) is not a function of x, and therefore it can be omitted. The term ρ(y, x) = -log {p (y|x)} is a modeling choice to expresses the probabilistic relationship between the desired image x and the measurements y. U (x) = -log {p (x)} brings in the influence of the statistical nature of the unknown into the formula. In (1.3), the first term, F (x) , is the data fidelity term, the second term; λU (x), is referred to as the regularization term, with U (x) being the energy function. The two terms are chosen to model the probabilities as energy functional. λ is the regularization parameter balancing these two terms. More information on this type of formulation will be provided in chapter 2. This approach of SR is in general known as

F (x) = k k=1 y k -H k x p p (1.4)
On the other hand, there has been some attention in choosing the 1 norm as the function ρ in (1.3) for image SR and restoration, where p = 1 in (1.3). The 2 norm corresponds to Gaussian distributed errors, the 1 norm mainly corresponds to the Laplacian error model. According to [START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF]), p = 1 results in a pixel-wise median, and p = 2 leads to a pixel-wise mean of all the measurements after motion compensation in the SR model. It has been proven that the 1 -norm fidelity can be more effective than the 2 -norm fidelity when the images contain non-Gaussian errors [START_REF] Cho | Handling outliers in non-blind image deconvolution[END_REF][START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF]).

For complicated types of noise model error, both the 1 norm and the 2013)). In cases with mixed error modes, the p norm function (1 p 2) is sometimes employed as the constraint function because of its convex property, and its applicability for the imaging model errors (Shen et al. (2015d).) When 1 p 2, it results in a weighted mean of the measurements. If the value of p is close to one, the solution is calculated with a more significant weight around the measures near the median value. In the case that the value of p is near two, the solution is approximated to the average value [START_REF] Farsiu | Fast and robust multiframe super resolution[END_REF]).

The regularization term plays a significant role in the regularized variational framework. As SR is a classical ill-posed inverse problem, regularization is therefore utilized to stabilize the inversion process [START_REF] Ng | A total variation regularization based superresolution reconstruction algorithm for digital video[END_REF]; [START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF]Sun et al. (2008)).

Regarding the Bayesian theorem, the regularization term represents the prior image modelling, providing the prior knowledge about the desired image [START_REF] Babacan | Variational bayesian super resolution[END_REF]; [START_REF] Beal | Variational algorithms for approximate Bayesian inference[END_REF]). Over the past ten years of dynamic development, there have been a large number of studies of various regularization for image restoration, and SR [START_REF] Bertaccini | An adaptive norm algorithm for image restoration[END_REF]; [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF]; [START_REF] Shen | Adaptive norm selection for regularized image restoration and super-resolution[END_REF]; [START_REF] Yuan | Multiframe super-resolution employing a spatially weighted total variation model[END_REF][START_REF] Yuan | Regional spatially adaptive total variation super-resolution with spatial information filtering and clustering[END_REF]).

In the early years, the smoothness of natural images was mainly studied, which results in the quadratic property of the regularizations [START_REF] Tikhonov | Solutions of ill-posed problems[END_REF]).

Tikhonov-based regularization is the representative smoothing constraint, whose energy function is usually defined as:

U (x) = Γx 2 2 (1.5)
where Γ is commonly an identity matrix or high-pass operator (e.g. a difference operator or a weighted Fourier operator), laplacian regularization is one of the famous regularizations used in SR and was developed from Tikhonov regularization by choosing the smoothing operator as the discrete 2-D operator.

The smoothing prior models are in contrast with the nature of images, where sharp details in images are always desirable for human beings in many applications, including remote sensing imaging, medical diagnosis and object recognition [START_REF] Greenspan | Super-resolution in medical imaging[END_REF]; [START_REF] Shen | Super-resolution reconstruction algorithm to modis remote sensing images[END_REF]; [START_REF] Zhang | A super-resolution reconstruction algorithm for surveillance images[END_REF]). Thus, 1 -norm based regularization are often a preference for their edge-preserving properties [START_REF] Goldstein | The split bregman method for l1-regularized problems[END_REF]; [START_REF] Pan | Efficient huber-markov edge-preserving image restoration[END_REF]). The Total Variation (TV) regularization was first proposed by Osher et al. [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]), based on the notion that an image is naturally "blocky" and discontinuous. The standard TV function is given as:

U (x) = ∂ ∂x 1 x + ∂ ∂x 2 x (1.6)
where ∂ ∂x 1 and ∂ ∂x 2 are the first-order image gradients in the horizontal and vertical directions respectively. Unlike the quadratic regularizations, edge information can be better preserved using TV regularization, with the 1 -norm to deal with the image information rather than the 2 -norm [START_REF] Bioucas-Dias | Total variation-based image deconvo-126 BIBLIOGRAPHY lution: a majorization-minimization approach[END_REF]; [START_REF] Zhu | An efficient primal-dual hybrid gradient algorithm for total variation image restoration[END_REF]).

Therefore, the TV prior model has been the most popular model for image processing in the last two decades and has been applied in fields such as image denoising, deblurring, segmentation, and SR [START_REF] Chan | Total variation blind deconvolution[END_REF]; [START_REF] Huang | A new total variation method for multiplicative noise removal[END_REF]).

I2S Super-resolution Framework

The These scene motions can occur due to the controlled motions in imaging systems, e.g., images acquired from orbiting satellites.

To implement SR in a real application at I2S, we have acquired the images for SR through hardware control by designing an imaging mechanism with hardware techniques.

In this method, the sensors can acquire images with known sub-pixel displacements. We capture four LR images with the translational shifts being (0, 0), (0, 0.5), (0.5, 0), and (0.5, 0.5). Our low-resolution images each have a half-pixel displacement of the scene. By doing so, we have achieved the goal of having low-resolution images, with each having new information which can be used to construct a HR image. The visual representation of our sub-pixel shift is shown in Figure (1.2). The resulting image will produce an image with up to 2 times higher resolution.

As a result of this resolution enhancement, we come across a negative effect on the final image. The final image visually does not bear the best quality and has the downside of a blurry effect on the image; this effect is the result of the fusion of the pixels from the four low-resolution images into one super image. This effect is displayed in simple formation in Figure (1.4). In this example, we only demonstrate two LR images where the second image is the 0.5-pixel displacement of the original image. In (1.4), picture (a)

is the original image and picture (b) is the image captured from the scene with 0.5-pixel displacement. The result shown in picture (c) shows that since we consider the 0.5-pixel movement of the sensor between two images, there are shared data and signals in which, after fusion, we get overlapped regions in the pixels. This fusion effect represents itself as a blurry effect.

A supplementary deconvolution algorithm is needed to improve the quality and accuracy to overcome this phenomenon. It is essential to attach a deconvolution algorithm after the super-resolved image to help it benefit from clear and visible quality.

The whole structure of our super-resolution scheme is depicted in Figure (1.5). We use four low-resolution images, where each image has 0.5-pixel displacement from the other. Next, we interpolate the pixel together to form a super-resolved image with up to 2 times higher resolution than the original image. Finally, we add a custom-designed deconvolution algorithm to overcome the blurry effect of the pixel interpolation. It is worth noting that the deconvolution algorithm has to be fast and accurate since the main application of this algorithm is for cameras mounted on the scanners of I2S to be used in practical application.

The following chapter will briefly review some of the famous optimization problems used for inverse problems, such as image denoising and image deconvolution. Later on, we will introduce the custom-designed and implemented deconvolution for our superresolution technique.

Structure of the thesis

In this chapter, we began by introducing the concept of super-resolution and its application in various areas. The necessity of having a super-resolution technique was discussed. Following on, a brief history and a review of the super-resolution was presented. It was shown that, in general, super-resolution is modelled as an inverse problem. Commonly, this inverse problem is tackled either in the frequency domain or spatial domain. Both methods have their advantage and disadvantage, and recently regularization method in the time domain has gained popularity. The super-resolution framework of I2S was demonstrated as the last topic. Despite being very effective in achieving a high-resolution image, this method lacked a significant setback in quality.

To solve this problem, we will introduce a deconvolution/denoising algorithm as a supplement in the next chapter.

In chapter 2, we will introduce inverse problems in signal and image processing.

We will describe the concept of an ill-posed inverse problem and the conditions for an inverse problem to become ill. Next, we will present image deconvolution as an inverse problem, and we will give a brief description of the point spread function or the kernel used in our algorithms. We will then provide background on sparsity and review the literature. Next, we will move on to the optimisation of inverse problems. We will present the definition of the proximal operator, where it is used to solve convex optimisation problems, and it will be the basis of our proposed algorithm. We will review three algorithms that are used commonly to solve optimisation problems. We will briefly describe the Douglas-Rashford splitting Algorithm, forward-backwards splitting, and ADMM algorithm.

Chapter 3 will provide our proposed image deconvolution algorithm. First, we will begin by giving an introduction of our motivation and why we prefer to use classical methods rather than new learning methods. Next, we will provide our problem formulation. We will review the variable splitting algorithm used in our proposed method and then demonstrate our proposed approach and general framework. Moving on, we will give an in-depth analysis of the parameters used in our image deconvolution algorithm and discuss their effect on the performance of the algorithm. At the end of the chapter, we will present our scientific experiments. We will evaluate the performance of image deconvolution and various images and measure the results with different image quality metrics.

Chapter 4 will provide information about the practical results of implementing the image deconvolution and super-resolution scheme in the scanners of I2S. We will shortly review the objectives of book scanning and the devices used for this purpose. Next, we will introduce MTF or FTM measurement, a standard measurement for image quality analysis in the industry. Moving on, we will review the implementation process of our algorithm from MATLAB to C++ programming language. The fifth section of chapter 4 will provide practical experiments and results on samples gathered from EAGLE and E -XTRA cameras of I2S. We will demonstrate the results in visual representation and also the image quality ratings using the MTF standard. Since our main objective is to develop a flexible algorithm, we will then examine the results of implementing the deconvolution algorithm on microscopic and Terahertz cameras of I2S.

In chapter 5, we will demonstrate the results obtained by applying our image deoncolution algoirhtm on Herschel observation maps. It will be shown that our algorithm provides better quality images with recently studied methods.

In chapter 6,we will demonstrate the second phase of our image deconvolution algorithm based on learning methods. First, we will review recent research and deep learning algorithms in various image processing applications such as medical diagnosis, detection and inverse problems. We will provide the concept of plug-and-play priors used in our image deconvolution algorithm. We will present our motivation and intuition of using such priors in our scheme. Our proposed plug-and-play approach uses Gaussian mixtures since they are good representatives of image priors. In the last section of this chapter, we will examine the experiments on various images captured with the scanner of I2S.

We will conclude the thesis in chapter 7, where we will summarize the main contributions and aspects. We will also talk briefly about the ongoing work extending our thesis.

Chapter 2: Inverse problems in signal and image processing

Before we give details of our frameworks and main contributions, we will discuss and review inverse problems in image and signal processing. Some background of optimization algorithms for inverse problems will be provided. We will introduce the proximal operator theory as the primary solution for our proposed method.

Introduction

This chapter will explain the mathematical formulation used in solving the problem of image deconvolution. It was shown in the previous chapter that our super-resolution problem could be modelled as an inverse problem. This chapter aims to formulate a deconvolution problem to tackle the blurring effect, which is encountered as an inverse problem.

In a simple form, inverse problems are the opposite of direct problems. In a direct problem, we intend to find an effect from the cause of phenomena, whereas in the inverse problem, we are faced with the effect, and we pursue to find the cause. The most popular and common practical case is when interpreting physical measurements of an unknown object of interest. In image processing, the direct problem is finding out how a given sharp photograph looks if the camera were incorrectly focused. The inverse problem of this example is well known as deblurring which is to find the sharp and clear image from a given blurry image.

In this example, the cause is the clear image, and the effect is the blurred image.

It is famously known that the Hubble telescope suffered a blurring problem after just starting to work where a deblurring had to be used as a correction algorithm.

In this thesis, we focus on the ill-posed inverse problems. In such circumstances, the inverse problem is much more difficult than solving the direct problem. To understand the ill-posed problem, we need to get familiar with the notion of a well-posed problem and its conditions, which was originally introduced by Jacques Hadamard as follows:

C1 : Existence. There has to be at least one solution.

C2: Uniqueness. There has to be at the most one solution.

C3: Stability. The solution depends continuously on the data.

We define H as the forward map, which is conceptually defined as H(cause)= effect.

The direct problem must be well-posed, which means H should be a well-defined, single-valued and continuous function. The inverse problem is ill-posed if H -1 does not exist or it is not continuous. Therefore, one of the three conditions fails for H -1 . Illposedness is related to interpretation tasks that are extremely sensitive to measurement and modelling errors. As a result, a successful inversion algorithm must be devised to be robust against measurement noise, computationally effective and mathematically justified by the appropriate analysis.

Image deconvolution as an inverse problem

Throughout this thesis, we denote vectors by lowercase letters, such as x and denote matrices by uppercase letters, such as X. We assume that vectors and matrices consist of real-valued elements. For instance, x ∈ R m represents a vector with m real-valued elements, and X ∈ R m×n represents a matrix with m rows and n columns that consists of real-valued elements. As introduced in chapter 1, in the case of discrete image deblurring, we consider the formulation below:

y = Hx + n (2.1)
In most cases, H is not very well known or is just an assumption. Also, in the case of deblurring, H does not depend on the data, and it is dependent on the acquisition system. Based on this fact, image deblurring (also known as image deconvolution) is considered an ill-posed inverse problem.

It is very tempting to try to solve the inverse problem of (2.1) by naive reconstruction:

x ≈ H -1 (y -n) (2.2)
However, since we are dealing with an ill-posed inverse problem, the approach (2.2) will fail. As described in section 1.3.3 regularization is needed for a successful and noise-robust solution of our linear inverse problem. We will explain the concept of regularization in the coming sections; however, before this, we will briefly explain the blur matrix H, which is also regarded as Point Spread Function (PSF).

Point spread function

The PSF kernel H models the intrinsic camera characteristics, such as diffraction when the light goes through a finite aperture, light averaging within the sensor itself, lens aberration, etc. The experimental procedure will minimize other blur sources like motion, atmospheric turbulence or defocus blur that may change from one snapshot to another. In this thesis, we do not seek to estimate these parameters. The diffraction kernel is determined by factors such as the shape, the size of the aperture, the focal length, and the wavelength of the monochromatic light.

Besides the diffraction kernel, other blur sources inherent in the camera's optical system are present in real-world cameras. These are mainly optical aberrations and anti-aliasing filters (that reduce aliasing but do not completely cancel these effects) introduced in the system prior to sampling [START_REF] Williams | Introduction to the optical transfer function[END_REF]; [START_REF] Zhao | Characteristic-analysis of optical low pass filter used in digital camera[END_REF]). The sampling process also introduces blur. With regards to the formulation of (2.1), image deconvolution is mainly divided into two categories: blind and non-blind deconvolution [START_REF] Campisi | Blind image deconvolution: theory and applications[END_REF]).

Blind image deconvolution seeks to estimate the proper image, assuming the blur is unknown. In this thesis, we consider our deconvolution to be non-blind. The point spread function may arise from a mathematical model of the blurring process or actual measurements of the point spread function, despite having the freedom of changing the PSF, i.e. H and being regarded as an input, for the time being; we estimate the PSF to be a 5×5 Gaussian normalized matrix as shown in Figure (2.2).

Dimensions of the computed or measured PSF array are often much smaller than N , size of the image, due to the fact that the spreading of the light takes place in a small region around the pixels. We can always conceptually augment such a small array with zeros to make the PSF array as large as needed in the equations above.

Background on sparsity

Parsimony is regarded as the representation of a phenomenon with the fewest elements possible. In signal processing, sparse representations intend to represent signals with as few as possible significant coefficients. This "parsimony" principle and the tendency of having "few elements", a.k.a. sparsity, is now a principal focus of attention in diverse signal/image processing, computer vision and machine learning applications. A matrix or a vector is said to be sparse if it has only a few (compared to the total size) non-zero entries.

In the case of a vector, its l 0 norm, x 0 , is simply the number of non-zero entries in the vector. From sparsity point of view, l 0 is the most appropriate norm, while l 2 says almost nothing. On the other hand, finding the sparsest solution is an optimization from the point of view of mathematical analysis.

Compressed sensing permits to sense signals directly with few samples going beyond the classical Nyquist rate. By using sparsity as a prior in an optimization framework, exact recovery can be achieved from a small set of linear non-adaptive measurements [START_REF] Baraniuk | Compressive sensing [lecture notes[END_REF]; [START_REF] Candès | Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information[END_REF]; [START_REF] Bibliography Donoho | Compressed sensing[END_REF]; [START_REF] Romberg | Imaging via compressive sampling[END_REF]).

Compressed sensing is mainly known for finding exact or approximate solutions for underdetermined linear systems of equations, which could not be solved using traditional linear algebra techniques. It showed that sampling under the Shannon-Nyquist rate is no longer impossible.

Sparsity is also used enormously in image restoration. Primary work of sparsitybased image restoration used thresholding wavelet coefficients. The motivation behind wavelet shrinkage is twofold. First, it is to be noted that the regularity of transient events can be observed across different scales -such an idea originated from Marr's vision theory and is at the basis of scale-space theory [START_REF] Lindeberg | Scale-space theory in computer vision[END_REF]). Wavelets have offered a principled method of decomposing signals in the scale space with specific desirable properties [START_REF] Champagnat | Deconvolution of sparse spike trains accounting for wavelet phase shifts and colored noise[END_REF][START_REF] Ruder | An overview of gradient descent optimization algorithms[END_REF]). Second, despite the difficulty of conducting signal-noise separation for a single sample or coefficient (microscopic level), the ensemble (macroscopic level) property of a signal is in most cases sufficient to distinguish itself from that of noise.

Later on modern sparsity-based image processing and restoration algorithms try to use various dictionary learning methods for compact representation.

As mentioned in the previous chapter, prior information in the formulation plays an important role in providing information about the reconstructed image; sparsity can be used extensively as prior information for regularized image restoration. A comprehensive study of sparse methods in image processing is done in [START_REF] Badri | Sparse and Scale-Invariant Methods in Image Processing[END_REF]).

Norms

Just in the phrases above, the term "norm" has been taken with somewhat excessive license. Let X be a real or complex vector space. A function . : Let us introduce the l p norms. Define, for x ∈ R n :

X → R + is called a norm if ∀x ∈ X:
x p = ( n i=1 |x i | p ) 1/p (2.3) if 1 ≤ p < ∞ this function is a norm. However, if 0 < p < 1 this function is a quasi-norm. Note that x 0 = lim p→0 x p . It is neither a norm or a quasi-norm.
To understand the relationship between sparsity and the l p norm, one can consider the case n = 2 and the various shapes of the l p unit ball in function of p (see Figure (2.3)):

if x 2 ≤ 1, the components of x can take all values uniformly between -1 and 1. if

x p ≤ 1 the components of x are constrained by the shape of the ball. The more p -→ 0, the stronger is the constraint. The less p, the more x is sparse.

Optimization of inverse problem

In this section, we will briefly explain some of the standard optimization methods for regularized image reconstruction. Before we continue to demonstrate the optimization methods, we have to note that in our work, we concentrate primarily on proximal algorithms since we desire to have fast and accurate results due to our practical goals.

As mentioned in previous sections of this thesis, many constrained convex optimization problems are in the following form:

Minimize H(x) = f (x) + g(x) (2.4)
where f is convex and differentiable, g is a convex function (differentiable or not),

x ∈ R n . We also mention some generalizations where g can be non-convex [START_REF] Parikh | Proximal algorithms[END_REF]). The formulation of the form (2.4) appears in many machine learning, dictionary learning and image processing applications.

Proximal operators are among the popular methods to solve the problem of the form (2.4). Various methods to solve problems of the form (2.4) have emerged in the literature based on proximal operators.

Definition 2.1

Definition: Let X = R n equipped with the standard norm 2 and f : X → R ∪ {+∞} a convex function. The proximal operator [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF]) (also known as the proximal mapping [START_REF] Beck | First-order methods in optimization[END_REF]) ([Ch. 6]) of f is the function ♣ prox f : X → X defined by:

prox f (y) = argmin x 1 2 x -y 2 X + f (x) (2.5)
Note that since f is supposed to be convex, the existence and the uniqueness of the minimum is guaranteed. Note that:

prox λf (y) = argmin x 1 2 x -y 2 X + λf (x) (2.6)
For a general reference on proximal operators see Figure (2.4). Using Equation (2.4), we start at the blue points, we aim to reach the minimum of f not globally. we keep on finding a minimum while staying inside a ball i.e. we stay close to y in relation to the constraint, blue point moves to the corresponding red points. This is continued until the blue points reach the minimum of f .

The main theorem about proximal operators comes from the minimum of the convex function f and the fixed points of its associated proximal operator [START_REF] Parikh | Proximal algorithms[END_REF]). The point x minimizes f if and only if :

x = prox f ( x) (2.7)
i.e., if x is a fixed point of prox f . The proximal operator can be computed analytically for many important f functions. The notion of a proximal operator is not restricted to convex functions. Indeed, it can been generalized for some proper (not necessarily convex) function f [START_REF] Urruty | Convex analysis and minimization algorithms[END_REF]). As a famous example in the popular case when f is the l 1 -norm, the solution reduces to the soft-thresholding [START_REF] Donoho | De-noising by soft-thresholding[END_REF])

that is given as follows: let

X = R and f (x) = x 1 . prox λf (y) = max(0, |y| -λ)sign(y) (2.8)
Many algorithms in the literature rely on proximal splitting, which entails introducing an intermediate variable and breaking the primary problem into subproblems corresponding to proximal operators. These algorithms are very popular in image processing because they deliver a high-quality result and scale well to large-scale issues while remaining efficient. It's also worth noting that the proximal operator is a pointwise operator, which means that the shrinkage is applied to each element separately, allowing proximal solvers to be implemented in parallel. In chapter 3, we will demonstrate the use of the proximal algorithm, which is the basis of our image deconvolution algorithm used in this thesis. The following part will briefly review some of the popular proximal-based algorithms for solving inverse problems.

Douglas-Rachford splitting algorithm

Considering the formulation of Equation ( 2.4) and with the assumption that f and g are closed convex functions, the Douglas-Rachford [START_REF] Combettes | A douglas-rachford splitting approach to nonsmooth convex variational signal recovery[END_REF])

iteration can be defined as follow:

Start at any y 0 and repeat for k = 0, 1, ...,:

x k+1 = prox f (y k ) y k+1 = y k + prox g (2x k+1 -y k ) -x k+1 (2.9)
This method is very handy when f and g have inexpensive proximal operators. The important point is that f , and g do not need to be smooth, and one only needs them to be "proximable". This algorithm was introduced in [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF]) as a generalization of an algorithm introduced by Douglas and Rachford in the case of quadratic minimization (which corresponds to the solution of a linear system).

Forward-backward splitting

In many situations and applications, the function g is neither differentiable nor even finite-valued. Faced with this challenge, problem Equation (2.4) cannot be minimized using normal and simple gradient-descent methods, which exist in many literature [START_REF] Ruder | An overview of gradient descent optimization algorithms[END_REF]). Forward-backward splitting uses the proximal operator, prox τ g (z), and finds a point close to the minimizer of g without straying too far from starting point z and is often referred to as a backward (or implicit) gradient-descent step with stepsize τ . If the proximal operator can be estimated, then it is feasible to solve (2.4) efficiently using the Forward-Backward Splitting (FBS) method (also known as the proximal gradient method). In simple terms, FBS can handle non-differentiable objectives and convex constraints while maintaining the simplicity of gradient-descent methods.

Due to the vast applications of FBS and its utility for sparse coding and regression, many variants of FBS have been developed to improve performance and ease of use. In its raw form, FBS requires the user to choose several convergence parameters that strongly affect the algorithm's performance and reliability. These include step sizes, stopping condition parameters, acceleration schemes, stability conditions, and initialization.

Forward-Backwards Splitting is a two-stage method that addresses each term in (2.6) separately. The FBS method is listed in Algorithm 1.

Algorithm 1: Forward-backward splitting while not converged do

1: xk+1 = x k -τ k ∇f (x k ) 2: x k+1 = prox g (x k+1 , τ k ) = argmin x τ k g(x) + 1 2 x -xk+1 2 3: end while ;
Algorithm 1 interchanges between forwarding gradient descent on f and backward gradient descent on g. Using a backward step for g can be advantageous in several ways. On the other hand, it is shown that step (2) always has a unique, well-defined solution [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]). The backward step has an important effect on the convergence of FBS.

The exact detail of the FBS algorithm can be found in numerous articles and books.

However, here we have just given a simple explanation of the algorithm in summary to provide a general idea.

Alternating Direction Method of Multipliers (ADMM)

ADMM is designed to solve the problem of the following type:

minimize f (x) + g(z) subject to Ax + Bz = c (2.10) with x ∈ R n and z ∈ R m , where A ∈ R p×n and c ∈ R p .
In most cases the assumption is that f and g are convex functions. The optimal value of the problem (2.10) will be denoted by:

p = inf {f (x) + g(z)|Ax + Bz = c} (2.11)
Note that the objective is separable into two sets of variables. To solve the problem, one considers the Augmented Lagrangian:

L ρ (x, z, y) = f (x) + g(z) + y T (Ax + Bz -c) + (ρ/2) Ax + Bx -c 2 2 (2.12)
where ρ is the step length. The original ADMM algorithm was proposed by Gabay and Mercier [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF]), and Glowinski, and Marrocco [START_REF] Glowinski | Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problèmes de dirichlet non linéaires[END_REF]). Many further findings in ADMM were discovered by Eckstein and Bertsekas [START_REF] Eckstein | On the douglas-rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF]). At iteration k, we minimize for x, then z, and finally, update y, keeping the other variables constant during each minimization.

This gives us the ADMM algorithm as below:

x k+1 = argmin x L p (x, z k , y k ) z k+1 = argmin x L p (x k+1 , z, y k ) y k+1 = y k + ρ(Ax k+1 + Bz k+1 -c) (2.13)
The constraint variables (A, B, and C) play an important role in the algorithm's effectiveness. A detail to consider is that the algorithm decouples the objective function on variables x and z and minimize them independently. It consists of an x-minimization step, z-minimization step, and a dual variable update. The dual variable update uses a step size parameter ρ. In ADMM, x and z are updated in an alternating or sequential fashion, accounting for alternating directions. The first two expressions in (2.13) can also be expressed in the proximal form.

Apart from the algorithm reviewed here, there also exists a popular algorithm socalled Alternating Minimization Algorithm [START_REF] Tao | Alternating direction algorithms for total variation deconvolution in image reconstruction[END_REF]), which is used for image reconstruction. We will mainly use this algorithm in our deconvolution algorithm.

This method is presented in depth in chapter 3, introducing our Tools and problem formulation.

Chapter conclusion

In this chapter, we began by introducing the fundamentals of modelling the image deconvolution problem as an inverse problem application and optimization of minimization problems. We explained the phenomenon of sparsity and its application on compressed sensing and image restoration. Later on in the chapter, we reviewed famous algorithms and optimization methods used for image reconstruction. The proximal operator was explained as an essential tool for optimization algorithms. We briefly examined three algorithms frequently used to solve inverse problems in various applications. In the next chapter, we will present our tools and problem formulation used for image deconvolution algorithms throughout our work. We will also introduce our work and approach. We will talk about our modelling and solution based on variable splitting and proximal operator to develop a general framework for a flexible and robust image deconvolution algorithm.

Chapter 3: Proposed method and experiments

Introduction

This chapter will explain the mathematical formulation and tools we have used for our image deconvolution algorithm. In our framework, we will use the proximal operator with variable splitting technique to reduce our algorithm's complexity and to have a fast algorithm.

In general, inverse problems have been tackled with either "classical methods of optimisation" or learning methods such as Convolutional Neural Networks(CNN). We would approach our deconvolution inverse problem using classical methods in our work, and this is due to some factors which make the classical methods favourable.

Despite the good performance of CNN or other learning algorithms, such algorithms are not always suitable in practice. These methods demand a large amount of training data to perform well, which may not always be available. Additionally, the performance of learning methods can be only guaranteed for a specific kind of data; changes in the nature of data may lead to loss of performance [START_REF] Mccann | Convolutional neural networks for inverse problems in imaging: A review[END_REF]). Learning methods such as CNN's work like a black box without tuning. Since the image quality will be measured with international standards at I2S, which will be explained in chapter 4, it is essential for I2S to have an algorithm tuned with parameters and does not act as a black box with no flexibility.

Since this thesis and its algorithms are for the practical purpose and implementation in I2S devices, it is essential to look for methods and algorithms that would best suit the practice. The primary usage of the algorithm will be for CopiBook scanners of I2S.

These scanners are designed for scanning books posters of various sizes for digitalisation.

Therefore, the algorithm must be as fast as possible since the scanner can not bear being slow to scan a book with many pages. Also, in our case, we do not have access to many is blurred by a kernel different from the ones in the training set. One of the primary key purposes of our deconvolution is to design an algorithm that can be flexible and can be used for various products, applications and circumstances at I2S. As shown in the coming chapters, thanks to our algorithm, I2S could use the algorithm not just for their scanners but also for their microscopes and terahertz camera. We have also applied our algorithm for astronomy images known as "observational maps", presented in a separate chapter.

Due to the explained factors, we will consider classical optimisation methods for our image deconvolution algorithm for the first stage of our implementation. One of the goals of our thesis is to develop a fast and low complex algorithm by preserving the accuracy for practical application in industry for products like connected image acquisition systems and embedded systems with limited computation capability. In the past decade, many papers have been published for non-blind image deconvolution either based on classical or learning methods without considering this practical limitation.

Most of these applications do not prove to be practical or feasible in computationally constrained conditions where complexity, accuracy and speed are essential. Our new algorithm, based on proximal operators, proves to be accurate and effective. It also considers a broad family of sparsity inducing functions which increases the flexibility of our algorithm.

Problem Formulation

As mentioned in chapter 2, we will consider Equation (3.2) as our formulation for image deconvolution modelling. In digital image processing, the general, discrete model for a linear degradation caused by blur and additive noise is formulated as:

y(k) = s∈Ω H H(k, s)x(s) + n(k) k = (k 1 , k 2 ) ∈ Ω (3.
y = Hx + n (3.2)
where H denotes a known convolution operator in matrix form and n additive Gaussian noise with zero mean and known variance σ 2 , for some of our experiments, we will consider H to be a 5×5 Gaussian matrix as shown on Figure (2.2). Typically, as discussed, these inverse problems are ill-posed; therefore, there is no exact solution, or the solution is not unique because H is not invertible or is extremely ill-conditioned. One of the effective methods of approaching this problem is using regularization methods (Liu et al. (2013)). In this case, we will consider the following optimization problem:

argmin x µ 2 ||Hx -y|| 2 2 + φ(Dx) (3.3)
In Equation (3.3), for simplicity we write φ(Dx) instead of:

φ(D x 1 x) + φ(D x 2 x),
which are respectively the matrix form of the first order gradients in both directions:

d x 1 = [1, -1] and d x 2 = [1, -1] T .
As discussed in depth in chapter 1, the first term of the equation measures the closeness of solution to the observation y, it is known as data fidelity term. The second term φ(u), the regularizer, enforces a prior knowledge about x into the solution. The parameter µ > 0 assigns a trade-off between the two assumptions.

Our thesis proposes a new fast image deconvolution algorithm using approximated proximal operators. It combines the idea of variable splitting with new sparsity inducing regularizers introduced in (Gholami and Hosseini ( 2011)). To the best of our knowledge, this has not been considered for image deconvolution applications. The motivation behind this is twofold: first, some functions have been shown to have better sparsity inducing performance relative to other regularizers such as p . Second, a quadratic penalty approach requires low computational complexity.

Tools and Proposed method

In this part, we will begin to explain the tools used for our optimization. First, we will explain variable splitting, a well-known method for optimization. Next, we will present the general framework that we have used for image deconvolution, and also, we will provide the proximal solution of the general framework.

Variable Splitting

In this section, we will give a demonstration of variable splitting, which has been utilized in several image processing applications. Here we consider an unconstrained optimization problem where the objective function consists of two functions:

x = argmin x {f (x) + g(x)} (3.4)
In the context of the thesis and our problem we assume the following:

f (x) = µ 2 ||Hx -y|| 2 2 g(x) = φ(Dx) (3.5)
In this point we will take the unconstrained optimization of (3.4) and we will split the variables into two variables x and u under the constraint that x = u. As a result the new problem is given by

(x, û) = argmin x=u {f (x) + g(u)} (3.6)
The question is why this method is favoured to the latter unconstrained problem?

The answer is that we can minimize two functions separately and use different techniques for each one of them. A penalty function is added to Equation (3.6) to attack and penalize the difference between x and u.

(x, û) = argmin x f (x) + g(v) + β 2 x -u 2 2 (3.7)
The constant parameter β controls the magnitude of the penalty function. It is important to note that the choice of β will affect the rate of convergence of the algorithm.

As the value of β increases, the constraint of x = u is more aggressively enforced.

However, as β becomes large, the problem becomes stiffer. The choice of β falls into the so-called Goldilocks Principle :

i) If β is too big then the algorithm becomes too stiff i.e ill conditioned.

ii) If β is too small, the algorithm will converge slowly.

The proper value of β is always a trade-off between (i) and (ii). In our work, as will be thoroughly explained later, β is constant β = 256. The variable splitting idea has been recently used in several image processing applications. The method was used in [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF]) to achieve a fast algorithm for TV-based image enhancement. It was also used in (Bioucas-Dias and Figueiredo ( 2008)) to ease the problem of involving compound (multiple) regularizers; i.e., instead of a single regularizer g(x) = φ(Dx) in (3.5), a linear combination of two (or more) regularizers g(x) = φ 1 (Dx)+φ 2 (Dx)+... are implemented.

Similar usage of variable splitting approach underlies the proposed split-Bregman methods [START_REF] Goldstein | The split bregman method for l1-regularized problems[END_REF]). On the other hand, instead of using a quadratic penalty technique, those methods attack the constrained optimization problem directly with a Bregman iterative algorithm [START_REF] Yin | Bregman iterative algorithms for \ell_1minimization with applications to compressed sensing[END_REF]). It has been proven that, when g is a linear function, i.e., g(u) = Gu, the Bregman iterative algorithm is equivalent to the augmented Lagrangian method [START_REF] Yin | Bregman iterative algorithms for \ell_1minimization with applications to compressed sensing[END_REF]).

Proposed Approach

As mentioned before we begin by considering the following model:

argmin x µ 2 ||Hx -y|| 2 2 + φ(Dx) (3.8)
By using variable splitting, we introduce a new variable, say u, to serve as the argument of φ(.) under the constraint that u = Dx. A quadratic function to penalize the difference, needs to be added. After applying the modification we consider the constrained problem:

argmin x,u µ 2 ||Hx -y|| 2 2 + φ(u) s.t u = Dx (3.9)
Given the constraint, it must be that u = Dx, so the solution to Equations (3.8) and (3.9) must be the same or very close. In order to enforce the constraint, we will add an additional term to the cost function that penalizes large differences between u and Dx. This then results in the unconstrained optimization problem given by argmin

x,u µ 2 ||Hx -y|| 2 2 + φ(u) + β 2 ||Dx -u|| 2 2 (3.10)
where β is a constant that controls the gain of the penalty term. It also enforces the difference between Dx and u to be as small as possible.

The problem formulated in the form of Equation (3.10) can be solved by an al-ternating minimization scheme, i.e. solving two sub-problems iteratively. It is easy to minimize the objective function in (3.10) with respect to either u or x.

(S1) argmin u φ(u) + β 2 ||Dx -u|| 2 2 (S2) argmin x µ 2 ||Hx -y|| 2 2 + β 2 ||Dx -u|| 2 2 (3.11)
The algorithm decouples the objective function on variables u and x and minimizes them independently. It consists of an u-minimization step and an x-minimization step (3.11). The solution of Equation ( S1) is calculated using shrinkage operators, which will be explained later on.

For a fixed u, (S2) is quadratic in x, and the minimizer x is calculated by taking the derivative of the Equation (S2) and setting it to zero to get:

(D T D + µ β H T H)x = D T u i + µ β H T y (3.12)
After splitting the discrete gradient operator we get:

(D T x 1 D x 1 + D T x 2 D x 2 + µ β H T H)x = D T x 1 u x 1 + D T x 2 u x 2 + µ β H T y (3.13)
Under the periodic boundary condition for u, D T x 1 D x 1 and D T x 2 D x 2 are all block circulant. In this case, the Hessian matrix on the left-hand side of (3.13) can be diagonalized by 2D discrete Fourier transform F. Using the convolution theorem of Fourier transforms, we can write:

x = F -1   F(D x 1 ) F(u x 1 ) + F(D x 2 ) F(u x 2 ) + ( µ β )F(H) F(y) F(D x 1 ) F(D x 1 ) + F(D x 2 ) F(D x 2 ) + ( µ β )F(H) F(H)   (3.14)
where F denotes the Fourier Transform, " ¯" denotes complex conjugate and " " component-wise multiplication, the division is computed component-wise. An advantage of this formulation is that Fast-Fourier will reduce computational complexity and simplify the algorithm.

In the following section, we will introduce the general regularization function which is used in this article and also we will show how we use the proximal operator as the solution of (S1).

General Framework

In this section, we will review some well-known sparsity-inducing regularizers and introduce a family of regularizers which have shown to provide better sparsity-

enhancing performance. The 1 norm, φ(u) = ||x|| 1 = i |u i |
, is probably the most popular regularizer. For this regularizer, (S1) has a closed-form solution: the soft threshold operator [START_REF] Donoho | De-noising by soft-thresholding[END_REF]). It has been shown that non-convex penalties can achieve better performance to promote sparsity [START_REF] Kochenderfer | Algorithms for optimization[END_REF]), e.g.

φ(u) = ||x|| p = i |u i | p (the p quasi-norm, 0 < p < 1)
. However, analytic solutions are available only for some values of p, i.e.1/2 and 2/3 [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF]).

Only numerical methods exist for other p values [START_REF] Nguyen | A computationally efficient superresolution image reconstruction algorithm[END_REF]).

As mentioned and discussed in chapter 2, closed-form solutions for (S1) using nonsmooth regularizers can be derived by the theoretical framework presented in [START_REF] Parikh | Proximal algorithms[END_REF]). Nevertheless, many regularizers do not have an analytic expression.

Due to the lack of closed-form solutions, some approximate methods have been studied for particular regularizers. As an example, a proximal first-order approximation for edge preservation has been studied in [START_REF] Badri | Fast edge-aware processing via first order proximal approximation[END_REF]).

In our thesis, we will take a broader approach. We will consider multiple sparsityinducing regularizers (general framework), including p . The proximal operators of these functions have been explained in detail in [START_REF] Gholami | A general framework for sparsity-based denoising and inversion[END_REF]). The general framework regularizers for p ∈ (0, 2] and q ∈ [-1, ∞) are formulated as:

ϕ p q (x) =      1 q (1 -(|x| p + 1) -q ) q = 0 log e (|x| p + 1) q = 0 (3.15)
Depending on the value of p and q, Equation ( 3.15) defines a potential function, each potential function has its own characteristic and can be used as prior information in second term of Equation (3.8). For instance, non-convex penalties can achieve better performance to promote sparsity [12], e.g.

φ(u) = ||x|| p = i |u i | p (the p quasi-norm, 0 < p < 1).
The potential function defined in Equation (3.15) for some different choices of p and q leads to some known potential functions as presented in Table 3.1.

p q Expression of ϕ p q (x) p -1

|x| p 1 1 |x| |x|+1 1 3 1 3 (1 -1 (|x|+1) 3 )
Table 3.1: Potential function and its corresponding functions according to the values of p and q.

Also, in Figure ( 3.1), we have depicted various forms of ϕ p q (x) presented in table 1. Each of the potential functions can be a prior information for any image denoising and deconvolution problem according to the criteria of the application being used. As mentioned before, it is the first time that a wide variety of image priors are being considered in an image deconvolution scheme. This allows us to use various functions according to practical circumstances, which is ideal for I2S scanners and datasets.

A&A proofs: manuscript no. debeam-deconvolution-denoise-arash Fig. 2: PSNR values after deconvolution using the regularizer ' p q pxq as a function of p and q. The region around p=1.35 and q=-1 provides the best performance. Closed-form solutions for (S1) using non-smooth regularizers (potential functions) can be derived by proximal operators defined in definition 2.1. The approximate proximal operators of (3.15) are provided as shrinkage/thresholding operators : The proposed general regularization and its approximation induce sparsity in the sense that several values, depending on the strength of regularization, will be exactly equal to zero.

prox 1 β ϕ p q (u) =        1 -p β |u| p-2 (|u| p +1) q+1 + u if |u| > η 1 β 0 otherwise
To assess the accuracy of the approximated proximal operators for the case of p , we carry out a comparison to see how close the proximal operators are to the exact numerical solution. Although the p norms are not differentiable, the precise solution of this particular function has been studied in detail in [START_REF] Lorenz | Non-convex variational denoising of images: Interpolation between hard and soft wavelet shrinkage[END_REF]). As it is demonstrated in Figure (3.2), The exact solution is represented in solid line with its corresponding approximation in dashed line. It is evident that the accuracy of proximal operators is high for all cases. It is worth mentioning that in some cases, the accuracy depends on the value of p: accuracy increases as p → 0. In general, this slight difference is so small that it can be neglected.

Proposed Algorithm

We now give the overall algorithm using the general framework for the sub-problem The algorithm can be controlled either by the number of iterations it needs or the cost function. It will be shown that thanks to our model of implementation, many iterations are not required for the algorithm to converge.

In many papers and articles, peak signal-to-noise ratio (PSNR) is used to measure the convergence of an algorithm. In this thesis, we would instead use Equation (3.3) itself to measure how effective the algorithm is in terms of deconvolution and also convergence.

f Herschel observation maps: deconvolution and sparse filtering for beam e↵ect and denoising.

t algorithm with various potential function configuration to obtain di↵erent results based on the corresponding regularizer. To illustrate the variation in performance, we in figure 2 show PSNR values for q P r´1, 3s and p P r0.1, 2s with each range divided into 64 equally spaced points. For this experiment, we have used one 1024 ˆ1024 sample with additive Gaussian noise ( "?). We can observe that the best ' p q pxq function lies around the region where p=1.35 and q=-1.

Article number, page 3 of 5

Cost(j) = µ 2 ||H xj -y|| 2 2 + φ(D xj ) (3.17)
At Iteration j + 1 the Cost(j + 1) will be calculated and compared with the cost at Iteration j: If Cost(j + 1) -Cost(j) then: the algorithm will continue. Else " xj will be returned as the result". This approach of defining the cost function is more effective since we are using the problem itself as means of measuring how close and effective the x is to the real x that we are trying to recover. In the rest of the thesis, we will only consider the approximated solution proposed in Equation (3.16).

Scientific Experiments

The main purpose of this algorithm is for the scanners of I2S in a real-world scenario and application. However, before we move into practical experiments and tests, which will be presented in the next chapter, we need to verify the capabilities of the algorithm using scientific experiments with scientific data. In this section, we will provide various experiments conducted on three images, and we will examine how the algorithm works in diverse conditions.

We begin by conducting a numerical experiment to quantify the loss in accuracy by using approximated proximal operators for image deconvolution/denoising application for the case ϕ(u) = u p . For comparison, the solution is computed by two methods: the exact solution, as proposed in [START_REF] Lorenz | Non-convex variational denoising of images: Interpolation between hard and soft wavelet shrinkage[END_REF]), and its approximation in Equation (3.16) for some p quasi-norms. Although [START_REF] Lorenz | Non-convex variational denoising of images: Interpolation between hard and soft wavelet shrinkage[END_REF]) provides the exact solution, this approach has very high computational complexity in comparison to approximated proximal operators, and the solutions take a few minutes.Whereas in our case, the approximated solutions are calculated within seconds. Since the difference cannot be assessed visually,

in

Notes on parameter β and µ algorithm

As mentioned in section 3.3.1, the parameter β plays a crucial role in the variable splitting technique. However, this role has not been entirely understood [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF]) and not many papers have been dealing with this matter in detail. To gain some perspective of the importance of the subject, here we provide some comments and some -1 (x) and ϕ 2 1 (x) ). This proves that for any φ(u), a constant β will provide a sufficient result. It is worth mentioning that an intermediate value of β can avoid numerical problems of the penalty method mentioned in section 2. A of [START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF]).

Algorithm 2 of (Wang et al. ( 2008))continuation scheme, i.e. gradually increasing β to large values : Algorithm 2 was proposed to tackle the difficulty arising due to the value of µ. In [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF]), µ was chosen according to the level of noise, which is not optimal in terms of peak signal-to-noise ratio (PSNR). This assumption led to small values of µ in the experiments and tests. In this context, it is necessary to have many iterations for the optimization problem to converge, and an algorithm that requires a high number of iterations is not efficient. To tackle this problem, the authors of [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF]) propose algorithm 2 to speed the algorithm. In algorithm 2, the value of β is multiplied by a factor at each iteration. In the next iteration, the optimization will be done with a larger β so that the algorithm can converge faster. We have discovered that this method may not be effective to accelerate the algorithm. Moreover, in higher iterations, β may become too large, i.e. β ≥ 2024, which may cause the algorithm to become too stiff. To tackle this problem, in our method, we choose µ according to optimal PSNR. In this case, the values of µ will have big values (see table 2 of next section).

Under this context, it is unnecessary to use the continuation scheme of algorithm 2 if µ is large. For this reason, our algorithm keeps β fixed. This reduces computational For more information about the algorithm please refer to sections 2.1 and 4 of Wang article. We have conducted an experiment to show our claim. We used the Cameraman image and deconvolution in this experiment using our proposed algorithm with large µ and constant β. We also performed deconvolution under the circumstances of algorithm 2 of [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF]) for the same image. At each iteration the value of β is increased with factor of κ .i.e β i = κ * β (i-1) . Three different values of κ were considered for this experiment. The results are presented on the right-hand side diagram of Figure (3.4).

Upon comparing the results of the two experiments, we can see that fixed β provides a level of accuracy comparable with that of a continuation scheme approach of algorithm 2. 1) , i: iteration value. We note that larger β values converge slower than small ones. Fixed β can provide comparable performance to dynamic if µ is large.

β i = κ * β (i-
Another note concerns the work done in [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF]). This paper considers an extension of the algorithm 2 of [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF]). It considers replacing l 1 with l p regularizers in the anisotropic case. Closed form solutions are applied for particular values of p (p = 1 2 and p = 2 3 ). For other values of p they propose a lookup table (LUT) which is inefficient and impractical. In this work [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF]), the following claims are made which do not match our observations: i) One iteration per value of β can be enough.

ii) Smaller values of β might require more iterations. To assess the correctness of the claims, we have performed the deconvolution on the Cameraman image with µ = 300 and for the case of ϕ(u) = u p . We repeated the experiment with five different values of fixed β. As illustrated in Figure (3.4) left, we can observe different results in contrast with the claims made in [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF]). Our observations provide the following results: i) Several iterations per value of beta are necessary.

ii) More iterations are needed when β is large instead of small.

In this section, we tried to give an analytical observation of the parameters of new and better. These parameters play an important role in optimising deconvolution and denoising problems. However, although they may seem not so important, these parameters are vital to be understood according to the application. The observations that we have provided in this section have not been discussed and analysed thoroughly in recent papers and research. We have tried to provide meaningful insight on these parameters to gain a comprehensive understanding of the role of such parameters in similar optimisation problems.

Experiments of deconvolution performance

In this part, we will perform numerous experiments on scientific data to observe the performance of the deconvolution algorithm. Our proposed algorithm uses various regularizers as prior information for image deconvolution application. It will be shown that we can find new functions and regularizers that outperform previous methods of the image to deconvolution. It will be the first time such new regularizers are considered for image denoising/deconvolution. We will also prove that our algorithm is flexible to be used in various applications under extreme conditions. To begin with, in the following section, we will discuss the measurement parameter that we have used to assess the quality of the images we have obtained through our proposed algorithm.

Image quality metrics

Here we will talk about three main image quality metrics used to assess image quality. Unfortunately, no general quality metric is commonly used. The performance of a quality metric is typically evaluated by its prediction accuracy, monotonicity and consistency. It is also anticipated to reflect the quality grades appointed by independent human viewers. Researchers and scientists are increasingly becoming involved in defining a well defined and accepted image quality metric.

Humans are ultimately the judges of image quality. Moreover, human beings have an inherent ability to determine what "good" or "poor" quality is. The human visual perception system is highly capable of identifying structural information from a scene and hence identifying the differences between the data extracted from a reference and a sample scene. Naturally, it is expected from a machine or computer to provide a similar assessment and perform better on tasks that involve differentiating between a sample and a reference image. For this purpose, various quality metrics have been studied and suggested as a tool for computers to assess the quality of signal and image as close as possible to human criteria. These approaches are primarily based on research and our current understanding of the Human Visual System (HVS).

Mean-Squared Error (MSE):

The MSE quantifies the global difference between a discrete signal and an original (reference) signal; in this sense, MSE is a full-reference metric. It is formulated by:

MSE = 1 N N -1 0 (x i -x i ) 2 (3.18)
where x i is the i-th pixel/sample of the original image/signal, x i is the i-th pixel/sample of the given image/signal and N is total number of pixel/sampled points. This approach is popularly used because of its simplicity. MSE bears a few drawbacks. It effectively considers only global/general variations and neglects localized errors which may happen in a certain region or area of an image. In addition, it is unable to detect structured or patterned errors. Structured/patterned errors include artefacts (an object or pattern that is not originally part of the original image) and distortions that possess a particular, not arbitrary pattern. More specifically, there is an interdependent relationship among adjacent pixels. MSE misses this type of error because the metric examines only absolute differences.

Peak Signal-to-Noise Ratio (PSNR):

Another commonly used and a close relative of the MSE metric is peak signal-tonoise ratio (PSNR). The PSNR is given by: PSNR = 10 log 10 L 2 M SE (3.19) In this formulation, L is the highest possible pixel value in the image, for example, 255 in the case of an 8-bit image. PSNR is a ratio that influences the quality of a signal's representation by comparing its maximum achievable strength to the power of corrupting noise. A higher PSNR value provides a higher image quality. PSNR is a metric that is frequently used to assess digital signal transmission quality. Since PSNR is a variation of the MSE, it still concentrates on pixel-by-pixel comparison and bears the same drawbacks of MSE. PSNR is also more sensitive to additive Gaussian noise.

Structural Similarity Index Measure (SSIM):

Natural image signals are highly structured and patterned: in these images, pixels exhibit strong dependencies; these dependencies carry essential information and data about the structure/pattern of the objects in the visual scene. SSIM was developed to provide a more reliable tool to assess the quality of various types of images. The SSIM is a famous quality metric used to measure the similarity between two images. It was developed by Wang et al. [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]) and is regarded to be in correspondence with the quality of the perception of the human visual system (HVS).

The SSIM is defined as:

SSIM(x, x) = [l(x, x)] α . [c(x, x)] δ . [s(x, x)] γ (3.20)
where α > 0, δ > 0, γ > 0 are parameters to adjust the relative importance of each component. The luminance comparison function is formulated as:

l(x, x) = 2µ x µ x + C 1 µ 2 x + µ 2 x + C 1 . (3.21)
The luminance comparison l function measures the closeness of the two images mean luminance (µ x and µ x). The constant C 1 is added to avoid instability when µ 2 x + µ 2

x is very close to zero.

The contrast comparison function is formulated as:

c(x, x) = 2σ x σ x + C 2 σ 2 x + σ 2 x + C 2 . (3.22)
The contrast comparison function c measures the closeness of the contrast of the two images. Here the contrast is measured by the standard deviation σ x and σ x. Here also C 2 is used to avoid singularity error.

The structure comparison function is formulated as: There are no precise rules for choosing between the SSIM or the PSNR in evaluating the image quality. Some studies have revealed that the MSE and so the PSNR perform poorly in discriminating structural content of images since different types of degradations applied to the same image can result in the same values of the MSE [START_REF] Wang | Mean squared error: Love it or leave it? a new look at signal fidelity measures[END_REF]). On the contrary, some studies have shown that the MSE, and consequently the PSNR, exhibit the best performance in assessing the quality of noisy images [START_REF] Avcibas | Statistical evaluation of image quality measures[END_REF]). and contrast can vary across a scene, we use the local luminance and contrast for our definition.

s(x, x) = σ xx + C 3 σ x σ x + C 3 . ( 3
The system diagram of the proposed quality assessment system is shown in Fig. 3. Suppose and are two nonnegative image signals, which have been aligned with each other (e.g., spatial patches extracted from each image). If we consider one of the signals to have perfect quality, then the similarity measure can serve as a quantitative measurement of the quality of the second signal. The system separates the task of similarity measurement into three comparisons: luminance, contrast and structure. First, the luminance of each signal is compared. Assuming discrete signals, this is estimated as the mean intensity (2) The luminance comparison function is then a function of and . Second, we remove the mean intensity from the signal. In discrete form, the resulting signal corresponds to the projection of vector onto the hyperplane defined by (3) We use the standard deviation (the square root of variance) as an estimate of the signal contrast. An unbiased estimate in discrete form is given by ( 4)

The contrast comparison is then the comparison of Finally, the three components are combined to yield an overall similarity measure (5) An important point is that the three components are relatively independent. For example, the change of luminance and/or contrast will not affect the structures of images.

In order to complete the definition of the similarity measure in (5), we need to define the three functions , , and , as well as the combination function . We also would like the similarity measure to satisfy the following conditions. 1) Symmetry: . 2) Boundedness:

. 3) Unique maximum:

if and only if (in discrete representations, for all ). For luminance comparison, we define [START_REF]2 Book Scanning and Scanners of I2S[END_REF] where the constant is included to avoid instability when is very close to zero. Specifically, we choose (7)

where is the dynamic range of the pixel values (255 for 8-bit grayscale images), and is a small constant. Similar considerations also apply to contrast comparison and structure comparison described later. Equation ( 6) is easily seen to obey the three properties listed above.

Equation ( 6) is also qualitatively consistent with Weber's law, which has been widely used to model light adaptation (also called luminance masking) in the HVS. According to Weber's Since there are no exact rules and regulations about how to assess the quality of the images in scientific papers, also to benefit from the positive aspects of the two image quality metrics; in this thesis, we will use both SSIM and PSNR matrix to evaluate the quality of the images obtained from our deconvolution algorithm. The following section will provide the results of image deconvolution experiments on Lena and Cameraman images to model complex images and signals. Also, we will consider a piecewise constant image to resemble the sparse image and signal.

Experiments

To begin with, we will present a notion of the dataset that we have used in our experiments. It is vital to understand that each image has its unique characteristics. In other words, if we take a look at the signal of each image, they would represent a unique property of their own. To better understand we have depicted the signal representation of three sample images without noise or blur in Figure The usage and application of this particular image will be described in the next chapter.

The images of Lena and the Cameraman contain a more complex signal. Since the Lena image has more structure and pattern, its signal is more complicated than the Cameraman. Figure scanners.The piecewise constant image characteristic is highly similar to text images.

The usage and application of this particular image will be described in the next chapter.

The images of Lena and the Cameraman contain a more complex signal. Since the Lena image has more structure and pattern, its signal is more complicated than the Cameraman. Figure To better assess the capability of the proposed algorithm, we perform deconvolution on our dataset using various regularizers.We report the PSNR and structural similarity index measure (SSIM) obtained with different regularizers after deconvolving for three images: Lena, Cameraman and the piecewise constant. The result of this experiment is presented in Table (3.3).As mentioned before, we perform deconvolution according To effectively observe the variation of PSNR for different types of signals and images, we analyzed a colour map to better understand this phenomenon. As mentioned at the beginning of this section, various images exhibit different properties and characteristics. Some photos tend to be sparse, some Hyper-Laplacian and others are complex signals. To illustrate this variation, we have used two 256 × 256 images: Lena and the piecewise constant image. We conduct image deconvolution for various values of p and q. As shown in Figure (3.8), PSNR values of deconvolution can be observed for q ∈ [-1, 3] and p ∈ [0.1, 1.5] with each range divided into 64 equally spaced points.

In other words, for each image, the deconvolution was done 4096 times, according to optimal µ, each time concerning a different ϕ p q (x) function. The additive Gaussian noise (σ = 0.02) is considered for this experiment. It can be seen that the pattern of PSNR varies from one image to the other. The observation also proves our theory that one general regularizer can not be used for extensive photos and images. The results show that the best ϕ p q (x) function for Lena lies around the region where p = 1.1 and q = 2.9, whereas the best function for Cameraman lies around the region where p = 0.9 and q = 3.

Color Map:

To benefit the appropriate ϕ p q (x) that fits a particular condition of image deconvolution/denoising, we propose a general scheme: The colour map can be acquired once according to the set of images or conditions that the user wants to use the algorithm. Once the appropriate values of p and q are found, the optimal ϕ p q (x) can be used for further experiments and usage. This scheme will guarantee that the images acquired bear the maximum achievable quality.

Visual Demonstration: Our CIFRE PhD aims to implement the algorithms within the scanners and cameras of I2S for the use of clients. In the end, the judge of the quality of scanned images will be the clients. They always prefer to have a product that visually looks good, is sharp and contains as much detail as possible. In consequence, one of the primary purposes of deconvolution is to provide a reconstructed image that can have the optimal quality perception of the human visual system (HVS); therefore, at any stage, aside from the image quality metrics, it is crucial to assess the results visually FIGURE 3.8 -PSNR values after deconvolution using the regularizer ' p q pxq as a function of p and q. The best ' p q pxq function for Lena lies around the region where p " 1.1 and q " 2.9. The best function for Cameraman lies around the region where p " 0.9 and q " 3.

to the results in table 2, One of the key results is that we see that it is possible to find a function among all ' p q pxq regularizers that can outperform regularizers such as l 1 [WYYZ08, ABDF10] and p norm [KF09, ZMZ `13]. A grid search was conducted to find the optimal values of p and q for ' p q pxq. This grid search will be explained shortly. Our experiments show that the variation in performance (PSNR) across values of p and q is similar for different images, convolution kernels and levels of noise. However, some minor variations are observed in each case.

Color Map :

To effectively observe the variation of PSNR for different types of signals and images, we analyzed a colour map to better understand this phenomenon. As mentioned at the beginning of this section, various images exhibit different properties and characteristics. 67 Figure 3.8: PSNR values after deconvolution using the regularizer ϕ p q (x) as a function of p and q. The best ϕ p q (x) function for Lena lies around the region where p = 1.1 and q = 2.9. The best function for Cameraman lies around the region where p = 0.9 and q = 3 with the human eye.

We used Lena and the piecewise constant image for a visual demonstration of deconvolution. We convolved the images with the kernel in previous experiments with additive Gaussian noise with σ = 0.04. The results of deconvolution are shown in Figure (3.10). For both sets of images, the top left is the original, the top right represents the blurred image, the bottom left is the result for the regularizer proposed in [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF]; [START_REF] Zuo | A generalized iterated shrinkage algorithm for non-convex sparse coding[END_REF]) with p : 3/4 and the bottom right is the regularizer of our method: ϕ 3 1.08 (x)( for Lena), ϕ 3 0.92 (x)( for piecewise constant). A significant enhancement can be visually assessed in both cases. The details are much apparent in our proposed regularizers, and the texture in the face of Lena is more visible in comparison to other regularizers. In the case of a piecewise constant, the lines, especially the diagonal lines, are clear and detectable. It will be presented that this will achieve a high TFM score.

It is vital for an algorithm to be flexible and to perform in various circumstances. For the deconvolution to be effective, our algorithm should survive suboptimal conditions.

Chapter conclusion

This chapter provided an extensive review and study of our proposed algorithm that we have used for our image deconvolution/denoising application. We first began by explaining the objectives and motivation for developing a fast and accurate image deconvolution algorithm. We then provided the formulation that we have used in our thesis. We proposed a method that relies on the variable splitting technique; this allows us to solve our optimization problem more effectively and efficiently. It also allows us to separate the problem into two main sub-problems, which can be solved alternatively. This proposition reduces the complexity of our algorithm. It will be shown in chapter 5 that it increases the modularity of our algorithm and allows us to run the algorithm using various parallel computing, which will make the algorithm very fast. We then proposed the general framework that we have used as the main regularizer term in our inverse problem. The potential functions that we have used in this thesis have never been used for any image denoising/deconvolution algorithm before. The potential function allows our algorithm to be flexible and to be able to be used in various applications. To reduce the complexity of our algorithm we have used the proximal operator as an approximation to the potential functions. It was shown that the approximations are very close to the exact solutions and, at the same time, are very fast. We gave an extensive analysis on the parameters of µ and β to provide a general view on the role of these parameters in solving the optimization problem by various experiments and providing the results with plots.

We also demonstrated new claims and observations, which will help better understand the parameters and the role for optimization. Later on, we evaluated three image quality metrics that are popularly used to evaluate the quality of restored images. All three metrics have their pros and cons in which we can benefit according to an application.

In our work to better evaluate our algorithm, we have used all three metrics to evaluate the performance of the image deconvolution algorithm. In the end, we demonstrated the results of various experiments of image deconvolution. The results showed the effectiveness of our proposed algorithm under various circumstances. Furthermore, our results proved to be solid both numerically and visually. It was proven that the potential function could outperform the classical regularizers, which have been extensively studied previously. It was also shown that our algorithm is resistant to changes and alterations in the parameters and conditions of the data. In the next chapter, we will move on to the actual implementation of the algorithm within the scanners and cameras of I2S. We will further provide simulations on real data in a real-world application.

Chapter 4: Practical Results

Introduction

In the previous chapter, we analysed the efficacy of our proposed algorithm with the scientific data. The results showed a good performance of the algorithm both visually and analytically. Previous results were based on scientific data and analysis. However, the main goal of this thesis is to be able to provide a well-defined algorithm for practical purposes. The main objective of I2S is to use our algorithm not just for the CopyBook OS scanner but also to adapt the algorithm for various products such as microscopes and other types of cameras. It is always challenging to take a developed algorithm or idea from a scientific domain and use it in a practical application. Researchers and engineers are constantly trying to build the gap between industry and science. In this regard, from the beginning of the development of scientific ideas, it is important to bear in mind the criteria and objective of a practical application so that from the beginning to the end of implementation, the bridge between science and industry could be built efficiently. Bearing this in mind, we have also aimed to develop our algorithm to be implemented within the scanners of I2S and other products. After verifying scientific data, it is essential to move on to the next step and evaluate the algorithm in real-world conditions. The algorithm may work within simulation software such as MATLAB, but it will face new challenges and obstacles when appropriately and fully implemented within a computer and faced with real-world data. In this chapter, we will begin by explaining the transformation of the algorithm from MATLAB to C++ software and the steps taken to optimise the algorithm. Next, we will provide the experiments by using the empirical data and images taken with the scanners of I2S. It will be shown that our deconvolution proved to be very effective in practice, and also it works for scanners.

We can also implement the algorithm with some adaptation and adjustments in other applications and devices at I2S.

Book Scanning and Scanners of I2S

We live in a digital world, and almost all data is available through digital files with digital devices such as mobile phones, laptops and tablets. One of the main areas that has gained much attention for digitalization is book scanning. Many organizations and companies seek to create a digital archive of their records and documents. For example, organizations such as a university or libraries are looking to preserve vital documents.

Book scanning or digitization is the procedure of converting physical books, magazines, maps, painted pictures into digital media such as images, electronic text, or electronic books (e-books) by using an image scanner with a special camera. Image scanners can be both manual and automated. In a conventional commercial image scanner, the book or the picture is placed on a flat surface, a light is transmitted to the book from above to produce enough luminance for the camera. A glass plate is extended on top of the book to the edge of the scanner, this helps to line up the book's spine easier, keep the book firm in its position and flatten the page for better acquisition. Pages may be turned by hand or by automated paper transfer devices. Once the book has been scanned, a dynamic quality check ensures that every single image is acquired to the highest possible standards. Various books and pages with different conditions are used in the scanning process. In most cases, digitalization deals with fragile, faded and old books which need special handling and care. In this case, it is widespread that specific passages and areas in the image can be hard to read. Therefore, in the background and post-processing, image processing and restoration algorithms, such as our proposed methods, are used to preserve or recover the missing information. Commercial book scanners differ from regular scanners; these book scanners are usually equipped with a high-quality digital camera with light sources on top. The camera is mounted on an extended frame or arm to provide easy access for the user or machine to flip through the pages. Figure (4.1) shows a sketch of a typical book scanner.

Our thesis is aimed at this type of scanner at I2S named CopiBook OS presented in Figure (4.2). This model is designed as A2 format book scanner with a 400×400 DPI resolution. The camera on this particular model is a 71 Mpixels commercially named Our proposed deconvolution algorithm will be implemented on the "EAGLE camera", and the super-resolution scheme with the deconvolution at its side will be integrated into the E-XTRA camera. In the following sections, we will provide the results of our method for images captured with both cameras. In the next section, we will summarise the standard image quality metrics widely used in industrial image processing applications and product verifications.

The Modulation Transfer Function (MTF or FTM)

In many industrial image processing applications, a measurement known as MTF is used to assess the quality of the images captured with various image capturing devices such as cameras. MTF measurements aim to estimate the resolution achieved by an optic system [START_REF] Burns | Ten tips for maintaining digital image quality[END_REF]). Resolution is a measure of how well spatial details are preserved. Measuring two factors as follow is required to define it,: 1) spatial detail and 2) preservation. These spatial detail and preservation parameters are not single measurements but rather a continuum of measurements; this is why a functional curve quantifying them, i.e., the MTF, can be depicted to measure the Simplifies some forms of system analysis bullet.

Compatible with several visual image quality descriptions.

The spatial frequency response (SFR) diagram is a curve that characterizes how an imaging system conserves the relative contrast (spatial detail) of increasing spatial frequency detail.

MTF is a metric quantifying the sharpness of the reconstructed image. The input variable in the horizontal axis of the SFR curve is spatial frequency, which increases to the right. Higher spatial frequencies are interpreted to more finely spaced details. The output response in the vertical axis is the preservation of contrast from an object to a digital image by camera or scanner.

with proposed ISO software solutions for reporting summary measures of camera resolution under ISO 12233 edition 2. 4 Though limiting resolution is a reasonable summary metric for objectively reporting spatial resolution, it does have limitations, particularly in predicting image sharpness. For instance, using a 10% SFR criterion, all of the images in Fig. 1 have the same limiting visual resolution. This is indicated by the loss of text visibility at the fourth text grouping from the top in each case. Notice, however, the remarkable differences in image sharpness between the three images. The rightmost clearly has the best image quality. The higher SFR values at all of the spatial frequencies in the companion graph predict this. This is followed by the middle and left most images with decreasing image sharpness, but equivalent limiting resolution.

Figure 1 illustrates the importance of focusing on the low to middle spatial frequency range for the measuring and predicting of perceived image sharpness and overall quality. This fact has not been lost on the image processing community and is the region where digital sharpening operations are generally beneficial, when applied in moderation. to an image that has been over-sharpened with an image processing algorithm. Though most of the character is readable, their degree of discernability is not the same. This difference is reflected in the shape of the MTF.

Sh

MTF is significantly popular in the industry since it is very practical. With a single measurement, it is able to measure parameters at the same time, and it is also very simple. It is one of the most important benchmarks for imaging performance in the scientific/industrial community because it is fundamental.

General scheme for MTF usage:

To use the MTF measurement effectively, we would like to provide a general guideline that we have used to measure the effectiveness of our algorithm using the MTF measurement. We have depicted the scenarios with the SFR curve measurement in Fig-

ure (4.4). The green curve shows a perfectly standard SFR curve that we would like to achieve after applying our image processing (denoising, deconvolution, interpolation).

If the resulting image bears low contrast and low visibility with blur, the curve would be similar to the red curve in Figure (4.4). In this case, for example, we need to apply the algorithm with more vigorous parameters to be able to move the red curve towards the green. The purple curve is the scenario when the resulting image has been sharpened too much; in this case, the curve needs to move towards the green curve with less intense parameters so that the image will have a standard resolution. The overshooting at the low frequencies of the purple curve is a good demonstration of either over a sharpened image or an artificial sharpening that is not desirable in industrial/scientific applications. We have used this simple but effective scheme in our experiments which will be presented in the coming sections.

(1)

(3) (4) If the sample image is a colour image, the SFR curve is depicted for four red, green, blue and grey channels. Apart from the SFR curve, a star rating system is also included in the MTF software. This system scores the image on a scale of one star to four stars. The sampling image will gain four stars if the SFR curve is very close to the standard perfect SFR curve compared to, and it will gain one star SFR for the curve is far from the usual criteria.

Practical Implementation

This section will begin to present the practical experiments with the images captured with I2S scanners. The images were captured with copybook OS scanners, as mentioned before. Two different data sets were gathered; the first sample of images was collected with the EAGLE camera. The images involved books, texts, newspapers and pictures. The same kind of images was later captured using the E-XTRA camera.

Before presenting the results of deconvolution and super-resolution, we will give a brief explanation of the programming language used for experiments.

Programming Languages

In this thesis, we have used MATLAB and C++ programming languages for our simulations. As shown in chapter 6, we have also used python for some parts of our simulations.

MATLAB programming language:

MATLAB (also regarded as matrix laboratory) is a fourth-generation high-level programming language. It has an interactive environment for numerical computation, visualization and programming. MATLAB programming language is more straightforward than most programming languages. It is regarded as a high-level language because it is closer to the human language than the computer or machine language.

Matlab is beneficial for matrix calculations. We have used MATLAB programming language for scientific experiments presented in chapter 3 since it is rapid and easy to use. It will help us speed the validation of our algorithm. Also, some of the functions were already defined in Matlab libraries. The programming of our algorithm was done in a modular way which makes it easy to understand and easy to use for further modifications. Although MATLAB is very effective for validating calculations and mathematics, it cannot be used for communicating with machines such as computers and microprocessors. Therefore, another programming language such as C or C++ has to be used to implement the algorithm in devices. We have used C++ as the basis of programming ask our algorithm to be implemented in the cameras and computers of I2S.

C++ programming language:

The main goal of our algorithm is for practical application in industry for products like connected image acquisition systems and embedded systems with limited computation capability. C++ is a great tool to help develop a solution that needs fast machine performance. C++ allows the user to program from the base level. This may seem very difficult and complex. However, it will provide the freedom to manipulate the algorithm commands for the machine quickly and efficiently. In C++, the user has to explain everything they intend to do and manipulate the source code; this provides a deeper understanding of how all the parts of the algorithm work. C++ is closer to the hardware level. Because of these reasons, it is used in many compilers as a backend programming language. Some of the benefits of C++ are listed below:

Portability.

Object-Oriented.

Memory Management. engineers gain an idea of the process and effects applied at each step to an image. As we mentioned, apart from analysing the results of image processing algorithms analytically, it is vital for the results to be acceptable for the human visual system. With the image watch plugin, we can easily observe the changes and variations of an image generally and pixel-wise view. Images can be seen in 4 channels explicitly; this helps analyse colour images. A sample screenshot of the Image watch software is presented in Figure (4.6).

ImGui:

Apart from developing the proper algorithm code in C++ with optimization, we have to build a user interface that will be compatible with the code, making it easier for the user to access the algorithm's parameters easily. The image watch software can only be used in debug mode, which is not very efficient when the user is only interested in viewing the results of the image deconvolution algorithm. As mentioned in chapter 3, our algorithm has three main parameters,µ, β and the potential function ϕ p q (x). After considering these factors, we have designed a user interface compatible with C++ called ImGui (also known as Dear ImGui). ImGui is a graphical user interface library for C++. It is designed to enable fast iterations of an algorithm or code and enable the programmers to create content creation tools and visualization tools (https://github.com/ocornut/imgui).

It is particularly suited for fullscreen and embedded applications. A demonstration of our image deconvolution algorithm user interface is given in Figure (4.7).

The user interface is divided into different sections. The values of µ and β are set in the parameter sections. We have a dedicated region to set the parameters of the H or PSF. The initial size, σ values and the shape of the PSF kernel can be selected in this part. Below the parameter selection region, some buttons are set to adjust the view of the image in the window underneath. The original image and the result are depicted at the user interface's bottom part. The user can access the resulting image at any point once the parameters are set. At each moment, by moving the mouse over a particular region of images, a zoomed comparison of the regions is shown as a small pop up window to better evaluate the result (the yellow region in Figure (4.7)). By making a left-click, the pop up can be fixed at any desired position on the page.

The main challenge of programming proposed algorithm:

Even though the OpenCV image processing library and MATLAB provide many functionalities for image processing algorithms, we have faced a particular problem in our programming. Our image deconvolution/denoising and super-resolution techniques are custom-designed with new functions and formulas. For the initial programming of our algorithm in MATLAB, although we benefited from functions already developed in MATLAB, many other Functions had to be written explicitly for our algorithms. On the other hand, we faced a significant challenge in transitioning the codes and functions of MATLAB to C++ language. Most of the functions that we have used in MATLAB or defined in MATLAB did not exist in the OpenCV library. We had to write the equivalent of functions from scratch in the C++ programming language for the transition. Although this proved to be a difficult task, in the end, we were able to provide a fully developed function and C++ for image deconvolution Denoising and super-resolution. In this process, once a new function was developed in C++, it had to be analytically compared to its counterpart in MATLAB to ensure that the results were the same and the function worked adequately. This challenge was one of the most time-consuming processes of implementing our algorithm to be practically used in the industry. Once the initial functional code was developed, the code was optimized and reviewed to make to code less complex and fast. The next part will demonstrate the practical experiments under developed C++ code using the dataset captured with the scanners.

Practical Experiments and Results

This section will demonstrate the results of the image deconvolution algorithm applied to the captured images from the CopiBook OS scanner. We will present the results of images captured using an EAGLE camera for the first part. We will compare the results visually, and we will also evaluate the results with MTF measurement. As discussed previously, scanners face different images from books to newspapers, from maps to pictures; therefore, we need to simulate the same conditions for our image to convolution algorithm. To ensure that our algorithm works for a variety of images, we have included different kinds of images with varying properties of the signal. Suppose the algorithm performs as well for all images. In that case, we can be confident that it will work for different signals, and it will be a robust algorithm to be used under various conditions and circumstances. Now we move on to the MTF measurement of the deconvolution results on the EA-GLE camera. As described before, we place a slant edge target next to our experimental images Figure (??). The captured image will be used as input for MTF software, and it detects the target automatically and provides the SFR curve for three RGB channels.

The scoring is provided in a separate table. The result of the MTF curves is shown in Figure (4.11). The result for the score is given in Figure (4.12).

The MTF curve results show that we are able to improve the last result of MTF for all channels. The Deconvolution curves are almost close to the ideal SFR curve presented in section 4.3. In the latest MTF measurement I2S was able to achieve a 3-star scoring, we are able to improve the scoring and achieve the maximum score, 4-star, with our deconvolution algorithm.

In summary, we have shown the results of our image the convolution algorithm applied to the images captured with the EAGLE camera of I2S. We conducted our experiment on a variety of image samples with different properties. The visual results showed that the images proved significantly; they had a better sharpness, contrast and visibility. The algorithm performed well for all types of data. We also evaluated our results with the well-known MTF measurement.

We were able to improve the latest MTF results of I2S in all four channels of an image. Previously I2S had achieved a three-star scoring for their MTF measurement.

With our algorithm, we have achieved a four-star rating, the maximum achievable score.

In the next section, we provide the results of our experiment on the super-resolved image according to the pipeline that we had shown in chapter 1. It will be shown that the deconvolution also works very well for the new images obtained with E-XTRA camera.

E-XTRA Camera Samples

Just like the EAGLE camera, the E-XTRA camera was calibrated according to the standards so that captured images would have the maximum quality possible. In capture process, we capture four LR images with the translational shifts being (0, 0), (0, 0.5), (0.5, 0), and (0.5, 0.5). According to the implementation explained in section 1. 4 and presented in Figure (1.3), the images are put together. In our work, we refer to the achieved high-resolution image as super-resolved image. After this step, the deconvolution algorithm is applied with its specific parameters to enhance the quality of the super-resolved image. We will also refer to the super-resolved image before deconvolution as a raw-super-resolved image.

Once again, due to the calibration system's confidentiality, the calibration method cannot be explained in this chapter. For the same reason, we could not disclose the parameters and the ϕ p q (x) applied on the super-resolved image. Only visual and MTF results are presented. The result of deconvolution on various samples of raw-superresolved images are shown in Figures (4.13) and (4.14). The images with stars represent the results of deconvolution.

In Figure (4.13), we have shown the results of images with sparse properties. As it can be observed, the space between the lines are well distinguished, and the numbers are more visible. For example, the parallel lines in the image with the number "8.0" have been explicitly recovered. In the raw-super-resolved image, the number of parallel lines could not be counted easily, whereas the quantity of the lines could be estimated easily. This is a significant enhancement of image quality.

In Figure (4.14), we have shown the results of images with complex properties. In all results, the patterns, the lines and the structure of images are appropriately restored. For instance, considering the face of the girl, the tiny spots on the skin are distinctly visible with detail. The pattern on the eye and eyebrows areas is noticeably visible compared to its raw counterpart. As another example, the boundaries on the structure of the leaf image is enhanced and extracted. Image contrast and sharpness have been significantly improved in all retrieved results, and the borders and edges are more detectable. The results in both figures presented above show that our proposed deconvolution is efficient for recovering raw-super-resolved images. Once again, the algorithm performs well for all types of samples proving the robustness of the algorithm in a different applications.

The sizes of the images captured with CopiBook scanners with E-XTRA camera are 20000 × 14192 pixels. These images are double the size (thanks to our super-resolution technique) of the images taken with the EAGLE camera. The Runtime of the experiment for images of 20000 × 14192 pixels for two different configurations is 10.5 seconds and The MTF curve for the raw-super-resolved image shows a very poor result. However, the result changes significantly after applying the deconvolution algorithm in 800 DPI.

After the competition, the curves are very close to the ideal desired MTF curve. The curve improves further once the image is considered in 600 DPI. We have achieved the best MTF curve possible for the new super-resolved image, and this achievement is the first in the industry and the applications of I2S. Furthermore, when we look at the MTF score in Figure (4.17), we have gained a four-star rating for both 800 and 600 DPI cases.

The amount after the noise has been reduced by 33% for 800 DPI and 66% for 600 DPI. These results show the significant effect of our deconvolution algorithm on the super-resolution scheme. We have been able to meet the industrial standards of image quality. We will now present the image quality and the resolution achieved through our super-resolution scheme compared to the images used at I2S before implementing the scheme. In Figure (4.18), image (a) shows the classic images of I2S before the super-resolution scheme. The same image is shown after the super-resolution scheme for 800 DPI,(b), 600 DPI,(c), cases. Their results speak for themselves; we can see an outstanding improvement in image quality. To be more precise, none of the numbers and the lines corresponding to them are visible in the classic image (a). Also, in the traditional image, (a), the quality is not so well, and each individual pixel is visible.

We cannot see any individual pixels in the results obtained through the super-resolution scheme (deconvolution included). The number itself can be seen clearly for both numbers of "7.1" and "8.0", and the parallel lines can be distinguished easily. We have to note that the space between the lines has been recovered ideally, which is representative of the accuracy of our algorithm. Background noise has also been reduced drastically, and artefacts have been removed.

Microscopic camera sample

Developing a robust and flexible deconvolution algorithm has always been a critical factor from the beginning of this thesis. This is why we have used the potential function The properties of medical images, such as microscopic images, are different from sparse or complex signals that we had analysed before; therefore, several tests with different conditions were conducted on different samples to find the optimal parameter value µ. We obtained the colour map for microscopic images as explained in section 3.5.2. We conducted image deconvolution for various values of p and q. PSNR values of deconvolution were observed for q ∈ [-1, 3] and p ∈ [0.1, 1.5] with each range divided into 32 equally spaced points. In other words, for each image, the deconvolution was done 1024 times, according to optimal µ, each time concerning a different ϕ p q (x) function. The additive Gaussian noise (σ = 0.04) is considered for this experiment. Usually, the noise is high for microscopic images. The parameters, the PSF kernel and exact ϕ p q (x) function for this experiment could not be disclosed because of confidentiality. The result of image deconvolution/denoising od different samples is displayed in Figure (4.20).

The results for both samples show a significant enhancement in image quality.

According to our medical collaborators for the ATLIS project, it is imperative for practitioners to be able to see the patterns in the samples clearly. The cells should be visible clearly for better diagnosis. Our proposed algorithm distinguishes the patterns with detail. It is also able to recover the smallest cells. Noise in both samples has been removed or reduced.

Terahertz imaging

Terahertz (THz) spectroscopy and imaging is a contactless, non-invasive technology emerging as a tool for analyzing material's characterization. It has been used for contact-free measurements of metals, semiconductors, 2D materials, and superconductors, pharmaceutical analysis, making it a popular research area for fundamental science, security, and medical applications.

Terahertz imaging has been developing very fast recently. Among various imaging

Commercialization of Super-resolution and Deconvolution

As shown in previous sections, the deconvolution algorithm and the super-resolution scheme showed promising effectiveness for various applications of I2S. For the first step, the deconvolution algorithm was considered a software upgrade for the EAGLE camera.

The results are validated with international standards of image quality and analysis, and after the validation, our deconvolution algorithm was commercially presented to customers under the commercial name "BOOSTPIC" . The algorithm was presented in conferences and workshops to clients and partners as a newly developed feature for CopiBook OS scanners. A screenshot of the commercial advertisement of the "BOOSTPIC" algorithm is presented in Figure (4.23). Up to the point of writing this thesis, the integration of the super-resolution scheme with deconvolution is at the final stage, and the new method will be mounted on the E-XTRA camera and later on, it will be commercially presented to customers and clients as a new product of I2S.

https://www.i2s.fr/fr/page/digisolution

Chapter conclusion

This chapter provided the practical results of implementing the deconvolution and the super-resolution scheme on real-world data. We first began by introducing the I2S book scanner currently being used at I2S, which was the subject of our scientific research. We are considering two different cameras off I2S in our study; the EAGLE camera and the E-XTRA camera. We briefly described how book scanners work, their applications, and why it is vital to have good quality image processing algorithms for such devices. We then provided the MTF as a standard metric for measuring the quality of images in various image processing applications and imaging devices for industrial. The scientific experiments of our algorithm were implemented in the MATLAB programming language. However, to implement the algorithms in computers and devices, we transferred the codes to the C++ language. This caused a particular challenge since most of the codes and functions we were using were not developed in the standard libraries of C++ and also OpenCV library. We used special libraries designed for OpenCV to build a friendly user interface to help the operator use the algorithm efficiently. In section 4.5, we provided our practical experiments on four different devices of I2S. We first began with the EAGle camera, which was currently being used in scanners of I2S. For this particular camera, we have only implemented the deconvolution algorithm. Visual results represented a significant enhancement in the quantity of the images captured. The FTM measurement gave a four-star rating for the classical EAGLE camera images deconvolution. Moving on, we provided the results of the super-resolution scheme on the images captured with the E-XTRA camera of I2S.

The proposed super-resolution scheme allowed the camera to produce images up to 800 DPI resolution. The visual assessment of the images showed significant achievement in quality. It was shown that we could produce images with more substantial details and accuracy with our super-resolution scheme. The MTF measurement presented a four-star grading for the 800 DPI and 600 DPI images; this is the first in the industry.

Our deconvolution, with parameter adaptation, can also be implemented on Microscopy and Terahertz cameras. For both cases, the algorithm represented a notable improvement in image quality. The proposed deconvolution algorithm has been commercialised for clients, and the super-resolution scheme is integrated and commercially presented as a new product. This chapter proved we have been able to achieve our main objective goals of developing an effective algorithm for various practical and industrial uses.

of the interstellar medium as accessible in the Herschel data, and study the quality of the singularity spectra obtained on filtered and de-beamed data. Since the inverse problem is the same as in chapter 3, we will exclude explaining the minimisation algorithm based on the approximal operator and the Fourier transform. Considering that the dataset is changed, it is essential to discover which values of p and q are optimal for this particular dataset.

Choice of the potential function

To identify the best potential function for Herschel image deconvolution (beam effect reduction), we conduct an experiment in which we evaluate the performance of all possible potential functions and measure their performance in terms of PSNR. In this experiment, we take a sample from simulation data and convolve the image with beam effect, H. Next, we use our algorithm with various potential function configurations to obtain different results based on the corresponding regularizer. To illustrate the variation in performance, we show in Fig. 5.1 PSNR values for q ∈ [-1, 3] and p ∈ [0.1, 2] with each range divided into 64 equally spaced points. For this experiment, we have used one 1024 × 1024 sample with additive Gaussian noise (σ = 0.005). We can observe that the best ϕ p q (x) function lies around the region where p = 1.35 and q = -1. It should be noted however, that such an "optimal" value of p = 1.35 is resulting from a PSNR computation using a fixed additive Gaussian noise and that, for real observation maps, there is an interval interval of acceptable p values, typically in the range [1. 3, 1.6] as will be shown in the experiments. For some maps, a value of p too close to a low value of p = 1.2 may result in oversmoothing small-scale gradients. Consequently, in the following sections, all our experiments on Herschel maps and simulation data are conducted setting q = -1 and p in the interval [1. 3, 1.6].

Herschel data and results

In this experiment we consider observation maps from the Herschel Gould Belt Survey ( [START_REF] André | From filamentary clouds to prestellar cores to the stellar imf: Initial highlights from the herschel gould belt survey[END_REF][START_REF] Pilbratt | Herschel space observatory-an esa facility for far-infrared and submillimetre 135 BIBLIOGRAPHY astronomy[END_REF]), acquired by the SPIRE instrument [START_REF] Griffin | The herschel-spire instrument and its in-flight performance[END_REF]) at high spatial resolution and dynamical range. We focus on the 250 µm

Processing

To evaluate the results obtained in noise reduction and de-beaming as contemplated in this study, we use the methodology presented in detail in [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]) from which we use the notations and to which we refer the reader. In particular, the statistics of turbulence in the ISM are obtained by the computation of the singularity spectrum in a microcanonical formulation as developped in the mentioned reference, and the complexity of the ISM is studied through the display of the singularity exponents, computed with the help of a local correlation measure; the core of the approach is based on the notion of reconstructibility of a complex signal which is shown to be computed locally though the hypothesis of the identity of two sets: the most singular manifiold and the most unpredictable manifold. The former is defined as the set of the most singular points in an observation map, and the latter is defined as the set of points which lead to a perfect reconstruction. If s is an observation map acquired at minimal spatial resolution r 0 = 2 -n , we define a local correlation measure H(µ n , x, r 0 ) according to formula (B.7) in [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]) with µ n referring to a discrete gradient measure, and x is a 2D point location in the domain of the observation map s. The singularity exponent at x is the quantity defined by Equation (B.5) of the previous reference. The singularity spectrum is then shown to be given by the relation: .3) with ρ n being the probability distribution function of the singularity exponents h(x) at resolution r 0 = 2 -n and ρ 0 is the most probable event i.e. the maximum of ρ n .

D(h) = 2 - log(ρ n (h)/ρ 0 ) log(r 0 ) . ( 5
The singularity spectrum h → D(h) appears in the limiting behavior of the probability distribution functions:

ρ n (h) ∼ c n r -D(h) 0 (n → +∞) (5.4)
with c n > 0. As long as the singularity spectrum is computed with enough precision, it is possible to distinguish different physical processes responsible for the turbulent statistics in the observation map s such as log-normal or log-Poisson processes. Besides, it is fundamental to notice that before the application of the multifractal formalism, scale invariance of the observational map must be checked, as described in section 4.2 of [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]). Such scale invariance has been checked for the observation maps studied in this work.

Musca

We start by showing, on the Musca observation map, the superiority of noise reduction methodology presented in this work over the L 1 -L 1 algorithm introduced in [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]). We refer to [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]) and the references therein for a presentation of the Herschel 250 µm observation map of the Musca filament. In Figure (5.2) we first recall, in the top panel, the map of the singularity exponents of the unfiltered observation map and, in the below panel, the resulting singularity exponents after having applied the noise reduction with p = 1.5, q = -1, λ = 0.1. Since these singularity exponents display the transitions within a signal, one could raise the question whether the filamentary coherent structures shown could equally be rendered by the gradient's norms of the Musca map. It turns out that the gradient's norms themselves possess a high dynamical range in such a way that the logarithm of the gradient's norms must be considered instead: if s is the signal of an observational map, log( ∇s (x)) can indeed be considered as a simple local correlation measure and consequently, an approximation of a singularity exponent at x = (x, y). In practice we compute the gradient of the observation map s at x = (x, y), ∇s(x) = ∂s ∂x (x), ∂s ∂y (x)

in Fourier space in the usual manner:

-ix ∂s ∂x = F -1 if 1 ∂F ∂f 1 (s) -iy ∂s ∂y = F -1 if 2 ∂F ∂f 2 (s)
(5.5)

with f = (f 1 , f 2 ) is the frequency vector, (x, y) denotes spatial coordinates, F is the Fourier transform, F -1 the inverse Fourier transform and i is the imaginary unit.

The log of the gradient's norm can then be considered as a local correlation measure by defining:

h(x) =          log ∇s(x) ∇s /log(l) : if ∇s(x) ∇s > ε log(ε) log(l) : else, (5.6)
with l = (l x l y ) -1 2 , l x , l y are the x and y lengths, in pixel units, of the array of values associated to the observational map s, ∇s is the average of the gradient's norms over the spatial domain of the observation map and ε is a threshold value, chosen here ε = 10 -30 .

In Figure (5.4) we display the resulting gradient's norms, as defined by Equation 5.6.

The resulting striations and coherent filamentary structures resemble a lot those of Figure 5.3 but it turns out that these values do not encode the same statistics of turbulence.

Indeed, if we use Equation 5.6 to compute the corresponding singularity spectrum as defined by Equation ( 5.3) on Musca raw and with noise reduction, we obtain the two curves shown in Figure (5.5), right. The resulting curves are too dissimilar to be labelled as correct singularity spectra. In the same figure, left, we show in blue and red the singularity spectra computed using the values of h(x) from formula (B.7) in [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF] and in green the curve computed with h(x) defined by 5.6 from the data of Musca with noise reduction (p = 1.5, q = -1, λ = 0.1). Consequently although the log of gradient's norms reveal interesting filamentary structures on observation maps with noise reduction, we must keep the values of the singularity exponents given formula (B.7) in [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF] for both visualization and the computation of correct statistics of turbulence.

Chapter conclusion

This chapter tried to extend the deconvolution algorithm to the astronomical data set.

We conducted primary experiments on Herschel observation maps to improve the quality and reduce the noise of the images for better observation. We found a new potential function that could provide better results than other related work. Although this work is in the preliminary stages, it has shown true potential for this type of image for better image analysis. More experiments and simulations are underway to further develop and adapt our image deconvolution/denoising algorithm for Herschel observation maps. for applications ranging from geophysical, microscopy, medical imaging and low light image processing. We will provide a brief overview of these evolving aspects.

Medical imaging:

The applications of Machine Learning (ML) has been increasing fast in the medical imaging field (

Computational photography:

The purpose of computational photography is to create visually engaging images that are as credible as possible to the scenes captured. These factors make deep learning a good candidate for computational photography reconstruction applications. For example, deep learning provides suitable solutions for extreme low-light imaging [START_REF] Chen | Learning to see in the dark[END_REF]). Deep learning also helps to estimate the depths of different objects in an observed scene from a picture [START_REF] Godard | Unsupervised monocular depth estimation with left-right consistency[END_REF]). Currently, deep learning

Plug-and-Play Priors

In this section, we talk about the idea of Plug-and-Play (PnP), a framework for the incorporation of advanced prior models. The principal goal of PnP is to allow the cost minimization models in classical minimization methods to be fused with learning models, such as DNNs, that are designed to generate outputs from inputs. First, we begin by giving the motivation and introduction to plug-and-play, then we will talk about how we use this method in our deconvolution algorithm.

Motivation and Intuition

As we have previously talked, inverse problems represent a very broad class of important problems in imaging in which the goal is to recover some unknown , x, from a set of noisy and degraded measurements, y. The classical formulation of this problems is:

y = Hx + n (6.1)
where y ∈ R m is the observed data, x ∈ R n is the (vectorized) underlying image to be inferred, H is the observation (blur) matrix, and n is assumed to be Gaussian zero-mean with known variance σ 2 . The maximum a posteriori (MAP) estimate (under prior p(x)) and regularized ML estimate are general methods for solving such inverse problems. Let p(y|x) be the conditional probability density function (pdf) of the measurements y given x, and p(x) be the pdf of the unknown, then the MAP estimate of x can be expressed in the form:

x = arg max x p(y|x)p(x) (6.2) x = argmin x {-log p(y|x) -log p(x)} (6.3)
which is written in the form:

x = argmin x {f (x) + h(x)} (6.4)
In regularized based approach inversion, f (x) can be any data fitting expression such as the potential function introduced in chapter 3, and h(x) can be any stabilizing function that penalizes irregular behaviour in the solution. This optimization framework for regularized approach has performed well. However, this approach is not modular. So the optimization of both terms must be done together. Image denoisers such as block-matching and 3D filtering and convolutional neural network denoisers have demonstrated promising improvements in image denoising performance.This is possible due to increasing complex image operations.

Question. Is there a way to fuse the traditional models of regularized inversion with modern denoising algorithms?

Plug-and-play (PnP) methods try to answer this question by providing a framework for fusing traditional iterative algorithms for imaging inverse problems with modern image denoisers, used in a black-box fashion.

Image denoising:

The special case of image denoising occurs when H = I. In this case the observation y consists of x plus additive Gaussian white noise. The MAP cost function is given by:

f (x, y) = argmin x 1 2 ||x -y|| 2 2 + φ(x) (6.5)
Recall definition of proximal operator from definition 2.1 assuming λ = 1, we get:

prox f (y) = argmin x 1 2 x -y 2 2 + f (x) (6.6)
When comparing Equation (6.6) with (6.5), the proximal operator could be seen as the MAP solution of a denoising problem. In this context, the argument of prox f is the noisy data; the noise is Gaussian i.i.d. with unit variance. The prior is p(x) ∝= exp(-f (x)). This interpretation is at the root of the method plug-and-play. Here, the denoiser acts as a regularizer for the reconstruction scheme and ensures good image reconstruction quality at each iteration of the algorithm. The final reconstructed image matches the measurements and satisfies the prior defined by the denoiser. In simple words, the key idea behind plug-and-play is to use denoisers, such as BM3D [START_REF] Dabov | Image denoising by sparse 3-d transformdomain collaborative filtering[END_REF]) or any other denoiser, in place of proximal operators in iterative traditional optimization algorithms. As mentioned before, this idea considers combining learning methods with classical methods simultaneously. One of the key advantages of this method is modularity. It will enable us to use any type of denoiser suitable for our application in the place of the proximal operator. This method is not superior to the proximal operator. It only provides the freedom of using a learned denoiser that, according to the application, may perform better than the proximal operator. Since it is modular, not many changes are necessary to be made for the classical algorithm. Any denoiser is plugged into the algorithm, and the whole algorithm is played so that the optimal result is obtained. This is how this particular method gets the name plug-andplay.

The idea of plug-and-play was first proposed in [START_REF] Venkatakrishnan | Plug-and-play priors for model 139 BIBLIOGRAPHY based reconstruction[END_REF]).

Since 2016 it has gained much attention in different applications. Most of the plugand-play methods are built on the classical ADMM algorithm [START_REF] Teodoro | Image restoration and reconstruction using targeted plug-and-play priors[END_REF]), however has never been considered in the classical alternating minimization method. We will use plug-and-play in our alternating minimization image deconvolution algorithm for the first time.

Proposed plug-and-play approach

In this part, we will demonstrate how we use the plug-and-play approach in our deconvolution scheme. Recall that for the deconvolution problem of Equation ( 6 We have multiplied and divided Equation (6.15) with p r (y) to have Equation (6.16).

p r (y) = N (y; 0, C r + σ 2 I). The Minimum Mean Square Error (MMSE) of x from an observation y can be defined as the posterior expectation:

x = E[x|y] = xp(x|y)dx (6.17)

We insert p(x|y) from Equation ( 6.16) into Equation (6.17):

x = 1 p(y) K r=1 ω r p r (y) x p(y|x)N (x; 0, C r ) p r (y) dx (6.18) Within Equation ( 6.18) we define xr as follow:

x = 1 p(y) K r=1 ω r p r (y) x p(y|x)N (x; 0, C r ) p r (y) dx xr (6.19) xr is the MMSE estimate of x from y, using a single Gaussian prior (the m-th component of the GMM).

Lemma 6.1. Linear MMSE estimator for linear observation process

♥

In observation of a linear process: y = Hx + n, where H is a known matrix and n is random noise vector with zero mean, E {n} = 0. with the assumption that cross-covariance matrix is zero, C XN = 0. We have the following:

E {x} = x E {y} = H x C XY = C X H T C Y = HC X H T + C N (6.20)
The expression for linear MMSE estimator of x is as follow:

x = C XY C -1 Y (yȳ) + x (6.21)

Using the information above we get:

x = C X H T HC X H T + C Z -1
y -H X + x (6.22)

For more information about Linear MMSE estimator please refer to ECE 830

Course Notes of John A. Gubner. We take Equation ( 6

Experiments

This section will present the results of image deconvolution based on the plug-andplay approach. We will use the dataset generated with I2S images captured with an EAGLE camera for the test.

Our proposed algorithm in chapter 3 will become as follow after the modifications:

The cost criteria are the same as in chapter 3. One of the advantages of this algorithm is modularity. We have shown that sub-problem (S1) can be interpreted as a denoiser. In this regard, step (2) of the algorithm can be replaced by any black box. The Algorithm 3: Image deconvolution using Plug-and-play priors Input: initialize x, H, β, i = 0, j = 0 Output: x 1: while not converged do 2: Compute u j with GMM denoiser according to (6.24); 3: Compute x j using u j according to (6.9) ; 4: Compute Cost ( j ) ; 5: j ← j + 1 6: end while 7: return x i ; 8: return Cost; number of iterations may differ according to the type and application of the denoiser being used. In a simple scheme, the algorithm performs denoising and deconvolution simultaneously. Most denoising is done with step 2 with the GMM denoiser, and step 3 performs deconvolution.

We now present the visual results of the results obtained with the methods described Still, the recovery is noticeable. We have also shown the zoomed version of the signals for the original, noisy and recovered image in Figure (6.4). We can see that the GMM based deconvolution/denoising algorithm has recovered the signal well; the black curve is close closer to the blue curve compared to the red curve. The restored image's signal representation shows that although the signal is recovered, some noise remains. This can be solved with more algorithm iterations or more training data. Test 2: For the second experiment, 25% of images were kept clean, for 25% we added noise with σ = 0.015, for 25% we added noise with σ = 0.025 and for the rest we added noise with σ = 0.035. We used the same complex test image as the previous experiment. For this test image we added with σ = 0.06. The result of denoising is Here we would like to provide some points according to the conditions where the deconvolution performance can be improved:

1-if the significance of deconvolution is not visible enough, a more substantial and bigger PSF can be chosen to help the deconvolution. This parameter shows its effect on the subproblem (S2), where we try to recover the x using the Fourier transform. 2-if the significance of deconvolution is not visible more iterations could be needed to allow the algorithm to perform the deconvolution. Please note this will add to the processing time of the algorithm, which is not always desirable. 3-large values of µ can help improve the deconvolution's performance. µ is used in sub-problem (S2); if a more significant value is chosen, more weight will be put to this sub-problem, enhancing the effect of the kernel and deconvolution.

In general, this is the first time that the idea of plug-and-play is used in alternating minimisation the scheme; therefore, more experiments and exercises are needed to define how we can approach denoising and deconvolution more efficiently.

Chapter conclusion

This chapter introduces the idea of plug-and-play priors for image denoising and deconvolution algorithms. Plug-and-play is a branch of recently developed learning methods in inverse problems. We briefly reviewed deep learning methods and algorithms such as CNNs and DNNs. Since 2015, many deep learning methods have become the focus of research and development in various applications such as medical imaging, computational photography, microscopy, and image processing. We gave our motivation and intuition of plug and play priors. We show that this algorithm can be used within the classical iterative methods. One of the key points of the plug and play is its modularity.

In a classical, alternative method, the proximal minimisation term could be replaced by a black box such as an image denoiser. In the next step, we gave the idea of the GMM based denoiser used in our algorithm. In the last section of this chapter, we

provided various examples and experiments of applying plug-and-play in our classical image deconvolution/denoising scheme provided in chapter 3.

The results show that the plug-and-play priors can perform reasonably well for denoising. However, the algorithm does not perform notably well for deconvolution compared to the classical results obtained in chapter 3. Moreover, supplementary experiments and tests should be carried out to precisely examine the role of parameters such as training data, the Kernel, µ and β in alternating minimisation scheme. Our thesis was the first time that plug-and play-priors were used in an alternating minimisation scheme. The results showed a promising potential for plug-and-play in classical approaches for deconvolution/denoising. This part of the thesis is still under integration to be implemented within the devices of I2S.

1. 1

 1 The concept of multi-frame super-resolution. The images on the left side represent the LR images of the same scene with sub-pixel displacement thus the HR image (the picture on the right side) can be acquired by fusing the complementary information with SR methods. . . . . . . . . 1.2 Visual demonstration of 0.5-pixel displacement of the scene captured . . 1.3 I2S implementation of super-resolution . . . . . . . . . . . . . . . . . 1.4 Reconstruction of the super-resolved image: by interlacing the original pixels and the shifted image . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Super-resolution scheme of this thesis with 4 RGB images. . . . . . . . 2.1 Demonstration of an Inverse problem application in image processing. . 2.2 3D demonstration of a sample PSF and its matrix value. . . . . . . . . . 2.3 Relationship between sparsity and lp norm . . . . . . . . . . . . . . . 2.4 Evaluating a proximal operator at various points. . . . . . . . . . . . . 3.1 Potential function for different values of p and q. . . . . . . . . . . . . 3.2 Top) Exact and approximate Proximal Operator of ϕ p -1 (x) = |x| p with β = 1. The approximation has great accuracy in comparison with exact solution. Only small variations can be observed around the threshold value. (Bottom) zoom at the threshold region. . . . . . . . . . . . . . 3.3 PSNR of images recovered in experiment 3.1 for different β, Lena (above) and Cameraman (below). The performance does not improve once β reaches a certain value. . . . . . . . . . . . . . . . . . . . . . . 3.4 (Left) Evolution of the cost function for deblurring Cameraman, µ = 300, p = 2/3. Different values of fixed β. (Right) Comparison between

3. 6

 6 Data set used for Experiments. . . . . . . . . . . . . . . . . . . . . . . 3.7 Two demonstration of signals: (a) Signal representation of original Images. (b) Representation of the original signals in the presence of blur and noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . . 4.1 Graphical demonstration of a typical manual book scanner scheme . . . 4.2 CopiBook OS scanner of I2S . . . . . . . . . . . . . . . . . . . . . . . 4.3 SFR plots of three images with different resolutions taken from (Williams and Burns (2008)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 The meaning of different SFR curves. . . . . . . . . . . . . . . . . . . 4.5 Slant Edge Target also known as QA62. . . . . . . . . . . . . . . . . . 4.6 Screenshot of the Image watch user interface. . . . . . . . . . . . . . . 4.7 Custom designed user interface of proposed deconvolution. . . . . . . . 4.8 Result of deconvolution on EAGLE camera sparse images. . . . . . . . 4.9 Result of deconvolution on EAGLE camera complex images. . . . . . . 4.10 Slant edge target as input image . . . . . . . . . . . . . . . . . . . . . . 4.11 MTF curve for EAGLE camera samples. . . . . . . . . . . . . . . . . . 4.12 Scoring table for MTF measurement for EAGLE camera . . . . . . . . ii LIST OF FIGURES 4.13 Result of deconvolution on raw-super-resolved images with sparse property. 80 4.14 Result of deconvolution on raw-super-resolved images with complex property. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.15 Slant edge target as input image for E-XTRA camera. . . . . . . . . . . 82 4.16 MTF curve for E-XTRA camera samples for two mode of 800 DPI and 600 DPI. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.17 Scoring table for MTF measurement for E-XTRA camera . . . . . . . . 84 4.18 Quality comparison of classical images of I2S and new super-resolved Image . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.19 ATLIS Modulo microscope for cell analysis . . . . . . . . . . . . . . . 86 4.20 Deconvolution/denoising result for two samples of microscopic imaging accquired with optimal ϕ p q (x) potential function. . . . . . . . . . . . . 87 4.21 TZcam terahertz camera of i2S . . . . . . . . . . . . . . . . . . . . . . 88 4.22 (a) THz image of a 1951 USAF resolution test chart at 0.97 THz acquired with TZcam and x0.25 magnification lens. (b) Deconvolved image (c) Plots of the Line Pairs pixel in grey values of the group-1 element 2 of the chart of original THz image(d) deconvolved image.(e) Sample industrial application at 2.5THz. (f) Deconvolved image of (e) . . . . . 89 4.23 Commercial representation of BOOSTPIC . . . . . . . . . . . . . . . . 90 5.1 PSNR values after deconvolution using the regularizer ϕ p q (x) as a function of p and q. The region around p=1.35 and q=-1 provides the best performance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.2 Singularity exponents of the Musca observation Herschel at 250 µm.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.1 Sample of training images for sparse data set . . . . . . . . . . . . . . . 114 6.2 Result of denoising using clean training images . . . . . . . . . . . . . 115 6.3 Result of denoising using mixture of noisy and clean training images . . 116 6.4 Result of denoising using mixture of noisy and clean training imageszoomed signal view . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.5 Sample of training images for complex data set . . . . . . . . . . . . . 117 6.6 Result of denoising using clean training complex images -zoomed view 117 6.7 Result of denoising using mixture of noisy and clean complex training images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 6.8 Result of deconvolution using a mixture of noisy and clean complex training images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.9 Signals of deconvolution using a mixture of noisy and clean complex training images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 List of Tables 1.1 Categories of SR algorithms -Classification of SR algorithms (SISO, MISO, and MIMO) and some potential applications for each class. . . . 3.1 Potential function and its corresponding functions according to the values of p and q. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Accuracy of approximation.The difference between approximated and exact solutions is negligible. . . . . . . . . . . . . . . . . . . . . . . .

3. 4

 4 PSNR for different kernel sizes for different regularizers. The proposed regularizer ϕ p q (x) provides a modest PSNR enhancement. The values in bold correspond to the best performance. . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: The concept of multi-frame super-resolution. The images on the left side represent the LR images of the same scene with sub-pixel displacement thus the HR image (the picture on the right side) can be acquired by fusing the complementary information with SR methods.

  (El-Khamy et al. (2006)), El-Khamy et al. recommended to first register multiple LR images in the wavelet domain, afterwards fuse the registered LR wavelet coefficients to attain a single image, and in the final stage to perform interpolation to have a higherresolution image.

  the general variational regularized SR framework. The formulation is equivalent to maximum likelihood (ML) estimation without the regularisation term. The MAP methods incorporate the prior constraints of the image and obtain the results by maximizing the cost function of the posterior probability. Regularization methods are popular because of their flexibility with edge-preserving priors and joint parameter estimation. Relatively, Bayesian estimation is used when the posterior probability distribution of the unknown parameters, instead of the specific parameters, is estimated.

2

 2 norm have their advantages and disadvantages. Some researchers have therefore employed improved techniques for the data fidelity term (Shen et al. (2015c); Suo et al. (2011); Yue et al. (2014b); Zeng and Yang (2013); Zhang et al. (

Figure 1 . 2 :

 12 Figure 1.2: Visual demonstration of 0.5-pixel displacement of the scene captured

Figure 1 . 4 :Figure 1 . 5 :

 1415 Figure 1.3: I2S implementation of super-resolution

Figure 2 . 1 :

 21 Figure 2.1: Demonstration of an Inverse problem application in image processing.

Figure 2 . 2 :

 22 Figure 2.2: 3D demonstration of a sample PSF and its matrix value.

  (a) x = 0 if and only if x = 0. (b) λx = |λ| x for all λ and x. (c) x + y ≤ x + y for all x and y (triangle inequality). If only (b) and (c) hold the function is called a semi-norm. And if only (a) and (b) hold the function is called a quasi-norm.

Figure 2 . 3 :

 23 Figure 2.3: Relationship between sparsity and lp norm

Figure 2 . 4 :

 24 Figure 2.4: Evaluating a proximal operator at various points.

Fig. 3 : 5 Figure 3 . 1 :

 3531 Fig. 3: Potential function for di↵erent values of p and q.

( 3 .

 3 16) Here (x) + denotes the positive part of x. As shown by this equation, it is computed by a simple closed-form expression for values bigger than η 1 β and set to zero if smaller.

(

  S1). As outlined in Algorithm 2 below, we minimize Equation (3.10) by solving the u and x sub-problems separately until the algorithm converges. We keep β constant (β= 256) for our academic experiments, making our algorithm simpler and reducing computational complexity. In step 2 of the algorithm, the solution of subproblem (S1) of Equation(3.11) is calculated using the approximate proximal operator of Equation(3.16) In step 3, the u obtained from step 2 is used to calculate the x of sub problem (S2) of Equation(3.11) with Equation(3.14). This process is repeated until the algorithm is converged.

Fig. 1 :

 1 Fig. 1: (Top) Exact and approximate Proximal Operator of ' p ´1pxq " |x| p with " 1. The approximation has great accuracy in comparison with exact solution. Only small variations can be observed around the threshold value. (Bottom) Zoom at the threshold region.

Figure 3 . 2 :

 32 Figure 3.2: Top) Exact and approximate Proximal Operator of ϕ p -1 (x) = |x| p with β = 1. The approximation has great accuracy in comparison with exact solution. Only small variations can be observed around the threshold value. (Bottom) zoom at the threshold region.

Algorithm 2 :

 2 Image deconvolution using Approximate Proximal Operators Input: initialize x, H, β, i = 0, j = 0 Output: x 1: while not converged do 2: Compute u j according to (S1) of Equation (3.11) for fixed x; 3: Compute x j using u j according to (S2) of Equation (3.11) ; 4: Compute Cost ( j ) ; 5: j ← j + 1 6: end while 7: return x i ; 8: return Cost; In this manner, we define the cost function of our algorithm in the following: Definition 3.1. definition of cost function ♣ Definition: In algorithm 2, for a given y ∈ R, H ∈ R and recovered x ∈ R at iteration j, the cost is defined as :

  scheme, the authors use large values of β for the algorithm to converge faster. However, as discussed before and also mentioned in[START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF]), large values of β may cause the problem to become ill-conditioned and create problems to assess the performance of deconvolution following the evolution of β. To verify the claim, we have conducted an experiment. In our test, we have used two images of Lena and the Cameraman. For this experiment we have used 3 different φ(u). Illustrated in Figure(3.3), it is evident that the quality of the reconstructed image is maintained once β has reached a certain value, i.e. β ∈ [2 7 , 2 8 ] for 3 different regularizers (i.e. ϕ 3/4

Fig. 6 :

 6 Fig.6: PSNR of images recovered in experiment 3.1 for different , Lena (above) and Cameraman (below). The performance does not improve once reaches a certain value.

Figure 3 . 3 :

 33 Figure 3.3: PSNR of images recovered in experiment 3.1 for different β, Lena (above) and Cameraman (below). The performance does not improve once β reaches a certain value.

Figure 3 . 4 :

 34 Figure 3.4: (Left) Evolution of the cost function for deblurring Cameraman, µ = 300, p = 2/3. Different values of fixed β. (Right) Comparison between fixed beta and dynamic beta.β i = κ * β (i-1) , i: iteration value. We note that larger β values converge slower than small ones. Fixed β can provide comparable performance to dynamic if µ is large.

  .23) The structure comparison function s measures the correlation coefficient between the two images x and x. σ xx is the covariance between the two images. The positive values of the SSIM index are within the range of 0 and 1 i.e SSIM ∈ [0, 1]. A value of 0 means no correlation between images, and 1 means that x = x.

Figure ( 3 . 5 )

 35 Figure(3.5) below shows the diagram flow of the structural similarity measurement system. Signal X and signal Y refer to the reference and sample images. Signal Y is the same as x of our formulation. For image quality assessment, it is usually recommended to apply the SSIM index locally rather than globally.

Fig. 3 .

 3 Fig. 3. Diagram of the structural similarity (SSIM) measurement system.

Figure 3 . 5 :

 35 Figure 3.5: Visual demonstration of the structural similarity (SSIM) measurement system (Wang et al. (2004)).

  (

FIGURE 3 . 6 -

 36 FIGURE 3.6 -Data set used for Experiments

Figure 3 . 6 :

 36 Figure 3.6: Data set used for Experiments.

  (3.7)(b) shows the same signals with blur and some added white Gaussian noise. The original signals show a different characteristic in the presence of noise and blur. For example, the piecewise constant signal behaves as a complex signal compared to its original version. Also, the signal for the Lena image has changed so much that it looks very different from its original counterpart.By looking at the blurred noises signals, it's apparent that the inverse problem of denoising/deconvolution is very ill-posed. This demonstration shows the importance of the prior information that has to be introduced to the optimization problem.From another point of view, the simple observation proves that regular classical regularizers are not sufficient for all image deconvolution/denoising applications. Consequently, the new potential function will provide new regularizers to better recover the signals following the prior information needed to recover its original property.

FIGURE 3

 3 FIGURE 3.7 -(a) Signal representation of original Images. (b) Representation of t original signals in the presence of blur and noise.

Figure 3 . 7 :

 37 Figure 3.7: Two demonstration of signals: (a) Signal representation of original Images. (b) Representation of the original signals in the presence of blur and noise.

3. 5 .

 5 EXPERIMENTS OF DECONVOLUTION PERFORMANCE a Lena. b piecewise constant.

Figure 4 . 1 :

 41 Figure 4.1: Graphical demonstration of a typical manual book scanner scheme

Figure 4 . 2 :

 42 Figure 4.2: CopiBook OS scanner of I2S

Figure 1 :

 1 Figure 1: SFR plots (top) and associated images demonstrating image sharpness and limiting resolution

Figure 4 . 3 :

 43 Figure 4.3: SFR plots of three images with different resolutions taken from (Williams and Burns (2008))

Figure 4 . 4 :

 44 Figure 4.4: The meaning of different SFR curves.

Figure 4 . 5 :

 45 Figure 4.5: Slant Edge Target also known as QA62.

Fast

  Microsoft Visual Studio IDE for our C++ software development. It is a very powerful IDE widely used for software development applications such as web services, mobile apps and software. It contains completion tools, compilers, and other features to facilitate the software development process. In our programming, We have benefited from the functions of the OpenCV library to develop our algorithm. OpenCV is an open-source computer vision library, Originally developed by Intel, publicly available for computer vision researchers and engineers . It is widely used for computer vision applications in Artificial Intelligence, Machine Learning algorithms, and any tasks that need image processing. Some of the main applications of OpenCV are as follow: Image processing Face detection / face recognition

Figure 4 . 6 :

 46 Figure 4.6: Screenshot of the Image watch user interface.

Figure 4 . 7 :

 47 Figure 4.7: Custom designed user interface of proposed deconvolution.

Figure 4 . 8 :

 48 Figure 4.8: Result of deconvolution on EAGLE camera sparse images.

Figure 4 . 9 :

 49 Figure 4.9: Result of deconvolution on EAGLE camera complex images.

Figure 4 .

 4 Figure 4.10: Slant edge target as input image

Figure 4 .

 4 Figure 4.11: MTF curve for EAGLE camera samples.

Figure 4 . 12 :

 412 Figure 4.12: Scoring table for MTF measurement for EAGLE camera

Figure 4 .

 4 Figure 4.13: Result of deconvolution on raw-super-resolved images with sparse property.

Figure 4 .

 4 Figure 4.14: Result of deconvolution on raw-super-resolved images with complex property.

Figure 4 .

 4 Figure 4.15: Slant edge target as input image for E-XTRA camera.

Figure 4 .

 4 Figure 4.16: MTF curve for E-XTRA camera samples for two mode of 800 DPI and 600 DPI.

Figure 4 .

 4 Figure 4.17: Scoring table for MTF measurement for E-XTRA camera

Figure 4 . 18 :

 418 Figure 4.18: Quality comparison of classical images of I2S and new super-resolved Image

Figure 4 .

 4 Figure 4.19: ATLIS Modulo microscope for cell analysis

FIGURE 4 .

 4 FIGURE 4.20 -Deconvolution/denoising result for two samples of microscopic imaging accquired with optimal ' p q pxq potential function. 103

Figure 4 .

 4 Figure 4.20: Deconvolution/denoising result for two samples of microscopic imaging accquired with optimal ϕ p q (x) potential function.

Figure 4 . 21 :

 421 Figure 4.21: TZcam terahertz camera of i2S

Figure 4 .

 4 Figure 4.22: (a) THz image of a 1951 USAF resolution test chart at 0.97 THz acquired with TZcam and x0.25 magnification lens. (b) Deconvolved image (c) Plots of the Line Pairs pixel in grey values of the group-1 element 2 of the chart of original THz image(d) deconvolved image.(e) Sample industrial application at 2.5THz. (f) Deconvolved image of (e)

Figure 4 .

 4 Figure 4.23: Commercial representation of BOOSTPIC

Figure 5 . 1 :

 51 Figure 5.1: PSNR values after deconvolution using the regularizer ϕ p q (x) as a function of p and q. The region around p=1.35 and q=-1 provides the best performance.

Figure 5 . 2 :

 52 Figure 5.2: Singularity exponents of the Musca observation Herschel at 250 µm. Top: Display of the singularity exponents of the raw, unfiltered map. Below: after noise reduction, p = 1.5, q = -1, λ = 0.1.

Figure 5 . 3 :

 53 Figure 5.3: Singularity exponents of the Musca observation Herschel at 250 µm, zoom on subregion, after noise reduction, p = 1.5, q = -1, λ = 0.1. Compare with (Yahia et al. (2021)), fig. 7.

Figure ( 5 . 3 )

 53 Figure(5.3) displays the map of the singularity exponents of a sub-region of the main filament. Comparing to[START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]), Fig.7, the extremely complex intrication of striations and sub-filaments at lower scale is now becoming even more evident. The values displayed in Figure(5.3) are the singularity exponents h(x) coming from formula (B.7) in[START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF] defined by a local correlation measure as explained in that reference.

Figure 5 . 4 :

 54 Figure 5.4: Visualization of the log of gradient's norms, Musca observation Herschel map, after noise reduction, p = 1.5, q = -1, λ = 0.1. Although the log of gradient's norms enhances the striations, the resulting spectrum, shown in Figure 5.5 is very different than the one obtained with the correlation measure.

Figure 5 . 5 :

 55 Figure 5.5:Left: comparison of obtained singularity spectra for Musca observation map: raw data, log of gradient's norms, correlation measure. The red and green curves are generated from the filtered observation map, the blue spectrum is generated from the raw map.The green spectrum has been computed from the values of h(x) given by Equation(5.6). The two other spectra (blue and red) make use of the values h(x) from formula (B.7) in[START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]. Right: comparison between the singularity spectrum of Musca unfiltered (blue) computed with h(x) given by Equation(5.6) and the singularity spectrum of Musca with noise reduction p = 1.5, q = -1, λ = 0.1 computed with h(x) given by Equation(5.6).

  new variable, u, to serve as the argument of φ(.) under the constraint that u = Dx. A quadratic function was added to penalize the difference: argmin problem formulated in the form of Equation (6.8) can be solved by an alternating Assume the additive Gaussian noise model for each image patch : (x; µ r , C r ) (6.12)In this formulation, ω 1 , ..., ω K are the mixing weights for GMM and N (.; µ, C)represents a Gaussian density of mean µ and covariance matrix C. We will assume that the components have zero-mean i.e µ r = 0 and the cross-covariance,C xn is also zero. The sum of two independent zero-mean Gaussians is a zero-mean Gaussian with covariance equal to the sum of the covariances. Under this condition and according to Equation (6.10) we have:p(y) = K r=1 ω r N x; µ r , C r + σ 2 I (6.13)Recall Bayes' theorem stated mathematically as the following equation: ' theorem with p(x) defined in Equation (6.12), we get :

  .22) of Lemma 6.1. In our case H = I, we have : xr = C r C r + σ 2 -1 y (6.23) Inserting GMM prior for y and p r (y) = N (y; 0, C r + σ 2 I) into Equation (6y; µ r , C r + σ 2 I xr K r=1 ω r N y; µ r , C r + σ 2 I (6.24)We use Equation (6.24) to obtain a clean x patch from noisy observations y. It can be written in packed form:y) = C r C r + σ 2 -1 y (6.26) ψ r (y) = ω r N y; µ r , C r + σ 2 I K r=1 ω r N y; µ r , C r + σ 2 I (6.27)In this formulation, ψ r (y) is the posterior probability that (noisy) patch y belongs to the r th component of the GMM and s r (y) is the MMSE estimate of that patch; if it was known that it had been generated by the r th component. It is evident that s r (y)and ψ r only depend on the noisy y. Although s r (y) is a linear function of y and the MMSE estimate x is a nonlinear function. After computing the MMSE patch estimates, x1 , ..., xN where N is the total number of patches, they are returned to their location and combined by averaging.

Test 1 :

 1 Figure (6.1) For this test, we used 16 images for training. These pictures had different orientations. For the first experiment, all training images were clean without noise and blur. So in this case H = I and n = 0. Test 1: We chose a piecewise constant image for the test image, and we added noise with σ = 0.03. The result of denoising is shown in Figure (6.2). We calculated the PSNR value of the noisy image, (a), and denoised image, (b), with the original image without noise. The PSNR value for the noisy image is 30.2 and for the denoised image is 36.4. The result shows that the denoiser performs reasonably well when the training

Figure 6 . 1 :Test 2 :

 612 Figure 6.1: Sample of training images for sparse data set

Figure 6 . 2 :

 62 Figure 6.2: Result of denoising using clean training images

Figure 6 . 3 :

 63 Figure 6.3: Result of denoising using mixture of noisy and clean training images

Figure 6 . 4 :Figure 6 . 5 :Figure 6 . 6 :

 646566 Figure 6.4: Result of denoising using mixture of noisy and clean training images -zoomed signal view

Figure 6 . 7 :Test 3 :

 673 Figure 6.7: Result of denoising using mixture of noisy and clean complex training images

Figure 6 . 9 :

 69 Figure 6.9: Signals of deconvolution using a mixture of noisy and clean complex training images

Table 1 . 1 :

 11 Categories of SR algorithms -Classification of SR algorithms (SISO, MISO, and MIMO) and some potential applications for each class.

multiple-input multiple-output (MIMO) super-resolution, commonly known as video-tovideo SR. A typical application of these techniques can be to enhance a video sequence captured by video cameras. A brief summary of the classification of the super-resolution techniques with potential applications is presented in Table

1

.1.

  The data fidelity term, ρ(y k -H k x), in the Equation (1.3) is used to constrain the residuals between the real LR images and the simulated ones obtained, and it is usually associated with the noise model. As an example, the 2 norm based linear least-squares

	term is widely used (Bouman and Sauer (1993); Molina et al. (2003); Shen et al. (2007);
	Song et al. (2010); Yuan et al. (2010b)), as p = 2 in (1.4). One of the key advantages
	of the 2 norm problem is that it is easy to solve and, in this regard, many efficient
	algorithms exist (Farsiu et al. (2004b); Shen et al. (2015b)). The result solved by the 2
	model would be optimal, as checked experimentally, when the model error is Gaussian
	distributed (Song et al. (2010)). Recalling the form of data fidelity in Equation (1.3), let
	us define a l p data fidelity term:

  images or training data for the Deep Neural Network (DNN) scheme. Because of the sensibility and practical nature of the scanners, any available database which is open source can not be used. Detailed training data needs to be gathered for scanners. From another point of view, the scanners are used in various applications at different locations.

	Therefore, we do not have exact information about the experimental data. Suppose we
	use precise training data for a special neural network. In that case, the algorithm will
	probably produce faulty results upon presenting the network with new data that is not
	compatible with the training data.
	Most DL-based algorithms are designed for precise usage and practical application;

these algorithms do not offer much flexibility. In other words, Convolutional Neural Networks (CNNs) have been proposed to learn a parametric model from a massive amount of data, i.e., hundreds of thousands of triplets of sharp images, blurry images and uniform blur kernel, that can be used as a learned prior

[START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF]

or to directly apply deblurring on a picture

[START_REF] Schuler | Learning to deblur[END_REF]

) but may fail if a test image

Table 3 .2: Accuracy

 3 table 3.2, we show the Mean Squared Error (MSE) values for various values of p. The MSEs considered here are computed: i) Between the approximated solution and the original image.

		Lena	Cameraman	Piecewise Constant
	p	Exact Approx Exact Approx Exact	Approx
	3/4	0.0024 0.0024 0.0035 0.0035 0.0014	0.0014
	1/2	0.0025 0.0025 0.0036 0.0036 0.0011	0.0011
	1/3	0.0027 0.0027 0.0041 0.0040 0.0016	0.0017

ii) Between the exact solution and the original image.

As the results show, the approximation provides a very close solution to the exact solution with less complexity and less time. The closeness of solutions supports the claim made previously: the use of approximated thresholding operators is almost indistinguishable from exact solutions. of approximation.The difference between approximated and exact solutions is negligible.

  deconvolution. In Proceedings of the Conference on Computer Vision and Pattern Recognition (CVPR), 2017. [28] Schuler CJ, Hirsch M, Harmeling S, Schölkopf B. Learning to deblur. IEEE transactions on pattern analysis and machine intelligence. 2015 Sep 23;38(7):1439-51. [29]Ongie G, Jalal A, Metzler CA, Baraniuk RG, Dimakis AG, Willett R. Deep learning techniques for inverse problems in imaging. IEEE Journal on Selected Areas in Information Theory. 2020 May 1;1(1):39-56. [43]A. Beck and M. Teboulle, "A fast iterative shrinkagethresholding algorithm for linear inverse problems," SIAM Journal on Imaging Sciences, vol. 2, pp. 183-202, 2009.[44]J.M. Bioucas-Dias and M.A.T. Figueiredo, "A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration," IEEE Transactions on ImageProcessing, vol. 16, pp. 2992-3004, 2007. Different values of fixed . (Right) Comparison between fixed beta and dynamic beta. i "  ˚ pi´1q , i: iteration value. We note that larger values converge slower than small ones. Fixed can provide comparable performance to dynamic if µ is large.

	2.6	¨10 4							2.4	¨10 4
	2.4 2.5						=1024 =512 =256 =125 =50		2.3						=1.3 =1.5 =2 =512
									2.2
	2.3
	2.2								2.1
	cost
	2.1								2
	2
									1.9
	1.9
									1.8
	1.8
	1.7	0	10	20	30	40	50	60	1.7	0	10	20	30	40	50	60
					iteration								iteration

Fig. 6: (Left) Evolution of the cost function for deblurring Cameraman, µ " 300, p " 2{3.

Table 3 . 4 :

 34 PSNR for different kernel sizes for different regularizers. The proposed regularizer ϕ p q (x) provides a modest PSNR enhancement. The values in bold correspond to the best performance.

https://www.i2s.fr/fr https://www.anrt.asso.fr/fr

http://opencv.org.

Acknowledgments

Comparison of PSNR and SSIM of reconstructed images using regularizers: l 1 ,l p and ϕ p q (x). The optimal value of µ is also reported. There exists one ϕ p q (x)) able to provide better performance in both scales.The values in bold correspond to the best performance.

to optimal PSNR and SSIM. Therefore, we have used grid search for each regularizer to find the optimal value of µ. It is important to note that the value of µ for each regularizer,ϕ p q (x), differs according to the quality measurement matric; this is because each quality matric considers different properties for scoring. The values of µ are reported accordingly in the table. The PSNR and SSIM values of the blurred image are also included. When compared to the PSNR and SSIM values of restored images, the algorithm's effectiveness is evident. In this experiment, we have used the kernel from (Zhang et al. (2012a)) and additive Gaussian noise with σ = 0.02. This kernel is particularly interesting since it is a very complex blur. This kernel destroys the information on the original image/signal extensively, therefore, making it harder for the algorithm to restore the data. This helps to evaluate ther robustness of the deconvolution. According to the results in table 2, One of the key results is that we see that it is possible to find a function among all ϕ p q (x) regularizers that can outperform regularizers such as l 1 [START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF]; [START_REF] Wang | A new alternating minimization algorithm for total variation image reconstruction[END_REF]) and l p norm [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF]; [START_REF] Zuo | A generalized iterated shrinkage algorithm for non-convex sparse coding[END_REF]). A grid search was conducted to find the optimal values of p and q for ϕ p q (x). This grid search will be explained shortly. Our experiments show that the variation in performance (PSNR) across values of p and q is similar for different images, convolution kernels and levels of noise. However, some minor variations are observed in each case. The scanner of I2S will be used in various unknown conditions; we have to test if the algorithm is resistant to changes. One of the main alterations occurs when the image's capturing conditions are varied; these factors include luminance, vibrations of the camera or scanner, and poor data. These alterations can be modelled by changing the properties of noise and the blur of the algorithm and examining how it performs.

To model such practical conditions and variations, we perform deconvolution with a uniform kernel of varying size; Table (3.4) reports the PSNR obtained by deblurring Lena. The result shows that as the conditions worsen, our proposed method and ϕ p q (x) function outperforms the other classical regularizers significantly. For example, in the case of the kernel size eleven by eleven, ϕ 1.1 3.2 (x) provides a considerable improvement. We see that for all kernel sizes, there is a ϕ p q (x) that offers better performance than previous regularizers. However, the optimal value of p and q may vary along with the size of the kernel. q (x) . For both images, an improvement can be assessed visually. The PSNR values for all images have been provided in Table 3.3

EAGLE Camera Samples

Before capturing the images, the EAGLE camera was calibrated according to the standards so that captured images would have the maximum quality possible. This is done to avoid any additional and unwanted parameters affecting the images, which may alter the exact performance evaluation of the deconvolution algorithm. Due to the confidentiality of the calibration system, the method of calibration cannot be explained in this chapter. Also, due to confidentiality, we could not disclose the parameters and the ϕ p q (x) used for the results. We will only present the visual results. Therefore, in the following parts, we will only show the screenshot of the green region of the Figure (4.7).

The result of deconvolution on various samples are shown on Figure (4.8) and (4.9).

The images with stars represent the results of deconvolution.

In Figure (4.8), we have shown the results of images with sparse properties. These images represent data such as book texts, newspaper texts, and simple figures. In region (1), the space between the lines are well distinguished, and the numbers are more visible.

In region (2), the text is more explicit and much easier to read than the original version.

The noise in both regions and other parts are reduced noticeably.

In Figure (4.9), we have shown the results of images with complex properties. In region (1), the patterns on the hair and face of the women are more visible. In region (2), the lines, structure of cloth and the curved lines have been well recovered. In region (3), small details of the hand such as dots and spots are better observed and visible.

In region (4), The handwritten text is clearer to read. The pattern and structure of the background paper have been recovered. In all retrieved results, the image contrast and sharpness have been significantly improved. The boundaries and edges between shapes and characters are more detectable. The results in both figures presented above show the effectiveness of the algorithm. The image quality has been enhanced for all types of data which proves the robustness of our proposed deconvolution algorithm.

The sizes of the images captured with CopiBook scanners with EAGLE camera are 10000 × 7096 pixels. These are huge images carrying a lot of information from the captured image. In Figures (4.8) and (4.9), we have shown zoomed regions of the large image because the effects of the image restoration can be better observed in Chapter 5: Results on Herschel observation maps This chapter will show the results of applying our deconvolution algorithm in Herschel observation maps. This is ongoing work, therefore. We will briefly review the results. The Herschel Space Observatory (aka Herschel) is an ESA facility available for the worldwide astronomical community . Herschel is opening a new opportunity to study how the universe has evolved into the universe we observe today. In this chapter, we will pursue to examine the effect of our new deconvolution algorithm on such image datasets and the improvements it may bring. 2021)). Indeed it is strongly believed that these filaments are involved in the process of star formation in the following manner: in a first step turbulence is thought to provide a complex intrication of filamentrary structures, while in a second step gravity takes over on the densest filaments [START_REF] André | Interstellar filaments and star formation[END_REF]). Moreover, it has been pointed out that the filaments all exhibit similar widths while only the densest ones contain prestellar cores (see [START_REF] André | Interstellar filaments and star formation[END_REF]; [START_REF] André | From filamentary networks to dense cores in molecular clouds: toward a new paradigm for star formation[END_REF]; [START_REF] Roy | How the power spectrum of dust continuum images may hide the presence of a characteristic filament width[END_REF] and Suri et al. ( 2019))

Introduction

Observations with

In [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]) it has been shown that to accurately access the multiscale properties of an Herschel observation map like the Musca SPIRE 250µm of the Gould Survey program, it is necessary to operate a reduction of the Gaussian noise present in an observation map; such noise comes from the cosmic infrared background (CIB) and https://www.cosmos.esa.int/web/herschel/home the cosmic microwave background, among other disturbances. This noise, of a Gaussian nature, tends to "log-normalize" the singularity spectrum of an observation map, thus concealing the statistics coming from the filamentary structures at small scales, and consequently making it almost useless to interpret a singularity spectrum in terms of underlying dynamics. In addition, the reduction of Gaussian noise makes it possible to visualize, via the geometric distribution of singularity exponents on an observation map, the extraordinary complexity of the spatial distribution of filamentary structures at small scales, i.e. structures which are responsible for the non-log-normal character of the obtained singularity spectra. In ( [START_REF] Yahia | Description of turbulent dynamics in the interstellar medium: multifractal-microcanonical analysis-i. application to herschel observations of the musca filament[END_REF]) the problem of noise reduction was solved by considering a sparse noise-reduction algorithm which tends to eliminate Gaussian noise while keeping the coherent gradient information at small scales. This was done by solving the L 1 -L 1 optimization problem. In other words, both F (x) and

In this chapter, we will slightly change the notation of the problem. However, the overall method of optimization and minimization are analogous to the Equations and method explained in chapter 3. In this chapter we will consider the following problem:

with: H matrix representation of a convolution operator H (PSF); if this convolution is periodic, H is then a (block) circulant matrix. In this chapter, the beam effect in the images is modeled in the matrix H. The goal is to recover s from the observational map y. In this case, we will consider the following optimization problem:

In the chapter, we will bring considerable improvements to the Gaussian noise reduction problem by considering other types of priors, which are more parameterizable and allow finer control over the preservation of filamentary structures. We will show on the one hand the spectacular improvements in visualization of these filamentary structures at different scales, thus exposing all the complexity of the geometric organization

Chapter 6: Plug-and-Play Priors

In this chapter, we will begin the second phase of our thesis. The first part involved the classical method of solving inverse problems such as deconvolution and denoising.

As mentioned before, other than classical methods, recently, much research has been conducted for learning methods such as CNNs. In the second part, we discuss a pipeline known as Plug-and-Play, which has recently gained much attention.

Introduction

The key objective of the first part of the thesis was to develop an efficient denoising algorithm and an effective super-resolution scheme to be implemented in cameras and products of I2S. Due to the limitations of available hardware, the main solution was to use the classical methods, which can be fast, less complex and easily used in practical understanding of plug-and-play and the benefits of this method. It will be shown that learning methods can be used and implemented within the technique and structure of a classical optimisation problem. Before we dig into plug-and-play, we will first provide a general overview of deep learning techniques for inverse problems and the advances deep learning methods have achieved in image processing. We will then present the idea of plug and play and consider using this technique in our proposed deconvolution algorithm.

Short review of Deep learning methods

Deep learning (LeCun et al. ( 2015)) approach has been known as a state-of-the-art tool in artificial intelligence with applications in face recognition [START_REF] Schroff | Facenet: A unified embedding for face recognition and clustering[END_REF]), machine translation [START_REF] Wu | Google's neural machine translation system: Bridging the gap between human and machine translation[END_REF]) and medical image analysis [START_REF] Gulshan | Development and validation of a deep learning algorithm for detection of diabetic retinopathy in retinal fundus photographs[END_REF]).

Researchers have reported improvements in applying CNNs

Microscopy:

In

Geophysical imaging:

Seismic inversion and imaging is the application of reconstructing the Earth's interior activities by modelling the physical reproduction of seismic waves. Comparing simulated synthetic measurements to actual acoustic recordings of reflected waves can be used to parametrize the models of ill-posed inverse problems. Some of these DL methods [START_REF] Hansen | Efficient monte carlo sampling of inverse problems using a neural network-based forward-applied to gpr crosshole traveltime inversion[END_REF]) include techniques that use generative models [START_REF] Mosser | Stochastic seismic waveform inversion using generative adversarial networks as a geological prior[END_REF]), which are constrained by partial differential equations [START_REF] Yang | Physics-informed generative adversarial networks for stochastic differential equations[END_REF]).

Other computational imaging applications:

Although minimization scheme, i.e. solving two sub-problems iteratively. It is easy to minimize the objective function in (6.8) with respect to either u or x as formulated in Equation (3.11). As discussed in the previous part, subproblem (S1) could be seen as the MAP solution of a denoising problem. Therefore we will replace this step with a GMM based denoiser: u will be calculated with the GMM denoiser. Among the many state-of-theart denoising methods, we opted for a GMM-based, patch-based denoiser proposed in [START_REF] Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF]). The reasons for this particular choice can be described in three factors:

Several authors have shown that a GMM is a good prior for clean image patches [START_REF] Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF]; [START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF]; [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]).

A GMM prior can be effectively learned from any external dataset of clean images [START_REF] Yu | Solving inverse problems with piecewise linear estimators: From gaussian mixture models to structured sparsity[END_REF]; Zoran and Weiss ( 2011)), or directly from noisy patches [START_REF] Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF]).

By learning the GMM prior from an external dataset of clean images, we can develop class-specific priors/denoisers [START_REF] Luo | Adaptive image denoising by targeted databases[END_REF]).

Sub-problem (S2) will be calculated as follow:

  (6.9)

Denoising and using Gaussian mixtures

This part will briefly describe the function of the GMM based denoiser that we have used in our plug-and-play implementation. In summary, we start by extracting all the patches from the training images. Next, we use the expectation-maximization algorithm to determine the parameters of the GMM, i.e. the GMM is learned using the expectation-maximization (EM) algorithm. Then, this GMM model is used as a prior to denoise the observed patches. Under this assumption, it can be done by computing the minimum mean squared error (MMSE) estimate in closed form [START_REF] Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF]).

In this subsection, unlike the other parts of the thesis, x and y they refer to an arbitrary individual patch. 6.9). The plug-and-play approach is able to recover the signal as close as possible to the original signal. Overshoots in the recovered signal have been reduced. The recovered signal shows that the plug-and-play approach performs well for denoising, but the ability to perform deconvolution is not as profound as denoising.

We have to know that the training process is time-consuming despite not using many training images. This is one of the main drawbacks of using learning methods. Also, as seen in the results, to make the algorithm more flexible, we need to extend the range of training images used, which will add to the time for training. The number of patches and GMM components also play a role in the processing time of the algorithm. The higher the number of these parameters, the more time it will take for the algorithm to perform deconvolution. Since we are using Fourier transform in our (S2) sub-problem, it was observed experimentally that our proposed algorithm is faster than other plug-and-play schemes based on ADMM. Due to confidentiality, the exact parameters and the timer processing cannot be disclosed. The plug-and-play scheme is still under experimental evaluation at I2S.

Notes on the performance:

The performance of the deconvolution depends on several items and parameters.

As we saw, the plug-and-play scheme works quite well for denoising because the training

Chapter 7: Conclusion and prospects

In this chapter, we conclude the thesis by summarizing our main contributions and results. We also highlight the ongoing works we are conducting as an extension.

Summary of the contributions

This thesis has presented the framework of developing a super-resolution scheme and image deconvolution algorithm for various cameras and products. The main objective of this thesis was to develop a practical algorithm that can be fast, accurate, and flexible. We presented the framework of developing a super-resolution scheme with the help of sensor displacement of the camera to achieve high-resolution images with up to 2 times the resolution of the original ones. To enhance the quality of the images, we proposed an image deconvolution algorithm with an alternating minimization scheme to be used aside and the super-resolution scheme. We considered classical optimization methods for our image deconvolution algorithm for the first stage of our implementation.

One of the goals of our thesis was to develop a fast and low complex algorithm by preserving the accuracy for practical application in the industry.

We can summarize our work in the following items:

We proposed the super-resolution scheme to capture low-resolution images, each image having half-pixel displacement. We assembled the low-resolution images to capture a higher resolution image, which will allow us to capture images of 800 DPI.

We proposed an image deconvolution/denoising algorithm based on variable splitting and proximal operators. Our proposed method proved to be fast, accurate and less complex, which made it ideal to implement in practical applications for products like connected image acquisition systems and embedded systems with limited computation capability.

For the first time, we propose the use of new potential functions to be used as prior information in imaging inverse problems. The scientific experiments proved that we could offer new functions that outperform the classical regularizers used in previous studies. Potential functions performed well for different kinds of images, such as sparse and complex. We used various image quality metrics to evaluate our scientific experiments, and all results showed an enhanced improvement in image quality.

We have successfully transitioned the algorithm to practical implementation. For the implementation of deconvolution, we have used C++ programming language.

We developed new libraries and functions which did not exist before for image deconvolution/denoising algorithms. We have successfully implemented the algorithm on two cameras of I2S used for book scanners.

The image deconvolution/denoising algorithm proved effective for other products and cameras of I2S. The algorithms were implemented for microscopic and Terahertz cameras, and the results showed improved quality and resolution in images captured with such cameras with parameter adaptation.

In an ongoing work, we also adopted the deconvolution/ denoising algorithm on the Herschel observation maps. With the new potential function, we have improved the denoising and deconvolution of observational maps.

For the second step of the thesis, we proposed plug-and-play priors for image inverse problems. This method provides the ability to use a learned denoiser within the structure of a classical optimization scheme. We suggested using plug-and-play priors in the alternating minimization scheme for the first time.

The denoising results were promising. However, the algorithm does not perform deconvolution as good as the classical method.

Future work and prospects

One of the main challenges we faced in our work consisted of modelling the PSF or H. The captured image goes through various steps and systems during the acquisition of images in image capturing devices such as cameras, microscopes, and scanners.

These systems include the interpolation, the lens, the A/D chip, and the digital image processing unit. Each system leaves its effect on the resulting image. These effects are 123 modelled with the PSF matrix. It is always difficult to define a general and global PSF for each image deconvolution application, especially in industry. To better model the PSF, it would be interesting to design a PSF estimation to be implemented before the deconvolution, which is specifically designed for the application. A lot of blind image deconvolution algorithms have been studied before. However, not many are considered for industrial use. In the user interface, the user can define the PSF before running the program. Would it be possible to improve the PSF selection using non-blind image deconvolution techniques?

We would like to use Expectation Propagation for further development of our proposed image deconvolution/denoising algorithm. Using this technique, we pursue to convert our proposed algorithm into a non-blind image deconvolution scheme.

The plug-and-play approach has been recently developed, and still, more evaluations are needed to understand this scheme. Our proposed plug-and-play can be further experimented with to improve the quality of image deconvolution. Various quantities of training data can be studied to further understand its effect on the resulting image.

BIBLIOGRAPHY Shen, H., Peng, L., Yue, L., Yuan, Q., and Zhang, L. (2015d). Adaptive norm selection for regularized image restoration and super-resolution. IEEE transactions on cybernetics, 46(6):1388-1399.

Shen, H., Zhang, L., Huang, B., and Li, P. (2007). A map approach for joint motion estimation, segmentation, and super resolution. IEEE Transactions on Image processing, 16(2):479-490.

Sinha, A., Lee, J., Li, S., and Barbastathis, G. (2017). Lensless computational imaging through deep learning. Optica, 4(9):1117-1125.

Song, H., Zhang, L., Wang, P., Zhang, K., and Li, X. ( 2010). An adaptive l 1-l 2 hybrid error model Suri, S., Sánchez-Monge, Á., Schilke, P., Clarke, S. D., Smith, R. J., Ossenkopf-Okada, V., Klessen, R., Padoan, P., Goldsmith, P., Arce, H. G., et al. (2019). The carma-nro orion survey-filamentary structure as seen in c18o emission. Astronomy & Astrophysics, 623:A142.