An alternative calculation method for the product of image fuzzy matrix and image vector under replicated boundary conditions

Abstract

This patent discloses an alternative calculation method for the product of image blur matrix and image vector under Replicated boundary conditions, including: 1) under Replicated boundary conditions, the product of image blur matrix and image vector (product one), and the product of transpose of image blur matrix and image vector (product two) are respectively converted into the sum of products of multiple block matrices with available structures and image vectors; 2) for each block matrix, a corresponding point spread function is constructed; 3) the convolution between the point spread function of the block matrix and the image or the corresponding boundary of the image is used to calculate each decomposed part in the matrix form of the two products instead; and 4) the two products are calculated. The disclosed calculation method can be integrated into large-scale image filtering and image restoration under Replicated boundary conditions to solve the problem that the product of image blur matrix and its transpose and image vector is difficult to calculate.

Description

Technical Field

[0001] The invention relates to image processing, and in particular to an alternative calculation of the product of an image blur matrix and its transposition with an image vector in linear image restoration under Replicated boundary conditions. Background Art

[0002] Image restoration is to solve the following degradation problem caused by linear shift invariance

[0003] f=Hx+η (1)

[0004] Where: Image blur matrix H∈R mn×mn By point spread function (PSF) h∈R p ×q The structure of the system is uniquely determined by the boundary conditions adopted; f∈R mn×1 From the degraded image F∈R m×n The vector obtained by row-first arrangement; η is the rearranged vector of additive noise.

[0005] In solving (1), the main computational task is to calculate the product of H and f (product 1), i.e. g1 = Hf, and the adjoint matrix of H * (In the real number domain, this is the transposed matrix H of H T ) and f (product 2), that is, g2 = H T f.

[0006] Generally, H and H T Both are extremely sparse and large matrices, making it difficult to directly calculate product 1 and product 2. In this regard, there are currently two alternative calculation methods for product 1 and product 2 under some boundary conditions: one is an alternative calculation method based on Kronecker products, and the other is an alternative calculation method based on diagonalizing the image blur matrix. Under deterministic boundary conditions, whether product 1 and product 2 can be calculated alternatively, and if so, which of the above methods is more effective, both depend on what boundary conditions are used and whether the point spread function is symmetric.

[0007] At present, the deterministic image boundary conditions include Zero boundary condition, Periodic boundary condition, Neumann boundary condition, Anti-Reflective boundary condition, Mean boundary condition and Replicated boundary condition. Under these conditions, the image blur matrix has a deterministic structure that is unrelated to the degraded image, and the structure is uniquely determined by the boundary conditions. This allows us to analyze these structures and apply them to the calculation of product-sum product and the design of the preset for formula (1).

[0008] The Replicated boundary condition is a new boundary condition, which was proposed by Chinese scholars Xu Zhou et al. in the paper "A boundary condition based deconvolution framework for image deblurring" in 2014. Xu Zhou et al. believe that this boundary condition takes into account the situation that the image pixels conform to the "heavy-tailed" distribution (that is, most pixels in the image have similar values ​​to their neighboring pixels). Therefore, it is a boundary condition that is more in line with the actual image and can be used in the restoration and processing of images with corresponding characteristics.

[0009] Xu Zhou et al. also proposed an image deblurring framework based on Replicated boundaries in the same paper: 1) Perform multiplicative decomposition on H, that is, convert H into a non-square biased blur matrix K∈R (m+2p1)(n+2q1)×mn With the non-square boundary condition matrix P∈R mn×(m+2p1)(n+2q1) The product of, that is, H = KP, thus transforming the image degradation problem f = Hx + ζ into f = KPx + ζ; 2) using the conjugate gradient algorithm to solve the linear equation P * K * y=P * K * KPx, to recover the unknown image x, and in the corresponding iterative solution, use two Fast Fourier Transformations (FFT) and one Inverse Fast Fourier Transformation to calculate the matrix form of product one and product two; 3) Determine the Replicated boundary conditions P and P * The structure of Px and P * x can be calculated, thus proposing an alternative calculation method for product 1 and product 2 under the Replicated boundary condition, which is also the only alternative calculation method under this boundary condition. The algorithm complexity of Xu Zhou et al.'s alternative calculation method for product 1 and product 2 is O(m×n×Nnz), where Nnz is the number of nonzero elements in the point spread function. Under other boundary conditions, whether the framework is feasible depends on the boundary matrices P and P *has been determined. Summary of the invention

