Hybrid CUDA acceleration of wavelet-based MRI image denoising

1 Introduction

The analysis of medical images using deep learning and neural networks has become a highly relevant and important research area, driving advancements in diagnosis and treatment across nearly all fields of medicine (Suganyadevi et al., 2022; Gao et al., 2025; Masoudi et al., 2021). These technologies enable early disease detection and support clinical decision-making in medical practice (Masoudi et al., 2021; Wang et al., 2025). The use of artificial intelligence has demonstrated promising results in medical image classification, segmentation, and quantitative analysis (Moreno López et al., 2021; Li et al., 2023).

Noise is one of the primary factors that degrade the quality of medical images and directly affect diagnostic accuracy (Kaur et al., 2021; Alcaín et al., 2021; Goyal et al., 2018). It leads to the loss of fine structural details, reduced contrast, and increased difficulty in visual interpretation. Therefore, there is a critical need for methods capable of effectively suppressing noise while preserving anatomical structures. In magnetic resonance imaging (MRI), common noise sources include thermal Gaussian noise, Rician noise, and artifacts caused by patient motion or magnetic field instability (Patil and Bhosale, 2022; Moeller et al., 2021). Studies have shown that preliminary noise reduction significantly improves classification accuracy; for instance, it increases from 68.74% to 73.30% for MRI and from 77.72% to 82.28% for CT (Masoudi et al., 2021).

Traditional filtering approaches, such as linear smoothing filters or Wiener filtering, can partially reduce noise but tend to blur anatomical boundaries (Patil and Bhosale, 2022; You et al., 2023). Median filters preserve edges more effectively; however, their performance deteriorates in low signal-to-noise ratio conditions. Modern methods, including Non-Local Means and anisotropic diffusion, provide higher reconstruction quality but are computationally expensive, which limits their applicability in real-time scenarios (Yuan, 2015; Cappabianco et al., 2019).

Wavelet-based filtering is widely used as an effective noise reduction technique because of its multiscale signal representation. It allows effective separation of noise from meaningful structures and preserves anatomical edges (Taassori, 2024; Patil and Bhosale, 2023). Wavelet transforms are also important for edge detection, as they separate high-frequency components (boundaries) from low-frequency background information (Parida and Bhoi, 2017). However, the widespread use of wavelet-based methods is limited by their high computational cost, which constrains practical deployment, particularly in latency-sensitive image-processing pipelines. Thus, in the context of this study, the main limitation of wavelet denoising is not methodological opacity or dependence on training data, but the computational cost of decomposition, thresholding, and reconstruction, which motivates the CUDA-based acceleration proposed in this work.

In parallel with classical and wavelet-based approaches, recent years have seen rapid progress in deep learning and hybrid architectures for medical image processing, including MRI (Zhang et al., 2023; Kang et al., 2024). CNN–Transformer hybrid models have demonstrated strong performance in medical image segmentation (Zhang et al., 2023), while lightweight CNN designs with pruning and quantization have been explored to reduce inference latency (Sriwiboon, 2025b). Learning-based MRI analysis and reconstruction methods have also shown strong performance in related tasks such as meningioma grading and 3D MRI super-resolution (Zhang et al., 2024; Kang et al., 2024), and transformer-based segmentation has been applied in clinical imaging studies (Sriwiboon and Phimphisan, 2025a). However, these methods belong to a data-driven design space. Their performance depends on the availability and representativeness of training data, the correspondence between training and target acquisition conditions, model adaptation procedures, and hardware-specific inference behavior. These requirements may be difficult to satisfy in settings where paired or sufficiently representative training data are unavailable, scanner protocols vary, or the denoising stage must remain transparent and reproducible.

For this reason, the present work focuses on wavelet-based denoising as a deterministic and interpretable alternative, as well as a possible preprocessing component for learning-based pipelines. In contrast to trained denoisers, the discrete wavelet transform does not require a training stage, paired reference data, or model retraining when the acquisition setting changes. Its multiresolution representation separates low-frequency anatomical structures from high-frequency detail and noise components, while the thresholding rule provides an explicit and mathematically controlled mechanism for noise suppression. In addition, the inverse transform preserves a well-defined reconstruction pathway, which is important for reproducibility and auditability in medical image-processing workflows. Therefore, the wavelet-based approach was selected in this study not as a replacement for deep learning in all medical imaging tasks, but as a transparent, training-data-independent denoising core whose computational bottleneck can be directly addressed through CUDA-level execution optimization.

Preprocessing remains relevant in medical imaging pipelines, including in combination with learning-based methods. For example, the study in (Mochurad, 2025) demonstrated that edge-detection-based preprocessing can improve segmentation quality by preserving fine structural details and anatomical boundaries. This observation supports the broader practical value of efficient image enhancement and denoising stages, although the present work focuses specifically on accelerating a classical wavelet-based denoising method rather than on modifying downstream learning models.

Processing speed is critically important for clinical systems, where decisions must be made in real time. Slow algorithms can cause delays that reduce the effectiveness of medical workflows and patient treatment (Bnou et al., 2020). Therefore, there is a need for parallel algorithms capable of processing large volumes of medical data with minimal latency. The use of graphics processing units (GPUs), field-programmable gate arrays (FPGAs), and application-specific integrated circuits (ASICs) makes it possible to accelerate image filtering and analysis processes (Wu et al., 2025).

Prior studies have also demonstrated the feasibility of parallel acceleration for computationally intensive sensor-processing and MRI-related image-analysis tasks, including LiDAR localization and brain MRI analysis (Mochurad and Kryvinska, 2021; Mochurad and Mochurad, 2025). These results support the view that GPU-oriented optimization can substantially reduce execution time in computationally intensive medical image-processing pipelines. In this context, wavelet-based denoising is a suitable target for implementation-level optimization because it combines established restoration behavior with a structured sequence of computational stages, including decomposition, thresholding, and reconstruction. In a straightforward GPU realization, these stages are typically executed as separate data-parallel kernels with intermediate results repeatedly written to and read from global memory. Such an organization limits concurrency between computational branches and increases synchronization overhead.

The objective of this study is to develop and evaluate a hybrid GPU-based method for accelerating a fixed wavelet-based MRI denoising pipeline based on the discrete wavelet transform (DWT) by combining data-level and task-level parallelization. The denoising core itself remains classical and is based on a single-level 2D discrete wavelet transform with adaptive soft-thresholding, whereas the main contribution of the work lies in the execution-level organization of this method on the GPU.

The contribution of this study includes the following elements:

• We developed a hybrid parallel GPU execution method for accelerating the complete DWT–thresholding–inverse DWT pipeline for MRI image denoising. In contrast to straightforward data-parallel implementations, the proposed method combines fine-grained data-level parallelism within CUDA kernels with task-level parallelism between the low- and high-frequency processing branches within a unified execution framework;
• We proposed an execution-level reorganization of the vertical filtering and inverse reconstruction stages based on dual-stream scheduling and fused reconstruction kernels. This organization reduces synchronization overhead and intermediate global-memory traffic and thereby improves the throughput of the wavelet denoising pipeline relative to a conventional sequential-stage data-parallel realization;
• We developed a memory-aware execution strategy integrated into the proposed hybrid method, including constant-memory placement of filter coefficients, symmetric boundary handling directly within kernels, and selective loop unrolling. This strategy improves memory-access efficiency and contributes to the stable runtime behavior of the computational stages;
• We performed an experimental validation of the proposed method against a reproduced data-parallel GPU baseline with sequential stage execution and showed that the proposed hybrid organization provides stable per-image runtime together with a 20%–30% throughput improvement while preserving the reconstruction behavior of the underlying wavelet denoising procedure.

The scientific novelty of the study lies in the development of a hybrid parallel method for GPU-based execution of wavelet MRI denoising, in which the standard DWT–thresholding–inverse DWT procedure is reorganized into a unified CUDA execution workflow. In this workflow, pixel-wise filtering is performed by data-parallel kernels, independent low- and high-frequency wavelet branches are scheduled on separate CUDA streams, inverse vertical reconstruction is also executed by independent CUDA streams, and the final horizontal reconstruction is performed using a fused kernel. Unlike existing GPU-based wavelet denoising implementations that primarily map decomposition, thresholding, and reconstruction onto consecutive data-parallel CUDA stages, the proposed method introduces a coordinated execution-level combination of data-level parallelism, task-level branch scheduling, dual-stream vertical-pass processing, dual-stream inverse reconstruction, fused horizontal reconstruction, constant-memory filter storage, and in-kernel symmetric boundary handling within one DWT–thresholding–inverse DWT workflow. This integrated organization distinguishes the proposed method from the reproduced data-parallel baseline of (Wu et al., 2025), because it reduces stage-level serialization, synchronization points, and intermediate global-memory traffic while preserving the mathematical form and reconstruction behavior of the underlying wavelet denoising model.

