Image Reconstruction and Evaluation: Applications on Micro-surfaces and Lenna Image Representation

: This article develops algorithms for the characterization and the visualization of micro-scale features by using a small number of sample points, and with a goal to mitigate for the measurement shortcomings which are often destructive or time consuming. The popular measurement techniques which are used in imaging micro-surfaces may include 3D stylus or interferometric profilometry and Scanning Electron Microscopy (SEM), where both generate can represent the surface characteristics in terms of 3D dimensional topology and greyscale image, respectively. Such images could be highly dense; therefore, traditional image processing techniques might be computationally expensive. We implement the algorithms in several case studies to rapidly examine the microscopic features of a Microelectromechanical System (MEMS) micro-surface, and then validate the results on a famous greyscale image; i.e. “Lenna” image. The contribution of this research include first, we develop local and global algorithm based on modified Thin Plate Spline (TPS) model to reconstruct high resolution images of the micro-surface’s topography, and its derivatives by using low resolution images. Second, we obtain a bending energy algorithm from our modified TPS model, and use it to filter out image defects. Finally, we develop a computationally efficient Windowing technique, which combines TPS and Linear Sequential Estimation (LSE), to enhance the visualization of images. The Windowing technique allows rapid image reconstruction based on the reduction of inverse problem.


Introduction
The measurement of Micro-Nano features requires specialized imaging tools such as optical microscopes, thermal cameras, SEM, contact or non-contact position sensors.SEM can image conductive and non-conductive surfaces at a range of magnification from Millimeter to Nano scale; however, imaging a large surface at high resolution is time-consuming, often destructive or impractical.For example, it is possible to capture a full view of a 3D microelectromechanical (MEMS) structure in figure (1a) by using SEM at 130X magnification; but the image quality at the interfacethe electrical and mechanical interconnect between the two structures -is poor.Thus, a 597X magnification of a smaller area was captured to reveal the micron features as shown in figure (1b).The problem becomes to stitch multiple images into a single image.Toward this challenge, one of our objectives is to develop an automatic and piece-wise algorithms that can reconstruct and enhance image visualization at minimum cost function.The reconstruction of images usually involves the interpolation of a large number of intermediate pixels from a smaller number of control points, which represents a digital elevation model (DEM) with the elevation being a quantity that represents the attribute interest.Typically, control points are sample locations of a captured or sampled pixel quantity at known (x,y) coordinates.Such quantities might be the pixel intensity matrix of an image or the sampled elevations of a micro-surface.DEMs carry vital information about the feature topography, and it can be used to characterize or visualize images by using reconstruction, feature extraction, and filtering enhancement-techniques [1][2][3].There is no single technique that is suitable for all types of image enhancement.For example, low pass and median filtering techniques alone are not appropriate for high quality reconstruction of image because they deteriorate its high frequency components [4].Therefore, researchers developed a wide range of techniques that target the image's components such as Markov random field (MRF), wavelet decomposition, Gaussian filter, multiquadratic biharmonic, Shepard's, -Function Interpolation and Thin Plate Spline (TPS) [5,6].TPS algorithm is used in the representation of 3D object from 2D images.TPS requires no assumption to be made with respect to the distribution or the location of the DEM.This method became popular in the medical imaging; because the reduction of the registration from global into local allows to analyze local changes [7,8].While many researches showed the application of TPS in the visualization of macroscopic objects, a little information is known about its capability to characterize or visualize Micro and Nano structures.Toward the end of this goal, the paper expands on the application of TPS algorithm, and it provides mechanisms for imagining MEMS structures: in sections 2, we derive mathematical models based on a modified TPS algorithm for the synthesis of the topological features of micro-surfaces.The interpolation algorithm is used to reconstruct images from few scattered data points.Then, we derive characterization algorithms to describe and extract image features in terms of topological slopes, curvatures, and bending components.We develop a linear sequential estimator (LSE) based on TPS algorithm for a rapid reconstruction of image.Furthermore, we describe enhancement techniques to improve the image visualization such as image stitching, defect removal and resolution improvement.In section 3, we conduct experiments to evaluate the performance of the proposed techniques with some of the popular methods.Finally, in section 4, we test the algorithms in two illustrative examples: micro-surface of a MEMS based on cloud of sampled points, and famous Lenna image.

