Image Restoration

Image restoration estimates a clean image from an observation that has been degraded in a known or partly known way: by noise, blur, missing pixels, a lower resolution, compression, haze, rain or snow. What separates it from image enhancement is the model. In the frame that image processing takes from Gonzalez and Woods, enhancement is judged subjectively, while restoration is objective because it inverts a mathematical or probabilistic model of the degradation. In our reading, every restoration method therefore has two parts, explicit or hidden in a network’s training data: a model of how the image was degraded, and a prior that chooses among the many clean images consistent with the observation.

The single-degradation tasks have their own articles: image denoising, image deblurring, super-resolution, image inpainting, image dehazing and, for under-exposed photographs, low-light image enhancement. This article treats what they share: degradation models as a family, the estimation problem they pose, the tasks without their own article (removing rain, snow and raindrops, and removing several kinds of degradation with one network), the general-purpose restoration networks that are evaluated on several tasks, and how results are compared. It builds on image noise, which supplies the noise models.

Problem Definition

A degradation model maps a clean image $x$ to an observation $y$. The family that the linked task articles use is linear with additive noise,

$$y = A x + n ,$$

with the image written as a vector, $A$ a linear operator and $n$ noise. The hub article writes the shift-invariant case as $g = h * f + \eta$; the family is wider than convolution:

Degradation	Operator $A$	Task article
Noise only	identity	image denoising
Blur	convolution with a point spread function	image deblurring
Missing pixels	a diagonal 0/1 mask	image inpainting
Lower resolution	blur, then subsampling	super-resolution
Compositions	products of the above	for example blur and decimation in super-resolution

Noise need not be additive and Gaussian: photon-limited images follow the Poisson–Gaussian model in image noise, and compression adds quantization error in a transform domain, which image compression describes. Real-ESRGAN’s second-order degradation, which ESRGAN covers, applies blur, resizing, noise and JPEG compression twice in sequence, a composition whose JPEG stages are not linear.

Weather degradations are not a linear operator applied to $x$ but layer models in which the clean image is mixed with other components. The haze imaging model, $\mathbf{I} = \mathbf{J}t + \mathbf{A}(1 - t)$, belongs to image dehazing. The rain model that Yang et al. (2017, Eq. 1) call widely used is the sum of a background layer $B$ and a rain-streak layer $S$, $O = B + S$, which makes deraining a separation of two signals. Yang, Tan, Feng, Liu, Guo and Yan added a binary map $R$ of where individual streaks are visible, $O = B + SR$ (Eq. 2), and then a model for heavy rain,

$$O = \alpha \Big( B + \sum_{t=1}^{s} S_t R \Big) + (1 - \alpha) A ,$$

with $s$ streak layers of different directions and the haze model’s transmission $\alpha$ and atmospheric light $A$ standing for rain accumulation, the veil that distant streaks form when they can no longer be seen individually (Yang et al., 2017, Eq. 3). Li, Tan and Cheong (2020, Eqs. 1–3) collect three such models from earlier work: heavy rain with streaks and veiling, $I = t(J + \sum_i R_i) + (1 - t)A$, after Li, Cheong and Tan (2019, not read); adherent raindrops on a window or lens, $I = (1 - M)J + K$, with a binary mask $M$ and the blurred imagery $K$ that the drops form, after Qian et al. (2018, not read); and falling snow, $I = zS + J(1 - z)$, with snowflakes $S$ and a binary snow mask $z$, after Liu et al.’s DesnowNet (not read; Li et al. cite both its 2017 preprint and its 2018 journal version). From the differences between these models they conclude that the problems do not share the same intrinsic characteristics, which they give as the reason why a method designed for one weather does not work on another.

Estimation. For the linear family with Gaussian noise of variance $\sigma^2$, the maximum a posteriori estimate maximizes $p(x \mid y) \propto p(y \mid x)\, p(x)$; taking the negative logarithm, it follows that

$$\hat x = \arg\min_x \; \frac{1}{2\sigma^2}\, \lVert A x - y \rVert^2 + R(x), \qquad R(x) = -\log p(x) ,$$

a data term set by the degradation model plus a regularizer set by the prior. A quadratic regularizer $R(x) = \tfrac{\lambda}{2\sigma^2} \lVert Dx \rVert^2$ on a difference operator $D$ makes the estimate the solution of the linear system $(A^{\mathsf T}A + \lambda D^{\mathsf T}D)\,\hat x = A^{\mathsf T} y$; for a circular blur it is the regularized inverse filter derived in image deblurring, and least squares treats the same penalty as generalized Tikhonov regularization. Total variation, sparsity in a wavelet basis and patch models give non-quadratic priors (total variation denoising, BM3D). Learned methods keep the structure and move the prior into data: Yang et al. (2017, Eq. 4) write their deraining problem as MAP estimation of $B$, $S$ and $R$ with priors on each, and note that a deep network learns those priors from training data and embeds them implicitly. Half-quadratic splitting with a learned denoiser in place of the prior step, which image deblurring describes for Zhang et al. (2017), keeps the operator $A$ explicit and learns only the prior.

