Remote sensing image color consistency optimization method based on multi-feature constraint area network adjustment

Abstract

This invention discloses a method for optimizing the color consistency of remote sensing images based on multi-feature constrained regional network adjustment. The method involves acquiring a set of remote sensing images to be stitched, first extracting overlapping regions between the images, constructing a multi-dimensional relative constraint that integrates HOG texture, pixel mean, and CDF distribution features, and introducing a global virtual control value and a regularization penalty term as absolute constraints to jointly form a system of linear equations. Then, the least squares method is used to solve the system of linear equations to obtain the correction parameters corresponding to each remote sensing image, thereby obtaining the corresponding quadratic correction function. The quadratic correction function is then used to perform color correction on the corresponding remote sensing image to obtain the corresponding global image. Finally, Wallis transform is used to perform pixel-by-pixel local enhancement on the global image to obtain the final corrected image.

Description

Technical Field

[0001] This invention belongs to the technical field of image processing, specifically relating to a method for optimizing the color consistency of remote sensing images based on multi-feature constrained regional network adjustment. Background Technology

[0002] With the continuous development of Earth observation technology, large-scale, high-resolution mosaics generated by stitching together multiple remote sensing images are playing an increasingly important role in fields such as geographic mapping, ecological monitoring, and disaster assessment. However, due to differences in factors such as solar altitude angle, atmospheric conditions, sensor viewing angle, and acquisition time, there are usually significant radiometric differences between the original remote sensing images. These radiometric differences are mainly manifested visually as inconsistencies in tone and significant differences in brightness, thus affecting the overall visual harmony and realism of the mosaic images. Therefore, effective color consistency optimization processing before image stitching is a key step in ensuring image quality and consistency.

[0003] Existing color correction methods can be broadly categorized into three types based on their scope and strategy: direct methods, path propagation methods, and global optimization methods. Direct methods (such as histogram matching and Wallis transform) map the tones of other images to a reference standard by selecting a reference image; however, the selection of the reference image lacks a unified standard and is prone to introducing color casts. Path propagation methods construct topological relationships between images (such as minimum spanning trees or shortest paths) to sequentially transmit color correction information. While this avoids dependence on a single reference, its sequential transmission mechanism easily leads to error accumulation and instability. To overcome the limitations of the above methods, global optimization methods have emerged. These methods construct a unified optimization model to simultaneously solve for the correction parameters of all images. Early works often used linear models for mapping, which struggled to handle complex nonlinear radiometric differences. Recent research focuses on improving local fidelity through more complex models or strategies, for example... To guide optimization, methods can be introduced such as spline curves, local enhancement, or the use of specific ground features (e.g., NDVI) and gradient information (e.g., variational constraints). For example, in 2020, Liu Shijie proposed a multi-image color equalization method based on regional network adjustment. This method first uses OpenMP multi-threading technology to accelerate the extraction of the Histogram of Oriented Gradients (HOG) features from the images, and then uses this to select high-confidence corresponding grids to establish robust image connectivity. Subsequently, a rigorous regional network adjustment model is constructed, including baseline image locking, virtual control point constraints, and special constraints for difficult areas such as water bodies, thereby preventing parameter drift while unifying the color baseline. Then, observations are adaptively weighted according to the quality of the corresponding grids, and the optimal color equalization parameters are solved through an iterative algorithm. Finally, these parameters are used to perform pixel-by-pixel grayscale transformation on the original images, effectively eliminating radiometric differences between multiple images and achieving uniform color across the entire survey area. The main drawbacks of this global optimization algorithm are as follows: 1. In the global optimization stage, the constraints mostly rely on single or low-order statistics (such as the mean). When faced with large-scale, nonlinear color differences caused by terrain undulations or complex features, the stability and adaptability of the model are insufficient.

[0004] 2. In existing global-local two-level strategies, the local enhancement step is often insufficient to fully recover the contrast loss and texture blur caused by global smoothing, resulting in unnatural transitions in the final mosaic. Summary of the Invention

[0005] The purpose of this invention is to provide a method for optimizing the color consistency of remote sensing images based on multi-feature constrained regional network adjustment. By combining the global optimization adjustment concept with the local adaptive enhancement algorithm, the problem of uneven color in image stitching is solved.

[0006] To achieve the above objectives, the present invention provides the following technical solution: A method for optimizing the color consistency of remote sensing images based on multi-feature constrained regional network adjustment includes the following steps: Step 1: Perform global multi-feature regional network adjustment The process involves acquiring a set of remote sensing images to be stitched together, first extracting overlapping regions between the images, constructing a three-dimensional relative constraint that integrates HOG texture, mean, and CDF distribution features, and then introducing a global virtual control value and a regularization penalty term as two-dimensional absolute constraints to form a system of linear equations. The least squares method is then used to solve the system of linear equations to obtain the correction parameters for each remote sensing image, thereby obtaining the corresponding quadratic correction function. Finally, the quadratic correction function is used to perform color correction on the corresponding remote sensing image to obtain the corresponding global image. Step 2: Perform local Wallis enhancement The Wallis transform is used to perform pixel-by-pixel local enhancement on the global image, thereby obtaining the final corrected image.

