CN117391955A - Convex set projection super-resolution reconstruction method based on multi-frame optical coherence tomography images - Google Patents

Abstract

The invention discloses a convex set projection super-resolution reconstruction method based on multi-frame optical coherence tomography, which mainly comprises two parts, wherein an optical coherence tomography low-resolution image sequence is firstly collected; screening multiple frames with the same section; performing angle correction on the data set; the second section includes: establishing an image degradation model; interpolating the reference frame to a target resolution; defining a closed convex set for each pixel of the high resolution image; calculating residual errors pixel by a degradation model; introducing all reference frame corrections frame by frame; and (5) performing correction circularly for a plurality of times. And reconstructing high-frequency information in the high-resolution image by fusing information in the multi-frame low-resolution image. The operation object of the reconstruction process is image blocks corresponding to a reference frame and a reference frame, different image blocks are registered according to corresponding sub-pixel displacement, firstly, displacement estimation and image reconstruction are carried out in a blocking mode, and then block images reconstructed by super resolution are spliced according to the original position, so that a whole super resolution image is obtained.

Description

Convex set projection super-resolution reconstruction method based on multi-frame optical coherence tomography images

Technical field

The invention belongs to optical medical image processing technology, and specifically relates to a method of constructing super-resolution images based on multi-frame optical coherence tomography scanning.

Background technique

Optical coherence tomography (Optical Coherence Tomography--OCT) has been widely used as a detection method in ophthalmology clinical diagnosis of retinal diseases due to its characteristics of non-invasive, non-contact and application in internal imaging of living bodies. The retinal B-scan image obtained by OCT provides clinicians with thickness and morphological information of each cell layer, which is of great significance for the detection and evaluation of diseases with obvious changes in retinal morphology.

Clinical diagnosis usually requires high-resolution OCT images to accurately reflect the complex pathological structure of the membrane layer. Although the micron-level resolution of OCT images has played an important role in clinical imaging, if the resolution of OCT images is further improved, it can provide more microstructure information of the film layer, which will help doctors achieve accurate diagnosis. Reveal the pathological manifestations of specific diseases and greatly promote the development of clinical medical diagnosis.

OCT system imaging mainly relies on the one-dimensional interference signal formed by the interference of backscattered light from different areas of the imaging background and the reference beam; an A-scan image is formed by demodulating the one-dimensional interference signal, and multiple A-scans formed by scanning in space The images are stitched together to form a two-dimensional B-scan image. Therefore, its axial resolution depends on the central wavelength and spectral bandwidth of the interference light source, and its lateral resolution is related to the numerical aperture of the focusing objective lens and the size of the light spot projected onto the surface of the object to be measured. Therefore, researchers have improved the imaging resolution of the OCT system by improving hardware configuration, including using ultra-broadband light sources and complex scanning beam focusing optical systems. However, this method is limited by the diffraction limit and laser technology. Improving the system hardware can only improve the resolution to a limited extent. Moreover, the complex optical system and light source increase the manufacturing cost of the OCT system, and it cannot be widely used in the market and medical front lines.

Compared with improving the hardware system, the method of improving the resolution of OCT images by reconstructing OCT images has been widely praised in recent years because it can improve image quality at a lower cost and is not limited by the diffraction limit. Existing super-resolution reconstruction methods can be divided into three categories: (1) super-resolution reconstruction based on a single image: (2) super-resolution reconstruction based on multiple images; (3) super-resolution reconstruction based on deep learning . Among them, the first type requires calculating the point spread function (PSF) based on the parameters of the optical system, and then using the calculated point spread function to perform a deconvolution operation to obtain a clear image. However, the original OCT data contains a large amount of speckle noise due to autocorrelation term interference during demodulation, and deconvolution of the noisy image will produce a large number of artifacts. Moreover, mature OCT systems are highly complex, optical parameters are difficult to obtain due to patent protection, and it is difficult to calculate point spread functions. Although the third category (super-resolution reconstruction based on deep learning) methods can efficiently generate higher-definition images, they rely on correct low-resolution and high-resolution mapping data sets, and their resolution is limited. The improvement factor is limited, the authenticity of the generated high-resolution images is reduced, and it cannot be applied in the field of medical images. Therefore, the second type of scheme proposed based on multi-image super-resolution reconstruction (retina OCT image) can reduce the impact of point spread function estimation errors based on single-image schemes on the reconstruction results. Furthermore, multi-image methods are robust to low-resolution image noise. When a certain low-resolution image contains noise or less information, other images can still provide effective information to reconstruct the lost details, thus improving the robustness of the reconstruction results.

Therefore, the advantage of the multi-image super-resolution reconstruction method is that it does not depend on the data set, can be extended to more application scenarios, and can meet the needs of OCT image super-resolution reconstruction in various scenarios.

A US patent (US5978109) discloses a super-resolution scanning optical device, which proposes a device for imaging light from a coherent light source in the form of fine spots on the conjugate surface, and a device for imaging the conjugate surface in the form of fine spots. Scan the device. The coherent light source consists of a first light source and a second light source with opposite phases. The side amplitude of the main image lobe of the first light source is offset by the amplitude of the main lobe of the second light source, thereby reducing the width of the main image lobe of the first light source. This achieves super-resolution below the diffraction limit without forming slit-like or annular openings. However, the disadvantage of this technical feature is that it cannot improve the lateral and axial resolution at the same time, and it requires the use of two light sources with strictly opposite phases, which is costly and technically difficult to implement.

A US patent (US5994690) discloses an improved OCT imaging system and proposes a method for estimating impulse response from the interferometer output interference signal. The impulse response can be obtained from cross-correlation and autocorrelation data. The autocorrelation power spectrum and cross-correlation power spectrum are obtained by performing Fourier transform on the autocorrelation data and cross-correlation data respectively. Taking the ratio of the cross power spectrum and the autopower spectrum, the transfer function of the linear translation invariant system of the interaction between OCT and tissue is obtained, and the light pulse response is obtained by performing inverse Fourier transform on the transfer function. Coherent demodulation and deconvolution technology are combined to resolve closely spaced reflection points in the sample and improve the axial resolution of the OCT system. However, the disadvantage of this technical feature is that the system and method cannot simultaneously improve the axial and lateral resolution of the OCT system, and its measurement or estimation of the power spectrum is affected by optical components, measurement electronic equipment, data acquisition systems and various noise sources. characteristics, which in turn affects the estimation and calculation of the transfer function, so the clarity of the reconstructed image cannot be guaranteed.

A Chinese patent (CN105976321) discloses an optical coherence tomography image super-resolution reconstruction method and device. It is proposed to divide the three-dimensional low-resolution OCT image into multiple similar frame groups in the time dimension, and divide each OCT image in the multiple similar frame groups into multiple film layers. Each frame of OCT image is divided into multiple overlapping image blocks, and multiple similar image blocks of each image block are determined within the film layer corresponding to each image block. Based on the determined multiple similar image blocks of each image block, the average image block of each image block is obtained. Through the mapping equation of the pre-constructed high-low resolution dictionary pair and the corresponding sparse coefficients, the average image block of each image block obtained is processed to obtain a high-resolution image of the three-dimensional OCT image. By constructing multiple sets of mutually matched high-low resolution OCT images as training samples, the accuracy of reconstructed super-resolution images is improved. However, the disadvantage of this technical feature is that its reconstruction accuracy relies on a large number of samples of high-low resolution images as training sets. This method is easily affected by individual differences, causing the accuracy to decrease, and is used in the process of calculating the average image block. The similar image blocks are the internal image blocks of the film layer. Although the reconstruction noise can be suppressed, the edge information will be weakened accordingly and the high-frequency components in the reconstructed image will be reduced.

A Chinese patent (CN116128728) discloses an OCT super-resolution reconstruction method and device based on equivariant learning and prior guidance. It is proposed that the first constraint function of the network to be trained is constructed based on the constructed OCT imaging model; the second constraint function of the network to be trained is constructed based on equivariant learning; and the third constraint function of the network to be trained is constructed based on the anatomical prior knowledge of the tissue being measured. Constraint function; use the collected noisy low-resolution images as a training set, and train the neural network to be trained based on the above three constraint functions; perform high signal-to-noise ratio super-operation on the OCT image to be reconstructed based on the trained network resolution reconstruction. Equivariant learning is used to achieve self-supervised super-resolution reconstruction, and the anatomical prior knowledge of the tested tissue is used to achieve self-supervised denoising, so that the network can simultaneously achieve both super-resolution and noise reduction functions for OCT image reconstruction. However, the disadvantage of this technical feature is that the reconstruction is based on the self-supervised deep learning method and prior knowledge of a single image. The image resolution is limited in the improvement factor, and the data authenticity of the reconstructed high-resolution image is reduced, and it is susceptible to erroneous data guidance. , causing the misdiagnosis rate to increase.

A Chinese patent (CN116630154) discloses a deconvolution super-resolution reconstruction method and device for OCT images. It is proposed to obtain low-resolution OCT images through Fourier transform as the input original data, construct an optimization function for sparse continuous prior deconvolution calculation, perform initial settings for reconstruction optimization, and then conduct iterative training and introduce intermediate variables for iterative calculations. . After completing the optimization iteration, the final deconvolution super-resolution reconstructed OCT image is output. This method can avoid artifacts caused by deconvolution iterative reconstruction and effectively improve the resolution of the reconstructed image. However, the disadvantage of this technical feature is that the point spread function (PSF) used in the process of constructing the optimization function is an estimated value, and the imaging principle of the OCT system is a scanning interference system, and the obtained signal is a one-dimensional signal, so it is not theoretically possible. There are two-dimensional PSFs available for reconstruction calculations. Inaccurate estimation of the point spread function may cause artifacts in the reconstructed image. And the deconvolution reconstruction based on a single image is limited by the amount of information, and the authenticity of the reconstructed high-resolution image decreases.

In summary, improving the resolution of the OCT system by improving hardware configuration has problems such as high cost and system complexity. Compared with improving the hardware system, improving the resolution of the OCT system through image reconstruction has more advantages, but the following problems limit its wide application in the clinical field: (1) The clinically used OCT system is highly integrated and has poor parameters. The differences are large and difficult to obtain, making it difficult and complex to calculate the point spread function; (2) B-scan images essentially do not have a two-dimensional point spread function that can be used for reconstruction calculations, so the point spread function is inaccurately estimated. It will cause artifacts in the reconstruction results, thereby reducing the image quality; (3) Reconstruction methods based on single images and deep learning rely on hardware system parameters and correct high-low resolution image pairs. Some missing information will cause the reconstructed image to be unrealistic. The accuracy decreases, causing the clinical misdiagnosis rate to increase; (4) For reconstruction methods based on single images and deep learning, the reconstruction resolution multiple is limited by the amount of data and model structure, and there is an upper limit, making it impossible to recover more detailed information.

