A method for global bias correction of an image for patch processing

Abstract

The application discloses a kind of global bias correction methods for block processing image, belong to computer image processing field.For the problem of global brightness inconsistency in existing block splicing technology, low efficiency of iterative solution, the application is by being blocked to input image and extracting adjacent block overlap area, the number of overlapping area pixels and average gray difference information are counted;Linear constraint matrix is constructed based on the bias difference relationship between adjacent blocks, and graph Laplace type matrix is formed from its normal equation;Regularization term is introduced to construct stable sparse linear equation set, solve the additive bias parameters of each block, and then carry out gray correction on each block.The method models the brightness relationship between blocks from a global perspective, which can effectively reduce the problem of global gray inconsistency, and is suitable for two-dimensional / three-dimensional image deep learning block reasoning splicing and correction scenarios where the multi-region acquisition image approximately exhibits additive brightness offset.

Description

Technical Field

[0001] This invention relates to the field of computer image processing, specifically to a block-based image global bias correction method, applicable to deep learning block inference stitching, and correction scenarios in multi-region acquired images where the bias is approximately additive. Background Technology

[0002] With the rapid improvement of graphics processing unit (GPU) computing power, deep learning-based image processing methods have been widely used in the field of computer vision. However, in real-world deployment environments with limited GPU memory, inference networks often struggle to process large images (such as pathological slide images, remote sensing images, and high-resolution microscopic images) in one go. To overcome the GPU memory bottleneck, the current mainstream approach is to divide large images into several sub-image blocks, input each sub-image into the network for inference, and then stitch the outputs of each sub-block together according to their spatial location to form a complete image.

[0003] During the block stitching process, due to potential systematic deviations in the acquisition, reconstruction, or network inference processes of each block, inconsistencies in grayscale often occur in the overlapping areas of adjacent blocks, resulting in seams or abrupt brightness changes in the stitched result. To address this issue, existing technologies mainly employ the following two approaches: The first type of approach is a stitching method based on local weighted fusion of overlapping regions. This method constructs a linear or Gaussian weight function based on the distance from the pixel to the boundary within the common overlapping area of ​​adjacent blocks, and then performs a weighted average of the pixel values ​​in the overlapping area to achieve local smoothing of grayscale transitions. For example, this weighted fusion strategy is used in OpenCV's stitching module and most deep learning block inference frameworks. However, this type of method essentially only performs local patching in the overlapping area and does not model and correct the overall brightness shift between blocks from a global perspective. When there are many blocks or systematic brightness drift exists between blocks, weighted fusion cannot eliminate global brightness inconsistencies, and the stitched result will still show obvious grayscale step distribution, seriously affecting the visual effect of the image and the accuracy of subsequent quantitative analysis.

[0004] The second category of approaches is global brightness correction methods based on iterative optimization. These methods model the brightness consistency between blocks as a global optimization problem, estimating the bias parameters of each block using iterative algorithms such as gradient descent. For example, some methods employ automatic differentiation techniques to calculate the gradient of the objective function with respect to the bias parameters of each block, and update the parameters through multiple iterations until convergence. While this type of method can achieve global brightness consistency to some extent, it requires multiple iterations, resulting in high computational complexity. Especially when the number of blocks is large (such as 3D volumetric data blocks or large-scale tiled images), the time consumed by iterative optimization increases significantly, making it difficult to meet the efficiency requirements of real-time processing or large-scale data production.

[0005] In addition, existing technologies have the following shortcomings: First, some methods are designed for specific two-dimensional image processes, and when extended to three-dimensional volume data or multi-dimensional image sequences, the adjacency relationship and overlapping area statistics need to be redesigned; Second, for the brightness inconsistency problem that is approximately manifested as additive brightness shift during multi-region image acquisition, existing methods lack a unified and efficient global correction mechanism; Third, although some methods based on histogram matching or color transfer can achieve global brightness adjustment, these methods use the entire image or entire layer as the processing unit, do not make full use of the adjacency relationship and overlapping area constraints between blocks, have limited correction accuracy, and are prone to introducing global tone distortion.

[0006] In summary, existing technologies for global brightness consistency correction in segmented image stitching generally suffer from shortcomings such as insufficient local repair capabilities, lack of global modeling, low iterative solution efficiency, difficulty in multidimensional expansion, and insufficient scene adaptability. There is an urgent need for an image correction method that can model the brightness relationship between segments from a global perspective and can be efficiently solved through sparse linear equations. Summary of the Invention

[0007] The purpose of this invention is to overcome the shortcomings of existing block image stitching techniques, such as poor global brightness consistency, low computational efficiency, and limited applicability, and to provide a global bias correction method for block processing. Specifically, this invention provides a global bias modeling method based on gray-level consistency constraints in overlapping regions. By constructing a constraint matrix of bias differences between adjacent blocks and forming a graphical Laplacian matrix from its normal equations, the brightness bias estimation problem of each block is transformed into a problem of solving a system of regularized sparse linear equations. This allows for unified modeling and correction of the overall brightness shift between blocks from a global perspective, effectively reducing global gray-level inconsistencies and brightness drift in the stitched image. This invention also provides a general correction method applicable to block processing of two-dimensional images and three-dimensional volumetric data, which can be embedded as an independent post-processing module in deep learning inference processes or multi-region image acquisition processing processes.

[0008] This invention provides a global image bias correction method for block-based image processing. Based on gray-level consistency constraints in overlapping regions, this method constructs a bias difference constraint matrix between adjacent block pairs. The normal equation of this constraint matrix forms a graphical Laplacian matrix, transforming the bias estimation problem into a regularized sparse linear equation system for solution. This achieves efficient and stable correction of the overall brightness shift of each block. The following section elaborates on the technical solution of this invention, including the method flow, algorithm principle, core parameters, and engineering implementation details.

[0009] The core technical feature of this invention lies in transforming the traditional iterative optimization problem of global brightness consistency into a problem of solving a sparse linear equation system based on a graph structure. By solving the regularized linear system, the bias parameters of each block are obtained, avoiding multiple rounds of gradient descent updates to the bias parameters. The overall processing flow of this method includes the following seven steps: Step S1: Divide the input image into blocks to obtain multiple sub-image block data; Step S2: Determine the adjacency relationship between each block and extract the overlapping area between adjacent blocks; Step S3: Count the number of pixels and grayscale differences in each overlapping region; Step S4: Construct a linear constraint matrix based on the bias difference relationship between adjacent blocks, and form a graph Laplace matrix from the normal equation of the linear constraint matrix. At the same time, construct a constraint vector based on the average gray level difference of the overlapping area. Step S5: Introduce regularization terms into the graph Laplace-type normal equations to construct a stable sparse linear system of equations; Step S6: Use the sparse linear equation solver to solve the regularized linear equation system and obtain the bias parameters of each block. Step S7: Perform grayscale correction on each block according to the obtained bias parameters.

