TWI639135B - Restoration method for blurred image - Google Patents

Abstract

The invention discloses a method for restoring a blurred image, which comprises the following steps: (1) Point spread function estimation: first performing a Fourier transform on the blurred image, and then taking the Fourier transformed value into its Cepstrum, and then from the inverted Calculate the Blur Angle and Blur Length of the Point Spread Function (PSF) in the spectral domain; (2) Deblur the single image: Spread the estimated blur angle and blur length to the point spread The function is normalized, and then the deconvolution of the single image is performed on the normalized point spread function; (3) the clear image register: the clear image is selected or the image is blurred by a single image. Clear Image Register; (4) Multiple image deblurring: first search for the corresponding feature points between the current image and the latest image stored in the clear image register, and use these feature points to derive the perspective transformation (Holeography Transformation) required to correct the corresponding pixels of the two images, and then calculate these two according to Temporal Information. Right image pixel weight (Weight) and compared to the weight of low substituted weights high pixel weight of the pixel way to generate restored pixel, and finally output a clear image of this recovery. .

Description

Fuzzy image restoration method

The invention relates to a method for restoring a blurred picture, in particular to a method for restoring a blurred image.

Due to the advancement of science and technology in recent years, computer vision technology has been widely used in important institutions or traffic monitoring systems for license plate identification to achieve a smart surveillance system for maintaining law and order and crime prevention. Motion blur can occur due to some external factors, but these fuzzy problems may also cause misjudgment of subsequent identification applications, so it is necessary to repair and compensate the blurred images.

In the field of image blur reduction, most of the images are deblurred (Image Deblurring), but in the single fuzzy reduction, there is a good reduction result for the artificial motion-Blurred Image, or constant velocity. The result of motion blur (Uniform Motion) is better, but it does not have better results in the restoration of non-uniform motion (Non-Uniform Motion) blurred image, and the calculation of fuzzy reduction is too large, but The requirements of related equipment and equipment are relatively increased, and the cost is also increased.

In view of the above-mentioned problems, it is an object of the present invention to provide a method for restoring a blurred image to solve the problems faced by the prior art.

Based on the above object, the present invention provides a method for restoring a blurred image, which comprises the steps of: performing a Fourier transform on a blurred image to generate a converted value, and then obtaining a inverse spectrum of the converted value, and then calculating a point spread from the cepstral domain. The fuzzy angle and the fuzzy length of the function; the point spread function is normalized according to the blur angle and the blur length, and then the normalized point spread function is subjected to the deconvolution of the single image to generate the restored image; the clear image or the restored image is selected. Clear image buffer; and searching for corresponding feature points between the current image and the restored image stored in the clear image register, using the feature points to derive the homography matrix required for the perspective transformation, and correcting the current image and The corresponding pixels of the image are restored, and the weights of the pixels in the current image and the restored image are calculated according to the time domain information and compared, and the low-weight pixels are replaced by high-weight pixels to generate the restored pixels, and another restored image is output.

Preferably, the clear image buffer stores the clear image in the image sequence when the image sequence is input, and the blurred image in the image sequence is restored to the restored image by the restoration method of the blurred image.

Preferably, the point spread function can be expressed by the following formula: Among them, L is the blur length, and θ is the blur angle.

Preferably, the cepstrum domain definition of the blurred image can be expressed by the following formula: C ( g ( x , y )) = F ^-1 {log| F ( g ( x , y ))|} where F represents a Fourier transform, F ^-1 represents an inverse Fourier transform, and g(x, y) is a blurred image.

Preferably, calculating the blur length and the blur angle of the point spread function may include the following steps: (1) inputting a cepstrum map to find a center point; (2) from a center point of 90°, 45°, and 0° Select the point of the largest cepstrum value as the search point; (3) Move to the search point and use a 5*5 mask to mask the center corresponding search point; (4) Cut the cepstrum in the mask The values are all added and the coordinates of the search point are stored; (5) determining whether the direction of the search point is offset, if otherwise, proceeding to step (6), if yes, proceeding to step (7); (6) calculating the distance between the search point and the center point Whether it is greater than the threshold, otherwise return to step (2) to repeat steps (2) to (5), and if so, proceed to step (7); (7) find the turning point of the falling and rising from the trajectory value distribution map; 8) Calculate the blur length and the blur angle according to the turning point and the center point; where the coordinates of the turning point are (x ₁ , y ₁ ), the coordinates of the center point are (x ₀ , y ₀ ), L is the blur length, and substituting (x ₁ , y ₁ ) and (x ₀ , y ₀ ) to obtain the fuzzy length are expressed by the following formula: R is the length of the base of the triangle formed by the turning point, the center point and the X-axis. Substituting x ₁ and x ₀ to obtain the length of the base is expressed by the following formula: R =| x ₁ - x ₀ |θ is the blur angle, Substituting the blur length and the length of the base to obtain the blur angle is expressed by the following formula:

As described above, the method for restoring a blurred image of the present invention can be applied to the restoration of a linear motion blur picture, and the proposed video deblurring method can be implemented in a common image monitoring system or a general moving vehicle. Device, which can be used to instantly blur the video through video streaming Yes, because the algorithm proposed in this case has a small computational complexity, it can be embedded in a general security surveillance camera system to enhance the added value of related products.

S11 to S14, S61 to S68‧‧‧ steps

FIG. 1 is a first flowchart of a method for restoring a blurred image of the present invention.

Figure 2 is a cepstrum of a constant velocity motion blur image.

Figure 3 is a cepstrum of a non-equal motion blurred image.

Figure 4 is a second flow chart of the method for restoring a blurred image of the present invention.

Figure 5 is a search trajectory map of the cepstrum.

Figure 6 is a 3D diagram of the cepstral domain.

Figure 7 is a search trajectory of the non-equal motion blur cepstogram.

Figure 8 is a trajectory numerical distribution map.

Figure 9 is a schematic diagram showing the results of single image restoration.

Figure 10 is a schematic diagram of the transformation between two images performed by the homography matrix.

Figure 11 is a schematic diagram showing the results of corner matching: (a) part is the t-1 picture; (b) part is the tth picture; (c) part is the corner point of the t-1 picture and the tth picture Match the result.

Figure 12 is a schematic diagram showing the results of perspective transformation.

Figure 13 is a schematic diagram showing the results of multiple image restoration processes.

In order to understand the features, contents, and advantages of the present invention, and the advantages thereof, the present invention will be described in conjunction with the drawings, and the description of the embodiments will be described in detail below. The use of instructions and supplementary instructions is not necessarily true to the actual ratio and accuracy of the invention. It is to be understood that the scope of the present invention is not to be construed as limiting the scope of the present invention.

The advantages and features of the present invention, as well as the technical methods of the present invention, are described in more detail with reference to the exemplary embodiments and the accompanying drawings, and the present invention may be implemented in various forms and should not be construed as limited thereby. The embodiments of the present invention, and the embodiments of the present invention are intended to provide a more complete and complete and complete disclosure of the scope of the present invention, and The scope of the patent application is defined.

Please refer to FIG. 1 , which is a first flowchart of a method for restoring a blurred image of the present invention. As shown, the blurred image restoration method 100 of the present invention comprises the following steps:

In step S11, the blurred image is subjected to Fourier transform to generate a converted value, and then the scrambled spectrum of the converted value is obtained, and then the blur angle and the blur length of the point spread function are calculated from the cepstral domain.

In step S12, the point spread function is normalized according to the blur angle and the blur length, and then the normalized point spread function is subjected to deconvolution of the single image to generate a restored image.

In step S13: a clear image or a restored image is stored in the clear image register.

In step S14: searching for corresponding feature points between the current image and the restored image stored in the clear image register, and deriving the homography matrix required for the perspective transformation by using the feature points, and correcting the current image and the restored image according to the feature point. Corresponding pixels are used to calculate and compare the weights of the pixels in the current image and the restored image according to the time domain information, and replace the low-weight pixels with high-weight pixels to generate restored pixels, and output another restored image.

Please refer to Figures 2 to 13 again; Figure 2 is a cepstrum image of a constant velocity motion blur image; Figure 3 is a cepstrum image of a non-equal motion blur image; Figure 4 is a blurred image of the present invention. Complex The second flowchart of the original method; the fifth graph is the search trajectory map of the cepstrum map; the sixth graph is the 3D map of the cepstrum domain; the seventh graph is the search trajectory map of the non-equal motion blur cepstogram Figure 8 is a trajectory numerical distribution map; Figure 9 is a schematic diagram of the results of single image restoration; Figure 10 is a schematic diagram of the transformation between two images for the homography matrix; Figure 11 is a corner matching The result is a schematic diagram: (a) part is the t-1 picture; (b) part is the tth picture; (c) part is the matching result of the corner point of the t-1 picture and the tth picture; the 12th picture is the perspective Schematic diagram of the result of the transformation; Figure 13 is a schematic diagram of the results of multiple image restoration processes. The above steps will be described in detail below in conjunction with the various drawings.

