An improved algorithm for diffractive optical element with high imaging quality

. An improved algorithm for diffractive optical element (DOE) with high imaging quality is proposed in this paper. The algorithm is designed based on amplitude division between signal and noise regions, further subdivides the noise region into two distinct parts. The image quality in the signal region will be effectively improved by employing a partition-constraint strategy, which imposes amplitude freedom on the ﬁ rst noise region while enforcing strict amplitude constraints on the second noise region. The principle of the algorithm, simulation analysis, and experimental results are presented. The simulation and experimental results demonstrate that the algorithm is feasible.


Introduction
Diffractive optical element (DOE) is widely used in beam shaping, optical encryption, optical communication, optical imaging, and other fields for its advantages of small size and flexible design [1][2][3][4][5].Research on design methods of DOE has academic and application value.Scholars have conducted extensive studies and achieved some research findings, such as Gerchberg-Saxton (GS) algorithm, inputoutput algorithm, simulated annealing algorithm, Yang-Gu algorithm, genetic algorithm, and other design methods [6][7][8][9][10][11].The GS algorithm plays a significant role in DOE design because of its features of fast optimization speed and high diffraction efficiency.However, it is well-known that the GS algorithm is limited by the selection of initial phase and tends to fall into local optimal solutions.At the same time, it is quite difficult to optimize the speckle noise by GS algorithm.To address these issues, scholars have done a lot of improvement and optimization researches, which can be mainly summarized into two categories: initial phase optimization and changing amplitude constraint methods.
In terms of initial phase optimization, Aagedal et al. proposed the use of pseudo-random or spherical phases as initial phase to reduce speckles [12].Pang et al. tried to control the quadratic phase with different parameters, eventually successfully shaping the Gaussian beam into a high-quality flat-top beam [13].Fischer and Sinzinger calculated the iterative initial phase by combining quadratic phase and error diffusion algorithm, and obtained a reconstructed image with good quality but slightly lower diffraction efficiency [14].Doskolovich et al. solved the initial phase distribution of specific patterns based on geometric methods, and obtained a piecewise smooth phase distribution through GS algorithm [15].For patterns with a uniformly distributed intensity, optimizing with an appropriate quadratic phase has obvious advantages over using a random phase.However, for patterns with more complex intensity distributions, it is invalid to use a quadratic phase as the initial value for obtain diffractive image with no artifacts.Instead, more complex methods, such as ray tracing and error diffusion, need to be employed to calculate a suitable initial phase.
Changing amplitude constraint method divides the target image into signal and noise regions.The zoning constraints strategy is considered to be a more universal GS optimization algorithm, and scholars have done a lot of research [16][17][18][19].The research results show that imposing different amplitude constraints on these two different regions can effectively reduce speckle noise in the signal window, thereby it will improve the quality of the reconstructed image.However, these optimization algorithms do not impose constraints on the noise region, using a random phase as the initial phase would result in a significant amount of speckle noise distributed in the noise region.In addition, the optical energy utilization will decrease with the enlargement of noise region.Therefore, these algorithms usually need to use quadratic phase as the initial phase [18].
However, the use of quadratic phase will introduce ringing and artifacts while suppressing speckle noise [20].
An improved algorithm based on GS is proposed in this paper, which can reduce speckle noise with no ringing and artifact.The algorithm can be illustrated as follow: Firstly, the target image plane is divided into three regions: signal region, first noise region, and second noise region.Secondly, imposes amplitude freedom on the first noise region while enforcing strict amplitude constraints on the second noise region.Thirdly, the amplitude in signal window is obtained by linearly superimposing the amplitude of the diffracted image plane and the amplitude of the target image.Finally, the optimal phase distribution is obtained while image quality matches three evaluation functions during the iteration process.During the iteration process, speckle noise is concentrated in the unconstrained first noise area, thereby improving the image quality within the signal area.Meanwhile, amplitude constraints in the second noise area reduce the intensity distribution of speckle noise within that area, effectively addressing the issue of reduced diffraction efficiency that arises as the zero-padding area increases.The relationship between image quality and design parameters are investigated by numerical simulation, which can give significant design recommendations for DOE with high quality.Furthermore, the DOE calculated by the proposed algorithm is fabricated by using our self-made Micro-Nano photolithography system.The high-quality diffractive image with no ringing and artifacts is obtained by optical reconstruction.