[0010] In the above process, steps S1 to S3 are the data preprocessing stage, steps S4 to S6 are the global modeling and analytical solution stage, and step S7 is the result output stage. The technical details and algorithm principles of each step are described in detail below.

[0011] Further, in step S1, the input image can be a two-dimensional grayscale image, a two-dimensional multi-channel color image, or three-dimensional volume data. For a two-dimensional image, let the size of the input image be H×W (height×width). Based on the memory capacity and network inference requirements, the image is divided into N_patch sub-image blocks. The block division strategy can adopt uniform grid partitioning or adaptive size partitioning. In the uniform grid partitioning mode, the size of a single sub-block is set to h×w, and the overlap widths of adjacent sub-blocks in the horizontal and vertical directions are o_x and o_y, respectively. Then, the number of blocks in the horizontal direction N_x and the number of blocks in the vertical direction N_y satisfy: N_x = ⌈(W - o_x) / (w - o_x)⌉; N_y = ⌈(H - o_y) / (h - o_y)⌉; Total number of patches N_patch = N_x × N_y.

[0012] To ensure the continuity of subsequent global modeling, the overlap widths o_x and o_y can be set according to network inference boundary effects, memory limitations, and stitching quality requirements. For 3D volumetric data, the block partitioning strategy can be extended to 3D uniform mesh partitioning, with the overlap thicknesses of adjacent sub-blocks in the x, y, and z directions being o_x, o_y, and o_z, respectively, and the total number of blocks N_patch = N_x × N_y × N_z.

[0013] Further, in step S2, the adjacency relationships between blocks are determined based on their spatial coordinates. For two-dimensional images, adjacency relationships can be constructed using adjacent blocks in the horizontal and vertical directions; for three-dimensional volume data, a six-neighborhood or twenty-six-neighborhood adjacency criterion can be used. Let 𝒩(i) be the neighborhood set of the i-th block, which contains the index j of all blocks that overlap with the i-th block. Adjacency relationships can be stored using a sparse adjacency list or an adjacent block list to reduce memory usage and improve the efficiency of subsequent matrix construction.

[0014] After determining the adjacency relationship, the pixel set of the overlapping region between adjacent blocks is extracted. Let V_i^j be the pixel set of the overlapping region between the i-th block and the j-th block in the local coordinate system of the i-th block, and V_j^i be the pixel set in the local coordinate system of the j-th block. The number of pixels in the overlapping region is denoted as N_ij, and its value is uniquely determined by the block spatial coordinates and the overlap width. For two-dimensional blocks, N_ij = o_x × h or N_ij = w × o_y or N_ij = o_x × o_y (diagonal overlap), depending on the adjacency direction; for three-dimensional blocks, N_ij is the number of voxels in the overlapping volume in the corresponding direction.

[0015] Furthermore, in step S3, to construct global consistency constraints, it is necessary to accurately statistically analyze the grayscale difference information of each overlapping region. This invention employs a method based on pixel-level grayscale difference mean statistics to improve robustness to noise and reduce the impact of outliers.

[0016] Specifically, for each pair of adjacent blocks (i, j), the grayscale difference between corresponding pixels within their overlapping region is calculated. Let the grayscale value of the nth pixel in V_i^j be v_i^j(n), and the grayscale value of the corresponding pixel in V_j^i be v_j^i(n). Then, the grayscale difference statistic R_ij of the overlapping region is defined as: R_ij = (1 / N_ij) × Σ_{n=1}^{N_ij} [v_j^i(n) - v_i^j(n)] The physical meaning of this statistic R_ij is: the average grayscale difference between the j-th block and the i-th block within the overlapping region. If R_ij > 0, it indicates that the overall brightness of the j-th block in the overlapping region is higher than that of the i-th block; if R_ij < 0, it indicates that the brightness is lower than that of the i-th block. This mean-based statistical method can compress the pixel or voxel differences in the overlapping region into a single adjacent block bias difference constraint.

[0017] In engineering implementation, to balance computational efficiency and memory usage, grayscale statistics of overlapping areas can be completed during the block inference stage or the post-processing stage. The statistical results can be stored in the form of an adjacent block pair list and a constraint vector: the adjacent block pair list records the two block indices corresponding to each constraint, and the constraint vector records the average grayscale difference R_ij of the corresponding overlapping areas.

[0018] Furthermore, in step S4, given that the core innovation of this invention lies in modeling the brightness inconsistency problem between blocks as a global optimization problem and designing an analytically solvable objective function, let the grayscale bias coefficient of the i-th block be b_i, which represents the offset amount for shifting and correcting the overall brightness of the block. After correction, the pixel grayscale value of the i-th block in the overlapping region V_i^j changes from v_i^j(n) to v_i^j(n) + b_i, and the pixel grayscale value of the j-th block in the overlapping region V_j^i changes from v_j^i(n) to v_j^i(n) + b_j.

[0019] To achieve global brightness consistency, this invention constructs the following global bias optimization objective function: L_global^b = Σ_{(i,j)∈ℰ} (1 / N_ij) × || (v_i^j + b_i) - (v_j^i +b_j) ||² (1) In Equation (1), L_global^b is the global bias optimization objective function; ℰ is the set of all adjacent block pairs; (i, j) is a pair of adjacent blocks; N_ij is the number of pixels in the overlapping area of ​​block i and block j; v_i^j is the pixel gray value of the i-th block corresponding to block j in the overlapping area; b_i is the bias parameter of the i-th block; ||·||² represents the squared L2 norm. The objective function (1) achieves global brightness consistency by minimizing the weighted sum of the corrected gray value differences of all adjacent block pairs in the overlapping area.

[0020] Furthermore, in step S4, the bias difference constraints of all adjacent block pairs are uniformly represented as a linear least squares problem. Let N_patch be the total number of blocks, K be the number of adjacent block pairs, then C is the linear constraint matrix of K × N_patch, b is the bias vector of N_patch × 1, and r is the constraint vector of K × 1.

[0021] This equation indicates that for each constraint b_i - b_j = R_ij, the expected bias difference between the i-th and j-th blocks is equal to the average gray-level difference R_ij of the overlapping regions. Since this system of equations is overdetermined and cannot precisely satisfy all constraints, the least squares method is used to solve for the optimal bias vector b, i.e., minimizing the sum of squares of the constraint residuals: The corresponding optimization problem is: b* = argmin_b L_global^b (2) In equation (2), b = [b_1, b_2, ..., b_{N_patch}]^T is the bias vector to be determined, and b* is the optimal bias parameter. This optimization problem is essentially an unconstrained quadratic optimization problem with respect to the bias variable b.