[0010] The technical problem to be solved by the present invention is to design a new indirect calculation method for product one and product two under the Replicated boundary condition. This method does not perform multiplication decomposition of the fuzzy matrix, but performs addition decomposition, so that the calculation accuracy of product one and product two is higher, and when the point spread function is of a regular size, the time consumption is equivalent to that of Xu Zhou's method.

[0011] In view of the above problems, the present invention provides an indirect calculation method for product one and product two under Replicated boundary conditions. The idea is: under Replicated boundary conditions, according to the structure of the image fuzzy matrix, first decompose product one and product two into the sum of multiple matrix and vector products, respectively, to ensure that the generated decomposition matrix has an available block structure; construct the convolution kernel of each block matrix by using the structure of each block matrix; use convolution instead to calculate the matrix form of the product of each decomposition matrix and the image vector; and obtain product one and product two by combining the calculation results.

[0012] To support the above ideas, in the present invention: Example 1 gives the specific steps of the alternative calculation method of product one and product two under the Replicated boundary condition; Example 2 converts product one and product two under the Replicated boundary condition into the sum of the products of multiple fuzzy matrices and the image vector respectively; Example 3 constructs the corresponding point spread for each decomposition matrix in the decomposition formula of product one and product two; Example 4 uses convolution based on the constructed point spread function to calculate the matrix form of the product of each decomposition matrix of product one and product two and the image vector; Example 5 illustrates the effectiveness of the convolution alternative method given by the present invention through experimental results.

[0013] The calculation method provided by the present invention provides an alternative calculation for large product one and product two under Replicated boundary conditions, which can be integrated into image filtering and image restoration under Replicated boundary conditions and has application value. The calculation accuracy of product one and product two provided by the present invention is relatively high. The time complexity of this algorithm mainly comes from the complexity of the convolution operation of product one and product two under zero boundary conditions, that is, O(m×n×p×q). Therefore, when the point spread function is a normal case, that is, a normal size and the element values ​​are random, the time consumption of this method is comparable to that of Xu Zhou's method, and when the point spread size is large and the number of non-zero elements is small, the patented method will be slower than Xu Zhou's method. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 The calculation process of the present invention is shown;

[0015] Figure 2 Shown is the grayscale image Saturno image used in the experiment, with a size of 512×512;

[0016] Figure 3 Shown is the grayscale image Cameraman image used in the experiment, with a size of 256×256;

[0017] Figure 4 The point spread function Saturnopsf used in the experiment is shown.

[0018] Figure 5 , Figure 6 and Figure 7 Shown are three randomly generated point spread functions PSFn1, PSFn2 and PSFn3 used in the image Cameraman in the experiment, with sizes of 53×53, 7×7, 17×13 and 31×37 respectively. DETAILED DESCRIPTION

[0019] In the following description, for the purpose of patent description, the specific contents of the embodiments described in the invention content are described in combination with the accompanying drawings. In order to express clearly and facilitate programming implementation, Matlab function forms are applied in some calculation formulas.

[0020] like Figure 1 As shown, the present invention gives the product g1 = Hf and the product g2 = H under the Replicated boundary condition. T An alternative method for calculating f includes the following steps:

[0021] (1) Input F and h, and rearrange F into an image vector f in row-first order;

[0022] (2) Decompose g1 and g2;

[0023] (3) constructing the point spread function of the fuzzy matrix involved in each decomposition part in the decomposition formula of g1 and g2, and based on this, convolving the point spread function with the image or the corresponding image boundary to calculate the matrix form of each decomposition part;

[0024] (4) Calculate the corresponding images of g1 and g2;

[0025] (5) The corresponding images of g1 and g2 are rearranged into vectors by row priority, and the obtained images are g1 and g2. g 2.

[0026] Second Embodiment

[0027] The present invention decomposes g1 into g1 = Hf = (H TT +H TL +HLT +H LL )f, the decomposition of g2 is g2=H T f=(H T TT +H T TL +H T LT +H T LL )f, in the two decomposition formulas: in each matrix, the superscript T represents the transpose of the matrix, the first subscript represents the inter-block structure of the current matrix, and the second subscript represents the intra-block structure of the current matrix. T and L represent the Toeplitz matrix and L-type matrix respectively. The specific structures of the T and L matrices are not given here, which will not affect the subsequent description of the calculation method.

[0028] Embodiment 3

[0029] The present invention provides a product-decomposition formula g1=Hf=(H TT +H TL +H LT +H LL )f and the product binary decomposition formula g2=H T f=(H T TT +H T TL +H T LT +H T LL )f, the construction method of the point spread function of each matrix is ​​as follows:

[0030] Let fliplr(.) represent the flipping of the matrix column direction, i.e. left-right, flipud(.) represent the flipping of the matrix row direction, i.e. up-down, conv2(.,.) represent the two-dimensional convolution, i:j represent all subscripts from i to j, V(:,:) in: indexes all rows and columns of the matrix V, and V(end,end) in: indexes the last row and last column of the matrix V; in addition, let (Note: We assume that the number of rows p and the number of columns q of the point spread function are both odd numbers, that is, p=2p1+1, q=2q1+1. Even if the size of the actual point spread function is not an odd number, it can be made into an odd number by appending one and only one zero row or / and one and only one zero column to the end of the row or column of the point spread function. The corresponding modified point spread function will not affect the calculation results of this patent). Then:

[0031] (1)H TT The point spread function is h1=h;

[0032] (2) Construction H TLThe point spread function has 2 kernels, corresponding to the non-zero elements in the first and last columns of each L-shaped sub-block. The specific method is and where i=1,...,q1 j=1,...,p;

[0033] (3) Construction H LT Point spread function, with 2 kernels, corresponding to H LT There are two non-zero column blocks at the beginning and end of and Where 1≤i≤p1;

[0034] (4) Construction H LL Point spread function, a total of 4 kernels, corresponding to H LL The first and last non-zero columns in the first and last two non-zero columns block (a total of 4 non-zero columns), the corresponding calculation formula is the same as H TL The core is related to and

[0035] (5)H T TT The point spread function is h2 = fliplr (flipud (h));

[0036] (6) Construction H T TL Point spread function, with 2 kernels, corresponding to the methods h 2_TL_1 =fliplr((h 1_TL_1 ) T ) and h 2_TL_2 =fliplr((h 1_TL_2 ) T );

[0037] (7) Construction H T LT The point spread function of the 2-kernel method is h 2_LT_1 =fliplr(flipud(h 1_LT_1 )) and h 2_LT_2 =fliplr(flipud(h 1_LT_2 ));

[0038] (8) Construction H T LL Point spread function, a total of 4 kernel methods are h 2_LL_1 =fliplr(flipud((h 1_LL_1 ) T )),h 2_LL_2 =fliplr(flipud((h 1_LL_2 ) T)),h 2_LL_3 =fliplr(flipud((h 1_LL_3 ) T )) and h 2_LL_4 =fliplr(flipud((h 1_LL_4 ) T )).

[0039] Fourth Embodiment

[0040] The present invention is based on convolution and provides a product-decomposition formula g1=Hf=(H TT +H TL +H LT +H LL )f and product g2 = H T f=(H T TT +H T TL +H T LT +H T LL )f, the calculation method of the product of each matrix and vector is as follows:

[0041] (1)H TT The matrix form G corresponding to f 1_TT , the size is m×n, and its calculation method is as follows: 1_TT =conv2(F,h,'same');

[0042] (2)H TL The matrix form of f is G 1_TL , the size is m×n, and the calculation method is: Calculate r 1_TL_1 =conv2(flipud((h 1_TL_1 ) T ),F(:,1)) and r 1_TL_2 =conv2(flipud((h 1_TL_2 ) T ),F(:,n)), then G 1_TL (:,1:q1)=r 1_TL_1 (p1+1:end-p1,:),G 1_TL (:,end-q1+1:end)=r 1_TL_2 (p1+1:end-p1,:),G 1_TL The pixel values ​​in other areas are 0;

[0043] (3)H LT The matrix form of f is G 1_LT , the size is m×n, and the calculation method is: Calculate r 1_LT_1 =conv2(h1_LT_1 ,F(1,:)) and r 1_LT_2 =conv2(h 1_LT_2 ,F(m,:)), then G 1_LT (1:p1,:)=r 1_LT_1 (:,q1+1:end-q1), G 1_LT (end-p1+1:end,:)=r 1_LT_2 (:,q1+1:end-q1), G 1_LT The pixel values ​​in other areas are 0;

