Wanqing Sun, Qinyou Deng, Ran Li^*
School of Optical-Electrical and Computer Engineering, University of Shanghai for Science and Technology, Shanghai, China.
DOI: 10.4236/jsip.2026.172002   PDF    HTML   XML   27 Downloads   166 Views   Citations

Abstract

To address the challenges of edge detection in square-wave fringe images corrupted by discrete spot noise during granular kinematics measurement, this paper proposes an edge detection framework based on gradient magnitude thresholding and applies it to the 3D reconstruction of granular beds. First, a multi-scale Perona-Malik anisotropic diffusion model is employed for denoising; by constructing an image pyramid and adaptively adjusting diffusion coefficients across scales, the method effectively suppresses spot noise while preserving fringe integrity. Second, the Otsu algorithm is utilized for automated threshold segmentation, where morphological structuring elements are adaptively refined according to fringe curvature. Finally, to circumvent the limitations of the traditional Canny algorithm—specifically edge fragmentation and hysteresis tracking failure—a gradient magnitude thresholding-based extraction method is selected. Numerical simulations on images with a curvature of 0.01 and 2000 noise points demonstrate that the proposed pipeline significantly outperforms direct Canny detection, increasing the Intersection over Union (IoU) from 0.51 to 0.71 and the F1-score from 0.67 to 0.82. Experimental validation on static granular beds yields a planar reconstruction error within one particle diameter (0.8 - 1.0 mm) and a slope error of 0.0015. By leveraging the structural priors of square-wave fringes, this method provides a reliable technical foundation for three-dimensional measurements in the field of particle kinematics.

Keywords

Granular Matter, Square-Wave Fringe, Edge Detection, Gradient Magnitude Thresholding, 3D Reconstruction

Share and Cite:

1. Introduction

Granular matter consists of systems composed of numerous discrete solid particles [1]-[2]. In the field of granular dynamics simulations, the Discrete Element Method (DEM) [3] is commonly employed to model particle motion. In experimental scenarios, Particle Image Velocimetry (PIV) [4] is typically used to reconstruct three-dimensional flow fields and validate simulation results. However, studies analyzing morphological evolution as an entry point remain scarce, whether in simulations or experimental settings. This is primarily attributed to the discrete and highly reflective nature of granular media used in experiments, which generates a large number of scattered light spots. While these high-density discrete spots are readily captured by PIV algorithms, they pose significant limitations for three-dimensional reconstruction. In recent years, structured light-based 3D reconstruction has gained considerable attention due to its high precision and ease of implementation, with applications spanning industrial inspection, biomedical measurement, autonomous driving, and other fields [5] [6].

The two most widely used structured-light algorithms are Fourier Transform Profilometry (FTP) [7] and Phase Measurement Profilometry (PMP) [8] [9]. Both methods rely on sinusoidal fringe patterns and recover 3D surface profiles through conventional feature extraction and phase-unwrapping algorithms [10]-[11]. In addition, Time-of-Flight (TOF) techniques [12] and laser-scanning methods [13] typically employ line lasers for data acquisition and are mainly applied to static surfaces with low reflectivity. In the present experiments, electroplated glass beads are used as the granular material. These beads exhibit inherently strong reflectivity, which is beneficial for PIV measurements but poses significant challenges for structured-light-based reconstruction. When sinusoidal fringe patterns are projected onto such highly reflective surfaces, phase ambiguity and abrupt phase discontinuities frequently occur. These effects cannot be effectively corrected by post-processing algorithms in either FTP or PMP [14]. Furthermore, while the granular bed consists of a large number of identical particles, stereo-vision-based methods [15] rely on feature matching between images captured by different cameras and therefore require additional artificial texture patterns with carefully designed complexity. Finally, line-laser projection also suffers from strong specular reflection. Although it is indeed an effective approach for 3D reconstruction along a single line, extending it to full-field surface reconstruction presents significant limitations. Dense laser lines lead to excessive reflections that prevent reliable image capture, whereas sparse laser lines fail to represent surface morphology accurately. In previous studies on single-point impacts [16], the impact-induced morphological evolution exhibited strong symmetry. In that work, Zheng proposed using point lasers to generate polarized laser curves projected onto the granular surface, forming clusters of red light spots from which a laser curve containing depth information of the deformed surface can be fitted. However, this method is only suitable for scenarios involving single-particle motion patterns. In more complex experimental scenarios in granular dynamics, the laser stripe cannot cover the entire observation region, and complex morphological analysis becomes inaccurate due to non-uniform fitting of the laser stripes. In these cases, line-laser techniques are sufficient to characterize the crater formation process accurately for symmetric impacts. However, dual-point impacts and other application scenarios generate more complex morphologies. Under such conditions, extending line-laser methods or sinusoidal fringe projection fails to provide accurate depth information.

