Space-time joint related single-photon counting three-dimensional imaging method

Abstract

The present invention discloses a spatio-temporal joint correlation single-photon counting three-dimensional imaging method. First, pixel-by-pixel preprocessing is performed on the image. According to the temporal correlation of signal photons, photon aggregation clusters are formed, photons within the screening threshold range are selected, and a small set of one aggregation cluster per pixel is constituted, thereby filtering out a large number of noise photons, compressing the data volume, and reducing the computational cost. Secondly, according to the neighborhood spatial correlation, the empty pixels in the entire image are filled to obtain a set of three-dimensional spatial point cloud data. Finally, a cost function of the first-order difference regularization terms in the horizontal and vertical directions of the neighborhood pixel abscissa is established, and the alternating direction multiplier method is used to solve the minimum value of the cost function to obtain an accurate depth estimate. The imaging effect of the present invention is superior to the traditional maximum likelihood estimation method and the time-correlated photon fast denoising method. The depth image obtained after processing has a small error and a high signal-to-noise ratio, improving the three-dimensional imaging quality.

Description

Technical Field

[0001] The present invention belongs to the field of data processing, and specifically relates to a spatio-temporal joint correlation single-photon counting three-dimensional imaging method. Background Art

[0002] Under the detection background of short integration time and high noise, long-distance photon counting lidar systems usually record large data volumes but extremely sparse data cubes. Signal photons are submerged in a large number of noise photons. The depth maps obtained by traditional image estimation methods have a large deviation from the true values, and it is difficult to accurately obtain the information of the detection target. Therefore, there is a challenge to estimate accurate depth images in the case of photon deficiency with less memory requirement, lower computational complexity, and faster imaging speed.

[0003] The alternating direction method of multipliers (ADMM) is a computational framework for solving optimization problems, applicable to solving distributed convex optimization problems of functions, especially statistical learning problems. The ADMM decomposes the large global problem into multiple smaller and easier-to-solve local subproblems through a decomposition and coordination process, and obtains the solution of the large global problem by coordinating the solutions of the subproblems.

[0004] Compared with traditional lidars, photon counting lidars using the principle of time of flight of photons can obtain more accurate depth information and reflectivity images of detection targets. In the military and civilian fields, the development of three-dimensional photon counting lidars has become mature, and the application scope is becoming more and more extensive. Therefore, while researching and developing three-dimensional lidar systems, researchers have proposed many effective methods to improve the quality of reconstructed images, which complement the hardware system.

[0005] The traditional Maximum Likelihood Estimation (MLE) algorithm is usually used in cases with a large amount of data or a high signal-to-noise ratio. However, in the case of short integration times or high background noise, the image imaging effect is not good. To solve the contradiction between the large measurement error caused by short integration times and the slow measurement speed caused by long integration times, Ahmed Kirmani et al. proposed a method of using the first detected photon for imaging in the case of low photon counts [Kirmani A, Venkatraman D, Shin D, et al. First-Photon Imaging[J]. Science, 2014, 343(6166):58-61.]; Donggeek Shin proposed a method applicable to array single-photon detectors, using a denoising method that combines the spatial and temporal characteristics of horizontal and vertical pixels [D Shin, Xu F, D Venkatraman, et al. Photon-efficient imaging with a single-photon camera[J]. Nature Communications, 2016, 7:12046.]. Feng et al. proposed a method for fast time-correlated photon denoising [Research on Key Technologies of Fast Scanning Photon Counting Lidar Imaging[D]. Nanjing University of Science and Technology, 2018.]. These methods have achieved good performance improvement by reducing the detection false alarm rate or increasing the image reconstruction accuracy rate. However, in the case of a large influence of detection environmental noise, the performance and efficiency will be limited. Summary of the Invention

[0006] The object of the present invention is to propose a spatio-temporal joint correlation single-photon counting three-dimensional imaging method.

[0007] The technical solution for realizing the object of the present invention is as follows: A spatio-temporal joint correlation single-photon counting three-dimensional imaging method, the steps are as follows:

