Laser speckle fundus image registration method and system

Abstract

The application relates to the technical field of image registration, and discloses a laser speckle fundus image registration method and system, which comprises the following steps: reading a laser speckle fundus video and discretizing the laser speckle fundus video into a sequence of original images, determining a reference frame and a floating frame, generating an enhanced fundus image after spherical projection correction, and obtaining a multi-scale blood vessel skeleton hierarchical expression; then constructing a local blood vessel tree signature vector based on the multi-scale blood vessel skeleton hierarchical expression, calculating a feature matching cost matrix; finally, constructing a topological graph by using a candidate feature point set, screening a set of certain matching point pairs, calculating a spatial geometric transformation matrix, outputting a high signal-to-noise ratio fundus image after floating frame registration is completed, and improving the feature point matching accuracy and the registration success rate of a multi-blood vessel fork point dense area, and eliminating the cascade registration failure phenomenon caused by false matching.

Description

Technical Field

[0001] This invention relates to the field of image registration technology, and more specifically, to a laser speckle fundus image registration method and system. Background Technology

[0002] Precise registration of fundus images is a core pre-treatment technique for the diagnosis of ophthalmic diseases and the analysis of fundus hemodynamics. Laser speckle fundus imaging technology, due to its ability to non-invasively acquire microcirculatory blood flow information in the fundus, has shown unique advantages in the early screening and monitoring of fundus diseases such as diabetic retinopathy and glaucoma. Because of factors such as physiological nystagmus in the subjects and optical distortion of the imaging equipment, the position of blood vessels in each frame of the laser speckle fundus video sequence may be offset. Registration techniques are needed to align the blood vessels in all frames with the reference frame to achieve multi-frame information fusion, thereby improving the image signal-to-noise ratio and providing reliable data support for subsequent medical analysis.

[0003] In the prior art, Chinese patent CN115409689B discloses a registration method and apparatus for multimodal retinal fundus images. It first acquires multimodal reference and target retinal fundus images, extracts blood vessels from the two types of images to obtain corresponding vascular maps, then determines correctly paired feature points through a screening operation to complete registration, performs affine transformation on the target images, and finally outputs the reference image and the transformed image to a display device in a checkerboard pattern, thus realizing the registration of multimodal fundus images. Chinese patent application CN116630377A discloses a dynamic registration quantification method and apparatus based on a vascular skeleton. It first establishes vascular standard values ​​for a dynamic vascular image sequence, then generates a reference skeleton and the skeleton of the target vascular image. Based on the reference skeleton, it constructs the geometric transformation relationship between the target skeleton and the reference skeleton, completes dynamic vascular image registration through rigid transformation, and then quantifies the dynamic changes of blood vessels and generates change curves, ensuring the accuracy of hemodynamic response assessment.

[0004] However, existing technologies still have significant limitations in the registration of laser speckle fundus images: the feature point selection in the existing technology CN115409689B does not construct multi-scale feature representations for the uneven blood vessel thickness in laser speckle images. The same blood vessel bifurcation point will exhibit different morphologies at different image processing scales. When there is a slight angular shift in the image, the relative width of thick and thin blood vessels will change, leading to a decrease in feature detection stability and causing local feature scale mismatch problems. The dynamic registration in the existing technology CN116630377A does not establish a global verification mechanism for blood vessel topology. In areas with dense multiple blood vessel bifurcation points, it is easy to confuse a bifurcation point in the reference frame with multiple neighboring candidate bifurcation points in the floating frame. Once an initial matching error occurs, it will trigger a cascading registration error, which will not only reduce the accuracy of feature point matching, but also lead to registration failure in areas with dense multiple blood vessel bifurcation points. It cannot provide a high-precision registration basis for multi-frame fusion of laser speckle fundus images, and ultimately affect the reliability of fundus blood flow information analysis. Summary of the Invention

[0005] The laser speckle fundus image data used in the implementation of this invention are all derived from medical image datasets with explicit authorization from the subjects, and the data collection process strictly complies with the relevant provisions of the Personal Information Protection Law. All subject identity information has been de-identified, retaining only the image data used for image registration and analysis. The implementation of this invention does not involve any unauthorized collection and processing of personal information.

[0006] To overcome the aforementioned deficiencies of existing technologies, this invention provides a laser speckle fundus image registration method and system. Through steps such as generating multi-scale hierarchical representations of the vascular skeleton, constructing local vascular tree signature vectors, filtering confirmed matching point pairs using a topological graph, and calculating the spatial geometric transformation matrix, laser speckle fundus image registration is achieved. This effectively improves the accuracy of feature point matching, increases the registration success rate in densely populated areas with multiple vascular bifurcation points, and completely eliminates cascade registration failures caused by incorrect matching. Ultimately, it outputs high signal-to-noise ratio fundus images, providing high-quality data support for subsequent fundus medical analysis.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] A laser speckle fundus image registration method includes:

[0009] Read the laser speckle fundus video and discretize it into a raw image sequence. Determine the reference frame and floating frame. Perform spherical projection correction on the reference frame and floating frame to generate enhanced fundus images. Process the enhanced fundus images to generate multi-scale vascular skeleton hierarchical expression.

[0010] Based on the multi-scale hierarchical representation of the vascular skeleton, topological key points are detected to form a candidate feature point set; the geometric properties of the local vascular branches of each topological key point in the candidate feature point set are tracked, a local vascular tree signature vector is constructed based on the geometric properties, and the feature matching cost matrix is ​​calculated based on the local vascular tree signature vector.

[0011] A baseline topology map and a floating topology map are constructed using a candidate feature point set. A set of sure-matching point pairs is selected from the baseline topology map and the floating topology map based on the feature matching cost matrix. A spatial geometric transformation matrix is ​​calculated based on the set of sure-matching point pairs. The floating frame is then registered based on the spatial geometric transformation matrix, and a high signal-to-noise ratio fundus image is output.

[0012] The method for determining the reference frame and the floating frame includes:

[0013] A center-periphery weighted scoring mechanism is applied to each image in the original image sequence to calculate a comprehensive quality score. The image with the highest comprehensive quality score is selected as the reference frame, and the remaining images in the original image sequence are used as floating frames to be registered.

[0014] The method for generating multi-scale hierarchical expressions of the vascular skeleton includes:

[0015] Three sets of Gaussian-matched filters, coarse-scale, medium-scale, and fine-scale, were set up to perform convolution and threshold segmentation on the enhanced fundus image, generating binary maps of main blood vessels, binary maps of the entire blood vessel network, and texture maps of microvessels.

[0016] Iterative erosion skeletonization was performed on the binary map of the main blood vessel, the binary map of the whole blood vessel network, and the microvascular texture map to obtain coarse-scale skeleton, medium-scale skeleton, and fine-scale skeleton. The coarse-scale skeleton, medium-scale skeleton, and fine-scale skeleton were combined to form a multi-scale hierarchical expression of the blood vessel skeleton.

[0017] The candidate feature point set includes a reference frame candidate feature point set and a floating frame candidate feature point set;

[0018] The method for forming a candidate feature point set includes:

[0019] Neighborhood connectivity analysis was performed on pixels on the mesoscale skeleton in the multi-scale vascular skeleton hierarchical representation to identify bifurcation points and intersection points. All bifurcation points and intersection points in the mesoscale skeleton were collected as topological key points. The topological key points of the baseline frame were aggregated into a candidate feature point set for the baseline frame, and the topological key points of the floating frame were aggregated into a candidate feature point set for the floating frame.

[0020] The method for tracking the geometric properties of local vascular branches at each topological key point in the candidate feature point set includes:

[0021] Centered on each topological keypoint in the candidate feature point set, the branch direction and branch length are traced along the mesoscale skeleton, and mapped back to the binary map of the whole vascular network to measure the branch width, thus obtaining geometric attributes including branch direction angle, branch length and branch width.

[0022] The method for constructing a local vascular tree signature vector based on geometric properties includes:

[0023] The angle between adjacent branches is calculated based on the branch direction angle in the geometric attributes, and the angle between adjacent branches is used as a rotation invariant angular feature. The branch length and branch width are normalized to obtain normalized length and normalized width, and the normalized length and normalized width are used as scale invariant geometric features. The rotation invariant angular features and scale invariant geometric features are combined to generate the local vascular tree signature vector of each topological key point.

[0024] The method for calculating the feature matching cost matrix based on the local vascular tree signature vector includes:

[0025] Topological keypoints in the candidate feature point set of the baseline frame are selected and paired with topological keypoints in the candidate feature point set of the floating frame to form multiple point pair combinations. The invariant distance between the local vascular tree signature vector of the topological keypoint in the candidate feature point set of the baseline frame and the local vascular tree signature vector of the topological keypoint in the candidate feature point set of the floating frame is calculated in each point pair combination. All point pair combinations are traversed to generate a feature matching cost matrix.

[0026] The method for constructing the baseline topology graph and the floating topology graph includes:

[0027] Each topological keypoint in the candidate feature point set of the baseline frame and the candidate feature point set of the floating frame is used as a vertex of the topological graph. Edges are established between the topological keypoints that are directly connected on the mesoscale skeleton to construct the baseline topological graph and the floating topological graph respectively.

[0028] The method for filtering the set of sure-matching point pairs includes:

[0029] Based on the feature matching cost matrix, a set of candidate matching point pairs is selected. For each candidate matching point pair in the set, a second-order proximity constraint verification is performed. Candidate matching point pairs that fail the constraint verification are eliminated to obtain a set of sure matching point pairs.

[0030] A laser speckle fundus image registration system is provided for implementing the aforementioned laser speckle fundus image registration method. The system comprises:

[0031] Multi-scale skeleton generation module: used to read laser speckle fundus video and discretize it into a raw image sequence, determine the reference frame and floating frame, perform spherical projection correction on the reference frame and floating frame to generate enhanced fundus image, process the enhanced fundus image, and generate multi-scale vascular skeleton hierarchical expression.

[0032] The signature vector construction module is used to detect topological key points based on multi-scale vascular skeleton hierarchical expression and form a candidate feature point set; track the geometric properties of local vascular branches of each topological key point in the candidate feature point set, construct local vascular tree signature vectors based on geometric properties, and calculate the feature matching cost matrix based on local vascular tree signature vectors.

[0033] Image registration module: Constructs a baseline topology map and a floating topology map using candidate feature point sets. Selects a set of sure-matching point pairs from the baseline topology map and the floating topology map based on the feature matching cost matrix. Calculates the spatial geometric transformation matrix based on the set of sure-matching point pairs. Registers the floating frames based on the spatial geometric transformation matrix and outputs a high signal-to-noise ratio fundus image.

