Relaxation-Based Radiometric Normalization for Multitemporal Cross-Sensor Satellite Images

Multitemporal cross-sensor imagery is fundamental for the monitoring of the Earth’s surface over time. However, these data often lack visual consistency because of variations in the atmospheric and surface conditions, making it challenging to compare and analyze images. Various image-normalization methods have been proposed to address this issue, such as histogram matching and linear regression using iteratively reweighted multivariate alteration detection (IR-MAD). However, these methods have limitations in their ability to maintain important features and their requirement of reference images, which may not be available or may not adequately represent the target images. To overcome these limitations, a relaxation-based algorithm for satellite-image normalization is proposed. The algorithm iteratively adjusts the radiometric values of images by updating the normalization parameters (slope (α) and intercept (β)) until a desired level of consistency is reached. This method was tested on multitemporal cross-sensor-image datasets and showed significant improvements in radiometric consistency compared to other methods. The proposed relaxation algorithm outperformed IR-MAD and the original images in reducing radiometric inconsistencies, maintaining important features, and improving the accuracy (MAE = 2.3; RMSE = 2.8) and consistency of the surface-reflectance values (R2 = 87.56%; Euclidean distance = 2.11; spectral angle mapper = 12.60).


Introduction
Multitemporal cross-sensor imagery has become increasingly common in remotesensing applications, as the data obtained allow the analysis of changes and trends over time [1,2]. However, one of the major challenges in the usage of multitemporal cross-sensor imagery is the removal of radiometric inconsistency among images, as shown in Figure 1, which can lead to errors in analysis and the misinterpretation of results [3]. Factors such as different sensor characteristics, topography, atmospheric conditions, and sun-sensor geometry can all contribute to inconsistencies among images [4,5]. To address this issue, researchers have developed various methods for image normalization. Image normalization involves the adjustment of the radiometric values of an image to make it more comparable to the other images in a dataset [6,7]. Histogram matching is a commonly used method for image normalization, in which the histogram of an image is matched to a reference histogram [8]. Another method is dark-bright target normalization, which manually selects dark and bright targets in the reference image and normalizes the target image [4,9]. One of the most widely used methods for image normalization is linear regression using iteratively reweighted multivariate alteration detection (IR-MAD) [10]. The application of IR-MAD involves the detection and correction of the differences between images using statistical techniques; the radiometric values of each image are adjusted iteratively until the differences between the images are minimized [11,12]. Although these methods have shown promising results, they also have limitations. For example, histogram matching can Sensors 2023, 23, 5150 2 of 23 be sensitive to outliers [13], while dark-bright target normalization requires the presence of a suitable target in the image. While effective, IR-MAD requires a reference image and may not perform well in complex datasets with atmospheric or surface conditions, seasonal changes, and other variations. 0 2 of 23 images using statistical techniques; the radiometric values of each image are adjusted iteratively until the differences between the images are minimized [11,12]. Although these methods have shown promising results, they also have limitations. For example, histogram matching can be sensitive to outliers [13], while dark-bright target normalization requires the presence of a suitable target in the image. While effective, IR-MAD requires a reference image and may not perform well in complex datasets with atmospheric or surface conditions, seasonal changes, and other variations. To address the limitations of current image-normalization methods, this study proposes a novel relaxation algorithm that improves image consistency both quantitatively and qualitatively without requiring a reference image. The relaxation algorithm iteratively adjusts the radiometric values of images using a relaxation parameter that controls the degree of smoothing. By applying an iterative process to reduce radiometric differences and adjust radiometric values, the proposed method aims to minimize errors in multitemporal cross-sensor images. The proposed method was tested on multitemporal cross-sensor-image datasets with complex features, such as seasonal changes, water, topography, desert, snow, and cloud cover, and demonstrated significant improvements in image consistency compared to other methods.

Satellite Data
Several multitemporal cross-sensor images acquired from Landsat 8 (LAND-SAT/LC08/C02/T1_L2) and Sentinel 2 (COPERNICUS/S2_SR) surface-reflectance products were used to evaluate the proposed relaxation method. These datasets were chosen based on diversity in terms of geographic features and atmospheric conditions. The details of the datasets are summarized in Table 1. Note that all spectral bands of the images except the thermal, cirrus, and panchromatic bands were used.  To address the limitations of current image-normalization methods, this study proposes a novel relaxation algorithm that improves image consistency both quantitatively and qualitatively without requiring a reference image. The relaxation algorithm iteratively adjusts the radiometric values of images using a relaxation parameter that controls the degree of smoothing. By applying an iterative process to reduce radiometric differences and adjust radiometric values, the proposed method aims to minimize errors in multitemporal cross-sensor images. The proposed method was tested on multitemporal cross-sensorimage datasets with complex features, such as seasonal changes, water, topography, desert, snow, and cloud cover, and demonstrated significant improvements in image consistency compared to other methods.