The practical significance of the proposed approach lies in enabling high-throughput wavelet-based MRI denoising with stable runtime and preserved reconstruction quality, which supports its use in latency-sensitive medical image-processing pipelines. The method was evaluated on 2D MRI slices under controlled Gaussian, Rician, and motion-related degradations, providing a reproducible proof-of-concept basis for computational benchmarking. However, the reported reconstruction metrics should not be interpreted as a direct estimate of clinical denoising performance on real volumetric or multi-scanner MRI acquisitions.

2 State-of-the-art

In the scientific literature, the main directions relevant to MRI image denoising include classical filtering techniques, wavelet-based methods, learning-based approaches, and parallel implementations of image-processing pipelines. Each of these categories is briefly discussed below.

2.1 Classical filtering methods

The most common classical denoising approaches include Gaussian, mean, median, Wiener, and bilateral filters (Taassori, 2024; Singh, 2016; Lew et al., 2024; Nitin and Gupta, 2022; Sun, 2023). Their advantages are simplicity of implementation and high processing speed, which makes them suitable for clinical use. However, their main limitation is reduced reconstruction quality: Gaussian and mean filters tend to blur anatomical edges (Singh, 2016; Wu et al., 2025), the median filter is effective for impulse noise but performs poorly for Gaussian noise in MRI (Wu et al., 2025), whereas Wiener and bilateral filters preserve edges more effectively but exhibit high computational cost, limiting scalability (Taassori, 2024; Sun, 2023).

2.2 Wavelet-based approaches

Wavelet filtering (DWT, stationary wavelet transform (SWT)) is an effective tool for multiscale signal analysis and preservation of anatomical structures (Taassori, 2024; Patil and Bhosale, 2023; Parida and Bhoi, 2017; Bnou et al., 2020; Ismael and Baykara, 2021; Sachin and Dharmpal, 2011). These methods excel at separating noise from useful image content, thereby maintaining structures and edges. At the same time, their practical use is often limited by the computational cost of sequential decomposition, thresholding, and reconstruction stages, especially when large image sets must be processed under latency constraints.

2.3 Deep learning methods

Deep learning methods have become an important direction in medical image denoising because they can learn complex spatial dependencies and image priors directly from data. CNN-based denoising models, residual networks, encoder–decoder architectures, U-Net-like models, attention-based networks, and CNN–Transformer hybrids are commonly used to suppress noise while preserving anatomical structures. These approaches can achieve high reconstruction quality in many reported settings, often reflected in improved PSNR and SSIM values, because they are able to model nonlinear relationships between noisy and clean image patterns and can adapt to complex noise distributions, including signal-dependent noise and structured artifacts (Ilesanmi and Ilesanmi, 2021; Kaur et al., 2023). They can also improve the accuracy of downstream tasks such as classification and segmentation by enhancing diagnostically relevant structures before automated analysis (Kaur et al., 2023; Mechria et al., 2022).

At the same time, learning-based denoising introduces several limitations that are important in medical imaging. These models usually require large and representative training datasets, and their performance may decrease when the target images differ from the training data because of scanner type, acquisition protocol, anatomy, pathology, or noise distribution. Supervised denoising methods may also require paired noisy and clean training data, while self-supervised or unpaired approaches still require sufficiently representative training distributions, which are not always available for clinical MRI. In addition, model training is computationally expensive, hyperparameter-dependent, and often requires dedicated GPU resources. Although inference can be fast after training, runtime behavior still depends on the target hardware, model size, and optimization framework. Another limitation is the reduced interpretability of deep models: the restoration process is less transparent than in transform-based methods, and in sensitive medical applications this may complicate reproducibility, auditability, and the assessment of possible texture smoothing or structure modification.

CNN-based frameworks have also demonstrated strong performance across broader medical image analysis tasks, including segmentation scenarios where modern architectures such as DeepLabv3+ are commonly used (Sriwiboon, 2025a). Image enhancement and preprocessing remain important even in CNN pipelines; for example, CLAHE-assisted CNN designs have been reported to improve the visibility of clinically relevant patterns in retinal imaging (Phimphisan and Sriwiboon, 2024). Beyond imaging, efficient CNN–Transformer hybrids have also been explored for other biomedical modalities, reflecting the broader trend toward compact and explainable models (Sriwiboon and Phimphisan, 2025b). Therefore, learning-based denoising methods represent a powerful and rapidly developing research direction, but they are methodologically different from the deterministic wavelet-based approach considered in this work. In the present study, they are discussed as an important related direction, while the main focus remains on accelerating a fixed, interpretable DWT–thresholding–inverse DWT denoising pipeline on the GPU.

2.4 Parallel algorithm implementations

The adoption of graphics processing units (GPUs) and CUDA technology has opened new opportunities for accelerating image processing pipelines. Prior work (Wu et al., 2025) demonstrated that parallel wavelet-based denoising algorithms can significantly reduce processing time. However, most existing solutions rely on synchronous data processing, which limits scalability (Wu et al., 2025). Hybrid strategies combining data-level and task-level parallelism have shown superior performance and better scalability in computational imaging and signal-processing pipelines, particularly when multiple processing stages must be executed concurrently. Prior studies demonstrate that integrating data-parallel execution with task-parallel scheduling improves throughput, reduces bottlenecks, and enables more efficient GPU and CPU utilization in medical imaging and related high-performance workflows (Kirimtat and Krejcar, 2024). However, the proposed work differs from (Wu et al., 2025) and related CUDA implementations of wavelet denoising in the level at which parallelism is organized (Wu et al., 2025). demonstrated the possibility of accelerating wavelet denoising by mapping wavelet operations onto CUDA kernels; however, such an implementation remains predominantly stage-oriented and data-parallel. In these approaches, decomposition, thresholding, and reconstruction are usually executed as consecutive GPU stages with synchronization between them. In contrast, the proposed method reorganizes the computational graph of the fixed DWT–thresholding–inverse DWT pipeline. The method introduces dual-stream scheduling for independent low- and high-frequency branches, concurrent execution of vertical filtering and inverse vertical reconstruction, fused horizontal reconstruction kernels, storage of filter coefficients in constant memory, and symmetric boundary handling directly inside CUDA kernels. Therefore, the proposed contribution consists in GPU-level reorganization of an established denoising procedure aimed at reducing synchronization overhead and intermediate global-memory traffic.

Based on the above analysis, the following conclusions can be drawn:

• classical methods are fast but provide relatively low reconstruction quality;
• wavelet-based approaches offer a trade-off between accuracy and performance, yet their sequential nature limits real-time use;
• deep neural networks can learn complex noise and image priors and often achieve high reconstruction quality, but their practical use depends on representative training data, model adaptation, interpretability requirements, and stable hardware-specific inference;
• parallel GPU-based methods offer real-time potential but require optimized parallel execution strategies.

Based on the limitations observed in prior studies, this work focuses on addressing the following research problem:

Despite significant progress in classical, wavelet-based, deep learning, and GPU-accelerated denoising approaches, existing solutions still face several critical limitations:

• Sequential wavelet implementations are too slow for large-scale clinical deployment;
• Deep neural networks require extensive training and computational resources, limiting real-time applicability;
• Existing GPU-based methods primarily rely on data-parallel strategies, which restrict scalability and underutilize GPU resources.

Therefore, the research problem addressed in this work is the development of a computationally efficient hybrid GPU-based method for accelerating a fixed wavelet MRI denoising pipeline while preserving the reconstruction behavior of the underlying denoising procedure. Unlike straightforward CUDA realizations based only on sequential stage execution of data-parallel kernels, the proposed approach reorganizes the DWT–thresholding–inverse DWT workflow through hybrid task-and-data parallel execution, dual-stream scheduling, and fused processing stages.

Accordingly, this study addresses the following research question:

RQ: Can a hybrid CUDA-based implementation of wavelet denoising for 2D brain MRI provide near-real-time slice-level throughput while preserving the reconstruction behavior of the underlying method and maintaining stable runtime across different controlled MRI noise and artifact models?

Primary endpoint:

• throughput ≥200 images/s;
• peak signal-to-noise ratio (PSNR) in the range of approximately 29–30 dB and structural similarity index (SSIM) around 0.90 under Gaussian noise in the evaluated experimental setting.

3 Materials and methods

3.1 Problem statement

Let I x , y denote a noisy 2D brain MRI slice defined on the discrete grid Ω = 1 , … , M × 1 , … , N, where I x , y is the intensity at pixel x , y. Let I ∼ x , y be the unknown noise-free image and I ^ x , y the restored estimate produced by the denoising operator.

We assume the presence of additive Gaussian noise Equation 1.where I ∼ : Ω → R – the true noise-free image. Here η x , y ∼ N 0 , σ 2 denotes additive Gaussian noise.

I x , y = I ∼ x , y + η x , y , (1)

The objective is to reconstruct an estimate I ^ of the true image I ∼ by minimizing the mean squared error with an additional regularization term that preserves structural details (see Equation 2):where λ > 0 is a regularization parameter and R u denotes the L ₁ regularization of wavelet coefficients to preserve anatomical structures.

I ^ = arg min u : Ω → R I − u 2 2 + λ R u , (2)

For the numerical implementation, the regularization term is expressed through the discrete wavelet transform (DWT) and thresholding, as shown in Equation 3.where W and W − 1 denote the forward and inverse DWT, respectively, and T τ is the soft-thresholding operator with threshold τ.

I ^ = W − 1 ∘ T τ ∘ W I , (3)

The regularization term can be written as Equation 4.

R u = W u 1 . (4)

3.2 Dataset and preprocessing

The experimental evaluation was performed using the Brain Tumor Classification (MRI) dataset (Hemant, 2021). The dataset contains two-dimensional brain MRI images corresponding to three tumor categories: meningioma, glioma, and pituitary tumor. In this study, all 3,064 available grayscale MRI slices were used for controlled denoising experiments. Each image was converted into a single-channel 512 × 512 image in float32 format, and intensity values were scaled to the 0–255 range for quantitative evaluation.

The dataset was selected because it provides a sufficiently large and reproducible set of 2D brain MRI slices for computational benchmarking of the proposed denoising pipeline. However, the dataset is originally intended for classification rather than for raw MRI reconstruction or clinical denoising validation. Therefore, paired noisy and noise-free clinical MRI acquisitions are not available. For this reason, the available images were treated as reference slices, and controlled degradations were generated synthetically to evaluate reconstruction behavior and runtime under reproducible conditions.

Three degradation models were considered: additive Gaussian noise, simulated Rician noise, and simplified motion-induced artifacts. Gaussian noise with standard deviation σ = 5–25 was added to model additive thermal and electronic noise components. Rician noise was simulated using the magnitude MRI model by adding two independent Gaussian components to the real and imaginary channels before magnitude reconstruction. Motion-related artifacts were approximated by rigid in-plane transformations, including random translations and rotations.

Rician noise was simulated using the standard magnitude MRI model. Since the available dataset contains real-valued grayscale magnitude images rather than raw complex MRI data, each reference slice I 0 x , y was treated as the real component of an underlying complex signal with zero imaginary component. Two independent Gaussian noise components were then added to the real and imaginary channels, and the final corrupted image was obtained by the magnitude operation given in Equation 5: where I 0 x , y denotes the real-valued reference grayscale MRI slice before simulated noise addition, and n 1 x , y and n 2 x , y are independent zero-mean Gaussian components n 1 , n 2 ∼ N 0 , σ R 2 representing noise in the real and imaginary channels of the complex MRI signal. In this simulation, the imaginary component of the reference image was assumed to be zero because raw complex-valued MRI data were not available. The noise standard deviation σ R = 5–25 in Table 1 was selected to match the noise-level range used in the Gaussian experiments. This simulation provides a controlled approximation of signal-dependent magnitude MRI noise and allows the method to be evaluated beyond purely additive Gaussian corruption.

Z R x , y = I 0 x , y + n 1 x , y + i n 2 x , y ,

I Ri с ian x , y = Z R x , y = I 0 x , y + n 1 x , y 2 + n 2 2 x , y . (5)

TABLE 1. Denoising performance of the proposed hybrid method under different MRI noise sources.

Motion-induced artifacts were generated by applying rigid in-plane transformations, including random horizontal and vertical shifts of up to ±15 pixels and rotations of up to ±10°, using linear interpolation to simulate structured spatial distortions. This procedure represents a simplified image-domain approximation of motion-related degradation. It does not reproduce the full physics of MRI motion artifacts, which may arise from k-space trajectory inconsistencies and can appear as ghosting, blurring, phase-related distortions, or spatially varying artifacts in clinical scans. Therefore, this simulation was used only as a controlled degradation model for methodological comparison on the available 2D dataset.

The selected degradation models were intended to approximate three different aspects of MRI noise formation under controlled and reproducible conditions. Additive Gaussian noise approximates thermal and electronic noise components that can appear in the acquisition chain. The Rician model provides a closer approximation to magnitude MRI data, where noise becomes signal-dependent after reconstruction from complex-valued measurements. The simulated motion artifacts approximate structured spatial degradation caused by patient movement, although they do not reproduce the full k-space acquisition process. Therefore, these models capture several important statistical and spatial characteristics of clinical MRI degradation, but they remain simplified approximations rather than complete clinical noise models. It should be emphasized that the present evaluation was designed as a controlled computational benchmark rather than as a full clinical validation study. The use of synthetic Gaussian noise, simulated Rician noise, and image-domain motion artifacts allowed the CPU implementation, the reproduced data-parallel GPU baseline, and the proposed hybrid CUDA implementation to be compared under identical and reproducible degradation conditions. However, real in vivo MRI acquisitions may contain spatially varying and scanner-dependent artifacts arising from coil sensitivity profiles, B0 and B1 field inhomogeneities, parallel imaging reconstruction, gradient nonlinearities, scanner-specific preprocessing, and patient-dependent motion patterns. These effects were not fully reproduced by the degradation models used in this study. Therefore, the reported PSNR, SSIM, and MSE values should be interpreted as controlled reconstruction indicators for the evaluated dataset and noise models, rather than as direct evidence of clinical denoising efficacy across scanners, protocols, and acquisition settings.

The same preprocessing protocol, image size, grayscale float32 representation, wavelet basis, thresholding rule, boundary handling strategy, timing protocol, and hardware environment were used for the CPU implementation, the reproduced data-parallel GPU baseline, and the proposed hybrid CUDA implementation. This ensured that the reported differences in runtime and throughput were caused by the proposed CUDA execution organization rather than by differences in input data preparation or denoising parameter tuning.

3.3 Discrete wavelet transform

The 1D discrete wavelet transform is defined as Equations 6, 7:

a k = ∑ i = 0 L − 1 h L i x 2 k − i , (6)

d k = ∑ i = 0 L − 1 h H i x 2 k − i . (7)

Here, h L i and h H i denote the low-pass and high-pass decomposition filter taps.

To perform denoising, a single-level 2D discrete wavelet transform (2D-DWT) with biorthogonal filters (bior 4.4) is applied. In this study, the biorthogonal 4.4 (bior4.4) wavelet was selected due to several properties that are advantageous for medical image denoising. First, bior4.4 provides near-symmetric analysis and synthesis filters, which helps reduce phase distortion and preserves anatomical edge geometry more accurately than asymmetric orthogonal wavelets such as Daubechies. Second, its compact support ensures efficient implementation with short filters, lowering computational cost and memory bandwidth requirements on the GPU. Third, biorthogonal wavelets offer linear-phase reconstruction, which is desirable for medical imaging because it avoids edge shifting and maintains the spatial alignment of fine structures. Compared with Haar wavelets, which are extremely compact but introduce block-like artifacts due to discontinuous basis functions, bior4.4 yields smoother approximations and better preserves soft-tissue gradients. These properties make the bior4.4 basis a suitable compromise between reconstruction fidelity, computational efficiency and stability in a real-time MRI denoising pipeline.

The image I is decomposed into four subbands:

• LL – low-frequency approximation,
• LH, HL, HH – high-frequency horizontal, vertical, and diagonal detail components.

The forward DWT is defined as Equations 8–11. where h L and h H are the low-pass and high-pass filter coefficients; L is the filter length; i , j are the coordinates of the subbands. In addition to the analysis filters h L k and h H k, the inverse discrete wavelet transform requires the corresponding reconstruction filters, denoted as h L r k and h H r k. These filters are used in the synthesis stage of the inverse DWT. The exact analysis and reconstruction filter coefficients used in this work are provided in Table 2.

L L i , j = ∑ k = 0 L − 1 ∑ l = 0 L − 1 h L k h L l I 2 i + k , 2 j + l , (8)

L H i , j = ∑ k = 0 L − 1 ∑ l = 0 L − 1 h L k h H l I 2 i + k , 2 j + l , (9)

H L i , j = ∑ k = 0 L − 1 ∑ l = 0 L − 1 h H k h L l I 2 i + k , 2 j + l , (10)

H H i , j = ∑ k = 0 L − 1 ∑ l = 0 L − 1 h H k h H l I 2 i + k , 2 j + l , (11)

TABLE 2. Exact biorthogonal (bior 4.4) filter coefficients used in all experiments: decomposition low-pass h L, decomposition high-pass h H reconstruction low-pass h L r, and reconstruction high-pass h H r.

Downsampling by a factor of 2 is applied in both directions. At image boundaries, symmetric reflection is used, as shown in Equation 12.

I − k = I k , I N + k = I N − k − 1 . (12)

3.4 Soft thresholding

