EVALUATION OF WAVELET AND NON-LOCAL MEAN DENOISING OF TERRESTRIAL LASER SCANNING DATA FOR SMALL-SCALE JOINT ROUGHNESS ESTIMATION

Terrestrial Laser Scanning (TLS) is a well-known remote sensing tool that enables precise 3D acquisition of surface morphology from distances of a few meters to a few kilometres. The morphological representations obtained are important in engineering geology and rock mechanics, where surface morphology details are of particular interest in rock stability problems and engineering construction. The actual size of the discernible surface detail depends on the instrument range error (noise effect) and effective data resolution (smoothing effect). Range error can be (partly) removed by applying a denoising method. Based on the positive results from previous studies, two denoising methods, namely 2D wavelet transform (WT) and non-local mean (NLM), are tested here, with the goal of obtaining roughness estimations that are suitable in the context of rock engineering practice. Both methods are applied in two variants: conventional Discrete WT (DWT) and Stationary WT (SWT), classic NLM (NLM) and probabilistic NLM (PNLM). The noise effect and denoising performance are studied in relation to the TLS effective data resolution. Analyses are performed on the reference data acquired by a highly precise Advanced TOpometric Sensor (ATOS) on a 20x30 cm rock joint sample. Roughness ratio is computed by comparing the noisy and denoised surfaces to the original ATOS surface. The roughness ratio indicates the success of all denoising methods. Besides, it shows that SWT oversmoothes the surface and the performance of the DWT, NLM and PNLM vary with the noise level and data resolution. The noise effect becomes less prominent when data resolution decreases.


INTRODUCTION
Terrestrial Laser Scanning (TLS) for remotely acquiring precise and detailed surface morphology is of broad interest in engineering geology and rock mechanics.Some of the applications include extracting 3D geological features, assessment of weathering and material classification using relative reflectance, and modelling bare earth under vegetation (Fowler et al., 2011).This research focuses on measuring the detailed surface morphology of rock joints (discontinuities), general referred to as "roughness".Rock joint roughness is a fundamental component of overall rock strength, which is key to evaluating rock mass stability.In order to properly introduce the morphology of rock joints into stability analyses, it is crucial to represent the roughness at engineering scale of interest, which depends on the areal extent, and in direction of expected movement.Traditional disc-and-compass or contour gauge measurements of joint morphology provide limited information on roughness due to scale restrictions (no scale-dependency) and reliance on 2D measurements (no direction-dependency).TLS offers an attractive alternative, since in a short time relatively precise and detailed 3D data of an in-situ large-scale discontinuity can be obtained.Sturzenegger and Stead (2009) showed that TLS technology and careful fieldwork allow the extraction of large-scale roughness profiles.However, the small-scale roughness estimation is limited by the TLS range error, which results in overestimated roughness (noise effect), and the effective data resolution, which limits the smallest observable roughness scale (smoothing effect).
As proved in some recent studies (Bitenc et al., 2015a;Khoshelham et al., 2011;Smigiel et al., 2013Smigiel et al., , 2011Smigiel et al., , 2008) ) range error can be successfully reduced by image denoising methods.The TLS denoised surfaces show details, which may otherwise be lost in noise.In our previous study (Bitenc et al., 2015a) denoising results of Discrete Wavelet Transform (DWT) were analysed, where different DWT methods as well as thresholding procedures were tested.Smigiel et al. (2011) compared the DWT and Non-Local Mean (NLM) denoising methods for a surface reconstruction of a medium (around one cubic meter) and small (around one cubic decimetre) archaeological objects.They advocated the use of NLM.The goal here is to compare the performance of DWT and NLM in terms of providing reliable joint surface roughness estimation.Besides, results of the best denoising method are analysed with respect to the TLS effective data resolution.
The paper is organized as follows: the main object of this research, the rock joint roughness, is defined and described in Section 2. Limitations of the chosen roughness measurement method TLS are explained in Section 3 and suggestions for improvement of noisy TLS data are given in Section 4. Section 5 shows results of the methodology applied on a rock joint sample, which lead to conclusions as written in Section 6.