[0007] Furthermore, for remote sensing image pairs with overlapping areas, when calculating the relative constraints of HOG textures, the overlapping areas are first uniformly divided into multiple grids, and the gradient direction histogram (HOG) and gradient distance corresponding to each grid are calculated. Then, using the gradient distance as a parameter, grids with the same name are selected. The grids with the same name are further selected according to the selection strategy. Finally, the difference in the correction value corresponding to the gray mean of the retained grids with the same name is used as the relative constraint of the HOG texture.

[0008] Furthermore, when filtering grids with the same name, the gradient distance corresponding to each grid is first calculated. Gradient distance Grids exceeding the set threshold are recorded as grids with the same name, and then the following filtering strategy is applied: a) Divide the corresponding grids within each overlapping region according to the gradient distance. Sort the data in ascending order, select the top 20% of the same-named grid cells, and retain those that satisfy the gradient distance. The same grid; b) Calculate the average gradient of each grid cell obtained in step a). , Only the average gradient is retained. >1 and The same grid with a value greater than 1, in which , These represent pairs of images from remote sensing. The average gradient corresponding to each image; c) Statistical analysis of the mean grayscale value of the corresponding grid obtained in step b) , And calculate the standard error corresponding to the difference. Remove values ​​with a difference greater than three times the standard error. The same name grid, in which , These represent pairs of images from remote sensing. The average grayscale value corresponding to each image.

[0009] Furthermore, for remote sensing image pairs with overlapping areas, when calculating the relative constraint of pixel mean, the pixel mean corresponding to the overlapping area is calculated separately, and the difference between the correction values ​​corresponding to the pixel mean of the two remote sensing images is used as the relative constraint of pixel mean. When calculating the relative constraints of CDF distribution features, the cumulative histogram CDF corresponding to the overlapping region is calculated separately. M quantiles are selected from each cumulative histogram CDF, and the difference between the correction values ​​corresponding to the pixel values ​​of each quantile in the two remote sensing images is used as the relative constraints of CDF distribution features.

[0010] Furthermore, when calculating the absolute constraints, the average brightness value of all pixels in each remote sensing image before correction is first calculated. Then, the average of all brightness values ​​is taken to obtain the global virtual control value. The difference between the correction value corresponding to the average brightness value of each remote sensing image and the global virtual control value is used as the global constraint. Calculate the cumulative histogram (CDF) of each remote sensing image before correction, sort them in ascending order, select the pixel samples in the top 30% distribution, and then use the difference between the corrected value and the original value of each pixel sample in each remote sensing image as the channel constraint.

[0011] Furthermore, the least squares method employs iterative reweighted least squares (IRWLS). Two remote sensing images with overlapping regions are denoted as remote sensing image pairs. When solving the linear equation system, the mean constraint weights corresponding to each remote sensing image pair are first calculated separately. Then, these mean constraint weights are arranged into a diagonal matrix as the mean constraint weight matrix. Similarly, the HOG grid constraint weight matrix, CDF quantile constraint weight matrix, virtual global control weight matrix, and channel differentiation weight matrix are calculated. Finally, the mean constraint weight matrix, HOG grid constraint weight matrix, CDF quantile constraint weight matrix, virtual global control weight matrix, and channel differentiation weight matrix are arranged into a diagonal matrix as the weight matrix for the least squares method.

