Microstructure Images Restoration of Metallic Materials Based upon KSVD and Smoothing Penalty Sparse Representation Approach

Microstructure images of metallic materials play a significant role in industrial applications. To address image degradation problem of metallic materials, a novel image restoration technique based on K-means singular value decomposition (KSVD) and smoothing penalty sparse representation (SPSR) algorithm is proposed in this work, the microstructure images of aluminum alloy 7075 (AA7075) material are used as examples. To begin with, to reflect the detail structure characteristics of the damaged image, the KSVD dictionary is introduced to substitute the traditional sparse transform basis (TSTB) for sparse representation. Then, due to the image restoration, modeling belongs to a highly underdetermined equation, and traditional sparse reconstruction methods may cause instability and obvious artifacts in the reconstructed images, especially reconstructed image with many smooth regions and the noise level is strong, thus the SPSR (here, q = 0.5) algorithm is designed to reconstruct the damaged image. The results of simulation and two practical cases demonstrate that the proposed method has superior performance compared with some state-of-the-art methods in terms of restoration performance factors and visual quality. Meanwhile, the grain size parameters and grain boundaries of microstructure image are discussed before and after they are restored by proposed method.


Introduction
Microstructure images of metallic materials play a significant role in industrial applications, such as metallography detection, surface topography measurement (STM), and micro-electro mechanical manufacturing (MEMM) [1][2][3], etc. However, the raw microstructure images are easily contaminated on their acquisition, storage, and transmission, which degrade the fidelity of the microstructure image. Further, in most cases, images re-cquisition is not possible in practice, as such damaged images are not suitable for subsequent processing [4,5].
Image restoration is an effective technique for recovering incomplete or damaged images approximate to the ideal images. Up to now, various image restoration strategies have been proposed, which can be classified into several categories: (i) Time domain analysis methods, such as adaptive filter denoising (AFD) [6,7]; (ii) Frequency domain analysis methods, such as wavelet-wavelet packet denoising (W-WPD) technique [8,9]; (iii) Data-driven approaches, such as partial differential equation (PDE) [10][11][12][13], and wavelet hidden Markov random field (WHMRF) [14]; (iv) Sparse representation (SR) techniques, such as redundant dictionary and non-convex penalty regularization [15][16][17][18], etc. Although the damaged images can be restored more precisely by the above methods, the drawbacks are also obvious. For example, the adaptive filter denoising method can achieve a good denoising effect under low noise level, but the denoising effect will be greatly reduced with the noise increases. The W-WPD can improve the quality of image, but some artifacts (e.g., ambiguity points) might be generated along with the reconstruction process. The main shortcomings, including optimal parameter selection and high computation cost, still remain unsolved in data-driven approaches. For the WHMRF method, it is difficult to flexibly integrate the spatial quantitative relation (between the damaged pixel point and its neighbor points) into the restoration model [19].
Today, SR techniques are based on the principle that an image can be sparsely represented in a time domain or frequency transform domain, where each image block could be reconstructed by sparse reconstructing algorithms. The core idea behind SR is how to represent the images more sparsely in the time/frequency transform domains, usually, the most common methods focus on the establishment of redundant multiscale transform and redundant dictionary. For example, in ref. [20], the curvelet transform is proposed to denoise the white noise, which exhibits higher perceptual quality than wavelet-based reconstructions. In ref. [21], the edge detection and fuzzy clustering algorithm are combined to preserve the edges of synthetic aperture radar (SAR) despeckling in translation-invariant second-generation band-wavelet transform (TIBT) domain. In ref. [22], to address the image denoising problem, the over-complete discrete cosine transform (DCT) dictionary, global trained dictionary, and K-means singular value decomposition (KSVD) dictionary, were proposed.
Another core point of the SR is reconstruction algorithm. Generally, image restoration is an inverse problem, in which the ideal images could be approximatively estimated from the noisy images. The difficulty knot of image restoration is that the inverse problems are often ill-posed (or non-deterministic polynomial-time hard, NP-hard). In traditional sparse representations, most of the methods are applied based on regularization-based technique, for example, matching pursuit (MP) [23], orthogonal matching pursuit (OMP) [24], and compressive sampling matching pursuit (CoSaMP) [25], some regularization models, such as total variation regularization function (TVRF) [26], sparse non-local regularization (SNLR) [27,28], etc. However, the traditional sparse representation methods may cause instability and obvious artifacts in the reconstructed images, especially for the restoration of microstructure image that include some smooth regions and multi-boundary and fine-textures, or when the noise level is strong.
To address the above issues in SR and exploit the spatial information of microstructure images, in this paper, a novel image restoration method based upon smoothing penalty sparse representation (SPSR) and adaptive over-complete KSVD dictionary is proposed for microstructure image of metallic materials, taking aluminum alloy 7075 (AA7075) material as an example. The nonconvex penalty regularization is introduced to address the image inverse problem mentioned above, and the success rate will be improved greatly in image reconstruction. Moreover, unlike the common procedure used in refs. [29,30], in which the sparse transform basis was used for measuring the image sparsity, this paper utilizes over-complete dictionary (e.g., trained KSVD dictionary) to promote the image sparsity under the given redundant system. The simulation and experimental results show the superiority of the proposed method compared with some state-of-the-art methods, such as wavelet packet method, the discrete cosine transformation (DCT) combined with OMP method, and the KSVD dictionary combined with OMP method, respectively. Meanwhile, the comparison results of grain parameters, such as grain diameter (mean), grain area (polygon), grain perimeter (ratio), etc., are also discussed in detail before and after processing based on the proposed method.
For the applicability of the proposed method, generally, a microstructure image is used that exhibits a lamellar (plate-like) structure or a structure that exhibits twinning, which can be divided into two kinds of regions: one is the image blocks containing lamellar (plate-like) structure and boundaries in structure regions, etc., and the other is the image blocks that distributed in smooth regions. Correspondingly, the optimization expression of the sparse representation method usually includes a likelihood part and prior knowledge part, furthermore, the external noise will play a critical role in smooth regions during image restoration, and its effect is relatively weak in structured regions, where there is a strong similarity between the two blocks, due to the pixel values being similar in smooth regions, and if the noise level is strong, the information characteristics of external noise in smooth regions might be regarded as structural information in sparse coefficients, wherein the classical optimization approach might remove the inveracious structure information that failed, leading to instability and obvious artifacts on the reconstructed image. For the proposed method, on the one hand, before the reconstruction algorithm is implemented, the information of lamellar (plate-like) structure and boundaries could be represented in over-complete KSVD trained dictionary, on the other hand, the prior knowledge part is improved by introducing a smoothing parameter, and the likelihood part is also improved via a regularization weight, therefore, the noise distributed in structure regions and smooth regions could be denoised adaptively, especially for the noise points hidden in the smooth regions.
The main contributions of this paper are summarized as follows: (1) The dictionary training algorithm, namely KSVD, is introduced, while the detailed structure information, such as lamellar (plate-like) structure and boundaries, could be represented accurately. (2) The smoothing parameter and regularization weight are designed based upon the traditional sparse representation method, and the noise distributed in structure regions and smooth regions could be denoised adaptively. (3) The grain parameters, such as grain diameter (mean), grain area (polygon), grain perimeter (ratio), etc., are discussed before and after they are processed by the proposed method, and the structural information that used for industrial applications, such as metallography detection, micro-electro mechanical manufacturing (MEMM) could be clearly isolated, which may open up a new field of application of material microstructure to industry.
The layout of this paper is organized as follows. Section 2 describes the algorithms of sparse representation and KSVD. Section 3 introduces smoothing penalty sparse representation (SPSR) algorithm, and its derivation and parameter selection, etc. In Section 4, the simulation and experimental results of the proposed method are presented with other approaches. Conclusions are shown in Section 5.