In the field of image preprocessing, traditional denoising algorithms exhibit limitations when processing high-density discrete light spots. Gaussian or median filtering often blurs edges, while nonlinear filters designed for salt-and-pepper noise tend to cause horizontal adhesion between stripes when processing bright spots formed by large-diameter particles [17]. The anisotropic diffusion model (P-M model) proposed by Perona and Malik provides a potential solution for balancing noise suppression and edge preservation [18]. However, the traditional P-M model suffers from a “staircase effect” in multi-scale feature representation, requiring integration with multi-resolution analysis (such as image pyramids) to adapt to complex granular image backgrounds. Furthermore, the accuracy of edge extraction directly determines the precision of square-wave stripe reconstruction. The classical Canny operator [19] is prone to edge fragmentation in low signal-to-noise ratio environments, and its fixed high-low thresholds struggle to adapt to granular experimental scenarios with uneven illumination.

To mitigate the effects of strong reflections and phase loss, this study employs square-wave stripe patterns with lower feature complexity. Compared to sinusoidal stripes, the intensity transitions of square-wave patterns enable more direct depth information extraction. Accordingly, given the limited sensitivity of high-speed cameras to dark regions, this paper proposes an edge detection framework based on gradient magnitude thresholding. First, a multi-scale Perona-Malik anisotropic diffusion model is constructed, utilizing image pyramids to adaptively adjust diffusion coefficients, effectively eliminating scattered light spots from particles while preserving edge integrity. Second, curvature-adaptive morphological structuring elements are introduced to optimize Otsu [20] segmentation in response to stripe curvature variations. Finally, an adaptive gradient magnitude threshold edge extraction algorithm is proposed, addressing the failure of traditional Canny algorithms under complex granular noise conditions. Through simulation experiments and static granular bed reconstruction validation, the proposed method provides a suitable image preprocessing pipeline for three-dimensional granular dynamics measurement.

The structure of this paper is organized as follows: In Section 2 elaborates on the image processing methodology. In Section 3, we propose a 3D reconstruction technical framework; Section 4 validates the effectiveness of the proposed method through experimental verification; Section 5 summarizes the work and discusses future research directions.

2. Related Work of Image Preprocessing

2.1. Anisotropic Diffusion-Based Denoising

In square-wave stripe image processing, traditional linear smoothing methods such as Gaussian filtering lead to edge blurring, and light spot clusters also generate artifacts, failing to meet the required edge extraction demands. This study employs the Perona-Malik anisotropic diffusion model, whose core equation is as follows:

$\frac{\partial I}{\partial t}= div\left(c\left(\left|\nabla I\right|\right)*\nabla I\right)$ (1)

where I represents image intensity, t is the time variable, and c is the edge-sensitive diffusion coefficient function. Tailored to the characteristics of square-wave stripes, this paper selects an exponential diffusion function:

$c\left(\left|\nabla I\right|\right)={\text{e}}^{-{\left(\frac{\left|\nabla I\right|}{K}\right)}^2}$ (2)

This function sharply attenuates the diffusion coefficient when the gradient value exceeds K, effectively preserving stripe boundary features.

To address the specific noise characteristics of square-wave stripes combined with discrete light spots, an adaptive mechanism with multi-scale parameters is designed. Figure 1 illustrates the multi-scale processing workflow:

Figure 1. Multi-scale processing flowchart.

Dynamic Adjustment of Gradient Threshold K:

$K_s=K_0\cdot \left(1+\alpha \cdot s\right)$ (3)

where s is the scale level (0 for the original scale), α is the scale growth factor, and $K_0$ is the base gradient threshold. The scale optimization function for the number of iterations n is shown as follows:

$n_s=\{\begin{array}{l}{n}_0,s=0\\ \lfloor n_0\cdot {\beta }^s\rfloor ,s>0\end{array}$ (4)

where β is the decay coefficient and $n_0$ is the base number of iterations. The number of pyramid levels was set to 3 based on empirical evaluation. We found that performance plateaued at 3 levels—further increasing the level count did not yield additional improvement, while reducing it below 3 resulted in degraded performance.

Edge-aware diffusion is a key technique in multi-scale image processing, addressing the detail loss or noise residue problems caused by the “one-size-fits-all” approach of traditional methods in edge and smooth regions. The mathematical expression for gradient magnitude calculation is:

$G\left(x,y\right)=\sqrt{{\left(\frac{\partial I}{\partial x}\right)}^2+{\left(\frac{\partial I}{\partial y}\right)}^2}$ (5)

Subsequently, the output is obtained through a fusion equation using an edge-aware weighting function, expressed as:

$I_{\text{fuse}}=w\ast I_{\text{fine}}+\left(1-w\right)\ast I_{\text{coarse}}$ (6)

where w is the edge weight, $I_{\text{fine}}$ is the diffusion result at the current scale, and $I_{\text{coarse}}$ is the upsampled result from the coarsest scale. This design ensures that fine-scale information is primarily used at stripe edges, while coarse-scale information is mainly used in light spot regions, achieving seamless transition through adaptive weighting.

In multi-scale denoising algorithms, α and β are two key parameters controlling behavior at different scales. This paper evaluates parameters and different denoising algorithms using PSNR (Peak Signal-to-Noise Ratio). PSNR is one of the classic objective quality evaluation metrics in digital image processing. Its core concept is to measure the degree of distortion by calculating the Mean Square Error (MSE) between the original and distorted images, then quantitatively evaluate through the ratio of the signal’s maximum power to the noise power. For grayscale images I and K of size m × n, the mathematical definition of MSE is as follows:

$\text{MSE}=\frac{1}{mn}{\displaystyle \sum }_{i=0}^{m-1}{\displaystyle \sum }_{j=0}^{n-1}{\left[I\left(i,j\right)-K\left(i,j\right)\right]}^2$ (7)

The PSNR calculation formula is shown as (8):

$\text{PSNR}=20\cdot {\log }_{10}\left(\frac{{\text{MAX}}_I^2}{\sqrt{\text{MSE}}}\right)$ (8)

where ${\text{MAX}}_I^2$ is the signal power and $\sqrt{\text{MSE}}$ is the noise power.

Compared to other traditional denoising algorithms, the evaluation metrics results for the denoising algorithm in this paper are shown in Table 1. The Structural Similarity Index (SSIM) has been added as an evaluation metric, which measures the similarity between two images. By comparing the luminance, contrast, and structure dimensions of uncompressed distortion-free images with distorted images to construct a mathematical model, it can better reflect human visual perception. As shown in Table 1, the proposed algorithm demonstrates good performance across different evaluation systems.

Table 1. Comparison of the performance of different denoising algorithms.

Method	PSNR (dB)	SSIM	Edge retention
Gaussian filtering	13.33	0.4438	0.1846
Bilateral filtering traditional	12.54	0.5509	0.1950
Anisotropic diffusion	12.12	0.5671	0.1967
Algorithm of this paper	12.94	0.5574	0.2166

2.2. Image Binarization and Morphological Processing

As a classic image binarization algorithm, the Otsu algorithm still possesses certain advantages. Otsu, as an algorithm for automatically determining the binarization segmentation threshold, was proposed by Japanese scholar Otsu in 1979. The algorithm selects the threshold that maximizes the between-class variance as the optimal segmentation threshold for the image by calculating all possible thresholds. Maximization of the between-class variance means maximization of the difference between foreground and background, thereby achieving effective image segmentation. As an unsupervised automatic threshold selection method, for images containing two classes (foreground and background), the optimal threshold t* satisfies:

$t^{*}={\text{argmax}}_t{\sigma }_B^2\left(t\right)$ (9)

where the between-class variance is defined as:

${\sigma }_B\left(t\right)= w_0\left(t\right)\cdot w_1\left(t\right)\cdot {\left[{\mu }_0\left(t\right)-{\mu }_1\left(t\right)\right]}^2$ (10)

where $w_0$ and $w_1$ are the weights of the background and foreground respectively, ${\mu }_0$ and ${\mu }_1$ are their respective gray-level means. The advantage of OTSU lies in its ability to automatically determine the threshold without manual setting, and it performs best for images with bimodal distributions.

After binarization processing, the stripe regions in the image are basically processed, but isolated noise points still exist. Therefore, this paper selects a symmetric structuring element $S_{\text{comb}}$ , whose shape is designed as an elongated rectangle with a size of 1 × 30 which was chosen to preserve vertical connectivity of the fringe patterns while suppressing horizontal connections. The vertical dimension was adjusted based on the curvature of the fringe contours observed in our experimental setup. Using this structuring element to perform morphological opening operation on the binary image, i.e., erosion followed by dilation, the operation process can be expressed as:

$B_{\text{clean}}=\left(B\ominus S_{\text{comb}}\right)\oplus S_{\text{comb}}$ (11)

After the above operations, the isolated white points in the image have been removed. To repair the edge jaggedness that may be caused by morphological operations, a small elliptical kernel is then selected for slight dilation processing, making the line edges exhibit natural smoothness.

2.3. Edge Extraction Strategy Based on Gradient Magnitude

As a classic edge detection method, the Canny algorithm performs excellently in many complex natural images, but it has obvious deficiencies when processing images with the following characteristics, such as square wave stripes: First, non-maximum suppression leads to edge breakage. The non-maximum suppression (NMS) of the Canny algorithm requires that edge pixels must be local maxima in the gradient direction at that point. However, the edges of square wave stripes are not ideal “peaks,” but rather have high gradient responses within several pixels’ width. NMS will forcibly retain only one pixel, resulting in edges that are too sparse or even broken. Second, the failure of double threshold hysteresis tracking. Canny’s double threshold strategy (high threshold + low threshold connection) assumes that edges are continuous. However, under light spot noise interference, stripe edges may be discontinuous, causing the hysteresis tracking algorithm to fail to connect broken edge segments. Third, sensitivity to noise. Even after Gaussian smoothing, Canny will still misidentify some high-gradient noise as edges.

In response to the limitations mentioned above, this paper proposes an edge extraction method based on direct thresholding of gradient magnitude. The core idea is: the edges of square wave stripes have extremely high gradient contrast, so only a global threshold needs to be set for the gradient magnitude to accurately segment edges, without the need for complex NMS and double threshold tracking.

The algorithm steps are as follows:

First, calculate the gradient magnitude of the image. Use the Sobel operator to calculate horizontal and vertical gradients respectively as Equation (12):

$G_x=I\ast S_x ,G_y=I\ast S_y$ (12)

where $S_x$ and $S_y$ are Sobel kernels, and * denotes convolution. Then calculate the gradient magnitude:

$G=\sqrt{G_x^2+G_y^2}$ (13)

The threshold selection is adaptive. By analyzing the distribution of gradient magnitudes, n% is chosen as the threshold:

$T_{grad}=\text{percentile}\left(G,n\right)$ (14)

In this simulation experiment, the 75% gradient threshold was selected for edge extraction. Assuming that edge pixels occupy a proportion p of the image, the optimal threshold corresponds to the (1 − p) ×100% percentile, which theoretically separates the edges from the background. In our experimental setup, edge pixels account for approximately 25% of each row, resulting in a 75% threshold. For other scenarios, the threshold should be adjusted according to the stripe density. Since the edge of each stripe typically spans 1 - 5 pixels per row, and subsequent 3D reconstruction requires selecting central points for row-wise depth calculation, the edge proportion can be estimated based on the reference plane geometry and stripe spacing. This estimated value subsequently determines the appropriate percentile threshold used for robust edge detection. While low-gradient pixels are mainly flat regions and noise. Finally, the edge image is obtained:

$E\left(x,y\right)=\{\begin{array}{l}255,\text{ if }G\left(x,y\right)>T_{grad}\\ 0,\text{ otherwise}\end{array}$ (15)

Compared to Canny, this method retains width information. For subsequent depth information calculation, the center pixel of the width can be selected to reduce errors. By eliminating NMS and hysteresis tracking, the time complexity is reduced to a certain extent. For the pipeline strategy provided in this paper, the overall complexity is O(T × M × N), where M and N are the length and width of the image.

3. Structural Light 3D Reconstruction Based on Square Wave Patterns

The main workflow of the 3D reconstruction algorithm employed in this study is divided into three stages: data collection, data processing, and error correction.

During data collection, the experimental setup is assembled and key parameters are measured. These parameters include the distance between the camera and the projector, the camera height, and the scaling factor between pixel coordinates and physical length units. Image frames without projected fringes are used for PIV analysis, while frames containing fringe deformation are used for 3D reconstruction.

The core of the data-processing stage is the extraction of depth information. This research method requires two images of overlaid fringes, namely the deformed fringe images of the reference plane and the target topography. The basic principle for depth acquisition is illustrated in Figure 2. One projected fringe is taken as an example to explain the method. Point C denotes the camera, point B denotes the projector, and point D represents a point on the inclined surface under measurement. Point A is defined as the intersection of the projector ray passing through point D with the reference plane. Point E is the intersection of the camera ray passing through point D with the reference plane. Here, H denotes the camera height, h denotes the height of the measured object, and d denotes the distance between the camera and the projector. The height of the measured point can be determined from a pair of similar triangles with known geometric parameters. The height h is given by:

$h=\frac{H\times L_{AE}}{d+L_{AE}}$ (16)

where $L_{AE}$ is derived from the displacement of the deformed fringe pattern and is expressed in pixel units.

Figure 2. Schematic diagram of 3D reconstruction algorithms: (a) Similar triangle method; (b) Causes of fringe errors.

The algorithm workflow is summarized as follows:

1) Use the edge detection process from Section 2 to obtain the edge information of the stripes.

2) During crater formation, particle splashing and the morphology of the impacting sphere can cause local blurring of the fringe boundaries during binarization. As crater formation is a gradual evolution process and the boundary-blurred region is relatively narrow, this region can be identified and compensated for through fitting function.

3) The pixel value of $L_{AE}$ is obtained by jointly processing the fringe pattern on the reference surface and the deformed fringe patterns in each frame. A ruler is placed within the camera’s field of view to determine the scaling factor ρ between pixel coordinates and physical lengths. If the camera position changes, this coefficient must be recalibrated. The actual length $L_{AE}$ is then calculated as $L_{AE}$ × ρ.

4) According to Equation (16), calculate the depth distance h from the reflection point D to the reference plane line by line.

The projection employed in this study consists of vertical square-wave fringes. By capturing images of multiple consecutive planes within the medium, the relevant error parameters can be determined. As illustrated in Figure 2(a), the extension of line CD coincides with the vertical line from the measured point D to the reference plane, forming segment DE. Consequently, the fringe $S_1$ containing point D is located directly beneath the camera, which yields the highest depth accuracy reconstruction based on the measured distance d.

However, the projector emits multiple parallel fringes. For fringes where the similar-triangle calculations are no longer rely on right triangles (i.e., non-orthogonal similar triangles), as intuitively illustrated in Figure 2(b), systematic discrepancies arise: on the left side of the high-speed camera, the experimentally measured $L_{A{E}^{\prime }}$ differs from $L_{AE}$ used in the similar-triangle calculation, with $L_{AE}$ being shorter than $L_{A{E}^{\prime }}$ ; On the right side of the high-speed camera, the experimentally measured $L_{I{J}^{\prime }}$ is likewise different from $L_{IJ}$ in the calculation, where $L_{IJ}$ is longer than $L_{I{J}^{\prime }}$ . Consequently, these deviations introduce reconstruction errors.

Since the fringe that provides the most accurate depth information based on the distance $d$ is the one located directly beneath the camera and perpendicular to the optical axis, the fringe $S_1$ on a plane with a known height and directly below the camera is selected as the reference. Moreover, because the projected vertical fringes are mutually parallel and perpendicular to the reference plane, this property can be exploited to determine the error parameters for each fringe and to perform error correction.

The error correction function is obtained through reconstruction using multiple planes with known heights. Each plane produces an inclined surface that contains error information. Using the known height, the corresponding height value on the inclined surface is selected as the reference line, from which the error parameters are calculated. Finally, by averaging the error functions derived from multiple planes, the overall accuracy can be further improved.

Due to the presence of concave and convex features on the measured surface, the projected fringes deform in opposite directions. The computed fringe displacements are signed values rather than absolute magnitudes, so the reconstructed depth values may be positive or negative depending on surface shape. Therefore, the error function is multiplied by the value to be corrected to ensure that the sign of the depth value is preserved after error correction, thereby avoiding incorrect morphological information in the final 3D reconstruction.

The depth correction procedure was performed using 10 reference planes positioned at known heights with a uniform spacing of 5 mm between adjacent planes. For each reference plane, the measured depth was compared against the ground truth to establish the correction function parameters, which were then averaged across all calibration measurements to improve robustness. Figure 3(a) shows a continuous plane with an actual measured height of 47 mm. The CCD camera is positioned near the coordinate origin. As the coordinate value increases, the position approaches the projector. Figure 3(b) presents the height information of the median row in Figure 3(a). The proposed form of the error correction function is given in (17). In this example, the parameters obtained for the median row are as follows: $a=9.671\times {10}^{-8}$ , $b=1.391\times {10}^{-4}$ and $c=0.8008$ .

$y=a{x}^2+bx+c$ (17)

4. Experiment Results

4.1. Image Preprocessing

To facilitate the evaluation of the rationality of the image processing algorithm pipeline in this paper, a computer is used to simulate and design square wave stripe images, and a large number of discrete light spots are added as noise.

Figure 3. Plane correction for 3D reconstruction: (a) 3D reconstruction of the 47 mm correction plane; (b) Height information of the median row in Figure 3(a); (c) Error correction function curve of Figure 3(b).

The algorithm model is evaluated by comparing the extracted stripe edges under noise-free and noisy conditions. This paper simulates three types of stripes with curvatures of 0, 0.001, and 0.0005 respectively. The discrete light spots are unipolar bright spot noise simulating optical scattering bright spots, with 500 points and 2000 light spots added respectively. And the size consists of random light spots with diameters ranging from 2 to 9 pixels. The image size is 400 × 400. Among them, Figure 4 shows the original image with 0.0005 curvature and the image after adding noise. The edge curves obtained directly using the Canny algorithm without processing are shown in Figure 5.

Figure 4. Fringe image with a curvature of 0.0005.

Figure 5. Unprocessed stripe edge image.

For the simulated noise image sets with three types of curvature, after processing through the algorithm pipeline in this paper, the edge processing results are shown in Figure 6.

Figure 6. Edge extraction results of high-noise images.

The evaluation results of IoU and F1 scores are shown in Table 2. The values in Table 2 are the averages of multiple measurements. Through the controlled variable method, comparisons are made with other denoising algorithms such as Gaussian filtering and edge processing methods such as the Canny algorithm, showing certain improvements in all cases. Therefore, this algorithm has certain advantages for image processing of discrete light spot clusters.

Table 2. Comparison of algorithm evaluation metrics.

Number of noise points	Method	Curvature	F1-score	IoU
500	Gaussian filtering	Vertical stripes	0.9732	0.9669
	Canny (50, 150)		0.9744	0.95
	Algorithm of this paper		0.9870	0.9744
2000	Gaussian filtering		0.9197	0.8513
	Canny (50, 150)		0.9713	0.9442
	Algorithm of this paper		0.9817	0.9641
500	Gaussian filtering	0.0005	0.8937	0.8078
	Canny (50, 150)		0.8557	0.7478
	Algorithm of this paper		0.9226	0.8564
2000	Gaussian filtering		0.7728	0.6298
	Canny (50, 150)		0.8021	0.6695
	Algorithm of this paper		0.8976	0.8142
500	Gaussian filtering	0.001	0.8382	0.7215
	Canny (50, 150)		0.7538	0.6049
	Algorithm of this paper		0.8757	0.7789
2000	Gaussian filtering		0.6985	0.5366
	Canny (50, 150)		0.6738	0.5080
	Algorithm of this paper		0.8219	0.7080

To better match actual conditions, the effect is set as Gaussian light spots, where the center has the highest grayscale value and the surrounding grayscale values gradually decrease. Under the condition of 0.001 curvature, 3000 noise light spots are set, and the effect is shown in Figure 7.

(a) Noisy original image (b) Edge processing of the original image

(c) Binarization processing result (d) Final result

Figure 7. Algorithm processing results.

The F1-Score is 0.8658 and the IoU index is 0.7634. When the curvature is 0.02, the F1-Score is 0.8014 and the IoU index is 0.6886.

For the images that need to be processed in the context of this paper, the original image and its edge processing results are shown in Figure 8. It can be intuitively seen that the application effect is good with high adaptability, which can provide a reliable data source and support for subsequent 3D reconstruction.

4.2. Experimental Verification of 3D Reconstruction

The particle experiment apparatus is described as follows: The imaging system consisted of a Revealer_5F04 camera and a LightCrafter 4500 projector. Images were captured at a resolution of 1720 × 236 pixels. Camera distortion parameters were obtained using Zhang’s calibration method. Notably, the camera and

(a) Original experimental image (b) Edge processing result

Figure 8. Application effect on actual images.

projector do not need to be mounted at the same height; through algebraic substitution in Equation (18), only the camera height needs to be measured. The scale factor was determined by capturing an additional image containing a calibrated ruler and computing the pixel-to-millimeter ratio from the ruler markings. Green electroplated glass beads with diameters ranging from 0.8 mm to 1.0 mm were used as the granular material, chosen because green exhibits high brightness in our experimental lighting conditions, requiring only high-contrast colored fringes for effective imaging. The projected fringe width was set to four times the particle diameter; narrower fringes tend to cause stripe adhesion and merging, while wider fringes result in loss of depth information during rapid impact motion.

Two verification experiments were designed for particle packing: planar height measurement of the granular bed and slope measurement of a flat inclined surface formed by the packed particles.

In the planar verification experiment, the actual height of the measured granular bed region was approximately 9.5 mm. The 3D reconstruction result is shown in Figure 9(a). As shown in Figure 9(b), the height ranged from 9.1 mm to 9.9 mm, with an error within one particle diameter (0.9 - 1.0 mm).

Figure 9. (a): Particle planar 3D reconstruction results; (b): The numerical range of height in this area.

In the inclined slope verification experiment, an inclined particle surface with a height of 32 mm and a width of 288 mm was constructed. The reference slope was $k_0=-0.1111$ . The reconstruction result is shown in Figure 10(a). By selecting an arbitrary row (for example, row 820), the function fitting result is given in (18), where $k_0=-0.1111$ and the absolute error is 0.0015.

$h=-0.1096x+7.415$ (18)

This paper proposes a simple and effective 3D reconstruction method that uses only square-wave fringes to reconstruct the surface of a granular medium. Previous approaches could only extract information from a single fringe to obtain 3D depth information for a cross-section [16], whereas the method proposed in this paper projects multiple fringes onto the surface of the granular medium and simultaneously extracts deformation data from all fringes, thereby capturing the 3D depth information at a given moment. The extracted data are then reconstructed, fitted, and optimized. In experiments involving multiple impacts that cause deformation of a granular bed, the proposed method can extract depth information at arbitrary spatial locations and at different moments in time, enabling subsequent experimental comparison and analysis. In static verification experiments, the slope error of the inclined surface was 0.0015, and the planar error was within one particle diameter (0.9 - 1.0 mm).

Figure 10. 3D reconstruction results of the inclined plane; (b): Comparison of the reconstruction results of row 820 in (a) with its fitting function.

This study validated the effectiveness of the proposed method under controlled static conditions. However, in practical dynamic granular flow scenarios, rapid motion may introduce motion blur, severe occlusion can compromise fringe extraction accuracy, and transparent or highly reflective particles may produce complex optical interference. The performance under these conditions requires further validation. Future work will incorporate high-speed cameras and adaptive exposure strategies to address dynamic scenarios.

5. Conclusion

In previous research in the field of particle kinematics, the strong reflective property of particles has been utilized for velocity analysis of kinematic fluids (PIV). However, this characteristic poses a challenge for classical algorithms in structured light 3D reconstruction. Moreover, since fluid analysis in the field of particle kinematics is holistic-based, the use of square wave stripes as active features is proposed to extract edge information and obtain 3D depth information through the similar triangle method. Among these, the image processing algorithm is a core component of this workflow. This paper, targeting the characteristics of particle images and within the acceptable error range in the field of particle science, proposes an advanced approach: first performing anisotropic diffusion denoising, then processing through Otsu and morphological operations to obtain smooth binarized image curves, and finally acquiring edge information through gradient magnitude. In simulated image experiments, this algorithm demonstrates advantages over other algorithms. Finally, through validation experiments with known particle beds, the operability and adaptability under actual conditions are confirmed, providing a 3D reconstruction analysis perspective and an effective measurement scheme for future research on more complex particle motion scenarios.

Funding

This work has been supported by the National Natural Science Foundation of China (12072200, 12372384), and the National Key Research and Development Program of China under Grant 2022YFF0607701.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Li, X., Li, B., Qie, L., Li, Y. and Wang, H. (2025) Dynamic Response of a Frame Structure Impacted by Debris Flow Containing Large Boulders. Structures, 74, Article 108543.[CrossRef]
[2]	Chen, J., Zhang, W., Cao, C., Yin, H., Wang, J., Li, W., et al. (2024) The Effect of the Check Dam on the Sediment Transport and Control in Debris Flow Events. Engineering Geology, 329, Article 107397.[CrossRef]
[3]	Cundall, P.A. (1971) A Computer Model for Simulating Progressive Large-Scale Movements in Blocky Rock Systems. Proceedings of the Symposium of the International Society for Rock Mechanics, Society for Rock Mechanics (ISRM), II-8.
[4]	Adrian, R.J. (1989) Scattering Particle Characteristics and Their Effect on Pulsed Laser Measurements of Fluid Flow: Speckle Velocimetry vs Particle Image Velocimetry. Applied Optics, 23, 1690-1691.[CrossRef] [PubMed]
[5]	Qu, Y., Huang, J. and Zhang, X. (2018) Rapid 3D Reconstruction for Image Sequence Acquired from UAV Camera. Sensors, 18, 225-244.[CrossRef] [PubMed]
[6]	Hu, T., Song, C., Zhuang, J. and Lyu, Y. (2025) Automatic Segmentation and Measurement System of 3D Point Cloud Images Based on RGB-D Camera for Rat Wounds. Biomedical Signal Processing and Control, 104, Article 107682.[CrossRef]
[7]	Lin, J.F. and Su, X. (1995) Two-Dimensional Fourier Transform Profilometry for the Automatic Measurement of Three-Dimensional Object Shapes. Optical Engineering, 34, Article 3297.[CrossRef]
[8]	He, Q., Ning, J., Liu, X. and Li, Q. (2025) Phase-Shifting Profilometry for 3D Shape Measurement of Moving Objects on Production Lines. Precision Engineering, 92, 30-38.[CrossRef]
[9]	Ciulla, C. (2025) Two-Dimensional Image Noise Removal and Reconstruction Using Discrete Fourier Transform, K-Space Filtering and Z-Space Filtering. Progress in Engineering Science, 2, Article 100056.[CrossRef]
[10]	Kraus, K. and Pfeifer, N. (1998) Determination of Terrain Models in Wooded Areas with Airborne Laser Scanner Data. ISPRS Journal of Photogrammetry and Remote Sensing, 53, 193-203.[CrossRef]
[11]	Jiang, Y., Wu, Z., Liu, Y., Yu, X. and Zhang, Q. (2023) High-Speed 3D Shape Measurement Using Efficient Moiré-Assisted Three-Frequency Heterodyne Phase Unwrapping Algorithm. Optics and Lasers in Engineering, 161, Article 107383.[CrossRef]
[12]	Sarbolandi, H., Lefloch, D. and Kolb, A. (2015) Kinect Range Sensing: Structured-Light versus Time-of-Flight Kinect. Computer Vision and Image Understanding, 139, 1-20.[CrossRef]
[13]	Scivetti, M., Pilolli, G.P., Corsalini, M., Lucchese, A. and Favia, G. (2007) Confocal Laser Scanning Microscopy of Human Cementocytes: Analysis of Three-Dimensional Image Reconstruction. Annals of Anatomy-Anatomischer Anzeiger, 189, 169-174.[CrossRef] [PubMed]
[14]	Zhang, S. (2018) High-Speed 3D Shape Measurement with Structured Light Methods: A Review. Optics and Lasers in Engineering, 106, 119-131.[CrossRef]
[15]	Li, D., Xu, L., Tang, X., Sun, S., Cai, X. and Zhang, P. (2017) 3D Imaging of Greenhouse Plants with an Inexpensive Binocular Stereo Vision System. Remote Sensing, 9, Article No. 508.[CrossRef]
[16]	Zheng, Y., Sun, W., Yang, H. and Li, R. (2024) Three-Dimensional Reconstruction of Impact Crater Morphology during the Impact Process. IEEE Access, 12, Article 163723.[CrossRef]
[17]	Zhang, Q., Huang, C., Yang, L. and Yang, Z. (2023) Salt and Pepper Noise Removal Method Based on Graph Signal Reconstruction. Digital Signal Processing, 135, Article 103941.[CrossRef]
[18]	Yesli, S. and Diaf, M. (2026) Three Unbiased Anisotropic Diffusion Filtering Models with Optimized Parameters for Rician Noise Removing in MR Images. Magnetic Resonance Imaging, 130, Article 110673.[CrossRef]
[19]	Emmanuel, A.E., Amusa, K.A., Erinosho, T.C. and Raji, M.T. (2026) Optimised Canny Edge Detection Algorithm for Medical Image Feature Mapping and Extraction. Engineering Applications of Artificial Intelligence, 167, Article 113763.[CrossRef]
[20]	Arino, A., Yamamoto, N., Sasaki, K., Kimura, R. and Aizawa, T. (2026) Quantitative Evaluation of Fatty Change in Rotator Cuff Muscles Using the Otsu Thresholding Method. Journal of Shoulder and Elbow Surgery, 35, 200-208.[CrossRef] [PubMed]