[0012] Furthermore, the mean constraint weights are calculated using the following formulas. HOG same-name grid constraint weights CDF quantile constraint weights Virtual global control weights Channel Differential Weights ; in, Indicates remote sensing image pairs ( The total number of pixels in the overlapping region of ) This represents the set of all remote sensing image pairs within the survey area; in, Indicates remote sensing image pairs ( The total number of overlapping grid cells with the same name in the region; Where M represents the total number of quantiles corresponding to the relative constraints used to calculate the CDF distribution characteristics; in, Indicates the first The variance of the global mean estimate for a remote sensing image. Indicates a hyperparameter; .

[0013] Furthermore, the three-dimensional relative constraints and two-dimensional absolute constraints corresponding to each remote sensing image are first combined into a five-dimensional constraint, which is then used as a linear equation system. All the linear equation systems are then constructed into a super-large linear equation system. Finally, the least squares method is used to solve the super-large linear equation system to obtain the correction parameters corresponding to each remote sensing image.

[0014] Furthermore, the mathematical model corresponding to the quadratic correction function is as follows: in, and These represent the color values ​​before and after correction, respectively. This represents the correction parameter to be determined.

[0015] Compared with the prior art, the beneficial effects of the present invention are: (1) Integrating multi-dimensional image features effectively overcomes model instability under complex terrain. Unlike the traditional method of color adjustment that only uses a single feature, this invention innovatively integrates HOG texture features, CDF features, and brightness mean constraints. This multi-dimensional feature constraint mechanism breaks the limitations of a single feature, comprehensively captures the radiation and structural information of the image, and significantly improves the stability and adaptability of the model when facing complex terrain.

[0016] (2) Constructing a global quadratic optimization model based on the idea of ​​regional network adjustment effectively alleviates the accumulation of errors in long paths. This invention designs a quadratic adjustment algorithm with multiple feature constraints, which breaks through the limitations of traditional linear models in processing large-scale images. Through joint calculation in the global scope, this strategy can more accurately fit the nonlinear radiation differences caused by complex lighting and atmospheric conditions, effectively alleviate the problem of error accumulation in the long path stitching process, and fundamentally improve the robustness and final accuracy of color correction.

[0017] (3) Adaptive Wallis filtering is introduced to achieve high-fidelity image detail compensation. Addressing the technical shortcomings of traditional global smoothing color processing, such as texture blurring and reduced contrast, this invention specifically introduces an adaptive Wallis filter for targeted compensation during the local optimization stage. This method accurately restores and enhances the local texture details of the image without compromising global color consistency, ensuring high fidelity and visual quality in the final output. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the overall framework of the present invention; Figure 2 and Figure 3 The diagram illustrates the optimization results of the consistency optimization method of the present invention on datasets Data1 and Data2, respectively. Detailed Implementation

[0019] To make the technical means, creative features, objectives and effects of this invention easier to understand, the following embodiments, in conjunction with the accompanying drawings, specifically illustrate the remote sensing image color consistency optimization method based on multi-feature constrained regional network adjustment of this invention. It should be noted that the description of these embodiments is for the purpose of helping to understand this invention, but does not constitute a limitation of this invention.

[0020] To address the characteristics of multi-source remote sensing data, this invention proposes a method for optimizing the color consistency of remote sensing images based on multi-feature constrained regional network adjustment. For multiple multispectral remote sensing images with significant color differences, a global multi-feature regional network adjustment is first performed. By extracting the correspondence of overlapping areas in the images, a multi-dimensional relative constraint integrating HOG texture, mean, and CDF distribution features is constructed. Virtual control points and regularization terms are introduced as absolute constraints. An adaptively weighted IRWLS algorithm is used to robustly solve the quadratic polynomial color uniformity model, thereby eliminating the main radiometric differences between images and unifying the global color tone. Subsequently, local Wallis enhancement is performed. The processing order is established based on breadth-first search, and the image variation coefficient guides adaptive content segmentation to avoid grid effects. After local statistical interpolation, a pixel-by-pixel Wallis transformation is performed to further optimize local contrast and texture details while maintaining global color consistency.

[0021] Specifically, the method of this invention includes two stages: The first stage is global multi-feature region network adjustment, which extracts the correspondence of overlapping regions, constructs relative constraints including HOG texture features, mean features, and CDF distribution features, and combines the absolute constraints of virtual control points and regularization terms, using adaptive weight design and robust IRWLS solution to obtain the parameters of the quadratic polynomial color balancing model; the second stage is local Wallis enhancement, which establishes the reference image order by constructing an adjacency graph through BFS, performs content adaptive segmentation based on the coefficient of variation calculation, and performs pixel-by-pixel Wallis transformation after local statistical interpolation. Therefore, the method of this invention can effectively eliminate color differences in large-scale image stitching while maintaining the local texture details and global color consistency of the image.

[0022] Phase 1: In the global phase, let the set of remote sensing images to be corrected, i.e., the set of remote sensing images to be stitched, be... For each image This invention defines a correction function to be estimated for each color channel. Furthermore, by constructing observation equations through "color correspondence / statistical feature correspondence" within overlapping areas of adjacent images, an overdetermined set of equations is formed for global solution, thereby avoiding the error accumulation and reference drift problems of local propagation methods. This idea is consistent with the mainstream paradigm of "color consistency processing based on global optimization," which treats the correction parameters of all images as global variables, obtains the global optimal solution by minimizing the energy function of the difference between adjacent corrected images, and finally corrects each image individually. 1. Quadratic Correction Model: The traditional linear model ( Although the calculation is relatively simple, it is difficult to simulate radiation differences when facing complex geographical environments or when there are large overall color differences in the image. To improve the model's fitting ability and avoid overly complex parameters, we independently choose a quadratic function model for each color channel: in, and These are the color values ​​before and after correction. These are the model correction parameters to be determined. The global objective is to solve for the parameter vectors of all images simultaneously in a single global solution.

[0023] We extend the concept of "consistent geometric positions of corresponding points" in regional network adjustment to "consistent multi-dimensional color features in overlapping areas." To construct a more robust observation equation, we introduce three complementary feature dimensions for constraint: (1) pixel mean, used to determine the overall brightness baseline; (2) HOG features, which are insensitive to illumination changes and used to match and constrain regions with similar gradient structures, especially in weakly textured regions where they are more reliable than pure brightness; (3) CDF quantiles, used to constrain the overall shape of color distribution and ensure consistency between contrast and dynamic range. The combination of these three features ensures the fidelity of color mapping from different levels.

[0024] For overlapping regions Any remote sensing image pair ( ), its corrected pixel mean and To minimize the difference, and to avoid strong correlation between RGB channels, this invention requires first converting the remote sensing image to the YCbCr color space, and then taking the pixel mean of each of the three channels. The calculation method for the three channels is the same, and the pixel mean of the overlapping area is taken. / Defined as: for overlapping regions The average of all valid pixel values ​​within the image is used to construct image pairs. The basic radiation reference constraint is defined by equation (1). According to equation (1), this constraint can be written as: Since simple mean constraints can only constrain the brightness information of an image, they cannot guarantee more color information. Therefore, to establish a finer and more robust correspondence within the overlapping region, we introduce Histogram of Oriented Gradients (HOG) features. HOG is a robust similarity measure that is relatively insensitive to color changes and has good invariance to local geometric and lighting changes. This invention introduces HOG as a local texture description, actively searching and matching homonymous grids with similar internal structures within the overlapping region. After establishing the grid, a more robust matching constraint is constructed by constraining the mean of the corresponding grids. Specifically: (1) Grid division and HOG description: a) Calculate image gradient information Calculate the overlapping region image pairs separately. The gradient information includes the gradient magnitude G(x,y) and the gradient direction α(x,y).

[0025] Gradient magnitude: Gradient direction: in, Represents pixels within the overlapping area The corresponding pixel value, Represents the coordinates of a pixel.

[0026] b) Histogram of Orientation Gradients in Gridded Images Divide the overlapping region into multiple non-overlapping grids, with each grid size set to [size missing]. For each pixel, divide the grid into 36 bins (10° per bin) with a 360° angular range. For each pixel within the grid, weighted summations are performed across the 36 bins based on its gradient direction α. ​​The weight is the gradient magnitude at that point. After processing all pixels within the grid using the above method, the gradient direction histogram of the current grid is obtained.