In response to the above problems, the present invention proposes a convex set projection super-resolution reconstruction method based on multi-frame OCT images, which does not require calculation of hardware system parameters and construction of corresponding high-low resolution data sets. Its reconstruction factor depends on the number of basic image frames, so it can theoretically break through the upper limit of the current super-resolution factor and obtain clear and realistic high-resolution images.

Contents of the invention

The purpose of this invention is to propose a convex set projection super-resolution reconstruction method based on multi-frame optical coherence tomography images, which reconstructs high-resolution images in high-resolution images by fusing different detailed information contained in multi-frame low-resolution images. frequency information to obtain clear, high-resolution images. The Chinese meaning of OCT (Optical Coherence Tomography) in the article is optical coherence tomography.

The solution proposed by the present invention includes the following steps:

(1) Collect OCT tomography low-resolution image sequences, change the scanning mode of the OCT system, and scan the same section in a short time to obtain multiple frames of B-scan images. Select one of the frames as the benchmark. Frames, other images are used as reference frames for calculating and reconstructing high-resolution images.

(2) Screen the multi-frame B-scan images of the same section, compare the reference frames with the selected reference frames, and select the reference frames with large PSNR (Peak Signal-to-Noise Ratio) and SSIM (Structural Similarity Index) values. Reference frames are fused for reconstruction. The selected base frames and filtered reference frames constitute the reconstructed low-resolution data set.

(3) Use the SIFT (Scale Invariant Feature Matching Algorithm) algorithm to perform angle correction on the above data set, adjust the direction of the reference frame to be consistent with the reference frame, and eliminate the angular influence in subsequent displacement estimation. All feature points and corresponding feature vectors in the base frame and reference frame are detected through the SIFT algorithm, and the feature points in the reference frame that match the feature points in the base frame are found by calculating the Euclidean distance of the feature vectors in the two images. Establish corresponding relationships. All matching feature points are screened, some feature points with the highest matching degree are selected, their coordinates are counted, and the average rotation angle of the two images is calculated. Correct the reference frame image according to the calculated rotation angle.

(4) Cut the base frame and reference frame into image blocks, and use the frequency domain registration method to calculate the displacement between the corresponding image blocks. In the previous step, the angle correction has been performed based on the SIFT algorithm, so the calculated displacement is ignored. Angle influence, it is considered that there is only displacement in the x and y directions, and the unit is the number of pixels.

Frequency domain registration principle:

The displacement of image f ₂ (x) in the spatial domain relative to the image f ₁ (u) corresponds to the phase shift of the image F ₂ (u) in the frequency domain relative to the image F ₁ (u), so the displacement parameter Δx can be calculated by calculating the phase difference ∠ (F ₂ (u)/F ₁ (u)) is obtained. The relative displacement of each image block is recorded, and convex set projection super-resolution image reconstruction is performed based on the sub-pixel displacement (Δx, Δy).

(5) The convex set projection image reconstruction process is mainly divided into the following six steps:

(5.1) Establish an image degradation model: First, establish an image acquisition model that connects the original high-resolution image and the low-resolution observation sequence, which is an iterative formula for the unknown high-resolution image, which can be expressed as:

(5.2) Reference frame interpolation to target resolution: Use linear interpolation method to interpolate the reference frame image, use the interpolated reference frame as the initial estimate of the high-resolution image, and use the prior information provided by the reference frame pixel by pixel on this basis. Perform pixel value correction.

(5.3) Define a closed convex set for each pixel of the high-resolution image: In the actual imaging process, noise is inevitably interfered, so the influence of additive noise is introduced based on the linear translation fuzzy constraint to make the noise The statistical properties are linked to the a priori boundary δ ₀ . If the additive noise is Gaussian distributed and its variance is σ _v , then δ ₀ should be equal to cσ _v (c≥0), determined by a certain statistical confidence. Define the following closed convex set for each pixel of the reconstructed high-resolution image:

in, Represents the residual difference between x(i ₁ , i ₂ ) and g(m ₁ , m ₂ ).

(5.4) Calculate the residual pixel by pixel according to the degradation model to determine whether the boundary limit is exceeded. The pixel value of the reference frame image after correcting the interpolation must limit the pixel value not to exceed the threshold: Calculate the reference frame and the simulated low-resolution downsampled by the point spread function The residual difference between rate images is corrected according to the following formula according to whether the residual value exceeds the boundary limit.

The corrected pixel value needs to be within the image threshold interval [0, 255]

(5.5) Introduce all reference frames frame by frame for correction: read in all reference frames in sequence, traverse all pixels in the reference frames, and calculate the residual. And the residual value is corrected on the reference frame that has been amplified by interpolation, and the correction position corresponds to:

(5.6) Loop multiple times for correction: Limited by the amplitude of the point spread function, one iteration cannot completely integrate the supplementary information in the reference frame into the high-resolution reference frame, so multiple iterations are needed to make different reference frames Different information of the same imaging object collected on the frame image is fully integrated.

The operation object of the above reconstruction process is the image blocks at the corresponding positions of the reference frame and the base frame. Different image blocks can be registered according to the corresponding sub-pixel displacements, and the pixel value correction position is more accurate. Therefore, displacement estimation and image reconstruction are performed in blocks, and then the super-resolution reconstructed block images are spliced according to their original positions to obtain the entire super-resolution image.

The characteristics and beneficial effects of the present invention are that the proposed convex set projection super-resolution reconstruction method based on multi-frame optical coherence tomography scanning mainly embodies the following four advantages:

(1) The horizontal and vertical resolution of the OCT system can be improved at the same time without the need to modify or upgrade the system hardware, thereby reducing the manufacturing cost and technical requirements of the high-resolution OCT system, and can break through the diffraction limit and obtain sub- Micron-level OCT pixel resolution.

