Compressed Sensing for STM imaging of defects and disorder
Brian E. Lerner, Anayeli Flores-Garibay, Benjamin J. Lawrie, Petro Maksymovych
CCompressed Sensing for STM imaging of defects and disorder ∗ Brian E. Lerner, Anayeli Flores-Garibay, Benjamin J. Lawrie, and Petro Maksymovych † Oak Ridge National Laboratory, 1 Bethel Valley Rd, Oak Ridge, TN 37831 (Dated: January 19, 2021)Compressed sensing (CS) is a valuable technique for reconstructing measurements in numerousdomains. CS has not yet gained widespread adoption in scanning tunneling microscopy (STM),despite potentially offering the advantages of lower acquisition time and enhanced tolerance tonoise. Here we applied a simple CS framework, using a weighted iterative thresholding algorithmfor CS reconstruction, to representative high-resolution STM images of superconducting surfacesand adsorbed molecules. We calculated reconstruction diagrams for a range of scanning patterns,sampling densities, and noise intensities, evaluating reconstruction quality for the whole imageand chosen defects. Overall we find that typical STM images can be satisfactorily reconstructeddown to 30% sampling - already a strong improvement. We furthermore outline limitations of thismethod, such as sampling pattern artifacts, which become particularly pronounced for images withintrinsic long-range disorder, and propose ways to mitigate some of them. Finally we investigatecompressibility of STM images as a measure of intrinsic noise in the image and a precursor to CSreconstruction, enabling a priori estimation of the effectiveness of CS reconstruction with minimalcomputational cost.
Keywords : Compressed sensing, scanning tunneling mi-croscopy
I. INTRODUCTION
Scanning tunneling microscopy (STM) and spec-troscopy (STS) have become indispensable techniques forelectronic, structural and magnetic characterization ofsurfaces with atomic resolution. STM has enabled inves-tigations of broken symmetry and vortex interactions insuperconductors [1, 2], enabled the band structure map-ping of quantum materials [3], and was used for the firstobservations of spatial LDOS modulations [4, 5].However, small tunneling currents limit the rate of cur-rent measurement to the millisecond timescale, so thatSTM measurements are characterized by comparativelylong measurement times [3]. This limitation becomes ap-parent in experiments that seek to probe extended sur-face areas, seek rare events such as low density defects,and want to strike a balance between high-resolutionmeasurement in real space and energy resolution. In suchcases, the ability to accurately reconstruct the underlyingperiodic and defect structure of nanoscale samples withreduced measurement time is highly desirable. ∗ This manuscript has been authored by UT-Battelle, LLC un-der Contract No. DE-AC05-00OR22725 with the U.S. Depart-ment of Energy. The United States Government retains and thepublisher, by accepting the article for publication, acknowledgesthat the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce thepublished form of this manuscript, or allow others to do so, forUnited States Government purposes. The Department of En-ergy will provide public access to these results of federally spon-sored research in accordance with the DOE Public Access Plan(http://energy.gov/downloads/doe-public-access-plan). † [email protected] Compressed sensing (CS) shows potential for meetingthis demand. CS is based on the notion that if a basis setcan be found where the signal is sparse (and as a corollarythe signal is compressible in that basis), accurate recon-struction is possible using fewer measurements than re-quired by the Shannon-Nyquist Sampling Theorem. CShas been successfully employed for diverse applicationsincluding radio interferometry [6], nuclear magnetic res-onance of protein structure [7, 8], recovery of correlationsof entangled photon pairs [9, 10], medical imaging [11, 12]and many more.An image is compressible by virtue of its sparsity ina transform domain. Most images in the natural worldhave a sparse frequency or wavelet representation, in-cluding those generated by scanning microscopies. In-deed, CS has been successfully implemented in scanningelectron [13], atomic force [14], and piezoresponse forcemicroscopy [15], and quasiparticle interference imagingby STS [3, 16]. However, a detailed understanding of thepotential of CS for STM has yet to be developed, partic-ularly with respect to imaging defects and disorder.In this paper, we explore the parameter space of asimple CS framework in the context of representativeSTM images from surfaces of superconductors and singlemolecule layers (introduced in section II). Our specific fo-cus is to emphasize the quality of reconstruction arounddefects and as a function of added noise. In sections IIIand IV, the basic methodology of CS is laid out, andthe framework is described. Using a soft weighted itera-tive thresholding (SWIT) algorithm of practical compu-tational complexity, we performed reconstructions acrossvariable noise perturbation intensities and sampling den-sities. These reconstructions are evaluated for structuralsimilarity index measure (SSIM) and mean squared error(MSE) and are used to calculate reconstruction diagramsin section V. Our results reveal that accurate reconstruc-tion can be obtained at sampling densities as low as 20-30% for images with both point and extended nanoscale a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n defects - i.e. with almost 5-fold compression. We alsonote artifacts arising in the reconstructions, and detailways of mitigating these deviations through proper algo-rithm configuration. To effectively apply CS in practice,it is very helpful to understand what types of images canbe effectively reconstructed. In V, we also characterizeour images using compressibility, finding compressibilityto be an effective measure of noise in the STM images,and a necessary, albeit not sufficient, criterion for effec-tive CS reconstruction. II. EXPERIMENTAL DATA
We applied CS to representative STM images of acleaved 100-surface of FeSe superconductor with Se va-cancy defects [17] (Fig. 1a) and two kinds of adsorbedmolecular layers - C on Ag(111) (Fig. 1c) and TCNQ(tetracyanoquinodimethane) on graphite (Fig. 1b). Eachof the sample images have a different size, lattice struc-ture, and point or extended defect. Moreover, as seenin Fig. 1d, the images represent three kinds of inten-sity distribution, centered on low values correspondingto the atomic lattice in the case of FeSe, a broader andmore uniform distribution in the case of TCNQ and adistinctly bimodal distribution for C , owing to a singleatomic step of the underlying substrate. III. CS BASICS
Sparsity regularization is a common approach to im-pose constraints on undefined optimization problems [18],which gave rise to CS methodology in the mid-2000s[19, 20]. CS is designed to reconstruct a signal x ∈ R n × from samples y ∈ R m × , where typically m (cid:28) n . Suc-cessful reconstruction is possible when x has a sparserepresentation α ∈ R n × , i.e. in some basis the num-ber of significant coefficients k in α is small compared to n . The CS algorithm computes α . Once obtained, x isrecovered using the basis transform Ψ ∈ R n × n : x = Ψ α (1)The sampling process has a matrix representation Φ ∈ R m × n constructed by stacking each measurement vector:Φ x = y (2)Substituting eq. 1 for x in eq. 2 and setting A = ΦΨ wehave: Aα = y (3)CS provides a solution α for this undetermined systemof equations by minimizing the sparsity of α under theconstraints of eq. 3, expressed as: min (cid:107) α (cid:107) (cid:96) s.t. Aα = y (4) FIG. 1. STM images of FeSe (a), TCNQ (b), and C (c),with representative defects magnified in each inset. (d) Thedistribution of normalized constant-current STM height foreach image. While this provides an exact solution, (cid:96) minimizationis a combinatorial optimization problem that is compu-tationally expensive, and intractably so for large signals[20]. Fortunately, the (cid:96) norm can be substituted to con-vert the problem into one of convex optimization, wherefor most inputs, α is recovered exactly [20]. IV. FRAMEWORK
The CS framework can utilize a variety of 1) samplingmatrices Φ, 2) transform matrices Ψ, and 3) optimiza-tion algorithms. Ψ should necessarily be chosen to ensuresparsity in the transform domain, but it should also beincoherent with Φ. The algorithm minimizes the sparsityin α while remaining correlated to the measurements y (eq. 3). In our reconstructions, we use Lissajous and ro-tated line trajectories for sampling patterns, the discretecosine transform (DCT), and a SWIT algorithm. Theelements of this framework, with special regard to theirapplicability for STM, are discussed in the following. A. Transform Matrix
STM images often exhibit a large amount of order andare generally smooth (i.e. differentiable in the absence ofnoise). As a result, the images lend themselves to sparsityin the DCT basis. The DCT transform matrix also hasthe advantage of being maximally incoherent with pointsampling matrices [21], and has a fast matrix implemen-tation [22]. This transform has been utilized in previ-
FIG. 2. DCTs of (a) FeSe, (b) TCNQ, and (c) C . (d)The intensity of the diagonal coefficients for each DCT, aswell as the DCT of an array of random Gaussian noise, whichdemonstrate varying sparsity levels. ous applications of CS [23–25], and has historically beenused for JPEG compression [26]. The discrete wavelettransform (DWT) is another commonly used dictionaryin compressed sensing, thought it works most efficientlywith dense sampling matrices with random entries likethose used for single-pixel imaging and is less incoherentthan DCT for point sampling matrices [22]. B. Sampling Matrix
When scanning a surface, it is conventional to use araster scan, where the probe traverses the sample in aseries of alternating lines, resulting in an evenly sampledgrid. The speed of the probe and the sampling frequencyare set based on the demands of the experiment. Whilethe design of the sampling matrix Φ in other CS applica-tions is often flexible (programmable with a spatial lightmodulator for optical CS applications, for instance), weare constrained to sampling along the continuous pathof the probe. Here, since we are concerned with thealgorithmic aspects of the reconstruction, we chose touse pre-existing STM images and resample them withsmooth Lissajous (Fig. 3d) and rotated line (Fig. 3a)patterns which make the methods more compatible withfast scanning. The sampling can furthermore be ran-domized along the sampling path, but we have not seena significant impact from such randomization.
FIG. 3. The path of the rotated line pattern is shown in (a),with simulated start and end points denoted by green andred circles. Despite sparse sampling of the image (b), decentreconstruction is achieved (c). The same process is also shownfor Lissajous (d-f). Reconstructions in this figure performedfor 20% sampling density and 100 iterations.
C. Optimization Algorithm
There are a variety of reconstruction algorithms thathave already been explored for other CS applications. Inthe convex optimization class, the (cid:96) norm is replacedby the (cid:96) norm. Greedy pursuit algorithms use an iter-ative approach where locally optimal decisions are madein each iteration. Iterative thresholding [27] is a type ofgreedy pursuit algorithm that has relatively low compu-tational complexity and is robust to noise. Due to thesebenefits, we employed a SWIT algorithm as successfullydemonstrated in [14]. The algorithm works as follows: α = 0 r = y for i in I : c = A T rα = η wst ( α + κ · c ) r = y − Aα i f (cid:107) r (cid:107) (cid:96) < (cid:15) (cid:107) y (cid:107) (cid:96) : break Initialization to α = 0 can be changed to an educatedguess and the stopping condition can be arbitrarily cho-sen, while the step size κ ensures convergence. The softweighted thresholding function η wst is implemented as: η wst = 1 w sgn ( x )( | wx | − t ) , | wx | − t > , | wx | − t ≤ t is customiz-able. Here, we set a fixed value on the number of nonzerocoefficients while initializing the algorithm. In each iter-ation, the coefficients are weighted as described above, t is adjusted to maintain the specified sparsity, and coeffi- FIG. 4. Reconstructed images for ten, five, and two-fold un-dersampling for TCNQ (a-c), C (d-f), and FeSe (g-i), withmagnified defects in insets. All reconstructions performed for100 iterations using the rotated line sampling pattern. cients below t are zeroed. By tuning the weights to modelexpected DCT dispersion, weighted iterative threshold-ing algorithms tend to outperform their non-weightedcounterparts [14]. Each of the reconstructions consti-tuting the reconstruction diagrams ran for 100 iterationsdue to computational considerations, though in our ex-periment we found that reconstruction tends to improveup to around 300 iterations–and sometimes many more–before plateauing. D. Quality Assessment
To understand the bounds of reconstruction, we eval-uated the SWIT algorithm while systematically varyingthe noise intensity δ and sampling density ρ . While iter-ative thresholding algorithms are noted for being noise-robust [28], little investigation has been carried out toconfirm this for reconstruction of STM images. In orderto test this, we generated 1 /f noise in Python and appliedit to pixels along the simulated measurement path so asto mimic varying noise levels during measurement. Thenoise perturbation scale for each image was normalized torange from 0.1–1 of the highest-peak FWHM in the im-age’s intensity histogram (Fig. 1d). We used FWHM asa measure of spread due to the approximately Gaussianshape of the distributions. We implemented rotated lineand Lissajous sampling patterns across ρ from 0.02–0.5.The patterns used here were generated using magni [29], a compressed sensing Python package for atomic forcemicroscopy.For each reconstruction in this δ – ρ parameter space,the quality of the reconstructed image was evaluatedfor SSIM and MSE. SSIM was calculated using scikit-image’s [30] default implementation, which is adaptedfrom [31]. The MSE is derived in the standard way,1 N (cid:88) ( χ − x ) (7)where N is the number of pixels, x is the reconstructedimage and χ is the base image.To perform these reconstructions, we build our CSframework with a DCT transform, due to the benefitsespoused in Sec. IV, in combination with the noted sam-pling patterns. We solve eq. 3 using a SWIT algorithmas described in the code block in Sec. IV, with κ = 0 . y − Ax ) and the 2-norm of y is less than a tolerance (cid:15) = 0 . η wst used in these reconstructionsare adopted from a Gaussian model of DCT structurein [14], which was used to successfully reconstruct AFMimages. V. RESULTS
Our first observation is that CS is generally very goodat reconstructing STM images even at a sampling den-
FIG. 5. Noise perturbation intensity ( δ ) vs. sampling density( ρ ) phase diagrams for reconstructions of TCNQ (a), C (b)and FeSe (c), with relevant reconstructions shown above andbelow the diagram for each sample. The parameters used forthe reconstructions in the top row are marked by green dotsin the respective diagrams; the bottom row parameters aremarked by yellow dots. sity as low as 20% of the original image. To ascertainthat this conclusion applies not only to spatial order inthe images, but also to defect sites, we have identified thelatter using state-pace methods for detection of protru-sions (using Laplacian of Gaussian filter), and then builtlocal masks of the defects, comparing reconstruction inthat local region. As seen in the insets of Fig. 4, sin-gle vacancies in FeSe and extended defects in the TCNQoverlayer (missing molecules) reconstruct well. At 50%sampling density, the reconstructed defects are indistin-guishable from their unsampled counterparts.The δ – ρ reconstruction diagrams (Fig. 5) demonstratethe method’s robustness to moderate 1 /f noise. All re-constructions have high SSIM above sampling density ρ ≈
30% which only begins to degrade at noise pertur-bations of 0.4 for TCNQ and 0.8 for FeSe. While high-noise distortions are apparent in the reconstructions ofTCNQ (Fig. 5a) and FeSe (Fig. 5c), the simplicity ofFeSe’s vacancies and the regularity of its lattice likelylead to smoother SSIM falloff at high noise. C , in starkcontrast, has a wholly noise-independent transition (Fig.5e). C also exhibits a sharp transition to higher SSIMat sampling density around 30%, which exceeds the tran-sition point of the other samples by 10-20%. Visual ex-amination of the reconstructions (Fig. 5h) reveals thepresence of sampling pattern artifacts at low SSIM whichdisappear after the transition line. The reasons for thisdeviation will be discussed below.Given that CS is predicated on the principle of com-pression, we explored the extent to which our CS resultscorrelate to image compressibility for typical STM im-ages as well as simulated arrays, one composed of pseudo-random Gaussian noise (Fig. 6b) and the other an or-dered lattice (Fig. 6c). We evaluated compressibility bytransforming each image and kept a compressed set of themost significant coefficients, setting the rest to 0 beforetransforming back to real space and evaluating the MSE.The pseudo-random noise image displays the highest er-ror across compression sizes, i.e. it is the most incom-pressible, while the ordered lattice is most compressible.STM images fall between these two extremes, as seen inFig. 6a. Intriguingly, there is a very significant differencebetween individual images, which actually goes againstthe trend that may be inferred from the visual inspectionof the original data in Fig. 1. C , not TCNQ or FeSe, isthe most compressible image, while FeSe is notably lesscompressible than either TCNQ or C . The difference incompressibility stems from the signal to noise ratio thatcharacterizes these images. To ascertain that this is thecase, in Fig. 6d, we plot compressibility of TCNQ as afunction of strength of added noise (measured as a frac-tion of the largest signal in the image). The compress-ibility curve very clearly traverses the range observed inFig. 6a, eventually becoming equivalent to noise. Wenote that all these images were all acquired on differentdays, with different physical tips and different instrumentconditions. The ability to “calibrate” the STM imagewith compressibility appears to be a valuable measure of the data quality and experimental results.We now show that the compressibility of an image gen-erally correlates with its CS performance. In Fig. 6e weplot the normalized CS reconstruction error vs the nor-malized DCT compression error as a function of noise,for three levels of data compression. For 5-fold com-pression (20% sampling), the correlation is reasonablygood, which confirms our notion. However, for smallerdensities, CS systematically produces higher error thanobtained by DCT compression, which reduces the corre-lation between the two techniques. We speculate thatpartly these deviations are due to CS being sensitiveto the compatibility of sampling and transform matri-ces with both each other and the image, as well as thealgorithm type and configuration.A striking disparity, however, appears for C , whichis the most compressible of the typical STM images (Fig.6a) but requires the highest sampling density to achievequality reconstruction. Interestingly, past the transitionline in both SSIM (Fig. 5) and MSE (Fig. 9) phase dia-grams. C generally has the highest SSIM, followed byTCNQ then FeSe. Resolving this puzzle depends on anunderstanding of how and when sampling pattern arti-facts appear, as their presence is the major cause of ρ dependence in the phase transition. We have found thatthis brand of artifact can be removed by properly con-figuring the SWIT algorithm. Small disturbances can be FIG. 6. The STM images, along with random Gaussian noise(b) and an ordered lattice (c), were transformed into the DCTbasis before being compressed and inverse transformed. TheMSE–normalized against the highest value of each curve– vs.the normalized compressed size, i.e. the compression ratioin the DCT domain, is shown for each image in (a). Thisprocedure is repeated for different levels of Gaussian noiseapplied to TCNQ pre-transform (d). (e) Compression errorvs. CS reconstruction error as a function of noise for varyingsampling/compression ratio.
FIG. 7. C reconstructions for Lissajous (b,e) and rotatedline (c,f) sampling patterns using a 1% threshold on the num-ber of non-zero coefficients and 300 iterations. The top re-constructions utilized a wide-variance DCT weight model (a)which was also used for the reconstructions in Fig. 5. Thoseon the bottom utilized a model with a severely limited vari-ance; the relevant low-frequency corner of this model is shownin the inset of (a). The diagonal of each model and sampleDCT is compared in (d). removed by increasing the number of iterations, but moreprominent artifacts require increased iterations and/orspecialized setup of the threshold function (eq. 5).In each iteration of the SWIT, the threshold functionweights each coefficient using a DCT model and basedon a specified threshold ratio, keeps a certain number ofcoefficients while setting the rest to 0. We show that set-ting the threshold ratio to 1%, instead of 5%, runningfor 300 iterations, and minimizing the variance in theweight model, the artifacts can be removed from C .Reconstruction with the Lissajous pattern was more re-sponsive (Fig. 7b) to the same DCT-model variance (Fig.7a) used for the phase diagrams, though interestingly therotated line reconstructions improved (Fig. 7f) only withseverely minimized variance (Fig. 7a inset).To determine the ideal thresholding function param-eters, we evaluated C and TCNQ for SSIM across arange of threshold ratios and variances (Fig. 8). We seethat SSIM falls off for TCNQ at low threshold ratios forall variances σ , and in the limit of low σ and thresh-old ratio– a trend consistent for both sampling patterns.This behavior is expected as reducing threshold ratio anddecreasing σ are both tantamount to applying a low-passfilter. Surprisingly, the filtering at low σ and thresh-old ratio produces distinctly higher SSIM for the defectcompared to the global image, though visual inspectionrevels intense lattice warping. The defect diagrams forboth samples show higher SSIM for rotated line than Lis-sajous, a difference especially stark for C . In contrastto TCNQ, which has similar trends in performance forboth patterns, C is quite different. For Lissajous, theSSIM falls off at at threshold ratios around 20% inde-pendently of σ . Rotated line maintains high SSIM across low σ for all thresholds, though a transition line devel-ops with increasing σ that exponentially falls to very lowthreshold ratios. At low threshold ratio, C is seem-ingly immune from SSIM degradation, though the defectdiagram has a slight dip at very low threshold. Visualinspection of reconstructions in this regime reveals heavyand unsatisfactory smoothing which retains a semblanceof the step defect and an accordingly high SSIM. For allsamples and patterns in Fig. 8 though, overlapping high-SSIM regions across global and defect diagrams reveal anoptimal parameter space for defect-lattice reconstructionand provide a proof-of-principle for effectively tuning thethresholding function parameters.To better understand C ’s sensitivity to sampling pat-tern, we refer back to its DCT (Fig. 2). Each DCT coef-ficient distribution features a cluster of high-magnitudecoefficients in the upper left-hand corner, i.e. for low fre-quencies. It is important to note that the spread is alsodependent on the image dimensions that dictate the fullextent of the DCT frequency range. TCNQ and FeSeexhibit denser low-frequency clusters and a scattering ofhigh-magnitude mini-clusters– features not present in thediffuse coefficient spread for C . The likely source ofthis spread is the multitude of randomized short-rangeorientations of individual C atoms (inset of Fig. 8).We postulate that the diffuse spread leads to complexfrequency-domain interactions with sampling patternsand the thresholding function, conditions that make itmore difficult to tune the algorithm’s parameters.In our studies, SSIM proved to be a faithful reconstruc- FIG. 8. SSIM evaluated for reconstructions of TCNQ (a,c)and C (b,d) across varying levels of σ (the width of thevariance in the DCT weight model) and threshold ratio (therelative number of of nonzero coefficients used by the opti-mization algorithm). The top and bottom rows respectivelycorrespond to reconstructions performed using Lissajous androtated line sampling patterns. All reconstructions performedwith sampling density ρ = 0 . FIG. 9. Noise perturbation intensity ( δ ) vs. sampling den-sity ( ρ ) MSE phase diagrams for reconstructions of TCNQ(a,d), C (b,e) and FeSe (c,f) for rotated line (top row) andLissajous (bottom row) sampling patterns with defect phasediagrams in the insets. All phase diagrams have been normal-ized to their respective maximum MSE. tion quality metric in terms of capturing the influenceof unwanted artifacts. Reconstructions were also eval-uated for MSE, another commonly used quality metric.MSE lacks SSIM’s useful universal scale, making cross-comparison of images and phase diagrams more difficult.Furthermore, MSE is not adept at capturing structuralartifacts [32], and this flaw is displayed in phase diagramscreated using the metric (Fig. 9). While they mod-erately resemble those for SSIM, these diagrams fail toproperly differentiate between good reconstructions andthose marred by artifacts. As a particularly harsh exam-ple, the poorly reconstructed image of TCNQ at noiseand sampling density equal to 0.1 gives a poor SSIM; the MSE, however, is given a median value. A reverseeffect occurs for FeSe at these parameters, but visual in-spection of the reconstruction yields long-scale structurelargely intact, seemingly further confirming SSIM’s util-ity. However, small-scale structure, i.e. the lattice anddefects, are perturbed, and MSE may be better for cap-turing such anomalies. VI. CONCLUSIONS