Ill-posedness and the null space. When $A$ has a non-trivial null space, as the mask and subsampling do, every image $x + v$ with $Av = 0$ explains the data exactly as well as $x$, and only the prior decides the component in the null space; linear algebra defines the null space, and singular value decomposition works through the null space of a small averaging operator. A blur whose transfer function has no zeros has a trivial null space but is ill-conditioned, so the data constrain its weakly transmitted frequencies only up to amplified noise. The example below measures how the error of one estimator splits between the null space and its complement.

Input and Output

• Input: one degraded image, sometimes with the parameters of the degradation (a noise level, a blur kernel, a mask, a scale factor). Multi-frame inputs, such as bursts (burst photography) and video, give additional observations of the same scene.
• Output: an estimate of the clean image at the size the model defines, the same size for denoising, deblurring, inpainting and weather removal, larger for super-resolution. Model-based methods may also return the estimated degradation (a kernel, a transmission map, a rain-streak layer or mask).

Variants

• Non-blind and blind. Non-blind methods are given $A$ and the noise level; blind methods estimate them or work without them. Blind deblurring is the classic case (image deblurring).
• One degradation or several. AirNet’s authors observe that almost all existing approaches could handle only one specific degradation (Li et al., 2022, §1). Methods for several degradations differ in what they assume at test time, and the name “all-in-one” covers three different settings:
  • AOD-Net, the “All-in-One Dehazing Network” (Li, Peng, Wang, Xu and Feng, 2017), handles one degradation, haze; “all-in-one” refers to estimating the transmission and the atmospheric light together in one variable, as image dehazing explains.
  • Li, Tan and Cheong (2020) remove several kinds of bad weather (rain with veiling, snow and adherent raindrops) with one network and one set of weights, using one encoder per weather type.
  • AirNet (Li, Liu, Hu, Wu, Lv and Peng, 2022) restores images with noise, rain or haze without being told which corruption, or at which level. Its authors state in a footnote that “unknown” here means unspecified, not unseen, that is, each test corruption is of a type present in training, and each image carries one degradation while the dataset contains several.
• Synthetic and real degradations. Supervised methods are trained on pairs made by applying a degradation model to clean images, and tested either on more synthetic pairs, where a reference exists, or on real degraded photographs, where none exists.

Challenges

• Model mismatch. A method is only as good as its degradation model where it is used. Fu et al. (2017) found that heavy rain makes images hazy, an effect absent from their synthetic training data, and applied a dehazing method before deraining; super-resolution describes how bicubic-trained networks fail on real images for the same reason.
• Separating degradation from content. Rain streaks, snow and noise look like image structure. Fu et al. (2017) note that existing methods have difficulty removing rain while preserving structure when an object’s orientation is similar to that of the streaks; Yang et al. (2017) list over-smoothing of non-rain regions as a limitation of earlier methods, because streaks and background texture overlap.
• Unknown and mixed degradations. Real images can carry several degradations at once and their strength varies across the image; AirNet’s authors give self-driving cars that meet rain and haze in sequence or together as their motivating case (Li et al., 2022, §1).
• No ground truth for real images. Large numbers of real degraded and clean pairs are hard to obtain (Fu et al., 2017, §2.3), so the quantitative results in the papers discussed here are on synthetic degradations, or on real noise with estimated references.
• Hallucination. Learned priors can invent plausible content in the null space of $A$, which matters where the image is evidence; image processing discusses the trade-off with classical methods.

Approaches

The classical approaches are model-based: filters and regularized inverses for denoising and deblurring, interpolation and back-projection for super-resolution, diffusion and patch copying for inpainting, and priors such as the dark channel for haze. Their task articles cover them. The learned methods below are networks trained on synthetic pairs, and the papers read here form a sequence: one network per degradation (DDN and JORDER, 2017), one network for several weather types (2020, 2022), general architectures trained separately for each task (MPRNet, 2021; Restormer, 2022), and one model for several unspecified corruptions (AirNet, 2022).

Deraining

Deep detail network (DDN). Fu, Huang, Zeng, Huang, Ding and Paisley (CVPR 2017) found that a convolutional network trained to map a rainy image $X$ directly to the clean image $Y$ underfit and gave a strong colour shift, which they attribute to the mapping covering the whole range of pixel values. They made two changes. The network predicts the negative residual $Y - X$, whose values are concentrated in a much narrower range, and adds it to the input. And it sees only the detail layer $X_{\text{detail}} = X - X_{\text{base}}$, where the base layer is a low-pass version of $X$ obtained with a guided filter of radius 15; removing the background leaves mostly rain streaks and object edges, and most detail-layer values are close to zero. The objective is

$$\mathcal{L} = \sum_{i=1}^{N} \big\lVert f(X_{i,\text{detail}}, W, b) + X_i - Y_i \big\rVert_F^2$$

