A small-size image-based source camera identification method

Abstract

The application discloses a source camera identification method based on small-size images, comprising the following steps: firstly, the image is divided into blocks, the real boundary is extended, then the double-tree complex wavelet is combined with the adaptive windowed wiener filter to obtain the noise residual of the image block, and the weight adjustment strategy is formulated according to the noise variance of the residual to obtain the weight matrix, and the reference PRNU noise of multiple cameras is calculated by combining the maximum likelihood estimation method; the test image is denoised by the double-tree complex wavelet combined with the adaptive windowed wiener filter, and the PRNU noise of the test image is extracted; the correlation between the PRNU noise of the to-be-tested image and the reference PRNU noise of all cameras is calculated, when the correlation with a certain reference camera is the highest, it is determined that the to-be-tested image comes from the camera. The application solves the problem that the noise interference of different regions of the small-size image is not fully considered in the existing camera PRNU estimation algorithm, and the variance weight is abnormal due to the traditional weighted average method, and effectively improves the small image source camera identification performance.

Description

Technical Field

[0001] This invention relates to the field of image forensics, and in particular to a source camera identification method based on small-sized images. Background Technology

[0002] Source detection using photoresponse non-uniform noise (PRNU) can identify not only different camera models but also different cameras of the same model. PRNU primarily arises from defects in the sensor manufacturing process and the varying light sensitivities of pixels due to silicon wafer inhomogeneities. Therefore, even sensors manufactured from the same silicon wafer will exhibit uncorrelated non-uniform noise. This characteristic makes PRNU a powerful fingerprint for identifying source devices and verifying image integrity.

[0003] Since the accuracy of source device recognition for large-sized images is already very high, if the accuracy of source device recognition for small-sized images can meet the actual needs, the required hardware storage space will be smaller, the algorithm's time efficiency will be higher, and it may be possible to realize source device recognition tasks on micro devices.

[0004] Zeng H, Wan Y, Deng K, et al. Source camera identification with Dual-Tree complex wavelet transform[J].IEEE Access,2020,8:18874-18883. This paper proposes using dual-tree complex wavelet transform and leveraging the local dependence of wavelet coefficients to shrink the wavelet coefficients, thereby further extracting PRNU from the image. Lawgaly A, Khelifi F. Sensor pattern noise estimation based on improved locally adaptive DCT filtering and weighted averaging for source camera identification and verification[J].IEEE Transactions on InformationForensics and Security,2016,12(2):392-404.(Lawgaly A, Khelifi F. Sensor pattern noise estimation based on improved locally adaptive DCT filtering and weighted averaging for source camera identification and verification[J].IEEE Transactions on InformationForensics and Security,2016,12(2):392-404.) proposes an improved locally adaptive discrete cosine transform filter and incorporates a weighted averaging technique that weights the noise based on the variance characteristics of the entire image.

[0005] Existing weighted averaging techniques fail to achieve minimum mean squared error because the degree of interference varies across different regions of the image, and simple weighting based on the variance of the noise residual does not yield reasonable weights. To address these issues, a more accurate method for estimating the camera reference unit (PRNU) is needed. Summary of the Invention

[0006] To address the technical problem this invention aims to solve, a source camera recognition method is provided that can achieve superior performance on small-sized images. This method can estimate more accurate light response non-uniform noise from small-sized images, effectively improving the accuracy of image source camera recognition, and also has strong robustness in recognizing highly compressed images.

[0007] To achieve the above objectives, the technical solution provided by this invention is as follows:

[0008] A source camera recognition method based on small-sized images includes:

[0009] S1, obtain the reference PRNU noise from multiple cameras and store these PRNU noises in the PRNU noise database;

[0010] S2, the test image is denoised using a dual-tree complex wavelet filter, and the PRNU noise of the test image is extracted by subtracting the denoised image from the test image.

[0011] S3. Calculate the correlation between the PRNU noise of the image under test and the reference PRNU noise of all cameras in the PRNU noise database. When the correlation with a certain reference camera is the highest, it is determined that the image under test comes from that camera.

[0012] Furthermore, in step S1, the specific steps for obtaining the reference PRNU noise of a certain camera are as follows:

[0013] A1. A set of images is captured by the camera, and a set of small-sized images are obtained by cropping the central part of the set of images as reference images for the camera; each reference image is divided into multiple image blocks, and the content of the reference image neighborhood is extended to the boundary of the image blocks;

[0014] A2, after processing the image blocks in step A1, apply dual-tree complex wavelet transform to each image block and then use a median estimator to estimate the noise variance of the wavelet subband. Then, apply the noise variance to windowed Wiener filtering to estimate the wavelet coefficients. Next, perform dual-tree complex wavelet inverse transform, then restore the image block boundary expansion and subtract the image block before denoising from the image block after denoising to obtain the noise residual of each image block.

[0015] A3: Weights are set according to the variance of interference noise on the image patch using a weighted average method, and the weights are adjusted for image patches with abnormal variance using a weight adjustment strategy; after performing weight adjustment operations on all image patches, a weight coefficient matrix is ​​obtained.

[0016] A4. The PRNU multiplicative factor K of the camera is obtained by combining the weight coefficient matrix with the weighted maximum likelihood estimation method.

[0017] A5. The multiplicative factor K is zero-mean processed, and then the zero-mean processed multiplicative factor K is multiplied with all reference images captured by the camera and then averaged to obtain the reference PRNU noise of the camera.

[0018] Furthermore, the specific process of using the median estimator to estimate the noise variance of the wavelet subband in step A2, and then applying the noise variance to windowed Wiener filtering to estimate the wavelet coefficients, is as follows:

[0019] After performing a two-dimensional dual-tree complex wavelet transform on each image block to obtain wavelet coefficients, the filtering estimation process can be expressed as equation (1):

[0020]

[0021] In the formula, W in W represents the wavelet coefficients before filtering. out (u,v) represents the wavelet coefficients at the filtered pixel (u,v), and the noise variance. and We can obtain the following from equations (2) and (3):

[0022]

[0023]

[0024] In equation (2), median() represents the median estimator, and N in equation (3) h It is a local window of size h×h centered at (u,v), where (i,j) represents a point within the local window, W temp This represents the wavelet coefficients of the high-frequency subband after the image undergoes dual-tree complex wavelet transform. Max() represents taking the maximum value between 0 and the variance estimate, while the min() function represents taking the minimum value of all window estimates.

[0025] Further, the step of setting weights based on the variance of interference noise on the image patch using a weighted average method includes:

[0026] Variance of image patch noise The calculation method is shown in equation (4), taking the image block in row p and column q of the i-th reference image as an example:

[0027]

[0028] P and Q represent the number of rows and columns of an image block into which a reference image is cut, respectively; This represents the number of pixels in an image patch that are rearranged in one direction. Indicates interference noise on image patches The components at each pixel As shown in equation (5), Then it means The mean;

[0029]

[0030] In equation (5), X i (p,q) represents the noise residual of the image patch. For X iThe mean of (p,q), where N is the number of reference images;

[0031] After obtaining the variance of the image patch interference noise, the weights of the image patches are obtained using the weighted average method of equation (6):

[0032]

[0033] In the formula, w k (p,q) represents the weight of the image patch at position (p,q) in the k-th reference image. Let represent the variance of the interference noise of the image patch at position (p,q) in the k-th image, and let i represent the i-th reference image.

[0034] Furthermore, the image blocks with aberrations are weighted using a weight adjustment strategy, including:

[0035] When the variance of the interference noise is greater than λ1 or less than λ2, the variance is adjusted. Specifically, σ in formula (6) is adjusted. ο 2 (p,q) is adjusted to be calculated as shown in equation (7). Where ο = i, k:

[0036]

[0037] In equation (7), λ1 and λ2 represent the upper and lower thresholds of the noise variance, respectively; r is the sensitivity factor. This represents the variance of the interference noise after being constrained using a weighting adjustment strategy.

[0038] Furthermore, the method of obtaining the PRNU multiplicative factor K of the camera using the weighted coefficient matrix combined with the weighted maximum likelihood estimation method includes:

[0039] The noise variance weighted maximum likelihood estimation algorithm is expressed as equation (8):

[0040]

[0041] X i Let I represent the noise residual of the i-th reference image. i Let w represent the i-th reference image. i This represents the weight coefficient matrix consisting of the weights of all image patches in the i-th reference image, where N is the number of reference images.

[0042] Furthermore, the zero-mean normalization process is as follows:

[0043] Subtract the column mean from each value in the column of the multiplicative factor K, and then subtract the row mean from each value in the row.

[0044] Furthermore, the calculation of the correlation between the PRNU noise of the image under test and the reference PRNU noise of all cameras in the PRNU noise database using peak correlation energy includes:

[0045] Peak-related energy PCE is expressed as:

[0046]

[0047] In the formula, MN is the product of dimensions of the non-uniform noise in the matched optical response, and NCC is... RQ (i,j) represents the two-dimensional cyclic cross-correlation between the camera's reference PRNU noise R and the test image's PRNU noise Q at point (i,j):

[0048]

[0049] Where m and n represent the two dimensions of R and Q, R(i,j) represents the value of R at point (i,j), and Q(k+i,l+j) represents the value of Q at point (k+i,l+j). The average value of R is... This represents the average value of Q.

[0050] Compared with the prior art, the present invention has the following technical features:

[0051] 1. True boundary extension solves the problem of boundary effects when performing wavelet decomposition-based denoising on image blocks. Furthermore, by extending the original image content, it avoids introducing obvious discontinuities along the image boundary after the extension, thus preventing the denoising effect from being affected.

[0052] 2. Weighted averaging techniques can provide estimates that are very close to the actual signal, but because the quality of PRNU extracted from different regions of the image varies, calculating weights based on the entire image is not ideal. Therefore, this invention obtains weights by analyzing block-level images.

[0053] 3. Using the variance of the interference noise instead of directly using the variance of the noise residual for weighting will result in a more accurate estimation of the camera reference PRNU (i.e., camera light response non-uniform noise).

[0054] 4. Considering the problem of unreasonable weighting when the variance of image patches is abnormal in previous methods, such as when the variance of interference noise in some image patches is close to 0, the weight will increase sharply as the variance decreases. This may cause the device reference PRNU in a certain area to be dominated by a small number of image patches, while the remaining image patches become meaningless in calculating the device reference PRNU; conversely, the role of image patches with excessively high variance in calculating the reference PRNU will be ignored. To address this shortcoming, a weight adjustment strategy is proposed for improvement.

[0055] 5. Using dual-tree complex wavelet transform can effectively separate high- and low-frequency information and can effectively extract PRNU noise from high-frequency coefficients. Compared with wavelet transform, dual-tree complex wavelet transform has good shift invariance and good direction selectivity. This method shows a significant performance improvement over DWT-based methods in regions around strong edges.

[0056] 6. The method of this invention solves the problems of insufficient consideration of noise interference in different regions of small-sized images and abnormal variance weights caused by traditional weighted averaging methods in existing camera PRNU estimation algorithms. It also proposes block-level true boundary extension to eliminate the inherent defects of filters. Therefore, in image source recognition, especially in the recognition of small-sized image sources, the method of this invention has higher recognition performance than the current state-of-the-art algorithms. Attached Figure Description

[0057] Figure 1 This is a flowchart illustrating the principle of the source camera tracing method based on small-sized images of the present invention.

[0058] Figure 2 The flowchart shows how to obtain the PRNU multiplicative factor K by combining the maximum likelihood estimation of the weight coefficient matrix. Detailed Implementation

[0059] This embodiment provides a source camera identification method based on small-sized images. First, the image is segmented and its true boundaries are expanded. Then, multi-scale dual-tree complex wavelet transform is applied to the image patches, followed by local adaptive threshold window Wiener filtering. After inverse transform and restoration of the expanded image, noise residuals are obtained. A weight adjustment strategy is formulated based on the noise variance of the residuals to obtain a weight matrix that reasonably represents the importance of image patches. Based on this weight matrix, maximum likelihood estimation is used to obtain the camera's photo response non-uniformity (PRNU), and a zero-mean method is used to suppress non-uniqueness noise. Finally, the correlation between the camera's PRNU noise and the PRNU of the natural image is calculated using the correlation energy peak, and source identification and detection are performed based on this correlation. This invention solves the problems of insufficient consideration of noise interference in different regions of small-sized images and the abnormal variance weights caused by traditional weighted averaging methods in existing camera PRNU estimation algorithms, effectively improving the performance of source camera identification for small images.

