CN109712112B - Aerial photography insulator image positioning method based on local features - Google Patents

Aerial photography insulator image positioning method based on local features Download PDF

Info

Publication number
CN109712112B
CN109712112B CN201811401151.2A CN201811401151A CN109712112B CN 109712112 B CN109712112 B CN 109712112B CN 201811401151 A CN201811401151 A CN 201811401151A CN 109712112 B CN109712112 B CN 109712112B
Authority
CN
China
Prior art keywords
image
point
scale
matrix
points
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811401151.2A
Other languages
Chinese (zh)
Other versions
CN109712112A (en
Inventor
赵俊梅
张利平
李晓
任一峰
白鑫
张灵菲
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
North University of China
Original Assignee
North University of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by North University of China filed Critical North University of China
Priority to CN201811401151.2A priority Critical patent/CN109712112B/en
Publication of CN109712112A publication Critical patent/CN109712112A/en
Application granted granted Critical
Publication of CN109712112B publication Critical patent/CN109712112B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Analysis (AREA)

Abstract

The invention belongs to the technical field of image processing in computer vision, and is applied to the power transmission line inspection in a power system. The positioning method of the aerial photography insulator image based on the local characteristics comprises the steps of collecting an acquired image shot along a power transmission line through an unmanned aerial vehicle, and carrying out graying processing on the acquired and reference insulator images; and (3) extracting the characteristic points of the grayed image on an MATLAB platform by using a Harris-Laplace detection method, extracting the characteristic points of the grayed image processed in the second step on the MATLAB platform by using an SURF algorithm, and performing image registration and insulator positioning on the MATLAB platform by using estimation geometric transformation. The invention fully exerts the advantages of rotation, scale and affine invariance of local features, and the registration and positioning of the local feature points have less number of feature points and reduced calculation amount.

Description

Aerial photography insulator image positioning method based on local features
Technical Field
The invention belongs to the technical field of image processing in computer vision, and is applied to the aspect of power transmission line inspection in a power system.
Background
The power transmission line is a vital component in a power grid system, is used as a main line of the power grid system, plays a decisive role in whether the whole power grid is reliable, long-term, safe and stable in operation, and is directly related to the healthy development of national economy in long-term and effective operation of the power grid system. Along with the implementation of transmission network construction engineering, transmission line sharply increases, and line maintenance work volume sharply increases, and the geographical environment is complicated, the climatic conditions circumstances such as changeable make traditional manual work patrol and examine the operation dangerous more, and efficiency and rate of accuracy are patrolled and examined not high. Compared with the traditional inspection mode, the helicopter inspection and unmanned aerial vehicle inspection mode has the advantages of high efficiency, flexible inspection mode, short image acquisition period, no influence of natural environment and the like, and gradually becomes the mainstream mode of power transmission line inspection. Among the power patrol targets, insulators are important inspection targets that are large in number and easily damaged in a power transmission network. In long-term operation, the insulator is easily damaged under the influence of severe weather such as strong wind, thunderstorm, ice mould, ice coating and the like, the normal operation of a power transmission network of a power system is further influenced, and a large-area power failure accident can be caused seriously. The aircraft inspection mode is to shoot and collect images along the power transmission line through the stable speed and the relatively fixed angle of helicopter or unmanned aerial vehicle, and can pass through various landforms such as mountains, rivers, grasslands, houses, arable land and the like in the shooting process, and also can be under various meteorological conditions (rain, snow, fog and the like), and the applicability and robustness of different backgrounds and noises formed by different landforms and different climatic environments to the detection and fault detection of the power transmission line are great challenges.
At present, insulator image detection methods are based on segmentation methods such as colors, textures and shapes, and the detection and positioning effects of the methods are good for the situations of simple background, high resolution and the like. However, most of the actual aerial insulator backgrounds are very complex, especially for long power transmission lines, the situations are different, and the detection and positioning effects are not good, so that an algorithm of local feature positioning is utilized. The insulator positioning method based on the local characteristic points can keep the characteristics of horizontal rotation, translation, scale change and illumination change of the image unchanged, has certain stability on target deformation, local shielding and noise interference, and is suitable for positioning insulators in complex aerial images. The method combines Harris-Laplace with SURF algorithm, and uses estimation geometric transformation to register and position the aerial insulator.
The Harris corner detection algorithm is based on image gray scale, and corner detection is completed by searching a maximum value point with drastically changed gray scale in a target image. The robustness of the algorithm is enhanced through Taylor series expansion, and the problem of registration of image rotation and a bevel edge, which cannot be well adapted by a Moravec operator, is solved.
The Harris algorithm has high resistance to geometric attacks such as rotation and translation, but has no scale invariance. Generally speaking, error detection can be effectively reduced in a larger scale, real feature points can be extracted, but positioning is not easy to be accurate; the feature points can be accurately positioned at a smaller scale, but the false detection rate is increased. The Schmid schmidt and Mikolajczyk schmidt propose an improvement method aiming at the problems, and propose a Harris-Laplace detection method, based on the principle that the feature points are detected under the larger scale at present, and then the real feature points are accurately positioned under the smaller scale, the real corner points can be more accurately positioned by building a scale space to detect the corner points. The method is to estimate the intensity of the feature by using an M matrix, regard the feature as an angular point after passing a threshold value, and then solve the feature scale. The characteristic points detected by the algorithm can resist common image processing such as compression, filtering and the like, and can keep good stability under the condition of large-scale scaling, and meanwhile, the characteristic scale proportion before and after scaling conversion is equal to the scaling proportion.
SURF feature points have stability to local or global disturbances, have scale invariance and rotation invariance, and are robust to noise, occlusion, affine transformation, illumination, and 3D view angle offset. The algorithm is also based on local features, which also has the advantages of the SIFT algorithm. When the SURF algorithm is used for detecting the characteristic points, a determinant and integral image method approximate to a Hessian matrix is used instead of a Gaussian difference algorithm (DoG algorithm) in the SIFT algorithm, so that the SURF algorithm can effectively reduce the calculation amount and effectively improve the speed. The SURF algorithm is represented by 64-dimensional vectors and Laplacian symbols, so that the registration speed can be further improved, and the registration times can be reduced. Therefore, the SURF algorithm is better than the SIFT algorithm in the aspects of the operation speed and the resolution performance and the registration of the feature points, and has better advantages.
The estimation geometric transformation algorithm is simple, the operation is good, and the robustness is good. The method is based on a random sampling consistency algorithm, the algorithm is an iterative algorithm for estimating parameters of a mathematical model, and the method is mainly characterized in that the accuracy of the parameters of the model is gradually improved along with the increase of the iteration times. The main idea is that most characteristics can be solved to meet the parameters of the mathematical model through a sampling and verification strategy.
Disclosure of Invention
The technical problem to be solved by the invention is as follows: aiming at the complex background of an aerial insulating sub-image, the method makes up the defects of detection and positioning based on color, shape and texture characteristics, fully utilizes a local characteristic detection algorithm to carry out registration and positioning, and improves the positioning accuracy.
The technical scheme adopted by the invention is as follows: the method for positioning aerial photography insulator images based on local features is carried out according to the following steps
Acquiring an acquired image shot along a power transmission line by an unmanned aerial vehicle, and carrying out graying processing on the acquired and reference insulator images;
secondly, extracting the characteristic points of the image after the graying processing on an MATLAB platform by utilizing a Harris-Laplace detection method, comprising the following steps
a1. Firstly, generating a scale space through convolution operation of Gaussian kernel functions with different scale parameters and an original image,
Figure BDA0001876344670000041
wherein L (x, y, sigma) is a scale space,
Figure BDA0001876344670000042
representing a convolution operation, I (x, y) is a grayscale image, g (x, y, σ) is a gaussian function,
Figure BDA0001876344670000043
a2. then, a scale-adaptive autocorrelation matrix M is determined, which is often used for feature extraction and local feature structure description of the image, and M is equal to mu (x, y, sigma) by using a representation method of second moment in Harris corner detection1D) Second moment in multi-scale:
Figure BDA0001876344670000044
in the formula, g (σ)1) The expression scale σ1The (x, y) represents the position of the pixel point in the image, and the (l), (x) represents the image after Gaussian smoothing, and the symbol
Figure BDA0001876344670000045
Represents the convolution, Lx(x,y,σD) And Ly(x,y,σD) Respectively, the use of Gaussian g (σ) for the imagesD) The result, σ, of the function in the x or y direction after smoothing1To an integral scale, σDIs a differential scale, σ1Variable, σ, representing the current dimension of the Harris cornerDA variable representing a change in the differential value in the vicinity of the corner;
a3. calculating a characteristic intensity function R of the image I (x, y) subjected to the graying processing at a pixel point (x, y), wherein R is det (u (x, y, sigma)1D))-k·trace2(u(x,y,σ1D) In the formula, R represents the angular point intensity, det represents a determinant of a matrix, trace represents a trace of the matrix, and k is a constant and is generally taken as (0.04-0.06). The position of the local maximum of the corresponding function R is the candidate feature point pkThe coordinates of (a).
a4. Calculating LoG function, judging each feature point by using LoG function, and looking at candidate feature point pkWhether the LoG function can obtain an extremum or not is as follows:
Figure BDA0001876344670000056
in the formula, Lxx(x,y,σ1) Is a second order differential of Gaussian
Figure BDA0001876344670000051
Convolution with image I at point (x, y), Lyy(x,y,σ1) Is a second order differential of Gaussian
Figure BDA0001876344670000052
Convolution with image I at point (x, y);
a5. verifying candidate feature point p of each scale space by using iterative methodkIf the result obtained by the LoG function operation is the regional extreme value in the whole scale space searching range, if not, the point is abandoned and the searching is continued until the characteristic point p with the maximum characteristic intensity R is foundmaxOtherwise, repeating steps a3 and a 4;
thirdly, extracting the characteristic points of the grayed image I (x, y) processed in the second step by using an SURF algorithm on an MATLAB platform, and comprising the following steps
b1. Detecting extreme points, wherein SURF selects points of extreme values obtained by a LoG function in an image scale space as candidate characteristic points, filtering an original image by using box filters with different sizes to form an image pyramid, performing extreme point detection on each layer in the image pyramid by using a Hessian matrix, and defining the Hessian matrix H (x ', sigma') with the scale of sigma 'at the points (x', y ') at one point (x', y ') in the image I' as follows:
Figure BDA0001876344670000053
where σ' is the scale in the image, Lx'x'(x ', y ', σ ') is a second-order differential of Gaussian
Figure BDA0001876344670000054
Convolution with the image I ' at the point (x ', y '), Lx'y'(x ', y ', σ ') is a second-order differential of Gaussian
Figure BDA0001876344670000055
Convolution with the image I ' at the point (x ', y '), Ly'y'(x ', y ', σ ') is highSecond order differential of the order
Figure BDA0001876344670000061
Convolution with the image I ' at point (x ', y ');
b2. positioning feature points, namely solving an extreme value of a scale image at (x ', y ', sigma ') according to a Hessian matrix, then firstly carrying out non-maximum inhibition in a3 x 3 three-dimensional field of the extreme value point, wherein 8 neighborhood pixels are arranged in the same layer, 9 neighborhood pixels are respectively arranged at a lower layer and an upper layer, and only the extreme value points which are larger or smaller than 26 fields around the previous scale, the next scale and the scale can be used as candidate feature points, and then carrying out interpolation operation in a scale space and an image space by using a Taylor series expansion formula to obtain the position and the scale value of the stable feature point;
b3. determining a principal direction, namely taking a characteristic point as a center, calculating Harr wavelet correspondences of points in a field with the radius of 6S (S is a scale value of the characteristic point) in the horizontal and vertical directions, giving Gaussian weight coefficients to the corresponding values, accumulating the responses in a range of 600 to form a new vector, traversing the whole circular region, and selecting the longest vector direction as the principal direction of the characteristic point;
b4. generating a feature descriptor, rotating coordinate axes to a main direction by taking a feature point as a center, wherein the main direction refers to a horizontal axis or a vertical axis, selecting a square area with the side length of 20S according to the main direction, dividing the window area into 4 multiplied by 4 subregions, calculating the wavelet correspondence in the 5 multiplied by 5S range, and calculating the d corresponding to Harr wavelets in the horizontal and vertical directions relative to the main directionx'And dy'The response value coefficients are given in the same way, then corresponding coefficients and absolute values of each sub-region are added to form a feature vector v (∑ dx ', ∑ dy', ∑ dx '|, Σ | dy' |), S represents a sampling step and is also a scale value of a feature point, so that a4 × (4 × 4) ═ 64-dimensional feature description vector is generated for each feature point, and vector normalization is performed, so that certain robustness is achieved for illumination;
step four, on an MATLAB platform, carrying out image registration by using estimation geometric transformation according to the following steps;
c1. the transformation matrix H is initialized to zero;
c2. setting the parameter count to be 0, starting random sampling, wherein the count is counting;
c3. setting the parameter K as the total random sampling number, if the count is less than K, performing the following operation
Increasing the count, namely, the count is equal to the count + 1;
sampling affine transformation method, randomly selecting 3 alignment points from two image reference images A and B acquisition images;
thirdly, calculating a transformation matrix and a matrix on the basis of the selected registration points
Figure BDA0001876344670000071
Point [ x, y ]]Projected point after affine
Figure BDA0001876344670000072
Fourthly, calculating distance measurement by adopting a random sample consensus (RANSAC) algorithm, wherein the calculation formula is as follows:
Figure BDA0001876344670000073
is a characteristic point of the image A,
Figure BDA0001876344670000074
is a characteristic point of the B-picture,
Figure BDA0001876344670000075
representing mapped points based on the matrix H,
Figure BDA0001876344670000076
representing the distance between the characteristic points, wherein t' is a threshold value, and Num is the number of the characteristic points;
the distance metric is determined and if it is less than the matrix H, the new matrix H' is substituted for H. Then, dynamically updating K, and if the registration point is mapped to the matrix H', exiting the sampling cycle;
c4. a fine transformation matrix H "is calculated if all feature points of the two images a and B are mapped to the matrix H'.
c5. Performing iterative refining
Firstly, displaying all characteristic points which can be mapped to a matrix H' as inline points;
secondly, calculating a conversion matrix H' by using the inline points;
if the distance measurement of the matrix H 'is smaller than the matrix H', the new matrix H 'replaces the previous matrix H', otherwise, the circulation is continued;
fifthly, positioning the insulator, namely after registration based on geometric transformation, displaying SURF characteristic points based on local characteristics between a reference insulator image and an actually acquired insulator image, then positioning the insulator sub-target, firstly obtaining the size of the reference image, converting the size into coordinates, and then converting the coordinate system (x, y) of the reference image into the coordinate system (u, y) of the registration image by using an affine geometric transformation method0,v0) And then, marking a rectangular frame on a target insulator image in the registration image to indicate that the positioning is successful, and segmenting the target insulator. And carrying out operations such as binarization, geometric feature extraction and the like on the segmented image so as to judge whether the image works normally.
The invention has the beneficial effects that: the method gives full play to the advantages of rotation, scale and affine invariance of local features, and if local feature points are registered and positioned, the number of the feature points is small, and the calculated amount is reduced; the registration value of the feature point is sensitive to position change, and the registration accuracy is improved; the extraction of the characteristic points can inhibit noise and better adapt to the conditions of gray scale change, local deformation, interference shielding and the like. The collected insulators are accurately positioned, the state of the insulators is conveniently analyzed, and the safe operation of a power grid is guaranteed.
Drawings
FIG. 1 is a flow chart of the present invention.
Detailed Description
As shown in FIG. 1, the method for positioning aerial insulating sub-images based on local features is carried out according to the following steps
Acquiring an acquired image shot along a power transmission line by an unmanned aerial vehicle, and carrying out graying processing on the acquired and reference insulator images;
step two, extracting characteristic points of the grayed image on an MATLAB platform by using a Harris-Laplace detection method, comprising the following steps
a1. Firstly, generating a scale space through convolution operation of Gaussian kernel functions with different scale parameters and an original image,
Figure BDA0001876344670000091
wherein L (x, y, sigma) is a scale space,
Figure BDA0001876344670000092
representing a convolution operation, I (x, y) is a grayscale image, g (x, y, σ) is a gaussian function,
Figure BDA0001876344670000093
x and y refer to the horizontal and vertical axis coordinates of the grayscale image, and σ denotes the scale in the image.
a2. Then, a scale-adaptive autocorrelation matrix M is determined, which is often used for feature extraction and local feature structure description of the image, and M is equal to mu (x, y, sigma) by using a representation method of second moment in Harris corner detection1D) Second moment in multi-scale:
Figure BDA0001876344670000094
in the formula, g (σ)1) Representation scale σ1The (x, y) represents the position of the pixel point in the image, and the (l), (x) represents the image after Gaussian smoothing, and the symbol
Figure BDA0001876344670000095
Represents the convolution, Lx(x,y,σD) And Ly(x,y,σD) Respectively, the use of Gaussian g (σ) for the imagesD) The result, σ, of the function in the x or y direction after smoothing1To an integral scale, σDIs a differential scale, σ1Variable, σ, representing the current dimension of the Harris cornerDA variable representing a variation in the differential value in the vicinity of the corner;
a3. calculating a characteristic intensity function R of the image I (x, y) subjected to the graying processing at a pixel point (x, y), wherein R is det (u (x, y, sigma)1D))-k·trace2(u(x,y,σ1D) In the formula, R represents angular point intensity, det represents a determinant of a matrix, trace represents a trace of the matrix, and k is a constant and is generally (0.04-0.06). The position of the local maximum of the corresponding function R is the candidate feature point pkThe coordinates of (a).
a4. Calculating LoG function, judging each feature point by using LoG function, and looking at candidate feature point pkWhether the LoG function can obtain an extreme value or not is as follows:
Figure BDA0001876344670000107
in the formula, Lxx(x,y,σ1) Is a second order differential of Gaussian
Figure BDA0001876344670000101
Convolution with image I at point (x, y), Lyy(x,y,σ1) Is a second order differential of Gaussian
Figure BDA0001876344670000102
Convolution with image I at point (x, y);
a5. verifying candidate feature point p of each scale space by using iterative methodkIf the result obtained by the LoG function operation is the regional extreme value in the whole scale space searching range, if not, the point is abandoned and the searching is continued until the characteristic point p with the maximum characteristic intensity R is foundmaxOtherwise, repeating steps a3 and a 4;
thirdly, extracting the characteristic points of the grayed image I (x, y) processed in the second step by using an SURF algorithm on an MATLAB platform, and comprising the following steps
b1. Detecting extreme points, selecting points with extreme values obtained by LoG function in image scale space as candidate characteristic points by SURF, and using different sizesThe block filter performs filtering processing on an original image to form an image pyramid, a Hessian matrix is used for carrying out extreme point detection on each layer in the image pyramid, and the Hessian matrix H (x ', sigma') with the scale sigma 'is defined as follows at a point (x', y ') in the image I':
Figure BDA0001876344670000103
where σ' is the scale in the image, Lx'x'(x ', y ', σ ') is a second-order differential of Gaussian
Figure BDA0001876344670000104
Convolution with the image I ' at point (x ', y '), Lx'y'(x ', y ', σ ') is a second-order differential of Gaussian
Figure BDA0001876344670000105
Convolution with the image I ' at point (x ', y '), Ly'y'(x ', y ', σ ') is a second-order differential of Gaussian
Figure BDA0001876344670000106
Convolution with the image I ' at point (x ', y ');
b2. positioning feature points, namely solving an extreme value of a scale image at (x ', y ', sigma ') according to a Hessian matrix, then firstly carrying out non-maximum inhibition in a3 x 3 three-dimensional field of the extreme value point, wherein 8 neighborhood pixels are arranged in the same layer, 9 neighborhood pixels are respectively arranged at a lower layer and an upper layer, and only the extreme value points which are larger or smaller than 26 fields around the previous scale, the next scale and the scale can be used as candidate feature points, and then carrying out interpolation operation in a scale space and an image space by using a Taylor series expansion formula to obtain the position and the scale value of the stable feature point;
b3. determining a principal direction, namely taking a characteristic point as a center, calculating Harr wavelet correspondences of points in a field with the radius of 6S (S is a scale value of the characteristic point) in the horizontal and vertical directions, giving Gaussian weight coefficients to the corresponding values, accumulating the responses in a range of 600 to form a new vector, traversing the whole circular region, and selecting the longest vector direction as the principal direction of the characteristic point;
b4. generating a feature descriptor, rotating a coordinate axis to a principal direction by taking a feature point as a center, wherein the principal direction refers to a horizontal axis or a vertical axis, selecting a square area with the side length of 20S according to the principal direction, dividing the window area into 4 multiplied by 4 sub-areas, calculating the wavelet correspondence in the 5 multiplied by 5S range, and calculating the d corresponding to Harr wavelets in the horizontal and vertical directions relative to the principal directionx'And dy'The response value coefficients are given in the same way, then corresponding coefficients and absolute values of each sub-region are added to form a feature vector v (∑ dx ', ∑ dy', ∑ dx '|, Σ | dy' |), S represents a sampling step and is also a scale value of a feature point, so that a4 × (4 × 4) ═ 64-dimensional feature description vector is generated for each feature point, and vector normalization is performed, so that certain robustness is achieved for illumination;
fourthly, on an MATLAB platform, carrying out image registration by using estimation geometric transformation according to the following steps;
(1) the transformation matrix H is initialized to zero;
(2) setting the parameter count to be 0, starting random sampling, wherein the count is counting;
(3) setting the parameter K as the total number of random samples, and if count is less than K, performing the following operation
Increasing the count, namely, count is equal to count + 1;
sampling affine transformation method, randomly selecting 3 alignment points from two image reference images A and B;
thirdly, calculating a transformation matrix and a matrix on the basis of the selected registration points
Figure BDA0001876344670000121
Point [ x, y ]]Projected point after affine
Figure BDA0001876344670000122
Fourthly, calculating distance measurement by adopting a random sample consensus (RANSAC) algorithm, wherein the calculation formula is as follows:
Figure BDA0001876344670000123
is a characteristic point of the image A,
Figure BDA0001876344670000124
is a characteristic point of the B-picture,
Figure BDA0001876344670000125
representing the mapped points based on the matrix H,
Figure BDA0001876344670000126
representing the distance between the characteristic points, wherein t' is a threshold value, and Num is the number of the characteristic points;
the distance metric is determined and if it is less than the matrix H, the new matrix H' is substituted for H. Then, dynamically updating K, and if the registration point is mapped to the matrix H', exiting the sampling cycle;
(4) a fine transformation matrix H "is calculated if all feature points of the two images a and B are mapped to the matrix H'.
(5) Performing iterative refining
Firstly, displaying all characteristic points which can be mapped to a matrix H' as inline points;
computing a conversion matrix H' ″ by using an inline point;
if the distance measurement of the matrix H 'is smaller than the matrix H', the new matrix H 'replaces the previous matrix H', otherwise, the circulation is continued;
fifthly, positioning the insulator, namely after registration based on geometric transformation, displaying SURF characteristic points based on local characteristics between a reference insulator image and an actually acquired insulator image, then positioning the insulator sub-target, firstly obtaining the size of the reference image, converting the size into coordinates, and then converting the coordinate system (x, y) of the reference image into the coordinate system (u, y) of the registration image by using an affine geometric transformation method0,v0) And then, marking a rectangular frame on a target insulator image in the registration image to indicate that the positioning is successful, segmenting the target insulator, and performing binarization and geometric feature extraction operations on the segmented image to further judge whether the image works normally.
It is to be noted that the present invention does not relate to insulator state recognition. And ending the positioning.

Claims (1)

1. The aerial photography insulator image positioning method based on the local features is characterized in that: the method comprises the following steps
Acquiring an acquired image shot along a power transmission line by an unmanned aerial vehicle, and carrying out graying processing on the acquired and reference insulator images;
secondly, extracting the characteristic points of the image after the graying processing on an MATLAB platform by utilizing a Harris-Laplace detection method, comprising the following steps
a1. Firstly, generating a scale space through convolution operation of Gaussian kernel functions with different scale parameters and an original image,
Figure FDA0003636382200000011
wherein L (x, y, sigma) is a scale space,
Figure FDA0003636382200000012
representing a convolution operation, I (x, y) is a grayscale image, g (x, y, σ) is a gaussian function,
Figure FDA0003636382200000013
a2. then, a scale-adaptive autocorrelation matrix M is determined, which is often used for feature extraction and local feature structure description of the image, and M is equal to mu (x, y, sigma) by using a representation method of second moment in Harris corner detection1D) Second moment in multi-scale:
Figure FDA0003636382200000014
in the formula, g (σ)1) The expression scale σ1(x, y) represents the position of the pixel point in the image, L (x) represents the image after Gaussian smoothing, and the symbol
Figure FDA0003636382200000015
Represents the convolution, Lx(x,y,σD) And Ly(x,y,σD) Respectively, the use of Gaussian g (σ) for the imagesD) The result of the function smoothing in the x and y directions, σ1To an integral scale, σDIs a differential scale, σ1Variable, σ, representing the current dimension of the Harris cornerDA variable representing a change in the differential value in the vicinity of the corner;
a3. calculating a characteristic intensity function R of the image I (x, y) subjected to the graying processing at a pixel point (x, y), wherein R is det (mu (x, y, sigma)1D))-k·trace2(μ(x,y,σ1D) In the formula, R represents the angular point intensity, det represents a determinant of a matrix, trace represents a trace of the matrix, and k is a constant and has a value of 0.04-0.06;
the position of the local maximum value of the function R is the candidate characteristic point pkThe coordinates of (a);
a4. calculating LoG function, judging each feature point by using LoG function, and looking at candidate feature point pkWhether the LoG function can obtain an extremum or not is as follows:
Figure FDA0003636382200000021
in the formula, Lxx(x,y,σ1) Is a second order differential of Gaussian
Figure FDA0003636382200000022
Convolution with image I at point (x, y), Lyy(x,y,σ1) Is a second order differential of Gaussian
Figure FDA0003636382200000023
Convolution with image I at point (x, y);
a5. verifying candidate feature point p of each scale space by using iterative methodkWhether the result obtained by LoG function operation is a regional extreme value in the whole search range of the scale space or not, and if not, whether the result is a regional extreme valueExtreme value, abandoning the point and continuing searching until finding the characteristic point p with the maximum characteristic intensity RmaxOtherwise, repeating steps a3 and a 4;
thirdly, extracting the characteristic points of the grayed image I (x, y) processed in the second step by using an SURF algorithm on an MATLAB platform, and comprising the following steps
b1. Detecting extreme points, namely selecting points with extreme values obtained by a LoG function in an image scale space as candidate characteristic points by SURF, filtering an original image by using box filters with different sizes to form an image pyramid, detecting the extreme points by using a Hessian matrix for each layer in the image pyramid, and defining the Hessian matrix H (x ', sigma') with the scale of sigma 'at the points (x', y ') in the image I' as follows:
Figure FDA0003636382200000031
where σ' is the scale in the image, Lx'x'(x ', y ', σ ') is a second-order differential of Gaussian
Figure FDA0003636382200000032
Convolution with the image I ' at point (x ', y '), Lx'y'(x ', y ', σ ') is a second-order differential of Gaussian
Figure FDA0003636382200000033
Convolution with the image I ' at point (x ', y '), Ly'y'(x ', y ', σ ') is a second-order differential of Gaussian
Figure FDA0003636382200000034
Convolution with the image I ' at point (x ', y ');
b2. positioning feature points, namely solving an extreme value of a scale image at (x ', y ', sigma ') according to a Hessian matrix, then firstly carrying out non-maximum inhibition in a3 x 3 three-dimensional field of the extreme value point, wherein 8 neighborhood pixels are arranged in the same layer, 9 neighborhood pixels are respectively arranged at a lower layer and an upper layer, and only the extreme value points which are larger or smaller than 26 fields around the previous scale, the next scale and the scale can be used as candidate feature points, and then carrying out interpolation operation in a scale space and an image space by using a Taylor series expansion formula to obtain the position and the scale value of the stable feature point;
b3. determining a principal direction, namely taking a characteristic point as a center, calculating Harr wavelet responses of points in the field with the radius of 6S in the horizontal and vertical directions, wherein S is a scale value of the characteristic point, giving Gaussian weight coefficients to the response values, accumulating the responses within the range of 60 degrees to form a new vector, traversing the whole circular region, and selecting the longest vector direction as the principal direction of the characteristic point;
b4. generating a characteristic descriptor, taking a characteristic point as a center, rotating a coordinate axis to a main direction, wherein the main direction refers to a horizontal axis or a vertical axis, selecting a square area with the side length of 20S according to the main direction, dividing the area into 4 multiplied by 4 sub-areas, calculating the wavelet response in the range of 5 multiplied by 5S, and calculating the d of Harr wavelet response in the horizontal direction and the vertical direction relative to the main directionx'And dy'The response value coefficients are given in the same way, then the response coefficients and the absolute values of the response coefficients of each sub-area are added to form a characteristic vector v (Σ dx ', ∑ dy', ∑ dx '|, Σ | dy' |), S represents a sampling step and is also a scale value of a characteristic point, so that a4 × (4 × 4) ═ 64-dimensional characteristic description vector is generated for each characteristic point, and then vector normalization is performed, so that certain robustness is achieved for illumination;
fourthly, on an MATLAB platform, carrying out image registration by using estimation geometric transformation according to the following steps;
c1. the transformation matrix H is initialized to zero;
c2. setting the parameter count to be 0, starting random sampling, wherein the count is counting;
c3. setting the parameter K as the total number of random samples, and if count is less than K, performing the following operation
Increasing the count, namely, the count is equal to the count + 1;
sampling affine transformation method, randomly selecting 3 alignment points from two image reference images A and B;
calculating transformation moment based on the selected registration pointMatrix, matrix
Figure FDA0003636382200000041
Point [ x, y ]]Projected point after affine
Figure FDA0003636382200000042
Fourthly, calculating distance measurement by adopting a random sample consensus (RANSAC) algorithm, wherein the calculation formula is as follows:
Figure FDA0003636382200000043
Figure FDA0003636382200000044
is a characteristic point of the image A,
Figure FDA0003636382200000045
is a characteristic point of the B-picture,
Figure FDA0003636382200000046
representing mapped points based on the matrix H,
Figure FDA0003636382200000047
representing the distance between the characteristic points, wherein t is a threshold value, and Num is the number of the characteristic points;
judging the distance measurement, and if the distance measurement is smaller than the matrix H, replacing H with a new matrix H'; then, dynamically updating K, and if the registration point is mapped to the matrix H', exiting the sampling cycle;
c4. computing a fine transformation matrix H "if all feature points of the two images a and B are mapped to matrix H';
c5. performing iterative refining
i. Displaying all feature points that can be mapped to the matrix H "as inline points;
ii. Calculating a transformation matrix H' "using the inlining points;
iii, if the distance measure of the matrix H ' "is smaller than the matrix H '", the new matrix H "" replaces the previous matrix H ' ", otherwise the cycle continues;
fifthly, positioning the insulator, namely after registration based on geometric transformation, displaying SURF characteristic points based on local characteristics between a reference insulator image and an actually acquired insulator image, then positioning the insulator sub-target, firstly obtaining the size of the reference image, converting the size into coordinates, and then converting the coordinate system (x, y) of the reference image into the coordinate system (u, y) of the registration image by using an affine geometric transformation method0,v0) And then, marking a rectangular frame on a target insulator image in the registration image to indicate that the positioning is successful, segmenting the target insulator, and performing binarization and geometric feature extraction operations on the segmented image to further judge whether the image works normally.
CN201811401151.2A 2018-11-22 2018-11-22 Aerial photography insulator image positioning method based on local features Active CN109712112B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811401151.2A CN109712112B (en) 2018-11-22 2018-11-22 Aerial photography insulator image positioning method based on local features

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811401151.2A CN109712112B (en) 2018-11-22 2018-11-22 Aerial photography insulator image positioning method based on local features

Publications (2)

Publication Number Publication Date
CN109712112A CN109712112A (en) 2019-05-03
CN109712112B true CN109712112B (en) 2022-06-24

Family

ID=66254933

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811401151.2A Active CN109712112B (en) 2018-11-22 2018-11-22 Aerial photography insulator image positioning method based on local features

Country Status (1)

Country Link
CN (1) CN109712112B (en)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110533635B (en) * 2019-07-30 2023-07-14 江苏理工学院 Soft package surface quality detection method based on machine vision
CN110852146B (en) * 2019-09-23 2023-05-16 合肥赛为智能有限公司 Unmanned aerial vehicle image feature point detection method
CN110764529B (en) * 2019-10-21 2020-07-21 安徽诺乐知识产权服务有限公司 Flight direction correction platform, method and storage medium based on target positioning big data
CN111754465B (en) * 2020-06-04 2023-06-09 四川大学 Insulator positioning and string dropping detection method
CN112013830B (en) * 2020-08-20 2024-01-30 中国电建集团贵州电力设计研究院有限公司 Accurate positioning method for inspection image detection defects of unmanned aerial vehicle of power transmission line
CN113156274B (en) * 2021-01-27 2024-05-28 南京工程学院 Non-contact detection system and method for degraded insulator based on unmanned aerial vehicle
CN115546568B (en) * 2022-12-01 2023-03-10 合肥中科类脑智能技术有限公司 Insulator defect detection method, system, equipment and storage medium

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102629330A (en) * 2012-02-29 2012-08-08 华南理工大学 Rapid and high-precision matching method of depth image and color image
CN102865859A (en) * 2012-09-21 2013-01-09 西北工业大学 Aviation sequence image position estimating method based on SURF (Speeded Up Robust Features)
CN102915540A (en) * 2012-10-10 2013-02-06 南京大学 Image matching method based on improved Harris-Laplace and scale invariant feature transform (SIFT) descriptor
CN103442180A (en) * 2013-08-27 2013-12-11 桂林电子科技大学 Binocular video splicing device based on SOPC and binocular video splicing method
CN104764748A (en) * 2015-05-04 2015-07-08 成都唐源电气有限责任公司 Method and system for positioning insulators and method and system for fault detection
CN107689058A (en) * 2017-09-01 2018-02-13 哈尔滨理工大学 A kind of image registration algorithm based on SURF feature extractions
CN108426566A (en) * 2018-02-28 2018-08-21 中国计量大学 A kind of method for positioning mobile robot based on multiple-camera

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8358840B2 (en) * 2007-07-16 2013-01-22 Alexander Bronstein Methods and systems for representation and matching of video content
CN103063166B (en) * 2013-01-05 2015-03-18 山西省电力公司大同供电分公司 Detection method and device for wind deflection angle of suspension type composite insulator chain
CN106940876A (en) * 2017-02-21 2017-07-11 华东师范大学 A kind of quick unmanned plane merging algorithm for images based on SURF

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102629330A (en) * 2012-02-29 2012-08-08 华南理工大学 Rapid and high-precision matching method of depth image and color image
CN102865859A (en) * 2012-09-21 2013-01-09 西北工业大学 Aviation sequence image position estimating method based on SURF (Speeded Up Robust Features)
CN102915540A (en) * 2012-10-10 2013-02-06 南京大学 Image matching method based on improved Harris-Laplace and scale invariant feature transform (SIFT) descriptor
CN103442180A (en) * 2013-08-27 2013-12-11 桂林电子科技大学 Binocular video splicing device based on SOPC and binocular video splicing method
CN104764748A (en) * 2015-05-04 2015-07-08 成都唐源电气有限责任公司 Method and system for positioning insulators and method and system for fault detection
CN107689058A (en) * 2017-09-01 2018-02-13 哈尔滨理工大学 A kind of image registration algorithm based on SURF feature extractions
CN108426566A (en) * 2018-02-28 2018-08-21 中国计量大学 A kind of method for positioning mobile robot based on multiple-camera

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Algorithm of Remote Sensing Image Matching Based on Corner-Point;Wang Changjie等;《2017 International Workshop on Remote Sensing with Intelligent Processing (RSIP)》;20170626;第1-4页 *
基于改进SURF 算法的人脸点云配准;郭昱等;《光学技术》;20180515;第44卷(第3期);第333-338页 *

Also Published As

Publication number Publication date
CN109712112A (en) 2019-05-03

Similar Documents

Publication Publication Date Title
CN109712112B (en) Aerial photography insulator image positioning method based on local features
CN109325935B (en) Power transmission line detection method based on unmanned aerial vehicle image
CN109949340A (en) Target scale adaptive tracking method based on OpenCV
CN109785379A (en) The measurement method and measuring system of a kind of symmetric objects size and weight
CN111428748A (en) Infrared image insulator recognition and detection method based on HOG characteristics and SVM
Li et al. Adaptive building edge detection by combining LiDAR data and aerial images
CN103697855B (en) A kind of hull horizontal attitude measuring method detected based on sea horizon
CN103295232B (en) Based on the SAR image registration method in straight line and region
CN111007531A (en) Road edge detection method based on laser point cloud data
CN104715474B (en) High resolution synthetic aperture radar linearity building object detecting method based on Based On Method of Labeling Watershed Algorithm
CN110428425B (en) Sea-land separation method of SAR image based on coastline vector data
CN103761731A (en) Small infrared aerial target detection method based on non-downsampling contourlet transformation
CN109523528B (en) Power transmission line extraction method based on unmanned aerial vehicle binocular vision SGC algorithm
CN109410248B (en) Flotation froth motion characteristic extraction method based on r-K algorithm
CN108257155B (en) Extended target stable tracking point extraction method based on local and global coupling
CN110245566B (en) Infrared target remote tracking method based on background features
CN110222661B (en) Feature extraction method for moving target identification and tracking
CN106530313A (en) Sea-sky line real-time detection method based on region segmentation
CN106156758B (en) A kind of tidal saltmarsh method in SAR seashore image
CN111242000A (en) Road edge detection method combining laser point cloud steering
CN103679740B (en) ROI (Region of Interest) extraction method of ground target of unmanned aerial vehicle
CN106886988A (en) A kind of linear goal detection method and system based on unmanned aerial vehicle remote sensing
CN109948629B (en) GIS equipment X-ray image fault detection method based on SIFT features
Tong et al. Transmission line extraction and recognition from natural complex background
CN109544609A (en) A kind of sidescan-sonar image matching process based on SIFT algorithm

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant