Retinal image registration method for myopia development

Abstract

The application faces the retinal image registration method for the development of myopia and belongs to the field of retinal image registration. The application firstly adopts a SuperPoint deep learning framework to detect and describe key points in a retinal image pair, then adopts a SuperGlue deep learning framework to match the key point set between the retinal image pair, and then estimates the initial pose of the test camera by using a PnP algorithm based on RANSAC. Finally, a particle swarm optimization algorithm PSO is used to optimize the three-dimensional space model, the three-dimensional space model increases the distortion model and the myopia development model, and the test image is registered and transformed accordingly. The method can provide the best registration accuracy in the case of myopia development or not.

Description

Technical Field

[0001] This invention belongs to the field of retinal image registration methods, specifically relating to a retinal image registration method for the development of myopia. Background Technology

[0002] Retinal image registration works on a pair of retinal images, with the goal of aligning the two images spatially. The Registration of Retinal Images Through Eye Modelling and Pose Estimation (REMPE) method, proposed by Hernandez-Matas et al., uses corresponding points between retinal image pairs as a medium to achieve registration by jointly estimating the relative pose between the cameras in the images and the ellipsoidal model parameters of the eye.

[0003] In the ideal process of acquiring fundus retinal images, the fundus camera typically has a fixed pose, with the eye centered in front of the camera. Its rotation can provide the camera with views of different retinal regions. In other words, the main reason why retinal images of the same person taken at different times need to be registered is the change in eye pose. REMPE takes this as its starting point and establishes the equivalent geometric structure between the pairs of retinal images to be registered.

[0004] The reference camera for acquiring the reference image is fixed, while the eye rotation and the pose of the test camera for acquiring the test image are defined relative to the reference camera. Registration of the reference and test images is achieved by simultaneously estimating the relative pose of the camera when acquiring the retinal image, as well as the shape and rotation of the eye. REMPE utilizes the point correspondence between the images for registration, which is accomplished by searching for the camera pose, eye shape, and pose that best interpret the 2D coordinate correspondence between the matching points. Its workflow is roughly as follows: First, keypoints are extracted and matched in the two retinal images to be registered. Then, an initial pose estimate is calculated using RANSAC and a spherical model. Subsequently, the pose is further optimized using a particle swarm optimization algorithm, while simultaneously estimating the semi-axis length and rotation angle of the ellipsoidal eye model. Finally, the optimal geometry is used for image registration transformation and output.

[0005] The shortcomings and deficiencies of the REMPE retinal image registration method during use are as follows: 1. On the FIRE retinal image registration dataset, which does not show any myopia progression, REMPE does not achieve the best registration accuracy, does not surpass VOTUS, and still has a certain gap with VOTUS in registration accuracy; 2. For retinal image pairs with obvious myopia progression, REMPE is difficult to achieve satisfactory registration results, and it cannot cope with the problem of the significant increase in the complexity of retinal image registration caused by the progression of myopia.

[0006] The main reason for these shortcomings is that for retinal image pairs showing signs of myopia development, it is necessary to superimpose the changes in the shape of the retinal surface caused by axial elongation and fundus curvature, which undoubtedly greatly increases the complexity of registration. However, existing retinal image registration methods (including the REMPE retinal image registration method) implicitly assume that the corresponding eyeballs or fundus retinal surfaces in the two images are identical in shape, and do not consider the increased registration difficulty and transformation model complexity caused by changes in the shape of the eyeballs or fundus retinal surfaces. Therefore, for retinal image pairs showing significant signs of myopia development, existing methods usually fail to achieve satisfactory registration results. Summary of the Invention

[0007] This invention addresses the problems of existing retinal image registration methods by proposing a retinal image registration method tailored to the development of myopia. The registration method of this invention incorporates a distortion model of the fundus camera during retinal imaging and constructs and introduces a novel myopia development model to approximately simulate the changes in axial elongation and fundus curvature caused by the development of myopia. This provides optimal registration accuracy regardless of whether myopia has developed.

[0008] To achieve the above objectives, the present invention provides a retinal image registration method for the development of myopia, comprising the following steps:

[0009] Step 1) Use the SuperPoint deep learning framework to detect and describe key points in retinal image pairs, and then use the SuperGlue deep learning framework to match key point sets between retinal image pairs;

[0010] Step 2) Estimate the initial pose of the test camera using the RANSAC-based PnP algorithm;

[0011] Step 3) Use the Particle Swarm Optimization (PSO) algorithm to optimize the three-dimensional spatial model, which includes a distortion model and a myopia development model, and then perform registration transformation on the test image accordingly.

[0012] Furthermore, the specific structure of the three-dimensional spatial model is as follows:

[0013] World coordinate system x w - y w - z w The origin is located at the center of the eyeball. w Reference camera C0 is used to acquire reference images. F The test camera c1 is used to acquire test images. F 1. The pose of the reference camera is fixed, while the relative pose between the test camera and the reference camera represents the change in eye pose when acquiring the reference image and the test image. (The reference eye...) E 0 and test eyeballs E 1 corresponds to the reference image and the test image, respectively, which simulate the changes in the shape of the eyeball as myopia develops;

[0014] Eyeball Model

[0015] An ellipsoid is used to approximate the shape of the eyeball, with its center fixed at the origin of the world coordinate system. The three orthogonal semi-axes of the ellipsoid are used. a , b , c ] around the world coordinate axis[ x w , y w , z w ] Rotation angle[ r a , r b , r c The symbol ] represents the eye pose; therefore, the final equation for the eye model is as follows:

[0016] ;

[0017] Where x is a 3D point on the surface of the eyeball model, and Q represents the rotation of the eyeball model; the eyeball model contains a total of 6 parameters: eyeball shape parameters [ a , b , c ] and attitude parameters[ r a , r b , r c ];

[0018] Camera model

[0019] The camera model consists of two parts: an intrinsic parameter matrix and an extrinsic parameter matrix. The intrinsic parameter matrix K depends on the camera itself, and its specific form is as follows:

[0020] ;

[0021] in, α x and α y This indicates the focal length in pixels. u 0 and v 0 represents the coordinates of the image center in the pixel coordinate system. γ This represents the skewness coefficient between the coordinate axes, ideally zero. α x and α y The final calculation formula is as follows:

[0022] ;

[0023] in, r This represents the radius of the retinal region in an image, expressed in pixels. l and k These represent the distance from the lens to the cornea and the field of view, respectively. Specific values ​​can be obtained from the specifications of the fundus camera. ρ This represents the radius of the initial spherical eyeball model, with a value of 12mm;

[0024] The extrinsic parameter matrix {R, t} reflects the position and orientation of the fundus camera in the world coordinate system, including the rotation matrix R and the translation vector t, as shown in the following formula:

[0025] ;

[0026] in,[ r θ , r ϕ , r ω The ] represents the angle of rotation of the camera coordinate system relative to the world coordinate system. t x , t y , t z ] indicates the position of the world coordinate system origin in the camera coordinate system; the unknown parameters in the camera model include the 6 parameters in the extrinsic parameter matrix, [ r θ , r ϕ , r ω ]and[ tx , t y , t z ];

[0027] Distortion model

[0028] During the acquisition of retinal images, the basic form of the distortion model is a fourth-order radial distortion model, with the specific formula as follows:

[0029] ;

[0030] in,[ x d , y d [] represents the distorted coordinates of a point (x, y) on the normalized image plane. k 1, k [2] represents the radial distortion coefficient, with different values ​​corresponding to different degrees of distortion. The parameters in the distortion model consist of two radial distortion coefficients. k 1, k 2];

[0031] Mapping from 3D retinal points to 2D image points

[0032] The specific mapping process of fundus camera imaging after introducing the distortion model is as follows:

[0033] Step 1) Using the extrinsic parameter matrix {R, t} of the fundus camera, the 3D retinal point P in the world coordinate system is... w = [ x w , y w , z w ] T Transform to the camera coordinate system to obtain P c = [ x c , y c , z c ] T ,Right now

[0034] ;

[0035] Step 2) P is converted using the following formula c Transforming to the normalized image plane, we get P1 = [ x , y , 1] T ;

[0036] ;

[0037] Step 3) In the normalized image plane, calculate the coordinates of point P after P1 has undergone fourth-order radial distortion using formula (5). 1_d = [ x d , y d , 1] T ;

[0038] Step 4) P is obtained through the intrinsic parameter matrix K of the fundus camera. 1_d After transforming from the physical image coordinate system to the pixel coordinate system, we obtain the 2D image point P. uv_d ,Right now

[0039] ;

[0040] Mapping from 2D image points to 3D retinal points

[0041] First, consider point P in the 2D image. uv_d Perform distortion correction to obtain a distortion-free 2D image point P. uv Then, the mapping of 3D retinal points is performed. The mapping process from 2D image points to 3D retinal points involves solving for the 2D image point P. uv Passing through the camera center c c The process of finding the intersection point of the ray formula (9) and the eyeball model formula (1) in the posterior hemisphere, where the intersection point is the corresponding 3D retinal point P. w ;

[0042] ;

[0043] Where λ is a coefficient to be solved, and P is the projection matrix of the camera;

[0044] Myopia Development Model

[0045] Two eyeballs are set up, corresponding to the reference image and the test image respectively. The a and b half-axis lengths and pose parameters of the two eyeballs are completely identical. Only the c half-axis length from the center of the eyeball to the fundus is different, and they can be optimized separately.

[0046] Fundus surface S 0 ( S 1) Points on p 1 ( p 2) The center point of the anterior surface of the eyeball p The ray formed by 0 is on another fundus surface S 1 ( S Intersection on 0) p 2 ( p 1) That is, the corresponding point after mapping.

[0047] Furthermore, the parameters of each component in the 3D spatial model are initialized as follows:

[0048] 3.1 The reference eyeball is set to have a radius that is not rotated relative to the world coordinate system. ρ = A sphere with a diameter of 12mm;

[0049] 3.2 Test the eyeball, setting the posterior axis length from the center of the eyeball to the fundus as... c = 12.5mm, and the other parameters are the same as the reference eyeball;

[0050] 3.3 The coefficients of the distortion model use the calibration results as initial values, i.e. k 1 = -0.5623, k 2 = 0.3317;

[0051] 3.4 The reference camera is fixed in the world coordinate system with coordinates c0 = [0, 0, -]. δ ] T = [0, 0, -( l + ρ )] T Furthermore, its orientation remains unchanged relative to the world coordinate system;

[0052] 3.5 The extrinsic parameter matrix {R1, t1} of the test camera was estimated by solving the PnP problem. The solution scheme used was the RANSAC algorithm. This algorithm estimates the camera pose using a set of 2D-3D corresponding points and the camera's projection matrix P. The set of 2D-3D corresponding points consists of 3D points mapped from the matched 2D keypoints on the reference image to a spherical fundus with a radius of 12mm by the reference camera, and the corresponding matched 2D keypoints on the test image.

[0053] Furthermore, the parameters that need to be optimized in the 3D spatial model include: eyeball shape [ a , b , c , c ʹ], eyeball posture [ r a , r b , r c ], distortion coefficient [ k 1, k 2], testing camera pose [ r θ , r ϕ , r ω ,t x , t y , t z ];Will F 0 and F Matched in 1 n Mapping the key points onto the same fundus surface yields... n Group corresponding point {q i , p i Their distance in 3D space is:

[0054] ;

[0055] Therefore, minimizing the distance between all corresponding 3D points is our goal in optimizing the parameters of the 3D spatial model. To improve the robustness of the algorithm to mismatched keypoints, we will... n Sort the 3D corresponding points in ascending order by distance, and only process the first few points. Summing the values:

[0056] ;

[0057] in, Indicates the first position after ascending order. j Distance between corresponding 3D points Let represent the sequence of parameters to be optimized in the 3D spatial model. The objective function shown in formula (11) is optimized using the particle swarm optimization algorithm.

[0058] Furthermore, the image registration and transformation data in step 3) is output in two ways: the first is the reference image. F 0 / Test Image F 1. The reference / test camera image is mapped onto the reference / test fundus, outputting the corresponding 3D point coordinates and RGB values ​​to generate a 3D model of the retinal surface; the second method is to use test images. F 1. The system can sequentially pass through the test camera, test fundus, reference fundus, and reference camera to output registered 2D floating-point coordinates and RGB values, which are then used to generate bilinear interpolation. F 1. The image after registration transformation.

[0059] Compared with existing technologies, the beneficial effects of this invention are as follows: The Retinal Image Registration for Myopia Development (RIRMD) method of this invention, based on REMPE, introduces a distortion model of the fundus camera during retinal imaging into the three-dimensional spatial model. Simultaneously, a novel myopia development model is further constructed and introduced into the three-dimensional spatial model to approximately simulate the changes in axial elongation and fundus curvature of the eyeball as myopia develops.

[0060] The method of this invention and several state-of-the-art retinal image registration methods were quantitatively evaluated and compared on public datasets with no myopia development and on our own constructed datasets with obvious myopia development. The results showed that: (1) on public datasets with no myopia development, the method of this invention achieved the best registration accuracy, registration stability and robustness; (2) on datasets with obvious myopia development, the method of this invention achieved the best registration accuracy with a very significant advantage. Attached Figure Description

[0061] Figure 1 A schematic diagram of the three-dimensional spatial model designed for this invention.

[0062] Figure 2 The imaging process of the fundus camera after the introduction of the distortion model.

[0063] Figure 3 This is a schematic diagram of the mapping scheme between corresponding points on two fundus surfaces.

[0064] Figure 4 This is a schematic flowchart of the registration method of the present invention.

[0065] Figure 5 Example of keypoint extraction and matching results.

[0066] Figure 6 This is an example of the ground truth point in the annotations for each image in the FIMD dataset.

[0067] Figure 7 The registration results of RIRMD and REMPE on the FIMD dataset. Detailed Implementation

[0068] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. To facilitate understanding of the retinal image registration method for myopia development according to the present invention, the three-dimensional spatial model constructed by the present invention will be introduced first.

[0069] This invention establishes a three-dimensional spatial model between pairs of retinal images showing the progression of myopia, based on retinal images acquired by a fundus camera and the changes in eyeball shape as myopia develops. For example... Figure 1 As shown, in this model, the world coordinate system ( x w - y w - z w The origin of ) is located at the center of the eyeball (c w The reference camera (c0) is used to acquire reference images. F The test camera (c1) is used to acquire test images. F 1), where the pose of the reference camera is fixed, and the relative pose between the test camera and the reference camera represents the change in eye pose when acquiring the reference image and the test image. Furthermore, the reference eye ( E 0) and test eyeballs ( E 1) Corresponding to the reference image and the test image respectively, they approximately simulate the changes in eyeball shape as myopia progresses. The following section will describe each component of the three-dimensional spatial model constructed above in more detail.

[0070] Eyeball Model

[0071] An ellipsoid is used to approximate the shape of a reference eyeball, with its center fixed at the origin of the world coordinate system. Furthermore, three orthogonal semi-axises of the ellipsoid are employed. a , b , c ] around the world coordinate axis[ x w , y w , z w ] Rotation angle[ r a , r b , r c The symbol ] represents the eye pose. Therefore, the final equation for the eye model is as follows:

[0072] ;

[0073] Where x is a 3D point on the surface of the eyeball model, and Q represents the rotation of the eyeball model. In summary, the eyeball model contains a total of 6 parameters: eyeball shape parameters [ a , b , c ] and attitude parameters[ r a , r b ,r c ].

[0074] Camera model

[0075] An ideal camera model consists of two parts: an intrinsic parameter matrix and an extrinsic parameter matrix. The intrinsic parameter matrix K depends on the camera itself and takes the following form:

[0076] ;

[0077] in, α x and α y This indicates the focal length in pixels. u 0 and v 0 represents the coordinates of the image center in the pixel coordinate system. γ This represents the tilt coefficient between the coordinate axes, ideally zero. It should be noted that we use the method proposed in REMPE to approximate the calculation. α x and α y The final calculation formula is shown here only. A more detailed description can be found in the reference (C. Hernandez-Matas, X. Zabulis and AA Argyros, "REMPE: Registration of Retinal Images Through Eye Modelling and Pose Estimation," in IEEE Journal of Biomedical and Health Informatics, vol. 24, no. 12, pp. 3362-3373, Dec. 2020, doi: 10.1109 / JBHI.2020.2984483.).

[0078] ;

[0079] in, r This represents the radius of the retinal region in an image, expressed in pixels. l and k These represent the distance from the lens to the cornea and the field of view, respectively. Specific values ​​can be obtained from the specifications of the fundus camera. ρ This represents the radius of the initial spherical eyeball model, with a value of 12mm.

[0080] The extrinsic parameter matrix {R, t} reflects the position and orientation of the fundus camera in the world coordinate system, including the rotation matrix R and the translation vector t, as shown in the following formula:

[0081] ;

[0082] in,[ r θ , r ϕ , r ω The ] represents the angle of rotation of the camera coordinate system relative to the world coordinate system. t x , t y , t z [] indicates the position of the origin of the world coordinate system in the camera coordinate system.

[0083] In summary, the unknown parameters in the camera model include the six parameters in the extrinsic parameter matrix. r θ , r ϕ , r ω ]and[ t x , t y , t z This is because the intrinsic parameter matrix is ​​equivalent to a known constant that can be easily obtained for a given retinal image.

[0084] Distortion model

[0085] Due to the concave shape of the retina and the convex shape of the camera lens, distortion often occurs during retinal image acquisition, which is not considered by the REMPE method. In this invention, we fabricated a ceramic checkerboard calibration plate with 20*16 1.0*1.0mm small squares and acquired multiple images of the calibration plate at different positions and orientations using a fundus camera. Then, using the calibration method in the reference (Z. Zhang, "A flexible new technique for camera calibration," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 22, no. 11, pp. 1330-1334, Nov. 2000, doi: 10.1109 / 34.888718.), we determined the basic form of the distortion model to be a fourth-order radial distortion model, with the specific formula as follows:

[0086] ;

[0087] in,[ x d , y d ] represents the coordinates of a point on the normalized image plane. x , y ) Generate distorted coordinate values, [ k 1, k [2] represents the coefficient of radial distortion, and different values ​​correspond to different degrees of distortion.

[0088] In summary, the parameters in the distortion model consist of two radial distortion coefficients. k 1, k 2).

[0089] Mapping from 3D retinal points to 2D image points

[0090] The process of mapping 3D retinal points to 2D image points is essentially the process of a fundus camera imaging the retina. The fundus camera imaging process after introducing a distortion model is as follows: Figure 2 As shown, the specific mapping process is as follows: (1) The 3D retinal point P in the world coordinate system is mapped using the extrinsic parameter matrix {R, t} of the fundus camera. w = [ x w , y w , z w ] T Transform to the camera coordinate system (CCS) to obtain P. c = [ x c , y c , z c ] T ,Right now

[0091] ;

[0092] (2) P is expressed by the following formula c Transforming to the normalized image plane in the Normalized Coordinate System (NCS), we obtain P1 = [ x , y , 1] T ;

[0093] ;

[0094] (3) In the normalized image plane, the coordinates of point P after the fourth-order radial distortion of P1 are calculated by formula (5). 1_d =[ x d , y d , 1] T ;

[0095] (4) P is obtained through the intrinsic parameter matrix K of the fundus camera. 1_d The 2D image point P is obtained by transforming from the Physical Image Coordinate System (PIS) to the Pixel Coordinate System (PCS). uv_d ,Right now

[0096] ;

[0097] Mapping from 2D image points to 3D retinal points

[0098] The mapping process from 2D image points to 3D retinal points can be understood as solving for 2D image point P. uv Passing through the camera center c c The process of the intersection of the formed ray (Formula (9)) with the posterior hemisphere of the eyeball model (Formula (1)) is the process of finding the intersection point, which is the corresponding 3D retinal point P. w It is important to note that after introducing the distortion model, it is necessary to first analyze the 2D image point P. uv_d Perform distortion correction to obtain a distortion-free 2D image point P. uv Only then can the subsequent 3D retinal point mapping work be carried out.

[0099] ;

[0100] Where λ is a coefficient to be solved, and P is the projection matrix of the camera.

[0101] Myopia Development Model

[0102] We approximate a model of myopia development by simulating changes in axial elongation and fundus curvature in three-dimensional space as the eyeball shape evolves with the development of myopia. For example... Figure 1 As shown, we set up two eyeballs, corresponding to the reference image and the test image respectively. The distance between these two eyeballs... a , b The semi-axial length and pose parameters are completely identical, except for the distance from the center of the eyeball to the fundus. c The different half-shaft lengths allow for individual optimization.

[0103] In addition, we have also set up a mapping scheme for corresponding points between two fundus surfaces, such as Figure 3 As shown, fundus surface S 0 ( S 1) Points on p 1 ( p 2) The center point of the anterior surface of the eyeball p The ray formed by 0 is on another fundus surface S 1 ( S Intersection on 0) p 2 ( p 1) That is, the corresponding point after mapping.

[0104] The three-dimensional spatial model constructed by the present invention has been described in detail above. Based on the above three-dimensional spatial model, the image registration method of the present invention will be described below.

[0105] The RIRMD method proposed in this patent achieves registration on a constructed 3D spatial model using corresponding points between retinal image pairs as a medium. The specific process of this framework is as follows: Figure 4 As shown, firstly, SuperPoint and SuperGlue are used to detect, describe, and match key points in the retinal image pairs; then, the initial pose of the test camera is estimated using the RANSAC-based PnP algorithm; finally, the 3D spatial model is optimized using the Particle Swarm Optimization (PSO) algorithm, and the test image is then registered accordingly.

[0106] Regarding the implementation process of the registration framework, this patented method adopts a similar process to REMPE, but there are several differences: 1. The methods used in the key point extraction and matching stages are different; 2. The method of this invention has more parameters in the initialization and optimization stages of the 3D spatial model, mainly due to the addition of parameters after introducing the distortion model and the myopia development model. These will be explained in more detail below.

[0107] Key point extraction and matching

[0108] In this invention, we employ the SuperPoint deep learning framework proposed by DeTone et al. (DeTone, Daniel et al. "SuperPoint: Self-Supervised Interest Point Detection and Description." 2018 IEEE / CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW) (2017): 337-33712.) to detect and describe key points in retinal images. Then, the SuperGlue deep learning framework proposed by Sarlin et al. (P.-E. Sarlin, D. DeTone, T. Malisiewicz and A. Rabinovich, "SuperGlue: Learning Feature Matching with Graph Neural Networks," 2020 IEEE / CVF Conference on Computer Vision and Pattern Recognition (CVPR), Seattle, WA, USA, 2020, pp. 4937-4946, doi: 10.1109 / CVPR42600.2020.00499.) is used to match keypoint sets between retinal image pairs. An example of keypoint extraction and matching results is given here, such as... Figure 5 As shown.

[0109] Initialization of 3D spatial model

[0110] In this patent, we initialize the parameters of each component in the 3D spatial model as follows: (1) The reference eyeball is set to a radius with no rotation relative to the world coordinate system. ρ = 12mm sphere; (2) For the test eyeball, the length of the posterior half axis from the center of the eyeball to the fundus is set to c = 12.5mm, and the other parameters are the same as the reference eyeball; (3) The coefficients of the distortion model use the calibration results as the initial values, that is k 1 = -0.5623, k 2 = 0.3317; (4) The reference camera is fixed in the world coordinate system with coordinates c0 = [0, 0, - δ ] T = [0, 0, -( l + ρ )]T , and its pose has no rotation relative to the world coordinate system; (5) The external parameter matrix {R1, t1} of the test camera is estimated by solving the PnP [5] problem. The solution scheme adopted is the RANSAC [5] algorithm. This algorithm estimates the pose of the camera by using a set of 2D-3D corresponding points and the projection matrix P of the camera. The set of 2D-3D corresponding points is composed of the 3D points on the spherical fundus with a radius of 12mm that are mapped by the reference camera to the 2D key points matched on the reference image and the corresponding 2D key points on the test image.

[0111] Optimization of 3D spatial models