Soft thresholding is applied to the high-frequency subbands to suppress noise while preserving significant structural details. The soft-thresholding operator is defined as in Equation 13:where the threshold is computed according to Equation 14:with α as the adaptation coefficient and β b as the subband-specific weight factor. The noise level σ is estimated according to Equation 15:where M H H is the median absolute deviation of the HH subband. In the single-level 2D DWT configuration used in this work, each detail subband is assigned its own threshold value. Three different thresholds τ L H, τ H L, τ H H are computed according to Equation 14, allowing the denoising strength to be adapted to the statistical properties of each subband. The explicit threshold definitions are: τ L H = α β L H σ, τ H L = α β H L σ, τ H H = α β H H σ. In all experiments, the subband weight factors are set to β L H = 0.8, β H L = 0.8, β H H = 1.2, so that the H H subband receives stronger thresholding, while the L H and H L subbands, which contain oriented anatomical edges, are denoised more conservatively. In this work, the soft-thresholding rule itself follows standard wavelet denoising formulations, while its subband parameterization is adapted to the considered brain MRI setting and integrated into the hybrid CUDA-based execution pipeline described in Section 3.6.

c ′ = s i g n c max c − τ , 0 , (13)

τ b = α σ β b , b ∈ L H , H L , H H , (14)

σ = M H H 0.6745 (15)

The denoising parameters were fixed before the comparative performance evaluation in order to isolate the effect of the proposed GPU execution strategy. The bior4.4 wavelet basis and single-level 2D DWT configuration were used in all experiments. This configuration was selected to preserve anatomical edge geometry while limiting intermediate memory allocation, synchronization overhead, and runtime variability for 512 × 512 MRI images.

The threshold values were computed adaptively from the median absolute deviation of the HH subband and then scaled using subband-specific factors. This design avoids manually selecting a separate threshold for each image or noise type. The same wavelet basis, decomposition level, thresholding rule, adaptation coefficient, and subband weights were used in all experiments, including Gaussian, Rician, and motion-related degradations, as well as in the CPU, reproduced data-parallel GPU, and proposed hybrid CUDA implementations. This ensured that the performance comparison reflected differences in execution organization rather than differences in denoising parameter tuning.

3.5 Inverse transform

In the inverse transform, vertical synthesis is performed first using the reconstruction low-pass and high-pass filters h L r and h H r. The vertical reconstruction is defined as Equations 16, 17.

L v i , j = ∑ k h L r k L L i − k , j + ∑ k h H r k L H i − k , j , (16)

H v i , j = ∑ k h L r k H L i − k , j + ∑ k h H r k H H i − k , j . (17)

This step produces the intermediate low-frequency and high-frequency components L v and H v, each of size H × W.

The horizontal reconstruction is then applied to obtain the denoised image I ^ (see Equation 18).

I ^ i , j = ∑ k h L r k L v i , j − k + ∑ k h H r k H v i , j − k . (18)

The same symmetric boundary extension rule described in Equation 12 is applied consistently in both the vertical and horizontal reconstruction stages.

3.6 Hybrid CUDA-based parallelization

To achieve high processing performance, a hybrid parallelization strategy is employed, combining both data-level and task-level parallelism on the GPU.

3.6.1 Data parallelism

Each pixel i , j is processed independently. The thread coordinates are defined as in Equation 19.where b x , b y are block indices, t x , t y are thread indices inside a block, and B x = B y = 16 was selected as the main configuration used in the experiments. Configurations of 8 × 8 and 32 × 32 were also evaluated.

i , j = b x · B x + t x , b y · B y + t y (19)

3.6.2 Task parallelism

Low- and high-frequency filtering is executed concurrently within a single computational process.

• − Horizontal filtering: both low- and high-pass filters are applied simultaneously to each row. Results are stored in separate buffers (Lo and Hi components). Filter coefficients are stored in constant GPU memory (four arrays, 10 elements each) to minimize latency and allow broadcast access by all threads.
• − Vertical filtering: performed by two parallel kernel calls on independent CUDA streams: one for Lo components (forming LL and LH), another for Hi components (forming HL and HH). Each kernel computes two outputs simultaneously with loop unrolling for reduced iteration count.

3.6.3 Inverse reconstruction

• − Vertical reconstruction processes the pairs (L L , L H) and (H L , H H) in parallel on independent CUDA streams using pitched memory allocation (cudaMallocPitch) for optimal coalesced access.
• − Horizontal reconstruction combines Lo and Hi components into the final image.

Synchronization between DWT stages is achieved using cudaStreamSynchronize (). Adaptive thresholding for L H , H L, and H H is performed in parallel, with threshold values dynamically computed based on noise variance and scale factors. To make the integration of data-level and task-level parallelism explicit, the proposed CUDA implementation is organized as a dependency graph of kernels rather than as an independent set of acceleration techniques. Data-level parallelism is used inside each CUDA kernel, where individual threads compute output samples of the corresponding convolution, thresholding, or reconstruction operation. Task-level parallelism is introduced only when the data dependencies allow independent execution of wavelet branches. In particular, after horizontal filtering, the low-frequency branch Lo and the high-frequency branch Hi are independent and are therefore assigned to different CUDA streams during the vertical DWT pass. The same principle is used during inverse vertical reconstruction, where the pairs (LL, LH) and (HL, HH) are processed concurrently. Table 3 summarizes the kernel-level organization of the proposed pipeline, including input and output buffers, stream assignment, and synchronization points.

TABLE 3. Kernel-level organization and data flow of the proposed hybrid CUDA pipeline.

The main synchronization points are placed after the forward vertical DWT streams and after the inverse vertical reconstruction streams. Thus, the Lo and Hi branches are not serialized during the vertical passes. Instead, the effective latency of these stages is governed by the slower branch and the stream synchronization overhead. The fused horizontal inverse kernel is launched only after both inverse vertical branches have completed. This execution order preserves the mathematical form of the single-level 2D DWT, adaptive soft-thresholding, and inverse DWT, while reducing stage-level serialization and intermediate global-memory traffic compared with a conventional sequential-stage CUDA implementation.

This hybrid strategy enables efficient pipeline execution and high-throughput processing on the evaluated hardware platform. The block diagram shown in Figure 1 illustrates the image processing pipeline on the GPU, highlighting the distribution of subtasks and the simultaneous use of data parallelism at all processing stages.

The block diagram (see Figure 1) illustrates the distribution of individual subtasks, such as horizontal and vertical convolution, inverse transform, and adaptive thresholding, while data parallelism is maintained at each processing stage, ensuring efficient utilization of GPU resources. To clarify the distinction between the proposed implementation and a conventional data-parallel CUDA baseline, the baseline executes the DWT, thresholding, and inverse DWT as separate sequential GPU stages. The proposed implementation preserves the same denoising operations but changes their execution schedule. In particular, independent low- and high-frequency branches are assigned to separate CUDA streams, vertical filtering and inverse vertical reconstruction are executed concurrently where data dependencies allow, and the final horizontal reconstruction is performed by a fused kernel. This organization reduces the number of synchronization points and decreases intermediate global-memory transfers without changing the mathematical form of the denoising operator.

The reduction in synchronization cost follows from the independence of the low- and high-frequency branches after horizontal filtering. In a sequential stage-oriented CUDA implementation, vertical filtering of the Lo and Hi branches and the corresponding inverse vertical reconstruction stages are executed one after another, with synchronization between consecutive kernels. In the proposed scheme, these independent branches are assigned to separate CUDA streams, so their execution can overlap where data dependencies allow. As a result, the effective latency of these stages is determined by the longer branch plus stream synchronization overhead, rather than by the sum of both branch runtimes.

The fused reconstruction kernel reduces memory-bandwidth pressure by avoiding an additional write-read cycle for intermediate reconstruction buffers. Instead of storing partially reconstructed low- and high-frequency components in global memory and launching a separate kernel to combine them, the proposed implementation performs filtering and accumulation within a single kernel. Together with constant-memory storage of the short bior4.4 filter coefficients and in-kernel symmetric boundary handling, this reduces intermediate global-memory traffic, kernel-launch overhead, and stage-level serialization.

Algorithm 1 presents the proposed hybrid CUDA-based wavelet denoising procedure, including parallel filtering, thresholding, and reconstruction steps.

Algorithm 1 Kernel-level hybrid CUDA workflow for wavelet-based MRI denoising.

In the single-level 2D DWT setting used in this work, the input image I has size H × W. The horizontal filtering step produces two intermediate arrays Lo and Hi of size H × W. After vertical filtering and downsampling by a factor of two in each direction, each subband LL, LH, HL, HH has size H 2 × W 2. The intermediate arrays used during the inverse transform, L v and H v, have size H × W, and the reconstructed image I ^ also has size H × W.

In all experiments, we use a two-dimensional CUDA block of blockDim = (16, 16, 1), and the grid dimensions are computed as gridDim_row = (ceil (W/16), ceil (H/16), 1) for horizontal filtering and gridDim_col = (ceil (W/16), ceil (H/16), 1) for vertical filtering and reconstruction kernels, where ceil (·) denotes rounding up to the nearest integer. Two CUDA streams are employed during inverse vertical reconstruction: one stream processes the pair (LL, LH’), and the second stream processes the pair (HL’, HH’). Filter coefficients are stored in constant GPU memory, and the same symmetric boundary extension as in Equation 12 is implemented inside each kernel.

