RU2523944C2 - Method of selecting extended linear objects on aerospace images - Google Patents

Abstract

FIELD: physics, computer engineering.

SUBSTANCE: invention relates to means of selecting linear objects on images. The method comprises processing an image using a finite impulse response (FIR) filter, which enables to determine points belonging to extended linear objects on the image from the positive response of the FIR filter, based on analysis of mean-square deviation values of brightness of points of the filter window when the filter window is turned relative to each analysed point on the image.

EFFECT: high accuracy of selecting extended linear objects on an image.

4 dwg

Description

The invention relates to computer technology and can be used to determine extended objects in satellite images when creating electronic maps for geographic information systems.

A known method for recognizing objects based on fractal wavelet analysis (US Pat. No. 2440608 RF, IPC G06K 9/00 (2006.01), 12/02/2010), which consists in applying the wavelet transform to the image, as a result of approximating the image array (A) a larger scale and four detailed image arrays for horizontal, vertical and diagonal details (H, V, D), applying this transformation iteratively to matrix A for each image level, thus obtaining a multi-level decomposition of the image into m sshtabam, 2-fold, allowing on the basis of the fractal dimension analysis and object recognition, regardless of the angle of observation and interference.

The disadvantages of this method are the computational complexity due to the use of wavelet transform and decomposition of the image in scale, as a result, low speed when using large-scale images, insufficient accuracy in recognizing extended linear objects.

Closest to the proposed method is a method of encoding information about geographic systems by images (Pat. No. 2374689 of the Russian Federation, IPC G06T 17/50 (2006.01), 11/27/2009), the essence of which is as follows: in the additional navigation window of the interface, coding automation programs are issued a document defining a set of coding rules for elements of a geographic system, which has functionalities corresponding to the types of elements of a geographic system and which automate all or part of the operations by code information about the elements of these types, and browser-based navigation tools that allow you to move from the description of the encoding rules of this type of elements of the geographical system to the description of the encoding rules of another type of elements of the geographical system structurally related to elements of this type due to the general standards of systems of this type. After that, a place is selected in the specified document that describes the encoding rules for objects of the type that has the next element of the geographic system, perform the appropriate actions to encode information about this element, then encode information about all or part of the associated system elements, using the automatic transition in the specified document to the encoding rules of neighboring elements using the indicated browser-like navigation tools.

The disadvantages of this method are the implementation complexity associated with the need to store many classifiers and attributes of objects, the probability of an error due to the need for the operator to analyze a large amount of input data.

An object of the invention is to increase the accuracy of highlighting extended linear objects in aerospace images.

The problem is solved in that in the known method for extracting extended linear objects, including obtaining an image in digital form, pre-processing the image and extracting extended linear objects, image processing using the FIR filter is added, which allows you to determine the points belonging to extended linear objects of the image by a positive response FIR filter, based on the analysis of the mean square deviation (RMS) of the brightness of the dots of the filter window when you rotate the filter window relative regarding each studied point of the image.

The invention can be used to highlight objects of the road network when creating electronic maps based on aerospace images for geographic information systems and meets the criterion of "industrial applicability".

The invention is illustrated by drawings, where Fig. 1 shows a diagram of a rotation of a filter window relative to a certain point, Fig. 2 is a drawing explaining a method for determining a brightness value of a point in a discrete space, Fig. 3 and Fig. 4 show a block diagram of a selection algorithm extended linear objects in aerospace images.

The method consists in determining image points with a positive response of the FIR filter, building on their basis an RMS map and a map of the angles of rotation of the FIR filter window, highlighting extended linear objects according to the obtained maps.

For each image point, the value of the response function of the FIR filter is determined. For this, the function of the filter response is studied at a point at each position of the filter window, rotated relative to the point with a step of rotation angle φ.

Figure 1 presents a diagram of the rotation of the filter window relative to the investigated point. The filter window is a straight line segment, the center of which is the studied image point. Thus, the width of the filter window is equal to one point, the length L of the filter window is set.

The angle α of rotation of the filter window is in the range from zero to one hundred eighty degrees, including the value of zero degrees: [0 °; 180 °). The number of possible positions of the filter window is equal to the number of rotation angles N _φ , into which the range [0 °; 180 °). The pitch angle φ is calculated as:

$$\varphi ={180}^{\circ }/N_{\varphi }\left(\mathrm{one}\right)$$

We study a segment equal to half the length of the filter window, due to the symmetry of the coordinates of the points of the filter window relative to the studied image point (center of coordinates), as well as the need for additional verification to exclude points that do not belong to the desired objects. Therefore, for one position of the filter window, two RMSE values are determined for each half of the filter window. The first half of the filter window includes points with a positive coordinate offset relative to the Y axis; the first RMSE value is calculated for them. The second half of the filter window includes points with a negative coordinate offset relative to the Y axis; for them, the second RMS value is calculated. Calculations are carried out for each position of the filter window determined by the angle of rotation relative to the investigated point. All first RMSE values for all investigated positions of the filter window form the first group of RMSE values, all second RMSE values for all investigated positions of the filter window form the second group of RMSE values.

The mathematical expectation M (1), the variance D (l) and the standard deviation σ (1) of the discrete random variable f (x, y) of the half window of the FIR filter of a certain position 1 are calculated according to the well-known formulas of probability theory [Kochetkov P.A. A short course in probability theory and mathematical statistics / P.A. Kochetkov. Textbook - Moscow, MGIU, 1999]:

$$M\left(l\right)=\frac{{\displaystyle \sum _i^l\sout{\int }\left(x_i,y_i\right)}}{L/2}\left(2\right)$$

where i is a certain point of half the filter window; L is the length of the filter window; f (x _i , y _i )

- the value of the brightness function at a certain point i of the half of the filter window.

$$D\left(l\right)=\frac{{\displaystyle \sum _i^l{\left[\sout{\int }\left(x_i,y_i\right)-M\left(l\right)\right]}^2}}{L/2}\left(3\right)$$

where i is a certain point of half the filter window; M (l) is the mathematical expectation of a discrete random variable f (x, y) half of the filter window of a certain position l; L is the length of the filter window; f (x _i , y _i ) is the value of the brightness function at a certain i-th point of half of the filter window.

$$\sigma \left(l\right)=\sqrt{D\left(l\right)}\left(\mathrm{four}\right)$$

where D (l) is the variance of the random variable f (x, y) of the filter window half of a certain position l of the filter window.

Among the found MSE values of all filter window positions, the minimum σ ₁ for the first group of MSE values and the minimum σ ₂ for the second group of MSE values are determined.

The found values of σ ₁ and σ ₂ must satisfy the condition:

$$\begin{array}{c}{\sigma }_{\mathrm{one}}\left(j_{\mathrm{one}}\right)<T_{\mathrm{one}},\\ {\sigma }_2\left(j_2\right)<T_{\mathrm{one}},\\ {\sigma }_{\mathrm{one}}\left(j_{\mathrm{one}}\right)\to 0\\ {\sigma }_2\left(j_2\right)\to 0\\ |\; |\; |{\sigma }_{\mathrm{one}}\left(j_{\mathrm{one}}\right)-{\sigma }_2\left(j_2\right)|\; |\; |\to 0\\ j_{\mathrm{one}}\approx {\mathrm{<figure-callout\; id="2"\; label="j"\; filenames="00000018.png"\; state="\{\{state\}\}">j</figure-callout>}}_2,\end{array}\left(5\right)$$

where j ₁ , j ₂ are the rotation angles of the filter window corresponding to the position of the filter window for the minimum RMS value σ ₁ and the minimum RMS value σ _2, respectively; T ₁ - threshold value, which is set taking into account the degree of overlap of the selected linear objects by other image objects.

According to condition (5), the values of σ ₁ and σ ₂ must have low slightly different values, not exceed a predetermined threshold T ₁ , the corresponding rotation angles should vary slightly.

An additional threshold check is introduced for the difference between the σ ₁ value and the mean standard deviation for the first group of mean-square deviations and the difference between the value of σ ₂ and the mean standard deviation for the second group of mean-square deviations:

$$\{\begin{array}{c}|\; |\; |{\tilde{\sigma }}_{\mathrm{one}}-{\sigma }_{\mathrm{one}}\left(j_{\mathrm{one}}\right)|\; |\; |>{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="00000018.png"\; state="\{\{state\}\}">T</figure-callout>}}_2,\\ |\; |\; |{\tilde{\sigma }}_2-{\sigma }_2\left(j_2\right)|\; |\; |>{\mathrm{<figure-callout\; id="2"\; label="T"\; filenames="00000018.png"\; state="\{\{state\}\}">T</figure-callout>}}_2,\end{array}\left(6\right)$$

- среднее значение СКО для положительного направления смещения координат точек при всех положениях окна фильтра; $${\tilde{\sigma }}_2$$

- среднее значение СКО для отрицательного направления смещения координат точек при всех положениях окна фильтра; Т₂ - пороговое значение, которое устанавливается для исключения точек объектов больших по площади, удовлетворяющих условию (5).Where $${\tilde{\sigma }}_{\mathrm{one}}$$

- average RMS value for the positive direction of the shift of the coordinates of the points at all positions of the filter window; $${\tilde{\sigma }}_2$$

- average RMS value for the negative direction of the shift of the coordinates of the points at all positions of the filter window; T ₂ is a threshold value that is set to exclude points of objects of large area that satisfy condition (5).

Since the filter window is a straight line segment, it is necessary to determine the position of the filter points in the nodes of the discrete grid. Figure 2 presents a figure explaining a method for determining the brightness of a point in discrete space.

To determine the coordinates of the points of the current position of the filter window, the formula is used, obtained by combining the straight line equation (for the current position of the filter window) and the circle that the filter window describes when rotating relative to the image point under study:

$$\{\begin{array}{c}y=\pm \frac{kh}{{\left(\mathrm{one}+k^2\right)}^{\mathrm{one}/2},}\\ x=\pm \frac{h}{{\left(\mathrm{one}+k^2\right)}^{\mathrm{one}/2},}\end{array}\left(7\right)$$

where h is the sequence number of the filter window point starting from the studied image point with the sequence number zero; k is the angular coefficient of the line.

As a result of calculations, the coordinates of the points of the current position of the filter window can take real values. To determine the brightness at the points of the filter window, it is necessary to make a transition to integer coordinates. By rounding up and down to integer values, four environment points are defined for each point (x _h , y _h ) of the filter window.

For each found j-th environment point with coordinates (x _j , y _j ) there is a distance:

$$\{\begin{array}{c}\Delta x=|\; |\; |x_h-x_j|\; |\; |,\\ \Delta y=|\; |\; |y_h-y_j|\; |\; |;\end{array}\left(8\right)$$

further, based on the degree of proximity of the environment points to the point with real coordinates, the value of the weight coefficient is calculated:

$$\omega \left(x_j;y_j\right)=\left(\mathrm{one}-\Delta x\right)\left(\mathrm{one}-\Delta y\right)\left(9\right)$$

The brightness of the point (x _h , y _h ) is determined as follows:

$${\sout{\int }}^{\prime }\left(x^{\prime },y^{\prime }\right)={\displaystyle \sum _j^x{\displaystyle \sum _j^y\omega \left(x_j,y_j\right)\sout{\int }\left(x_j,y_j\right)}}\left(10\right)$$

In the case when the angle of rotation of the filter window is ninety degrees, it is impossible to determine the angular coefficient of the straight line k, since the tangent of the angle of ninety degrees is not defined. But since the location of the filter window in this case is vertical, the coordinates of the points and the brightness values in them are determined without additional calculations.

The total brightness of the point is determined according to the brightness of the environment points, taking into account their weight coefficients.

The response function of the FIR filter at the point under study gives a positive response if condition (5) and condition (6) are satisfied simultaneously for this point. For points with a positive response of the FIR filter response function, a map of the found minimum SDE values and a map of the rotation angles corresponding to these values are constructed. Based on the data obtained, after additional processing based on combining the found points with a positive response of the FIR filter and suppressing noise, the desired extended linear objects are constructed.

The map of the found minimum values of the standard deviation is described by the formula:

$$Q=F_Q\left(\sout{\int }\left(x,y\right)\right),\left(x,y\right)\in Q,\left(\mathrm{eleven}\right)$$

where Q is the standard deviation, at points with a nonzero response of the FIR filter function; F _Q (.) Is the response function of the FIR filter; f (x, y) is the brightness value of the point with coordinates (x, y).

Thus, the standard deviation map is a map of the original image, on which the calculated minimum standard deviations for the image points defined as belonging to extended linear objects are marked in accordance with the coordinates of these points.

The map of the rotation angles of the filter window is described by the formula:

$$\Psi =F_{\Psi }\left(Q\right)\left(12\right)$$

where ψ is a map of rotation angle values; F _ψ (.) Is the function of the angle of rotation of the filter window.

Thus, the map of rotation angles is a map of the original image on which the values of the rotation angles of the filter window are marked, corresponding to the positions of the filter window with the minimum standard deviation, for image points defined as belonging to extended linear objects, in accordance with the coordinates of these points.

The algorithm for highlighting extended linear objects in aerospace images is presented in the form of a block diagram in figure 3 and figure 4.

In block 301 (figure 3) the image is input and pre-processed.

Image is understood as a matrix of pixel brightness values of a digital image. As source images, visual images of the processed Earth remote sensing data presented in digital format can be used. Further image processing is controlled by the controller.

Image pre-processing is performed automatically. As the image preprocessing, the median filter is applied, which allows noise suppression with preservation of brightness differences (contours) without distorting stepwise and sawtooth image functions. Next, the transition to block 302 (figure 3).

In block 302 (Fig. 3), the initial image processing parameters are input. Based on the image parameters and its brightness characteristics, the user sets the following initial parameters: the value of the filter window length L, the number of rotation angles N _φ , the threshold value T ₁ , the threshold value T ₂ .

The filter window length L is set according to the image scale of the aerospace image. The larger the image scale, the larger the value of L.

The value of the number of rotation angles N _φ affects the amount of time required for the operation of the algorithm and the noise level of the map of the rotation angles of the filter window. The larger the value of N _φ , the more time is required for the algorithm to work and the noisier is the map of the angles of rotation of the filter window, however, too low a value of N _φ can lead to the loss of useful information.

The threshold value T ₁ limits the minimum acceptable value of the standard deviation satisfying the search conditions and is set taking into account the degree of overlap of the selected linear objects with other image objects to reduce errors in the search for points.

The threshold value T _{2 of the} difference between the minimum standard deviation and the average value of standard deviation for all possible positions of the filter window of any direction of displacement is set to exclude points of areal objects (fields, lakes, houses, etc.) that are not related to the desired extended linear objects, but satisfy the specified verification conditions (5).

The value of the parameters T ₁ and T _{2 is} also affected by the value of the overall level of brightness and contrast of the image. Next, the transition to block 303 (figure 3).

At block 303 (FIG. 3), verification of completion of image processing is performed. If the processing is not completed, the next image point to be examined is selected in block 304 (Fig. 3). If the processing is completed, based on the data obtained, the construction of the RMS map and the map of the rotation angles of the filter window in block 306 are performed (Fig. 3).

In block 304 (FIG. 3), the coordinates of the next image point under investigation are determined. Next, the transition to block 305 (figure 3).

In block 305 (Fig. 3), the initial value of the angle of rotation of the filter window is set to zero degrees. Next, the transition to block 401 (figure 4).

In block 401 (Fig. 4), the value of the angle α of rotation of the filter window is checked, which should be in the range [0 °; 180 °). If the value of the angle α does not exceed the values allowed by the specified range, there is a further check of the value of the angle α in block 402 (Fig. 4). If the value of the angle α exceeds the permissible values of the specified range, the transition to block 408 (Fig. 4) determines the minimum standard deviation σ ₁ among the standard deviations of the first group.

In block 402 (FIG. 4), it is checked whether the angle α of rotation of the filter window is equal to ninety degrees. If the angle α is not equal to ninety degrees, there is a transition to block 403 (Fig. 4), where according to formula (10) in the discrete space the brightness values of points with a positive coordinate offset relative to the Y axis are determined. If the angle α is ninety degrees, use the formula ( 10) for calculating the brightness of the points it becomes impossible, since it is impossible to calculate the value of the angular coefficient k, since the angle tangent of ninety degrees is not defined. However, the location of the filter window in this case is vertical, the coordinates of the points and the brightness values in them are determined according to the image. In this case, a transition is made to block 405 (Fig. 4) for calculating the RMS value at the current position of the filter window for the first RMS group.

In block 403 (Fig. 4), for all points with a positive coordinate offset relative to the Y axis, at the current position of the filter window, the brightness values in the discrete space are calculated according to formula (10). All found brightness values of points are saved. Next, the transition to block 404 (figure 4).

In block 404 (Fig. 4), for all points with a negative coordinate offset relative to the Y axis, at the current position of the filter window, the brightness values in the discrete space are calculated according to formula (10). All found brightness values of points are saved. Next, the transition to block 405 (figure 4).

In block 405 (Fig. 4), the RMS value is calculated at the current position of the filter window for the first RMS group. RMSD is calculated by the formula (4). Moreover, if the angle α is equal to ninety degrees for the current position of the filter window, the brightness values of image points with a positive coordinate offset relative to the Y axis corresponding to the position of the filter window are used to calculate the standard deviation. If the angle α is not equal to ninety degrees, the brightness values of the points calculated in block 403 are used to calculate the RMSE (FIG. 4). The found MSE value is stored in the first MSE group. For the found RMS value, the corresponding value of the rotation angle α of the current position of the filter window is stored. Next, the transition to block 406 (figure 4).

In block 406 (Fig. 4), the RMS value is calculated at the current position of the filter window for the second RMS group. RMSD is calculated by the formula (4). Moreover, if the angle α is equal to ninety degrees for the current position of the filter window, the brightness of the image points with a negative coordinate offset relative to the Y axis corresponding to the position of the filter window are used to calculate the standard deviation. If the angle α is not equal to ninety degrees, the brightness values of the points calculated in block 404 are used to calculate the RMSE (FIG. 4). The found MSE value is stored in the second MSE group. For the found RMS value, the corresponding value of the rotation angle α of the current position of the filter window is stored. Next, the transition to block 407 (figure 4).

In block 407 (Fig. 4), the value of the rotation angle α of the filter window increases by a step of the rotation angle φ. The pitch of the rotation angle φ is calculated by the formula (1). Next, the transition to block 401 (Fig. 4), where the new value of the rotation angle α of the filter window is checked.

In block 408 (Fig. 4), the minimum standard deviation σ ₁ is determined among the values of the standard deviation of the first group. The found value of σ _{1 is} preserved. Next, the transition to block 409 (figure 4).

In block 409 (FIG. 4), the angle of rotation σ _{1 of} the filter window corresponding to the minimum RMS value σ ₁ found in block 408 (FIG. 4) among the RMS values of the first group is determined. The value of α _{1 is} saved. Next, the transition to block 410 (figure 4).

In block 410 (Fig. 4), the minimum standard deviation σ ₂ is determined among the values of the standard deviation of the second group. The found value of σ _{2 is} retained. Next, the transition to block 411 (figure 4).

In block 411 (FIG. 4), the angle of rotation α _{2 of} the filter window corresponding to the minimum RMS value σ ₂ found in block 410 (FIG. 4) among the RMS values of the second group is determined. The value of α _{2 is} stored. Next, the transition to block 412 (figure 4).

In block 412 (Fig. 4), the fulfillment of condition (5) is checked for the values of σ ₁ found in block 408 (Fig. 4), in block 410 (Fig. 4), the values of σ ₂ and the corresponding values of the angle α ₁ defined in the block 409 (FIG. 4) and an angle α ₂ defined in block 411 (FIG. 4). If condition (5) is fulfilled, there is a transition to block 413 (Fig. 4), in which the average value of standard deviation σ ₁ _ _cp for the first group of standard deviations is calculated. If condition (5) is not fulfilled, it is considered that the point under investigation does not belong to extended linear objects, there is a transition to block 303 (Fig. 3), in which the image processing is completed.

In block 413 (figure 4), the average RMSE σ ₁ _ _Wed for the first group of RMS is calculated. Next, the transition to block 414 (figure 4).

In block 414 (Fig. 4), the mean value of the standard deviation σ ₂ _ _sr for the second group of standard deviations is calculated. Next, the transition to block 415 (figure 4).

In block 415 (FIG. 4), the fulfillment of condition (6) is checked for the value of σ ₁ _ _Wed calculated in block 413 (FIG. 4), calculated in the block in block 414 (FIG. 4) of the value of σ ₂ _ _Wed calculated in block 408 (FIG. 4) of the value of σ ₁ calculated in block 410 (FIG. 4) of the value of σ ₂ . If condition (6) is fulfilled, there is a transition to block 416 (Fig. 4), in which the smaller of the values of σ ₁ and σ _{2 is} determined. If condition (6) is not fulfilled, it is considered that the point under investigation does not belong to extended linear objects, there is a transition to block 303 (Fig. 3), in which the image processing is completed.

In block 416 (figure 4), the smaller of the values of σ ₁ and σ _{2 is} determined. Next, the transition to block 417 (figure 4).

In block 417 (figure 4), the angle of rotation of the filter window is determined corresponding to the smaller of the values of σ ₁ and σ ₂ . Next, the transition to block 418 (figure 4).

In block 418 (Fig. 4) for the studied point, in accordance with the coordinates of this point in the image, the MSE value found in block 416 (Fig. 4) is saved and the rotation angle value found in block 417 (Fig. 4) for the MSE card and rotation angle maps, respectively. The studied point is considered to belong to an extended linear object, its coordinates are saved. Next, the transition to block 303 (figure 3), in which the completion of image processing is checked.

In block 306 (Fig. 3), a RMS map and a map of the rotation angles of the filter window are constructed. After processing the entire image, on the basis of all the data stored in block 418 (Fig. 4), the RMS map, formula (11) and a map of the rotation angles of the filter window, formula (12) are built. Thus, information on the minimum RMSE value and the filter window position at which this RMS was found for points that, according to the algorithm, are considered to belong to extended linear objects, is stored on these maps. Next, the transition to block 307 (figure 3).

In block 307 (Fig. 3), additional processing of the RMS map and the map of the rotation angles of the filter window, as well as the selection of extended linear objects in the image, takes place. Additional processing of the RMS map and the map of the angles of rotation of the filter window is carried out by combining the found points of extended linear objects, taking into account the value of the minimum RMS of the map RMS and the corresponding position of the filter window of the map of the angles of rotation of the filter window of each point. Based on the results of additional processing of the standard deviation map and the map of rotation angles, a diagram of selected extended linear image objects is constructed. According to this scheme, the selection of extended linear objects in the original image.

This invention improves the accuracy of highlighting extended linear objects in aerospace images due to the developed method based on the use of an FIR filter, which allows you to determine the points belonging to extended linear objects of the image by the positive result of the response function of the FIR filter, based on the analysis of the MSE of the brightness of the filter window dots when you rotate the filter window relative to each point in the image.

Claims (1)

The method for extracting extended linear objects in aerospace images, including digitally obtaining an image, preprocessing an image and extracting extended linear objects, is characterized in that it contains image processing using an FIR filter that allows you to determine the points belonging to extended linear image objects by a positive response FIR filter, based on the analysis of the mean square deviation of the brightness of the dots of the filter window dots when the window f filter relative to each studied point of the image.

Images