[0008] Step 1, using a long-distance photon counting lidar system, collecting point cloud data, and preprocessing the image data;

[0009] Step 2, for the preprocessed data, adaptively obtain the photon flight time of each pixel point to further filter out noise;

[0010] Step 3, complete empty pixels according to neighborhood information;

[0011] Step 4, add the depth information obtained in Step 3 to the total variation regularization term to construct a cost function, and use the alternating direction multiplier method to iteratively update to obtain the optimal solution as the depth information estimation value.

[0012] Preferably, the preprocessing of the image data includes: setting a gating threshold according to prior information and filtering out obvious noise points in the image.

[0013] Preferably, for the preprocessed data, the specific method for adaptively obtaining the photon flight time of each pixel point and further filtering noise is as follows:

[0014] At the pixel point (i, j), sort the flight times of the photons, group L photons into a photon unit, sort the differences of the flight times of adjacent photons in the n-neighborhood, calculate the deviation Δi of adjacent photon units, and find the minimum value of the deviation Δi min of the corresponding photon unit:

[0015]

[0016] where is the flight time of each photon unit at the pixel point (i, j);

[0017] Use the flight time of the unit at the position l corresponding to the deviation Δi of adjacent photon units min as the initial estimated depth t of the target depth initial :

[0018]

[0019] For the photon flight time, determine whether the next group of photons is a signal photon. If the difference between the flight time of the next group of signal photons and the initial estimated depth is less than the threshold T p :

[0020]

[0021] then this type of photon unit is called a signal photon unit, and record the flight times of the photons that meet the conditions in the set . Select a threshold K. If the number of photons in the set at the pixel point (i, j) reaches the calibration value K, or the difference between the photon flight time and the initial depth information is greater than the threshold T p range, at this time the determination of this pixel point ends, and proceed to the photon determination of the next pixel point.

[0022] Preferably, the specific method for filling in empty pixels according to neighborhood information is as follows: Process the image pixel by pixel. At the pixel point (i, j), fill in the empty pixel points with the mean value according to the neighborhood information of its surrounding 8 points.

[0023] Preferably, for a pixel (i, j), if the number of marked non-empty pixels among the eight nearest neighboring pixels (i-1, j-1), (i-1, j), (i-1, j+1), (i, j-1), (i, j+1), (i+1, j-1), (i+1, j), (i+1, j+1) is not less than 3, it is denoted as a small empty pixel; otherwise, it is denoted as a large hole pixel. The small pixel is filled with the mean value of 8-neighborhood non-empty pixels, and the large empty pixel is regarded as an empty pixel.

[0024] Preferably, the specific method for iteratively updating to obtain the optimal solution using the alternating direction method of multipliers is as follows:

[0025] At pixel (i, j), the photon count follows a Poisson distribution, and joint likelihood estimation is performed:

[0026]

[0027] where s i,j,t = r i,j g0(t - t i,j ) + b i,j + d i,j , s i,j,t represents all the photon signals received by the detector, y i,j,t is the lidar observation model, t i,j ≥0 represents the distance between the detector and the target surface, i.e., the depth information, r i,j ≥0 represents the intensity of the target, b i,j ≥0 represents the background and dark photon energy levels, g0 represents the impulse response of the system, t, r, b are N×1 vectors, where N = N c ×N r , T is the time bin, and the cost function is obtained by minimizing the negative log-likelihood function:

[0028]

[0029] where const represents a constant that does not affect the cost function estimation. Assuming the background noise b i,j = 0 and the instrument impulse response is approximately Gaussian distributed where c1 is a constant, σ is a hyperparameter of the system, and is a constant. The depth and intensity estimates of the image are obtained according to the maximum likelihood estimation method, and the cost function is simplified to:

[0030]

[0031] where, Denote the constant term that does not affect the cost function estimation. The total variation regularization term is adopted to improve image estimation. According to the neighborhood spatial correlation, the cost function C after total variation can be obtained TV (t, r):

[0032] C TV (t, r) = L(t, r) + τ1TV(t) + τ2TV(r)