with $f$ a ResNet-style network of 26 layers with 16 filters of size 3 × 3 (Fu et al., 2017, Eq. 4, §3.1). They synthesized 14,000 rainy/clean pairs from 1,000 clean images, each rendered with 14 streak orientations and magnitudes, trained on 9,100 and tested on the remaining 4,900; on those, they report an SSIM of 0.90 ± 0.05 against 0.87 ± 0.07 for a Gaussian-mixture layer-prior method and 0.83 ± 0.09 for discriminative sparse coding (their Table 1). Applied to images with noise and JPEG artifacts, the same structure removed those too (their §3.8), the authors’ argument being that it suits corruptions that are high-frequency.

JORDER. Yang et al. (CVPR 2017) built a multi-task network on their rain model. A contextualized dilated network extracts features through three paths of two 3 × 3 dilated convolutions each, with receptive fields of 5 × 5, 9 × 9 and 13 × 13; from these features the network predicts the streak map $R$, then the streak layer $S$, then the background $B$, each from the previous predictions, trained with losses on all three. The recurrent version, JORDER-R, applies the network repeatedly and removes rain part by part, each pass working on the previous pass’s output. For rain accumulation they trained a separate network on synthetic veils and found the order streak removal, accumulation removal, streak removal better than the other orders they tried, because removing the veil makes faint streaks that the first pass missed more visible (§5.2, Fig. 8). They also introduced the Rain100L (one streak type) and Rain100H (five streak directions) test sets, synthesized from BSD200 images (§6). On Rain100H, measured on the luminance channel, they report a PSNR of 23.45 dB and SSIM 0.7490 for JORDER-R against 22.15 dB for JORDER and 15.66 dB for discriminative sparse coding (their Table 2).

Desnowing and adverse weather

Snow appears in single-weather work such as DesnowNet (Liu et al., 2018, not read), which Li, Tan and Cheong (2020) cite for the snow model and the Snow100K dataset, and in methods that remove several kinds of weather at once.

All in One Bad Weather Removal. Li, Tan and Cheong (CVPR 2020) handle rain with veiling, snow and adherent raindrops in one network with one set of weights. The generator has one encoder per weather type and a shared decoder. Between them, a neural architecture search chooses how to combine the encoders’ features from a set of operations that includes, besides separable and dilated convolutions, four designed after the weather models: a deveiling operation based on the haze model, a decomposition into low and high frequencies by guided filters, a residue operation, and self-attention in place of the raindrop attention maps that earlier methods learned from masks. A discriminator in the manner of a generative adversarial network judges whether an output is real and also classifies its degradation type, and the classification loss for one type updates only that type’s encoder (§3). Training samples 9,000 images per epoch from each of Outdoor-Rain, Snow100K and an augmented Raindrop set (§4.1). The authors report PSNR/SSIM of 24.71/0.898 on the rain-and-fog Test1 set against 21.56/0.855 for the heavy-rain method HRGAN, 28.33/0.8820 on Snow100K-L against 27.17/0.8983 for DesnowNet, and 31.12/0.9268 on Raindrop, where their PSNR is below that of AttentGAN (31.57) and of Quan et al. (31.44) (their Tables 1–3). The text says correctly that the method does not have the best raindrop PSNR, but Table 2 marks Quan et al.’s 31.44 as the best PSNR although AttentGAN’s 31.57 is higher, and the text cites Fig. 9 for the snow comparison, which is in Fig. 8.

TransWeather. Valanarasu, Yasarla and Patel (CVPR 2022) replaced the multiple encoders with a single transformer encoder and decoder. The encoder adds intra-patch transformer blocks that attend among sub-patches of each patch, on the argument that rain streaks, raindrops and snowflakes are small and several can fall inside one patch. The decoder uses learnable weather-type queries, in the manner of the object queries of the detection transformer DETR, which attend to the encoder’s features to produce a task feature vector. Training uses a smooth $L_1$ loss and a perceptual loss on VGG16 features, on the same training distribution as Li et al., which they call “All-Weather”: 9,000 Snow100K images, 1,069 Raindrop images and 9,000 Outdoor-Rain images (§4.1). With one model for all three tests they report 31.05/0.9509 on Test1, 33.78/0.9287 on Snow100K-L and 34.55/0.9502 on RainDrop, against 24.71, 28.33 and 31.12 dB for Li et al.’s network in the same tables, and against task-specific baselines that were trained separately on each dataset, including MPRNet (21.90 dB on Test1, after training on Outdoor-Rain) and SwinIR (their Tables 1–3). The paper’s queries are 512 learned embeddings, and visualized attention maps of some queries respond more to fog, rain or snow (§5).

The paper has several inconsistencies. Its Eq. 3 writes the snow model as $I = M \odot S + M \odot (1 - z)$, which differs from Li et al.’s $I = zS + J(1 - z)$ and leaves out the clean image. Its Table 1 lists “RESCAN + Dehaze (DHF)” twice, where Li et al.’s Table 1 has one de-hazing-first and one de-raining-first row with the same numbers. Its Table 3 gives AttentGAN 30.55 dB on RainDrop, where Li et al.’s Table 2 gives 31.57 dB. And it cites DDMSNet once as [73], its entry in the reference list, and once as [72], which is a dehazing network by Zhang et al.

