ANALYSIS OF CLOUDINESS BY SEGMENTATION AND MONITORING OF SATELLITE MAP IMAGES

The algorithm of the dynamic threshold segmentation of images using clipping plane in a three-dimensional XYZ image space is proposed. To build the clipping plane of the dynamic threshold the precession and nutation angles as the base threshold values are found. The developed algorithm is applied to the satellite map images to get cloudiness intensity. The satellite map images are transformed by segmentation and inversion. The segmented and inverted images are scanned to receive the distributed cumulative histograms. By the help of so-called cloudiness meter the statistical data is processed for calculation and monitoring of cloudiness in Ukraine. The formulas to create an image of the distributed cumulative histogram are considered. Formulas to reconstruct images of the rotated satellite map images are proposed. The satellite weather map images were taken from the Wunderground services. The clustering algorithm is used to classify the regions of Ukraine by cloudiness intensity, which were created distributed cumulative images. The clustering algorithm is based on the agglomerative procedure.


INTRODUCTION
The analysis of satellite maps allows us to determine the location of thunderclouds, hurricanes, fronts, fogs, clouds.Weather phenomena change in time and space.Therefore, sequences of maps can be used for their analysis, which allows predicting further behavior and movement of phenomena within a short period of time [1][2][3][4].
Satellite images are infrared images transmitted from geostationary satellites.Satellite images show the distribution of clouds.The horizontal resolution of digital images is variable and depends on the position of the satellite.In the central part of the Meteosat map, one pixel corresponds to a square of five kilometers, and with a smaller angle of observation, the square is increased to eight kilometers.On the maps, the colors indicate the height of the cloud mass: white, blue, brown and provide information about the height of the peaks of the cloud.
In order to monitor and analyze the clouds in Ukraine, we need to have access to a service that regularly provides an image of the surface, as well as software for further processing of this data.Two components to process images are the segmentation and clustering algorithms.According to the survey in [5], we can divide the first ones into two classes: those based on the finding of the threshold of intensity, as well as those in which the division is carried out by the allocation of regions of the image with certain properties.In the first one, the thresholds of intensity are determined on the basis of histograms.Among them, there are algorithms for searching, entropy [6] and minimum error algorithms [7,8] and others.The method of separating the graph into parts [9,10] represents the methods of the second class.Most algorithms are quite cumbersome, especially those that use patterns of graphs or are based on statistical calculations.In addition, relatively simple traditional threshold segmentation algorithms are not suitable for segmentation of images with significant defects in lighting or exposure.
Many papers present a survey of the latest image segmentation techniques using clustering [11][12][13].
Fuzzy, C-Means, K-Means clustering are the most wide-spread clustering approaches to image segmentation because of their robust characteristics for data classification.The representation of images by an adequate number of clusters (objects in the image) can better present their contents.The process of clustering is used in pattern recognition to roll up the characteristics of images in their search [14][15][16].
The paper [17] proposes multilevel image thresholding for image segmentation using several P-metaheuristic algorithms, including whale optimization algorithm, grey wolf optimizer, teaching-learning-based optimization algorithm and some others, which have been presented recently.A new method of multilevel thresholds for image segmentation using grey wolf optimizer is proposed in [18].This metaheuristic algorithm is applied to multilevel threshold problem using Kapur's entropy and Otsu's between class variance functions.Based on the mentioned approaches the geoinformation technology of cloudiness analysis is being created.Some examples are considered in electronic resources and publications [19][20][21].
This paper is based on the research provided by the grant support of the State Fund for Fundamental Research (project N 33651).

SEGMENTATION BY INCLINED PLANE
We consider the method of adaptive segmentation of images, in which the segmentation plane is constructed at an angle to the plane OX , OY (Fig. 1).To construct the plane of the cut, it is necessary to have the coordinates of the three points through which it will pass.

Figure 1 -Segmented image by inclined plane
Let three points in space be given: O x , y ,z ,   , .

A=y z z +y z z +y z z B=z y y +z y y +z y y
C=x y y +x y y +x y y D= x y z y z +x y z yz +x y z y z The arbitrary point z of the found plane is calculated by the formula: where w is a value of the coordinate along the axis OX ; h is along the axis OY .Now it is possible to take over all pixels of the image to determine if the intensity of this pixel is less or more than the plane value at a specific point.
, then the pixel replenishes the list of the light segment.
As it can be seen in the formula, the input data for the division algorithm are the coordinates of three points in the space of the intensity of the image pixels.In fact, these points should set the local thresholds for segmentation.The task is complex as from the point of view of algorithms for determining the threshold, and as for analyzing the surface of the image to select the points themselves and their coordinates.One possible method is to find the thresholds in the regions of the image (the choice of region forms is a separate task) and the linear approximation of the values of the thresholds of the plane 0 = D + Cz + By + Ax .The simplified procedure for determining the initial points for plotting is to select three points out of four possible, namely: at the corners of the rectangle of the image:

CALCULATION OF ROTATION ANGLES
We divide image A into four equal parts (regions) 11 A (Fig. 2a), and let us assume that their thresholds of segmentation are 11 T , 12 T , 21 T , 22 T (Fig. 2b).We find the mean values of the thresholds for segmentation for the left-right and the upper-lower orientation (two of them are given in Fig. 3).Also, we find the total threshold value as the average value of the thresholds of all parts of the image: We find angles α, β of the rotation.To do this, consider triangles   From here we find the angles: Now, let consider the plane of the parallel plane XOY and the normal to it (the vector [0,0,1] is directed upwards (on the axis of intensity ( ) Z I ).This plane can be represented as an equation: where are the coordinates of the normal vector of this plane.The equation of a plane parallel to the plane XOY with a normal vector [0,0,1] is written as follows: where D is some constant.
To obtain a dynamic threshold, we need to turn the normal [0,0,1] so that the equation based on it is a plane that reflects the dynamic transition of the threshold value.To do this, we use the values of the angles that were calculated before.Using the data angles ( , )    , we can obtain matrices of affine transformations (rotations) relative to axis OY for angle α ; relative to the axis for angle β .
Since we will perform both rotations of the normal vector[0,0,1] , we find the composition of affine transformations: cos sin sin sin cos 0 cos sin sin cos sin cos cos After finding the composition of the transformations, we apply it to the normal vector [0,0,1] : 0 cos sin sin sin cos 0 sin cos 0 0 cos sin 0 sin . 1 sin cos sin cos cos 1 cos cos The plane equation, which will reflect the dynamic threshold, will be the following: The final formula for calculating the segmentation threshold ( 2) is reduced to the form: sin cos sin cos cos where, α , β are the angles of rotation of the plane relative to axis OY (nutation) and OX (precession), T is the value of the total threshold of segmentation, w is the difference between running point x and half width, h is the difference between running point y and half height.The constructed plane will have the form as it can be seen in Fig. 5b.For comparison, the value of the constant threshold is given in Fig. 5a.

CLOUDINESS IMAGES TRANSFORMATION
The satellite maps of Ukraine with cloudiness are being obtained from a free resource [3] every three hours.They have cloudiness on the neighboring territories and blue colors of the neighboring seas.The developed algorithm was applied for three satellite images without and with cloudiness in visible spectrum.The original images and segmented ones are given in Table 1.

Original images Segmented images
Cumulative histograms of images are given in Fig. 6.In this figure we can see three charts: for the original images.Neighboring sees and countries have influence on histograms.

Figure 6 -Cumulative histograms of satellite map images
To exclude the influence of the additional objects on the result of the analysis, we cut and segment all the objects outside Ukraine as well as outside cloudiness.On the map, the territory of Ukraine is presented as a closed region outlined by a black contour.Therefore, we use a well-known algorithm of image filling in four directions with a stack [16].The color of filling comes from the input image or is provided by users.We choose the black color as having no influence on cloudiness intensity.The segmented and black color images are given in Table 2.

Table 2. Segmented and colored satellite map images
Cumulative histograms of segmented images are given in Fig. 7.We can see that cloudiness pixels are within the intensity interval from 120 to 200 .These pixels are the object of analysis.

DISTRIBUTED CUMULATIVE HISTOGRAM FOR ROTATED IMAGES
We present the technique helping to perform cloudiness monitoring.At first we calculate two distributed histograms as two sets of ( ) N M ordinary histograms (for every column and row of the image pixel matrix): ( ) { ( )}, 0, 255, 1, The distributed histogram shows frequency of pixel intensity values in columns ) (c V i and in rows ) (r V j .In the image histograms, the OX axis corresponds to the grey level intensities ( 0 255  ) in N columns ( M rows) and the OY axis corresponds to the frequency of these intensities.
Then we calculate two distributed cumulative histograms as sets of frequency sums: 0 ( ) ( ), 0, 255 , 1, where ( )  For further processing of the DCH we present it by a flat 2D image on the plane OX , OI is a top view on the three-dimensional distributed histogram along the OI axis as it is seen in Fig. 8.In the new image, each value of the pixel intensity corresponds to the pixel frequency in the columns or rows given by the DCH in Fig. 8: ( ) 255 ( ) / , 0, 255, ( ) 0 ( ) 255 ( ) / , 0, 255, ( ) 0 For three satellite images from Table 1, three DCH are presented in Fig. 9.In Fig. 9 it is difficult to distinguish pixels responsible for different objects.To make them to be more visible we fill the closed regions of white and black colors by grey color using the flood-fill algorithm [18].So, we get an informative part of the distributed cumulative histogram ( IDCH ).In Fig. 10  Histograms in Fig. 10 have different sizes and intensity.Our task is to measure these properties.Moreover, we want to get similar information from different points of viewing.Depending on viewing angle distributed histograms of one map image demonstrate different cloudiness for different coordinates being in a sum as a constant value.Then, we can compare cloudiness features got from different rotation angles and determine directions of the greatest cloudiness intensity.So, we organize socalled cloudiness meter.But instead of moving a camera around the image (Fig. 11) we rotate the sample image inside the circle.
The calculation and processing begins with the selection and loading of the image to be analyzed.The value of the RGB image components is averaged, so the pixel component value R G B    is calculated by the formula: , (24)   where , ,  The user has the opportunity to specify the parameters of the processing, including its mode, the number of rows or columns taken into account for each histogram, the size of the image and the angle of its rotation.
Rotation of the original image is carried out by applying combinations of affine transformations and changing the size of the canvas image.The rotation is automatically applied to displaying the original image in the program window.
Dimensions of a new canvas for an image after rotation can be determined using the formulas: where d is the diagonal length of the original image,  is the angle between the diagonal and the base of the original image,  is the angle of image rotation, , A B   are the width and height of the canvas of the rotated image.
To perform affine transformations to the original image, the following transformation matrices are used: . (26)   Passed through each pixel of the image, the pixel intensity value is stored in the array corresponding to its row (or column in the case of column analysis).The length of each array is determined by the dimension of the pixel color component.For this paper it is equal to 256 .When counting, transparent pixels ( 0   ) are thrown away.The process of passing can be represented by the following expressions: , (27)   where [256]    In the case that each histogram consists of several rows (columns) of the original image, the passage process will have a slightly different look: The resulting values of rows (or columns) are accumulated in arrays of values for cumulative histograms constructing.Each of the cumulative histograms can be represented in the program window, with the slider allowing to select the row (column) for which the cumulative histogram will be built: M M is the cumulative histogram, constructed in rows and in columns of the original image, respectively.
On the basis of the constructed cumulative histograms, a new image is formed.In this case, each cumulative histogram determines one row (column) of the new image.In case of processing the rows, the height of the new image will coincide with the height of the output image, and the width is determined by the number of elements in the cumulative histogram.For column processing, the height and width values will be reversed.The value of each pixel of this row is formed on the basis of the corresponding element value of the cumulative histogram.Correspondingly, for rows processing the new image will have the size of the B N  pixels, and the pixel values will be determined by the formula: . (30)   When processing the columns, a new image will have a size N A  , and the pixel value formula will look like the following: , (31)   where , , , (32)   where , A B   are the new image sizes.All the previous formulas are implemented in software allowing rotating the map image to build the distributed cumulative histograms (for example, in Fig. 12).

CLOUDINESS INTENSITY MEASUREMENT
We consider methods for measurement of cloudiness intensity in different directions and estimated directions with maximal cloudiness.
A simple method is based on a formula to calculate a mean value of pixel intensity in columns and rows of the image pixels matrix: where , i j b is pixel intensity in i -row and j -column ), W and H represents the number of columns and the number of rows respectively.
To exclude the influence of large values of white pixels we invert images of distributed cumulative histograms (Fig. 9, 10) and use the inverted images (Fig. 13) for calculations.Then we apply a number of pixels n as a width of clouudiness corridor.With it for every chart we seek a coordinate interval having a maximal intensity value by the following formula: 0 ( ) ( ) max ( ), 0, 179 where W is a width of an image, n is chosen from a query for the needed territory band.Formula ( 35) is being applied to all distributed cumulative histograms got from rotated images (Fig. 12).For ideal case, we have 180 histograms because other 180 histograms repeat previous ones.To reduce calculation in practical case we accept a number of experiments to be smaller.Projecting the corridor on the transformed satellite map image we get the cropped part of the image having the coordinate m (Fig. 15).The coordinate m is different for all rotated images.We chose the cropped part with the greatest value of ) ( m i x  and apply for them the formulas (33), (34).New chart in Fig. 16 demonstrates how cloudiness is distributed within the cropped part.The second measurement technique is clustering.We cover an image of the distributed cumulative histogram by a grid and its rectangles become the object for measurement.
We expose an example of the DCH in Fig. 17.Black color of zero intensity does not affect the intensity value of the rectangles covering the image.So, on the clustered matrix we measure intensity by a number of indexed rectangles and intensity value of each rectangle.
z .Here xw and yh are the width and length of the image, the height (intensity) coordinate z is unknown.Let us consider the method of finding thresholds in regions of an image of a uniform rectangular shape and linear approximation of the values of their thresholds by

Figure 2 -
Figure 2 -Division of the image into equal parts (a), and threshold values in a three-dimensional space (b)

Figure 3 -
Figure 3 -Threshold values for segmentation for the width (a) and height (b) measure of the image (width x or height y ), and,

Figure 4 -
Figure 4 -Angles of rotation: (a) relative to the axis OY (nutation), (b) OX (precession) -based on the found threshold values

Figure 5 -
Figure 5 -The static threshold (a) and the adaptive threshold of segmentation (b)

Figure 7 -
Figure 7 -Cumulative histograms of segmented and colored map images from Table 2 frequency in column and row, N , M are the numbers of columns and rows.A schematic example of the distributed histogram is given in Fig.8.Colors of the top planes are proportional to heights of parallelepipeds.

Figure 8 -
Figure 8 -Distributed cumulative histogram of abstract image

Figure 10 -
Figure 10 -Informative distributed cumulative histograms of three cloudiness images: (a) -without clouds, (b) -middle clouds, (c) -big clouds are the color component values of the new image pixels, , A B are the width and height of the original image.

Figure 11 -
Figure 11 -Scanning the sample image to get cumulative histograms.
for storing histograms data, constructed in rows and in the columns of the original image, respectively, ij R is the value of the color component of the original image pixel, ij  is the value of the transparency of the original image pixel.

Figure 12 -
Figure 12 -Examples of the image positions to scan them from OX plane by different angles and calculate the DCH

Figure 13 -Figure 14 -
Figure 13 -Inverted distributed cumulative histogram for image by 0 degrees Now we apply a formula for mean value of pixel intensity in rows of image matrix (4) to the inverted distributed histogram of cloudiness got for different angles of observation: 0, 90 and 45,135 degrees.As a result we have charts in Fig. 14.

Figure 17 K
Figure 17 -Distributed cumulative histogram for image by 90 degrees In Fig.18, one can see two clustered matrixes of the DCH shown in Fig.13and Fig.17 . Lower parts of matrixes demonstrate measured intensity of cloudiness from one distributed cumulative histogram.The 179 histograms are being clustered to get data for all cloud distribution.

Figure 18 -
Figure 18 -Clustered matrix of distributed cumulative histogram