Statistical Analysis

While analyzing randomly rough surfaces we often need a statistical approach to determine some set of representative quantities. Within Gwyddion, there are several ways of doing this. In this section we will explain the various statistical tools and modules offered in Gwyddion, and also present the basic equations which were used to develop the algorithms they utilize.

Scanning probe microscopy data are usually represented as a two-dimensional data field of size N×M, where N and M are the number of rows and columns of the data field, respectively. The real area of the field is denoted as Lx×Ly where Lx and Ly are the dimensions along the respective axes. The sampling interval (distance between two adjacent points within the scan) is denoted Δ. We assume that the sampling interval is the same in both the x and y direction and that the surface height at a given point (xy) can be described by a random function ξ(xy) that has given statistical properties.

Note that the AFM data are usually collected as line scans along the x axis that are concatenated together to form the two-dimensional image. Therefore, the scanning speed in the x direction is considerably higher than the scanning speed in the y direction. As a result, the statistical properties of AFM data are usually collected along the x profiles as these are less affected by low frequency noise and thermal drift of the sample.

Statistical Quantities Tool

Statistical quantities include basic properties of the height values distribution, including its variance, skewness and kurtosis. The quantities accessible within Gwyddion by means of the Statistical Quantities tool are as follows:

  • Mean value, minimum, maximum and median.
  • RMS value of the height irregularities: this quantity is computed from data variance.
  • Grain-wise RMS value which differs from the ordinary RMS only if masking is used. The mean value is then determined for each grain (contiguous part of mask or inverted mask, depending on masking type) separately and the variance is then calculated from these per-grain mean values.
  • Ra value of the height irregularities: this quantity is similar to RMS value with the only difference in exponent (power) within the data variance sum. As for the RMS this exponent is q = 2, the Ra value is computed with exponent q = 1 and absolute values of the data (zero mean).
  • Height distribution skewness: computed from 3rd central moment of data values.
  • Height distribution kurtosis: computed from 4th central moment of data values.
  • Projected surface area and surface area: computed by simple triangulation.
  • Mean inclination of facets in area: computed by averaging normalized facet direction vectors.
  • Variation, which is calculated as the integral of the absolute value of the local gradient.
  • Estimated differential entropy of value distribution, calculated from the histogram as explained in the description of Entropy function. It is displayed in natural units of information (nats).

Tip

By default, the Statistical Quantities tool will display figures based on the entire image. If you would like to analyze a certain region within the image, simply click and drag a rectangle around it. The tool window will update with new numbers based on this new region. If you want you see the statistical quantities for the entire image again, just click once within the data window and the tool will reset.

More precisely, RMS (σ), skewness (γ1), and kurtosis (γ2) are computed from central moments of i-th order μi according to the following formulas:

Note that Gwyddion computes the excess kurtosis. Add 3 to obtain the area texture parameter Sku.

The surface area is estimated by the following method. Let zi for i = 1, 2, 3, 4 denote values in four neighbour points (pixel centres), and hx and hy pixel dimensions along corresponding axes. If an additional point is placed in the centre of the rectangle which corresponds to the common corner of the four pixels (using the mean value of the pixels), four triangles are formed and the surface area can be approximated by summing their areas. This leads to the following formulas for the area of one triangle (top) and the surface area of one pixel (bottom):

The method is now well-defined for inner pixels of the region. Each value participates on eight triangles, two with each of the four neighbour values. Half of each of these triangles lies in one pixel, the other half in the other pixel. By counting in the area that lies inside each pixel, the total area is defined also for grains and masked areas. It remains to define it for boundary pixels of the whole data field. We do this by virtually extending the data field with a copy of the border row of pixels on each side for the purpose of surface area calculation, thus making all pixels of interest inner.

Surface area calculation triangulation scheme (left). Application of the triangulation scheme to a three-pixel masked area (right), e.g. a grain. The small circles represent pixel-center vertices zi, thin dashed lines stand for pixel boundaries while thick lines symbolize the triangulation. The surface area estimate equals to the area covered by the mask (grey) in this scheme.

Statistical Functions Tool

