the variational principle of leveling

Abstract

This work describes and practically implements the method of leveling illumination in images of various origins. Insufficient or uneven lighting can result in loss of detail, low contrast, or incorrect color reproduction. Lighting that is too bright can result in areas where detail is lost due to excessive brightness. Also, during the formation and analysis of various medical images, the problem of insufficient illumination of the image arises, that is, the image looks dark, details are difficult to see, especially in the shadows. The practical part of the work contains a description of the mathematical model of the transition from the variational problem (in the optimization form) to the differential form. A comparative analysis of existing computational methods of mathematical physics was also conducted. The working time and the number of iterations for obtaining results that form images with uniform illumination for various images were experimentally determined.

Keywords

• eol>Poisson's equation
• variational problem formulation
• computational methods
• elliptic second order equation
• equalization of illumination1

-

illumination in images

Anastasiia

Roman Yarovoy3†

1. Introduction

The problem of image illumination refers to how light affects the quality and perception of an image. If the lighting is insufficient or uneven, it can result in loss of detail, low contrast, or incorrect color reproduction. Lighting that is too bright can result in areas where detail is lost due to excessive brightness. Also, during the formation and analysis of various medical images, the problem of insufficient illumination of the image arises, that is, the image looks dark, details are difficult to see, especially in the shadows. One of the special problems is the uneven illumination of the object under study. Such a problem complicates the correct perception of this object.

Mathematical apparatus of varying complexity is used to analyze unevenly lit images and correct the above-mentioned lighting problems. The methods of leveling lighting in images use a matrix filtering apparatus (based on the principles of convolution), which is known from the course of linear algebra [1], [2]. This device is simple from the point of view of software implementation, and the processor time required to implement the implemented approach is relatively small. The Laplace and Fourier transform apparatus is more complicated than the previous one, but it gives a better visual result of correcting uneven lighting in the image at the output [1], [3]. The processor time for image processing with this mathematical apparatus is significantly greater than the time for image processing with matrix filters. The most wellknown mathematical apparatus for solving the problem of processing unevenly lit images is the variational principle, which involves the transition from the optimization formulation of the minimization of a quadratic functional of a special form to finding the solution of the Poisson equation with the appropriate boundary conditions that correspond to the given problem. To date, there are many methods of solving the main classes of problems of mathematical physics, as well as optimization problems of mathematical physics, where functional constraints are equations describing a specific process [3]–[5]. The main computational methods for finding the numerical solution of the Poisson equation are the finite difference method and the finite element method [2], [5], [6]. It should be noted that the finite element method has an advantage over the finite element method for solving the problem of equalizing uneven illumination in the case of examining an image area that has a very complex geometric shape. This paper will describe the image analysis method based on the variational principle and the finite difference method.

2. Problem formulation

2.1. Technical formulation of the problem

The technical task of the work consists in conducting an analysis of unevenly lit images based on the field of the given field of gradients of these images, to obtain a visually improved image. To conduct an analysis of classical computational methods of elliptic equations of mathematical physics of the second order, which are the basis of the variational approach of Poisson image processing. 2.2. Mathematical formulation of the problem: transition from the variational formulation of the problem to the differential formulation To correct unevenly distributed lighting in the image, you need to preserve all boundary elements of the objects that are in the image under study. In order to distinguish one object from another, it is necessary to have data about the gradient of the entire image, which is presented in the form of two matrices for each channel, if the image is full-color and is presented in digital form in the RGB palette [7]. The Poisson approach involves saving this information on the basis of such a quadratic functional I 1 (u), which must be minimized:

❑ ❑ I 1 (u )=∬ (u'x−v'x )2 ds +∬ (u'y−v'y )2 ds → min ,

G G where u – the searched image, which is the result of editing by saving, moreover u'x=∂ u / ∂ x , u'y=∂ u / ∂ y , v – input image for which field gradients are known, presented as matrices based on formulas v'x=∂ v / ∂ x , v'y=∂ v / ∂ y , G – the computational domain that defines the image itself.

It is obvious that in order to equalize the unevenly distributed illumination in the image, it is necessary to minimize another functional, which has the form:

❑ I 2 (u )=∬ (u−u )2 ds → min ,

G where u – average value of illumination. The functional I 2 (u ) consists in minimizing the variance. A reduction in dispersion implies a much more uniform distribution of the brightness of the pixels of the entire image. In this case, we have the functional of the variational formulation of the problem of leveling uneven illumination in the image given by and:

I (u )=α1 I 1 (u )+ α2 I 2 (u )=¿ ❑ ❑ ❑ ¿ α1(∬ (u'x−v'x )2 ds +∬ (u'y−v'y )2 ds)+α2∬ (u−u)2 ds → min ,

G G G

α1+ α2=1.

In functionality I (u ) parameters α1 and α2 є weighting factors that determine the strength of influence of functionals I 1 (u ) and I 2 (u ) respectively. It is obvious that the last functional I (u ) can also be written with one parameter:

I (u )=α1(1⋅ I 1 (u )+ α2 ⋅ I 2 (u )) I 1 (u )+ λ⋅ I 2 (u )=¿ α1 ❑ ❑ ❑ ¿∬ (u'x−v'x )2 ds +∬ (u'y−v'y )2 ds + λ∬ (u−u )2 ds → min ,

G G G

The last equivalence transformation “~” is true given that the extremal of the initial functional and the resulting functional are the same function. It should also be noted that λ=α2/ α1 , α1 ≠ 0. It is obvious that α1 ≠ 0, because at α1=0 all information about the initial image (field of image gradients) is lost. This means that at α1=0 the task loses any meaning.

In order to close the mathematical formulation of the variational problem, it is also necessary to mention the boundary conditions. Taking into account the fact that the boundary of the image objects is not tracked at the boundary of the entire image, and also that the minimization of variance for I 2 is better achieved with unfixed boundaries, then

| = | ∂ u ∂ u ∂ x ¿ ∂ x ¿= ∂∂uy|Up= ∂∂uy|Down= ∂∂un|Г=0 ,¿ where Г – border of the investigated image.

Such a variational problem can be easily reduced to the differential formulation of the same problem by applying the Euler-Lagrange equation, which has the form: ∂ S (u , u'x , u'y) − ∂