[0034] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0035] This invention effectively resolves scale ambiguities of feature points with varying vessel thicknesses by generating a multi-scale hierarchical representation of the vascular skeleton, laying a stable multi-scale geometric foundation for subsequent feature extraction. The local vascular tree signature vector constructed based on this hierarchical representation accurately encodes the local topology and geometric attributes of feature points, achieving preliminary accurate matching of feature points when combined with a feature matching cost matrix. Furthermore, by constructing a topology graph and filtering a set of confirmed matching points, the topological ambiguity problem in densely populated areas of multi-vascular bifurcation points is resolved, ensuring high purity of the matching point pairs. Finally, registration is completed through a spatial geometric transformation matrix, outputting a high signal-to-noise ratio image. This improves the accuracy of feature point matching and the registration effectiveness in densely populated areas of multi-vascular bifurcation points, fundamentally eliminating the cascade registration failure problem caused by incorrect matching, and providing high-quality registration data support for subsequent medical analysis of laser speckle fundus images. Attached Figure Description

[0036] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0037] Figure 1 A flowchart of a laser speckle fundus image registration method provided in an embodiment of the present invention;

[0038] Figure 2A flowchart illustrating the center-periphery weighted scoring mechanism provided in this embodiment of the invention;

[0039] Figure 3 This is a schematic diagram of the multi-scale vascular skeleton hierarchy provided in an embodiment of the present invention;

[0040] Figure 4 This is a schematic diagram illustrating the principle of neighborhood connectivity analysis and topological key point identification provided in an embodiment of the present invention;

[0041] Figure 5 A functional block diagram of a laser speckle fundus image registration system provided in an embodiment of the present invention. Detailed Implementation

[0042] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0043] Example 1

[0044] Please see Figure 1 As shown, this embodiment provides a laser speckle fundus image registration method, including:

[0045] Step S10: Read the laser speckle fundus video and discretize it into an original image sequence. Determine the reference frame and floating frame. Perform spherical projection correction on the reference frame and floating frame to generate an enhanced fundus image. Process the enhanced fundus image to generate a multi-scale vascular skeleton hierarchical expression.

[0046] Further, step S10 includes:

[0047] Step S11: Read the laser speckle fundus video and discretize it into an original image sequence. Apply a center-periphery weighted scoring mechanism to each image in the original image sequence to calculate a comprehensive quality score. Select the image with the highest comprehensive quality score as the reference frame and use the remaining images in the original image sequence as floating frames to be registered.

[0048] Specifically, the laser speckle fundus video is dynamic fundus image data acquired using laser speckle imaging technology. Its duration is approximately T seconds, with T typically ranging from 3 to 8 seconds. The video frame rate is F frames per second, with F typically ranging from 30 to 80 frames per second. The resulting discretized original image sequence contains N grayscale fundus images, where N is the product of T and F. The fundus images in the original image sequence are grayscale images. Each grayscale fundus image exhibits the inherent low signal-to-noise ratio characteristic of laser speckle imaging. Furthermore, the unavoidable physiological tremors of the examined eye during acquisition result in motion blur in some images, decreased sharpness of blood vessel edges, and complete loss of information in localized areas due to eyelash occlusion caused by blinking. The optical characteristics of fundus imaging concentrate the light source energy in the central region of the image, resulting in sufficient brightness and high blood vessel contrast in the central region, while the peripheral regions far from the center show significant brightness attenuation and blurred blood vessel details.

[0049] See Figure 2 The implementation process of the center-periphery weighted scoring mechanism is as follows: Each grayscale image of the fundus is divided into a central region and a peripheral region. The central region is defined as a circular region with the geometric center of the image as the center and the radius of the circle obtained by multiplying the length of the short side of the image by the radius coefficient Rc of the central region as the radius. The method of determining Rc is to calibrate it according to the optical parameters of the laser speckle imaging device and the typical position range of the optic disc. For example, the range of Rc values ​​is 0.2 to 0.4. The peripheral region is defined as the annular region outside the central region.

[0050] The gray-level histogram entropy and edge sharpness score are calculated for the central region and for the surrounding region, respectively. The gray-level histogram entropy is calculated using the information entropy formula, which statistically represents the gray-level values ​​of pixels within the region as a histogram. The negative logarithmic weighted sum of the probabilities of each gray level is then calculated. A higher entropy value indicates a more uniform gray-level distribution and richer image information. The edge sharpness score is calculated using a gradient magnitude statistical method. The Sobel operator is applied to the pixels within the region to extract the horizontal and vertical gradients. The mean of the gradient magnitudes is used as the edge sharpness score. A higher score indicates sharper edge transitions and higher image clarity. The grayscale histogram entropy and edge sharpness score are normalized. The grayscale histogram entropy is normalized by dividing the original entropy by the theoretical maximum entropy. For example, for an eight-bit grayscale image, the theoretical maximum entropy is 8 bits, and the normalized entropy ranges from 0 to 1. The edge sharpness score is normalized by dividing the original score by the current image's global maximum gradient magnitude, and the normalized edge sharpness ranges from 0 to 1. The formula for calculating the overall quality score Q is Q=Wc×(H'c+E'c)+Wp×(H'p+E'p), where H'c is the gray-level histogram entropy value of the normalized central region, E'c is the edge sharpness score of the normalized central region, H'p is the gray-level histogram entropy value of the normalized peripheral region, E'p is the edge sharpness score of the normalized peripheral region, Wc is the weight of the central region, and Wp is the weight of the peripheral region. Wc>Wp because the light source is concentrated in the center, making the image quality of the central region contribute more to the extraction of blood vessel features, while the peripheral region has a relatively smaller impact on the overall registration effect due to insufficient brightness and large distortion. For example, Wc is 0.7 and Wp is 0.3. The center-periphery weighted scoring mechanism, rather than uniform weighted scoring, is adopted to match the quality assessment results with the optical energy distribution characteristics of laser speckle imaging. Images with high vascular clarity and low distortion in the central region receive higher scores, while low quality in the peripheral region due to insufficient illumination does not excessively lower the overall score. The selected reference frame thus has higher vascular recognizability in the optic disc and its surrounding specific radius region, providing a high-quality geometric basis for subsequent vascular skeleton extraction and feature point detection, and avoiding the decrease in registration accuracy caused by poor reference frame quality.

[0051] The reference frame is the image with the highest overall quality score in the original image sequence, serving as the spatial reference coordinate system for registration of all floating frames. Floating frames are the remaining images in the original image sequence excluding the reference frame, and they need to be aligned with the reference frame through spatial transformation. Selecting the frame with the highest overall quality score, rather than a frame in the middle of the sequence or a frame at a specific location, as the reference frame ensures that the reference frame has optimal vessel clarity and minimal motion blur in the central region. This results in higher geometric stability for the vascular skeleton and topological keypoints extracted subsequently from the reference frame, smaller measurement errors in feature point positions and branch attributes, and makes the quality of the candidate feature point set from the reference frame superior to that extracted from any frame, thereby improving the reliability of feature matching.

[0052] Step S12: Construct a virtual hemispherical retinal model, use the back projection algorithm to map the reference frame and floating frame onto the surface of the virtual hemispherical retinal model and re-unfold it, perform contrast-limited histogram equalization on the unfolded image to obtain the enhanced reference frame and enhanced floating frame respectively, and use the enhanced reference frame and enhanced floating frame as the enhanced fundus image.

[0053] The construction process of the virtual hemispherical retinal model includes: Anatomically, the retina of the eyeball is an approximately spherical crown-shaped curved surface. When the laser speckle imaging system projects this surface onto a planar image sensor, perspective distortion occurs. This manifests as the spacing and length of blood vessels in the central region of the image being close to the actual anatomical dimensions, while the spacing and length of blood vessels in the peripheral regions far from the center are compressed. The virtual hemispherical retinal model is a mathematical surface established based on the average axial length of the human eye and the radius of curvature of the retina. The implementation process of the back projection algorithm is as follows: The optical center O is located based on the imaging characteristics of the highest brightness at the image center. center O center The coordinates are obtained by calculating the centroid of the image brightness distribution; with O center Establish a mapping relationship between the coordinates of the planar image and the surface coordinates of the virtual hemispherical retina model, with the projection origin as the origin; for any pixel P in the planar image... plane According to its relationship with O center The distance is calculated to the latitude and longitude coordinates P on the virtual hemispherical retinal model surface corresponding to the point. sphere ; For P sphere The corrected planar coordinates P are obtained by applying equidistant cylindrical projection. correctedThis restores the compressed spacing and length of blood vessels in the surrounding areas. Using spherical projection correction instead of directly using the original image aims to eliminate perspective distortion caused by the bowl-shaped structure of the eyeball, ensuring that the geometric dimensions of blood vessels in each region of the image maintain a consistent proportional relationship with the actual anatomical dimensions. This results in the extracted vascular skeleton having uniform scale characteristics in space, avoiding scale ambiguity caused by the same blood vessel exhibiting different widths in the image center and surrounding areas. This allows the multi-scale filter in step S13 to stably detect blood vessels of the same diameter across the entire image.

[0054] The implementation process of contrast-limited histogram equalization is as follows: The corrected planar image is divided into several non-overlapping sub-blocks. For each sub-block, a gray-level histogram is calculated independently and histogram equalization is performed to stretch the gray-level distribution within the sub-block to the full range. To avoid excessive amplification of noise during histogram equalization, a contrast limit threshold Cclip is set. When the frequency of a certain gray level in the sub-block histogram exceeds Cclip, the excess is evenly distributed to other gray levels. Cclip is determined by calibrating according to the typical noise level of the laser speckle image. An exemplary Cclip value range is 2 to 4. The equalization results of adjacent sub-blocks are smoothly transitioned using bilinear interpolation to eliminate gray-level jumps at the sub-block boundaries. The use of contrast-limited histogram equalization instead of global histogram equalization aims to enhance the contrast between blood vessels and the background in local areas while suppressing noise amplification. This results in sharper blood vessel edges and a more significant grayscale difference between blood vessels and the background in the enhanced fundus image. This provides a higher signal-to-noise ratio input for subsequent Gaussian-matched filter blood vessel detection, increasing the distinguishability between blood vessel and background response values ​​and reducing the false positive and false negative rates during threshold segmentation. The enhanced reference frame is the image obtained after performing spherical projection correction and contrast-limited histogram equalization on the reference frame. The enhanced floating frame is the image obtained after performing the same processing on the floating frame. The enhanced reference frame and the enhanced floating frame are collectively referred to as the enhanced fundus image.

[0055] Step S13: Set three sets of Gaussian-matched filters for coarse scale, medium scale and fine scale, and perform convolution and threshold segmentation on the enhanced fundus image respectively to generate binary map of main blood vessel, binary map of whole blood vessel network and microvascular texture map;

[0056] The Gaussian-matched filter is a directional filter specifically designed for detecting linear structures. Its design principle is based on the prior knowledge that the gray-level distribution of blood vessel cross-sections approximates a Gaussian shape. The filter kernel function exhibits a Gaussian distribution perpendicular to the blood vessel's orientation and a uniform distribution parallel to it. It produces a high response value when convolved with the gray-level profile of the blood vessel and a low response value when convolved with the background texture. The sigma values ​​of the three Gaussian-matched filters are set as σc, σm, and σf, corresponding to coarse, medium, and fine scales, respectively. σc is determined based on the typical diameter range of the main retinal vessels, obtained through population statistics in ophthalmic imaging. An exemplary σc value ranges from 6 to 10 pixels, corresponding to a main vessel diameter range of 15 to 25 pixels. σm is determined based on the typical diameter range of retinal branch vessels, with an exemplary σm value range of 3 to 5 pixels, corresponding to a branch vessel diameter range of 7 to 14 pixels. The method for determining σf is based on the typical diameter range of retinal microvessels. For example, the σf value ranges from 1 to 2 pixels, corresponding to a microvessel diameter range of 1 to 6 pixels. Each filter group contains multiple kernel functions in different directions, with N directions in total. dir The method for determining N is based on the isotropic requirement of blood vessel orientation, for example, N. dir The values ​​are set in 12 directions, with adjacent directions spaced 15 degrees apart. Using three sets of filters with different sigma values ​​instead of a single-scale filter is particularly necessary for laser speckle fundus images: the inherent low signal-to-noise ratio and speckle noise in laser speckle imaging blur vessel edges and severely reduce contrast. If a single-scale filter is used, the coarse-scale filter, while detecting the main vessel, will completely lose information about fine vessels, resulting in an incomplete vascular network topology; the fine-scale filter, while detecting microvessels, will misidentify the edges of the main vessel as multiple parallel vessels, generating numerous false bifurcation points; the meso-scale filter, under the interference of laser speckle noise, will cause breaks at vessel boundaries, disrupting vascular connectivity. The multi-scale hierarchical extraction strategy allows vessels of different diameters to be detected independently at their respective optimal scale levels. The filter parameters at each scale level are optimized for vessels within a specific diameter range and the corresponding noise level. The large kernel size of the coarse-scale filter smooths laser speckle noise while maintaining the integrity of the main vessel; the fine-scale filter focuses on microvessels without being affected by the edges of the main vessel; and the meso-scale filter achieves a balance between noise suppression and detail preservation. This multi-scale hierarchical extraction is not simply a matter of information overlay, but a necessary prerequisite for ensuring the correctness of vascular network topology under high-noise conditions caused by laser speckle.

[0057] The generation process of the main vessel binary map is as follows: The kernel functions in each direction of the coarse-scale filter bank are convolved with the enhanced fundus image. For each pixel location, the maximum value among all directional response values ​​is taken as the vessel response value at that location, forming a coarse-scale response map; a high threshold T is set. high T high The method for determining the threshold is based on the statistical distribution of response values ​​from the coarse-scale response map. The high percentile of the response value distribution is taken as the threshold; for example, the 95th percentile is used. Response values ​​exceeding T are considered thresholds. high The pixels in the first position are labeled as blood vessel pixels, and the remaining pixels are labeled as background pixels, forming an initial binary image. A morphological closing operation is then performed on the initial binary image to fill the voids inside the blood vessels. The structuring element of the closing operation is a circle with a radius of R. close R close The value is set based on half the width of the main blood vessel; after closing the operation, a binary image of the main blood vessel is obtained. The generation process of the binary image of the whole blood vessel network is as follows: a mesoscale filter bank is used to perform convolution operation with the enhanced fundus image, and the maximum directional response value is taken to form a mesoscale response map; a mesoscale threshold T is set. mid T mid The method for determining this is also based on the percentiles of the response value distribution; for example, the percentile is the 85th percentile. Response values ​​exceeding T... mid The pixels are labeled as blood vessel pixels, forming a binary map of the entire blood vessel network, which includes main vessels and branch vessels. The generation process of the microvascular texture map is as follows: a fine-scale filter bank is used to perform convolution operation with the enhanced fundus image, and the maximum response value in the direction is taken to form a fine-scale response map; a local adaptive thresholding method is used to segment the fine-scale response map. The local adaptive threshold is calculated by taking 50%-70% of the average response value in the neighborhood of each pixel position as the threshold; pixels with response values ​​exceeding the local adaptive threshold are labeled as blood vessel pixels, forming a microvascular texture map. Different threshold strategies are used to segment three scale levels to match the signal-to-noise ratio differences of blood vessels at each scale. The coarse-scale main vessels have a high signal-to-noise ratio and are suitable for using a high fixed threshold, while the fine-scale microvascular vessels have a low signal-to-noise ratio and are greatly affected by the local background and are suitable for using an adaptive threshold. This makes the binary maps of the three scale levels have low false detection and false negative rates. The main vessel binary map does not contain isolated points that are falsely detected due to noise, and the microvascular texture map does not lose real blood vessels due to excessively high fixed thresholds.

[0058] Step S14: Perform iterative erosion skeletonization processing on the main vessel binary map, the whole vessel network binary map, and the micro-vessel texture map respectively to obtain coarse-scale skeleton, medium-scale skeleton, and fine-scale skeleton. Combine the coarse-scale skeleton, medium-scale skeleton, and fine-scale skeleton to form a multi-scale vascular skeleton hierarchical expression.

[0059] The skeletonization process is as follows: Skeletonization is a morphological image processing operation that peels the target region in a binary image layer by layer down to the center line, which is only a single pixel wide, while maintaining the topological connectivity of the target. Iterative erosion skeletonization is performed on the main blood vessel binary image. The iterative erosion uses an 8-connected structuring element. Each iteration deletes boundary pixels that meet the deletion criteria: the deletion of the pixel will not cause a break or hole in the target region. The iteration continues until no pixels meet the deletion criteria, resulting in a coarse-scale skeleton. The same skeletonization process is then performed on the full blood vessel network binary image and the micro-vessel texture image, respectively, to obtain the medium-scale skeleton and the fine-scale skeleton. (See also...) Figure 3 The coarse-scale skeleton only contains the centerline of the main blood vessel, with fewer but stable bifurcation and intersection points. The meso-scale skeleton contains the centerlines of both the main and branch blood vessels, with a moderate number of bifurcation and intersection points covering the main vascular network around the optic disc. The fine-scale skeleton contains the centerlines of all detectable blood vessels, with numerous bifurcation and intersection points, but some locations are affected by noise. Iterative erosion skeletonization, rather than morphological refinement, is used to strictly maintain the topological connectivity of the vascular network. The three branches at the bifurcation points remain connected to the same pixel after skeletonization, ensuring that the topological key points after skeletonization accurately correspond to the anatomical bifurcation positions of the original blood vessels. This provides accurate topological input for the neighborhood connectivity analysis in the subsequent step S20, avoiding the loss of bifurcation points or the generation of pseudo-bifurcation points due to skeleton breakage.

[0060] The multi-scale hierarchical representation of the vascular skeleton is combined as follows: coarse-scale, meso-scale, and fine-scale skeletons are stored as three independent levels. These three levels share the same image coordinate system but contain vascular centerline information of varying densities. The high noise level and low contrast of laser speckle images pose a serious risk of topological errors to single-scale skeletonization. If single-scale skeleton extraction is performed only on the original noisy image, laser speckle noise will cause instability in vascular edge detection, resulting in numerous breakpoints after skeletonization. Real vascular bifurcations are incorrectly identified as multiple independent endpoints. Simultaneously, noise-induced pseudo-vascular responses will generate false bifurcation points, severely disrupting the vascular network's topology. Multi-scale hierarchical representation, by independently extracting skeletons at different noise suppression levels, ensures that the topology of each scale level remains correct under the corresponding signal-to-noise ratio: the coarse-scale skeleton ensures the integrity of the main vessel topology through strong noise suppression; the meso-scale skeleton maintains the connectivity of branch vessels under moderate noise suppression; and the fine-scale skeleton, although affected by noise, retains microvascular information. This hierarchical strategy is not a simple superposition of multi-scale information but a necessary prerequisite for ensuring topological correctness and reliable feature point extraction in the high-noise environment of laser speckle. The hierarchical inclusion relationships between different levels further validate the consistency of the topological structure, providing a reliable geometric basis for subsequent feature extraction and matching.

[0061] The multi-scale vascular skeleton hierarchical representation output in step S10 provides the geometric basis for the topological keypoint detection in step S20. Without the multi-scale hierarchical extraction in step S10, step S20 would face the dilemma of simultaneously detecting both coarse and fine vascular bifurcation points on a single-scale skeleton. Coarse vascular bifurcation points, due to their large vessel width, would generate multiple pseudo bifurcation points after skeletonization, while fine vascular bifurcation points, due to their low signal-to-noise ratio, would result in breaks after skeletonization. The mixture of these two types of bifurcation points would lead to a large number of unstable feature points in the candidate feature point set, resulting in a large number of mismatched candidates in the feature matching cost matrix of step S24. The spherical projection correction in step S10 provides the coordinate basis after distortion correction for the spatial geometric transformation matrix calculation in step S30. Without the spherical projection correction in step S12, the transformation matrix calculated in step S30 would need to simultaneously fit rigid body transformation and perspective distortion. The increased degrees of freedom of the transformation model would lead to a decrease in fitting stability under a limited number of matching point pairs, and the registration accuracy would decrease accordingly. The quality of the reference frame selected in step S10 determines the quality of the candidate feature point set of the reference frame extracted in step S20. If the quality of the reference frame is poor, the positional error of the topological key points in the candidate feature point set of the reference frame will be large, and the measurement error of the branch geometric attributes will be large. This will cause the local blood vessel tree signature vector generated in step S23 to fail to accurately describe the local topological structure of the feature points, making it difficult for the graph matching consistency verification in step S32 to distinguish between correct matching and structural mismatch.

[0062] In step S10, the center-periphery weighted scoring mechanism makes the reference frame selection biased towards images with clear blood vessels in the central region. The clear blood vessel edges make the peak value of the Gaussian-matched filter response value sharper. The sharp peak value of the response value makes the blood vessel boundary position after threshold segmentation more accurate. The accurate blood vessel boundary position makes the centerline position after skeletonization closer to the anatomical center of the blood vessel. The skeleton close to the anatomical center makes the bifurcation point position detected by neighborhood connectivity analysis more stable. The stable bifurcation point position makes the measurement error of the angle feature and length feature of the local blood vessel tree signature vector smaller. The signature vector with smaller measurement error makes the distance value of the correctly matched point pair in the feature matching cost matrix smaller and the distance value of the mismatched point pair relatively larger. The increased distance difference between the correct match and the mismatch makes it easier to distinguish between the two types of matching in the graph matching consistency verification in step S32, thereby improving the accuracy of feature point matching. Spherical projection correction makes the geometric size ratio of blood vessels in different regions of the image uniform. The uniform size ratio ensures that blood vessels of the same diameter are detected by filters at the same scale level across the entire image. Blood vessels detected at the same scale level have consistent topological properties after skeletonization. The consistent topological properties make the local blood vessel tree signature vectors of corresponding blood vessel bifurcation points in the reference frame and the floating frame more similar. The high similarity of the signature vectors makes the distance value of the correct matching point pair in the feature matching cost matrix smaller, thereby improving the registration success rate of densely populated areas with multiple blood vessel bifurcation points. Multi-scale vascular skeleton hierarchical representation separates coarse and fine vessel bifurcation points into different scale levels. The separated bifurcation points have stable local topological structures within their respective scale levels. The stable local topological structures ensure that the neighbor relationships in the topology graph constructed in step S31 accurately reflect the physical connections of the vascular network. The accurate neighbor relationships enable the second-order proximity constraint verification in step S32 to effectively identify structural mismatches. After the structural mismatches are eliminated, the purity of the set of sure-match point pairs is improved. The pure set of sure-match point pairs results in a smaller fitting residual for the spatial geometric transformation matrix calculated in step S33. The transformation matrix with a smaller fitting residual results in a higher degree of vascular overlap between the registered floating frame and the reference frame. Thus, the cascade registration failure phenomenon caused by mismatches is eliminated.

[0063] Step S20: Construct a local vascular tree signature vector based on the multi-scale vascular skeleton hierarchical expression, and calculate the feature matching cost matrix based on the local vascular tree signature vector; the candidate feature point set includes the baseline frame candidate feature point set and the floating frame candidate feature point set;

[0064] Further, step S20 includes:

[0065] Step S21: Perform neighborhood connectivity analysis on the pixels on the mesoscale skeleton in the multi-scale vascular skeleton hierarchical expression, identify bifurcation points and intersection points, collect all bifurcation points and intersection points in the mesoscale skeleton as topological key points, gather the topological key points of the reference frame into a reference frame candidate feature point set, and gather the topological key points of the floating frame into a floating frame candidate feature point set.

[0066] Specifically, step S20 takes the multi-scale vascular skeleton hierarchical expression generated in step S10 as input, extracts topologically significant feature points from it, and performs geometric encoding on each feature point to form a feature descriptor that can be used for cross-frame matching. Topological key points refer to the location points in the vascular network with special connection structures, including two categories: vascular bifurcation points and vascular intersection points. A vascular bifurcation point is the location where one vascular branch into two or more vascular branches, and a vascular intersection point is the location where two independent vascular branches intersect each other in space. The reason for selecting the mesoscale skeleton instead of the coarse-scale skeleton or the fine-scale skeleton as the feature extraction object is as follows: the coarse-scale skeleton only contains the centerline of the main vascular branch, and the number of topological key points is too small, making it difficult to provide sufficient constraints for the calculation of the spatial transformation matrix; the fine-scale skeleton contains the centerlines of all detectable vascular branches, and the number of topological key points is too large, with some positions affected by noise, resulting in pseudo-bifurcation points. Introducing a large number of unstable feature points will increase the computational burden of subsequent matching and the risk of mismatch; the mesoscale skeleton contains the centerlines of the main vascular branch and branch vascular branches, and the number of topological key points is moderate and the positions are stable, achieving a balance between feature point richness and detection stability.

[0067] See Figure 4 The implementation process of neighborhood connectivity analysis is as follows: Traverse each skeleton pixel on the mesoscale skeleton, and check the number of skeleton pixels (connection points) connected to it within its eight-neighbor window for each skeleton pixel. The eight-neighbor area refers to the eight adjacent pixel positions within a three-by-three pixel window centered on the current pixel, excluding the center point. This includes four orthogonal neighborhoods (top, bottom, left, right) and four diagonal neighborhoods (top left, top right, bottom left, bottom right). The criterion for determining connection points is: if a neighboring pixel belongs to the mesoscale skeleton and is spatially directly adjacent to the current pixel, it is determined to be a connection point. Let N be the number of connection points for the current skeleton pixel. connect When N connect When N equals 1, the current pixel is located at the end of the blood vessel, which is the endpoint; when N equals 1, the current pixel is located at the end of the blood vessel, which is the endpoint. connect When N equals 2, the current pixel is located in the middle of the blood vessel and is a normal skeleton point; when N equals 2, the current pixel is located in the middle of the blood vessel and is a normal skeleton point. connect When N equals 3, the current pixel connects to three blood vessel branches and is determined to be a bifurcation point; when N equals 3, the current pixel connects to three blood vessel branches and is determined to be a bifurcation point. connectWhen the value is 4 or greater, the current pixel connects to four or more blood vessel branches and is identified as an intersection point. The reason for using eight-neighborhood connectivity analysis instead of four-neighborhood analysis is that the orientation of the vascular skeleton in the image is arbitrary, and diagonal blood vessel connections are also common. Examining only four-neighborhood connections would miss diagonally connected branches, leading to incomplete bifurcation point identification or misjudgment of the number of connection points. The reason for collecting bifurcation points and intersection points as topological keypoints is that bifurcation points and intersection points have unique topological positions in the vascular network, and their local connectivity structure does not change due to image translation, rotation, or slight scaling, naturally possessing geometric stability as registration feature points. All topological keypoints extracted from the mesoscale skeleton corresponding to the enhanced reference frame constitute the reference frame candidate feature point set, and all topological keypoints extracted from the mesoscale skeleton corresponding to the enhanced floating frame constitute the floating frame candidate feature point set.

[0068] Step S22: Track the geometric properties of local vascular branches of each topological key point in the candidate feature point set: With each topological key point in the candidate feature point set as the center, track the branch direction and branch length along the mesoscale skeleton, and map back to the binary map of the whole vascular network to measure the branch width, so as to obtain the geometric properties including the branch direction angle, branch length and branch width.

[0069] Step S22 tracks the geometric properties of local vascular branches, including three dimensions: branch direction angle, branch length, and branch width. Local vascular branches refer to several vascular skeleton segments extending from a topological keypoint. Bifurcations naturally connect three branches, and intersections naturally connect four branches. The tracking operation starts from the topological keypoint and proceeds pixel-by-pixel along the mesoscale skeleton in each connection direction, recording the pixel coordinate sequence along the tracking path. The measurement process for the branch direction angle is as follows: Centered on the topological keypoint p, track along the i*th branch to a distance p equal to a preset step size L. step Let the position be q. i* L step The step size is set based on the curvature variation characteristics of blood vessel branches near the bifurcation point. Too short a step size will cause the orientation angle measurement to be excessively affected by single-pixel noise, while too long a step size will cause the orientation angle to deviate from the true direction of the branch at the bifurcation point. An example is L... step The value ranges from 8 pixels to 15 pixels; calculate the distance from point p to point q. i* The angle between the line connecting the i-th branches and the positive direction of the horizontal axis of the image is denoted as the branch direction angle θ of the i-th branch. i* The angle range is from 0 degrees to 360 degrees. The reason for measuring the branch direction angle relative to the horizontal axis of the image rather than relative to any arbitrary reference line is that the horizontal axis of the image has a consistent definition in all frames of the image, making the measurement results of the branch direction angle comparable and providing a unified angular benchmark for subsequent calculation of the angle between adjacent branches.

[0070] The process of measuring branch length is as follows: Starting from the topological key point p, continue tracing along the i*th branch until the next topological key point or blood vessel endpoint is encountered, and record the endpoint as p. i*,end ; Statistics from p to p i*,end The tracking path passes through all skeleton pixels. The sum of the distances between adjacent pixels is calculated as the branch arc length. The distance between adjacent pixels in the orthogonal direction is one pixel unit, and the distance between adjacent pixels in the diagonal direction is... Each pixel unit; the calculated arc length is denoted as L, the branch length of the i-th branch. i* The reason for using arc length measurement instead of straight-line distance between the start and end points is that blood vessels do not follow a strictly straight path on the retina. The arc length of curved blood vessels can more accurately reflect the actual anatomical length of the blood vessels. Straight-line distance will underestimate the length of curved blood vessels, causing blood vessels with different degrees of curvature to be incorrectly judged as similar in length characteristics.

[0071] The branch width measurement process is as follows: Since the mesoscale skeleton retains only the centerline with a single pixel width after skeletonization, the blood vessel width information is lost in the skeleton and needs to be mapped back to the full blood vessel network binary image generated in step S13 for measurement. The coordinates of the topological key point p are located in the full blood vessel network binary image. Several sampling positions are selected along the i*th branch. An exemplary number of sampling positions is 3 to 5, and the sampling positions are evenly distributed from the branch start point to the preset step size L. step Within the range; at each sampling location, measure the width of the blood vessel region in the binary image of the entire blood vessel network along a direction perpendicular to the branch direction. The measurement method is to search from the sampling location to both sides in a vertical direction until the background pixel is encountered. The sum of the search distances on both sides is the width of the blood vessel cross section at that sampling location. Take the average value of the blood vessel cross section widths at all sampling locations, and record the average value as the branch width w of the i-th branch. i* The rationale for using multi-location sampling and averaging instead of single-point measurement is that the width of blood vessels varies slightly along the vessel's direction. Single-point measurement results are easily affected by local variations in vessel thickness. Multi-location sampling and averaging can obtain a more stable width estimate and reduce the interference of random measurement errors on subsequent feature matching.

[0072] Step S23: Construct a local vascular tree signature vector based on geometric attributes: Calculate the angle between adjacent branches based on the branch direction angle in the geometric attributes, and use the angle between adjacent branches as a rotation-invariant angular feature; Normalize the branch length and branch width to obtain normalized length and normalized width, and use the normalized length and normalized width as scale-invariant geometric features; Combine the rotation-invariant angular features with the scale-invariant geometric features to generate a local vascular tree signature vector for each topological key point.

[0073] Step S23 constructs the local vascular tree signature vector. The local vascular tree signature vector is a compact mathematical description of the local topological structure of the topological keypoint and has invariance to rigid body transformations of the image. Rigid body transformation refers to spatial transformations that maintain the shape and size of an object, including translation, rotation, and mirroring. The displacement and rotation between adjacent frames in the original image sequence caused by nystagmus belong to the category of rigid body transformations. Referring to Table 1, the construction process of the rotation invariant angle feature is as follows: Sort each branch of the topological keypoint p in ascending order of branch direction angle. Taking the bifurcation point as an example, let the sorted branch direction angle sequence be θ1, θ2, θ3, where θ1 is the first branch direction angle, θ2 is the second branch direction angle, and θ3 is the third branch direction angle. Calculate the angle between adjacent branches: the angle between the first and second branches δθ12 = θ2 - θ1, the angle between the third and second branches δθ23 = θ3 - θ2, and the angle between the third and first branches δθ31 = 360° - θ3 + θ1. The sum of the included angles of adjacent branches is always equal to 360°, satisfying the closure constraint condition. The basis for using the included angles of adjacent branches, rather than the branch direction angles themselves, as the rotation-invariant angular feature is that when the image rotates, all branch direction angles will increase or decrease by the same rotation angle synchronously, and the individual branch direction angles do not have stability as they change with the rotation angle; the included angles of adjacent branches are obtained by subtracting the two branch direction angles, and the rotation angle is eliminated in the subtraction process, so the included angles of adjacent branches do not change with the image rotation, and have rotation invariance.

[0074] The process of constructing scale-invariant geometric features is as follows: The branch lengths are normalized by dividing the length of each branch by the maximum value L of all branch lengths within the local region. max The normalized length L' is obtained. i* L' i* =L i* / L max L max The value is the maximum of the branch lengths of all connected branches of the current topological key point p; the branch width is normalized by dividing the branch width of each branch by the average width w of all branches in the local region. avg The normalized width w' is obtained. i* w' i* =w i* / w avg w avgThe value is the arithmetic mean of the branch widths of all connected branches at the current topological keypoint p. Normalization is used because when the image is scaled, all branch lengths and widths are magnified or reduced proportionally, and the normalized relative ratio is unaffected by the scaling ratio, exhibiting scale invariance. The selection of the maximum length and average width of the local region as the normalization benchmark is based on the following: the maximum length corresponds to the main vessels in the local vascular network at the current bifurcation point, and its length exhibits high consistency across different frames; the average width corresponds to the overall thickness of the local vascular network at the current bifurcation point, which can offset the impact of individual branch width measurement errors on the normalization result.

[0075] The local vascular tree signature vector is combined as follows: rotation-invariant angular features and scale-invariant geometric features are arranged in a fixed order to form the local vascular tree signature vector Sp of feature point p. Taking a bifurcation point as an example, Sp is expressed as a nine-dimensional vector composed of δθ12, δθ23, δθ31, L'1, L'2, L'3, w'1, w'2, and w'3. The first three components are the angles between adjacent branches, the middle three components are the normalized lengths, and the last three components are the normalized widths. If the topological keypoint is a crosspoint connecting four branches, the local vascular tree signature vector is expanded to a twelve-dimensional vector, containing the angles between four adjacent branches, four normalized lengths, and four normalized widths. For each topological keypoint in the baseline frame candidate feature point set and the floating frame candidate feature point set, a corresponding local vascular tree signature vector is generated in the above manner.

[0076] Table 1. Physical meaning and calculation source of each component of the local vascular tree signature vector.

[0077]

[0078] Step S24: Select topological key points from the candidate feature point set of the reference frame and pair them with topological key points from the candidate feature point set of the floating frame to form multiple point pair combinations; calculate the invariant distance between the local vascular tree signature vector of the topological key point in the candidate feature point set of the reference frame and the local vascular tree signature vector of the topological key point in the candidate feature point set of the floating frame in each point pair combination, and traverse all point pair combinations to generate a feature matching cost matrix.

[0079] Step S24 calculates the feature matching cost matrix, which is a two-dimensional matrix recording the similarity measures of all possible matching pairs between the baseline frame candidate feature point set and the floating frame candidate feature point set. Let the baseline frame candidate feature point set contain N... ref N topological key points, the floating frame candidate feature point set contains N float If there are N topological key points, then the dimension of the feature matching cost matrix M is N. ref row × N floatLet the matrix element M[i,j] represent the invariant distance between the i-th topological keypoint in the candidate feature point set of the baseline frame and the j-th topological keypoint in the candidate feature point set of the floating frame. Let S be the local vascular tree signature vector of the topological keypoints in the candidate feature point set of the baseline frame. ref The local vascular tree signature vector of the topological key points in the floating frame candidate feature point set is S. float Angle weighting coefficient W θ The length weighting coefficient is W. L The width weighting coefficient is W. w Based on W θ W L and W w Calculate S ref With S float The weighted Euclidean distance is used as an invariant distance. The specific calculation method for the weighted Euclidean distance is: angle weight coefficient W... θ × Sum of squared angle differences + Length weighting coefficient W L × Sum of squared length differences + Width weighting coefficient W w × Sum of the squares of the width differences, then calculate the square root.

[0080] The sum of squared angle differences is calculated as follows: Topological keypoints in the candidate feature point set of the reference frame are defined as reference frame topological keypoints; topological keypoints in the candidate feature point set of the floating frame are defined as floating frame topological keypoints; and the local vascular tree signature vector S of the reference frame topological keypoints is calculated as follows: ref The angular components and the local vascular tree signature vector S of the topological keypoints in the floating frame float The corresponding angular components are subtracted one by one, and the squares of the differences are summed. The calculation method for the sum of squares of length and width differences is similar to that for the sum of squares of angular differences; the squares of the differences between the normalized length and normalized width components are summed separately. Weighting coefficient W θ W L W w To balance the contributions of angular, length, and width features in distance calculation, the weighting coefficients are determined based on the measurement stability of each feature in laser speckle fundus images: angular features have high measurement stability and are assigned a larger weight; length features have moderate measurement stability and are assigned a medium weight; width features have slightly lower stability due to boundary extraction errors in the binary image of the whole vascular network and are assigned a smaller weight. An exemplary weight configuration is W. θ Values: 0.5, W L Values: 0.3, W w The value is 0.2, and the sum of the three weights equals 1.

[0081] The process of generating the feature matching cost matrix is ​​as follows: traverse each topological keypoint p in the candidate feature point set of the baseline frame.ref,i For each p ref,i Traverse each topological keypoint p in the candidate feature point set of the floating frame float,j Calculate p ref,i Local vascular tree signature vector S ref,i With p float,j Local vascular tree signature vector S float,j The invariant distance D(S) between ref,i ,S float,j The calculation results are stored in matrix element M[i,j]. After traversal, each row of the feature matching cost matrix M corresponds to a topological key point in the candidate feature point set of the baseline frame, and each column corresponds to a topological key point in the candidate feature point set of the floating frame. The smaller the value of the matrix element, the more similar the corresponding point pair is in the local topological structure, and the higher the matching probability. If the topological key point in the candidate feature point set of the baseline frame is a bifurcation point and the topological key point in the candidate feature point set of the floating frame is an intersection point, the dimensions of their local vascular tree signature vectors are different. The signature vector of the bifurcation point is 9-dimensional (3 branches), and the signature vector of the intersection point is 12-dimensional (4 branches). The dimension mismatch indicates that the two are of different types. At this time, the invariant distance is set to a maximum value, indicating that the point pair has no matching possibility, thus avoiding the incorrect matching of the bifurcation point and the intersection point.

[0082] In step S20, neighborhood connectivity analysis precisely locates bifurcation and intersection points in the mesoscale skeleton, ensuring that all topological key points in the candidate feature point set have clear topological meaning and eliminating interference from ordinary skeleton points and endpoints in the matching process. Geometric attribute tracking measures branch direction angles, branch lengths, and branch widths along the mesoscale skeleton, providing each topological key point with multidimensional geometric information describing its local vascular tree structure, avoiding insufficient feature discrimination due to relying solely on position coordinates. Rotation-invariant angular features convert branch direction angles into angles between adjacent branches, eliminating the influence of image rotation on angular features and preventing rotational displacement caused by fremitus between the reference frame and the floating frame from interfering with feature matching. Scale-invariant geometric features locally normalize branch lengths and widths, eliminating the influence of image scaling on geometric features and preventing differences in vessel size caused by changes in imaging distance from interfering with feature matching. The local vascular tree signature vector combines rotation-invariant angular features and scale-invariant geometric features into a compact feature descriptor, ensuring that the local topological structure of each topological key point is fully encoded in a fixed-length vector, facilitating subsequent vector distance calculation and comparison. The feature matching cost matrix records the invariant distances of all possible matching point pairs, enabling step S30 to quickly filter candidate matching point pairs with similar local topological structures based on the distance values, thus reducing the computational scale of global topological consistency verification. The local vascular tree signature vector constructed based on relative relationships makes the feature descriptor inherently robust to rigid body transformations caused by nystagmus, allowing direct feature matching with the reference frame without pre-alignment of floating frames, simplifying the registration process. The combination of multi-dimensional geometric attributes gives the local vascular tree signature vector sufficient discriminative power. Even in densely vascular regions with multiple adjacent bifurcation points, the signature vectors of each bifurcation point exhibit different feature patterns due to differences in the combination of branch angle, branch length, and branch width, effectively alleviating the difficulty of identifying ambiguities at multiple vascular bifurcation points. The weighted Euclidean distance formula balances the contribution of various features through weight coefficients, allowing the more stable angle features to play a dominant role in the matching decision, reducing the negative impact of the less stable width features on the matching results, and further improving the accuracy of feature matching. Step S20 transforms the problem of inaccurate matching of a single feature point in laser speckle fundus image registration under the interference of nystagmus, rotation and scaling into a feature descriptor matching problem with rigid body transformation invariance. By utilizing the relative invariance of the local vascular tree structure, the accuracy of feature point matching is improved, providing high-quality candidate matching point pairs for graph matching consistency verification and spatial geometric transformation matrix calculation in step S30.

[0083] Step S30: Construct a baseline topology map and a floating topology map using the candidate feature point set; select a set of sure-matching point pairs from the baseline topology map and the floating topology map according to the feature matching cost matrix; calculate the spatial geometric transformation matrix according to the set of sure-matching point pairs; register the floating frame according to the spatial geometric transformation matrix; and output a high signal-to-noise ratio fundus image.

[0084] Further, step S30 includes:

[0085] Step S31: Take each topological key point in the candidate feature point set of the reference frame and the candidate feature point set of the floating frame as a vertex of the topological graph, establish edges between the topological key points directly connected on the mesoscale skeleton, and construct the reference topological graph and the floating topological graph respectively.

[0086] Specifically, step S30 takes the candidate feature point set and feature matching cost matrix output in step S20 as input, encodes the physical connection structure of the vascular network into a topological graph using graph theory methods, and uses the neighbor consistency relationship in the topological graph to globally verify the candidate matching point pairs, eliminating false matches caused by local feature similarity. Finally, based on the clean set of sure-matching point pairs, the spatial geometric transformation matrix is ​​calculated and image registration and fusion are completed. The topological graph is a basic data structure in graph theory, consisting of a vertex set and an edge set. Vertices represent node elements in the graph, and edges represent the connection relationships between nodes. In the laser speckle fundus image registration scenario, the topological graph is used to represent the physical connection structure between topological key points in the vascular network, with topological key points as vertices and vascular skeleton segments between adjacent topological key points as edges.

[0087] The construction process of the baseline topology graph is as follows: A candidate feature point set of the baseline frame is extracted, containing all bifurcation and intersection points extracted from the mesoscale skeleton corresponding to the enhanced baseline frame. Each topological keypoint in the candidate feature point set of the baseline frame is used as a vertex in the baseline topology graph. The vertex attributes include the position coordinates of the topological keypoint in the image coordinate system and the local vascular tree signature vector. Each pair of topological keypoints on the mesoscale skeleton is traversed to determine whether there is a directly connected vascular skeleton segment between the two topological keypoints. The criterion for direct connection is that tracing along the mesoscale skeleton from one topological keypoint can reach another topological keypoint without passing through other topological keypoints on the tracing path. If two topological keypoints satisfy the direct connection criterion, an edge is established for the corresponding two vertices in the baseline topology graph. The edge attributes include the arc length of the vascular skeleton segment connecting the two topological keypoints. The construction process of the floating topology graph is the same as that of the baseline topology graph: a candidate feature point set of the floating frame is extracted, and the floating topology graph is generated according to the above vertex and edge construction rules. The rationale for using a topological graph to represent the vascular network structure instead of just a list of feature point coordinates is that the topological graph not only preserves the geometric location information of the topological key points, but also records the physical topological structure of the vascular network through the connection relationships of the edges. Whether two topological key points are adjacent and the attributes of adjacent edges provide structural constraints for subsequent graph matching consistency verification, which is an information dimension that a coordinate list cannot provide.

[0088] The rationale for using topological keypoints as vertices, rather than selecting arbitrary pixels on the skeleton, is as follows: Topological keypoints have unique topological positions in the vascular network, their number of neighbors is fixed, and their neighbor relationships are stable. The topological graph constructed with topological keypoints as vertices has a clear graph structure, facilitating subsequent neighbor consistency verification. If arbitrary pixels on the skeleton are used as vertices, the graph will contain a large number of ordinary skeleton points with only two neighbors. These points have no topological significance and their neighbor relationships are singular, failing to provide effective structural constraints for mismatch removal. The rationale for establishing edges between directly connected topological keypoints, rather than connecting all topological keypoints pairwise, is as follows: Directly connected edges correspond to real vascular skeleton segments, and the attributes of the edges are directly related to the anatomical length of the blood vessels, possessing physical meaning. If all topological keypoints are connected pairwise, the edges will contain a large number of virtual connections spanning multiple blood vessel segments, and the attributes of the edges will lose physical meaning, failing to reflect the true topological structure of the vascular network.

[0089] Step S32: Select a set of candidate matching point pairs based on the feature matching cost matrix, perform second-order proximity constraint verification on each candidate matching point pair in the set, and remove candidate matching point pairs that fail the constraint verification to obtain a set of sure matching point pairs;

[0090] Step S32 performs graph matching consistency verification, which is a process of verifying the structural constraints of candidate matching point pairs using the neighbor relationships in the topological graph. The selection process for the candidate matching point pair set is as follows: From the feature matching cost matrix M generated in step S24, for each topological key point p in the candidate feature point set of the baseline frame... ref,i Find the topological key point p in the floating frame candidate feature point set that minimizes the matrix element M[i,j]. float,j , will point to p ref,i With p float,j As a candidate matching point pair; set a distance threshold D match D match The determination method is based on the typical variation range of the local vascular tree signature vector between different frames at the same topological key point, exemplified by D. match The value is taken as 15%-25% of the average value of each component of the local vascular tree signature vector; if the value of M[i,j] is greater than the distance threshold D match Then determine the candidate matching point pair p ref,i With p float,j If the local topological structure difference is too large, the point pair is excluded from the candidate matching point pair set; for each topological key point p in the floating frame candidate feature point set... float,j Perform the same operation to find the topological keypoint p in the candidate feature point set of the reference frame that minimizes the matrix element M[i,j]. ref,i If M[i,j] is less than or equal to the distance threshold Dmatch and the point pair p ref,i With p float,j If both the bidirectional nearest neighbor condition and the condition are met, then the pair of points is added to the candidate matching pair set. The bidirectional nearest neighbor condition refers to p ref,i It is p float,j The nearest neighbor in the candidate feature point set of the reference frame and p float,j It is p ref,i The nearest neighbor in the candidate feature point set of the floating frame. The reason for using the bidirectional nearest neighbor condition to screen candidate matching point pairs instead of only the unidirectional nearest neighbor is that: the unidirectional nearest neighbor may produce many-to-one matching, that is, multiple topological key points in the candidate feature point set of the floating frame will take the same topological key point in the reference frame as their nearest neighbor, which violates the principle of uniqueness of registration matching; the bidirectional nearest neighbor condition requires that the matching relationship is valid in both directions, effectively eliminating many-to-one matching and improving the initial quality of the candidate matching point pair set.

[0091] The implementation process of the second-order proximity constraint verification is as follows: Suppose there exists a candidate matching point pair in the candidate matching point pair set, and vertex A in the base topology graph and vertex A' in the floating topology graph constitute this matching pair; Find all neighboring vertices of vertex A in the base topology graph. A neighboring vertex is defined as a vertex that is directly connected to A through an edge. Let the set of neighboring vertices of A be NB(A)={B,C,D}. The number of neighbors is equal to the number of branches of the topological key point represented by A. The number of neighbors of a branching point is three, and the number of neighbors of a cross point is four; Find all neighboring vertices of vertex A' in the floating topology graph. Let the set of neighboring vertices of A' be NB(A')={B',C',D'}; If A and A' are a correct match, then each neighboring vertex of A should have a corresponding matching object in the candidate matching point pair set, and this matching object should be located in the set of neighboring vertices of A'. For example, if B and B' form a candidate matching pair in the candidate matching point pair set, and B belongs to NB(A) and B' belongs to NB(A'), then the matching relationship between B and B' is topologically consistent with the matching relationship between A and A'; conversely, if the matching object of B in the candidate matching point pair set is B'', and B'' does not belong to NB(A'), then the matching relationship between B and B'' contradicts the matching relationship between A and A' in terms of topological structure, and there is at least one incorrect match between A and A' or B and B''.

[0092] The verification criterion for neighbor matching consistency is as follows: For candidate matching point pairs A and A', count the number of vertices in the neighbor vertex set NB(A) of A that satisfy the neighbor matching condition. The neighbor matching condition is that the neighbor vertex has a matching object in the candidate matching point pair set and the matching object belongs to the neighbor vertex set NB(A') of A'. Let the number of vertices satisfying the neighbor matching condition be N. consistent Let the total number of neighboring vertices of A be N. total Calculate the neighbor consistency ratio R consistent =N consistent / N total Set a consistency threshold T consistent T consistent The method for determining T is based on the typical number of neighbors at a vascular network bifurcation point and the allowable tolerance for neighbor matching failures. An example is T... consistent The value range is from 0.5 to 0.8; refer to Table 2, if R consistent Greater than or equal to T consistent If R consistent Less than T consistent If the neighbor structure of candidate matching point pair A and A' cannot correspond, it is considered a structural mismatch, and it is removed from the candidate matching point pair set.

[0093] Table 2. Decision logic for verifying second-order proximity constraints:

[0094]

[0095] The iterative execution process of the second-order proximity constraint verification is as follows: Perform the above neighbor consistency verification on each candidate matching point pair in the candidate matching point pair set, and mark the candidate matching point pairs that fail the verification as to be removed; after completing one round of verification, remove all candidate matching point pairs to be removed from the candidate matching point pair set; since the removal of candidate matching point pairs may change the neighbor matching status of other candidate matching point pairs, the neighbor consistency verification needs to be performed again on the updated candidate matching point pair set; iteratively execute the verification and removal operations until there are no new point pairs to be removed in the candidate matching point pair set in two consecutive rounds of verification, at which point the candidate matching point pair set reaches a stable state; output the candidate matching point pair set in the stable state as the sure matching point pair set. The rationale for using iterative verification instead of single-round verification is as follows: In single-round verification, the existence of an incorrect matching point pair may cause its neighbor's correct matching point pair to show false positive consistency due to the incorrect match. After the incorrect matching point pair is removed, the consistency verification result of the correct matching point pair may change. Iterative verification gradually removes incorrect matches, so that the consistency verification of the remaining matching point pairs gradually approaches the true state, avoiding the interference of incorrect matches on the verification result of correct matches.

[0096] The rationale for employing second-order proximity constraint verification, rather than relying solely on the distance between local vascular tree signature vectors for matching and filtering, is as follows: Local vascular tree signature vectors only describe the local geometry surrounding a single topological keypoint. In densely populated vascular regions, multiple bifurcation points with similar local structures may exist, and the distances between their signature vectors are close. Relying solely on distance filtering makes it difficult to distinguish the correct match for these similar bifurcation points. Second-order proximity constraints introduce the structural relationship between topological keypoints and their neighbors. Even if two bifurcation points have similar local structures, the configurations of their neighboring bifurcation points are usually different. Verifying the consistency of neighbor matching can effectively distinguish bifurcation points that are locally similar but have different global positions. The core of second-order proximity constraint verification lies in transforming the physical connectivity of the vascular network into a structured constraint for matching verification: This method fully utilizes the essential property of the vascular network as a connected graph. Each bifurcation point does not exist in isolation but forms a definite connection relationship with other bifurcation points through blood vessels. This connection relationship maintains topological invariance across different frame images. By verifying the consistency of neighbor relationships between matching point pairs, the second-order proximity constraint extends the single-point matching problem into a local subgraph matching problem, with the constraint strength increasing exponentially: if a bifurcation point has 3 neighbors, then the matching consistency of all 4 points must be satisfied simultaneously, reducing the probability of incorrect matching. The implementation of the second-order proximity constraint verification improves the registration success rate in densely bifurcated regions: in the optic disc periphery region with densely bifurcated vessels, traditional methods frequently fail because they cannot distinguish similar bifurcation points, while the second-order proximity constraint accurately identifies the unique identity of each bifurcation point through the differences in neighbor relationships; even when laser speckle noise causes local feature instability, as long as the physical connectivity of the vessels remains intact, the second-order proximity constraint can still reliably filter out correct matches, making the registration algorithm more robust to noise and local feature degradation.

[0097] Step S33: Calculate the spatial geometric transformation matrix using the set of sure matching point pairs, and apply the spatial geometric transformation matrix to map the enhanced floating frame to the enhanced reference frame to obtain the registered floating frame; superimpose and fuse the registered floating frame and the enhanced reference frame, and after iteratively processing all enhanced floating frames, take the average value to output a high signal-to-noise ratio fundus image.

[0098] Step S33 calculates the spatial geometric transformation matrix using the set of confirmed matching point pairs and performs image registration fusion. The spatial geometric transformation matrix describes the mapping relationship from the coordinate system of the floating frame image to the coordinate system of the reference frame image, transforming the coordinates of each pixel in the floating frame into its corresponding position in the coordinate system of the reference frame. The calculation of the spatial geometric transformation matrix adopts a two-level transformation model: affine transformation as the global coarse registration model, and thin-plate spline interpolation as the local fine deformation model.

[0099] The calculation process of affine transformation is as follows: Affine transformation is a geometric transformation that preserves straight lines and parallelism. It can characterize the combined effects of translation, rotation, scaling, and shearing. The transformation matrix is ​​a 2x3 matrix containing six parameters to be solved. Assume the set of sure-matching point pairs contains K pairs of matching points. Each pair consists of the coordinates of the topological keypoints in the reference frame and the topological keypoints in the floating frame. The affine transformation parameters are solved using the least squares method, constructing an overdetermined system of equations. Each pair of matching points in the system provides two constraint equations, corresponding to the mapping relationship between the x and y coordinates, respectively. The six parameters of the affine transformation matrix are obtained by solving the least squares solution of the overdetermined system of equations. The least squares solution minimizes the sum of squares of the coordinate mapping residuals of all matching point pairs. The minimum size requirement for the set of sure-matching point pairs is three pairs of non-collinear matching point pairs. Three pairs of matching point pairs provide exactly six constraint equations, matching the six parameters of the affine transformation. When the number of matching point pairs exceeds three, the system of equations becomes overdetermined, and the least squares solution can, to some extent, offset the coordinate measurement errors of individual matching point pairs.

[0100] The specific combination of the two-stage transformation is as follows: First, an affine transformation matrix is ​​used to transform the coordinates of the topological keypoints of each sure-matching point pair in the floating frame candidate feature point set to the coordinate system of the reference frame, obtaining the coordinates of the floating frame feature points after affine transformation. The residual vector between the coordinates of the floating frame feature points after affine transformation and the corresponding coordinates of the topological keypoints of the reference frame is calculated. This residual vector reflects the local nonlinear deformation components that the affine transformation cannot fit. Using the coordinates of the floating frame feature points after affine transformation as the control points and the residual vector as the deformation at the control points, a thin-plate spline interpolation model is constructed. Thin-plate spline interpolation is a non-rigid deformation model that can fit local nonlinear deformations while maintaining overall smoothness. It is suitable for correcting local deformations caused by residual distortion of the ocular bowl structure in fundus images. The mathematical form of thin-plate spline interpolation includes an affine transformation part and a radial basis function part. In the two-stage transformation model, thin-plate spline interpolation is only responsible for fitting the residual deformation, and the parameters of its affine transformation part are set to zero to avoid repetition with the first-stage affine transformation. The radial basis function corresponds to the local deformation components, and each sure-matching point pair contributes a radial basis function center. The thin-plate spline interpolation parameters are solved by constructing a system of linear equations. The dimension of the equation system is related to the number of sure-matching point pairs, and the solution process involves matrix inversion or solving the linear equation system. The final spatial geometric transformation is a combination of affine transformation and thin-plate spline deformation compensation: for any pixel coordinate P in the enhanced floating frame... float First, affine transformation is applied to obtain the intermediate coordinates P. affine Then query the thin plate spline interpolation model in P affine Deformation compensation amount δP at the location tps The final transformation result is P final =P affine+δP tps The rationale for combining affine transformation and thin-plate spline interpolation, rather than using either transformation model alone, is as follows: Affine transformation alone can only fit global rigid body transformations and uniform scaling, and cannot compensate for local deformations; the registration accuracy is limited by the completeness of distortion correction. When using thin-plate spline interpolation alone, if it is certain that the number of matching point pairs is small or unevenly distributed, the deformation estimation in the sparse region at the center of the radial basis function is unstable. The two-stage transformation model first eliminates global displacement and rotation through affine transformation, and then compensates for local deformations through thin-plate spline interpolation, thus leveraging the advantages of each. Furthermore, the two-stage separate calculation avoids the numerical instability problem that occurs when thin-plate splines simultaneously fit large-scale rigid body transformations and small-amplitude local deformations.

[0101] The process of generating the registered floating frame is as follows: Apply a spatial geometric transformation matrix to each pixel coordinate in the enhanced floating frame to calculate its corresponding position in the coordinate system of the enhanced reference frame; Since the transformed coordinates are usually non-integer, bilinear interpolation is used to obtain the gray value of the corresponding position. The bilinear interpolation is calculated by weighting the gray values ​​of the four integer coordinate pixels around the target position. The weights are related to the distance from the target position to each integer coordinate; Assign the interpolated gray value to the corresponding pixel position in the registered floating frame. After traversing all pixel coordinates, the complete registered floating frame is obtained.

[0102] The overlay and fusion process is as follows: the registered floating frame and the enhanced reference frame are overlaid pixel-level in the same coordinate system; the overlay uses a weighted average method, and the gray value of the enhanced reference frame at position (x,y) is set to I. ref (x,y), the grayscale value of the registered floating frame at position (x,y) is I. reg (x,y), the merged grayscale value is I fused (x,y)=W ref ×I ref (x,y)+W reg ×I reg (x,y), where W ref To enhance the fusion weights of the reference frame, W reg For the fusion weights of the registered floating frames, W ref With W reg The sum equals one; higher weights are assigned to vascular regions. The determination of vascular regions is based on the union of the three scale binary images generated in step S13, that is, performing a logical OR operation on the main vascular binary image, the whole vascular network binary image, and the microvascular texture image. The position marked as a vascular pixel in any binary image is determined as a vascular region. If the position (x,y) is marked as a vascular pixel in the union of vascular regions, the fusion weight W of the floating frame after registration at that position is increased. reg The increase is set according to the need for vascular enhancement; for example, for non-vascular areas: Wref =0.6, W reg =0.4, for vascular regions: W ref =0.4, W reg =0.6. The reason for using a higher weight for the vascular region is that the goal of registration is to align the blood vessels with the reference frame. The registration accuracy of the vascular region directly affects the reliability of subsequent blood flow analysis. Increasing the weight of the vascular region can enhance the accumulation effect of vascular information during the fusion process.

[0103] The process of iteratively processing all floating frames is as follows: Each frame in the original image sequence, except for the reference frame, is treated as a floating frame and processed sequentially; steps S12 to S33 are repeated for each floating frame to generate the corresponding registered floating frame; all registered floating frames and the enhanced reference frame are accumulated in the same coordinate system, and the pixel-level average of the accumulated result is taken to obtain a high signal-to-noise ratio fundus image. The principle of improving the signal-to-noise ratio through multi-frame averaging is that the noise generated by laser speckle imaging is random between different frames. After multiple frames are superimposed, the noise components cancel each other out, and the signal components are superimposed on each other. The signal-to-noise ratio increases approximately proportionally with the square root of the number of superimposed frames. The reason for using post-registration averaging instead of direct averaging is that nystagmus causes the position of blood vessels to shift between different frames. If direct averaging is performed without registration, the edges of blood vessels will be blurred due to non-alignment, and the contrast of blood vessels will decrease. After registration, the positions of blood vessels in all frames are aligned with the reference frame, and the averaging operation will not blur the edges of blood vessels, thus maintaining or even enhancing the clarity of blood vessels.

[0104] The mesoscale skeleton in the multi-scale vascular skeleton hierarchy representation generated in step S10 provides the basis for edge connections in the topology graph construction in step S31. The single-pixel width characteristic of the mesoscale skeleton makes the determination of whether there is a direct connection between two topological key points simple and clear. The binary image of the main blood vessel generated in step S10 provides the basis for vascular region determination in the weighted fusion in step S33, allowing the fusion weights to be set differently according to the location of the blood vessels. The edge structure of the topology graph directly inherits the topological connectivity of the mesoscale skeleton, and the edge attributes correspond precisely to the anatomical blood vessel length, enabling the second-order proximity constraint verification to be based on the real vascular network structure, thus ensuring the reliability of the verification results. The signature vector constructed based on the local topology structure in step S20 already reflects the local similarity of a single topological key point. Step S32 further verifies the structural consistency between the topological key point and its neighbors. The dual constraints of local similarity and global consistency form a hierarchical mismatch elimination mechanism. Topological key points with similar local features but different global positions are effectively identified and eliminated in the second layer of constraints.

[0105] In step S30, the topology graph construction explicitly encodes the physical connection structure of the vascular network as edge relationships in the graph, allowing direct querying of neighbor information for each topological key point, thus providing a structured input for second-order proximity constraint verification. Second-order proximity constraint verification utilizes the consistency of neighbor matching to globally verify candidate matching point pairs. Point pairs with similar local features but mismatched neighbor structures are identified as structural mismatches and eliminated, making the purity of the set of sure-matching point pairs higher than that of the candidate matching point pairs selected solely based on local features. This increased purity of the set of sure-matching point pairs means that the input to the spatial geometric transformation matrix contains little or no mismatching point pairs, the least-squares estimate of the affine transformation parameters is closer to the true value, the radial basis function center position of the thin-plate spline interpolation is accurate, and the transformation model can accurately fit the spatial mapping relationship between the base frame and the floating frame. The accurate transformation model ensures that the vessel positions in the registered floating frames are highly aligned with the reference frames. When multiple frames are superimposed and averaged, the vessel edges are not blurred due to misalignment. Noise components cancel each other out while vessel information is superimposed, resulting in a high signal-to-noise ratio fundus image with high vessel contrast and low noise level.

[0106] By incorporating the physical connectivity structure of the vascular network into the matching verification process, even in densely populated vascular regions with multiple bifurcation points exhibiting similar local features, the differences in neighbor configurations at each bifurcation point can serve as a distinguishing factor. This effectively separates correctly matched and incorrectly matched point pairs, significantly improving the registration success rate in densely populated multi-vascular bifurcation regions. The iterative second-order proximity constraint verification possesses self-correcting capabilities. A small number of initially retained incorrectly matched point pairs are gradually eliminated in subsequent iterations due to changes in neighbor matching relationships, ensuring that the set of confirmed matched point pairs reaches a high purity state after iterative stabilization. The synergy of the two-level transformation models allows for effective fitting of both global rigid body transformations and local nonlinear deformations, freeing registration accuracy from the expressive power of a single transformation model. Differential fusion weights based on the main vessel binary map enable more comprehensive information accumulation in the vascular region across multiple frames, resulting in clearer and sharper vascular structures in high signal-to-noise ratio fundus images, providing high-quality input data for subsequent hemodynamic analysis. Step S30 transforms the problems of mismatch due to similar local features and difficulty in identifying multiple vascular bifurcation points in laser speckle fundus image registration into a consistency verification problem based on graph topology. By utilizing the constraint of the neighbor relationship of the vascular network, the accuracy of feature point matching and the registration success rate of dense areas of multiple vascular bifurcation points are both improved.

[0107] Example 2

[0108] This embodiment, based on Embodiment 1, provides a laser speckle fundus image registration system, such as... Figure 5 As shown, it includes:

[0109] Multi-scale skeleton generation module: used to read laser speckle fundus video and discretize it into a raw image sequence, determine the reference frame and floating frame, perform spherical projection correction on the reference frame and floating frame to generate enhanced fundus image, process the enhanced fundus image, and generate multi-scale vascular skeleton hierarchical expression.

[0110] The signature vector construction module is used to detect topological key points based on multi-scale vascular skeleton hierarchical expression and form a candidate feature point set; track the geometric properties of local vascular branches of each topological key point in the candidate feature point set, construct local vascular tree signature vectors based on geometric properties, and calculate the feature matching cost matrix based on local vascular tree signature vectors.

[0111] Image registration module: Constructs a baseline topology map and a floating topology map using candidate feature point sets. Selects a set of sure-matching point pairs from the baseline topology map and the floating topology map based on the feature matching cost matrix. Calculates the spatial geometric transformation matrix based on the set of sure-matching point pairs. Registers the floating frames based on the spatial geometric transformation matrix and outputs a high signal-to-noise ratio fundus image.

[0112] Furthermore, in the multi-scale skeleton generation module, the method for determining the reference frame and the floating frame includes:

[0113] A center-periphery weighted scoring mechanism is applied to each image in the original image sequence to calculate a comprehensive quality score. The image with the highest comprehensive quality score is selected as the reference frame, and the remaining images in the original image sequence are used as floating frames to be registered.

[0114] Furthermore, the method for generating multi-scale vascular skeleton hierarchical expressions includes:

[0115] Three sets of Gaussian-matched filters (coarse-scale, meso-scale, and fine-scale) are set up to perform convolution and thresholding on enhanced fundus images to generate binary maps of main vessels, full vascular network, and microvascular textures. Iterative erosion skeletonization is then performed on the main vessel binary map, full vascular network binary map, and microvascular texture map to obtain coarse-scale skeleton, meso-scale skeleton, and fine-scale skeleton. The coarse-scale skeleton, meso-scale skeleton, and fine-scale skeleton are combined to form a multi-scale hierarchical representation of the vascular skeleton.

[0116] Furthermore, in the signature vector construction module, the method for calculating the feature matching cost matrix based on the local vascular tree signature vector includes:

[0117] Topological keypoints in the candidate feature point set of the baseline frame are selected and paired with topological keypoints in the candidate feature point set of the floating frame to form multiple point pair combinations. The invariant distance between the local vascular tree signature vector of the topological keypoint in the candidate feature point set of the baseline frame and the local vascular tree signature vector of the topological keypoint in the candidate feature point set of the floating frame is calculated in each point pair combination. All point pair combinations are traversed to generate a feature matching cost matrix.

[0118] Furthermore, in the image registration module, the method for constructing the baseline topology map and the floating topology map includes:

[0119] Each topological keypoint in the candidate feature point set of the baseline frame and the candidate feature point set of the floating frame is used as a vertex of the topological graph. Edges are established between the topological keypoints that are directly connected on the mesoscale skeleton to construct the baseline topological graph and the floating topological graph respectively.

[0120] In the image registration module, the method for filtering the set of sure-matching point pairs includes:

[0121] Based on the feature matching cost matrix, a set of candidate matching point pairs is selected. For each candidate matching point pair in the set, a second-order proximity constraint verification is performed. Candidate matching point pairs that fail the constraint verification are eliminated to obtain a set of sure matching point pairs.

[0122] The methods and systems of this application may be implemented in many ways. For example, they may be implemented by software, hardware, firmware, or any combination of software, hardware, and firmware. The above-described order of steps for the method is for illustrative purposes only, and the steps of the method of this application are not limited to the order specifically described above, unless otherwise specifically stated.

[0123] In addition, the parts of the technical solutions provided in the embodiments of this application that are consistent with the implementation principles of the corresponding technical solutions in the prior art have not been described in detail, so as to avoid excessive elaboration.

[0124] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the invention. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A laser speckle fundus image registration method, characterized in that, The method includes: Read the laser speckle fundus video and discretize it into a raw image sequence. Determine the reference frame and floating frame. Perform spherical projection correction on the reference frame and floating frame to generate enhanced fundus images. Process the enhanced fundus images to generate multi-scale vascular skeleton hierarchical expression. The method for generating a multi-scale vascular skeleton hierarchical representation includes: setting three sets of Gaussian-matched filters at coarse, medium, and fine scales, performing convolution and threshold segmentation on the enhanced fundus image respectively to generate a main vessel binary map, a full vascular network binary map, and a microvascular texture map; performing iterative erosion skeletonization processing on the main vessel binary map, the full vascular network binary map, and the microvascular texture map respectively to obtain a coarse-scale skeleton, a medium-scale skeleton, and a fine-scale skeleton; and combining the coarse-scale skeleton, the medium-scale skeleton, and the fine-scale skeleton to form a multi-scale vascular skeleton hierarchical representation. Topological key points are detected in the mesoscale skeleton of the multi-scale vascular skeleton hierarchical representation to form a candidate feature point set; the geometric properties of the local vascular branches of each topological key point in the candidate feature point set are tracked, a local vascular tree signature vector is constructed based on the geometric properties, and the feature matching cost matrix is ​​calculated based on the local vascular tree signature vector. A baseline topology map and a floating topology map are constructed using a candidate feature point set. A set of sure-matching point pairs is selected from the baseline topology map and the floating topology map based on the feature matching cost matrix. A spatial geometric transformation matrix is ​​calculated based on the set of sure-matching point pairs. The floating frame is registered based on the spatial geometric transformation matrix. The registered floating frame and the baseline frame are differentially weighted and fused with the vascular region determined by the union of the main vessel binary map, the whole vessel network binary map and the microvessel texture map to output a high signal-to-noise ratio fundus image.

2. The laser speckle fundus image registration method according to claim 1, characterized in that, The method for determining the reference frame and the floating frame includes: A center-periphery weighted scoring mechanism is applied to each image in the original image sequence to calculate a comprehensive quality score. The image with the highest comprehensive quality score is selected as the reference frame, and the remaining images in the original image sequence are used as floating frames to be registered.

3. The laser speckle fundus image registration method according to claim 2, characterized in that, The candidate feature point set includes a reference frame candidate feature point set and a floating frame candidate feature point set; The method for forming a candidate feature point set includes: Neighborhood connectivity analysis was performed on pixels on the mesoscale skeleton in the multi-scale vascular skeleton hierarchical representation to identify bifurcation points and intersection points. All bifurcation points and intersection points in the mesoscale skeleton were collected as topological key points. The topological key points of the baseline frame were aggregated into a candidate feature point set for the baseline frame, and the topological key points of the floating frame were aggregated into a candidate feature point set for the floating frame.

4. The laser speckle fundus image registration method according to claim 3, characterized in that, The method for tracking the geometric properties of local vascular branches at each topological key point in the candidate feature point set includes: Centered on each topological keypoint in the candidate feature point set, the branch direction and branch length are traced along the mesoscale skeleton, and mapped back to the binary map of the whole vascular network to measure the branch width, thus obtaining geometric attributes including branch direction angle, branch length and branch width.

5. The laser speckle fundus image registration method according to claim 4, characterized in that, The method for constructing a local vascular tree signature vector based on geometric properties includes: The angle between adjacent branches is calculated based on the branch direction angle in the geometric attributes, and the angle between adjacent branches is used as a rotation invariant angular feature. The branch length and branch width are normalized to obtain normalized length and normalized width, and the normalized length and normalized width are used as scale invariant geometric features. The rotation invariant angular features and scale invariant geometric features are combined to generate the local vascular tree signature vector of each topological key point.

6. The laser speckle fundus image registration method according to claim 5, characterized in that, The method for calculating the feature matching cost matrix based on the local vascular tree signature vector includes: Topological keypoints in the candidate feature point set of the baseline frame are selected and paired with topological keypoints in the candidate feature point set of the floating frame to form multiple point pair combinations. The invariant distance between the local vascular tree signature vector of the topological keypoint in the candidate feature point set of the baseline frame and the local vascular tree signature vector of the topological keypoint in the candidate feature point set of the floating frame is calculated in each point pair combination. All point pair combinations are traversed to generate a feature matching cost matrix.

7. The laser speckle fundus image registration method according to claim 6, characterized in that, The method for constructing the baseline topology graph and the floating topology graph includes: Each topological keypoint in the candidate feature point set of the baseline frame and the candidate feature point set of the floating frame is used as a vertex of the topological graph. Edges are established between the topological keypoints that are directly connected on the mesoscale skeleton to construct the baseline topological graph and the floating topological graph respectively.

8. The laser speckle fundus image registration method according to claim 7, characterized in that, The method for filtering the set of sure-matching point pairs includes: Based on the feature matching cost matrix, a set of candidate matching point pairs is selected. For each candidate matching point pair in the set, a second-order proximity constraint verification is performed. Candidate matching point pairs that fail the constraint verification are eliminated to obtain a set of sure matching point pairs.

9. A laser speckle fundus image registration system, used to implement the laser speckle fundus image registration method according to any one of claims 1-8, characterized in that, The system includes: Multi-scale skeleton generation module: used to read laser speckle fundus video and discretize it into a raw image sequence, determine the reference frame and floating frame, perform spherical projection correction on the reference frame and floating frame to generate enhanced fundus image, process the enhanced fundus image, and generate multi-scale vascular skeleton hierarchical expression. The signature vector construction module is used to detect topological key points based on multi-scale vascular skeleton hierarchical expression and form a candidate feature point set; track the geometric properties of local vascular branches of each topological key point in the candidate feature point set, construct local vascular tree signature vectors based on geometric properties, and calculate the feature matching cost matrix based on local vascular tree signature vectors. Image registration module: Constructs a baseline topology map and a floating topology map using candidate feature point sets. Selects a set of sure-matching point pairs from the baseline topology map and the floating topology map based on the feature matching cost matrix. Calculates the spatial geometric transformation matrix based on the set of sure-matching point pairs. Registers the floating frames based on the spatial geometric transformation matrix and outputs a high signal-to-noise ratio fundus image.

Images