Microlens array spherical aberration correction method based on high-order linear ordinary differential equation

Abstract

The invention relates to a microlens array spherical aberration correction method based on a high-order linear ordinary differential equation. The problem that an existing microlens array spherical aberration correction method is poor in microlens array spherical aberration correction effect is solved. The method comprises the following steps: firstly, establishing a high-order linear ordinary differential equation used for describing a single-lens light propagation path in a micro-lens array, and introducing a non-constant coefficient function; based on a set initial condition and the obtained pixel initial coordinates of the element image array, solving a simultaneous equation of the light field space geometrical relationship and a high-order linear ordinary differential equation by utilizing a formulation reduced-order iterative numerical solution algorithm to obtain an aberration matrix of the element image array; and feeding back the aberration matrix to a spatial function of the element image array mapped to the light field, and changing the pixel horizontal coordinate and vertical coordinate positions of the corresponding element image according to the distance from the pixel point to the optical axis of the corresponding lens to complete the spherical aberration correction of the microlens array. According to the invention, the spherical aberration correction capability of the integrated imaging system can be improved, and the stereo imaging quality is improved.

Description

Technical Field

[0001] The present invention relates to the field of integral imaging technology, and particularly to a method for correcting spherical aberration of a microlens array based on a higher-order linear ordinary differential equation. Background Art

[0002] With the development of science and technology, three-dimensional display technology plays an increasingly important role in many fields, and more fields also hope to apply this three-dimensional scene reproduction technology to practice, especially in industries such as medical treatment, military, and entertainment. Integral imaging is a display technology that can view continuous parallax and full-color three-dimensional stereoscopic images without wearing any auxiliary devices, abbreviated as integral imaging, and has received extensive attention from research scholars in various countries since its proposal. The principle of integral imaging is to use a specific lens array and a display screen to achieve naked-eye stereoscopic display by multi-angle convergent imaging of multiple homologous points. The platform is easy to build and the stereoscopic effect is remarkable.

[0003] Integral imaging mainly has two links: the light field acquisition process and the light field reconstruction process. Light field acquisition is to obtain multiple viewpoint information of a scene by using a real or virtual microlens array or camera array, that is, an elemental image array (EIA), and then restore the stereoscopic light field through the elemental image array and a lens array with the same parameters. In an integral imaging system implemented using circular lenses, the lens arrays used in the light field acquisition process and the light field reconstruction process are exactly the same. As a phase modulation device, the lens only modulates the propagation direction of the light beam and has a high light transmittance. Therefore, the design of the microlens array and the generation of the elemental image array are the key hubs of stereoscopic display.

[0004] An ideal optical system (also called a Gaussian optical system) can form a perfect image. However, in an actual optical system, due to non-ideal factors such as manufacturing materials and processes, the difference between the image formed by the actual optical system and the image formed in the paraxial region is called aberration. Aberration objectively exists in each actual optical system, resulting in the fact that the light beam from any object point cannot finally converge into a focused perfect image point, but forms a divergent diffusion spot. There are a total of 7 types of aberrations, among which spherical aberration, coma aberration, astigmatism, field curvature, and distortion are monochromatic aberrations, and axial chromatic aberration and lateral chromatic aberration are chromatic aberrations caused by light of different wavelengths.

[0005] Precisely controlling light using precision optical devices is the key to restoring three-dimensional scenes in three-dimensional light field display systems. However, actual three-dimensional light field display systems are not ideal optical systems. The inherent aberrations of lenses lead to a decline in the quality of reconstructed three-dimensional images, thus greatly limiting the application and implementation of three-dimensional light field display technology. In an integral imaging three-dimensional light field display system, lenses with a large field of view can improve the viewing angle of the three-dimensional light field display system. However, lenses with a large field of view usually introduce relatively severe aberrations, and severe aberrations will result in a reduction in the viewing angle of the three-dimensional light field display system and the overall blurring of the reconstructed images at the edge viewpoints.

[0006] Optical optimization is one of the powerful means to alleviate the adverse effects brought by lens aberrations. In "Depth-enhanced head-mounted light field displays based on integral imaging" (Wang X, Hua H. Optics Letters, 2021, 46(5): 985-988.), Wang et al. used a structure of a single-center lens array coupled with an optical fiber bundle to eliminate most of the off-axis aberrations. However, most of the above optical optimization methods may lead to an increase in system complexity, volume, or manufacturing difficulty. In "360° light field 3D display system based on a triplet lenses array and holographic functional screen" (Gao X, Sang X, Yu X, et al. Chinese Optics Letters, 2017, 15(12): 121201.), Gao et al. designed a triplet lens to balance the aberrations carried by a single lens, and increased the viewing angle of the three-dimensional light field display system to 70 degrees. In "A desktop light field display system with uniform viewpoint distribution" (Xu B, Yu XB, Gao X, et al. Chinese Journal of Liquid Crystals and Displays, 2022, 37(5): 573-580.), Xu et al. designed an aspherical lens to improve the problems of crosstalk between viewpoints and improper perspective relationship caused by lens aberrations. Inevitably, the use of multiple groups of lenses will increase the system complexity, and at the same time, it also poses more stringent requirements on the lens manufacturing process. Prefiltering is another means to mitigate the adverse effects of aberrations in the light field system. In "Wavefront aberration correction for integral imaging with the pre-filtering function array" (Zhang W, Sang X, Gao X, et al. Optics Express, 2018, 26(21): 27064-27075.), Zhang et al. proposed a method of using Wiener filter prefiltering to improve the three-dimensional light field display image. However, when using this method at a large viewing angle with severe aberrations, the resulting image will produce a "ringing effect", which will also seriously reduce the image quality. Therefore, this method is only applicable to three-dimensional light field display systems with a small viewing angle. In recent years, with the booming development of computer technology, machine learning technology has been applied to the optimization research of imaging quality in optical systems.In "Optical-aberrations-corrected light field re-projection for high-quality plenoptic imaging" (Chen Y, Jin X, Xiong B. Optics Express, 2020, 28(3): 3057-3072.), Chen et al. proposed a light field reprojection method that uses machine learning to correct the aberrations of a single-lens light field camera.

[0007] Integral imaging is a three-dimensional imaging technology using a lens array. The design of the lens array directly affects the imaging quality. Especially when the lens has aberrations, the imaging effect will be distorted. Therefore, designing a suitable lens array to reduce aberrations is one of the key issues. Aberrations objectively exist in the optical system, and existing technologies can only suppress aberrations within an acceptable range through various optimization means. Most of the current methods for suppressing aberrations are based on the primary aberration theory and use Zernike polynomials to describe the optical structure of the lens. Therefore, the calculated initial optical structure only considers the suppression of primary aberrations. In fact, there are still a large number of higher-order aberrations that need to be corrected. For example, the series expansion of spherical aberration is shown in Equation (1):

[0008] δU′ = a 1 L 2 + a 2 L 4 + a 3 L 6 +... (1)

[0009] where L is the lens aperture, the first term is the primary aberration, and the other terms are higher-order aberrations. a 1、 a 2、 a 3 … are coefficients. It can be seen from the above formula that the larger the aperture of the lens, the more serious the spherical aberration.

[0010] Various studies usually increase the display viewing angle by increasing the system complexity (including the optical structure complexity and the system assembly complexity). However, increasing the system complexity will bring errors to the assembly process, ultimately resulting in the display effect not meeting the design standards, and the increase in operation and process difficulty increases the processing time and cost of the display system.

[0011] In actual integrated imaging engineering calculations, the aberration correction of microlenses is expressed in the form of high-order ordinary differential equations, which can generally only be solved numerically. Traditional numerical solution methods such as the finite difference method or the Runge-Kutta method have great limitations: low accuracy and severely limited solvable range. When the size of the microlens is extremely small, the light propagation path is very sensitive to slight changes in the shape of the lens surface, making it difficult to obtain the exact path of the light. Summary of the invention

[0012] Aiming at the problem that the traditional numerical solution method used in the existing microlens array spherical aberration correction has low precision and severely limited solvable range, and cannot accurately describe the focusing behavior of the microlens array light, resulting in poor microlens array spherical aberration correction effect, the present invention proposes a microlens array spherical aberration correction method based on high-order linear ordinary differential equations.

[0013] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0014] A method for correcting spherical aberration of a microlens array based on a high-order linear ordinary differential equation comprises the following steps:

[0015] Step 1: Establish a high-order linear ordinary differential equation for describing the light propagation path of a single lens in a microlens array, introduce a non-constant coefficient function into the high-order linear ordinary differential equation, and set the initial conditions of the high-order linear ordinary differential equation;

[0016] Step 2: Obtain the initial coordinates of each pixel in the element image array, and based on the initial conditions and the initial coordinates, use a formulated reduced-order iterative numerical solution algorithm to solve the light field spatial geometric relationship and the high-order linear ordinary differential equation to obtain the aberration matrix of the element image array, wherein the aberration matrix includes the longitudinal offset and the lateral offset corresponding to each pixel, wherein the formulated reduced-order iterative numerical solution algorithm includes the following steps:

[0017] Step 2.1: After dividing the solution interval into steps, calculate the approximate value of the highest-order derivative at each division point;

[0018] Step 2.2: Using a numerical integration formula or Taylor formula, the low-order derivative values ​​of each segmentation point are uniformly expressed by the highest-order derivative value of the segmentation point, and the high-order linear ordinary differential equation is substituted to calculate the highest-order derivative value of each segmentation point;

[0019] Step 2.3: In order from high order to low order, the highest order derivative value of each segmentation point is used to iteratively obtain the numerical solution of each order derivative of all segmentation points, i.e., the aberration offset, and all the aberration offsets constitute the aberration matrix;

[0020] Step 3: Feed back the aberration matrix into the spatial function of the light field mapped by the elemental image array, and change the horizontal and vertical coordinates of the pixels of the corresponding elemental image according to the distance from the pixel point to the optical axis of the corresponding lens, so as to complete the spherical aberration correction of the microlens array.

[0021] The beneficial effects of the present invention are as follows:

[0022] (1) The present invention can obtain the offsets of the rays of all nodes in the specified solution interval relative to the optical axis, rather than just the offsets of the interval endpoints. Therefore, there are many controllable structural parameters for spherical aberration, with high precision, and the calculation precision required for 3D display can be arbitrarily controlled by adjusting the step size, thereby controlling the calculation time. Moreover, the approximate calculation part of the algorithm is basically not affected by the equation coefficient function terms, and the number of sensitive terms affecting the precision can be minimized to the greatest extent;

[0023] (2) The solution interval of the present invention is not affected by the coefficient function terms of the high-order differential equation of spherical aberration, and is only limited by the floating-point operation upper limit of the computer system. The calculation process is completely formulated and is easy to be programmed and implemented, having strong engineering application advantages. Description of the Drawings

[0024] Figure 1 is a flowchart of a method for correcting spherical aberration of a microlens array based on a high-order linear ordinary differential equation according to an embodiment of the present invention;

[0025] Figure 2 is a schematic diagram of the correction effect of the method for correcting spherical aberration of the microlens array of the present invention;

[0026] Figure 3 is a flowchart of a formula-based reduced-order iterative numerical solution algorithm;

[0027] Figure 4 is a schematic diagram of a simplified analysis of spherical aberration of a circular lens;

[0028] Figure 5 is a program interface and the control parameters included;

[0029] Figure 6 are the (a) numerical solutions, (b) exact solutions, (c) absolute errors, and (d) relative errors of y(x) and its derivative values of each order near x = 7;

[0030] Figure 7 is the calculation result of the professional mathematical software Maple;

[0031] Figure 8 are the (a) numerical solutions, (b) exact solutions, (c) absolute errors, and (d) relative errors of y(x) and its derivative values of each order near x = 10;

[0032] Figure 9Comparison chart of calculated value and accurate value;

[0033] Figure 10 Element image array before correction;

[0034] Figure 11 Element image array after correction;

[0035] Figure 12 Light field image generated from the element image array before correction;

[0036] Figure 13 Light field image generated from the element image array after correction. Detailed implementation manners

[0037] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and preferred embodiments.

[0038] Since each lens unit in the lens array is exactly the same, only one lens unit is analyzed in this embodiment. When light passes through the lens, because the actual lens is a non-ideal spherical surface, the light rays far from the optical axis will not be focused at the same focal point as the paraxial light rays. The light rays emitted from the object point will form a diffusion spot after passing through the lens system, that is, the actual image point. The difference between the actual image point and the ideal image point is called spherical aberration. Severe spherical aberration makes it difficult for the human eye to form stereoscopic vision in the peripheral field of view, thus greatly reducing the effective viewing angle and imaging clarity of the system. The most influential factor on the clarity quality of stereoscopic imaging is spherical aberration. The larger the spherical aberration, the larger the diffusion spot, which ultimately leads to problems such as blurred image quality.

[0039] Therefore, the present invention takes spherical aberration as an example to elaborate on the aberration correction method for the microlens array. It should be noted that the present invention is also applicable to the high-order calculation of other types of aberrations.

[0040] See Figure 1 and Figure 2 , this embodiment provides a spherical aberration correction method for a microlens array based on a high-order linear ordinary differential equation, specifically including the following steps 1 to 3.

[0041] Step 1: First, establish a high-order linear ordinary differential equation for describing the light ray propagation path of a single lens in the microlens array, introduce a non-constant coefficient function into the high-order linear ordinary differential equation, and also set the initial conditions of the high-order linear ordinary differential equation.

[0042] Since the arrangement and surface shape of each microlens in the integral imaging system will affect the light ray focusing path, which will introduce some non-linear effects, a high-order linear ordinary differential equation can be used to describe the light ray propagation path of a single lens, and a non-constant coefficient function is introduced to consider the correction of aberrations such as spherical aberration. The general mathematical model of the high-order linear ordinary differential equation is shown in Equation (2):

[0043]

[0044] Initial conditions:

[0045] Among them, y(x) is the offset of the light ray relative to the optical axis, that is, the relative spherical aberration displacement or the light ray offset of the lens array, which is used to describe the deviation of the imaging light ray; x is the distance that the light ray propagates along the optical axis, x ∈ [a, b], and a and b are the two endpoints of the solution interval; n is the order of differentiation; f(x) represents the inhomogeneous term of the external influence of the system, which describes the additional deviation caused by external factors. Among them, the coefficient functions a k (x) (k = 0, 1,... n - 1) are related to the path of the light ray, the properties of the lens material, and the surface shape, etc. For example, the coefficient function a 3 (x) reflects the third-order derivative effect caused by the non-uniform surface curvature of the microlens, and the coefficient function a 2 (x) represents the refractive index change when the light ray passes through different lens surfaces, and the coefficient function a 1 (x) is related to the lens surface shape and array arrangement, and the coefficient function a 0 (x) is related to the lens geometry and material properties. The specific function forms will be given in the application examples below.

[0046] Step 2: Obtain the initial coordinates of each pixel in the element image array, and based on the set initial conditions and the obtained initial coordinates, use the formula-based reduced-order iterative numerical solution algorithm to solve the simultaneous equations of the optical field spatial geometric relationship and the high-order linear ordinary differential equation, and obtain the aberration matrix of the element image array, where the aberration matrix includes the longitudinal offset and the lateral offset corresponding to each pixel in the element image array.

[0047] Traditional numerical solutions lack high-precision calculation methods, so the high-order model of spherical aberration is generally simplified to within the second order, and then each order derivative term in the discrete equation (or system of equations) is eliminated to obtain a discrete sequence of solutions. Aiming at the problem that high-order linear ordinary differential equations are not easy to solve and have large errors, the present invention proposes a formula-based reduced-order iterative numerical solution algorithm for lens spherical aberration correction. After the solution interval is segmented by step size, the discrete numerical solution of the highest-order derivative at each segmentation point is first obtained, and then the discrete numerical solution of the lower-order derivative at this point is recursively obtained from the obtained high-order solution until the target solution is obtained. The specific steps are as follows:

[0048] Step 2.1: After the solution interval is segmented by step size, calculate the approximate value of the highest-order derivative at each segmentation point;

[0049] Step 2.2: Use the numerical integration formula or Taylor formula to express the low-order derivative values at each segmentation point in terms of the highest-order derivative value at that segmentation point, and substitute them into the high-order linear ordinary differential equation to calculate the highest-order derivative value at each segmentation point;

[0050] Step 2.3: In the order from high order to low order, use the highest-order derivative values at each segmentation point to iteratively find the numerical solutions of the derivatives of all orders at all segmentation points, i.e., the aberration offsets. All aberration offsets form an aberration matrix.

[0051] As Figure 3 shown, the formula-based reduced-order iterative numerical solution algorithm is specifically as follows:

[0052] Divide the solution interval [a, b] into a = x 0 < x 1 < … < x m = b, with step size h = Δx i = x i - x i-1 , where i = 0, 1, 2, 3, …, m. For ease of calculation, it is required that the step size h is a fixed value and does not change with i.

[0053] For the case with initial conditions of , it is only necessary to divide the solution interval [a, b] into a = x m < … < x 1 < x 0 = b, and the step size h = Δx i = x i - x i-1 , where i = 0, 1, 2, 3, …, m.

[0054] Substitute the initial conditions into formula (2) to obtain formula (3):

[0055]

[0056] where represents the simplification of y (n) (x 0 ).

[0057] Starting from i = 1, y (n-r) (x i-1 ) is regarded as a known term, r = 1, 2, 3, …, n - 1, n. Select a certain numerical integration formula or Taylor formula to find and then obtain formula (4) from S(r, i).

[0058] y (n-r) (x i ) = y (n-r) (x i-1 ) + S(r, i) = φ r [y(n) (x i )] (4)

[0059] Among them, φ r [y (n) (x i )] is the numerical calculation formula of y (n-r) (x i ). Selecting different calculation methods for S(r, i) will result in different φ r [y (n) (x i )].

[0060] That is to say, using y (n) (x i ) to represent y (n-1) (x i ), and then using y (n-1) (x i ) to represent y (n-2) (x i ), so that in turn, y (n) (x i ) can be used to represent y (n-2) (x i ), y (n-3) (x i ), …, y′(x i ), y(x i )(note that this process does not involve equation (2) itself). Substituting all the obtained y (n-r) (x i ) into equation (2), at this time, equation (2) only contains one unknown y (n) (x i ), so as to find the approximate value of y (n) (x i ), denoted as The formula is as follows:

[0061]

[0062] Then use Substitute it back into equation (4) to successively obtain the approximate value of y (n-r) (x i ) Denoted as From this, the discrete numerical solutions of each order of the high-order linear ordinary differential equation (i.e., equation (2)) with respect to y(x) based on the step size h can be obtained That is, the highest-order derivative value at each segmentation point is obtained.

[0063] The approximate calculation part of this method has nothing to do with the ordinary differential equation itself and only uses Therefore, the calculation accuracy is for the coefficient term ak (x) is insensitive and the segmentation point x can be obtained simultaneously i Approximate values of each order derivative of y at y(x i ) = y (0) (x i ), i = 0, 1, 2, 3, …, m.

[0064] Next, different numerical integrations or Taylor formulas are used for the numerical integration term S(r, i) in Equation (4), and at the segmentation point x i (i = 1, 2, 3,..., m), the order reduction iterative calculation of y(x i ) is performed, and the corresponding different numerical solution formulas for the higher-order linear ordinary differential equation can be obtained.

[0065] The calculation of S(r, i) is implemented in the form of numerical integration. The following takes the trapezoidal numerical integration and Taylor formula methods as examples to give the formulaic numerical solution formulas.

[0066] (1) Trapezoidal numerical integration formula:

