A Vignetting Correction Method for Remote Sensing Images Based on Low-Rank Modeling and Polynomial Fitting

1. Introduction

Vignetting is a common radiometric nonuniformity in imaging systems, typically manifested as a gradual attenuation of image brightness from the center toward the edges. For remote sensing images, this spatial radiometric nonuniformity not only degrades visual quality but also affects radiometric calibration, land surface parameter retrieval, target detection, and consistency in multi-scene image mosaicking [1,2]. In wide-swath imaging, multi-detector mosaicking, and thermal infrared remote sensing, vignetting, striping, and response drift are often coupled, making image nonuniformity more pronounced [8,9]. However, in the absence of vignetting-free ground-truth images or dedicated calibration data, vignetting correction usually has to be performed under non-reference conditions. In such cases, existing methods can only rely on single-image statistical priors or multi-image shared constraints to estimate the vignette field, making them more susceptible to interference from scene content, texture structures, striping, and response drift, which in turn affects the stability of vignette recovery and correction accuracy. Therefore, developing a stable, generalizable, and remote-sensing-oriented vignetting correction method remains of clear research value.

Image vignetting correction methods can be roughly divided into three categories. The first category consists of methods based on physical calibration or radiometric calibration, which recover the vignette distribution through exposure modeling, camera response estimation, flat-field calibration, or on-orbit calibration [8,9,10,11]. Such methods usually have good physical interpretability, but they often require additional calibration data, stable imaging conditions, or rigorous experimental procedures, resulting in high deployment costs in practical remote sensing tasks. The second category comprises methods based on single-image statistical priors, such as texture segmentation, gradient distribution symmetry, or radial gradient constraints, to estimate the vignette field [3,4]. These methods reduce dependence on external calibration data, but their estimation results are easily disturbed by scene content when image content is complex, texture structures are strong, or imaging geometry is non-ideal. The third category extends radiometric nonuniformity correction to multi-view or multi-frame scenarios and restores radiometric consistency through multi-image information [12,13,14], rather than explicitly extracting the vignette field. Existing studies mainly focus on single-image vignette fitting, physical calibration, or relative radiometric correction, while explicit modeling of a vignette structure shared by multiple images remains relatively rare, especially in remote sensing scenarios.

In imaging, vignetting is essentially a low-frequency and smooth spatial radiometric field. Therefore, polynomial models and their variants have long been important tools for vignette modeling. Classical methods usually employ two-dimensional polynomials, radial polynomials, local polynomials, or deformable radial polynomials to approximate the vignette distribution [5,6,7]. These methods offer simple models, few parameters, easy implementation, and low computational cost, making them particularly suitable for describing smooth, low-frequency brightness attenuation fields. For vignetting in remote sensing images, polynomial fitting can effectively recover large-scale brightness variation trends and is therefore highly practical in engineering applications. Nevertheless, polynomial models also have clear limitations. First, fitting from a single image is easily affected by surface textures and radiometric structures. Second, low-order polynomials are difficult to accurately describe complex off-center and asymmetric vignetting. Third, when multi-frame shared constraints are absent, the fitting result often exhibits large instability. In other words, polynomial fitting is effective in representing smooth low-frequency structures, but it is not well suited to separating the vignette field from complex image content.

Statistical analysis of long-term experimental data from space optical systems shows that, when the space environment and observation angle vary within a limited range, the vignette field usually changes slowly over a short time scale and is unlikely to undergo significant variations. Its spatial distribution and attenuation amplitude generally exhibit good stability. Therefore, it can be approximately regarded as satisfying the short-term stability assumption.

The worst-corner vignette degree is used for quantitative evaluation.

$$V_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{n}\mathrm{e}\mathrm{r},\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{s}\mathrm{t}}=1-{\eta }_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{n}\mathrm{e}\mathrm{r},\mathrm{m}\mathrm{i}\mathrm{n}}$$

(1)

