Multi-particle 3D shape reconstruction method based on projection consistency

Abstract

A multi-particle 3D shape reconstruction method based on projection consistency is proposed. This method uses a three-angle laser interferometric particle imaging system to acquire defocused images from three viewpoints. The defocused speckle image is obtained by cropping the central M×N pixel range, and a phase retrieval algorithm is used for reconstruction to obtain the 2D particle shape. A grayscale threshold of 70 or higher is set as the target contour region, which is then binarized and smoothed. Four candidate results are projected into 3D space using Matlab, generating 3D structures in the intersecting regions. These structures are then reprojected back to each viewpoint. The error S between the reprojected image and the original input projection is calculated, and the image with the smallest S is selected as the final reconstruction result. This method effectively eliminates trivial blur interference in the phase retrieval process through a systematic iterative optimization and error evaluation mechanism, achieving high-precision reconstruction of multi-particle 3D shapes.

Description

Technical Field

[0001] This invention belongs to the field of image processing and three-dimensional reconstruction in optical measurement, and specifically relates to a multi-particle three-dimensional shape reconstruction method based on projection consistency. Background Technology

[0002] Irregularly shaped particles are widely present in the natural environment, daily life, and engineering practice. Measuring the size, shape, velocity, concentration, and three-dimensional position of these particles in a particle distribution field is of significant scientific importance for studying environmental climate change and assessing engineering production safety. Based on Mie scattering theory, interferometric particle imaging technology obtains information related to particle parameters by analyzing and processing interferometric images of irregular particles in the frequency domain. In fields such as meteorology, combustion, nuclear industry, and cosmetic testing, a large number of irregular particles exist in production activities. The shape parameters of these particles affect their scattered light distribution; therefore, the study of particle shape is of great and far-reaching significance in the measurement of irregular particles.

[0003] Currently, there are many methods for measuring irregular particles based on laser interferometry theory. However, the Fourier phase retrieval methods used to reconstruct the two-dimensional particle projection contours all suffer from trivial ambiguity, which may lead to significant errors in three-dimensional reconstruction. This problem arises because the Fourier phase retrieval problem has multiple solutions. Different signals have the same Fourier transform intensity, meaning that the solution for recovering phase information from the Fourier amplitude is not unique. The following three transforms and their linear combinations have the same Fourier transform amplitude: global phase shift, spatial shift, and conjugate transpose. This type of solution does not change the overall contour of the phase retrieval solution; it only shifts or rotates the original image by 180°. Global phase shift does not change the intensity of the original image; the resulting ambiguity is called trivial ambiguity.

[0004] Patent CN116698708A proposes a method for reconstructing the shape of irregular particles based on a phase retrieval algorithm. This method utilizes a laser interferometric particle imaging system to acquire defocused speckle images of irregular particles, and then reconstructs the two-dimensional shape of the measured particle from the defocused speckle image using a phase retrieval algorithm. This method features a simple acquisition system and good particle shape reconstruction results, but it is limited to two-dimensional shape measurement, and its reconstruction performance may be poor when measuring multiple particles.

[0005] Patent CN117890285A proposes a method for three-dimensional measurement of irregular particles based on dual-angle defocused interferometric images. This method acquires two-dimensional shape information of particles from defocused interferometric images and then reconstructs the three-dimensional contours of particles by combining the shape information from two angles. This method has a simple acquisition system and achieves good three-dimensional shape reconstruction results, but it does not consider the trivial blurring problem in two-dimensional reconstruction, and the reconstruction target is only a single particle.

[0006] The aforementioned patents all measure the shape of irregular particles using a laser interferometric imaging system, but they do not consider the trivial fuzziness problem that may exist in the two-dimensional phase retrieval algorithm. Furthermore, the reconstruction target is only a single particle. In three-dimensional reconstruction, if the input two-dimensional projection image may have a 180° rotation error due to phase fuzziness, the existing methods cannot reliably eliminate these interferences and reconstruct the correct three-dimensional model. Summary of the Invention

[0007] The purpose of this invention is to address the problems of existing reconstruction methods not considering trivial fuzziness in two-dimensional phase retrieval algorithms and the fact that the reconstruction target is only a single particle. This invention provides a multi-particle 3D shape reconstruction method based on projection consistency. This method can directly achieve accurate 3D shape reconstruction of multi-particle systems based on defocus speckle maps acquired by laser interferometric imaging systems. By effectively eliminating trivial fuzziness interference in the phase retrieval process, it provides technical support for high-precision 3D topography measurement of complex multi-particle fields.

[0008] This invention uses a hybrid input-output (HIO) phase retrieval algorithm with wavelet sparse constraints to obtain the two-dimensional shape information of particles from defocused interferometric images. It eliminates incorrect three-dimensional reconstruction caused by phase ambiguity during phase retrieval in two-dimensional reconstruction by using a projection consistency method. Furthermore, it reconstructs the three-dimensional contour of particles from multiple angles based on the shape information from a back-projection algorithm. This method is applicable to the three-dimensional shape reconstruction of single particles and multiple particles of various irregular particles with sub-millimeter size.

[0009] Technical solution of the present invention

[0010] The steps of the multi-particle 3D shape reconstruction method based on projection consistency proposed in this invention are as follows:

[0011] Step 1. Data Acquisition and Preprocessing: Using a three-angle laser interferometry particle imaging system, at the same defocus distance at the imaging end, the scattering defocus images I1, I2 and I3 of the multi-particle under test are acquired at the three scattering angles α, β and γ. The images within the M×N pixel range of the center of the defocus images I1, I2 and I3 are cropped to obtain the defocus speckle map O, which is saved as O1, O2 and O3.

[0012] Furthermore, the three scattering angles α, β, and γ are mutually perpendicular.

[0013] Furthermore, the shear speckle region is preferably square, with region sizes M and N being equal and both set as N1, and N1 being ≥256.

[0014] Step 2. Field Initialization: Initialize the complex scattering field distribution to be recovered, make an initial guess about the shape of the multi-particle sample, and set it as follows:

[0015]

[0016] It is a random amplitude distribution. It is a random phase distribution, where i is the imaginary unit.

[0017] Step 3. Iterative phase recovery: For the defocus speckle map O obtained at each scattering angle, the complex scattering field is iteratively updated by combining forward and backward propagation with amplitude constraints (such as wavelet sparsity constraints) and hybrid input-output (HIO) algorithms until convergence, and the two-dimensional multi-particle shape reconstruction results L1, L2, and L3 at each scattering angle are obtained.

[0018] Furthermore, the specific steps of iterative phase recovery include:

[0019] Step 3-1. Forward propagation; calculate the frequency domain estimate at the k-th iteration. .

[0020]

[0021] It is the Fast Fourier Transform. This is the surface estimate of the object in the k-th iteration, where k = 1, 2, 3, ... , To maximize the number of iterations, the initial guess set in step 2 is used in the first iteration. .

[0022] Step 3-2. Amplitude replacement; replace the frequency domain estimate with the square root of the defocus speckle map O from Step 1. The strength, and retain The phase is used to obtain the frequency domain estimate after amplitude constraint. The defocus speckle patterns O are O1, O2, and O3, respectively.

[0023]

[0024] Step 3-3. Backpropagation; Perform inverse Fast Fourier transform on the frequency domain estimate after amplitude constraint to calculate the intermediate surface estimate. ,

[0025]

[0026] Steps 3-4. Wavelet sparse update; estimate the intermediate surface according to the following formula. Apply wavelet sparsity and update.

[0027]

[0028] in, and These are wavelet transform and inverse wavelet transform, respectively. It is a hard threshold function:

[0029]

[0030] Where c represents the wavelet coefficient variable, and the threshold of the hard threshold function. Selected by the following method, The coefficient matrices of the multiple sub-bands decomposed by wavelet transform are flattened into a one-dimensional vector, and the total energy is calculated.

[0031]

[0032] in, This represents the total number of all wavelet coefficients. Next, all coefficients are sorted in descending order of their absolute values ​​to obtain the sequence. And calculate its cumulative energy:

[0033]

[0034] Find the smallest index position This allows the accumulated energy to meet the requirements. And set the magnitude of the position coefficient as the threshold τ.

[0035] Furthermore, the wavelet basis of the wavelet transform is a family of functions that satisfy orthogonality or bioorthogonality, are compactly supported or non-compactly supported, and have vanishing moments of orders 2-10. The introduction of wavelet sparsity constraints, especially under low sampling rate conditions (e.g., when the focal distance is small or the object size is large), can significantly improve the quality of the two-dimensional reconstruction results.

[0036] Steps 3-5: Shrink the package support domain update; estimate after wavelet sparse step. Smoothing is performed using a Gaussian kernel function with a standard deviation of 1-2, as shown in the following formula:

[0037]

[0038] Operator It's a convolution operator, the standard deviation σ has a value of 1-2, then... The threshold t is used to binarize the domain into the support domain S(x,y).

[0039]

[0040] The threshold t is calculated by the formula t = 0.2max(A(x, y)).

[0041] Steps 3-6. Update using the Hybrid Input / Output Phase Recovery (HIO) algorithm; then smooth the operation... Update using the following HIO formula, the updated result is... As input for the next iteration,

[0042]

[0043] Among them, feedback parameters Set it between 0.75 and 0.95.

[0044] Step 3-7. Repeat steps 3-1 to 3-6 until the maximum number of iterations is reached. Output the two-dimensional multi-particle shape reconstruction results L1, L2, and L3 of the multi-particle under test in the three directions of α, β, and γ.

[0045] Furthermore, the maximum number of iterations Set to 500-2000 times.

[0046] Step 4. Morphological processing: Morphological processing is performed on the two-dimensional multi-particle shape reconstruction results L1, L2, and L3 obtained by iterative phase recovery to make the edges smoother.

[0047] Furthermore, the morphological processing is a closed operation, and the structuring element size is 3×3 to 9×9.

[0048] Step 5. 3D Reconstruction and Optimization: Based on the 2D shape reconstruction results of each scattering angle after Step 4, candidate 3D reconstruction models are generated, evaluated by calculating projection error, and the optimal model is selected as the final 3D reconstruction result H.

[0049] Furthermore, the 3D reconstruction and optimization steps are as follows:

[0050] Step 5-1. Generate candidate 3D reconstruction models; In 3D reconstruction, consider that each 2D projection has two rotation states of 0° or 180°. The three viewpoints combine to form a total of 8 candidate arrangement schemes. Due to the symmetry of global rotation, the actual independent candidate combinations are 4 of them. Therefore, back-projection 3D reconstruction is performed on the 4 candidate combinations to obtain 4 candidate 3D models.

[0051] Step 5-2. Calculate the reprojection error; reproject each candidate 3D model along the given three observation directions α, β, γ to generate a reprojected image K(i, j); calculate the error S between the reprojected image K(i, j) and the original projection J(i, j) using the following formula:

[0052]

[0053] The original projection J(i, j) refers to the projection map of the set of 2D reconstructed images used to generate the current candidate 3D reconstruction model in each viewpoint direction (α, β, γ). The cubic function applies a non-linear penalty to the matching error. This design can significantly amplify the contribution of larger errors, making it easier to distinguish the comprehensive error differences between different objects, thereby optimizing the selection or evaluation of the target structure. , , Let K(i, j) be the reprojection image of the candidate 3D model under the acquisition viewpoints in the three directions of α, β, and γ, and let J(i, j) be the mean square error between the original projection and the reprojection image K(i, j). The formula for calculating J(i, j) is as follows:

[0054]

[0055] Step 5-3. Select the optimal solution; since a correct 3D reconstruction should present consistent structural features from different viewpoints, its reprojection error should be minimized. Therefore, the candidate model with the smallest reprojection error S is selected as the final 3D reconstruction result H.

[0056] Advantages and beneficial effects of the present invention:

[0057] This invention proposes a 3D reconstruction method to eliminate reconstruction ambiguities caused by trivial fuzziness in multi-particle 2D phase retrieval. Based on the principle of projection consistency, this method effectively eliminates the interference of trivial fuzziness during phase retrieval on 3D reconstruction. Its correction algorithm is computationally efficient and has a simple process, significantly reducing computational complexity while ensuring reconstruction accuracy. This method can achieve accurate reconstruction of the 3D shape of multiple particles using the defocus speckle map obtained by a laser interferometric imaging system, without the need for additional illumination paths, resulting in a simple and reliable system structure. This method can reconstruct the 3D morphology of each particle in a multi-particle system relatively well, providing an effective technical means for high-precision 3D measurement of complex multi-particle fields. Attached Figure Description

[0058] Figure 1 This is a flowchart illustrating the reconstruction process of the multi-particle 3D shape reconstruction method based on projection consistency in this embodiment of the invention.

[0059] Figure 2 This is an illustration of the trivial fuzziness present in the phase retrieval algorithm of this invention, wherein: (a) is a simulated interference speckle pattern of a PBY-shaped simulated particle, and (b) and (c) are possible reconstruction results of the speckle.

[0060] Figure 3These are the effect diagrams of the two-dimensional particle shape reconstruction method based on defocus speckle sampling at different sampling rates in this invention. Among them, (a) is the reconstruction result of the HIO algorithm with wavelet sparsity constraints, and (b) is the reconstruction result of the HIO algorithm without wavelet sparsity constraints. It can be seen that the former has better reconstruction performance in the case of larger objects.

[0061] Figure 4 This is the particle three-dimensional reconstruction process in this invention. (a), (b), (c), and (d) represent all possible reconstruction results, respectively.

[0062] Figure 5 These are the effect diagrams of multi-particle 3D reconstruction in this invention. (a) is the simulated multi-particle 3D model, and (b) is the 3D model reconstructed by the method described above.

[0063] Figure 6 This is a schematic diagram of an embodiment of the three-angle laser interferometric particle imaging system of the present invention, wherein: 1 is a laser, 2 is a beam expanding and collimating system, 3, 5, and 7 are the first, second, and third imaging lenses, respectively, 4, 6, and 8 are the first, second, and third CCD cameras, respectively, 9 is the particle being measured, α, β, and γ are the scattering angles of the three acquisition angles, and the three viewing angles are perpendicular to each other. Detailed Implementation

[0064] The flowchart of the multi-particle 3D shape reconstruction method based on projection consistency in this embodiment of the invention is as follows: Figure 1 As shown, the specific steps are as follows:

[0065] Step 1: Data Acquisition and Preprocessing

[0066] (1). Construct a three-angle laser interferometric particle imaging system (see Figure 6 At the imaging end, defocused images are acquired at three scattering angles: α, β, and γ. The acquired defocused images contain shape information of the particles at these three angles. The three scattering angles α, β, and γ are mutually perpendicular.

[0067] (2). At the same defocus distance, the defocus images at the above three scattering angles are acquired and recorded as I1, I2 and I3.

[0068] (3). The images within the center N1×N1 pixel range of the defocused images I1, I2 and I3 are cropped respectively to obtain the defocused speckle map O, and saved as O1, O2 and O3 respectively.

[0069] Furthermore, the size N1 of the shear speckle region is set to ≥256, and in this embodiment it is set to 512.

[0070] Step 2, Field Initialization:

[0071] (4). Make an initial guess about the particle shape and set it as follows:

[0072]

[0073] It is a random amplitude distribution. It is a random phase distribution, where i is the imaginary unit.

[0074] Step 3, Iterative Phase Recovery:

[0075] (5). Forward propagation. Calculate the frequency domain estimate at the k-th iteration. ,

[0076]

[0077] It is the Fast Fourier Transform. This is the surface estimate of the object in the k-th iteration, where k = 1, 2, 3, ... , To maximize the number of iterations, the initial guess set in step 2 is used in the first iteration. .

[0078] (6) Amplitude replacement. Replace the frequency domain estimate with the square root of the defocus speckle map O from step 1. The strength, and retain The phase is used to obtain the frequency domain estimate after amplitude constraint. The defocus speckle patterns O are O1, O2, and O3, respectively.

[0079]

[0080] (7) Backpropagation. Perform an inverse fast Fourier transform on the frequency domain estimate after amplitude constraint to calculate the intermediate surface estimate. ,

[0081]

[0082] (8) Wavelet sparse update. The intermediate surface is estimated according to the following formula. Apply wavelet sparsity and update.

[0083]

[0084] in, and These are wavelet transform and inverse wavelet transform, respectively. It is a hard threshold function:

[0085]

[0086] c represents the wavelet coefficient variable, and the threshold of the hard thresholding function. Selected by the following method, The coefficient matrices of the multiple sub-bands decomposed by wavelet transform are flattened into a one-dimensional vector, and the total energy is calculated.

[0087]

[0088] in, This represents the total number of all wavelet coefficients. Next, all coefficients are sorted in descending order of their absolute values ​​to obtain the sequence. And calculate its cumulative energy:

[0089]

[0090] Find the smallest index position This allows the accumulated energy to meet the requirements. And set the magnitude of the position coefficient as the threshold τ.

[0091] Furthermore, the wavelet basis of the wavelet transform is a family of functions that satisfy orthogonality or bioorthogonality, are compactly supported or non-compactly supported, and have a certain number (2-10th order) of vanishing moments. The introduction of wavelet sparsity constraints, especially under low sampling rate conditions (e.g., when the focal distance is small or the object size is large), can significantly improve the quality of the two-dimensional reconstruction results.

[0092] (9) Shrinking the wrapper support domain update. Multi-particle shape estimation after wavelet sparse step. Perform a Gaussian kernel smoothing operation with a standard deviation of 1-2. The formula is as follows:

[0093]

[0094] Operator It is a convolution operator, standard deviation -2, then The threshold t is used to binarize the domain into the support region S(x, y).

[0095]

[0096] The threshold t is calculated by the formula t = 0.2max(A(x, y)).

[0097] (10) Update using the Hybrid Input / Output Phase Recovery (HIO) algorithm. After smoothing... Update using the following HIO formula to obtain the input for the next iteration. ,

[0098]

[0099] Among them, feedback parameters The value is set between 0.75 and 0.95; in this embodiment, it is set to 0.9.

[0100] (11). Repeat steps (5)-(10) in step 3 until the maximum number of iterations is reached. Output the reconstruction results L1, L2, and L3 of the two-dimensional multi-particle shape of the multi-particle under test in the three directions of α, β, and γ.

[0101] Furthermore, the maximum number of iterations The frequency is set to 500-2000 times; in this embodiment, it is set to 800 times.

[0102] Step 4, Morphological processing:

[0103] (12). Based on the reconstructed shapes L1, L2, and L3 of the multi-particles in the three directions of α, β, and γ obtained in step 3, morphological processing is performed to make the edges smoother.

[0104] Furthermore, the morphological processing is a closing operation, and the size of the structuring element is 3×3 to 9×9, which is 5×5 in this embodiment.

[0105] Step 5, 3D Reconstruction and Optimization:

[0106] (13) Generate candidate 3D reconstruction models; In 3D reconstruction, considering that each 2D projection has two rotation states of 0° or 180°, the three-view combinations constitute a total of 8 candidate arrangement schemes. Due to the symmetry of global rotation, the actual independent candidate combinations are 4 of them. Therefore, back-projection 3D reconstruction is performed on the 4 candidate combinations to obtain 4 candidate 3D models.

[0107] (14). Calculate the reprojection error. Reproject each candidate 3D model along the three given observation directions α, β, γ to generate a reprojected image K(i, j); calculate the error S between the reprojected image K(i, j) and the original projection J(i, j) using the following formula:

[0108]

[0109] The original projection J(i, j) is the projection map of the set of two-dimensional reconstructed images used to generate the current candidate three-dimensional reconstruction model in the α, β, and γ viewpoint directions. , , Let K(i, j) represent the mean square error between the reprojected image K(i, j) of the candidate 3D model and the original projected image J(i, j) in the three directions α, β, and γ, respectively. The calculation formula is as follows:

[0110]

[0111] (15) Select the optimal solution and choose the candidate model with the smallest reprojection error S as the final three-dimensional reconstruction result H.

[0112] Example 1: Three-Angle Laser Interferometric Particle Imaging System

[0113] The three-angle laser interferometry particle imaging system used in this invention is as follows: Figure 6 As shown, 1 is the laser, 2 is the beam expander and collimator system, 3, 5, and 7 are the first, second, and third imaging lenses, respectively, 4, 6, and 8 are the first, second, and third CCD cameras, respectively, and 9 is the particle being measured.

[0114] Laser 1 is a semiconductor laser with a wavelength of 532nm. The emitted laser beam is collimated by a microscope objective with a magnification of 10×, a pinhole filter with a aperture width of 10μm, and a collimating lens. The scattering angles α and β are both 45°, and γ is 90°. The three acquisition angles are perpendicular to each other. The imaging lenses that collect the scattered light from the particles are Nikon 50mm f / 1.4D fixed-focus lenses 3, 5, and 7. The detectors 4, 6, and 8 are CCDs with a pixel size of 5.86μm, an effective pixel count of 1080×1920, and a frame rate of 165fps. The system magnification is 2.27, the object distance is 70mm, the image distance is 208mm, and the defocus distance is 14mm. The frontal dimensions of the irregular particle to be measured are 900μm×1028μm, and the side thickness is 83μm.

[0115] Example 2: Multi-particle 3D shape reconstruction

[0116] A 512×512 pixel region was cropped from the defocus speckle images obtained from three viewpoints. Phase retrieval was used for reconstruction to obtain the two-dimensional particle shapes. A grayscale threshold of 70 or higher was set as the target contour region, which was then binarized and segmented. The target contour region was smoothed. The four candidate results were projected into 3D space in Matlab, generating 3D structures in the intersecting regions. These structures were then reprojected onto each viewpoint, and the error S between the reprojected image and the original input projection was calculated. The image with the smallest S was selected as the final reconstruction result.

[0117] First, the three-angle laser interferometry particle imaging system built in Example 1 ( Figure 6 In the three mutually perpendicular directions α, β, and γ, defocused images are acquired respectively, and the central 512×512 pixel area is cropped to obtain speckle data O1-O3 (step 1 in the specific implementation).

[0118] In the phase retrieval stage (step 3), a hybrid iterative algorithm is adopted: forward propagation is performed through Fourier transform, and the shape estimation is updated by inverse transform after replacing the amplitude with measured data; in particular, the introduction of db4 wavelet sparsity constraints (step 3, (8)) significantly improves the reconstruction quality under low sampling rates, such as Figure 3 As shown, the fuzzy results without wavelet sparsity are compared ( Figure 3 (a) After adding the outline, the reconstructed contour is significantly clearer. Figure 3 (b)). To resolve the 180° rotation ambiguity in phase recovery ( Figure 2 ), and perform three-dimensional backprojection on the four projection combinations respectively (step (13) in step 5), and use the reprojection error Selecting the optimal model: such as Figure 4 In the two-particle case shown, the one with the smallest error is... Figure 4 (b) was selected as the final reconstruction result. Figure 5 The effectiveness of this method on various two-particle samples is further demonstrated. Through a systematic iterative optimization and error evaluation mechanism, this method achieves high-precision reconstruction of multi-particle 3D shapes.

Claims

1. A multi-particle 3D shape reconstruction method based on projection consistency, comprising: Step 1. Data Acquisition and Preprocessing: Using a three-angle laser interferometry particle imaging system, at the same defocus distance at the imaging end, scattering defocus images of the multi-particle under test are acquired at three scattering angles α, β, and γ. The image within the M × N pixel range of the center of the defocus image is cropped to obtain a defocus speckle map O, which is saved as O1, O2, and O3. The three scattering angles α, β, and γ are mutually perpendicular. Step 2. Field Initialization: Initialize the complex scattering field distribution to be recovered, make an initial guess about the shape of the multi-particle to be measured, and set it as... ; Step 3. Iterative phase recovery: For the defocused speckle map O obtained at each scattering angle, the complex scattering field is iteratively updated by combining forward and backward propagation with amplitude constraints and a hybrid input-output phase recovery algorithm until convergence, and the two-dimensional multi-particle shape reconstruction results L1, L2, and L3 at each scattering angle are obtained. Step 4. Morphological processing: Morphological processing is performed on the two-dimensional multi-particle shape reconstruction results L1, L2, and L3 obtained by iterative phase recovery to make the edges smoother; Step 5. 3D Reconstruction and Optimization: Based on the 2D shape reconstruction results of each scattering angle after Step 4, candidate 3D reconstruction models are generated, evaluated by calculating projection error, and the optimal model is selected as the final 3D reconstruction result H.

2. The method according to claim 1, characterized in that, The shear speckle region in step 1 is a square region with equal M and N sizes, both set to N1, and N1 ≥ 256.

3. The method according to claim 1, characterized in that, The initial guess of the shape of the multi-particle to be tested in step 2 is: It is a random amplitude distribution. It is a random phase distribution, where i is the imaginary unit.

4. The method according to claim 1, characterized in that, The iterative phase recovery described in step 3 includes: Step 3-1. Forward propagation; calculate the frequency domain estimate at the k-th iteration. , It is the Fast Fourier Transform, f k (x, y) is the surface estimate of the object in the k-th iteration, where k = 1, 2, 3, ... , To maximize the number of iterations, the initial guess set in step 2 is used in the first iteration. ; Step 3-2. Amplitude replacement; replace the frequency domain estimate with the square root of the defocus speckle map O from Step 1. The strength, and retain The phase is used to obtain the frequency domain estimate after amplitude constraint. The defocus speckle patterns O are O1, O2, and O3, respectively. ； Step 3-3. Backpropagation; Perform inverse Fast Fourier transform on the frequency domain estimate after amplitude constraint to calculate the intermediate surface estimate. , Steps 3-4. Wavelet sparse update; estimate the intermediate surface according to the following formula. Apply wavelet sparsity and update. in, and These are wavelet transform and inverse wavelet transform, respectively. It is a hard threshold function; Steps 3-5: Shrink the package support domain update; estimate after wavelet sparse step. Smoothing is performed using a Gaussian kernel function with a standard deviation of 1-2, as shown in the following formula: Operator It is a convolution operator, with a standard deviation σ ranging from 1 to 2, then... The threshold t is used to binarize the domain into the support region S(x, y). The threshold t is calculated by the formula t = 0.2max(A(x, y)); Steps 3-6. Update using the Hybrid Input / Output Phase Recovery (HIO) algorithm; then smooth the operation... Update using the following HIO formula, the updated result is... As input for the next iteration, in, For feedback parameters; Step 3-7. Repeat steps 3-1 to 3-6 until the maximum number of iterations is reached. Output the two-dimensional multi-particle shape reconstruction results L1, L2, and L3 of the multi-particle under test in the three directions of α, β, and γ.

5. The method according to claim 4, characterized in that, The hard threshold function described in steps 3-4 is as follows: Where c represents the wavelet coefficient variable, and the threshold of the hard threshold function. Selected by the following method, The coefficient matrices of the multiple sub-bands decomposed by wavelet transform are flattened into a one-dimensional vector, and the total energy is calculated. in, This is the total number of all wavelet coefficients; next, all coefficients are sorted in descending order of their absolute values ​​to obtain the sequence. And calculate its cumulative energy: Find the smallest index position This allows the accumulated energy to meet the requirements. And set the magnitude of the position coefficient as the threshold τ.

6. The method according to claim 4, characterized in that, The wavelet basis of the wavelet transform is a family of functions that satisfy orthogonality or biorthogonality, are compactly supported or non-compactly supported, and have vanishing moments of order 2 to 10.

7. The method according to claim 4, characterized in that, The feedback parameters Set to between 0.75 and 0.95; maximum number of iterations. Set to 500~2000 times.

8. The method according to claim 1, characterized in that, The morphological processing described in step 4 is a closing operation, and the structuring element size is 3×3 to 9×9.

9. The method according to claim 1 or 2, characterized in that, The three-dimensional reconstruction and optimization steps described in step 5 are as follows: Step 5-1. Generate candidate 3D reconstruction models; In 3D reconstruction, consider that each 2D projection has two rotation states of 0° or 180°. The three viewpoints combine to form a total of 8 candidate arrangement schemes. Due to the symmetry of global rotation, the actual independent candidate combinations are 4 of them. Therefore, back-projection 3D reconstruction is performed on the 4 candidate combinations to obtain 4 candidate 3D models. Step 5-2. Calculate the reprojection error; reproject each candidate 3D model along the given three observation directions α, β, γ to generate a reprojected image K(i, j); calculate the error S between the reprojected image K(i, j) and the original projection J(i, j) using the following formula: The original projection J(i, j) is the projection map of the 2D reconstructed image used to generate the current candidate 3D reconstruction model in each viewpoint direction (α, β, γ). , , Let K(i, j) represent the mean square error between the reprojected image K(i, j) of the candidate 3D model and the original projection J(i, j) in the three directions α, β, and γ, respectively. The calculation formula is as follows: ； Step 5-3. Select the optimal solution; select the candidate model with the smallest reprojection error S as the final 3D reconstruction result H.

Images