Fusión de imágenes satelitales usando la transformada de Brovey y enriquecimiento espectral sobre computación heterogénea CPU/GPU Satellite-image fusion using Brovey transform and spectral richness calibration on heterogeneous computing CPU/GPUs

Objetivo: Acelerar el proceso de fusión de imágenes satelitales mediante un modelo de procesamiento paralelo masivo basado en computación heterogénea (CPU/GPU). Metodología: Para cumplir este propósito, se realizó una implementación paralela-heterogénea de la transformada de Brovey. La implementación se evaluó con un conjunto de imágenes de 4 tamaños diferentes y se comparó frente a una implementación secuencial en CPU. Resultados: Se evidenció un Speed-up de hasta 532X en el proceso de fusión de una imagen de 8192 píxeles. En cuanto a calidad de la imagen obtenida, al obtener el coeficiente de correlación entre la imagen fusionada y la pancromática, se obtuvo un promedio de detalle espacial por banda del 0.9714 en un espacio de color (R,G,B). La implementación se encuentra disponible en https://github.com/Parall-UD/ParallelBrovey-Transform. Conclusiones: Se demostró que la programación paralela masiva en plataformas heterogéneas CPU/GPU representa una solución efectiva a la alta exigencia computacional del proceso de fusión de imágenes satelitales sin afectar la calidad de la imagen obtenida, comparada con soluciones secuenciales.

Revista Investigación e Innovación en Ingenierías, vol. 9, 2°. 1, pp. 7-25, Enero-julio 2021. DOI: https://doi.org/10.17081/invinno.9.2.3961 Introducción Image processing is aimed at improving image quality so that relevant information can be detected, boosting pattern and feature recognition [1]. When taking an image as the starting point for data processing, it is possible to improve image quality by applying several processes, e.g. pixel-by-pixel operations, segmentation, brightness adjustments, contrast calibration, equalization, filtering, and image fusion. In [2], a literature survey gathered a set of techniques addressing image processing at a pixel level for teledetection. These techniques were based on Principal Component Analysis (PCA), Intensity-Hue-Saturation (IHS), the Brovey transform, and the Wavelet transform among others. The studies conclude that image processing methods have a significant impact on the quality of the resulting (processed) images.
Image fusion combines relevant information from a given set of images to produce a single image that exhibits higher quality when compared to any of the original images. For satellite image fusion it is necessary to provide a Panchromatic image, which has high spatial resolution, and a Multispectral image, which is has considerable spectral information [3]. In [4], a study was conducted employing seven image-fusion techniques. The effectiveness of the techniques was assessed using both quantitative and qualitative metrics. Based on a series of statistical metrics, the Wavelet transformed was found to be the most effective technique, followed by the Brovey transform. Moreover, in [5], the Fast Haar Wavelet Transform (FHWT) was identified as an excellent technique for satellite image fusion, achieving the best performance in several metrics such as the correlation coefficient, the Relative average spectral error, the ERGAS metric, and the Universal Quality index Qu.
The capacity of parallel computing using Graphic Processing Units (GPUs) has raised an increasing interest among researchers and engineers studying high-performance computing [6]. For example, in [7], a search for electromagnetic emission sources was proposed based on interferometric data obtained from spectralline radio-astronomy using parallel programing techniques; a 3.2 times faster processing was achieved when compared to the CPU process alone. Other studies [8] proposed an image-fusion algorithm employing convolution-based neural networks by means of a parallel processing architecture; results were compared to those obtained with traditional methods, showing an improvement in the network's training time and satisfactory quality of the fusion process. When conducting satellite image fusion, the goal is to cover large areas of land. To this end, large-scale images are used, which significantly increases the processing time when employing a serial architecture. For example, in [9], assumptions are made about accelerating processing techniques, such as super-resolution (SR), by using GPU parallel processing. This paper presents preliminary results on satellite image fusion employing the Brovey transform along with spectral richness calibration using a CPU/GPU heterogeneous computing architecture. To guarantee correlation and quality, both spectral and spatial, several statistical metrics are assessed, such as the correlation coefficient, the Universal Image Q index, as well as the ERGAS and Bias tests [10]. The purpose of employing this method using a parallel processing architecture is to reduce the time necessary to complete the satellite image fusion stage together with the calibration of spectral richness. The type of image processing introduced herein is conducted using Python over the CUDA parallel processing architecture. This paper is organized as follows. Section 2 presents some background and related work. Section 3 introduces the proposed methodology by parts; namely the general approach, the method, and the assessment procedure. Moreover, results are analyzed with regard to running times, visual comparison, and statistical metrics. Finally, conclusions and future works are presented.

Materials and methods
Image processing refers to a variety of techniques that maximize the usefulness of the information contained in an image for a particular target context [11]. In [12], image processing techniques were applied to assess the performance of a convolution-based neural network when classifying images in the context of augmented reality applications. The results were significant in terms of feature extraction and pattern recognition improvements, namely the image processing techniques increased classification accuracy by 17%, which shows the importance of applying image processing techniques.
Image fusion is commonly defined as the process of combining information, extracted from two or more images of the same scene, to consolidate a single compound image of higher quality, which provides improved visual experience and better informatic processing [13]. This improvement occurs at different levels, namely a pixel level, a functional level and a strategic level. Satellite image fusion is commonly performed between a panchromatic image (PAN) and a multispectral image (MUL), which can be seen as a structural, functional system. The structural content is provided by the PAN, since it brings high spatial resolution, whereas the functional content is determined by the MUL due to its spectral richness, which resides in the bands [14]. As a result, when conducting satellite image fusion, spatial information is properly added to the multispectral image, increasing the spatial resolution while maintaining spectral richness.
To conduct image fusion, some methods are available, such as the HIS transform, the main component transform, the Wavelet transform, and the Brovey transform, among others.
For this study, image fusion based on the Brovey transform is employed; this method decomposes the pixel space of the multispectral image into color and luminance components. Then the method computes a mathematical synthesis to obtain colored images and high-resolution images. Each band of colored images is multiplied by the color-band integrated proportion of the color field of the data in high resolution, resulting in an automatic re-sampling process that restores each band to a high-resolution size [15], as shown in (1).
(1) where R, G, B correspond to the bands of the multispectral image, namely, red, green and blue, and PAN represents the data from the panchromatic image; Bandi is every newly generated band.
When conducting image fusion, typically, it is necessary to compare the original image with the processed image so that the intended goal can be tested. In order to evaluate the usefulness of fused images, compared to the original multispectral images, quantitative tests exist, such as correlation, entropy, bias and the universal image Q index test [10].
Another subject of study in the present paper is heterogeneous computation. At present, the development of processing systems has focused on producing devices with simultaneous running capabilities using two approaches. The first approach corresponds to the design of multi-core CPUs; the second approach focuses on many-thread processing systems such as GPUs (Graphic Processing Units), which optimize the performance by running parallelizable processes [16]. Heterogeneous computation refers to systems that use a variety of computation units such as CPUs, GPUs, DSPs, FPGAs [8]. For the present study, a CPU/GPU architecture is defined as the heterogeneous processing system.
Parallel computing resources include a multi-core CPU as well as a multi-core GPU. In general, the CPU prepares and transfers the data, whereas the GPU is in charge of the arithmetic operations [17]. In [18], a framework for cooperative heterogeneous computation was proposed to allow an efficient use of the cores of a host CPU as well as the CUDA cores intended for the GPU. As a result, the proposed framework attained the goals of parallel running.
Based on the previous description, two main topics are identified for the present study; namely image processing using satellite image fusion techniques based on the Brovey transform, and heterogeneous processing using a CPU/GPU architecture. By combing these two topics in [19], a method for satellite image fusion was proposed based on a deep multichannel model. The experimental stage of this model involved a CPU/GPU architecture. The proposed model significantly improved the quality of the fusion process when compared to other methods in the literature. In [20], real-time parallel algorithms were proposed based on multi-GPUs intended for the processing of medical images. One of the proposed algorithms resulted in image processing 265 times faster than its serial CPU counterpart; another algorithm showed an improvement of 667 times faster processing compared to the host CPU implementation.
In this context, the present study is aimed at implementing the Brovey transform as the satellite image fusion method along with a spectral calibration stage. This proposal is conducted using a CPU/GPU heterogeneous computing architecture so that running times can be compared to their serial processing counterparts. The results are also compared in terms of visual evaluation and statistical index assessment.
This section introduces the stages of this study by presenting a general approach, the implementation of the Brovey transform, and the assessment criteria.

General Approach
process of satellite image fusion using the Brovey transform with spectral calibration. The implementation of such a method for image fusion is conducted using both serial processing (CPU only) and parallel processing (CUDA-driven CPU/GPU heterogeneous computing architecture). In order to implement the Brovey transform, Python code is used by importing several modules such as numpy for serial processing, and pycuda for parallel processing. For the evaluation stage, the following items are defined: serial and parallel running times (i.e. execution times), resulting images after applying the Brovey transform, and the set of statistical indices such as correlation coefficient, Universal Image Q Index, ERGAS and Bias (Figure 1).

Implementation of the Brovey transform
The proposed method consists in applying the Brovey transform using its typical structure. Figure 2 presents the schematic representation of this method. The method begins with two input images, namely the multispectral and panchromatic images. The first image is decomposed into its RGB bands, then the sum of the bands is carried out. Subsequently, the Brovey transform is applied by taking each of the previously decomposed bands for division into the sum of the bands, then multiplying this result by 3 (scalar) and also by data from the panchromatic image. This process results in three new bands, namely Rn, Gn and Bn. From this point onwards, a spectral calibration stage is implemented. To this end, each of the newly generated bands is processed to find their corresponding maximum (max) and minimum (min) values. The min value is subtracted from each band and the resulting data is multiplied by 255 for subsequent normalization (division) by the difference between the max and min values. Once the bands with spectral calibration are obtained, a process for concatenation of the bands is conducted, which yields the new image that holds the calibrated spectral richness from the multispectral image together with the spatial resolution from the panchromatic image. For the serial implementation of the presented method, a library called skimage was used so as to load .TIF images of any depth. Additionally, the numpy library is employed to search the maximum and minimum values, which are necessary for the spectral calibration stage. This implementation is composed mainly of for cycles, which enclose each of the steps required to complete both the Brovey transform and the spectral calibration stage.
To carry out the implementation of the Brovey transform along with the spectral calibration stage in a parallel fashion, the CPU is used as a multi-core device whereas the GPU is used as a many-core device over the a CUDA framework. Moreover, the pycuda library is required to apply this framework based on Python code. As the algorithm runs, libraries such as gpuarray, skcuda and ElementwiseKernel nget involved to allow the computation of the spectral calibration stage. Figure 3 shows the resulting CPU/GPU interaction.

Assessment
The proposed implementation was carried out using the following equipment: Intel (R) Xeon (R) CPU E-52697 v3 @ 2.60GHZ with NVIDIA Tesla K80 graphics card. Also, the images used for the assessment of the proposal were obtained from satellites IKONOS and LANDSAT. Moreover, the assessment criteria to measure processing performance correspond to the running time, a visual analysis, and the correlation between the original and the processed image. To assess correlation, the following statistical indices were computed; correlation coefficient, universal image quality index Q, ERGAS, Bias; also, spatial gain is measured.

Correlation coefficient
This index defines the correlation level between fused images and original images. This index is considered, as reliable and so has been widely used in the literature. The value of the correlation coefficient is defined in the interval [-1,1]; a value of 0 indicates no correlation, whereas a value of 1 indicates that the two images are equal. A value of -1 indicates that images are the exact opposite of one another. Equation (2) presents the formal description of this index [10]. (2)

Universal Image Quality Index Qu
This index describes the quality of processed images. It is considered as a robust index since it is based on average values and considers the standard deviation of the data. This index is defined in the interval (0,1) [10]. Equation (3) describes the function to compute this index.
This section shows the results obtained in the present study in terms of three assessment criteria, namely the serial and parallel running time (execution time), the visual impact of the images obtained after applying the Brovey transform, and the values of the various statistical indices chosen.

Running time (execution time)
Running times are presented for each type of computational architecture, also considering the size of the images. Table 1 shows the average running time in milliseconds obtained from five trials per size.   In order to show the graphical behavior obtained from the Brovey transform applied on CPU/GPU heterogeneous computing, Figure 5 presents a detailed curve of this behavior.

Resulting images
This section presents the images obtained after applying the satellite image fusion process using the Brovey transform method. The method takes the Multispectral and Panchromatic images as input and completes the process by calibrating spectral richness. For every image below, a specific geographic point was selected so that proper zoom-in frames clearly show that the resulting image captures the spectral richness of the multispectral input image as well as the spatial detail of the panchromatic input image. Figure 6 shows the resulting image after applying the Brovey transform over images with 1024 pixels. In this case, a frame containing an aircraft was selected as the specific geographic point so that the spectral richness and the spatial detail of the image can be observed.     Figure 9 shows the resulting image after applying the Brovey transform on a 8192-pixel image. The geographic point displays an industrial site located on the outskirts of Bogota. Here also, spectral richness and spatial detail can be observed.   Table 2 shows the values obtained after computing the correlation coefficient, the Universal Image Q Index as well as the ERGAS and Bias indices. These indices correspond to the processed imaged when compared to the multispectral input image; therefore, this set of results are concerned with spectral richness. The same assessment process was applied to all the image sizes used in the study.  Thus, these results focus on assessing the spatial detail. As in the previous case, indices were computed for all image sizes in the study.  Table 4 shows the correlation coefficients obtained when comparing the multispectral image and the panchromatic image so that the spatial gain can be computed.

Spatial gain
To compute the spatial gain associated to each of the bands (between the multispectral panchromatic image and the processed image), the data stored in the CC band-to-panchromatic column of Table 2 is used.
A band-to-band subtraction of the CC band-to-panchromatic values is applied to the original image using the data summarized in Table 4. Table 5 shows the results.

Results analysis
Using the results presented above, the following analysis is threefold. First, the analysis focuses on runningtime performance, namely serial and parallel processing performance. Second, a visual comparison of results is provided. Finally, the analysis focuses on the data gathered from the set of statistical metrics.
Regarding the running-time analysis, Table 1 and Figure 4 indicate that both serial processing and parallel processing exhibit similar functions, namely increasing functions of the size of images. However, significant differences in time can be observed when comparing the performance of the two types of computing architecture. A detailed look at Figure 5 (the function close-up for parallel processing) confirms the increasing function of image size. However, such a function is only polynomial increasing with small factors, which indicates that significantly lower running times can be achieved with parallel computing compared to the serial computing of the processing algorithms. Table 6 shows how many times faster the parallel algorithm is compared to its serial counterpart.

Source: Own elaboration
Based on the results in Table 6 Regarding the analysis focused on visual comparison, especial attention is devoted to the previously introduced geographical points. Fig.10 shows a visual comparison for the 1024-pixel image. Figure 10C shows how the fused image exhibits high spatial detail together with spectral richness. In this resulting image, the object can be easily recognized (the aircraft). Such an easily recognizable object appears a lot blurrier in the first image. Figure 10A corresponds to the original multispectral image, whereas Figure10B corresponds to the original panchromatic image.   (Figure 12C), allow the visualization of particular spatial details that are not distinguishable in Figure 12A, which corresponds to the original multispectral image; Figure 12B shows the original panchromatic image also for comparison purposes. Finally, the visual comparison for the case of 8192-pixel images is shown in Figure 13. The geographic point chosen for this set of images corresponds to the industrial site located in the outskirts of Bogota. Figure   13C exhibits great spatial detail together with color discrimination of objects due to spectral richness. Figure 13A corresponds to the original multispectral image and Figure 13B corresponds to the original panchromatic image  In addition, the ERGAS coefficient has an average value of 15.63%, which corresponds to a high percentage, sice its ideal value should also be near zero. Finally, the spatial gain indicator of Table 7, based on the value of a correlation coefficient, shows the average spatial gain for each of the RGB bands.

Conclusions
The analysis presented leads to conclude that, when evaluating image processing performance in terms of running time, the implementation of the proposed algorithm using the CPU/GPU heterogeneous computing architecture exhibits significantly shorter processing times. The reductions in processing time remain significant despite observing a degree-3 polynomial time-increasing function on the size of images for both the serial and parallel processing scenarios, at work [21] they conclude that as the image size increases, time increases significantly in CPU and relatively low in GPU.. Improvements can be dimensioned by observing specific results obtained with the serial and parallel architectures, namely, for the largest image size, processing time reductions led to 532.53 times faster processing, as you can see in the work done in [22,23,24]. For this particular study, when employing the Brovey transform along with a spectral richness calibration stage over heterogeneous processing, the size of the image has negligible impact on running times. This can be observed in the reduced time differences recorded in the CPU/GPU scenarios for different image sizes.
The analysis derived from visual comparison between the original multispectral image, the panchromatic image, and the processed image indicates that the observable changes in images are significant, especially when comparing the resulting image with the input multispectral image. This is achieved by incorporating better spatial definition, leading to images in which objects can be clearly identified, also preserving and calibrating the spectral richness from the multispectral input image. Moreover, this level of spatial detail combined with a high degree of spectral richness provides a broad and detailed source of information that can be used in other application studies such as those on object detection, object recognition and object prediction.
By gathering the values of several statistical indices used in both spatial and spectral assessment, it can be concluded that, for the proposed implementation, the correlation coefficient applied on a band-to-band basis shows no significant differences when changing the size of images. The same is true for the correlation coefficient in the case of the band-to-panchromatic comparison. Therefore, the size of the image has no direct impact on the correlation between the original images and the processed images. However, when comparing these two correlation indices, band-to-band vs band-to-panchromatic, it can be observed that, on average, the CC of the band-to-panchromatic case is higher than the band-to-band CC. Hence, the