Suboptimal threshold estimation for detection of point-like objects in radar images

There are many methods for detection of point-like features in gray-leveled bitmap images. The problem of defining a threshold for acceptance or rejection of the results is usually neglected or left to experts. In this paper, a novel method of estimating suboptimal detection threshold values is proposed. It is based on overlapping the results of two or three different methods parametrized with respective thresholds. The quality functions (of two or three variables), whose global extrema (maximum) approximately correspond to the suboptimal levels of thresholds for the used methods, are defined. This method was applied to a series of the bitmaps generated by a radar sensor and by simulated bitmaps.


Introduction
To alleviate the problem of detection and identification of potentially problematic structures in maritime traffic, the Automatic Identification System (AIS) has been introduced into service 15 years ago. However, not all objects at sea have an AIS transceiver; moreover, some deliberately avoid using it, trying to squeeze through a network of radars. There is a great probability that some of these vehicles have been designed to be difficult for detection and localization in space by any sensor (e.g., visual, radar, sonar, or infrared sensor systems). The speed of these vehicles can be below 5 knots, making them hard to detect even with radars with moving target indication (MTI, Doppler) processing. Merchant ships are equipped with classic maritime radars, which usually do not include any kind of MTI processing. This kind of system is commonly used in maritime traffic control centers. The radar operator has to adjust various parameters: gain, brilliance, anti-clutter sea, anti-clutter rain, and tune (manually or in automatic mode) without any guarantee that these objects will be detected.
Many papers deal with image enhancement methods in the form of clutter-reduction techniques. They can be *Correspondence: dzoran1@gmail.com † Equal contributor 1 Horizon Systems -workshop for HW/SW development, Pirot, Serbia, Knjaza Milosa 2 ulaz 2 stan 27, 18300 Pirot, Serbias Full list of author information is available at the end of the article divided in three main categories: based on image processing, based on statistical signal processing [1,2] and constant false alarm rate (CFAR) systems [3,4], and based on artificial neural networks (ANNs), like in [5,6], which are relatively new. The clutter-reduction methods usually delete objects with small reflections and with high fluctuations. Several classical methods of processing and detecting objects in radar images are presented in [7].
The approach in this paper is an attempt to deal with a raw, nonfiltered radar image in order to do detection of problematic objects. Enhanced detection of small, pointlike objects in gray-leveled bitmaps from radar images is a focus of this research.
Reflection of electromagnetic energy from small objects is very low and has large fluctuations in the intensity and distribution around the object. This phenomenon leads to a situation where small objects are hard to be detected by classical methods. It is likely that such objects would not be noticed or detected (i.e., declared as a plot).
Methods based on image processing techniques are used in this paper. The problem of selecting point-like objects in bitmaps belongs to the class of fundamental problems in computer vision systems. In [8], it is pointed at several categories of approaches to this problem. One category is based on the methods of gradient analysis, i.e., changes in the intensities of the resolution elements (pixels) are analyzed. The second category includes methods based on search procedures based on the template matching.
In a series of papers, Kenney et al. [9,10] performed a detailed comparative analysis of these methods. In [10], axiomatic approach to several algorithms of detection is shown. The original proposal was made to compare the methods for detecting characteristic points (corners). Methods like Harris-Stephens, Forstner, Shi-Tomasi, and Rohr are generalized so that the data on position can be in one, two, or N dimensions, and intensity can be from one to M dimensions.
The research in [11] is an example of using normalized cross-correlation (NCC) methods of identification and registration of images. Brunelli and Poggio [12] compared NCC with other methods of searching for patterns in the image and identified the method of localization of characteristic points in the tasks of human face recognition. An example of recent work that deals with the methods of template matching is the work of Lamberti et al. [13].
One of the practical examples of its implementation is given inside the open computer library sponsored by Intel under the name Open Computer Vision (OPEN-CV). Bouguet [14] in his detailed technical description used the Shi and Tomasi method to detect point-like features and follow up these features by using the method of pyramidal implementation based on the work of Lukas and Kanade [15].
All of the abovementioned methods in some part utilize a threshold associated with a decision to accept the result or to refuse. In many studies in the field of characteristic point detection (corners), the threshold question is ignored or it is taken as a predefined constant.
By doing the original experiments described in the papers [16][17][18][19], two methods are used: normalized crosscorrelation and the Shi-Tomasi method. By experimenting with these two methods, it was observed (by a human operator) that the threshold values should have to be changed periodically for obtaining good detection results. The threshold values depend on atmospheric phenomena, sea conditions, and so on. It has been observed that when the thresholds were fine adjusted manually by the experienced operator, the detections of the first and second methods were close as regards the position.
This paper is an attempt to propose a new method which is based on comparing the results of different detection methods.
The present research deals with three different methods where the resulting detections depend on separate thresholds. The results of all three methods are compared. The overlapping quality functions are defined. The maximum value of those functions correspond to new automatic threshold determination. By doing so, need for human operator intervention is eliminated.
The proposed method is tested by a bitmap produced by a maritime radar located on the coast. This radar is working in amplitude/noncoherent mode. In order to do better testing of this method, simulations are made. A certain number of locations are selected to test this method. Around these positions, the intensity values of neighbor pixels are simulated to look like they originate from highly fluctuating objects. By changing the statistical parameters of the simulation, it is possible to test the new proposed method.
The paper is organized as follows: Section 2 presents the methods used in this paper. The Harris and Stephens method is introduced, and the threshold expressed in percent is defined. In a similar way, the Shi-Tomasi method and the appropriate threshold are introduced. After a brief introduction, the normalized cross-correlation method, which belongs to a different category than the first two methods, is presented. Similarly to the first two methods, the threshold is defined. Section 3 explains how the new overlapping quality measure is defined based on counting detections resulting from the application of three different methods. The quality measure gives positive numbers, the bigger values corresponding to better overlapping. In order to find the global maximum of the quality function, two approaches are made. One is brutal force, calculating the results for all discrete threshold values and a more sophisticated one using the well-known hill climbing algorithm. Section 4.1 presents the results of this novel method applied to the recorded radar image of a radar surveillance station. To deeply test this procedure, the simulation model described in Section 4.3 is made. The final conclusions are given in Section 5.

The methods
This section gives a brief overview of the three known methods for detection of characteristic points in an image. These are the Harris-Stephens and Shi-Tomasi, as representatives of the auto-correlation method, and the normalized cross-correlation, which belongs to another category of methods.
The example of a radar bitmap image, which is within the focus of this research, is given in Figure 1. This bitmap is a part of the maritime radar screen. The observed range is 18 nautical miles from the coastal radar station position. The radar screen is north up oriented. The values are 8bit gray leveled and correspond to the power level of the amplitude of the returned signal which is reflected from the space under surveillance. The pixel's X and Y coordinates represent the east and north positive locations from the location of the radar sensor. The sample cell is 75 m and 450 cells to full range. The original bitmap is 900 × 900 pixels, and in this paper, only some parts are shown. The call of each of the methods produces a result that is a matrix of the same dimensions as the matrix representing the starting bitmap. Auto-correlation methods at locations corresponding to the corners give the results that have higher values.
In practice, there is no procedure which defines the threshold, when to accept or reject the result. The threshold is usually selected from case to case relying on the opinion of an expert. There is no clear boundary. The situation is similar when using the normalized crosscorrelation method. Auto-correlation methods produce a similar resulting matrix. The results are invariant with respect to the changes in light intensity, rotation, translation, affine transformations, and scaling. Featuring normalized cross-correlation is invariant with respect to changes in the intensity of the image, but not to the other image transformations. The introduced template matrix is an attempt to ensure the invariance to the other transformations.
The basic idea exploited in this paper is that the most significant (good) results of all applied methods have best overlap only for certain levels of thresholds. Thresholds for each of the methods and functions of overlapping quality are defined.
If there is a maximum of the quality function, there is the best match of the results of the applied methods.
In other words, the optimal thresholds are obtained for all three methods.

Harris-Stephens auto-correlation method
This auto-correlation method is proposed in [20]. Revisiting the research of Moravec [21], a function designed to detect both edges and corners as a linear combination of determinant and trace of μ squared is proposed. Starting from a bitmap image represented by matrix A(x, y), matrix G(p) is calculated as: observing conditions: Matrix G(p) can be represented in the form: The measure of corner (or edge) quality is calculated according to [20] as: where k is the empirical constant. Threshold T h for the Harris-Stephens auto-correlation method can be defined as: where T H represents the relative threshold expressed in percents. It is required for Equation 6 to find the minimal and maximal values of R(p), labeled as R min and R max , respectively. When threshold is applied to R(p), values greater than T h are taken as detections, and values of R(p) less than T h are assigned a zero value. After application of the threshold matrix, R is the subject of spatial filtering, where R is divided into submatrices and only the local maximums are kept, while all other values are set to zero. The total number of local maximums corresponds to the number of detected point-like objects. From these locations, a set of detected features O H is formed.

Shi-Tomasi auto-correlation method
The second analyzed and modified method for detection of point-like object detection is proposed in [22]. In essence, the Shi-Tomasi method is similar to Harris-Stephens, where, from the bitmap image represented with matrix A(x, y), matrix G(p x , p y ) is calculated by (1). Matrix G(p) can be represented in the same way as in (4).
The main difference with regard to the Harris-Stephens method is the calculation of matrix G eigenvalues (7), followed by the smaller value of λ selection (8), which is assigned to all points of matrix L(p): Standard implementations of this procedure require the sorting of L(p) values in ascending order, followed by selection of the first N j locations, where N j is a previously defined constant. Sorting of a large number of points is a time-consuming operation, and it can significantly affect the time of execution of this method.
As an efficient alternative to using the first N j -sorted elements, a relative threshold T S is suggested. The relative threshold represents a ratio in percents compared to the extremal values of L(p), and the appropriate absolute threshold value is: Values of matrix L(p) smaller than T s are assigned a zero value, and then, spatial filtering is performed on the modified matrix L(p), resulting in its division into submatrices whose local maximum locations correspond to the point-like object locations. From these locations, a set of detected features O S is formed.

Normalized cross-correlation-based method
The third method used for detecting and locating pointlike objects is the search pattern method based on normalized cross-correlation. A matrix of dimensions (2ω x + 1) × (2ω y + 1) is used as a template for this search. Further, with ω x = ω y = 3, the template bitmap is a square 7 × 7 matrix, and it is used as a point-like object template Equation 10.
This type of matrix was chosen because dimensions of the objects that look like points are about six to seven pixels. The selected mean value of gray intensity is 127, between the minimum 0 and maximum 255 of each pixel. Distribution of the pixels corresponds roughly to a circle. This was intended to bridge the fact that normalized cross-correlation is not invariant to rotation, translation, affine transformations, and scaling. Martin and Crowley [23] compared the correlation techniques and emphasized this fact: 127 127 127 0 0  0 127 127 127 127 127 0  0 127 127 127 127 127 0  0 127 127 127 127 127 0  0 0 127 127 127 0 The initial task is to solve the problem by minimizing the sum of square differences (MSSD) between the image in bitmap form represented by matrix A(x, y) and the template matrix T in the well-known form given by Equation 11. The main goal is to find locations in matrix A(x, y) with minimum ε(p x , p y ): Dimensions of the bitmap A(x, y) are N x and N y , and the region of ε(p x , p y ) values is computed over ω x < p x < (N x − ω x ), ω y < p y < (N y − ω y ). A good match between pattern T and image A(x, y) is achieved at the positions with the smallest value of ε(p x , p y ). Another known form for ε(p x , p y ) computation is given by Equation 12, known as the normalized cross-correlation. This form is independent of intensity fluctuation and is given as: where A(p x , p y ) is the mean window value of size (2ω x + 1, 2ω y + 1) , and centered at (p x , p y ), while T is the mean value of pattern T. To unify symbols with the previous section relation, K (p) = K (p x , p y ) = (p x , p y ) is made.
The outcome of Equation 12 is in the range K (p) ∈ (−1, 1). The maximum negative values of K (p) correspond to a good match with the original pattern, while large positive values implicate matching with the inverted pattern. Threshold T C is introduced in a similar way like in the above two methods, and relative threshold T c is defined as: The assumption is that the best overlapping results of the comparison of an element of O H , O S , and O C lead to approximately optimally selected thresholds. To determine the quality (the goodness of overlap), Euclidean distance between elements of sets O H , O S , and O C is calculated. When the distance is less than some preassigned value, it can be considered that there is a good fit.
In this paper, the approach is chosen whereby the main matrix is divided into uniformly spaced submatrices where detections are counted. The selected dimensions of submatrices are 12 × 12. The number of detections inside each of the submatrices for the respective method is labeled by hm, sm, and cm. If any of these numbers is larger than zero then the number of detected objects hu, su, and cu is incremented for all methods (Equation 14). When the resulting detections come from two methods in every one of the submatrices, the two methods overlap (Equation 15). Equation 16 shows how the resulting detections are counted for all three methods: The number of successful Harris-Stephens vs Shi-Tomasi methods overlapping is denoted by hsu; scu stands for the overlapping Shi-Tomasi vs normalized crosscorrelation (CC); and the Harris-Stephens vs normalized CC successful hits is marked by hcu Equation 15. The number hscu represents successful overlapping of all three methods (Equation 16).
The overlapping quality for each combination of these methods is obtained by Equations 17 to 20 where Q hs , Q sc , Q hc , and Q hsc , respectively, represent the quality measure for overlapping Harris-Stephens vs Shi-Tomasi, Shi-Tomasi vs normalized CC, Harris-Stephens vs normalized CC, and Harris-Stephens vs Shi-Tomasi vs normalized CC.
The intervals for thresholds T H , T S , and T C can be changed within the values for which the number of detection is greater than zero and less than the maximum number of filtered results. For example, if the image has dimension 900 × 900 and spatial filter window is 9 × 9, the thresholds make sense if the number of detections is from 0 to 100 × 100 = 10, 000. In the case that both methods give the maximum number of overlapping detections, the proposed method does not give useful results.
If the two thresholds are varied from the minimum to maximum meaningful values, one can form a surface which may have one or more local peaks, but only one of them is global and represents the solution for selection of optimal thresholds.

Hill climbing optimization
'Hill climbing' is a mathematical optimization technique for the iterative searching of local or global extreme value of functions of one or more variables. These procedures belong to the widely used methods in artificial intelligence systems [24]. For the experiment in this work, a brutal force method is used to locate the global maximum for the case where two variables are considered. The 'Hill climbing' procedure is developed following the recommendations from [25] and [24], for finding the global extreme in the case of the function of three variables. In this way, automatic finding of suboptimal values of thresholds T H , T S , and T C was made. The main idea about this suboptimal solution is to choose random values of thresholds and then to search for the maximum in the neighborhood. This procedure is repeated for a limited number of times. The number of necessary calculations done in this way is greatly reduced compared to the brutal force method. There is a possibility that this algorithm ends up in a local extremum. The hill climbing method belongs to a special area of research and, in these experiments, is used as is.

The experiments
The idea about the proposed procedure was borne through the solving of the tasks of detecting objects at sea by radar which produces bitmap images. On these bitmapbased images, small objects are shown like small points. Implementing the method proposed by Shi/Tomasi has given satisfactory results, but parameter T S in this method should be determined with the help of experts in interpreting radar image. As an alternative, the normalized cross-correlation method could be implemented in which experts also have to adjust the threshold T C . The same situation is with the Harris-Stephens method. By one radar system, the first two methods were implemented to work in parallel. In the first experiment thresholds of the first and second methods according to the procedures described in Sections 2.1 and 2.2 were left to the operator by choice. When the thresholds were varied, it was observed that for certain values of the thresholds the results of both methods came to a good overlap. For small values of the thresholds, it has been observed that the number of detections increases, but the new detections did not overlap well. A large number of newly detected locations cannot overlap -they diverge. In other words, besides overlapped detections, the number of those which cannot overlap also increases.
As a measure of goodness of the overlap, the criteria In Equations 17 to 20 are defined. In order to test this method, an attempt was made to develop a simulation model. In this newly created test model, it is possible to set up new point-like objects with different statistical parameters which can be varied in a controlled way. By doing this, the whole newly proposed model can be deeply tested.

Real marine radar image
The proposed procedure was tested by the radar images with locations of detections as shown in Figure 1. The same figure with marked regions is shown in Figure 2. Two regions are selected in Figure 2. A large circle at the right of the image marks the region with well-visible objects. These objects may be detected by classical methods like sliding window. A large rectangle covers the area where one small boat shows reflections with high fluctuation. Figure 3 shows six successive bitmap images (successive radar scans) inside the area of the big rectangle.
The brutal force method (i.e., calculation of the quality functions for all values of the thresholds) to find the global extremes of Equations 17, 18, and 19 over the bitmap shown in Figure 1 gives the results presented by the surf diagrams in Figures 4, 5, and 6. These diagrams clearly   Figure 2 presented as a sequence of the same region registered during a 10-min period.
show the existence of existing global extrema. Figure 7 shows objects detected by all three methods with lowthreshold values. The results of different methods are marked by different colors. The main principle used in this study is that when reducing the threshold, one can observe scattering detections of different methods. Finding the optimal thresholds described in Section 3 results in the image with detections shown in Figure 8. It may be noted that the structure from the center of the rectangle of Figure 2 is detected by all three methods.
To simplify the proposed method, the choice is made to use integer values for the thresholds. In this example, meaningful ranges (where the number of detections is greater than zero and less than the maximum, i.e., one per each submatrix filter) are chosen. Thresholds T H , T S , and T C in this example may have 72, 52, and 88 discrete values, respectively. In order to find the global maximum of Equation 20, it is necessary to combine all detections 329,472 times.
By using the hill climbing method, the suboptimal solution is tested. The computing time has now been drastically reduced, as mentioned in Section 3.1. In this case, random numbers with uniform distribution between the minimal and maximal threshold values are selected for T H , T S , and T C , and the results for Q hsc are compared. If a new value of Q hsc (T H , T S , T C ) is greater than the previous value, the search is continued to find if there is any higher value of Q hsc in the neighborhood. The whole procedure is repeated thousands of times, which decreases possibility of ending up at one of the local extrema. In the experiments made while preparing this paper, typical number of attempts was around 15,000 to obtain an approximate location of a global extremum. Figure 8 shows a situation where the thresholds are nearly optimal. Circles of different colors which symbolize the results obtained by separate methods show a relatively good overlapping.

Discussion of the results
The resulting surf diagrams show the existence of global extrema. The diagram in Figure 4 shows that the extremum is noticed for thresholds T C ≈ 32 i T H ≈ 8.   Figure 6 shows the location of maximum overlapping of quality functions for thresholds T S ≈ 27 − 30 and T C ≈ 30 − 35. From these diagrams, it can be concluded that there is an approximate range of thresholds for each method where the defined quality functions give maximum.
Equation 20 was introduced to calculate the overlapping quality for all three methods. The maximum was found by  using the hill climbing method. Like in the previous brutal force methods, similar threshold values are obtained.
In Figure 8, the resulting detection locations are marked by different colors. Thresholds for all three methods are calculated with the hill climbing method. The object with a big fluctuation located at the upper left part of the image was detected by all three methods. There is a chance for misinterpretation when there is only one object. To overcome this possibility, a series of simulations is made.

Simulation
Simulation is made in order to test the proposed procedure. The mean value and standard deviation are calculated for the bitmap shown in Figure 1. The mean value of 45 and standard deviation of 6 are calculated. Then, 50 new objects are placed over the original bitmap. It is assumed that their intensities and densities have big fluctuations. Simulated objects were generated using a pseudorandom generator with normal distribution supposing that real objects are small, with weak radar reflexive surface and that they have echoes which look similar to additive white Gaussian noise (AWGN) with some mean value and variance.
Over Figure 1, by using the pseudorandom number generator with Gaussian distribution, a group of pixels, of elliptical shape, with bigger axis tangentially directed to the sensor position is generated. Length of the major axis was chosen to be seven pixels and the small axis a length of four pixels. Values of the pixel intensity in the groups with the simulated object are created by varying the mean value with a standard deviation of 50. If there are values produced with the pseudorandom generator with Gaussian distribution that are less than 0, then they are made equal to 0, and if greater than 255, they are limited to 255. The positions of the simulated objects are assigned to the set Os (defined previously). Part of the whole bitmap is shown Figure 8 Detection with suboptimal threshold. The image shows detections over Figure 1, with nearly optimal values of thresholds -maximum overlap. Green circles are detection results when Harris-Stephens method is used; blue circles are detection results when Shi-Tomasi method is used; and red circles are detection results when normalized cross-correlation method is used.