Definition
Rock joint roughness refers to local departures of the joint surface from planarity.For large-scale discontinuities, the joint roughness consists of large-scale (waviness or primary roughness) and small-scale (unevenness or secondary roughness) components (ISRM, 1978).Waviness can be described as a large-amplitude and low-frequency signal.Priest (1993) defined it as "surface irregularities with a wavelength greater than about 10 cm".Unevenness represents a smallamplitude and high-frequency signal.It covers finer scales of 5 to 10 cm and is superimposed on the waviness.Roughness changes with the direction (anisotropy) and should be parameterized in a direction of the most probable shear failure.When the shear direction is not known a-priori, roughness should be measured in all possible directions.

Parameterization
The angular threshold method, hereafter referred to as the Grasselli parameter (GP) after Grasselli (2001), was chosen to parameterize the roughness.The reasons are that it considers 3D surface and describes roughness anisotropy, and appears to be less sensitive to TLS measurement random errors as discclinometer method (Bitenc et al., 2015b) and fractals (Khoshelham et al., 2011).
The Grasselli parameter is based on experimental data that can identify potential rock contact areas during direct shear testing of artificial rock joints.Highly accurate and detailed ATOS (Advanced TOpometric Sensor) measurements were used to reconstruct (triangulate) the surface of the rock joint before and after the shearing.Based on joint surface damage patterns, it was found that only those areas of the joint surface that face the shear direction and are steeper than a threshold inclination θ * , which depends on the normal load, provide shear resistance.Regression analysis of the sum of those areas results in an empirical equation for the roughness, which is defined as: Where θ * is the maximum apparent dip angle of the surface in the shear (analysis) direction and C is an empirical fitting parameter calculated via a non-linear least-squares regression.More details on the roughness parameter development and calculation can be found in (Bitenc et al., 2015a;Grasselli, 2006).

TLS LIMITATIONS
When estimating surface details that are close to limits of the TLS instrument capabilities, the laser scanner range error and the effective data resolution need to be considered.In the following, range error and resolution are discussed in the context of the Rigel VZ-400 laser scanner used in our research.

Range error
In general, TLS data accuracy depends on several factors: (i) imprecision of laser scanner mechanism, (ii) geometric properties of scanned surface (scanning geometry), (iii) physical properties of scanned surface material, and (iv) environmental (atmospheric) conditions (Soudarissanane et al., 2011).
Resulting measurement errors are composed of systematic and random errors.Systematic errors are usually removed by a proper calibration procedure (Lichti, 2007).The remaining random errors are attributed mainly to range error and are therefore referred to as range noise.In the official specifications the VZ400 range error (positioning precision) is 3 mm (Riegl, 2015).As commonly holds and as shown in (Vezočnik, 2011) the actual noise is smaller.His experiment showed that for scanning distances up to 65 m and incidence angles up to 60ᵒ the noise reaches 2.2 mm.

Effective data resolution
TLS data resolution refers to the ability to resolve two objects on adjacent sight lines and depends on sampling interval (nominal point spacing) and laser beamwidth (footprint size) (Lichti and Jamtsho, 2006).In theory a very high data resolution can be achieved by repeating the scans.But the actual resolution would be lower, since the laser light diverges with the distance travelled and illuminates certain area on the surface.The range measured is the average distance to this illuminated area.As a result, at longer distances and thus larger beam-widths the surface details are averaged (smoothed) out.
The resulting (effective) TLS data resolution was studied in (Lichti and Jamtsho, 2006;Pesci et al., 2011).Based on the concept of Average Modulation Transfer Function (AMTF), which combines the sampling interval and the size of the beam footprint, Lichti and Jamtsho (2006) defined a new resolution measure called Effective Instantaneous Field Of View (EIFOV).Mills (2015) calculated the EIFOV for the VZ400.

DENOISING METHODS
Denoising methods are a promising tool to minimize the effect of TLS range error on surface roughness estimation.Once the complexity of a 3D randomly scattered point cloud is reduced by gridding to a regular 2.5D surface, a wide range of existing image processing algorithms can be used.An overview of image denoising methods and further references can be found in (Buades et al., 2005;Smigiel et al., 2011;Zhang et al. 2014).In this research the two methods (DWT and NLM), together with their two variants, are tested for a reliable roughness estimation.
In the following, the relevant details of the methods are shortly described.