1. Point spread function estimation

The point spread function can be regarded as a pulse function in an optical system. In mathematics, the point source (input) can be represented by a point pulse function, and the output light field distribution can be called an impulse response (Impulse Response), and this represents the pulse received by the image. response. If the imaging system produces an inverted image, the image plane coordinate axis can simply be reversed from the object plane coordinate axis. Calculating the image plane fold integral is a simple process without distortion.

1.1, image conversion

Most of the diffusion function estimation methods are processed in the spatial domain. The main way is to estimate the local feature points such as points, edges and lines. In these methods, if you want to estimate the point spread function, you must The type of blur is known in advance, and the accuracy of these methods will decrease when the blur is more serious. Therefore, the system is based on the frequency domain method to calculate the length ( L ) and angle (θ) of the blur.

Among them, the point spread function of linear motion blur is approximated by equation (1).

Among them, L is the blur length, and θ is the blur angle.

1.2, cepstral domain conversion

The cepstrum domain was first proposed in the field of one-dimensional signal processing and was later introduced into image processing. The cepstrum domain of the blurred image g ( x, y ) is defined as equation (2), where F and F ^-1 respectively represent Fourier Transform, inverse Fourier transform, from which the cepstrum domain of the image is the logarithm of the original image power spectrum and then inverse Fourier transform.

C ( g ( x , y ))= F ^-1 {log| F ( g ( x , y ))|} (2)

1.3, calculate the blur length and angle

In the experiment of the present invention, it is found that in the case of constant velocity motion blur, there are two distinct black points in the cepstrum diagram, as shown by the white circle in Fig. 2, and we only need to calculate between the center point and the point. The distance is also the blur length of the image, and then the angle θ formed by the point and the X axis is calculated, and this θ is also the blur angle; but if it is non-equal motion blur, there will be no obvious black spots, as shown in Fig. 5. As shown, the present invention therefore proposes an algorithm that calculates the blur length and angle whether it is a constant velocity or a non-equal motion blur. As shown in FIG. 4, the algorithm flow includes the following steps:

In step S61: input a cepstrum map to find a center point.

In step S62, the point of the largest cepstrum value is selected from the three directions of 90°, 45°, and 0° of the center point to use it as a search point.

In step S63: moving to the search point and using a 5*5 mask, the mask center corresponds to the search point.

In step S64: all the cepstrum values in the mask are added and the coordinates of the search point are stored.

In step S65, it is judged whether or not the direction of the search point is shifted. If not, the process proceeds to step S66, and if yes, the process proceeds to step S67.

In step S66, it is calculated whether the distance between the search point and the center point is greater than the threshold value. Otherwise, the process returns to step S62 to repeat steps S62 to S65, and if yes, the process proceeds to step S67.

In step S67, the turning point of the falling and then rising is found from the trajectory value distribution map.

In step S68, the blur length and the blur angle are calculated according to the turning point and the center point.

As stated above, assuming that a cepstrum diagram of non-equal motion blur is input, the first step is to find the center point, and the second step is selected from the three directions of 90°, 45°, and 0° of the center point. The point of the largest cepstrum value is used as the search point, the third step is moved to the search point, then a 5*5 mask is used, the mask center corresponds to the search point, and the fourth step adds all the cepstrum values in the mask. And storing the search point coordinates, the fifth step is to determine whether the search point direction is offset, if not offset, calculate whether the distance between the search point and the center point is greater than the threshold, if otherwise, repeat the above steps in the second step, if it is in the first The search is stopped when the five-step offset occurs, as shown by the trace in the white circle in Figure 5.