[0044] (4)H LL The matrix form of f is G 1_LL , the size is m×n, and the calculation method is: Calculate r 1_LL_1 =conv2((h 1_LL_1 ) T ,F(1,1)),r 1_LL_2 =conv2((h 1_LL_2 ) T ,F(1,n)),r 1_LL_3 =conv2((h 1_LL_3 ) T ,F(m,1)) and r 1_LL_4 =conv2((h 1_LL_4 ) T ,F(m,n)), then G 1_LL (1:p1,1:q1)=r 1_LL_1 , G 1_LL (1:p1,end-q1+1:end)=r 1_LL_2 , G 1_LL (end-p1+1:end,1:q1)=r 1_LL_3 , G 1_LL (end-p1+1:end,end-q1+1:end)=r 1_LL_4 , G 1_LL The pixel values ​​in other areas are 0;

[0045] (5)H T TT The matrix form of f is G 2_TT , the size is m×n, and the calculation method is,G 2_TT =conv2(F,fliplr(flipud(h)),'same');

[0046] (6)H T TL The matrix form of f is G 2_TL , the size is m×n, and the calculation method is: Calculate r 2_TL_1 =conv2(h2_TL_1 ,F(:,1:q1)),r 2_TL_2 =conv2(h 2_TL_2 ,F(:,end-q1+1:end)), then G 2_TL (:,1)=r 2_TL_1 (p1+1:end-p1,q1), G 2_TL (:,n)=r 2_TL_2 (p1+1:end-p1,q1), G 2_TL The pixel values ​​in other areas are 0;

[0047] (7)H T LT The matrix form of f is G 2_LT , the size is m×n, and the calculation method is: Calculate r 2_LT_1 =conv2(h 2_LT_1 ,F(1:p1,:)),r 2_LT_2 =conv2(h 2_LT_2 ,F(end-p1+1:end,:)), then G 2_LT (1,:)=r 2_LT_1 (p1,q1+1:end-q1), G 2_LT (m,:)=r 2_LT_2 (p1,q1+1:end-q1), G 2_LT The pixel values ​​in other areas are 0;

[0048] (8)H T LL The matrix form of f is G 2_LL , the size is m×n, and the calculation method is: Calculate r 2_LL_1 =conv2(h 2_LL_1 ,F(1:p1,1:q1)),r 2_LL_2 =conv2(h 2_LL_2 ,F(1:p1,end-q1+1:end)),r 2_LL_3 =conv2(h 2_LL_3 ,F(end-p1+1:end,1:q1)),r 2_LL_4 =conv2(h 2_LL_4 ,F(end-p1+1:end,end-q1+1:end)), then G 2_LL (1,1)=r 2_LL_1 (p1,q1),G 2_LL (1,n)=r 2_LL_2 (p1,q1),G 2_LL (m,1)=r 2_LL_3 (p1,q1),G 2_LL(m,n)=r 2_LL_4 (p1,q1),G 2_LL The pixel values ​​in other areas are 0.

[0049] Fifth embodiment

[0050] Under the RPBC boundary condition, the matrix result of the product one is given as G1=G 1_TT +G 1_TL +G 1_LT +G 1_LL , the matrix result of product two is G2=G 2_TT +G 2_TL +G 2_LT +G 2_LL The calculation methods of the matrices on the right sides of the two equations have been given in Example 4.

[0051] Sixth Embodiment

[0052] The present invention uses an experiment to verify the effectiveness of the method given in the above embodiment from two aspects: calculation error and calculation time. The experiment is completed using Matlab2014a on a ThinkPad notebook with a CPU of 2.0GHz and a memory of 4G. The experimental materials are: the experimental image F includes a Saturno image with a size of 512×512 and a Cammerman image with a size of 256×256, respectively. Figure 2 and Figure 3 For experimental image 1, the point spread function h uses Saturnopsf with a size of 53×53, and its energy is mainly concentrated in the range of h(26:30,25:29), see Figure 4 ; For experimental image 2, the point spread function h uses three randomly generated matrices with sizes of 7×7, 17×13 and 31×37, see Figure 5 , Figure 6 and Figure 7 .

[0053] Tables 1 and 2 below compare the root mean square error (Rmse) and the time (in seconds) obtained by the method of the present invention and the method of Xuzhou et al. in the process of calculating product one and product two. The calculation method of Rmse is in It is the alternative calculation result of product one or product two using the method of this patent or XuZhou method, and g is the real product one or product two.

[0054] Table 1 Experimental results of product 1