All-in-one restoration for unknown corruption

AirNet (Li, Liu, Hu, Wu, Lv and Peng, CVPR 2022) treats noise, rain and haze as one problem without telling the network which one it faces. A contrastive-based degradation encoder learns a representation of the degradation by contrastive learning: two crops of the same degraded image are a positive pair, because they share a degradation, and patches from other images are negatives. The representation is taken from the encoder’s first layer, a spatial tensor rather than a vector. A degradation-guided restoration network of five groups of five blocks, each with two modules, uses it to modulate deformable convolutions, whose sampling offsets it helps predict, and spatial feature transforms, which scale and shift features. The loss is an $L_1$ reconstruction loss plus the contrastive loss (§3).

The paper compares two settings with the same baselines retrained: one model per degradation (“one-by-one”) and one model for all three (“all-in-one”), over BSD68 at $\sigma = 15, 25, 50$, Rain100L (200 training and 100 test pairs) and SOTS from RESIDE (§4.1). In the all-in-one setting AirNet has the highest average, 31.20 dB PSNR and 0.9095 SSIM, against 30.17/0.8990 for MPRNet retrained the same way (Table 4). In the one-by-one setting, MPRNet has the higher average (31.98 against 30.38 dB), with 38.26 against 34.90 dB on Rain100L and 28.21 against 23.18 dB on SOTS, while AirNet is ahead on BSD68 at all three noise levels; the text summarizes Table 4 as AirNet being superior to all the baselines in most cases. The authors leave open how it would handle blur and snow (§6). The repository README describes checkpoints for each single task and for the all-in-one setting; the repository has no licence file (as of November 2026), so its code was not read.

General restoration backbones

Between single-task networks and single models for many degradations lie architectures designed for restoration in general and trained separately for each task.

MPRNet. Zamir, Arora, Khan, Hayat, Khan, Yang and Shao (CVPR 2021) restore in three stages. The first two are encoder–decoders based on U-Net with channel attention blocks, which capture context; the last stage, an original-resolution subnetwork without downsampling, keeps fine detail. The input is split into four non-overlapping patches for the first stage, two for the second and the whole image for the third. A supervised attention module between stages produces an intermediate restored image, supervised by the ground truth, and computes from it attention masks that reweight the features passed on; cross-stage feature fusion carries multi-scale features from earlier stages to later ones. Each stage predicts a residual added to the degraded input, and the loss sums over the three stages a Charbonnier loss and an edge loss on the Laplacian of the images, weighted 0.05 (§3). Separate models are trained for deraining, deblurring and denoising, 40, 96 and 80 channels wide respectively (§4.2). For deraining they trained on 13,712 pairs from several datasets and tested on five sets, Test100, Rain100H, Rain100L, Test2800 and Test1200, computing PSNR and SSIM on the Y channel of YCbCr; MPRNet’s average was 32.73 dB against 30.75 dB for the previous best method in their table, MSPFN (Table 2).

Rajaei, Rajaei and Damavandi (2023) analysed the released MPRNet models for IPOL, in an MLBriefs article whose source code IPOL notes has not been reviewed. On four GoPro test images MPRNet’s deblurring PSNR was higher than SRN’s and MIMO-UNet’s (§3.1, Fig. 4). For denoising they added Gaussian noise to images unlike the SIDD training data: NAFNet scored higher on all four images at $\sigma = 10$, MPRNet on all four at $\sigma = 40$ (Figs. 5–7). Their conclusion that MPRNet’s performance “seems inferior” to the other methods used does not match their own deblurring figures. They also write that MPRNet was trained on Rain100L, while MPRNet’s Table 1 lists 0 training samples from Rain100L. The MPRNet repository is under an Academic Public License for non-commercial use (as of November 2026); only the licence was read.

Restormer. Zamir, Arora, Khan, Hayat, Khan and Yang (CVPR 2022) designed a transformer for full-resolution restoration. Standard self-attention costs grow quadratically with the number of pixels, so earlier restoration transformers such as SwinIR (described in super-resolution) and Uformer attend within 8 × 8 windows. Restormer’s multi-Dconv head transposed attention instead computes attention across feature channels: queries, keys and values are formed by 1 × 1 point-wise and 3 × 3 depth-wise convolutions, and their product is a $C \times C$ attention map, so the cost is linear in the number of pixels. Its feed-forward network is gated, the element-wise product of two branches one of which passes through a GELU. The network is a four-level encoder–decoder with layer normalization, it predicts a residual added to the input, and it is trained with an $L_1$ loss and progressive learning, from 128 × 128 patches in large batches to 384 × 384 patches in small ones (§3, §4). On the same five deraining test sets, with rows for earlier methods identical to MPRNet’s Table 2, Restormer reports an average of 33.96 dB, 1.05 dB above the best earlier entry, SPAIR (32.91 dB), and above MPRNet’s 32.73 dB (Restormer, Table 1). It also reports results for motion deblurring, defocus deblurring and denoising, each from a separately trained model. The repository is under the MIT licence (as of November 2026).