One-dimensional statistical functions can be accessed by using the Statistical Functions tool. Within the tool window, you can select which function to evaluate using the selection box on the left labeled Output Type. The graph preview will update automatically. You can select in which direction to evaluate (horizontal or vertical), but as stated above, we recommend using the fast scanning axis direction. You can also select which interpolation method to use. When you are finished, click Apply to close the tool window and output a new graph window containing the statistical data.

Tip

Similar to the Statistical Quantities tool, this tool evaluates for the entire image by default, but you can select a sub-region to analyze if you wish.

Height and Angle Distribution Functions

The simplest statistical functions are the height and slope distribution functions. These can be computed as non-cumulative (i.e. densities) or cumulative. These functions are computed as normalized histograms of the height or slope (obtained as derivatives in the selected direction – horizontal or vertical) values. In other words, the quantity on the abscissa in angle distribution is the tangent of the angle, not the angle itself.

The normalization of the densities ρ(p) (where p is the corresponding quantity, height or slope) is such that

Evidently, the scale of the values is then independent on the number of data points and the number of histogram buckets. The cumulative distributions are integrals of the densities and they have values from interval [0, 1].

First-Order vs. Second-Order Quantities

The height and slope distribution quantities belong to the first-order statistical quantities, describing only the statistical properties of the individual points. However, for the complete description of the surface properties it is necessary to study higher order functions. Usually, second-order statistical quantities observing mutual relationship of two points on the surface are employed. These functions are namely the autocorrelation function, the height-height correlation function, and the power spectral density function. A description of each of these follows:

Autocorrelation Function

The autocorrelation function is given by

where z1 and z2 are the values of heights at points (x1, y1), (x2, y2); furthermore, τx =  x1 −  x2 and τy =  y1 −  y2. The function w(z1, z2, τx, τy) denotes the two-dimensional probability density of the random function ξ(xy) corresponding to points (x1, y1), (x2, y2), and the distance between these points τ.

From the discrete AFM data one can evaluate this function as

where m = τxx, n = τyy. The function can thus be evaluated in a discrete set of values of τx and τy separated by the sampling intervals Δx and Δy, respectively. The two-dimensional autocorrelation function can be calculated with Data ProcessStatistics2D Autocorrelation. Note that unlike in the one-dimensional case, it is assumed that suitable levelling has been already performed and the zero level set meaningfully. The 2D Autocorrelation function performs no levelling on its own. This gives you more control and is useful in some cases. Nevertheless, if you are more interested in the best contrast in the ACF than in defined quantitative values, you typically want to use Mean Zero Value beforehand.

For AFM measurements, we usually evaluate the one-dimensional autocorrelation function based only on profiles along the fast scanning axis. It can therefore be evaluated from the discrete AFM data values as

The one-dimensional autocorrelation function is often assumed to have the form of a Gaussian, i.e. it can be given by the following relation

where σ denotes the root mean square deviation of the heights and T denotes the autocorrelation length.

For the exponential autocorrelation function we have the following relation

Autocorrelation function obtained for simulated Gaussian randomly rough surface (i.e. with a Gaussian autocorrelation function) with σ ≈ 20 nm and T ≈ 300 nm.

We can also introduce the radial ACF Gr(τ), i.e. angularly averaged two-dimensional ACF, which of course contains the same information as the one-dimensional ACF for isotropic surfaces:

Note

For optical measurements (e. g. spectroscopic reflectometry, ellipsometry) the Gaussian autocorrelation function is usually expected to be in good agreement with the surface properties. However, some articles related with surface growth and oxidation usually assume that the exponential form is closer to the reality.

Height-Height Correlation Function

The difference between the height-height correlation function and the autocorrelation function is very small. As with the autocorrelation function, we sum the multiplication of two different values. For the autocorrelation function, these values represented the different distances between points. For the height-height correlation function, we instead use the power of difference between the points.

For AFM measurements, we usually evaluate the one-dimensional height-height correlation function based only on profiles along the fast scanning axis. It can therefore be evaluated from the discrete AFM data values as

where m = τxx. The function thus can be evaluated in a discrete set of values of τx separated by the sampling interval Δx.

The one-dimensional height-height correlation function is often assumed to be Gaussian, i.e. given by the following relation

where σ denotes the root mean square deviation of the heights and T denotes the autocorrelation length.

For the exponential height-height correlation function we have the following relation

In the following figure the height-height correlation function obtained for a simulated Gaussian surface is plotted. It is fitted using the formula shown above. The resulting values of σ and T obtained by fitting the HHCF are practically the same as for the ACF.