Satellite Data
Several multitemporal cross-sensor images acquired from Landsat 8 (LANDSAT/LC08/ C02/T1_L2) and Sentinel 2 (COPERNICUS/S2_SR) surface-reflectance products were used to evaluate the proposed relaxation method. These datasets were chosen based on diversity in terms of geographic features and atmospheric conditions. The details of the datasets are summarized in Table 1. Note that all spectral bands of the images except the thermal, cirrus, and panchromatic bands were used.
Dataset #1 consists of surface-reflectance images from Landsat 8 and Sentinel 2 sensors over Tainan City, Taiwan, acquired from January to March 2020. The main objective of this dataset is to analyze the impact of normalization on urban areas and airport images, which could provide valuable information on land-cover and land-use patterns for use in further studies, such as urbanization studies [14,15]. Dataset #2 comprises surface-reflectance images of Mut, Egypt, one of the ancient cities located in the Western Desert. The city experiences a hot and arid climate, with little rainfall. The landscape is dominated by sand dunes and rocky mountains. This dataset focuses on a specific part of the city, and its objective is to analyze the impact of normalization on desert images, which poses significant challenges due to atmospheric aerosols and bright-reflective surfaces [16,17]. Dataset #3 Sensors 2023, 23, 5150 3 of 23 includes surface-reflectance images of Yakutsk, Russia, a city in the Sakha Republic known for its extremely cold subarctic winters. These images were acquired between February and March 2021, at the end of the winter season. The primary objective of this dataset is to investigate the impact of normalization on snowy images, which can be challenging to process due to their high reflectance, which often results in loss of detail and information in bright areas [18]. Datasets #4 and #5 are composed of surface-reflectance images of urban areas in Manaus, Brazil, and Dubbo, Australia, respectively. Both datasets contain cloudy images, but the amount of cloud pixels differs between the two. Dataset #4 has fewer cloud pixels than Dataset #5. The objective of these datasets is to study the impact of normalization on images with cloud pixels. Clouds can present challenging issues in image processing as they can obscure the underlying surface features and vary in thickness and height, resulting in different levels of reflectance [19,20]. Dataset #6 comprises surface-reflectance images from Legrena, a small coastal city located in southeastern Greece. The city's landscape is characterized by rocky cliffs, valleys, and beaches, making it a unique dataset with diverse geographical features. The dataset was selected for its water bodies, the topography effects caused by cliffs, and shadow effects on the acquired images. Processing of images of water bodies is challenging due to their significantly low surface-reflectance values compared with other inland surfaces [21]. Water is influenced by various factors, such as atmospheric conditions, water depth, and the sun angle, which can result in varying levels of surface reflectance and color [22]. Moreover, the topography of a region can lead to variations in surface reflectance due to changes in slope, aspect, and terrain ruggedness [23,24]. Additionally, shadow and shading can affect images, making accurate analysis difficult [25]. The objective of this dataset is to analyze the impact of normalization on these areas. Dataset #7 is the final dataset, and it consists of surface-reflectance images over Nashville City, Tennessee, USA. The images were acquired during the transition from summer to autumn season. The objective of this dataset is to identify the impact of normalization on seasonally affected pixels. Seasonally affected pixels are pixels in an image that change due to seasonal environmental changes, such as variations in temperature, humidity, and atmospheric conditions [26,27].

Overview of Relaxation Method
The proposed relaxation method consists of five steps, as illustrated in Figure 2. The first step involves selecting images from two satellite datasets, followed by image preprocessing in the second step. In the third step, pseudo-invariant feature (PIF) extraction is carried out using IR-MAD between paired images. Image normalization and optimization are performed using a regression model and satisfactory accuracy assessment in the next step. The third and fourth steps are repeated iteratively to improve the surface-reflectance is carried out using IR-MAD between paired images. Image normalization and optimization are performed using a regression model and satisfactory accuracy assessment in the next step. The third and fourth steps are repeated iteratively to improve the surface-reflectance consistency and enhance image visualization. Finally, the normalized images are aligned locally and globally to create a global normalized image cube. Each of these steps is described in detail below.

Image Selection and Preprocessing
Metadata filtering is an initial step in data processing, in which unsuitable images are filtered out. The parameters in this filtering process involve metadata of images, such as sensor type, acquisition date, and geographical location. Dates are filtered to obtain images acquired during specific seasons or events, while geographic filtering is useful in selection of images of a particular area of interest. Row-and path-based filtering are important for selecting images with similar viewing angles or creating mosaics. The objective of this filtering is to ensure that the selected images meet the specific analysis requirements and yield dependable results.
In this study, image pre-processing was used to prepare the satellite images for downstream analysis. The initial step in image pre-processing is image stacking, where Sentinel 2 images are combined with Landsat 8 images to create a seamless composite image [28,29]. Next, regions of interest (ROIs) are selected to ensure that the images have the same area and viewpoint of objects. This step is essential to eliminate any distortions caused by differences in sensor viewing angles and ensure that the images are correctly aligned [30]. Image resampling is then applied, whereby the resolution of Sentinel 2 images is reduced to a 30-m spatial resolution. This step is necessary to ensure compatibility with Landsat 8 images, which have a lower spatial resolution. Finally, the pre-processing stage is concluded with a pixel-aligned image cube, in which the images are aligned based on their pixels.

PIF Extraction
The next step in the image-processing pipeline is the extraction of pseudo-invariant features (PIFs). These PIFs refer to stable features on the Earth's surface that remain relatively unchanged over time, and they can be used to normalize remote-sensing data [31]. Examples of PIFs include urban areas, roads, and bare soil, which tend to have consistent surface-reflectance values over time. Although PIFs are not completely invariant and can be affected by seasonal or weather-related changes, they remain relatively stable and can be used to account for variations in atmospheric and environmental conditions during image normalization [32]. By incorporating PIFs into the normalization process, the accuracy and comparability of remote-sensing data can be improved, making it easier to analyze and interpret changes over time.
In this study, we performed iteratively reweighted multivariate alteration detection (IR-MAD) to select the PIFs from paired images. The IR-MAD approach is an efficient method to extract time-invariant features for radiometric normalization that builds upon the traditional multivariate alteration detection (MAD) algorithm [33,34]. As shown in Figure 2, the input dataset for PIF extraction is an image cube of multi-temporal and cross-sensor images denoted as S : {S 1 , S 2 , . . . , S n }, where n refers to the number of images. This image cube is the result of the image-pre-processing stage, in which the images from the two sensors are stacked, resampled and pixel aligned, so they are compatible with each other for further processing. The MAD is a method used to identify changes between two multispectral images, which utilizes traditional canonical correlation analysis (CCA) to model the linear combination of two multispectral images based on their order of correlation [35]. The differences between ordered pairs are referred to as MAD variates, which are represented by Equation (1). These MAD variates illustrate the variance of the two multispectral images (S 1 and S n ) based on their total k-spectral bands, ranked from the highest to the lowest. The eigenvectors a and b, along with their corresponding eigenvalues, are given in Equation (2). In this equation, C 11 and C nn represent the variance matrix of a single set variable (S 1 and S n , respectively), and the covariance between them is represented by C 1n and C n1 .  When the difference image follows a multivariate normal distribution, the sum of squared MAD variates (Equation (3)) can be shown to follow a chi-square distribution X 2 with f degrees of freedom equal to the number of spectral bands, expressed as Equation (4).
where σ MAD i represents the standard deviations of MAD variates and {MAD 1 , · · · , MAD k } is defined as The resulting distribution can be used to identify statistically significant anomalies in the difference image. This is achieved by calculating the probability that the sum of squared MAD variates will exceed a certain threshold within the chi-square distribution [33,36]. The pixels that satisfy the following probability conditions are selected as PIFs candidates.
where Pr(no change) is used to select the PIFs in kernel space. The selection of PIFs is formulated as where t represents the fixed threshold from percentiles in the X 2 distribution, usually larger than 0.9 to mask out the water pixels, cloud pixels, and shadow pixels in the image [33,37]. Next, an iterative reweighting scheme is performed to improve the detection of weak anomalies by using ω j as a weight factor for selected pixel-j. The reweighting scheme involves estimating weights based on the PIF selection (ω j ) in the previous iteration, which helps to identify weak anomalies that may have been missed in the initial weighting scheme (ω j = 1) [38]. This weight enters the calculation of mean, variance, and covariances (N is the number of pixels) of S 1 and S n , respectively.
for the mean value of S 1 and S n , and for the variance and covariance between S 1 and S n . In IR-MAD, iterations are performed until a convergence criterion is reached, which is usually determined by a maximum number of iterations or a minimum change in the MAD score between iterations. The MAD score measures the magnitude of change in pixel values between two images and is updated in each iteration. Once the MAD score reaches a stable value, the iteration is stopped, and PIFs are extracted.

Image Normalization and Optimization
After the PIFs are selected and extracted, the following step is image normalization. In this study, we applied the regression model (Equation (9)) to transform the radiometric condition of the target image (image S n ) into the radiometric condition of the reference image (image S 1 ) [39][40][41]. This regression utilizes PIFs from previous stage to remove the pixels that are constantly changing, such as clouds, water, and even vegetation in some cases.
where α n→1 and β n→1 are the slope and intercept of image S n to image S 1 , which is obtained from Equation (10), below. In this equation, σ(PIF S1 ) and σ(PIF Sn ) denote the standarddeviation values from PIF images of S 1 and S n . Next, PIF S1 and PIF Sn represent the mean values of PIF images of S 1 and S n respectively.
Furthermore, we propose an iterative optimization algorithm to obtain consistent surface-reflectance values of normalized images. The proposed relaxation algorithm minimizes a real value function or an error function by constructing a sequence of iterations [42,43]. The algorithm operates by starting with an initial value of normalization parameters (slope (α) and intercept (β) variables) from the previous normalization result and repeatedly updates them until the desired level of convergence is obtained. The proposed algorithm gradually aligns the radiometric values of the target image with those of any reference image, which results in a more consistent set of images.
In this study, two different networks were used to perform relaxation for image normalization, as illustrated in Figure 3. These two networks were used to align the radiometric conditions of all images with each other without relying on a reference image. Figure 3a shows a network in which each of the n-images in the dataset is connected to its two neighboring images (n − 1 and n + 1), and the first image is connected to the last image in the dataset. The network in this figure is called a ring network. In contrast, Figure 3b displays more links between its images than Figure 3a, with each image in this network connected to all the other images. This network is called a fully connected network. The inputs to these networks are multitemporal cross-sensor images, and the outputs are the normalized images. Hence, two different error functions were used in this study to match the conditions in Figure 3, expressed in Equations (11) and (12).
In the previous equations, A = α 1→2 ; α 1→3 ; . . . ; α (n−1)→n and B = {β 1→2 ; β 1→3 ; . . . ; β (n−1)→n } were the set of slope and intercept components of the normalization. These sets (A, B) are updated through a relaxation-iteration process until the minimum error value is obtained. Next, O : {O 1 , O 2 , . . . , O n } represents the PIFs masked images formulated in Equations (13) and (14). These are masked images that are generated using common PIFs, which are selected PIFs that are shared among multiple images and can be obtained by combining multiple PIF-selection results. (14) In this study, two different networks were used to perform relaxation malization, as illustrated in Figure 3. These two networks were used to a metric conditions of all images with each other without relying on a referen ure 3a shows a network in which each of the n-images in the dataset is co two neighboring images (n − 1 and n + 1), and the first image is connected to in the dataset. The network in this figure is called a ring network. In cont displays more links between its images than Figure 3a, with each image i connected to all the other images. This network is called a fully connected inputs to these networks are multitemporal cross-sensor images, and the o normalized images. Hence, two different error functions were used in this s the conditions in Figure 3, expressed in Equations (11) and (12).  The (11) and (12) are error functions in this study, with the objective of obtaining the minimum total error value of each image against other images in the dataset. Specifically, Equation (11) is the error function for the relaxation algorithm, which is carried out by the ring network ( Figure 3a). This equation calculates the sum error from each PIF-masked image with respect to their two neighboring images. Equation (12) is the error function for the relaxation algorithm, which is performed by using a fully connected network (Figure 3b). This equation calculates the total error from each PIF-masked image against all its neighboring images.
The name of relaxation is derived from the fact that it gradually "relaxes" the constraints on the normalization components, allowing them to gradually approach a solution that satisfies all the constraints simultaneously. This process is repeated until the error function is minimized to an acceptable level, at which point the final set of values can be considered consistent.