[0022] Construct the graphical Laplace-type normal equations from the constraint matrix C and the constraint vector r: First, by taking the partial derivative with respect to each bias variable b_k (k = 1, 2, ..., N_patch) and setting it to zero, we obtain the following system of derivative equations: { ∇{b_1} L_global^b = 0 ∇{b_2} L_global^b = 0 ⋮ ∇{b_{N_patch}} L_global^b = 0} (3) In equation (3), ∇{b_k} represents the partial derivative operator with respect to the bias variable b_k, and L_global^b is the global bias optimization objective function defined in equation (1). This system of equations represents obtaining N_patch simultaneous equations by taking the partial derivatives with respect to all bias variables and setting them to zero.

[0023] Taking the k-th block as an example, its partial derivative equation is: ∇_{b_k} L_global^b = 0 (4) In equation (4), ∇_{b_k} represents the partial derivative operator with respect to the bias variable b_k, and L_global^b is the global bias optimization objective function defined in equation (1).

[0024] Considering the symmetry of the constraint terms between adjacent blocks in the objective function (1) (i.e., the block pairs (i, j) and (j, i) describe the same overlapping region), the gradient needs to be expanded during the differentiation process. Specifically, only the adjacency constraint terms directly related to b_k are retained: ∇{b_k} [ Σ{j∈𝒩(k)} (1 / N_kj) × || V_k^j + b_k - (V_j^k + b_j) ||_2^2 ] = 0 (5) In equation (5), V_k^j is the pixel grayscale vector of the k-th block in the overlapping region V_k^j, V_j^k is the pixel grayscale vector of the j-th block in the overlapping region V_j^k, b_k is the bias parameter of the k-th block, b_j is the bias parameter of the j-th block, and ||·||_2^2 represents the squared L2 norm. This equation represents taking the partial derivative of the adjacency constraint term of the objective function (1) with respect to b_k and setting it to zero.

[0025] Take the chain rule derivative of the L2 norm squared term in equation (5). Let the number of pixels in the overlapping region be N_kj, then the dimensions of vectors V_k^j and V_j^k are both N_kj × 1. The derivative of the L2 norm squared term with respect to b_k is: ∇_{b_k} || V_k^j + b_k - (V_j^k + b_j) ||2^2 = 2 × [V_k^j + b_k -(V_j^k + b_j)]^T · 1{N_kj} Where 1_{N_kj} is an all-1 vector of N_kj × 1. Substituting the above result into equation (5) and considering the normalization coefficient 1 / N_kj, we get: Σ_{j∈𝒩(k)} (2 / N_kj) × [V_k^j + b_k - (V_j^k + b_j)]^T · 1_{N_kj} = 0 (6) Note that [V_k^j]^T · 1_{N_kj} = Σ_{n=1}^{N_kj} v_k^j(n), which is the sum of the gray values ​​of the pixels in the overlapping region. Transform the vector inner product in equation (6) into a scalar summation form: Σ_{j∈𝒩(k)} (2 / N_kj) × Σ_{n=1}^{N_kj} [v_k^j(n) + b_k - v_j^k(n) - b_j] = 0 Canceling the constant factor 2 and separating the known and unknown quantities: Σ_{j∈𝒩(k)} [ (Σ_{n=1}^{N_kj} [v_k^j(n) - v_j^k(n)]) / N_kj + b_k- b_j ] = 0 (7) In equation (7), the first term (Σ[v_k^j - v_j^k]) / N_kj is the grayscale difference R_kj statistically obtained in step S3. Further rearranging equation (7), separating variables and rearranging terms: Σ_{j∈𝒩(k)} b_k - Σ_{j∈𝒩(k)} b_j = Σ_{j∈𝒩(k)} R_kj Let n_k^{nbr} = |𝒩(k)| be the number of neighbors of the k-th block (i.e., the number of blocks that overlap with this block), then the above formula can be written as: n_k^{nbr} × b_k - Σ_{j∈𝒩(k)} b_j = - Σ_{j∈𝒩(k)} (Σ_{n=1}^{N_kj} [v_k^j(n) - v_j^k(n)]) / N_kj (8) The physical meaning of equation (8) is that the bias parameter b_k of the k-th block is linearly constrained by the bias parameters b_j of all its adjacent blocks and the gray-scale difference observation of the overlapping region. This equation transforms the solution of each bias variable into a linear constraint relationship with respect to its neighborhood.

[0026] Furthermore, in step S4, the bias difference constraints of all adjacent block pairs can be uniformly represented as a linear least squares problem: C × b ≈ r; Further construct the normal equation: C^TC × b = C^T r (9) In equation (9), C is the linear constraint matrix of adjacent blocks of K × N_patch, and K is the number of adjacent blocks; C^TC is the graph Laplace matrix of N_patch × N_patch, which describes the adjacency relationship between blocks; b is the bias vector of N_patch × 1; and r is the constraint vector of K × 1, which represents the observation obtained by statistical analysis of the average gray difference of the overlapping area.