Modified TPS Algorithm
The TPS is a non-stochastic algorithm that was first developed to minimize the bending energy of thin-plate structures.This characteristic could be advantageous for smoothing the deformations due to large image pixilation or small surface roughness [9,10].This is because TPS has minimal bending energy properties over the entire object while it maintains the correlation with the DEMs [11].TPS algorithm uses radial basis function (RBFs) such as logarithmic RBFs and multiquadratic function that emphasizes on the local characteristics of the object, i.e., the reconstruction is largely unaffected at locations far from the RBFs.There is no unique procedure to define the RBF's; therefore, several basis functions have been developed and used in attempt to reconstruct objects from the measured DEMs [5,8,12].In this paper, we want to obtain a continuous, smooth, and differentiable model of an image or micro-surface that passes through its DEM.DEM's can be constructed with different patterns of data points as shown in figure .(2),depending on the target surface, collection procedures and techniques employed.Toward this goal, we utilize the biharmonic partial differential equation of a thin plate under point loading [3,12,13].We write where equations ( 1) and ( 2) correspond to the logarithmic RBFs of ( ) = ( ) and the boundary conditions, respectively.The radial distance is ( − ) + ( − ) . is the height or the pixel intensity level which is measured at ( , ) coordinate.a's and F's are the coefficients of a nominal plane and the weighting factors of micro-surface, respectively.= 16 / is the elastic constant at the interpolated point.and are the spring and rigidity constants, respectively.
The interpolation function ( , ) is divided into two parts [10]: (1) a linear part + + , and (2) the sum of the independent principal warp ∑ ( ) + ∑ ( ) + .In this form, the logarithmic function is not defined when value is zero.Therefore, we add a small dimensionless quantity relative to to guarantee that equation ( 1) is continuously differentiable function.This yield a modified TPS algorithm The above equation suggests that interpolated area could be either sharpened when is large or smoothen if it is small.The a's and F's coefficients in the boundary equation ( 2) are unknowns, and they can be estimated by using DEM.We describe DEM by any arbitrary N control points whose Cartesian locations and magnitude values are known.We can write an N system of linear equations by substituting the DEM in equation (3).The boundary conditions provide additional three linear equations.The system of linear equations can be written in matrix form where we define the TPS coefficient vector by and the measurement vector The measurement model matrix H becomes The inverse of H always exists except at singularity; which is either caused by dependent measurements, i.e. replication, or if the H matrix is constructed from three collinear samples.Therefore, if H is a positive definite and symmetric matrix, then we can estimate the coefficient vector from The estimation algorithmic equation (3) becomes a continuous and differentiable model whose coefficients are found from the information that is measured for the micro-surface.

