This paper introduces a novel approach to quantitative contrast-enhanced spectral mammography (CESM) using triple energy K-edge imaging and a two-pass material decomposition method. The primary objective is to accurately determine local projected iodine concentration by leveraging the iodine K-edge in presence of glandular and adipose breast tissue. Addressing the ill-posed problem of established dual-energy methods, our approach accounts for three base materials with three non-redundant measurements. The proposed method involves two key steps. First, we estimate the local total thickness D via three material decomposition with the use of an iterative beam hardening correction, followed by strong denoising. Second, we use D as a constraint computing the recombined iodine image in classical dual-energy technique. In addition based on noise propagation the dose distribution and tube voltage for each image was optimized with respect to minimize standard deviation. Simulation studies using a numerical CIRS dual-energy phantom demonstrate thickness independent effective background cancellation, accurate reconstruction of relative projected iodine concentrations at moderate image noise levels. Sensitivity to remaining inaccuracies was simulated by adding scatter and deviation of used spectra for reconstruction and beam hardening correction. Furthermore, measurements were performed in laboratory system to test the feasibility. In conclusion, the proposed method provides a robust theoretical framework for accurately reconstructing projected iodine concentration images in CESM. For real measurements the method requires a highly accurate estimation of x-ray spectrum and knowledge of attenuation values.

Radiologists have preferred visual impressions or 'styles' of X-ray images that are manually adjusted to their needs to support their diagnostic performance. In this work, we propose an automatic and interpretable X-ray style transfer by introducing a trainable version of the Local Laplacian Filter (LLF) [1]. From the shape of the LLF's optimized remap function, the characteristics of the style transfer can be inferred and reliability of the algorithm can be ensured. Moreover, we enable the LLF to capture complex X-ray style features by replacing the remap function with a Multi-Layer Perceptron (MLP) and adding a trainable normalization layer. We demonstrate the effectiveness of the proposed method by transforming unprocessed mammographic X-ray images into images that match the style of target mammograms and achieve a Structural Similarity Index (SSIM) of 0.94 compared to 0.82 of the baseline LLF style transfer method from [2].

Deep learning-based image analysis offers great potential in clinical practice. However, it faces mainly two challenges: scarcity of large-scale annotated clinical data for training and susceptibility to adversarial data in inference. As an example, an artificial intelligence (AI) system could check patient positioning, by segmenting and evaluating relative positions of anatomical structures in medical images. Nevertheless, data to train such AI system might be highly imbalanced with mostly well-positioned images being available. Thus, we propose the use of synthetic X-ray images and annotation masks forward projected from 3D photon-counting CT volumes to create realistic non-optimally positioned X-ray images for training. An open-source model (TotalSegmentator) was used to annotate the clavicles in 3D CT volumes. We evaluated model robustness with respect to the internal (simulated) patient rotation α on real-data-trained models and real&synthetic-data-trained models. Our results showed that real&synthetic- data-trained models have Dice score percentage improvements of 3% to 15% across different α groups compared to the real-data-trained model. Therefore, we demonstrated that synthetic data could be supplementary used to train and enrich heavily underrepresented conditions to increase model robustness.

The existence of metallic implants in projection images for cone-beam computed tomography (CBCT) introduces undesired artifacts which degrade the quality of reconstructed images. In order to reduce metal artifacts, projection in-painting is an essential step in many metal artifact reduction algorithms. In this work, a hybrid network combining the shift window (Swin) vision transformer (ViT) and a convolutional neural network is proposed as a baseline network for the inpainting task. To incorporate metal information for the Swin ViT-based encoder, metal-conscious self-embedding and neighborhood-embedding methods are investigated. Both methods have improved the performance of the baseline network. Furthermore, by choosing appropriate window size, the model with neighborhood-embedding could achieve the lowest mean absolute error of 0.079 in metal regions and the highest peak signal-to-noise ratio of 42.346 in CBCT projections. At the end, the efficiency of metal-conscious embedding on both simulated and real cadaver CBCT data has been demonstrated, where the inpainting capability of the baseline network has been enhanced.