The Review of Sparse Representation (SR)
Generally, the image degradation process can be described by a degradation matrix R as follows: where y is the damaged image, x is the latent ideal image and n is external noise. Usually, image x is assumed to be K-sparse, i.e., which having K non-zero components (the pixel value of the measurement matrix in the damaged area is 0, in contrast, undamaged area is 1). According to the framework of compressed sensing (CS) [31,32], a discrete 1-D signal with limited length x 0 that can be represented in the domain ϕ = [ϕ 1 , ϕ 2 , . . . , ϕ N ], where α is the transformation coefficients of signal x in domain ϕ. Usually, the commonly used sparse transformation base ϕ includes discrete wavelet transform (DWT) basis, discrete cosine transform (DCT) basis, Gabor basis and curvelet basis, etc. [33]. Further, if the 1-D discrete signal is expanded into 2-D image, i.e., x 0 → x , Equation (1) can now be rewritten as follows: where A = Rϕ denotes sensing matrix. Here, Candes et al. [32,34,35] showed that the Equation (3) can only be solved when the sensing matrix A satisfies the conditions of incoherence or restricted isometry property (RIP), where the RIP is given by where parameter δ controls the level of data discrepancy based on an estimation of the noise variance. Usually, the common measurement matrix includes Gaussian random matrix (GRM), local Fourier matrix (LFM) and local Hadamard matrix (LHM), etc.