Characterization of MEMS Surfaces by using Modified TPS algorithm
There are two popular types of bi-dimensional (2D to 2D ½ ) representations of microscale surfaces: SEM, and surface profilometery .SEM produces image of a sample by scanning it with a focused beam of electrons that contains information about the sample's surface topology.The produced image contains information about the planner dimension of the surface, ( , ), but it does not directly provide depth information.The surface profilometery approach uses contact or noncontact micron measurement sensing instrumentation which is coordinated in a multi-axis microresolution motion stage.The produced image is typically a 2D ½ microsurface constructed form a cloud of points of ( , , ).For each of them, the same image-based analysis can be performed [14].MEMS surfaces exhibit topographical components that range from Nano, to Micro, to Meso scale size.The non-stationary and the multi-scale nature of a MEMS surface could limit its description which is typically provided by majority of the topographical characterization method employed [4,13,15].Such surface characterization algorithms tend to provide parameterized functions that are strongly scale dependent; however, one scale parameter cannot provide sufficient information about the surface in all scales.This is because a scale parameter could conflict with multi-scale parameters as in the Permengham parameters that are used to describe 2½D micro-surfaces.These parameters were designed to characterize particular morphological features such as shape, elongation, shape complexity, and surface directionality [16].To overcome such shortcoming, we define the microsurface characteristics of MEMS surfaces by using the topological definitions that are commonly used in the description of geological landscapes; particularly, we obtain models from the derivatives of equation (3) to characterize MEMS surfaces.For example, the incremental change in micro-surface height as a point moves in lateral or horizontal directions is caused by a topography that contains peaks -hills -and valleys -dips.The angle switching between hills and valleys corresponds to changes in the continuity of micro-surface slopes.This corresponds to finding the slopes of a microsurface in the x and y directions at the interpolated points; therefore, the slopes are the first partial derivative of the reconstructed micro-surface in x and y at sample k are The second partial derivative of the reconstructed micro-surface represents the slope variation of the micro-surface, or in the geological terms: the dip variation of the geological unit [5].The surface curvatures in x, y and xy or yx are We should point out that the second partial derivatives of the reconstructed micro-surface at the reconstructed points could be used for finite-element microstructural analysis [10].
The strike of a dipping layer is commonly used to describe the geological landscapes, and it is defined by the intersection of that layer with an imaginary horizontal plane; forming a line of equal elevations [6].The strike is calculated by setting the x and y derivatives of equation (3) to zero.This gives This slope is the tangent of the angle between the positive x-axis and the strike line.The strike angle, δ, is the angle measured clockwise from the positive y-axis to the strike line, and it is a measure of the topographical contour slopes ⁄ .The strike angle, estimated in degree, becomes A vector which is normal to the micro-surface at the reconstructed point is defined by the azimuth, , and the inclination -or dip angle θ in degrees.The direction of the true dip angle is always perpendicular to the strike; therefore, azimuth measures the normality of a topographical micro-surface, and it is measured clockwise from the positive y-axis -or the north-such that the geological strike is equal to α+90 o .The azimuth is given by By definition, dips at the reconstructed points depict positive inclination that is orthogonal to the strike line; therefore, the inclination or dip angle θ must be positive and is given by The final characteristic that we obtain in this section describes the energy stored due to the deformations of micro-surface.The deformation index of a thin plate is subjected to slight bending, and it is derived for micro-surface based on the entire distributed bending energies, i.e., = ∬ + + [17].Applying this definition to equation (11)(12)(13), we can define the following bending energy index at reconstructed point The constant of proportionality, , has unit of Ј/m 2 , and it is used as a qualitative scaling factor.The TPS function, z, in equation ( 3) minimizes the nonnegative quantity of Iz; therefore, a planar micro-surface will have zero bending energy at areas where curvatures have high energy.A reconstructed topography of the bending energy is used as a qualitative measure of the geometry of a micro-surface at a location.For example, it views the negative features on a MEMS surface image such as a micromachined surface [13].

Image Processing Based on TPS algorithms
The classical image reconstruction techniques can be categorized into: intensity-based method, contour-based method, and shape-based method [18].The intensity-based method requires information about the pixel intensity values and their location within the image.We found previously that 2D micro-surface can be reconstructed from cloud of points whose location and height values are known in the Cartesian coordinate.It is possible to utilize TPS algorithm in the characterization and reconstruction of intensity-based images by replacing the z-value with intensity value.In this section, we will further develop techniques based on TPS algorithm for several image-processing enhancement, including denoising and rapid reconstruction.

TPS Based Image-Reconstruction Method
The quality of an image tends to have better visualization when: (1) the image is reconstructed at high resolution, and (2) the original pixel intensity values of the image; or the control points, is inherited in the reconstructed image.Reconstruction of an image involves the interpolation of regularly spaced intermediate pixels.Therefore, we select TPS with logarithmic RBFs algorithm, to interpolate the value of a new pixel within image region.Because the TPS coefficients in equation ( 8) is solved by using the control points, the modified TPS algorithm, which is equation ( 3), passes through these control points; therefore, the value of an interpolated pixel depends on the selected neighborhood that contains the control points.There are two types of neighborhoods: 1.A spatial neighborhood: this technique reconstructs a region within an image by using the information that is available within it, and without knowledge from other neighboring regions.The method can be applied either globally for the entire image, or locally for a selected region within the image.

Spatio-to-temporal neighborhood: this technique reconstructs an image by sequentially
updating the information from multiple neighborhoods.For example, an image could be reconstructed and warped from several regions.

