Robust Algorithm for Large-Scale Gaussian Patterns Localization

Efficient accurate Gaussian localization is an important topic in many applications, e.g. localization based super-resolution microscopy and image scanning microscopy, which requires large-scale Gaussian patterns localization for accurate super-resolution image reconstruction. Existing Gaussian localization methods usually require high signal-to-noise image and the existing standard fitting algorithm usually requires manually inputting a good initial value for all parameters, which could be not convenient to use and difficult to guarantee high robustness for large-scale Gaussian localizations with a computer. It would be even more challenge to detect all the Gaussian patterns with high-dynamic-range of amplitudes, as well as to estimate a good initial value for all parameters for efficient Gaussian fitting and guarantee high robustness of the localization algorithm for low signal-to-noise ratio image data with strong background. In this paper, we propose an efficient Gaussian patterns detection technique and a robust Gaussian fitting method for accurate Gaussian fitting without initial estimation. In our technique, a fast Pearson correlation algorithm is proposed to improve the efficiency of the calculation of normalized cross correlation for large scale object detection with template matching. By introducing blind background estimation, a modified iterative least-squares Gaussian fitting algorithm without initials estimation is proposed for robust Gaussian fitting with noisy data with strong background. The simulation shows that the performance of the proposed detection technique is high for low SNR image and an efficiency improvement of 27% can be achieved; the proposed Gaussian fitting algorithm is capable of calculating all parameters without initial estimation, and the resulting fitting accuracy is very close to exiting standard methods, which indicates that image signal-to-noise ratio higher than 10dB is required to obtain subpixel accuracy.