Purpose:Digital breast tomosynthesis (DBT) has been introduced more than a decade ago. Studies have shown higher breast cancer detection rates and lower recall rates, and it has become an established imaging method in diagnostic settings. However, full-field digital mammography (FFDM) remains the most common imaging modality for screening in many countries, as it delivers high-resolution planar images of the breast. To combine the advantages of DBT with the faster acquisition and the unique in-plane resolution capabilities known from FFDM, a system concept was developed for application in screening and diagnosis. Approach:The concept comprises an X-ray tube with adaptive focal spot position based on the flying focal spot (FFS) technology and optimized X-ray spectra. This is combined with innovative algorithmic concepts for tomosynthesis reconstruction and synthetic mammograms (SMs). Results:An X-ray tube with FFS was incorporated into a DBT system that performs 50-deg wide tomosynthesis scans with 25 projections in 4.85 s. Laboratory evaluations demonstrated significant improvements in the effective modular transfer function (eMTF). The improved eMTF as well as the effectiveness of the algorithmic concepts is shown in images from a clinical evaluation study. Conclusions:The DBT system concept enables high spatial resolution at short acquisition times. This leads to improved microcalcification visibility, reduced risk of motion artifacts, and shorter breast compression times. It shifts the in-plane resolution of DBT into the high-resolution range of FFDM. The presented technology leap might be a key contributor to facilitating the paradigm shift of replacing FFDM with DBT plus SM.

Wide-angle digital breast tomosynthesis (DBT) is well known to offer benefits in mass perceptibility compared to narrow-angle DBT due to reduced anatomical overlap. Regarding the perceptibility of micro-calcifications the situation is somehow inverted. On the one hand this can be related to effects during data acquisition and their impact on the system MTF. On the other hand there is a wider spread of calcifications in depth direction in narrow-angle DBT, which distributes calcifications over more slices. This is equivalent to an inherent thicker slice for high spatial frequencies. In this work we want to assume an equivalent quality of raw data and only focus on the effects of different acquisition angles in the reconstruction. We propose an algorithm which creates so-called hybrid thick DBT slices and optimizes the visualization of calcifications while preserving the high mass perceptibility of thin wide-angle DBT slices. The algorithm is purely based on filtered backprojection (FBP) and can be implemented in an efficient manner. For validation simulation studies using the VICTRE (FDA) pipeline are performed. Our results indicate that hybrid thick-slices in wide-angle DBT enable to successfully solve the contrarian imaging tasks of high mass and high calcification perception within one imaging setup.

Denoising algorithms are sensitive to the noise level and noise power spectrum of the input image and their ability to adapt to this. In the worst-case, image structures can be accidentally removed or even added. This holds up for analytical image filters but even more for deep learning-based denoising algorithms due to their high parameter space and their data-driven nature. We propose to use the knowledge about the noise distribution of the image at hand to limit the influence and ability of denoising algorithms to a known and plausible range. Specifically, we can use the physical knowledge of X-ray radiography by considering the Poisson noise distribution and the noise power spectrum of the detector. Through this approach, we can limit the change of the acquired signal by the denoising algorithm to the expected noise range, and therefore prevent the removal or hallucination of small relevant structures. The presented method allows to use denoising algorithms and especially deep learning-based methods in a controlled and safe fashion in medical x-ray imaging.

The existence of metallic implants in projection images for cone-beam computed tomography (CBCT) introduces undesired artifacts which degrade the quality of reconstructed images. In order to reduce metal artifacts, projection inpainting is an essential step in many metal artifact reduction algorithms. In this work, a hybrid network combining the shift window (Swin) vision transformer (ViT) and a convolutional neural network is proposed as a baseline network for the inpainting task. To incorporate metal information for the Swin ViT-based encoder, metal-conscious self-embedding and neighborhoodembedding methods are investigated [1]. Both methods have improved the performance of the baseline network. Furthermore, by choosing appropriate window size, the model with neighborhood-embedding could achieve the lowest mean absolute error of 0.079 in metal regions and the highest peak signal-to-noise ratio of 42.346 in CBCT projections. At the end, the efficiency of metal-conscious embedding on both simulated and real cadaver CBCT data has been demonstrated, where the inpainting capability of the baseline network has been enhanced.

Digital breast tomosynthesis ( DBT) enables significantly higher cancer detection rates compared to full-field digital mammography (FFDM) without compromising the recall rate. However, regarding microcalcification assessment established tomosynthesis system concepts still tend to be inferior to FFDM. To further boost the clinical role of DBT in breast cancer screening and diagnosis, a system concept was developed that enables fast wide-angle DBT with the unique in-plane resolution capabilities known from FFDM. The concept comprises a novel X-ray tube concept that incorporates an adaptive focal spot position, fast flat-panel detector technology, and innovative algorithmic concepts for image reconstruction. We have built a DBT system that provides tomosynthesis image stacks and synthetic mammograms from 50 degrees tomosynthesis scans realized in less than five seconds. In this contribution, we motivate the design of the system concept, present a physics characterization of its imaging performance, and outline the algorithmic concepts used for image processing. We conclude with illustrating the potential clinical impact by means of clinical case examples from first evaluations in Europe.

BACKGROUND:Due to the high attenuation of metals, severe artifacts occur in cone beam computed tomography (CBCT). The metal segmentation in CBCT projections usually serves as a prerequisite for metal artifact reduction (MAR) algorithms. PURPOSE:The occurrence of truncation caused by the limited detector size leads to the incomplete acquisition of metal masks from the threshold-based method in CBCT volume. Therefore, segmenting metal directly in CBCT projections is pursued in this work. METHODS:Since the generation of high quality clinical training data is a constant challenge, this study proposes to generate simulated digital radiographs (data I) based on real CT data combined with self-designed computer aided design (CAD) implants. In addition to the simulated projections generated from 3D volumes, 2D x-ray images combined with projections of implants serve as the complementary data set (data II) to improve the network performance. In this work, SwinConvUNet consisting of shift window (Swin) vision transformers (ViTs) with patch merging as encoder is proposed for metal segmentation. RESULTS:The model's performance is evaluated on accurately labeled test datasets obtained from cadaver scans as well as the unlabeled clinical projections. When trained on the data I only, the convolutional neural network (CNN) encoder-based networks UNet and TransUNet achieve only limited performance on the cadaver test data, with an average dice score of 0.821 and 0.850. After using both data II and data I during training, the average dice scores for the two models increase to 0.906 and 0.919, respectively. By replacing the CNN encoder with Swin transformer, the proposed SwinConvUNet reaches an average dice score of 0.933 for cadaver projections when only trained on the data I. Furthermore, SwinConvUNet has the largest average dice score of 0.953 for cadaver projections when trained on the combined data set. CONCLUSIONS:Our experiments quantitatively demonstrate the effectiveness of the combination of the projections simulated under two pathways for network training. Besides, the proposed SwinConvUNet trained on the simulated projections performs state-of-the-art, robust metal segmentation as demonstrated on experiments on cadaver and clinical data sets. With the accurate segmentations from the proposed model, MAR can be conducted even for highly truncated CBCT scans.

In this work, we propose a U-Net-based super-resolution neural network, SRU-Net, to create emulated high spatial resolution (eHR) CT images from low spatial resolution (LR) CT images. As resolution could be defined by the modulation transfer function in CT reconstruction, we propose the novel approach based on CT reconstruction kernels to create realistic multi-detector CT (MDCT) synthetic LR images from high-resolution cone-beam CT (CBCT) scans. Keeping a constant sampling grid size of 0.20 × 0.20mm2, we reconstruct two types of MDCT-like LR images and one corresponding HR image from the same CBCT raw data and train two models respectively. We validated the performance of the trained models on unseen LR CBCT images. We then applied the trained network to MDCT images. Mean squared error, structural similarity index measures and peak signal-to-noise ratio of two models show significant improvements (p < 0.001) in the eHR images.

In several image acquisition and processing steps of X-ray radiography, knowledge of the existence of metal implants and their exact position is highly beneficial (e.g. dose regulation, image contrast adjustment). Another application which would benefit from an accurate metal segmentation is cone beam computed tomography (CBCT) which is based on 2D X-ray projections. Due to the high attenuation of metals, severe artifacts occur in the 3D X-ray acquisitions. The metal segmentation in CBCT projections usually serves as a prerequisite for metal artifact avoidance and reduction algorithms. Since the generation of high quality clinical training is a constant challenge, this study proposes to generate simulated X-ray images based on CT data sets combined with self-designed computer aided design (CAD) implants and make use of convolutional neural network (CNN) and vision transformer (ViT) for metal segmentation. Model test is performed on accurately labeled X-ray test datasets obtained from specimen scans. The CNN encoder-based network like U-Net has limited performance on cadaver test data with an average dice score below 0.30, while the metal segmentation transformer with dual decoder (MST-DD) shows high robustness and generalization on the segmentation task, with an average dice score of 0.90. Our study indicates that the CAD model-based data generation has high flexibility and could be a way to overcome the problem of shortage in clinical data sampling and labelling. Furthermore, the MST-DD approach generates a more reliable neural network in case of training on simulated data.