where${\eta }_{\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{r}\mathrm{n}\mathrm{e}\mathrm{r},\mathrm{m}\mathrm{i}\mathrm{n}}=\min \left({\displaystyle \frac{I_{\mathrm{L}\mathrm{T}}}{I_c}},{\displaystyle \frac{I_{\mathrm{R}\mathrm{T}}}{I_c}},{\displaystyle \frac{I_{\mathrm{L}\mathrm{B}}}{I_c}},{\displaystyle \frac{I_{\mathrm{R}\mathrm{B}}}{I_c}}\right)$,$I_c={\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}}_{\left(x,y\right)\in {\Omega }_c}\widehat{I}\left(x,y\right)$,$I_{\mathrm{L}\mathrm{T}}={\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}}_{\left(x,y\right)\in {\Omega }_{\mathrm{L}\mathrm{T}}}\widehat{I}\left(x,y\right)$,$I_{\mathrm{R}\mathrm{T}}={\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}}_{\left(x,y\right)\in {\Omega }_{\mathrm{R}\mathrm{T}}}\widehat{I}\left(x,y\right)$,$I_{\mathrm{L}\mathrm{B}}={\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}}_{\left(x,y\right)\in {\Omega }_{\mathrm{L}\mathrm{B}}}\widehat{I}\left(x,y\right)$,$I_{\mathrm{R}\mathrm{B}}={\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}}_{\left(x,y\right)\in {\Omega }_{\mathrm{R}\mathrm{B}}}\widehat{I}\left(x,y\right)$.

Here, $I_c$ denotes the median low-frequency radiance of the central region, $I_{\mathrm{L}\mathrm{T}}$ denotes the median low-frequency radiance of the upper-left region, $I_{\mathrm{R}\mathrm{T}}$ denotes the median low-frequency radiance of the upper-right region, $I_{\mathrm{L}\mathrm{B}}$ denotes the median low-frequency radiance of the lower-left region, $I_{\mathrm{R}\mathrm{B}}$ denotes the median low-frequency radiance of the lower-right region. $\widehat{I}\left(x,y\right)$ denotes the low-frequency radiance field obtained by Gaussian low-pass filtering.

As shown in Figure 1, the vignette field usually satisfies the short-term stability assumption under similar observation directions and spatiotemporal imaging conditions; that is, its spatial distribution and attenuation characteristics remain approximately unchanged within a short time window. However, when the spatiotemporal observation environment changes significantly, the vignette field may become nonstationary, and its attenuation degree and spatial morphology may vary markedly or even exhibit abrupt changes. It should also be noted that, from the perspectives of radiometric field transmission and sensor response mechanisms, vignetting is mainly a multiplicative spatial attenuation process, whereas additive terms usually only represent background bias, response drift, or local additional disturbances. Based on this understanding, this paper models the vignette field as a multiplicative-dominant common component and treats additive effects as secondary disturbances.

Studies on low-rank decomposition have shown that, when an observation matrix can be decomposed into a low-rank background and sparse disturbances, stable recovery can be achieved through convex optimization or approximate decomposition [15,16,17]. This idea has since been widely used in background modeling, shadow correction, illumination field separation, and hyperspectral anomaly detection [18,19,20]. For vignette correction, if multiple frames share similar vignette distributions, while different scene contents, local anomalies, and noise mainly appear as non-shared disturbances, then organizing multi-frame data into a matrix and applying low-rank decomposition may make it possible to extract the common vignette structure. However, low-rank decomposition itself is more focused on structure extraction and does not directly guarantee that the recovered result satisfies the physical property of vignetting as a smooth low-frequency field. Therefore, using low-rank decomposition alone often fails to obtain an ideal vignette field. Modeling and solving the shared vignette field across multiple frames under the same imaging condition remain relatively uncommon, especially in remote sensing scenarios. Meanwhile, existing studies typically emphasize parametric fitting, with limited systematic analysis of the respective roles, complementarity, and joint influence of low-rank modeling and polynomial fitting on final vignette estimation and correction performance.

In addition to affecting ground radiometric consistency, vignetting also directly influences the accuracy of star-based radiometric calibration. Stars have been used as stable calibration sources for on-orbit radiometric calibration and imaging performance assessment of remote sensing satellites, and a key prerequisite is that the same star should exhibit a consistent radiometric response under different observation conditions and at different field positions [21,22]. However, vignetting introduces field-dependent multiplicative attenuation across the image plane, causing the point-spread energy distribution and the background-subtracted DN value of a star to vary with field position. This reduces the consistency of stellar signal extraction and further affects the stability of calibration results [23,24]. For low-irradiance stellar targets in particular, background subtraction and signal-to-noise control are critical to photometric accuracy. If vignetting is not effectively suppressed, the field-dependent response variation it introduces will be directly translated into calibration error. Existing studies have shown that star-based calibration can achieve approximately 2% calibration accuracy under reasonable signal-to-noise conditions [22], which implies that vignette correction should minimize field-dependent variations in stellar energy as much as possible to improve star-based calibration accuracy.

This paper formulates remote sensing image vignette correction as a vignette-field extraction problem. Specifically, a multi-frame data matrix is first constructed in the logarithmic domain, and low-rank decomposition is used to estimate the common vignette field. Polynomial fitting is then introduced to further recover a smooth vignette field, thereby improving the physical plausibility and stability of the result. Finally, a correction gain is constructed in the original intensity domain to achieve image vignette correction. The proposed method accounts for both vignette-field extraction capability and low-frequency field modeling capability.

The main contributions of this paper are as follows: (1) based on the short-term stability of vignetting, a multi-frame remote sensing image vignette-correction framework is proposed, in which vignette-field extraction is formulated as a low-rank decomposition problem; (2) the correction method combines low-rank modeling and polynomial fitting, simultaneously exploiting multi-frame image diversity and the smooth low-frequency prior of the vignette field; and (3) experiments verify the effectiveness of the proposed method under various vignette-disturbed scenarios, and the results demonstrate that it achieves the best overall performance in terms of MAE, MAD, CenterMAE, and EdgeMAE.

2. Materials and Methods

Based on the systematic analysis in the previous section, a remote sensing image vignette correction method based on low-rank modeling and polynomial fitting is proposed. The overall method jointly considers vignette-field extraction capability and smooth-field modeling capability.

Figure 2. Algorithms.

Figure 2. Algorithms.

2.1. Datasets

The experimental data used in this study consist of three parts: real remote sensing images, simulated vignetted images, and star images. These three types of data serve different purposes in the evaluation. The real remote sensing images are used to demonstrate the applicability of the proposed method under practical imaging conditions, the simulated vignetted images are used for controlled quantitative analysis, and the star images are used to further assess the influence of vignette correction on radiometric consistency and calibration accuracy.

The real remote sensing images were collected from SDGSAT-1, Landsat 8/9, and ASTER, covering different imaging conditions and representative surface scenes. These data provide diverse image content and radiometric characteristics for evaluating correction performance under practical conditions. The simulated vignetted images were generated by superimposing vignette fields extracted from real detector radiometric images onto remote sensing radiance images. In this way, the simulated data preserve measured detector-related vignette characteristics while allowing controlled evaluation of correction performance.

The remote sensing images and the simulated vignetted images were all of size 2048 × 2048 pixels, whereas the star images were of size 1024 × 1024 pixels. As listed in Table 1, the star dataset is accompanied by auxiliary information such as stellar magnitude, right ascension, and declination, which is used to support star identification and subsequent calibration-related analysis.

2.2. Correction Method

2.2.1. Logarithmic Domain Transformation

Considering that vignetting usually acts on an image as a multiplicative factor in the imaging model, the original image I can be expressed as:

$$I=R·V+N$$

(2)

where $R$ is the ideal image, $V$ is the spatial vignette field, and$N$denotes the low-amplitude background noise.

Taking the logarithm of both sides yields the following approximate representation.

$$\log I\iff \log R+\log V$$

(3)

2.2.2. Data Matrix Construction

Assume that there are $n$ images in total, and that each image has a size of $H\times W=m$. Each image is vectorized column-wise into a column vector of length, and all images are concatenated to obtain the observation matrix $D_0$.

Since the gray-value range of the original images is large, the matrix is first normalized as $D$ to reduce dimensional effects and improve numerical stability.

2.2.3. Vignetting Field-Sparse Anomaly Separation for Multi-Frame Images

This paper models the data matrix using a low-rank-sparse decomposition model, i.e.,

$$D=L+S$$

(4)

where $L$ represents the vignette field shared by multiple images, and $S$ represents the non-shared residual term relative to the vignette field, including differences in scene content, local textures, noise, and disturbance terms.

The reason for modeling $L$ as a low-rank term is that, under the short-term stability assumption of the vignette field, different images share the same or approximately the same vignette distribution, and their vectorized columns are therefore strongly correlated. In the ideal case, the shared vignette component can be represented as $\mathit{L}=v{\mathbf{1}}^T$, corresponding to a rank-1 structure. It can therefore be approximately regarded as low-rank in practice. On the other hand, scene content and local disturbances appear as non-shared variations relative to the shared vignette component and can be further approximated as a sparse residual term S in the low-rank-sparse decomposition framework.