[0033] where τ1 and τ2 are regularization parameters, where and are the first-order differences of the abscissa and ordinate at pixel n. To minimize the cost function, the cost function is minimized by the alternating direction method of multipliers. The convex optimization problem of the alternating multiplier method is as follows:

[0034]

[0035] where the joint estimation of depth and intensity is set as a parameter H (j) is a selection matrix used to select the depth and intensity data in the joint estimation;

[0036] Here, assume a full-rank matrix M as:

[0037]

[0038] where j = 1, 2, 3, 4, and the cost function is split to obtain:

[0039] H (1) = [0 N , K]

[0040] H (2) = [K, 0 N

[0041] g3(u (3) = τ1TV(z1) + τ2TV(z2), H (3) = I 2N

[0042] H (4) = [I N , 0 N

[0043] where g1, g2, g3, and g4 are the intensity, depth, TV regularization, and unnecessary constant terms corresponding to the images after the cost function is split, respectively, where u (3) 、u (4) are the solutions for each update, r​​i ML , Calculated by the maximum likelihood estimation method, I N , K N is the identity matrix, and the alternating direction method of multipliers is used to iteratively update each term:

[0044]

[0045]

[0046]

[0047] Among them, s is the variable value of the obtained minimum value, and k is the number of iterative updates, that is, when the updated is less than the threshold Threshold, the iterative loop is exited to obtain the depth estimate of the target.

[0048] Compared with the prior art, the significant advantages of the present invention are as follows: the present invention utilizes the time correlation between detection data and the spatial correlation of the target to be measured to reduce the error in reconstructing the depth image of the target. According to the different characteristics of signal photons and noise photons in the detector response in the time domain, the required signal photons are extracted, the noise photons are separated, and then based on the Poisson probability distribution model of single photons and the target spatial correlation, the depth image of the target is reconstructed through the extracted signal photons; the imaging effect of the present invention is better than that of the traditional maximum likelihood estimation method and the time-correlated photon fast denoising method, the error of the obtained depth image after processing is small, the signal-to-noise ratio is high, and the three-dimensional imaging quality is improved.

[0049] The present invention will be further described in detail below with reference to the accompanying drawings. Description of the Drawings

[0050] Figure 1 is a schematic diagram of the overall process of the method for adaptive image denoising and reconstruction based on photon counting lidar.

[0051] Figure 2 is a schematic diagram of the logic of the three-dimensional information perception system of the long-distance photon counting lidar supporting the experimental method in the present invention.

[0052] Figure 3 Schematic diagram of the lidar system supporting the experimental method in the present invention.

[0053] Figure 4 is a schematic diagram of the process of the present invention for obtaining the target depth using time correlation.

[0054] Figure 5 is a schematic diagram of the process of the present invention for adding the depth information to the total variation regularization term to form a cost function and using the alternating direction method of multipliers to iteratively update to obtain the optimal solution.

[0055] Figure 6 For verifying the experimental results of a spatio-temporal joint correlation single-photon counting three-dimensional imaging method of the present invention, Figure 6 (a) is an experimental target scene diagram for verifying the present invention, Figure 6 (b) is an original point cloud diagram of the experimental target for verifying the present invention, Figure 6 (c) is a depth map processed manually and used as a reference value for comparison of different methods, Figure 6 (d) is a target depth image processed by an imaging method of a traditional maximum likelihood estimation algorithm when the sampling integration time is 20 ms, Figure 6 (e) is a target depth image processed by an imaging method of time-domain denoising proposed by Feng when the sampling integration time is 20 ms, Figure 6 (f) is a target depth image processed by the imaging method proposed by the present invention when the sampling integration time is 20 ms. Detailed implementation manners

[0056] Some technical details of the present invention will be further described in detail below. The present invention will be further described in detail with reference to the accompanying drawings.

[0057] As Figure 1 shown, a spatio-temporal joint correlation single-photon counting three-dimensional imaging method uses a long-distance photon counting lidar system to collect a set of point cloud data; according to the collected point cloud data, preprocess the image data, apply a windowing method, set the signal photon flight time gating to filter a large amount of noise information; according to the preprocessed data, adaptively obtain the photon flight time of each pixel point to further filter noise; according to the obtained photon flight time, complete the empty pixels according to the neighborhood information; construct a cost function by adding a total variation regularization term to the obtained depth information, and use the alternating direction method of multipliers to iteratively update to obtain the optimal solution, and this solution is used as the final depth information estimation value. The specific steps are as follows:

[0058] Step 1: Set the pixel size N Figure 2 ×N r ×N c of the collected image, parameters such as the integration time of each pixel through the lidar system schematic diagram as r ×N c ×T shown, collect the point cloud data required for the experiment N

[0059] Step 2: Based on the point cloud data obtained after preliminary denoising in Step 1, perform the next processing, in combination with the flowchart in Figure 4. In the case of long-distance detection, the amount of noise photons received by the system is much larger than the amount of signal photons. However, signal photons tend to cluster together more than noise photons. Thus, a method for determining signal photons is to find the largest cluster of photons. Process the data pixel by pixel. At pixel point (i,j), first sort the time of flight (ToF) of the photons, and form photon units with n photons. Assume the number of photons in this pixel is L. Sort the differences between the ToF of adjacent photons in the n-neighborhood, calculate the deviation Δi between adjacent photon units, and find the minimum value Δi of Δi min of the corresponding photon unit:

[0060]

[0061] where is the time of flight of each photon unit at pixel point (i,j), and the ToF of the unit at position l corresponding to the deviation Δi between adjacent photon units is used as the initial estimated depth t of the target depth min : initial :

[0062]

[0063] For the time of flight of the photons, determine whether the next group of photons is a signal photon. If the difference between the time of flight of the next group of signal photons and the initial estimated depth is less than the threshold T p :

[0064]

[0065] then this type of photon unit is called a signal photon unit, and the time of flight of the photons that meet the conditions is recorded in the set . Select the threshold K. If the number of photons in the set at pixel point (i,j) reaches the calibration value K, or the difference between the time of flight of the photons and the initial depth information is greater than the threshold T p range, at this time the determination of this pixel point ends, and the photon determination of the next pixel point is entered.

[0066] Step 3: In the case of low integration time, the number of photons acquired by the detector decreases. When using a very low acquisition time, the number of photons received by the sensor drops sharply. Under these conditions, many empty pixels are generated. Without additional information, it is difficult to perform depth estimation. According to the data obtained in Step 2, due to the presence of missing pixels, based on spatial correlation, the neighborhood information of pixel points is used to fill the empty pixels. For a pixel (i,j), if the number of non-empty pixels marked among the eight nearest surrounding pixels (i-1,j-1), (i-1,j), (i-1,j+1), (i,j-1), (i,j+1), (i+1,j-1), (i+1,j), (i+1,j+1) is not less than 3, it is denoted as a small empty pixel; otherwise, it is denoted as a large hole pixel. Different interpolation templates are matched for different empty pixels for filling. Small pixel points are filled with the mean value of 8-neighborhood non-empty pixels, and large empty pixel points are regarded as empty pixels.

[0067]

[0068] Step 4: Based on the ToF data obtained in Step 3, perform iterative update of the depth information of the image, combined with Figure 5 the flowchart of, and through the alternating direction multiplier method, obtain relatively accurate depth information of the target.

[0069] LiDAR observation model y i,j,t , where (i,j) ∈ {1,…,N r} × {1,…,N c}, representing the photon count at pixel (i,j) within the time period t, following a Poisson distribution:

[0070] y i,j,t ~Poisson(s i,j,t )

[0071] where s i,j,t =r i,j g0(t - t i,j ) + b i,j + d i,j , t i,j ≥0 represents the distance between the detector and the target surface, representing the depth information, r i,j ≥0 represents the intensity of the target, b i,j ≥0 represents a constant, representing the background and dark photon energy levels, and g0 represents the impulse response of the system. Assuming the independence between observed pixels leads to the joint likelihood:

[0072]

[0073] where t, r, b are N×1 vectors, where N = N r N c , and T is the time bin.

[0074] Obtained by minimizing the negative log-likelihood function:

[0075]

[0076] where const represents a constant that does not affect the estimation of the cost function. For simplicity of the cost function, the classical estimation method assumes no background noise (b i,j = 0), and assumes that the instrument impulse response is approximately Gaussian-distributed: where c1 is a constant, σ is a hyperparameter of the system, and is a constant, obtained according to MLE, and Based on this setting, the cost function is (after removing unnecessary constants):

[0077]

[0078] where, represents a constant term that does not affect the estimation of the cost function.

[0079] In the case of extremely low photons, since the maximum likelihood estimation method estimates and with poor quality, the estimated image is thus subjected to regularized image restoration processing, which is explained by considering the following cost function:

[0080] C(t,r) = L(t,r) + φ(t,r)

[0081] where φ(t,r) is the regularization term. Further, the total variation regularization term is adopted to improve the image estimation. According to the neighborhood spatial correlation, the following cost function can be obtained:

[0082] C TV (t,r) = L(t,r) + τ1TV(t) + τ2TV(r)

[0083] where τ1 and τ2 are regularization parameters, where and are the first-order differences of the abscissa and ordinate at pixel n.

[0084] The cost function is minimized by the alternating direction method of multipliers. The convex optimization problem of the alternating multiplier method is:

[0085]

[0086] where a parameter z is set for the joint estimation of depth and intensity, H (j)It is a selection matrix used to select data on depth and intensity in joint estimation. Assume a full-rank matrix:

[0087]

[0088] where j = 1, 2, 3, 4, the cost function can be split to obtain:

[0089] H (1) = [0 N , K]

[0090] H (2) = [K, 0 N

[0091] g3(u (3) ) = τ1TV(z1) + τ2TV(z2), H (3) = I 2N

[0092] H (4) = [I N , 0 N

[0093] where g1, g2, g3, g4 are respectively the intensity, depth, TV regularization, and unnecessary constant terms corresponding to the images after splitting the cost function, where u (3) , u (4) are respectively the solutions for each update, r i ML , are obtained by the maximum likelihood estimation method, I N , K N is the identity matrix. The alternating direction method of multipliers is used to iteratively update each term:

[0094]

[0095]

[0096]

[0097] where s is the variable value for which the minimum is obtained, and k is the number of iterative updates. It is equal to when the updated is less than the threshold Threshold, the iterative loop is exited, and the depth estimate of the target is obtained.​​

Claims

1. A spatio-temporal joint correlation single-photon counting three-dimensional imaging method, characterized in that, The steps are as follows: Step 1: Use a long-distance photon-counting lidar system to collect point cloud data and preprocess the image data; Step 2: For the preprocessed data, adaptively obtain the photon flight time of each pixel point to further filter out noise; Step 3: Complete the empty pixels according to the neighborhood information; Step 4: Add the depth information obtained in Step 3 to the total variation regularization term to construct a cost function, and use the alternating direction method of multipliers to iteratively update to obtain the optimal solution as the depth information estimate value.

2. The spatio-temporal joint correlation single-photon counting three-dimensional imaging method according to claim 1, wherein Preprocessing the image data includes: setting the gating threshold according to prior information and filtering out obvious noise points in the image.

3. The spatio-temporal joint correlation single-photon counting three-dimensional imaging method according to claim 1, characterized in that, The specific method for adaptively obtaining the photon flight time of each pixel point for the preprocessed data and further filtering out noise is as follows: At the pixel point (i, j), sort the flight times of photons, group L photons to form a photon unit, perform a difference sorting on the flight times of adjacent photons in the n-neighborhood, calculate the deviation Δi between adjacent photon units, and find the photon unit corresponding to the minimum deviation Δi min when: In the formula, is the flight time of each photon unit at the pixel point (i, j); From the deviation Δi of adjacent photon units min The flight time of the unit at the corresponding position l is used as the initial estimated depth t of the target depth initial : For the photon flight time, determine whether the next group of photons is signal photons. If the difference between the flight time of the next group of signal photons and the initial estimated depth is less than the threshold T p : Then this type of photon unit is called a signal photon unit, and the photon flight time that meets the conditions is recorded in the set . Select a threshold K. If the number of photons in the set at pixel point (i, j) reaches the calibration value K, or the difference between the photon flight time and the initial depth information is greater than the threshold T p range, at this time the determination of this pixel point ends, and the photon determination of the next pixel point enters.

4. The spatio-temporal joint correlation single-photon counting three-dimensional imaging method according to claim 1, wherein The specific method for completing the empty pixels according to the neighborhood information is as follows: Process the image pixel by pixel. At the pixel point (i, j), fill the empty pixel point with the mean value according to the neighborhood information of its surrounding 8 points.

5. The spatio-temporal joint correlation single-photon counting three-dimensional imaging method according to claim 1 or 4, characterized in that, If for the pixel (i, j), the number of marked non-empty pixel points among its 8 nearest surrounding pixel points (i - 1, j - 1), (i - 1, j), (i - 1, j + 1), (i, j - 1), (i, j + 1), (i + 1, j - 1), (i + 1, j), (i + 1, j + 1) is not less than 3, it is denoted as a small empty pixel point; otherwise, it is denoted as a large hole point. The small pixel points are filled with the mean value of the 8-neighborhood non-empty pixels, and the large empty pixel points are regarded as empty pixel points.

6. The spatio-temporal joint correlation single-photon counting three-dimensional imaging method according to claim 1, wherein The specific method for using the alternating direction method of multipliers to iteratively update to obtain the optimal solution is as follows: At the pixel (i, j), the photon count follows a Poisson distribution, and a joint likelihood estimate is performed: where s i,j,t = r i,j g0(t - t i,j ) + b i,j + d i,j , s i,j,t represents all the photon signals received by the detector, y i,j,t is the lidar observation model, t i,j ≥ 0 represents the distance between the detector and the target surface, i.e., the depth information, r i,j ≥ 0 represents the intensity of the target, b i,j ≥ 0 represents the background and dark photon energy levels, g0 represents the impulse response of the system, t, r, b are N×1 vectors, where N = N c × N r , T is the time bin, and the cost function is obtained by minimizing the negative log-likelihood function: where const represents a constant that does not affect the cost function estimation, assuming that the background and dark photon energy level b i,j = 0, and the instrument impulse response is approximately Gaussian distributed where c1 is a constant, σ is a system hyperparameter, and is a constant, and the depth and intensity estimates of the image are obtained according to the maximum likelihood estimation method and the cost function is simplified to: Among them, represents a constant term that does not affect the cost function estimation. The total variation regularization term is used to improve image estimation. According to the neighborhood spatial correlation, the cost function C after total variation can be obtained TV (t, r): C TV (t, r) = L(t, r) + τ1TV(t) + τ2TV(r) where τ1 and τ2 are regularization parameters, where and are the first-order differences of the abscissa and ordinate at pixel n. To minimize the cost function, the cost function is minimized by the alternating direction method of multipliers. The convex optimization problem of the alternating multiplier method is: Among them, the joint estimation of depth and intensity is set as a parameter H( j ) is a selection matrix used to select the data of depth and intensity in the joint estimation; Here, assume a full-rank matrix M as: where j = 1, 2, 3, 4, and the cost function is split to obtain: where \(g_1\), \(g_2\), \(g_3\), \(g_4\) are the intensity, depth, TV regularization, and unnecessary constant term of the corresponding images after the cost function is split, respectively. Among them u( 3 ), u( 4 ) are the solutions for each update, r i ML , are calculated by the maximum likelihood estimation method, I N , K N is the identity matrix, and the alternating direction method of multipliers is used to iteratively update each term: Among them, s is the variable value for which the minimum value is obtained, and k is the number of iterative updates. That is, when the updated is less than the threshold Threshold, the iterative loop is exited to obtain the depth estimation of the target.

Images