The proposed method combines data-parallel execution, where each pixel is processed by a dedicated CUDA thread, with task-parallel execution, in which low- and high-frequency branches as well as the LL and LH subbands and the HL and HH subbands are processed in parallel CUDA streams. In this context, task-level parallelism refers specifically to the concurrent execution of independent wavelet branches on separate CUDA streams, whereas data-level parallelism refers to pixel-wise convolution and thresholding inside CUDA kernels. The use of dual CUDA streams reduces synchronization overhead, while constant-memory caching and loop unrolling improve memory-access efficiency.

4 Modeling and results

4.1 Evaluation metrics and experimental setup

To evaluate the denoising quality, the PSNR and SSIM metrics are used. PSNR is defined as Equation 20.where RMSE is the root mean square error between the reference and denoised images. In this work, float32 images are scaled to the 0–255 range. SSIM measures structural similarity between images based on luminance, contrast, and structural information. Both metrics are computed using the scikit-image library. In addition to these reconstruction metrics, runtime, throughput, and uncertainty estimates are reported in order to characterize both the computational and reconstruction behavior of the proposed method under controlled degradations.

P S N R = 20 · ⁡ log 10 255 R M S E , (20)

Consistent with the research question (RQ), the primary endpoint of the study was high slice-level denoising throughput together with preserved reconstruction quality under Gaussian noise in the evaluated experimental setting. Specifically, the primary endpoint was defined as achieving a processing speed of at least 200 images per second together with PSNR in the range 29–30 dB and SSIM around 0.90 on the Brain Tumor MRI dataset.

Secondary endpoints were: stability of per-image runtime when varying the dataset size; robustness of PSNR and SSIM under Rician noise and motion-induced artifacts; and relative speedup over a CPU implementation and over an existing GPU data-parallel wavelet baseline.

Experiments were conducted in Google Colab using an Intel Xeon CPU and an NVIDIA Tesla T4 GPU (compute capability 7.5), with CUDA 12.5 and OpenCV 4.5.4. The GPU implementation was developed using CUDA C/C++ and the CUDA Runtime API. GPU kernels were used for horizontal and vertical wavelet filtering, adaptive thresholding, and inverse reconstruction. CUDA streams were used for task-level scheduling, while cudaEvent timing and std::chrono::high_resolution_clock were used for kernel-level and end-to-end runtime measurements, respectively. OpenCV was used for image loading, preprocessing, and format conversion, whereas scikit-image was used for PSNR and SSIM computation. The implementation details are summarized in Algorithm 1, which provides reproducible pseudocode of the complete DWT–thresholding–inverse DWT workflow. To ensure a fair baseline comparison, the data-parallel GPU wavelet denoising baseline described by (Wu et al., 2025) was reproduced as a stage-oriented CUDA implementation of the same DWT–thresholding–inverse DWT pipeline. The reproduced baseline used the same input images, image size, grayscale float32 representation, bior4.4 wavelet filters, single-level 2D DWT configuration, adaptive soft-thresholding rule, symmetric boundary handling, preprocessing steps, timing protocol, and hardware environment as the proposed method. Where implementation-level details were not explicitly specified in the reference study, the baseline was implemented conservatively as a sequential stage-oriented data-parallel CUDA workflow, in which decomposition, thresholding, and inverse reconstruction were executed as separate GPU kernels with synchronization between stages, without dual-stream branch scheduling or fused reconstruction. This ensured that the comparison isolated the effect of the proposed execution-level reorganization rather than differences in denoising parameters, input data, or hardware conditions. The detailed dataset description and preprocessing protocol are provided in Section 3.2. In the present section, we report the experimental setup, evaluation metrics, degradation scenarios, runtime measurements, and comparative results obtained using the dataset and preprocessing procedure described above.

In addition to Gaussian noise, two other MRI degradation sources, Rician noise and motion-induced artifacts, were also evaluated to assess the robustness of the proposed method under different noise conditions.

All three noise models were evaluated on the full Brain Tumor Classification (MRI) dataset (Hemant, 2021), ensuring consistent and comparable testing conditions across noise types. The average PSNR and SSIM values computed over all images and noise levels are reported in Table 1. To complement PSNR and SSIM, the mean squared error (MSE) was also computed as a direct pixel-wise distortion measure between the reference and denoised images. MSE provides an additional quantitative estimate of reconstruction error and helps interpret PSNR values in terms of absolute intensity-domain deviations. Lower MSE values indicate lower residual distortion after denoising.

As shown in Table 1, the proposed method exhibits consistent reconstruction behavior across all three MRI noise categories, although the differences between noise types are more pronounced than in purely additive noise settings. Gaussian noise yields the highest reconstruction quality, with PSNR exceeding 30 dB and SSIM above 0.92. Rician noise leads to a moderate reduction in performance due to its signal-dependent nature, lowering PSNR by approximately 2.4 dB relative to the Gaussian case. The additional MSE values confirm the same tendency as PSNR: Gaussian noise produces the lowest reconstruction error, whereas Rician noise and motion-related artifacts lead to larger residual deviations. The motion case has the highest MSE because rigid spatial distortions introduce structured pixel-wise differences that are more severe than purely additive noise. At the same time, SSIM remains relatively stable, indicating that part of the anatomical structure is preserved despite the increased pixel-wise error.

For Rician noise and motion-induced artifacts, the results remain lower than for the Gaussian case but stable across repeated trials, indicating that the proposed implementation preserves predictable reconstruction behavior under multiple controlled degradation types. The uncertainty analysis (Table 4; Figure 2) further supports this observation by showing narrow confidence intervals and compact per-case distributions for both PSNR and SSIM.

TABLE 4. PSNR and SSIM mean ± std and 95% confidence intervals across 5 runs and all subjects.

These results indicate that the proposed system remains computationally applicable when the input is affected by different controlled degradation types, because the same DWT–thresholding–inverse DWT workflow and the same CUDA execution schedule are used in all cases. However, the denoising effectiveness is degradation-dependent. The method performs best for additive Gaussian noise, remains stable but less accurate for Rician noise, and shows the largest residual error for motion-related artifacts because these artifacts introduce structured spatial distortions rather than only high-frequency random noise.

Execution time of GPU operations was measured using two approaches: CUDA events for kernel-level timing and std::chrono::high_resolution_clock for end-to-end runtime, including cudaDeviceSynchronize (). Table 5 reports the latency of the core processing stages, whereas the end-to-end per-image runtime and throughput values used for the main performance analysis are reported separately in Tables 4, 6. Results were saved in CSV files for further analysis. The algorithm is optimized for GPUs with compute capability ≥6.0, ensuring compatibility with modern architectures.

TABLE 5. Timing breakdown of the core processing stages of the proposed GPU implementation on Tesla T4.

TABLE 6. Comparison of total execution time and performance gain for different parallel approaches.

To improve the auditability of the runtime evaluation, we provide a decomposition of the total latency into its main computational components. The pipeline consists of CPU preprocessing, host-to-device data transfer, GPU kernel execution and device-to-host data transfer. All measurements were performed on the same hardware used in the experiments, namely, an Intel Xeon CPU and an NVIDIA Tesla T4 GPU. The timing breakdown is summarized in Table 5, which presents the average latency of the core processing stages for a 512 × 512 input image on the evaluated hardware platform.

As shown in Table 5, GPU kernel execution accounts for most of the latency associated with the core processing stages, while data-transfer overheads remain relatively small due to the single-level architecture and the fixed filter size. This breakdown complements the end-to-end runtime analysis reported later in Tables 4, 6 and provides a transparent view of the low-level execution profile within the tested configuration.

These timing measurements indicate that the runtime reduction is mainly associated with GPU-side execution rather than host-device transfer overhead. The H2D and D2H transfer times are small compared with the complete end-to-end runtime, whereas the kernel execution stage is the dominant part of the core GPU processing path. Therefore, the main optimization target is the organization of GPU kernels and intermediate memory access. The comparison with the reproduced data-parallel GPU implementation further supports this interpretation: the proposed method reduces the per-image runtime from 5.96 ms to 4.65 ms for the full dataset in end-to-end mode, which corresponds to a reduction of approximately 22%. This improvement is consistent with the expected effect of overlapping independent wavelet branches and reducing intermediate global-memory traffic through fused reconstruction. To clarify how the proposed execution strategies affect synchronization costs and memory-bandwidth usage, Table 7 summarizes the execution-level role of each CUDA optimization used in the proposed implementation.

TABLE 7. Execution-level interpretation of the proposed CUDA optimizations.

To provide quantitative support for the claimed memory-aware execution behavior, an additional CUDA profiling analysis was performed on the NVIDIA Tesla T4 GPU for 512 × 512 MRI slices. The profiling focused on kernel-level execution time, achieved occupancy, global memory throughput, and the relative contribution of each kernel group to the total GPU processing time. The measurements were obtained using CUDA event timing and profiling counters collected over repeated runs after warm-up iterations. The same input size, bior4.4 filters, block configuration, boundary handling rule, and thresholding procedure as in the main experiments were used. The results are summarized in Table 8.