Discrete wavelet transform
Denoising by discrete wavelet transform (DWT) has its origin in research of Donoho (1995).Wavelet Transforms belong to transforming domain filtering methods.DWTs decompose the surface into different scales (levels) with different space and frequency resolution by translating and dilating a single function, the mother wavelet.Since noise is characterized by high frequency fluctuations, it is likely that thresholding high frequency components (details) of DWT reduces noise and preserves low frequency components that present general trend.

Mother wavelet and level:
For a multi-level DWT, the decision has to be made about mother wavelet and number of levels.In our research the most general Daubechies wavelet (db) was chosen as the best match with the surface shape, and the decomposition level was determined by the surface size.
Considering those criteria, our experience indicates that the choice of wavelet (e.g.db3 or db4) and number of levels (e.g. 3 or 4) influence on the shape of components, but have negligible effects on denoising results.Therefore, the effects of wavelets and number of levels are not further studied here.
4.1.2Thresholding: Proper threshold selection is a key to eliminate the noise while preserving surface details.Different threshold methods have been proposed; e.g.fixed form thresholding by Donoho and Johnstone (1994), and penalized thresholding by Birgé and Massart (1997).The threshold can be applied globally (one value for all levels) or locally (one value for each level) and in a soft or hard thresholding mode.In the hard mode, coefficients that are smaller than a threshold are set to zero.In the soft mode, additional coefficients that are above the threshold are reduced for the threshold Therefore, soft thresholding results in a smoother profile and hard thresholding introduces steps.Based on our previous study (Bitenc et al., 2015a) the penalised low global hard threshold is applied.

Non-Local Mean
Non-local Mean is a widely used data-adaptive image denoising method introduced by Buades et al. (2005).Algorithms belong to pixel domain methods, where a non-local spatial filter is used to correct a noisy image.NLM is based on similarities of pixel neighbourhoods, assuming there is some redundancy i.e. selfsimilarity within an image.NLM algorithm corrects the noisy image rather than separates the noise (oscillatory) from true image (smooth), as is the case of DWT.In recent years many variants of NLM have been developed as for example Probabilistic NLM (PNLM) (Wu et al., 2013).

NLM versus PNLM:
The weight function of the classic NLM has a disadvantage of assigning non-zero weights to dissimilar patches.The problem is pointed out in (Duval et al., 2011), but no clear reasoning exist for this inadequacy.Wu et al. (2013) developed a new probabilistic weight function and showed that PNLM denoised results outperform the NLM.Besides, the PNLM works for any types of noise and the NLM assumes the Gaussian noise with zero mean and unknown variance.For those reasons and because the authors provide the Matlab code, the PNLM is tested in this research for TLS data denoising.

Data
To study the sensitivity of Grasselli roughness parameter (GP) to the noisy TLS data having certain resolution, higher precision and much denser ATOS data were used as a reference.A rock joint formed in fossiliferous limestone having dimensions of 20×30 cm (Figure 2, top) was imaged in the laboratory at approximately 0.5 m distance by ATOS I measurement system (Capture3D, 2013).On average, point density was 15 points per square millimetre (Figure 2, bottom).