The reason for the deflection is because the value at the falling point is mostly negative, so it will deflect along the positive value around the search. The depression can be clearly seen from the 3D image in Fig. 6, but Non-constant motion blur usually does not have a significant drop point, so there will be no deflection when searching for the fifth step, as shown in Figure 7, then enter the sixth step to determine if the distance between the search point and the center point exceeds The threshold is stopped, and the seventh step is passed through the points stored in each mask, and the turning point that turns to rise when the value gradually decreases can be calculated. It can be clearly seen from the numerical distribution map of the trajectory in FIG. Finally, calculate the distance between the turning point and the center point, assuming that the turning point is ( x ₁ , y ₁ ), the center point is ( x ₀ , y ₀ ), and then calculate the two points by using the Euclidean distance formula (3). The distance between them is the blur length L , and the two points and the x-axis can be regarded as one triangle, L is the oblique side and the bottom side is R. The length of the base can be calculated by the formula (4), and then the formula is passed. 5) The inverse cosine formula can be used to calculate the angle.

R =| x ₁ - x ₀ | (4)

2, single image to blur

In this step, it is mainly divided into two parts, namely, the point spread function normalization and the single image deconvolution. Since the point spread function (PSF) is unknown before the restoration, the iterative operation requires an initialization seed point. Therefore, the initial point spread function can be obtained by the above-mentioned calculated blur length L and angle θ , and then calculated by the Richardson-Lucy iterative restoration algorithm, and the Richardson-Lucy calculation is performed by the formula (6). In the formula (6), B represents a blurred image, the point spread function in accordance with the result of the sharp image I _m average of the different positions, the CM image I _m is the total number of clarity. The update formula of Lucy-Richardson algorithm for I _m can be derived into formula (7), where I ^{t +1} represents the updated image in the iterative operation, and I ^t is the image before the update. Representing the Convolution Operation, and the Richardson-Lucy algorithm assumes that the image noise is based on Poisson Distribution, and the method of the present invention assumes that the image noise is a Gaussian distribution. (7) can be expressed using equation (8).

2.1, point spread function normalization

According to the above formula, to deconvolve the blurred image, the fuzzy length and angle obtained by the point spread function estimation can be regarded as a two-dimensional matrix, and in order for the Richardson-Lucy method to perform the folding operation smoothly, The formula (9) is used to normalize the values in the two-dimensional matrix between 0 and 1, and the sum is 1, where the denominator is the sum of the values of the elements in the matrix.

2.2 single image deconvolution

The fuzzy image is deconvolved by the formula (8). Since the Richardson-Lucy method iterative algorithm starts with I ^{t +1} unknown, the original blurred image is first initialized to I ^t , and the image is iterated through the above formula and restored to clear. The image, as shown in Figure 9, is the result of the restoration.

3, clear image register

When inputting the image sequence, it will first determine whether it is a blurred image. If the image is a clear image, it will be stored in the clear image buffer and output a clear image. If the image is blurred and the blur is large, then The single image will be deblurred first, and the result will be stored in the clear image register. The scratchpad will store K clear images or restored images, and will be shaved when a new screen enters. The oldest picture, in addition to before the image enters the scratchpad, there will be a counter for the cumulative action. When the value in the counter before entering the register exceeds K, it means in the current register. The pictures have not been restored by a single image, so the single image will be forced to be restored and then stored in the clear image register, thus preventing the accumulation of blurred picture errors.

4, multiple images to blur

In the actual image sequence, it usually contains clear images and blurred images. If each blurred image is deblurred by using a single image, it will be time consuming and cost-effective to be used for instant application, so the method of the present invention utilizes The concept of Temporal information to achieve instant video repair, mainly including image correction (Alignment) and video restoration. Firstly, the K image in the image buffer includes the original clear image and the restored image. Then, in order to obtain the data of the temporal neighboring domain, the feature points of the current image and the front image are searched for calculating the perspective transformation (Perspective Transformation). The Homography Matrix is used to obtain the transformation between the two images. After the perspective transformation, the pixels between the two images can be aligned and aligned; after the corrected image is obtained, the K -alignment is used. After the image image, the time neighborhood weighting calculation is performed, so that each pixel of the current image obtains a new calculation result.

4.1, image correction

