CN106447708A - OCT eye fundus image data registration method - Google Patents

Abstract

The present invention selects the OCT fundus image, uses the Canny edge detection method to extract the retinal edge of the OCT fundus image, and collects these edges in the format of point cloud data; secondly, adopts the method of spatial grid division to extract the feature points of the point cloud data; Thirdly, apply the singular value decomposition algorithm to calculate the transformation matrix between the point clouds to be registered, and eliminate the obvious position error; finally, use the improved iterative closest point algorithm for precise registration, and apply the obtained rotation matrix and translation matrix to the original On the OCT fundus image, the final result is obtained. When dealing with dense point clouds with a large amount of data, the algorithm of the present invention has obvious advantages in terms of time complexity and registration accuracy. In most cases, the algorithm of the present invention improves the efficiency of the traditional iterative closest point algorithm by 70%. The present invention is not only effective for OCT fundus image registration and splicing, but also creates a fundus retinal image with a larger field of view.

Description

A registration method of OCT fundus image data

technical field

The present application relates to image registration technology, in particular, to a method for registration of three-dimensional fundus image data using OCT imaging.

Background technique

Image registration technology usually refers to finding a series of spatial transformations for an image so that it has the same spatial position as the feature information corresponding to another image or multiple images. For the fundus data of OCT three-dimensional imaging, the result of registration should make the layers of the retina in different fundus images aligned with each other without obvious faults.

Currently, image registration techniques are mainly divided into two categories: feature-based image registration and mutual information-based image registration. The feature-based image registration method first extracts the feature information of the image, and then uses these features as a model for registration. The result of feature extraction is a set of data containing image features and a description of these data. Each data is represented by a set of attributes. Further descriptions of attributes include the orientation of edges and the size of radians. These features constitute the local features of the image, and there are interrelationships among the local features, such as geometric relations, radiometric relations, topological relations and so on. Global features can be described by the relationship between these local features. Usually, registration based on local features is mostly based on points, lines or edges, while the registration criterion of global features is a method of registration based on the relationship between local features. Registration methods based on mutual information usually refer to image registration methods with maximum mutual information. Image registration based on mutual information uses the generalized distance between the joint probability distribution of two images and the probability distribution when they are completely independent to estimate mutual information, and it is used as a measure of multimodal medical image registration. When two images based on a common scene achieve the best registration, the gray level mutual information of their corresponding pixels should be the largest. Because the registration based on mutual information is sensitive to noise, in general, it is necessary to preprocess the image by filtering and segmentation, and then perform sampling, transformation, interpolation, optimization and other processes to achieve the purpose of registration.

There are many documents on two-dimensional medical image registration. Kratika Sharma and Ajay Goyal (Sharma K, Goyal A. Classification based survey of image registration methods. International Conference on Computing, Communications & Networking Technologies, 2013, 1-7) divided the steps of image registration methods into It consists of four steps, which are spatial relationship method, scaling method, pyramid method and constant descriptor wavelet method. Mei-sen Pan and other scholars (Pan M, Jiang J, Rong Q, et al. Amodified medical image registration. Multimedia Tools & Applications, 2013, 70: 1585-1615) proposed a B-spline gradient operator method based on edge detection . This method is simple to implement and has low computational load and good registration accuracy. And for two completely overlapping images with noise, Lucian Ciobanu and Luís (Ciobanu L, L. Iterative filtering of SIFT keypoint matches for multi-view registration in Distributed Video Coding. Multimedia Tools & Applications, 2011, 55:557-578) uses an iterative strategy to reconstruct the relationship between point pairs and eliminate noise points. This strategy results in a significant reduction in the overall number of noise points while maintaining the original high-speed accurate matching. To make image registration algorithms more accurate and robust, a novel idea is to model the entire image distribution. Scholars such as Shihui Ying (Ying S, Wu G, Wang Q, et al. Groupwise Registration via Graph Shrinkage on the Image Manifold. IEEE Conference on Computer Vision & Pattern Recognition, 2013, 2323-2330) first introduced this concept, and matched image features The standard procedure can be summarized as a dynamic shrinkage problem. In addition to the above methods, some other feature-based registration algorithms can also get better results.

However, optical coherence tomography fundus data is three-dimensional volume data composed of multiple two-dimensional optical coherence tomography images superimposed. The above methods will face the limitation of operating memory and computing time when processing such data. Therefore, in order to register 3D fundus data, a scheme to obtain optical coherence tomography volume data with a large field of view is to use the montage method (Li Y, Gregori G, Lam B L, et al. Automatic montage of SD-OCT datasets. Optics express, 2011, 19:26239-26248). This method uses resampling, interpolation, and cross-correlation to piece together complete optical coherence tomography volume data using vascular ridges as features of interest. This montage approach stitches scattered, partially overlapping OCT images into a single 3D image with a large field of view. However, this method fails to complete the registration when the vascular ridges, which are features of interest, are blurred in fundus images. In addition, some scholars use existing tools and platforms to generate volumetric data with a larger field of view. Meng Lu (MengL.Acceleration method of 3D medical images registration based on computed device architecture.Bio-medical materials and engineering,2014,24:1109-1116) developed an accelerated registration algorithm based on the computing device architecture provided by NVIDIA , the algorithm can improve the performance of 3D medical image registration, accelerate the calculation speed, and is suitable for processing large-scale data. In addition, scholars such as Stephan Preibisch (Preibisch S, Saalfeld S, Tomancak P. Globally optimal stitching of tiled 3D microscopic image acquisitions. Bioinformatics, 2009, 25:1463-1465) also implemented a stitching plug-in on the imageJ platform, which can be used without In the case of prior knowledge, the scattered small-block data can be reconstructed into a whole. In addition to the above-mentioned stitching plug-ins, other types of stitching tools have also begun to be applied. However, due to the bottleneck of the optical coherence tomography equipment itself and the involuntary movement of the human eye during the scanning process, the obtained optical coherence tomography fundus data will be doped with some small non-rigid transformations (Chen Guolin. Non-rigid medical image Research and Implementation of Registration Method [D]. Master Thesis of Nanjing University of Science and Technology, 2009). Therefore, the above method has certain limitations when dealing with such special clinical ophthalmology optical coherence tomography images.

In summary, for the fundus data of OCT three-dimensional imaging, how to quickly and accurately create an algorithm for fundus images with a larger field of view, so as to provide help for clinicians in the diagnosis and treatment of ophthalmic diseases, is an urgent need to solve in the existing technology technical problems.

Contents of the invention

The object of the present invention is to propose a registration method for OCT fundus image data to quickly and accurately create fundus retinal image data with a larger field of view.

For reaching this purpose, the present invention adopts following technical scheme:

A method for registering OCT fundus image data, comprising the steps of:

Step 1. Use the Canny edge detection method to perform denoising and edge detection on the image:

Grayscale the original image, filter the image, then calculate the gradient magnitude of the image, perform non-maximum suppression on the gradient magnitude, and then use the double threshold algorithm to detect and connect the edges;

Step 2. Extraction and visualization of image point cloud data:

The image processed in step 1 is superimposed into a 3D volume data in the original order, and each edge point is extracted into a point in the 3D point cloud according to its spatial position, and the 3D coordinates of the point cloud points are superimposed correspondingly The spatial coordinates of the obtained three-dimensional volume data;

Step 3. Extraction of point cloud data boundary feature points:

Assign the points in the point cloud to different spatial grids according to their spatial coordinates, then find all the spatial grids belonging to the boundary of the point cloud, and finally extract the point cloud points in the boundary grid as the boundary of the point cloud data Feature points;

Step 4. Use the singular value decomposition algorithm to complete the initial registration of point cloud data:

Let P represent the original set, Q represent the control set, see formula 3, define the target matrix as follows:

C _P and C _Q are the centroids of the original set P and the control set Q respectively, M represents the number of point clouds in the point cloud data set, P _i and Q _i represent the i-th point in the original set P and the control set Q respectively, and The target matrix E adopts singular value decomposition (Singular Value Decomposition, SVD), then E=UDV ^T . where the columns of U are the eigenvectors of the EE ^T matrix, and the columns of V are the eigenvectors of the E ^T E matrix. V ^T is the transpose matrix of V and D is a diagonal matrix, D=diag(d _i ), where d _i is the singular value of E, let

Then the rotation matrix R=UBV ^T , the translation matrix T=C _Q -RC _P , apply the obtained matrices R and T to the original set P to eliminate the large displacement that may exist between the original set P and the control set Q under the initial conditions error;

Step 5. Use the improved iterative closest point algorithm to complete the precise registration of point cloud data:

Step 5.1 Initialization. For given two point cloud sets P and Q, specify a convergence threshold τ,

Step 5.2 Use Formula 5 to calculate the weight corresponding to each point in the point cloud data, compare the threshold ε and exclude points with lower weights,

Q _B is the corresponding point of the point PA in the set Q, Dis(PA, Q _{B )} represents the Euclidean distance between _{PA and Q B , and} Dis _MAX represents the maximum value of the Euclidean distance between point pairs, and the Euclidean distance adopts the formula 6 is calculated to get,

Step 5.3 Iterate the following steps until the minimum root mean square error of Equation 7 converges to the given threshold τ:

Step 5.3.1 calculates the Euclidean distance between point clouds in sets P and Q according to formula 6,

Among them, ω _x , ω _y , ω _z respectively represent the weight factors of M-estimation in each coordinate direction, (x _A , y _A , z _A ), (x _B , y _B , z _B ) are the midpoints of the set P The spatial coordinates of point B in A and Q,

Step 5.3.2 For the set P that excludes points with lower weights, find the point with the closest Euclidean distance in the set Q as the corresponding point and store it in the nearest point set,

Step 5.3.3 Use formula 7 to calculate the rotation matrix R and translation matrix T between the calculation set P and the nearest point set using the least square method,

Step 5.3.4 Apply the rotation matrix R and the translation matrix T to the set P to obtain a new set, use formula 7 to calculate whether the minimum root mean square error converges to the given threshold τ, if so, end the operation, otherwise use the step 5.3 Perform iterative calculations.

Further, in step 1, the grayscale formula of the original image is Gray=0.299R+0.587G+0.114B, the image filtering adopts Gaussian filtering, and the first-order finite difference is used to calculate the partial derivative of the image pixel matrix with respect to the horizontal and vertical directions .

Further, in step 1, double thresholds are used to detect and connect edges. The function of high threshold is to suppress edges, and all gradient magnitudes higher than this threshold can be regarded as edges. The function of low threshold is to connect edges, and all Points whose gradient magnitude is higher than a low threshold are regarded as edges and connected to become the result of the final edge detection.

Further, in step 3, the directed bounding box is used as the minimum bounding box of the point cloud data to divide the point cloud data points into different spatial grids,

Decompose the point cloud data into several spatial grids according to the same volume, and the size of each spatial grid is defined as:

In the formula, L represents the number of point cloud points in the point cloud data, and V represents the volume of the smallest bounding box. V/L is the reciprocal of the point cloud density, representing the average size of the space occupied by each point cloud point in the point cloud data. Let the initial size S _grid of the spatial grid be K times the reciprocal of the point cloud density, and divide the minimum bounding box of the point cloud data into several spatial grids according to the size of S _grid .

Further, use the boundary seed grid algorithm to find all the boundary space grids, and extract the point cloud points contained in this boundary space grid as the boundary feature points of the point cloud data, and divide the space grids into two types, space and real grid , the spatial grids that do not contain any point cloud points are empty spaces, and the other spatial grids are solid grids. Use spatial coordinates (x, y, z) to represent the spatial position of a certain grid, and use the function f to represent the type of this spatial grid. If a certain spatial grid is a real grid, f(x, y, z) = 1, Otherwise f(x,y,z)=0, use formula 2 to determine whether a certain spatial grid is the boundary spatial grid:

When U(x,y,z)≤1, it means that at most one pair of the upper, lower, left and right, front and rear grids in the six neighborhoods of this spatial grid is a real grid, then this spatial grid is a boundary grid .

Further, in step three, 8≤K≤24.

Further, in step five, the convergence threshold τ=0.2, and the comparison threshold ε is 0.2≤ε≤0.4.

Further, in step five, the weight formula of M-estimation in each coordinate direction is

v is the standardized residual value in each coordinate direction, and c is a constant.

Further, c=1.345.

The present invention selects the OCT fundus image, uses the Canny edge detection method to extract the retinal edge of the OCT fundus image, and collects these edges in the format of point cloud data; secondly, adopts the method of spatial grid division to extract the feature points of the point cloud data; Thirdly, apply the singular value decomposition algorithm to calculate the transformation matrix between the point clouds to be registered, and eliminate the obvious position error; finally, use the improved iterative closest point algorithm for precise registration, and apply the obtained rotation matrix and translation matrix to the original On the OCT fundus image, the final result is obtained.

When dealing with dense point clouds with a large amount of data, the algorithm of the present invention has obvious advantages in terms of time complexity and registration accuracy. In most cases, the algorithm of the invention improves the efficiency of the traditional iterative closest point algorithm by 70%. The invention is not only effective for registration and splicing of OCT fundus images, but also creates a fundus retinal image with a large visual field with high accuracy.

Description of drawings

Fig. 1 is the flow chart of the OCT fundus image data registration method according to the specific embodiment of the present invention;

Fig. 2 is the state diagram of the boundary grid adopting boundary grid seed algorithm according to the specific embodiment of the present invention;

Fig. 3 is a comparison diagram of point cloud data registration process according to a specific embodiment of the present invention;

Fig. 4 is according to the specific embodiment of the present invention the enlarged comparison diagram of the partial results of point cloud data overlapping area;

Fig. 5 is a position comparison diagram of the transformation matrix acting on the original set and the control set after registration according to a specific embodiment of the present invention;

Fig. 6 is a comparison diagram of registration results for "Stanford Rabbit" point cloud data according to a specific embodiment of the present invention;

Fig. 7 is a comparison diagram of the time consumption of the improved algorithm according to the present invention and the traditional iterative algorithm.

detailed description

The present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described here are only used to explain the present invention, but not to limit the present invention. In addition, it should be noted that, for the convenience of description, only some structures related to the present invention are shown in the drawings but not all structures.

The principle of the present invention is: use the Canny edge detection method to denoise the medical fundus image, and use the obtained edge as the feature of the image, extract these features into point cloud data, and use the method of the minimum bounding box to extract the point cloud data. The boundary feature points of , use the boundary feature points to complete the initial registration by singular value decomposition, and accurately register the point cloud data by improving the iterative closest point method.

Referring to Fig. 1, it shows the flow chart of OCT fundus image data registration method according to the present invention, specifically comprises the following steps:

Step 1. Use the Canny edge detection method to perform denoising and edge detection on the image:

The edge of the image refers to the critical area where the image changes sharply from a certain gray value to another gray level, and it usually separates two areas where the brightness changes significantly. The edge of the image carries most of the information of the image and is an important part of image feature extraction. The first step includes the following sub-steps:

(1) Grayscale the original image: In practical application scenarios, most images are in color format. The image in color format needs to be grayscaled first. A more common grayscale method is to perform weighted average of the sampling values of the three RGB channels of the image, that is, Gray=(R+G+B)/3. For medical ophthalmology images, considering the physiological characteristics of the human eye, the gray scale formula in the present invention is Gray=0.299R+0.587G+0.114B.

(2) Filter the image, and then calculate the gradient magnitude of the image:

Image filtering is an important process in image processing. It suppresses and eliminates the influence of noise on image processing while retaining the original features of the image as much as possible, thereby improving the accuracy and reliability of image processing. Common filtering methods mainly include mean filtering, median filtering, Gaussian filtering and bilateral filtering. The Canny edge detection algorithm in the present invention adopts Gaussian filtering.

The edge of the image is generally in the area where the gray value changes greatly, and the first-order finite difference can be used to calculate the partial derivative of the image pixel matrix with respect to the horizontal and vertical directions. For a point with a larger partial derivative, it means that the gray level of the area where the point is located has a large change, and it may be an edge point of the image.

(3) Perform non-maximum suppression on the gradient magnitude, and then use a double threshold algorithm to detect and connect edges:

From the conclusion of the low-order derivative, it can be seen that the point with the maximum gradient amplitude is not necessarily the point with the largest change in the gray value of the region, and it does not mean that this point must be an edge point, so the Canny edge detection algorithm in the present invention has carried out non- Maximum suppression. That is to search for the eight neighbors of the maximum gradient amplitude point, keep the point with the largest gradient amplitude value, and omit the other points to ensure that there is no more than one edge point in any small area.

The Canny edge detection algorithm uses a double threshold to detect and connect edges. The role of a high threshold is to suppress edges, and all gradient magnitudes above this threshold are considered edges. Therefore, the higher the high threshold, the fewer edges are detected. The role of the low threshold is to connect the edges. Due to the inhibition of the high threshold, many detected edge points cannot be connected into an edge. At this time, the Canny edge detection algorithm regards all points whose gradient magnitude is higher than the low threshold among the edge points as edges and connects them to become the result of the final edge detection.

Step 2. Extraction and visualization of image point cloud data

Point cloud data consists of a set of points in three-dimensional space, and is usually used to describe the three-dimensional spatial structure of an object. Each point cloud point contains at least a set of three-dimensional coordinates representing its own spatial position, and some point clouds also have high-dimensional coordinates, which are used to represent information such as color or light reflection intensity. After the original fundus optical coherence tomography image is processed in step 1, a series of images including the edges of each retinal layer are obtained.

Superimpose these images into a 3D volume data in the original order, and extract each edge point into a point in the 3D point cloud according to its spatial position, and the 3D coordinates of the point cloud points correspond to the 3D volume obtained by superposition The spatial coordinates of the data.

Step 3. Extraction of point cloud data boundary feature points

In this step, the points in the point cloud are first assigned to different spatial grids according to their spatial coordinates, then all the spatial grids belonging to the point cloud boundary are found, and finally the point cloud points in the boundary grid are extracted As the boundary feature points of point cloud data.

Furthermore, to divide the point cloud data points into different spatial grids, it is first necessary to find the smallest bounding box that can contain the entire point cloud, and then divide the smallest bounding box into several spatial grids, so that any point cloud data Points are contained in a certain spatial grid.

The present invention uses a directional bounding box as the smallest bounding box of point cloud data, and the directional bounding box is a bounding box method with any direction, which can contain the whole object more compactly. After obtaining the minimum bounding box of the point cloud data, the present invention decomposes the point cloud data into several spatial grids according to the same volume, and the size of each spatial grid is defined as:

In the formula, L represents the number of point cloud points in the point cloud data, and V represents the volume of the smallest bounding box. V/L is the reciprocal of the point cloud density, representing the average size of the space occupied by each point cloud point in the point cloud data. Let the initial size S _grid of the spatial grid be K times the reciprocal of the point cloud density, and divide the minimum bounding box of the point cloud data into several spatial grids according to the size of S _grid . K in Formula 1 is a variable parameter. For the intensive point cloud data processed by the present invention, experiments have proved that when K is assigned a value of 8-24, the size of the spatial grid can contain enough point cloud data points.

Furthermore, the present invention uses a boundary seed grid algorithm to find all boundary space grids, and extracts point cloud points contained in the boundary space grids as boundary feature points of point cloud data. In order to extract the boundary spatial grid, the spatial grid can be divided into two types, blank and real grid, the spatial grid that does not contain any point cloud points is blank, and the other spatial grids are real grid. Use spatial coordinates (x, y, z) to represent the spatial position of a certain grid, and use the function f to represent the type of this spatial grid. If a certain spatial grid is a real grid, f(x, y, z) = 1, Otherwise f(x,y,z)=0. The coordinates of the six neighbor grids adjacent to this grid in space are (x-1, y, z), (x+1, y, z), (x, y-1, z), (x, y+1, z), (x, y, z-1), (x, y, z+1). Use formula 2 to determine whether a certain spatial grid is the boundary spatial grid:

U(x,y,z) is the sum of three products, let the coordinates of the middle grid be (x,y,z), then (x-1,y,z), (x+1,y,z) , (x,y-1,z), (x,y+1,z), (x,y,z-1), (x,y,z+1) represent the six-neighborhood grid of this grid location. f(x-1,y,z)·f(x+1,y,z)=1 if and only if the upper and lower positions of the grid are real cells, f(x,y-1,z)·f(x ,y+1,z)=1 if and only if the left and right positions of the grid are real cells, f(x,y,z-1)·f(x,y,z+1)=1 if and only if the grid The front and rear positions of the grid are all real grids. When U(x,y,z)≤1, it means that at most one pair of the upper, lower, left and right, front and rear grids in the six neighborhoods of this spatial grid is a real grid, then this spatial grid is a boundary grid .

Referring to FIG. 2 , it is a state diagram of a boundary grid using a boundary grid seed algorithm according to a specific embodiment of the present invention. The left figure is a plane boundary, the middle figure is a line boundary, and the right figure is a point boundary.

Step 4. Use the singular value decomposition algorithm to complete the initial registration of point cloud data

For the boundary feature points of the above point cloud data, the singular value decomposition algorithm is used to calculate the rotation matrix and translation matrix between the registration sets to complete the initial registration of the point cloud data.

Let P represent the original set, Q represent the control set, see formula 3, define the target matrix as follows:

C _P and C _Q are the centroids of the original set P and the control set Q respectively, M represents the number of point clouds in the point cloud data set, P _i and Q _i represent the i-th point in the original set P and the control set Q respectively, and The target matrix E adopts singular value decomposition (Singular Value Decomposition, SVD), then E=UDV ^T . where the columns of U are the eigenvectors of the EE ^T matrix, and the columns of V are the eigenvectors of the E ^T E matrix. V ^T is the transpose matrix of V and D is a diagonal matrix, D=diag(d _i ), where d _i is the singular value of E, that is, the square root of the eigenvalue of the E'E matrix, so that

Then the rotation matrix R=UBV ^T , the translation matrix T=C _Q -RC _P , apply the obtained matrices R and T to the original set P to eliminate the large displacement that may exist between the original set P and the control set Q under the initial conditions An error, eg, P'=RP+T, provides a good state for the precise registration step.

Step 5. Use the improved iterative closest point algorithm to complete the precise registration of point cloud data

The main idea of the traditional iterative closest point algorithm is to find out all the corresponding points of P and Q according to some geometric features, and use these corresponding points as the objects of registration. Then construct the objective equation representing the degree of registration, and calculate the rotation matrix and translation matrix in the current situation by finding the optimal solution of the objective equation. Finally, this process is iterated continuously until the objective equation satisfies a certain set threshold.

The traditional iterative closest point algorithm uses the point-to-point distance as a geometric feature to find the corresponding point. For two points P _i (x _p ,y _p ,z _p ) and Q _i (x _q ,y _q ,z _q ) in sets P and Q, the Euclidean distance between them is expressed as

For any point in the set P, use the above formula to calculate the distance from each point in Q to this point, and select the point with the smallest Euclidean distance as the corresponding point of this point in the set Q. Then find the rotation matrix R and the translation matrix T, and apply them to P _i , then the position of the obtained point is RP _i +T, and use the least squares method to construct the objective equation

Among them, N represents the number of points in the point cloud, and E represents the sum of the squares of the distances from each point in the set P to the corresponding point in the set Q after transformation.

It can be seen from the above objective equation that when E is the smallest, it means that the relative position of the corresponding point of this iteration is the closest. Therefore, the rotation matrix and translation matrix that make E the smallest are the optimal solutions for this iteration. For the solution of the rotation matrix and translation matrix, the method of separating translation and rotation is adopted, that is, the initial value estimation of the translation matrix T is performed first, and the centers of gravity of P and Q are first calculated as

and

Then the estimated value of translation between sets P and Q is p _c -q _c , and the objective equation is simplified as

In this way, the rotation matrix R in this iteration process is obtained, and the translation matrix T is calculated by using T=Q-RP, and the above process is repeated until the optimal solution E of the objective equation converges to a given threshold.

For point cloud data with better initial conditions, the iterative closest point algorithm can obtain more accurate registration results. However, the traditional iterative closest point algorithm also has some shortcomings.

First of all, the iterative closest point algorithm assumes that two point cloud sets are inclusive, that is, one point cloud set is a subset of the other point cloud set, which is usually difficult to satisfy.

Secondly, in the process of selecting corresponding point pairs, for any point in the point cloud data set, the algorithm will calculate the Euclidean distance between this point and all points in another set. Assuming that the two point cloud sets each have N _P and N _Q points, then the time complexity of this step is O(N _P ×N _Q ). This makes the algorithm spend a lot of time calculating the Euclidean distance of the corresponding point pairs, and the calculation cost is very high. In addition, the algorithm defaults the point pair with the closest Euclidean distance as the corresponding point. Due to the existence of a certain amount of noise points, a certain error occurs after this step, which makes the algorithm fall into a local minimum.

Therefore, the present invention proposes an improved iterative closest point algorithm: in the traditional iterative closest point algorithm, all points in the point cloud collection are given the same weight, so all points will participate in the calculation process of the Euclidean distance between points , which is the bottleneck of the algorithm. The algorithm of the invention assigns different weights to different points, and the farther the Euclidean distance between point pairs is, the smaller their weights are. Suppose P _A is a point in the set P, then the weight of P _A is expressed as

Q _B is the corresponding point of point PA in set Q, and Dis(PA , Q _{B )} represents the Euclidean distance between _PA and _{Q B.} Dis _MAX represents the maximum value of the Euclidean distance between point pairs. For a given threshold ε, all corresponding point pairs with a weight lower than ε will be discarded and will not participate in the iterative calculation process. Threshold ε is a variable parameter, which is used to trade off the time complexity of registration and the accuracy of registration. If ε is reduced, it means that more corresponding point pairs will participate in the iterative calculation process, making the registration result more accurate, but at the same time it also increases the number of calculations and increases the time complexity of the algorithm. Preferably, ε is 0.2-0.4.

Because the point cloud data may have noise, these noises will affect the calculation of the center of gravity of the point cloud and the calculation of the iterative nearest point objective equation, thereby reducing the registration accuracy. The present invention introduces M-estimation to eliminate the influence of noise points on the experimental results. The basic idea of M-estimation robust regression is to determine the weight of each point according to the size of the regression residual to achieve the purpose of robustness. In order to reduce the influence of abnormal points, different weights can be assigned to different points, that is, a larger weight is given to a point with a small residual error, and a smaller weight is given to a point with a larger residual error. In the present invention, the weight formula of M-estimation in each coordinate direction is

The weight factor for a point is defined as

v is the standardized residual value in each coordinate direction, c is a constant, generally 1.345. Therefore, when v belongs to the interval (-c,c), the M-estimator degenerates into a classical least squares estimate. And when the residual v is greater than c, the weight factor w _v decreases as the residual increases. The calculation method of the Euclidean distance between the middle point A of the corresponding set P and the middle point B of Q is shown in formula 6:

Among them, ω _x , ω _y , ω _z respectively represent the weight factors of M-estimation in each coordinate direction, (x _A , y _A , z _A ), (x _B , y _B , z _B ) are points A and B respectively The space coordinates of , the iterative equation at this time is as follows:

In summary, the improved iterative closest point algorithm is as follows:

Step 5.1 Initialization. For given two point cloud sets P and Q, specify a convergence threshold τ, preferably, τ is 0.2;

Step 5.2 Use formula 5 to calculate the weight corresponding to each point in the point cloud data, compare the threshold ε and exclude points with lower weight, preferably, ε is 0.2 to 0.4;

Step 5.3 Iterate the following steps until the minimum root mean square error of Equation 7 converges to the given threshold τ:

Step 5.3.1 calculates the Euclidean distance between point clouds in sets P and Q according to formula 6,

Step 5.3.2 For the set P that excludes points with lower weights, find the point with the closest Euclidean distance in the set Q as the corresponding point and store it in the nearest point set,

Step 5.3.3 Use formula 7 to calculate the rotation matrix R and translation matrix T between the calculation set P and the nearest point set using the least square method,

Step 5.3.4 Apply the rotation matrix R and the translation matrix T to the set P to obtain a new set, use formula 7 to calculate whether the minimum root mean square error converges to the given threshold τ, if so, end the operation, otherwise use the step 5.3 Perform iterative calculations.

Therefore, the present invention selects the OCT fundus image, uses the Canny edge detection method to extract the retinal edge of the OCT fundus image, and collects these edges in the form of point cloud data; secondly, adopts the method of spatial grid division to extract the features of the point cloud data point; again, apply the singular value decomposition algorithm to calculate the transformation matrix between the point clouds to be registered, and eliminate the obvious position error; finally, use the improved iterative closest point algorithm for precise registration, and apply the obtained rotation matrix and translation matrix to On the original OCT fundus image, the final result is obtained.

Example 1:

In the following, the present invention will be further described in combination with experiments, and the present invention will be compared with the classic iterative closest point method to verify its accuracy and robustness.

The computer configuration used in the experiment is Intel Core E7500 dual-core CPU, the main frequency is 2.93GHz, the memory is 1GB×2 DDR2, the operating system is Microsoft Windows Win7 Sp1 flagship version, and the algorithm implementation platform is MicrosoftVisual Studio 2010.

In order to evaluate the registration performance of the algorithm, the algorithm of the present invention uses two fundus optical tomography images as experimental data, which are adjacent parts of the human retina structure. The overlap of the two datasets is approximately 75 x 500 x 375 voxels. Figure 3 shows the registration process of point cloud data from four different perspectives. The left area of the picture shows the initial positions of the two point cloud sets (the red points represent the original set, and the green points represent the control set), the right area shows the real-time registration results during the iteration process, and the outside of the right point cloud A contour represents the minimal directed bounding box of a point cloud.

After the above iterative process ends, the final registration result is obtained. In order to make the experimental results more intuitive, Figure 4 shows a partial enlarged view of the registration results. As shown in Figure 4, the left area of the image shows the initial positions of the original set and the control set, and it can be observed that there is a large misalignment between them. The image on the left shows the positional relationship of the overlapping areas of two data sets before registration, and the image on the right shows the registration result of the algorithm in this paper. The result processed by the algorithm of the present invention is displayed on the right side of the image, and the image shows that most of the misaligned points have been relatively accurately registered.

After the registration process, we apply the obtained series of transformation matrices to the original set and the control set, and render them with imageJ. Figure 5 visualizes the original data and the resulting data after registration. Figure 5 renders four fundus optical coherence tomography image volume data. The above two images are the above-mentioned original set and control set respectively, the latter two are the experimental results of the present invention, the lower left image is the side view of the registration result, the lower right is the top view of the registration result, and the concave part in the figure is the human eye macular center. The traditional fundus volume data registration algorithm often leads to the phenomenon of cliffs in the registered retinal inner nucleus layer, photoreceptor cell layer and pigment epithelium layer. However, from the results shown in Figure 5, there is no obvious retinal layer fault in the registration results, and no clear mosaic traces can be observed in the overlapping areas of the data. This optical coherence tomography with a larger field of view Body data may be of some help to clinicians in the prevention and diagnosis of ophthalmic diseases.

In the present invention, Figs. 3-5 are only used to show the effects obtained by the present invention, but the effects of the present invention are not only expressed depending on Figs. 3-5.

Example 2

In addition, the present invention uses the registration error to evaluate the registration accuracy of the algorithm, and uses the time consumption length to evaluate the registration efficiency of the algorithm. The registration error represents the percentage of the number of points that failed to match corresponding points in the registration process to the total number of point cloud points, and it has the following calculation form

where Success(P _A ,Q _B ) is defined as

In the above formula, N represents the number of corresponding point pairs involved in the calculation. (P _A , Q _B ) represents a pair of corresponding points. Success(P _A , Q _B ) is used to represent the registration result of the corresponding point pair (P _A , Q _B ). If the Euclidean distance between the corresponding point pairs after registration is less than the given threshold δ, it means that the registration of the corresponding point pair is successful, and Success(P _A , Q _B )=1. _The value of δ in the algorithm of the present invention is 0.15 times of the Euclidean distance before accurate registration (PA, QB ₎ , that is, the corresponding point pair whose Euclidean distance is less than 15% of the original distance after registration is regarded as successful registration. Through statistics and calculation, Table 1 shows the comparison between the algorithm of the present invention and the traditional iterative closest point algorithm proposed by Besl in terms of time consumption and registration error.

Table 1 The comparison of the experimental results of the algorithm of the present invention and the traditional iterative closest point algorithm

As an experimental comparison, the present invention conducts experiments on some open source point cloud datasets. Figure 6 shows the registration results of the point cloud data of the public library "Stanford Rabbit" using the algorithm of the present invention. The left picture shows the initial position of the point cloud data, and the right picture shows the registered point cloud data. This point cloud data has a total of 9731 point cloud points. It takes 11.023 seconds to use the traditional iterative closest point algorithm, and the registration error is 0.001130. Using the improved iterative closest point algorithm takes 3.994 seconds and the registration error is 0.000794. By improving the traditional iterative closest point algorithm, the present invention significantly improves the registration of point cloud data in terms of calculation time and registration accuracy.

In order to compare the improved algorithm and the traditional algorithm more intuitively, the present invention conducts comparative experiments on point cloud data of different sizes, and uses a line graph to show the difference in calculation time of the algorithm. In Table 2, the first data set with a point cloud data volume of 9731 is "Stanford Rabbit", and the second is "Dragon" in the Stanford 3D scanning warehouse (http://graphics.stanford.edu/data/3Dscanrep/) , and the rest are point cloud data in the actual project. Table 2 and Fig. 7 below show the comparison of registration time between the algorithm of the present invention and the traditional iterative closest point algorithm for different data sets.

Table 2 Point cloud data collections of different sizes, the calculation time comparison between the algorithm of the present invention and the traditional iterative closest point algorithm

As shown in Table 2 and Figure 7, when dealing with dense point clouds with a large amount of data, the algorithm of the present invention has obvious advantages in terms of time complexity and registration accuracy. In most cases, the algorithm of the invention improves the efficiency of the traditional iterative closest point algorithm by 70%.

The above experimental results verify the effectiveness of the present invention for OCT fundus image registration and splicing, and at the same time verify the accuracy and robustness of the present invention for creating retinal images with a larger field of view.

The above content is a further detailed description of the present invention in conjunction with specific preferred embodiments. It cannot be determined that the specific embodiments of the present invention are limited thereto. Under the present invention, several simple deduction or substitutions can also be made, all of which should be considered as belonging to the protection scope of the present invention determined by the submitted claims.

Claims (9)

1. A method for registration of OCT fundus image data, comprising the steps of: Step 1. Use the Canny edge detection method to perform denoising and edge detection on the image: Grayscale the original image, filter the image, then calculate the gradient magnitude of the image, perform non-maximum suppression on the gradient magnitude, and then use the double threshold algorithm to detect and connect the edges; Step 2. Extraction and visualization of image point cloud data: The image processed in step 1 is superimposed into a 3D volume data in the original order, and each edge point is extracted into a point in the 3D point cloud according to its spatial position, and the 3D coordinates of the point cloud points are superimposed correspondingly The spatial coordinates of the obtained three-dimensional volume data; Step 3. Extraction of point cloud data boundary feature points: Assign the points in the point cloud to different spatial grids according to their spatial coordinates, then find all the spatial grids belonging to the boundary of the point cloud, and finally extract the point cloud points in the boundary grid as the boundary of the point cloud data Feature points; Step 4. Use the singular value decomposition algorithm to complete the initial registration of point cloud data: Let P represent the original set, Q represent the control set, see formula 3, define the target matrix as follows: C _P and C _Q are the centroids of the original set P and the control set Q respectively, M represents the number of point clouds in the point cloud data set, P _i and Q _i represent the i-th point in the original set P and the control set Q respectively, and The target matrix E adopts singular value decomposition (Singular Value Decomposition, SVD), then E=UDV ^T . where the columns of U are the eigenvectors of the EE ^T matrix, and the columns of V are the eigenvectors of the E ^T E matrix. V ^T is the transpose matrix of V and D is a diagonal matrix, D=diag(d _i ), where d _i is the singular value of E, let Then the rotation matrix R=UBV ^T , the translation matrix T=C _Q -RC _P , apply the obtained matrices R and T to the original set P to eliminate the large displacement that may exist between the original set P and the control set Q under the initial conditions error; Step 5. Use the improved iterative closest point algorithm to complete the precise registration of point cloud data: Step 5.1 Initialization. For given two point cloud sets P and Q, specify a convergence threshold τ, Step 5.2 Use Formula 5 to calculate the weight corresponding to each point in the point cloud data, compare the threshold ε and exclude points with lower weights, Q _B is the corresponding point of the point PA in the set Q, Dis(PA, Q _{B )} represents the Euclidean distance between _{PA and Q B , and} Dis _MAX represents the maximum value of the Euclidean distance between point pairs, and the Euclidean distance adopts the formula 6 calculated to get, Step 5.3 Iterate the following steps until the minimum root mean square error of Equation 6 converges to the given threshold τ: Step 5.3.1 calculates the Euclidean distance between point clouds in sets P and Q according to formula 6, Among them, ω _x , ω _y , ω _z respectively represent the weight factors of M-estimation in each coordinate direction, (x _A , y _A , z _A ), (x _B , y _B , z _B ) are the midpoints of the set P The spatial coordinates of point B in A and Q, Step 5.3.2 For the set P that excludes points with lower weights, find the point with the closest Euclidean distance in the set Q as the corresponding point and store it in the nearest point set, Step 5.3.3 Use formula 7 to calculate the rotation matrix R and translation matrix T between the calculation set P and the nearest point set using the least square method, Step 5.3.4 Apply the rotation matrix R and the translation matrix T to the set P to obtain a new set, use formula 7 to calculate whether the minimum root mean square error converges to the given threshold τ, if so, end the operation, otherwise use the step 5.3 Perform iterative calculations. 2. the OCT fundus image data registration method according to claim 1, is characterized in that: In step 1, the grayscale formula of the original image is Gray=0.299R+0.587G+0.114B, the image filtering adopts Gaussian filtering, and the first-order finite difference is used to calculate the partial derivative of the image pixel matrix with respect to the horizontal and vertical directions. 3. OCT fundus image data registration method according to claim 2, is characterized in that: In step 1, double thresholds are used to detect and connect edges. The function of high threshold is to suppress edges, and all gradient magnitudes higher than this threshold can be regarded as edges. The function of low threshold is to connect edges, and all gradient magnitudes Points above a low threshold are considered edges and concatenated to become the final edge detection result. 4. OCT fundus image data registration method according to claim 1, is characterized in that: In step three, the directed bounding box is used as the minimum bounding box of the point cloud data to divide the point cloud data points into different spatial grids, Decompose the point cloud data into several spatial grids according to the same volume, and the size of each spatial grid is defined as: In the formula, L represents the number of point cloud points in the point cloud data, and V represents the volume of the smallest bounding box. V/L is the reciprocal of the point cloud density, representing the average size of the space occupied by each point cloud point in the point cloud data. Let the initial size S _grid of the spatial grid be K times the reciprocal of the point cloud density, and divide the minimum bounding box of the point cloud data into several spatial grids according to the size of S _grid . 5. the OCT fundus image data registration method according to claim 4, is characterized in that: Use the boundary seed grid algorithm to find all the boundary space grids, and extract the point cloud points contained in this boundary space grid as the boundary feature points of the point cloud data, and divide the space grids into two types, blank and real grids, which do not contain The spatial grid of any point cloud point is a blank space, and the other spatial grids are solid grids. Use spatial coordinates (x, y, z) to represent the spatial position of a certain grid, and use the function f to represent the type of this spatial grid. If a certain spatial grid is a real grid, f(x, y, z) = 1, Otherwise f(x,y,z)=0, use formula 2 to determine whether a certain spatial grid is the boundary spatial grid: When U(x,y,z)≤1, it means that at most one pair of the upper, lower, left and right, front and rear grids in the six neighborhoods of this spatial grid is a real grid, then this spatial grid is a boundary grid . 6. the OCT fundus image data registration method according to claim 4 or 5, is characterized in that: In step three, 8≤K≤24. 7. OCT fundus image data registration method according to claim 1, is characterized in that: In step five, the convergence threshold τ=0.2, and the comparison threshold ε is 0.2≤ε≤0.4. 8. OCT fundus image data registration method according to claim 7, is characterized in that: In step five, the weight formula of M-estimation in each coordinate direction is $\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; =</span>\left\{\begin{array}{cc}\mathrm{<span\; class="notranslate">\; 1</span>}& <span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&le;</span>\mathrm{<span\; class="notranslate">\; c</span>}\\ \frac{\mathrm{<span\; class="notranslate">\; c</span>}}{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; |</span>}& <span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; v</span>}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; ></span>\mathrm{<span\; class="notranslate">\; c</span>}\end{array}\right.$ v is the standardized residual value in each coordinate direction, and c is a constant. 9. OCT fundus image data registration method according to claim 8, is characterized in that: c=1.345.