Experimental workflow
The ATOS point cloud (APC) was first interpolated into a grid of 0.5 mm.This grid represents noiseless and very detailed reference surface of the sample.The original APC was too detailed to enter the Grasselli roughness computation.Further, TLS data were simulated by gridding APC into grid sizes from 3 mm to 25 mm corresponding to the Riegl VZ400 effective data resolution (EIFOV) at distances from approximately 7 m to 65 m, considering the results in (Mills, 2015).Additionally, grid sizes of 30, 40 and 50 mm were added to cover also larger distances.Each of those 11 grids (original grids) was corrupted with five levels of noise (σ) ranging from 0.5 to 2.5 mm (noisy grids).
The noisy grids were denoised by DWT, SWT, NLM and PNLM.DWT and SWT transforms were executed on three levels using a general purpose Daubechies wavelet db3 (i.e.db3 has three vanishing moments).Penalised low threshold is computed for alpha of 2. Parameters for NLM and PNLM are taken from (Wu et al., 2013), where patch size is 7 and search window 21 pixels.The weights are computed using the known noise level σ. (P)NLMs were not computed for grids larger than 10 mm, since too few image pixels were included in the "image".For DWT and SWT the surfaces of coarser resolution were extended to match the size of the 3 mm grid.The extension does not influence results of reconstructed image (Strang and Nguyen, 1997).
For the reference and each original, noisy and denoised grid the Grasselli roughness parameter (GP) was computed.The analysis direction changes clockwise from 0ᵒ (+Y-axis direction) to 355ᵒ in 5ᵒ steps (see Figure 2, top).Comparing GPs of the original grids (GP o ) to the noisy (GP n ) and denoised grids (GP d ) the Grasselli roughness Parameter ratio (GP ratio) is computed as: (2)