Height-height correlation function obtained for simulated Gaussian randomly rough surface with σ ≈ 20 nm and T ≈ 300 nm.

Power Spectral Density Function

The two-dimensional power spectral density function can be written in terms of the Fourier transform of the autocorrelation function

Similarly to the autocorrelation function, we also usually evaluate the one-dimensional power spectral density function which is given by the equation

This function can be evaluated by means of the Fast Fourier Transform as follows:

where Pj(Kx) is the Fourier coefficient of the j-th row, i.e.

If we choose the Gaussian ACF, the corresponding Gaussian relation for the PSDF is

For the surface with exponential ACF we have

In the following figure the resulting PSDF and its fit for the same surface as used in the ACF and HHCF fitting are plotted. We can see that the function can be again fitted by Gaussian PSDF. The resulting values of σ and T were practically same as those from the HHCF and ACF fit.

PSDF obtained for data simulated with Gaussian autocorrelation function.

We can also introduce the radial PSDF Wr(K), i.e. angularly integrated two-dimensional PSDF, which of course contains the same information as the one-dimensional PSDF for isotropic rough surfaces:

For a surface with Gaussian ACF this function is expressed as

while for exponential-ACF surface as

Tip

Within Gwyddion you can fit all statistical functions presented here by their Gaussian and exponential forms. To do this, fist click Apply within the Statistical Functions tool window. This will create a new graph window. With this new window selected, click on GraphFit Graph.

Minkowski Functionals

The Minkowski functionals are used to describe global geometric characteristics of structures. Two-dimensional discrete variants of volume V, surface S and connectivity (Euler-Poincaré Characteristic) χ are calculated according to the following formulas:

Here N denotes the total number of pixels, Nwhite denotes the number of white pixels, that is pixels above the threshold. Pixels below the threshold are referred to as black. Symbol Nbound denotes the number of white-black pixel boundaries. Finally, Cwhite and Cblack denote the number of continuous sets of white and black pixels respectively.

For an image with continuous set of values the functionals are parametrized by the height threshold value ϑ that divides white pixels from black, that is they can be viewed as functions of this parameter. And these functions V(ϑ), S(ϑ) and χ(ϑ) are plotted.

Range

The range distribution is a plot of the growth of value range depending on the lateral distance. It is always a non-decreasing function.

For each pixel in the image and each lateral distance, it is possible to calculate the minimum and maximum from values that do not lie father than the given lateral distance from the given pixel. The local range is the difference between the maximum and minimum (and in two dimensions can be visualised using the Rank presentation function). Averaging the local ranges for all image pixels gives the range curve.

Area Scale Graph

The area scale graph shows the ratio of developed surface area and projected area minus one, as the function of scale at which the surface area is measured. Because of the subtraction of unity, which would be the ratio for a completely flat surface, the quantity is called excess area.

Similar to the correlation functions, the area excess can be in principle defined as a directional quantity. However, the function displayed in the tool assumes an isotropic surface – the horizontal or vertical orientation that can be chosen there only determines the primary direction used in its calculation.

Masked Statistical Functions

Most of the statistical functions support masking. In other words, they can be calculated for arbitrarily shaped image regions. For some of them, such as the height or angle distribution, no further explanation is needed: only image pixels covered (or not covered) by the mask are included in the computation. However, the meaning of the function is less clear for correlation functions and spectral densities.

We will illustrate it for the ACF (the simplest case). The full description can be found in the literature [1]. The discrete ACF formula can be rewritten

where Ωm is the set of image pixels inside the image region for which the pixel m columns to the right also lies inside the region. If we calculate the ACF for a rectangular region, for example the entire image, nothing changes. However, this formula expresses the function meaningfully for regions of any shape. The only missing part is how to perform the computation efficiently.

For this we define ck,l as the mask of valid pixels, i.e. it is equal to 1 for pixels we want to include and 0 for pixels we want to exclude. It then follows that

Furthermore, if image data are premultiplied by ck,l the ACF sum over the irregular region can be replaced by the usual ACF sum over the entire image. Both are efficiently computed using FFT.