Purpose This paper studies spatial resolution that is achievable with a fast slot‐scanning tomosynthesis approach for orthopedic examinations. Hereby, we use parallel scanning motion implemented in a twin robotic x‐ray system. Methods We have measured and analyzed the modulation transfer function (MTF) for various combinations of scanning speed, x‐ray tube voltage, pulse length, the nominal focal spot size as well as source‐to‐object distances. Moreover, we present a theoretical model which describes the system in terms of the MTF. The system was equipped with newly developed linear trajectory prototypes for slot scanning. The acquired images form the basis for a small‐angle tomosynthesis reconstruction. In total, three different scanning speeds (27, 14, 8 cm/s), pulse lengths (1, 2, 4 ms), tube voltages (80, 100, 120 kV), two nominal focal spot sizes (0.6, 1.0), and three source‐to‐object distances (950, 1050, 1150 mm) were investigated. To determine the resolution capabilities, we measured the MTF for the given parameter space. The results were then used to design a filter that yields a desired resolution in the reconstructed image. In addition, we also measured the noise power spectrum (NPS) to show the influence of the aforementioned filters on the noise distribution. Results We have shown that the presented model is in good agreement with the performed measurements. Scanning speed and pulse width have an impact on the MTF in the scanning direction. Up to a travel distance of 0.3 mm during an x‐ray pulse, an isotropic resolution can be achieved. Longer pulse width or higher scanning speed cause anisotropic resolution. Moreover, it is shown that none of the investigated parameters have an influence on the MTF perpendicular to the scanning direction (slot direction). The 10% MTF value ranges between 9 and 18 lp/cm in the scanning direction and about 18 lp/cm in slot direction. Tube voltage, nominal focal spot size as well as the source‐to‐object distance showed no major impact on the system MTF. In terms of the anisotropic resolution capabilities, we have shown that limiting the resolution in the slot direction to obtain isotropic resolution is possible yet at the cost of an inhomogeneous noise pattern. However, maintaining the resolution in slot direction will provide a better edge response and a more homogeneous noise texture at the cost of an inhomogeneous image resolution. Conclusions We have demonstrated the feasibility of the slot‐scanning technique using a twin‐robotic x‐ray system. Even the fastest scanning mode (27 cm/s) yields image resolution on a level that is sufficient for typical orthopedic examinations in terms of musculoskeletal (MSK) measurements. Moreover, it could be shown that the application of specifically designed target MTFs on two‐dimensional x‐ray images is feasible.

The estimation of patient dose using Monte Carlo (MC) simulations based on the available patient CT images is limited to the length of the scan. Software tools for dose estimation based on standard computational phantoms overcome this problem; however, they are limited with respect to taking individual patient anatomy into account. The purpose of this study was to generate whole-body patient models in order to take scattered radiation and over-scanning effects into account. Thorax examinations were performed on three physical anthropomorphic phantoms at tube voltages of 80 kV and 120 kV; absorbed dose was measured using thermoluminescence dosimeters (TLD). Whole-body voxel models were built as a combination of the acquired CT images appended by data taken from widely used anthropomorphic voxel phantoms.MC simulations were performed both for the CT image volumes alone and for the whole-body models. Measured and calculated dose distributions were compared for each TLD chip position; additionally, organ doses were determined. MC simulations based only on CT data underestimated dose by 8% -15% on average depending on patient size with highest underestimation values of 37% for the adult phantom at the caudal border of the image volume. The use of whole-body models substantially reduced these errors; measured and simulated results consistently agreed to better than 10%.This study demonstrates that combined whole-body models can provide three-dimensional dose distributions with improved accuracy. Using the presented concept should be of high interest for research studies which demand high accuracy, e. g. for dose optimization efforts. (C) 2014 Associazione Italiana di Fisica Medica. Published by Elsevier Ltd. All rights reserved.