[0055]

[0056] Table 2 Experimental results of product 2

[0057]

[0058] It can be seen from Table 1 and Table 2 above that in terms of calculation accuracy, the Rmse of product one and product two calculated by the method of the present invention is smaller and has higher accuracy.

[0059] The first three rows of Tables 1 and 2 above show that when the point spread function is the more conventional case represented by the first three, the time consumption of this method is comparable to that of Xu Zhou et al.; the last row of the two tables shows that when the point spread function becomes larger, especially when there are fewer non-zero elements in the point spread function, the time consumption of this method is greater than that of Xu Zhou et al., mainly because of the time consumption of G 1_TT and G 2_TT The convolution calculation.

Claims

1. A method for image processing under Replicated boundary conditions, for image F∈R m×n and point spread function h∈R p×q The method is based on the Replicated boundary condition, and the image blur matrix H∈R is composed of the elements of h. mn×mn The row-first rearranged image vector f∈R with F mn×1 The product of the two, that is, the product g1 = Hf, and the transpose of H T The product of g2 and f, that is, the product of g2 = H T f is calculated instead using the following steps: (1) Input F and h; (2) Perform additive decomposition on the fuzzy matrix and its transpose in the calculation formulas of g1 and g2. Specifically, decompose H into H TT , H TL , H LT and H LL , H T Decomposed into H T TT , H T TL , H T LT and H T LL In each matrix, the first subscript represents the inter-block structure of the matrix, the second subscript represents the intra-block structure of the matrix, the subscript T represents the Toeplitz matrix, the subscript L represents the L-type matrix, and the superscript T represents the transpose of the matrix; (3) Decompose g1 and g2. Specifically, 18 convolution kernels are constructed for each decomposed part, 16 of which are for H TL , H LT , H LL , H T TL , H T LT and H T LL The fuzzy matrix construction is only applicable under the Replicated boundary condition. These convolution kernels are used as point spread functions and convoluted with F, the four borders of F, and the four vertex pixels of F to calculate the non-zero area of ​​the matrix corresponding to each decomposed part in g1 and g2; (4) Calculate the corresponding images of g1 and g2; (5) The corresponding images of g1 and g2 are rearranged into vectors in row priority order to obtain g1 and g2.

2. The image processing method according to claim 1, characterized in that The pair g1 and g2 are decomposed, and the specific decomposition method is: The decomposition of g1 is g1=Hf=(H TT +H TL +H LT +H LL )f, the decomposition of g2 is g2=H T f=(H T TT +H T TL +H T LT +H T LL )f.

3. The image processing method according to claim 2, characterized in that For the matrices corresponding to the decomposition parts in the decomposition formulas of g1 and g2, the point spread function is constructed as follows: Use fliplr(X) to flip X left to right, flipud(X) to flip X upside down, and V(:,:) in : means to get all rows or columns of matrix V. and but: (1)H TT The point spread function is h1=h; (2) Construction H TL Point spread function, with 2 kernels, is and Where i=1,...,q1j=1,...,p; (3) Construction H LT Point spread function, with 2 kernels, is and Where 1≤i≤p1; (4) Construction H LL Point spread function, with 4 kernels, the method is and Where i=1,...,p1; (5)H T TT The point spread function is h2 = fliplr (flipud (h)); (6) Construction H T TL Point spread function, with 2 kernels, method h 2_TL_1 =fliplr((h 1_TL_1 ) T ) and h 2_TL_2 =fliplr((h 1_TL_2 ) T ); (7) Construction H T LT Point spread function, with 2 kernels, method h 2_LT_1 =fliplr(flipud(h 1_LT_1 )) and h 2_LT_2 =fliplr(flipud(h 1_LT_2 )); (8) Construction H T LL Point spread function, 4 kernels in total, method h 2_LL_1 =fliplr(flipud((h 1_LL_1 ) T )),h 2_LL_2 =fliplr(flipud((h 1_LL_2 ) T )),h 2_LL_3 =fliplr(flipud((h 1_LL_3 ) T )) and h 2_LL_4 =fliplr(flipud((h 1_LL_4 ) T )).

