Phase recovery method based on implicit regularization and BDCA algorithm

Abstract

The present invention discloses a phase recovery method based on implicit regularization and algorithm BDCA. According to the algorithm BDCA, a minimal solution of a phase recovery model based on implicit regularization is solved, and a phase-recovered image u is output. The calculation formula of the phase recovery model based on implicit regularization is as follows: #imgabs0#wherein: u is a clean image, u∈[0,1], λ is a regularization parameter, b is an observed degraded image, #imgabs1# is a Fourier transform operator, #imgabs2# represents a complex domain space with dimensions n and m respectively, control parameters ∈>0, <·,·> represents an inner product; #imgabs3# is a Huber regularization function, #imgabs4# is a linear operator, and the Huber regularization function is #imgabs5#. The present invention effectively performs phase recovery when the image has high noise.

Description

Technical Field

[0001] The invention relates to a phase recovery method based on implicit regularization and algorithm BDCA, and belongs to the technical field of computational imaging. Background Art

[0002] In optics, the information of high-frequency light waves can be easily recorded by optical devices. However, it is quite difficult to obtain the phase information of these lights because the diffraction plate can only capture the absolute value of the Fourier transform sample. In addition, during the acquisition process, the amplitude information obtained is always affected by severe noise, which makes the problem more complicated. Since the phase of the image contains rich information, it is necessary to reconstruct it.

[0003] Assume that the image is degraded by Gaussian white noise ζ with mean 0, The general phase recovery problem can be expressed as:

[0004]

[0005] Among them: unknown clean image is the image we observed, is a linear operator in complex Euclidean space, Ω is a discrete lattice, and n = n1 × n2 is the image size. Model (1) is obviously a non-convex and non-smooth optimization problem. We can clearly observe that if there is no prior information, the model does not have a unique solution.

[0006] The focus of phase recovery is to consider how to recover images or signals from Fourier transform amplitudes with noise. In order to solve this problem, existing phase recovery algorithms often use prior knowledge to recover images. For example, Gerchberg and Saxton first proposed the error-reduction algorithm, which reconstructs the phase through direct transformation between Fourier and real domains. Candès, Li, and Soltanolkotabi proposed a Wirtinger flow algorithm to solve the phase recovery problem of non-convex functions through theoretical analysis. Wang et al. proposed a method with a computational complexity of k. 2 The nlogn sparse truncated amplitude flow algorithm is used to reconstruct an n-dimensional sparse k signal (k<<n). Vaswani, Nayer and Eldar use truncated Wirtinger flow and alternating minimization techniques to deal with low-rank phase recovery problems. There are many other phase recovery methods, such as phaselift method, differential mapping iterative method, hybrid input-output algorithm, etc.

[0007] Although the above method can reconstruct the phase through the amplitude information of Fourier transform, it may not be applicable to images with high noise. Moreover, some existing algorithms will make the non-convex model converge to some local minimum values. Therefore, we need to consider finding the global optimal solution of the phase recovery task without changing the non-convex nature of the non-convex problem. In this case, it is necessary to design a fast, effective and convergent algorithm for non-convex problems. Summary of the invention

[0008] Objective: To overcome the deficiencies in the prior art, the present invention provides a phase recovery method based on implicit regularization and algorithm BDCA.

[0009] Technical solution: To solve the above technical problems, the technical solution adopted by the present invention is:

[0010] A phase recovery method based on implicit regularization and BDCA algorithm includes the following steps:

[0011] According to the algorithm BDCA, the minimal solution of the phase recovery model based on implicit regularization is solved and the phase recovered image u is output.

[0012] The calculation formula of the phase recovery model based on implicit regularization is as follows:

[0013]

[0014] Where: u is the clean image, u∈[0,1], λ is a regularization parameter, b is the observed degraded image, is the Fourier transform operator, They represent the complex domain space with dimensions n and m respectively, the control parameter ∈>0, and <·,·> represents the inner product. is a continuously differentiable Huber canonical function, is a linear operator, and the specific Huber canonical function is

[0015] As a preferred solution, the minimal solution of the phase recovery model based on implicit regularization is solved according to the BDCA algorithm, and the phase-recovered image u is output, including the following steps:

[0016] Step 1: Initialize parameter u 0 , v 0 =u 0 , Dual variable initialization parameter q 0 =0,w 0 =0,d 0= 0, backward search step parameters α,ξ,r1,r2,r3,λ,β,ρ,∈, set the maximum number of inner iterations Intmax, the maximum number of outer iterations Outmax;

[0017] Step 2: When k = 1:Outmax, calculate h using equation (3) k ;

[0018] Step 3: When m = 1:Intmax, update v by equations (8), (10), (15)-(17), and (20) respectively. m+1 ,p m+1 ,z m +1 ,u m+1 Subproblem, get v m+1 ,p m+1 ,z m+1 ,u m+1 ,u m ;

