Measurement of three-dimensional shrinkage deformations and volumes of stabilised soft clay during drying with Structure from Motion photogrammetry

The paper applies the Structure from Motion photogrammetry technique to measure three-dimensional shrinkage deformations of a cylindrical specimen during drying. A clay–binder mixture was statically compacted, and its shrinkage behaviour was observed for four weeks with a newly established photogrammetric system. The complete surface information including XYZ coordinates and RGB values was then reconstructed and formed into dense point cloud data, with which the volume change and shrinkage deformations are computed. In volume change determination, the reconstructed dense point cloud achieves a fine reproduction of the specimen surface, leading to satisfactory accuracy. The computed volume deviation was within 495.3 mm3 (0.24%) and 4170 mm3 (1.83%) for the dummy cylinder and clay specimen, respectively. Both values are significantly smaller than the volume deviations computed with a classical approach that involves artificially creating markers on the specimen surface. In the computation of shrinkage deformations, the paper proposes a novel method that detects the fast point feature histogram of the point clouds and then matches the detected descriptors, leading to an improvement in three-dimensional deformation determination. The compacted sample shrinks, and the computed axial deformations increase quasi-linearly with height. However, due to the noticeable tilt that occurred during drying, the calculated radial deformations are scattered at the same height and appear misleading. Circumferential deformations and strains are therefore adopted to represent the true material behaviour.


Introduction
Uneven drying shrinkage of clays often leads to the creation of desiccation cracks [33,47].Such cracks affect the performance of earth structures and result in practical problems including landslides and erosions.Mapping the existing shrinkage and cracks quickly and reliably is of practical importance.Moreover, to advance the theory related to clay cracking, quantifying the local deformation is vital, especially in an automatic and three-dimensional (3D) manner.Experimental investigations into the topic (i.e.drying shrinkage of compacted clays) have been extensively performed by, among others, Fleureau, Taibi and their colleagues [11-13, 19, 22, 39, 48], Marinho [29,30], Birle et al. [3], Romero et al. [34,35], Alonso et al. [1], Kodikara [18], Liu et al. [27] and Wang et al. [45].These studies aimed mainly at establishing the relationship between the physical properties of soil and suction (e.g.shrinkage-swelling, saturation-desaturation and soilwater retention curves).However, insufficient reports are available to address the full-field 3D drying shrinkage of compacted clays, which has been deemed a long-standing issue for the serviceability of dams, airfields and pavements [41].
The goal of this paper is to advance automatic deformation measurement methods to monitor the drying shrinkage of clays.The proposed method uses Structure from Motion (SfM) photogrammetry in a geotechnical laboratory.Shrinkage deformations and volumetric changes are quantified with SfM photogrammetry for a compacted clay-binder mixture along the drying path.Unlike previous studies [1, 3, 11-13, 18, 22, 27, 29, 30, 34, 35, 39, 45, 48] that aimed at establishing shrinkageswelling curves, the present research focuses on the computation of full-field 3D shrinkage deformations and volumetric changes.To achieve this goal, the 3D surface of the compacted clay is reconstructed.A novel method is introduced to detect and match feature points in the timeelapsed point clouds and compute local deformations without any predefined markers.The advance is valuable because there is no need to create delicate markers on the specimen or the latex membrane.Instead, the unique features of the specimen provide direct indices of data correspondences and avoid potential errors, e.g.resulting from specimen-membrane slippage [24].

Structure from Motion photogrammetry
Measurements of local deformations of soils were advanced in geotechnical laboratories from local strain gauges, including linear variable differential transformer [7], Hall effect semiconductor [5,6] and local deformation transducer [15], to recent micro-computed tomography [8], digital image correlation [49] and fibre Bragg gratings [20].The advent of new sensing and imaging techniques has expanded the geotechnical testing environment and object while improving measurement accuracy and speed.The close-range stereophotogrammetry, applied in triaxial tests by Salazar, Coffman and their colleagues [37,38], Zhang, Li and their colleagues [25,26,50] and more recently by Nishimura [31] and Fang et al. [9,10], characterised the 3D deformation of soil specimens under mechanical loading.
The reported stereophotogrammetry fixed multiple compact-type digital cameras around the specimen, took 2D photographs and then reconstructed the partial or complete 3D surface of the specimen.Despite the barrier (e.g.optical ray refraction) caused by the triaxial cell wall and water within the cell, the achieved accuracy and precision reached the order of 10 -3 mm for specimens of 75 mm in diameter and 150 mm in height [31].
Extending the above research, Li et al. [21] used closerange SfM photogrammetry for measurements of full-field radial deformations of compacted clays during uniaxial compression tests.SfM is a photogrammetric range imaging technique for reconstructing the 3D structure from a set of 2D images taken by a motion camera [43].Instead of using multiple cameras, SfM photogrammetry requires one camera and takes photographs at multiple viewing angles to generate the complete 3D surface of a structure.Compared to stereophotogrammetry, SfM photogrammetry exhibits more flexibility in terms of the camera position (no need to be fixed) and may include more image pixels in 3D construction software for high-detail reconstruction.This technique is used in the current study.
3 Material, testing procedure and data analysis

Malmi clay and InfraStabi80 binder
The studied material is a mixture of intact soft clay and a hydraulic binder.The clay samples were taken from a homogenous layer at depths between 5 and 5.5 m at Malmi Airport, Helsinki, with clay properties given in Fig. 1.Clay Fig. 1 Profile of the geotechnical index properties of Malmi clay (w: gravimetric water content; q: bulk density; LOI: loss on ignition; and s u : undrained shear strength) was mixed with InfraStabi80 binder produced by Ecolan Oy (Finland).This binder has the lowest carbon emission per strength among four sustainable binders recently investigated by Li et al. [21] since it consists of 80% fly ashes and 20% quick cement.The fly ashes were taken as an industrial by-product in the emission calculations, leading to the lowest carbon footprint value of the binder being 140 kg CO 2 -eq/tn in comparison with that of CEMIII cement, 470 kg CO 2 -eq/tn.
3.2 Preparation of the dummy specimen and compacted clay-binder mixture Before testing soil specimens, a dummy specimen was prepared by firstly painting a metallic cylinder with red primer, then creating an array of 150 white circular spots and finally drawing a small blue dot on top of each white spot (Fig. 2a).The dummy, with a nominal diameter of 50 mm and height of 100 mm, was tested in trial tests to examine the errors of the established photogrammetric system.Three soil specimens were also prepared by mixing intact Malmi clay with the binder (100 kg/m 3 ) and compacting the clay-binder mixture in a 50-mm-diameter PVC tube.Axial stress was applied by incrementally mounting 1, 4, 5 and 10 kg weights onto a rigid piston.The maximum consolidation stress (100 kPa) was maintained for 72 h in the 6 °C cold room.After consolidation, the densities were 1.46, 1.44 and 1.46 g/cm 3 for specimens no. 1, 2 and 3, respectively, whereas the corresponding water contents were 45.3%, 44.5% and 45.0%, respectively.Compared to clays often compacted by dynamic efforts (e.g.Proctor compaction), the surface of the specimen in the present study contained considerably larger voids.The discontinuous and uneven surfaces represented realistic features of clays under low in situ consolidation stress but resulted in challenges during marker preparation.An array of 140 circular spots was generated by spraying white primer onto a thin perforated latex membrane.A blue dot was carefully drawn onto each white spot using a fine marker pen (Fig. 2b).Due to the difference in colour, the soil surface, paints and markers were segmented by RGB thresholding.

Structure from Motion photogrammetry studio, camera settings and testing protocol
The prepared specimen was glued onto a metallic pedestal that had been screwed on a rotary table (Fig. 3).A 0.9 9 0.6 9 0.4 m 3 (inner length by width by height) polystyrene chamber covered the specimen.In the chamber, five coded circular control points were glued onto the external wall of the pedestal to facilitate the computation of the camera position.During the scan, a remotely controlled stepwise motor drove the rotary table, with the rotation speed and pause time set in the system program.Luminosity in the front side of the specimen was balanced by adjusting the light intensity of the two symmetrically distributed LED lighting sources.The background was black to reduce refraction and enhance the specimenbackground contrast.
The photographs used for the photogrammetry were taken with a Canon EOS 1D X Mark II camera equipped with a Tamron SP 90 mm fixed lens.The camera was screwed onto a rigid and height-adjustable frame.The camera settings (see Table 1) were as follows: ISO value 200, aperture F/32 and shutter speed 0.25 s.The combination of these parameters ensured a depth of field of 60.1 mm, larger than the specimen diameter (50 mm).The photographs (RAW format, 5472 9 3648 pixels) were taken at two angles of elevation (0 and 30 degrees) with each taking 120 photographs, i.e. one photograph was taken when the specimen azimuth angle changed by three degrees.The rotary system automatically controlled the specimen rotation (0.5 s) and photographing (1.5 s) to reduce the motion blur since no physical contact was allowed during photographing.Each scanning process took about 8 min, during which the shrinkage of the specimen was less than 0.1% assuming the maximum volumetric strain of 15%/day at the initial stage.The evaporation slowed down after the first day and the corresponding volume change was much less than 0.1%, which was considered negligible in this study.To validate the established SfM photogrammetric studio setup and methodology, the prepared dummy was scanned twice.
Following consolidation, the clay specimen was placed in another room with approximately 30% relative humidity at 24 °C and photographed at 0, 17, 96 and 672 h.After 672 h, the final diameter and height of the specimen were measured with a calliper.Since the specimen surface was uneven, each measurement was repeated five times and the average value was taken.The specimen was then submerged in non-wetting oil to determine its volume.The method may take into account some non-visible pores connected to the surface that could be penetrated by the liquid.To mitigate the errors, readings were taken within five seconds after submersion.Yet, the measured volume is smaller than that from a perfect reconstruction of the specimen surface.
Finally, the sample was dried at the temperature of 105 °C for 72 h to determine its water content [16].In addition to specimen #1 used for photogrammetry, two extra specimens were simultaneously prepared for the measurement of water evaporation (specimen #2) and small-strain stiffness variation (specimen #3).Those specimens were kept in the same room and tested at the same time as specimen #1.To measure the water evaporation, specimen #2 was placed on top of an analytical balance, with its weight, diameter and height measured daily in the first week and weekly afterwards.Specimen #3 was mounted on a triaxial pedestal and covered by a loading cap.The variation in small-strain stiffness (i.e.shear modulus) was measured weekly by exciting a 1 kHz shear wave from the specimen bottom and capturing the signal from the top via a pair of bender elements.To ensure stable contacts and representative results, the specimen and bender elements were kept in place during the whole drying process, with each measurement repeated three times.

3D reconstruction and displacement computation
RealityCapture software performed 3D reconstruction of the specimen from 2D photographs.It uses a sequence of algorithms, including image matching, feature extraction, triangulation and texturing.3D reconstruction of SfM photogrammetry starts with finding the camera positions (Fig. 4a), i.e. camera motion matrices (see Appendix for formulations).The next step is to generate a sparse point cloud by a trigonometric algorithm based on the control points, camera motion matrices and lens data (Fig. 4b).Then, a dense point cloud is generated by processing all the image pixels (Fig. 4c).Finally, the software maps the textures and colours, adopted from photographs, on the dense point cloud (Fig. 4d).
The textured model was exported into an XYZ-format file that contained an N 9 6 matrix, with N & 7 9 10 7 , i.e. *7 million points having XYZ coordinates and RGB values.The models must be registered before further calculation.Figure 5 shows an example of registering the model at T = 672 h (aligning cloud, A n93 matrix) with respect to the model at T = 0 h (reference cloud, R m93 matrix).A 3D affine transformation matrix was defined by comparing the five coded markers that had been glued on the pedestal (control points R 0 -R 4 and A 0 -A 4 in Fig. 5).The transformation matrix contains scaling, rotation and translation information.It was determined by enforcing the control points of T = 672 h to match the control points of T = 0 h with the least square errors.Finally, the defined transformation matrix was applied to the aligning cloud (T = 672 h), leading to a newly aligned cloud (A 0 n93 matrix) that can be compared with the reference cloud for more quantitative analysis.MATLAB script analysed the data and plotted the figures in 3D space.It also computed specimen deformations and strains based on the XYZ coordinates in the reference and compare datasets, as schematically shown in Fig. 6a, where d xyz , d xy , d z and d h denote three-dimensional, radial, axial and circumferential deformations, respectively, and e r , e a and e h denote the radial, axial and circumferential strains, respectively.All the deformations and strains have a positive sign in the calculations because the soil experienced shrinkage, which was normally assumed positive according to the sign convention in soil mechanics.
Furthermore, the MATLAB script estimated specimen volumes based on the volume approximation by the convex hull (MATLAB syntax 'Convhull') and the sum of the elementary tetrahedron volumes [25].Only the latter method is presented here due to its better accuracy than the 'Convhull'.As shown in Fig. 6b, the volume computation consists of the generation of an enclosed triangular surface mesh, inserting an arbitrary point within the point cloud, connecting this inserted point with each point on the

Structure from Motion systematic errors
Systematic measurement errors resulting from, among others, uncertainties in 3D reconstruction, marker identification and registration were examined by scanning the dummy cylinder twice (Fig. 7a, b).After comparing the geometric positions of the markers in each scan (Fig. 7c), the maximum deviation between the two scans was computed for the marker-to-marker (M2M) Euclidean distance.The maximum deviation is 95 lm, corresponding to systematic errors within three pixels.Such an error is an order of magnitude larger than that by Nishimura [31], who combined stereophotogrammetry and digital image correlation techniques to reach subpixel accuracy.In the X, Y and Z directions, the computed deviations are typically within 70 lm (Fig. 7d).
The volumes of the cylinder were computed using dense point clouds and markers.As synthesised in Table 2, volume estimations by dense point clouds are approximately fifty times more accurate than those by markers in this study.The measurement errors are 495.3 and 466.1 mm 3 (relative errors of 0.24% and 0.23%) for the former cases and -22266.5 and -22249.4mm 3 (11.01 and 11.00%) for the latter.The large error obtained from the marker measurements is due to the small number of markers on the specimen (150), leading to insufficient accuracy of the triangularisation.This effect can be theoretically quantified for the cylindric dummy whose geometry is regular.In the cross section, the triangular area is sin(pi/30) 9 cos(pi/ 30) 9 R 2 = 0.203R 2 , whereas the realistic fan-shaped area is pi/15 9 R 2 = 0.209R 2 .In the axial direction, the 150 targets covered about 91 mm of the 100-mm-height dummy.By considering both axial and radial underestimations, the markers should theoretically underestimate the volume by 11.6%, which is reasonably close to the values obtained in experiments, with errors of 11.01% and 11.0% for measurements #1 and #2, respectively.For more complex geometry (e.g.objects with irregular shapes or surface defections), establishing the mathematical correction becomes challenging.Improving the marker coverage by increasing the number of markers would likely improve the accuracy of the photogrammetric measurements of the specimen volume.According to previous experimental results by Li et al. (22), the shrinkage of compacted Cha ˆtaignier clay (liquid limit w L = 71 and similar to that of Malmi clay) was practically negligible when the suction increased from approximately 10 to 100 MPa.In this study, the suction oscillation between 121 and 244 MPa should have little impact on the volume change because major shrinkage has been theoretically finished at smaller suction values.Thus, the suction oscillation is deemed not to affect the SfM photogrammetric scan.The initial water content of the specimens (w i & 45%) corresponds to a suction of 0.81 MPa (measured by the WP4C dewpoint potentiometer), meaning that the specimens followed the drying path, and their water content decreased gradually to reach an equilibrium with the environment in terms of suction values.As the binding reaction consumes water and the water evaporates, the specimens shrink and thus generate axial and radial shrinkage strains.In Fig. 8b, the water content decreased rapidly in the first three days due to (i) chemical reactions during hydration and (ii) the large difference in suction between the specimen and the ambient environment.Water evaporation slowed down from the third day and practically stopped after the 8th day.Similarly, the height and diameter of the specimen decreased mostly in the first three days and remained constant afterwards.The  Volume calculated from the buoyancy force is considered the reference value ultimate axial strain was 4.9% and 4.7% for specimens #1 and #2, respectively, whereas the radial strain was 6.3% and 5.7%, respectively.The difference in axial and radial strains (i.e. e a \ e r ) was unexpected but seemed to result from multiple factors including anisotropy induced during sample preparation and evaporation restraint at the specimen bottom that was glued onto the pedestal.Figure 8c, d shows the results of bender element tests performed on specimen #3.Since the received signals have clear peaks, no correction is applied [44], and the peak-topeak time was used to calculate the shear wave velocity and small-strain stiffness.The specimen stiffness steadily rises after the 8 th day, while the water content might be constant (based on measurements of specimen #2).The continuous gain in stiffness in this period is difficult to explain, but one hypothesis is that it is a result of the pozzolanic reaction that typically lasts weeks for stabilised soft clays [17].

Drying shrinkage computed from markers
MATLAB script plotted the dense point clouds in 3D (Fig. 9a).The plots show realistic surface features, including coordinates and colours of the specimen at different drying times.To ensure a more meaningful analysis, the 140 marked blue dots were extracted and are separately plotted in Fig. 9b  fitting the H T0 :d z plots with least square regression (i.e.continuous lines) and are equal to the inverse of these slopes (i.e.Dd z /DH).In the 3rd diagram of Fig. 10a, radial components d xy show significant discrepancies at the same H T0 .These discrepant points appear to distribute symmetrically along the average lines, which implies the occurrence of sample tilt during drying.To increase the readability of the results, one can use circumferential deformation and strain (denoted as d h and e h , respectively), as suggested by Nishimura [31].Compared to the average radial strain, the circumferential strain is closer to the calliper measurement, i.e. dashed lines in Fig. 10b.
The mentioned tilt is highlighted by projecting the 3D deformation vectors of Fig. 9b onto the XY plane.For clarity, the projected 2D arrows are simplified by circular points that correspond to arrowheads (Fig. 11).At T = 17 h, the d x :d y plot of fourteen markers at H T0- & 10 mm (dots connected by the thin line) is not centred at the origin, which implies that shrinkage deformations are non-uniform, orientation-dependent and larger in the northwest than in the southeast.Such a phenomenon typically occurs when the specimen is heterogeneous.As the drying time T increases from 17 to 96 h, the specimen continues to shrink, and the corresponding d x :d y plots expand, with a preference toward the northwest direction.
Volume estimations of specimen #1 were achieved using the meshes in Fig. 12.The computed volume at T = 672 h is 144578.4mm 3 , remarkably smaller than those calculated from the calliper and from the buoyancy force measurements, i.e. 189099.5 and 181341.9mm 3 , respectively.Such an underestimation is not surprising mainly because the markers fail to cover the whole surface of the specimen.In addition, as mentioned earlier, 140 tracked markers are insufficient to generate dense meshes to represent the complex surface condition of the compacted specimen.To improve the volume estimation, the coverage of markers must be complete, and the mesh density should be larger.

Drying shrinkage computed from dense point clouds
140 tracked markers were used to compute local shrinkage deformations and strains, as shown in Fig. 9b.The method has been well established in geotechnical laboratory testing (e.g.[25,31]).The resulting deformation field is reasonably accurate, but the computed volume change is not, mainly because a limited number of markers does not replicate the entire 3D surface accurately enough.In addition, the creation of artificial markers, either on the latex membrane or on the specimen, is a cumbersome process.Computing the deformations and volumes without any artificial markers becomes attractive and necessary.This subsection, therefore, introduces developments based on the constructed dense point clouds.Figure 13a-c shows the computed C2C distances with downsampling (DS) factors of 10000, 100 and 1, respectively.As the point cloud size increases, the computed C2C distances become finer.With a downsampling factor of 10000, the reference and compare clouds are too sparse to ensure a realistic nearest neighbourhood, and the histogram charts of the computed C2C distances in Fig. 13a are significantly different in comparison with those in Fig. 13b, c.For the specimen at T = 672 h, the C2C distances of the top surface (approximately 5.0 mm) are close to the global measurement by calliper (i.e.4.9 mm from Fig. 10).
Figure 14 compares the local deformations computed with markers (diamonds and dashed lines) and those based on the C2C algorithm (circles and continuous lines).The comparison shows poor correlations between the two methods in terms of 3D deformations d xyz and vertical deformations d z (Fig. 14), even though the radial  In addition, the density of the point clouds affects the computed C2C distances.To make a quantitative analysis, the above calculated C2C distances were thinned by averaging the X, Y and Z coordinates of the points within each height increment (e.g. 1 mm).The profiles of d xyz , d xy and d z in Fig. 14 are therefore reduced to 110 points for downsampling factors of 1, 100 and 10000 (DS1, DS100 and DS10000, respectively).Considering the densest point cloud DS1 as the reference, one may calculate the computational errors for DS100 and DS10000.As shown in Fig. 15, deformations computed from the DS100 point cloud deviate slightly from those computed from the raw cloud, i.e. -0.1 \ (d DS100 -d) \ 0.2 mm.The deviations increase by a factor of 20 when using the DS10000 dataset, i.e. -2 \ (d DS10000 -d) \ 4 mm.For future 3D deformation computations, the well-known C2C algorithm must be used with caution, and a careful check of data correspondences and point cloud density is highly suggested.

A new method to match data and compute shrinkage deformations
An alternative to improve the C2C algorithm regarding the deformation computation is to establish a more realistic correspondence between the reference and compare point clouds.A novel idea is to establish data correspondence based on feature detection and matching.Due to the distinctive features on the specimen surface (e.g.speckles, blobs, corners and edges), data correspondence can be theoretically established by detecting and mapping these unique features.The idea has been developed in the computer vision community, where sophisticated feature descriptors, e.g.Scale-Invariant Feature Transform (SIFT) [28] and Binary Robust Independent Elementary Features (BRIEF) [4], allow accurate local feature detection for 2D images.Regarding the 3D point clouds of this study, an open-source 3D feature descriptor, named Fast Point Feature Histogram (FPFH) [36], is applied to detect features.FPFH descriptors compute and return an N-by-33 matrix, with N & 7.69 million for the reference cloud (T = 0 h) and 7.23 million for the compare cloud (T = 672 h).FPFH feature detection of the two complete clouds requires a 7.69 9 7.23 trillion (approximately 556 TB) array, which exceeds the maximum RAM by an order of three.The issue can be solved by using high-performance computing clusters or downsampling the point clouds.The latter strategy was adopted in this paper, and the dense point clouds were downsampled using a box grid filter with dimensions of 0.5 9 0.5 9 0.5 mm 3 .The grid average method (taking the grid centroid to represent each grid) is used to downsample the clouds to approximately 110 and 104 thousand points for the specimen at T = 0 h and T = 672 h, respectively (Fig. 16a).The feature detection and matching of the downsampled clouds require approximately 114 GB RAM and take approximately four minutes.As shown in Fig. 16b, approximately 3500 features are detected and matched.These features contain significant outliers that are first filtered by enforcing the height requirement, i.e. the height decreases during drying for the same feature point.Following the first filtration, approximately 2200 feature points remained and output realistic data correspondences (Fig. 16c).These filtered features still contain outliers that can be used only when the users tolerate a certain amount of inaccuracy.For higher demands, the filtered features should be further one-by-one inspected to avoid any false data mapping.In practice, this can be achieved by either developing specific algorithms or visual inspection.Figure 16d shows an example of approximately 500 filtered feature pairs that were confirmed by visual inspection.As shown in Fig. 16e, the deformations calculated from the feature-to-feature distance work well in comparison with the marker-to-marker deformations.Because the raw datasets were downsampled, the effect of the grid dimension needs to be assessed in the future as the computing force increases.

Computation of shrinkage volumes
Computing the shrinkage volume requires at first a consistent definition of volume since different techniques measure different types of volume.The buoyancy force discounted all the voids that are submissible to oil, i.e. voids inter-connected to the surface.It thus measures the 'inner' volume of the specimen.The dense point cloud neglected the voids that are visible to the camera (e.g.nearsurface voids) and counted all the other voids.The 140 markers counted all the visible and invisible voids, but the problem was that these markers did not fully cover the soil surface.The calliper method measured the realistic external volume considering both near-surface and inner voids, although there were errors related to the diameter and height measurements.In this study, the 'volume' is defined as the 'envelope' of the soil mass, and the reference is taken from the calliper measurement.Theoretically, the buoyancy force and marker techniques will underestimate the reference volume because the former measures the inner volume and the latter has insufficient coverage.
Figure 17 shows the volume computation based on dense point clouds with *7 million points per cloud.The calculated volume values, together with those measured by buoyancy force, calliper and markers, are synthesised in Table 3.As expected, the buoyancy force and marker techniques underestimate the volume by 4-5% and 18-31%, respectively.Volume calculations based on dense point clouds agree well with the physical measurements, with errors of 4170 and 2964 mm 3 (1.83 and 1.59%), respectively.

Discussion
The current study investigated the feasibility of SfM photogrammetry in measuring local shrinkage deformation and volumetric changes.The results demonstrate that the applied SfM photogrammetry achieves an accuracy of approximately three pixels (95 lm) for the studied specimen and camera settings.The volume estimation errors are within 495.3 mm 3 (0.24%) for the dummy and 3170 mm 3 (1.83%) for the compacted clay.The latter is likely overestimating the error of the proposed method because the reference volume measured by the calliper has measurement errors for the studied specimen.By applying the FPFH algorithm, the authors detected features of the reconstructed 3D point clouds and established data correspondence without predefined markers.The method seems useful for future applications of SfM photogrammetry, especially in full-field deformation computations.
Compared with classic close-range stereophotogrammetry, SfM photogrammetry requires fewer cameras (one is enough), and the camera position is not necessarily fixed during photography.These features allow SfM photogrammetry to be more applicable in large-scale field tests with the advantage of a lower budget.The measuring accuracy of SfM photogrammetry mainly depends on the pixel size of the images and on 3D reconstruction algorithms.From the theoretical point of view, the achievable accuracy should be close for stereophotogrammetry and SfM photogrammetry because both image-based techniques have similar underlying principles for deformation computations.However, one potential inconvenience of SfM photogrammetry is that the specimen must be static or quasi-static (i.e.negligible shape change) during photographing, which is not necessarily required in stereophotogrammetry.This implies that loadings should be paused during SfM photographing.
SfM photogrammetry is also comparable with another image-based technique (e.g.stereo-digital image correlation) that has gained increasing attention in geotechnical laboratory testing [2,40,46].One difference between SfM photogrammetry and stereo-digital image correlation is the reconstructed objects.SfM photogrammetry reconstructs the full-field XYZ coordinates and RGB colours of each pixel of the photographs, whereas the latter mainly reconstructs the recognised textures, e.g.soil surface features, artificially painted speckles, and so on.In stereodigital image correlation, deformation computation is based on the subsets mapped in reference and compares 3D models.Like stereophotogrammetry, stereo-digital image correlation basically requires multiple cameras to be installed around the object.Therefore, the same limitations   of stereophotogrammetry apply to stereo-digital image correlation.Besides, as Li et al. [24] discussed, the application of stereo-digital image correlation in triaxial testing is questionable due to the nonlinear distortion induced by refraction and computational inefficiency caused by the lack of an affordable post-processing program.
To improve the SfM photogrammetry 3D reconstruction accuracy, one can consider taking another set of closerange photographs in addition to the global photographs that cover the whole specimen and pedestal (Fig. 3).This can be done by reducing the subject distance or by combining the SfM and close-range stereo photogrammetry, as that introduced by Li et al. [23].Additionally, feature detection and matching processes can be optimised with better automation and fewer outliers.Algorithms regarding feature detection for 3D point clouds are scarce.The applied FPFH descriptors detect features with significant outliers that must be filtered and eliminated with particular care.In addition, for point clouds of approximately seven million points, the open-source feature-matching process requires tremendous computer RAM that exceeds the limit capacity of most current workstations.These two inconveniences need to be considered in future applications.Regarding the first issue (higher accuracy and automation in data detection and matching), future developers are suggested to take advantage of recent advances in stereo vision (e.g.functional maps by Ovsjanikov et al. [32]).Finally, SfM photogrammetry should be examined by more field tests.Future tests on large-scale objects, especially in field conditions, will bring complementary feedback.

Conclusion
While 3D image-based techniques have gained popularity in experimental mechanics, their integration into the geotechnical laboratory and field testing still faces practical challenges related to cost, specimen size and post-processing algorithms.Current research applies an economical 3D vision technique (i.e.SfM photogrammetry) in a laboratory-scale drying test.By scanning the compacted specimen with SfM photogrammetry, 3D models were reconstructed in the form of dense point clouds that allow higher accuracy (error within 1.83%) in shrinkage volume computation than the artificially marked target points.The paper introduced an open-source algorithm to detect and map feature points, based on which local deformations and strains are computed without any predefined markers.The presented SfM photogrammetry is affordable and can be used in nearly every geotechnical laboratory.Future developments of SfM photogrammetry into geotechnical testing may focus on data mapping algorithms and field tests.

Appendix: Basic formulation of SfM photogrammetry
SfM photogrammetry is a technique for recovering threedimensional scene structure and camera motion from a sequence of images.Formulation of SfM photogrammetry was initiated in the computer vision community with the orthogonal assumption [42] and served as the backbone of current SfM reconstruction algorithms.The basic principles of SfM photogrammetry, including the SfM problem, formulation and solution, are summarised below.
Supposing a cylindrical soil specimen with one tracked feature point P in the world coordinate system ð x; ŷ; ẑÞ, one image was taken with the orthogonal camera projection (Fig. 18a).In the image coordinate system ð î; ĵ; kÞ, the projection point Q has the coordinate (u, v).The SfM problem is defined as computing the 3D coordinate of the tracked feature point P(x, y, z) in the world coordinate system from its projection Q(u, v) in the image coordinate systems.
The above SfM problem can be expressed as Eq. ( 2), where P is a 3 9 1 matrix denoting the 3D coordinates of the feature point P in the world coordinate system, and C is a 3 9 1 matrix denoting the 3D coordinates of the camera centre in the world coordinate system.Geometrically, the situation in Eq. 2 can be simplified by moving the world coordinate system origin to the centroid of the soil specimen and therefore the camera origin to the specimen centre (Fig. 18b).By centroid subtraction, Eq. 2 can be rearranged into a new matrix equation (Eq.3), in which u and v are known.In practice, multiple image frames (e.g.F) are taken, and each has multiple feature points (e.g.P p with p = 1, 2, … N) (Fig. 18c).Equation (3) was expanded to Eq. ( 4), where W denotes the observation matrix, M denotes the camera motion matrix, and S denotes the scene structure matrix.The SfM problem now becomes a mathematic question, i.e. how can matrices S and M be solved from observation matrix W?
Fig. 18 Schematic drawing of object generation from world coordinate system to image coordinate system u 1;1 u 1;2 . . .u 1;N u 2;1 u 2;2 . . .u The observation matrix W is a 2F 9 N matrix, whose rank is in fact smaller than or equal to 3 according to the rank theorem.By applying the singular value decomposition and enforcing rank constraints, W can be simplified and rewritten in the following form: where U and V T are the sub-matrices (the first three columns and rows) of U and V T , respectively.The 3 9 3 diagonal matrix R can be factorised into R 1/2 9 R 1/2 and then multiplied by a unit matrix Q 9 Q -1 : The combination of Eqs. ( 4) and (6) gives the following expression: where Q can be solved by enforcing the orthonormality constraints: Since Q has nine variables, three images are required to solve Q.As Q is solved, the camera motion matrix M and the structure matrix S can be determined according to Eq. ( 6).

Fig. 4
Fig. 4 Sequences of 3D reconstruction: a camera position reconstruction; b sparse point cloud generation; c dense point cloud generation; and d dense point cloud with textures and colours

Fig. 6
Fig. 6 Basic principle of computing a deformations and b volumes [22] (d xyz : 3D deformation; d xy : radial deformation; d z : axial deformation; d h : circumferential deformation; e r : radial strain; e a : axial strain; e h : circumferential strain; and V: volume)

Fig. 7
Fig. 7 SfM reconstruction and measurement errors for the dummy: a 3D reconstruction of scan #1; b 3D reconstruction of scan #2; c deviations based on marker-to-marker Euclidean distance; and d histogram distribution of measurement errors . The tracked markers at T = 0 h are considered the reference dataset, based on which 3D displacement vectors d xyz were computed.Compared to the specimen diameter and height, 3D displacement vectors indicated by the arrows in Fig. 9b are relatively small and are scaled by a factor of 2 for better visualisation.The 3D displacement vectors d xyz in Fig. 9b are decomposed into vertical and radial components, which are plotted against the specimen height at T = 0 h, i.e.H T0 .In the 2nd diagram of Fig. 10a, the axial components d z increase quasi-linearly from 0 to 3.4, 4.4 and 5.0 mm as H T0 rises from 10 to 104 mm for T = 17, 96 and 672 h, respectively.Global axial strains can be defined by linearly

Fig. 8
Fig. 8 Variation of a ambience temperature and relative humidity, b axial strain, radial strain and water content, c shear wave and d small-strain stiffness of the clay specimen during drying

Fig. 9
Fig. 9 3D deformations of the clay specimen during drying: a reconstructed point clouds plotted in XYZ space and b extracted 140 markers with the corresponding 3D displacement vectors (vector length scaled by a factor of 2)

Fig. 10
Fig. 10 Decomposed deformations and strains at different specimen heights: a height plotted against deformations and b height plotted against strains

Fig. 11
Fig. 11 Plot of the Y component of 3D deformation against the X component for the specimen at T = 17, 96 and 672 h

Fig. 13 1 Fig. 14
Fig. 13 Computation of deformation based on C2C algorithm with the downsampling factor of a 10000; b 100 and c 1

Fig. 15
Fig. 15 Effect of the downsampling factor on the computed deformations

Fig. 16
Fig. 16 Computation of deformation based on feature detection and matching technique: a 3D reconstruction of the downsampled specimens at T = 0 and T = 672 h; b feature extraction and matching; c first outlier filtration; d second outlier filtration; and e comparison of deformations calculated from makers and features

Fig. 17
Fig. 17 Volume estimations from dense point cloud for the specimen at T = 0, 17, 96 and 672 h Fig. 17 Volume estimations from dense point cloud for the specimen at T = 0, 17, 96 and 672 h

Table 1
Parameters of the camera and lens

Table 2
Comparison of volume measurements by physical and photogrammetric methods for the dummy (Vol.: Volume; E: error; and RE: Relative error)