[0027] The elements of matrix A are defined as follows: A[i][j] = { n_i^{nbr}, i = j { -1, j ∈ 𝒩(i) { 0, j ∉ 𝒩(i) (10) In equation (10), the diagonal element A[i][i] equals the number of adjacency constraints involved in the i-th block, and the off-diagonal element A[i][j] takes the value of -1 or the negative value of the corresponding weight when there is an adjacency constraint between block i and block j, otherwise it is 0. This matrix can be obtained by C^TC, which is equivalent to the graph Laplace matrix of the block adjacency graph. For a two-dimensional four-neighborhood block layout, each row of the matrix contains at most 5 non-zero elements; for a three-dimensional six-neighborhood block layout, each row of the matrix contains at most 7 non-zero elements.

[0028] The elements of the constraint vector r are defined as follows: C^T r is the right-hand side term of the normal equation, which can be obtained by matrix multiplication of the constraint vector r between adjacent blocks; if written in node form, its k-th element is equivalent to the signed sum of the average gray-level differences of the overlapping regions adjacent to the k-th block.

[0029] r[k] = Σ_{j∈𝒩(k)} R_kj (11) In equation (11), r[k] is the sum of the average grayscale differences between the k-th block and all its neighboring blocks within the overlapping region. The construction of this constraint vector only requires traversing the adjacency table once, with a computational complexity of O(N_patch × n_{avg}), where n_{avg} is the average number of neighborhoods of the block. For regular grid blocks, n_{avg} is a constant.

[0030] In the engineering implementation, a K × N_patch sparse constraint matrix C is first constructed using the COO format. For each adjacent block pair (i, j), 1 is written in the i-th column of the corresponding row, and -1 is written in the j-th column. Then, C is converted into a sparse format suitable for the solver, and C^TC and C^Tr are calculated.

[0031] Furthermore, in step S5, since the graph Laplace matrix has a global translational degree of freedom, and may exhibit ill-conditioned or near-singular problems when the block adjacency relationship is weak or the number of constraints is insufficient, a regularization term is introduced for improvement: (C^TC + αI) × b = C^T r (12) In equation (12), I is the identity matrix of N_patch × N_patch, and α is a preset regularization coefficient. In one embodiment, α can be 1e^{-6}; in other embodiments, α can be adjusted according to the dynamic range of the image, the number of blocks, and the matrix condition number.

[0032] The physical meaning of the regularization term αI is: to introduce an L2 regularization constraint on the bias parameter b in the bias estimation, to suppress the excessively large overall correction amount, and to improve the condition number of the graph Laplace matrix, thereby improving the stability of the numerical solution.

[0033] Further, in step S6, the regularized linear equation system (12) is solved using a linear algebra solver to obtain the bias vector b. Since the matrix (C^TC + αI) is a symmetric positive definite sparse matrix, the present invention preferably employs the following two solution strategies: In a possible implementation, a sparse linear equation solver is used to solve the regularized system of linear equations. Specifically, C^TC and C^Tr are first constructed, and then (C^TC + αI)b = C^Tr is solved to obtain the bias vector b.

[0034] When the number of blocks is large, iterative solution methods or preprocessing solution methods suitable for sparse positive definite linear systems can be used as optional implementation methods; the specific solver can be selected according to the engineering language, matrix size and memory constraints.

[0035] In this way, the bias estimation is transformed into solving a sparse linear system, avoiding hundreds or thousands of gradient descent parameter updates to the original objective function.

[0036] After solving, the bias parameters b_k (k = 1, 2, ..., N_patch) of each block are obtained. In a possible implementation, the bias vector can be centered, i.e., b ← b - mean(b), so that the average value of the biases of all blocks is 0, so as to avoid introducing an overall brightness shift.

[0037] Further, in step S7, grayscale correction is performed on each block based on the obtained bias parameters. The correction method adopts an additive model, that is, the bias parameter b_k is uniformly added to the grayscale values ​​of all pixels in the k-th block: I_k^{corrected}(x, y) = I_k^{original}(x, y) + b_k Where I_k^{original}(x, y) is the pixel grayscale value of the k-th block at coordinates (x, y) before correction, and I_k^{corrected}(x, y) is the pixel grayscale value after correction. For multi-channel color images, this bias correction can be applied independently to each color channel (e.g., applying the same or independent bias parameters to the RGB channels respectively), or it can be applied only to the luminance channel (e.g., the Y channel in the YCbCr color space) to maintain the chromaticity information.

[0038] During the calibration process, the corrected grayscale values ​​need to be dynamically truncated to prevent overflow. I_k^{corrected}(x, y) = I_k^{original}(x, y) + b_k; When a fixed bit depth image needs to be output, dynamic range truncation or quantization can be further performed in the output stage.

[0039] Where I_k^{original}(x, y) is the pixel gray value of the k-th block at coordinates (x, y) before correction, b_k is the bias parameter of the k-th block, and I_k^{corrected}(x, y) is the pixel gray value after correction.

[0040] After the independent correction of each block is completed, the corrected blocks can be stitched together into a complete image according to their original spatial positions. Since global offset correction has reduced the additive brightness shift between blocks from an overall perspective, the grayscale transition in the overlapping areas of the stitched image is more natural. To further improve the visual effect, post-processing methods such as linear weighted fusion can be used in conjunction with the original stitching process.

[0041] The technical solution of this invention has good scalability and can be adapted to various imaging scenarios and data dimensions: (1) Two-dimensional image expansion: For deep learning block inference scenarios, after the neural network completes the inference of each sub-block, the inference result blocks are directly extracted as input, and the global bias correction process from step S2 to step S7 is executed. This correction module can be embedded as an independent post-processing unit into the existing inference framework without modifying the network structure and training process.

[0042] (2) 3D volumetric data expansion: For 3D medical images or microscopic volumetric data, the block strategy is expanded to a 3D uniform grid. The adjacency relationship can adopt the six-neighbor or twenty-six-neighbor criteria, and the overlapping area is a 3D voxel block. The objective function and the regularized linear equation system remain unchanged. Only the pixel index needs to be expanded to the voxel index, and the number of pixels in the overlapping area needs to be expanded to the number of voxels.

[0043] (3) Multi-region image acquisition correction extension: For images acquired from different fields of view or different regions, each acquisition region can be regarded as a block. When there is spatial overlap between regions and the brightness difference is approximately additive, the same global bias correction model is constructed after statistically analyzing the gray difference of the overlapping regions.

[0044] (4) Application of pathological slides and remote sensing image mosaicking: For multi-region acquisition scenarios such as pathological slide scanning and remote sensing image mosaicking, each acquired tile is regarded as a block, and the overlapping area between tiles is pre-determined by the scanner or mosaicking algorithm. This method can be directly applied to such scenarios to achieve global brightness equalization of large-scale tile images.

[0045] In summary, the technical solution of this invention performs global modeling through adjacent block constraint matrices and graph Laplacian normal equations, and achieves additive brightness bias correction of block images by solving a regularized sparse linear system. It has technical features such as global consistency, computational efficiency, numerical stability and multidimensional data adaptability.

[0046] The beneficial effects achieved by this invention are as follows: 1. Achieve global brightness consistency correction, effectively reducing splicing seams and grayscale drift caused by additive bias. This invention constructs an adjacent block pair constraint matrix and forms a graph Laplacian matrix from its normal equation, modeling the adjacency relationship between block images as a global graph structure. It utilizes the grayscale difference information of all overlapping regions to construct a unified set of linear constraint equations, thereby jointly estimating and correcting the overall brightness shift of each block from a global perspective.

[0047] 2. Solving linear equations with high computational efficiency and strong numerical stability. This invention transforms the traditional bias estimation problem, which requires multiple iterations of optimization, into a problem of solving a regularized sparse linear equation system through analytical modeling. Compared with iterative optimization methods based on gradient descent or automatic differentiation backpropagation, this invention avoids multiple rounds of gradient descent updates to the bias parameters and can improve solution efficiency by utilizing a sparse matrix structure.

[0048] 3. High versatility and wide applicability, supporting multi-dimensional data and various imaging scenarios. The technical framework of this invention is based on the adjacency relationship and overlapping area statistics between blocks, and does not depend on a specific network structure. It can be applied to two-dimensional image block stitching and three-dimensional volume data block processing. For the problem of approximately additive brightness shift in multi-region acquired images, the same global bias correction framework can also be used. Attached Figure Description

[0049] Figure 1 This is an overall flowchart of the image global offset correction method for block processing described in this invention; Detailed Implementation

[0050] The term "embodiment" used herein, as an example, is not necessarily to be construed as superior to or better than other embodiments. Performance tests in these embodiments of the invention, unless otherwise specified, employ conventional testing methods in the art. It should be understood that the terminology used herein is merely for describing particular implementations and is not intended to limit the scope of the disclosure.

[0051] Unless otherwise stated, the technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains; other experimental methods and techniques not specifically mentioned herein refer to experimental methods and techniques commonly used by one of ordinary skill in the art.

[0052] To better illustrate the content of this invention, numerous specific details are provided in the following detailed embodiments. Those skilled in the art should understand that the invention can be practiced even without certain specific details. In the embodiments, some methods, means, instruments, and devices well-known to those skilled in the art are not described in detail, in order to highlight the main points of the invention.

[0053] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All applications utilizing the concept of the present invention are protected.

[0054] This embodiment uses deep learning inference of digital pathological slide images as an application scenario to explain in detail the specific implementation process of the present invention in the processing of large-size two-dimensional images.

[0055] The input pathological slide image is a 16-bit grayscale image of 8192×8192 pixels, with grayscale values ​​ranging from 0 to 65535. Due to GPU memory limitations, the deep learning segmentation network cannot process the entire image at once and must divide it into several sub-blocks for separate inference. A uniform grid partitioning strategy is adopted, with the size of a single sub-block set to 2048×2048 pixels, and the overlap width of adjacent sub-blocks in both the horizontal and vertical directions being 256 pixels (i.e., an overlap rate of 12.5%).

[0056] Based on the formula for calculating the number of blocks: N_x = ⌈(8192 - 256) / (2048 - 256)⌉ = ⌈7936 / 1792⌉ = 5 N_y = ⌈(8192 - 256) / (2048 - 256)⌉ = 5 Total number of patches N_patch = N_x × N_y = 25.

[0057] The spatial coordinates of each block are determined by its row and column indices (p, q), where p = 0, 1, 2, 3, 4 are row indices, and q = 0, 1, 2, 3, 4 are column indices. The coordinates of the top-left corner of the k-th block (k = p × N_x + q) in the entire image are: x_k = q × (2048 - 256) = q × 1792 y_k = p × (2048 - 256) = p × 1792 The block size is uniformly 2048×2048 pixels, and the parts that exceed the image boundary are processed with zero padding or mirror padding.

[0058] A four-neighbor adjacency criterion is adopted. For internal non-boundary blocks (row index p ∈ {1,2,3}, column index q ∈ {1,2,3}), their neighborhood includes blocks in the four directions of top, bottom, left, and right; for boundary blocks, the number of neighborhoods is reduced accordingly. For example, the 0th block located in the top left corner (p=0, q=0) only has a right neighbor (block 1) and a bottom neighbor (block 5), with a neighborhood count of n_0^{nbr} = 2; the 12th block located in the center (p=2, q=2) has four neighbors: top (block 7), bottom (block 17), left (block 11), and right (block 13), with a neighborhood count of n_12^{nbr} = 4.

[0059] For horizontally adjacent block pairs (i, j), the overlapping region is a rectangular strip with a width of 256 pixels and a height of 2048 pixels, with the number of pixels N_ij = 256 × 2048 = 524288. For vertically adjacent block pairs, the overlapping region is also a rectangular strip with a width of 2048 pixels and a height of 256 pixels, with the same number of pixels (524288).

[0060] The extraction of overlapping pixel sets is achieved through array slicing. Let the i-th block be the current block, and its right neighboring block be j. Then the overlapping region V_i^j of the i-th block is the rightmost 256 columns of pixels in block i, i.e., the array slice [:, -256:]; the overlapping region V_j^i of the j-th block is the leftmost 256 columns of pixels in block j, i.e., the array slice [:, :256]. The extraction of overlapping regions in the vertical direction is similar, using row-direction slices [-256:, :] and [:256, :].

[0061] After the deep learning network completes the inference for each block, the inference output block (rather than the original input block) is extracted for grayscale statistics. Assume that the network output is a probability map or feature map with the same size as the input, and the grayscale value is also 16 bits deep.

[0062] For each pair of adjacent blocks (i, j), calculate the average grayscale difference of corresponding pixels within the overlapping region. Taking the horizontal adjacent block pair (0, 1) as an example: R_01 = (1 / 524288) × Σ_{m=0}^{2047} Σ_{n=0}^{255} [v_1^0(m, n) -v_0^1(m, n)] Where v_0^{1}(m, n) is the grayscale value of the pixel in the m-th row and n-th column of the overlapping region of block 0, and v_1^{0}(m, n) is the grayscale value of the corresponding pixel in the overlapping region of block 1. This summation operation is implemented through vectorized matrix operations, which can be efficiently performed using scientific computing libraries such as NumPy.

[0063] The statistical results are stored in the form of a list of adjacent blocks and constraint vectors. In this embodiment, there are 40 undirected adjacent blocks (5×4=20 pairs in the horizontal direction and 5×4=20 pairs in the vertical direction). Therefore, C is a 40×25 sparse constraint matrix and r is a 40×1 constraint vector.

[0064] Construct a 40×25 sparse linear constraint matrix C. For each pair of adjacent blocks (i, j), set C[k][i] = 1, C[k][j] = -1, and the remaining elements to 0 in the corresponding constraint row. Then construct a 25×25 graph Laplace matrix C^TC.

[0065] Initialize the diagonal element A[k][k] = 0 (k = 0, 1, ..., 24).

[0066] Iterate through all adjacent block pairs (i, j): Execute A[i][i] += 1, A[j][j] += 1 (accumulate the number of neighbors); Execute A[i][j] = -1, A[j][i] = -1 (mark the adjacency relationship).

[0067] Taking the 0th block as an example, it is adjacent to block 1 (right neighbor) and block 5 (lower neighbor), so A[0][0] = 2, A[0][1] = -1, A[0][5] = -1, and the rest of the off-diagonal elements are 0.

[0068] Construct a constraint vector r, where the k-th element is the average gray level difference R_ij between the corresponding adjacent blocks.

[0069] r[k] = Σ_{j∈𝒩(k)} R_kj Taking the 0th block as an example: Taking the 0th block as an example, the constraint formed with its right neighbor block 1 can be written as b_0 - b_1 = R_01, and the constraint formed with its lower neighbor block 5 can be written as b_0 - b_5 = R_05.

[0070] Where R_01 is the average grayscale difference between the overlapping areas of block 0 and block 1, and R_05 is the average difference between block 0 and block 5.

[0071] Matrix C is stored using a sparse format such as COO or CSC. In this embodiment, the number of non-zero elements in C is 2 × 40 = 80; the 25 × 25 graph Laplace matrix obtained from C^TC can continue to be stored using a sparse format.

[0072] Introduce a regularization coefficient α = 1e^{-6}. Construct the regularization matrix: A_reg = C^TC + αI Where I is a 25×25 identity matrix. After regularization, A_reg = C^TC + αI, and together with the right-hand side C^Tr, they form a regularized linear system of equations.

[0073] The sparse linear equation solver is used to solve the regularized normal equations.

[0074] (1) Perform Cholesky decomposition on A_reg: A_reg = L × L^T, where L is a lower triangular sparse matrix; (2) Solve the forward equation L × y = r to obtain the intermediate vector y; (3) Solve the backward equation L^T × b = y to obtain the bias vector b.

[0075] The solution process can call scipy.sparse.linalg.spsolve or other equivalent sparse linear equation solvers. In this embodiment, the equation to be solved is (C^TC + αI)b = C^T r.

[0076] The obtained bias vector b contains 25 elements, each corresponding to a grayscale bias parameter for one of the 25 blocks. After solving, the bias vector can be centered, i.e., b ← b - mean(b), to make the bias mean 0.

[0077] Grayscale correction is performed on each block based on the obtained bias parameters. Taking block 12 as an example, if its bias parameter b_12 = -1200, then additive correction is performed on all pixels in this block, that is, I_{12}^{corrected} = I_{12}^{original} + b_12, which is equivalent to subtracting 1200 from the whole.

[0078] I_{12}^{corrected}(x, y) = I_{12}^{original}(x, y) + b_12 After all 25 blocks have been corrected, they can be stitched together into a complete image according to their original spatial coordinates. Because global offset correction reduces the additive brightness shift between blocks, the grayscale transition in the overlapping areas of the stitched image is more natural; if necessary, weighted fusion can be performed in conjunction with the original stitching process.

[0079] This embodiment uses a stack of three-dimensional fluorescence images acquired by a confocal microscope as an application scenario to illustrate the extended implementation process of the present invention in the block processing of three-dimensional volume data.

[0080] The input 3D volume data is 1024×1024×512 voxels (x×y×z) with 16-bit grayscale depth. Due to the high memory consumption of 3D convolutional networks, the volume data needs to be divided into 3D sub-blocks for separate inference. A uniform grid is used, with sub-block sizes of 256×256×128 voxels. The overlap thickness of adjacent sub-blocks in the x, y, and z directions is 32, 32, and 16 voxels, respectively.

[0081] Calculation of the number of blocks: N_x = ⌈(1024 - 32) / (256 - 32)⌉ = ⌈992 / 224⌉ = 5 N_y = ⌈(1024 - 32) / (256 - 32)⌉ = 5 N_z = ⌈(512 - 16) / (128 - 16)⌉ = ⌈496 / 112⌉ = 5 The total number of patches N_patch = 5 × 5 × 5 = 125.

[0082] The six-neighbor adjacency criterion is adopted (face contact is considered adjacency). For internal voxel blocks, the number of neighbors is 6; for face boundary blocks, the number of neighbors is 5; for edge boundary blocks, the number of neighbors is 4; and for corner boundary blocks, the number of neighbors is 3.

[0083] The overlapping region in the x-direction has a size of 32×256×128 voxels, with a voxel count N_ij = 32×256×128 = 1048576; the same applies to the y-direction. The overlapping region in the z-direction has a size of 256×256×16 voxels, with a voxel count N_ij = 256×256×16 = 1048576.

[0084] Overlapping region extraction uses three-dimensional array slicing. For example, the overlapping region of block i and its right neighbor block j in the x direction: V_i^j is the rightmost 32 voxels of block i, and the slicing operation is [:, :, -32:]; V_j^i is the leftmost 32 voxels of block j, and the slicing operation is [:, :, :32].

[0085] For each pair of 3D adjacent blocks, calculate the average gray-level difference of the corresponding voxels within the overlapping region. Taking block pair (i, j) as an example: R_ij = (1 / N_ij) × Σ [v_j^i(x, y, z) - v_i^j(x, y, z)] The statistical process is implemented using three-dimensional vectorized operations, and the adjacency statistics of 125 blocks can be completed in a few seconds.

[0086] For a 5×5×5 3D block mesh, the total number of blocks is 125. If undirected adjacent block pairs with six neighbors are used, there are 300 adjacency constraints, and C is 300×125; if bidirectional six-neighbor traversal is used in the program, there are 600 directed constraints, and C is 600×125. Subsequently, a 125×125 graph Laplace matrix is ​​constructed using C^TC.

[0087] The constraint vector r is constructed in the same way as in the two-dimensional scene, with each term being the average gray level difference R_ij of the corresponding adjacent three-dimensional block overlapping voxel regions.

[0088] By introducing a regularization coefficient α = 1e^{-6}, we construct A_reg = C^TC + αI, and use a sparse linear equation solver to solve (C^TC + αI)b = C^T r, thus obtaining the bias parameters for 125 blocks.

[0089] Voxel-level grayscale correction is performed on each 3D sub-block based on the bias parameters: V_k^{corrected}(x, y, z) = V_k^{original}(x, y, z) + b_k; When a fixed bit depth data needs to be output, truncation or quantization can be performed during the output stage.

[0090] The 125 corrected sub-blocks were stitched together into a complete dataset according to their three-dimensional spatial coordinates. The additive brightness shift at the boundaries between adjacent sub-blocks in the stitched 3D image was reduced, which is beneficial for subsequent 3D cell segmentation and quantitative analysis.

[0091] This embodiment uses high-throughput fluorescence microscopy multi-field scanning imaging as an application scenario to illustrate the specific implementation process of the present invention in addressing the problem of inconsistent brightness in multiple fields of view, which is approximately manifested as additive brightness shift, in images acquired from multiple regions.

[0092] Imagine using an inverted fluorescence microscope to scan a 6 mm × 4 mm tissue section with a 20× objective lens. The single field-of-view image size is 2048 × 2048 pixels, and the pixel resolution is 0.5 μm / pixel. To cover the entire tissue area, the microscope stage moves in a serpentine scanning trajectory, acquiring a total of M × N = 6 × 4 = 24 field-of-view images (Tiles). Each field-of-view image is a 16-bit single-channel fluorescence image (e.g., the DAPI nuclear staining channel).

[0093] Adjacent fields of view maintain a 10% overlap in both the horizontal (X-axis) and vertical (Y-axis) directions, i.e., an overlap width of 205 pixels (approximately 2048 × 10%). The total number of patches is N_patch = 24, and each patch is numbered k = 0, 1, ..., 23 according to the scanning order. Its corresponding stage row and column indices are (p, q), where p = 0, 1, 2, 3 are row indices, and q = 0, 1, 2, 3, 4, 5 are column indices.

[0094] During the acquisition process, due to the slight attenuation of LED light source intensity over time (thermal drift during long-term scanning), fine-tuning of camera exposure time between different batches, and differences in fluorescence signal intensity caused by uneven tissue section thickness, there is a systematic brightness shift between different fields of view. Let the original image of the k-th field of view be I_k^{original}, with a size of 2048×2048.

[0095] The four-neighbor adjacency criterion is adopted to directly determine the adjacency relationship based on the physical scanning coordinates of the stage, without the need for additional image registration. For internal fields of view (row index p ∈ {1,2}, column index q ∈ {1,2,3,4}), its neighborhood includes the fields of view in the four directions of top, bottom, left, and right; for boundary fields of view, the number of neighborhoods is reduced accordingly. For example, the 0th field of view (p=0, q=0) located at the scanning start point has only a right neighbor (field of view 1) and a lower neighbor (field of view 6), with a neighborhood count of n_0^{nbr} = 2; the 14th field of view (p=2, q=2) located at the center has four neighborhoods, with a neighborhood count of n_14^{nbr} = 4.

[0096] For horizontally adjacent view pairs (i, j), the overlapping region is a rectangular strip with a width of 205 pixels and a height of 2048 pixels, with a pixel count N_ij = 205 × 2048 = 419840. For vertically adjacent view pairs, the overlapping region is also a rectangular strip with a width of 2048 pixels and a height of 205 pixels, with the same pixel count of 419840.

[0097] The pixel set of the overlapping region is extracted through precise stage position mapping. Let the step distance of the microscope stage in the X direction be d_x = 2048 - 205 = 1843 pixels, and the step distance in the Y direction be d_y = 1843 pixels. The global coordinates of the upper left corner of the k-th field of view are (q × d_x, p × d_y). Based on the global coordinates, the common overlapping region between adjacent fields of view can be directly cropped, avoiding interference from image registration errors on grayscale statistics.

[0098] Specifically, for horizontal adjacency pairs (i, j), the overlapping region V_i^j of view i is its rightmost 205 columns of pixels, i.e., the array slice [:, -205:]; the overlapping region V_j^i of view j is its leftmost 205 columns of pixels, i.e., [:, :205]. For vertical adjacency pairs, row-direction slices [-205:, :] and [:205, :] are used.

[0099] Before calculating the grayscale difference in the overlapping area, preprocessing operations such as background subtraction, outlier filtering, or effective pixel selection can be performed on the overlapping area; alternatively, the average grayscale difference can be directly calculated for the corresponding pixels in the overlapping area.

[0100] Background subtraction, outlier filtering, and effective pixel selection are optional preprocessing steps that can be performed depending on the fluorescence image signal-to-noise ratio, background non-uniformity, and actual application requirements.

[0101] For each pair of adjacent views (i, j), calculate the mean grayscale difference of the effective pixels: R_ij = (1 / N_ij) × Σ [v_j^i(x, y) - v_i^j(x, y)] Where Ω represents the set of valid pixels that meet the signal strength screening criteria, and N'_ij represents the number of valid pixels. This statistic reflects the systematic brightness shift of field of view i relative to field of view j within the overlapping region.

[0102] Taking field of view 0 and field of view 1 (horizontally adjacent) as an example, assuming the average effective pixel grayscale difference R_01 = 42.3 (16-bit grayscale), it indicates that the overall fluorescence signal of field of view 0 is about 42 grayscale levels higher than that of field of view 1. This difference may be caused by a slight decrease in light source intensity or local tissue thickness differences during the acquisition time interval between the two fields of view.

[0103] The adjacency statistics for all 24 fields of view can be stored as a list of adjacent field of view pairs and a constraint vector. For a 6-column × 4-row scan grid, there are 4 × (6-1) = 20 horizontal adjacency pairs and (4-1) × 6 = 18 vertical adjacency pairs, for a total of 38 undirected adjacency constraints.

[0104] Construct a 38×24 sparse linear constraint matrix C. For each pair of adjacent views (i, j), set C[k][i] = 1, C[k][j] = -1, and the remaining elements to 0 in the corresponding constraint row; then construct a 24×24 graph Laplace matrix C^TC.

[0105] Initialize the diagonal element A[k][k] = 0 (k = 0, 1, ..., 23).

[0106] Iterate through all 38 adjacent view pairs (i, j): Execute A[i][i] += 1, A[j][j] += 1; Execute A[i][j] = -1, A[j][i] = -1.

[0107] Taking field of view 0 as an example, it is adjacent to field of view 1 (right) and field of view 6 (bottom), so A[0][0] = 2, A[0][1] = -1, A[0][6] = -1. Taking field of view 14 (center) as an example, it is adjacent to field of view 8 (top), 20 (bottom), 13 (left), and 15 (right), so A

[14]

[14] = 4, corresponding to -1 for the four non-diagonal elements.

[0108] Matrix C is stored in a sparse format such as COO, CSR or CSC, and C^TC is a 24×24 graph Laplace sparse matrix.

[0109] The k-th element of the constraint vector r is defined as the average gray-scale difference R_ij of the overlapping region corresponding to the k-th adjacent view constraint.

[0110] r[k] = Σ_{j∈𝒩(k)} R_kj Taking the 0th field of view as an example: r[0] = (R_01 + R_06) Where R_01 represents the difference from field of view 1, and R_06 represents the difference from field of view 6.

[0111] In fluorescence microscopy images, the fluorescence signal intensity may differ significantly between different fields of view (e.g., between tissue edges and the center), potentially resulting in a large condition number for the Laplacian matrix A. By introducing a regularization coefficient α = 1e^{-5}, a regularization matrix is ​​constructed: A_reg = C^TC + αI After regularization, the matrix A_reg = C^TC + αI can improve the numerical stability of the linear equation system.

[0112] The regularized system of equations (C^TC + αI)b = C^T r is solved using a sparse linear equation solver.

[0113] The solution yields 24 bias parameters, each corresponding to a global additive brightness offset for one of the 24 fields of view. After solving, the bias vector can be centered, i.e., b ← b - mean(b), so that the average value of the biases for all fields of view is 0.

[0114] Global grayscale correction is performed on the fluorescence images of each field of view based on the obtained bias parameters. Taking field of view 0 as an example, if its bias parameter b_0 = 85.2, then additive correction is performed on all pixels in this field of view: I_0^{corrected}(x, y) = I_0^{original}(x, y) + b_0; When a 16-bit unsigned integer is required as the output, truncation or quantization can be performed during the output stage.

[0115] Since the bias parameter is a floating-point number, the correction process uses floating-point operations to maintain precision, and the final result is quantized into a 16-bit unsigned integer.

[0116] Because global offset correction reduces the additive brightness shift between different fields of view, the grayscale transition at the boundaries of the stitched panoramic fluorescence image is more natural. Overlapping areas can be combined with simple linear weighted fusion or other existing stitching processes to further reduce the impact of mechanical positioning errors or local noise.

[0117] The corrected panoramic image can be directly used for subsequent cell nucleus segmentation, density statistics, and morphological analysis, avoiding the difficulty of adaptive segmentation threshold and quantitative analysis bias caused by inconsistent brightness.

[0118] In summary, the present invention overcomes the shortcomings of the prior art and achieves the aforementioned technical effects through the following technical means: First, in response to the shortcomings of existing technologies that "rely only on local weighted fusion and are difficult to correct global brightness drift", this invention constructs a global adjacency relationship model through adjacent block constraint matrices and graph Laplace normal equations to achieve global bias correction.

[0119] Second, in response to the shortcomings of existing technologies that "rely on iterative optimization for solution, have high computational complexity and low efficiency", this invention improves the solution efficiency by using analytical modeling and regularized sparse linear equation system solution to avoid multiple rounds of gradient descent updates to the bias parameters.

[0120] Third, in response to the shortcomings of existing technologies that are "difficult to apply to multi-dimensional images and different imaging conditions", this invention achieves unified correction of additive brightness bias in two-dimensional images, three-dimensional volume data and multi-region acquired images through unified graph structure modeling and scalable adjacency definition.

[0121] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

[0122] The technical solutions and technical details disclosed in the embodiments of this invention are merely illustrative of the inventive concept of this invention and do not constitute a limitation on the technical solutions of this invention. Any conventional changes, substitutions or combinations made to the technical details disclosed in the embodiments of this invention have the same inventive concept as this invention and are within the protection scope of the claims of this invention.

Claims

1. A method for global image offset correction in block processing, characterized in that, Includes the following steps: S1: The input image is divided into blocks to obtain multiple sub-image block data; S2: Determine the adjacency relationship between each block and extract the overlapping area between adjacent blocks; S3: Statistical analysis of pixel count and grayscale difference information in each overlapping region; S4: Construct a linear constraint matrix based on the bias difference relationship between adjacent blocks, and form a graph Laplace matrix from the normal equation of the linear constraint matrix. At the same time, construct a constraint vector based on the average gray difference of the overlapping area. S5: Introduce regularization terms to the graph Laplace-type normal equations to construct a stable sparse linear system of equations; S6: Use a sparse linear equation solver to solve the regularized system of linear equations and obtain the bias parameters for each block. S7: Perform grayscale correction on each block according to the obtained bias parameters.

2. The method according to claim 1, characterized in that, The statistical analysis of the number of pixels and grayscale differences in each overlapping region includes: For each pair of adjacent blocks, calculate the average gray-level difference of corresponding pixels or voxels within their overlapping regions to obtain the average gray-level difference constraint R_ij: R_ij = (1 / N_ij) × Σ_{n=1}^{N_ij} [v_j^i(n) - v_i^j(n)] Where N_ij is the number of pixels or voxels in the overlapping area between the i-th block and the j-th block, v_i^j(n) is the gray value of the n-th pixel or voxel in the overlapping area of ​​the i-th block, v_j^i(n) is the gray value of the corresponding pixel or voxel in the overlapping area of ​​the j-th block, and R_ij represents the average gray value difference between the j-th block and the i-th block in the overlapping area.

3. The method according to claim 1, characterized in that, The construction of a linear constraint matrix based on the bias difference relationship between adjacent blocks, and the formation of a graphical Laplace matrix from its normal equation, includes: Construct a sparse linear constraint matrix C of size K × N_patch, where K is the number of adjacent block pairs and N_patch is the total number of blocks. For the k-th adjacent block constraint, if its corresponding adjacent block pairs are i and j, then C[k][i] = 1, C[k][j] = -1, and the remaining elements are 0. The graph Laplace matrix describing the block adjacency relationship can be obtained from C^TC. A[i][j] = n_i^{nbr}, when i = j; A[i][j] = -1, when j ∈ 𝒩(i); A[i][j] = 0, when j ∉ 𝒩(i); Where n_i^{nbr} is the number of neighborhoods of the i-th block, and 𝒩(i) is the neighborhood set of the i-th block.

4. The method according to claim 1, characterized in that, The step of constructing a constraint vector based on the average gray-level difference of the overlapping region includes: Construct a K × 1 constraint vector r, where K is the number of adjacent block pairs; for the k-th adjacent block constraint, if its corresponding adjacent block pairs are i and j, then r[k] = R_ij; r[k] = R_ij = (1 / N_ij) × Σ_{n=1}^{N_ij} [v_j^i(n) - v_i^j(n)] Where R_ij represents the average grayscale difference between the j-th block and the i-th block in the overlapping region, N_ij is the number of pixels or voxels in the overlapping region between the i-th block and the j-th block, v_i^j(n) is the grayscale value of the n-th pixel or voxel in the overlapping region of the i-th block, and v_j^i(n) is the grayscale value of the corresponding pixel or voxel in the overlapping region of the j-th block.

5. The method according to claim 1, characterized in that, The introduction of a regularization term into the graph Laplace-type normal equation includes: Introducing the regularization coefficient α and the identity matrix I, we construct a regularized system of linear equations: (C^TC + αI) × b = C^T r Where α ranges from 1e^{-6} to 1e^{-3}, C is the linear constraint matrix of adjacent blocks, b is the bias vector to be determined, and r is the constraint vector composed of the average gray difference of the overlapping areas.

6. The method according to claim 1, characterized in that, The method of solving linear equation systems using a linear algebra solver includes: A sparse linear equation solver is used to solve the regularized normal equations to obtain the bias vector; the solver can be selected from sparse direct method, sparse LU decomposition or other equivalent sparse linear system solution methods according to the engineering implementation. When the number of blocks is large, iterative solution methods or preprocessing solution methods suitable for sparse positive definite linear systems can be used as optional implementation methods.

7. The method according to claim 1, characterized in that, The grayscale correction process for each block based on the calculated bias parameters includes: For each block, the corresponding bias parameter is uniformly added to the grayscale values ​​of all pixels. I_k^{corrected}(x, y) = I_k^{original}(x, y) + b_k Where I_k^{original}(x, y) is the pixel gray value of the k-th block at coordinates (x, y) before correction, b_k is the bias parameter of the k-th block, and I_k^{corrected}(x, y) is the pixel gray value after correction.

8. The method according to claim 1, characterized in that, The input image is either a two-dimensional image or three-dimensional volume data. When the input image is a two-dimensional image, the block processing adopts a two-dimensional uniform grid division, and the adjacency relationship adopts the four-neighbor or eight-neighbor determination criterion. When the input image is three-dimensional volume data, the block processing adopts three-dimensional uniform grid division, the adjacency relationship adopts the six-neighbor or twenty-six-neighbor determination criterion, and the overlapping area is a three-dimensional voxel block.

9. The method according to claim 1, characterized in that, The method is applicable to image stitching after deep learning block inference, brightness correction of multi-region acquired images, or global bias correction of high-throughput fluorescence microscopy multi-field scanning imaging.

10. The method according to claim 1, characterized in that, Before calculating the number of pixels and grayscale difference information in each overlapping region, the method further includes: Before calculating the number of pixels and grayscale differences in each overlapping region, at least one of the following preprocessing operations is performed on the overlapping region: background subtraction, outlier filtering, or effective pixel selection.

Images