I. INTRODUCTION
Gaussian fitting problem is an important topic for many applications, such as localization problem in Image Scanning Microscopy [1], [2], correlation curve and surface fitting in Fluorescence Correlation Spectroscopy [3] and the analysis of atomic resolution images [4], etc. By accurate Gaussian patterns localization and fitting, many information about the studied object can be obtained, for example, the amplitude and sigma indicate the energy and the width of a pulse, and the coordinate indicates the exact position of a molecule.
So far, the neural network and deep learning based localization technique have been developed for localization The associate editor coordinating the review of this manuscript and approving it for publication was Shiqi Wang. microscopy [5], which has been demonstrated powerful. The advantage of deep-learning method is that no prior information is needed, which makes it a straightforward technique without any parameter tuning. However, large scale data set is required for training the network to guarantee good performance. Besides, accuracy fitting is hard to achieve, as an exact solution that is in continuous space is required for a fitting problem. In addition, acquiring high-quality data set could also be a challenge, as the ground-truths for experimental data set usually is not available in practice. These make it not easy to use in practice. In contrast, computer vision (CV) based detection method [6], or signal processing based detection techniques [7] are much easier to implement, with capacity of achieving high efficiency and accuracy, which makes it very suitable to large-scale Gaussian patterns localization applications. We here focus on high-performance signal processing methods, i.e. correlation based detection and fitting based localization method, which has been demonstrated reliable and accurate, and therefore widely used in field of image processing.
Basically, Gaussian patterns localization contains two steps: Gaussian patterns detection and Gaussian fitting. In the first step, the Gaussian patterns detection is implemented by performing correlation between the measured images with an appropriate Gaussian template, and then candidates are determined by peaks-finding from the calculated correlation map. In practice, the template matching is implemented by the normalized cross correlation (NCC) [7]. However, the pixelby-pixel template matching could be very time-consuming. Some acceleration versions of NCC has been proposed, such as sum-table scheme acceleration technique [8], and the hardware acceleration scheme based on field programmable gate array (FPGA) [9].
Based on the detecting result, Gaussian fitting is applied to obtain all parameters of the Gaussian model. So far, there has been several basic fitting methods for Gaussian fitting problem. The existing Gaussian fitting methods mainly are polynomialized Gaussian fitting technique and accurate Gaussian model with non-linear optimization based fitting technique. The polynomialized fitting technique is based on the fact that the Gaussian function contains a quadratic function in the power term of the model function, which therefore can be converted into to a polynomial function by applying the logarithm operator to both side of the model function. By rearranging, the resulting polynomial function is converted to a linearized system. In this way, the coefficients of the polynomial can be estimated by solving the linear function, then all parameters of the Gaussian model function is restored from these estimated coefficients. The basic Gaussian fitting technique of this type is Caruanas' algorithm [10], which is very computationally efficient, however its accuracy can be affected by noise seriously. A modified version of Caruanas' algorithm has been proposed [11] for improving the antinoise performance of fitting algorithm and achieving better fitting accuracy. This algorithm introduced Taylor expanding to analyze the polynomialized function, which indicates that noise could be magnified by the logarithm operation and that weighting the residual error to reduce the impact of noise is necessary; an iteration form is also proposed to reduce the impact of noise for data with long range of Gaussian attenuation, where noise-effect dominates. This algorithm has been demonstrated much more efficient and robust to noise. However background term is ignored in its Gaussian fitting model, i.e. the background of the measured data in application is not taken into consideration at all. The performance of the fitting algorithm could be affected seriously if the data for fitting contains large noise with strong background. Therefore, baseline correction is required to improve the fitting precision [12], otherwise, its accuracy could be affected by noise and the background of measured data. The background in experiment, however, always exists in reality and it is not easy to precisely estimate due to measurement environment by algorithm. For this reason, this type of methods may not suitable to accurate fitting application, such as surface fitting and localization. In practice, the non-linear optimization based fitting method seems to be more often to use, which establishes an accurate fitting model and solve the fitting problem with a non-linear optimization algorithm, such as Least-Square (LSQ) optimization method [13] and Maximum Likelihood (ML) estimator [14]. There has been many existing nonlinear optimization solvers are available for solving 2D Gaussian fitting problem, take Matlab optimization toolbox for example. In practice, most of existing solvers require appropriate initial value for all parameter for generating good fitting result. Even though the non-linear optimization based fitting method is powerful, it might be not convenient to use in practice. In contrast, a solver without preestimating the initial value for all parameters would be much more promising for applications, especially for developing a robust and easy-to-implement algorithm for large scale robust and accurate Gaussian fitting.
In this article, we propose a robust and easy-to-use Gaussian patterns localization technique. Our contribution includes introducing an efficient Gaussian patterns detection technique based on NCC with FFT acceleration and a robust Gaussian fitting algorithm without parameters initial-value estimation. The rest parts of this paper are organized as follow: In part 2, we show the proposed Gaussian patterns detection technique and demonstrate its high suitability for fast large scale Gaussian patterns localization problem. In part 3, we analyze the LSQ based Gaussian fitting method and show a modified Gaussian fitting algorithm without initial estimation. In part 4, we verify the proposed Gaussian localization algorithm and study their performance by simulation. Finally, we summarize our work in part 5.

II. HIGH PERFORMANCE GAUSSIAN PATTERNS DETECTION
In application, robust object detection is an important problem, especially for object detection with noisy data. To address this problem, we here employ correlation-based object detection method and show an efficient detection algorithm for large-scale Gaussian patterns detection.
Since the model is known, the optimal method to detect an object would be by correlation. In practice, Pearson correlation is widely used technique for calculating correlation between two signals [15]. Mathematically, the cross correlation between vector α and β can be represented by where α, β is the average of all elements of α and β respectively. The advantage of the Pearson correlation is that it does not dependent on background as the baseline of the data sets has been removed as formula (1) shows, which makes it very suitable for analyzing data with strong background.