[0112] The parameters that need to be optimized in the 3D spatial model constructed in this invention include: (1) eyeball shape [ a , b , c , c ʹ];(2)Eyeball posture[ r a , r b , r c ];(3)Distortion coefficient[ k 1, k 2]; (4) Test camera pose [ r θ , r ϕ , r ω , t x , t y , t z ].like Figure 1 As shown, F 0 and F Matched in 1 n Mapping the key points onto the same fundus surface yields... n Group corresponding point {q i , p i Their distance in 3D space is:

[0113] ;

[0114] Therefore, minimizing the distance between all corresponding 3D points is our objective in optimizing the parameters of the 3D spatial model. To improve the algorithm's robustness to mismatched keypoints, we will... n Sort the 3D corresponding points in ascending order by distance, and only process the first few points. Summing the values:

[0115] ;

[0116] in, Indicates the first position after ascending order. j Distance between corresponding 3D points This represents the sequence of parameters to be optimized in the 3D spatial model.

[0117] Furthermore, this invention employs the Particle Swarm Optimization (PSO) algorithm to optimize the objective function shown in equation (11). This algorithm utilizes… p Each particle searches for candidate solutions in a given search space. g Subsequent iterations improve the design to achieve optimization.

[0118] Data output

[0119] After obtaining the optimal 3D spatial model, the data can be output in two ways. First, the reference image... F 0 (Test image) F 1) The reference camera (test camera) can be mapped onto a reference fundus (test fundus), outputting the corresponding 3D point coordinates and RGB values, thereby generating a 3D model of the fundus retinal surface. It is worth noting that this patent generates 3D retinal surfaces with different curvatures from two images showing the development of myopia, reflecting the changes in the shape of the fundus retina caused by the development of myopia. Second, the test image... F 1. The system can sequentially pass through the test camera, test fundus, reference fundus, and reference camera to output registered 2D floating-point coordinates and RGB values, which are then used to generate bilinear interpolation. F 1. The image after registration transformation.

[0120] The implementation process of the 3D spatial model and registration framework constructed in this invention has been described in detail above. To verify the effectiveness of this invention, a quantitative evaluation and comparison of the proposed method and several state-of-the-art methods were performed on a retinal image registration dataset with and without myopia development.

[0121] Experimental setup

[0122] Our experiments were conducted on a computer with the following specifications: Intel(R) Core(TM) i7-8700 CPU @ 3.20GHz; 16.0 GB RAM; Windows 10 Professional operating system. Furthermore, the retinal image registration method proposed in this invention, geared towards the development of myopia, was implemented in Python using OpenCV and PyTorch.

[0123] After debugging the SuperPoint and SuperGlue algorithms, the present invention finally uses the following configuration to extract and match key points in retinal images: (1) For the SuperPoint algorithm, the confidence threshold of the key point detector is set to 0.001, and the other parameters are set to default values; (2) For the SuperGlue algorithm, the loaded weight file is of type "outdoor", and the other parameters are set to default values.

[0124] For the hyperparameters of the PSO algorithm used in this invention, unless otherwise specified, we use the following configuration for experiments: (1) number of particles p = 10000, number of iterations g = 300; (2) The search space and maximum search step size of all parameters to be optimized are shown in Table 1; (3) Individual learning factor c 1 = 1.6, group learning factor c 2 = 1.8, inertia weight w = 0.5.

[0125] Table 1. Search space and maximum search step size for the parameters to be optimized.

[0126]

[0127] Dataset

[0128] We conducted experiments on two datasets. The first is the publicly available Fundus Image Registration (FIRE) dataset without myopia development. The second is our own self-constructed Fundus Image Myopia Development (FIMD) dataset with myopia development.

[0129] The FIRE dataset consists of three classes of retinal image pairs, the features and number of which are shown in Table 2. Each image pair is labeled with 10 pairs of ground truth correspondences to quantitatively evaluate the registration accuracy.

[0130] Table 2 Categories of FIRE datasets

[0131]

[0132] The FIMD dataset contains 70 pairs of retinal images showing clear myopia development. Each image pair is labeled with 12 widely distributed and uniformly distributed ground truth points (an example is given). Figure 6 As shown in the figure, it is used to quantitatively evaluate the registration accuracy.

[0133] Quantitative evaluation indicators

[0134] The registration error of each image pair is calculated. n c (FIRE:) n c = 10; FIMD: n c = 12) The average registration error corresponding to the ground truth points is obtained, in pixels. Then, the registration success rate of all image pairs in the statistical dataset is calculated under different error thresholds (from 0 to 25 pixels). (If the registration error of an image pair is less than the given error threshold, then the image pair is considered successfully registered. The percentage of successfully registered image pairs in the statistical dataset is the registration success rate under that error threshold.) Finally, the normalized area under the curve (AUC) is calculated based on the obtained registration success rate curve, which is the quantitative evaluation index used in the experimental evaluation of this invention. This index reflects the overall registration success rate of a method on a dataset, with a value between 0 and 1, and a larger value is better.

[0135] Experimental results on the FIRE dataset

[0136] We quantitatively compared the registration performance of our proposed method RIRMD with several other existing methods on the FIRE dataset without myopia development, including GDB-ICP proposed by Yang et al. in 2007, SURF-PIIFD-RPM proposed by Wang et al. in 2015, GFEMR proposed by Wang et al. in 2019, VOTUS proposed by Motta et al. in 2019, and REMPE proposed by Hernandez-Matas et al. in 2020.

[0137] The experimental results are shown in the table below. It can be seen that the method proposed in this invention is superior to all comparative methods overall, including VOTUS, which REMPE failed to surpass. Furthermore, compared to the most competitive REMPE and VOTUS, the method proposed in this invention... S , P , A The method proposed in this invention achieves optimal or near-optimal (with only a small difference from the optimal result) results on all three types of data, which indicates that the method has stronger stability and robustness.

[0138] In summary, the RIRMD method proposed in this patent achieves the best registration results in retinal image registration tasks where there is no myopia development.