It is important to note that the amount of information available in the data about the ACF for given m depends on Ωm. Unlike for entire rectangular images, it does not have to decrease monotonically with increasing m and can vary almost arbitrarily. There can even be holes, i.e. horizontal distances m for which the region contains no pair of pixels exactly this apart. If this occurs Gwyddion replaces the missing values using linear interpolation.

Row/Column Statistics Tool

The Row/Column Statistics tool calculates numeric characteristics of each row or column and plots them as a function of its position. This makes it kind of complementary to Statistical Functions tool. Available quantities include:

  • Mean value, minimum, maximum and median.
  • RMS value of the height irregularities computed from data variance Rq.
  • Skewness and kurtosis of the height distribution.
  • Surface line length. It is estimated as the total length of the straight segments joining data values in the row (column).
  • Overall slope, i.e. the tangent of the mean line fitted through the row (column).
  • Tangent of β0. This is a characteristics of the steepnes of local slopes, closely related to the behaviour of autocorrelation and height-height correlation functions at zero. For discrete values it is calculated as follows:
  • Standard roughness parameters Ra, Rz, Rt.

In addition to the graph displaying the values for individual rows/columns, the mean value and standard deviation of the selected quantity is calculated from the set of individual row/column values. It must be emphasised that the standard deviation of the selected quantity represents the dispersion of values for individual rows/columns and cannot be interpreted as an error of the analogous two-dimensional quantity.

Two-Dimensional Slope Statistics

Several functions in Data ProcessStatistics operate on two-dimensional slope (derivative) statistics.

Slope Distribution calculates a plain two-dimensional distribution of derivatives, that is the horizontal and vertical coordinate on the resulting data field is the horizontal and vertical derivative, respectively. The slopes can be calculated as central derivatives (one-side on the borders of the image) or, if Use local plane fitting is enabled, by fitting a local plane through the neighbourhood of each point and using its gradient.

Slope Distribution can also plot summary graphs representing one-dimensional distributions of quantities related to the local slopes and facet inclination angles given by the following formula:

Three different plots are available:

  • Inclination (θ), the distribution of the inclination angle ϑ from the horizontal plane. Of course, representing the slope as an angle requires the value and the dimensions to be the same physical quantity.
  • Inclination (gradient) is similar to the ϑ plot, except the distribution of the derivative v is plotted instead of the angle.
  • Inclination (φ) visualises the integral of v2 for each direction φ in the horizontal plane. This means it is not a plain distribution of φ because areas with larger slopes contribute more significantly than flat areas.

Angle Distribution function is a visualization aid that does not calculate a distribution in the strict sense. For each derivative v the circle of points satisfying

is drawn. The number of points on the circle is given by Number of steps.

Lattice measurement

Data ProcessStatisticsMeasure Lattice

Parameters of regular two-dimensional structures can be obtained by analysing data transformed into a form that captures their regularity. The ACF and PSDF are particularly well-suited for this because they exhibit peaks corresponding to the Bravais lattice vectors for the periodic pattern.

Peaks in the two-dimensional ACF correspond to lattice vectors in the direct space (and their integer multiples and linear combinations). Peaks in the two-dimensional PSDF correspond to lattice vectors of the inverse lattice. The matrices formed by the lattice vectors in direct and frequency space are transposed inverses of each other. With a suitable transformation either can be used to measure the lattice.

The Gwyddion function Measure Lattice can utilise either ACF or PSDF for the measurement. Since regular lattices with large periods correspond to peaks in the ACF image that are far apart but peaks in the PSDF that are close, it is preferable to use the ACF in this case to improve the resolution. Conversely, the PSDF image is more suitable if the lattice periods are small. In the dialogue you can switch freely between direct and frequency space representations and the selected lattice is transformed as necessary.

If scan lines contain lots of high-frequency noise the horizontal ACF (middle row in the ACF image) can be also quite noisy. In such case it is better to interpolate it from the surrounding rows by enabling the option Interpolate horizontal ACF.

The button Estimate attempts to find automatically the lattice for the image. When it does not select the vectors as you would like you can adjust the vectors in the image manually. The selection can be done either in a manner similar to Affine Distortion (this corresponds to choosing Show lattice as Lattice), or only the two vectors can be shown and adjusted (Vectors). Either way, when you choose approximately the peak positions, pressing the Refine button adjusts the positions of the vectors to improve the match with the peaks. Button Reset can be used to restore the vectors to a sane initial state.