Iterative reconstruction (IR) methods have recently re-emerged in transmission x-ray computed tomography (CT). They were successfully used in the early years of CT, but given up when the amount of measured data increased because of the higher computational demands of IR compared to analytical methods. The availability of large computational capacities in normal workstations and the ongoing efforts towards lower doses in CT have changed the situation; IR has become a hot topic for all major vendors of clinical CT systems in the past 5 years. This review strives to provide information on IR methods and aims at interested physicists and physicians already active in the field of CT. We give an overview on the terminology used and an introduction to the most important algorithmic concepts including references for further reading. As a practical example, details on a model-based iterative reconstruction algorithm implemented on a modern graphics adapter (GPU) are presented, followed by application examples for several dedicated CT scanners in order to demonstrate the performance and potential of iterative reconstruction methods. Finally, some general thoughts regarding the advantages and disadvantages of IR methods as well as open points for research in this field are discussed.

PURPOSE:Monte Carlo (MC) simulation is an established technique for dose calculation in diagnostic radiology. The major drawback is its high computational demand, which limits the possibility of usage in real-time applications. The aim of this study was to develop fast on-site computed tomography (CT) specific MC dose calculations by using a graphics processing unit (GPU) cluster.METHODS:GPUs are powerful systems which are especially suited to problems that can be expressed as data-parallel computations. In MC simulations, each photon track is independent of the others; each launched photon can be mapped to one thread on the GPU, thousands of threads are executed in parallel in order to achieve high performance. For further acceleration, the authors considered multiple GPUs. The total computation was divided into different parts which can be calculated in parallel on multiple devices. The GPU cluster is an MC calculation server which is connected to the CT scanner and computes 3D dose distributions on-site immediately after image reconstruction. To estimate the performance gain, the authors benchmarked dose calculation times on a 2.6 GHz Intel Xeon 5430 Quad core workstation equipped with two NVIDIA GeForce GTX 285 cards. The on-site calculation concept was demonstrated for clinical and preclinical datasets on CT scanners (multislice CT, flat-detector CT, and micro-CT) with varying geometry, spectra, and filtration. To validate the GPU-based MC algorithm, the authors measured dose values on a 64-slice CT system using calibrated ionization chambers and thermoluminesence dosimeters (TLDs) which were placed inside standard cylindrical polymethyl methacrylate (PMMA) phantoms.RESULTS:The dose values and profiles obtained by GPU-based MC simulations were in the expected good agreement with computed tomography dose index (CTDI) measurements and reference TLD profiles with differences being less than 5%. For 10(9) photon histories simulated in a 256 × 256 × 12 voxel thorax dataset with voxel size of 1.36 × 1.36 × 3.00 mm(3), calculation times of about 70 and 24 min were necessary with single-core and multiple-core central processing unit (CPU) solutions, respectively. Using GPUs, the same MC calculations were performed in 1.27 min (single card) and 0.65 min (two cards) without a loss in quality. Simulations were thus speeded up by factors up to 55 and 36 compared to single-core and multiple-core CPU, respectively. The performance scaled nearly linearly with the number of GPUs. Tests confirmed that the proposed GPU-based MC tool can be easily adapted to different types of CT scanners and used as service providers for fast on-site dose calculations.CONCLUSIONS:The Monte Carlo software package provides fast on-site calculation of 3D dose distributions in the CT suite which makes it a practical tool for any type of CT-specific application.

OBJECTIVE:Mammography, today's standard imaging approach, has deficits with respect to the superimposition of anatomical structures. Dedicated CT of the breast so far indicated that it can provide superior soft-tissue imaging, but that it still has significant limitations with respect to spatial resolution and dose. We have assessed novel dedicated breast CT technology.METHODS:Based on simulations and measurements we developed novel technology which uses direct-conversion CdTe material and photon-counting electronics with 100 μm detector element size for close to 100% dose efficiency. We assessed the potential for the imaging of microcalcifications of 100 to 200 μm diameter and soft-tissue lesions of 1 to 5 mm diameter by simulations at dose levels between 1 and 6 mGy.RESULTS:Microcalcifications of 150 μm and soft-tissue lesions of 2 mm diameter were found to be clearly detectable at an average glandular dose of 3 mGy. Separate displays are required for high-resolution microcalcification and for low-resolution soft-tissue analysis. Total CT data acquisition time will be below 10 s.CONCLUSION:Dedicated breast CT may eventually provide comprehensive diagnostic assessment of microcalcifications and soft-tissue structures at dose levels equivalent to or below those of two-view screening mammography.