The continuous image matching is corrected by the perspective transformation method. As shown in Fig. 10, a homography matrix H can be calculated through two consecutive images to perform image conversion. In Fig. 10, the t -1 image is H through the screen. The shooting angle becomes the same as the t-th sheet. The transformation processing may include panning, rotation, and zooming. Then, using Shi-Tomasi corners to detect the corner points in the image, the corner point is used as the feature point, and the most subsequent one should be The characteristic points of the matrix, as shown in Fig. 11, the perspective transformation of the image through the homography matrix, the main purpose of which is to correct the image in the clear image buffer to the same position as the current blurred image to be restored, so The image weight value calculation in the subsequent video repair processing can be performed.

In image imaging, basically the camera follows a perspective transformation model of small hole imaging, assuming that any point in the three-dimensional space x = [ X , Y , Z ] ^T , its corresponding point m = [ u , v ] ^T , where The homogeneous coordinates are x' = [ X , Y , Z , 1 ] ^T and m ' = [ u , v , 1] ^T . According to the perspective transformation model, the three-dimensional space point x and its corresponding point m should satisfy the formula (10). Where λ is a constant, f _x and f _y are the camera focal lengths, s is the inclination coefficient of the corresponding point coordinates, depending on the angle between the X and Y axes of the corresponding point coordinate system, R is a 3×3 rotation matrix, and T is 3 A translation vector of ×1. If the XY plane of the three-dimensional coordinate system is overlapped with the two-dimensional plane, the Z-value of the three-dimensional coordinates can be represented as 0. If the ith column element of the rotation matrix R is represented by r _i , the equation (10) can be rewritten as the equation (11). Where x' represents the homogeneous coordinate of any point x on the two-dimensional plane, then x' = [ X , Y , 1] ^T , so that equation (12) can be generated, where H is between the three-dimensional plane and the corresponding two-dimensional plane The homography matrix, because of the λ constant, transforms the mapping between 3D and 2D into a set of homography matrices. Four or more known feature points are required to obtain the homography matrix.

As described above, when establishing a homography matrix, at least four sets of feature points are required, and each set of feature points ( X, Y ) and its corresponding feature points ( u, v ) are brought into the formula (12) λm '= Hx ' Four sets of formulas (13) are generated, wherein each set of corresponding points can provide a linear equation of two homography matrices H , and four sets of formulas (13) are used to calculate the values of h ₁₁ ~ h _{33 in} the perspective transformation matrix H. Find the perspective transformation matrix H. If the parameter λ is set to 0, the above function can calculate its homography using the least square method (Least Squares). However, not all points can adapt to this transmission transformation. If the outlier is calculated, it will cause a great deviation in its homography estimate. Therefore, the RANSAC (RANdom SAmple Consensus) method is used to deal with this situation, and a set of four points is randomly selected in the feature point group in the picture, and the set and the least square method are used to estimate the homography matrix H , by minimizing the error. The best function matching of the square sum finding data satisfies the best function matching of the current transform matrix, that is, the sum of the squares of the error between the data and the actual data is the smallest, and returns the consistent concentrated data, and finally calculates the mass ratio of the homography matrix. , the best set of the set as the value of the homography matrix. As shown in FIG. 12 by the perspective transform result, the figure for the first t H ^{t -2} -2 screen homography matrix, H ^{t -1} t -1 for the first picture homography, t-visible After the -2 picture and the t -1th picture are calculated via their respective homography matrices, the picture can be rotated to the same level as the t-th picture.

4.2, video repair

Image correction is more convenient for pixel-to-pixel matching processing, which makes it easier to match the values of pixels. This makes multi-image restoration more efficient. By matching pixels to pixels, you can compare the weight values between pixels. The weight value is calculated as formula (14), l _{i, x} is to find the pixel point to replace the current blurred image, use the CK clear image in the image register to calculate the weight value, and finally select the pixel with the highest weight value. To replace the blurred image pixel value, first calculate the weighted average of the surrounding pixels for each sequential pixel q of the image CB _{ti,q in} the register, as shown in equation (15), where I _t,p is to be replaced. Currently, the restored pixel value of the blurred picture pixel p , t represents the time serial number of the picture frame, where the calculation is performed on the CK picture in the temporary register, where W is the sum of the weights w ( t, p, ti, q ), and w (t, p, ti, q ) weight value is calculated as shown in formula _{(16), H (CB ti} , q) of CB _{ti, q} via fluoroscopy image after the H process transitions, the current through the calculation of H (CB _{ti, q)} Fuzzy frame (B _{t, p)} for the Euclidean distance between the exponential function calculation, so it is Weight values normalized between [0,1], if this both _{H (CB ti, q),} (B t, p) the farther the distance, the exponentiation result is close to zero, which means that this The point does not have any reference significance. Here, the constant σ = 0.5 is set, which mainly adjusts the Euclidean distance calculated by H ( CB _{ti, q} ) and ( B _t,p ) to make the distance larger. After the exponential operation approaches 0, the pixel value with the largest weight value is selected to replace the current blurred picture pixel value. As shown in Fig. 13, after the image restoration processing, the left part of the figure is the original blurred image, and the right part is the restored result, and the edge portion of the content of the restored image becomes clearer.