TABLE 8. CUDA profiling results for the proposed hybrid implementation on Tesla T4.

The sum of the mean kernel times in Table 8 is 0.250 ms/image, which is consistent with the core GPU kernel execution time reported in Table 5. Occupancy and memory throughput are reported at the kernel-group level because these metrics characterize individual kernel behavior and are not directly additive across heterogeneous CUDA kernels.

The profiling results indicate that the vertical DWT filtering and inverse vertical reconstruction stages are the dominant GPU-side components, accounting for 24.4% and 23.2% of the total kernel execution time, respectively. This is consistent with the memory-access pattern of vertical filtering, where access to neighboring samples is less favorable than in row-wise filtering. The achieved occupancy values remain between 68% and 86% across the kernel groups, indicating that the implementation provides sufficient thread-level parallelism without severe underutilization of the GPU. The highest global memory throughput is observed for the dual-stream vertical DWT stage, confirming that this part of the pipeline is mainly memory-bandwidth sensitive.

The fused horizontal inverse reconstruction kernel reduces the need for an additional intermediate global-memory write-read cycle and accounts for 16.4% of the GPU kernel time. This supports the design choice of combining final filtering and accumulation in one kernel. Overall, the profiling results show that the runtime improvement of the proposed implementation is associated not only with parallel branch scheduling but also with reduced intermediate memory traffic and more compact reconstruction execution. The remaining bottlenecks are primarily related to vertical memory access and inverse vertical reconstruction, which suggests that future optimization may focus on shared-memory tiling, improved memory coalescing for column-wise operations, and architecture-specific tuning of block dimensions.

To evaluate the effectiveness of the proposed hybrid approach, experimental studies were conducted. All experiments were repeated five times to reduce variability, and the reported values represent the mean execution time with a standard deviation below 10% for CPU runs and below 5% for GPU runs. To quantify uncertainty in reconstruction quality, each experiment was repeated five times with independent noise realizations. For every noise level and every subject, PSNR and SSIM were computed as mean ± standard deviation. The 95% confidence intervals were computed using C I = x ¯ ± 1.96 s n s, where x ¯ is the mean PSNR and SSIM over all evaluated samples, s is the corresponding standard deviation, and n is the number of independent observations used to estimate the statistics. These values are reported in Table 4. Per-image distributions for all subjects are visualized using boxplots (Figure 2), providing insight into inter-subject variability and robustness across anatomical cases.

Figure 2 shows per-subject PSNR and SSIM distributions for Gaussian, Rician, and motion-induced noise. Boxplots summarize median performance, interquartile ranges, and outliers, demonstrating consistent behavior across subjects and confirming that the method remains stable across the evaluated anatomical cases.

The speedup factor S p was computed as Equation 21, where T s is the end-to-end CPU runtime and T p is the end-to-end runtime of the proposed hybrid GPU method.

S p = T s T p , (21)

The results are presented in Table 9.

TABLE 9. Comparison of total execution time and GPU speedup factor for different dataset sizes.

The obtained results demonstrate that the GPU implementation significantly outperforms the CPU version in terms of processing speed. For a single image, a maximum acceleration factor of 124× was achieved, while for larger batches the speedup stabilizes at approximately 68–70×. Importantly, this acceleration does not alter the reconstruction produced by the underlying denoising procedure. This is confirmed by the PSNR and SSIM values for both CPU and GPU implementations presented in Table 10, which show nearly identical results.

TABLE 10. Comparison of reconstructed image quality metrics for CPU and GPU implementations.

A visual demonstration of the GPU implementation is presented using an example image from the Brain Tumor Classification (MRI) dataset before and after applying the proposed method (Figure 3). The images illustrate a direct comparison between CPU and GPU processing on the same data under identical visualization settings.

To support the visual comparison between the CPU and GPU outputs, the corresponding CPU–GPU difference heatmap is presented in Figure 4.

The visual assessment complements the quantitative PSNR, SSIM, and MSE results. In Figure 3, the CPU and GPU denoising results are visually almost identical, indicating that the GPU implementation accelerates the computation without altering the reconstruction behavior of the underlying wavelet denoising procedure. Figure 4 further supports this observation by showing that the CPU–GPU differences are confined to negligible floating-point deviations. All qualitative results were generated using identical window and level settings to ensure consistent visual comparison.

To compare with existing approaches, the data-parallel algorithm described in (Wu et al., 2025) was reproduced within the same experimental environment. The comparison was performed as a controlled within-study reimplementation rather than as a direct comparison with absolute runtimes reported in the original paper, because runtime values depend strongly on hardware, software stack, memory configuration, and implementation details. Both the reproduced baseline and the proposed method were tested on identical image sets using the same wavelet basis, thresholding rule, boundary handling, preprocessing pipeline, data type, timing protocol, and hardware platform. The baseline reconstruction output was also checked against the CPU implementation to ensure that the reproduced data-parallel GPU version preserved the same denoising behavior up to negligible floating-point differences. The comparison results with the proposed hybrid method are presented in Table 6. The performance gain S A over the data-parallel GPU baseline was computed as in Equation 22, where T D P denotes the runtime of the reproduced data-parallel GPU implementation and T H Y B denotes the runtime of the proposed hybrid method.

S A = T D P − T H Y B T D P × 100 % . (22)

As seen in Table 6, the hybrid approach achieves lower execution time across all test sets, with the performance gain ranging from 4.17% (for 5 images) to 24.68% (for 100 images).

Figure 5 shows the dependence of total processing time on the number of images:

• CPU (blue curve) compared to hybrid GPU approach (orange curve).
• hybrid GPU approach (orange curve) compared to data-parallel GPU algorithm (green curve).
• As shown in Figure 5, with an increasing number of images, the processing time of the CPU implementation grows linearly and significantly exceeds the execution time of the GPU implementation using the hybrid approach. Furthermore, the hybrid method provides a lower execution time compared to the GPU data-parallel implementation from (Wu et al., 2025). For a more detailed analysis, Table 11 presents the comparison of the per-image processing time depending on the total number of processed images for different algorithm versions.
• The results in Table 11 show that, when evaluated in end-to-end mode, the processing time per image remains almost constant for the GPU implementations: approximately 4.4–4.6 m for the hybrid method and 4.8–6.2 m for the data-parallel method, regardless of the dataset size.
• Overall, the experimental results show that the proposed hybrid implementation provides substantial acceleration relative to the CPU version and lower latency than the reproduced data-parallel GPU baseline, while preserving the denoising output of the underlying wavelet formulation and maintaining stable reconstruction quality across the evaluated conditions.

TABLE 11. Comparison of per-image processing time for different algorithm implementations across varying dataset sizes.

These experiments evaluate scalability with respect to dataset size on the tested Tesla T4 platform. The results show that the per-image runtime of the proposed hybrid implementation remains nearly stable when the number of processed images increases from a single image to the full dataset of 3,064 images. Therefore, the scalability claim in this study refers to dataset-size scalability and stable end-to-end runtime behavior under the evaluated hardware configuration. Empirical validation across multiple GPU architectures was not performed in the present work and is therefore not claimed as an experimentally demonstrated result. Since the implementation relies on standard CUDA mechanisms, including CUDA streams, constant memory, and fused kernels, it is expected to be portable to modern CUDA-capable GPUs; however, quantitative multi-architecture benchmarking is left for future work.

4.2 Sensitivity analysis of wavelet basis and threshold parameters

To examine the possible dependence of reconstruction quality and runtime behavior on the selected wavelet and threshold configuration, an additional sensitivity analysis was performed. This analysis was not used to tune the proposed method for the main comparison with the CPU implementation and the reproduced GPU baseline. Instead, it was designed to assess whether the selected bior4.4 configuration provides a stable balance between reconstruction accuracy, structural preservation, and computational efficiency.

Three groups of parameters were varied: the wavelet basis, the threshold scaling coefficient, and the effective filter length. The tested wavelet bases included Haar, Daubechies-4 (db4), Symlet-4 (sym4), Coiflet-2 (coif2), and biorthogonal 4.4 (bior4.4). These bases were selected because they differ in filter length, symmetry properties, compactness of support, and reconstruction characteristics. The threshold scaling coefficient α was varied in the range 0.6–1.4 with a step of 0.2, while the subband weights βLH, βHL, and βHH were kept fixed in order to isolate the effect of global threshold strength. For all configurations, the same input images, noise levels, preprocessing steps, boundary handling rule, CUDA block size, timing protocol, and hardware platform were used. PSNR, SSIM, MSE, per-image runtime, and throughput were recorded for each configuration.

The sensitivity analysis was performed under Gaussian noise because it provides a controlled degradation model and allows direct comparison of reconstruction metrics across parameter settings. The purpose of this analysis was not to identify a universally optimal wavelet basis for all MRI acquisition conditions, but to verify that the selected bior4.4 basis and the adopted threshold range provide stable reconstruction behavior and a consistent execution profile for the considered CUDA-based implementation. The results of the wavelet-basis analysis are summarized in Table 12.

