A coaxial digital holographic particle field three-dimensional reconstruction method combining structure tensor and guided filtering
By combining structural tensors and guided filtering, the problems of insufficient focusing feature representation and noise interference in the three-dimensional reconstruction of coaxial digital holographic particle fields are solved, and high-precision and robust reconstruction of particle fields is achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- GUANGDONG UNIV OF TECH
- Filing Date
- 2026-04-01
- Publication Date
- 2026-06-30
AI Technical Summary
Existing coaxial digital holographic particle field 3D reconstruction methods have insufficient ability to focus feature representation during the fusion of multi-distance reconstruction results, and are easily affected by noise and background disturbances, resulting in unstable particle recognition and localization, and problems such as diffraction ring artifact interference, particle adhesion, and repeated detection.
By combining structural tensor and guided filtering methods, spatial consistency fusion of multi-distance reconstruction results is achieved. A focus response map is constructed using structural tensor and guided filtering is performed. Combined with pseudo-loop gating and self-focusing search, robust extraction and precise depth localization of particle candidate regions are realized.
It improves the accuracy and stability of particle 3D reconstruction, reduces diffraction artifact interference, enhances particle target response, and ensures the accuracy of particle depth estimation and the integrity of particle coordinate and size measurement.
Smart Images
Figure CN122312906A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of three-dimensional reconstruction of coaxial holographic particle fields, and more specifically, relates to a method for three-dimensional reconstruction of coaxial digital holographic particle fields that combines structural tensors and guided filtering. Background Technology
[0002] Digital holography is an imaging and measurement technique that reconstructs a target field by using a hologram formed by the interference of object light and reference light, combined with numerical diffraction calculations. Due to its non-contact, full-field recording, and ability to recover three-dimensional information, digital holography has been widely applied in particle measurement, flow field diagnosis, spray analysis, and three-dimensional particle localization. Among these, coaxial digital holography, with its compact optical path structure and simple system implementation, can acquire holographic images containing particle spatial distribution information under single-exposure conditions, thus possessing high application value in three-dimensional particle field reconstruction.
[0003] Existing coaxial digital holographic particle field 3D reconstruction methods typically rely on multi-distance numerical reconstruction to obtain reconstruction results at different axial positions, and combine this with a focusing evaluation function to determine the axial position and spatial distribution of the particles. However, in practical applications, due to the large number of particles, their dense distribution, significant differences in particle size, and the susceptibility of the imaging process to noise and background disturbances, each reconstruction layer often suffers from problems such as defocusing blur, background undulations, and diffraction artifacts. Especially under coaxial imaging conditions, particles easily generate obvious diffraction ring structures, which form a strong response at non-focused locations, thus interfering with the identification and localization of real particles.
[0004] Furthermore, existing methods often employ focusing evaluation based on single features such as grayscale, gradient, variance, or Laplacian, independently judging each reconstructed layer. While these methods are relatively simple to implement, they are prone to problems such as unstable focusing evaluation, large depth estimation errors, and discontinuous extended depth-of-field results when particle edges are unclear, local textures are weak, or noise interference is strong. Meanwhile, during particle candidate region extraction, relying solely on single-threshold segmentation, local extremum detection, or conventional morphological processing can easily lead to false detections, missed detections, and duplicate detections when particles are clustered, adjacent particles are closely spaced, or diffraction rings are prominent, thus affecting the accuracy of particle 3D coordinate reconstruction and size parameter measurement.
[0005] Therefore, it is still necessary to propose a method for three-dimensional reconstruction of coaxial digital holographic particle fields, so as to balance the ability to characterize focusing features and spatial consistency in the process of fusing multi-distance reconstruction results, and effectively suppress diffraction artifacts, improve the separation effect of adjacent particles, and reduce duplicate detection and missed detection in the particle detection and measurement process, thereby improving the accuracy and robustness of particle field three-dimensional reconstruction. Summary of the Invention
[0006] The purpose of this invention is to overcome the shortcomings of the prior art and provide a coaxial digital holographic particle field three-dimensional reconstruction method that combines structural tensor and guided filtering.
[0007] The technical solution of the present invention to solve the above-mentioned technical problems is:
[0008] A method for three-dimensional reconstruction of coaxial digital holographic particle fields combining structure tensor and guided filtering includes the following steps:
[0009] S1. Acquire a coaxial digital holographic intensity image under the coaxial digital holographic imaging optical path, and perform normalization preprocessing on the coaxial digital holographic intensity image to obtain the holographic image A to be reconstructed;
[0010] S2. At multiple preset reconstruction distances along the imaging optical axis, Fresnel diffraction convolution reconstruction is performed on the holographic image A to be reconstructed to obtain the reconstruction intensity image corresponding to each preset reconstruction distance, and a three-dimensional reconstruction intensity stack B is constructed according to the axial order of the preset reconstruction distances.
[0011] S3. Bandpass enhancement is performed on the reconstruction intensity images of each layer in the 3D reconstruction intensity stack B. Based on the bandpass-enhanced reconstruction intensity images of each layer, the focus response map of each layer corresponding to the preset reconstruction distance is constructed through the structure tensor. Then, the guide map generated by the maximum intensity projection of the 3D reconstruction intensity stack B along the preset reconstruction distance dimension is used as the guide image to perform guide filtering on the focus response map of the corresponding preset reconstruction distance layer to obtain a 3D weighted prior stack with spatial consistency and edge preservation characteristics.
[0012] S4. The three-dimensional weight prior stack is axially normalized pixel by pixel in the preset reconstruction distance dimension to obtain the fusion weight of each pixel in multiple preset reconstruction distances; for each pixel, the soft depth map C is calculated by axial weighted average according to the distribution of its fusion weight in the preset reconstruction distance, and the reconstruction intensity value at the preset reconstruction distance corresponding to the maximum fusion weight of each pixel is used as the pixel value corresponding to the extended depth image to generate the extended depth image D.
[0013] S5. Based on the extended depth-of-field image D and the fusion weight of each pixel, construct a candidate score map, and obtain the particle candidate connected domain and the corresponding candidate center set E through double threshold morphological reconstruction, distance transformation, h-maxima transformation, region maximum value label extraction, adaptive label merging and label-controlled watershed segmentation.
[0014] S6. For each candidate center in the candidate center set E, initialize its axial depth position according to the depth value corresponding to the coordinates of the candidate center in the soft depth map C, and adaptively determine the local analysis region by combining the scale information of the corresponding particle candidate connected domain in step S5; then, within the local analysis region, perform pseudo-ring gating based on the response difference between the central region and the ring region to eliminate pseudo-target responses caused by the diffraction ring structure.
[0015] S7. Within the axial neighborhood of the initialized axial depth, the axial gradient variance of the reconstructed image in the local analysis region is used as the focusing evaluation function to determine the precise axial depth of the particle through self-focusing search. Then, using the axial gray variance features of the local reconstruction intensity stack near the precise axial depth of the particle and the gray features of the focusing layer corresponding to the precise axial depth of the particle, binary clustering is performed on the particle target in the local analysis region to obtain the particle mask, and the particle size parameters are calculated based on the particle mask.
[0016] S8. After deep refinement and size measurement, perform deduplication processing based on the spatial proximity of candidate centers and the overlap relationship of particle masks; at the same time, perform omission detection processing in the remaining areas of the image not covered by the detected particle masks, and finally output the three-dimensional spatial coordinates and size parameters F of all valid particles.
[0017] Preferably, in step S1, under the coaxial digital holographic imaging optical path, the image acquisition camera at the end of the optical path acquires the coaxial digital holographic intensity image formed by the interference of the object light and the reference light; the acquired coaxial digital holographic intensity image is subjected to 0-1 linear normalization preprocessing to reduce the influence of light intensity distribution differences and background light disturbances under different acquisition conditions on the subsequent numerical reconstruction results, and a standardized holographic image A to be reconstructed is obtained.
[0018] Preferably, in step S2, multiple preset reconstruction distances are sampled at uniform intervals within a preset axial range along the optical axis of the coaxial digital holographic imaging optical path; the holographic image A to be reconstructed is converted to the frequency domain by fast Fourier transform, and the frequency domain data is multiplied by the propagation kernel spectrum corresponding to each preset reconstruction distance. After inverse fast Fourier transform, the complex amplitude reconstruction result of the corresponding preset reconstruction distance is obtained. The reconstructed intensity image is obtained by taking the modulus and squaring, and stacked to generate a three-dimensional reconstruction intensity stack B.
[0019] Preferably, in step S2, the propagation kernel is a chirplet wavelet kernel, or a diffraction propagation kernel that satisfies the Fresnel diffraction integral condition and has the same diffraction propagation calculation accuracy as the chirplet wavelet kernel. The core parameters of the propagation kernel are jointly determined by the working laser wavelength of the imaging system, the physical pixel size of the image sensor, the corresponding preset reconstruction distance, and the Gaussian window function parameters of the kernel function.
[0020] Preferably, in step S3, the structure tensor serves as a focus evaluation index to assess the sharpness of the reconstructed intensity images corresponding to each preset reconstruction distance in the 3D reconstruction intensity stack B, while suppressing the interference of unidirectional edge structures on the focus evaluation results. The specific steps for calculating the structure tensor and constructing the focus response map are as follows:
[0021] For the reconstructed intensity images of each layer in the 3D reconstruction intensity stack B after bandpass enhancement, Gaussian smoothing is first performed using a Gaussian smoothing kernel of a preset first size to suppress image noise and obtain the denoised reconstructed intensity image.
[0022] Using the Sobel operator, the gray-level gradient in the x and y directions is calculated for each pixel in the reconstructed intensity image after noise reduction.
[0023] Then, for each pixel, construct the outer product matrix of its gray-level gradient in the x-direction and gray-level gradient in the y-direction;
[0024] For the gradient outer product matrix corresponding to each pixel, a Gaussian weighted average operation is performed within a Gaussian window of a preset second size to obtain the structure tensor corresponding to the pixel at the current preset reconstruction distance.
[0025] The structure tensor contains four matrix elements: gradient energy in the x-direction, gradient energy in the y-direction, and two symmetric gradient-related terms. The gradient energy in the x-direction is calculated by averaging the squares of the gray-level gradients in the x-direction and the gray-level gradients in the y-direction, respectively, and the gradient-related terms are calculated by averaging the product of the gray-level gradients in the x-direction and the gray-level gradients in the y-direction.
[0026] Based on the structure tensor corresponding to each pixel, the focus response value of that pixel is calculated. Based on the focus response values of all pixels, a focus response map of the layer corresponding to the current preset reconstruction distance is generated. The focus response value is calculated using the trace of the structure tensor, which is the sum of the gradient energy in the x-direction and the gradient energy in the y-direction of the corresponding pixel. The magnitude of the focus response value is positively correlated with the local focus degree of the corresponding pixel at the current preset reconstruction distance.
[0027] Preferably, in step S3, the generation rule of the three-dimensional weight prior stack is as follows:
[0028] For the 3D reconstruction intensity stack B, along the preset reconstruction distance dimension, take the maximum reconstruction intensity value of each pixel under all preset reconstruction distances to generate a 2D guide map;
[0029] The focus response map of each preset reconstruction distance corresponding to the layer, constructed based on the structure tensor, is used as the filtering input for the corresponding layer;
[0030] Using the generated two-dimensional guide image as the guide image, guided filtering is performed on the input to be filtered corresponding to each preset reconstruction distance layer to obtain the filtered output results with spatial consistency and edge preservation characteristics of each layer. The filtered output results of each layer are stacked in axial order according to the preset reconstruction distance to obtain a three-dimensional weighted prior stack.
[0031] Preferably, in step S4, the steps for generating the fusion weights, soft depth map C, and extended depth-of-field image D are as follows:
[0032] A non-negative constraint is applied to the weight prior values in the three-dimensional weight prior stack. After setting the negative values to zero, axial normalization is performed pixel by pixel along the preset reconstruction distance dimension to obtain the fusion weight of each pixel at each preset reconstruction distance.
[0033] For each pixel, the fusion weight at each preset reconstruction distance is used as the weighting coefficient. The axial position values at the corresponding preset reconstruction distance are weighted and summed to obtain the soft depth value of the pixel. After traversing all pixels, a soft depth map C is generated.
[0034] For each pixel, the preset reconstruction distance layer corresponding to the maximum fusion weight is selected, and the reconstruction intensity value of the corresponding pixel in the layer is taken as the pixel value of the extended depth image. After traversing all pixels, the extended depth image D is synthesized.
[0035] Preferably, in step S5, the construction steps of the candidate score map are as follows:
[0036] The extended depth-of-field image D generated in step S4 is subjected to texture enhancement processing to strengthen the edge and texture features of particle targets in the image, and the texture enhancement result is obtained.
[0037] Multiple sets of ring kernel templates of different scales are preset. The texture enhancement result is convolved and matched with the ring kernel templates of each scale respectively. The matching response at each scale is calculated. The absolute values of the matching responses at all scales are summed to obtain the total multi-scale ring kernel matching response.
[0038] Combining the fusion weights obtained in step S4, the total response of the multi-scale ring kernel matching is subjected to pixel-by-pixel weighted modulation to finally obtain the candidate score map.
[0039] Preferably, in step S6, the execution rules for the pseudo-loop gating are as follows:
[0040] Within the adaptively determined local analysis region, for each candidate center in the candidate center set E, the central region and the outer ring region are divided with the candidate center as the center and the preset particle feature size as the radius.
[0041] Calculate the average response intensity of the central region and the annular region in the candidate score map. If the average response intensity of the annular region is greater than that of the central region, the candidate center is determined to be a pseudo-target response caused by the diffraction ring and is removed.
[0042] Preferably, in step S7, the steps for generating the particle mask and calculating its size parameters are as follows:
[0043] For a particle target that has completed self-focusing depth refinement, within the local analysis area, for each pixel, based on the reconstruction intensity values of all axial reconstruction layers of the local reconstruction intensity stack near the particle's precise depth, the axial gray-scale variance is calculated to obtain the axial gray-scale variance feature of the pixel.
[0044] The axial gray-level variance feature of each pixel is combined with the gray-level feature of the pixel in the focusing layer corresponding to the precise depth of the particle to form a two-dimensional segmentation feature vector for particle segmentation.
[0045] Based on the segmentation feature vectors of all pixels, the local analysis region is segmented by binary clustering of particle targets and background to obtain an initial particle mask; morphological post-processing is performed on the initial particle mask to obtain the final particle mask.
[0046] Based on the particle mask, the number of pixels in the particle target area is counted, the particle area is calculated, and the equivalent diameter of the particle's equal-area circle is calculated based on the particle area.
[0047] Preferably, in step S8, the rules for deduplication and omission detection are as follows:
[0048] Deduplication: If the spatial distance between two candidate centers is less than the preset minimum particle spacing, and the corresponding particle mask overlap is greater than the preset overlap threshold, then retain the particle target with larger particle size and higher focus response, and remove the other duplicate target.
[0049] Missed detection processing: In the remaining areas of the image not covered by the detected particle mask, candidate score map construction and watershed segmentation are re-executed to supplement the detection of missed particle targets.
[0050] Compared with the prior art, the present invention has the following advantages:
[0051] 1. This invention aims to utilize the effective characterization capability of structural tensors for particle focusing features and the edge-preserving properties of guided filtering to construct a spatial consistency fusion method for multi-distance reconstruction results. This method enhances particle target features while suppressing artifacts and background interference. Combined with candidate region localization, pseudo-loop suppression, depth refinement, particle size measurement, and deduplication and supplementation, it achieves high-precision reconstruction of particle three-dimensional spatial coordinates and size parameters.
[0052] 2. This invention constructs a focus response based on a structure tensor on the multi-distance numerical reconstruction results, and combines guided filtering to generate weighted priors with spatial consistency and edge-preserving properties, thereby forming fusion weights and achieving effective fusion of reconstruction information from multiple preset reconstruction distances. Compared to existing methods that independently evaluate each reconstruction layer based solely on single grayscale, gradient, or variance information, this invention can more effectively characterize the focusing features of particles, enhance particle target response, suppress background gradual variation components and artifact interference, thereby improving the quality of extended depth-of-field images and the stability and accuracy of particle depth estimation.
[0053] 3. This invention constructs a candidate score map based on texture enhancement results from extended depth-of-field images, multi-scale ring kernel matching responses, and fusion weights. Combined with a label-controlled watershed segmentation framework, it achieves robust extraction of particle candidate regions and effective separation of adjacent or adherent particles. Compared to existing methods that rely solely on single-threshold segmentation or local extremum detection, this invention reduces the impact of diffraction ring spurious responses and particle adhesion on detection results, improving the accuracy of particle candidate center localization.
[0054] 4. This invention achieves precise reconstruction of particle axial depth by performing pseudo-loop gating within the local analysis region of candidate particles and constructing a self-focusing metric based on gradient variance in the axial depth neighborhood. Combined with hierarchical sampling search to determine the reconstruction position corresponding to the extreme value of the focusing metric, this invention effectively suppresses pseudo-target responses caused by diffraction ring structures, reduces the risk of depth misjudgment, and improves the reliability of precise particle depth localization.
[0055] 5. This invention utilizes the axial gray-level variance features of the locally reconstructed intensity stack near the precise depth of particles and the gray-level features of the focusing layer to construct segmentation features, achieving effective differentiation between particle targets and the background. Combined with a particle mask, it completes the measurement of particle area and equivalent diameter. Simultaneously, through deduplication processing based on the spatial proximity relationship of candidate centers and the overlap relationship of particle masks, and supplementary detection processing for uncovered remaining areas, it reduces the phenomenon of duplicate labeling and missed detection of the same particle. Compared with existing technologies, this invention can simultaneously improve the accuracy of particle 3D positioning and size measurement, enhancing the integrity and robustness of the coaxial digital holographic particle field 3D reconstruction results. Attached Figure Description
[0056] Figure 1 This is a flowchart of the three-dimensional reconstruction method for coaxial digital holographic particle fields that combines structural tensors and guided filtering according to the present invention.
[0057] Figure 2 This is an algorithmic framework diagram of the coaxial digital holographic particle field three-dimensional reconstruction method combining structural tensor and guided filtering according to the present invention. Detailed Implementation
[0058] The present invention will be further described in detail below with reference to the embodiments and accompanying drawings, but the embodiments of the present invention are not limited thereto.
[0059] See Figures 1-2 The present invention provides a method for three-dimensional reconstruction of coaxial digital holographic particle fields by combining structural tensors and guided filtering, which includes the following steps:
[0060] S1. Acquire a coaxial digital holographic intensity image under the coaxial digital holographic imaging optical path, and perform normalization preprocessing on the coaxial digital holographic intensity image to reduce the influence of light intensity distribution differences under different acquisition conditions on subsequent reconstruction results, and obtain the holographic image to be reconstructed A.
[0061] S2. At multiple preset reconstruction distances along the imaging optical axis, Fresnel diffraction convolution reconstruction is performed on the holographic image A to be reconstructed to obtain the reconstruction intensity image corresponding to each preset reconstruction distance, and a three-dimensional reconstruction intensity stack B is constructed according to the axial order of the preset reconstruction distances.
[0062] S3. Bandpass enhancement is applied to the reconstruction intensity images of each layer in the 3D reconstruction intensity stack B to suppress diffraction rings and background gradual components, and to enhance particle-scale texture. Based on the bandpass-enhanced reconstruction intensity images of each layer, a focus response map of the corresponding layer at each preset reconstruction distance is constructed using the structure tensor to characterize the local focus of each pixel at the current preset reconstruction distance. Then, the guide map generated by the maximum intensity projection of the 3D reconstruction intensity stack B along the preset reconstruction distance dimension is used as the guide image to perform guide filtering on the focus response map of the corresponding preset reconstruction distance layer to preserve edges, thereby obtaining a 3D weighted prior stack with spatial consistency and edge preservation characteristics.
[0063] S4. The three-dimensional weight prior stack is axially normalized pixel by pixel in the preset reconstruction distance dimension to obtain the fusion weight of each pixel in multiple preset reconstruction distances; for each pixel, according to the distribution of its fusion weight in the preset reconstruction distance, the soft depth map C is calculated by axial weighted averaging, and the reconstruction intensity value at the preset reconstruction distance corresponding to the maximum fusion weight of each pixel is used as the pixel value corresponding to the extended depth image to generate the extended depth image D.
[0064] S5. Construct a candidate score map based on the extended depth-of-field image D and the fusion weight of each pixel; extract particle candidate regions using dual-threshold morphological reconstruction; then perform distance transformation on the particle candidate regions and h-maxima transformation on the distance transformation results to suppress low-significance local maxima; subsequently, combine region maximum value label extraction, adaptive label merging and label-controlled watershed segmentation to obtain the candidate center set E, so as to achieve the initial separation and localization of adjacent particles or adherent particles;
[0065] S6. For each candidate center in the candidate center set E, initialize its axial depth position according to the depth value corresponding to the coordinates of the candidate center in the soft depth map C, and adaptively determine the local analysis region by combining the scale information of the corresponding particle candidate connected domain in step S5; then, within the local analysis region, perform pseudo-ring gating based on the response difference between the central region and the ring region to eliminate pseudo-target responses caused by the diffraction ring structure.
[0066] S7. Within the axial neighborhood corresponding to the initial axial depth, construct a self-focusing metric using the axial gradient variance of the reconstructed image within the local analysis region; then, use a hierarchical search strategy to traverse the reconstruction positions within the neighborhood, retrieve the reconstruction coordinates corresponding to the extreme value of the focusing metric, and determine the coordinates as the precise axial depth of the corresponding particle; after determining the precise axial depth of the particle, use the axial gray variance features of the local reconstruction intensity stack near the precise axial depth of the particle and the gray features of the focusing layer corresponding to the precise axial depth of the particle to perform binary clustering segmentation of the particle targets within the local analysis region, obtain the particle mask, and calculate the particle size parameters based on the particle mask;
[0067] S8. After depth refinement and size measurement, perform deduplication processing based on the spatial proximity of candidate centers and the overlap of particle masks to eliminate duplicate markers of the same particle in adjacent reconstruction layers or adjacent candidate regions; at the same time, perform omission detection processing in the remaining areas of the image not covered by the detected particle mask, and finally output the three-dimensional spatial coordinates and size parameters F of all effective particles.
[0068] See Figures 1-2 In step S1, under the coaxial digital holographic imaging optical path, the image acquisition camera at the end of the optical path acquires the coaxial digital holographic intensity image formed by the interference of the object light and the reference light; the acquired coaxial digital holographic intensity image is subjected to 0-1 linear normalization preprocessing to reduce the influence of light intensity distribution differences and background light disturbances under different acquisition conditions on the subsequent numerical reconstruction results, and a standardized holographic image A to be reconstructed is obtained.
[0069] See Figures 1-2 In step S2, multiple preset reconstruction distances Sampling is performed at uniform intervals within a preset axial range along the optical axis of the coaxial digital holographic imaging optical path. The holographic image A to be reconstructed is converted to the frequency domain using a Fast Fourier Transform. The frequency domain data is multiplied by the propagation kernel spectrum corresponding to each preset reconstruction distance, and then subjected to an Inverse Fast Fourier Transform to obtain the complex amplitude reconstruction result for the corresponding preset reconstruction distance. The squared value is then used to obtain the reconstructed intensity image. These images are stacked in axial order according to the preset reconstruction distances to generate a three-dimensional reconstructed intensity stack B. The k-th layer reconstructed intensity image of the three-dimensional reconstructed intensity stack B is denoted as... In addition, the propagation kernel adopts the chirplet wavelet kernel, or a diffraction propagation kernel that satisfies the Fresnel diffraction integral condition and has the same diffraction propagation calculation accuracy as the chirplet wavelet kernel. The core parameters of the propagation kernel are jointly determined by the working laser wavelength of the imaging system, the physical pixel size of the image sensor, the corresponding preset reconstruction distance, and the Gaussian window function parameters of the kernel function.
[0070] See Figures 1-2 In step S3, the structure tensor serves as a focus evaluation index, used to assess the clarity of the reconstructed intensity images corresponding to each preset reconstruction distance in the 3D reconstruction intensity stack B, while suppressing the interference of unidirectional edge structures such as diffraction rings on the focus evaluation results. The specific steps for calculating the structure tensor and constructing the focus response map are as follows:
[0071] For the reconstructed intensity images of each layer in the 3D reconstruction intensity stack B after bandpass enhancement processing , First, Gaussian smoothing is performed using a Gaussian smoothing kernel of a preset first size to suppress image noise, resulting in a denoised reconstructed intensity image. Then, using the Sobel operator, the gray-level gradients in the x and y directions are calculated for each pixel in the denoised reconstructed intensity image. Next, for each pixel, an outer product matrix of its x-direction and y-direction gray-level gradients is constructed. For each pixel's gradient outer product matrix, a Gaussian weighted average operation is performed within a Gaussian window of a preset second size to obtain the structure tensor corresponding to that pixel at the current preset reconstruction distance. The structure tensor contains four matrix elements: x-direction gradient energy, y-direction gradient energy, and two symmetric gradient correlation terms. The x-direction gradient energy is the result of Gaussian weighted averaging of the square of the x-direction gray-level gradient, the y-direction gradient energy is the result of Gaussian weighted averaging of the square of the y-direction gray-level gradient, and the gradient correlation terms are the result of Gaussian weighted averaging of the product of the x-direction and y-direction gray-level gradients.
[0072] In this embodiment, the expression for the structure tensor is:
[0073] ;
[0074] in, This represents the weighted average operation within the Gaussian window. For the structure tensor at different depths, for Directional gradient energy, for Directional gradient energy, This is a gradient-related term, reflecting the directional consistency of the local structure.
[0075] Based on the structure tensor corresponding to each pixel, the focus response value of that pixel is calculated. Based on the focus response values of all pixels, a focus response map of the layer corresponding to the current preset reconstruction distance is generated. The focus response value is calculated using the trace of the structure tensor, which is the sum of the gradient energy in the x-direction and the gradient energy in the y-direction of the corresponding pixel. The magnitude of the focus response value is positively correlated with the local focus degree of the corresponding pixel at the current preset reconstruction distance.
[0076] See Figures 1-2 In step S3, the generation rules for the three-dimensional weight prior stack are as follows:
[0077] For the 3D reconstruction intensity stack B, along the preset reconstruction distance dimension, the maximum reconstruction intensity value of each pixel at all preset reconstruction distances is taken to generate a 2D guide map; wherein, the 2D guide map can be specifically represented as:
[0078] ;
[0079] in, It is a two-dimensional guide diagram; Indicates that the 3D reconstruction intensity stack B is at the 3D level. The reconstruction intensity value corresponding to each preset reconstruction distance;
[0080] The focus response map of each preset reconstruction distance corresponding to the layer, constructed based on the structure tensor, is used as the filtering input for the corresponding layer;
[0081] Using the generated 2D guide image as the guide image, guided filtering is performed on the input to be filtered corresponding to each preset reconstruction distance layer, resulting in filtered outputs with spatial consistency and edge preservation characteristics for each layer, i.e.:
[0082] ;
[0083] in, Indicates the first Focus response map corresponding to a preset reconstruction distance Indicates For a two-dimensional guide diagram, Perform guided filtering. Indicates the first The weighted prior graph corresponding to each reconstruction distance;
[0084] Finally, the filter outputs of each layer are stacked in axial order according to the preset reconstruction distance to obtain the three-dimensional weighted prior stack.
[0085] See Figures 1-2 In step S4, the steps for generating the fusion weights, soft depth map C, and extended depth-of-field image D are as follows:
[0086] A non-negative constraint is applied to the weight prior values in the three-dimensional weight prior stack. After setting the negative values to zero, axial normalization is performed pixel by pixel along the preset reconstruction distance dimension to obtain the fusion weight of each pixel at each preset reconstruction distance. This fusion weight is obtained by applying a non-negative constraint to the weight prior and normalizing it pixel by pixel along the reconstruction distance, and can be specifically expressed as follows:
[0087] ;
[0088] in, Indicates the first Weighted prior maps corresponding to preset reconstruction distances This represents the normalized fusion weights. This indicates the total number of preset reconstruction distances;
[0089] For each pixel, its fusion weight at each preset reconstruction distance is used as a weighting coefficient. The axial position values at the corresponding preset reconstruction distance are then summed using weighted averages to obtain the soft depth value of that pixel. After traversing all pixels, a soft depth map C is generated. This soft depth map C can be represented as:
[0090] ;
[0091] in, Indicates the first A preset reconstruction distance, Represents pixels The soft depth value at that location;
[0092] For each pixel, the preset reconstruction distance layer corresponding to the maximum fusion weight is selected, and the reconstruction intensity value of the corresponding pixel in that layer is taken as the pixel value of the extended depth image. After traversing all pixels, the extended depth image D is synthesized, which can be represented as:
[0093] ;
[0094] in, .
[0095] See Figures 1-2 In step S5, the construction of the candidate score map and the acquisition of the particle candidate center set E are as follows:
[0096] The extended depth-of-field image D generated in step S4 is subjected to texture enhancement processing to strengthen the edge and texture features of the particle target; multiple sets of ring kernel templates of different scales are preset, and the texture enhancement results are convolved and matched with the ring kernel templates of each scale respectively to calculate the matching response at each scale. The absolute values of the matching responses at all scales are accumulated to obtain the total multi-scale ring kernel matching response.
[0097] Combining the fusion weights obtained in step S4, the total response of the multi-scale ring kernel matching is weighted and modulated pixel by pixel to finally obtain the candidate score map, the mathematical expression of which is:
[0098] ;
[0099] in, For candidate score graphs, To extend the depth map, To integrate weights, To extend the depth-of-field map texture enhancement result, For the first Ring core templates at various scales For scale sets, This is a convolution operation;
[0100] Figure 2 The “Score(maskedbybestW)” in the figure illustrates the response distribution characteristics of the candidate score map after the fusion weight modulation.
[0101] Based on this, candidate regions of particles are extracted through dual-threshold morphological reconstruction; then, distance transformation is performed on the candidate regions, and h-maxima transformation is used to suppress local maxima with low significance to generate marker points. These marker points are then used to perform marker-controlled watershed segmentation to obtain the candidate center set E, thereby achieving the separation and localization of adjacent or adherent particles.
[0102] See Figures 1-2 In step S6, pseudo-ring gating is achieved by comparing the response differences between the central region and the ring region within the candidate center's neighborhood, in order to eliminate pseudo-target responses caused by the diffraction ring structure. Specifically:
[0103] For each candidate center in the candidate center set E, the central region and the outer ring region are divided with the candidate center as the center and the preset particle feature size as the radius. The average response intensity of the central region and the ring region in the candidate score map are calculated respectively. If the average response intensity of the ring region is greater than that of the central region, the candidate center is determined to be a false target response caused by the diffraction ring and is removed.
[0104] After the false target removal is completed, the self-focusing depth refinement process is entered. Specifically, the self-focusing metric adopts the gradient variance metric based on the gradient operator in the local analysis region, and hierarchical sampling search is used in the axial depth neighborhood to determine the reconstruction position corresponding to the extreme value of the focusing metric, which is used as the precise depth of the particle.
[0105] See Figures 1-2 In step S7, the process of particle mask generation and size parameter calculation is as follows:
[0106] For a particle target that has completed self-focusing depth refinement, within the local analysis area, for each pixel, based on the reconstruction intensity values of all axial reconstruction layers of the local reconstruction intensity stack near the particle's precise depth, the axial gray-scale variance is calculated to obtain the axial gray-scale variance feature of the pixel.
[0107] ;
[0108] in, To reconstruct the number of floors, For the first The axial depth position corresponding to the reconstructed layer. For pixels Average reconstruction intensity on each reconstruction layer;
[0109] The aforementioned axial grayscale variance features can highlight the differences in axial variation of particles near the focal point and suppress the weak response of the background; at the same time, they can extract the grayscale values of pixels at the precise depth of the particles corresponding to the focal layer. As the grayscale feature of the focusing layer; the axial grayscale variance feature of each pixel is combined with the grayscale feature of the focusing layer corresponding to the precise depth of the particle to form a two-dimensional segmentation feature vector for particle segmentation.
[0110] Based on the two-dimensional segmentation feature vectors of all pixels, the local analysis region is segmented by binary clustering of particle targets and background to obtain an initial particle mask; morphological post-processing is performed on the initial particle mask to obtain the final particle mask.
[0111] Based on the particle mask, the number of pixels in the particle target area is counted, the particle area is calculated, and the equivalent diameter of the particle's equal-area circle is calculated based on the particle area to complete the quantization of particle size.
[0112] See Figures 1-2 In step S8, the particle results that have completed depth refinement and size measurement are subjected to deduplication and omission detection processing in sequence. Deduplication is achieved by relying on the spatial proximity relationship of candidate centers and the overlap relationship of particle masks. Omission detection processing is only performed in the remaining areas of the image that are not covered by the detected particle mask. Finally, the three-dimensional spatial coordinates and size parameters F of all effective particles are output. Figure 2 This shows the three-dimensional distribution of particles and the superposition effect of detection results;
[0113] The specific execution rules for deduplication and omission detection are as follows:
[0114] Deduplication: If the spatial distance between two candidate centers is less than the preset minimum particle spacing, and the overlap of the corresponding particle masks is greater than the preset overlap threshold, then the particle target with the larger size and higher focus response is retained, and duplicate targets are removed.
[0115] Missed detection processing: In the remaining areas of the image not covered by the detected particle mask, the candidate score map is reconstructed and the watershed segmentation is performed to detect the missed particle targets.
[0116] The above are preferred embodiments of the present invention, but the embodiments of the present invention are not limited to the above content. Any changes, modifications, substitutions, combinations, or simplifications made without departing from the spirit and principle of the present invention shall be considered equivalent substitutions and shall be included within the protection scope of the present invention.
Claims
1. A coaxial digital holographic particle field three-dimensional reconstruction method combining structure tensor and guided filtering, characterized by, Includes the following steps: S1. Acquire a coaxial digital holographic intensity image under the coaxial digital holographic imaging optical path, and perform normalization preprocessing on the coaxial digital holographic intensity image to obtain the holographic image A to be reconstructed; S2. At multiple preset reconstruction distances along the imaging optical axis, Fresnel diffraction convolution reconstruction is performed on the holographic image A to be reconstructed to obtain the reconstruction intensity image corresponding to each preset reconstruction distance, and a three-dimensional reconstruction intensity stack B is constructed according to the axial order of the preset reconstruction distances. S3. Bandpass enhancement is performed on the reconstruction intensity images of each layer in the 3D reconstruction intensity stack B. Based on the bandpass-enhanced reconstruction intensity images of each layer, the focus response map of each layer corresponding to the preset reconstruction distance is constructed through the structure tensor. Then, the guide map generated by the maximum intensity projection of the 3D reconstruction intensity stack B along the preset reconstruction distance dimension is used as the guide image to perform guide filtering on the focus response map of the corresponding preset reconstruction distance layer to obtain a 3D weighted prior stack with spatial consistency and edge preservation characteristics. S4. Perform axial normalization on the three-dimensional weight prior stack pixel by pixel in the preset reconstruction distance dimension to obtain the fusion weight of each pixel in multiple preset reconstruction distances. For each pixel, a soft depth map C is calculated by axial weighted average based on the distribution of its fusion weight over a preset reconstruction distance. The reconstruction intensity value at the preset reconstruction distance corresponding to the maximum fusion weight of each pixel is used as the pixel value corresponding to the extended depth image to generate an extended depth image D. S5. Based on the extended depth-of-field image D and the fusion weight of each pixel, construct a candidate score map, and extract the particle candidate connected components and the corresponding candidate center set E through image segmentation, label extraction and merging processing; S6. For each candidate center in the candidate center set E, initialize its axial depth position according to the depth value corresponding to the coordinates of the candidate center in the soft depth map C, and adaptively determine the local analysis region by combining the scale information of the corresponding particle candidate connected domain in step S5; then, within the local analysis region, perform pseudo-ring gating based on the response difference between the central region and the ring region to eliminate pseudo-target responses caused by the diffraction ring structure. S7. Within the axial neighborhood of the initialized axial depth, the axial gradient variance of the reconstructed image in the local analysis region is used as the focusing evaluation function to determine the precise axial depth of the particle through self-focusing search. Then, using the axial gray variance features of the local reconstruction intensity stack near the precise axial depth of the particle and the gray features of the focusing layer corresponding to the precise axial depth of the particle, binary clustering is performed on the particle target in the local analysis region to obtain the particle mask, and the particle size parameters are calculated based on the particle mask. S8. After deep refinement and size measurement, perform deduplication processing based on the spatial proximity of candidate centers and the overlap relationship of particle masks; at the same time, perform omission detection processing in the remaining areas of the image not covered by the detected particle masks, and finally output the three-dimensional spatial coordinates and size parameters F of all valid particles.
2. The method according to claim 1, wherein, In step S1, under the coaxial digital holographic imaging optical path, the image acquisition camera at the end of the optical path acquires the coaxial digital holographic intensity image formed by the interference of the object light and the reference light; the acquired coaxial digital holographic intensity image is subjected to 0-1 linear normalization preprocessing to reduce the influence of light intensity distribution differences and background light disturbances under different acquisition conditions on the subsequent numerical reconstruction results, and a standardized holographic image A to be reconstructed is obtained.
3. The method for three-dimensional reconstruction of coaxial digital holographic particle fields combining structural tensor and guided filtering according to claim 2, characterized in that, In step S2, multiple preset reconstruction distances are sampled at uniform intervals within a preset axial range along the optical axis of the coaxial digital holographic imaging optical path; the holographic image A to be reconstructed is converted to the frequency domain by fast Fourier transform, and the frequency domain data is multiplied by the propagation kernel spectrum corresponding to each preset reconstruction distance. After inverse fast Fourier transform, the complex amplitude reconstruction result of the corresponding preset reconstruction distance is obtained. The modulus is squared to obtain the reconstruction intensity image, and the images are stacked to generate a three-dimensional reconstruction intensity stack B.
4. The method for three-dimensional reconstruction of coaxial digital holographic particle fields combining structural tensor and guided filtering according to claim 3, characterized in that, In step S2, the propagation kernel is a chirplet wavelet kernel, or a diffraction propagation kernel that satisfies the Fresnel diffraction integral condition and has the same diffraction propagation calculation accuracy as the chirplet wavelet kernel. The core parameters of the propagation kernel are jointly determined by the working laser wavelength of the imaging system, the physical pixel size of the image sensor, the corresponding preset reconstruction distance, and the Gaussian window function parameters of the kernel function.
5. The coaxial digital holographic particle field three-dimensional reconstruction method combining structural tensor and guided filtering according to claim 4, characterized in that, In step S3, the structure tensor serves as a focus evaluation index, used to assess the sharpness of the reconstructed intensity images corresponding to each preset reconstruction distance in the 3D reconstruction intensity stack B, while suppressing the interference of unidirectional edge structures on the focus evaluation results. The specific steps for calculating the structure tensor and constructing the focus response map are as follows: For the reconstructed intensity images of each layer in the 3D reconstruction intensity stack B after bandpass enhancement, Gaussian smoothing is first performed using a Gaussian smoothing kernel of a preset first size to suppress image noise and obtain the denoised reconstructed intensity image. Using the Sobel operator, the gray-level gradient in the x and y directions is calculated for each pixel in the reconstructed intensity image after noise reduction. Then, for each pixel, construct the outer product matrix of its gray-level gradient in the x-direction and gray-level gradient in the y-direction; For the gradient outer product matrix corresponding to each pixel, a Gaussian weighted average operation is performed within a Gaussian window of a preset second size to obtain the structure tensor corresponding to the pixel at the current preset reconstruction distance. The structure tensor contains four matrix elements: gradient energy in the x-direction, gradient energy in the y-direction, and two symmetric gradient-related terms. The gradient energy in the x-direction is calculated by averaging the squares of the gray-level gradients in the x-direction and the gray-level gradients in the y-direction, respectively, and the gradient-related terms are calculated by averaging the product of the gray-level gradients in the x-direction and the gray-level gradients in the y-direction. Based on the structure tensor corresponding to each pixel, the focus response value of that pixel is calculated. Based on the focus response values of all pixels, a focus response map of the layer corresponding to the current preset reconstruction distance is generated. The focus response value is calculated using the trace of the structure tensor, which is the sum of the gradient energy in the x-direction and the gradient energy in the y-direction of the corresponding pixel. The magnitude of the focus response value is positively correlated with the local focus degree of the corresponding pixel at the current preset reconstruction distance.
6. The method for three-dimensional reconstruction of coaxial digital holographic particle fields combining structural tensor and guided filtering according to claim 5, characterized in that, In step S3, the generation rules for the three-dimensional weight prior stack are as follows: For the 3D reconstruction intensity stack B, along the preset reconstruction distance dimension, take the maximum reconstruction intensity value of each pixel under all preset reconstruction distances to generate a 2D guide map; The focus response map of each preset reconstruction distance corresponding to the layer, constructed based on the structure tensor, is used as the filtering input for the corresponding layer; Using the generated two-dimensional guide image as the guide image, guided filtering is performed on the input to be filtered corresponding to each preset reconstruction distance layer to obtain the filtered output results with spatial consistency and edge preservation characteristics of each layer. The filtered output results of each layer are stacked in axial order according to the preset reconstruction distance to obtain a three-dimensional weighted prior stack.
7. The method for three-dimensional reconstruction of coaxial digital holographic particle fields combining structural tensor and guided filtering according to claim 6, characterized in that, In step S4, the steps for generating the fusion weights, soft depth map C, and extended depth-of-field image D are as follows: A non-negative constraint is applied to the weight prior values in the three-dimensional weight prior stack. After setting the negative values to zero, axial normalization is performed pixel by pixel along the preset reconstruction distance dimension to obtain the fusion weight of each pixel at each preset reconstruction distance. For each pixel, the fusion weight at each preset reconstruction distance is used as the weighting coefficient. The axial position values at the corresponding preset reconstruction distance are weighted and summed to obtain the soft depth value of the pixel. After traversing all pixels, a soft depth map C is generated. For each pixel, the preset reconstruction distance layer corresponding to the maximum fusion weight is selected, and the reconstruction intensity value of the corresponding pixel in the layer is taken as the pixel value of the extended depth image. After traversing all pixels, the extended depth image D is synthesized.
8. The method for three-dimensional reconstruction of coaxial digital holographic particle fields combining structural tensor and guided filtering according to claim 7, characterized in that, In step S5, the construction steps of the candidate score map are as follows: The extended depth-of-field image D generated in step S4 is subjected to texture enhancement processing to strengthen the edge and texture features of particle targets in the image, and the texture enhancement result is obtained. Multiple sets of ring kernel templates of different scales are preset. The texture enhancement result is convolved and matched with the ring kernel templates of each scale respectively. The matching response at each scale is calculated. The absolute values of the matching responses at all scales are summed to obtain the total multi-scale ring kernel matching response. Combining the fusion weights obtained in step S4, the total response of the multi-scale ring kernel matching is subjected to pixel-by-pixel weighted modulation to finally obtain the candidate score map.
9. The method for three-dimensional reconstruction of coaxial digital holographic particle fields combining structural tensor and guided filtering according to claim 8, characterized in that, In step S6, the execution rules for the pseudo-loop gating are as follows: Within the adaptively determined local analysis region, for each candidate center in the candidate center set E, the central region and the outer ring region are divided with the center as the center and the preset particle feature size as the radius. Calculate the average response intensity of the central region and the annular region in the candidate score map. If the average response intensity of the annular region is greater than that of the central region, the candidate center is determined to be a pseudo-target response caused by the diffraction ring and is removed.
10. The method for three-dimensional reconstruction of coaxial digital holographic particle fields combining structural tensor and guided filtering according to claim 9, characterized in that, In step S7, the steps for generating the particle mask and calculating its size parameters are as follows: For a particle target that has completed self-focusing depth refinement, within the local analysis area, for each pixel, based on the reconstruction intensity values of all axial reconstruction layers of the local reconstruction intensity stack near the particle's precise depth, the axial gray-scale variance is calculated to obtain the axial gray-scale variance feature of the pixel. The axial gray-level variance feature of each pixel is combined with the gray-level feature of the pixel in the focusing layer corresponding to the precise depth of the particle to form a two-dimensional segmentation feature vector for particle segmentation. Based on the segmentation feature vectors of all pixels, the local analysis region is segmented by binary clustering of particle targets and background to obtain an initial particle mask; morphological post-processing is performed on the initial particle mask to obtain the final particle mask. Based on the particle mask, the number of pixels in the particle target area is counted, the particle area is calculated, and the equivalent diameter of the particle's equal-area circle is calculated based on the particle area.