4. The image processing method according to claim 3, characterized in that For each decomposition part in the decomposition formula of g1 and g2, the specific calculation method in matrix form is: Use conv2(.,.) to represent a two-dimensional convolution, use i:j to represent all integers from i to j, use end in V(end,end) to index the last row and last column of the corresponding matrix V, and use G 1_TT , G 1_TL , G 1_LT , G 1_LL , G 2_TT , G 2_TL , G 2_LT and G 2_LL , respectively represent H TT f、H TL f、H LT f、H LL f、H T TT f、H T TL f、H T LT f and H T LL The matrix form of the matrix-vector product of f is: (1) G 1_TT The calculation method is G 1_TT =conv2(F,h1,'same'); (2) G 1_TL The calculation method is: Calculate r 1_TL_1 =conv2(flipud((h 1_TL_1 ) T ),F(:,1)) and r 1_TL_2 =conv2(flipud((h 1_TL_2 ) T ),F(:,n)), then G 1_TL (:,1:q1)=r 1_TL_1 (p1+1:end-p1,:),G 1_TL (:,end-q1+1:end)=r 1_TL_2 (p1+1:end-p1,:),G 1_TL The pixel values ​​in other areas are 0; (3)G 1_LT The calculation method is: Calculate r 1_LT_1 =conv2(h 1_LT_1 ,F(1,:)) and r 1_LT_2 =conv2(h 1_LT_2 ,F(m,:)), then G 1_LT (1:p1,:)=r 1_LT_1 (:,q1+1:end-q1), G 1_LT (end-p1+1:end,:)=r 1_LT_2 (:,q1+1:end-q1), G 1_LT The pixel values ​​in other areas are 0; (4)G 1_LL The calculation method is: Calculate r 1_LL_1 =conv2((h 1_LL_1 ) T ,F(1,1)),r 1_LL_2 =conv2((h 1_LL_2 ) T ,F(1,n)),r 1_LL_3 =conv2((h 1_LL_3 ) T ,F(m,1)) and r 1_LL_4 =conv2((h 1_LL_4 ) T ,F(m,n)), then G 1_LL (1:p1,1:q1)=r 1_LL_1 , G 1_LL (1:p1,end-q1+1:end)=r 1_LL_2 , G 1_LL (end-p1+1:end,1:q1)=r 1_LL_3 , G 1_LL (end-p1+1:end,end-q1+1:end)=r 1_LL_4 , G 1_LL The pixel values ​​in other areas are 0; (5)G 2_TT The calculation method is G 2_TT =conv2(F,h2,'same'); (6)G 2_TL The calculation method is: Calculate r 2_TL_1 =conv2(h 2_TL_1 ,F(:,1:q1)),r 2_TL_2 =conv2(h 2_TL_2 ,F(:,end-q1+1:end)), then G 2_TL (:,1)=r 2_TL_1 (p1+1:end-p1,q1), G 2_TL (:,n)=r 2_TL_2 (p1+1:end-p1,q1), G 2_TL The pixel values ​​in other areas are 0; (7)G 2_LT The calculation method is: Calculate r 2_LT_1 =conv2(h 2_LT_1 ,F(1:p1,:)),r 2_LT_2 =conv2(h 2_LT_2 ,F(end-p1+1:end,:)), then G 2_LT (1,:)=r 2_LT_1 (p1,q1+1:end-q1), G 2_LT (m,:)=r 2_LT_2 (p1,q1+1:end-q1), G 2_LT The pixel values ​​in other areas are 0; (8)G 2_LL The calculation method is: Calculate r 2_LL_1 =conv2(h 2_LL_1 ,F(1:p1,1:q1)),r 2_LL_2 =conv2(h 2_LL_2 ,F(1:p1,end-q1+1:end)),r 2_LL_3 =conv2(h 2_LL_3 ,F(end-p1+1:end,1:q1)),r 2_LL_4 =conv2(h 2_LL_4 ,F(end-p1+1:end,end-q1+1:end)), then G 2_LL (1,1)=r 2_LL_1 (p1,q1),G 2_LL (1,n)=r 2_LL_2 (p1,q1),G 2_LL (m,1)=r 2_LL_3 (p1,q1),G 2_LL (m,n)=r 2_LL_4 (p1,q1),G 2_LL The pixel values ​​in other areas are 0.

5. The image processing method according to claim 4, characterized in that Calculate the corresponding images G1 and G2 of g1 and g2. The specific calculation method is: Calculate the corresponding image G1 of g1 = G 1_TT +G 1_TL +G 1_LT +G 1_LL The corresponding image of g2 is G2=G 2_TT +G 2_TL +G 2_LT +G 2_LL .

Images