Where $v\in {\mathbb{R}}^{m\times 1}$denotes the estimated common background vector, and$\mathbf{1}\in {\mathbb{R}}^{m\times n}$ denotes an all-ones vector. This form means that each column of matrix $L$ is identical; that is, all images share the same vignette field. Therefore,

$$\underset{v,\mathit{s}}{\min }\lambda {‖\mathit{S}‖}_{1}s.t.\mathit{D}=v{\mathbf{1}}^T+\mathit{S}$$

(5)

Where $\lambda$ is the regularization parameter of the sparse term. In this paper, $\lambda ={\displaystyle \frac{1}{\sqrt{m}}}$ is used to balance the vignette field and image textures.

To solve the above problem, the Augmented Lagrangian Method (ALM) is adopted for iterative optimization, and the augmented Lagrangian function is constructed as:

$$\mathcal{L}\left(v,S,Y,\mu \right)=\lambda {‖\mathit{S}‖}_1+〈Y,D-v{\mathbf{1}}^T-S〉+{\displaystyle \frac{\mu }{2}}{‖D-v{\mathbf{1}}^T-S‖}_F^2$$

(6)

where $Y$ is the Lagrange multiplier, $\mu >0$ is the penalty parameter, $〈·,·〉$ denotes the matrix inner product, and ${‖·‖}_F$ denotes the Frobenius norm.

During initialization, we first set:

$$Y_0=D$$

(7)

Then, it is normalized using the joint dual norm,

$$Y_0\leftarrow {\displaystyle \frac{D}{max\left({‖D‖}_2,{‖D‖}_{\mathrm{\infty }}/\lambda \right)}}$$

(8)

Where ${‖D‖}_2$ denotes the spectral norm of the matrix, and ${‖D‖}_{\infty }=\underset{i,j}{\max }\left|D_{ij}\right|$ denotes the maximum absolute value among all matrix elements. This initialization keeps the multiplier variable within an appropriate numerical scale and improves the stability of subsequent iterations. The penalty parameter is initialized as:

$${\mu }_0={\displaystyle \frac{1.25}{{‖D‖}_2}}$$

(9)

and its upper bound is set as:

$${\mu }_{max}={10}^7{\mu }_0$$

(10)

At the k-th iteration, $S^k$,$Y^k$ and ${\mu }^k$ are first fixed, and the background vector $v$ is updated. Define the intermediate variable,

$${\tilde{\mathcal{L}}}^k=D-S^k+{\displaystyle \frac{1}{{\mu }^k}}{Y}^k$$

(11)

Since the low-rank term is constrained to a column-shared structure $v{\mathbf{1}}^T$, the update of $v$ can be written as the following least-squares problem.

$$v^{k+1}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{v}{\min }{‖{\tilde{\mathcal{L}}}^k-v{\mathbf{1}}^T‖}_F^2$$

(12)

The closed-form solution ${\tilde{\mathcal{L}}}^k$ is the row-wise mean along the column direction, i.e.,

$$v_i^{k+1}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{j=1}^n}{\tilde{L}}_{ij}^k,i=\mathrm{1,2},\dots ,m$$

(13)

The low-rank term is then constructed as:

$$L^{k+1}=v^{k+1}{\mathbf{1}}^T$$

(14)

After obtaining $L^{k+1}$, fix $L^{k+1}$, $Y^k$, and ${\mu }^k$, and update the sparse term $S$. Define:

$${\tilde{S}}^k=D-L^{k+1}+{\displaystyle \frac{1}{{\mu }^k}}{Y}^k$$

(15)

The subproblem of $S$ can be written as:

$$S^{k+1}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{S}{\min }\lambda {‖\mathit{S}‖}_1+{\displaystyle \frac{{\mu }^k}{2}}{‖S-{\tilde{S}}^k‖}_F^2$$

(16)

Its closed-form solution is given by the soft-thresholding operator,

𝑆𝑘+1=𝑆𝜆/𝜇𝑘𝑆𝑘

(17)

where:

$$S_{\tau }\left(x\right)=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(x\right)\mathrm{m}\mathrm{a}\mathrm{x}\left(\left|x\right|-\tau ,0\right)$$