TABLE 12. Sensitivity of denoising performance to wavelet basis and effective filter length under Gaussian noise.

As shown in Table 12, the wavelet basis affects both reconstruction metrics and runtime behavior. The Haar basis provides the lowest per-image runtime because of its shortest filter, but it yields lower PSNR and SSIM values, indicating weaker structural preservation. Wavelets with longer and smoother filters improve reconstruction metrics, although they require more convolution operations and therefore increase per-image runtime. Among the tested configurations, bior4.4 provides the most stable balance between PSNR, SSIM, MSE, and throughput. It achieves the highest PSNR and SSIM while maintaining throughput above 200 images/s, which is consistent with the primary endpoint of the study.

To further assess the role of threshold strength, an additional analysis was performed for the bior4.4 basis by varying the global threshold scaling coefficient α while keeping the subband-specific weights fixed. The results are presented in Table 13.

TABLE 13. Sensitivity of reconstruction metrics to the threshold scaling coefficient α for bior4.4 under Gaussian noise.

As shown in Table 13, changing α mainly affects reconstruction metrics, whereas runtime remains nearly unchanged. This is expected because threshold scaling changes the coefficient-shrinkage strength but does not alter the number of CUDA kernels, memory transfers, filtering operations, or synchronization points. Lower α values preserve more high-frequency components but leave more residual noise, whereas higher α values suppress noise more aggressively and may smooth weak anatomical details. The best reconstruction behavior was observed at α = 1.0, which corresponds to the fixed threshold configuration used in the main experiments. Overall, the sensitivity analysis confirms that the selected bior4.4 configuration is not an isolated parameter choice, but a stable setting that preserves reconstruction quality while maintaining the execution profile required for high-throughput CUDA-based denoising.

4.3 Discussion

The experimental results provide a positive answer to the research question formulated in Section 1. The primary endpoint, which was defined as achieving high slice-level throughput together with preserved reconstruction quality under Gaussian noise, was met: the proposed hybrid method processes 215–220 images per second while reaching PSNR in the range of 29–30 dB and SSIM around 0.90 on the Brain Tumor MRI dataset. At the same time, the GPU implementation substantially outperforms the CPU version in terms of processing speed: for a single image, the maximum acceleration achieved is 124×, and for larger datasets the speedup stabilizes at 68–70×, with negligible loss of precision compared to the CPU baseline.

The observed acceleration has practical implications for clinical image-processing workflows. The achieved throughput of 215–220 slices per second and the stable per-image runtime of approximately 4.4–4.6 ms indicate that the proposed method can be used as a lightweight preprocessing stage before visualization, segmentation, classification, or quantitative analysis modules in latency-sensitive MRI pipelines. In practical terms, such acceleration may reduce waiting time during batch preprocessing of MRI slices and allow denoising to be performed close to image acquisition or before downstream analysis without becoming the dominant computational bottleneck. Since the proposed method preserves the classical wavelet denoising operator and does not require model training, it can be integrated as an implementation-level module into existing CPU and GPU medical image-processing workflows. However, its current clinical applicability should be interpreted at the 2D slice-processing level; routine integration into volumetric clinical MRI workflows would require additional validation on in vivo data, Digital Imaging and Communications in Medicine (DICOM) and Picture Archiving and Communication System (PACS)-compatible processing, and extension to inter-slice or 3D denoising. The absence of direct validation on multi-scanner in vivo MRI data is an important limitation of the present study. The current experiments demonstrate that the proposed CUDA execution strategy accelerates a fixed wavelet-based denoising pipeline while preserving its reconstruction behavior under controlled degradation conditions. They do not, however, establish scanner-independent clinical generalizability. A complete clinical validation would require evaluation on raw or minimally preprocessed in vivo MRI data acquired using different scanner vendors, field strengths, pulse sequences, reconstruction pipelines, and anatomical protocols. Such validation should include both quantitative and qualitative criteria: runtime and throughput, no-reference image quality measures when clean references are unavailable, radiologist- or expert-based assessment of anatomical detail preservation, evaluation of spatially varying residual noise, and the effect of denoising on downstream segmentation or classification tasks. Therefore, the present results should be interpreted as computational evidence of feasibility and efficiency, while multi-scanner in vivo validation remains a necessary next step before routine clinical deployment.

To complement the direct experimental evaluation, a direct comparison with baseline methods implemented on the same dataset was conducted (Table 14). This provides an objective reference point for evaluating the measured runtime benefit and reconstruction behavior of the proposed approach.

TABLE 14. Direct performance comparison with baseline methods on the Brain Tumor MRI dataset.

As shown in Table 14, the proposed hybrid method achieves a 20%–30% runtime improvement over the reproduced data-parallel GPU baseline while maintaining PSNR of 29.8 dB and SSIM of 0.90 in the evaluated setting.

The runtime advantage of the proposed method should be interpreted as a result of more efficient GPU execution organization. Both the reproduced data-parallel GPU baseline and the proposed method are based on the same general paradigm of wavelet denoising. The difference lies in how the computational stages are scheduled and how intermediate data are processed on the GPU. The proposed dual-stream and fused-kernel organization reduces stage-level serialization and intermediate global-memory traffic, which explains the lower per-image processing time reported in Table 11. At the same time, the PSNR and SSIM values remain within the expected range because the denoising operator itself is not replaced by a different restoration model.

The speed–accuracy trade-off should be interpreted at two levels. At the execution level, the proposed CUDA organization improves runtime without changing the mathematical denoising operator. Therefore, the reduction in processing time is not obtained by simplifying the reconstruction model, but by more efficient scheduling and memory handling on the GPU. This is supported by the close PSNR and SSIM values between the CPU and GPU implementations and by the negligible CPU–GPU difference shown in the heatmap. At the algorithmic level, reconstruction accuracy depends on the selected wavelet basis and threshold parameters. The sensitivity analysis in Tables 12, 13 shows that this dependence is systematic and interpretable. Shorter filters reduce per-image runtime but may provide weaker structural preservation, whereas smoother and near-symmetric wavelet bases improve PSNR and SSIM at a moderate increase in runtime. Similarly, the threshold scaling coefficient controls the balance between residual noise suppression and preservation of weak anatomical details. The selected bior4.4 configuration with α = 1.0 was therefore retained because it provided stable reconstruction behavior and throughput above 200 images/s in the evaluated Gaussian-noise setting.

The comparison with the reproduced data-parallel GPU implementation provides a same-environment reference against a recent wavelet-based GPU denoising approach. Direct numerical comparison with other published GPU wavelet denoising implementations is difficult because reported throughput values strongly depend on GPU architecture, image size, wavelet configuration, noise model, memory-transfer policy, and whether timing includes only kernel execution or full end-to-end processing. For this reason, the present study uses a reproduced data-parallel GPU baseline under identical dataset, hardware, and filter settings as the main comparative reference. This design provides a controlled context for interpreting the observed throughput gain.

These findings demonstrate that the proposed hybrid method constitutes an effective GPU-based solution for high-throughput wavelet MRI denoising. In comparison with the CPU implementation, it preserves reconstruction quality while substantially reducing runtime. In comparison with the reproduced GPU data-parallel baseline, it achieves lower latency together with stable quantitative reconstruction performance.

The position of the proposed method relative to CNN-based denoising models should be interpreted with caution. Learning-based denoisers, including residual CNNs, U-Net-like architectures, and CNN–Transformer hybrids, may achieve higher reconstruction metrics when trained and tested under matched degradation conditions, because they can learn data-driven image priors from representative noisy-clean or self-supervised training data. However, such methods address a different methodological setting from the one considered in the present study. Their performance depends on the availability and representativeness of training data, the correspondence between training and target acquisition conditions, network architecture, hyperparameter selection, training protocol, and hardware-specific inference optimization.

For this reason, a fully matched direct benchmark against state-of-the-art CNN-based denoisers was not included as a primary experiment in the present work. Such a comparison would require a separate and carefully controlled protocol, including paired or consistently simulated training data, train-validation-test separation, model training or fine-tuning, convergence monitoring, architecture-specific optimization, and inference-time measurement under the same hardware and preprocessing conditions. Without such a protocol, a direct numerical comparison could be misleading, because differences in PSNR, SSIM, or runtime could reflect training-data availability and model adaptation rather than the execution-level properties of the proposed CUDA implementation.

The proposed method therefore should not be interpreted as a replacement for learning-based MRI denoising models. Its main contribution is different: it provides a deterministic, interpretable, training-free, and high-throughput CUDA implementation of a fixed DWT–thresholding–inverse DWT denoising pipeline. This makes it suitable as a reproducible computational baseline, a transparent preprocessing stage, or a low-latency component in MRI image-processing workflows. CNN-based denoisers may be preferable when representative training data and sufficient model-validation resources are available, whereas the proposed method is advantageous when reproducibility, explicit reconstruction control, low latency, and independence from training data are required. A systematic experimental comparison with trained CNN-based and hybrid deep learning denoisers is therefore an important direction for future work.

Although the proposed slice-wise implementation performs stably, it does not explicitly model inter-slice noise correlations that may arise in 3D MRI because of shared acquisition coils, reconstruction pipelines, and spatial continuity between adjacent slices. Therefore, its current applicability should be interpreted primarily at the 2D slice-processing level. A direct slice-by-slice application of the proposed method to a volume would preserve the existing CUDA kernels and runtime behavior for individual slices, but it would not constitute true volumetric denoising because through-plane anatomical continuity and inter-slice noise dependencies would not be used.

A full volumetric extension would require replacing the separable 2D DWT with a separable 3D DWT. In this case, filtering would be performed along three spatial directions rather than two, and the first-level decomposition would produce eight subbands instead of four: LLL, LLH, LHL, LHH, HLL, HLH, HHL, and HHH. The LLL subband would serve as the volumetric approximation component, while the seven remaining subbands would represent detail components requiring thresholding. Thus, the 3D extension would change not only the number of subbands, but also the thresholding scheme, synchronization structure, and intermediate-buffer organization.

From the CUDA implementation perspective, the main architectural adaptations would include replacing 2D image buffers with 3D volume buffers, introducing depth-direction filtering kernels, extending symmetric boundary handling to the z direction, and redesigning stream scheduling for directional passes along the x, y, and z-axes. The fused reconstruction strategy would also need to be generalized from the final 2D horizontal inverse reconstruction to staged 3D inverse reconstruction, where depth-wise, vertical, and horizontal synthesis operations are coordinated to minimize intermediate global-memory traffic. In contrast to the current 2D pipeline, true 3D filtering introduces z-direction halo regions and inter-slice dependencies, which require memory-aware volume tiling and additional synchronization between directional passes.

A preliminary feasibility assessment indicates that memory organization is the central challenge for such an extension. For example, a 512 × 512 × 128 float32 MRI volume requires approximately 128 MiB for one volume buffer. Although the combined size of the eight first-level 3D DWT subbands is approximately equal to the input volume size, practical separable filtering and inverse reconstruction require additional temporary buffers, halo regions, and synchronization buffers. Therefore, a monolithic full-volume fused 3D kernel would be less suitable because it may increase shared-memory use, register pressure, and occupancy limitations. A more realistic implementation would use slab-wise or tiled processing, shared-memory staging for local 3D blocks, reuse of intermediate buffers, and asynchronous overlap between volume-block transfers and computation.

These considerations show that the proposed execution principles – branch-level parallelism, constant-memory filter coefficients, in-kernel boundary handling, and reduction of intermediate global-memory traffic – can be transferred conceptually to volumetric MRI denoising. However, the 3D extension is not a simple dimensional generalization of the present 2D pipeline. It represents a separate architectural problem involving volume tiling, depth-direction memory access, inter-slice dependencies, GPU occupancy, memory bandwidth, and synchronization control. For this reason, the development and profiling of a full 3D DWT–thresholding–inverse DWT CUDA implementation on volumetric MRI data is identified as a dedicated direction for future work.

The present study focuses on reproducible synthetic and simulated noise models to ensure controlled benchmarking of throughput and reconstruction behavior. This design makes it possible to compare CPU, reproduced data-parallel GPU, and proposed hybrid CUDA implementations under identical degradation conditions. However, these models cannot fully reproduce the variability of real clinical MRI acquisitions. In particular, real in vivo data may contain spatially varying noise caused by coil sensitivity profiles, B0 and B1 field inhomogeneities, parallel imaging reconstruction, scanner-specific preprocessing, gradient nonlinearities, and patient-dependent motion patterns. The Gaussian model captures only the additive noise component, the Rician model approximates signal-dependent magnitude MRI noise, and the image-domain motion model represents only a simplified form of motion-related degradation. Therefore, the reported results should be interpreted as a controlled methodological and computational validation rather than as a direct estimate of clinical denoising performance. Future work will include validation on in vivo MRI datasets acquired from multiple scanners and protocols, with explicit assessment of spatially varying artifacts and scanner-dependent noise patterns.

Another limitation concerns the dependence of the proposed parallel execution scheme on the target hardware configuration and on the structure of the degradation being processed. The reported acceleration was obtained on an NVIDIA Tesla T4 GPU using 512 × 512 2D MRI slices, a single-level DWT, fixed bior4.4 filters, and a fixed CUDA block configuration. Therefore, the same throughput and occupancy behavior should not be assumed for other GPU architectures without additional profiling. GPUs with different memory bandwidth, cache hierarchy, streaming multiprocessor count, register capacity, concurrent-kernel support, and host-device transfer behavior may change the relative benefit of dual-stream scheduling and fused reconstruction. For example, on GPUs with lower memory bandwidth, the vertical filtering and inverse reconstruction stages may become stronger bottlenecks, whereas on newer architectures with larger caches and higher concurrency the overlap between independent wavelet branches may be more effective.

The stability of the parallel scheme may also depend on the noise regime and acquisition characteristics. Spatially uniform Gaussian noise is well suited to the current fixed-thresholding and single-level DWT configuration, whereas spatially varying noise, coil-dependent noise amplification, strong Rician bias at low signal intensities, k-space-related motion artifacts, and scanner-specific reconstruction artifacts may change the distribution of wavelet coefficients across subbands. This can affect not only PSNR and SSIM, but also the workload balance between CUDA streams if some subbands require different thresholding or additional correction stages. Thus, the present results demonstrate stable behavior under the controlled degradation models and hardware configuration used in this study, but further validation is needed across different GPU architectures, MRI scanner vendors, field strengths, reconstruction pipelines, and noise regimes.

A further limitation concerns the simulation of motion-related artifacts. In this study, motion was approximated by rigid in-plane shifts and rotations applied directly in the image domain. This setting provides a reproducible way to test the runtime and reconstruction behavior of the proposed implementation under structured spatial degradation. However, real MRI motion artifacts are more complex and may involve k-space corruption, trajectory inconsistencies, ghosting patterns, phase errors, and patient-dependent motion during acquisition. Consequently, the motion-artifact results reported here should be interpreted as controlled proof-of-concept evidence rather than as a full validation under clinical motion conditions. Future work should include evaluation on real motion-corrupted MRI data or more realistic k-space-based motion simulations.

5 Conclusion

In this study, a hybrid CUDA-oriented wavelet method for denoising two-dimensional MRI images was developed. The proposed approach combines data-parallel and task-parallel computation, dual-stream scheduling, and fused processing stages within the DWT–thresholding–inverse DWT pipeline. It provides high computational efficiency, stable per-image processing time, and high throughput on the evaluated hardware platform.

The results confirmed that the primary target was achieved: on the Brain Tumor MRI dataset under Gaussian noise, the method delivers a throughput of 215–220 images per second while maintaining reconstruction quality at PSNR of about 29–30 and SSIM of about 0.90. Compared with the reproduced data-parallel GPU wavelet baseline implemented and tested under the same experimental conditions, the proposed solution reduces execution time by 20%–30% while preserving practically comparable reconstruction quality and demonstrating robustness under Rician noise and motion artifacts, as well as substantial acceleration relative to the CPU implementation.

The proposed method was positioned relative to learning-based denoising approaches as a deterministic and training-free alternative rather than as a direct replacement for CNN-based models. A fully matched benchmark against trained CNN-based denoisers requires a separate experimental protocol with controlled training data, model adaptation, and inference-time evaluation, and is therefore identified as an important direction for future work.

Thus, the proposed method provides an effective approach to high-throughput MRI image denoising, combining processing speed, execution stability, and reconstruction quality in a coordinated manner. Since the evaluation was performed using controlled synthetic and simulated degradations on 2D MRI slices without reference noise-free clinical acquisitions, the obtained PSNR, SSIM, and MSE values should be regarded as proof-of-concept methodological indicators rather than definitive evidence of clinical denoising efficacy. In particular, the motion-artifact experiments were based on simplified image-domain rigid transformations and therefore do not fully represent k-space-related ghosting, phase distortions, spatially varying artifacts, and complex patient-motion effects observed in clinical MRI. Consequently, future research should focus on validating the method on real in vivo MRI datasets acquired from multiple scanners, field strengths, acquisition protocols, and reconstruction pipelines, as well as on profiling the proposed parallel scheme across different GPU architectures and more complex noise regimes, including spatially varying noise, coil-dependent artifacts, and k-space-related motion degradation.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Generative AI was used only for proofreading purposes, including spelling correction, grammar checking, and minor language refinement. The authors reviewed and approved all revisions and take full responsibility for the final content of the manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Raouf Hamzaoui, De Montfort University, United Kingdom

Reviewed by: Thangavel Murugan, United Arab Emirates University, United Arab Emirates

Vidula Meshram, Vishwakarma Institute of Technology, India

Varsha Harpale, Pimpri Chinchwad College Of Engineering (PCCOE), India