[0067] From it can be obtained that:

[0068]

[0069] By mathematical induction, Equation (6) is successively obtained:

[0070]

[0071] Among them:

[0072]

[0073] Substitute into Equation (2), and solve to obtain Equation (8):

[0074]

[0075] Then substitute the obtained as y (n) (x i ) into Equation (6), and then the respective values of

[0076] The following gives the pseudo-code of the formulaic order reduction iterative numerical solution algorithm in the programming process: Step 1:

[0077] Calculate y (n) (x 0 ) according to Equation (3).

[0078] Step 2:

[0079] i = 0; while i ≤ n

[0080] {

[0081] i = i + 1;

[0082] Let r = 1, 2, 3,..., n, and calculate d(r) according to Equation (7);

[0083] Calculate according to Equation (8)

[0084] Let r = 1, 2, 3,..., n, substitute y (n) (x i ) into Equation (6), and output successively

[0085] }

[0086] (2) Taylor's formula:

[0087] Perform the Taylor expansion of y(x) up to the (r = 0, 1, 2, 3,..., n - 1)-th order derivative of y at the point x i (x (r) (x i ) at the point x i-1 up to the (n - 1)-th order, as shown in Equation (9):

[0088]

[0089] Take According to Equation (9), obtain Equation (10):

[0090]

[0091] Substitute into Equation (5), and solve to obtain

[0092] where

[0093] a r (x i ) represents the value of the coefficient function at the point x i .

[0094] Then substitute the obtained into Equation (10), and calculate respectively to obtain

[0095] In summary, for the formulaic reduced-order iterative numerical solution algorithm adopted by the present invention, after the solution interval is segmented according to a specified step size based on the accuracy requirement, the approximate value of the highest-order derivative at each segmentation point is calculated. Secondly, using the specified numerical integration formula or Taylor formula, the values of each lower-order derivative at each segmentation point are uniformly expressed in terms of the highest-order derivative value at the segmentation point, and substituted into the original equation to calculate the highest-order derivative value at the segmentation point. Finally, in the order from high order to low order, the numerical solutions of the derivatives of all orders (including the 0th-order derivative) at all segmentation points are iteratively obtained using the highest-order derivative value at each point, that is, the spherical aberration offset. Since the number of terms involved in the approximate calculation is greatly reduced, this algorithm can significantly improve the calculation accuracy, expand the solution space, and the form of its solution is consistent for different numerical integration processing methods, with a completely formulaic solution formula, which is easy to implement through programming.

[0096] Spherical aberration calculation should originally be analyzed on the two-dimensional plane y(x r ,x t ) of the lens. Because the lens array for integral imaging uses circular lenses and the imaging system is rotationally symmetric, that is, the analysis of spherical aberration in the vertical and horizontal directions centered on the optical axis is equivalent, and only the object points in the meridional plane of a single lens need to be considered, and the imaging in other directions will be the same. Therefore, the above spherical aberration analysis simplifies the two-dimensional lens plane y(x r ,x t ) to the analysis including the optical axis and the vertical plane, that is, to solve y(x r )(abbreviated as y(x)), as Figure 4 shown. The solution process of the horizontal plane y(x t ) is exactly the same and will not be elaborated here.

[0097] Step 3: Feed back the aberration matrix to the spatial function of the light field mapped by the elemental image array, and change the pixel abscissa and ordinate positions of the corresponding elemental image according to the distance from the pixel point to the optical axis of the corresponding lens to complete the spherical aberration correction of the microlens array.

[0098] In this step, when the calculated spherical aberration value y(x r ,x t ) is fed back to the spatial function of the light field mapped by the elemental image array, according to the distance from the object point (the pixel point in the elemental image) to the optical axis of the corresponding lens, the pixel abscissa and ordinate positions of some elemental images are respectively changed, and the transformation relationship from the original pixel distribution to the spherical aberration-free pixel distribution can be obtained.

[0099] A point δ k (x r ,x t ,z) in the light field corresponds to a point δ k ′(x r ′,xt ′), (x r ′, x t ′) are the pixel coordinates of point δ k ′ on the element image display surface. Spherical aberration causes the imaging point of the light field to shift, resulting in the focused light becoming a light spot. Compensate the obtained spherical aberration value into the light field imaging point coordinates, and use the geometric relationship of the light path from the homologous light rays on the element image array passing through the lens array to the light field to feedback to the pixel points in the element image, and obtain a certain pixel point (x r , x t ) of the corrected element image and the mapping relationship with the light field, as shown in Equations (11) and (12).

[0100]

[0101] Among them, z = z(x r , x t ) is the depth of the light field at point (x r , x t ); is the optical center of the lens; P I is the diameter of the pixel in the reconstructed three-dimensional light field; P D is the pixel diameter on the display plane where the element image array is loaded, that is, the element image pixel size; w is the distance between the element image display surface and the lens array.

[0102] By feedback the target solution to the element image array according to the light information propagation path between the element image array and the lens array, without changing the parameters of the lens array, an element image array containing spherical aberration correction information is obtained, improving the spherical aberration correction ability of the integral imaging system, and further realizing the improvement of the three-dimensional imaging quality.

[0103] The following further explains the formula-based reduced-order iterative numerical solution algorithm:

[0104] The models targeted by classical algorithms (such as the Runge-Kutta method) are generally first-order problems For higher-order problems, they need to be transformed into first-order models for solution, and its iterative form is y i ′ = Φ(x i-1 , x i , y i-1 ), and is realized by By making an approximate calculation for . This integral requires the joint participation of each term of the equation and is sensitive to each coefficient term a k (x) of the higher-order linear ordinary differential equation.

[0105] For the higher-order spherical aberration solution problem in the form of Equation (2) of the present invention, only by making an approximation for Make an approximate calculation as shown in formula (4), then substitute formula (4) into formula (5) and solve the algebraic equation to obtain Finally, the Substituting equation (4) into equation (4) yields Therefore, the introduction of the approximate calculation part is irrelevant to the parameters of the high-order model (the rounding effect of its calculation is not considered here), and the relevant terms are only Calculation accuracy for coefficient term a k (x) and the right-hand term f(x) are both insensitive. This can maximize the analytical nature of the calculation process and reduce errors.

[0106] According to the requirements of engineering applications, the formula (4) can also be Select other numerical methods such as Simpson numerical integration, rectangular numerical integration algorithm, and substitute formula (4) into the original equation, and the corresponding numerical solution consistent with the structure of the aforementioned trapezoidal numerical integration algorithm can also be obtained. Different formula schemes have slightly different computational complexities, and there are also certain differences in the corresponding solution accuracy. Comparing the two methods selected in this embodiment, it is found that the computational complexity of using the trapezoidal numerical integration algorithm to calculate formula (8) is lower and the effect is better. This formula solution is also an important factor in the generalization and ease of programming of the present invention in engineering applications.

[0107] The present invention can obtain the offset of the light of all nodes in the specified solution interval relative to the optical axis, rather than just the offset of the interval endpoints. Therefore, the controllable structural parameters of spherical aberration are many and the accuracy is high. The calculation accuracy required for three-dimensional display can be arbitrarily controlled by adjusting the step size, thereby controlling the calculation time. The approximate calculation part of the algorithm is basically not affected by the coefficient function terms of the equation, and the number of sensitive terms that affect the accuracy can be minimized. In addition, the interval that the present invention can solve is not affected by the coefficient function terms of the high-order differential equation of spherical aberration, and is only limited by the floating-point operation upper limit of the computer system. The calculation process is completely formulated and easy to program and implement, with strong engineering application advantages.

[0108] The following only uses the "trapezoidal numerical integration algorithm" as an example for error analysis. Other methods are similar and will not be described in detail.

[0109] Let the step length h = Δx i =x i -x i-1 (In actual calculations, the step size is usually fixed and does not change with i). According to formula (4), it can be deduced that in the trapezoidal numerical integration algorithm, the single-step error is shown in formula (13).

[0110]

[0111] Among them, η is between x i -1 and xi at some point between.

[0112] Although the error is at the same level in form as that of the second-order Runge-Kutta method, since there is no definite absolute error and relative error for the numerical solution of differential equations, the order of error alone cannot fully explain the actual error level. The formulaic reduced-order iterative numerical solution algorithm of the present invention is only sensitive to h, y(x i-1 ), y (n-r) (x i-1 ). Like all numerical algorithms for differential equations with unilateral initial conditions, its error will generate transmission and progression as the solution interval extends and the number of steps increases. However, through extensive comparison with various classical algorithms, its accuracy and the length of the solvable interval are far higher than those of classical algorithms.

[0113] To verify the estimation performance of the present invention, a numerical solution calculation is carried out here considering an ordinary differential equation with a significantly monotonically increasing form. An executable general program of the formulaic reduced-order iterative numerical solution algorithm is written based on Delphi. Operating environment: Intel(R) Core(TM) i3-4010U@1.70GHz, memory is 4G, and the operating system is Win7. At the same time, the professional mathematical software Maple is used for comparison to perform the same numerical solution calculation on the selected complex equation.

[0114] The function used for testing is y = e 3x , and the corresponding equation is as follows:

[0115]

[0116] Input the initial value into the designed general program, and set to solve in the interval [0, 10]. The program interface and the included control parameters are as Figure 5 shown.

[0117] It can be seen that even for a solution function with exponential growth such as y = e 3x , when x = 7, the approximate value of y(7) calculated by the formulaic reduced-order iterative numerical solution algorithm is 1318811778.0449, and the error from the accurate value 1318811778.0419 is only 0.0029 (higher precision and significant figures are used in the internal calculation of the program, so it is not 0.0030), and the relative error is approximately 0( Figure 6 ). The absolute errors of its first to fourth derivative values are all less than 0.1, and the relative errors are basically 0. At the same time, the approximate value of y(7) obtained by the professional mathematical software Maple for the same calculation is 1318815734.43911, and the error from the accurate value 1318811778.0419 is 3956.39721, and the relative error is 0.000003( Figure 7)。The error is much larger than that of the formulated reduced-order iterative numerical solution algorithm, and the absolute errors of its derivative values of each order are also much larger than those of the formulated reduced-order iterative numerical solution algorithm.

[0118] The results of this example also show that the solvable intervals of the formulated reduced-order iterative numerical solution algorithm can far exceed the solution range of Maple. The maximum solvable interval calculated by Maple is [0, 7.5230122], while the corresponding calculation interval of the formulated reduced-order iterative numerical solution algorithm is only limited by the floating-point operation upper limit of the program running system. Only the numerical solutions of y(10), y(30), and y(100) are shown here ( Figure 7 ). The same conclusion is obtained for other examples and will not be elaborated here. From the perspective of computational time complexity, the calculation speed of Maple is better than that of the formulated reduced-order iterative numerical solution algorithm, but Maple does not give the approximate values of the functions and their derivative values of each segmentation point. However, the operation time of the designed program includes the output time-consuming for providing a large number of intermediate results, and it is not intuitive to compare.

[0119] The approximate value obtained by the formulated reduced-order iterative numerical solution algorithm for y(10) at the point x = 10 is 10686442522172.6943, and the error from the accurate value of 10686442522148.8067 is only 23.8876, and the relative error is less than 10 -12 , almost 0 ( Figure 8 ). It can also be found through Figure 9 that the relative error of the formulated reduced-order iterative numerical solution algorithm remains approximately 0 in the interval [0, 10]. The absolute errors of its derivative values of the 1st to 4th orders are 47.638, 70.836, 3.759, and 658.437 respectively, and the relative errors are almost all 0 ( Figure 8 ). Therefore, the numerical solution of the high-order equation targeted by the formulated reduced-order iterative numerical solution algorithm has high precision and can effectively ensure the accuracy of pixel movement during the rearrangement of element images.

[0120] To verify the improvement effect of spherical aberration correction on three-dimensional imaging, a specific application example is given below. In this application example, the resolution of the display screen is 1920×1080, which is composed of a 70×80 lens array circular aperture microlens. The lens diameter is 1 mm, and the focal length f is 2.5 mm. The distance w between the lens and the element image display surface is 3.15 mm, and the pixel diameter P D = 478 PPI (0.04725 mm) on the element image display surface. The depth range is 16.75 mm, and the pixel diameter P I = 0.246 mm on the central depth plane. The resolution of each element image is 18×18. The propagation path of light is described by a fourth-order linear ordinary differential equation, specifically expressed as Equation (15).

[0121]

[0122] Among them, the coefficient a 3 (x), a 2 (x), a 1 (x), a 0 (x), and f(x) are all functions of the variable x, reflecting the geometric characteristics, material properties of the integral imaging lens array, and the behavior of light on the surfaces of different microlenses.

[0123] a 3 (x) reflects the third-order derivative effect caused by the non-uniform surface curvature of the microlens. The surface curvature of the lens varies with position, and a typical curvature distribution can be expressed as: a 3 (x) = 1 / R 0 (1 + αx 2 ), usually it is a function of the position x, where R 0 = 150 μm is the initial curvature radius of the lens, and α = 0.0008 is a constant describing the curvature change.

[0124] a 2 (x) describes the refractive index change when light passes through different lens surfaces. In the integral imaging lens array, the refractive index n(x) of the gradient refractive index material varies with position, affecting the second-order derivative of light propagation. The refractive index of the integral imaging microlens array usually shows a parabolic distribution, so there is where β = 0.0003 is the refractive index change coefficient.

[0125] a 1 (x) is related to the lens surface shape and array arrangement. Since the lens array in integral imaging is actually an aspherical lens, its surface height z(r) is usually expressed by Taylor expansion as where r is the lens curvature, A 4 and A 6 are aspherical coefficients used to control the high-order effects. Here, take A 4 = -8×10 -7 , A 6 = 6×10 -10 . Therefore, a 1 (x) is expressed as a 1 (x) = A 4 x 3 +A 6 x 5 .

[0126] a 0 (x) represents the influence of the lens geometry and material properties on the light shift, and can be expressed as:

[0127] where R(x) is the radius of curvature of the lens, and k = 0.02 is the coefficient for correcting spherical aberration.

[0128] f(x) represents the external input to the system, such as the influence of the light source position, additional compensating optical elements, or the external environment (temperature, pressure, etc.) on the system, and is expressed as f(x) = ζx 2 , where ζ = 0.001.

[0129] Using the above initial conditions and the initial coordinates of each pixel of the element image array, the simultaneous solution of the spatial geometric relationship of the light field and the fourth-order linear ordinary differential equation is carried out to obtain the longitudinal and lateral offsets corresponding to 1,808,800 pixels, that is, the aberration matrix. The element image array before correction is as Figure 10 shown, and the element image array after correction is as Figure 11 shown, but the element image cannot intuitively reflect the improvement of the light field quality by the fusion of aberration correction information. Through Figure 12 and Figure 13 comparison, it can be clearly found that, when the resolution remains unchanged (the parameters of the light field remain unchanged), the visual clarity of the reconstructed light field by fusing aberration compensation information is better than that of the original light field.

[0130] When the present invention establishes the mathematical mapping relationship between the spatial data of the light field, it rearranges the element image array by using the spherical aberration correction information. According to the initial conditions, the original coordinates of each pixel on the display surface of the element image, and the differential equation, the actual aberration matrix of each pixel of the light field is calculated, and then through the optical path geometric relationship from the element image array to the lens array and then to the light field, the new light field depth matrix [z(x, y)] with fused aberration compensation information m×n is used to calculate the element image pixel coordinates corresponding to the three-dimensional pixels. Using this element image array, without changing the system hardware parameters such as the lens array, the image quality of the 3D light field image can be enhanced, and while ensuring the three-dimensional depth, the visual clarity of the light field can be improved.

[0131] The technical features of the above-described embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the technical features in the above-described embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as the scope described in this specification.

[0132] The above-described embodiments only represent several implementation manners of the present invention, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the invention patent. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several deformations and improvements can still be made, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the invention patent should be subject to the appended claims.

Claims

1. A method for correcting spherical aberration of a microlens array based on a high-order linear ordinary differential equation, characterized in that: The following steps are involved: Step 1: Establish a high-order linear ordinary differential equation for describing the light propagation path of a single lens in a microlens array, introduce a non-constant coefficient function into the high-order linear ordinary differential equation, and set the initial conditions of the high-order linear ordinary differential equation; Step 2: Obtain the initial coordinates of each pixel in the element image array, and based on the initial conditions and the initial coordinates, use a formulated reduced-order iterative numerical solution algorithm to solve the light field spatial geometric relationship and the high-order linear ordinary differential equation to obtain the aberration matrix of the element image array, wherein the aberration matrix includes the longitudinal offset and the lateral offset corresponding to each pixel, wherein the formulated reduced-order iterative numerical solution algorithm includes the following steps: Step 2.1: After dividing the solution interval into steps, calculate the approximate value of the highest-order derivative at each division point; Step 2.2: Using a numerical integration formula or Taylor formula, the low-order derivative values ​​of each segmentation point are uniformly expressed by the highest-order derivative value of the segmentation point, and the high-order linear ordinary differential equation is substituted to calculate the highest-order derivative value of each segmentation point; Step 2.3: In order from high order to low order, the highest order derivative value of each segmentation point is used to iteratively obtain the numerical solution of each order derivative of all segmentation points, i.e., the aberration offset, and all the aberration offsets constitute the aberration matrix; Step 3: Feedback the aberration matrix to the spatial function of mapping the element image array to the light field, and change the horizontal and vertical coordinates of the pixels of the corresponding element image according to the distance from the pixel point to the optical axis of the corresponding lens to complete the spherical aberration correction of the microlens array.

2. The method for correcting spherical aberration of a microlens array based on a high-order linear ordinary differential equation according to claim 1, characterized in that: The general mathematical model of the high-order linear ordinary differential equation is: Initial conditions: Among them, y(x) represents the offset of the light relative to the optical axis; x is the distance the light propagates along the optical axis, x∈[a,b], a and b are the two endpoints of the solution interval; f(x) represents the non-homogeneous term of the external influence of the system; a k (x) is the coefficient function, k=0,1,...n-1.

3. The method for correcting spherical aberration of a microlens array based on a high-order linear ordinary differential equation according to claim 2, characterized in that: In step 2.2, the process of solving the highest-order derivative value using the trapezoidal numerical integration formula includes: Divide the solution interval [a,b] into a=x0 <x1<…<x m =b, step length h = Δx i =x i -x i-1 , i=0,1,2,3,…,m, the step length h is a fixed value and does not change with i; Depend on We can get: Mathematical induction leads to formula (6): in: Will Substituting into equation (2), we can obtain equation (8): Then the desired As y (n) (x i ) into formula (2), we can respectively obtain r=1,2,3,...,n.

4. The method for correcting spherical aberration of a microlens array based on a high-order linear ordinary differential equation according to claim 3, characterized in that: When the initial condition is When , the solution interval [a, b] is divided into a=x m <… <x1<x0=b。 5. The microlens array spherical aberration correction method based on high-order linear ordinary differential equations according to claim 1 or 2, characterized in that: In step 2.2, the process of using Taylor's formula to solve the highest-order derivative value includes: Put y(x) at x i The r-order derivative y of a point (r=0,1,2,3,...,n-1) (r) (x i ) in x i-1 The point is expanded to the n-1th order Taylor expansion, as shown in formula (9): Pick According to formula (9), we can get formula (10): Will Substituting into equation (5), we get in, Then the desired Substituting into equation (10), we can find 6. The method for correcting spherical aberration of a microlens array based on a high-order linear ordinary differential equation according to claim 1 or 2, characterized in that: The numerical integration formula adopts Simpson numerical integration formula or rectangular numerical integration formula.

Images