[0019] Update the multiplier d m+1 =d m +r1(z m+1 -Au m+1 );

[0020] w m+1 =w m +r2(u m+1 -v m+1 );

[0021]

[0022] when When the threshold is reached, the inner iteration is terminated, and y k =u m+1 ,u k =u m ;

[0023] Step 4: Set the backward search direction d k =y k -u k ;

[0024] Step 5: When d k = 0, terminate the external iteration, let u k+1 =u k , output u k+1 When d k ≠0, go to step 6;

[0025] Step 6: Select any parameters And make When F(y k +β k dk )>F(y k )-αβ k 2 ||d k || 2 Execute β when the condition is met k =ξβ k , F(*) represents formula (1).

[0026] Step 7: Let u k+1 =y k +β k d k , when u k+1 -u k = 0, terminate the external iteration and output u k+1 ; when u k+1 -u k ≠0, go to step 2.

[0027] Output the iteration result u from step 5 or step 7 k+1 Then it is the final phase recovery result u.

[0028] As a preferred solution: the calculation formula of formula (3) is as follows:

[0029]

[0030] in, represents the subdifferential, A T represents the conjugate transpose of A, h k ,u k is the estimated value of the kth iteration, and ρ is a positive parameter that controls the degree of convexity of the model.

[0031] As a preferred solution, the calculation formula of formula (8) is as follows:

[0032]

[0033] Among them, u m 、w m Represent the estimated values ​​of u, w at the mth iteration, v m+1 represents the estimated value of v at the m+1th iteration. r2 is the Lagrangian parameter.

[0034] As a preferred solution, the calculation formula of formula (10) is as follows:

[0035]

[0036] represents the estimated value of the m+1th iteration, and r3 is the Lagrangian parameter.

[0037] As a preferred solution, the calculation formulas of formulas (15), (16) and (17) are as follows:

[0038]

[0039]

[0040]

[0041] represents the estimated value of the m+1th iteration, δ m represents the estimated value of the mth iteration, r1 is the Lagrange parameter, and sign(*) is the sign function.

[0042] As a preferred solution, the calculation formula of formula (20) is as follows:

[0043]

[0044] in, represent The complex conjugate of is a linear operator, represent The conjugate transpose of represents the real part, represents the imaginary part, The prime numbers of integers s1 and s2 are n1 and n2, and t1 and t2 are the positions of the 2D signals. Is an imaginary unit.

[0045] As a preferred solution, the threshold is set to 10 -5 .

[0046] Beneficial effects: The present invention discloses a phase recovery method based on implicit regularization and algorithm BDCA, and proposes a phase recovery model based on implicit regularization, which can effectively perform phase recovery when the image has high noise. Due to the non-convex nature of phase recovery, we consider the non-convex minimization model of this task in the solution process. And considering the adverse effects of noise in the phase recovery task, we use regularization to smooth the noise, so an implicit regularization term is constructed. According to the decomposability of the phase recovery model, we use the algorithm BDCA with convergence guarantee to effectively solve the model. Secondly, in order to verify the robustness and flexibility of the model, two simple and effective total variation regularization and tight frame regularization are used to smooth the noise. Numerical experimental results show that the model we proposed is flexible and efficient in performing phase recovery tasks for degraded images with high noise. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 is the phase recovery result when the noise intensity σ=0.

[0048] Figure 2 is the phase recovery result when the noise intensity σ=10.

[0049] Figure 3 This is the phase recovery result when the noise intensity is σ=30.

[0050] Figure 4 It is the energy functional curve of the phase recovery model based on implicit regularization. DETAILED DESCRIPTION

[0051] The present invention will be further described below in conjunction with specific embodiments.

[0052] A phase recovery method based on implicit regularization and BDCA algorithm includes the following steps:

[0053] Step 1: Establish a phase retrieval model based on implicit regularization

[0054] Assuming that the clean image u is a non-negative bounded real value, u∈[0,1], in order to remove the noise in the degraded image b as much as possible, we use regularization to smooth b, then the phase recovery model based on implicit regularization is:

[0055]

[0056] Where: λ is a positive parameter and b is the observed degraded image. is the Fourier transform operator. is an implicit regularization that is highly effective in removing noise. It can be a gradient operator Wavelet Operator Low-rank operators and so on.

[0057] Huber canonical function Φ ∈ (u) The form is as follows

[0058]

[0059] where ∈>0, and is continuously differentiable.

[0060] Step 2: Use the BDCA algorithm to solve the minimal solution of phase recovery based on implicit regularization

[0061] According to the decomposable characteristics of the model, the phase recovery model based on implicit regularization is F(u)=G(u)-H(u):

[0062]

[0063] Among them, the coefficient ρ is a positive parameter that controls the degree of strong convexity of the model, and G(u) and H(u) are both strongly convex terms.

[0064] Using the BDCA algorithm, we obtain the following formula

[0065]

[0066] in represents the subdifferential, A T represents the conjugate transpose of A, h k ,u k is the estimated value of the kth iteration, u k+1 is the estimated value of the k+1th iteration, and ρ is a positive parameter that controls the degree of convexity of the model.

[0067] Considering u k+1 The difficulty of solving the subproblem is that we use the splitting algorithm ADMM to solve the subproblem. After introducing the auxiliary variable, u k+1 The sub-problem is written as the following constraint form expressed as formula (5):

[0068]

[0069] Where z = (z0, z1, z2), BOX constraint p=(p i ,p i+n ) T , p i ,p i+n represent The components of are in two different directions, similar to the components of the gradient operator in the horizontal and vertical directions. n is the number of pixels in image u, and z0, z1, z2 represent the three channels of z. Define b: = (b0, b1, b2), where b0, b1, b2 are the three channels of image b.

[0070] Formula (5) can be expressed as Lagrangian function as formula (6):

[0071]

[0072] Among them, w, q, d are Lagrange multipliers, r1, r2, r3 are Lagrange parameters, and d = (d0, d1, d2) is defined, where d0, d1, d2 are the three channels of parameter d.

[0073] According to the ADMM algorithm, the variables of formula (6) are split and each sub-problem w, p, z, u is solved separately.

[0074] (1) The v subproblem can be expressed as

[0075]

[0076] The closed-form solution can be obtained as

[0077]

[0078] v m+1 Represents the estimated value of v at the m+1th iteration. m 、w m They represent the estimated values ​​of u and w at the mth iteration respectively.

[0079] (2) The p-subproblem can be expressed as

[0080]

[0081] Since the above formula is splittable, it can be written as (p i ,p i+n ), we define So the solution is

[0082]

[0083] Represents the estimated value of the m+1th iteration.

[0084] (3) To solve the z subproblem, we need to solve each of its components separately. Let’s take solving z0 as an example:

[0085]

[0086] Decomposing the minimization problem

[0087]

[0088] There are two sub-problems, namely |z0| and sign(z0).

[0089] When z0≠0 Otherwise sign(0) = c, is an arbitrary constant of unit length. To minimize |z0|, we have

[0090]

[0091] The closed-form solution is

[0092]

[0093] make Then the closed-form solution can be reconstructed as

[0094]

[0095] represents the estimated value of the m+1th iteration, δ m Represents the estimated value at the mth iteration.

[0096] Similarly, z1 and z2 are solved according to this process. Therefore, we get

[0097]

[0098]

[0099] Represents the estimated value of the m+1th iteration, and sign(*) is the sign function. .

[0100] (4) Solve the u subproblem

[0101]

[0102] By calculating ||z-Au|| 2 / 2, we get

[0103]

[0104] stands for Fourier transform, stands for inverse Fourier transform, represents the complex conjugate of z, represents the real part, represents the imaginary part, where the amplitude modulation

[0105] 0≤t j ≤n j -1,j=1,2, the prime numbers of integers s1,s2 are n1,n2, t1,t2 are the positions of the 2D signal, Is an imaginary unit.

[0106]

[0107] in According to the above step-by-step solution process, we summarize the process as Algorithm 1.

[0108]

[0109]

[0110] Iterating the regular framework model (1) according to the iteration steps in Algorithm 1 will eventually output the iteration result in step 5 or step 7, u k+1 Then it is the final phase recovery result u.

[0111] To verify the effectiveness of the phase recovery algorithm based on implicit regularization and BDCA in this chapter, we selected the images in the test set Set18 as the test images of the phase recovery algorithm. And we consider the phase recovery model based on implicit regularization in the operator is the gradient operator and wavelet operator The recovery results in two cases. Operator uses the gradient operator When replacing, the model is recorded as TV model; the implicit regularization Wavelet operator When replacing, the model is recorded as TF model.

[0112] In practical applications, the image amplitude observation |Au| may contain different levels of noise due to interference from uncontrollable factors. In order to test the stability of the phase recovery algorithm under different levels of noise, we use the parameter σ to control the noise level in the simulation experiment. First, in the case of no noise σ = 0, we test the stability of our phase recovery algorithm. Figure 1 As shown in the figure, (a) is a noise-free image, and (b) is the amplitude information of image (a). We can see that (b) lacks the structural information and contour information of the target in (a), so it is difficult to perform the phase recovery task. According to the experimental results of the TV model and TF model in (c) and (d), the model based on implicit regularization has better phase recovery results, overcomes the step effect of total variation regularization, and utilizes the sparsity advantage brought by the wavelet operator.

[0113] In order to demonstrate the denoising ability of the phase recovery method based on implicit regularization and algorithm BDCA, this paper selected two test images for experiments under the noise intensity of Gaussian white noise σ=10, 20, and the size of the two images was fixed to 500×500. The experimental results are shown in Figure 2 and Figure 3 As shown in Figure 2, the phase recovery algorithm in this chapter can maintain image texture and details and has good denoising ability. Under different noise intensities, the restoration results of our method are closer to the original image and the restored details are richer.

[0114] For non-convex functions, we demonstrate the convergence of the algorithm through experiments. Figure 4 We found that the energy functional of the objective equation (1) is monotonically decreasing. The curve decreases rapidly without violent oscillations. The experiment converges after about 20 iterations. These experimental results verify that our algorithm has good convergence and reflect the robustness of the strategy.

[0115] This chapter proposes a phase recovery method based on implicit regularization and the BDCA algorithm, using the prior information of the image itself as an implicit regularization term. According to the splitting characteristics of the non-convex model, this chapter uses the composite splitting algorithm BDCA+ADMM for fast solution. The original problem is decomposed into multiple easy-to-solve sub-problems and solved independently, which speeds up the convergence of the algorithm. Experimental results show that our phase recovery algorithm is robust to degraded images with different noise levels, flexibly handles phase recovery tasks with different regularizations, and has good convergence performance.

[0116] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principle of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A phase recovery method based on implicit regularization and BDCA algorithm, characterized by: The steps include: According to the algorithm BDCA, the minimal solution of the phase recovery model based on implicit regularization is solved, and the phase recovered image u is output; The calculation formula of the phase recovery model based on implicit regularization is as follows: Where: u is the clean image, u∈[0,1], λ is a regularization parameter, b is the observed degraded image, is the Fourier transform operator, Represent the complex domain space with dimensions n and m respectively, control parameter ∈>0, and <·,·> represents the inner product; is the Huber regularization function, is a linear operator, and the Huber canonical function is According to the algorithm BDCA, the minimal solution of the phase recovery model based on implicit regularization is solved and the phase recovery image u is output, which includes the following steps: Step 1: Initialize parameters Dual variable initialization parameter q 0 =0,w 0 =0,d 0 = 0, backward search step parameters α,ξ,r1,r2,r3,λ,β,ρ,∈, set the maximum number of inner iterations Intmax, the maximum number of outer iterations Outmax; Step 2: When k = 1:Outmax, calculate h using equation (3) k ; Step 3: When m = 1:Intmax, update v by equations (8), (10), (15)-(17), and (20) respectively. m+1 ,p m+1 ,z m+1 ,u m +1 Subproblem, get v m+1 ,p m+1 ,z m+1 ,u m+1 ,u m ; Update the multiplier d m+1 =d m +r1(z m+1 -Au m+1 ); In m+1 =in m +r2(u m+1 -v m+1 ); q m+1 =q m +r3(p m+1 -Bu m+1 ); when When , the inner iteration terminates, and let y k =u m+1 ,u k =u m ; Step 4: Set the backward search direction d k =y k -u k ; Step 5: When d k = 0, terminate the external iteration, let u k+1 =u k , output u k+1 When d k ≠0, go to step 6; Step 6: Select any parameters And make When F(y k +β k d k )>F(y k )-αβ k 2 ||d k || 2 Execute β when the condition is met k =ξβ k , F(*) represents formula (1); Step 7: Let u k+1 =y k +β k d k , when u k+1 -u k = 0, terminate the external iteration and output u k+1 ; when u k+1 -u k ≠0, go to step 2; Output the iteration result u from step 5 or step 7 k+1 Then it is the final phase recovery result u; The calculation formula of formula (3) is as follows: in, represents the subdifferential, A T represents the conjugate transpose of A, h k ,u k is the estimated value of the kth iteration, ρ is a positive parameter that controls the degree of strong convexity of the model; The calculation formula of formula (8) is as follows: Among them, u m 、w m Represent the estimated values ​​of u, w at the mth iteration, v m+1 represents the estimated value of v at the m+1th iteration; r2 is the Lagrangian parameter; The calculation formula of formula (10) is as follows: represents the estimated value of the m+1th iteration, r3 is the Lagrangian parameter; The calculation formulas of formulas (15), (16) and (17) are as follows: represents the estimated value of the m+1th iteration, δ m represents the estimated value of the mth iteration, r1 is the Lagrange parameter, sign(*) is the sign function; The calculation formula of formula (20) is as follows: in, represent The complex conjugate of , B is a linear operator, represent The conjugate transpose of represents the real part, represents the imaginary part, The prime numbers of integers s1 and s2 are n1 and n2, and t1 and t2 are the positions of the 2D signals. Is an imaginary unit.

2. The phase recovery method based on implicit regularization and BDCA algorithm according to claim 1, characterized in that: Said is continuously differentiable.

Images