Figure 9
The original image with added simulations. The result of the procedure using optimal threshold selection over Figure 1 with added simulated objects. Green circles are detection results when Harris-Stephens method is used; blue circles are detection results when Shi-Tomasi method is used; and red circles are detection results when normalized cross-correlation method is used. Centers of black circles represent positions of the simulated objects.
in Figure 9. In this figure, the bigger black circles represent locations of the simulated objects, and the smaller colored circles represent detections produced by different methods.
Amplitude variations of the simulated objects are so specific that these locations are sometimes hard to be noticed by the human eye. They are marked with big black circles like in Figure 9. Figure 10 summarizes the results of ten successive tests generating 50 objects and a varying mean value from 40 to 120 with a standard deviation of 50 (it is considered to be a big fluctuation). Values bigger than 120 lead to almost 100% detection. Then point-like objects are possible to be detected by using a constant threshold. The average line shows a tendency for better hitrate as the mean value increases. It can be seen that even near the mean value of the whole picture, it is possible to detect many simulated objects.
When the simulated object intensity mean value is 45, there are about 45% detection hits. It can be noticed Figure 10 Simulation of ten successive tests. The hitrate diagram shows the percentage of detection by applying thresholds determined by the proposed procedure to the image with simulated objects. Standard deviation is 50 and mean value is varied from 40 to 120. that the percentage of hits grows when the mean value increases. In the region around 70, there is a fall of the percentage of hits. It is possible that coastal objects produce this effect, or imprecision of the random generator, or maybe the hill climbing method ended in a region of local extrema.
The results of detection in the real and simulated scenarios confirmed the efficiency of the proposed automatic threshold adjustment.
Simulation of real radar signals is complex. Image processing methods are usually tested by a predefined set of test images. The simulation is made in order to test this method. It shows that the maxima of our quality functions exist even for objects which are hard to be noticed by the human/operator's eye.

Conclusions
The task of selecting a threshold level in a number of methods is done manually and is based on the opinions of individual experts in particular fields of interests. This paper is an attempt to automate the process of decision making on the threshold levels.
The starting point is the assumption that each of the methods for a certain threshold of acceptance of the results gives detections that are valid. By lowering the threshold of acceptance, the method starts to detect details that do not correspond to the actual small objects. It has been observed that the detection of almost insignificant objects by all three methods do not match. By using the proposed procedure and function of the quality of overlap, it was possible to detect this moment, and the results are suboptimal threshold values. This paper combined the detections of the methods of Shi-Tomasi, Harris-Stephens, and normalized cross-correlation. It is believed that, in a similar way, other methods can be combined and that they will produce similar overlap.
In practice, it is proven by the mentioned experiments that in the case of coastal radar stations this approach provides satisfactory results. The proposed procedure does not work if the image has the maximum number of objects obtained by filtering. There is a need for existence of a part of the picture in which the results of all methods diverge. The application of hill climbing optimization used in this method makes the results close to real time. On the other hand, theoretical comparative analysis of these three methods is not a trivial process. This proposed overlapping procedure gave satisfactory results. Testing by simulations of this proposed method by a radar image yielded good results. The eventual implementation of the proposed method in the FPGA environment to work faster is one possibility. This method has to be tested in a real dense navigational area with a lot of small boats equipped with AIS transceivers. This will have true positions obtained by AIS and the measured data with high fluctuations produced by radar. By this test, the proposed model will be subjected to further proofs.