TPS Based Denoising Method
Image denoising is the process of filtering out impulsive defects, and it typically requires algorithms for detection and removal of defects, and an algorithm that reconstruct the removed defects.Roman et al. [19] introduced a TPS method for denoising of old paintings with the additive noise modeled as random variable with zero-mean Gaussian distribution, and with the goal to achieve an acceptable trade-off between the noise removal and smoothing of edges.With the goal to decrease unwanted smoothing effect of TPS denoising, the proposed method by Roman uses a weighted TPS approximation of image data computed sequentially over small overlapping patches, and centered at each pixel.In this paper, however, we describe procedures of a modified TPS image denoising that can: first, detect large defects, second recover lost information, and finally obtain a smooth reconstruction of the image.The method is described as follow: 1. Construct bending energy image; which is a greyscale intensity image of the defected image, and it is computed by using the bending energy index -equation ( 18) -in the global coordinate.The impulsive defects are located in the constructed bending energy image because they are artifacts that tend to have black or white level with strong edges at their spatial boundaries [20]; therefore, a white pixel in the bending energy image corresponds to a black pixel in the original image, and it refers to a local maximum.A black pixel indicates a local minimum, and it corresponds to a white pixel in the original image.2. Obtain the location of the defects from step (1), and then remove them to create a defectfree DEM; i.e., the original pixel intensity values are free of defect.3. Calculate the TPS coefficients from the updated DEM -from step (2) -by using equation (8).There are two ways to calculate the TPS coefficients: • Global: select all the available information from the updated DEM.
• Local: select a neighborhood for each defect.4. Interpolate the intensity value of the removed data by using equation (3). 5.If required, reconstruct high-resolution image by using the method described in section 4.1.

