Image Feature Extraction: Spot Detection

1. What is a spot?

A spot usually refers to a area with a different color and gray level from the surrounding area. In actual maps, there are often many such spots. For example, a tree is a spot, a lawn is a spot, and a house can also be a spot. Because a spot represents a region, it has better stability and better anti-noise ability than a simple corner, so it plays a very important role in image registration.

At the same time, sometimes the spots in images are also areas of our concern, such as in the medical and biological fields, we need to extract the location or number of spot with special significance from some X-ray photos or cell microscopic photos.

For example, the plane in the middle sky, the flower disk of the sunflower, and two spots in the X-ray tomography.

In the field of vision, the main idea of spot detection is to detect areas in the image that are larger or smaller than the pixel gray value around it. There are two methods to achieve this:

1. A derivative method is called a differential detector;
2. Watershed Algorithm Based on Local extreme values.

Here we will focus on the first method, mainly to detect log spots. The simpleblobdetector dot detection operator in opencv implements the second method. Here we will also introduce its interface usage.

2. Basic principles of log spot detection 2.1

Using the Laplace of Gaussian (log) operator to detect image spots is a very common method for two-dimensional Gaussian Functions:

$ G (x, y; \ sigma) =\ frac {1} {2 \ pi \ Sigma ^ 2} exp (-\ frac {x ^ 2 + y ^ 2} {2 \ Sigma ^ 2}) $

Its Laplace transformation is:

$ \ Delta ^ 2G = \ frac {\ partial ^ 2g} {\ partial x ^ 2} + \ frac {\ partial ^ 2g} {\ partial y ^ 2} $

The normalized Gaussian Rapp is transformed:

$ \ Delta ^ 2 _ {norm} = \ Sigma ^ 2 \ Delta ^ 2G = \ Sigma ^ 2 (\ frac {\ partial ^ 2g} {\ partial x ^ 2} + \ frac {\ partial ^ 2g} {\ partial y ^ 2 }) =-\ frac {1} {2 \ pi \ Sigma ^ 2} [1-\ frac {x ^ 2 + y ^ 2} {\ Sigma ^ 2}] \ cdot exp (-\ frac {x ^ 2 + y ^ 2} {2 \ Sigma ^ 2 }) $

The normalization algorithm sub-shows a circular symmetry function on two-dimensional images, as shown in. We can use this operator to detect the spots in the image, and by changing the value of $ \ Sigma $, we can detect two-dimensional spots of different sizes.

2.2 log principles

In fact, it is more intuitive to explain why the LOG operator can detect the spots in the image:

Convolution between an image and a two-dimensional function is actually to obtain the similarity between the image and this function. Likewise, the convolution of the image and Gaussian Laplace function is actually to obtain the similarity between the image and Gaussian Laplace function. When the spot size in the image approaches the same shape as the Gaussian Laplace function, the Laplace response of the image reaches the maximum.

From the perspective of probability, assume that the source image is like a density function of location-related random variable X, while log is the density function of random variable Y, then the density distribution function of the random variable X + Y is the convolution form of the two functions (for the theory of this part, refer to this blog article on Probability and Statistics ). If you want X + Y to reach the maximum value, it is best to keep x consistent with Y, that is, when x rises, y also rises, and when X is the largest, Y is also the largest.

So how was the LOG Operator conceived?

In fact, we know that Laplace can be used to detect local extreme points in an image, but it is sensitive to noise. So before we perform Laplace convolution on the image, we use a Gaussian low-pass filter to perform convolution on the image, the goal is to remove noise points in the image. This process can be described as follows:

First, the image $ f (x, y) $ is filtered by Gaussian Kernel with the variance of $ \ Sigma $ to remove noise in the image.

$ L (X, Y; \ sigma) = f (x, y) * g (x, y; \ sigma) $

Then the Laplace image of the image is:

$ \ Delta ^ 2 = \ frac {\ partial ^ 2L} {\ partial x ^ 2} + \ frac {\ partial ^ 2L} {\ partial y ^ 2} $

In fact, there are the following equations:

$ \ Delta ^ 2 [g (x, y) * f (x, y)] = \ Delta ^ 2 [g (x, y)] * f (x, y) $

Therefore, we can first obtain the Laplace operator of the Gaussian Kernel and then perform convolution on the image. That is, the steps described at the beginning.

2.3 Implementation of log Operators

Mat feat: gethogkernel (size & ksize, double sigma)
{
Mat kernel (ksize, cv_64f );
Point centpoint = point (ksize. Width-1)/2, (ksize. Height-1)/2 ));
// First calculate Gaussian
For (INT I = 0; I <kernel. Rows; I ++)
{
Double * pdata = kernel. PTR <double> (I );
For (Int J = 0; j <kernel. Cols; j ++)
{
Double Param =-(I-centpoint. y) * (I-centpoint. y) + (J-centpoint. x) * (J-centpoint. x)/(2 * Sigma * sigma );
Pdata [J] = exp (PARAM );
}
}
Double maxvalue;
Minmaxloc (kernel, null, & maxvalue );
For (INT I = 0; I <kernel. Rows; I ++)
{
Double * pdata = kernel. PTR <double> (I );
For (Int J = 0; j <kernel. Cols; j ++)
{
If (pdata [J] <EPS * maxvalue)
{
Pdata [J] = 0;
}
}
}

Double sumkernel = sum (kernel) [0];
If (sumkernel! = 0)
{
Kernel = kernel/sumkernel;
}
// Now calculate Laplacian
For (INT I = 0; I <kernel. Rows; I ++)
{
Double * pdata = kernel. PTR <double> (I );
For (Int J = 0; j <kernel. Cols; j ++)
{
Double Addition = (I-centpoint. y) * (I-centpoint. y) + (J-centpoint. x) * (J-centpoint. x)-2 * Sigma * sigma)/(Sigma * sigma );
Pdata [J] * = addition;
}
}
// Make the filter sum to zero
Sumkernel = sum (kernel) [0];
Kernel-= (sumkernel/(ksize. Width * ksize. Height ));

Return kernel;
}

2.4 Multi-scale detection

We noticed that when $ \ Sigma $ has a certain scale, only the corresponding radius of the spot can be detected. What is the radius of the spot detected, we can evaluate and export the normalized two-dimensional Laplace Gaussian OPERATOR:

The normalized Gaussian Laplace function is:

$ \ Delta ^ 2 _ {norm} =-\ frac {1} {2 \ pi \ Sigma ^ 2} [1-\ frac {x ^ 2 + y ^ 2} {\ Sigma ^ 2}] \ cdot exp (-\ frac {x ^ 2 + y ^ 2} {2 \ Sigma ^ 2 }) $

Finding the extreme value of $ \ Delta ^ 2 _ {norm} $ is equivalent to obtaining the following formula:

$ \ Frac {\ partial (\ Delta ^ 2 _ {norm}) }{\ partial \ Sigma} = 0 $

Get:

$ (X ^ 2 + y ^ 2-2 \ Sigma ^ 2) \ cdot exp (-\ frac {(x ^ 2 + y ^ 2 )} {2 \ Sigma ^ 2}) $

$ R ^ 2-2 \ Sigma ^ 2 = 0 $

For spots in the image, the Gaussian Laplace response value reaches the maximum when the scale $ \ Sigma = r/\ SQRT {2} $. Similarly, if the circular spot in the image is black and white, the Gaussian Laplace response value of the circular spot is minimized at $ \ Sigma = r/\ SQRT {2} $. The scale $ \ Sigma $ value when the Gaussian Laplace response reaches the peak value is called the feature scale.

In the case of multiple scales, the points that reach the maximum (or minimum) in space and scale are what we expect. For two-dimensional images $ I (x, y) $, calculate the discrete Laplace response value of the image at different scales, and then check each point in the location space; if the Laplace response value of this point is smaller than or less than the value of the other 26 cubic space fields (9 + 8 + 9), then this point is the detected image spot.

3. opencv for Spot Detection

The class for detecting blobs in opencv is simpleblobdetector. The definition of this class in opencv is as follows:

class SimpleBlobDetector : public FeatureDetector{public:struct Params{    Params();    float thresholdStep;    float minThreshold;    float maxThreshold;    size_t minRepeatability;    float minDistBetweenBlobs;    bool filterByColor;    uchar blobColor;    bool filterByArea;    float minArea, maxArea;    bool filterByCircularity;    float minCircularity, maxCircularity;    bool filterByInertia;    float minInertiaRatio, maxInertiaRatio;    bool filterByConvexity;    float minConvexity, maxConvexity;};SimpleBlobDetector(const SimpleBlobDetector::Params &parameters = SimpleBlobDetector::Params());protected:    ...};

The algorithm steps are as follows:

1. For the [minthreshold, maxthreshold) interval, take thresholdstep as the interval and perform multiple binarization.
2. For each binary image, use findcontours () to extract the Connected Domain and calculate the center of each connected domain.
3. Put all the centers obtained based on 2 together. Some very close points [How much is controlled by themindistbetweenblobs to be close to] are classified as a group, corresponding to a bolb feature ..
4. Estimate the final blob features and corresponding radius from those points obtained from 3 and return them with key points.

At the same time, this method supports feature extraction. There are a total of five options, which are not described here. By default, the Blob feature of the black circle is extracted. The following is an example

Int main (INT argc, char ** argv)
{
Mat image = imread (argv [1]);
Vector <keypoint> keypoints;
Simpleblobdetector: Params;

Simpleblobdetector blobdetect (Params );
Blobdetect. Create ("simpleblob ");
Blobdetect. Detect (image, keypoints );
Cout <keypoints. Size () <Endl;
Drawkeypoints (image, keypoints, image, scalar (255, 0, 0 ));

Namedwindow ("blobs ");
Imshow ("blobs", image );
Waitkey ();
Return 0;
}

In general, opencv's spot detection effect is not good, but it is not as effective as LOG Operator detection in some images.

4. Extended reading

A kernel similar to the LOG filter core is a Gaussian differential dog filter core, which is defined:

$ D (x, y, \ sigma) = (g (x, y, k \ sigma)-g (x, y, \ sigma) * I (X, y) = L (X, Y, k \ sigma)-L (X, Y, \ sigma) $

$ K $ is the proportional factor between two adjacent scales.

Dog can be seen as an approximation of log, but it is more efficient than log.

The Differential Operators described above have good effects in the near-circle spot detection, but these detection operators are limited to only detecting circular spots and cannot estimate the direction of the spots, because log operators are symmetric in the center. If we define a deformation of a two-dimensional Gaussian Kernel and note that it has different variance in the X and Y directions, then this operator can be used to detect vertices with directions.

$ G (x, y) = \ mathcal {A} \ cdot exp (-[A (ax ^ 2 + 2bxy + cy ^ 2)]) $

$ A =\frac {cos ^ 2 \ Theta} {2 \ Sigma ^ 2_x} + \ frac {sin ^ 2 \ Theta} {2 \ Sigma ^ 2_y} $

$ B =-\ frac {sin2 \ Theta} {2 \ Sigma ^ 2_x} + \ frac {sin2 \ Theta} {4 \ Sigma ^ 2_y} $

$ C =\ frac {sin ^ 2 \ Theta} {2 \ Sigma ^ 2_x} + \ frac {cos ^ 2 \ Theta} {2 \ sigma_y ^ 2} $

$ \ Mathcal {A} $ is a compliance factor.

5. References

1. Modern Digital Image-detailed explanation of Processing Technology Improvement and Application Cases

2. Local immutable features and descriptions of images

3. Lindeberg, T. Feature Detection with automatic scale Selection

4. Hui Kong. A generalized Laplacian of Gaussian filter for BLOB detection and its applications.

5. opencv2 marathon 20th lap -- principles and implementation of Blob Feature Detection

Image Feature Extraction: Spot Detection