Using an Improved Regularization Method and Rigid Transformation for Super-Resolution Applied on MRI Data

1. Introduction

Super-resolution (SR) has gained significant attention in recent years, particularly in medical imaging applications, where the resolution of acquired images is often limited by hardware constraints, time limitations, and patient comfort considerations. Traditional medical imaging modalities such as Magnetic Resonance Imaging (MRI) and Computed Tomography (CT) produce images at a resolution that can restrict the level of detail observable for diagnostic purposes. Increasing the resolution of these images through hardware improvements is often costly and impractical. As a result, computational techniques like SR have emerged as a powerful alternative, allowing high-resolution (HR) images to be reconstructed from low-resolution (LR) inputs without the need for expensive hardware [1,2,3].

The principle behind SR methods is to overcome limitations by leveraging redundant information from multiple LR images or sequences, often involving complex algorithms like regularization methods and machine learning models [4]. Various approaches, from classic interpolation methods to more advanced neural network-based models, have been employed to enhance the quality of medical images in terms of spatial resolution, signal-to-noise ratio (SNR), and edge preservation [5].

In medical imaging, SR is particularly valuable because it enhances the quality of images used in diagnostic processes. For instance, MRI scans are used to assess various medical conditions, and improving their resolution can lead to more accurate diagnoses. Using SR techniques like Wiener filter regularization [1] or edge-preserving high-frequency regularization [3] allows for better visual quality in images without increasing acquisition costs or hardware requirements. Furthermore, neural network-based SR methods [2] have demonstrated promising results in improving image quality with reduced computational time, making them suitable for real-time applications in clinical settings.

The concept of SR in medical imaging has evolved significantly over the years, with various methodologies proposed to tackle the challenges of resolution enhancement. One of the early applications of SR to MRI was proposed by Peled et al. [6], where an Iterative-Back-Projection (IBP) method was used to enhance MRI images of human white matter fiber tracts. While this method showed some promise, it was limited by the use of synthetic image data, which does not fully capture the complexities of real-world medical imaging. Subsequent work by Scheffler [7] addressed this limitation by highlighting the importance of utilizing original image data for more reliable SR reconstruction.

More recently, the integration of machine learning techniques into SR models has shown great promise. For example, a method combining iterative regularization with feed-forward neural networks was proposed by Babu et al. [2], yielding improved results over previous methods due to its capability to handle noise and produce clearer, higher-resolution images. This method demonstrates the potential of neural networks to enhance SR models by reducing computational complexity while maintaining high image quality.

Bayesian methods have also been a major area of exploration in SR research [8,9,10,11,12]. Aguena et al. [1] introduced a Bayesian approach to MRI SR, which employed a Wiener filter to regularize the iterative solution. This method achieved notable improvements in both noise reduction and edge preservation. Similarly, Ben-Ezra et al. [4] proposed a regularized SR framework for brain MRI, incorporating domain-specific knowledge to improve the quality of SR reconstructions. Their approach outperformed traditional maximum a posteriori (MAP) estimators in terms of both edge clarity and overall image quality.

Moreover, Ahmadi and Salari [3] proposed a high-frequency regularization technique that combines edge-preserving methods with traditional SR models. Their approach allows for enhanced edge definition in MRI images without the need for image segmentation, offering a computationally efficient solution suitable for clinical applications.

The Accelerated Proximal Gradient Method (APGM), as outlined in [13], is a well-known optimization technique commonly used for solving inverse problems in imaging, including SR. APGM accelerates the convergence of proximal gradient methods, which are widely adopted for SR tasks involving regularization. Its primary strength lies in its speed, as it converges more quickly than traditional gradient methods, making it suitable for large-scale imaging problems. However, APGM’s effectiveness is heavily dependent on the choice of regularizer, which influences how well the method can balance smoothness and sharpness in the reconstructed image. Poorly chosen regularizers can introduce artifacts or excessively smooth the image. While APGM is flexible and powerful, it requires careful tuning to achieve optimal results, especially when handling high-frequency details.

Block-matching and 3D filtering (BM3D), a renowned denoising algorithm discussed in [14], uses a collaborative filtering approach to reduce noise while preserving image structures. For super-resolution tasks, BM3D can act as a regularizer that effectively manages noise without compromising edges and textures. Its block-matching mechanism compares similar patches in the image, applying 3D filtering to reduce noise in these matched blocks. Although BM3D excels at preserving textures and fine details in natural images, its computational complexity can be high, particularly when dealing with large images or complex noise patterns. Additionally, the block-matching process may struggle in scenarios where image structures do not align well with the blocks, leading to potential loss of detail in areas with intricate textures.

Gu et al. [15] enhanced the BM3D algorithm by introducing weighted nuclear norm minimization, which improved its performance in image denoising. This modification further solidified BM3D as a versatile and widely adopted tool in image processing. Although BM3D is primarily an image denoising algorithm rather than a typical regularization technique, denoising often relies on regularization to reduce noise and improve image quality. BM3D utilizes collaborative filtering and 3D transform-domain techniques to achieve denoising, making it more aligned with advanced signal processing than conventional regularization methods.

Total variation (TV) regularization, a widely used technique in inverse problems, aims to promote sparsity in image gradients, leading to smoother regions while preserving sharp edges. TV regularization is known for its simplicity and its ability to retain edge information, making it a popular choice in SR tasks. However, TV regularization often suffers from the staircasing effect, where smooth regions of the image appear blocky or exhibit artificial edges. Over-regularization can further result in a loss of fine details, which limits the technique’s applicability in images with rich textures or high-frequency content [16].

Rapid and Accurate Image Super Resolution (RAISR), introduced in [17], is a learning-based SR method that is both computationally efficient and fast. It works by learning filters that are adaptive to local image features, such as gradients and edge orientations. Unlike deep learning-based methods that often require significant computational resources, RAISR is lightweight and quick, making it an attractive option for real-time applications. Despite its efficiency, RAISR tends to fall short when compared to more advanced SR methods like deep neural networks in terms of recovering high-frequency details. Its performance is highly dependent on the quality of the learned filters, and it may struggle with images that have complex structures or varying noise levels.

PPPV1, as described in [18], is a video super-resolution method, based on the Plug-and-Play (PnP) framework. This method iteratively refines images, ensuring the gradual recovery of fine details over multiple iterations. A key feature of PPPV1 is its reliance on a denoising module, originally based on DnCNN (Denoising Convolutional Neural Network), to remove noise from the image during each step of the reconstruction process. While DnCNN is effective at reducing noise, it sometimes introduces oversmoothing, especially in high-frequency regions where texture and fine details are critical. The proposed improvement to this method involves replacing DnCNN with a custom prior for denoising. This change allows for more control over detail preservation and texture recovery, potentially reducing the risk of oversmoothing. A well-designed custom prior can provide better balance between noise suppression and sharpness, which could lead to more accurate and visually appealing results, particularly in areas with intricate patterns or high-frequency details.

In addition to neural network-based and Bayesian approaches, convex optimization methods have been explored. Kawamura et al. [5] applied convex optimization techniques to MRI SR, producing state-of-the-art results by carefully balancing noise suppression and detail preservation.

Overall, the field of SR in medical imaging is rapidly advancing, with numerous approaches showing great potential in improving diagnostic imaging and reducing the need for high-cost imaging hardware. The next phase of research will likely focus on integrating these various techniques into more robust, real-time systems suitable for clinical environments.

In this paper, we propose an improved PPP regularization method for MRI super-resolution, utilizing an effective prior specifically designed for denoising and handling motion between frames. Building upon the foundation of our previous method (PPPV1 [18]), our approach incorporates an innovative denoiser that significantly enhances performance. Unlike traditional methods that primarily focus on denoising, our method integrates these advances into an MRI super-resolution framework.

2. Materials and Methods

The acquisition model we are assuming is:

$$\begin{array}{c}\hfill \mathbf{y}=A\mathbf{x}+\epsilon ,\end{array}$$

(1)

where:

• $\mathbf{y}$ is the full set of low resolution (LR) frames, described as $\mathbf{y}={[{{\mathbf{y}}_1}^T,{{\mathbf{y}}_2}^T,\dots ,{{\mathbf{y}}_p}^T]}^T$, where ${\mathbf{y}}_k,k=1,2,...,p$ are the p LR images. Each observed LR image is of size $N_1\times N_2$. Let the kth LR image be denoted in lexicographic notation as ${\mathbf{y}}_k={[y_{k,1},y_{k,2},\dots ,y_{k,M}]}^T$, for $k=1,2,...,p$ and $M=N_1{N}_2$.
• $\mathbf{x}$ is the desired high resolution (HR) image, of size $L_1{N}_1\times L_2{N}_2$, written in lexicographical notation as the vector $\mathbf{x}={[x_1,x_2,...,x_N]}^T$, where $N=L_1{N}_1{L}_2{N}_2$ and $L_1$ and $L_2$ represent the up-sampling factors in the horizontal and vertical directions, respectively.
• $\epsilon ={[{\epsilon }_1,{\epsilon }_2,...,{\epsilon }_p]}^T$, where ${\epsilon }_k$ is the noise vector for frame k and contains independent zero-mean Gaussian random variables.
• $A={[A_1,A_2,...,A_p]}^T$ is the degradation matrix which performs the operations of blur, rigid transformation and subsampling.

Assuming that each LR image is corrupted by additive noise, we can then represent the observation model as [19]:

$$\begin{array}{c}\hfill {\mathbf{y}}_k=A_k\mathbf{x}+{\epsilon }_{k}for1\le k\le p\end{array}$$

(2)

where

$$A_k=S{B}_k{M}_k.$$

(3)

$M_k$ is a matrix of size $L_1{N}_1{L}_2{N}_2\times L_1{N}_1{L}_2{N}_2$ that performs the rigid transformation, $B_k$ represents a $L_1{N}_1{L}_2{N}_2\times L_1{N}_1{L}_2{N}_2$ blur matrix, and S is a $N_1{N}_2\times L_1{N}_1{L}_2{N}_2$ subsampling matrix. In our case $B_k=I$, since we assumed no added blur on video frames.

The goal is to find the estimate $\widehat{\mathbf{x}}$ of the HR image $\mathbf{x}$ from the p LR images $y_k$ by minimizing the cost function

$$\widehat{\mathbf{x}}=arg\underset{\mathbf{x}\in {\mathbb{R}}^N}{min}f\left(\mathbf{x}\right)withf\left(\mathbf{x}\right)=g\left(\mathbf{x}\right)+h\left(\mathbf{x}\right),$$

(4)

where $g\left(\mathbf{x}\right)={\sum }_{k=1}^p{\displaystyle \frac{1}{2}}{\parallel A_k\mathbf{x}-{\mathbf{y}}_k\parallel }_2^2$ is the “fidelity to the data” term, and $h\left(\mathbf{x}\right)$ is the regularization term, which offers some prior knowledge about $\mathbf{x}$. In this study, we adopt the Plug-and-Play Priors approach, in which the ADMM algorithm is modified so that the proximal the proximal operator related to $h\left(\mathbf{x}\right)$ is replaced by a denoiser that solves the problem of Eq. (5). The denoiser used is based on the work by Chantas et al. [20].

The following outlines the algorithm we propose:

1.

The first step of our algorithm is to evaluate the term $M_k$ from the Equation (3), by using rigid registration. Rigid registration, also known as rigid body registration or rigid transformation, is a fundamental technique in medical image processing and computer vision. It is used to align two images by performing translations and rotations while preserving the shape and size of the structures within the images [22].

In a 2D plane, a rigid transformation can be represented using a $3\times 3$ matrix, often referred to as the transformation matrix. For example, a 2D translation can be represented as [23]:

$$T=\left[\begin{array}{ccc}1& 0& t_x\\ 0& 1& t_y\\ 0& 0& 1\end{array}\right]$$

Rotation and reflection matrices can also be formulated similarly. The result of the rigid transformation is represented as an affine transformation matrix. This matrix captures the translation and rotation parameters applied to the original image [23].

We assume that one of the LR images, ${\mathbf{y}}_{mid}$ (typically the middle one), is produced from the HR image $\mathbf{x}$, by applying only downsampling, without transformation. Thus, $M_{mid}=I$. Rigid transformation is calculated between ${\mathbf{y}}_{mid}$ and the rest of the LR images. Following that, we get $M_k$ for the remaining $p-1$ images.

2.

The subsequent phase is centered on employing the PnP-ADMM technique. We execute the PnP-ADMM, adhering to the procedure outlined in Algorithm 1 until reaching convergence, in order to minimize the problem described by Eq. (4). The initial HR image guess, ${\mathbf{x}}^0$, is generated from ${\mathbf{y}}_{mid}$ using the pseudo-inverse of ${\mathbf{A}}_{mid}$. Here, D represents the denoising operator, introduced and discussed in Section 2.1, and g is formulated as $g\left(\mathbf{x}\right)={\sum }_{k=1}^p{\displaystyle \frac{1}{2}}{\parallel A_k\mathbf{x}-{\mathbf{y}}_k\parallel }_2^2$.

Algorithm 1 PnP-ADMM [24]

1: ${\mathbf{u}}^0=\mathbf{0},{\mathbf{x}}^0,and\gamma >0$ 2: for$k=1,2,...,t$ do 3:    ${\mathbf{z}}^k\leftarrow pro{x}_{\gamma g}({\mathbf{x}}^{k-1}-{\mathbf{u}}^{k-1})$ 4:    ${\mathbf{x}}^k\leftarrow D({\mathbf{z}}^k+{\mathbf{u}}^{k-1})$ 5:    ${\mathbf{u}}^k\leftarrow {\mathbf{u}}^{k-1}+({\mathbf{z}}^k-{\mathbf{x}}^k)$ 6: end for 7: return  ${\mathbf{x}}^t$

We next explain the modification made to the standard ADMM algorithm to obtain PnP-ADMM. Line 4 or the standard ADMM is ${\mathbf{x}}^k\leftarrow pro{x}_{\beta h}({\mathbf{z}}^k+{\mathbf{u}}^{k-1})$. In the PnP-ADMM, the proximal operator is replaced by a denoiser D that solves the problem

$$\begin{array}{c}\hfill \mathbf{z}={\mathbf{x}}_0+\mathbf{w}, \mathrm{where} {\mathbf{x}}_0\sim p,\mathbf{w}\sim N(0;\beta I).\end{array}$$

(5)

It can be shown that the Maximum A Posteriori (MAP) estimator $\widehat{{\mathbf{x}}_0}$ of ${\mathbf{x}}_0$ is the proximal operator:

$$\widehat{{\mathbf{x}}_0}=pro{x}_{\beta h}\left(\mathbf{z}\right)=arg\underset{\mathbf{x}\in {\mathbb{R}}^N}{min}\{{\displaystyle \frac{1}{2}}{\parallel \mathbf{x}-\mathbf{z}\parallel }_2^2+\beta h\left(\mathbf{x}\right)\},$$

(6)

for $h\left(\mathbf{x}\right)=-log\left(p\right(\mathbf{x}\left)\right)$.

2.1. The Denoising Algorithm

In this section, we describe the algorithm we use to implement the denoising step of Eq. (6). The algorithm is a simplification of that proposed in [20], it is formulated in a probabilistic (Variational Bayes) context and utilizes an effective prior distribution, which we describe in short next.

2.1.1. The Prior Distribution

The prior distribution we employ for the denoising step was proposed in [20] for single image Super-Resolution, and it is of the form:

$$p\left(\mathbf{x}\right)\propto \prod _{w\in \mathsf{\Omega }}\left(\sum _{\delta \in \mathcal{D}}{\left(1+{\displaystyle \frac{\lambda }{\nu }}{ϵ}_{w,\delta }\left(\mathbf{x}\right)\right)}^{-{\displaystyle \frac{\nu +1}{2}}}\right),$$

(7)

where $\lambda ,\nu$ are the real-positive distribution parameters and ${ϵ}_{w,\delta }$ is a similarity measure between two patches each of center pixel w and $w+\delta$. The above distribution is produced after integrating out the hidden variables of the prior in [20]. However, this form in never explicitly used (it is not necessary) in the optimization algorithm. We show it here in this form for simplicity of presentation. Indeed, $h\left(\mathbf{x}\right)$ enables us to interpret the prior in a deterministic context, analogous with the penalty function imposed on the video frames, see equation (6).

We introduce a similarity measure between two image patches, denoted as ${\mathcal{N}}_w$ and ${\mathcal{N}}_{w^{\prime }}$, where $\mathbf{x}\left(w\right)$ and $\mathbf{x}\left(w^{\prime }\right)$ represent the central pixel of the first and second patch, respectively.

The complete set of pixel coordinates is represented by $\mathsf{\Omega }=\{1,\cdots ,N\}$. Furthermore, we define $\delta$ as the integer displacement between the center pixels of the two patches, such that $w^{\prime }=w+\delta$. For measuring similarity, we employ a weighted Euclidean norm, represented by ${ϵ}_{w,\delta }$, to quantify the difference between ${\mathcal{N}}_w$ and ${\mathcal{N}}_{w^{\prime }}$ (or ${\mathcal{N}}_{w+\delta }$) as follows:

$${ϵ}_{w,\delta }=\sum _{i\in \mathsf{\Omega }}{\mathbf{v}}_{\delta }^2\left(i\right){\gamma }_w\left(i\right),$$

(8)

where ${\mathbf{v}}_{\delta }$ is defined by: ${\mathbf{v}}_{\delta }={\mathbf{Q}}_{\delta }\mathbf{x}$ and ${\mathbf{v}}_{\delta }^2$ indicates the vector obtained by squaring each element of ${\mathbf{v}}_{\delta }$. ${\mathbf{Q}}_{\delta }$ represents the difference operator, an $N\times N$ matrix, such that the i-th component of ${\mathbf{Q}}_{\delta }\mathbf{x}$ equals $\mathbf{x}\left(i\right)-\mathbf{x}\left(i^{\prime }\right)$ for all $i,i^{\prime }\in \mathsf{\Omega }$ with $i^{\prime }-i=\delta$. The matrix ${\mathbf{G}}_w$ is an $N\times N$ diagonal matrix, where its diagonal elements corresponding to the pixels in ${\mathcal{N}}_w$ are the only non-zero values, specifically, ${\mathbf{G}}_w(i,i)=0$ for all i not in ${\mathcal{N}}_w$. Lastly, we denote by ${\gamma }_w$ the $N\times 1$ vector with elements the weights of the weighted norm: the closer to the central pixel of the patches the larger the weight value.

The norm defined by (8) retains its value even if the summation (8) runs over only the subset ${\mathcal{N}}_w\subset \mathsf{\Omega }$ instead of $\mathsf{\Omega }$, since ${\gamma }_w\left(i\right)=0$ for $i\notin {\mathcal{N}}_w$. However, we use the full summation range over $\mathsf{\Omega }$ for enabling fast computations with the Fast Fourier Transform, as explained next.

The distance between the patch ${\mathcal{N}}_{w=1}$ and an arbitrary patch ${\mathcal{N}}_{w^{\prime }}$, $w^{\prime }\in \mathsf{\Omega }$, is $\delta =w-w^{\prime }=1-w^{\prime }$. Given that the image patches correspond to ${\mathbf{g}}_1$ and ${\mathbf{g}}_{w^{\prime }}$, it is:

$${\mathbf{g}}_{w^{\prime }}\left(i\right)={\mathbf{g}}_{w=1}\left(i-\delta \right)={\mathbf{g}}_1\left(i+1-w\right),\forall i\in \mathsf{\Omega }.$$

(9)

As we can see, each ${\gamma }_{w^{\prime }}$, is a circularly shifted by $w^{\prime }$ version of ${\mathbf{g}}_1\equiv \mathbf{g}$ (denoted simply by $\mathbf{g}$ from now on). The formula (8) for calculating ${ϵ}_{w,\delta }$, expressed in terms of $\gamma$, is:

$${ϵ}_{w,\delta }=\sum _{i\in \mathsf{\Omega }}{\mathbf{v}}_{\delta }^2\left(i\right){\mathbf{g}}_w\left(i\right)=\sum _{i\in \mathsf{\Omega }}{\mathbf{v}}_{\delta }^2\left(i\right)\mathbf{g}(i+1-w).$$

(10)

Clearly, the values of ${ϵ}_{w,\delta }$ for all w’s, are the result of the correlation between ${\mathbf{v}}_{\delta }^2$ and $\mathbf{g}$, since the indices of ${\mathbf{v}}^2$ and $\mathbf{g}$ always differ by the constant $1-w$. To calculate the correlation required for the super-resolution technique discussed in the following section, we use the Fast Fourier Transform (FFT). This approach decreases the computational complexity of the algorithm from $O\left(N^2\right)$, typical for correlation calculations, to $O(NlogN)$, which is the complexity for multiplication in the DFT (Discrete Fourier Transform) domain.

2.1.2. Denoising in PnP-ADMM

Next, we describe the algorithm we employ in the PnP-ADMM context of Algorithm 1, and specifically for the denoising step (line 4). The algorithm we employ, as a denoising sub-problem of the general super-resolution algorithm (Algorithm 1), is in essence a special case of the VBPS algorithm in [20], where there is no blurring nor decimation. Mathematically speaking, this means that the imaging operator $\mathbf{DH}$ is the $N\times N$ identity matrix $\mathbf{I}$, as shown in line 8 of Algorithm 2.

Algorithm 1: Variational Bayes Patch Similarity Denoising

Input: Noisy image ${\mathbf{z}}^k+{\mathbf{u}}^{k-1}$. Output: Denoised image ${\mathbf{x}}^k$. Initialization: Image initial estimate: Set ${\alpha }_{\mathrm{new}}=\alpha /2$, where $\alpha$ is the regularization parameter obtained from [21]. Then, set ${\mathbf{m}}^{\left(0\right)}={\mathbf{x}}_{\mathrm{Stat}}$, where ${\mathbf{x}}_{\mathrm{Stat}}$ is the super-resolved image obtained after setting $\alpha ={\alpha }_{\mathrm{new}}$. Parameter selection: Set $t=0$, and $\beta ={N/\parallel \mathbf{x}-\mathbf{z}\parallel }_2^2$, $\lambda ={10}^3{\alpha }_{\mathrm{new}}$, $\nu =7$, $rmax=280$, $\mathrm{MAXITER}=25$ and err $={10}^{-7}$. 1: while$\parallel {\mathbf{m}}^{\left(t\right)}-{\mathbf{m}}^{(t-1)}{\parallel }_2^2/N>err$ AND $t<\mathrm{MAXITER}$ do 2:    for every $\delta$ in $\mathcal{D}$ do 3:      ${\mathbf{v}}_{\delta }\leftarrow {\mathbf{Q}}_{\delta }{\mathbf{m}}^{\left(t\right)}$ 4:      for every w in $\mathsf{\Omega }$ do 5:         Calculate the expectations of the following model’s random variables: $${\langle a_{w,\delta }\rangle }_{\left(t\right)}={\displaystyle \frac{1+\nu }{\lambda {\widehat{ϵ}}_{w,\delta }+\nu }},$$ $${\langle z_{w,\delta }\rangle }_{\left(t\right)}={\displaystyle \frac{e^{-\frac{\lambda }{2}{\langle a_{w,\delta }\rangle }_{\left(t\right)}{\widehat{e}}_{w,\delta }-\frac{\nu }{2}log{\langle a_{w,\delta }\rangle }_{\left(t\right)}}}{{\sum }_{\delta }{e}^{-\frac{\lambda }{2}{\langle a_{w,{\delta }^{\prime }}\rangle }_{\left(t\right)}{\widehat{e}}_{w,{\delta }^{\prime }}-\frac{\nu }{2}log{\langle a_{w,{\delta }^{\prime }}\rangle }_{\left(t\right)}}}},$$ where $\widehat{ϵ}$ is the $ϵ$ in (8), calculated with the image estimation provided in the previous iteration $t-1$, 6:         calculate ${\mathbf{b}}_{\delta }^{\left(t\right)}\left(w\right)={\langle a_{w,\delta }\rangle }_{\left(t\right)}{\langle z_{w,\delta }\rangle }_{\left(t\right)}$, for all w and $\delta$, 7:         set ${\Lambda }_{\delta }^{\left(t\right)}=\mathrm{diag}\{{\mathbf{b}}_{\delta }^{\left(t\right)}\ast \mathbf{g}\}$ (convolution), 8:         $t\leftarrow t+1$ 9:         Obtain ${\mathbf{m}}^{\left(t\right)}$ by solving the linear system $\left(\beta \mathbf{I}+\lambda {\sum }_{\delta }{\mathbf{Q}}_{\delta }^T{\Lambda }_{\delta }^{\left(t\right)}{\mathbf{Q}}_{\delta }\right){\mathbf{m}}^{\left(t\right)}=\beta \mathbf{y}$ with the Conjugated Gradients algorithm. 10:      end for 11:    end for 12: end while 13: T=t; ${\mathbf{x}}^k={\mathbf{m}}^{\left(T\right)}$.

More specifically, the imaging model assumed for the denoising step is a simplified form of Eq. (2.1) in [20], because it is now $\mathbf{DH}=\mathbf{I}$ (i.e., no blur/decimation, hence it is just the identity matrix). Also, in this form, ${\mathbf{z}}^k+{\mathbf{u}}^{k-1}$ has the role of the “noisy image” and ${\mathbf{x}}^k$ is the uncorrupted one, meant to be estimated by the denoising algorithm.

In parallel with imaging model, we assume the imaging model, i.e., the prior distribution introduced above and given by Eq. (5). This is in essence the prior distribution for the uncorrupted image to be estimated via the denoising procedure. This means that the Algorithm 2 is the result of the adoption of both the imaging model mentioned above and the prior (5) for $\mathbf{x}$. Lastly, note that the denoising Algorithm 2 selects automatically, in the initialization step, the noise variance $\beta$, among other parameters.

We implemented our method in SCICO [25], which is an open source library for computational imaging that includes implementations of several algorithms.

To evaluate our method, the widely-used publicly available dataset named the cancer image archive (TCIA) [26] was used, in order to compare our results to the previously proposed method. Specifically, we conducted experiments using a dataset of LR brain MRI images and a corresponding HR reference dataset.

3. Results

Our method with the effective prior achieved notable improvements in image quality, as demonstrated by Figure 1 and Figure 2.

To objectively evaluate the effectiveness of our improved technique, we calculated the PSNR and conducted comparisons with both alternative approaches and enhanced versions of our own method. Specifically, we compared against PPPV1 [18], APGM (accelerated proximal gradient method) [13], BM3D (Block-matching and 3D filtering) [16], Total Variation [14], RAISR (Rapid and Accurate Image Super Resolution) [17] and MIRNetv2 [27], as well as with the pseudo-inverse and the denoised pseudo-inverse images. The difference between PPPV1 and the currently proposed method is that now we use a custom prior instead of DnCNN for the denoising, while we use the same rigid transformation. The outcomes, detailed in Table 1, unequivocally demonstrate that our method surpasses others in delivering higher image quality.

The Wilcoxon signed-rank test was used to compare the PSNR values of the proposed method with the respective values for PPP V1, Pseudo-inverse, Denoised Pseudoinverse, APGM, BM3D and TV methods. The results obtained with those statistical tests are shown in Figure 3 and indicated statistically significant differences between the PPP and the other six methods, since no per-slice data was available for RAISR and MIRNetv2.

Considering the perceptual quality of the frames, it is obvious from Table 2 that our method gives the best results, outperforming the other methods on all datasets. This outcome demonstrates the robustness and effectiveness of our method in enhancing the natural quality of super-resolved videos for this specific dataset.

4. Discussion

In summary, the results presented in this study highlight the superior performance of the proposed method in the field of video super-resolution. This method consistently outperforms state-of-the-art techniques, as demonstrated by the substantial PSNR gains observed on the datasets used for evaluation. The following key takeaways can be drawn:

• The experimental results demonstrate the superiority of our approach over existing techniques, underscoring its potential for clinical applications in neuroimaging.
• The practical implications of our results suggest that our method holds great promise for applications where MRI slices quality enhancement is paramount.
• Computational efficiency is another significant advantage of our method. Unlike Deep Neural Network-based methods, our approach does not rely on neural networks and requires no training, making it faster and less resource-intensive.

These findings make a strong case for the adoption of our method in MRI enhancement and upscaling tasks. We believe that the approach we suggest has the potential to contribute significantly to the field of video super-resolution and benefit a wide range of applications.