Joint roughness profiling using photogrammetry

We propose an automated camera setup for photogrammetric roughness analysis in the laboratory environment. The developed fast and low-cost automation setup can be very useful for tedious and laborsome manual field logging practices. The photographs are processed in MATLAB to obtain disparity maps. Coding routines for stereo photogrammetry and digital measurements are written in MATLAB. Secondly, 6 effecting factors (projecting an image onto core face, depth of field, brightness, camera-to-object to baseline distance ratio, projected image size, and occlusion) influencing noise in roughness depth maps computed by employing stereo photogrammetry are investigated. After deciding the best values that allow the lowest amount of noise, depth maps of 6 core faces are computed. Using the 3D point cloud generated, roughness profile measurements are made. Then, 8 profile measurements are made for each core face, both manually and digitally. The accuracy of the disparity maps has been verified by comparing 48 joint roughness coefficient (JRC) measurements made manually using a profile gauge. It was proved that surface roughness can be measured very fast in millimetric accuracy with an average root mean square error (RMSE) of 3.50 and mean absolute error (MAE) of 3.02 by the help of the proposed set-up and calibration.


Introduction
Joint roughness is an important parameter for the characterization of the shear strength of rock discontinuities. The joint roughness coefficient (JRC) is a widely used two-dimensional (2D) index, which is determined focusing on the rock joint profile through visual inspection, back-calculation, and empirical-formula estimation (Barton 1973;Tse and Cruden 1979). It is a 2D indicator that is extensively used to represent the roughness of a given rock joint profile. It is recommended by the International Society of Rock Mechanics (ISRM) for roughness estimation (ISRM 1978). The appearance of the discontinuity surface is compared visually with the profiles shown and the JRC value corresponding to the profile ( Fig. 1) which most closely matches that of the discontinuity surface is chosen. In the case of small-scale laboratory specimens, the scale of the surface roughness will be approximately the same as that of the profiles illustrated. However, in the field the length of the surface of interest may be several meters or even tens of meters and the JRC value must be estimated for the full scale surface. Barton and Choubey (1997) proposed 10 typical profiles ( Fig. 1) for estimating JRC by visual comparison. To reduce the bias in assigning JRC, Tse and Cruden (1979) developed two parameters that are strongly correlated to the 10 typical profiles, namely Z 2 and SF, which are based on the root mean square and the mean square of the first derivative of the profile, respectively. A series of statistical parameters have been developed to calculate JRC (Yu and Vayssade 1991;Yang et al. 2001;Tatone and Grasselli 2010;Jang et al. 2014;Li and Zhang 2015).
Surface roughness measurement methods can be grouped as contact and non-contact methods. Contact methods include profilograph (Du et al. 2009) and profile gauge (Barton and Choubey 1997). Non-contact methods can be listed as shadow profilometer (Franklin and Maerz 1988), laser imaging (Mah et al. 2012(Mah et al. , 2016, photogrammetry (Jessell et al. 1995;Cravero et al. 2001;Lee and Ahn 2004;Unal et al. 2004;Bistacchi and Griffith 2011;Sturzenegger and Stead 2009;Kim et al. 2015;Vollgger and Cruden 2016), an X-ray computed tomography (Diaz et al. 2017). Laser imaging produce higher accuracy and resolution images than photogrammetry for long range; however, photogrammetry has significant advantages in terms of economic feasibility, portability, and legibility of color image information and convenience. High-resolution images are required to obtain reliable JRC values from photogrammetry models (Haneberg 2007;Baker et al. 2008;Poropat 2008;2009;Tannant 2015;Ozturk and Saricam 2018;Saricam and Ozturk 2018;Zhang et al. 2018;Azarafza et al. 2019;2021). The controlling factors regarding the photogrammetric 3D models have been widely discussed by many studies (Fraser 1984;Poropat 2009;Dai and Lu 2010;Fooladgar et al. 2013;Dai et al.2014). The influencing factors which are categorized as camera and planning factors were discussed in detail by Kim et al. (2015). They proposed JRC error models for both laboratory and field photogrammetric surveys (Kim et al. 2015;.
Studies have been performed stereo photogrammetry using variety of algorithms by the help of commercially available software packages. These packages have procedural similarities, and their processes generally require specific positions and orientations of cameras as well as a pre-defined set of ground control points to achieve accurate 3D models.
In this study, first we propose an automated camera setup for photogrammetric roughness analysis in the laboratory environment. In this set up, the cores are placed in front of a camera slider, such that lens of the camera is normal to the  (Barton and Choubey 1997) fracture surfaces, with a pre-defined distance between them. A stationary projector is placed in a way that it projects a feature-rich image that fully covers the fracture surfaces of the cores. The camera slider is automated by using a Raspberry Pi and a stepper motor to achieve precise movements which can be very useful for field logging practices. Coding routines for controlling the motors are written in Python 3. Since two photographs are needed for stereo photogrammetry, the automated camera is configured such that the camera, which is calibrated in advance, takes two photographs, which are a few millimeters apart from each other, of each fracture surface. The photographs are then processed in MATLAB to obtain disparity maps. Coding routines for stereo photogrammetry and digital measurements are written in MATLAB (The MathWorks 2017).
Secondly, 6 affecting factors (projecting an image onto core face, depth of field, brightness, camera-to-object to baseline distance ratio, projected image size, and occlusion) influencing noise in roughness depth maps computed by employing stereo photogrammetry are investigated. After deciding the best values that allow the lowest amount of noise, depth maps of 6 core faces are computed. Using the 3D point cloud generated, roughness profile measurements are made. Then, 8 profile measurements are made for each core face, both manually and digitally. The accuracy of the disparity maps has been verified by comparing 48 JRC measurements made manually using a profile gauge.

Material and methods
Photogrammetry is a non-contact method used to obtain 3D geometry of the objects (Linder 2009). The object's image is taken from several different positions and angles. Given that the positions and parameters of the cameras are known, by finding the same pixel in different images, location of the pixel in the 3D space can be calculated.
In stereo photogrammetry, only two photographs of the object are taken. The principle of stereo photogrammetry is shown in Fig. 2. Because the position of object point P is different in left and right image planes (PL and PR), with already-known intrinsic and extrinsic parameters of the cameras, location of P in the 3D space can be calculated.
The equipment used in this study contains a digital camera (Canon EOS 7D Mark II with Canon EF-S 18-135 mm IS STM), a projector (NEC NP215 DLP), a tripod, a tripod head, a camera slider (Jieyang 120 cm Carbon Slider), a Raspberry Pi (Raspberry Pi 3 Model B +), a stepper motor (Nema 17), a stepper motor driver (A4988), a servo motor, a camera remote (Canon RC-6), AA batteries, a battery holder, a profile gauge, a breadboard, a checkerboard pattern, timing belt, and pulleys. The setup can be observed in Fig. 3 and Fig. 4.
The camera is mounted on top of the camera slider using the tripod head, which is mounted to the tripod. The camera slider is automated using the stepper motor, the timing belt, and the pulleys. The Raspberry Pi is used to control the position of the camera and to take pictures. While the stepper motor is used to move the camera sideways on the slider, the servo motor is used to press the button on the camera remote in order to take pictures. The stepper motor is driven by the A4988 driver. The breadboard is used to construct the circuits and to power the motors. Coding routines for controlling the motors are written in Python 3.
The profile gauge is utilized to manually measure roughness of core faces. The checkboard pattern is used for stereo camera calibration. Coding routines for stereo photogrammetry and digital measurements are written in MATLAB (The MathWorks 2017).

Theory and calculations
To make accurate measurements using stereo photogrammetry, camera calibration should be done precisely. With camera calibration, parameters of the camera's image sensor and its lens are estimated. The parameters are divided into two groups as intrinsic and extrinsic. The extrinsic parameters are related to the position of the camera in the world. They are composed of a rotation and a translation matrix. The intrinsic parameters contain the optical center, the focal length, and the skew coefficient. After estimating these parameters, lens distortions can be fixed, as well as measuring the size of the objects in world units.
In this study, camera calibration is done by MAT-LAB's Stereo Camera Calibrator App. The app requires a checkerboard pattern to estimate the camera parameters. Therefore, a checkboard pattern is photographed several times with the camera being on the left and the right positions. Since the setup contains only one camera which is moved from one position to another,  (Poropat 2006) it is crucial to move the camera precisely every time. A slight difference in position would make making precise measurements impossible. Although there is only one camera, there are two pictures taken to make measurements using stereo photogrammetry. Hence, from now on, the camera is referred to as left and right cameras, for the sake of simplicity. The left camera is the camera used to take a picture from the left side of the core. Likewise, the right camera is the camera used to take a picture from the right side.
Several photographs are excluded from the calibration to reduce the reprojection errors below 1 pixel. After the stereo camera parameters are estimated, two pictures are taken using the left and the right camera.
Next, using the stereo camera parameters, stereo image rectification is done. Stereo image rectification undistorts the images and projects them onto a common image plane such that the corresponding pixels are located in the same row coordinates.
Then, disparity map is computed using the rectified stereo images. The disparity map contains the information of how much a pixel is moved from one image to another, in pixels. Using the disparity map and the stereo camera parameters, a 3D point cloud is computed. The 3D point cloud contains the metric depth information of each pixel. Using the point cloud, profile measurements are made. Flowchart of this process can be observed in Fig. 5.

Computing JRC from digitized profiles
JRC computation is performed using the regression equation established by Tse and Cruden (1979), which uses Z 2 (the root mean square) coefficient (Myers 1962): where M is the number of intervals, D x is the interval size, and y i is the ith interval. This parameter changes with the interval sizes. In this study, the interval size is used as 1 mm. Therefore, the following equation is employed for JRC calculations for sampling size of 1 mm (Yu and Vayssade 1991):

Improving quality of depth maps
In photogrammetry, to find the corresponding points, the images should contain features. For example, 3D model of a white planar surface cannot be easily created using photogrammetry, since it does not contain enough features. In addition to this, photograph quality also affects the amount of noise in the disparity map. In this section, parameters causing noise are investigated. These parameters include projecting an image onto the core face, camera settings such as depth of field and brightness, ratio of baseline distance to camera-to-object distance, size and sharpness of the projected image, and occlusion. Since improper determination of each of these parameters can lead to data corruption in photogrammetry, accurate determination of these parameters is crucial. In the next sections, each of these parameters are explained in detail.

Effect of projecting an image onto core face
Since finding corresponding points in the two images requires existence of features, regions that do not have pixels with different colors cannot be matched with their corresponding point. Hence, in order to increase the number of features in the single-color regions, an image is projected onto the core surfaces. Figure 6 shows the effect of projecting an image on the noise. In this figure, each row contains the left image of the core face and the millimetric depth map. The depth map at the top right is retrieved by using the photos taken without projecting an image, while the one at the bottom right is computed for the photos of the core whose face is projected with an image. As can be seen from Fig. 6, image projection greatly reduces the noise.

Effect of depth of field
In photography, depth of field (DoF) refers to the zone where objects appear sharp. The larger the DoF, the greater the zone in which the objects appear sharp. However, focusing on just the object of interest with a small DoF results in sharper pixels compared to covering a larger zone.
To understand the effect of DoF on noise, an experiment with 6 different DoF values is conducted. Because changing DoF affects the brightness of the image, ISO settings are adjusted to keep the brightness levels of the images similar.
The DoF values used in the experiment are f5, f10, f16, f20, f25, and f32. Among these values, f32 resulted in the lowest amount of noise.

Effect of brightness
Brightness of the photograph affects the performance of finding corresponding pixels. To find the best setting for brightness, an experiment with 7 different ISO values is conducted. The ISO values are 320, 500, 1250, 1600, 2000, 2500, and 5000. The lowest amount of noise is achieved with ISO 1600.

Effect of camera-to-object to baseline distance ratio
Baseline distance is the distance between two cameras. The baseline distance (B) is kept as 60 mm for all experiments. Camera-to-object distance (D) is the distance between the camera and the closest point of the core face to the camera.
To observe the effects of D/B ratio on the noise, an experiment with 4 different D values is conducted. The distance values are 350 mm, 400 mm, 450 mm, and 500 mm, which result in B/D ratios of 5.83, 6.67, 7.5, and 8.33, respectively. The ratios 6.67 and 7.5 show the lowest noise, where the distances are 400 mm and 450 mm, respectively.

Effect of projected image size
Because image projection reduces the noise greatly, size of the projected image is also analyzed to obtain the lowest amount of noise.
Comparing the lowest zoom level to the highest zoom level, the highest zoom level has more noise. This is because

Effect of occlusion
Occlusion occurs when a point is hidden by another point such that the camera cannot see it. To test the effect of occlusion on noise, 3 photographs of a core face are taken by rotating the core face in a way that it results in occlusion. As the occlusion decreases, the noise decreases. Hence, the placement of the core in front of the camera has an important role in reducing noise.

Digitizing manual measurements
Manual roughness measurements are made by using the profile gauge. A total of 48 roughness measurements from 6 different core faces were taken. These manual measurements were digitized to find JRC values. Digitization was performed in MATLAB (The MathWorks 2017). This was done by photographing the manual profile measurements and diagnosing the profiles found in these photographs by image processing methods.
First, each received profile is recorded as a different image. For each profile image, the images in the RGB color space are converted to grayscale using the "rgb2gray" function. Then, by the thresholding method, white and white pixels in the image are masked. A threshold value of 180 was used. Thus, pixels with a value between 0 and 180 have been removed and only dark pixels are left. By the help of the "imclose" function, the gaps in the obtained profile mask are closed to prevent unwanted ruptures in the profile. Then, using the "bwboundaries" function, all boundaries in the resulting mask were found. With the assumption that the Fig. 6 Effect of projecting an image on noise (colored camera-to-object distance depth legend in mm) longest boundary belongs to the profile, the longest boundary is selected from the boundaries and the inside of the border is filled with the "imfill" function. The "bwmorph" function was then used to reduce the filled-in boundary to a thickness of 1 pixel. The shortest path between the left-most and right-most pixels of the skeleton was found in order to remove the branching from the resulting skeleton. Using the "bwdistgeodesic" function to find the shortest path, the distances from the left-most pixels to all pixels and distances from all right pixels to all pixels are calculated. These distances were then collected and the lowest distance values between the obtained distances were masked. These values Fig. 7 Images from different angles of the reconstructed core surface from the three-dimensional point cloud indicate the shortest path to the right-most pixel to the rightmost pixel by following the pixels on the skeleton. This gives the value of the manual profile in pixels. While the manual measurements were being photographed, 10-cm long lines were drawn on the paper. Using these lines, it is calculated manually how many pixels a millimeter corresponds to. The profile from the pixel unit obtained by the image processing method is converted to the millimeter unit using this manually calculated value. Thus, the manual measurement is digitized.

Depth mapping and determination of the core
First of all, two pictures of the core are rectified by using the stereo camera parameters obtained after the camera calibration, so that the depth map can be generated. Rectified photographs are converted to gray color scales using the "rgb2gray" function of MATLAB. Then, by using the disparity function, two images are obtained as a map of contradiction. The offset map calculates how many pixels the same point is between two photos. The distortion of the points close to the camera is high, and the distortion of the photos away from the camera is low. This can be explained by the same effect, i.e., the parallax effect, as the moon is always in a fixed position when the car moves, but the trees on the side of the road move too quickly. After obtaining the map, a three-dimensional point cloud is created by using the "reconstructScene" function, with the help of the contradiction map and the stereo camera parameters.
The determination of the core in the depth map obtained starts with the conversion of the depth map to a binary mask. This mask is obtained by specifying values greater than zero as "1" in the depth map and all other values as "0." Then, using the imclose function on this mask, the morphological shutdown is applied to clear the gaps on the surface of the core due to the noise in the depth map. The mask obtained by this process is applied morphologically using the imopen function. The mask obtained is a mask that can roughly mask the location of the core in the depth map. In order to more precisely determine the pixels representing the core, the edges of the depth map are found using the edge function, with the Canny edge detection algorithm. The image showing the edges is masked by a mask which can be found in the previous stage, which can roughly find the surface of the core, and only the edges below that mask are left. This masked edge image before the "imclose" function with the morphological shutdown, followed by the process of closing the holes with the function "imfill." The depth map is masked using the mask obtained as a result of these processes. The pixels obtained as a result of this masking process are pixels that represent the core in the depth map more precisely. Since these pixels may contain some noise as well as core pixels, the mask containing these pixels is morphologically decomposed with a small structural element and this mask is cleared of noises. The mask obtained from this process is a mask that removes the core from the depth map. Figure 7 was obtained by digital reconstruction of the point cloud obtained after filtering the depth map with the core mask. This image shows the core surface obtained from stereo photogrammetry in four different angles. Figure 8 shows a roughness measurement from the depth map of the core surface obtained by stereo photogrammetry. In this way, the measurements taken along the red line shown are shown in the graph below. In this graph, the Y-axis shows the depth of the surface in millimeters, and the X-axis shows its distance from the pixel-type point to the starting point.

Digital measurement of roughness
Since there are noises in the depth map, the obtained profile values are not always as in Fig. 8. The profile values include noise values such as "NaN," "Inf," "-Inf." Therefore, before the JRC calculation is made, all of the noise values in the profile are replaced by "NaN." Then, the "fillmissing" function is replaced with the average of values near the "NaN" values using the "movmedian" method. Thus, the profile values are cleared of noises. This applies to all axis values in the received profile, i.e., the X, Y, and Z axis values.
When performing the JRC calculation, the depth value taken is important. In this study, the measurement range is 1 mm. That is, the JRC is calculated by one of the values taken every 1 mm from the profile. For example, a value of 101 is taken from a 100-mm long profile. Even if the diameter of the cores is 64 mm, the length of the profile may vary, as the surfaces may be inclined. For this reason, when selecting the values to be used for the JRC calculation from the digitally acquired profile, the distance from one point to the other should be as close to 1 mm as possible. In order to achieve this, the distance from one point to another using the X and Y values of the profile points was calculated by the Euclidean method. After the measurement is taken from one point, it is ensured that the other point to be taken from the measurement is the distance to the previous point which is the closest to 1 mm. Thus, even if the surface of the core is inclined, measurement is always taken at 1-mm intervals.
After obtaining a profile from the depth map, the JRC value of this profile was calculated and the lowest JRC value that the profile could take. The JRC calculation is sensitive to the rotation of the received profile, as it is aware of the depth values between the two adjacent points. Figure 9 shows the effect of the rotation angle of the profile on the JRC. In this way, all profiles are rotated versions of the same profile. The JRC value to be obtained from this profile is 2.15, while the JRC value taken in the 15° rotated version is 13.22.
In order to minimize the errors received due to rotation, the digitally measured profile is rotated at various angles to calculate more than one JRC for one profile and the lowest JRC value between the values is assigned as the JRC value of the profile. The angle at which the profile was rotated was found using the convex hull of the profile. For each line in the convex envelope, the profile is rotated so that the line is the X-axis of the profile. Thus, from the profile, the The profiles were rotated using the rotation matrix. The angle of inclination of the convex envelope line to be taken is calculated and the negative of this angle is used as the angle of rotation. The rotation matrix is given below.
(4) = cos −sin sin cos the rotated profile was obtained by multiplying the X and Y values of the original profile by the rotation matrix.

Estimating JRC
To test the accuracy of the stereo photogrammetry in measuring roughness of core faces, 3D point clouds for 6 core faces are generated. For each face, 8 profile measurements are made, both digitally and manually. The cores can be observed in Fig. 10. Black lines on the faces of each core represent the profile lines. Cores are numbered the same as Fig. 10 Six core faces used in the laboratory tests. Front, side, and top views of the cores are in the first, second, and third column, respectively their row number. The first row contains core face 1, and the last row shows core face 6.
Photographs of the cores are taken by projecting an image onto their faces. For each pair of photographs, the baseline distance is 60 mm. For each photograph, the camera settings are configured as ISO 1600, shutter 40, f32, and focal length 50 mm. The camera-to-object distance is 400 mm.
The digital measurements are made by following the black profile lines drawn onto the core faces. The measurement method can be observed in Fig. 11.
Manual and digital JRC measurements for 6 core faces can be observed in Fig. 12. A total of 48 JRC values were used to analyze the results. The photogrammetric JRC measurements mainly overestimated the JRC values. These overestimated values were attributed to the distorted waviness in the profiles (Fig. 13).

Error analysis and discussion
The root mean squared error (RMSE) and the mean absolute error (MAE) are commonly employed to measure the accuracy of the predicted values and the observed values. As a natural measure of average error magnitude, the advantages of MAE have also been reported (Willmott and Matsuura 2005;Chai and Draxler 2014;Kim et al. 2016).
Errors, shown in Table 1, are calculated by finding normalized JRC (NJRC) (Eq. 5) dividing digitally measured JRC to manually measured JRC. RMSE (Eq. 6) values is also calculated for each core face.
where JRC m is the manually measured JRC value and JRC p is the JRC values obtained from the photogrammetry models. NJRC is a normalized value which indicates the accuracy of the photogrammetric JRC values compared to manual measurements. Even though the noise is tried to be reduced to minimum, some amount of noise exists in the metric depth maps. The noisy pixels are assigned a value using MATLAB's "fillmissing" function. In other words, for noisy pixels, instead of the actual depth, a predicted value is used. Hence, this causes some amount of discrepancy between manual and digital measurements.
The correct value of the roughness depends largely on camera calibration. Therefore, the camera calibration should be carried out very precisely, without causing any errors in the measurement. Larger sensor sizes and longer focal lengths also can achieve high-resolution images which might decrease the noise.
Another issue to note is that the camera slider system should not move due to camera movements or any other factor. A movement that occurs during the camera movement causes a high level of noise in the map of inconsistency and negatively affects the measurement results. If the movement of the camera slider system cannot be prevented, instead of using a single camera, pairs of cameras produced for stereophotogrammetry consisting of two pre-calibrated cameras can be used.

Conclusions
Roughness profiles of core faces were measured by employing close range stereo photogrammetry. In order to improve the quality of the depth maps, the effects of projecting an image onto the core faces, depth of field, brightness, camerato-object to baseline distance ratio, projected image size, and occlusion were investigated. The lowest amount of noise for baseline distance of 60 mm was achieved with 400 mm of camera-to-object distance, ISO 1600, shutter 40, f32, and focal length of 50 mm, with core face projected with an image.
Accuracy of the method was tested by comparing JRC values of 48 manually and digitally measured roughness profiles retrieved from 6 core faces. It was proved that surface roughness can be measured in millimetric accuracy with an average RMSE and 3.50 MAE of 3.02 by the help of the proposed set-up and calibration.
3D image data are recorded fast, handled efficiently, and can be used to further characterize the discontinuity compared to manual measurements which are labor intensive and limits the number of measurements that can be acquired. Since limited number of measurements can be achieved manually, this can cause biasing of results, especially where anisotropy might be observed in roughness which might lead to inaccurate estimates of joint surface frictional strength. The stereo photogrammetry analysis method set-up and calibration described here which can be very useful for field logging practices provides a means of generating 3D roughness map which generates many more data points than the time consuming sparse manual measurements.
Funding The authors would like to thank the Scientific and Technological Research Council of Turkey, TUBITAK, Grant No:116M692, for financial support.

Data availability
The datasets and the code used or analyzed during the current study are available from the corresponding author on request.