[0060] See Figure 1 The specific implementation steps of this invention are as follows:

[0061] S1: Calculate the reference PRNU noise for multiple cameras and store these PRNU noises in the PRNU noise database.

[0062] In step S1, the specific steps for obtaining the reference PRNU noise of a certain camera are as follows:

[0063] A1. A set of images is captured by the camera, and a set of small-sized images is obtained by cropping the central portion of the set of images as reference images for the camera. Each reference image is divided into image blocks of 16×16 pixels or 32×32 pixels, and the content of the neighborhood of the reference image is extended to the boundary of the image block. In this embodiment, the small-sized images refer to images of 512×512 pixels or smaller.

[0064] This step can improve the accuracy of wavelet coefficient assignment in wavelet transform-based denoising filters and eliminate boundary effects caused by filter defects.

[0065] A2. After processing each small image block in step A1, perform dual-tree complex wavelet transform on each image block and then use a median estimator to estimate the noise variance of the wavelet subband. Then, apply the noise variance to windowed Wiener filtering to estimate the wavelet coefficients. Next, perform dual-tree complex wavelet inverse transform to achieve the effect of denoising the image block. After that, restore the image block boundary extension and subtract the image block before denoising from the image block after denoising to obtain the noise residual of each image block.

[0066] Specifically, the process of using a median estimator to estimate the noise variance of the wavelet subband in step A2, and then applying the noise variance to windowed Wiener filtering to estimate the wavelet coefficients, is as follows:

[0067] After performing a two-dimensional dual-tree complex wavelet transform on each image block to obtain wavelet coefficients, the filtering estimation process can be expressed as equation (1):

[0068]

[0069] In the formula, W in W represents the wavelet coefficients before filtering. out (u,v) represents the wavelet coefficients at the filtered pixel (u,v), and the noise variance. and We can obtain the following from equations (2) and (3):

[0070]

[0071]

[0072] In equation (2), median() represents the median estimator, and N in equation (3) h It is a local window of size h×h centered at (u,v), where (i,j) represents a point within the local window, W tempThis represents the wavelet coefficients of the high-frequency subband after the image undergoes a dual-tree complex wavelet transform. Max() indicates taking the maximum value between 0 and the variance estimate, while min() indicates taking the minimum value of all window estimates. Then, the estimated wavelet coefficients are subjected to a two-dimensional dual-tree complex wavelet inverse transform to obtain the denoised image patch.

[0073] A3 sets weights based on the variance of interference noise on the image patch using a weighted average method, and adjusts the assigned weights for image patches with abnormal variance using a weight adjustment strategy; after performing weight adjustment operations on all image patches, the weight coefficient matrix w is obtained.

[0074] Specifically, the weights are set using a weighted average method based on the variance of interference noise on the image patch, and the assigned weights are adjusted for image patches with abnormal variance using a weight adjustment strategy, as follows:

[0075] Variance of image patch noise The calculation method is shown in equation (4), taking the image block in row p and column q of the i-th reference image as an example:

[0076]

[0077] P and Q represent the number of rows and columns of an image block into which a reference image is cut, respectively; The number of pixels on an image patch rearranged in one direction (i.e., vertical or horizontal). Indicates interference noise on image patches The components at each pixel As shown in equation (5), Then it means The mean.

[0078]

[0079] In equation (5), X i (p,q) represents the noise residual of the image patch. For the average X i (p,q) thus obtains the approximate PRNU factor of the image region (the approximate PRNU factor of the position is obtained by averaging the image patches at the same position of the noise residuals of a set of small images), and N is the number of reference images used to calculate the camera reference PRNU noise.

[0080] After obtaining the variance of the image patch interference noise, the weights of the image patches can be obtained using the weighted average method of equation (6):

[0081]

[0082] In the formula, w k(p,q) represents the weight of the image patch at position (p,q) in the k-th reference image. Let represent the variance of the interference noise of the image patch at position (p,q) in the k-th image, and let i represent the i-th reference image.

[0083] The following problem arises when obtaining image patch weights according to equation (6):