LSE -TPS Based Warping Algorithm
In this section, we develop a LSE-TPS warping algorithm; which is a method for reconstruction of an image from multiple images.The goal is to warp multiple images into a single image, such as a torn image, while maintaining their spatial boundaries at the interface smooth.The spatial boundaries could be reconstructed at the same image resolution by using the method in subsubsection 2.3.1,where a "good" warping would be dependent on a single neighborhood.
Our objective becomes to estimate the TPS coefficients from multiple neighborhoods; particularly we want to add weighting factors and directionality among neighborhoods.This way it becomes possible to emphasize the reconstruction on selective neighborhoods.In this section, we refer neighborhood by a patch.Toward the end of this goal, we develop a Linear Sequential Estimation (LSE) algorithm that updates the TPS coefficients between successive patches.
Suppose there is a sequentially arranged patches.The application of LSE requires equal sizes of radial basis; i.e., each patch has n number of control points.Upon merging two successive patches -( ) and ( + 1) -into a single measurement model, we can write = + ; where = , = , and = . ( This measurement model relates the states vector to the measurement matrix H, and the measurement error, , which is an additive Gaussian noise.In a static image, the ′s are set to zero; therefore, the H matrix becomes symmetric; where the elements in the main diagonal are equal to zero, and the TPS coefficient in vector are not cross-correlated.The covariance of the measurement error is stationary in a single image, and it could be represented by a random diagonal variance matrix of zero mean and known standard deviation.In snapshot imaging, such as SEM, the measurement errors and the process weight matrix, , are neither correlated forward nor backward in time.Also, we assume that and are not cross correlated within an image i.e. = 0.
Our goal is to minimize the propagation of the measurement errors when the estimated coefficients are updated from one patch of data to another -window to window.We suggest that there exist an optimal coefficient vector X such that the Jacobin cost function, = , is minimum; i.e.
where the process weight matrix is diagonal, and it is given by and corresponds to patches ( ) and ( + 1), respectively.The size of each process matrix is ( + 3) × ( + 3) .The optimization of the Jacobian cost function with respect to produces minimum error.This is because the Hessian, which is the gradient of H, is positive definite; therefore, the optimization process can be written into the following covariance recursion process form Where equations (23 and 25) are the information matrices, equation ( 24) is the Kalman estimator gain matrix, and equation ( 26) is the Kalman update equation.These recursion equations estimate the TPS coefficient of patch ( + 1) by using its observation measurement and the TPS coefficient of patch ( ) ; therefore, a patch could be reconstructed based on its current information and the information available from the successive patch.The LSE-TPS method requires a-priori estimate of the coefficients.We can initialize the recursion process at k=0 by

LSE -TPS Based Rapid Reconstruction Method
The solution of the radial basis functions requires matrix inversion, which can be problematic for large numbers of control points [5,6,17].Equation ( 3) is dependent on the arithmetic operations that are required to compute the inverse of H; therefore, a large number of control points increases the number of computations.This mainly causes inaccuracy in the estimation of the TPS coefficient due to the accumulation of round-off errors.To improve the computational efficiency of the reconstruction and the visualization, we aim to reduce the size of the matrix H. Toward this goal; we develop the following Windowing technique: 1. Divide the image into L number of patches; each has equal number of control points.For example, figure (3a) contains 4-patches.2. Reconstruct each patch by using its own control points.3. Warp the patches at their spatial boundaries -for example, regions Ω and in figure 3.a-by selecting a neighboring set of control points for each boundary.Then, reconstruct the boundaries by using adequate method from subsubsection 2.3.1, or the LSE-TPS algorithm that we described in subsubsection 2.3.3.

Discussion
In this section, we test the performance of the modified TPS algorithms and compare it with popular interpolation methods used in image processing.

Performance Studies of the Modified TPS and Windwoing Techniques
We study the perfromance of the Modified TPS for a grayscale image by evalauting its noise-tosignal ratio (SNR) and the accurcy of the interpolation algorithim, which is presented in subsubsection 2.3.1.We define the accurcy by compairing the pixel value of the interpolated location with the true value of the pixel.To conduct such experiment, we extracted a 100 × 100 Pixel     We define the accumulative error between the true pixel value and the interpolated pixel at a location ( , ) by where is the number of samples, and the true values of ̃ vector are obtained from original image of100x100 Pixel 2 shown in figure (4.b).The total error based on the modified TPS with ɛ=10 -6 is calculated at column 37 and it is equals to 363.5187.It should be noted that the intermediate points suffers from small deviation from the true value as shown in figure (5) and as noted by non-zero values of accumulative error .The estimated value of a pixel at the DEM location is equal to the true value, i.e. coincident.This is shown in figure (6) which plots three overlaid surfaces in three dimension: (1) the triangulation of true pixel of image in figure (4.b), (2) the interpolation of the entire image whose pixels are evaluated at the location of each true pixel by using the Modified TPS method described in subsubsection 2.3.1, and (3) the intermediate pixels which are interpolated along column 37.
The effect of correction factor ɛ in the accuracy of the modified TPS algorithm in grayscale image is studied for a range of values with results shown in table 1.The method stays accurate for small values of ɛ whose values are comparable to the maximum pixel intensity value in grayscale image, but the method fails when ɛ becomes comparably large.
We utilize MATLAB built-in functions to compare the accuracy of the Modified TPS with a few popular interpolation methods used in image processing toolbox.
= interp2( , , , , , 'Method') function returns matrix containing elements corresponding to the elements of and and determined by interpolation within the two-dimensional function specified by matrices , , and .The 'Method' that were tested are 'nearest' or nearest neighbor interpolation, 'linear' or Linear interpolation (default), 'spline' or cubic spline interpolation, and 'cubic' or cubic interpolation, as long as data is uniformly-spaced.Otherwise, this method is the same as 'spline'.We use these methods to interpolate for the target values and then compare it with the modified TPS by using the cost function .The simulation results are summarized in table (2), and show that the modified TPS is on average 770% more accurate than other listed methods.
Another measure used to check on the quality of the interpolation method effectiveness is the relationship between SNR and the interpolation resolution.We designed a case study where a true image with 50 × 50 Pixel 2 is used to reconstruct 100x100 Pixel 2 image.The SNR values of the entire true images -50 × 50 Pixel 2 and 100 × 100 Pixel 2 -are 14.9 and 15.1, respectively.Then, we compare the values of the SNR which is obtained for each method to reconstruct 100 × 100 Pixel 2 image.table (2) shows that the nearest method has the lowest SNR among other methods and it is sufficiently equal to the SNR for the true image of 50x50 Pixel 2 .This is due to the locality nature of the "nearest" method which emphasize on the information available around the region of interpolated pixel.However, the modified TPS has highest SNR value of 16.4 which could be explained by the global nature of the applied TPS method where all DEMs in the 50 × 50 Pixel2 is weighted in the interpolation.The last performance measure of the modified TPS method is to understand how the SNR value varies with the number of the interpolated pixels in a given image.We use the source image of 50x50 Pixel 2 in figure (4.b) to increase the image resolution by using modified TPS method, and with ɛ = 10 - 6 .The interpolation is carried out to interpolate the intermediate points between pixels, and with step up resolutions as indicated in table 3. It can be concluded that the value of the calculated SNR increases as the resolution of the reconstructed image increase, however, to understand how such change is related to the size of the reconstructed image relative to the true image, we define SNR error percentage by = ( − )/ × 100%.Where is the signal-noise ratio of the reconstructed image, and is the signal-noise-ratio of the true image used in the reconstruction.The simulation results in figure (7) shows exponential convergence of .I.e. the signal-to-noise ratio of a reconstructed image has maximum limit where it reaches saturation at certain pixel resolution.In this case study, the pixel resolution that corresponds to the saturation limit could be obtained from the following suggested fit-model: * = ( + )/( + ), where is the pixel intensity expressed by ( ), and , , is a set of undetermined coefficients.The fit analysis is carried out by using MATLAB "fittool" GUI function.The coefficients , , are evaluated and equal to 12.91, −101, −7.425 .The "Goodness of fit" is evaluated by using the regression coefficient, R 2 , which is equal to 0.9999.The final value of the fit-model, * , can be found by using "limit theory", and it gives * ( → ∞) = 12.91%, i.e. ( → ∞) = 16.82.
The pixel intensity that corresponds to 1% within the range of the final value of is ~ 381 ×381 Pixel 2 .This explains that the SNR value of the reconstructed image would not change significantly at high resolution.On the other hand, the computational efficiency of the windowing technique, which was discussed in subsubsection 2.3.4,could be evaluated by comparing it with the global TPS method in terms of the number of computations required to invert H matrix.For example, the number of the computations required to invert H by using "Householder transformation algorithm" is 3/2n 3 +O(n 2 ) [21], where n is the total number of control points within an image.The Householder is a stable and efficient algorithm that is typically used for the inversion of large symmetrical matrices [19[22]; however, when an image is divided into L patches, the number of the computations becomes 3L/2(n/L) 3 +O (n 2 ); therefore, the number of the computation decreases when the number of patches increases.For example, suppose an image contains 1024 pixels.The number of computation required to inverse H matrix with L=1 and L=24 are 1.6 10 9 and 2.8 10 6 , respectively.

Previous Studies
A comparison between the bending energy methods is compared with a few popular methods in table 4. In addition, table 5 provided a comparison between the developed reconstructions techniques with a few popular fundamental methods used in image restoration.

SIFT
The method uses corner detection algorithm, image matching, affine deformation, viewpoint change, noise, and illumination changes [23].It shows stability in scale and rotation, robustness to localization error, and reduces the time and computational complexity of images [24].

Harris corner detector
The method is based on the local auto-correlation function of a signal particularly to detect gray scale image location that has large gradients in all directions at a predetermined scale.The method has good consistency on natural imagery [25].

SUSAN
The method is based on non-linear filtering, and it associate each pixel of the image a small area of similar brightness from neighboring pixels.It is used as edge finder [26].

SURF
The method is based on multiscale space theory and feature detection by using Hessian matrix, which has good performance and accuracy [23], but needs further improvement when the rotation is large [24].
This Study Our method uses density information of the extracted features, which is incorporated to regularize the overall energy function to reduce the impact from outliers.It is based on corner detection using local or global properties, which is an effective tool to implement the registration of local image parts, especially when the key points are evenly distributed across the whole image.

Spatial Filtering Method
Typically, smoothing filters are used for blurring and for noise reduction.Sharpening of the image is done by spatial differentiation which enhances edges and other discontinuities and de-emphasizes areas with slowly varying intensities [7].Example of spatial methods include Sobel Method [27] , Prewitt Method [28], and Laplacian method [29].

Genetic Algorithm
The method is a stochastic optimization algorithm that is based on four steps: selection, crossover, mutation and fitness function [30].A genetic algorithm needs less information to solve any problem as compared to conventional optimization methods.

Neural Network Method
Restoration is carried out by parameter evaluation and image construction procedure using learning BP neural network to determine the best image transform function [31].Neural Network is used to remove noise, motion blur, out-of-focus blur, and distortion caused by lowresolution images [32].

Fuzzy Method
Used to handle uncertainty to enhance contrast of the images and to remove unclearness in edges, boundaries, regions and features of image [33].This study We described modified TPS methods, which are based on the original pixel intensity values to evaluate model parameter based on logarithmic RBF's of bending surface of minimum deformation.The model can restore selected areas or entire image, and it is suitable to increase the image resolution, and restore lost information.We combined TPS and LSE in a windowing technique to improve the computation time of reconstruction.

Application Results
In this section, we implement techniques presented in section 2 into two applications.Software codes were developed in this study by using MATLAB program [34].

Characterization of Micro-surfaces Topology
In this illustrative example, we characterize the topology of a micro-surface that contains roughness and waviness components.The height variation of a MEMS surface in figure (8a) is measured at different location ( , ) , and then it is registered in DEM.A 3D MEMS profiler, "Veeco" [35] is used to scan the surface at spatial pixel resolution of 0.193 , and with a total number of control points of 479 × 639 Pixel 2 .Figure (8b) shows a 3D visualization of the pad surface plotted by using all the measured control points.The characterization of such data-rich model is computationally inefficient; therefore, a sufficiently small number of control points would be preferred.
A smaller data set of 10 × 10 Pixel 2 is uniformly sampled from the DEM at a lateral pixel resolution of 0.965μ .Figure (9a) shows a linear triangulation of 100 control points covering an area of 92.5 × 123.3 µm 2 .Equation ( 8) is applied to extract the TPS coefficients, and then the micro-surface is constructed at a smaller pixel size of 0.5μ by using equation (3). Figure (9b) shows a continuous and smooth micro-surface that passes through all the 100 control points.The variation of the slopes in the x and y coordinates are obtained from equations ( 9) and (10), and then they are plotted at 0.5μ resolution, as shown in figures (9c) and (9d), respectively.Similarly, the curvature variation in the x and y coordinates are constructed by using equations (11) and (12), with the results shown in figures (9e) and (9f), respectively.These later surfaces measure the topography fluctuation; where a spiky point implies that there is a saddle point close to right angle.
The strike angle in equation ( 15) is a measure of the topographical contour slopes.Figure (9g) shows a reconstruction of the strike surface; i.e., non-uniform fluctuation of strike angles within ∓90 range.The dips surface in figure (9h) is constructed from equation (17), and it shows positive inclination orthogonal to the strike line.The bending energy in figure (9i) is constructed from equation (18), and it shows the bent variation in the surface due to curvatures.The nodes, which are located on the infinite part of the extrapolated micro-surface, have negligible bending energies; therefore, the energy index can be used to detect the sudden changes within micro-surface such as holes, groves, and standing features in MEMS structures.

Visual Enhancement of Grayscale Lenna Image
SEM imaging suffers from the following shortcomings: (1) higher magnification leads to smaller field of view; therefore warping multiple images is needed for larger view, (2) high magnification often causes some damages to the specimen due to high electron field; therefore, there might be a need to filter out defects, and (3) high magnification could be non-feasible; therefore, there might be a need to reconstruct a higher resolution image from the lower magnification image.Because our algorithms can be applied to any greyscale SEM image, we select a standard picture from the image processing literature that is analogues to SEM images, and the goal is to guide the readers with the comparison.In this section, a noise-free Lena image [36], which is shown in figure ( This later image will be used throughout this section.An impulsive defect in form of black-box is added to figure (10b) as shown in figure (12a).The bending energy equation ( 18) is applied to detect the edges of the defect.Figure (12b) is the reconstructed bending energy image, and it shows the border of the defect cornered by white pixels; while the pixels away from defected area are registered black.We further implement the TPS methods to enhance the visualization of a defected image.A defect is added to a poor resolution image of 26 × 26 Pixel 2 as shown in figure (13a).The defect is detected by using the bending energy algorithm.The control points, which are corresponding to the defect, are removed from the image; leading to a DEM defect-free.This new DEM is used in equation (3) to reconstruct a higher resolution image of 51 × 51 Pixel 2 by using the global TPS algorithm, and the result is shown in figure (13b).Table 4. provides a comparison between the bending energy method with a few popular methods.
We investigate the windowing techniques for a 51 × 51 Pixel 2 image.A quick comparison between global and local TPS methods showed that the number of computations required to inverse H matrix of a 26 × 26 Pixel 2 are 198, 977 and 37,500 respectively; however, the boundaries between the square windows are not smooth, and they become visually sharp when the resolution increases.This is observed in figures (10b) and (10c).To mitigate such shortcoming, the rest of this section studies the application of warping by using the procedures that we suggested in sections 4.   In the last warping method, we use the LSE-TPS to construct the missing column 25 by using the information from both column 24 and 26.In this method, the local TPS coefficients of each column 24 and 26 are computed separately.A new TPS coefficients for column 25 are estimated from equations (23)(24)(25)(26)(27), then column 25 is reconstructed by using equation (3).The warping quality depends on the weighting matrix and the direction of propagation.For example, column 25 in figures (17a) and (17b) is reconstructed by propagating the coefficients from left-right and right-left, respectively.7

Conclusions
In this paper, we developed TPS and LSE techniques for the reconstruction and characterization of micro-surfaces and greyscale images.The techniques included: (1) reconstruction of an unknown micro-surface from a set of small data sample points, (2) restoration of bad samples, (3) filtering and warping of greyscale images, and (4) enhancement of surface visualization and characterization of MEMS structures.Because the techniques are obtained from the modified ordinary TPS surface, the source data, which is used to reconstruct micro-surfaces and images, are preserved.The reconstruction of high resolution images were conducted by using local or global logarithmic-based TPS method.Local TPS restoration is dependent on a subset of data close to the point of evaluation; whereas a global method is dependent on all data.The benefit of global reconstruction of microsurfaces is capitalized in its simplicity; however, the method is slow and inaccurate for large size of source data.On the other hand, we found that the windowing and warping techniques, when combined, could overcome the pitfall of the global method.Finally, we showed that the bending energy algorithm detects fault in greyscale image, and locates bent micro-surface on MEMS structures.

Figure 1 .
SEM of a 3D interconnected MEMS structures: (a) The low magnification image shows the entire structure; (b) The high magnification image shows the area of interest.

Figure 2 .
Figure 2. Sampling patterns: DEM's distributed in x and y plane.
Figure (3b) illustrate warping of multiple patches into a single image.
2 image from figure (1.a), then reduced the resolution of the extracted image by half.This is done by removing the intermediate pixels from the 100 × 100 Pixel 2 image.The results are shown in figure (4.a & 4.b).

Figure 4 .
Figure 4. Gray scale images extracted from SEM picture in figure (1.b):(a) Resolution is reduced down to 100x100 Pixel 2 ; (b) Resolution is reduced down to 50 × 50 Pixel 2 .

Figure 5 .
Figure 5.Comparison between the true values and the interpolated pixels which are evaluated between the intermediate points along column 37 of the image in figure (4.b).The true values are obtained from the DEM in column 74 of the image in figure (4.a).We use the global modified TPS method to double the resolution of the 50 × 50 Pixel 2 image in figure (4.b).The pixel intensity is based on grayscale image with a maximum value of 256.The intermediate points, which bisect the mid distances between every two adjacent pixels in column 37 of figure (4.b), is interpolated using equation (3).

( 28 )Figure 6 .
Figure 6.Three dimensional surface of the interpolated points superimposed on the true image surface.

Figure 7 .
Figure 7. Semilog plot of the SNR error percentage, , and the size of the reconstructed image based on modified TPS method.

Figure 14 .
Figure (14a) is divided into 5x5 square windows.The modified TPS algorithm is then implemented at the local coordinates of each window with the goal to increase the overall image resolution.The overall image resolution has increased from 100 × 100 Pixel 2 to 250 × 250 Pixel 2 as shown in figures (14b) and (140c), respectively.Reconstruction by using local TPS method: (a) A 50 × 50 Pixel 2 image is divided into 25 square windows; (b) Reconstructed image at resolution of 100 × 100 Pixel 2 ; (c) Reconstructed image at resolution of 250 × 250 Pixel 2 .
Figure (15a) shows an image torn into four pieces, where the white area refers to the missing control points whose spatial locations (x,y) are known.The goal is to warp the pieces into a single image by using global TPS, local TPS and LSE-TPS methods.First, in the global TPS method, the missing pixels are evaluated by using the control points from all the four pieces together.The missing pixels are constructed at same image resolution, and the result is shown in figure (15b).Second, the local TPS is used to warp the torn image in figure (16a).The control points, which are located in the left or right part of the torn image, are used separately in the reconstruction of the missing pixels that are located in column 25.Figures (16b) and (16c) are reconstructed by using the information in the left part, and right parts, respectively.

Table 1 .
Effect of correction coefficient in the accuracy of the modified TPS method.

Table 2 .
Effect of correction factor in the accuracy of modified TPS method.
**Reconstructed to resolution of 100x100 Pixel2 by using DEM size of 50x50 Pixel2.ϯ Evaluated along column 37. Ϙ Evaluated for the entire reconstructed image with a resolution of 100x100 Pixel2

Table 3 .
SNR relative to interpolation size in a modified TPS method.
* Interpolation by using DEM size of 50 × 50 Pixel2.ωEvaluated for the entire reconstructed image.

Table 4 .
Comparison with methods used in edge detection.

Table 5 .
Comparison with Fundamental methods used in image restoration.