In our reading, the line from DDN to Restormer moves the prior from the model into the data while keeping a trace of the model: DDN and JORDER build the rain model into the network’s outputs, MPRNet and Restormer predict a residual added to the input, and the all-in-one networks spend part of their capacity on recognizing which degradation is present.

Datasets and Benchmarks

The benchmarks used by the papers above are mostly synthetic, made by applying a degradation model to clean images:

• Noise: BSD68 and Urban100 with Gaussian noise at fixed levels (AirNet tests on both), and the real-noise sets DND and SIDD; image denoising describes BSD68, DND and SIDD.
• Blur: GoPro, made by averaging frames of high-speed video, and HIDE and RealBlur; image deblurring covers the first.
• Rain: Rain100L and Rain100H (Yang et al., 2017), Test100, Test1200 and Test2800, collected for evaluation in MPRNet’s Table 1. Test2800 is the test part of Rain14000, which that table attributes to Fu et al. (2017) and splits into 11,200 training and 2,800 test pairs; Fu et al. themselves trained on 9,100 of their 14,000 pairs and tested on 4,900.
• Weather mixtures: Outdoor-Rain and its Test1 set (rain with fog), Snow100K with its Snow100K-L test set, and the real Raindrop pairs, combined as “All-Weather” by TransWeather.
• Haze: RESIDE, with its Outdoor Training Set and SOTS test set, described in image dehazing.

Numbers are comparable only within one protocol. JORDER-R’s 23.45 dB on Rain100H (JORDER’s Table 2) and MPRNet’s 30.41 dB on Rain100H (MPRNet’s Table 2) come from different training data, and the table rows they sit in were produced under different protocols. AirNet’s Rain100L numbers come from training on Rain100L’s own 200 pairs (or all three tasks), not on MPRNet’s 13,712.

Evaluation Metrics

With a reference image, results are reported as PSNR and SSIM; deraining papers compute both on the luminance or Y channel (Yang et al., 2017, §6; MPRNet, §4.3; Restormer, §4.1). MPRNet also converts the PSNR and SSIM differences to relative reductions in RMSE and in DSSIM, $(1 - \text{SSIM})/2$, against the best method. Without a reference, as for real rainy photographs, papers show images: Fu et al. (2017) state that no standard, high-quality reference-free measure existed and give only qualitative results for real data. Image quality assessment covers no-reference measures, and LPIPS a learned perceptual distance.

Example: One Estimator, Four Degradations

The example takes the quadratic MAP estimate

$$\hat x = \arg\min_x \lVert A x - y \rVert^2 + \lambda \lVert D x \rVert^2 ,$$

with $D$ the periodic forward differences, so that $D^{\mathsf T} D$ is the negative of the five-point Laplacian, and solves $(A^{\mathsf T}A + \lambda D^{\mathsf T}D)\hat x = A^{\mathsf T}y$ by conjugate gradients. It feeds the same solver either the right operator $A$ or the identity, the model of a method that only denoises, for four degradations of a 256 × 256 grayscale crop of scikit-image’s bundled chelsea image:

• noise: Gaussian noise of standard deviation 0.02 (on a 0–1 scale), $A = I$;
• blur: a circular Gaussian blur with $\sigma_b = 2$ pixels (mode="wrap"), then noise;
• mask: half the pixels observed at random, noise on the observed ones, the others set to 0;
• down2: the same blur, then every second pixel in each direction, then noise.

All operators are circular, and the data are generated with exactly the operators the right model uses, so the right model is exact; boundaries, which matter for real photographs, are treated in image deblurring. For down2 the identity model needs the observation on the fine grid, and the code prints both pixel replication and linear interpolation between the samples. The regularization weight is printed at $\lambda = 0.003$ and $0.03$.

Two measures are printed for five seeds (the noise and the mask change with the seed): PSNR with peak 1, and for the mask and down2, the share of the right model’s squared error $e = \hat x - x$ that lies in the null space of $A$. For the mask the null-space component is the error on the unobserved pixels. For down2 it is $e - A^{\mathsf T}(AA^{\mathsf T})^{-1}Ae$; $AA^{\mathsf T}$ is circulant on the coarse grid, with eigenvalues the average of $|\hat h|^2$ over the four fine-grid frequencies that alias to each coarse one, so the projection takes two FFTs. The code also checks its components: the conjugate-gradient solution for the blur against the closed-form Fourier solution of the same energy; that on unobserved mask pixels, where the normal equations reduce to $\lambda (D^{\mathsf T}D\hat x)_i = 0$, the solution is discrete-harmonic; that the down2 projector is idempotent; the adjoint of the down2 operator; and that the solver converged in every run.

import os
os.environ["OMP_NUM_THREADS"] = "1"; os.environ["OPENBLAS_NUM_THREADS"] = "1"; os.environ["MKL_NUM_THREADS"] = "1"
import numpy as np
from scipy import ndimage
from scipy.sparse.linalg import LinearOperator, cg
from skimage import data, color