[0084] When the interference noise variance of some image patches is close to 0, the weights increase sharply as the variance decreases. This may cause the device reference PRNU in a certain area to be dominated by a small number of image patches, rendering the remaining image patches meaningless in the calculation of the reference PRNU. Conversely, the role of image patches with excessively high variance in the calculation of the reference PRNU will be ignored. To address this issue, when the noise variance is abnormal, the sensitivity of the weights to changes in noise variance should be reduced.

[0085] Therefore, this invention performs threshold judgment on the interference noise variance. When the interference noise variance is greater than λ1 or less than λ2, the variance is adjusted, thereby reducing the sensitivity of the weight to abnormal interference noise variance; σ in formula (6) ο 2 (p,q) is adjusted to be calculated as shown in equation (7). Where ο = i, k:

[0086]

[0087] In equation (7), λ1 = 3.6 and λ2 = 0.4, which represent the upper and lower limits of the noise variance, respectively; r is the sensitivity factor, r = 3, which determines the sensitivity of the weights to changes in variance; σ ο 2 (p,q) represents the variance of the originally calculated interference noise. This represents the variance of the interference noise after being constrained using a weighting adjustment strategy.

[0088] A4, such as Figure 2 As shown, the PRNU multiplicative factor K of the camera is obtained by using the weighted coefficient matrix w combined with the weighted maximum likelihood estimation method.

[0089] The noise variance weighted maximum likelihood estimation algorithm can be expressed as equation (8):

[0090]

[0091] X i I represents the noise residual of the i-th reference image (the noise residual of the entire image obtained after merging the noise residuals of image patches). i Let w represent the i-th reference image. iThis represents the weight coefficient matrix, which consists of the weights of all image blocks in the i-th reference image obtained through steps A1 to A3 above.

[0092] A5. The multiplicative factor K is zero-mean processed to further suppress unnecessary artifacts caused by CFA interpolation and other equipment processing operations. Then, the zero-mean processed multiplicative factor K is multiplied with all reference images captured by the camera and then averaged to obtain the reference PRNU noise of the camera.

[0093] The zero-mean process involves subtracting the column mean from each value in the column of the multiplicative factor K, and then subtracting the row mean from each value in the row.

[0094] S2, the test image is denoised using a dual-tree complex wavelet filter, and the PRNU noise of the test image is extracted by subtracting the denoised image from the test image.

[0095] S3. Use peak-to-correlation energy (PCE) to calculate the correlation between the PRNU noise of the image under test and the reference PRNU noise of all cameras in the PRNU noise database. When the correlation is the highest, it is determined that the image under test comes from that camera; otherwise, the image under test does not come from that camera.

[0096] In this step, the peak-related energy is expressed as:

[0097]

[0098] In the formula, MN is the product of dimensions of the non-uniform noise in the matched optical response, and NCC is... RQ (i,j) represents the two-dimensional cyclic cross-correlation between the camera's reference PRNU noise R and the test image's PRNU noise Q at point (i,j):

[0099]

[0100] Where m and n represent the two dimensions of R and Q, R(i,j) represents the value of R at point (i,j), and Q(k+i,l+j) represents the value of Q at point (k+i,l+j). The average value of R is... The average value of Q, while NCC RQ (0,0) is:

[0101]

[0102] The above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.

Claims