( ∂ u ∂ x ∂ S (u , u'x , u'y) ∂ u'x

In this case, we will get: ¿ 0+((u'y−v'y)2)'u' +0=2(u'y−v'y) .

y

Let's substitute the obtained intermediate results into the Euler-Lagrange equation. We will have: 2 λ (u−u)− ∂ (2(u'x−v'x))− ∂ (2(u'y−v'y))=¿

∂ x ∂ y ¿−2( ∂ (u'x−v'x)+

∂ x After simplification, we get the equation: ∂ (u'y−v'y)− λ (u−u)).

∂ y u'x'x+u'y'y− λu=v'x'x+ v'y'y− λ u .

The last equation is Poisson's equation. Taking into account the boundary conditions (derivatives along the normal are equal to zero), it is obvious that if the solution of the equation is some function u¿ ( x , y ), then the function will also be a solution u¿ ( x , y )+C, where C – an arbitrary constant. That is, we will prove:

(u¿+C )'x'x+(u¿+C )'y'y− λ (u¿+C )=v xx+ v yy− λ (u¿+C ) based on

It is obvious that

(u¿)'x'x+(u¿)'y'y− λ u¿=v xx+ v yy− λ u¿ . (u¿+C )'x'x+(u¿+C )'y'y− λ (u¿+C )−(v xx+ v yy− λ (u¿+C ))=¿ ¿ (u¿)' ' xx+(u¿)'y'y− λ u¿− λC −v xx−v yy + λ u¿+ λC =¿ (u¿)'x'x+(u¿)'y'y− λ u¿+ λ u¿−v xx−v yy=0

That is, this means that to solve the problem, you can put u=0 to fix the only solution of the problem, which will be easily pulled under the limits of the permissible pixel intensity values (for real numbers - from 0 to 1, or for integers numbers from 0 to 255), for example, according to the classic formula: where

λu− Δ u=− Δ f , λ>0 Δ u= ∂∂2xu2 + ∂∂2yu2 , Δ f = ∂∂2xf2 + ∂∂2yf2 ,

In ( 1 ), we assume that we have a known function f ( x , y ). This means that Δf is easily found by classical well-known difference schemes. Boundary conditions for equation ( 1 ) are equality of zero derivatives.

Suppose we need to write the difference equation for ( 1 ) [8]–[10]. Then we will use central difference schemes for any node ( xi , y j) , i=1 , N x+2 , j=1 , N y +2 . Steps in spatial coordinates x and y let's put equal hx , h y . Then: u= u−umin або u=[255⋅ u−umin ].

umax−umin umax−umin

Next, computational experiments will be conducted for the boundary value problem for the Poisson equation, which are known in classical mathematical physics.

3. Computational experiments

3.1. A numerical method for solving the problem

We have a mathematical model that contains a well-known linear equation of mathematical physics of the 2nd order [8], [9]:

| ∂2 u ∂ x2 (xi, y j)

| ∂2 u ∂ y2 (xi, y j) ≈ ≈ | ∂2 f ∂ x2 (xi, y j)

| ∂2 f ∂ y2 (xi, y j) ≈ ≈ , ( 1 ) ( 2 ) ( 3 ) u ( xi−1 , y j)−2 u ( xi , y j)+u ( xi+1 , y j) u ( xi , y j−1)−2 u ( xi , y j)+u ( xi , y j+1) . f ( xi−1 , y j)−2 f ( xi , y j)+ f ( xi+1 , y j) f ( xi , y j−1)−2 f ( xi , y j)+ f ( xi , y j+1) .

(hx)2 (h y )2 (hx)2 (h y )2

Similarly, we write the second derivatives in ( 1 ) for the known function f ( x , y ) in the righthand side of equation ( 1 ).

Let's enter the notation ui , j=u ( xi , y j) , f i , j=f ( xi , y j) , i=1 , N x+1 , j=1 , N y +1 . Then ( 2 ) and ( 3 ) will have the form: ∂2 u| ∂ x2 ( xi , y j) ∂2 f | ∂ x2 ( xi , y j) ≈ ui−1 , j−2hu2xi , j +ui+1 , j , ∂∂2yu2|( xi , y j) ≈ f i−1 , j−2hf2xi , j + f i+1 , j , ∂∂2yf2|( xi , y j) i=2 , N x +1 , j=2 , N y +1 .

≈ ui , j−1−2 ui , j +ui , j+1 ,

h2y ≈ f i , j−1−2 f i , j + f i , j+1 , h2y Using ( 4 ), we finally write down) the difference scheme for ( 1 ): λ ui , j−(

ui−1 , j−2hu2xi , j +ui+1 , j + ui , j−1−2hu2yi , j +ui , j+1 )=¿ ¿−( f i−1 , j−2hf2xi , j + f i+1 , j + f i , j−1−2hf2yi , j + f i , j+1 ),

i=2 , N x +1 , j=2 , N y +1 .

Let's transform equation ( 5 ), collecting all the coefficients near the points ( xi , y j ) , ( xi−1 , y j ) , ( xi+1 , y j ) , ( xi , y j−1) and ( xi , y j+1), and we will see the 5-point difference scheme for equation ( 1 ) in the new notation.

ui−1 , j−2 ui , j +ui+1 , j + ui , j−1−2 ui , j +ui , j+1 − λ ui , j=¿

h2x h2y ¿ f i−1 , j−2 f i , j + f i+1 , j + f i , j−1−2 f i , j + f i , j+1 , i=2 , N x +1 , j=2 , N y +1 .

h2x h2y h12x ui−1 , j + h12x ui+1 , j + h12y ui , j−1+ h12y ui , j+1−( h22x + h22y + λ)ui , j=¿ ¿ f i−1 , j−2 f i , j + f i+1 , j + f i , j−1−2 f i , j + f i , j+1 , i=2 , N x +1 , j=2 , N y +1 .

h2x h2y We will introduce new notations in the difference equation ( 6 ):       coefficient near ui−1 , j: A X i , j=1 / h2x; coefficient near ui+1 , j: C X i , j=1 / h2x; coefficient near ui , j−1: A Y i , j=1 / h2y; coefficient near ui , j+1: C Y i , j=1 / h2y; coefficient near ui , j: Bi , j=2 / h2x +2 / h2y + λ; coefficient of the right side of the equation Di , j:

Di , j= f i−1 , j−2 f i , j + f i+1 , j + f i , j−1−2 f i , j + f i , j+1 .

h2x h2y ( 4 ) ( 5 ) ( 6 ) In this case, ( 6 ) will turn into the following differential equation:

A X i , j ui−1, j+C X i , j ui+1, j+ A Y i , j ui , j−1+C X i , j ui , j+1−Bi , j ui , j= Di , j , i=2 , N x +1 , j=2 , N y +1 . ( 7 )

Equation ( 7 ) contains the unknown function u (ui , j) at each point of the calculation area ( xi , y j) , i=2 , N x +1 , j=2 , N y +1.

For ( 7 ), you can apply any known computational method of linear algebra, or use multi-grid methods, which will significantly speed up the calculations, given that the number of unknowns in the calculation area is large (pixels of the image). Also, it should be noted that taking into account that λ >0, then the difference equation ( 7 ) has a diagonal advantage, which consists in the fact that in each row the module of each diagonal element of the SLAR is not less than the sum of the modules of the other elements of the row [8]–[10].

Now, to close the issue of the model, consider the consideration of the zero derivative at all boundaries of the calculation area [10].

Let our calculation area be in a rectangle [ xa ; xb ] × [ ya ; yb ] . To begin with, consider the zero derivative on the left boundary: ∂ u| =0.

∂ x (xa; y) Based on this, we need to adjust equation ( 7 ) on the left boundary: Given that

A X 2, j u1, j+C X 2, j u3 , j+ A Y 2, j u2, j−1+C X 2, j u2, j+1−¿

−B2, j u2, j= D2, j , j=2 , N y +1 .

| ∂ u ∂ x (xa; y) ≈ u2, j−u1, j =0 , u1, j=u2, j , j=2 , N y +1 .

hx Substitute the obtained boundary condition ( 9 ) on the left into ( 8 ): Finally, let's transform ( 10 ) into the form:

A X 2, j u2, j+C X 2, j u3 , j+ A Y 2, j u2, j−1+¿ +C X 2, j u2, j+1−B2, j u2, j= D2, j , j=2 , N y +1 .

C X 2, j u3 , j+ A Y 2, j u2, j−1+C X 2, j u2, j+1−¿ −( B2, j− A X 2, j) u2, j= D2, j , j=2 , N y +1 .

In equation (11), it should be noted that the coefficients are adjusted:

B2, j=B2, j− A X 2, j , A X 2, j=0 , j=2 , N y +1 .

Now we need to correct equation ( 7 ) on the right boundary: ( 8 ) ( 9 ) ( 10 ) (11) ∂∂ux|(xb; y) ≈ uNx+2, j−uNx+1, j =0,uNx+1, j=uNx+2, j , j=2, N y+1.

hx

XNx+1, juNx, j+C XNx+1, juNx+1, j+ A Y Nx+1, juNx+1, j−1+¿ +C XNx+1, juNx+1, j+1−BNx+1, juNx+1, j=DNx+1, j , j=2, N y+1.

Substitute into (12) the obtained boundary condition (13) on the right: Substitute into (16) the obtained boundary condition (17) on the right: (12) (13) (14) (15) (16) (17) (18) (19) Finally, let's transform (14) into the form:

XNx+1, juNx, j+ A Y Nx+1, juNx+1, j−1+C XNx+1, juNx+1, j+1−¿ −(BNx+1, j−C XNx+1, j)uNx+1, j=DNx+1, j , j=2, N y+1.

In equation (15), it should be noted that the coefficients are adjusted:

BNx+1, j=BNx+1, j−C XNx+1, j ,C XNx+1, j=0 , j=2, N y+1.

We adjust equation ( 7 ) at the lower limit:

A Xi,2ui−1,2+C Xi,2ui+1,2+ A Y i,2ui,1+C Xi,2ui,3−¿ −Bi,2ui,2=Di,2,i=2, N x+1.

≈ ui,2−ui,1=0,ui,1=ui,2,i=2, N x+1.

∂∂uy|(x; ya) hy A Xi,2ui−1,2+C Xi,2ui+1,2+ A Y i,2ui,2+C Xi,2ui,3−¿

−Bi,2ui,2=Di,2,i=2, N x+1.

A Xi,2ui−1,2+C Xi,2ui+1,2+C Xi,2ui,3−¿ −(Bi,2−A Y i,2)ui,2=Di,2,i=2, N x+1.

Finally, let's transform (18) into the form: In equation (19), it should be noted that the coefficients are adjusted:

Bi,2=Bi,2−A Y i,2, A Y i,2=0,i=2, N x+1.

We adjust equation ( 7 ) at the lower limit:

A X i , N y+1 ui−1, N y+1+C X i , N y+1 ui+1, N y+1+ A Y i , N y+1 ui , N y+C X i , N y+1 ui , N y+2−¿

−Bi , N y+1 ui , N y+1= Di , N y+1 , i=2 , N x +1 . Substitute into (20) the obtained boundary condition (21) on the right: Finally, let's transform (22) into the form:

A X i , N y+1 ui−1, N y+1+C X i , N y+1 ui+1, N y+1+ A Y i , N y+1 ui , N y+¿ +C Y i , N y+1 ui , N y+1−Bi , N y+1 ui , N y+1= Di , N y+1 , i=2 , N x +1 .

A X i , N y+1 ui−1, N y+1+C X i , N y+1 ui+1, N y+1+ A Y i , N y+1 ui , N y−¿

−( Bi , N y+1−C Y i , N y+1) ui , N y+1= Di , N y+1 , i=2 , N x +1 .

In equation (23), it should be noted that the coefficients are adjusted:

Bi , N y+1=Bi , N y+1−C Y i , N y+1 , C Y i , N y+1=0 , i=2 , N x +1 .

For a rectangular area, the values at the corners of this area also raise questions. Here, too, everything is simple. As an option, you can make "stubs" as the arithmetic mean of two known nearest neighboring limit values: u1,1= u1,2+2 u2,1 , uN x+2, N y+2= uN x+1, N y+2+2 uN x+2, N y+1 , u1, N y+2= u1, N y+1+u2, N y+2 , uN x+2,1= uN x+1,1+uN x+2,2 .

2 2

It should be noted that the boundary nodes of the region (24) do not take part in the calculation procedures, since the 5-point difference scheme does not pass through these nodes.

3.2. Testing the numerical method on various images

The resulting algorithm was tested on 30 different images. The timing results for the first six images are shown in Tables 1–4. It should be noted that the block sweep algorithm performed better than all others for solving the discrete analog for the set task of equalizing the uneven distribution of pixel intensity in the image. The algorithm of the Scipy library of the Python programming language also performed well [11]–[13]. It should also be noted that computational methods of linear algebra more effectively solve systems of linear algebraic equations when the matrices of this system have a diagonal advantage [14], [15]. The greater the difference between the modulus of the diagonal element and the sum of the modules of the other elements of the row (and for each row), the faster the iterative method or algorithm converges. (20) (21) (22) (23) (24) 1 2 3 4 5 6 1 2 3 4 5 6 216 450 216 249 192 341 216 450 216 249 192 341 197 297 218 186 256 404 197 297 218 186 256 404

Next, we test the same algorithm based on the already known approximate algorithm of longitudinal-transverse running (the results are shown in Table 3).

We will do the same for the analytical block sweep algorithm, which makes it possible to obtain an exact solution to the problem, taking into account the fact that the discrete analogue of the Poisson equation has a clear sparse and block structure. (results are shown in Table 4). Table 4 Time characteristics of finding a solution to the problem using the matrix run method for six different images

Image VPeirxteiclsal HoPriixzeolnstal PTioxtealls

The visual results of the program are shown below.

As can be seen from Figure 1, all images have a better pixel intensity distribution as a result of applying the algorithm.

4. Conclusions

The practical result of this work is a comparative analysis of unevenly lit images based on the field of the given field of gradients of these images, to obtain a visually improved image.

In the work, mathematical dependencies are obtained for constructing an approximate solution to the set problem of equalizing the uneven distribution of pixel intensity. An analysis of classical computational methods of second-order elliptic equations of mathematical physics, which are the basis of the variational approach of Poisson image processing, was carried out. Seidel's method performed the worst in terms of time, the block run method using the built-in Python library Scipy gave the best result in terms of time.

In order to obtain better time indicators for a specific accuracy, it is best to use multi-grid approximate methods, which make it possible to significantly reduce the number of iterations, which means to reduce the time to execute the program.

It should also be noted that the speed of iterative methods and algorithms will be higher if calculations are performed in the normal system (from 0 to 1) and not in the classic system (from 0 to 255). [11] Azimzadeh P., Forsyth P.A. (2016) Weakly Chained Matrices, Policy Iteration, and Impulse