SEEDS = [0, 1, 2, 3, 4]
LAMS = [0.003, 0.03]              # fixed at planning, not tuned
SIG_B, SIG_N = 2.0, 0.02          # Gaussian blur (circular), noise standard deviation

x = color.rgb2gray(data.chelsea())[22:278, 97:353]   # 256 x 256 crop, values in [0, 1]
N = x.shape[0]

def blur(z): return ndimage.gaussian_filter(z, SIG_B, mode="wrap")
def lap(z):                       # D^T D for periodic forward differences
    return 4 * z - np.roll(z, 1, 0) - np.roll(z, -1, 0) - np.roll(z, 1, 1) - np.roll(z, -1, 1)
def down(z): return blur(z)[::2, ::2]
def down_T(y):
    z = np.zeros((N, N)); z[::2, ::2] = y; return blur(z)
def up_rep(y): return np.kron(y, np.ones((2, 2)))          # pixel replication
def up_lin(y):                                              # linear interpolation on the sample grid
    for ax in (0, 1):
        z = np.repeat(y, 2, axis=ax)
        sl = [slice(None)] * 2; sl[ax] = slice(1, None, 2)
        z[tuple(sl)] = 0.5 * (y + np.roll(y, -1, axis=ax))
        y = z
    return y
ident = lambda z: z

infos = []
def solve(A, AT, y, lam):
    """Quadratic MAP estimate: argmin |A x - y|^2 + lam |D x|^2, by conjugate gradients."""
    op = LinearOperator((N * N, N * N), dtype=float,
                        matvec=lambda v: (AT(A(v.reshape(N, N))) + lam * lap(v.reshape(N, N))).ravel())
    s, info = cg(op, AT(y).ravel(), rtol=1e-10, maxiter=5000)
    infos.append(info)
    return s.reshape(N, N)

psnr = lambda xh: 10 * np.log10(1.0 / np.mean((xh - x) ** 2))

delta = np.zeros((N, N)); delta[0, 0] = 1
H = np.fft.fft2(blur(delta))                                 # transfer function of the blur
w = 2 * np.pi * np.fft.fftfreq(N)
L = (2 - 2 * np.cos(w))[:, None] + (2 - 2 * np.cos(w))[None, :]   # eigenvalues of D^T D
H2 = np.abs(H) ** 2; M2 = N // 2
G = 0.25 * (H2[:M2, :M2] + H2[M2:, :M2] + H2[:M2, M2:] + H2[M2:, M2:])  # eigenvalues of A A^T (down2)
def range_part_down2(e):          # orthogonal projection onto the row space of A = S B
    return down_T(np.real(np.fft.ifft2(np.fft.fft2(down(e)) / G)))

res = {}                          # (degradation, lam, quantity) -> list over seeds
checks = {}
r0 = np.random.default_rng(100)   # component checks on random arrays
u, vc = r0.normal(size=(N, N)), r0.normal(size=(M2, M2))
checks["blur transfer function H(0) - 1"] = abs(H[0, 0].real - 1)
checks["down2 adjoint |<Au,v> - <u,A^T v>| / |<Au,v>|"] = abs(np.sum(down(u) * vc) - np.sum(u * down_T(vc))) / abs(np.sum(down(u) * vc))
checks["up_lin and up_rep keep the samples: max error"] = max(np.abs(up_lin(vc)[::2, ::2] - vc).max(), np.abs(up_rep(vc)[::2, ::2] - vc).max())
add = lambda k, v: res.setdefault(k, []).append(v)
for seed in SEEDS:
    rng = np.random.default_rng(seed)
    mask = rng.random((N, N)) < 0.5
    noise = rng.normal(0, SIG_N, (N, N))
    add(("data", 0, "observed fraction"), mask.mean()); add(("data", 0, "noise std"), noise.std())
    ops = {"noise": (ident, ident, x + noise),
           "blur":  (blur, blur, blur(x) + noise),
           "mask":  (lambda z: z * mask, lambda z: z * mask, mask * (x + noise)),
           "down2": (down, down_T, down(x) + noise[::2, ::2])}
    for name, (A, AT, y) in ops.items():
        for lam in LAMS:
            xr = solve(A, AT, y, lam)
            add((name, lam, "right"), psnr(xr))
            ids = {"identity": y} if name != "down2" else {"identity, replicate": up_rep(y), "identity, linear": up_lin(y)}
            for lab, yy in ids.items():
                xi = solve(ident, ident, yy, lam)
                add((name, lam, lab), psnr(xi)); add((name, lam, "gain vs " + lab), psnr(xr) - psnr(xi))
                if name == "noise": add(("noise", lam, "max|right - identity|"), np.abs(xr - xi).max())
            e = xr - x
            if name == "mask":
                add((name, lam, "null-space share"), np.sum(e[~mask] ** 2) / np.sum(e ** 2))
                if seed == 0: checks["mask lam %.3f: max|(D^T D x)_i| on unobserved" % lam] = np.abs(lap(xr)[~mask]).max()
            if name == "down2":
                r = range_part_down2(e)
                add((name, lam, "null-space share"), np.sum((e - r) ** 2) / np.sum(e ** 2))
                if seed == 0: checks["down2 lam %.3f: |P(Pe) - Pe| / |e|" % lam] = np.linalg.norm(range_part_down2(r) - r) / np.linalg.norm(e)
            if name == "blur" and seed == 0:
                xf = np.real(np.fft.ifft2(np.conj(H) * np.fft.fft2(y) / (H2 + lam * L)))
                checks["blur lam %.3f: max|CG - closed-form FFT|" % lam] = np.abs(xr - xf).max()

