Applied Science, Faculty of
Non UBC
Electrical and Computer Engineering, Department of
DSpace
Computation 11 (1): 7 (2023)
Haider, Hassaan
Shah, Jawad Ali
Kadir, Kushsairy
Khan, Najeeb
In the Compressed Sensing (CS) framework, the underdetermined system of linear equation (USLE) can have infinitely many possible solutions. However, we intend to find the sparsest possible solution, which is l₀-norm minimization. However, finding an l₀ norm solution out of infinitely many possible solutions is NP-hard problem that becomes non-convex optimization problem. It has been a practically proven fact that l₀ norm penalty can be adequately estimated by l₁ norm, which recasts a non-convex minimization problem to a convex problem. However, l₁ norm non-differentiable and gradient-based minimization algorithms are not applicable, due to this very reason there is a need to approximate l₁ norm by its smooth approximation. Iterative shrinkage algorithms provide an efficient method to numerically minimize l₁-regularized least square optimization problem. These algorithms are required to induce sparsity in their solutions to meet the CS recovery requirement. In this research article, we have developed a novel recovery method that uses hyperbolic tangent function to recover undersampled signal/images in CS framework. In our work, l₁ norm and soft thresholding are both approximated with the hyperbolic tangent functions. We have also proposed the criteria to tune optimization parameters to get optimal results. The error bounds for the proposed l₁ norm approximation are evaluated. To evaluate performance of our proposed method, we have utilized a dataset comprised of 1-D sparse signal, compressively sampled MR image and cardiac cine MRI. The MRI is an important imaging modality for assessing cardiac vascular function. It provides the ejection fraction and cardiac output of the heart. However, this advantage comes at the cost of a slow acquisition process. Hence, it is essential to speed up the acquisition process to take the full benefits of cardiac cine MRI. Numerical results based on performance metrics, such as Structural Similarity (SSIM), Peak Signal to Noise Ratio (PSNR) and Root Mean Square Error (RMSE) show that the proposed tangent hyperbolic based CS recovery offers a much better performance as compared to the traditional Iterative Soft Thresholding (IST) recovery methods.
Citation: Haider, H.; Shah, J.A.;Kadir, K.; Khan, N. SparseReconstruction Using HyperbolicTangent as Smooth l1-NormApproximation. Computation 2023, 11,7. https://doi.org/10.3390/computation11010007Academic Editors: Anando Sen andDemos T. TsahalisReceived: 11 October 2022Revised: 29 November 2022Accepted: 16 December 2022Published: 4 January 2023Copyright: © 2023 by the authors.Licensee MDPI, Basel, Switzerland.This article is an open access articledistributed under the terms andconditions of the Creative CommonsAttribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).computationArticleSparse Reconstruction Using Hyperbolic Tangent as Smoothl1-Norm ApproximationHassaan Haider 1 , Jawad Ali Shah 1,*, Kushsairy Kadir 2 and Najeeb Khan 31 Department of Electrical and Computer Engineering, Faculty of Engineering and Technology,International Islamic University, Islamabad 44000, Pakistan2 Electrical Section, British Malaysian Institute, Universiti Kuala Lumpur, Gombak 53100, Malaysia3 Department of Electrical and Computer Engineering, University of British Columbia,Vancouver, BC V6T 1Z4, Canada* Correspondence: jawad.shah@iiu.edu.pk; Tel.: +92-333-5274448Abstract: In the Compressed Sensing (CS) framework, the underdetermined system of linear equation(USLE) can have infinitely many possible solutions. However, we intend to find the sparsest possiblesolution, which is l0-norm minimization. However, finding an l0 norm solution out of infinitely manypossible solutions is NP-hard problem that becomes non-convex optimization problem. It has been apractically proven fact that l0 norm penalty can be adequately estimated by l1 norm, which recasts anon-convex minimization problem to a convex problem. However, l1 norm non-differentiable andgradient-based minimization algorithms are not applicable, due to this very reason there is a needto approximate l1 norm by its smooth approximation. Iterative shrinkage algorithms provide anefficient method to numerically minimize l1-regularized least square optimization problem. Thesealgorithms are required to induce sparsity in their solutions to meet the CS recovery requirement.In this research article, we have developed a novel recovery method that uses hyperbolic tangentfunction to recover undersampled signal/images in CS framework. In our work, l1 norm and softthresholding are both approximated with the hyperbolic tangent functions. We have also proposedthe criteria to tune optimization parameters to get optimal results. The error bounds for the proposedl1 norm approximation are evaluated. To evaluate performance of our proposed method, we haveutilized a dataset comprised of 1-D sparse signal, compressively sampled MR image and cardiac cineMRI. The MRI is an important imaging modality for assessing cardiac vascular function. It providesthe ejection fraction and cardiac output of the heart. However, this advantage comes at the cost of aslow acquisition process. Hence, it is essential to speed up the acquisition process to take the fullbenefits of cardiac cine MRI. Numerical results based on performance metrics, such as StructuralSimilarity (SSIM), Peak Signal to Noise Ratio (PSNR) and Root Mean Square Error (RMSE) show thatthe proposed tangent hyperbolic based CS recovery offers a much better performance as compared tothe traditional Iterative Soft Thresholding (IST) recovery methods.Keywords: compressed sensing; holography cardiac cine MRI; l1-norm smooth approximations;hyperbolic tangent function; soft thresholding1. IntroductionCompressed Sensing (CS) exploits the sparsity of signals in a certain domain tofind a near-optimal solution to the underdetermined system of linear equations. In CS,the sampling of signals depends on the information rate rather than its bandwidth. CStechnique facilitates simultaneous acquisition and compression of compressible or sparsesignals that potentially reduce the acquisition time. The CS is a data acquisition methodthat allows for the reconstruction of a signal from very few measurements if the signal istransformed in a sparsifying domain, and these measurements are highly incoherent withrespect to its sparsifying transform. Unfortunately, most of the reconstruction techniquesof compressively sampled signals are computationally expensive and non-linear [1–3].Computation 2023, 11, 7. https://doi.org/10.3390/computation11010007 https://www.mdpi.com/journal/computationComputation 2023, 11, 7 2 of 21CS has recently been used to reconstruct under-sampled biomedical images by ex-ploiting the sparsity of biomedical images in the sparsifying domain. The Fourier-encodednature of the MR image scanning process and the existence of suitable sparsifying trans-form domains, i.e., Wavelets, Contourlets, total variation, etc., make the MRI a potentiallysuitable application of CS [4]. Incoherent sampling, which is another important require-ment of CS, can be accomplished with the variable density k-space sampling method tointroduce noise-like random aliasing artefacts during the MR image recovery. Variabledensity k-space under-sampling pattern samples with high density from the center of thek-space that contains maximum energy of the MR images and undersamples the outerk-space region with lower density to efficiently reduce the MR image scanning time [5].Cardiac cine magnetic resonance imaging (MRI) is an emerging medical imagingmodality to evaluate the growth of Cardiac-Vascular Disease (CVD). It is useful in evalu-ating the cardiac wall thickness and motion in CVD patients [6–8]. Further, cardiac cineMRI aids in performing the quantitative study of ejection fraction and cardiac output of theheart. The ejection fraction is the percentage of blood that is ejected out of the ventricleswith each contraction. This amount is used to determine heart failures and other typesof heart diseases [9]. Cardiac output measures the amount of blood pumped by the heartper minute. However, these advantages are limited by the lengthy acquisition process ofcardiac cine MRI that requires multiple breaths-holds of the patient and extended patientengagement in MRI scanner. Therefore, it is essential to accelerate the image acquisition incardiac cine MRI by using fast pulse sequences and/or by reducing the number of samplestaken during data acquisition [10,11]. As the former approach is inherently limited bydifferent constraints, much research interest is moved to the latter approach. CS can beapplied successfully to the cardiac cine MRI, where sparsity is exploited in the temporaldimension [4,12]. However, improving speed and efficiency of CS recovery methods isan active area of interest for researchers working in medical imaging especially MRI. Thekey conditions for the CS framework to work are sparsity, non-linear reconstruction andincoherent undersampling.In MR imaging, sparsity can be accomplished by transforming the image in its sparserepresentation. To fulfil the condition of incoherent sampling in MR imaging, variousundersampling patterns can be utilized, such as radial lines sampling and variable densitysampling [4]. Non-linear reconstruction numerical techniques involve l1-norm regulariza-tion in order to find sparse solutions to the least-squares optimization problem. l2-normbased regularization provides the linear and simplest solution to under-determined systemproblem. However, it minimizes the energy of the error and distributes it over all solutionset that results in non-sparse solution that does not fits in CS framework. Similarly lp-norm(1 < p < ∞) based regularization, as value of p starts growing it tends to penalize only thelargest parameter, such as max function, and some bad parameters may hide under thelargest parameters, which results in less-sparse solution. For this very reason l1-norm is thepreferred regularization, as it promotes sparsity that perfectly fits in CS framework [1–3].However, l1-norm penalty is non-differentiable, so applying efficient optimization methodsthat involve derivative are not feasible. Therefore, various methods have been proposedto resolve the l1-norm regularization problem. The IST based recovery methods havesuccessfully been utilized to efficiently reconstruct images from under-sampled data in theCS framework [5,13]. An iterative hard thresholding-based recovery method is proposedfor the compressed sensing problem [14]. However, this algorithm has limited performanceas compared to the soft thresholding-based methods. Random filters for compressivesampling [15], Bregman iterative algorithms for compressed sensing [16], and a weighted l1minimization recovery algorithm [17] are proposed to solve the compressed sensing recov-ery problem. Lately, smooth l1-norm penalty based sparse signal reconstruction methodwas evolved for approximation of l1-norm that uses a hyperbolic tangent function [18]. Theresearch shows that this technique can be used for the reconstruction of undersampled MRimages from fewer acquired samples, which allows fast imaging without compromisingspatial resolution. Jawad et al. proposed that the wavelet thresholding can be implementedComputation 2023, 11, 7 3 of 21using hyperbolic tangent function. It was explained that the differentiable hyperbolicfunction provides a much more accurate recovery than IST techniques [19–21]. It wasexperimentally shown that the hyperbolic tangent function performs a much improvedsignal recovery as compared to the hard threshold, soft threshold and Garrote thresholdfunctions [22,23].In this paper, CS technique is applied to reduce MRI scanning time, CS exploits thesparsity of MRI in Fourier domain that enables us to take fewer samples without compro-mising on the quality of image recovered from undersampled MRI. We propose a noveland more efficient CS recovery algorithm based on the hyperbolic tangent function forapproximating l1-norm and shrinkage operation for accurate recovery of compressivelysampled sparse signals, MR images and the cardiac cine MRI. We introduce smooth ap-proximation of l1-norm, the hyperbolic tangent function, where steepest descent algorithmis applicable for minimization of objective function. The error bounds for the proposedl1-norm penalty are presented in this paper. We have used the soft thresholding techniquebased on the hyperbolic tangent function that is inspired by the maximum a posteriori(MAP) noise estimator. In this work, we also recommend the efficient criteria for the tuningparameters. Performance analysis of the proposed method is shown using simulations; torecover random 1-D sparse signal, 2-D MR image and clinical cardiac cine MRI. Severalquantitative performance measures are used apart from qualitative depiction, i.e., MeanSquare Error (MSE), Root Mean Squared Error (RMSE), Signal to Noise Ratio (SNR), PeakSignal to Noise Ratio (PSNR), Improved Signal to Noise Ratio (ISNR), correlation, fitness,and Structural Similarity (SSIM) in order to prove the supremacy of proposed method overexisting recovery techniques.2. Materials and MethodsReconstruction of undersampled signal through CS is an optimization problem, whichpromotes sparsity in our solution by minimizing the l1–norm.2.1. Proposed MethodLet zeRn be the signal in a vector form and yeCm be the undersampled measurements.Then, the CS recovery function is written as:f (z) =12‖y−ΦΨHz‖22 + λ‖z‖1 (1)whereΦ is the sampling domain of the signal z, whereasΨ represents sparsifying transform.The tuning parameter λ in Equation (1) provides an important trade-off parameter betweenfidelity and sparsity. The performance of our algorithm is dependent on proper thresholdlevel selection. We have employed the fixed value expression, depending on the signaldimensions and its noise variance [13].λ = σv√2 ln(n) (2)where σv is the noise standard variance and n is the length of the sparse signal.Since the hyperbolic tangent function has properties, such as non-convex, odd, smoothanalytical bounded function that is monotonically increasing, the slope of the function atthe origin can be tuned to any desired value [21]. So, our proposed approximation for thel1 norm in Equation (1) is defined as:‖z‖1 ∼=n∑i=1zitanh(γzi) (3)Since the hyperbolic tangent function is used as a smooth and differentiable approxi-mation to l1-norm. Therefore, the value of γ is taken quite high to make it closer to l1-norm,as shown in Figure 1. It is also providing the benefits of the smoothness and differentiability.Equation (1) can now be written as:Computation 2023, 11, 7 4 of 21f (z) =12‖y−ΦΨHz‖22 + λ∑ni=1 zitanh(γzi) (4)Computation 2023, 10, x FOR PEER REVIEW 5 of 21 Figure 1. 𝒍𝟏 norm approximation using the hyperbolic tangent function for different values of γ = (1, 4, 6, and 10). As the value of gamma continues to increase and the approximation is closer to the actual 𝒍𝟏 norm, however, it is less smooth. The proposed technique gives us the flexibility to choose between the level of smoothness and accuracy. Unlike the 𝑙ଵ norm, we can apply the unconstraint optimization techniques, where gradient needs to be calculated and the proposed approximation is twice differentiable and the 1st and 2nd order gradients of the proposed 𝑙ଵ norm are shown in Equations (10) and (11), respectively. ∇(‖𝒛‖) ≈ 𝑡𝑎𝑛ℎ(𝛾𝒛) − 𝛾𝒛 (𝑡𝑎𝑛ℎ(𝛾𝒛)ଶ − 1) (10)∇ଶ(‖𝒛‖) ≈ 2𝛾(𝛾𝒛 𝑡𝑎𝑛ℎ(𝛾𝒛) − 1)(𝑡𝑎𝑛ℎ(𝛾𝒛)ଶ − 1) (11)As the value of 𝛾 approaches infinity the error between ‖𝐳‖ଵ and ‖𝒛‖ఊ approaches zero. We here propose the simple lemma to determine the error bounds for ‖𝒛‖ and ‖𝒛‖ఊ. Lemma 1. The proposed smooth function of 𝑙ଵ norm 𝑓(𝑧) = 𝒛 𝑡𝑎𝑛ℎ (𝛾𝒛) fulfils the sufficient and necessary convexity condition in the interval 𝒛 ∈ [−1,1] as its derivative 𝑓′(𝒛) defined by Equation (9) is monotonically non-decreasing and its second derivative 𝑓′′(𝒛) defined by Equation (10) is nonnegative for 0 < 𝛾 ≤ 1; Lemma 2. 𝑙ଵ norm approximation error bounds for any 𝒛 ∈ ℝ and 𝛾 > 0. ห‖𝒛‖ଵ − ‖𝒛‖ఊห ≤12𝛾 (12)Proof. Let us consider two cases, first case for 𝑧 > 0, 𝑝(𝒛, 𝛾) − (𝒛)ା = 𝒛ଶ (1 + 𝑡𝑎𝑛ℎ(𝛾𝑧)) − 𝒛 (13)=𝒛2 (1 + 𝑡𝑎𝑛ℎ(𝛾𝑧)) − 𝒛 = 𝒛2 (tanh(𝛾𝒛) − 1) Now, we can find the maximum value of tanh(𝛾𝑧) to find the upper bound for Equation (13). As we know that the maximum value of tanh(𝛾𝑥) is 1, so we can write: -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 100.10.20.30.40.50.60.70.80.91zL 1 norm approximations L1 Normγ=1γ=4γ=6γ=10Figure 1. l1 norm approximation using the hyperbolic tangent function for different values ofγ = (1, 4, 6, and 10). As the value of gamma continues to increase and the approximation is closero the actual l1 norm, howev r, it is less smooth. Th proposed techniq e gives us the flexibility tochoose b twe n the level of moothness nd accuracy.For steepest descent algorithm, vector differentiation can not be used. Hence, it can berewritten as Equation (4) in lement form to find the artial derivative. Let A = ΦΨH , thenthe element-wise equation is defined as:f (z) =12∑i(Az− y)i(Az− y)i + λzitanh(γzi) (5)Let A = φΨH , then partial derivative of Equation (4) in element form is formulated as:∂ f (z)∂zl=∑ijAijAilzj −∑iyiAil + λ(tanh(γzl) + zlγ(1− tanh2(γzl)))(6)Hence, the steepest descent algorithm for lth update is:(∆z)l = −η∂ f (z)∂zl(7)Equation (7) is used to find a solution using the steepest descent algorithm.2.2. Error Bounds for Proposed Smooth l1-NormThe error bounds for the proposed smooth l1 norm approximation defined by Equation (2)are derived in this section [24]. The l1 norm approximation is proposed based on the fol-lowing two principles.1. |z| = (z)+ + (−z)+, where (z)+ = max{z, 0} is the plus function;2. This plus function can be smoothly approximated as:(z)+ ≈ p(z,γ) =12[z+ z.tanh(γz)] (8)Computation 2023, 11, 7 5 of 21From Equation (8), we can write a smooth approximation of l1 norm:‖z‖1 = (z)+ + (−z)+ ≈ p(z,γ) + p(−z,γ)= z2 [1+ tanh(γz)]− z2 [1+ tanh(−γz)]= z2 tanh(γz) +z2 tanh(γz)= z tanh(γz)= ‖z‖γ(9)Equation (9) represents the γ approximation of the l1 norm, as shown in Figure 1.Unlike the l1 norm, we can apply the unconstraint optimization techniques, where gradientneeds to be calculated and the proposed approximation is twice differentiable and the 1st and2nd order gradients of the proposed l1 norm are shown in Equations (10) and (11), respectively.∇(‖z‖) ≈ tanh(γz)− γz(tanh(γz)2 − 1)(10)∇2(‖z‖) ≈ 2γ(γz tanh(γz)− 1)(tanh(γz)2 − 1)(11)As the value of γ approaches infinity the error between ‖z‖1 and ‖z‖γ approacheszero. We here propose the simple lemma to determine the error bounds for ‖z‖ and ‖z‖γ.Lemma 1. The proposed smooth function of l1 norm f (z) = z tanh(γz) fulfils the suffi-cient and necessary convexity condition in the interval z ∈ [−1, 1] as its derivative f ′(z) de-fined by Equation (9) is monotonically non-decreasing and its second derivative f ′′(z) defined byEquation (10) is nonnegative for 0 < γ ≤ 1;Lemma 2. l1 norm approximation error bounds for any z ∈ R andγ > 0.∣∣∣‖z‖1 − ‖z‖γ∣∣∣ ≤ 12γ (12)Proof. Let us consider two cases, first case for z > 0,p(z,γ)− (z)+ = z2 (1+ tanh(γz))− z= z2 (1+ tanh(γz))− z= z2 (tanh(γz)− 1)(13)Now, we can find the maximum value of tanh(γz) to find the upper bound forEquation (13). As we know that the maximum value of tanh(γx) is 1, so we can write:maximaztanh(γz) =eγz − e−γzeγz + e−γz = 1 (14)Using Equation (14), the relationship between γ and z can be easily derived as:z =12γ(15)By inserting the value of z from Equation (15) in Equation (14)p(z,γ)− (z)+ ≤14γ(16)Figure 2 shows the proposed method obeys error bounds define by Equation (16).Computation 2023, 11, 7 6 of 21Computation 2023, 10, x FOR PEER REVIEW 6 of 21 𝑚𝑎𝑥𝑖𝑚𝑎௭ 𝑡𝑎𝑛ℎ (𝛾𝒛) =𝑒ఊ𝒛 − 𝑒ିఊ𝒛𝑒ఊ𝒛 + 𝑒ିఊ𝒛 = 1 (14)Using Equation (14), the relationship between 𝛾 and 𝑧 can be easily derived as: 𝑧 = 12𝛾 (15)By inserting the value of 𝒛 from Equation (15) in Equation (14) 𝑝(𝒛, 𝛾) − (𝒛)ା ≤14𝛾 (16)□ Figure 2 shows the proposed method obeys error bounds define by Equation (16). Figure 2. 𝒍𝟏 norm approximation error bounds for 𝒛 > 𝟎, the green line shows the upper bound proved mathematically in Equation (16), whereas the dotted red line shows the actual error between the proposed 𝒍𝟏 norm smooth approximation and actual non-differentiable 𝒍𝟏 norm. The error is maximum at approximately zero and approaches zero as 𝜸 → ∞. Figure 3 shows graphically that our proposed approximation function obeys the bounds defined by Equation (19). F re 2. l1 norm approximation err r bounds for z > 0, the green line shows the upper boundroved mathematically in Equation (16), whereas the dotted red line shows the actual error betweenthe proposed l1 norm s ooth approximation and actual non-differentiable l1 norm. The error ismaximum at approximately zero and approaches zero as γ→ ∞ .Figure 3 shows graphically that our proposed approximation function obeys thebounds defined by Equation (19).For z ≤ 0,0 ≤ p(z,γ)− (z)+ = p(z,γ) ≤ p(0,γ)= z2 (tanh(γz)− 1) ≤ 0= 14γ(17)As p is the monotonically increasing function. Hence, from Equations (17) and (18),p(z,γ) will dominate (z)+, so|p(z,γ)− (z)+| ≤ 14γ (18)From Equation (8), we can insert ‖z‖ = (z)+ + (−z)+∣∣∣‖z‖1 − ‖z‖γ∣∣∣ = ∣∣p(z,γ) + p(−z,γ)− ((z)+ + (−z)+)∣∣≤ |p(z,γ)− (z)+|+∣∣p(−z,γ)− (−z)+∣∣≤ 14γ + 14γ = 12γ(19)Figure 3 shows the error bounds versus error in the smooth approximation.Let us define ‖z‖(1,γ) as a smooth approximation to the l1 norm function ‖z‖1 for avector z ∈ Rn as:‖z‖(1,γ) = ∑ni ‖zi‖γ∣∣∣‖z‖(1,γ) − ‖z‖1∣∣∣ ≤ 2n 14γ = n2γ (20)Hence, we can conclude that:Computation 2023, 11, 7 7 of 21limγ→∞‖z‖(1,γ) = ‖z‖1 ∀z ∈ Rn (21)Let L : Rn → R by any continuous cost function and defined by f (z) = L(z) + ‖z‖1and fγ(z) = L(z) + ‖z‖(1,γ). If we define z =argminzf (z) and zγ =argminzfγ(z). Bydefinition of f and fγ and from Equation (20), it can be concluded thatlimγ→∞ fγ(z) = f (z) ∀z ∈ Rn (22)Computation 2023, 10, x FOR PEER REVIEW 7 of 21 Figure 3. 𝒍𝟏 norm approximation error bounds for 𝒛 < 𝟎, the green line shows the upper bound proved mathematically in Equation (19), whereas the dotted red line shows the actual error between proposed 𝒍𝟏 norm smooth approximation and actual non-differentiable 𝒍𝟏 norm. The error is maximum at approximately zero and approaches zero as 𝜸 → ∞. For 𝑧 ≤ 0, 0 ≤ 𝑝(𝑧, 𝛾) − (𝑧)ା = 𝑝(𝑧, 𝛾) ≤ 𝑝(0, 𝛾) = 𝑧2 (𝑡𝑎𝑛ℎ(𝛾𝑧) − 1) ≤ 0 = 14𝛾 (17) As 𝑝 is the monotonically increasing function. Hence, from Equations (17) and (18), 𝑝(𝑧, 𝛾) will dominate (𝑧)ା, so |𝑝(𝑧, 𝛾) − (𝑧)ା| ≤14𝛾 (17)From Equation (8), we can insert ‖𝒛‖ = (𝑧)ା + (−𝑧)ା ห‖𝑧‖ଵ − ‖𝑧‖ఊห = |𝑝(𝑧, 𝛾) + 𝑝(−𝑧, 𝛾) − ((𝑧)ା + (−𝑧)ା)| (1) ≤ |𝑝(𝑧, 𝛾) − (𝑧)ା| + |𝑝(−𝑧, 𝛾) − (−𝑧)ା| ≤ 14𝛾 +14𝛾 =12𝛾 Figure 3 shows the error bounds versus error in the smooth approximation. Let us define ‖𝑧‖(ଵ,ఊ) as a smooth approximation to the 𝑙ଵ norm function ‖𝑧‖ଵ for a vector 𝑧 ∈ ℝ௡ as: ‖𝑧‖(ଵ,ఊ) = ෍ ‖𝑧௜‖ఊ௡௜ ห‖𝑧‖(ଵ,ఊ) − ‖𝑧‖ଵห ≤ 2𝑛14𝛾 =𝑛2𝛾 (2) Hence, we can conclude that: igure 3. l1 norm approximati n error bounds f z < 0, the green line shows the upper bound provedmathematically in Equation (19), whereas the dotted red line shows the actual error between proposedl1 norm smooth approximation and actual non-differentiable l1 norm. The error is maximum atapproximately zero and approaches zero as γ→ ∞ .In addition, it is a known fact that f (z) ≤ f (z)∀z. In particular f (z) ≤ f (zγ), then:f (z) ≤ f (zγ) = L(zγ) + ‖zγ‖1= L(zγ) + ‖zγ‖1 + ‖zγ‖(1,γ − ‖zγ‖(1,γ)=(L(zγ) + ‖zγ‖(1,γ))+(‖zγ‖1 − ‖zγ‖(1,γ))= fγ(zγ) +(‖zγ‖1 − ‖zγ‖( ,γ)) (23)This implies that f (z)− fγ(zγ) ≥ − n2γ from Equation (21), similarly (z)− fγ(zγ) ≤ n2γ ,hence proved that limγ→∞ fγ(z) = f (z).It can be further stated that:| f (zγ)− f (z)| = | f (zγ)− f (z)− fγ(zγ) + fγ(zγ)|≤ | f (zγ)− fγ(zγ)|+ | fγ(zγ)− f (z)| (24)Hence, it proved that limγ→∞ f (zγ) = f (z). Moreover, if L is strictly convex, it can beeasily proven that: limγ→∞zγ = z.Computation 2023, 11, 7 8 of 213. The Map Estimator and Proposed Thresholding MechanismConventionally l1-norm minimization has inherent soft thresholding [13]. However,when we approximate the l1-norm by the hyperbolic tangent function, thresholding is notdone implicitly. The hard thresholding operator proposed in [14] can be defined by thefollowing equation.Sβ(z) ={z |z| > β0 otherwise(25)We have used a new thresholding function based on the tangent hyperbolic function.Therefore, β is an important parameter in controlling the under-sampling noise, whichhas Gaussian distribution [21]. To obtain the optimum value of β, the thresholding de-pends upon the undersampling noise. Therefore, the following data-driven thresholdingparameter β is used [13,25].β =σ2vσz(26)With σz is the standard deviation of sparse signal and σv the standard deviation of Gaussian-like noise produced due to under-sampling.To enhance the performance under different scenarios, different mathematical thresh-olding operators could be found in the literature [26–28]. The main idea in this approach ismapping the values nearer to the origin to zero and those that are further away from theorigin are shrunk towards zero.The basic denoising technique aims to find the estimate of original image or signalfrom its perturbed set of observations, as shown in Equation (27).y = z+ v (27)where y ∈ Rn is the noisy image, z ∈ Rn is the original signal and v is the zero-meanGaussian noise with probability distribution function (pdf) given by:pv(θ) =1√2piσ2vexp(‖θ‖222σ2v)(28)By taking the Wavelet transform of Equation (27), we get:q = s+ v (29)where q = Ψy and s = Ψz represent the sparsifying domain for noisy image and theoriginal image, respectively. As Wavelet transform is the linear operator, therefore thezero-mean Gaussian noise v after transformation will not change. The MAP estimation ofrandom vector s is given by:sˆ = maxs∈Rnp(s|q) (30)By using Bayes’ rule, one can ignore p(q) as it is independent of s, MAP estimator canbe written as:sˆ = maxs∈Rnp(q|s)ps(s) (31)The problem defined in Equation (31) can be further simplified by taking p(q|s) = pv(q− s) :sˆ = maxs[pv(q− s)]ps(s)= maxs[lnpv(q− s) + lnps(s)]= maxs[ln{1√2piσ2vexp(− ‖q−s‖222σ2v)}+ lnps(s)]= maxs[ln{(1√2piσ2v)nexp(− ‖q−s‖222σ2v)}+ lnps(s)]= maxs[− ‖q−s‖222σ2v+ f (s)](32)Computation 2023, 11, 7 9 of 21where f (s) = ln ps(s). By differentiating the argument of Equation (32) w.r.t. s andequating the result to zero, we can calculate the MAP estimator for Wavelet coefficients ofthe noise-free image as:(qi − sˆi)σ2v+ f ′(sˆi) = 0, 1 ≤ i ≤ n (33)The pdf of biomedical images are more peaked at the center than Gaussian, so Lapla-cian can better estimate the distribution of Wavelet domain coefficients, i.e.,ps(si) =1√2σvexp(√2σv|si|)(34)gives f ′(sˆi) = −√2σ sig(sˆi). Solving Equation (33) will result inqi = sˆi +√2 sig(sˆi) (35)Let β =√2σ2v and solve Equation (35) for sˆi to formulate the nonlinear shrinkage:sˆi = Sβ(q) = max{|q| − β, 0}.sig(q) (36)Equation (36) can be further elaborated as:Sβ(q) ∼={sgn(q)(|q| − β) |q| > β0 otherwise(37)In this paper, novel thresholding approach has been proposed, which is used onthe hyperbolic tangent, as the hyperbolic tangent function slope can be adjusted fromthe origin and it is a bounded function that makes it an suitable surrogate function forsoft thresholding. Hence, hyperbolic tangent based soft thresholding can be describedmathematically by following equation:Sβ(q) ∼={cz{tanh(α(|q| − β))} |q| > β0 otherwise(38)where β is a thresholding parameter and parameter α is used to control the shape of thehyperbolic tangent function. If α is closer to zero, Equation (38) approximately changesinto the soft thresholding function. When α approaches ∞, Equation (38) changes to thehard thresholding function, as shown in Figure 4. Our proposed Algorithm 1 starts as asoft thresholding function and smoothly changes to a hard threshold at higher iterations.The proposed approximation of this soft thresholding results in a better reconstruction ascompared to the conventional soft thresholding method [21] as illustrated in Figure 4.Algorithm 1. Proposed Algorithms.Inputs:Sensing matrix Fu, measurement vector yeCm, parameters γ, λ and β,Output:A k-sparse vector xˆ ∈ RnInitialization: Initialize x0, Index i = 0Step-1 (Sparse Representation): zi = ΨxiStep-1 (Gradient Computation): Find ∇ f (zi) using Equation (5)Step-2 (Solution Update): Compute the update using Equations (6) and (7).Step-3 (Shrinkage): Estimate Solution using Equation (38), i.e., zˆi+1 = Sβ(zi+1)Step-4 (Repeat): If stopping criterion is not met, i = i + 1 & go to step 1Output: xˆ = ΨH zˆiComputation 2023, 11, 7 10 of 21Computation 2023, 10, x FOR PEER REVIEW 10 of 21 thresholding. Hence, hyperbolic tangent based soft thresholding can be described mathe-matically by following equation: 𝑆ఉ(𝑞) ≅ ൜𝑐𝑧{𝑡𝑎𝑛ℎ൫𝛼(|𝑞| − 𝛽)൯ } |𝑞| > 𝛽0 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 (19)where 𝛽 is a thresholding parameter and parameter 𝛼 is used to control the shape of the hyperbolic tangent function. If 𝛼 is closer to zero, Equation (38) approximately changes into the soft thresholding function. When 𝛼 approaches ∞, Equation (38) changes to the hard thresholding function, as shown in Figure 4. Our proposed Algorithm 1 starts as a soft thresholding function and smoothly changes to a hard threshold at higher iterations. The proposed approximation of this soft thresholding results in a better reconstruction as compared to the conventional soft thresholding method [21] as illustrated in Figure 4. Figure 4. Hyperbolic tangent function based thresholding for alpha 𝛼 = (2, 4, 8, 16), the value of 𝛼 determines the slope of soft thresholding. The proposed method gives us the flexibility to shape the curves using 𝛼 depending upon its application. Algorithm 1. Proposed Algorithms. Inputs: Sensing matrix ℱ௨, measurement vector 𝑦𝜖ℂ௠, parameters 𝛾, 𝜆 and 𝛽, Output: A k-sparse vector 𝐱ො ∈ 𝑅௡ Initialization: Initialize 𝐱଴, Index 𝑖 = 0 Step-1 (Sparse Representation): 𝒛𝒊 = 𝚿𝒙𝒊 Step-1 (Gradient Computation): Find ∇𝑓(𝑧௜) using Equation (5) Step-2 (Solution Update): Compute the update using Equations (6) and (7). Step-3 (Shrinkage): Estimate Solution using Equation (38), i.e., 𝒛ො𝒊ା𝟏 =𝑆ఉ(𝐳௜ାଵ) Step-4 (Repeat): If stopping criterion is not met, 𝑖 = 𝑖 + 1 & go to step 1 Output: xො = Ψୌ?̂?௜ Figure 4. yperbolic tangent function based thresholding for alpha α = (2, 4, 8, 16), the value of αdetermines the slope of soft thresholding. The proposed method gives us the flexibility to shape thecurves using α depending upon its application.4. Results and DiscussionsIn order to gauge the performance of proposed algorithm, we have applied ouralgorithm to 1-D sparse signal, Co pressively sampled MR image and Cine Cardiac MRI.MRXCAT simulator is used to evaluate the proficiency of recovery algorithms in the field ofCardiac MRI. We have evaluated the performance of our proposed technique quantitativelyand qualitatively. The performance measures that are used in this research article are:pictorial depiction of under-sampling artefacts, Structural Similarity (SSIM), Peak Signal toNoise Ratio (PSNR) and Root Mean Square Error (RMSE).4.1. 1-D Sparse Signal RecoveryThe proposed algorithm is applied for the recovery of the 1-D sparse signal recovery,where the random sparse signal of length n = 512 is created in MATLAB and the supportfor the sparse signal was generated randomly with K = 85 non-zero elements. The randomsparse signal is compressively sampled using a random measurement matrix A ∈ R256×512with only m = 256 measurements.Figure 5 shows the fitness achieved by the proposed algorithm and soft-thresholdingmethod. The proposed method achieved faster convergence as compared to soft threshold-ing. Figure 6 shows sparsity effect on successful recovery achieved by the soft thresholdingand proposed algorithm. The proposed algorithm performs much better even with ahigher sparsity level as compared to the soft thresholding technique. Similarly, in Figure 7,the proposed method recovered the sparse signal with great accuracy, whereas the softthresholding technique failed to accurately recover the sparse signal. The accuracy of theproposed technique was also measured against performance measures, such as SNR, MSEand correlation, as shown in Table 1. The proposed algorithm performed much better thanthe soft thresholding method against all these performance measures. The time comparisonfor proposed algorithm is 1.57 s as compared to 1.34 s by conventional soft thresholding.Computation 2023, 11, 7 11 of 21Computation 2023, 10, x FOR PEER REVIEW 11 of 21 4. Results and Discussions In order to gauge the performance of proposed algorithm, we have applied our algo-rithm to 1-D sparse signal, Compressively sampled MR image and Cine Cardiac MRI. MRXCAT simulator is used to evaluate the proficiency of recovery algorithms in the field of Cardiac MRI. We have evaluated the performance of our proposed technique quantita-tively and qualitatively. The performance measures that are used in this research article are: pictorial depiction of under-sampling artefacts, Structural Similarity (SSIM), Peak Sig-nal to Noise Ratio (PSNR) and Root Mean Square Error (RMSE). 4.1. 1-D Sparse Signal Recovery The proposed algorithm is applied for the recovery of the 1-D sparse signal recovery, where the random sparse signal of length n = 512 is created in MATLAB and the support for the sparse signal was generated randomly with K = 85 non-zero elements. The random sparse signal is compressively sampled using a random measurement matrix 𝐴 ∈ 𝑅ଶହ଺௫ହଵଶ with only m = 256 measurements. Figure 5 shows the fitness achieved by the proposed algorithm and soft-thresholding method. The proposed method achieved faster convergence as compared to soft thresh-olding. Figure 6 shows sparsity effect on successful recovery achieved by the soft thresh-olding and proposed algorithm. The proposed algorithm performs much better even with a higher sparsity level as compared to the soft thresholding technique. Similarly, in Figure 7, the proposed method recovered the sparse signal with great accuracy, whereas the soft thresholding technique failed to accurately recover the sparse signal. The accuracy of the proposed technique was also measured against performance measures, such as SNR, MSE and correlation, as shown in Table 1. The proposed algorithm performed much better than the soft thresholding method against all these performance measures. The time compari-son for proposed algorithm is 1.57 s as compared to 1.34 s by conventional soft thresholding. Figure 5. Fitness achieved by the soft thresholding and proposed algorithm. The proposed algo-rithm converges rapidly as compared to the soft thresholding technique. Figure 5. Fitness achieved by the soft thresholding and proposed algorithm. The proposed algorithmconverges rapidly as compared to the soft thresholding technique.Computation 2023, 10, x FOR PEER REVIEW 12 of 21 Figure 6. Sparsity effect on successful recovery achieved by the soft thresholding and proposed al-gorithm. The proposed algorithm performs much better even with higher sparsity level as com-pared to the soft thresholding technique. Figure 7. (a) The recovered sparse signal from the proposed algorithm; (b) The recovered sparse signal from soft thresholding. Table 1. Performance comparison of different sparsity transforms using mean squared error in the transform domain. Temporal FFT performs better in cardiac cine MRI. Performance Metrics Soft Thresholding Proposed Algorithm MSE 1.00 × 10−2 1.61 × 10−4 Fitness 0.8664 0.0224 SNR 12.6712 30.6259 Correlation 0.9787 0.9995 Figure 6. Sparsity effect on successful recovery achieved by the soft thresholding and proposedalgorithm. The proposed algorithm performs much better even with higher sparsity level as comparedto the soft thresholding technique.Computation 2023, 11, 7 12 of 21Computation 2023, 10, x FOR PEER REVIEW 12 of 21 Figure 6. Sparsity effect on successful recovery achieved by the soft thresholding and proposed al-gorithm. The proposed algorithm performs much better even with higher sparsity level as com-pared to the soft thresholding technique. Figure 7. (a) The recovered sparse signal from the proposed algorithm; (b) The recovered sparse signal from soft thresholding. Table 1. Performance comparison of different sparsity transforms using mean squared error in the transform domain. Temporal FFT performs better in cardiac cine MRI. Performance Metrics Soft Thresholding Proposed Algorithm MSE 1.00 × 10−2 1.61 × 10−4 Fitness 0.8664 0.0224 SNR 12.6712 30.6259 Correlation 0.9787 0.9995 i . (a) The recovered sparse signal from the proposed algorithm; (b) Th recovered sparsesignal fro soft thresholding.a le 1. erfor a ce co ariso of iffere t s arsity tra sfor s si g ea s are error i t etransform domain. Temporal FFT performs better in cardiac cine MRI.Performance Metrics Soft Thresholding Proposed AlgorithmMSE 1.00 × 10−2 1.61 × 10−4Fitness 0.8664 0.0224SNR 12.6712 30.6259Correlation 0.9787 0.99954.2. 2-D Compressively Sampled MR Image RecoveryThe random sampling at the CS image acquisition produces incoherent and noise-likeartefacts in its sparsifying domain. In case of MR imaging or similar Fourier domainencoded biomedical imaging, where the MR image is in the spatial domain, the linearreconstruction (where, missing Fourier data points are replaced by zero and the resultantimage inverse Fourier transform is taken) produces artifacts similar to additive Gaus-sian noise. The type of noise produced by subsampling is governed by undersamplingpatterns [29]. In order to recover an image, the compressed sensing recovery essentiallybecomes an image denoising problem. Using this analogy of CS encoding and noisy image,the first step in recovering the original image is to estimate a noise, this is achieved bymaximum a posteriori (MAP) estimator. The proposed algorithm is also implemented torecover a 2-D Compressively sampled real human brain MR image of size 256 × 256. Thehuman brain MR image is a fully sampled scanned image by a 1.5 Tesla GE-HDxt-MRIscanner with Gradient Echo (GE) sequence and 8 channels head coils with the specifications,i.e., TE = 10 msec, flip angle = 90◦, bandwidth = 31.25 KHz, slice thickness=3 mm, TR = 55,and image dimensions = 256 × 256, at St. Mary’s Hospital, London, UK. This MR image iscompressively sampled by taking only 25% samples in k-space.Figure 8 shows the performance of the proposed method with respect to StructuralSimilarity (SSIM). The proposed method achieved much better SSIM as compared tosoft thresholding. The Peak Signal to Noise Ratio (PSNR) accomplished by proposedmethod is shown in Figure 9. Figure 10 shows the (a) Original 2D Brain MR Image,Computation 2023, 11, 7 13 of 21(b) Conventional Soft Thresholding based recovered 2D MRI, (c) 2D Brain MR Imagerecovered from undersampled image, (d) Difference of original and soft thresholding image,(e) Difference of proposed recovery method image with original image. The Difference isscaled up by 1000 in order to enhance its visibility. Table 2 shows the proposed methodhas outperformed soft thresholding method in terms of PSNR and SSIM. Table 3 showsthe performance of proposed algorithm and soft thresholding in terms of Mean SquareError (MSE), Improved Signal to Noise Ratio (ISNR), Correlation, SSIM, SNR, PSNR afterthe 15 iterations of soft thresholding and the proposed algorithm. The results show thatproposed method achieved much better results in as compared to soft thresholding.Computation 2023, 10, x FOR PEER REVIEW 13 of 21 4.2. 2-D Compressively Sampled MR Image Recovery The random sampling at the CS image acquisition produces incoherent and noise-like artefacts in its sparsifying domain. In case of MR imaging or similar Fourier domain encoded biomedical imaging, where the MR image is in the spatial domain, the linear reconstruction (where, missing Fourier data points are replaced by zero and the resultant image inverse Fourier transform is taken) produces artifacts similar to additive Gaussian noise. The type of noise produced by subsampling is governed by undersampling patterns [29]. In order to recover an image, the compressed sensing recovery essentially becomes an image denoising problem. Using this analogy of CS encoding and noisy image, the first step in recovering the original image is to estimate a noise, this is achieved by maximum a posteriori (MAP) estimator. The proposed algorithm is also implemented to recover a 2-D Compressively sampled real human brain MR image of size 256 × 256. The human brain MR image is a fully sampled scanned image by a 1.5 Tesla GE-HDxt-MRI scanner with Gradient Echo (GE) sequence and 8 channels head coils with the specifications, i.e., TE = 10 msec, flip angle = 90°, bandwidth = 31.25 KHz, slice thickness=3 mm, TR = 55, and image dimensions = 256 × 256, at St. Mary’s Hospital, London, UK. This MR image is compres-sively sampled by taking only 25% samples in k-space. Figure 8 shows the performance of the proposed method with respect to Structural Similarity (SSIM). The proposed method achieved much better SSIM as compared to soft thresholding. The Peak Signal to Noise Ratio (PSNR) accomplished by proposed method is shown in Figure 9. Figure 10 shows the (a) Original 2D Brain MR Image, (b) Conven-tional Soft Thresholding based recovered 2D MRI, (c) 2D Brain MR Image recovered from undersampled image, (d) Difference of original and soft thresholding image, (e) Differ-ence of proposed recovery method image with original image. The Difference is scaled up by 1000 in order to enhance its visibility. Table 2 shows the proposed method has outper-formed soft thresholding method in terms of PSNR and SSIM . Table 3 shows the perfor-mance of proposed algorithm and soft thresholding in terms of Mean Square Error (MSE), Improved Signal to Noise Ratio (ISNR), Correlation, SSIM, SNR, PSNR after the 15 itera-tions of soft thresholding and the proposed algorithm. The results show that proposed method achieved much better results in as compared to soft thresholding. Figure 8. Structural Similarity of proposed and soft thresholding algorithm for recovery of compres-sively sampled MR image against each iteration. i r . tr t r l i il rit f r s s ft t r s l i l rit f r r r f c r s-sively sa pled R i age against each iteration.Computation 2023, 10, x FOR PEER REVIEW 14 of 21 Figure 9. Correlation of proposed and soft thresholding algorithm of recovered compressively sam-ple MR image. Figure 10. (a) Original 2D Brain MR Image, (b) Conventional Soft Thresholding based recovered 2D MRI, (c) 2D Brain MR Image recovered from undersampled image, (d) Difference of original and soft thresholding image, (e) Difference of proposed recovery method image with original image. The Difference is scaled up by the factor of 1000 in order to enhance its visibility. . se and soft thresholding algorithm of recovered compressively sampleMR image.Computation 2023, 11, 7 14 of 21Computation 2023, 10, x FOR PEER REVIEW 14 of 21 Figure 9. Correlation of proposed and soft thresholding algorithm of recovered compressively sam-ple MR image. Figure 10. (a) Original 2D Brain MR Image, (b) Conventional Soft Thresholding based recovered 2D MRI, (c) 2D Brain MR Image recovered from undersampled image, (d) Difference of original and soft thresholding image, (e) Difference of proposed recovery method image with original image. The Difference is scaled up by the factor of 1000 in order to enhance its visibility. Figure 10. (a) Original 2D Brain MR Image, (b) Conventional Soft Thresholding based recovered 2DMRI, (c) 2D Brain MR Image recovered from undersampled image, (d) Difference of original and softthresholding image, (e) Difference of proposed recovery method image with original image. TheDifference is scaled up by the factor of 1000 in order to enhance its visibility.Table 2. Performance comparison of conventional soft thresholding and proposed method withdifferent compression levels, i.e., 5% to 50% of subsampling of the original 2-D MR image. Theseresults show that the proposed method achieves better results in terms of SSIM and PSNR at varyingcompression ratios.Compression RatioSoft Thresholding Proposed AlgorithmSSIM PSNR SSIM PSNR5 % 0.6843 75.9056 0.7048 76.160910% 0.7786 78.9320 0.8175 79.658020% 0.8994 82.0316 0.8472 83.762830% 0.9407 87.3535 0.9790 91.162040% 0.9724 91.2540 0.9920 96.128150% 0.9884 95.4245 0.9955 99.5496Table 3. Performance comparison of different sparsity transforms using mean squared error in thetransform domain. Temporal FFT performs better in cardiac Cine MRI.Performance Metrics Soft Thresholding Proposed AlgorithmMSE 1.38 × 10−4 0.73 × 10−4PSNR 86.7195 89.4497ISNR 28.3832 31.1135SSIM 0.9346 0.9711SNR 26.0298 28.7491Correlation 0.9980 0.99894.3. Cardiac Cine Magnetic Resonance Imaging RecoveryThe proposed algorithm is applied to MRXCAT, which produces breath-held under-sampled cardiac cine MR image data. For MRXCAT, the following parameters were set:Computation 2023, 11, 7 15 of 21recovery matrix size: 256 × 256 with 24 cardiac phases, with an image resolution set as1 × 1 × 1 mm3, TR = 3 ms, TE = 1.5 ms. Five different acceleration rates R = (2, 4, 8,12, 20) were used to assess the performance of the proposed method. For in vivo data,the following parameters were used: reconstruction matrix size: 256 × 256, 25 cardiacphases, with FOV of 375 mm. TE = 1 ms, TR = 3 ms and flip angle = 600. Five accelerationrates are R = (2, 4, 8, 12, 20) are used to evaluate the performance of the proposedmethod. The reconstructed images are matched with the fully sampled original generatedcardiac cine MRI as shown in Figure 11. All images are recovered in MATLAB by theproposed algorithm.Computation 2023, 10, x FOR PEER REVIEW 15 of 21 Table 2. Performance comparison of conventional soft thresholding and proposed method with dif-ferent compression levels, i.e., 5% to 50% of subsampling of the original 2-D MR image. These results show that the proposed method achieves better results in terms of SSIM and PSNR at varying com-pression ratios. Compression Ratio Soft Thresholding Proposed Algorithm SSIM PSNR SSIM PSNR 5 % 0.6843 75.9056 0.7048 76.1609 10% 0.7786 78.9320 0.8175 79.6580 20% 0.8994 82.0316 0.8472 83.7628 30% 0.9407 87.3535 0.9790 91.1620 40% 0.9724 91.2540 0.9920 96.1281 50% 0.9884 95.4245 0.9955 99.5496 Table 3. Performance comparison of different sparsity transforms using mean squared error in the transform domain. Temporal FFT performs better in cardiac Cine MRI. Performance Metrics Soft Thresholding Proposed Algorithm MSE 1.38 × 10−4 0.73 × 10−4 PSNR 86.7195 89.4497 ISNR 28.3832 31.1135 SSIM 0.9346 0.9711 SNR 26.0298 28.7491 Correlation 0.9980 0.9989 4.3. Cardiac Cine Magnetic Resonance Imaging Recovery The proposed algorithm is applied to MRXCAT, which produces breath-held under-sampled cardiac cine MR image data. For MRXCAT, the following parameters were set: recovery matrix size: 256 × 256 with 24 cardiac phases, with an image resolution set as 1 × 1 × 1 mm3, TR = 3 ms, TE = 1.5 ms. Five different acceleration rates R = (2, 4, 8, 12, 20) were used to assess the performance of the proposed method. For in vivo data, the following parameters were used: construction matrix size: 256 × 256, 25 cardiac phases, with FOV of 375 mm. TE = 1 ms, TR = 3 ms and flip angle = 600. Five acceleration rates are 𝑅 = (2, 4, 8, 12, 20) are us d to evaluate the performance of the propos d m th d. The recon-structed images are matched with the fully sampled original generated cardiac cine MRI as shown in Figure 11. All images are recovered in MATLAB by the proposed algorithm. (a) (b) (c) Figure 11. (a) Short axis cardiac cine MRI with completely sampled diastolic frame. (b) Sparsifying transform of cine cardiac MRI diastolic frame with temporal Fourier transform (Ψ), which results in sparse representation, (c) Another sparse representation of cardiac cine MR image (diastolic frame) using total variation transform (Ψ). To gauge the efficiency of proposed algorithm, we use MRXCAT simulator software. It is designed for the analysis of reconstruction algorithms performance in the area of car-diac cine MRI. MRXCAT simulator is used to evaluate the proficiency of recovery algo-rithms in the field of Cardiac MRI. We have evaluated the performance of our proposed Figure 11. (a) Short axis cardiac cine MRI with completely sa led diastolic frame. (b) Sparsifyingtransform of cine cardiac MRI diastolic frame with temporal Fourier transform (Ψ), which results insparse representation, (c) Another sparse representation of cardiac cine MR image (diastolic frame)using total variation transform (Ψ).To gauge the efficiency of proposed algorithm, we use MRXCAT simulator software. Itis designed for the analysis of reconstruction algorithms performance in the area of cardiaccine MRI. MRXCAT simulator is used to evaluate the proficiency of recovery algorithms inthe field of Cardiac MRI. We have evaluated the performance of our proposed techniquequantitatively and qualitatively. The performance measures that are used in this researcharticle are: pictorial depiction of under-sampling artefacts, Structural Similarity (SSIM),Peak Sign l to Noise Ratio (PSNR) nd Root Mean Square Error (RMSE).To evaluate the performance of proposed recovery technique qualitatively, we haveexperimentally depicted the recovered diastolic and systolic frames using acceleration ratesof (R = 2, 4, 8, 12, 20). The quantitative assessment of the proposed algorithm is done usingRMSE, PSNR and SSIM. Comparison between proposed algorithm and traditional softthresholding technique is also performed. Figure 11 depicts the proficiency of the proposedalgorithm at various acceleration rates while comparing with the soft thresholding. The firstcolumn shows the diastolic frame at different acceleration rates of cine cardiac MR imageand the second column represents the systolic frame of cine MRI. The top row depicts theresults of traditional IST algorithm, while the bottom row depicts proposed method results.Table 4 shows the performance comparison of different sparsity transform using meansquared error in the transform domain. It measures the average of the error squares betweenthe reconstructed and the acquired coefficients in the sparse domain. The proposed tangenthyperbolic based approximation performs well in temporal FFT as compared to the othersparse transforms. In particular, at higher acceleration rates, the tangent hyperbolic tangentbased proposed technique shows much improved recovery of CS images.Computation 2023, 11, 7 16 of 21Table 4. Performance comparison of different sparsity transforms using mean squared error in thetransform domain. Temporal FFT performs better in cardiac cine MRI.Acceleration Rates Spatial Domain Total Variation Temporal FFT2 0.1096 0.1123 0.07284 0.2321 0.1849 0.08488 0.2810 0.2438 0.094812 0.3533 0.2684 0.104320 0.4756 0.2982 0.1150Figure 12 shows simulated data where (a) compares proposed method at bottom rowwith IST at top row with acceleration rate of 2. The arrow in (a) depicts the very minutepresence of artefacts, while (b) depicts the performance of proposed method at accelerationrate equal to 4. The arrow in (b) depicts the presence of artefacts (c) shows the results ofboth algorithms with acceleration rate set at 8. The artefacts due to subsampling becomegradually more visible in IST results as highlighted by arrow mark (d) depicts the resultswhen acceleration rate is set at 12. Both techniques depicts the artefacts, however theseartefacts are visible in the IST as mentioned by the white arrow in the figure (e) shows verymuch degraded image quality of IST, while comparing it with the proposed method. Thesubsampling artefacts dominate the traditional IST result when acceleration rate R is setat 20.Computation 2023, 10, x FOR PEER REVIEW 17 of 21 20). We have comp red our method with the tr ditional iterative soft th eshol ing tech-nique. Figure 14 depi ts the effici ncy of our pr posed technique at various cceleration rates as compared to IST algorithm. Red line sh ws the results of our method at the dif-ferent acceleration rates. Figure 15 shows the PSNR of soft thresholding method and un-der-sampled data are shown with the green and blue lines, respectively. Figure 12. Simulated data (a) Compares proposed method at bottom row with IST at top row with acceleration rate of 2. The arrow in (a) depicts the very minute presence of artifacts. Here (b) depicts the performance of proposed method at acceleration rate equal to 4. The arrow in (b) depicts the presence of artefacts (c) shows the results of both algorithms with acceleration rate set at 8. (d) de-picts the results when acceleration rate is set at 12. These artefacts are visible in the IST as mentioned by the white arrow in the figure (e) shows very much degraded image quality of IST, while com-paring it with the proposed method. . Figure 13. This figure depicts the efficiency of proposed algorithm by means of SSIM index. As acceleration rate increases, the SSIM of proposed algorithm degrades slowly while comparing it with IST algorithm. 2 4 6 8 10 12 14 16 18 200.40.450.50.550.60.650.70.750.80.850.9Acceleration RateStructural Similarity Index ProposedSoft ThresholdingUnder-sampled ImageFigure 12. Simulated data (a) Compares proposed method at bottom row with IST at top row withacceleration rate of 2. The arrow in (a) depicts the very minute presence of artifacts. Here (b) depictsthe performance of proposed method at acceleration rate equal to 4. The arrow in (b) depicts thepresence of artefacts (c) shows the results of both algorithms with acceleration rate set at 8. (d) depictsthe results when acceleration rate is set at 12. These artefacts are visible in the IST as mentioned bythe white arrow in the figure (e) shows very much degraded image quality of IST, while comparing itwith the proposed method.Computation 2023, 11, 7 17 of 21To evaluate the recovered images quantitatively, we have used SSIM to compare theproposed technique with the IST technique. Figure 13 depicts the SSIM of our proposedalgorithm, iterative soft thresholding (IST) and undersampled images. The quality ofundersampled images is visibly quite poor. The efficiency of our proposed technique andIST technique is almost similar at low acceleration rates. However, the visible quality of thesoft thresholding based recovered images decreases as acceleration rates are increased, ascompared to the proposed method.Computation 2023, 10, x FOR PEER REVIEW 17 of 21 20). We have compared our method with the traditional iterative soft thresholding tech-nique. Figure 14 depicts the efficiency of our proposed technique at various acceleration rates as compared to IST algorithm. Red line shows the results of our method at the dif-ferent acceleration rates. Figure 15 shows the PSNR of soft thresholding method and un-der-sampled data are shown with the green and blue lines, respectively. Figure 12. Simulated data (a) Compares proposed method at bottom row with IST at top row with acceleration rate of 2. The arrow in (a) depicts the very minute presence of artifacts. Here (b) depicts the performance of proposed method at acceleration rate equal to 4. The arrow in (b) depicts the presence of artefacts (c) shows the results of both algorithms with acceleration rate set at 8. (d) de-picts the results when acceleration rate is set at 12. These artefacts are visible in the IST as mentioned by the white arrow in the figure (e) shows very much degraded image quality of IST, while com-paring it with the proposed method. . Figure 13. This figure depicts the efficiency of proposed algorithm by means of SSIM index. As acceleration rate increases, the SSIM of proposed algorithm degrades slowly while comparing it with IST algorithm. 2 4 6 8 10 12 14 16 18 200.40.450.50.550.60.650.70.750.80.850.9Acceleration RateStructural Similarity Index ProposedSoft ThresholdingUnder-sampled ImageFigure 13. This figure depicts the efficiency of proposed algorithm by means of SSIM index. Asacc leration rate increases, I of proposed algorithm degrades slowly while comparing it wi hIST algorithm.Figure 14 shows the requisite iterations for the image reconstruction in both themethods. The proposed method solves the problem in six iterations, while the IST recov-ery technique takes ten iterations to reach the optimal solution. In this result, the dataconsistency is used to show the performance of proposed method.Computation 2023, 10, x FOR PEER REVIEW 18 of 21 Figure 14. Comparison of cardiac cine MRI recovery at number of iterations; it can be seen from results that proposed method converges to an optimal solution in lesser iterations as compared to traditional thresholding. Table 5 elaborates the comparison of performance using root mean squared error (RMSE). It measures the error squares average between the recovered samples and the actual samples. The proposed method performance is much superior as compared to the IST technique. In particular, while operating at higher acceleration rates, the tangent hy-perbolic method performs much better in recovering the images. However, the efficiency of traditional IST algorithm degrades at higher acceleration rates. Table 5. Comparison of proposed method with conventional IST algorithm with RMSE. Proposed method performance is much better as acceleration rates are increased. Acceleration Rates Undersampled Image Iterative Soft Thresholding Proposed Method Simulated Data 2 0.081 0.0365 0.0353 4 0.1218 0.0472 0.0372 8 0.1498 0.0702 0.0419 12 0.1583 0.0775 0.0485 20 0.1782 0.0941 0.0606 In vivo Data 2 0.085 0.0099 0.0056 4 0.106 0.0241 0.0172 8 0.1170 0.0495 0.0206 12 0.120 0.0567 0.0338 20 0.1398 0.0585 0.0551 To evaluate the performance of reconstructed images qualitatively, using in vivo data, we have shown a comparison between the proposed method and the soft threshold-ing method in Figure 16. We have used five acceleration rates R = (2, 4, 8, 12, and 20) to show the comparison between the proposed method and the soft thresholding method. In this figure, the performance of the proposed technique and IST algorithm is similar at lower acceleration rates. However, the blurring artefacts in IST based recovered images are more prominent at higher acceleration rates as compared to the proposed algorithm as indicated by the white arrows. Figure 14. Comparison of cardiac cine MRI recovery at number of iterations; it can be s en fromresults that proposed method converges to an optimal solution in le ser iterations as compared totraditional thresholding.To valuate the efficiency of our recovered algorithm quanti atively, the sults a eshown the reconstruction results using PSNR at various accelerati n rates (R = 2, 4, 8, 12, 20).We have compared our method with th traditional iterative soft thresholding technique.Figure 14 depicts the efficiency of our proposed technique at various acceleratio ratesComputation 2023, 11, 7 18 of 21as compared to IST algorithm. Red line shows the results of our method at the differentacceleration rates. Figure 15 shows the PSNR of soft thresholding method and under-sampled data are shown with the green and blue lines, respectively.Computation 2023, 10, x FOR PEER REVIEW 19 of 21 Figure 15. The performance using the peak signal-to-noise ratio (PSNR). PSNR of our method is better at all acceleration rates as compared to the soft thresholding method. Figure 16. In real vivo data (a) Compares proposed method at bottom row with IST at top row with acceleration rate of 2. The arrow in (a) shows very minute artifacts. (b) depicts the performance of proposed method at acceleration rate equal to 4. The arrow in (b) depicts the presence of artefacts (c) shows the results of both algorithms with acceleration rate set at 8. The artefacts due to subsam-pling become gradually more visible in IST results as highlighted by arrow mark (d) depicts the results when acceleration rate is set at 12. Both techniques depicts the artefacts, however these arte-facts are visible in the IST as mentioned by the white arrow in the figure (e) show very much de-graded image quality of IST, while comparing it with the proposed method. The subsampling arte-facts dominate the traditional IST result when acceleration rate R is set at 20. 5. Conclusions 2 4 6 8 10 12 14 16 18 2081012141618202224Acceleration RatePSNR ProposedUnder-sampled ImageSoft ThresholdingFigure 15. The performance using the peak signal-to-noise ratio (PSNR). PSNR of our method isbetter at all acceleration rates as compared to the soft thresholding method.Table 5 elaborates the comparison of performance using root mean squared error(RMSE). It measures the error squares average between the recovered samples and theactual samples. The proposed method performance is much superior as compared tothe IST technique. In particular, while operating at higher acceleration rates, the tangenthyperbolic method performs much better in recovering the images. However, the efficiencyof traditional IST algorithm degrades at higher acceleration rates.Table 5. Comparison of proposed method with conventional IST algorithm with RMSE. Proposedmethod performance is much better as acceleration rates are increased.AccelerationRatesUndersampledImageIterative SoftThresholdingProposedMethodSimulated Data2 0.081 0.0365 0.03534 0.1218 0.0472 0.03728 0.1498 0.0702 0.041912 0.1583 0.0775 0.048520 0.1782 0.0941 0.0606In vivo Data2 0.085 0.0099 0.00564 0.106 0.0241 0.01728 0.1170 0.0495 0.020612 0.120 0.0567 0.033820 0.1398 0.0585 0.0551To evaluate the performance of reconstructed images qualitatively, using in vivo data,we have shown a comparison between the proposed method and the soft thresholdingmethod in Figure 16. We have used five acceleration rates R = (2, 4, 8, 12, and 20) to showthe comparison between the proposed m thod an the soft thresholding method. In thisfigure, the performance f the proposed technique and IST algorithm is similar at loweracceleration rates. However, the blurring artefacts in IST based recovered images are moreprominent at higher acceleration rates as compared to the proposed algorithm as indicatedby the white arrows.Computation 2023, 11, 7 19 of 21 Computation 2023, 10, x. https://doi.org/10.3390/xxxxx www.mdpi.com/journal/computation Figure 16. In real vivo data (a) Compares proposed method at bottom row with IST at top row withacceleration rate of 2. The arrow in (a) shows very minute artifacts. (b) depicts the performanceof proposed method at acceleration rate equal to 4. The arrow in (b) depicts the presence of arte-facts (c) shows the results of both algorithms with acceleration rate set at 8. The artefacts due tosubsampling become gradually more visible in IST results as highlighted by arrow mark (d) depictsthe results when acceleration rate is set at 12. Both techniques depicts the artefacts, however theseartefacts are visible in the IST as mentioned by the white arrow in the figure (e) show very muchdegraded image quality of IST, while comparing it with the proposed method. The subsamplingartefacts dominate the traditional IST result when acceleration rate R is set at 20.5. ConclusionsIn this paper, the novel CS recovery algorithm is proposed for compressively sampledsparse signals and biomedical images. The proposed method is applied in the 1-D sparsesignal, 2-D real human brain MRI and cardiac cine MRI. In our proposed algorithm, wehave introduced a hyperbolic tangent smooth approximation of non-differentiable l1-normand shrinkage. The experimental results quantitative analysis based on SSIM, PSNR,RMSE of recovered sparse signal and MR images have outperformed the conventional ISTalgorithm. The qualitative observations show significant improvement in the proposedmethod, especially at higher acceleration rates on Cine Cardiac MR images. In future, thisresearch work can be further enhanced to incorporate machine learning techniques usingComputation 2023, 11, 7 20 of 21large compressively sampled MRI datasets to restore accurate images, and the proposedmethod can be implemented on CS MRI scanners to reduce patient anxiety.Author Contributions: Conceptualization, H.H. and J.A.S.; methodology; software, H.H.; valida-tion, K.K., N.K. and J.A.S.; formal analysis, H.H.; investigation, H.H. and J.A.S.; resources, J.A.S.;writing—original draft preparation, H.H.; writing—review and editing, J.A.S. and N.K.; visualization,H.H. and J.A.S.; supervision, J.A.S.; project administration, J.A.S. and K.K.; funding acquisition, K.K.All authors have read and agreed to the published version of the manuscript.Funding: This research was partially funded by British Malaysia Intitute, Universiti of Kaul Lumpur,including publication of this research article.Data Availability Statement: 1-D dataset was generated using MATLAB. 2-D MR image was ac-quired from St. Mary’s Hospital, London, England. While breath-held under-sampled cardiac cinedataset is generated using MRXCAT software.Acknowledgments: We acknowledge the support given by University of Kaula Lumpur, BritishMalaysian Institute for funding the research publications. We also acknowledge Higher EducationCommission Pakistan for providing support in terms of hardware for experiments.Conflicts of Interest: The authors declare no conflict of interest.References1. Donoho, D.L. Compressed sensing. IEEE Trans. Inf. Theory 2006, 52, 1289–1306. [CrossRef]2. Candes, E.J.; Wakin, M.B. An Introduction To Compressive Sampling. IEEE Signal Process. Mag. 2008, 25, 21–30. [CrossRef]3. Candes, E.J.; Tao, T. Near-Optimal Signal Recovery From Random Projections: Universal Encoding Strategies? IEEE Trans. Inf.Theory 2006, 52, 5406–5425. [CrossRef]4. Lustig, M.; Donoho, D.; Pauly, J.M. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magn. Reson. Med.2007, 58, 1182–1195. [CrossRef]5. Lustig, M.; Donoho, D.L.; Santos, J.M.; Pauly, J.M. Compressed Sensing MRI. IEEE Signal Process. Mag. 2008, 25, 72–82. [CrossRef]6. Auger, D.A.; Bilchick, K.C.; Gonzalez, J.A.; Cui, S.X.; Holmes, J.W.; Kramer, C.M.; Salerno, M.; Epstein, F.H. Imaging left-ventricular mechanical activation in heart failure patients using cine DENSE MRI: Validation and implications for cardiacresynchronization therapy. J. Magn. Reson. Imaging 2017, 46, 887–896. [CrossRef]7. Liu, J.; Feng, L.; Shen, H.-W.; Zhu, C.; Wang, Y.; Mukai, K.; Brooks, G.C.; Ordovas, K.; Saloner, D. Highly-accelerated self-gatedfree-breathing 3D cardiac cine MRI: Validation in assessment of left ventricular function. Magn. Reson. Mater. Phys. Biol. Med.2017, 30, 337–346. [CrossRef]8. van Amerom, J.F.P.; Lloyd, D.F.A.; Price, A.N.; Kuklisova Murgasova, M.; Aljabar, P.; Malik, S.J.; Lohezic, M.; Rutherford, M.A.;Pushparajah, K.; Razavi, R.; et al. Fetal cardiac cine imaging using highly accelerated dynamic MRI with retrospective motioncorrection and outlier rejection. Magn. Reson. Med. 2018, 79, 327–338. [CrossRef]9. Paulus, W.J.; Tschöpe, C.; Sanderson, J.E.; Rusconi, C.; Flachskampf, F.A.; Rademakers, F.E.; Marino, P.; Smiseth, O.A.; DeKeulenaer, G.; Leite-Moreira, A.F.; et al. How to diagnose diastolic heart failure: A consensus statement on the diagnosis of heartfailure with normal left ventricular ejection fraction by the Heart Failure and Echocardiography Associations of the EuropeanSociety of Cardiology. Eur. Heart J. 2007, 28, 2539–2550. [CrossRef]10. Yerly, J.; Gubian, D.; Knebel, J.-F.; Schenk, A.; Chaptinel, J.; Ginami, G.; Stuber, M. A phantom study to determine the theoreticalaccuracy and precision of radial MRI to measure cross-sectional area differences for the application of coronary endothelialfunction assessment. Magn. Reson. Med. 2018, 79, 108–120. [CrossRef]11. Ahmed, A.H.; Qureshi, I.M.; Shah, J.A.; Zaheer, M. Motion correction based reconstruction method for compressively sampledcardiac MR imaging. Magn. Reson. Imaging 2017, 36, 159–166. [CrossRef]12. Gamper, U.; Boesiger, P.; Kozerke, S. Compressed sensing in dynamic MRI. Magn. Reson. Med. 2008, 59, 365–373. [CrossRef][PubMed]13. Donoho, D.L. De-noising by soft-thresholding. IEEE Trans. Inf. Theory 1995, 41, 613–627. [CrossRef]14. Blumensath, T.; Davies, M.E. Iterative hard thresholding for compressed sensing. Appl. Comput. Harmon. Anal. 2009, 27, 265–274.[CrossRef]15. Tropp, J.A.; Wakin, M.B.; Duarte, M.F.; Baron, D.; Baraniuk, R.G. Random Filters for Compressive Sampling and Reconstruction.In Proceedings of the 2006 IEEE International Conference on Acoustics Speech and Signal Processing Proceedings, Toulouse,France, 14–19 May 2006.16. Yin, W.; Osher, S.; Goldfarb, D.; Darbon, J. Bregman Iterative Algorithms for $\ell_1$-Minimization with Applications toCompressed Sensing. SIAM J. Imaging Sci. 2008, 1, 143–168. [CrossRef]17. Khajehnejad, M.A.; Xu, W.; Avestimehr, A.S.; Hassibi, B. Weighted `<inf>1</inf>minimization for sparse recovery with priorinformation. In Proceedings of the 2009 IEEE International Symposium on Information Theory, Seoul, Korea, 28 June–3 July 2009;pp. 483–487.Computation 2023, 11, 7 21 of 2118. Shah, J.; Qureshi, I.; Omer, H.; Khaliq, A. A modified POCS-based reconstruction method for compressively sampled MR imaging.Int. J. Imaging Syst. Technol. 2014, 24, 203–207. [CrossRef]19. Shah, J.A.; Qureshi, I.M.; Omer, H.; Khaliq, A.A.; Deng, Y. Compressively sampled magnetic resonance image reconstructionusing separable surrogate functional method. Concepts Magn. Reson. Part A 2014, 43, 157–165. [CrossRef]20. Bilal, M.; Shah, J.A.; Qureshi, I.M.; Kadir, K. Respiratory Motion Correction for Compressively Sampled Free Breathing CardiacMRI Using Smooth l(1)-Norm Approximation. Int. J. Biomed. Imaging 2018, 2018, 7803067. [CrossRef]21. Shah, J.; Qureshi, I.M.; Deng, Y.; Kadir, K. Reconstruction of Sparse Signals and Compressively Sampled Images Based on Smoothl1-Norm Approximation. J. Signal Process. Syst. 2017, 88, 333–344. [CrossRef]22. He, C.; Xing, J.; Li, J.; Yang, Q.; Wang, R. A New Wavelet Thresholding Function Based on Hyperbolic Tangent Function. Math.Probl. Eng. 2015, 2015, 528656. [CrossRef]23. Lu, J.-y.; Lin, H.; Ye, D.; Zhang, Y.-s. A New Wavelet Threshold Function and Denoising Application. Math. Probl. Eng. 2016,2016, 3195492. [CrossRef]24. Schmidt, M.; Fung, G.; Rosales, R. Fast Optimization Methods for L1 Regularization: A Comparative Study and Two NewApproaches. In Proceedings of the Machine Learning: ECML 2007, Berlin/Heidelberg, Germany, 17–21 September 2007; pp.286–297.25. Chang, S.G.; Bin, Y.; Vetterli, M. Spatially adaptive wavelet thresholding with context modeling for image denoising. IEEE Trans.Image Process. 2000, 9, 1522–1531. [CrossRef]26. Sendur, L.; Selesnick, I.W. Bivariate shrinkage functions for wavelet-based denoising exploiting interscale dependency. IEEETrans. Signal Process. 2002, 50, 2744–2756. [CrossRef]27. Nasri, M.; Nezamabadi-Pour, H. Image denoising in the wavelet domain using a new adaptive thresholding function. Neurocom-puting 2009, 72, 1012–1025. [CrossRef]28. Prinosil, J.; Smekal, Z.; Bartusek, K. Wavelet Thresholding Techniques in MRI Domain. In Proceedings of the 2010 InternationalConference on Biosciences, Online, 7–13 March 2010; pp. 58–63.29. Jin, J.; Yang, B.; Liang, K.; Wang, X. General image denoising framework based on compressive sensing theory. Comput. Graph.2014, 38, 382–391. [CrossRef]Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individualauthor(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury topeople or property resulting from any ideas, methods, instructions or products referred to in the content.
Article
eng
Reviewed
Vancouver : University of British Columbia Library
Multidisciplinary Digital Publishing Institute
CC BY 4.0
Faculty
Researcher
compressed sensing
holography cardiac cine MRI
l1-norm smooth approximations
hyperbolic tangent function
soft thresholding
Sparse Reconstruction Using Hyperbolic Tangent as Smooth l1-Norm Approximation
Text