As described above, the method for restoring a blurred image of the present invention can be applied to the restoration of a linear motion blur picture, and the proposed video deblurring method can be implemented in a common image monitoring system or a general moving vehicle. The device can perform the function of instant deblurring through video streaming. Since the algorithm proposed in this case has a small computational complexity, it can be embedded in a general security surveillance photography system to enhance the added value of related products.

The embodiments described above are merely illustrative of the technical spirit and the features of the present invention, and the objects of the present invention can be understood by those skilled in the art, and the scope of the present invention cannot be limited thereto. That is, the equivalent variations or modifications made by the spirit of the present invention should still be included in the scope of the present invention.

Claims (5)

A method for restoring a blurred image includes the following steps: performing Fourier transform on a blurred image in an image sequence to generate a converted value, then obtaining a inverse spectrum of the converted value, and calculating a point spread function from the cepstral domain Fuzzy angle and fuzzy length; normalizing the point spread function according to the blur angle and the blur length, and then performing a deconvolution of the single image on the normalized point spread function to generate a restored image; selecting the image sequence Clear image or the restored image is stored in the clear image register; and the corresponding feature point between the current image and the restored image stored in the clear image register is used, and the feature point is used to derive the requirement for the perspective transformation And correcting the current image and the corresponding pixel of the restored image according to the time domain information, and calculating the weights of the pixels in the current image and the restored image according to the time domain information, and comparing the low weight pixels with the high weight pixels Generating a restored pixel and outputting another restored image; wherein the current image is temporarily stored in the clear image The clarity of the image or the restored image is next to be treated of the blurred image. The method for restoring a blurred image according to the first aspect of the invention, wherein the clear image register stores the clear image in the image sequence when the image sequence is input, and the blurred image in the image sequence The restored image is restored by the restored method of the blurred image. The method for restoring a blurred image according to claim 1, wherein the point spread function is expressed by the following formula: Among them, L is the blur length, and θ is the blur angle. The method for restoring a blurred image according to claim 3, wherein the cepstrum domain definition of the blurred image is expressed by the following formula: C ( g ( x , y ))= F ^-1 {log| F ( g ( x , y ))|} where F represents a Fourier transform, F ^-1 represents an inverse Fourier transform, and g(x, y) is a blurred image. The method for restoring a blurred image according to claim 4, wherein the calculating the blur length of the point spread function and the blur angle comprises the following steps: (1) inputting a cepstrum map to find a center point; (2) Selecting the point of the largest cepstral value from the three directions of the center point of 90°, 45°, and 0° as the search point; (3) moving to the search point and using a 5*5 cover a cover, the center of the mask corresponds to the search point; (4) all the cepstrum values in the mask are added and the coordinates of the search point are stored; (5) determining whether the direction of the search point is offset, if otherwise entering Step (6), if yes, proceed to step (7); (6) calculate whether the distance between the search point and the center point is greater than a threshold, if otherwise, return to the step (2) to repeat the step (2) to the step ( 5), if yes, enter the step (7); (7) find the turning point of the falling and rising from the trajectory value distribution map; and (8) calculate the blur length and the blur angle according to the turning point and the center point; the turning point of coordinates (x _1, y _1), the center point of the coordinates _{(x 0, y 0),} L for the length fuzzy Substituting (x _1, y ₁₎ and (x _0, y ₀₎ to obtain the blur length lines represented by the following formula: R is the length of the base of the triangle formed by the turning point, the center point and the X axis, and substituting x ₁ and x ₀ to obtain the length of the base is expressed by the following formula: R =| x ₁ - x ₀ |θ The blur angle is substituted into the blur length and the length of the base to obtain the blur angle expressed by the following formula: