=Paper=
{{Paper
|id=Vol-2210/paper46
|storemode=property
|title=Investigation of hyperspectral image pixel signatures by the empirical mode decomposition method
|pdfUrl=https://ceur-ws.org/Vol-2210/paper46.pdf
|volume=Vol-2210
|authors=Pavel Pahomov,Alexander Borusyak,Vadim Turlapov
}}
==Investigation of hyperspectral image pixel signatures by the empirical mode decomposition method==
https://ceur-ws.org/Vol-2210/paper46.pdf
Investigation of hyperspectral image pixel signatures by the
empirical mode decomposition method
P A Pakhomov1, A V Borusyak1 and V E Turlapov1
1
National Research Lobachevsky State University of Nizhny Novgorod, Gagarin Avenue 23,
Nizhny Novgorod, Russia, 603950
Abstract. The signature of hyperspectral image (HSI) pixels and their decomposition into
empirical modes (EM) and low-frequency residuals are investigated. On the basis of estimates
related to the EM-decomposition method, the possibility of switching from a 2-byte
representation of the values of the HIS-signature to a 1-byte one is examined using the example
of the Moffett Field from the AVIRIS spectrometer. It is revealed that the localization of the
minimum window sizes for the first EM is correlated with the localization of the significant
influence of the atmosphere; the first low-frequency residues have a fairly high correlation
coefficient with the signature and the first 2 of them and their EM are most interesting for use;
50 of the 224 HIS-channels are noisy and can be excluded from consideration; EM with
practically no loss of accuracy can be reduced to a 1-byte representation. The management of
the classification capabilities of signatures by changing the threshold value of the correlation
coefficient with the sample, as well as the application of the 1st and 2nd low-frequency
residues in place of the signature, was studied. Classification capabilities of signatures in a 1-
byte representation are almost equivalent to a 2-byte one, which makes it possible to put a
signature with 1-byte representation as the object of compression. For the wavelet
decomposition of the HSI data array, in combination with a 1-byte representation, a
nearlossless compression ratio of 6.65 is obtained.
1. Introduction
Hyperspectral image (HSI) compression is the central task in the processing of such images. The first
technological stage, where this task was of vital importance, was the transmission of Earth remote
sensing (ERS) data to Earth, since it was not practicable to store and fully process such large amounts
of data onboard a spacecraft/aircraft. Thus, in [1] published in 1997, the algorithm of Context-Based
Adaptive Lossless Image Coding (CALIC) was proposed. This algorithm was two-dimensional, i.e. it
compressed images channel-by-channel and, therefore, it was called 2D-CALIC. By that time, it had
been shown that the classified context adaptive prediction implemented in the adaptive selection of
adaptive predictors (ASAP) and adaptive combination of adaptive predictors (ACAP) methods yielded
good results. The ACAP method was equipped with three-dimensional predictors obtained by learning
based on land data. Both methods were applied in the well-known Airborne Visible/Infrared Imaging
Spectrometer (AVIRIS), where a data density of 6-5 bit/pixel was achieved. Another popular approach
is based on the development of nonlinear predictors. This was implemented in the LOCO-I algorithm,
standardized as JPEG-LS [2], and in the above-mentioned 2D-CALIC [1], which was soon upgraded
by the authors to 3D and was called 3D-CALIC [3]. In the 2004 publication [4], an optimized version
of the 3D-CALIC method, called M-CALIC (M-multiband), was proposed. It was based on the multi-
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018)
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
channel nature of context-based prediction. The transition from the lossless compression algorithm to
the near-lossless one in the M-CALIC method was performed in accordance with [5] by applying the
sampling of the values with a step (2 δ +1), where δ is the permissible rounding error for pixel values.
The choice of the preferable type of compression very much depends on the applications where the
HSI thus stored will be used. The papers on the interpretation of hyperspectral aerospace
measurements for diagnosing the state of natural and technogenic objects, coastal waters, and crop
areas [6], [7], [8] underline the importance of a qualitative improvement in the channel resolution in
the HSI for precise classification of the type and state of the objects being observed, and note that
these problems cannot be solved by means of multispectral images alone. At the same time, the high
correlation of neighboring channels cannot be seen as a reason for excluding one of these channels as
uninformative because the entire small fraction of the differing information is necessary to accurately
determine the object of interest or its state.
One of the first approaches in the development of compression algorithms for HSI is to optimize
the number of spectral channels while preserving the information value of hyperspectral imaging. This
approach is reflected both in domestic [9] and in foreign sources [10], [11]. By the early 2000s, two
main lines of research were defined to address the problem of optimizing the number of HSI channels
[10]: 1) Feature Extraction; 2) Band Selection. With the possible exception of [12], these two lines of
research have in common the possibility of detecting the most informative channels at the initial stage
by the projective optimization method. In this method, the optimal approach is to project the initial
HSI X of the size n × N (n is the number of channels, N is the number of pixels per channel) into a
new HSI Y of a smaller size m × N (m < n), which maximizes the projection index J=J(Y), where
Y AT X , аnd A is an n × m matrix. Usually, the value of the index is estimated through the spectral
data variance by the principal component method [9], [11], and the matrix A is formed from the
columns (eigenvectors) of the covariance matrix of the original data. Sometimes, after applying the
principal component method, the method of independent components, well-proven on non-Gaussian
distributions, can be applied [9],[13].
The divergence of information J between two pixels of the HSI was chosen as a criterion in [14] for
the synthesis of images reflecting the spatial distribution of the amplitude ratios of like pixels in
different channels. Also of interest is the integral estimate for the amount of information proposed in
[9] taking into account the signal-to-noise ratio.
In a series of works of 2008-2013 [15], a three-stage algorithm for lossless HSI compression was
proposed, including:
1) Taking into account the relationship between HSI channels by calculating the correlation,
constructing a linear predictor of the next channel value from the previous one, and forming arrays
of deviations of the value of the next channel from the predicted one.
2) Forming an auxiliary data structure for storing unique pairs of groups of element values in a 1-
byte representation, as well as pointers to these groups.
3) Compression of the data obtained after the transformations by means of a standard entropy
algorithm by processing the generated auxiliary data structures.
An approximately 40% gain in compression was achieved compared to JPEG-LS. It was
established that about 45% of the gain was obtained by using a 1-byte representation of deviations
from the predicted values; 15-26% of the gain was obtained due to the reordering of the channels
during compression.
In a series of papers published in 2009-2016 [16],[17],[18],[19] a method of hierarchical lossless
compression for the purpose of storing HSI was constructed. The following requirements for the
method were formulated: 1) the possibility of compression of multi-channel images; 2) quick access to
fragments of compressed images at various scales; 3) low computational complexity of
decompression; 4) strict error control; 5) high efficiency in "no-error" mode and with small errors; 6)
use of interchannel dependencies; 7) quick access to the specified components of compressed images
at various scales; 8) compression of 16-bit images.
The following dependencies on the channel number are investigated and used for compression: 1)
the correlation coefficient with the next channel; 2) the channel average E; 3) the difference between
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 353
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
the minimum and maximum in the channel; 4) channel variance D. A method of "common reference
components" for channel groups in combination with the "sliding approximation" method within the
group is proposed, which permits to independently compress groups of neighboring channels. The best
achieved value of lossless compression ratio was about 3.6; the recommended number of hierarchy
levels was 4. The nature and the range of value changes in the channel X i max ( xki ) min ( xki )
k k
were shown depending on the channel number i. Thus, for the image "Urban and Mixed Environment"
(SpecTIR spectrometer), the maximum value of the changes X was about 32000, while for the
"Cuprite-1" image (AVIRIS) - about 12000. In this case, the change curves X i , E ( X i ) , D( X i )
look more informative for object classification than the correlation ratio curve.
The results of the research into the methods of lossless HSI compression as of 2013 were fixed in
the relevant standards [20], [21]. In particular, the values of the lossless compression ratio of the order
of 4-5 were stated as very high figures for HSI.
The research performed during the last decade into the role of noise in the compression of HSI,
including its role in near-lossless compression and compression with losses, is covered in a series of
publications [22], [23], [24]. It was found that: 1) the noise in the hyperspectral images is signal-
dependent and is almost uncorrelated spatially; 2) the noise parameters in adjacent channels are
generally close, although the overall range of variation of the signal-to-noise ratio over hyperspectral
images is very wide; 3) in almost all hyperspectral images the quality of approximately 80-85 percent
of the total number of channels is close to ideal, that is, the PSNR is close to or above 35 dB.
In this situation, it would be of interest to explore some universal methods capable of analyzing the
signatures of both individual pixels and their integral forms, in terms of classification with respect to
object and background or for detection the presence of noise, as the classification of the object and the
background or for the presence of noise detection.
2. Description of the empirical mode decomposition method
In our opinion, one can use for the above purposes the method of empirical mode decomposition
(EMD). The empirical mode decomposition method was published in 1998 [25], then it was adapted in
2008 for images [26] and is currently applied for a number of image processing tasks [27]. In
particular, it can also be applied for hyperspectral images to decompose the signature of an image
pixel, similar to a signal, into intrinsic modes having a space/time-varying periodicity. Relying on this
property, we want to isolate the highest-frequency part of the signature, which by its nature can
represent noise, and also allocate channels (channel groups) that determine individual attributes of the
object class specified by the sample pixel.
The decomposition into empirical modes is based on the following assumptions for the signal: 1)
the signal has at least two extrema: one maximum and one minimum; 2) the characteristic time scale is
determined by the interval between the extrema; 3) if the data is completely devoid of extrema (trend),
but may contain inflection points, then the signal can be divided into parts in order to reveal the
extrema.
The time interval between successive extrema is taken as the determination of the time scale for the
intrinsic oscillation mode, since it not only provides a much higher resolution, but it can also be
applied to the data with a non-zero mean (positive or negative values, without zero crossings). In our
case, the axis of the channel numbers plays the role of the time axis.
The mode extraction method, called sifting, is described as follows [25]. The decomposition
method uses envelopes built on local maxima and minima separately. For this purpose, the local
extrema of the f(t) signal are identified, and all maxima are interpolated by the cubic spline line as the
upper envelope U(t). In the same way, the lower spline envelope L(t) is constructed on the minima.
The mean between the upper and lower envelope R(t)=(U(t)+L(t))/2 receives the low-frequency
residue status and is used for further transformations as f(t), and the difference of the functions f(t) and
R(t) receives the status of the first empirical mode 1 (t ) .
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 354
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
The empirical mode is a function with the following properties: 1) the number of the function
extrema on the interval considered differs by no more than one from the number of zero crossings of
this function; 2) the half-sum of the upper and lower envelope of the function is close to zero.
As a result of the decomposition of the signal f(t) into empirical modes, we get:
N
t)
f( r
(t)
i(
t), (1)
i
1
where i (t ) are empirical modes, r(t) is the trend residue. The first modes contain high-frequency
components of the signal, and the last modes and the residue are low-frequency components.
Fast decomposition algorithm
We use the fast adaptive decomposition method proposed in [26]. Since we are considering a
signature, the signal is one-dimensional: it depends only on the channel number. Accordingly, the
algorithm is reduced to a one-dimensional version, and it is also simplified in comparison with [26]
and [27] while preserving the idea of the original method [25]. Namely, the operations of constructing
envelopes (upper and lower), and then, calculation of the points of the mean curve R(i) between them,
are replaced by smoothing (averaging) over a symmetric window of width w.
Algorithm steps:
1. Assign the signature of the k-th pixel as the signal f (i) ,
2. Initialize the size w of the processing window with a value of 3 (w =
3), the number of the empirical mode is q=0.
3. Calculation of low-frequency residue
1 iw/ 2
R(i) f ( j) ,
w j i w / 2
(2)
the values corresponding to the positions of the window that go beyond the
limits of the signature are taken to be equal to the edge value of the
signature.
4. Construct an empirical mode:
q q 1, q (i) f (i) R(i) (3)
5. Find and enter into the arrays pU and pL all the points of local
extrema (minima and maxima) of the current empirical mode q
within the window W of width w with the center in the current channel that
satisfy the conditions:
for the array pU (i) ( j ), j Ww (i) , (4)
for the array pL (i) ( j ), j Ww (i) , (5)
where Ww (i ) is a window of width w with the center in the channel i =
Figure 1. General view of 1,…,n.
the Moffett Field HIS. 6. If the number of maxima in pU or the number of minima in pL is less
than 2, then R cannot be further decomposed and the process terminates.
7. For each local maximum, find the distance dmax along the channel axis to another nearest
maximum, for each local minimum, to another minimum d min and then, take the smallest of them:
d = min (dmin ,dmax ) , and update the window size w 2d / 2 1 (6)
8. Set f (i) = R(i) and perform steps 3-7.
Figure 3 shows an example of signal decomposition into empirical modes.
3. Experimental research
The initial data are freely available hyperspectral images obtained with the use of the Airborne
Visible/Infrared Imaging Spectrometer (AVIRIS). Figure 1 shows the general view of the Moffett
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 355
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
Field HSI in three channels: 29th - R = 638.2nm; 19th - G = 540.6nm; 11th - B = 462.8nm. This
sample HSI shows various types of terrain, including water bodies, green vegetation and urban areas.
As it was mentioned above, when constructing and analyzing compression algorithms for almost any
data, the question arises about the presence of noisy data to which noise was added during the
registration process or those data that were distorted as a result of transmission through the
communication channels. The volume of such data can be significant and it can affect the resulting
compression ratios. In some cases, information noise can be filtered in one way or another by data
processing algorithms. Therefore, it is quite important to remove this noise before the compression of
the original image. The chosen test sample is stored in the ENvironment for Visualizing Images format
(ENVI) and has 2 bytes of information per pixel in each channel. Its width is 752 pixels, and its height
is 1924 pixels. In total, the chosen HSI has 224 channels (about 10 nm wide each), which cover the
wavelength range from 0.365 μm to 2.497 μm. To investigate the correlation of images of adjacent
HSI channels, a linear Pearson correlation coefficient was calculated between the current channel and
the next one (223 pairs):
rXY ( X X )(Y Y ) . (7)
( X X ) 2
(Y Y ) 2
The result of the calculation is shown in Table 1.
Table 1. Distribution of the values of the Pearson correlation coefficient for adjacent HSI channels.
Range Number of channels In percentages, %
0,9999-0,99999 32 14.34
0,999-0,9999 91 40.80
0,99-0,999 48 21.52
0,9-0,99 40 17.93
0,3-0,9 12 5.83
As can be seen from Table 1, over 94% of all channels have the correlation coefficient exceeding
0.9, and over 76% of channels have the correlation coefficient exceeding 0.99. This suggests that the
neighboring channels are usually very similar to each other, despite their modulation by the influence
of the atmosphere. Channels in which the correlation coefficient with the adjacent channel is lower
than 0.99, mainly fall within the wavelength ranges (near 1.4 µm, 99-128 channels of our HIS, and 1.9
µm, 153-166 channels) of almost complete absorption of light by a mixture of water vapor and carbon
dioxide. Approximately the same high absorption there is in the band above 2.5µm (channels 218-
224). In this HSI it is about 50 channels (or 22% of this HSI), the contents of which are almost
impossible to use because of the high noise level. If necessary, these channels can be deleted without
loss to solve problems of classification of objects of the image.
3.1. An example of a decomposition of a signature into empirical modes and low-frequency residuals.
Investigation of the possibilities
As standard for HSI storage is used 2 bytes per pixel per channel. The range
of values of this type of data is significantly wider than that of universally
used formats for storing visual information. One of an important reason for
the introduction of the 2-byte accuracy of HSI is the high accuracy of modern
spectrometers, which is achieved through special calibration methods [28].
At the same time, a fairly authoritative source [29] claims that already
because of the influence of the atmosphere on the measurements of the
spectrometer, it is difficult to talk about achieving accuracy higher than 2-
5%. This makes relevant the question of how many bits are needed to store
one pixel in the HSI channel to exclude losses for future classification of
HSI-objects, and to determine their state.
Figure 2. The pixel
The application of the method of empirical modes will be considered
chosen for constructing
using the water signature as an example. The pixel of the water sample is
empirical modes.
marked in red in Figure 2. Figure 4 shows an important for further analysis
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 356
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
part of the result of the decomposition of the water signature into empirical modes and low-frequency
residuals (R). The complete decomposition of the signature contains 15 empirical modes and a trend.
Table 2 shows a number of dependencies on the number of the empirical mode (EM) for the first
eight modes: 1) for the window size w; 2) for the correlation coefficient of the original signature with
its corresponding low-frequency residue; 3) for maximum deviation of the empirical mode from zero
(Max) in absolute units of a 2-byte representation and in percentage to the maximum value of this
signature and the maximum value throughout the HSI. Values in percent are given for comparison
with an estimate of the achievable accuracy of HSI in 2-5% given by D.Landgrebe (1999).
a)
200
100
0
b) c)
-100
100
100
0
0
-100
d) e)
-100
Figure 3. Decomposition of the signature into empirical modes and low-frequency residuals: a –
signature (f), first (R1) and second (R2) low-frequency residuals; b-e – empirical modes 1-4.
The values of the Pearson correlation coefficient between signature and low-frequency residues are
quite high: about from 0.999 to 0.99 (see table 2). Specifically, the significance of these values for the
classification of HSI objects will be considered below, but on the whole this indicates that the
signature in the correlation estimates can be replaced by at least its first low-frequency residues
corresponding to the required accuracy.
Table 2. The change in the window size and in the maximum of the empirical mode (EM) with the
its number up to 8 steps, the correlation coefficients of the signature S with its low-frequency
residuals R1-R8.
Empirical Mode (EM) Number 1 2 3 4 5 6 7 8
Window size w 3 3 3 3 3 3 7 10
Correlation coefficient rSR 0,9991 0,9987 0,9982 0,9979 0,9976 0,9973 0,9950 0,9906
Max of the EM 306 99 47 36 26 21 267 215
In [%] to Max of the signature 9.5 3.1 1.5 1.1 0.8 0.7 8.3 6.7
In [%] to Max of HSI 4.1 1.3 0.6 0.5 0.3 0.3 3.6 2.9
In the first 6 steps of the decomposition, the algorithm automatically selects a minimal blur window
with a size of 3 channels. This may indicate the presence of noise, at least in part of the channels of
our HSI. There are certain successes in correcting the influence of the atmosphere on the hyperspectral
image [30], [31], but for channels where the signatures are close to zero, the correction is practically
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 357
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
useless. From the point of view of the compression of the HSI, it is advantageous to detect noisy
channels and remove them from consideration. Let's consider, whether assignment of the minimum
size of the window can serve as the detector of noisiness of the channel.
Figure 4. Distance (d) between adjacent extrema of the 1-st empirical mode: maxima (points with
fill) and minima (point without fill), via the minimum of which the window size w is determined. The
horizontal axis is the channel number.
To do this, consider the distances between the current and next extremums (between two adjacent
maxima, dmax, or minima, dmin), by which and according (6) the window size w is determined along
the axis of the channel numbers. The zone of noisy channels should
give a smaller local window size than the local window size in the
"clean" channels zone. The minimum possible value of d, for which
w = 3, is 2. Figure 4 shows the distances between extrema for the
first empirical mode.
As can be seen from the graph, the presence of the values d = 2 is
not limited to intervals 99-128, 153-166, although in these intervals
they are the most. The appearance of d = 2 beyond the limits of these
intervals correlates well with the windows of the influence of the
atmosphere, but do not lead to a completely noisy channel. As a
result, via the values of d=2, we can only select the channels of
suspects for noise and channels that are substantially modulated by
the influence of the atmosphere.
3.2. Classification capabilities of signatures and their low-frequency
residues on the basis of values of the correlation coefficient
As is known, each unique object of HSI has an individual signature.
The difference between image objects is manifested in the difference
in the behavior of signatures in some sequence of channels.
Let us consider the possibilities of classifying objects of the HSI
via the signature of the object sample, using for this purpose different
values of the Pearson correlation coefficient. Let us also consider
how the transition to a 1-byte representation of the signature or the
Figure 5. Samples of use of its low-frequency residuals instead of the signature will affect
objects for classifying. the classification capabilities. Let's select several pixels-samples of
objects belonging to different objects (Figure 5). Sample 1 represents vegetation (tree crowns); sample
2 - concrete (from the airport runway); sample 3 - soil sample (Soil1) in the airport area; Samples 4
and 5 are water samples (Water1, Water2), but sample 4 looks like a darker one; sample 6 - soil
sample (Soil2) from the field. Figure 7 shows the signatures that are constructed as sample averages
from signatures that have a correlation coefficient with a sample signature of at least 0.999. We
calculate the correlation coefficient of the selected signature-samples with the signatures of all other
pixels of the HSI. We shall establish the thresholds to be investigated for classifying an object based
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 358
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
on the Pearson correlation coefficient with the sample signature: {0.9; 0.99; 0.999; 0.9999}. To the
desired class we will include all the signatures that have a correlation coefficient with a sample greater
than or equal to a given threshold. In addition to the binary classification (object / background) of the
2-byte representation of the HSI, a binary classification was performed on the first and second low-
frequency residue, as well as on the signature in a 1-byte representation. Then the difference in the
classification results was calculated.
Figure 6. Signatures of samples of objects marked for classifying objects in Figure 5.
Examples of the classification of the initial HSI for different thresholds for the correlation
coefficient are shown in Figure 7. The results of using in classifying the first or second low-frequency
residue instead of the signature sample are shown in Table 3 (only for the 0.9 threshold). Comparative
characteristics of the classification results for all signature samples in a 2-byte and 1-byte is given in
Table 4 (only for the 0.999 threshold). These data show that the classification by the value of the
correlation coefficient with signature sample can be used as a classification tool. In a number of cases,
the threshold value of the correlation coefficient can play a generalizing role, as shown in Figure 7 (a
and c). The threshold (T) of the correlation coefficient equal to 0.9 for the signature of sample 1
(Vegetation) made it possible to highlight almost the entire forest park and urban area, excluding
building and road surface, and gave a classification mask that is practically complementary to the
classification at T = 0.9 for sample 5 (Water2). The threshold T = 0.9 for sample 5 (Water2) captured
besides water also the fields adjacent to the coastline, road surfaces, urban buildings. For the threshold
T = 0.99 (Figure 7, b), we see in the Vegetation class only a very small number of objects. And at T =
0.999 - only 6 pixels corresponding to the sample. For sample 5, Water2, the threshold T = 0.99
(Figure 7, d) gives already almost the entire surface of the water of the gulf and the flowing rivers, as
well as the flood areas adjacent to the shoreline. The threshold T = 0.999 (Figure 7, e, f) already
completely separates the objects specified by different samples of water 4 and 5. The value of the
threshold at 0.999 in all cases makes it possible to classify objects with signature shapes almost
completely similar to the sample signature. The threshold value of 0.9999 guarantees the selection of
objects completely identical to the sample or with the similar signature shapes, but allowing deviations
in values not exceeding 2-5%. Undoubtedly, the classification problem for the selected 6 samples can
be solved on the basis of estimates [16]: X i , E ( X i ) , D( X i ) , where i is the channel number.
Moreover, in the application of these estimations, which are intervalwise, one, two or three intervals
will suffice, in which the signatures are clearly ranked by the value of the estimates.
Let us further consider the results of using, instead of the sample, the signature of its first or second
low-frequency residue (see Table 3). Visually at T = 0.9 and T = 0.99, the difference in classification
by signature and low-frequency residuals is observed as some smoothing and narrowing along the
boundary of objects, increasing with increasing threshold. At T = 0.99, the detection density of the
object area can also decrease (loss from 3% density, to Water1, and up to 30% - to Vegetation). At T =
0.999, significant losses occur in the object.
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 359
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
a b c d e f
Figure 7. Examples of classification by object samples at different threshold values (T) of the
correlation coefficient: a) sample 1 – Vegetation, T=0.9; b) sample 1 – Vegetation, T=0.99; c) sample
5 – Water2, T=0.9; d) sample 5 – Water2, T=0.99; e) sample 5 – Water2, T=0.999; f) sample 4 –
Water1, T=0.999.
Table 3. Classification of objects by the threshold of the correlation coefficient 0.9 for signatures of
samples S in a 2-byte representation and their low-frequency residues R1 and R2.
Sample Pixels in Class Mismatched pixels
N Object coordinates via R1 via R2
via S via R1 via R2
(x;y) pcs % pcs %
1 Vegetation 442; 1332 566744 558819 557542 8169 1.4 9846 1.7
2 Concrete 447; 1275 1018521 1014854 1015257 15303 1.5 21956 2.2
3 Soil1 426; 1351 120603 116374 115699 4229 3.5 4904 4.1
4 Water1 209; 470 767970 759931 756722 8061 1.0 11310 1.5
5 Water2 315; 775 873712 864921 861536 8797 1.0 12206 1.4
6 Soil2 540; 630 294865 283321 280976 11544 3.9 13889 4.7
3.3. Ability to switch from a 2-byte representation of signatures to a 1-byte
The possibility of switching from a 2-byte representation of signatures to a 1-byte one was evaluated
on the basis of a visual comparison of the results of a binary classification and a quantitative
evaluation of the difference of the classification masks. A 1-byte representation is obtained from a 2-
byte division of double byte values by 32 with rounding. In all cases of visual comparison, the
difference was not noticeable. The difference in masks in all cases remained within 7.5% (the worst
case for the Water2 sample, see Table 4). The unmatched pixels in all cases were distributed over the
area of the classified object in proportion to the density of the object.
Correlation coefficients between the signatures in the 2-byte and 1-byte representations for the
samples: Vegetation, Concrete, Water1, Water2, Soil1, Soil2. Their values were 0.99997, 0.99998,
0.99994, 0.99995, 0.99998, 0.99999 respectively. The correlation coefficient decreases with increasing
divisor. For example, for the Water2 sample, when the divisor is changed in the sequence {16, 32,
64}, the correlation coefficient with the 2-byte signature, while remaining high, changes as follows:
0.99998; 0.99995; 0.99980.
This suggests that it is possible to switch to calculations with a 1-byte representation, if necessary,
saving the divisor for backwards compatibility with a 2-byte representation.
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 360
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
Table 4. Classification of objects at the correlation threshold of 0.999 for sample signatures in a 2-
byte and 1-byte representation.
Sample Pixels in class Pixels in class Mismatched pixels
N Object
coordinates (x;y) (2-byte) (1-byte) pcs %
1 Vegetation 442; 1332 6 6 0 0
2 Concrete 447; 1275 430 429 7 1.6
3 Soil1 426; 1351 920 904 20 2.2
4 Water1 209; 470 60972 56418 4564 7.5
5 Water2 315; 775 27176 26453 801 2.9
6 Soil2 540; 630 831 826 5 0.6
3.4. Improving the performance of calculations in the HSI classification and evaluating the possibility
of the signature-based HSI compression
In order to improve the performance of the classification, an algorithm is implemented on the graphics
processor. The working time of the classification algorithm based on the correlation coefficient with
the sample signature is defined as the average (for 10 passes) total classification time for all 6 selected
signatures. Instead of the sample signature, its low-frequency residue can be used. The measured
classification time for the central (CPU) and graphics (GPU) processors is shown in Table 5.
Table 5. Performance achieved on CPU and GPU for HSI classification via value of correlation
coefficient.
Hardware Number of Core Clock frequency, MHz Time (sec)
CPU: AMD Phenom II X4 955 4 2500 257. 07
GPU: NVIDIA GeForce GTX 970 1664 (CUDA) 1152 0.31
From Table 5, it can be seen that the use of a graphics processor makes it possible to carry out the
classification procedure 836 times (3 orders of magnitude) faster than on a central processor. In this
case, we can talk about the operation of the algorithm in real time, which makes its use in desktop data
processing applications quite realistic. The performance for the same CPU of the empirical mode
decomposition was 0.91 seconds per thousand signatures.
In the interest of evaluating the capabilities of signature-based compression, an attempt is made to
increase the number of zero or close to zero values in a compressible array by applying wavelet
decomposition. Based on the length of the signature array, 7 steps of the Haar transformation are used
for the signature wavelet decomposition of the HSI in the 1-byte representations. After that, the source
2-byte HSI and the result of the Haar transformation were archived by ZIP (without losses) with the
maximum compression ratio (Table 6). The compression ratio for 1-byte signatures was 6.65. Since
we replaced the double-byte representation with a 1-byte representation, we will treat this compression
as nearlossless.
Table 6. Compression of the entire HSI data array.
File size, MB Without wavelet decomposition 1-byte signature with wavelet decomposition
Source file 648,18 324,09
Archived file 421,12 97,48
4. Conclusion
The signature of hyperspectral image (HSI) pixels and their decomposition into empirical modes (EM)
and low-frequency residuals are investigated. On the basis of estimates related to the EM-
decomposition method, the possibility of switching from a 2-byte representation of the values of the
HIS-signature to a 1-byte one is examined using the example of the Moffett Field from the AVIRIS
spectrometer.
It is revealed that: 1) the localization of the minimum window sizes for the first EM is correlated
with the localization of the significant influence of the atmosphere; 2) the first low-frequency residues
have a fairly high correlation coefficient with the signature; 3) the greatest interest for the
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 361
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
decomposition of the signature is represented by one or two first EM and the corresponding low-
frequency residue; 4) the 1st and 2nd modes on the significant part of the channel axis are close to
zero and can be reduced to a 1-byte representation practically without loss of accuracy; 5) 50 of the
224 HIS-channels are noisy and can be excluded from consideration.
The management of the classification capabilities of signatures by changing the threshold value of
the correlation coefficient with the sample, as well as the application of the 1st and 2nd low-frequency
residues in place of the signature, was studied. Classification capabilities of signatures in a 1-byte
representation are almost equivalent to a 2-byte one, which makes it possible to put a signature with 1-
byte representation (as eight senior digits) as the object of compression.
The classification procedure is implemented on the GPU, which accelerated its execution more
than 800 times, to fractions of a second.
Classification according to the samples using [16] on the basis of estimates X i , E ( X i ) , D( X i ) ,
where i is the channel number, can additionally reduce the number of channels necessary for
classification. If the estimates are applied in a series of channel intervals, one or three intervals may be
sufficient, in which the signatures are clearly ranked by the values of the estimates.
For the wavelet decomposition of the HSI data array, in combination with a 1-byte representation, a
nearlossless compression ratio of 6.65 is obtained.
For the future works it is interesting to investigate also the approaches based on the machine
learning methods like publications [19],[33],[34],[35].
5. References
[1] Wu X and Memon N 1997 Context-based, adaptive, lossless image coding IEEE Trans.
Commun. Apr. 45 437-444
[2] Weinberger M J, Seroussi G and Sapiro G 2000 The LOCO-I lossless image compression
algorithm: Principles and standardization into JPEG-LS IEEE Trans. Image Processing 9 1309-
1324
[3] Wu X and Memon N 2000 Context-based lossless interband compression - Extending CALIC
IEEE Trans. Image Processing 9 994-1001
[4] Magli E, Olmo G and Quacchio E 2004 Optimized onboard lossless and near-lossless
compression of hyperspectral data using CALIC IEEE Geoscience and Remote Sensing Letters
1(1) 21-25
[5] Wu X, Memon N and Sayood K 1995 A Context-Based, Adaptive, Lossless/Nearly-Lossless
Coding Scheme for Сontinuous-Tone Images ISO/IEC JTC 1/SC 29/WC 1 202
[6] Kozoderov V V, Kondranin T V, Kazancev O, Bobylev V I, Scherbakov M V, Borzyak V V,
Dmitriev E V, Egorov V D, Kamencev V P, Belyakov A Yu and Loginov S B 2009 Processing
and interpretation of hyperspectral aerospace measurement data for remote diagnostics of
natural and technogenic objects Earth Research from Space 2 36-54
[7] Chaban L N, Vecheruk G V and Т S 2009 Gavrilova. Investigation of the possibilities of land
cover classification using hyperspectral imagery in thematic processing packages of remote
sensing data Proceedings of MIPT 1(3) 171-180
[8] Chaban L N, Vecheruk G V, Kondranin T V, Kudryavtsev S V and Nikolenko A A 2012
Modeling and thematic processing of images identical to video data of ERS hyperspectral
equipment under development Current Problems of Earth Remote Sensing from Space 9(2) 111-
121
[9] Popov M A and Stankevich S A 2006 Methods for optimizing the number of spectral channels
in problems of ERS data processing and analysis Contemporary Problems of Earth Remote
Sensing from Space 1 106-112
[10] Jiménez L and Landgrebe D A 1999 Supervised Classification in High Dimensional Space:
Geometrical, Statistical, and Asymptotical Properties of Multivariate Data IEEE Transactions
on Systems, Man, And Cybernetics-Part C: Applications and Reviews 28(1) 39-53
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 362
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
[11] Arzuaga-Cruz E, Jimenez-Rodriguez L O and Velez-Reyes M 2003 Unsupervised Feature
Extraction and Band Subset Selection Techniques Based on Relative Entropy Criteria for
Hyperspectral Data Analysis Proc. SPIE 5093 462-473
[12] Lee C and Langrebe D A 1993 Feature Extraction Based on Decision Boundaries IEEE Trans.
on Pattern Analysis and Machine Intelligence 4(15) 388-400
[13] Myasnikov E V 2017 Hyperspectral image segmentation using dimensionality reduction and
classical segmentation approaches Computer Optics 41(4) 564-572 DOI: 10.18287/2412-6179-
2017-41-4-564-572
[14] Nakariyakul S and Casasent D 2004 Hyperspectral feature selection and fusion for detection of
chicken skin tumors Proc. SPIE 5271 128-139 DOI: 10.1117/12.517443
[15] Zamyatin A V and Sarinova A Zh 2013 An algorithm for compressing hyperspectral aerospace
images with the account of inter-band correlation Applied Informatics 5(47) 35-42
[16] Gashnikov M V and Glumov N I 2014 Hierarchical compression for hyperspectral image
storage Computer Optics 38(3) 482-488
[17] Gashnikov M V and Glumov N I 2016 Onboard processing of hyperspectral data in the remote
sensing systems based on hierarchical compression Computer Optics 40(4) 543-551 DOI:
10.18287/2412-6179-2016-40-4-543-551
[18] Gashnikov M V, Glumov N I, Kuznetsov A V, Mitekin V A, Myasnikov V V and Sergeev V V
2016 Hyperspectral remote sensing data compression and protection Computer Optics 40(5)
689-712 DOI: 10.18287/2412-6179-2016-40-5-689-712
[19] Vorobiova N S, Sergeyev V V and Chernov AV 2016 Information technology of early crop
identification by using satellite images Computer Optics 40(6) 929-938 DOI: 10.18287/2412-
6179-2016-40-6-929-938
[20] Auge E, S´anchez J E, Kiely A, Blanes I and Serra-Sagrista J 2013 Performance impact of
parameter tuning on the CCSDS-123 lossless multi- and hyperspectral image compression
standard Journal of Applied Remote Sensing 7(1) 16
[21] Blanes I, Magli E and Serra-Sagrisat J 2014 A tutorial on Image Compression on Optical Space
Imaging Systems IEEE Geoscience and Remote Sensing Magazine 2(3) 8-26
[22] Abramov S K and Lukin V V 2015 Problems of automating the processing of hyperspectral
remote sensing images Aerospace Engineering and Technology 6(123) 101-110
[23] Abramova V V, Abramov S K and Lukin V V 2015 Multistage Iterative Method for Blind
Evaluation of Mixed Noise Characteristics on Images ITS 6(1) 8-14
[24] Lukin V, Abramov S, Ponomarenko N, Uss M, Zriakhov M, Vozel B, Chehdi K and Astola J
2011 Methods and automatic procedures for processing images based on blind evaluation of
noise type and characteristics SPIE Journal of Applied Remote Sensing 5(1) 27 DOI:
10.1117/1.3539768
[25] Huang N E, Shen Z, Long S, Wu M C, Shih H H, Zheng Q, Yen N-C, Tung C C and Liu H H
1998 The Empirical Mode Decomposition and Hilbert Spectrum for nonlinear and non-
stationary time series analysis Proc. R. Soc. London A 454 903-995
[26] Bhuiyan S M A, Adhami R R and Khan J F 2008 A novel approach of fast and adaptive
bidimensional empirical mode decomposition IEEE International Conference on Acoustics,
Speech and Signal Processing 1313-1316
[27] Guryanov F and Krylov A 2017 Fast medical image registration using bidirectional empirical
mode decomposition Signal Processing: Image Communication 1-6
[28] Introduction to Hyperspectral Imaging: Tutorial MicroImages Inc. (Access mode:
http://www.microimages.com/documentation/Tutorials/hyprspec.pdf)
[29] Podlipnov V V and Skidanov R V 2017 Calibration of an imaging hyperspectrometer Computer
Optics 41(6) 869-874
[30] Landgrebe D 1999 Information Extraction Principles and Methods for Multispectral and
Hyperspectral Image Data Information Processing for Remote Sensing 3-38
[31] Denisova A Y and Myasnikov V V 2016 Atmospheric correction of hyperspectral images using
small volume of the verified data Computer Optics 40(4) 526-534 DOI: 10.18287/2412-6179-
2016-40-1-526-534
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 363
Image Processing and Earth Remote Sensing
P A Pakhomov, A V Borusyak and V E Turlapov
[32] Denisova A Y, Juravel Y N and Myasnikov V V 2016 Estimation of parameters of a linear
spectral mixture for hyperspectral images with atmospheric distortions Computer Optics 40(3)
380-387 DOI: 10.18287/2412-6179-2016-40-3-380-387
[33] Sirota A A and Dryuchenko M A 2015 Generalized image compression algorithms for
arbitrarily-shaped fragments and their implementation using artificial neural networks Computer
Optics 39(5) 751-761 DOI: 10.18287/0134-2452-2015-39-5-751-761
[34] Sergeev V V and Yuzkiv R R 2016 A parametric model for the autocorrelation function of space
hyperspectral data Computer Optics 40(3) 416-421 DOI: 10.18287/0134-2452-2016-40-3-416-
421
[35] Kopenkov V N and Myasnikov V V 2016 Development of an algorithm for automatic
construction of a computational procedure of local image processing, based on the hierarchical
regression Computer Optics 40(5) 713-720 DOI: 10.18287/2412-6179-2016-40-5-713-720
Acknowledgements
This work was supported by a grant from the Russian Science Foundation (project No. 16-11-00068).
IV International Conference on "Information Technology and Nanotechnology" (ITNT-2018) 364