[0027] (2) Gradient distance and same-name grid selection: To measure the similarity of grid content, gradient distance is used. An evaluation should be conducted. Generally speaking, The smaller the value, the more similar the grid content. When GD is less than the threshold T, it is judged as a "same-name grid". Constraints are only established on the selected set of same-name grids (a consistency constraint is applied to the mean of the same-name grids), thereby avoiding incorrect matching caused by weak textures / shadows / variable regions and improving the reliability and robustness of the constraints.

[0028] In the formula, Indicates the grid number, This indicates the direction bin number, which is 36 here; This represents the gradient direction histogram of the k-th grid with the same name.

[0029] Due to the diverse characteristics of ground features, texture differences can exist, such as interference from cloudy or foggy conditions. Simultaneously, weakly textured regions with low gradient sizes (such as lakes or cloud areas) make gradient similarity susceptible to noise and unreliability, making it difficult to guarantee that ground features in overlapping areas are consistent, which can affect subsequent adjustment. Therefore, in practical applications, we screen these homonymous grids using the following strategy: a) Divide the grids of the same name according to gradient distance Sort in ascending order, take the top 20%, and retain those that meet the requirements. The same grid; b) Calculate the average gradient of the grid with the same name. , Only retain >1 and A grid with a gradient greater than 1 is defined as follows: the average gradient represents the average gradient magnitude of all pixels within the grid. Since the average gradient value is less than 1, it indicates that the gradient change in this region is relatively slow. Generally, this region is a cloud, a water body, or another region with indistinct gradient characteristics. Such regions have low reliability and stability and should be removed before processing. c) Even after following the above steps for filtering, the results may still contain local outliers. To ensure the consistency of the filtered grid content, further removal is necessary. First, calculate the average gray value of each remaining grid with the same name. , To find the mean square error of the difference, the formula is as follows: in, This represents the number of grid cells with the same name remaining after the first two removal steps. Let Δ be the mean error, and Δ be the difference in the mean gray values ​​of the remaining grid cells with the same name. According to error theory, the difference Δ is greater than three times the mean error. The probability of random errors is only 3%. Therefore, when the amount of data is limited, applying a standard error of three times the difference between the grayscale mean and the remaining grid cells with the same name can further eliminate local outliers caused by mismatches.

[0030] (3) By implementing mean consistency constraints on these rigorously selected local regions with similar content, the reliability and physical rationality of the constraints are greatly enhanced, thereby ensuring that the model can still obtain accurate relative radiation relationships in complex surface environments: in, , These represent the average gray values ​​of the corresponding grid cells with the same name in the overlapping area after strict screening.

[0031] The two constraints mentioned above only constrain brightness and image gradient information, insufficient to guarantee complete consistency in tonal distribution. To ensure consistency in contrast and dynamic range for the corrected image, we introduce a concept based on the cumulative distribution function (CDF, ​​or cumulative histogram of image values). Specifically, for each image pair, we first calculate the CDF within the overlapping region, then use the equal probability quantiles corresponding to the two CDFs as the color correspondence. We select M quantiles (e.g., 0.1, 0.2, ..., 0.9) uniformly within the interval; here, we default M to 16. Let... and Image pairs ( In overlapping areas The pixel value corresponding to the m-th quantile. and These are all corrected pixel values, and the pixel values ​​at all 16 quantiles need to be included in the constraint calculation. Our constraint requires that these values ​​remain equal after correction, so that their corrected CDF distributions are similar. The above constraint can be described as: 2. In the aforementioned global adjustment model, if only relative constraints between images (such as overlap mean constraints, HOG corresponding grid constraints, CDF quantile constraints, etc.) are used, two typical problems will arise: when there are weakly connected regions in the survey area (low overlap rate, sparse connecting edges, and many edge images), the normal equations are prone to rank deficiency or ill-conditioned conditions; even if the equations are solvable, the lack of an absolute reference may lead to free drift in the overall solution space, manifested as "overall color drift". To improve the stability of the solution and enhance fidelity, this invention introduces two types of absolute constraints for stabilization: virtual global control and regularization penalty terms.

[0032] 1) We introduce a virtual global control as an absolute constraint to provide a uniform radiometric reference for the entire block. This constraint enforces a uniform mean luminance across all image-corrected images. They all converge to a virtual global value determined by the average initial brightness of all images. .

[0033] Virtual global value The determination method is as follows: For all original images, i.e., each remote sensing image before correction, first calculate the average brightness value corresponding to all pixels, which is the arithmetic mean of the Y channel components in the YCbCr space, denoted as . Then, these values ​​are averaged to obtain the virtual global value: This constraint not only effectively regularizes the correction function to suppress color drift, but more importantly, it ensures the solvability of the equation system, which is key to obtaining a stable and reliable solution. We can obtain the following equation: Represents the original image The arithmetic mean of the |Y channel components of all pixels is the average brightness. This indicates that the mean luminance is adjusted using a quadratic correction function. Perform the correction and obtain the correction value; 2) To prevent the correction function from overfitting in non-overlapping regions and to ensure radiometric fidelity, this invention introduces a regularization penalty term. We abandon uniform sampling, which cannot represent the dominant color tone of the image, and instead adopt an adaptive sampling strategy based on histograms: from each original image... The cumulative histogram (CDF) distribution takes the high-density regions in the top 30%, which represent the pixel values ​​of the dominant color tone. By forcing the corrected values ​​of these samples to approximate their original values. This effectively preserves the core color characteristics of the image. Taking a set of samples is to constrain the quadratic function simultaneously at multiple grayscale positions of the main color tone, avoiding local satisfaction and overall drift caused by using only a single point. Multiple samples are equivalent to applying bandwidth constraints to the main color tone range, making the correction smoother and more stable.

[0034] Finally, we combine the absolute control (6) and (7) regularization constraints with the relative constraints (2), (4), and (5) to form a large, overdetermined system of linear equations for global solution. The simplified form is: in, All parameters to be determined The resulting long vector is the correction parameter abc in the above formulas. It is a coefficient matrix. It is a vector of observations (or constant terms). It is the error vector, and its ultimate goal is 0.

[0035] 3. Robust Solution with Iterative Re-weighted Least Squares (IRWLS) Because this invention incorporates multiple constraints with different properties, a reasonable weighting mechanism is needed to avoid bias in the solution towards a particular type of constraint and ensure the acquisition of the overall optimal solution. This mechanism aims to harmonize the influence of each constraint, preventing the model solution process from favoring a specific constraint, thereby ensuring the acquisition of the overall optimal solution. This invention employs the following weighting design: 1) Mean constraint weight: We calculate the number of pixels between overlapping regions, therefore, we weight this constraint term. Set as the area of ​​the overlapping region This is the normalized percentage of the total area of ​​all overlapping regions, expressed in pixels. The specific calculation method for this weight is defined as follows: in, Indicates image pair ( The total number of pixels in the overlapping area. This is the set of all overlapping image pairs within the survey area. The specific format is as follows: 2) HOG Corresponding Grid Constraint Weights: When constructing the correspondence between corresponding HOG grids, considering the geometric and radiometric inconsistencies between adjacent images, and to objectively quantify the uncertainty of this matching relationship, we will assign weights to a pair of images (…). The total number of overlapping grid cells with the same name between ) Defined as the "connection strength" between the two. Accordingly, we introduce a specific weight for the constraints of each overlapping region. The weight is set as the normalized ratio of the number of grid cells with the same name in the area to the total number of all grid cells with the same name in the survey area. The formula for calculating the weight is as follows.

[0036] in, This represents the set of all overlapping image pairs within the survey area. The specific matrix form is as follows: 3) CDF Quantile Constraint Weights: For the CDF constraint, we assign equal weights to the M quantile matching terms. This ensures that each sampling point on the CDF curve contributes equally to the alignment of the overall distribution. Therefore, the weight of this constraint term is 1 / M. The mathematical expression of this constraint term is: The specific matrix form is as follows: 4) Virtual Global Control Weights: For virtual global control constraints, to differentiate the constraint strength based on the confidence level of the global statistics for each image, we designed a weight related to the statistical variance. Following common practice, the weight of this constraint is set as the inverse of the variance of the global mean estimate, and adjusted by a factor. Scaling is performed. Global weights of images The calculation is as follows: in, It is the first The variance of the global mean estimate for the image. As a hyperparameter, it is used to adjust the relative strength of the global constraint term in the overall objective function. The specific matrix form is as follows: 5) Channel Differentiation Weights: For regularization constraints, the human eye is more sensitive to changes in luminance (L channel) than to changes in chroma (a, b channels). Based on this prior knowledge, we assigned differentiated weights to different channels, applying a higher weight to the regularization of the luminance channel. A lower weight was applied to the chroma channel. The specific calculations are as follows: The specific matrix form is as follows, with corresponding weight matrices used for different channels: For the different weighting designs mentioned above, the six diagonal matrices are combined into a total weight matrix (all values ​​except for the diagonal within the block are 0). Therefore, the overall weight matrix... It can be described as: When processing large-scale remote sensing imagery, overlapping areas inevitably contain anomalous features such as clouds, cloud shadows, and water reflections, which can create erroneous constraints (i.e., gross errors). Traditional least squares estimation aims to minimize the L2 norm of the residuals, a characteristic that makes it extremely sensitive to gross errors; a single gross error can cause a significant deviation from the entire solution. To obtain robust solutions, we introduce the Iterative Reweighted Least Squares (IRWLS) framework.

[0037] By using this weight matrix Integrating into the IRWLS framework, in optimization theory, the weighted least squares method aims to solve for minimizing the quadratic residual form. Correction parameters The closed-form solution for x is: IRWLS solves the least squares problem through an iterative process. In each iteration, the algorithm first calculates the residuals of all constraint equations based on the currently estimated parameters. Subsequently, it assigns a new weight to each residual using a robust weighting function (Huber's M-estimator). The weighting function is designed to assign a very small weight to constraints with large residuals, and a larger weight to "normal" constraints with small residuals. These new weights are then used to update the weighted least squares problem and initiate the next iteration. Essentially, this strategy seamlessly integrates parameter estimation and outlier diagnosis within a single iterative loop. By dynamically "removing" unreliable observations (large residual terms) from the cost function, IRWLS converges to a solution dominated solely by high-quality observations, resulting in highly accurate and robust global color correction parameters.

[0038] Phase 2: Local Wallis Enhancement This stage aims to further optimize local details and eliminate residual color differences. Through local enhancement strategies, the final image maintains the seamless macroscopic color consistency achieved in the first stage, while its microscopic texture details and local contrast are effectively restored and enhanced in the second stage.

[0039] in, and These are points in the target image (the final image obtained) and the globally corrected image (the image after the first stage of correction). Pixel value at that location, and These represent the standard deviations of pixel values ​​in the target image and the globally corrected image, respectively. and These are the mean pixel values ​​of the target image and the globally corrected image, respectively. In this method, the globally corrected image is first divided into uniform blocks, and the standard deviation and mean of each block are calculated. To avoid block artifacts, these statistics are used to generate smooth, spatially varying transform parameters for each pixel through bilinear interpolation. Finally, each pixel is subjected to an independent Wallis transform based on the interpolated parameters for its location. Through this local enhancement strategy, the final image maintains seamless macroscopic color transitions while effectively restoring its microscopic texture details and local contrast.

[0040] The specific process is as follows: 1) Reference image and processing order determination (BFS) Each image to be processed is treated as a node in a graph, and an adjacency graph is constructed based on the spatial adjacency / overlap relationships between images. Breadth-First Search (BFS) is then performed on this graph to obtain a node visit sequence, which is used as the image processing order. Simultaneously, for each image in the sequence, its predecessor node on the BFS visit path is taken as the "current reference image." To establish the "initial reference image," the most suitable image for use as the global baseline is selected from all images to be processed. A typical selection method is to calculate the image with the largest average gradient for each image as the original reference. The statistics of the "current reference image" and the "original reference image" are then fused at a certain ratio to obtain the target reference mean and reference standard deviation for that image, which are used for overall alignment during subsequent enhancement.

[0041] 2) Content-adaptive block calculation For each image to be processed, global statistics (mean and standard deviation) are calculated, and the coefficient of variation (CV) is obtained from this to characterize the strength of image content / brightness variations. Based on the image CV, the block size and row / column block division are adaptively determined.

[0042] 3) Calculate local statistics For the previously segmented image blocks, calculate the mean and standard deviation within each block, and further calculate the mean and standard deviation of the block's corner points (corner points are obtained by weighting the contributions of multiple surrounding blocks; if a block's statistics are abnormal or zero, it can be excluded from corner point estimation); the local mean of the pixel is obtained by interpolating the statistics of adjacent corner points / blocks using bilinear interpolation. With local standard deviation .

[0043] 4) Pixel-by-pixel Wallis correction Finally, a pixel-by-pixel Wallis transform is performed on each pixel to transform the original pixel value. Map to new pixel values ​​using the method of "subtracting the local mean, scaling by the standard deviation, and adding the reference mean". , that is, use Controlling local contrast stretching or compression, using By controlling the brightness shift, a processing effect can be achieved that maintains global style consistency while locally adaptively enhancing details and suppressing uneven lighting. The entire process advances multiple images one by one in BFS order, so that the reference relationship propagates along the adjacency graph, ultimately resulting in an enhanced result with more consistent brightness / contrast and clearer details in large-scale stitching scenes.

[0044] To quantitatively evaluate the effectiveness of the method of the present invention, we use three indicators to quantitatively evaluate the quality of the color correction results.

[0045] Structural Similarity (SSIM): SSIM measures the similarity between the processed image and the original image in three dimensions: structure, contrast, and brightness. The closer the value is to 1, the less image distortion and the higher the color fidelity.

[0046] and Each represents one of two images. , and These represent the three components: brightness, contrast, and structure.

[0047] Gradient Loss (GL): GL is used to calculate the gradient changes in the processed image. The smaller the GL value, the less the color correction process damages the texture structure of the original image.

[0048] in It is the source image. This indicates the corrected result. Calculate the difference between pixel orientation maps. yes The number of pixels.

[0049] One-dimensional entropy (OEI): OEI represents the amount of information contained in the grayscale distribution of an image. The larger the OEI value, the richer the grayscale levels contained in the image and the stronger the detail representation.

[0050] Where pk represents the frequency of pixel values ​​k. The larger the OEI, the greater the amount of information, and the richer the grayscale of the image.

[0051] Experimental environment: All experiments were conducted under Windows 10 environment. CPU: Intel Core i9-10920X, GPU: NVIDIA GeForce RTX 2080, Memory size: 64GB, Programming language: C++11, OpenCV version: 3.4.5.

[0052] 4.3 Technical Effects Compared with current advanced color optimization algorithms, the method of this invention demonstrates superior performance. Figure 2 and Figure 3 The table below shows the before-and-after optimization results of the algorithm on different datasets, demonstrating good color optimization performance. Quantitative metrics are shown in the table below. It is important to note that the schemes and arrangements of this application shown in the exemplary embodiments are merely exemplary. Although only a few embodiments are described in detail in this disclosure, those who consult this disclosure will readily understand that many modifications are possible (e.g., variations in various parameter values ​​(temperature, power, humidity, etc.), installation arrangements, names, colors, logical orders, etc.) without substantially departing from the novel teachings and advantages of the subject matter described in this application. Therefore, all such modifications are also included within the scope of the invention, and the order or sequence of any process or method steps may be changed or rearranged according to alternative embodiments. In the claims, any "apparatus plus function" clause is intended to cover the structure described herein for performing the function, and not only structurally equivalent but also equivalent in structure. Other substitutions, modifications, alterations, and omissions may be made in the design, operation, and arrangement of the exemplary embodiments without departing from the scope of the invention. Therefore, the invention is not limited to the particular embodiments but extends to a variety of modifications that still fall within the scope of the appended claims.

[0053] Furthermore, in order to provide a concise description of exemplary embodiments, not all features of actual embodiments (i.e., those features that are not relevant to the best mode of carrying out the invention as currently considered, or those features that are not relevant to implementing the invention) may be omitted.

[0054] It should be understood that numerous specific implementation decisions can be made during the development of any practical implementation, such as in any engineering or design project. Such development efforts may be complex and time-consuming, but for those skilled in the art who benefit from this disclosure, the development effort will be a routine work of design, manufacturing, and production without requiring much experimentation.

[0055] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for optimizing the color consistency of remote sensing images based on multi-feature constrained regional network adjustment, characterized in that... Includes the following steps: Step 1: Perform global multi-feature regional network adjustment The process involves acquiring a set of remote sensing images to be stitched together, first extracting overlapping regions between the images, constructing a three-dimensional relative constraint that integrates HOG texture, mean, and CDF distribution features, and then introducing a global virtual control value and a regularization penalty term as two-dimensional absolute constraints to form a system of linear equations. The least squares method is then used to solve the system of linear equations to obtain the correction parameters for each remote sensing image, thereby obtaining the corresponding quadratic correction function. Finally, the quadratic correction function is used to perform color correction on the corresponding remote sensing image to obtain the corresponding global image. Step 2: Perform local Wallis enhancement The Wallis transform is used to perform pixel-by-pixel local enhancement on the global image, thereby obtaining the final corrected image.