The lattice vectors are always displayed as direct-space vectors, with components, total lengths and directions. The angle shown as φ is the angle between the two vectors.

Facet Analysis

Data ProcessStatisticsFacet Analysis

Facet analysis enables to interactively study orientations of facets occuring in the data and mark facets of specific orientations on the image. The left view displays data with preview of marked facets. The right smaller view, called facet view below, displays two-dimensional slope distribution.

The centre of facet view always correspond to zero inclination (horizontal facets), slope in x-direction increases towards left and right border and slope in y-direction increases towards top and bottom borders. The exact coordinate system is a bit complex and it adapts to the range of slopes in the particular data displayed.

Facet plane size controls the size (radius) of plane locally fitted in each point to determine the local inclination. The special value 0 stands for no plane fitting, the local inclination is determined from symmetric x and y derivatives in each point. The choice of neighbourhood size is crucial for meaningful results: it must be smaller than the features one is interested in to avoid their smoothing, on the other hand it has to be large enough to suppress noise present in the image.

Illustration of the influence of fitted plane size on the distribution of a scan of a delaminated DLC surface with considerable fine noise. One can see the distribution is completely obscured by the noise at small plane sizes. The neighbourhood sizes are: (a) 0, (b) 2, (c) 4, (d) 7. The angle and false color mappings are full-scale for each particular image, i.e. they vary among them.

Both facet view and data view allow to select a point with mouse and read corresponding facet normal inclination value ϑ and direction φ under Normal. When you select a point on data view, the facet view selection is updated to show inclination in this point.

Button Find Maximum sets facet view selection to slope distribution maximum (the initial selection position).

Button Mark updates the mask of areas with slope similar to the selected slope. More precisely, of areas with slope within Tolerance from the selected slope. The facet view then displays the set of slopes corresponding to marked points (note the set of selected slopes may not look circular on facet view, but this is only due to selected projection). Average inclination of all points in selected range of slopes is displayed under Mean Normal.

Entropy

Data ProcessStatisticsEntropy

The Statistical Quantities tool can display an estimate of the differential entropy of the value distribution. This function can, in addition, estimate the differential entropy of slope distribution and also helps seeing how it is calculated.

The Shannon differential entropy for a probability density function can be expressed

where X is the domain of the variable x. For instance for the height distribution x represents the surface height and X is the entire real axis. For slope distribution x is a two-component vector consisting of the derivatives along the axes, and X is correspondingly the plane.

There are many more or less sophisticated entropy estimation methods. Gwyddion uses a relativaly simple histogram based method in which the above formula is approximated with

where pi and wi are the estimated probability density, i.e. normalized histogram value, and width of the i-th histogram bin.

Of course, the estimated entropy value depends on the choice of bin width. With the exception of uniform distribution, for which the estimate does not depend on the bin size. From this follows that for reasonable distributions a suitable bin width is such that the estimated entropy does not change when the width changes somewhat. This is the criterion used for choosing a suitable bin width.

In practice this means the entropy is estimated over a large range of bin widths (typically many orders of magnitude), obtaining the scaling curve displayed in the graph on the right side of the dialogue. The entropy is then estimated by finding the inflexion point of the curve. If there is no such point, which can happend for instance when the distribution is too close to a sum of δ-functions, a very large negative value is reported as the entropy.

The entropies are displayed as unitless since they are logarithmic quantities. For absolute comparison with other calculations, it should be noted that in Gwyddion they are always calculated from quantities in base SI units (e.g. metres as opposed to for instance nanometres) and natural logarithms are used. The values are thus in natural units of information (nats).

The absolute entropy depends on the absolute width of the distribution. For instance, two Gaussian distributions of heights with different RMS values have different entropies. This is useful when we want to compare the absolute narowness of structures in the height distribution. On the other hand, it can be useful to compare them in a manner independent on the overall scale. For this we can utilise the fact that the Gaussian distribution has the maximum entropy among all distributions with the same RMS value. The difference between this maximum and the estimated entropy is displayed as entropy deficit. By definition, it is always positive (apart from numerical issues).

References

[1] D. Nečas and P. Klapetek: One-dimensional autocorrelation and power spectrum density functions of irregular regions. Ultramicroscopy 124 (2013) 13-19.