Control. SIAM J. No. 54 ( 3 ), 1341–1364. [12] Chen Y., Wan J., Lin J. (2018) Monotone mixed finite difference scheme for Monge–Ampère equation. J. Sci. Comput. 76, 1839–1867. [13] Chen Y., Wan J. (2017) Multigrid methods for convergent mixed finite difference scheme for

Monge-Ampère equation. Comput. Visual. Sci., 1–15. [14] Chen Y., Wan J. (2018) Numerical method for image registration model based on optimal mass transport. Inverse Problems Imaging 12( 2 ), 401–432. [15] Chow S.N., Li W., Zhou H. (2019) A discrete Schrödinger equation via optimal transport on graphs. Journal of Functional Analysis 276, 2440–2469.

References

[1] Ganesh R. Naik , Wellington P. ( 2024 ) Biomedical Signal Processing. A Modern Approach . ISBN 10: 103206191X, ISBN 13: 9781032061917

[2] Xue-Cheng T. , Suhua W. , Haiguang L. ( 2018 ) Mathematical Methods in Image Processing and Inverse Problems: IPIP 2018, Beijing, China, April 21-24 (Springer Proceedings in Mathematics & Statistics , 360 ) 1st ed. 2021 . DOI: https://doi.org/ 10.1007/ 978 -981-16-2701-9.

[3] Zhou H. , Qin R. , Liu Z. , Qian Y. , Ju X. ( 2022 ) Optimizing Performance of Image Processing Algorithms on GPUs . Proceeding of 2021 International Conference on Wireless Communications, Networking and Applications. WCNA 2021. Lecture Notes in Electrical Engineering . Springer, Singapore. DOI: https://doi.org/10.1007/ 978 -981-19-2456-9_ 95 .

[4] Mitsunori G. ( 2019 ) Image Processing Using ImageJ . Japanese Journal of Radiological Technology 75, no. 7 . DOI: http://dx.doi.org/10.6009/jjrt.2019 _jsrt_75.7 .688.

[5] Sivakumar S. , Saminathan S. , Ranjana R. , Mohan M. , Pareek P. ( 2023 ) Malware Detection Using the Machine Learning Based Modified Partial Swarm Optimization Approach , International Conference on Applied Intelligence and Sustainable Computing (ICAISC) , Dharwad, India, 2023 , pp. 1 - 5 . DOD: 10 .1109/ICAISC58445. 2023 . 10199796 .

[6] Lande J. , Maheswari M. , Kamali S. M. , Dineshkumar R. ( 2023 ) Hybrid Optimization Based Long Short-Term Memory for Android Malware Detection . International Conference on Evolutionary Algorithms and Soft Computing Techniques (EASCT) , Bengaluru, India, 2023 , pp. 1 - 5 . DOD: 10 .1109/EASCT59475. 2023 . 10393408 .

[7] Kotian P. , Sonkusare R. ( 2021 ) Detection of Malware in Cloud Environment using Deep Neural Network , 6th International Conference for Convergence in Technology, Maharashtra, India, 2021 , pp. 1 - 5 . DOI: 10 .1109/I2CT51068. 2021 . 9417901 .

[8] Khaidurov V. , Zaporozhets A. , Tsiupii T. ( 2022 ) Creation of High-Speed Methods for Solving Mathematical Models of Inverse Problems of Heat Power Engineering . Springer, Systems Decision and Control in Energy III , vol. 399 , 41 - 74 . DOI: https://doi.org/10.1007/978-3- 030 -87675- 3 . ISSN: 2198 - 4182 .

[9] Khaidurov V. , Yailymov B. , Shelestov A. ( 2023 ) Mathematical Model for Determining the Geometric Location of the Environmental Pollutant Based on Sensor Data . Proceedings of the IEEE International Conference on Intelligent Data Acquisition and Advanced Computing Systems: Technology and Applications , IDAACS, pp. 703 - 707 .

[10] Khaidurov V. , Tsiupii T. , Zhovnovach T. ( 2021 ) Modelling of ultrasonic testing and diagnostics of materials by application of inverse problems . CEUR Workshop Proceedings , 2021 , 3039, Pp. 1 - 5 .