2. The method for optimizing color consistency of remote sensing images based on multi-feature constrained regional network adjustment according to claim 1, characterized in that: For remote sensing image pairs with overlapping areas, when calculating the relative constraints of HOG textures, the overlapping areas are first uniformly divided into multiple grids, and the gradient direction histogram (HOG) and gradient distance corresponding to each grid are calculated. Then, using the gradient distance as a parameter, grids with the same name are selected. The grids with the same name are further selected according to the selection strategy. Finally, the difference in the correction value corresponding to the gray mean of the retained grids with the same name is used as the relative constraint of the HOG texture.

3. The method for optimizing color consistency of remote sensing images based on multi-feature constrained regional network adjustment according to claim 2, characterized in that: When filtering grids with the same name, first calculate the gradient distance corresponding to each grid. Gradient distance Grids exceeding the set threshold are recorded as grids with the same name, and then the following filtering strategy is applied: a) Divide the corresponding grids within each overlapping region according to the gradient distance. Sort the data in ascending order, select the top 20% of the same-named grid cells, and retain those that satisfy the gradient distance. The same grid; b) Calculate the average gradient of each grid cell obtained in step a). , Only the average gradient is retained. >1 and The same grid with a value greater than 1, in which , These represent pairs of images from remote sensing. The average gradient corresponding to each image; c) Statistical analysis of the mean grayscale value of the corresponding grid obtained in step b) , And calculate the standard error corresponding to the difference. Remove values ​​with a difference greater than three times the standard error. The same name grid, in which , These represent pairs of images from remote sensing. The average grayscale value corresponding to each image.

4. The method for optimizing color consistency of remote sensing images based on multi-feature constrained regional network adjustment according to claim 1, characterized in that: For remote sensing image pairs with overlapping areas, when calculating the relative constraint of pixel mean, the pixel mean corresponding to the overlapping area is calculated separately, and the difference between the correction values ​​corresponding to the pixel mean of the two remote sensing images is used as the relative constraint of pixel mean. When calculating the relative constraints of CDF distribution features, the cumulative histogram CDF corresponding to the overlapping region is calculated separately. M quantiles are selected from each cumulative histogram CDF, and the difference between the correction values ​​corresponding to the pixel values ​​of each quantile in the two remote sensing images is used as the relative constraints of CDF distribution features.

5. The method for optimizing color consistency of remote sensing images based on multi-feature constrained regional network adjustment according to claim 1, characterized in that: When calculating the absolute constraints, first calculate the average brightness value of all pixels in each remote sensing image before correction, then average all the average brightness values ​​to obtain the global virtual control value, and then use the difference between the correction value corresponding to the average brightness value of each remote sensing image and the global virtual control value as the global constraint. Calculate the cumulative histogram (CDF) of each remote sensing image before correction, sort them in ascending order, select the pixel samples in the top 30% distribution, and then use the difference between the corrected value and the original value of each pixel sample in each remote sensing image as the channel constraint.

6. The method for optimizing color consistency of remote sensing images based on multi-feature constrained regional network adjustment according to claim 1, characterized in that: The least squares method employs iterative reweighted least squares (IRWLS). Two remote sensing images with overlapping regions are denoted as remote sensing image pairs. When solving the linear equation system, the mean constraint weights corresponding to each remote sensing image pair are first calculated separately. These mean constraint weights are then arranged into a diagonal matrix as the mean constraint weight matrix. Similarly, the HOG grid constraint weight matrix, CDF quantile constraint weight matrix, virtual global control weight matrix, and channel differentiation weight matrix are calculated. Finally, these matrixes are arranged into a diagonal matrix as the weight matrix for the least squares method.

7. The method for optimizing color consistency of remote sensing images based on multi-feature constrained regional network adjustment according to claim 6, characterized in that: Calculate the mean constraint weights using the following formulas. HOG same-name grid constraint weights CDF quantile constraint weights Virtual global control weights Channel Differential Weights ; in, Indicates remote sensing image pairs ( The total number of pixels in the overlapping region of ) This represents the set of all remote sensing image pairs within the survey area; in, Indicates remote sensing image pairs ( The total number of overlapping grid cells with the same name in the region; Where M represents the total number of quantiles corresponding to the relative constraints used to calculate the CDF distribution characteristics; in, Indicates the first The variance of the global mean estimate for a remote sensing image. Indicates a hyperparameter; 。 8. The method for optimizing color consistency of remote sensing images based on multi-feature constrained regional network adjustment according to claim 6, characterized in that: First, the three-dimensional relative constraints and two-dimensional absolute constraints corresponding to each remote sensing image are combined into a five-dimensional constraint, which is then used as a linear equation system. Next, all the linear equation systems are constructed into a super-large linear equation system. Finally, the least squares method is used to solve the super-large linear equation system to obtain the correction parameters corresponding to each remote sensing image.

9. The method for optimizing color consistency of remote sensing images based on multi-feature constrained regional network adjustment according to claim 1, characterized in that: The mathematical model corresponding to the quadratic correction function is as follows: in, and These represent the color values ​​before and after correction, respectively. This represents the correction parameter to be determined.

Images