KSVD Algorithm
To enhance the sparsity of the damaged image, in this section, the adaptive over-complete dictionary is employed to substitute the traditional sparse transform basis, that is, the image blocks will be represented by trained KSVD method [36], i.e., A = RD. It should be noted that the matrix A may not satisfy the RIP criterion, due to the over-complete dictionary D; to address this issue, the mutual coherence technique [37] is applied to substitute the matrix A, i.e., where α i and α j are the i-th and j-th column in matrix A, respectively. The adaptive over-complete KSVD dictionary method is designed as follows: Step (1). Dictionary initialization. Suppose that D ∈ R n×k (k << n) is the over-complete dictionary, for example, over-complete discrete cosine transform (DCT) dictionary.
Step (2). Dictionary presentation. Suppose that D ∈ R n×k (k << n) is the over-complete DCT dictionary, Y = {y i } N i = 1 is a N-training samples set and X = {x i } N i = 1 is a solution vectors set of Y, hence, min where T 0 is the maximum value of nonzero vector in sparse coefficients and · F is Frobenius norm.
Step (3). Dictionary updated. Assume that d k refer to the k-th column of the pending update dictionary, then Equation (6) equivalent to where E k represents pre-computed error matrix, x k T represents the k-th row of X which actually gives the sparse coefficients or weights of d k . For the purpose of singular value decomposition (SVD), four parameters are defined as follows: where Ω k is a matrix with size N × |ω k |. The parameter Y R k includes a subset of the samples that are currently using d k element, E R k represents error columns that correspond to samples that using d k element, respectively. Returning to the objective Equation (7), it can be rewritten as Step (4). Singular value decomposition. The E k is decomposed by SVD method, i.e., E k = U∆V T , ∧ d k is the first column updated by d k in U, then the second column will be updated successively until the last column, otherwise, go back to step (3), and the new trained dictionary ∧ D will be finally obtained.

