Boxin Zuo

Meixia Geng

Xiangyun Hu

Earthquake Science

English

Crossref

RESEARCH PAPERMultiscale anisotropic diffusion for ringing artifact suppressionin geophysical deconvolution dataBoxin Zuo . Xiangyun Hu . Meixia GengReceived: 21 September 2015 / Accepted: 6 July 2016 / Published online: 19 August 2016 The Author(s) 2016. This article is published with open access at Springerlink.comAbstract Ringing artifact degradations always appear inthe deconvolution of geophysical data. To address thisproblem, we propose a postprocessing approach to suppressringing artifacts that uses a novel anisotropic diffusionbased on a stationary wavelet transform (SWT) algorithm.In this paper, we discuss the ringing artifact suppressionproblem and analyze the characteristics of the deconvolu-tion ringing artifact. The deconvolution data containingringing artifacts are decomposed into different SWT sub-bands for analysis, and a new multiscale adaptive aniso-tropic filter is developed to suppress these degradations.Finally, we demonstrate the performance of the proposedmethod and describe the experiments in detail.Keywords Deconvolution  Ringing artifacts  Anisotropicdiffusion  Stationary wavelet transform algorithm Multiscale1 IntroductionDeconvolution algorithms have been used to process datafrom a wide range of geophysical techniques, including forgravity, magnetic, and seismic data. Zuo and Hu (2012)also introduced a deconvolution algorithm for processingmagnetotelluric data. However, the ringing artifacts pro-duced by these deconvolution algorithms are rarely repor-ted. Most of the research on ringing artifacts focuses on theenhancement of image compression techniques (Peronaand Malik 1990; Weickert 1998; Pollak et al. 2000;Lysaker et al. 2003; Ding et al. 2005). Popovici andWithers (2007) proposed using a bilateral filter to suppressringing artifacts from compressed images. The discretewavelet transform algorithm was first used by Kim et al.(2009) to analyze ringing artifacts in multilevel edge mapsfor compressed images. These algorithms perform well forthe suppression of ringing artifacts in compressed images,but are not suitable for processing deconvolution data fromgeophysical measurements.The ringing artifacts that occur in the deconvolution ofgeophysical data are not the same as the artifacts in com-pressed images. In general, the ringing artifacts in decon-volution algorithms present as stronger oscillations than theartifacts in compressed images (Fig. 1).The dotted lines L1 and L3 represent the same dataselected from one row of observational data. Line L2 is theresult of compressing L1 using the JPEG 2000 codedalgorithm, while L4 is the result of deconvoluting L3.Geophysical deconvolution algorithms introduce ringingartifacts as a result of errors in the observational data, andthey commonly generate intense oscillations in areas ofdata discontinuity, e.g., edges and noise. They also oftenappear as false anomalous edges. The characteristics of theringing artifacts seen after the deconvolution of geophysi-cal data are very different from those of the ringing arti-facts in compressed images. Therefore, it is not possible todirectly suppress the ringing artifacts resulting from thedeconvolution of geophysical data using traditional diffu-sion algorithms for compression images. Based on thisanalysis, we propose a novel method to suppress ringingartifacts arising from the deconvolution of geophysicaldata.The discrete wavelet transform is a valuable signalanalysis tool and can highlight specific subband features(Arivazhagan et al. 2007) through the decomposition ofB. Zuo  X. Hu (&)  M. GengChina University of Geosciences, Wuhan 430074, Chinae-mail: xyhu@163.com123Earthq Sci (2016) 29(4):215–220DOI 10.1007/s11589-016-0161-4multiple scalars in complex data. Discrete wavelet trans-form is commonly applied for edge detections and con-trast enhancements in various applications (Huang et al.2002; Nakashizuka et al. 2004; Hatami et al. 2005).Pesquet et al. (1996) suggested using the stationarywavelet transform (SWT) algorithm to make the subbandcoefficients time invariant and to improve its power. Inthe SWT method, the coefficients have the same length asthe original data at every level. It uses a stationarytransform instead of traditional wavelet downsampling byadding an upsampling procedure on every scale coeffi-cient (Gao et al. 2008). The redundancy of the SWTalgorithm facilitates the detection of weak anomalies inthe data, especially the detection of the false edges ofringing artifacts.A nonlinear anisotropic diffusion method was proposedby Perona and Malik (1990) and has been widely studied ina number of scientific fields over the past two decades(Black et al. 1998; Pollak et al. 2000; Chen and Hsu 2005;Ding et al. 2005). The advantage of this method is that itcan adaptively smooth noise while simultaneously pre-serving the actual model boundaries. Therefore, we used itto analyze ringing artifacts in the decomposed subbands ofthe SWT algorithm. In summary, we propose an aniso-tropic diffusion method to suppress the ringing artifacts inthe results of the deconvolution of geophysical data basedon SWT analysis.The background to the SWT algorithm is brieflyreviewed and then the proposed ringing artifact anisotropicfilter based on the SWT algorithm is discussed. Wedescribe our proposed multiresolution system and thenpresent our experimental results and discussion beforeproviding conclusions on our work.2 Background to the SWT algorithmAssuming H and G are the coefficients of the wavelet filters,then H[j](H[0] = H, j = 1, 2,…, n) and G[j](G[0] = G, j = 1,2, …, n) can be obtained by inserting a zero between everyadjacent pair of elements. If a0 is an initial signal sequence,the signal stationary wavelet (Pesquet et al. 1996) isajþ1 ¼ H½jajbjþ1 ¼ G½jbj ;ð1Þwhere aj?1 represents the approximation coefficient of theSWT algorithm and bj?1 is the detail. We used aj(e1,…, ej)and bj(e1, …, ej) to denote the coefficients obtained fromthe SWT algorithm at level j/jðajþ1; bjþ1Þ ¼ Hajþ1ðejþ1Þ þ Gbjþ1ðejþ1Þ; ðejþ1 ¼ 0Þ;ð2Þujðajþ1; bjþ1Þ ¼ Hajþ1ðejþ1Þ þ Gbjþ1ðejþ1Þ; ðejþ1 ¼ 1Þ;ð3Þwhere * represents a conjugate operator. The inverse SWT isajðejÞ ¼ 12½/jðajþ1; bjþ1Þ þ ujðajþ1; bjþ1Þ: ð4Þ3 Ringing artifact anisotropic diffusion filter basedon the SWT algorithmThe basic mechanism behind our approach can be analyzedby comparing the detailed SWT coefficients of thedeconvolution results with the observed ringing artifacts(Fig. 2). LL2 represents the output from the low-pass filterin the horizontal and vertical directions, whereas HLi areFig. 1 a Ringing artifacts in a compressed image (Marziliano and Dufaux 2004). b Ringing artifacts produced during deconvolution ofgeophysical data (Zuo et al. 2009)216 Earthq Sci (2016) 29(4):215–220123the outputs from the high-pass filter in the vertical andhorizontal directions. All the high-frequency waveletcoefficients are magnified by a coefficient of 2 in this paperso that they can be more easily seen.In Fig. 2d, the low-frequency wavelet coefficient LL2does not appear in any of the ringing artifacts, so it does notneed to be filtered. However, Fig. 2e, f shows high levelsof ringing artifacts in the decomposed high-frequencycoefficients. In addition, the contrast in these coefficientsincreases as the level of decomposition k increases. A newanisotropic filter is suggested by this analysis, which workson the decomposed SWT coefficients as follows.For a high-frequency wavelet coefficient Hk(x, y) decom-posed by the SWT algorithm and with a finite power and nozero values over the data support value X, k is the level ofdecomposition of the SWT algorithm k (k = 0, 1, 2,…, n).The output Hk(x, y, t) can be expressed according to thefollowing partial differential equation:oot Hkðx; y; tÞ ¼ r  ½cðqÞrHkðx; y; tÞHkðx; y; 0Þ ¼ Hkðx; yÞ; ðoHkðx; y; tÞ=on~ÞjoX ¼ 0;ð5Þwhere qX denotes the border of X and on~ is the outernormal of qX. A new diffusivity function c(q, k) is definedin Eq. (6):cðq; kÞ ¼ ð1þﬃﬃkpÞ=f1þ ½q2ðx; y; t; kÞ  q20ðt; kÞ=½1þ q20ðt; kÞÞg; ð6Þwhere q(x, y, t, k) is the instantaneous coefficientdetermined byqðx; y; t; kÞ ¼ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃð1=2ÞðjrHk=HkjÞ2  ð1=16Þðr2Hk=HkÞ2½1þ ð1=4Þðr2Hk=HkÞ2s:ð7Þq(x, y, t, k) serves as the boundary detector for the data. Theringing artifact scale function q0(t, k) is defined to controlthe degree of smoothing behavior as defined in Eq. (8):q0ðt; kÞ ¼ ðvarðHkðx; y; tÞÞÞ1=2=meanðHkðx; y; tÞÞ; ð8Þwhere varðÞ and meanðÞ are the mean of the intensityvariance and the mean of the directions of the variants inthe high-frequency wavelet coefficients decomposed bythe SWT algorithm over a ringing artifact area,respectively.The proposed anisotropic filter has some characteristicsthat are not shared with conventional anisotropic diffusion.In previous research on anisotropic filters, c(q) was equal toa Gaussian filter Kr as Kr ¼ 12pr2 expð jxj2þjyj22r2 Þ (WeickertFig. 2 a Actual model. b Observation data. c Deconvolution results showing ringing artifacts. d Low-frequency wavelet coefficients LL2.e High-frequency wavelet coefficients HL1. f High-frequency wavelet coefficients HL2Earthq Sci (2016) 29(4):215–220 2171231998) and q0() was defined as a domain function (Yu andActon 2002). In this paper, Eqs. (6) and (8) define c(q, k),q(x, y, t, k), and q0(t, k) in the wavelet domain without theassumption of a Gaussian distribution. Moreover, c(q, k),q(x, y, t, k), and q0(t, k) are related to the SWT level ofdecomposition k.As shown in Fig. 2, as the SWT level of decomposi-tion k is increased, the contrast in the ringing artifactincreases dramatically. The function c(q, k) then presentsa high value and the proposed filter shows higher diffu-sion behavior as the SWT level of decomposition k isincreased.For a fixed level k, according to Eq. (6), the filter showsencouraging anisotropic diffusion in high-contrast edgeregions where q(x, y, t, k) is larger than q0(t, k). Otherwise,it acts as an isotropic filter over the ringing and non-edgeregions.Using the proposed anisotropic filter, the filter results forthe decomposition of the high-frequency wavelet coeffi-cients in Fig. 2d, f by the SWT algorithm are shown inFig. 3.4 Multiresolution anisotropic diffusionThe proposed multiresolution ringing artifact anisotropicdiffusion method is shown as a flow chart in Fig. 4. First,the deconvolution result is decomposed by the SWTalgorithm, then the decomposed coefficients in everydetailed subband are filtered by the proposed anisotropicfilter. Finally, the inverse SWT algorithm is applied to thefiltered results and the final result is the output.5 Results and discussionTo test the proposed anisotropic diffusion algorithm, thealgorithm was applied to the synthetic model shown inFig. 2a. The model was fitted by a Gaussian kernelK ¼ eðs21þs22Þ=ð2r2Þ=PP eðs21þs22Þ=ð2r2Þwith a support sizes1 = s2 = 6, r = 2. The observation data are shown inFig. 2b. A standard deconvolution algorithm was applied tothe data with various values of point spread function sup-port size bias (2, 4, 6, 8). The results of using our algorithmfor deconvolution and diffusion are shown in Fig. 5.As Fig. 5a shows, the results of the proposed algorithmgive a slightly low-pass effect. This phenomenon is due tothe smoothing behavior of anisotropic diffusion. It is alsoclear that the discontinuities in the data are preservedduring the diffusion procedure. Although the numericalrange of the geophysical data and the common images arenot the same, the characteristics of these ringing artifacts inthese data are similar. Commonly, the textures of imagesare very complex, and its deconvolution result will containcomplicated ringing artifact. So we discuss the suppressionof the ringing artifact in images in the experiment sectionas the example. The geophysical data processing is similarto this procedure.The effectiveness of the algorithm in suppressing ring-ing artifacts was evaluated by the mean absolute error(MAE) relative to the original model (Fig. 3a), as definedin Eq. (5):MAE ¼ 1MNXmi¼1Xmj¼1jIði; j; tÞ  Iorij: ð5ÞThe results of the proposed SWT method based onmultiresolution anisotropic diffusion are shown in Table 1.By comparing the results with those of the bilateral filter, itis easy to see that the numerical precision of the proposedSWT-based anisotropic diffusion results is higher than thatfor the bilateral filter, particularly when the ringing arti-facts are strong. The proposed diffusion results also have abetter visual quality. The bilateral filter results (Fig. 5c)seem to show fewer ringing artifacts and less blur than theproposed diffusion result (Fig. 5b). Based on the numericalstatistical results in Table 1, it is easy to see that the dataobtained using conventional anisotropic diffusion (Fig. 5d)are oversmoothed during filtering.The conventional diffusion and the bilateral filter did notpreserve the weak edges. This behavior leads to a decreaseFig. 3 a Filtering result for Fig. 2e. b Filtering result for Fig. 2fRinging artifact data SWTAnisotropicfilterInverseSWTOutputdataFig. 4 Flow chart for SWT anisotropic diffusion treatment of ringing artifacts218 Earthq Sci (2016) 29(4):215–220123in the actual model details, and some false edges appear inthe results. In contrast, the proposed method avoids thisbehavior. Under the conditions of strong ringing artifacts(e.g., Fig. 5a, 4th image), the proposed anisotropic diffu-sion filter gives a better perceptual quality and numericalprecision than the other filters. The MATLAB (version2014b)-based programs were executed using an Intel i5-2500 k CPU, 3.30 GHz, 8.0 GB RAM system and theexecution time for a 256 9 256 size model was\1 s.6 ConclusionsThis paper has introduced a diffusion algorithm that sup-presses ringing artifacts for use in the postprocessing ofdeconvolution data from geophysical applications. Theproposed anisotropic diffusion filter was able to simulta-neously suppress the ringing artifacts and preserve themodel boundaries. Experiments showed that this algorithmsignificantly enhanced both the visual and numerical staticFig. 5 Diffusion results for deconvolution. Comparison of results obtained by the proposed algorithm with those from other related algorithms.a Deconvolution results obtained with point spread function support size estimated bias values of 2, 4, 6, and 8, respectively. b Diffusion resultsobtained by the proposed algorithm. c Results from a filtered by a bilateral filter. d Results from a filtered by conventional anisotropic diffusionEarthq Sci (2016) 29(4):215–220 219123performance of the deconvolution results. This algorithmalso worked well for the suppression of high-level ringingartifacts.Open Access This article is distributed under the terms of theCreative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use,distribution, and reproduction in any medium, provided you giveappropriate credit to the original author(s) and the source, provide alink to the Creative Commons license, and indicate if changes weremade.ReferencesArivazhagan S, Deivalakshmi S, Kannan K, Gajbhiye BN, MuralidharC, Lukose SN, Subramanian MP (2007) Multi-resolution systemfor artifact removal and edge enhancement in computerizedtomography images. Pattern Recogn Lett 28:1769–1780Black JM, Sapiro G, Marimont DH, Heeger D (1998) Robustanisotropic diffusion. IEEE Trans Image Process 7(3):421–432Chen B, Hsu EW (2005) Noise removal in magnetic resonancediffusion tensor imaging. Magn Reson Med 54:393–407Ding Z, Gore JG, Anderson AW (2005) Reduction of noise indiffusion tensor images using anisotropic smoothing. MagnReson Med 53:485–490Gao Q, Zhao Y, Lu Y (2008) Despeckling SAR images usingstationary wavelet transform combining with directional filterbanks. Appl Math Comput 205:517–524Hatami S, Hosseini R, Kamarei M, Ahmadi H (2005) Wavelet basedfingerprint image enhancement. In: IEEE international sympo-sium on circuits systems, vol 5, Kobe, Japan, May 2005,pp 4610–4613Huang MY, Tseng DC, Liu MSC (2002) Wavelet image enhancementbased on Teager energy operator. In: 16th internationalconference on pattern recognition, vol 2, Quebec, Canada,August 11–15, pp 993–996Kim S, Jun S, Lee E, Shin J, Paik J (2009) Ringing artifact removal indigital restored images using multi-resolution edge map. Int JSignal Process Image Process Pattern Recognit 2(4):95–106Lysaker M, Lundervold A, Tai X-C (2003) Noise removal usingfourth-order partial differential equation with applications tomedical magnetic resonance images in space and time. IEEETrans Image Process 12(12):1579–1590Marziliano P, Dufaux F (2004) Perceptual blur and ringing metrics:application to JPEG2000. Signal Process Image Commun19:163–172Nakashizuka M, Aoki K, Nitta T (2004) A simple edge-weightedimage enhancement filter using wavelet scale products. In: 47thIEEE international midwest symposium on circuits and systems,Hiroshima, Japan, July 25–28, pp 285–288Perona P, Malik J (1990) Scale-space and edge detection usinganisotropic diffusion. IEEE Trans Pattern Anal Mach Intell12(7):629–639Pesquet JC, Krim H, Carfantan H (1996) Time-invariant orthonormalwavelet representations. IEEE Trans Signal Process 44:1964–1970Pollak I, Willsky AS, Krim H (2000) Image segmentation and edgeenhancement with stabilized inverse diffusion equations. IEEETrans Signal Process 9(2):256–266Popovici I, Withers WD (2007) Locating edges and removing ringingartifacts in JPEG images by frequency-domain analysis. IEEETrans Image Process 16:1470–1474Weickert J (1998) Anisotropic diffusion in image processing, vol 2.Teubner, Stuttgart, pp 184–189Yu YJ, Acton ST (2002) Speckle reducing anisotropic diffusion.IEEE Trans Image Process 11(11):1260–1270Zuo BX, Hu XY (2012) Geophysical model enhancement techniquebased on blind deconvolution. Comput Geosci 49:170–181Zuo BX, Ming DL, Tian JW (2009) Perceptual ringing metric toevaluate the quality of images restored using blind deconvolu-tion algorithms. Opt Eng 48(3):1–9Table 1 Experimental and numerical statistical resultsFigure 3a Figure 5a(1st–4th)Figure 5b (1st–4th),proposed methodFigure 5c (1st–4th),bilateral filterFigure 5d (1st–4th),conventional anisotropicfilterMean absolute error 0.0735 0.0522 0.0612 0.05780.1549 0.0698 0.1245 0.08850.2219 0.0749 0.1923 0.12160.1172 0.0382 0.1063 0.0633Increasing precision (%) N/A 28.97 16.73 21.36N/A 54.93 19.63 42.87N/A 66.24 13.34 45.20N/A 67.41 9.30 45.99Peak signal-to-noise ratio 24.7306 28.5947 25.6992 26.544519.5355 26.0079 20.3793 22.525316.9621 25.3502 17.5304 20.071915.4078 23.3484 15.8285 18.6624Structural similarity 0.6549 0.7073 0.6737 0.67190.3569 0.5804 0.4847 0.56390.2272 0.5516 0.3118 0.45780.1544 0.5185 0.2046 0.3816220 Earthq Sci (2016) 29(4):215–220123

Multiscale anisotropic diffusion for ringing artifact suppression in geophysical deconvolution data

https://link.springer.com/content/pdf/10.1007%2Fs11589-016-0161-4.pdf

Multiscale anisotropic diffusion for ringing artifact suppression in geophysical deconvolution data

Abstract

Similar works

Full text

Available Versions

Crossref