Image-Cube Creations and Alignments
The creation of an image cube involves stacking multiple images acquired at different times to form a three-dimensional data structure. The image cube in this study is the result of the normalization process and represents a time series of images with consistent radiometric values [44,45]. It is created by selecting the best result of normalization from several iterations and conducting a satisfactory accuracy assessment.
After the normalized images are selected, they are aligned to create a local cube with the times of acquisition in ascending order, along with their specific locations. This alignment is critical to ensure that the images are in the correct spatial orientation relative to each other, allowing accurate comparison and analysis over time [46]. After the local cube is created, a coordinate-reference-system transformation is performed to make a global image cube that contains local image cubes. The global-coordinate-reference system used for the image cube is typically WGS 84, which is a commonly used global reference system in remote sensing. This transformation is necessary to ensure that images are accurately georeferenced, enabling accurate spatial analysis and interpretation [47]. The resulting image cube can be used to detect and analyze changes over time, such as land-use and land-cover changes, vegetation dynamics, and urban growth.

Qualitative Assessment
During the visual assessment stage, we compared five out of seven datasets separately, which consisted of images with seasonal changes (United States of America), as well as water and topography (Greece), desert (Egypt), snow (Russia), and cloud (Brazil) features. The remaining datasets (Taiwan and Australia) are shown in Appendix A. This step allowed us to evaluate the effectiveness of our normalization algorithm in preserving visual quality while reducing spectral variability.

Seasonal Features
The image dataset for the United States of America image ( Figure 4) presented a challenging task due to the seasonal changes that affected the urban areas. However, all the normalization methods, including our proposed relaxation algorithm, performed well in reducing the image inconsistencies. However, our approach stood out as it improved the image consistency and maintained the seasonally affected pixels effectively. In addition, the seasonal transitions were smoother in our proposed relaxation algorithm than in the IR-MAD using a ring network. These results highlight the effectiveness of our approach in handling complex image datasets. The image dataset for the United States of America image ( Figure 4) presented a challenging task due to the seasonal changes that affected the urban areas. However, all the normalization methods, including our proposed relaxation algorithm, performed well in reducing the image inconsistencies. However, our approach stood out as it improved the image consistency and maintained the seasonally affected pixels effectively. In addition, the seasonal transitions were smoother in our proposed relaxation algorithm than in the IR-MAD using a ring network. These results highlight the effectiveness of our approach in handling complex image datasets.  Figure 5 presents the Greece image dataset, which focuses on areas with water and topographical features. Overall, normalization was performed effectively on this dataset. However, image inconsistency still occurred. The presence of water pixels could be one of the factors contributing to this inconsistency. Water has a high reflectance in blue and green spectral bands, but low reflectance in red and near-infrared bands, resulting in color variations that can vary significantly, depending on the spectral bands used in the image [48]. Additionally, water surfaces can be affected by changes in lighting conditions and atmospheric effects, which can further complicate the normalization process. Moreover, the effect of topography on normalization can also pose challenges in maintaining consistent image features across different terrain types. In areas with significant topographies, there can be variations in lighting conditions, shading, and perspective, which can affect  Figure 5 presents the Greece image dataset, which focuses on areas with water and topographical features. Overall, normalization was performed effectively on this dataset. However, image inconsistency still occurred. The presence of water pixels could be one of the factors contributing to this inconsistency. Water has a high reflectance in blue and green spectral bands, but low reflectance in red and near-infrared bands, resulting in color variations that can vary significantly, depending on the spectral bands used in the image [48]. Additionally, water surfaces can be affected by changes in lighting conditions and atmospheric effects, which can further complicate the normalization process. Moreover, the effect of topography on normalization can also pose challenges in maintaining consistent image features across different terrain types. In areas with significant topographies, there can be variations in lighting conditions, shading, and perspective, which can affect the appearance of images [49,50]. This can result in losses of contrast and detail in areas with high elevation, or the overemphasis of features in areas with low elevation. Our proposed relaxation algorithm and IR-MAD using a fully connected network may present difficulties in maintaining topographical features during normalization, which can lead to decreased contrast in images.  Figure 6 presents the Egypt image datasets. From this dataset, we found that our proposed relaxation algorithm with IR-MAD using the fully connected network showed a similar level of improvement in image consistency. However, the improvement using IR-MAD with the ring network was only partial and only applied to some of the images. As shown in the figure, we observed the presence of bright pixels in the middles of the images. These pixels may have been caused by different factors, including sensor noise, saturation, and atmospheric effects [51]. In addition, aerosol effects, such as dust and sand particles, which are common in desert images, may also have contributed to the appearance of these bright pixels [52]. Another possibility is that the normalization process introduced artifacts, which led to the bright pixels.  Figure 6 presents the Egypt image datasets. From this dataset, we found that our proposed relaxation algorithm with IR-MAD using the fully connected network showed a similar level of improvement in image consistency. However, the improvement using IR-MAD with the ring network was only partial and only applied to some of the images. As shown in the figure, we observed the presence of bright pixels in the middles of the images. These pixels may have been caused by different factors, including sensor noise, saturation, and atmospheric effects [51]. In addition, aerosol effects, such as dust and sand particles, which are common in desert images, may also have contributed to the appearance of these bright pixels [52]. Another possibility is that the normalization process introduced artifacts, which led to the bright pixels. Figure 7 shows the Russia image dataset, which consisted of snowy images. However, it was difficult to assess the effectiveness of the normalization methods, including IR-MAD and our proposed relaxation algorithm, due to the dominance of snow pixels in the images. Therefore, the normalization methods applied to these images may not have exerted a significant impact on their appearance. This made it difficult to visually assess any improvements in the consistency of the images after the normalization. Nevertheless, we observed a significant change in the second image of the dataset, which became more similar to its neighboring images after the normalization process.

Snow Features
a similar level of improvement in image consistency. However, the improvement using IR-MAD with the ring network was only partial and only applied to some of the images. As shown in the figure, we observed the presence of bright pixels in the middles of the images. These pixels may have been caused by different factors, including sensor noise, saturation, and atmospheric effects [51]. In addition, aerosol effects, such as dust and sand particles, which are common in desert images, may also have contributed to the appearance of these bright pixels [52]. Another possibility is that the normalization process introduced artifacts, which led to the bright pixels.   Figure 7 shows the Russia image dataset, which consisted of snowy images. However, it was difficult to assess the effectiveness of the normalization methods, including IR-MAD and our proposed relaxation algorithm, due to the dominance of snow pixels in the images. Therefore, the normalization methods applied to these images may not have exerted a significant impact on their appearance. This made it difficult to visually assess any improvements in the consistency of the images after the normalization. Nevertheless, we observed a significant change in the second image of the dataset, which became more similar to its neighboring images after the normalization process.  Figure 8 presents the Brazil image dataset, which focused on urban areas but also included images with cloud presence. We found that despite the presence of clouds, image normalization was still effectively conducted using both IR-MAD and our proposed re-  Figure 8 presents the Brazil image dataset, which focused on urban areas but also included images with cloud presence. We found that despite the presence of clouds, image normalization was still effectively conducted using both IR-MAD and our proposed relaxation algorithm. While the appearance of objects behind the clouds improved after the normalization, the clouds themselves remained relatively unchanged in appearance. However, we did not observe any noticeable artifacts or distortions in the cloud features resulting from the normalization process. Overall, the normalization methods, especially our proposed relaxation algorithm, showed promising results in reducing image inconsistencies and improving visual consistency.

Quantitative Assessment
The quantitative assessment of the image-normalization process involves several methods, including spectral comparison, temporal comparison, and statistical comparison. Spectral comparison involves comparing the reflectance values of different wavelengths to identify differences in spectral response. Temporal comparison involves comparing the surface-reflectance values of the same area over time to identify changes or trends. Statistical comparison involves calculating various statistical parameters, such as the mean and correlation coefficient, to assess the accuracy and consistency of normalization results. Together, these methods provide a comprehensive assessment of the effectiveness of the image-normalization process.

Spectral Comparisons
Spectral comparison involves comparing the histograms of surface-reflectance values of the same image before and after normalization. This comparison provides insights into the distribution of surface-reflectance values across different spectral bands and how they change due to normalization [53]. Significant changes in the shape or distribution of the histogram may indicate problems in the normalization process or highlight differences in the surface-reflectance values of an image before and after normalization. Figure 9 dis-

Quantitative Assessment
The quantitative assessment of the image-normalization process involves several methods, including spectral comparison, temporal comparison, and statistical comparison. Spectral comparison involves comparing the reflectance values of different wavelengths to identify differences in spectral response. Temporal comparison involves comparing the surface-reflectance values of the same area over time to identify changes or trends. Statistical comparison involves calculating various statistical parameters, such as the mean and correlation coefficient, to assess the accuracy and consistency of normalization results. Together, these methods provide a comprehensive assessment of the effectiveness of the image-normalization process.

Spectral Comparisons
Spectral comparison involves comparing the histograms of surface-reflectance values of the same image before and after normalization. This comparison provides insights into the distribution of surface-reflectance values across different spectral bands and how they change due to normalization [53]. Significant changes in the shape or distribution of the histogram may indicate problems in the normalization process or highlight differences in the surface-reflectance values of an image before and after normalization. Figure 9 displays histograms for images from each dataset that compare the original images (before normalization) and the IR-MAD and relaxation-algorithm images (after normalization). The histograms of the IR-MAD and proposed relaxation algorithm showed similar distributions to the original images, indicating that both normalization methods preserved the original distribution of the surface-reflectance values in the images. This suggests that both methods can successfully normalize images by maintaining the features and details in images during the normalization process [54].
ors 2023, 23, 5150 13 of Figure 9. Histogram comparisons of the images before normalization and after normalization.

Temporal Comparisons
Temporal comparison involves analyzing changes in surface-reflectance values ov time. Figure 10 shows the temporal trend of surface reflectance between the original i ages, those obtained through IR-MAD, and those from the relaxation algorithm. The s face-reflectance value in the figure is the average value of the surface-reflectance of a giv area over a time period. This value was calculated by adding up all the individual refl tance values of each pixel in the area and dividing the total by the number of pixels. Bas on Figure 10, it is evident that the design of the network plays a crucial role in improvi the consistency of surface-reflectance values between multitemporal cross-sensor imag Our proposed relaxation algorithm and the IR-MAD using a fully connected netwo aligned well with the original images in that they followed the same trend, with relativ consistent surface-reflectance values over time, as depicted by the black and orange lin in Figure 10. On other hand, the proposed relaxation algorithm and the IR-MAD usin

Temporal Comparisons
Temporal comparison involves analyzing changes in surface-reflectance values over time. Figure 10 shows the temporal trend of surface reflectance between the original images, those obtained through IR-MAD, and those from the relaxation algorithm. The surfacereflectance value in the figure is the average value of the surface-reflectance of a given area over a time period. This value was calculated by adding up all the individual reflectance values of each pixel in the area and dividing the total by the number of pixels. Based on Figure 10, it is evident that the design of the network plays a crucial role in improving the consistency of surface-reflectance values between multitemporal cross-sensor images. Our proposed relaxation algorithm and the IR-MAD using a fully connected network aligned well with the original images in that they followed the same trend, with relatively consistent surface-reflectance values over time, as depicted by the black and orange lines in Figure 10.
On other hand, the proposed relaxation algorithm and the IR-MAD using a ring network demonstrated a similar pattern, with the IR-MAD presenting a surface-reflectance value that contradicted the original, and our method attempted to reduce the difference but remained opposite to the original. The large disparity between the normalized images may have been due to various factors, such as the complexity of the scene, or the atmospheric conditions.

Loss-Value Measurements
Loss-value measurements assess the effectiveness of a normalization algorithm by evaluating the change in the total error value throughout the iteration process. This measurement helps to determine whether an algorithm achieves a minimum error value. Ideally, the objective value should decrease with each iteration, indicating that the algorithm is making progress towards minimizing the differences between normalized images [43,55].
The fluctuations in the loss-values in this study are presented in Figure 11. The initial objective value in iteration-0 is the original loss-value obtained by Equation (11) for the ring network and Equation (12) for the fully connected network. Iteration-1 represents the IR-MAD's loss value, while the proposed relaxation's loss-value is selected from the minimum value of iterations 2 to 10. Based on Figure 11, the overall normalization using the ring network and the fully connected network presented the same trends through the iterations. This figure shows that our proposed relaxation method using a fully connected network needs between one and five additional iterations to achieve a minimum lossvalue, while the proposed relaxation algorithm using the ring network needs nine additional iterations. In addition, the normalizations using the ring network had large fluctuations compared with the normalizations using the fully connected networks, which indicates that the normalization using the ring network may have been less efficient.

Loss-Value Measurements
Loss-value measurements assess the effectiveness of a normalization algorithm by evaluating the change in the total error value throughout the iteration process. This measurement helps to determine whether an algorithm achieves a minimum error value. Ideally, the objective value should decrease with each iteration, indicating that the algorithm is making progress towards minimizing the differences between normalized images [43,55].
The fluctuations in the loss-values in this study are presented in Figure 11. The initial objective value in iteration-0 is the original loss-value obtained by Equation (11) for the ring network and Equation (12) for the fully connected network. Iteration-1 represents the IR-MAD's loss value, while the proposed relaxation's loss-value is selected from the minimum value of iterations 2 to 10. Based on Figure 11, the overall normalization using the ring network and the fully connected network presented the same trends through the iterations. This figure shows that our proposed relaxation method using a fully connected network needs between one and five additional iterations to achieve a minimum loss-value, while the proposed relaxation algorithm using the ring network needs nine additional iterations. In addition, the normalizations using the ring network had large fluctuations compared with the normalizations using the fully connected networks, which indicates that the normalization using the ring network may have been less efficient. The specific loss-values in this study are listed in Table 2. According to the table, our proposed relaxations outperformed the IR-MAD algorithm in terms of reducing the error value. In particular, the relaxation algorithm achieved a much lower loss value in all locations compared to the IR-MAD algorithm. This suggests that the relaxation algorithm is more effective in minimizing the differences between normalized images and, thus, in improving the consistency of surface-reflectance values between multitemporal cross-sensor images. Table 2. Comparisons of loss-values of original images and normalized images, which represent the total error between the images in each dataset. Bold font indicates the lower value of loss-values.

Accuracy Measurements
The accuracy measurement in this study involved evaluating the error between the normalized images using metrics such as the mean absolute error (MAE) and the root mean square error (RMSE). These metrics provide an indication of how close the normalized images are to each other and identify any differences or discrepancies between them [56,57]. Smaller MAE and RMSE values indicate better alignment and similarity between normalized images. The formula of the MAE and RMSE between two normalized images are written in Equations (15) and (16). The specific loss-values in this study are listed in Table 2. According to the table, our proposed relaxations outperformed the IR-MAD algorithm in terms of reducing the error value. In particular, the relaxation algorithm achieved a much lower loss value in all locations compared to the IR-MAD algorithm. This suggests that the relaxation algorithm is more effective in minimizing the differences between normalized images and, thus, in improving the consistency of surface-reflectance values between multitemporal crosssensor images. Table 2. Comparisons of loss-values of original images and normalized images, which represent the total error between the images in each dataset. Bold font indicates the lower value of loss-values.

Accuracy Measurements
The accuracy measurement in this study involved evaluating the error between the normalized images using metrics such as the mean absolute error (MAE) and the root mean square error (RMSE). These metrics provide an indication of how close the normalized images are to each other and identify any differences or discrepancies between them [56,57]. Smaller MAE and RMSE values indicate better alignment and similarity between normalized images. The formula of the MAE and RMSE between two normalized images are written in Equations (15) and (16).
where P i is the surface-reflectance value in pixel i of the normalized image, M i is the corresponding surface-reflectance value in pixel i of the other normalized image, and N is the total number of pixels in the image. Table 3 presents the accuracy measurements from the original images (before normalization), IR-MAD, and proposed relaxation (both after normalization). Based on the table, the MAE and RMSE values for relaxation were consistently lower than those for the IR-MAD and the original images. This suggests that our proposed relaxation algorithm is a more effective normalization method for reducing image discrepancies and errors.

Correlation Measurements
Correlation measurements are used to evaluate the similarity between images before normalization and after normalization. The correlation coefficient ranges from −1 to 1, with 1 indicating a perfect positive correlation, 0 indicating no correlation, and −1 indicating a perfect negative correlation. A correlation coefficient close to 1 suggests that the normalized images are highly similar to the original images, while a correlation coefficient close to 0 indicates a low similarity between the images [58]. In this stage, the correlation coefficient is calculated for each spectral band for the entire image, formulated in Equation (17).
where k represents the number of bands, P j and M j are the surface-reflectance values of the j-th band in two images, P and M, P j and M j are the means value of two images, and σP and σM are the standard deviations of P and M, respectively. Table 4 shows the correlation measurements for the original images, the IR-MAD, and the relaxation algorithm for the image normalization. According to the table, the IR-MAD achieved the highest correlation coefficient in some of the experiments. However, the overall assessment shows that our proposed relaxation method achieved the highest correlation coefficient among the three methods (original, IR-MAD, and relaxation). This indicates that our relaxation algorithm is able to effectively maintain image similarity across different spectral bands before and after normalization, leading to a more accurate and consistent image representation.

Spectral Distance Measurements
Spectral distance measurements are used to evaluate similarity between normalized images [4]. Generally, there are two methods for distance measurements: Euclidean distance (ED) and spectral angle mapper (SAM). The Euclidean distance is a measure of the distance between two points in a multidimensional space [59]. In the context of image normalization, the points represent the surface-reflectance values of each pixel in the images. A shorter Euclidean distance indicates a greater degree of similarity between two pixels. The Euclidean distance between two pixels is calculated using the following formula: where P i and M i are the surface reflectance of pixel i in the two normalized images, P and M. On other hand, the spectral angle mapper is a measure of the angle between the spectral vectors of two pixels [60]. A smaller spectral angle indicates a greater degree of similarity between two pixels. The spectral angle between two pixels is calculated using the following formula: Table 5 shows the spectral distance measurements of the Euclidean distance and spectral angle mapper for the original images, the IR-MAD images, and the relaxation-algorithm images. The purpose of comparing the spectral distance before and after normalization is to determine improvements in image similarity. In some of the experiments, the IR-MAD achieved better spectral distance measurements. However, the proposed relaxation algorithm outperformed the IR-MAD and the original images, with consistently lower spectral distance values. This suggests that our relaxation algorithm is more effective in reducing the differences between images, resulting in a more similar representation of images after normalization. Based on the visual assessment, our proposed relaxation algorithm was found to outperform the IR-MAD and the original images in terms of improving the image consistency. Despite the complexity of some of the datasets, including the presence of clouds, water, and seasonal changes, our approach effectively reduced the image inconsistencies while maintaining the important features. However, it was observed that our relaxation algorithm was not able to preserve topographic features as well as the IR-MAD using the ring network. Therefore, the statistical assessment results suggest that our proposed relaxation produced better assessments than the IR-MAD and the original images. There could be several possible reasons behind the superior performance of our proposed relaxation algorithm r.
Firstly, our proposed relaxation algorithm incorporates a looping structure that repeats the normalization process until a desired level of image similarity is achieved. This allows the IR-MAD to continually refine the normalization parameters based on the difference between the current normalized image and its paired image. As a result, the algorithm can adjust the normalization parameters to better match the paired image, leading to a more accurate and consistent normalization. This iterative approach enables the relaxation algorithm to effectively reduce discrepancies between images, resulting in better assessments compared to the original images and the IR-MAD. Moreover, our proposed relaxation can ensure that normalization is applied consistently across all images. This is because the normalization parameters may vary slightly between images. Therefore, the iteration process helps to ensure the normalization is applied consistently to all images, which can improve the overall image similarity.
Based on the experimental results, our proposed relaxation algorithm demonstrated high effectiveness in improving the image consistency for various satellite-image datasets. However, one of its limitations is that it may not preserve certain features in images. For example, in some cases, topographic features may be lost during the normalization process, which can result in a decrease in contrast in the normalized images. Additionally, while our approach was able to maintain seasonal affected pixels, there may still have been some inconsistencies due to the complex features of the dataset, such as clouds and seasonal changes. Furthermore, we found that the improvements in image consistency for snowy images were insignificant. This may have been due to the nature of snow-covered surfaces, which are already uniform in color and texture, making it difficult for normalization methods to exert a significant impact. Therefore, it is important to carefully assess the results of the normalization process and ensure that the key features of the images are preserved while considering the limitations of the relaxation process.

Conclusions
Based on the findings of this study, the proposed relaxation algorithm demonstrated its effectiveness in enhancing the image consistency in multitemporal cross-sensor-image datasets. The combination of qualitative and quantitative assessments revealed that our algorithm surpassed both the IR-MAD method and the original images in reducing the image inconsistencies, preserving the crucial features, and enhancing the accuracy and consistency of the surface-reflectance values. These results indicate that the relaxation algorithm is a robust solution for addressing the challenges posed by complex datasets, including those with seasonal changes, water, urban areas, and cloud cover.
In addition to the algorithm's performance, the network design employed in this study deserves attention. Two different network configurations, namely a ring network and a fully connected network, were utilized for the relaxation process in image normalization. The ring network connected each image in the dataset to its neighboring images, while the fully connected network established links between all the images. The use of these network designs allowed the comprehensive alignment of the radiometric conditions across all the images without relying on a reference image.
The comparison between the two network configurations demonstrated the significant influence of the network design on the algorithm's effectiveness. The normalization process using the fully connected network was faster than that using the ring network. In addition, the fully connected network produced superior normalization results to those obtained with the ring network. This suggests that the fully connected network, as the network design of choice, not only ensures faster processing, but also ensures better normalization outcomes.
However, it is important to acknowledge that our relaxation algorithm exhibits limitations in preserving topographic features compared to the IR-MAD method that utilizes the ring network. To overcome this drawback, future studies could explore potential strategies to combine the strengths of these different methods, aiming to achieve further improvements in normalization outcomes for multitemporal cross-sensor images. By leveraging the advantages of both approaches, researchers can potentially enhance the preservation of topographic features while maintaining the overall effectiveness of the relaxation algorithm.
Overall, this study emphasizes the significance of employing effective image-normalization methods to enhance the accuracy and consistency of multitemporal cross-sensor-image datasets. The proposed relaxation algorithm stands as a valuable tool for researchers and practitioners working with these types of datasets, offering enhanced capabilities for reliable analysis. Moreover, further research could investigate the use of the relaxation algorithm in different remote-sensing applications to assess its potential in diverse contexts, expanding its utility beyond the scope of this study. By exploring its adaptability and performance in various scenarios, we can gain deeper insights into the versatility and broader applications of the relaxation algorithm in the field of remote sensing.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The following figures (Figures A1 and A2) present the results of the normalization of the urban-features images in Taiwan and the cloud and seasonal features in Australia, respectively, using the proposed relaxation algorithm and the IR-MAD algorithm with fully connected and ring networks. A visual comparison of the figures reveals that the proposed relaxation algorithm produced normalized images that were visually comparable to those obtained using the IR-MAD algorithm with the fully connected networks. However, the normalized images from the IR-MAD algorithm with the ring network show minor changes compared to the original images, indicating that the inconsistency was not effectively reduced. Therefore, the proposed relaxation algorithm outperformed the IR-MAD algorithm in terms of reducing the image inconsistency. and performance in various scenarios, we can gain deeper insights into the versatility and broader applications of the relaxation algorithm in the field of remote sensing.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The following figures (Figures A1 and A2) present the results of the normalization of the urban-features images in Taiwan and the cloud and seasonal features in Australia, respectively, using the proposed relaxation algorithm and the IR-MAD algorithm with fully connected and ring networks. A visual comparison of the figures reveals that the proposed relaxation algorithm produced normalized images that were visually comparable to those obtained using the IR-MAD algorithm with the fully connected networks. However, the normalized images from the IR-MAD algorithm with the ring network show minor changes compared to the original images, indicating that the inconsistency was not effectively reduced. Therefore, the proposed relaxation algorithm outperformed the IR-MAD algorithm in terms of reducing the image inconsistency.