print("Checks (seed 0):")
for k, v in checks.items(): print(f"  {k:48s} {v:.1e}")
print(f"  {'all CG runs converged':48s} {all(i == 0 for i in infos)}")
fmt = lambda v, p: f"{min(v):.{p}f}" if len(v) == 1 or min(v) == max(v) else f"{min(v):.{p}f} to {max(v):.{p}f}"
print(f"\nObserved fraction of the mask: {fmt(res[('data', 0, 'observed fraction')], 4)};"
      f" noise std: {fmt(res[('data', 0, 'noise std')], 4)}")
print(f"\nPSNR in dB, min to max over seeds {SEEDS}")
for (name, lam, q), v in res.items():
    if name == "data": continue
    p = 4 if "share" in q else 2
    if "max|" in q: print(f"  {name:6s} lam {lam:<6} {q:28s} {max(v):.1e}"); continue
    print(f"  {name:6s} lam {lam:<6} {q:28s} {fmt(v, p)}")

Output (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2; about 2 seconds):

Checks (seed 0):
  blur transfer function H(0) - 1                  0.0e+00
  down2 adjoint |<Au,v> - <u,A^T v>| / |<Au,v>|    2.0e-16
  up_lin and up_rep keep the samples: max error    0.0e+00
  blur lam 0.003: max|CG - closed-form FFT|        6.6e-09
  blur lam 0.030: max|CG - closed-form FFT|        1.4e-09
  mask lam 0.003: max|(D^T D x)_i| on unobserved   2.4e-07
  mask lam 0.030: max|(D^T D x)_i| on unobserved   2.0e-08
  down2 lam 0.003: |P(Pe) - Pe| / |e|              2.9e-14
  down2 lam 0.030: |P(Pe) - Pe| / |e|              3.2e-14
  all CG runs converged                            True

Observed fraction of the mask: 0.4982 to 0.5008; noise std: 0.0200 to 0.0201

PSNR in dB, min to max over seeds [0, 1, 2, 3, 4]
  noise  lam 0.003  right                        34.06 to 34.09
  noise  lam 0.003  identity                     34.06 to 34.09
  noise  lam 0.003  gain vs identity             0.00
  noise  lam 0.003  max|right - identity|        0.0e+00
  noise  lam 0.03   right                        34.78 to 34.82
  noise  lam 0.03   identity                     34.78 to 34.82
  noise  lam 0.03   gain vs identity             0.00
  noise  lam 0.03   max|right - identity|        0.0e+00
  blur   lam 0.003  right                        25.46 to 25.49
  blur   lam 0.003  identity                     26.68 to 26.70
  blur   lam 0.003  gain vs identity             -1.24 to -1.19
  blur   lam 0.03   right                        28.41 to 28.44
  blur   lam 0.03   identity                     26.80 to 26.82
  blur   lam 0.03   gain vs identity             1.60 to 1.64
  mask   lam 0.003  right                        30.98 to 31.06
  mask   lam 0.003  identity                     9.85 to 9.88
  mask   lam 0.003  gain vs identity             21.10 to 21.19
  mask   lam 0.003  null-space share             0.7495 to 0.7548
  mask   lam 0.03   right                        31.06 to 31.14
  mask   lam 0.03   identity                     10.25 to 10.28
  mask   lam 0.03   gain vs identity             20.79 to 20.87
  mask   lam 0.03   null-space share             0.7783 to 0.7831
  down2  lam 0.003  right                        24.72 to 24.87
  down2  lam 0.003  identity, replicate          26.11 to 26.15
  down2  lam 0.003  gain vs identity, replicate  -1.43 to -1.26
  down2  lam 0.003  identity, linear             26.79 to 26.83
  down2  lam 0.003  gain vs identity, linear     -2.10 to -1.96
  down2  lam 0.003  null-space share             0.1670 to 0.1730
  down2  lam 0.03   right                        27.67 to 27.74
  down2  lam 0.03   identity, replicate          26.19 to 26.23
  down2  lam 0.03   gain vs identity, replicate  1.45 to 1.52
  down2  lam 0.03   identity, linear             26.80 to 26.84
  down2  lam 0.03   gain vs identity, linear     0.86 to 0.89
  down2  lam 0.03   null-space share             0.3296 to 0.3348

Every check held: the conjugate-gradient and Fourier solutions agree to within $10^{-8}$, the mask solution is harmonic on the unobserved pixels to the solver’s tolerance, the projector is idempotent to rounding, every run converged, and the mask fraction and noise level are as specified. With noise alone the two models are the same model, and the outputs agree exactly.