Noise effect
Noise effect, which is defined as the surface roughness increase with the noise, is shown in Figure 3 for the reference grid, and original and noisy grids of 5 mm.GP parameter is sensitive to noise.For the data resolution of 5 mm and the data noise of 1.5 mm (yellow curve), the GP is on average overestimated for around 73 % (GP ratio).Noise effect varies with analysis direction and has distinct maximums at around 45° and 225° and minimums at 135° and 315° (sinusoidal pattern").This pattern is connected to the ratio of grid size and noise level.It diminishes with the larger grid sizes and smaller noise.Comparing the GP of reference ATOS surface (dark red curve) with the GP o (light blue curve), the smoothing effect can be observed.

Denoising performance
Performance of the four denoising methods applied on the noisy grids of 5 mm is shown as error plot in Figure 4. Lines show the mean and the error bars the standard deviation of GPs of all analysis directions.
The error plots show that SWT results in oversmoothed surface, DWT performs best for noise levels below 2 mm and the NLM takes over in case TLS data are noisier.Same pattern of denoising performance is observed when the four methods are applied on grid sizes up to 10 mm.In Figure 5 the GP ratio of DWT denoised grids of 5 mm is shown with respect to analysis direction.As observed already in Figure 4, the DWT performs less successfully in case of 2.5 mm noise.Denoised surfaces do not show the "sinusoidal" pattern anymore (compare to Figure 3).
Figure 5. GP ratio for the denoised grids of 5 mm.

Smoothing effect
The smoothing effect, which is defined as the roughness decrease with decreasing TLS data resolution (increasing grid size or EIFOV), is presented in Figure 6.With increasing the EIFOV from 0.5 mm to 50 mm the surface is gradually smoothed out.The smoothing effect is larger in directions of small-scale roughness; in case of our sample in directions clockwise from approximately 270° to 10°.Higher roughness for larger EIFOV shows in which direction a largerscale features exist; compare the GP peaks for different EIFOV in Figure 6 with analysis directions indicated on the sample in Figure 2, top.

Synthesis
The noise and smoothing effects are presented together in Figure 7 and 8.As shown already in Figure 3, the noise has a significant influence on GP; e.g. if the EIFOV is 5 the GP increases with the five noise levels (from 0.5 mm to 2.5 mm) for 11%, 37%, 73% 101% and 135%, respectively.But the noise effect drops very fast when the EIFOV increase; e.g. if the EIFOV is 10 mm the roughness increases with the five noise levels for 5%, 23%, 45%, 64% and 85% This means, that at longer ranges, when the EIFOV increases, the noise effect gets less important.
The error bars in Figure 7 show the standard deviation of GP computed for all analysis directions.The standard deviation increases with the noise level, which means that the GP varies more with the analysis direction for higher noise levels.The standard deviation mostly increases also with the EIFOV, which indicates higher anisotropy of the sample's surface at largerscales.
Figure 8. GP ratio of noisy surfaces versus EIFOV.
The denoising performance of the DWT with the EIFOV is shown in Figure 9. Comparing the curves with the ones in Figure 8, a great improvement of GP values can be observed, especially for smaller EIFOVs and noise levels; e.g. if the EIFOV is 5 mm the GP ratio for the five noise levels (from 0.5 mm to 2.5 mm) is -3%, -3%, 1%, -2% and 28%.The graphs do not show a clear trend, which means that DWT denoising does not depend on the data resolution.However, the results for large grid sizes might not be reliable, since these contain too few grid cells (pixels).

CONCLUSIONS
In this research the usability of noisy TLS data with limited data resolution for joint surface roughness estimation is investigated.
The TLS data as acquired with Riegl VZ400 are simulated by gridding reference ATOS data and by adding noise to the grid points.In order to reduce the noise effect, the performance of wavelet transform and Non-Local Mean denoising methods is analysed.The following conclusions are obtained:  Grasselli roughness parameter is very sensitive to surface topography representation, which makes it a good tool to analyse limitations of TLS data.
 All denoising methods improve noisy surfaces (GP ratio approaches to one) and thus the capabilities of TLS for modelling smaller scale roughness.Roughness of DWT denoised surfaces is mostly underestimated.
 Further testing of DWT, NLM and other denoising methods could be beneficial, but the noise effect decreases rapidly with the resolution.Thus an improvement of denoised surface would not have much impact on the roughness.
 EIFOV controls the discernible roughness scale.Therefore it is important to know the EIFOV for the scanning distance and to adjust the resolution of raw TLS data (e.g. by gridding or wavelet transform), so it corresponds to EIFOV.
 TLS effective data resolution cannot be improved by data processing.The footprint size of a scanner depends on the scanning geometry (range and incidence angle).Therefore it is advised to scan a surface as close as possible, optimally in perpendicular direction.
Future work will involve an investigation of correlation between effective data resolution and roughness scale.Further, the EIFOV increase with distances more than 100 m should be determined, in order to be able to estimate in-situ roughness of remote rock joints.
4.1.3DWT versus SWT:A disadvantage of conventional decimated DWT is that it is not shift invariant (because of downsampling).This means that the DWT of the translated and original surfaces are not the same.Since shift-invariance is important for many applications such as change detection, denoising and pattern recognition, a new type of DWT was developed, namely non-decimated or Undecimated Wavelet Transform (UDWT).The principles of DWT and UDWT are compared in Figure1for the 1D case and three levels.DWT downsamples the signal by two, in contrast to UDWT, where the signal is left unchanged but filters are upsampled by two on each consecutive level(Fugal, 2009).

Figure 1 .
Figure 1.Downsampling the signal in case of DWT (left) and upsampling the wavelet in case of UDWT (right) (Fugal, 2009).In our research the UDWT algorithm developed by Coifman and Donoho (1995) is used, for which a Matlab implementation is available.It averages some slightly different DWTs to define Stationary Wavelet Transform (SWT).
4.2.1 NLM parameters: NLM algorithm includes three parameters, namely similarity (N) and search (S) window, and filtering parameter (h).The similarity window or central patch is a square neighbourhood from which the weight is defined as a function of Euclidian distance, i.e. locations close to query pixel to be denoised get higher value.The search window limits the area of an image within which the similar neighbourhoods are searched for.The smaller the S the shorter is the computational time, but more local the filtering.The filtering parameter controls the decay of the weights as a function of the Euclidean distance.Success of NLM denoising depends mostly on the weights assigned to the noisy pixels neighbouring the pixel to be denoised.Parameter values are data dependent and difficult to tune.For example,Ville and Kocher (2009) suggest SUREbased optimization of parameters selection.In our research we took general values for the parameters as given in(Wu et al., 2013)

Figure 2 .
Figure 2. The sample of 20×30 cm rock joint formed in fossiliferous limestone with indicted analysis direction (top) and sample's ATOS point cloud; coordinates are given in meters (bottom).

Figure 3 .
Figure 3. Changes of the noise effect with analysis direction for the noisy 5 mm grids.Reference surface and original grid of 5 mm are added for comparison.

Figure 4 .
Figure 4. Comparing performance of the four denoising methods on the noisy grids of 5 mm.