CN105405146A - Feature density clustering and normal distribution transformation based side-scan sonar registration method - Google Patents

Abstract

The invention relates to the technical field of image registration, and particularly relates to a feature density clustering and normal distribution transformation based side-scan sonar registration method. The method mainly solves the problem that current method can not register sidescan sonar image rapidly and effectively. The method mainly comprises the following steps: firstly, using a mean value filter to carry out filtering on a to-be-registered sidescan sonar image; calculating image gradients based on the filtering; after thresholding a gradient image and getting all feature points, using a density-based spatial clustering of applications with noise (DBSCAN) to carry out clustering on pixel points of the gradient image, so as to obtain image feature points; and finally, using a normal distribution transformation (NDT) algorithm to optimize feature points, so as to obtain a registration relationship of two images. According to the side-scan sonar registration method, the problem that fewer image feature points cause registration failure in the feature-based method is overcome. Registration failure caused by large rotation angle of the two images is prevented. The method can be widely applied to fields of seabed topography and landform probe, AUV autonomous navigation, image processing and so on.

Description

The side scan sonar method for registering of feature based Density Clustering and normal distribution transform

Technical field

The present invention relates to image registration techniques field, particularly the side scan sonar method for registering of feature based Density Clustering and normal distribution transform.

Background technology

Image registration is the technology be transformed at the image that diverse location obtains by sensor under the same coordinate system.The Information Availability obtained by image registration, in submarine navigation device perception surrounding environment, effectively can improve the autonomous perception of submarine navigation device.Underwater environment is because its low visibility, unglazed photograph, meeting are to features such as underwater lighting system generation scatterings, and traditional optical sensor is used for underwater picture acquisition and there is many challenges.Because side scan sonar does not affect by these optical considerations, use it to carry out underwater picture acquisition and more and more come into one's own, side scan sonar registration technology is also corresponding to grow up.

At present, the method for sonar image registration mainly contains two kinds: the method for feature based and the method based on region.The method of feature based is used for the situation that the more and image of unique point in image is successive frame, such as Harris Corner Detection Algorithm (NegahdaripourS, Onprocessingandregistrationofforward-scanacousticvideoim agery, ComputerandRobotVision, 2005) and SIFT (ScaleInvariantFeatureTransform) in conjunction with the algorithm (NegahdaripourS of RANSAC (RANdomSAmpleConsensus), Dynamicsceneanalysisandmosaicingofbenthichabitatsbyfsson arimaging-issuesandcomplexities, InProceedingsofOCEANS, 2011).But these methods feature count less or that the match is successful characteristic number less time registration result poor.In order to the method that the solves feature based Problems existing when registration sonar image, the method for registering based on region is used to sonar image registration.Based on the method for registering in region using entire image as handling object, as phase correlation method.But this method is only applicable to the little situation of rotation angle.

AUV (AutonomousUnderwaterVehicle) carries side scan sonar and extracts ambient condition information, overall topography and geomorphology is under water obtained by registration sonar image, improve AUV to the perception of environment, and how fast and effeciently registration sidescan-sonar image is theme urgently to be resolved hurrily at present.

Summary of the invention

The present invention is directed to deficiency of the prior art, propose the side scan sonar method for registering of a kind of feature based Density Clustering and normal distribution transform.First mean filter is used to carry out filtering to sidescan-sonar image subject to registration, the gradient of computed image on this basis, and use density-based spatial clustering method DBSCAN (Density-BasedSpatialClusteringofApplicationswithNoise) to carry out cluster to the pixel in gradient image to obtain image characteristic point, finally use normal distribution transform NDT algorithm (NormalDistributedTransform) to optimize unique point, obtain the registration relation of two width images.This technical method can be applicable to the fields such as seafloor topography detection, AUV independent navigation and image procossing.

The side scan sonar method for registering of feature based Density Clustering and normal distribution transform specifically comprises the following steps:

Step one, the original sidescan-sonar image of median filter process.

There is noise in original sidescan-sonar image, median filter first will be used sonar image denoising.

Step 2, compute gradient magnitude image.

Sobel operator (Sobeloperator, Sobel Operator) is used to calculate gradient magnitude image.

Step 3, thresholding gradient image obtains whole unique point.

Step 4, uses DBSCAN cluster to obtain image characteristic point.

DBSCAN is a kind of density-based spatial clustering algorithm, this algorithm by the Region dividing with sufficient density for bunch, and find in the noisy spatial database of tool arbitrary shape bunch, it bunch will to be defined as the maximum set of the point that density is connected.The object of DBSCAN algorithm is to filter density regions, finds consistency sample point.

Step 5, uses NDT algorithm (NormalDistributedTransform, normal distribution transform) to obtain the registration result of two width images.

Two dimensional character point is converted to the probability distribution that zonal cooling can be micro-by NDT algorithm.When having the image of cluster feature point to carry out NDT process to a width, first image is divided into the grating image of formed objects, to the grid cell comprising unique point, carrying out zonal cooling with probability density form to it can be micro-, and represents the probability distribution of unique point in each grid cell with normal distribution:

$p\left(x\right)\∝\exp (-\frac{{(x-q)}^T{\mathrm{\Σ}}^{-1}(x-q)}{2})$

Wherein, x is the two-dimensional coordinate of i point in a grid cell, and q is the mean value of the characteristic point position comprised in grid cell, namely $q=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{x}_i;$ Σ is the covariance matrix of these points, $\mathrm{\Σ}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}(x_i-q){(x_i-q)}^T.$

The spatial transform relation T of two width images can be expressed as:

$T:\left(\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right)=\left(\begin{array}{cc}cos\mathrm{\θ}& -sin\mathrm{\θ}\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$

Wherein, $\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$ Be the translation relation of two width images, θ is the rotation relationship of two width images.

Wherein, the concrete cluster process of step 4 is as follows:

A. scan whole gradient image, find any one core point, this core point is expanded.The method expanded finds the data point be connected from all density of this core point.

B. travel through all core points in this core point neighborhood, find the point be connected with these pixel point density, until there is no the pixel that can expand.Finally be clustered into bunch boundary node be all non-core data point.

C. after be exactly rescan gradient image (search out before not comprising bunch in any pixel), find not by the core point of cluster, repeat step above again, this core point is expanded until data centralization does not have new core point.The data point that data centralization is not included in any bunch just forms abnormity point.Last data set is exactly the image characteristic point obtained after cluster.

In step 5, the object of registration is exactly these conversion parameters in order to find between two width sonar images.Specifically comprise the steps:

A. the NDT of a width sonar image is calculated;

B. initialization translation parameters and rotation parameter (use 0 or odometer data carry out initialization);

C. for another width sonar image, transformed in the coordinate system of piece image according to initiation parameter;

D. the normal distribution of each point in sonar image after calculating coordinate change;

E. the normal distribution sum by evaluating each point calculates the mark of translation parameters and rotation parameter;

F. use Hessian matrix method to be optimized these marks, calculate new estimates of parameters;

G. step c is returned, until meet convergent requirement.

Parameter $\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$ Can be expressed as with the mark of θ:

$score\left(p\right)=\underset{i}{\Σ}\exp (-\frac{{(x_i^{\′}-q_i)}^T{\mathrm{\Σ}}_i^{-1}(x_i^{\′}-q_i)}{2})$

Wherein, p=(t _x, t _y, θ) ^t, x _i' be by the some x in the second width image according to coordinate conversion parameter _itransform to coordinate in the first width image coordinate system.Q _iwith x ' respectively _iaverage and covariance.Be optimized mark in above-mentioned e step, make mark score (p) maximum exactly, namely-score (p) is minimum.The detailed operation of Newton iterative is as follows:

In order to make function f minimum, to solve an equation in each iteration:

HΔp＝-g

The Δ p calculated is added with parameter current, obtains new parameter

p←p+Δp

P 0 value initialization, so p=Δ p now.

G is the transposition gradient of f:

$g_i=\frac{\∂f}{\∂p_i}$

H is the Hessian matrix of f:

$H_{ij}=\frac{\∂f}{\∂p_i\∂p_j}$

Definition: q=x ' _i-q _i, then one of-score (p) can be written as:

$s=-\exp \frac{-q^T{\mathrm{\Σ}}^{-1}q}{2}$

The then Grad of s for:

$\begin{array}{c}{\tilde{g}}_i=-\frac{\∂s}{\∂p_i}=-\frac{\∂s}{\∂p}\frac{\∂q}{\∂p_i}\\ =q^T{\Σ}^{-1}\frac{\∂q}{\∂p_i}\exp \frac{-q^T{\Σ}^{-1}q}{2}\end{array}$

Each in Hessian matrix can be expressed as:

$\begin{array}{c}{\tilde{H}}_{ij}=-\frac{\∂s}{\∂p_j\∂p_j}=-\exp \frac{-q^T{\Σ}^{-1}q}{2}\\ \left(\right(-q^T{\Σ}^{-1}\frac{\∂q}{\∂p_i}\left)\right(-q^T{\Σ}^{-1}\frac{\∂q}{\∂p_j})+\\ (-q^T{\Σ}^{-1}\frac{\∂q}{\∂p_i\∂p_j})+(-\frac{\∂q^T}{\∂p_j}{\Σ}^{-1}\frac{\∂q}{\∂p_i}))\end{array}$

The Hessian matrix calculated is optimized mark score (p), obtains new estimates of parameters.New estimates of parameters substitutes into step of registration c, until Δ p meets convergent requirement.

The invention has the beneficial effects as follows: the difficulty overcoming the registration failure because image characteristic point is few in feature based point methods, do not exist when the two width image rotation angle problem of registration failure greatly, the fields such as seafloor topography detection, AUV independent navigation and image procossing can be widely used in yet.

Accompanying drawing explanation

Fig. 1 is schematic flow sheet of the present invention.

Fig. 2 is the results contrast figure of medium filtering image and original image.

Fig. 3 is that Sobel operator calculates gradient magnitude image.

Fig. 4 is that thresholding gradient image obtains whole unique point.

Fig. 5 carries out cluster for using DBSCAN to image characteristic point.

Fig. 6 is registration result image.

Embodiment

Below in conjunction with specific embodiments and the drawings, the present invention will be further described.

As shown in Figure 1, the side scan sonar method for registering of feature based Density Clustering and normal distribution transform specifically comprises the following steps:

Step one, as shown in Figure 2, the original sidescan-sonar image of median filter process:

There is noise in original sidescan-sonar image, median filter first will be used sonar image denoising.Original image is as shown in Fig. 2 (a) (b), and filtered result is as shown in Fig. 2 (c) (d).

Step 2, compute gradient magnitude image:

Use Sobel operator (Sobeloperator, Sobel Operator) to calculate gradient magnitude image, result as shown in Figure 3.

Step 3, as shown in Figure 4, thresholding gradient image obtains whole unique point.

Step 4, uses DBSCAN cluster to obtain image characteristic point.

DBSCAN is a kind of density-based spatial clustering algorithm, this algorithm by the Region dividing with sufficient density for bunch, and find in the noisy spatial database of tool arbitrary shape bunch, it bunch will to be defined as the maximum set of the point that density is connected.The object of DBSCAN algorithm is to filter density regions, finds consistency sample point.Concrete cluster process is as follows:

A. scan whole gradient image, find any one core point, this core point is expanded.The method expanded finds the data point be connected from all density of this core point.

B. travel through all core points in this core point neighborhood, find the point be connected with these pixel point density, until there is no the pixel that can expand.Finally be clustered into bunch boundary node be all non-core data point.

C. after be exactly rescan gradient image (search out before not comprising bunch in any pixel), find not by the core point of cluster, repeat step above again, this core point is expanded until data centralization does not have new core point.The data point that data centralization is not included in any bunch just forms abnormity point.Last data set is exactly the image characteristic point obtained after cluster.Fig. 5 is the experimental result using DBSCAN image characteristic point to be carried out to cluster, by cluster, removes the unique point of mistake.

Step 5, uses NDT algorithm (NormalDistributedTransform, normal distribution transform) to obtain the registration result of two width images.

Two dimensional character point is converted to the probability distribution that zonal cooling can be micro-by NDT algorithm.When having the image of cluster feature point to carry out NDT process to a width, first image is divided into the grating image of formed objects, to the grid cell comprising unique point, carrying out zonal cooling with probability density form to it can be micro-, and represents the probability distribution of unique point in each grid cell with normal distribution:

$p\left(x\right)\∝\exp (-\frac{{(x-q)}^T{\mathrm{\Σ}}^{-1}(x-q)}{2})$

Wherein, x is the two-dimensional coordinate of i point in a grid cell, and q is the mean value of the characteristic point position comprised in grid cell, namely $q=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{x}_i;$ Σ is the covariance matrix of these points, $\mathrm{\Σ}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}(x_i-q){(x_i-q)}^T.$

The spatial transform relation T of two width images can be expressed as:

$T:\left(\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right)=\left(\begin{array}{cc}cos\mathrm{\θ}& -sin\mathrm{\θ}\\ sin\mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$

Wherein, $\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$ Be the translation relation of two width images, θ is the rotation relationship of two width images.

The object of registration is exactly these conversion parameters in order to find between two width sonar images.Specifically comprise the steps:

A. the NDT of a width sonar image is calculated;

B. initialization translation parameters and rotation parameter (use 0 or odometer data carry out initialization);

C. for another width sonar image, transformed in the coordinate system of piece image according to initiation parameter;

D. the normal distribution of each point in sonar image after calculating coordinate change;

E. the normal distribution sum by evaluating each point calculates the mark of translation parameters and rotation parameter;

F. use Hessian matrix method to be optimized these marks, calculate new estimates of parameters;

G. step c is returned, until meet convergent requirement.

Parameter $\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$ Can be expressed as with the mark of θ:

$score\left(p\right)=\underset{i}{\Σ}\exp (-\frac{{(x_i^{\′}-q_i)}^T{\mathrm{\Σ}}_i^{-1}(x_i^{\′}-q_i)}{2})$

Wherein, p=(t _x, t _y, θ) ^t, x _i' be by the some x in the second width image according to coordinate conversion parameter _itransform to coordinate in the first width image coordinate system.Q _iwith x ' respectively _iaverage and covariance.Be optimized mark in above-mentioned e step, make mark score (p) maximum exactly, namely-score (p) is minimum.The detailed operation of Newton iterative is as follows:

In order to make function f minimum, to solve an equation in each iteration:

HΔp＝-g

The Δ p calculated is added with parameter current, obtains new parameter

p←p+Δp

P 0 value initialization, so p=Δ p now.

G is the transposition gradient of f:

$g_i=\frac{\∂f}{\∂p_i}$

H is the Hessian matrix of f:

$H_{ij}=\frac{\∂f}{\∂p_i\∂p_j}$

Definition: q=x ' _i-q _i, then one of-score (p) can be written as:

$s=-\exp \frac{-q^T{\mathrm{\Σ}}^{-1}q}{2}$

The then Grad of s for:

$\begin{array}{c}{\tilde{g}}_i=-\frac{\∂s}{\∂p_j}=-\frac{\∂s}{\∂q}\frac{\∂q}{\∂p_i}\\ =q^T{\Σ}^{-1}\frac{\∂q}{\∂p_i}\exp \frac{-q^T{\Σ}^{-1}q}{2}\end{array}$

Each in Hessian matrix can be expressed as:

$\begin{array}{c}{\tilde{H}}_{ij}=-\frac{\∂s}{\∂p_i\∂p_j}=-\exp \frac{-q^T{\Σ}^{-1}q}{2}\\ \left(\right(-q^T{\Σ}^{-1}\frac{\∂q}{\∂p_i}\left)\right(-q^T{\Σ}^{-1}\frac{\∂q}{\∂p_j})+\\ (-q^T{\Σ}^{-1}\frac{\∂q}{\∂p_i\∂p_j})+(-\frac{\∂q^T}{\∂p_j}{\Σ}^{-1}\frac{\∂q}{\∂p_i}))\end{array}$

The Hessian matrix calculated is optimized mark score (p), obtains new estimates of parameters.New estimates of parameters substitutes into step of registration c, until Δ p meets convergent requirement.Registration result as shown in Figure 6.

By reference to the accompanying drawings embodiments of the invention are elaborated above, but the present invention is not limited to above-described embodiment, in the ken that those of ordinary skill in the art possess, the various changes made under the prerequisite not departing from present inventive concept, all should belong to patent covering scope of the present invention.

Claims (3)

1. a side scan sonar method for registering for feature based Density Clustering and normal distribution transform, comprises the following steps:

Step one, the original sidescan-sonar image of median filter process: have noise in original sidescan-sonar image, first will use median filter to sonar image denoising;

Step 2, compute gradient magnitude image: use Sobel operator to calculate gradient magnitude image;

Step 3, thresholding gradient image obtains whole unique point;

Step 4, uses DBSCAN cluster to obtain image characteristic point: by the Region dividing with sufficient density for bunch, and find in the noisy spatial database of tool arbitrary shape bunch, it bunch will to be defined as the maximum set of the point that density is connected; The object of DBSCAN algorithm is to filter density regions, finds consistency sample point;

Step 5, uses NDT algorithm to obtain the registration result of two width images:

Two dimensional character point is converted to the probability distribution that zonal cooling can be micro-by NDT algorithm; When having the image of cluster feature point to carry out NDT process to a width, first image is divided into the grating image of formed objects, to the grid cell comprising unique point, carrying out zonal cooling with probability density form to it can be micro-, and represents the probability distribution of unique point in each grid cell with normal distribution:

$p\left(x\right)\∝\exp (-\frac{{(x-q)}^T{\mathrm{\Σ}}^{-1}(x-q)}{2})$

Wherein, x is the two-dimensional coordinate of i point in a grid cell, and q is the mean value of the characteristic point position comprised in grid cell, namely $q=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{x}_i;$ Σ is the covariance matrix of these points, $\mathrm{\Σ}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}(x_i-q){(x_i-q)}^T;$

The spatial transform relation T of two width images can be expressed as:

$T:\left(\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right)=\left(\begin{array}{cc}cos\mathrm{\θ}& -sin\mathrm{\θ}\\ sin\mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)+\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$

Wherein, $\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$ Be the translation relation of two width images, θ is the rotation relationship of two width images.

2. the side scan sonar method for registering of feature based Density Clustering according to claim 1 and normal distribution transform, is characterized in that, in step 4, the concrete cluster process of algorithm is as follows:

A. scan whole gradient image, find any one core point, this core point is expanded; The method expanded finds the data point be connected from all density of this core point;

B. travel through all core points in this core point neighborhood, find the point be connected with these pixel point density, until there is no the pixel that can expand; Finally be clustered into bunch boundary node be all non-core data point;

C. rescan gradient image, scans content search out before not comprising bunch in any pixel, find not by the core point of cluster, then repeat step above, this core point is expanded until data centralization does not have new core point; The data point that data centralization is not included in any bunch just forms abnormity point; Last data set is exactly the image characteristic point obtained after cluster.

3. the side scan sonar method for registering of feature based Density Clustering according to claim 1 and normal distribution transform, is characterized in that, in step 5, the object of registration is exactly these conversion parameters in order to find between two width sonar images; Specifically comprise the steps:

A. the NDT of a width sonar image is calculated;

B. initialization translation parameters and rotation parameter;

C. for another width sonar image, transformed in the coordinate system of piece image according to initiation parameter;

D. the normal distribution of each point in sonar image after calculating coordinate change;

E. the normal distribution sum by evaluating each point calculates the mark of translation parameters and rotation parameter;

F. use Hessian matrix method to be optimized these marks, calculate new estimates of parameters;

G. step c is returned, until meet convergent requirement;

Parameter $\left(\begin{array}{c}{t}_x\\ t_y\end{array}\right)$ Can be expressed as with the mark of θ:

$score\left(p\right)=\underset{i}{\Σ}\exp (-\frac{{(x_i^{\′}-q_i)}^T{\mathrm{\Σ}}_i^{-1}(x_i^{\′}-q_i)}{2})$

Wherein, p=(t _x, t _y, θ) ^t, x ' _iby the some x in the second width image according to coordinate conversion parameter _itransform to coordinate in the first width image coordinate system; q _iwith x ' respectively _iaverage and covariance; Be optimized mark in above-mentioned e step, make mark score (p) maximum exactly, namely-score (p) is minimum; The detailed operation of Newton iterative is as follows:

In order to make function f minimum, to solve an equation in each iteration:

HΔp＝-g

The Δ p calculated is added with parameter current, obtains new parameter

p←p+Δp

P 0 value initialization, so p=Δ p now;

G is the transposition gradient of f:

$g_i=\frac{\∂f}{\∂p_i}$

H is the Hessian matrix of f:

$H_{ij}=\frac{\∂f}{\∂p_i\∂p_j}$

Definition: q=x ' _i-q _i, then one of-score (p) can be written as:

$s=-\exp \frac{-q^T{\mathrm{\Σ}}^{-1}q}{2}$

The then Grad of s for:

$\begin{array}{c}{\tilde{g}}_i=-\frac{\∂s}{\∂p_i}=-\frac{\∂s}{\∂q}\frac{\∂q}{\∂p_i}\\ =q^T{\mathrm{\Σ}}^{-1}\frac{\∂q}{\∂p_i}\exp \frac{-q^T{\mathrm{\Σ}}^{-1}q}{2}\end{array}$

Each in Hessian matrix can be expressed as:

$\begin{array}{c}{\tilde{H}}_{ij}=-\frac{\∂s}{\∂p_i\∂p_j}=-\exp \frac{-q^T{\mathrm{\Σ}}^{-1}q}{2}\\ (\left(-q^T{\mathrm{\Σ}}^{-1}\frac{\∂q}{\∂p_i}\right)\left(-q^T{\mathrm{\Σ}}^{-1}\frac{\∂q}{\∂p_j}\right)+\\ \left(-q^T{\mathrm{\Σ}}^{-1}\frac{\∂q}{\∂p_i\∂p_j}\right)+\left(-\frac{\∂q^T}{\∂p_j}{\mathrm{\Σ}}^{-1}\frac{\∂q}{\∂p_i}\right))\end{array}$

The Hessian matrix calculated is optimized mark score (p), obtains new estimates of parameters; New estimates of parameters substitutes into step of registration c, until Δ p meets convergent requirement.