(18)

The reconstruction residual at the current iteration is computed as:

$$Z^{k+1}=D-L^{k+1}-S^{k+1}$$

(19)

and the Lagrange multiplier is updated as:

$$Y^{k+1}=Y^k+{\mu }^k{Z}^{k+1}$$

(20)

The penalty parameter is increased according to:

$${\mu }^{k+1}=\mathrm{m}\mathrm{i}\mathrm{n}\left(\rho {\mu }^k,{\mu }_{max}\right)$$

(21)

where $\rho =1.1$. Compared with a larger growth factor, this setting avoids excessively rapid growth of the penalty parameter, which may otherwise cause the iterations to be prematurely locked and thus hinder sufficient convergence of the background term.

To determine whether the algorithm has converged, the normalized residual is defined as:

$${Err}^{k+1}={\displaystyle \frac{{‖Z^{k+1}‖}_F}{{‖D‖}_F}}$$

(22)

To monitor the stability of the common background vector $v$, its relative change rate is further computed as:

$$∆v^{k+1}={\displaystyle \frac{{‖v^{k+1}-v^k‖}_2}{{‖v^k‖}_2+{10}^{-6}}}$$

(23)

Although low-rank decomposition can extract the common-structure vignette field from multiple images, the constraint in (5) does not impose smoothness. For vignetting, the true vignette field usually appears as a continuous, low-frequency, and smoothly varying spatial radiometric attenuation. The result obtained by low-rank decomposition alone may still contain local fluctuations, residual textures, or non-ideal structures. Therefore, polynomial fitting is further introduced after low-rank decomposition to obtain a vignette field that better matches practical imaging conditions.

As shown in Figure 3, the vignette field obtained after low-rank decomposition mainly characterizes the common structural component in multiple images. However, the result may still be affected by local residual textures and non-smooth disturbances, and is therefore insufficient to fully represent the vignette distribution of the actual imaging system. After further polynomial fitting, the obtained vignette field better conforms to the physical characteristics of the true vignette field.

2.2.4. Two-Dimensional High-Order Polynomial Fitting

We define the obtained vignette mask as$\mathcal{l}$, and subtract the global mean of the logarithmic observation matrix to remove the overall bias,

$${\mathcal{l}}_{\mathrm{a}\mathrm{d}\mathrm{j}}=\mathcal{l}-\overline{D_{\mathrm{l}\mathrm{o}\mathrm{g}}}$$

(24)

It is then reshaped into a two-dimensional image, defined as $\Gamma \in {\mathbb{R}}^{H\times W}$, which serves as the initial estimate of the two-dimensional vignette field.

Since the directly decomposed result $\Gamma$ may still contain local undulations and numerical fluctuations, a two-dimensional high-order polynomial fitting is further performed to obtain a smooth, continuous, and analytically representable vignette field.

Grid coordinates are constructed and normalized, with the coordinates set as: $x=\mathrm{1,2},\dots ,W$,$y=\mathrm{1,2},\dots ,H$

$$X^{\prime }={\displaystyle \frac{X-\overline{X}}{W/2}}, Y^{\prime }={\displaystyle \frac{Y-\overline{Y}}{H/2}}$$

(25)

where $\overline{X}$ and $\overline{Y}$ are the means of the horizontal and vertical coordinates, respectively.

The fitted surface is defined as:

$$P\left(X^{\prime },Y^{\prime }\right)=\sum _{p+q\le \mathrm{N}}{a}_{pq}{\left(X^{\prime }\right)}^p{\left(Y^{\prime }\right)}^q$$

(26)

where $a_{pq}$ denotes the polynomial coefficient to be solved, the indices satisfy $p,q\ge 0$, and $\mathrm{N}$ is the polynomial order.

The total number of coefficients in the two-dimensional $N$polynomial is:

$$N_c={\displaystyle \frac{\left(N+1\left)\right(N+2\right)}{2}}$$

(27)

The two-dimensional image $\Gamma$ is vectorized column-wise into a column vector.

$$\mathcal{z}=\mathrm{v}\mathrm{e}\mathrm{c}\left(\Gamma \right)\in {\mathbb{R}}^{m\times 1}$$

(28)

All basis functions ${\left(X^{\prime }\right)}^p{\left(Y^{\prime }\right)}^q$ are expanded term by term and concatenated to construct the design matrix.

$$A=\left[\begin{array}{ccccccc}1& X_1^{\prime }& Y_1^{\prime }& {\left(X_1^{\prime }\right)}^2& X_1^{\prime }{Y}_1^{\prime }& {\left(Y_1^{\prime }\right)}^2& \dots \\ 1& X_2^{\prime }& Y_2^{\prime }& {\left(X_2^{\prime }\right)}^2& X_2^{\prime }{Y}_2^{\prime }& {\left(Y_2^{\prime }\right)}^2& \dots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \\ 1& X_m^{\prime }& Y_m^{\prime }& {\left(X_m^{\prime }\right)}^2& X_m^{\prime }{Y}_m^{\prime }& {\left(Y_m^{\prime }\right)}^2& \dots \end{array}\right]$$

(29)

The coefficient vector is denoted by:

$$p={\left[a_{00},a_{10},a_{01},a_{20},a_{11},a_{02},\dots \right]}^{\mathrm{T}}$$

(30)

The fitting model can then be written as:

$$Ap\approx \mathcal{z}$$

(31)

The coefficients are solved by least squares:

$$p=\arg \underset{p}{\min }{‖Ap-\mathcal{z}‖}_2^2$$

(32)

After fitting, the smooth surface ${\Gamma }_{smooth}$ is reconstructed.

Finally, the minimum value is subtracted so that the vignette field is zero-referenced, yielding the final smooth vignette field.

$${\Gamma }_{\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{a}\mathrm{l}}={\Gamma }_{\mathrm{s}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}}-\mathrm{m}\mathrm{i}\mathrm{n}\left({\Gamma }_{\mathrm{s}\mathrm{m}\mathrm{o}\mathrm{o}\mathrm{t}\mathrm{h}}\right)$$

(33)

3. Results

3.1. Quantitative Evaluation

For quantitative analysis of the final extracted results, four metrics are used: mean absolute error (MAE), maximum absolute deviation (MAD), center-region mean absolute error (CenterMAE), and edge-region mean absolute error (EdgeMAE).

Let the corrected image be $\widehat{I}$, the vignette-free reference image be $I_{gt}$, and the image size be$H\times W$. Denote the total number of pixels by $N=H\times W$. The pixel-wise error is defined as:

$$e\left(i,j\right)=\widehat{I}\left(i,j\right)-I_{gt}\left(i,j\right)$$

(34)

Thus, MAE, MAD, CenterMAE, and EdgeMAE are expressed as follows.

$$MAE={\displaystyle \frac{\frac{1}{N}{\sum }_{i=1}^H{\sum }_{j=1}^W\left|e\left(i,j\right)\right|}{L}}$$

(35)

$$MAD={\displaystyle \frac{\underset{i,j}{\max }\left|e\left(i,j\right)\right|}{L}}$$

(36)

where $L$ is the upper limit of the dynamic range; it is 4095 for 12-bit images and 65535 for 16-bit images.

$$CerterMAE={\displaystyle \frac{\left(\frac{1}{\left|{\Omega }_c\right|}{\sum }_{\left(i,j\right)\in {\Omega }_c}\left|e\left(i,j\right)\right|\right)}{L}}$$

(37)

$$EdgeMAE={\displaystyle \frac{\left(\frac{1}{\left|{\Omega }_e\right|}{\sum }_{\left(i,j\right)\in {\Omega }_e}\left|e\left(i,j\right)\right|\right)}{L}}$$

(38)

where $r\left(i,j\right)=\sqrt{{\left({\displaystyle \frac{j-x_c}{W/2}}\right)}^2+{\left({\displaystyle \frac{i-y_c}{H/2}}\right)}^2}$denotes the normalized radial distance, the center region is ${\Omega }_c=\left\{\left(i,j\right) | r\left(i,j\right)\le 0.3\right\}$, and the edge region is${\Omega }_e=\left\{\left(i,j\right) | 0.3\le r\left(i,j\right)\le 1.0\right\}$.

For star images, the quantitative evaluation in this study is conducted from two aspects. First, the direct purpose of vignette correction is to reduce spatial nonuniformity in image-plane radiometric response. If the consistency of stellar responses across the image plane is improved after correction, then the spatial nonuniformity caused by vignetting can be considered effectively suppressed. The image-plane nonuniformity metric is therefore used for this purpose.

$$UR={\displaystyle \frac{1}{\overline{R}}}\sqrt{{\displaystyle \frac{1}{H\times W}}{\displaystyle \sum _{i=1}^H}{\displaystyle \sum _{j=1}^W}{\left[R\left(i,j\right)-\overline{R}\right]}^2}\times 100\%$$

(39)

Second, vignette correction is intended to restore response uniformity rather than alter the intrinsic radiometric information of the target itself. Therefore, after local background subtraction, the net stellar DN value should remain as stable as possible. If the correction process significantly changes the background-subtracted stellar net signal, it indicates that additional radiometric bias has been introduced while removing vignetting. Based on these considerations, the correction performance of different methods is quantitatively evaluated in terms of both nonuniformity reduction and preservation of net stellar DN values.

$$D_{net}=D_{star}-{\overline{D}}_{bg}$$

(40)

where $D_{net}$ denotes the net stellar DN value, $D_{star}$ denotes the DN response of the target region, and ${\overline{D}}_{bg}$ denotes the DN value of the local background region.

3.2. Measured Results and Analysis

As shown in Figure 4 and Figure 5, the vignette fields recovered by the three methods and their final correction results are further visualized. The correction performance of each method is analyzed and summarized based on the measured results.

For the four datasets, the results of the three methods in terms of MAE, MAD, CenterMAE, and EdgeMAE are shown in Figure 6. The proposed full method achieves the best overall performance, with MAE, MAD, CenterMAE, and EdgeMAE values of 0.482%, 3.646%, 0.138%, and 0.519%, respectively. Compared with the low-rank-only method, these four metrics are reduced by 22.9%, 32.9%, 71.7%, and 19.9%, respectively. Compared with the polynomial-only method, the reductions are 89.9%, 85.4%, 94.6%, and 89.9%, respectively. These results indicate that low-rank modeling can effectively extract the common vignette structure, while polynomial fitting further improves the smoothness and correction accuracy of vignette estimation.

As shown in Table 2, the image-plane nonuniformity of the four star samples are significantly reduced after vignette correction. Specifically, the nonuniformity metric decreases from 1.39%, 1.82%, 1.90%, and 1.92% before correction to 0.59%, 0.70%, 0.77%, and 0.80% after correction, respectively. This indicates that the response consistency at different field positions is substantially improved after correction. The result demonstrates that the proposed method can effectively suppress the spatial nonuniformity caused by vignetting and enhance the stability of stellar responses.

Meanwhile, the background-subtracted stellar DN values change only slightly before and after correction. This indicates that the correction process mainly restores image-plane response consistency without introducing significant bias into the net radiometric signal of the stars themselves. The visual results in Figure 7 further show that obvious spatial brightness nonuniformity is present before correction, whereas the background becomes more uniform after correction and the response consistency around the stellar targets is improved. Taken together, these results indicate that the proposed method reduces image-plane nonuniformity while preserving the stability of net stellar DN values, which is beneficial for subsequent star-based calibration analysis.

Figure 7. Stellar images before and after correction.

Figure 7. Stellar images before and after correction.

Table 2. Quantitative Assessment of Stellar Image Processing.

No.	Stellar net DN before correction	Image non-uniformity before correction (%)	Stellar net DN after correction	Image non-uniformity after correction (%)
1	2935	1.39	2932	0.59
2	2027	1.82	2027	0.70
3	2907	1.90	2881	0.77
4	2639	1.92	2625	0.80

Table 2. Quantitative Assessment of Stellar Image Processing.

No.	Stellar net DN before correction	Image non-uniformity before correction (%)	Stellar net DN after correction	Image non-uniformity after correction (%)
1	2935	1.39	2932	0.59
2	2027	1.82	2027	0.70
3	2907	1.90	2881	0.77
4	2639	1.92	2625	0.80

4. Discussion

The results show that the proposed full method achieves the best overall performance, indicating that low-rank modeling and polynomial fitting contribute in complementary ways to vignette correction. Specifically, the full method yields MAE, MAD, CenterMAE, and EdgeMAE values of 0.482%, 3.646%, 0.138%, and 0.519%, respectively. Compared with the low-rank-only method, these four metrics are reduced by 22.9%, 32.9%, 71.7%, and 19.9%, respectively; compared with the polynomial-only method, the reductions are 89.9%, 85.4%, 94.6%, and 89.9%, respectively. These results indicate that low-rank modeling can effectively extract the common vignette structure, whereas polynomial fitting further improves the smoothness of vignette estimation and the final correction accuracy.

From a mechanistic perspective, low-rank modeling effectively extracts the vignette component shared by multiple images acquired under similar conditions, whereas polynomial fitting further enforces the smooth low-frequency property expected for a physically plausible vignette field. Therefore, the former mainly captures shared structure, while the latter refines its spatial representation.

Compared with the low-rank-only method, the results of the full method suggest that low-rank decomposition alone is not sufficient to recover a physically realistic vignette field. Although low-rank modeling can isolate the common component across multiple images, the estimated field may still contain local fluctuations and non-ideal structures. These residual components are not fully consistent with the smooth spatial attenuation typically associated with vignetting. Polynomial fitting therefore plays an important refinement role by improving the continuity and physical plausibility of the recovered field.

The polynomial-only method, on the other hand, lacks an explicit shared-structure constraint and is therefore more sensitive to scene-dependent radiometric variations and local image content. As a result, smooth fitting alone may not reliably separate the vignette component from non-shared radiometric variations. This explains its significantly larger errors in the experiments. In contrast, the full method combines shared low-rank modeling with low-frequency polynomial fitting, thereby leveraging both shared-component extraction and smooth-field modeling to achieve higher correction accuracy.

The current framework relies on the assumption that the vignette field remains approximately stable within a short time period under similar imaging conditions. When this assumption is satisfied, multiple images can jointly constrain a common vignette field and improve estimation stability. However, if the vignette distribution changes significantly over time or across acquisition conditions, the shared low-rank model may no longer correspond to a single stable field, and the recovered result may become a compromise estimate among multiple vignette patterns. In addition, the present formulation is based on a relatively simple low-rank structure and a smooth polynomial model, which may be insufficient for more complex cases involving adaptive rank changes or locally varying high-frequency radiometric nonuniformities.

For star images, the correction effect is further analyzed from two perspectives: the suppression of image-plane nonuniformity and the preservation of net stellar signals. As shown in Table X and Figure 7, the image-plane nonuniformity of the four stellar samples is clearly reduced after correction, decreasing from 1.39, 1.82, 1.90, and 1.92 to 0.59, 0.70, 0.77, and 0.80, respectively. This indicates that response consistency at different field positions is significantly improved, demonstrating that the proposed method effectively suppresses the spatial radiometric nonuniformity caused by vignetting.

At the same time, the background-subtracted net stellar DN values remain largely stable, changing from 2935, 2027, 2907, and 2639 to 2932, 2027, 2881, and 2625, respectively. Compared with the pronounced improvement in the nonuniformity metric, the net DN values change only slightly, indicating that the correction mainly improves response consistency without introducing obvious bias into the intrinsic stellar radiometric signal. The visual comparison in Figure 7 also shows that the background distribution becomes more uniform after correction and the local radiometric response around the stellar targets becomes more stable. This further indicates that the proposed method effectively reduces image-plane nonuniformity while preserving net stellar radiometric information, which is beneficial for subsequent star-based radiometric calibration analysis.

Overall, the results show that combining shared low-rank modeling with polynomial fitting is an effective vignette-correction strategy. Low-rank modeling is used to identify the common vignette-related component, whereas polynomial fitting further improves its smoothness and stability, leading to better correction performance than either component alone and providing an effective solution for remote sensing image vignette correction.

5. Conclusions

This paper proposes a vignette correction method for remote sensing images based on low-rank modeling and polynomial fitting. By exploiting the shared vignette field among multiple images, the method constructs a data matrix in the logarithmic domain, extracts the common vignette component through low-rank decomposition, and further recovers a smooth vignette field using polynomial fitting. The final correction is performed in the original intensity domain. Overall, the proposed method provides a practical solution that balances common vignette-field extraction and smooth-field modeling.

The current method is mainly based on the assumption that the shared vignette field can be described by a low-rank model. Under more complex imaging conditions, however, the vignette field may be difficult to characterize using a fixed rank, especially a structure close to rank 1. Future work will therefore focus on adaptive rank selection or variable-rank modeling to enhance the expressive power of the model and improve robustness and generalization under more challenging scenarios.