1. A small size image based source camera identification method, characterized in that, The method comprises the following steps: S1, obtaining reference PRNU noises of multiple cameras and storing the reference PRNU noises in a PRNU noise database; S2, performing denoising on a test image by using a dual-tree complex wavelet filter, and subtracting the denoised image from the test image to obtain PRNU noise of the test image; S3, calculating the correlation between the PRNU noise of the test image and the reference PRNU noises of all cameras in the PRNU noise database, and determining that the test image comes from a camera with the highest correlation. In step S1, the specific steps for obtaining the reference PRNU noise of a camera are as follows: A1, capturing a group of images by using the camera, and obtaining a group of small-size images by intercepting the center part of the group of images as reference images of the camera; dividing each reference image into multiple image blocks, and expanding the content of the reference image neighborhood to the boundaries of the image blocks; A2, estimating the noise variance of the wavelet subband by using a median estimator after the dual-tree complex wavelet transform of the image blocks processed in step A1, and then applying the noise variance to a windowed Wiener filter to estimate the wavelet coefficients, followed by inverse dual-tree complex wavelet transform, then restoring the image block boundary expansion and subtracting the denoised image block from the image block before denoising to obtain the noise residual of each image block; A3, setting the weight according to the variance of the interference noise on the image block by using a weighted average method, and adjusting the weight given to the image block with abnormal variance by using a weight adjustment strategy; obtaining a weight coefficient matrix after the weight adjustment operation on all image blocks; A4, obtaining the PRNU multiplicative factor K of the camera by using the weight coefficient matrix and a weighted maximum likelihood estimation method; A5, performing zero-mean processing on the multiplicative factor K, and then multiplying the multiplicative factor K after the zero-mean processing with all the reference images captured by the camera and then averaging to obtain the reference PRNU noise of the camera; The weight is set according to the variance of the interference noise on the image block by using a weighted average method, which comprises: variance of the image block interference noise The calculation is shown in equation (4) using the image block of the ith reference image p rows q columns as an example: P and Q represent the number of rows and columns of the reference image cut into image blocks, respectively; the number of pixels on the image block rearranged in one direction, representing the interference noise on the image block the components of each pixel, as shown in equation (5), then represents the mean value of In formula (5), X i (p, q) represents the noise residual of the image block, is X i the mean of (p, q), and N is the number of reference images; After obtaining the variance of the interference noise of the image block, the weight of the image block is obtained by using the weighted average method of formula (6): where w k (p, q) denotes the weight of the image block at position (p, q) of the kth reference image, denotes the variance of the interference noise of the image block at position (p, q) of the kth image, i denotes the ith reference image; The weight given to the image block with abnormal variance is adjusted by using a weight adjustment strategy, which comprises: When the interference noise variance is greater than λ1 or less than λ2, the variance is adjusted, specifically σ o 2 (p, q) is adjusted to be calculated as formula (7) where o = i, k: In formula (7), λ1 and λ2 represent upper and lower threshold values of noise variance, respectively; r is a sensitive factor, represents interference noise variance after restriction by the weight adjustment strategy.

2. The small-size image based source camera identification method of claim 1, wherein, For the specific process of estimating the noise variance of the wavelet subband by using a median estimator and then applying the noise variance to a windowed Wiener filter to estimate the wavelet coefficients in step A2, the filtering estimation process can be represented by formula (1) after two-dimensional dual-tree complex wavelet transform is performed on each image block to obtain wavelet coefficients: The PRNU multiplicative factor K of the camera is obtained by using the weight coefficient matrix and a weighted maximum likelihood estimation method, which comprises: where W in denotes the wavelet coefficients before filtering, W out (u, v) denotes the wavelet coefficients at pixel point (u, v) after filtering, and σ2 and are obtained from equations (2), (3). In formula (2), median() represents a median estimator, N h is a local window with (u, v) as the center point and h x h as the size, (i, j) represents a point in the local window, W temp represents wavelet coefficients of a high frequency subband after a dual-tree complex wavelet transform of the image, Max() represents taking the maximum value of 0 and the variance estimate, and the min() function represents taking the minimum value of all window estimates. 3.The small-size image-based source camera identification method of claim 1, wherein, The noise variance weighted maximum likelihood estimation algorithm is represented by formula (8): The process of zero-mean processing is as follows: X i denotes the noise residual of the i-th reference image, I i denotes the i-th reference image, w i denotes a weight coefficient matrix composed of weights of all image blocks of the i-th reference image, N is the number of reference images. 4.The small-size image-based source camera identification method of claim 1, wherein, Subtract the average value of the column from each value in the column of the multiplicative factor K, and then subtract the average value of the row from each value in the row. The correlation between the PRNU noise of the test image and the reference PRNU noises of all cameras in the PRNU noise database is calculated by using peak correlation energy, which comprises: 5.The small-size image-based source camera identification method of claim 1, wherein, The peak correlation energy PCE is represented by formula (9): ​ where MN is the product of the dimensions of the matched optical response non-uniform noise, NCC RQ (i,j) is the two-dimensional cross-correlation of the reference PRNU noise R of the camera with the PRNU noise Q of the test image at point (i,j): where m and n represent the two dimensions of R and Q, R(i,j) represents the value of R at point (i,j), and Q(k+i,l+j) represents the value of Q at point (k+i,l+j), is the average value of R, is the average value of Q.

Images