Smoothing Penalty Sparse Representation (SPSR)
It should be noted that y = Aα + n belongs to a highly underdetermined equation (HUE) [38,39], there are infinite set of solutions. The problem of image reconstruction by sparse representation under residual error constraint can be calculated by where c is a threshold of the residual error. Moreover, the prior knowledge of the original image is usually utilized to regularize the solution under residual error constraint is expressed as where λ is regularization weight and ζ(x) regularization term. From the perspective of Bayesian estimation, the Aα − y 2 2 and ζ(x) can be viewed as the likelihood part and prior knowledge part, respectively. Therefore, the ζ(x) prior knowledge plays a significant role in image restoration based on sparse representation. For a microstructure image, it can also be divided into two types: the first is the image blocks containing details, boundaries, and singular points, and the second is the image blocks located in smooth regions. For the former, the details, boundaries, and singular points that determine the similarity between two blocks, and the influence of external noise, is relatively weak. However, in smooth regions, there is a strong similarity between the two blocks due to the pixel values being similar; here the influence of external noise will play a critical role in image restoration. If the noise level is strong, the information of noise in smooth regions is regarded as structural information in sparse coefficients. Meanwhile, the classical optimization approaches and regularization approaches cannot remove the false structural information contained in the sparse coefficients, and the traditional methods may cause instability and obvious artifacts in the reconstructed images.
represents the rearrangement of absolute values of r(α (k+1) ) in the decreasing order, and r(α) s+1 is the (s + 1)th component value of r(α). Note that, if ε k+1 = 0, choose α (k+1) to be an approximation of sparse solution and stop this iteration. (7) Let k = k + 1, and return to Step (4) to continue. End For the SPSR method, the following theorem summarizes the results for 0 < q ≤ 1, thus, we have the following theorem which can prove the above proposed algorithm. Theorem 1. Error estimation theorem [40]. Suppose that x 0 is s-sparse signal which satisfies Ax 0 = b. The smooth parameter ε k → ε * with k → ∞ . Matrix A satisfies the RIP of order 2s with δ 2s < 1, when ε * > 0, the sequence {x (k) } has at least one convergent subsequence. Suppose that the limit ε k = ε * is a local optimal solution for Equation (12), we have, where δ s (x ε * ) 2 is the approximate error of x ε * , which satisfies δ s (x ε * ) 2 = inf y 2,0 ≤s x ε * − y 2 . For the special case, when ε * = 0, there must exist a convergent subsequence converging to point x 0 , it satisfies, where C 1 , C 2 are C 3 are independent positive constants. To prove the Theorem 1, the following two lemmas are required.
Proof. According to arithmetic-geometric mean inequality [41], i.e., This completes the proof of Lemma 1. Furthermore, where C 4 is an independent positive constant.
Proof. We first compute the following formula, According to [ Using Lemma 1 and substituting Equation (20) to Equation (19), we have From the result of Equation (21), the Equation (17) can be calculated. It should be noted from Equation (17) that the L q (x (k) , ε k , λ) is monotonically decreasing sequence, hence, for all k ≥ 1 and 1 ≤ i ≤ n, there exists a positive number β which satisfies , and the Lemma 2 is proved conclusively.
Herein, combining the above inequalities in Lemma 1 and Lemma 2, the Theorem 1 can be proved ultimately, for simplicity, the detailed proof process was omitted. In the next section, the choice of regular operator q will be discussed in detail via a simulation experiment.
For the choice of regular operator q (0 < q ≤ 1), let q varying among {0.1, 0.5, 0.7, 1}. Firstly, the matrix A is designed by rand-function rand(64, 256) in MATLAB, and the initialization signal x 0 has t non-zero narrow impulses that subject to the standard Gaussian distribution (SGD), the locations of non-zeros are uniformly and randomly generated, and t, varying among {8, 10, 12, . . . , 32}. The penalty parameter λ = 10 −6 is small enough to approximately enforce Ax = Ax 0 and δ = 0.09, which is measured over 100 times in terms the perfect reconstruction. Taking the SPSR algorithm iterative 1000 times, if the recovery error satisfies x r − x 0 2 / x 0 2 ≤ 10 −3 , the iteration is stopped, where x r stands for a recovered vector. The recovery algorithm is SPSR method (q ∈ {0.1, 0.5, 0.7, 1}), the recovery success rate curve of different. q with sparsity is shown in Figure 1. From Figure 1, it can be seen that q = 0.1, q = 0.5 performed better than q = 0.7 and much better than q = 1. Moreover, q = 0.5 gives a higher success than q = 0.1 slightly. We emphasize that our results do not counter the intuition that a smaller q should recover more sparse vectors. Generally, a smaller q value makes the minimizing functional more non-convex, but more difficult to solve. In addition, in this algorithm, it has been discovered that if smoothing parameter ε decreased slowly, the performance of q = 0.1 reach further improved. However, the running time with q = 0.1 also became much longer. For example, in this simulation, the execution time of parameter q = 0.5 is 4.4550 s, and execution time of q = 0.7 and q = 1 are 5.7359 s and 8.6348 s, respectively, but the execution time of q = 0.1 is 9.4643 s, it should be noted that the execution time becoming longer when parameter q is selected smaller.

Case Verifications and Disscussion
In order to quantitatively calculate restoration effect, three restoration performance standards (RPSs) are employed for comparison, i.e., the normalized mean-square error (NMSE), peak signal-tonoise ratio (PSNR), and structural similarity (SSIM), respectively. It should be noted that the PSNR reflects the approximation degree between the recovered image and original image, and the SSIM reflects the similarity degree between the recovered image and original image. Three RPSs are

Case Verifications and Disscussion
In order to quantitatively calculate restoration effect, three restoration performance standards (RPSs) are employed for comparison, i.e., the normalized mean-square error (NMSE), peak signal-to-noise ratio (PSNR), and structural similarity (SSIM), respectively. It should be noted that the PSNR reflects the approximation degree between the recovered image and original image, and the SSIM reflects the similarity degree between the recovered image and original image. Three RPSs are defined as where x is an original image and ∧ x is a recovered image.
where µ is the mean of image, σ is variance or covariance, c 1 and c 2 is the small constant which can enhance the stability of calculation results, respectively. In the following experiments, the microstructure image is divided into 8 × 8 blocks. The DCT dictionary is used as initial dictionary, and KSVD is used as dictionary training algorithm, then the damaged microstructure images are applied as training samples. The KSVD dictionary training parameters setting is shown in Table 1.

Simulation Case of Aluminum Alloy 7075 Material
To verify the superiority of proposed SPSR (q = 0.5) approach, an ideal microstructure orientation image of aluminum alloy 7075 material is generated by cellular automaton simulation (CAS) method [42,43], see Figure 2, Figure 2a is the ideal color image, and Figure 2b is the grayscale image, respectively. The artificial Gaussian noises are added in the ideal microstructure image, and Figure 3 presents the visual quality of the reconstructed images based on wavelet packet and SPSR (here, q = 0.5) algorithms when the noise level is fixed on 0.1, 0.5, 0.8, and 1.0, respectively. It can be seen from Figure 3 that the proposed SPSR (q = 0.5) method clearly provides a significant improvement to the fidelity of the resultant image (see the 4th column). The wavelet packet method is also applied to process the noisy orientation image (see the 3rd column), and provides some improvement, but not nearly as dramatic as that achieved with the SPSR method. More importantly, as we can see, the proposed SPSR method achieves a significantly better visual quality than wavelet packet method in smooth regions.
from Figure 3 that the proposed SPSR (q = 0.5) method clearly provides a significant improvement to the fidelity of the resultant image (see the 4th column). The wavelet packet method is also applied to process the noisy orientation image (see the 3rd column), and provides some improvement, but not nearly as dramatic as that achieved with the SPSR method. More importantly, as we can see, the proposed SPSR method achieves a significantly better visual quality than wavelet packet method in smooth regions.  from Figure 3 that the proposed SPSR (q = 0.5) method clearly provides a significant improvement to the fidelity of the resultant image (see the 4th column). The wavelet packet method is also applied to process the noisy orientation image (see the 3rd column), and provides some improvement, but not nearly as dramatic as that achieved with the SPSR method. More importantly, as we can see, the proposed SPSR method achieves a significantly better visual quality than wavelet packet method in smooth regions.    Table 2 shows quantitative results of the restoration performance standards (RPSs) obtained for the images shown in Figure 3. It can be see that the wavelet-packet method does provide a moderate improvement particularly at the higher noise levels as evidenced by Figure 3. Further, as the images show as well, the SPSR (q = 0.5) method provides much more significant improvement in RPSs.
It should be noted that adding more noise to the orientation image lead to much higher NMSE, lower PSNR and SSIM, and the data did not follow the expected trend, likely due to so many missing points. In the next two sections, two kinds of damaged aluminum alloy 7075 (AA7075) orientation images are investigated. The first damaged image is contaminated, due to the charged particles and dust that exists in the electron back-scattered diffraction (EBSD) system. The second damaged image is a low-pixel (or low-resolution) image that is contaminated due to the thicker contamination membrane residues on the AA7075 sample surface. The orientation images were acquired by the Oxford Instruments AZtecHKL EBSD system. Orientation images were collected and recorded at 114 × 114 pixel. The Hough was run on the patterns after compression to 96 × 96 pixels with a 9 × 9 convolution mask, 1 • theta step size, and searching for a maximum of 10 peaks.

Experimental Case 1
The initial DCT dictionary of the first damaged image and the over-complete dictionary that was trained by KSVD algorithm are shown in Figure 4a,b, respectively. It is obvious that the structure of KSVD dictionary is more abundant than the DCT dictionary, the reason is that the KSVD dictionary fully reflects the detailed structure characteristics of the image, equivalently, the KSVD dictionary can precisely represent the sparse characterization of block sub-image.  Table 2 shows quantitative results of the restoration performance standards (RPSs) obtained for the images shown in Figure 3. It can be see that the wavelet-packet method does provide a moderate improvement particularly at the higher noise levels as evidenced by Figure 3. Further, as the images show as well, the SPSR (q = 0.5) method provides much more significant improvement in RPSs. It should be noted that adding more noise to the orientation image lead to much higher NMSE, lower PSNR and SSIM, and the data did not follow the expected trend, likely due to so many missing points. In the next two sections, two kinds of damaged aluminum alloy 7075 (AA7075) orientation images are investigated. The first damaged image is contaminated, due to the charged particles and dust that exists in the electron back-scattered diffraction (EBSD) system. The second damaged image is a low-pixel (or low-resolution) image that is contaminated due to the thicker contamination membrane residues on the AA7075 sample surface. The orientation images were acquired by the Oxford Instruments AZtecHKL EBSD system. Orientation images were collected and recorded at 114 × 114 pixel. The Hough was run on the patterns after compression to 96 × 96 pixels with a 9 × 9 convolution mask, 1° theta step size, and searching for a maximum of 10 peaks.

Experimental Case 1
The initial DCT dictionary of the first damaged image and the over-complete dictionary that was trained by KSVD algorithm are shown in Figure 4a,b, respectively. It is obvious that the structure of KSVD dictionary is more abundant than the DCT dictionary, the reason is that the KSVD dictionary fully reflects the detailed structure characteristics of the image, equivalently, the KSVD dictionary can precisely represent the sparse characterization of block sub-image.    From Figure 5c-f, it can be observed that all the reconstruction methods are able to reasonably recover the damaged image.
The quantitative RPSs between the actual and recovered image based on the above methods are summarized in Table 3. In this work, due to the ideal image being unknown, the above RPS results are designed based on damaged and recovered image, therefore, it should be noted that the larger the value of NMSE, the smaller the values of PSNR and SSIM, and the more accurate of restoration effect.
Form Table 3, it is clear that the proposed approach achieves better restoration result than other methods. In the smooth regions, through the proposed approach, one can find that the clarity and fidelity of image grain-boundary are enhanced significantly, and the influence of the speckle noise is reduced greatly.   The quantitative RPSs between the actual and recovered image based on the above methods are summarized in Table 3. In this work, due to the ideal image being unknown, the above RPS results are designed based on damaged and recovered image, therefore, it should be noted that the larger the value of NMSE, the smaller the values of PSNR and SSIM, and the more accurate of restoration effect.
Form Table 3, it is clear that the proposed approach achieves better restoration result than other methods. In the smooth regions, through the proposed approach, one can find that the clarity and fidelity of image grain-boundary are enhanced significantly, and the influence of the speckle noise is reduced greatly.
Furthermore, to explore the grain distribution and its statistical result before and after restoration, the maximum entropy threshold segmentation (METS) algorithm [44] is introduced in Figure 5b,f, respectively. The segmentation effect and grain number tags created by image pro-plus software [45] are shown in Figures 6 and 7, respectively. It is obvious that breakpoints of segmentation curves in Figure 7a are less than Figure 6a, and the number of grains in Figure 7b is more than Figure 6b. Furthermore, to explore the grain distribution and its statistical result before and after restoration, the maximum entropy threshold segmentation (METS) algorithm [44] is introduced in Figure 5b,f, respectively. The segmentation effect and grain number tags created by image pro-plus software [45] are shown in Figures 6 and 7, respectively. It is obvious that breakpoints of segmentation curves in Figure 7a are less than Figure 6a, and the number of grains in Figure 7b is more than Figure 6b.   Table 4 summarizes the statistical results of grain size, which consists of grain diameter (mean), grain area (polygon) and grain perimeter (ratio). It can be seen that the sample number is 64 before being restored, and it has increased rapidly to 232 after being restored by the proposed method; this is because more grains are calculated, due to the sealing character of the segmentation curves in Figure 7b.  Furthermore, to explore the grain distribution and its statistical result before and after restoration, the maximum entropy threshold segmentation (METS) algorithm [44] is introduced in Figure 5b,f, respectively. The segmentation effect and grain number tags created by image pro-plus software [45] are shown in Figures 6 and 7, respectively. It is obvious that breakpoints of segmentation curves in Figure 7a are less than Figure 6a, and the number of grains in Figure 7b is more than Figure 6b.   Table 4 summarizes the statistical results of grain size, which consists of grain diameter (mean), grain area (polygon) and grain perimeter (ratio). It can be seen that the sample number is 64 before being restored, and it has increased rapidly to 232 after being restored by the proposed method; this is because more grains are calculated, due to the sealing character of the segmentation curves in Figure 7b.   Table 4 summarizes the statistical results of grain size, which consists of grain diameter (mean), grain area (polygon) and grain perimeter (ratio). It can be seen that the sample number is 64 before being restored, and it has increased rapidly to 232 after being restored by the proposed method; this is because more grains are calculated, due to the sealing character of the segmentation curves in Figure 7b.  Figure 8 illustrates the histogram distribution of grain diameter, and Figure 9 shows the relationship between the grain diameter and grain area before and after being restored, respectively. From Figure 8a, it should be noted that the grain number is not counted when the grain size larger than 10 microns, because most of the grain boundaries are broken in the original damaged image, and only some continuous grains will be counted by the image pro-plus software; on the contrary, the grain can be calculated easily based upon the proposed method. From the results shown in Figure 9, the correlation value is 0.5537 in Figure 9a, and Figure 9b is 0.8086. From the scatter diagram of grain diameter (mean) after being restored by proposed method, it can be seen that the minimum value of the grain diameter (mean) is 2.868 µm, and grain area is 10.04 µm 2 , and most of the grains with diameters below 10 µm can be detected. It is proven that the grain distribution in Figure 6b is more dispersed than Figure 7b. Overall, it concluded from the experiments that, by both visual comparison and statistical assessment, the proposed method shows better restoration performance compared with some state-of-the-art methods.
Materials 2018, 11, x FOR PEER REVIEW 14 of 20 Figure 8 illustrates the histogram distribution of grain diameter, and Figure 9 shows the relationship between the grain diameter and grain area before and after being restored, respectively. From Figure 8a, it should be noted that the grain number is not counted when the grain size larger than 10 microns, because most of the grain boundaries are broken in the original damaged image, and only some continuous grains will be counted by the image pro-plus software; on the contrary, the grain can be calculated easily based upon the proposed method. From the results shown in Figure 9, the correlation value is 0.5537 in Figure 9a, and Figure 9b is 0.8086. From the scatter diagram of grain diameter (mean) after being restored by proposed method, it can be seen that the minimum value of the grain diameter (mean) is 2.868 μm, and grain area is 10.04 μm 2 , and most of the grains with diameters below 10 μm can be detected. It is proven that the grain distribution in Figure 6b is more dispersed than Figure 7b. Overall, it concluded from the experiments that, by both visual comparison and statistical assessment, the proposed method shows better restoration performance compared with some state-of-the-art methods.   Figure 8 illustrates the histogram distribution of grain diameter, and Figure 9 shows the relationship between the grain diameter and grain area before and after being restored, respectively. From Figure 8a, it should be noted that the grain number is not counted when the grain size larger than 10 microns, because most of the grain boundaries are broken in the original damaged image, and only some continuous grains will be counted by the image pro-plus software; on the contrary, the grain can be calculated easily based upon the proposed method. From the results shown in Figure 9, the correlation value is 0.5537 in Figure 9a, and Figure 9b is 0.8086. From the scatter diagram of grain diameter (mean) after being restored by proposed method, it can be seen that the minimum value of the grain diameter (mean) is 2.868 μm, and grain area is 10.04 μm 2 , and most of the grains with diameters below 10 μm can be detected. It is proven that the grain distribution in Figure 6b is more dispersed than Figure 7b. Overall, it concluded from the experiments that, by both visual comparison and statistical assessment, the proposed method shows better restoration performance compared with some state-of-the-art methods.

Experimental Case 2
For the low pixel (or low resolution) image, usually, in the engineering application, it is difficult to test the metal macro-mechanical property based on the existing grain boundary and dislocation orientation. Similar to the above experimental steps, the KSVD algorithm is applied to the DCT initial dictionary, and the training parameter setting coincided with Table 1. The DCT dictionary and trained KSVD dictionary are presented in Figure 10a and Figure 10b, respectively.   Table 5. Unfortunately, from Figure 11c and Table 5, it is obvious that the grain boundaries in recovered images are more unclear than the damaged image, based on the wavelet packet method. As shown in the close-up views of Figure 11d,e, KSVD + OMP leads to better restoration of edges and fewer artifacts than DCT + OMP, that is, because the KSVD dictionary contains more detailed structural information than the DCT dictionary. More importantly, the proposed algorithm still outperforms the other three methods; as demonstrated in Figure 11f, the number of black artifacts is reduced dramatically by the proposed method, and the grain boundaries of multi-regions become more legible.

Experimental Case 2
For the low pixel (or low resolution) image, usually, in the engineering application, it is difficult to test the metal macro-mechanical property based on the existing grain boundary and dislocation orientation. Similar to the above experimental steps, the KSVD algorithm is applied to the DCT initial dictionary, and the training parameter setting coincided with Table 1. The DCT dictionary and trained KSVD dictionary are presented in Figures 10a and 10b, respectively.

Experimental Case 2
For the low pixel (or low resolution) image, usually, in the engineering application, it is difficult to test the metal macro-mechanical property based on the existing grain boundary and dislocation orientation. Similar to the above experimental steps, the KSVD algorithm is applied to the DCT initial dictionary, and the training parameter setting coincided with Table 1. The DCT dictionary and trained KSVD dictionary are presented in Figure 10a and Figure 10b, respectively.   Table 5. Unfortunately, from Figure 11c and Table 5, it is obvious that the grain boundaries in recovered images are more unclear than the damaged image, based on the wavelet packet method. As shown in the close-up views of Figure 11d,e, KSVD + OMP leads to better restoration of edges and fewer artifacts than DCT + OMP, that is, because the KSVD dictionary contains more detailed structural information than the DCT dictionary. More importantly, the proposed algorithm still outperforms the other three methods; as demonstrated in Figure 11f, the number of black artifacts is reduced dramatically by the proposed method, and the grain boundaries of multi-regions become more legible.   Table 5. Unfortunately, from Figure 11c and Table 5, it is obvious that the grain boundaries in recovered images are more unclear than the damaged image, based on the wavelet packet method. As shown in the close-up views of Figure 11d,e, KSVD + OMP leads to better restoration of edges and fewer artifacts than DCT + OMP, that is, because the KSVD dictionary contains more detailed structural information than the DCT dictionary. More importantly, the proposed algorithm still outperforms the other three methods; as demonstrated in Figure 11f, the number of black artifacts is reduced dramatically by the proposed method, and the grain boundaries of multi-regions become more legible.  Furthermore, in order to discuss the restoration effect on grain boundary, the image segmentation and edge detection methods are applied on Figure 11b,f. Figure 12 shows the comparison of grain boundaries for the second damaged image. As illustrated in Figure 12a, the grain boundaries with small size (see the red circle) are still relatively obscure, and the grain number cannot be estimated accurately. After recovery via the proposed algorithm, the grain boundaries and grain number have a greater improvement, and are increasing. Compared with Figure 12a, it should be noted from Figure 12b that a larger number of continuous lamellar (plate-like) structures are recovered, and some breakpoints are reconnected by the proposed method, image segmentation, and edge detection method, meanwhile, the grain number will be increased via the image pro-plus software. Due to that the ideal microstructure image is unknown, comparing Figures 11 and 12, some data that may not show in ideal microstructure image might be created (thus over-correcting), this phenomenon might be improved by adjusting parameter q. Generally, for the SPSR method, the reconstruction effect and sharpness degree will be further improved when the regular operator q that is selected is much smaller, such as q = 0.1. However, the running time with q = 0.1 also became much longer. In this work, the regular operator q is selected to 0.5, as one can find in Table 5, that the running time of the proposed method (3371.2670 s) is much longer than wavelet packet (0.4969) and other sparse representation methods (DCT + OMP is 165.4376 s, and KSVD + OMP is 2928.2968 s) due to dictionary training and its iteration operations, thus, the faster calculation method based on adaptive selection  Furthermore, in order to discuss the restoration effect on grain boundary, the image segmentation and edge detection methods are applied on Figure 11b,f. Figure 12 shows the comparison of grain boundaries for the second damaged image. As illustrated in Figure 12a, the grain boundaries with small size (see the red circle) are still relatively obscure, and the grain number cannot be estimated accurately. After recovery via the proposed algorithm, the grain boundaries and grain number have a greater improvement, and are increasing. Compared with Figure 12a, it should be noted from Figure 12b that a larger number of continuous lamellar (plate-like) structures are recovered, and some breakpoints are reconnected by the proposed method, image segmentation, and edge detection method, meanwhile, the grain number will be increased via the image pro-plus software. Due to that the ideal microstructure image is unknown, comparing Figures 11 and 12, some data that may not show in ideal microstructure image might be created (thus over-correcting), this phenomenon might be improved by adjusting parameter q. Generally, for the SPSR method, the reconstruction effect and sharpness degree will be further improved when the regular operator q that is selected is much smaller, such as q = 0.1. However, the running time with q = 0.1 also became much longer. In this work, the regular operator q is selected to 0.5, as one can find in Table 5, that the running time of the proposed method (3371.2670 s) is much longer than wavelet packet (0.4969) and other sparse representation methods (DCT + OMP is 165.4376 s, and KSVD + OMP is 2928.2968 s) due to dictionary training and its iteration operations, thus, the faster calculation method based on adaptive selection method of regular operator q will be explored in future studies for eliminating the drawback of over-correcting. method of regular operator q will be explored in future studies for eliminating the drawback of overcorrecting.

Conclusions
To address the image restoration problem, a new image reconstruction technique for microstructure image based on KSVD and smoothing penalty sparse representation (SPSR) algorithm is proposed in this paper. In image sparse representation, traditional orthogonal basis functions are replaced by trained KSVD dictionary, and the trained KSVD dictionary represents the sparse characterization of block sub-image probably due to the traditional sparse representation methods that may cause instability and obvious artifacts in the reconstructed image, especially for the restoration of microstructure images, including some smooth regions, or when the noise level is strong. The proposed algorithm can overcome the above issues, which improves the reconstruction accuracy significantly. Moreover, the damaged microstructure image usually brings statistics missing in the analysis of microstructure grain parameters, and the proposed method can effectively address this drawback, which has a high engineering application value in metal materials, manufacturing, and microstructure fields.
Although the proposed method improves the reconstruction quality significantly, it still needs future improvements, where the complexity level and computational time of the proposed approach is rather high, due to dictionary training and its iteration operations. It is suggested that faster calculation methods will be explored in future studies.
Author Contributions: Conceptualization, formal analysis, investigation and writing the paper was done by Qing Li. Review and suggestion were provided by Steven Y. Liang. All authors have read and approved the final manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Conclusions
To address the image restoration problem, a new image reconstruction technique for microstructure image based on KSVD and smoothing penalty sparse representation (SPSR) algorithm is proposed in this paper. In image sparse representation, traditional orthogonal basis functions are replaced by trained KSVD dictionary, and the trained KSVD dictionary represents the sparse characterization of block sub-image probably due to the traditional sparse representation methods that may cause instability and obvious artifacts in the reconstructed image, especially for the restoration of microstructure images, including some smooth regions, or when the noise level is strong. The proposed algorithm can overcome the above issues, which improves the reconstruction accuracy significantly. Moreover, the damaged microstructure image usually brings statistics missing in the analysis of microstructure grain parameters, and the proposed method can effectively address this drawback, which has a high engineering application value in metal materials, manufacturing, and microstructure fields.
Although the proposed method improves the reconstruction quality significantly, it still needs future improvements, where the complexity level and computational time of the proposed approach is rather high, due to dictionary training and its iteration operations. It is suggested that faster calculation methods will be explored in future studies.