[0139] Table 3 Experimental results on the FIRE dataset

[0140]

[0141] Experimental results on the FIMD dataset

[0142] We quantitatively compared the registration performance of RIRMD, GDB-ICP, SURF-PIIFD-RPM, GFEMR, and REMPE on the FIMD dataset with myopia occurrence and development (this experiment could not be conducted because VOTUS does not have publicly available implementation code or executable program).

[0143] The experimental results are shown in the table below. It can be seen that the method proposed in this invention is significantly superior to all comparative methods, and even compared to the most competitive REMPE, the method proposed in this invention still has a significant advantage. Therefore, the method proposed in this invention, RIRMD, achieves the best registration effect with a significant advantage in retinal image registration tasks involving the occurrence and progression of myopia.

[0144] Table 4 Experimental results on the FIMD dataset

[0145]

[0146] Qualitative comparison results on the FIMD dataset

[0147] We qualitatively compared the registration results of RIRMD and the most competitive REMPE on the FIMD dataset, such as... Figure 7 As shown, the RIRMD method proposed in this invention can perform more accurate registration of retinal image pairs with significant myopia development, whether in the thicker main vessel region or the thinner branch vessel region.

[0148] In summary, the Retinal Image Registration for Myopia Development (RIRD) method proposed in this invention further extends and improves upon the work of REMPE, achieving accurate registration of retinal image pairs with significant myopia development. Furthermore, experiments have verified that the method proposed in this invention exhibits higher accuracy, stability, and robustness compared to state-of-the-art methods in retinal image registration tasks with and without myopia development.

[0149] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A retinal image registration method for the development of myopia, characterized in that: Follow these steps: Step 1) Use the SuperPoint deep learning framework to detect and describe key points in the reference image and the test image respectively, and then use the SuperGlue deep learning framework to match the key point sets between the reference image and the test image. Step 2) The matched 2D key points on the reference image are mapped onto a spherical fundus with a radius of 12mm by the reference camera to obtain 3D points. These 3D points and the matched 2D key points on the test image form a set of 2D-3D corresponding points. Then, the initial pose of the test camera is estimated by using the PnP algorithm based on RANSAC. Step 3) With the goal of minimizing the Euclidean distance between the matched keypoints in 3D space, the particle swarm optimization algorithm (PSO) is used to optimize the initialized 3D spatial model. The 3D spatial model includes a distortion model and a myopia development model. The optimization parameters of the two models include eye shape parameters, eye pose parameters, fourth-order radial distortion coefficients, and test camera pose parameters. Based on these parameters, the test image is registered and transformed.

2. The retinal image registration method for myopia development according to claim 1, characterized in that: The specific structure of the three-dimensional spatial model is as follows: World coordinate system x w - y w - z w The origin is located at the center of the eyeball. w Reference camera C0 is used to acquire reference images. F 0, while test camera c1 is used to acquire test images. F 1. The pose of the reference camera is fixed, while the relative pose between the test camera and the reference camera represents the change in eye pose when acquiring the reference image and the test image. (The reference eye...) E 0 and test eyeballs E 1 corresponds to the reference image and the test image, respectively, which simulate the changes in the shape of the eyeball as myopia develops; Eyeball Model An ellipsoid is used to approximate the shape of the eyeball, with its center fixed at the origin of the world coordinate system. The three orthogonal semi-axes of the ellipsoid are used. a , b , c ] around the world coordinate axis[ x w , y w , z w ] Rotation angle[ r a , r b , r c The ] represents the eye pose, and the final equation for the eye model is as follows: ； Where x is a 3D point on the surface of the eyeball model, and Q represents the rotation of the eyeball model; the eyeball model contains a total of 6 parameters: eyeball shape parameters [ a , b , c ] and attitude parameters[ r a , r b , r c ]; Camera model The camera model consists of two parts: an intrinsic parameter matrix and an extrinsic parameter matrix. The intrinsic parameter matrix K depends on the camera itself, and its specific form is as follows: ； in, α x and α y This indicates the focal length in pixels. u 0 and v 0 represents the coordinates of the image center in the pixel coordinate system. γ Indicates the skewness coefficient between the coordinate axes. α x and α y The calculation formula is as follows: ； in, r This represents the radius of the retinal region in an image, expressed in pixels. l and k These represent the distance from the lens to the cornea and the field of view, respectively. Specific values ​​are obtained from the specifications of the fundus camera. ρ This represents the radius of the initial eyeball model, with a value of 12mm; The extrinsic parameter matrix {R, t} reflects the position and orientation of the fundus camera in the world coordinate system, including the rotation matrix R and the translation vector t, as shown in the following formula: ； in,[ r θ , r ϕ , r ω The ] represents the angle of rotation of the camera coordinate system relative to the world coordinate system. t x , t y , t z ] indicates the position of the world coordinate system origin in the camera coordinate system; the unknown parameters in the camera model include the 6 parameters in the extrinsic parameter matrix, [ r θ , r ϕ , r ω ]and[ t x , t y , t z ]; Distortion model During the acquisition of retinal images, the basic form of the distortion model is a fourth-order radial distortion model, with the specific formula as follows: ； in,[ x d , y d [] represents the distorted coordinates of a point (x, y) on the normalized image plane. k 1, k [2] represents the radial distortion coefficient, with different values ​​corresponding to different degrees of distortion. The parameters in the distortion model consist of two radial distortion coefficients. k 1, k 2]; Mapping from 3D retinal points to 2D image points The specific mapping process of fundus camera imaging after introducing the distortion model is as follows: Step 1) Using the extrinsic parameter matrix {R, t} of the fundus camera, the 3D retinal point P in the world coordinate system is... w = [ x w , y w , z w ] T Transform to the camera coordinate system to obtain P c = [ x c , y c , z c ] T ,Right now ； Step 2) Use formula (7) to convert P c Transforming to the normalized image plane, we get P1 = [ x , y , 1] T ; ； Step 3) In the normalized image plane, calculate the coordinates of point P after P1 has undergone fourth-order radial distortion using formula (5). 1_d =[ x d , y d , 1] T ; Step 4) P is obtained through the intrinsic parameter matrix K of the fundus camera. 1_d After transforming from the physical image coordinate system to the pixel coordinate system, we obtain the 2D image point P. uv_d ,Right now ； Mapping from 2D image points to 3D retinal points First, consider point P in the 2D image. uv_d Perform distortion correction to obtain a distortion-free 2D image point P. uv Then, the mapping of 3D retinal points is performed. The mapping process from 2D image points to 3D retinal points involves solving for the 2D image point P. uv Passing through the camera center c c The process of finding the intersection point of the ray formula (9) and the eyeball model formula (1) in the posterior hemisphere, where the intersection point is the corresponding 3D retinal point P. w ; ； Where λ is a coefficient to be solved, and P is the projection matrix of the camera; Myopia Development Model Two eyeballs are set up, corresponding to the reference image and the test image respectively. The a and b half-axis lengths and pose parameters of the two eyeballs are completely identical. Only the c half-axis length from the center of the eyeball to the fundus is different, and they can be optimized separately. Fundus surface S 0 ( S 1) Points on p 1 ( p 2) The center point of the anterior surface of the eyeball p The ray formed by 0 is on another fundus surface S 1 ( S Intersection on 0) p 2 ( p 1) That is, the corresponding point after mapping.

3. The retinal image registration method for myopia development according to claim 2, characterized in that: The parameters of each component in the 3D spatial model are initialized as follows: 3.1 The reference eyeball is set to have a radius that is not rotated relative to the world coordinate system. ρ A sphere with a diameter of 12mm; 3.2 Test the eyeball, setting the posterior axis length from the center of the eyeball to the fundus as... c = 12.5mm, all other parameters are the same as the reference eyeball; 3.3 The coefficients of the distortion model use the calibration results as initial values, i.e. k 1 = -0.5623, k 2 = 0.3317; 3.4 The reference camera is fixed in the world coordinate system with coordinates c0 = [0, 0, -]. δ ] T = [0, 0, -( l + ρ )] T Furthermore, its orientation remains unchanged relative to the world coordinate system; 3.5 The extrinsic parameter matrix {R1, t1} of the test camera is estimated by solving the PnP problem. The solution scheme adopted is the RANSAC algorithm. This algorithm estimates the camera pose using a set of 2D-3D corresponding points and the camera's projection matrix P. The set of 2D-3D corresponding points consists of 3D points mapped from the matched 2D keypoints on the reference image to a spherical fundus with a radius of 12mm by the reference camera, and the corresponding matched 2D keypoints on the test image.

4. The retinal image registration method for myopia development according to claim 3, characterized in that: The parameters that need to be optimized in the 3D spatial model include: eye shape [ a , b , c , c ʹ], eyeball posture [ r a , r b , r c ], distortion coefficient [ k 1, k 2], testing camera pose [ r θ , r ϕ , r ω , t x , t y , t z ];Will F 0 and F Matched in 1 n Mapping the key points onto the same fundus surface yields... n Group corresponding point {q i , p i Their distance in 3D space is: ； Therefore, minimizing the distance between all corresponding 3D points is our goal in optimizing the parameters of the 3D spatial model. To improve the robustness of the algorithm to mismatched keypoints, we will... n Sort the 3D corresponding points in ascending order by distance, and only process the first few points. Summing the values: ； in, Indicates the first digit after ascending order. j Distance between corresponding 3D points Let represent the sequence of parameters to be optimized in the 3D spatial model. The objective function shown in formula (11) is optimized using the particle swarm optimization algorithm.

5. The retinal image registration method for myopia development according to claim 1, characterized in that: In step 3), the image registration and transformation data is output in two ways: the first is the reference image. F 0 / Test Image F 1. The reference / test camera image is mapped onto the reference / test fundus, outputting the corresponding 3D point coordinates and RGB values ​​to generate a 3D model of the retinal surface; the second method is to use test images. F 1. The system can sequentially pass through the test camera, test fundus, reference fundus, and reference camera to output registered 2D floating-point coordinates and RGB values, which are then used to generate bilinear interpolation. F 1. The image after registration transformation.

Images