(2) Super-resolution reconstruction technology based on multi-frame images is different from deconvolution technology based on single images. It does not have high dependence on the point spread function and does not require accurate calculation of the point spread function. It avoids the difficulty of obtaining optical scanning system parameters due to the high integration level of the OCT system, and reduces the complexity of calculations. High-resolution images can be reconstructed simply by changing the scanning mode of the OCT system, so this method can be applied to all parameter models of OCT equipment.

(3) The resolution multiple that can be improved depends on the number of frames in the reconstructed data set. Generally speaking, the square of the reconstruction multiple should not exceed the number of frames in the data set (k ² ≤ n). Therefore, it is different from the ultra-high resolution based on "deep learning". "Resolution reconstruction technology" is limited by model parameters, and can achieve super-resolution reconstruction at any multiple if the data set supports it.

(4) This method performs reconstruction based on the prior information provided by multiple frames of low-resolution images. The amount of information of a single imaging object is richer than "super-resolution reconstruction based on single image deconvolution" and "deep learning", so the reconstruction The image realism is higher than the other two methods. It makes up for the shortcomings of super-resolution technology based on single image deconvolution because the OCT system itself does not have a two-dimensional space point spread function, so the point spread function is inaccurately estimated. It avoids the need for a large number of correct high-low resolution images for training sets based on deep learning-based super-resolution technology. It can achieve high-efficiency reconstruction of real high-resolution OCT images.

The method proposed by the present invention can realize high-efficiency reconstruction of real high-resolution OCT images and reduce the manufacturing cost of high-resolution OCT equipment. It does not require the support of a large number of data sets, nor does it require in-depth knowledge of the parameters of the OCT optical scanning system, thus reducing the computational cost and complexity, and is applicable to all parameters and models of OCT systems. More importantly, this method can theoretically break through the diffraction limit and achieve submicron-level OCT scanning accuracy.

Description of the drawings

Figure 1 is a flow chart of calculation and reconstruction of the present invention.

Figure 2(a) shows the peak signal-to-noise ratio of the reference frame and the reference frame.

Figure 2(b) shows the structural similarity between the reference frame and the base frame.

Figure 3 shows the image displacement of the reference frame relative to the base frame.

Figure 4 is an image of the SIFT algorithm feature point detection and matching.

Figure 5 shows the reconstructed high-resolution image and the comparison interpolated high-resolution image.

Detailed ways

In order to further clarify the method of the present invention, the technical solution of the present invention is described in detail below with reference to the accompanying drawings and through specific embodiments.

The convex set projection super-resolution reconstruction method based on multi-frame optical coherence tomography images, the technical solution is to reconstruct the high-frequency information in the high-resolution image by fusing the information contained in the multi-frame low-resolution image to obtain clear high resolution image. The reconstruction method includes two parts: motion estimation and convex set projection image reconstruction. Motion estimation includes 4 steps, and convex set projection image reconstruction includes 6 steps.

Motion estimation:

(1) Collect OCT tomography low-resolution image sequences to change the scanning mode of the optical coherence tomography system so that it can quickly scan the same section to obtain multiple frames of B-scan images. Select one of the images as the reference frame and the other images. As a reference frame, used to calculate the reconstructed high-resolution image.

(2) Screen the multi-frame B-scan images of the same section, compare the reference frames with the selected benchmark frames, and select the reference frames with large PSNR and SSIM values to be fused with the benchmark frames for reconstruction. The selected benchmarks frames and filtered reference frames constitute the reconstructed low-resolution data set.

(3) Use the SIFT algorithm to perform angle correction on the data set obtained in step 2, adjust the direction of the reference frame to be consistent with the reference frame, and eliminate the angular influence in subsequent displacement estimation. Use the SIFT algorithm to detect all feature points and corresponding feature vectors in the base frame and the reference frame, and calculate the Euclidean distance of the feature vectors in the two images to find the feature points in the reference frame that match the feature points in the base frame, and establish Correspondence. Filter all matching feature points, select the feature point with the highest matching degree, and calculate its coordinates. Calculate the average rotation angle of the two images, and correct the reference frame image according to the calculated rotation angle.

(4) Cut the base frame and reference frame into image blocks, and use the frequency domain registration method proposed by Vandewalle et al. to calculate the displacement between the corresponding image blocks. In step 3, the angle correction has been performed based on the SIFT algorithm, so the calculated The angular influence that exists in the displacement is considered to only exist in the x and y directions, and is measured in the number of pixels.

Vandewalle frequency domain registration principle calculation formula:

The displacement of image f ₂ (x) relative to image f ₁ (x) in the spatial domain corresponds to the phase shift of image F ₂ (u) relative to image F ₁ (u) in the frequency domain. The displacement parameter Δx can be obtained by calculating the phase difference ∠(F ₂ (u)/F ₁ (u)). The relative displacement of each image block is recorded, and convex set projection super-resolution image reconstruction is performed based on the sub-pixel displacement (Δx, Δy).

After the above is completed, the high-resolution image is reconstructed using the convex set projection method. The convex set projection method assumes that the unknown image is an element in a suitable Hilbert space, and each low-resolution reference frame is used as a precursor to the unknown image. Based on empirical knowledge, the restriction generates a closed convex set containing the solution in the Hilbert space, and then introduces the restriction of the amplitude boundary, thereby deriving an iterative formula for solving the unknown image, iterating from the initial estimate, and calculating the super-resolution image.

Convex set projection image reconstruction includes the following six steps:

(5.1) Establish an image degradation model: First, establish an image acquisition model that links the original high-resolution image and the low-resolution observation sequence. This model is an iterative formula for the unknown high-resolution image, which can be expressed as:

Among them, g _l (m ₁ , m ₂ ) is the l-th low-resolution image, x is the original high-resolution image, h _l represents the spatial point spread function (PSF), which is called the degradation function, and eta _l represents the addition function. sex noise. m ₁ × m ₂ is the pixel size of the low-resolution image, n ₁ × n ₂ is the pixel size of the high-resolution image, i ₁ × i ₂ is the number of local pixels in the high-resolution image, and its size is related to the spatial point spread function h _l consistent.

(5.2) Interpolate the reference frame to the target resolution: use linear interpolation method to interpolate the reference frame image, use the interpolated reference frame as the initial estimate of the high-resolution image, and use the a priori information provided by the reference frame (low The information contained in the resolution image) performs pixel value correction on a pixel-by-pixel basis.

(5.3) Define a closed convex set for each pixel of the high-resolution image: In the actual imaging process, noise is inevitably interfered, so the influence of additive noise is introduced based on the linear translation fuzzy constraint to make the noise The statistical properties are linked to the a priori boundary δ ₀ . If the additive noise is Gaussian distributed and its variance is σ _v , then δ ₀ is equal to cσ _v (c≥0), determined by a certain statistical confidence. Define the following closed convex set for each pixel of the reconstructed high-resolution image: Among them,/> Represents the residual difference between x(i ₁ , i ₂ ) and g(m ₁ , m ₂ ). Among them, g(m ₁ , m ₂ ) represents a single pixel in the low-resolution reference frame; x(i ₁ , i ₂ ) represents the regional pixel corresponding to g(m ₁ , m ₂ ) in the high-resolution initial estimation image . Its regional scale is [M ₁ , M ₂ ], which represents the template size of the spatial point spread function h (m ₁ , m ₂ ; i ₁ , i ₂ ). When x(i ₁ , i ₂ ) represents a real high-resolution image pixel, r ^(y) (m ₁ , m ₂ ) is consistent with the distribution of noise η _l (m ₁ , m ₂ ). The definition of a closed convex set represents the absolute value of the difference between the value of the reference frame image at the pixel point (m ₁ , m ₂ ) and the simulated imaging process at that point, which is limited to the preset boundary condition δ ₀ =cσ within _v .

(5.4) Calculate the residual pixel by pixel according to the degradation model to determine whether the boundary limit is exceeded. The pixel value of the reference frame image after correcting the interpolation must limit the pixel value not to exceed the threshold: Calculate the reference frame and the simulated low-resolution downsampled by the point spread function The residual difference between rate images (that is, the low-resolution image after the high-resolution initial estimation image is downsampled by the point spread function). Depending on whether the residual value exceeds the boundary limit, the gray value of the pixel is corrected according to the following formula.

The corrected pixel value needs to be within the image threshold interval [0, 255],

(5.5) Introduce all reference frames frame by frame for correction: read all reference frames in sequence, traverse all pixels in the reference frame, calculate the residual, and correct the residual value on the reference frame that has been amplified by interpolation, and correct the position correspondence for:

Where (m ₁ , m ₂ ) is the pixel coordinate position of the low-resolution reference frame, (Δx, Δy) is the sub-pixel displacement estimated by frequency domain registration, (i ₁ , i ₁ ) is the interpolated high-resolution reference frame The initial rate estimates the image pixel coordinates, k is the interpolation multiple, and the rounding function is in square brackets. After all reference frames are corrected according to the above rules, an iteration is completed;

(5.6) Loop multiple times for correction: Limited by the amplitude of the point spread function, one iteration cannot completely integrate the supplementary information in the reference frame into the high-resolution reference frame, so multiple iterations are needed to make different reference frames Different information of the same imaging object collected on the frame image is fully integrated.

The operation object of the above reconstruction process is the image blocks at the corresponding positions of the reference frame and the base frame. Different image blocks are registered according to the corresponding sub-pixel displacements, so that the pixel value correction position is more accurate. Therefore, the displacement estimation and image reconstruction are performed in blocks first. Then the super-resolution reconstructed block images are spliced according to their original positions to obtain the entire super-resolution image.

In order to obtain a low-resolution data set for reconstruction, the present invention changes the scanning mode of the OCT device and scans the same section of the retina within 1 second to obtain 25 B-scan images with an image size of 256×200pix. Since the imaging object of OCT is the retina of the living human eye with motion changes and is affected by random shot noise, the low-resolution image sequence in the time series contains different light and dark boundary information of the same imaging target imaged at different positions. This is also The key to being able to perform super-resolution reconstruction. However, it is necessary to compare the structural similarity of low-resolution image sequences. This is because the retina of the living human eye is affected by eye movements, and there will be structural deformation caused by local pulling of the retina. The low-resolution images taken at this time can no longer be regarded as a sequence of images imaging the same imaging object. Therefore, the obtained low-resolution images need to be screened according to structural similarity.

During screening, two indicators, Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index (SSIM), were used to compare the distortion and similarity of image sequences. According to Figure 2, it can be seen that the change trends of the two indicators are consistent. The peak signal-to-noise ratio is chosen for measurement in order to remove low-resolution images that are excessively affected by noise. The structural similarity index is selected for measurement in order to remove low-resolution images that are affected by eye movements and cause changes in the imaging target. The first frame of low-resolution image is selected as the base frame and the filtered 2 to 20 frames of images are used as reference frames to form the reconstructed low-resolution data set. At the same time, the SIFT algorithm is used to detect the overall displacement of the first 20 frames of images. As shown in Figure 3, there is a sub-pixel displacement relationship between the low-resolution reference frames numbered 2, 3, 6, 7, 8, 10, 11, 12, 13, 16, and 18 and the reference frame. Other reference frames do not. There is a direct sub-pixel displacement relationship, so the accuracy of imaging information for the same imaging object decreases. Therefore, in the reconstruction calculation, the low-resolution images numbered above are given priority for reconstruction.

The SIFT algorithm is used to perform angle correction on the above data set, and the direction of the reference frame is adjusted to be consistent with the reference frame to eliminate the angular influence in subsequent displacement estimation. The feature point parameters detected by the SIFT algorithm include a 128-dimensional feature point description vector, feature point coordinates, scale factor and main direction of rotation. As shown in Figure 4, 693 feature points were detected in the base frame and 700 feature points were detected in the reference frame. By calculating the Euclidean distance of the feature vectors in the two images, we find the feature points in the reference frame that match the feature points in the reference frame, and establish a corresponding relationship. All matching feature points are screened, some feature points with the highest matching degree are selected, statistics are made on the difference in the main directions of rotation, and the average rotation angle of the two images is calculated. Rotate the reference frame image according to the calculated average rotation angle for correction. Eliminate the angular influence in the process of estimating displacement using the frequency domain method.

Next, the base frame and reference frame are cropped into image blocks, and the displacement between the corresponding image blocks is calculated using the frequency domain registration method proposed by Vandewalle et al. The reason why the frequency domain method is used here for displacement estimation instead of directly based on the feature points detected by SIFT is that although the SIFT algorithm is resistant to noise interference, the internal points of the retinal OCT image layer structure are affected by additive noise and are easily used as feature points. detection, resulting in false matching, and the accuracy is lower than that of the frequency domain registration method based on image blocks. The reason for cropping the block image for matching is that multiple sub-pixel displacement vectors can be combined on the entire image to perform pixel value correction, which can add information to the correction position more accurately. Therefore, the low-resolution image is cut into 16 horizontal and vertical 4×4 image blocks according to the size, and frequency domain registration is performed between corresponding blocks. A single reference frame can use 16 sub-pixel displacement vectors to locally fine-tune the reference frame. In order to obtain clear high-resolution reconstructed images.

First, a mapping model between the original high-resolution image and the low-resolution reference frame image is established, which is the iterative formula for the unknown high-resolution image:

The reference frame is interpolated to the target resolution as an initial estimate of the high-resolution image. On this basis, the prior information provided by the reference frame is used to perform pixel value correction and fine-tuning pixel by pixel. The interpolation target multiple depends on the number of reference frame images for reconstruction, and follows k ² ≤ n (k is the interpolation magnification factor, n is the number of reference frame images), as shown in Figure 5. The number of reference frames used for 4x reconstruction is no less than 16 frames. , so 20 reference frame images are used, and the number of reference frames used for 3x reconstruction is no less than 9 frames, so 12 reference frame images are used.

Define a closed convex set for each pixel of the high-resolution image:

Here,/> Represents the residual difference between x(i ₁ , i ₂ ) and g(m ₁ , m ₂ ). Among them, g(m ₁ , m ₂ ) represents a single pixel in the low-resolution reference frame, and x(i ₁ , i ₂ ) represents the regional pixel corresponding to g(m ₁ , m ₂ ) in the high-resolution initial estimation image. , whose regional scale is [7, 7], which is the template size of the spatial point spread function h (m ₁ , m ₂ ; i ₁ , i ₂ ). The relationship between the scale of the spatial point spread function and the multiple of the reconstruction magnification resolution is M ₁ =M ₂ =2k-1, and the spatial point spread function is in the form of a two-dimensional Gaussian function. Inaccurate modeling of the spatial point spread function will lead to inaccurate estimation of pixel correction values, but the multiple iterations of the convex set projection algorithm can make up for this flaw. It mainly relies on the accuracy of a large amount of prior information, so it is very suitable for use with There is no two-dimensional point spread function in the OCT scanning system itself. It is considered here that the distribution of noise η _l (m ₁ , m ₂ ) conforms to the normal distribution N (0, 1), so the definition of a closed convex set means that the value of the reference frame image at the pixel point (m ₁ , m ₂ ) is consistent with the simulated imaging The absolute value of the difference between the values of the process at this point is limited to the preset boundary condition δ ₀ =σ _v =1 (c=1).

Next, calculate the residual between the reference frame and the simulated low-resolution image downsampled by the point spread function (i.e., the low-resolution image after the high-resolution initial estimate image is downsampled by the point spread function). According to the residual value Whether it exceeds the boundary limit, the gray value of the pixel is corrected according to the following formula, and the corrected pixel value is limited to the image threshold interval [0, 255].

Read all reference frames in sequence, traverse all pixels in the reference frames, and calculate the residuals. And the residual value is corrected on the reference frame that has been amplified by interpolation, and the correction position corresponds to:

Where (m ₁ , m ₂ ) is the pixel coordinate position of the low-resolution reference frame, (Δx, Δy) is the sub-pixel displacement estimated by frequency domain registration, (i ₁ , i ₂ ) is the interpolated high-resolution reference frame The initial rate estimates the image pixel coordinates, and k is the interpolation multiple. After all reference frames are corrected according to the above rules, an iteration is completed. Multiple iterations are performed to fully integrate different information of the same imaging object collected on different reference frame images.

The operation object of the above reconstruction process is the image blocks at the corresponding positions of the reference frame and the base frame. After the displacement estimation and image reconstruction are completed based on the single block image, the high-resolution block images that have been super-resolution reconstructed need to be spliced according to the original position. Obtain clear, realistic and detailed information-rich super-resolution images of the entire image.

Detailed description of the results obtained by the super-resolution reconstruction method:

Referring to Figure 1, this is the computational reconstruction process of the present invention: first, collect multiple frames of the same cross-section OCT B-scan image sequence and select one of the frames as the reference frame. After that, angle correction and displacement estimation are performed. The SIFT algorithm is used to extract feature points and feature vectors in the base frame and reference frame, and the Euclidean distance of the feature vectors is calculated. The minimum Euclidean distance is the matching point. According to the coordinates of the matching points in the base frame and the reference frame, the average rotation angle of the image is calculated. Rotate the reference frame according to the average rotation angle so that it is consistent with the direction of the base frame. Next, the image is cropped into image blocks, and the frequency domain method is used to estimate the sub-pixel displacement between the base frame and the corresponding image blocks of the reference frame for image reconstruction.

In the process of convex set projection image reconstruction, the reference frame image is first interpolated and enlarged to the target reconstruction multiple, and then the reference frame image is read in sequence and subjected to high-pass filtering. Only the high-frequency information of the low-resolution reference frame is fused to the final reconstructed high-resolution image to sharpen the features in the image. Calculate the spatial point spread function (PSF), and use the point spread function to downsample the interpolated and amplified reference frame. Calculate the residual between the downsampled initial frame and other reference frames, and correct some pixel values of the initial frame based on the residual value to achieve the purpose of fine-tuning to make the image clear, and limit the adjusted pixel value to the range of 0 to 255. All reference frames are read in sequentially to perform the above calculations. When all reference frames complete the above calculation process, an iterative process is completed. After multiple iterations, the reconstruction is completed, and the final reconstructed high-resolution image is output. The objects of the above reconstruction are the image blocks corresponding to the base frame and the reference frame, and the final entire image is formed by splicing the reconstructed image blocks according to their original positions.

See Figure 2(a) and Figure 2(b), which shows the peak signal-to-noise ratio and structural similarity of the reference frame compared with the baseline frame. The first frame of image collected is used as the reference frame, and the other frames are numbered sequentially as reference frames. It can be seen that starting from the 21st frame of the image, the peak signal-to-noise ratio and structural similarity are significantly lower than those of the previous 20 frames. The lower the image structural similarity, the less information about the same imaging object contained in the low-resolution image. The reconstructed image is susceptible to artifacts caused by erroneous information, so only the first 20 frames of images are used for reconstruction.

Refer to Figure 3, which shows the sub-pixel displacement estimation of the reference frame relative to the base frame. Taking the reference frame image of image 1 as the origin, the sub-pixel displacement estimation of images from frames 2 to 25 is shown in Figure 3. Draw a circle with a single pixel displacement as the radius, and there is a sub-pixel displacement relationship between the image represented by the serial number in the circle and the reference frame. Although the image represented by the serial number outside the circle differs from the reference frame by more than one pixel distance, its displacement includes sub-pixel distance, so the newly added pixels of the interpolation can still be corrected during the reconstruction process.

See Figure 4, which shows the SIFT algorithm feature point detection and matching. Among them: a) is the distribution of feature points detected in the reference frame; b) is the distribution of feature points detected in a reference frame; c) is the matching of feature points in the two frames of images. Record the coordinate positions of the matched feature points in the two frames of images, calculate the rotation angle and average it as the angle correction value to rotate the reference frame to the same direction as the reference frame.

See Figure 5, which shows a high-resolution image reconstructed by the convex set projection super-resolution reconstruction method based on multi-frame OCT images. Among them: a) is the base frame image; b) is the base frame image enlarged 4 times by bicubic interpolation; c) is the image reconstructed based on this method with a resolution increased 4 times, which is based on 20 frames of low-resolution images. The reconstruction result after 10 iterations; d) is the image reconstructed based on this method with a resolution increased by 3 times, which is the reconstruction result based on 12 frames of low-resolution images. Among them, the low-resolution images are filtered, and the 12 images with the closest displacement distance and the highest structural similarity to the reference frame are selected for reconstruction. It can be seen from the reconstruction results that the high-resolution image is enlarged compared to bicubic interpolation.

The high-resolution image with the same magnification reconstructed by the method of the present invention is clearer, the visual effect is significantly improved, the contrast is enhanced, the boundaries are obvious, and the high-frequency information is richer. The 3x reconstruction result is based on a higher quality reference frame, so some details are significantly enhanced and consistent with the reference frame, making the reconstruction result more realistic.

Claims (1)