Four outcomes were planned.

• The right model never loses to the identity model. This failed at $\lambda = 0.003$ for the blur (gain −1.24 to −1.19 dB) and for down2 (−1.43 to −1.26 dB against replication, −2.10 to −1.96 dB against linear interpolation). It held at $\lambda = 0.03$ for the blur (+1.60 to +1.64 dB) and for down2 (+1.45 to +1.52 and +0.86 to +0.89 dB), at both weights for the mask, and for noise, where the gain is exactly zero. In our reading, the reason is noise amplification, a reading the example does not test. A deconvolving model divides by the blur’s transfer function where the penalty allows it. For the continuous Gaussian, whose transfer function $\exp(-2\pi^2\sigma_b^2 f^2)$ the sampled kernel approximates, $\hat h = 0.169$ at $f = 0.15$ cycles per pixel along one axis, where the eigenvalue of $D^{\mathsf T}D$ is $\hat l = 2 - 2\cos(2\pi f) = 0.824$. The right model’s filter $\hat h / (\hat h^2 + \lambda \hat l)$ then multiplies noise at that frequency by 5.4 at $\lambda = 0.003$ and 3.2 at $\lambda = 0.03$, where the identity model’s $1/(1 + \lambda \hat l)$ passes it at about 1.
• The gain is largest for the mask. This held at both weights: 21.10 to 21.19 dB and 20.79 to 20.87 dB, against at most 0.00 and 1.64 dB for the other degradations. It is largely built into the set-up, since the identity model treats the unobserved half of the image as black, and the example illustrates it rather than tests it.
• More than half of the mask’s error lies in its null space. This held: 0.7495 to 0.7548 at $\lambda = 0.003$ and 0.7783 to 0.7831 at $\lambda = 0.03$, a share the plan had derived before the run from the harmonic fill in the limit of small $\lambda$, so the example illustrates that derivation rather than tests it. The unobserved pixels are half the image and carry about three quarters of the error. The planning research had expected the share to be “≈ all of the error”; the output shows 0.75 to 0.78, not all of it. In our reading the rest is noise that the observed pixels keep.
• More than half of down2’s error lies in its null space. This failed: 0.1670 to 0.1730 at $\lambda = 0.003$ and 0.3296 to 0.3348 at $\lambda = 0.03$, although the null space has three quarters of the dimensions. In our reading, again untested: as $\lambda$ shrinks, the estimate comes to fit the data exactly, $A\hat x = y$, so that $Ae = n$ and the row-space part of the error tends to $A^{\mathsf T}(AA^{\mathsf T})^{-1}n$, the noise divided by eigenvalues of $AA^{\mathsf T}$ that the blur makes small at high frequencies. A large null space does not mean that most of the error is there.

What the example shows is the article’s opening claim in miniature: one prior and one solver give results from about 25 to 35 dB across the four degradations, and which operator the solver is given matters most for the mask, where the identity model takes the missing pixels for black.

The set-up was fixed before the first run, and some of its choices are open. The two values of $\lambda$ come from the planning research, where they were used on a different image, and are printed both rather than tuned. The upsampler for the identity model on down2 is a modelling choice the problem leaves open, so both are printed. In planning, the image was changed from scikit-image’s camera, which the deblurring examples use, to chelsea, and noise was restricted to the observed pixels of the mask after a planning script had also added it to the unobserved ones. The script was not changed after its first run. The example is distinct from the iteration counts that optimization compares on a similar quadratic, from the pseudoinverse of an averaging operator in singular value decomposition, and from the harmonic fill of a hole in image inpainting.

Applications

• Photography: denoising, deblurring and super-resolution in camera pipelines and photo editors; removal of compression artifacts from stored images (image compression).
• Outdoor vision: removing rain, snow, haze and raindrops before detection or segmentation in autonomous driving, video surveillance and drone vision. Fu et al. (2017) and TransWeather’s authors name outdoor vision systems, surveillance and navigation as motivations.
• Scientific and medical imaging: deconvolution in microscopy and astronomy, and denoising of low-dose images, where restored images are measurements and invented detail is a risk (image deblurring, image denoising).
• Training data: degradation models generate training pairs for learned restoration and super-resolution (ESRGAN).

Open Problems

• Real degradations. Networks trained on synthetic rain, snow or noise meet different statistics in real images, where reference images are rarely available for measuring the gap.
• Unknown and combined degradations. AirNet’s authors report that training on more degradations made denoising harder, while deraining and dehazing did not follow the same pattern, call for further study, and leave blur and snow open (Li et al., 2022, §4.4, §6).
• Comparable evaluation. Training sets, test splits, colour channels and retraining protocols differ between papers, so results are comparable only within one paper’s table, and the papers above contain several internal inconsistencies.
• Fidelity and invented detail. Learned priors fill the null space with plausible content; in our reading, the papers discussed here do not report how much of an output is measured and how much is inferred.