Metallic implants are responsible for various artifacts in flat-detector computed tomography visible as streaks and dark areas in the reconstructed volumetric images. In this paper a novel method for a fast reduction of these metal artifacts is presented using a three-step correction procedure to approximate the missing parts of the raw data. In addition to image quality aspects, this paper deals with the problem of high correction latencies by proposing a reconstruction and correction framework, that utilizes the massive computational power of graphics processing units (GPUs). An initial volume is reconstructed, followed by a 3-dimensional metal voxel segmentation algorithm. These metal voxels allow us to identify metal-influenced detector elements by using a simplified geometric forward projection. Consequently, these areas are corrected using a 3D interpolation scheme in the raw data domain, followed by a second reconstruction. This volume is then segmented into three materials with respect to bone structures using a threshold-based algorithm. A forward projection of the obtained tissueclass model substitutes missing or corrupted attenuation values for each detector element affected by metal and is followed by a final reconstruction. The entire process including the initial reconstruction, takes less than a minute (5123 volume with 496 projections of size 1240x960) and offers significant improvements of image quality. The method was evaluated with data from two FD-CT C-arm systems (Artis Zee and Artis Zeego, Siemens Healthcare, Forchheim, Germany).

Metallic implants generate streak-like artifacts in flat-detector computed tomography (FD-CT) reconstructed volumetric images. This study presents a novel method for reducing these disturbing artifacts by inserting discarded information into the original rawdata using a three-step correction procedure and working directly with each detector element. Computation times are minimized by completely implementing the correction process on graphics processing units (GPUs). First, the original volume is corrected using a three-dimensional interpolation scheme in the rawdata domain, followed by a second reconstruction. This metal artifact-reduced volume is then segmented into three materials, i.e. air, soft-tissue and bone, using a threshold-based algorithm. Subsequently, a forward projection of the obtained tissue-class model substitutes the missing or corrupted attenuation values directly for each flat detector element that contains attenuation values corresponding to metal parts, followed by a final reconstruction. Experiments using tissue-equivalent phantoms showed a significant reduction of metal artifacts (deviations of CT values after correction compared to measurements without metallic inserts reduced typically to below 20 HU, differences in image noise to below 5 HU) caused by the implants and no significant resolution losses even in areas close to the inserts. To cover a variety of different cases, cadaver measurements and clinical images in the knee, head and spine region were used to investigate the effectiveness and applicability of our method. A comparison to a three-dimensional interpolation correction showed that the new approach outperformed interpolation schemes. Correction times are minimized, and initial and corrected images are made available at almost the same time (12.7 s for the initial reconstruction, 46.2 s for the final corrected image compared to 114.1 s and 355.1 s on central processing units (CPUs)).

Median filtering is a commonly used technique in smoothing and denoising applications. Based on the vector programming model of modern commodity graphics processing units (GPUs), which directly support for minmax operations, compare and select as fundamental instructions, we implemented the branchless vectorized median (BVM) filter proposed in reference [1] using NVIDIA's compute unified device architecture (CUDA). The BVM filter keeps track of a sorted array from which values are deleted and to which new values are inserted. Although it is of O(M2) computational complexity while other sort algorithms are of O(M ln M) computational complexity, at least for typical data, it may outperform other implementations. The mainly reason is that this algorithm is branchless, and it makes use of data-level parallelism thereby its runtime is data- independent and highly predictable. We describe some important criteria such as the memory layout for a fast accessing scheme and discuss the bottlenecks in the branchless vectorized median computation. We provide performance benchmarks in comparsion to other implementations, a median filter on GPUs based on comparing a pivot value to all values, and the same branchless vectorized median implementation on CPUs. The comparison uses constant data, linear data, and random data. The runtime of BVM is independent of the data and shows a factor up to 4.6 faster than the pivot median filter. Although the performance of CUDA-based BVM filter is roughly 25% slower compared to a CPU-based (8 cores) routine, it is still a cheaper solution for many applications. We also present some factors such as the array size, the number of arrays and the filter size that influence the total BVM performance. An application of median filter for ring artifacts reduction will be demonstrated. The processing time is up to 3.7 times faster than the optimized CPU-based (four cores) routine.