1. The convex set projection super-resolution reconstruction method based on multi-frame optical coherence tomography images is characterized by reconstructing the high-frequency information in the high-resolution image by fusing the information contained in the multi-frame low-resolution image. To obtain clear, high-resolution images, the reconstruction method includes the following steps: (1) Collect optical coherence tomography low-resolution image sequences to change the scanning mode of the optical coherence tomography system so that it can quickly scan the same section to obtain multiple frames of B-ultrasound scanning images, and select one of the frames as the reference frame. , other images are used as reference frames for calculating and reconstructing high-resolution images; (2) Screen the multi-frame B-ultrasound scan images of the same section, compare the reference frames with the selected reference frames, and select the reference frame with the largest peak signal-to-noise ratio and structural similarity index value to fuse with the reference frame for reconstruction. , the selected base frame and the filtered reference frame constitute the reconstructed low-resolution data set; (3) Use the scale-invariant feature matching algorithm to perform angle correction on the data set obtained in step 2, adjust the direction of the reference frame to be consistent with the reference frame, and eliminate the angular influence in subsequent displacement estimation. Through scale-invariant feature matching The algorithm detects all feature points and corresponding feature vectors in the base frame and the reference frame. By calculating the Euclidean distance of the feature vectors in the two images, it finds the feature points in the reference frame that match the feature points in the base frame and establishes a corresponding relationship. , screen all matching feature points, select the feature point with the highest matching degree, count its coordinates, calculate the average rotation angle of the two images, and correct the reference frame image according to the calculated rotation angle; (4) Cut the base frame and reference frame into image blocks, and use the frequency domain registration method to calculate the displacement between the corresponding image blocks. In step 3, the angle correction has been performed based on the scale-invariant feature matching algorithm, so the calculated The angular influence in the displacement is considered to only exist in the x and y directions, in units of the number of pixels; Frequency domain registration principle calculation formula: The displacement of image f ₂ (x) in the spatial domain relative to the image f ₁ (x) corresponds to the phase shift of the image F ₂ (u) in the frequency domain relative to the image F ₁ (u). The displacement parameter Δx can be calculated by calculating the phase difference ∠ (F ₂ (u)/F ₁ (u)) is obtained, the relative displacement of each image block is recorded, and convex set projection super-resolution image reconstruction is performed based on the sub-pixel displacement (Δx, Δy); (5) The convex set projection image reconstruction process is mainly divided into the following six steps: (5.1) Establish an image degradation model: First, establish an image acquisition model that links the original high-resolution image and the low-resolution observation sequence. This model is an iterative formula for the unknown high-resolution image, which can be expressed as: Among them, g _l (m ₁ , m ₂ ) is the Frame low-resolution image, x(i ₁ , i ₂ ) is the original high-resolution image, h _l (m ₁ , m ₂ ; i ₁ , i ₂ ) represents the spatial point spread function (PSF), also known as drop Prime function, η _l (m ₁ , m ₂ ) represents additive noise; m ₁ × m ₂ is the pixel size of the low-resolution image, i ₁ × i ₂ is the number of local pixel points in the high-resolution image, and its size is related to the spatial point The diffusion function h _l is consistent; (5.2) Reference frame interpolation to target resolution: Use linear interpolation method to interpolate the reference frame image, use the interpolated reference frame as the initial estimate of the high-resolution image, and use the prior information provided by the reference frame pixel by pixel on this basis. Perform pixel value correction; (5.3) Define a closed convex set for each pixel of the high-resolution image: In the actual imaging process, noise is inevitably interfered, so the influence of additive noise is introduced based on the linear translation fuzzy constraint to make the noise The statistical properties are related to the a priori boundary δ _0. If the additive noise is Gaussian distributed and its variance is σ _v , then δ ₀ is equal to cσ _v (c≥0), which is determined by a certain statistical confidence. Each pixel of the reconstructed high-resolution image defines the following closed convex set: in, Represents the residual between x(i ₁ , i ₂ ) and g(m ₁ , m ₂ ), where g(m ₁ , m ₂ ) represents a single pixel in the low-resolution reference frame; x(i ₁ , i ₂ ) represents the regional pixel corresponding to g (m ₁ , m ₂ ) in the high-resolution initial estimation image, its regional scale is [M ₁ , M ₂ ], and represents the spatial point spread function h (m ₁ , m _{2 )} ; i ₁ , i ₂ ) template size, when x (i ₁ , i ₂ ) represents the real high-resolution image pixel point, r ^(y) (m ₁ , m ₂ ) and noise η _l (m ₁ , m ₂ ) distribution, the definition of a closed convex set represents the absolute value of the difference between the value of the reference frame image at the pixel point ( _mi , m ₂ ) and the simulated imaging process, the value at this point, which is limited to the preset The boundary condition δ ₀ =cσ _v is within; (5.4) Calculate the residual pixel by pixel according to the degradation model to determine whether the boundary limit is exceeded. The pixel value of the reference frame image after correcting the interpolation must limit the pixel value not to exceed the threshold: Calculate the reference frame and the simulated low-resolution downsampled by the point spread function The residual difference between rate images. Depending on whether the residual value exceeds the boundary limit, the gray value of the pixel is corrected according to the following formula: The corrected pixel value needs to be within the image threshold interval [0, 255], (5.5) Introduce all reference frames frame by frame for correction: read all reference frames in sequence, traverse all pixels in the reference frame, calculate the residual, and correct the residual value on the reference frame that has been amplified by interpolation, and correct the position correspondence for: Where (m ₁ , m ₂ ) is the pixel coordinate position of the low-resolution reference frame, (Δx, Δy) is the sub-pixel displacement estimated by frequency domain registration, (i ₁ , i ₂ ) is the interpolated high-resolution reference frame The initial rate estimates the image pixel coordinates, k is the interpolation multiple, and the rounding function is in square brackets. After all reference frames are corrected according to the above rules, an iteration is completed; (5.6) Loop multiple times for correction: Limited by the amplitude of the point spread function, one iteration cannot completely integrate the supplementary information in the reference frame into the high-resolution reference frame, so multiple iterations are needed to make different reference frames Different information of the same imaging object collected on the frame image is fully integrated. The operation object of the above reconstruction process is the image blocks at the corresponding positions of the reference frame and the base frame. Different image blocks are registered according to the corresponding sub-pixel displacements, so that the pixel value correction position is more accurate. Therefore, the displacement estimation and image reconstruction are performed in blocks first. Then the super-resolution reconstructed block images are spliced according to their original positions to obtain the entire super-resolution image.