Method
As shown in Figure 1, the target image is divided into three regions.The central region is the signal region (S), the checkered region outside the signal region is the first noise region (N 1 ), and the remaining all-black region is the second noise region (N 2 ).
According to the Nyquist sampling theorem [21], calculate the minimum required number of sampling points.The target image is zero-padded to M Â N, serving as the actual target light intensity distribution I t during optimization.Based on this, the target amplitude distribution A t is calculated from I t .The formula for calculating the number of sampling points is given below: where k is the operating wavelength of DOE, d is the distance between the imaging plane and the DOE, pix h is the sampling interval on the DOE plane, pix o is the sampling interval on the diffraction image plane.Before starting iteration, a random phase u 0 is multiplied to A t and propagated to the DOE plane by an inverse Fourier transform.The phase u k ð Þ h on the DOE plane is taken as the initial phase for iterative optimization.The flowchart of the algorithm is shown in Figure 2.
The specific steps are as follows: The complex amplitude distribution on the back surface of the DOE after the incident plane wave passes through the DOE phase modulation is given by: The symbol "k" in the formula represents the iteration number.In the first iteration, the initial phase is set to u Then, the light field emitted from the DOE will propagate to the image plane by Fourier transform, the complex amplitude distribution on the image plane is as follows: The intensity distribution on the image plane can be calculated by: where Ã is the conjugate of the diffracted light field distribution.The evaluation of the diffraction light field distribution is conducted using three assessment functions: diffraction efficiency (g), structural similarity index measure (SSIM), and peak signal-to-noise ratio (PSNR).The expressions for these functions are as expressed by: The parameters I o and I t respectively represent the intensity distribution on the diffraction image plane and the target intensity distribution.When the sampling points are reduced to the size of the signal window, the calculated results represent the SSIM and PSNR within the signal window.In order to obtain a high-quality reconstructed image, amplitude constraints are applied to the partitioned blocks on the image plane.The amplitude in the signal region is constrained by linearly superimposing the amplitude of the diffracted image plane and the amplitude of the target image.The amplitude in the first noise region remains unchanged, and the amplitude in the second noise region is constrained to 0: o ðu; vÞ is the corrected amplitude distribution on the diffracted plane, a is the parameter used to adjust the weight.The complex amplitude of the corrected diffracted image plane is: After the amplitude constraint, the complex amplitude of the light field on the DOE plane is obtained by using inverse Fourier transform on the complex amplitude of the corrected diffracted image plane: Then, the phase information on the DOE plane is retained and multiplied by the input plane wave light field to obtain the corrected complex amplitude of the light field on the DOE plane: Finally, return to the previous propagation process and perform the next iteration until the set number of iterations is reached or the evaluation function meets the requirements, to eventually obtain the required phase distribution on the DOE plane.

Simulation analysis
The diffraction theory indicates that larger diffraction angle corresponds to higher spatial frequency.When the first noise area is enlarged, the diffraction angle of the speckle noise distributed in this area also increases, inevitably leading to an increase in high-frequency information on the DOE, which will increase the difficulty of fabricate DOE If the first noise area is too small, it offers limited amplitude freedom, and the improvement in the quality of the diffraction image is not significant.Therefore, it is necessary to find an appropriate size for the first noise area that can both enhance the quality of the reconstructed image and not overly increase the complexity of DOE fabrication.Additionally, the weight value a, as a crucial parameter controlling the constraints within the signal window, also significantly affects the image quality within the signal window.The relationship between reconstructed image and parameter values have been explored in this section.Firstly, to find an appropriate weight value a, we use the images "Baboon", "Barbara", "Cameraman", "Goldhill" and "Peppers" as target images in the simulation.The signal window is set to 512 Â 512 pixels, and the image is zeropadded to a size of 2000 Â 2000.During the iterations, the amplitudes outside the signal region are constrained to 0, while keeping the initial constant.The weight value a starts from 0.5 and changes in increments of 0.1, up to 2.5.We compared the PSNR of the diffraction image in the signal window for different weight values and diffraction efficiency g.For each target, 10 different initial phases are used, and the average PSNR and g are calculated.The simulation analysis results are shown in Figure 3. Figure 3a presents five target images for simulation analysis, the relationship between PSNR, g and weight value a are respectively given by Figures 3b and 3c.
From Figure 3b, the PSNR will improve with a increases from 0.5 to 2, which can reach a peak near 2. In addition, the PSNR starts to decrease rapidly as a continues to increase to greater than 2. In Figure 3c, as a gradually increases from 0.5 to 2.1, the diffraction efficiency slowly decreases but still remains above 93%.However, when a exceeds 2.1, the diffraction efficiency experiences a significant decline.Through simulation analysis, we determined that a equal to 2 can close to the best performance.
Secondly, for a reasonable division of the noise region, we designed another simulation using "Baboon" as the object image.The signal window is set to 500 Â 500 pixels, zero-filled to 2000 Â 2000 pixels.While keeping the zerofilling region size constant, the width of the first noise region started from 500 pixels and increased by 50 pixels each time   the entire region is filled.By comparing the diffraction image quality under different window sizes, we aimed to find the appropriate region division.In the simulation, 20 different initial phases are used to calculate the PSNR and SSIM within the signal window under the conditions of the same initial phase but different sizes of the first noise region.The average values are then calculated.The simulation analysis results are shown in Figure 4.
By observing Figure 4, it can be seen that PSNR rapidly increases with the enlargement of the first noise area, reaching its peak when expanded to 2.6 times the length of the signal window area.After that, there is a decrease of 2-13 dB as the first noise area continues to expand.The trend of SSIM is similar to PSNR, with the difference that SSIM converges after Ln/l reaches 2.5.Considering both PSNR and SSIM, when Ln/l = 2.5, which means the length of the first noise region is 2.5 times that of the signal region, it results in the best performance.
Based on the analysis from the two simulations, we have determined that the optimal parameters for the algorithm are a weight value a is equal to 2, and the first noise region should be set to 2.5 times the width of the signal window.To verify the feasibility of the above findings, the Baboon image is used for testing.The simulation platform uses Windows 11 operating system, Core I7-12700 (2.10 GHz), NVIDIA GeForce RTX 3070 Ti.The parameter settings are shown in Table 1.
After zero-padding, the image has dimensions of 2000 Â 2000, with the signal region window size set at 512 Â 512.The size of the first noise region is defined as 2.5 times that of the signal window, resulting in dimensions of 1280 Â 1280, while the remaining area is designated as the second noise region.A comparison is conducted between the performance of the proposed algorithm and GS algorithm.After 100 iterations, the phase distribution is obtained while the PSNR in the signal region reaches the optimal state.Figures 5a and 5b respectively represent the images reconstructed by the proposed algorithm and GS algorithm within the signal window.Comparing Figures 5a and 5b, due to the large amount of speckle noise, the image quality reconstructed by GS algorithm decreased a lot.In contrast, the proposed algorithm can reconstruct the required light field distribution better, significant improvements in contrast and clarity.
A comparison of PSNR and SSIM is illustrated in Figures 6a and 6b.The proposed algorithm converges more rapidly than GS algorithm as observed from Figures 6a and  6b.The proposed algorithm can achieve satisfactory results after approximately 13 iterations and reach convergence after around 20 iterations.However, GS algorithm converges in approximately 30 iterations.The convergence speed of the proposed algorithm is 1.5 times faster than GS algorithm.The evaluation parameters for the reconstructed images are presented in Table 2.
The parameters PSNR, SSIM and g in Table 2 are represent the peak signal-to-noise ratio, structural similarity and diffraction efficiency of the entire diffraction plane respectively.The PSNR S and SSIM S represent the peak signal-to-noise ratio and structural similarity of the reproduced images in the signal region.The parameters PSNRs and SSIMs are increased by 23.094 dB and 0.362, which represents an improvement of over 127% and 57% compared to the GS algorithm, indicating a substantial improvement in image quality.There is a 3% difference of PSNR between the two algorithms.In terms of optical energy utilization, the proposed algorithm exhibits 21.5% decrease compared to GS algorithm.These can be attribu-ted to the presence of some speckle noise in the first noise region.The results show that the proposed algorithm greatly improves the reconstructed image quality in the signal region by sacrificing some intensity in the signal region.

Experiment and results
To further verify the feasibility of the proposed algorithm, we use a self-made high-resolution Micro-Nano photolithography system based on DMD parallel exposure to fabricate the calculated DOE and conduct optical reconstruction.An 8th order quantized DOE can already achieve 94.9% diffraction efficiency.Excessive quantization orders do not contribute to a significantly improvement in diffraction efficiency, but instead increases the difficulty of DOE fabrication.Therefore, in the experiments conducted in this paper, the phase is quantized to an 8th order.Key parameters of the Micro-Nano photolithography system and the photolithography plate used for fabricating DOE are shown in Table 3.
The fabrication process of the DOE is illustrated in Figure 7a.The 8th-order quantized phase distribution of the target image is calculated using the proposed algorithm.obtained phase distribution is then exposed onto a photoresist plate by using our self-made Micro-Nano photolithography system, which is designed based on DMD parallel exposure.The DOE is obtained after the process of exposure and development.Finally, the DOE is directly illuminated by laser beam as Figure 7b.
In the experiment, the characters "DOE" is considered as the image target.The calculation parameters are shown in Table 4. Figures 8a, 8b, and 8c are respectively presented the target image, micrograph of the DOE taken by laser microscope (Keyence VK-X250), and the reconstructed image.In the microscopic image of the DOE in Figure 8b, different colors represent different depths, gradually increasing from red to blue.
As predicted in the previous simulations, the algorithm proposed in this paper produces quite high-quality optical reproductions without ringing and artifacts.Meanwhile, through the amplitude constraint in different regions, the diffraction image in the signal area has high quality and uniformity.With an 8-level phase quantization, the PSNR within the signal area is 17.08 dB, and the SSIM is 0.8031.
In the first noise area, where the amplitude lacks constraints, there is some speckle noise, while in the second noise area, almost no speckle noise can be observed.

Conclusion
This paper introduces an improved algorithm for diffractive optical element based on amplitude subdivision and weighted correction method, which can significantly improve the imaging quality with no ringing and artifacts.These advantages are expected to address the challenges of low imaging quality in DOE, as well as the difficulty in balancing zero-padding size and diffraction efficiency during the phase optimization process.Leveraging the benefits of the proposed algorithm, including the absence of conjugate images and the presence of zero-order light, with its quite high diffraction efficiency, we anticipate that the proposed algorithm can facilitate the practical application of DOE, such as in beam shaping and optical imaging, among others.

Fig. 3 .
Fig. 3.The relationship between image quality and weight parameter.

Fig. 4 .
Fig. 4. The variety rule of image quality with the changing of window sizes.

Fig. 5 .
Fig. 5. Simulation analysis.(a) The reconstructed image within the signal window using the proposed algorithm (b) The reconstructed image within the signal window using GS algorithm.

Fig. 6 .
Fig. 6.(a) PSNR in the signal window of the diffraction image relative to the number of iterations, (b) SSIM in the signal window of the diffraction image relative to the number of iterations.

Table 2 .
Evaluation metrics for optimization results of different algorithms.

Table 3 .
Parameters of the Micro-Nano photolithography system and photolithography plate.