VOLUME 9, 2021
For image correlation calculation, the above correlation is performed to a template T , and a sub-image of an image I centered at position (m, n), and it has to goes through all pixels to calculate a complete correlation map. Apparently, the drawback of this approach is that the computing is timeconsuming for large scale image processing. To avoid this problem, we here introduce a fast algorithm to implement Pearson correlation between two images. Actually, the formula (1) can be rewritten by Assume that the size of the sub-image is M W × N W , which is identical to that of the template and W is the window function defining the support of a sub-image, i.e.
which are corresponding to the sum, average and the L2-norm of template T respectively; let which are corresponding to the sum, sum-of-squares and the average of the sub-image, respectively; let τ = 1 M W ·N W , the cross correlation between a template and an image then can be represented by As the correlation operation can be calculated by the fast Fourier transformation (FFT) [16], T ⊗ I , S 1 and S 2 can be calculated efficiently. Therefore, the proposed technique has advantage in large scale image processing. The correlation detection technique tries to find a pattern similar to the template in morphology, without being sensitive to the amplitude. This, however, could result in over sensitive to the Gaussian artifacts induced by noise. To avoid this drawback, we introduce an energy-based metric to improve the detection accuracy. In fact, the term S 2 − τ S 2 1 , which calculates the variance of all pixels inside the window W , already shows local-regional energy distribution in the measured image. Theoretically, the peaks of energy distribution take place at the centers of the Gaussian patterns. Since the artifact usually contains much lower energy than a true Gaussian pattern, a threshold can be introduced, i.e. Th( S 2 − τ S 2 1 ), to remove the candidates with low energy. Formula (3) then is modified by where Determine an appropriate threshold value is essential to guarantee the performance of detection. We here use the energy of the template as the threshold, e.g. τ = ν 2 . As the template is normalized, the threshold τ = (A · ν) 2 , where A denotes the amplitude of the Gaussian pattern with the lowest energy that can be detected, can be used instead. The detail of implementing this algorithm is shown in Algorithm 1.
Apparently, a Gaussian pattern probably takes place where the correlation value in the correlation map is high. As the correlation value is between 0 and 1, it is very easy to detect all potential Gaussian patterns by thresholding technique and get their rough positions by finding the peaks in the calculated correlation map by Algorithm 1.

III. LSQ BASED GAUSSIAN FITTING WITHOUT INITIAL ESTIMATION
Original LSQ Gaussian fitting has been demonstrated a much more efficient technique for Gaussian fitting problem. However, the intrinsic drawback is that the background is not included in the model function, which impedes its widespread application in practice and makes it only suitable to data without background. We here show that it is indeed also possible to be applied to data set with strong background, by introducing an outer-loop iteration based on existing highperformance LSQ algorithm and blind background estimation process.
Mathematically, a 2D Gaussian model function can be represented by where A denotes the amplitude, and (a, b) is the coordinate of the Gaussian pattern, σ x and σ y are the sizes of the Gaussian pattern in each direction. Take the logarithm on both side of equation (1), which yield Linearized the equation (2) with respect with the coefficients of Gaussian model function, we can get where y, z), and write a cost function based on least-squares, we have Apparently, there is a drawback in the above cost function, i.e. the logarithm is nonlinear transformation, which enhanced pixels with low intensity; due to the nonlinearity of the logarithm function, the pixels with lower intensity could be affected by the noise seriously. A good way to avoid this problem would be introducing a weight to depressed the noise amplification effect [11], the cost function then is rewritten by (9) namely, Set the gradient of the above cost function over all parameters and a linear system can be obtained by It shows that M actually is a symmetric matrix, therefore, only half of the elements in M are needed to calculate. In practice, the Gaussian model function itself or its power function can be used as the weight function, i.e.
where G k−1 is the predicted object function calculated with the estimated parameters in last iteration and the γ is a parameter controlling the penalty degree, so that the above linear system can be solved iteratively. The detail of this iterative version of LSQ (iLSQ) is shown in Algorithm 3. The advantage of iLSQ Gaussian fitting algorithm is that no initial estimation is needed, which makes implementing Gaussian pattern localization very convenient in practice, while high robustness to noise can be achieved. However, as mentioned above, the drawback is that its accuracy can be affected by a background, which actually is very common in experiment. As background is not taken into consideration in the model function, it cannot be used to data with strong background. To avoid this problem, we here introduce a simple algorithm, which is capable of estimating all parameters based on iLSQ without initial estimation.
The algorithm starts with estimating the background by taking the minimum pixel of the input image. Then it removes the background from the measured image and estimates the parameters except the background by the iLSQ method with the background-removed image. With the calculated parameters, the Gaussian model function without background is estimated. The background then is estimated by taking the average of the difference between the estimated model function and the measured image. This process is repeated until the algorithm converge. By this way, the algorithm does not depend on the initial value of all parameters anymore. This makes it suitable for noisy data with strong background, to which estimating the initial value from the measured data for each parameter could be challenge for general non-linear optimization based Gaussian fitting methods. In comparison with general fitting techniques, the proposed algorithm is much easier. The detail of the implementation of the proposed algorithm is shown in Algorithm 2.

IV. RESULTS
We verify the proposed techniques by simulation. First, we simulate large-scale focal light-spots localization, then test the performance of the proposed Gaussian fitting algorithm and make a comparison with existing techniques. The configuration of the computer used is: Quad-core CPU (i7-6700HQ, 2.6GHz) with 16GB RAM. The algorithm coding and testing environment is based on Matlab 2018b platform.
In the simulation on large-scale Gaussian patterns detection, the multi-focal measurement [17] for super-resolution confocal microscopy are simulated by generating the homogenous Gaussian patterns. In order to simulate a data set close to the real case in experiment, we choose some typical values of the parameters of the Gaussian pattern and VOLUME 9, 2021 noise level. The background is set to 200, the sigma value for each Gaussian pattern is set to 2, and the amplitudes are generated randomly, which is set within the range between [100,200], which is similar to a case with moderate signalto-noise ratio (SNR) in most of general experiments. Poisson noise is applied to the simulated data to simulate camera shot noise. The noisy image is just as Figure 1 (a).
In order to obtain high-quality correlation map, the noisy image are de-noised by the median filter. The correlation map then is calculated by (4) with the de-noised image. To detect all Gaussian patterns with high accuracy, a threshold value, e.g. 0.6, is set to guarantee only the strong correlation peaks are selected for detection.
The correlation map is calculated by (4) with the normalized Gaussian template, which however requires an appropriate sigma value for high-quality correlation map. In our algorithm, the sigma value of the Gaussian template is estimated by scanning approach starting from σ = 1, until the peak-value of the average of all pixels higher than 0.5 in the correlation map is reached. With the estimated sigma value, a correlation map is calculated for fine estimation.
Then 100 candidates with highest correlation value are used to estimate all parameters by the proposed Gaussian fitting algorithm, by which the background and the sigma value of the candidates are obtained for accurate detection. The value of A used to calculate the energy-based metric threshold τ is set to 20. To make a comparison with existing computer vision (CV) detection method, the same data set is also processed by general OpenCV detection method, which binarizes the image by hard thresholding and then calculates the center of each Gaussian pattern by connectedCompo-nentsWithStats() function. This technique is widely used to small objects detection in application. As the OpenCV technique actually is an amplitude-based detection method, which means that a proper threshold value estimation is always required. In our simulation, the threshold value for OpenCV detection function is simply set to a value that shows the best result.
The simulation results in Figure 1 show that both techniques are able to detect all Gaussian patterns successfully with high SNR image. To demonstrate the performance of the proposed technique, the calculated result of the correlation Algorithm 2 Gaussian Fitting Algorithm input: image data I for fitting, threshold value , maximum iterations k max ; begin k = 1; Bg = min (I ) ; while (true) The 2D iLSQ Algorithm input: image data I for fitting, threshold value , maximum iterations k max ; begin k = 1; w = I ; while (true) get M by (17); get K by solving the linear equation in (16); break; end k = k + 1; end end Output: θ k map and the detection result are shown. Figure 1 (c) shows the correlation map filtered by the proposed energy-based thresholding and Figure 1 (d) shows the filtered correlation map with correlation-value thresholding. The peaks in correlation map, which correspond to the centers of the Gaussian patterns, are high and their amplitudes are quite homogenous, even though the large dynamic of amplitudes of the Gaussian patterns and the strong background. It shows that the filtered correlation map is very clean even though the data set is noisy and the correlation peaks clearly shows the rough positions of all Gaussian patterns. As these surfaces of correlation map are convex and so homophonous, all Gaussian patterns can be detected by peak finding method easily. In comparison with amplitude-based technique, the proposed correlation based detection is not sensitive to noise so much.
We then further demonstrate the performance of the proposed detection technique with low SNR image data. In simulation, we generate homogenous Gaussian patterns with strong background with lower amplitudes to simulate the low SNR focal light-spot images. The amplitudes of the Gaussian patterns are set within range [30 100], and the background is set to 200. Poisson noise is applied to the image data as well. The simulated image is just as Figure 2 (a) shows. It shows that the background is so strong that many Gaussian patterns are hard to be recognized, even from the de-noised image. The simulation results are just as Figure 2 (b)-(d) shows. Figure 2 (b) shows that the OpenCV detection method cannot guarantee all Gaussian patterns are detected for low SNR data set, even with optimal threshold-value. In contrast, Figure 2 (c) shows that the correlation peaks are still quite homogenous, and all Gaussian patterns can be detected successfully. This indicates that correlation based detection technique with thresholding is much easier and more robust than the amplitude-based detection method, which usually involves complex threshold-value estimation.
Simulations show that the noise causes a lot of artifacts in calculated correlation map, which no doubt could affect the detecting accuracy seriously. By introducing energy-based metric, the artifacts is removed almost completely. Just as an appropriate sigma value is important for high-quality correlation map, the energy-based metric required an appropriate threshold value is essential for accurate detection to noisy data. Simulations show that taking a value of 1.5 folds of the standard variance of noise for the value of the thresholding amplitude A is a good choice. Simulation also shows that the proposed algorithm can be more efficient than existing standard method. It takes 33ms to calculate a correlation map, while the Matlab normxcorr2 function, i.e. the abovementioned fast NCC method, takes 42ms. As the proposed correlation technique is totally based on FFT, it can be more efficient than the standard method.
To show the performance of the proposed algorithm to noise, we test the fitting accuracy to different SNR level of simulated Gaussian patterns. The simulation result are just as Figure 4 shows, in which some results by existing fitting techniques are also shown for comparison. It shows that the root-mean-square error (RMSE) of all Gaussian fitting techniques decreased with increasing the SNR, but the resulting accuracies are different. The RMSE by iLSQ algorithm without background removal is much higher than that with background removal, which indicates that background is an important factor that could affect the accuracy of iLSQ algorithm seriously. The performance of the proposed algorithm is quite close to the iLSQ with background removal in accuracy. For low SNR, the accuracy of iLSQ with background removal and the proposed algorithm is lower than Matlab fit function, but it is not so much and still acceptable; for high SNR, all  algorithms work quite well and achieve similar accuracy. It also shows that an image SNR higher than 10dB is required for sub-pixel accuracy.   It shows that iLSQ is the most efficient, while the proposed algorithm takes more CPU time to calculate all parameters, but still much more efficient than the Matlab fitter. Simulation shows that the proposed algorithm takes about 10 iterations to converge, for which it takes about 10-fold CPU time of general iLSQ algorithm; this actually is the cost for optimal solution searching by the algorithm without prior information regarding the parameters. However, it makes the algorithm quite convenient to use in practice, since the initial-value estimation for all parameters is avoided.

V. CONCLUSION
In this paper, we show a simple algorithm for robust largescale Gaussian localization. In the proposed technique, a high-performance efficient detection technique based on NCC with FFT acceleration is proposed to localize all Gaussian patterns roughly and a modified iLSQ algorithm without parameters initial-value estimation is proposed to calculate their accurate positions. Simulations show that largescale Gaussian patterns can be detected efficiently by the proposed correlation based recognition algorithm with high robustness to noise. The proposed Gaussian fitting technique is capable of estimating all parameters of Gaussian model including background without dedicated initial estimation, which makes it much more convenient to use in application than existing non-linear optimization based fitting methods. These results show high performance and potentiality for the application of large-scale robust Gaussian patterns localization in the proposed technique. In future work, developing parallel computing package for the proposed algorithms for large-scale Gaussian patterns localization is necessary.
SHUN QIN received the M.Sc. degree in electronic science and technology from the Beijing Institute of Technology, Beijing, China, in 2014, and the Ph.D. degree in physics from the University of Göttingen, Göttingen, Germany, in 2019. His research interests include signal/image processing, advanced computing, optical imaging, inverse problem, super-resolution microscopy, simulation, and digital system design.