A blood vessel image enhancement method based on multi-scale Gaussian filtering

Abstract

A kind of blood vessel image enhancement method based on multiscale Gaussian filter, comprising the following steps: step one: input endoscopic blood vessel image;Step two: the input blood vessel image is separated into RGB three channel planes;Step three: using multiscale Gaussian filter and bilateral filter to filter G component to obtain detail layer and base layer;Step four: detail layer is multiplied by a gain coefficient, then merged with base layer to obtain two enhancement results;Step five: two enhancement results are converted to HIS space, I component is carried out 4-layer wavelet decomposition to obtain 1 low frequency and 12 high frequency components, and the rule of high frequency weighted average low frequency modulus value maximum is used to fuse to obtain enhanced I component, then merged with HS, converted to RGB space to obtain the final enhancement result.The present application can carry out different adaptive enhancement to different levels of blood vessels to obtain more abundant blood vessel images, which can assist doctors in early diagnosis of cancer.

Description

Technical Field

[0001] This invention belongs to the field of auxiliary medical diagnosis and image enhancement, and relates to a method for enhancing blood vessel images. Background Technology

[0002] After investigating the specific tissue information of cancerous areas, it was found that cancerous areas usually lead to changes in tissue composition. Studies have shown that there are significant differences in hemoglobin content between normal and cancerous tissues. Hemoglobin is abundant in the blood, and the blood vessels in malignant tumor tissues are more abundant, with relatively disordered blood vessel arrangement, increased tortuosity, uneven diameter, and more branches. Furthermore, there is a significant difference in blood oxygen saturation between normal and cancerous tissues. Therefore, significantly enhancing vascular information in suspected cancerous areas can help doctors make better diagnoses. Early-stage cancers often have no obvious clinical symptoms or signs, and endoscopically they mainly show minor changes in mucosal structure. Therefore, enhancing vascular features to amplify these diagnostic features and visually observable shape characteristics can help doctors shorten the diagnosis time for early-stage cancers, thus preventing the need for optimal treatment.

[0003] Currently, the mainstream proprietary technologies for enhancing endoscopic vascular images include NBI (narrow band imaging), FICE (flexible spectral imaging color enhancement), and I-Scan. All three aim to enhance vascular features for real-time processing. NBI is a hardware image enhancement technology that utilizes narrow band illumination, while FICE and I-Scan are both electronic staining techniques. FICE reconstructs a color image by calculating and selecting the reflection intensity at several specific wavelengths to enhance the target of interest. It requires rigorous calibration of the endoscope system before use, and the color of the enhanced image differs significantly from the true color. I-Scan technology not only allows for the selection of different wavelength combinations for image display but also incorporates surface enhancement and contrast enhancement methods, although its algorithm is more complex. Summary of the Invention

[0004] To overcome the shortcomings of existing technologies, this invention provides a vascular image enhancement method based on multi-scale Gaussian filtering, which can perform different adaptive enhancements on blood vessels at different levels to obtain richer vascular images, thereby assisting doctors in the early diagnosis of cancer.

[0005] The technical solution adopted by this invention to solve its technical problem is:

[0006] A method for enhancing blood vessel images based on multi-scale Gaussian filtering includes the following steps:

[0007] Step 1: Input endoscopic vascular images;

[0008] Step 2: Separate the input blood vessel image into three channel planes: RGB.

[0009] Step 3: Filter the G component using a multi-scale Gaussian filter and a bilateral filter to obtain the detail layer and the base layer;

[0010] Step 4: Multiply the detail layer by a gain factor based on the result of Step 3, and then merge it with the base layer to obtain two enhancement results;

[0011] Step 5: Convert the two enhancement results from Step 4 to the HIS space, perform 4-level wavelet decomposition on the I component to obtain 1 low-frequency and 12 high-frequency components, and use the rule of maximizing the modulus value of the low-frequency component by weighting the high-frequency component to obtain the enhanced I component. Then merge it with the HIS space and convert it to the RGB space to obtain the final enhancement result.

[0012] Furthermore, in step three, the G component is enhanced using a multi-scale Gaussian filter, as shown in the following formula:

[0013] f(x,y)=b(x,y)+d(x,y)

[0014] Where b(x,y) is the base layer, d(x,y) is the detail layer, and f(x,y) is the input graph.

[0015] A two-dimensional Gaussian kernel function g(x,y) with mean 0 and variance σ is:

[0016]

[0017] How the grassroots receive it:

[0018] b(x,y)=f(x,y)⊙g(x,y)

[0019] ⊙ represents convolution;

[0020] The above formula is used to perform multi-scale augmentation calculations by substituting the values ​​w1 = 0.5, w2 = 0.5, and w3 = 0.25 for the three scales.

[0021] b1(x,y)=f(x,y)⊙g1(x,y); b2(x,y)=f(x,y)⊙g2(x,y); b3(x,y)

[0022] = f(x,y)⊙g3(x,y)

[0023] d1(x,y)=f(x,y)-b1(x,y); d2(x,y)=b1(x,y)-b2(x,y); d3(x,y)

[0024] = b2(x,y) - b3(x,y)

[0025] D(x,y)=(1-w1*sgn(d1(x,y)))*d1(x,y)+w2*d2(x,y)+w3

[0026] *d3(x,y)

[0027] Where bi(x,y) is the i-th base layer, di(x,y) is the i-th detail layer, and D(x,y) is the final enhanced detail layer; (x,y) refers to the pixel point at (x,y) coordinates, i = 1, 2, 3.

[0028] The enhanced detail layer D(x,y) is then added to the original image to obtain a new enhanced image output1.

[0029] Furthermore, in step four, the G component of the original image is enhanced by bilateral filtering. The G component of the original image and the bilateral filter are convolved to obtain the base layer. Then, the base layer is subtracted from the original image to obtain the detail layer. Finally, the detail layer is multiplied by a gain coefficient to obtain the enhanced detail layer, which is then merged with the base layer to obtain the enhanced image output2, as shown in the following formula:

[0030] b(x,y)=f(x,y)⊙Q(x,y,s,l)

[0031] d(x,y)=f(x,y)-b(x,y)

[0032] z(x,y)=b(x,y)+k*d(x,y)

[0033] Where Q(x,y,s,l) is the convolution kernel of the bilateral filter, k is the gain coefficient, b(x,y) is the base layer, d(x,y) is the detail layer, and z(x,y) is the final enhancement result; s and l are the detail parameters of the bilateral filter, which control the weight of the range filter and the spatial domain sensitivity of the edge-preserving property, respectively.

[0034] Furthermore, in step five, output1 and output2 are converted to the HIS space, and their I components are decomposed into 4-level wavelet decomposition to obtain 1 low-frequency and 12 high-frequency components. These components are then fused according to the rule of maximizing the modulus of the low-frequency components and weighted averaging the high-frequency components to obtain the merged I component. This merged I component is then combined with the original HS component. Finally, the result is converted from the HIS space to the RGB space to obtain the final enhanced image output, as shown in the following formula:

[0035] Continuous wavelet transform formula:

[0036]

[0037] Where a is the translation parameter, b is the scaling parameter, t is the time parameter in the frequency domain, and x(t) is the function for transforming the image to the frequency domain;

[0038] The ψ(t) used is the Halm wavelet function, and the corresponding scaling function is:

[0039]

[0040] Formula for inverse continuous wavelet transform:

[0041]

[0042] ψ1(t) is the dual function of ψ(t) and the dual function must satisfy the following:

[0043]

[0044] Decomposing output1 and output2 yields two sets of wavelet decomposition coefficients as follows:

[0045]

[0046] When j=1, the original image is input on the right, and the high-frequency and low-frequency components are obtained:

[0047]

[0048] The high-frequency components are:

[0049]

[0050] A, B, and C represent the high-frequency components in the horizontal, vertical, and diagonal directions, respectively.

[0051] The low-frequency components are:

[0052]

[0053] This yields high-frequency and low-frequency components. The high-frequency components are weighted and averaged, while the low-frequency components are fused by taking the largest modulus. Finally, the enhanced output is obtained through inverse transformation.

[0054] The technical concept of this invention is as follows: After analyzing the hypoxia and anoxic absorption spectra of proteins, the g-component of the image is specifically processed using three-scale Gaussian filtering based on Retinex theory to obtain three different scale detail layers. These layers are then added together with different weights to obtain multi-scale detail layer enhancement. Next, bilateral filtering is used to decompose the image into detail layers and base layers. Each detail layer is multiplied by a gain coefficient and then merged into two enhanced images with different detail information. Then, the I-component of the HIS of the two enhanced images is decomposed using wavelet decomposition at four levels to obtain one low-frequency and twelve high-frequency coefficients. Finally, the high-frequency coefficients are weighted and averaged, while the low-frequency coefficients are maximized to obtain a wavelet fusion image. Finally, the HIS is converted to an RGB image. In summary, subjectively speaking, this method, observed by the naked eye, significantly enhances and better displays blood vessels and some easily overlooked microvessels while maintaining the original image's color accuracy, greatly enriching vascular information. Furthermore, the enhancement results outperform NBI, FICE, and I-SCAN in various objective evaluation metrics.

[0055] The vascular enhancement method of this invention mainly combines the detailed information obtained from multi-scale Gaussian filtering with the vascular boundary information enhanced by bilateral filtering to obtain more information about the blood vessels. This fusion enriches the vascular information. Visual illustrations using oral vascular images acquired via laboratory endoscopy, along with comparisons of enhanced images from three mainstream and widely used existing technologies (NBI, FICE, and I-SCAN) with the enhanced images of this invention using three evaluation metrics (entropy, signal-to-noise ratio, and average gradient), demonstrate that this invention outperforms the existing technologies.

[0056] The main benefits of this invention are that it can adaptively enhance blood vessels at different levels to obtain richer vascular images, thereby assisting doctors in the early diagnosis of cancer. Attached Figure Description

[0057] Figure 1 This is a flowchart of a blood vessel image enhancement method based on multi-scale Gaussian filtering.

[0058] Figure 2 This is an image of oral blood vessels acquired using a laboratory endoscope 1.

[0059] Figure 3 The image shows oral blood vessels captured using a laboratory endoscope. Detailed Implementation

[0060] The present invention will now be further described with reference to the accompanying drawings.

[0061] Reference Figures 1-3A method for enhancing blood vessel images based on multi-scale Gaussian filtering includes the following steps:

[0062] Step 1: Input endoscopic vascular images;

[0063] Step 2: Separate the input blood vessel image into three channel planes: RGB.

[0064] Step 3: Filter the G component using a multi-scale Gaussian filter and a bilateral filter to obtain the detail layer and the base layer;

[0065] In step three, the G component is enhanced using a multi-scale Gaussian filter, as shown in the following formula:

[0066] f(x,y)=b(x,y)+d(x,y)

[0067] Where b(x,y) is the base layer, d(x,y) is the detail layer, and f(x,y) is the input graph.

[0068] A two-dimensional Gaussian kernel function g(x,y) with mean 0 and variance σ is:

[0069]

[0070] How the grassroots receive it:

[0071] b(x,y)=f(x,y)⊙g(x,y)

[0072] ⊙ represents convolution;

[0073] Substituting the above formula with the three scales w1 = 0.5, w2 = 0.5, and w3 = 0.25, we can perform multi-scale augmentation calculations.

[0074] b1(x,y)=f(x,y)⊙g1(x,y); b2(x,y)=f(x,y)⊙g2(x,y); b3(x,y)

[0075] = f(x,y)⊙g3(x,y)

[0076] d1(x,y)=f(x,y)-b1(x,y); d2(x,y)=b1(x,y)-b2(x,y); d3(x,y)

[0077] = b2(x,y) - b3(x,y)

[0078] D(x,y)=(1-w1*sgn(d1(x,y)))*d1(x,y)+w2*d2(x,y)+w3

[0079] *d3(x,y)

[0080] Where bi(x,y) is the i-th base layer, di(x,y) is the i-th detail layer, and D(x,y) is the final enhanced detail layer; (x,y) refers to the pixel point at (x,y) coordinates, i = 1, 2, 3.

[0081] The enhanced detail layer D(x,y) is then added to the original image to obtain a new enhanced image output1.

[0082] Step 4: Multiply the detail layer by a gain factor based on the result of Step 3, and then merge it with the base layer to obtain two enhancement results;

[0083] In step four, the G component of the original image is enhanced by bilateral filtering. The G component of the original image and the bilateral filter are convolved to obtain the base layer. Then, the base layer is subtracted from the original image to obtain the detail layer. Finally, the detail layer is multiplied by a gain coefficient to obtain the enhanced detail layer, which is then merged with the base layer to obtain the enhanced image output2, as shown in the following formula:

[0084] b(x,y)=f(x,y)⊙Q(x,y,s,l)

[0085] d(x,y)=f(x,y)-b(x,y)

[0086] z(x,y)=b(x,y)+k*d(x,y)

[0087] Where Q(x,y,s,l) is the convolution kernel of the bilateral filter, k is the gain coefficient, b(x,y) is the base layer, d(x,y) is the detail layer, and z(x,y) is the final enhancement result; s and l are the detail parameters of the bilateral filter, which control the weight of the range filter and the spatial domain sensitivity of the edge-preserving property, respectively.

[0088] Step 5: Convert the two enhancement results from Step 4 to the HIS space, perform 4-level wavelet decomposition on the I component to obtain 1 low-frequency and 12 high-frequency components, and use the rule of maximizing the modulus value of the low-frequency component by weighting the high-frequency component to obtain the enhanced I component. Then merge it with the HIS space and convert it to the RGB space to obtain the final enhancement result.

[0089] In step five, output1 and output2 are converted to the HIS space, and their I components are decomposed into 4-level wavelet decomposition to obtain 1 low-frequency and 12 high-frequency components. The I components are then fused according to the rule of taking the maximum modulus value of the low frequency and taking the weighted average value of the high frequency to obtain the merged I component. This merged I component is then combined with the original HS, and finally converted from the HIS space to the RGB space to obtain the final enhanced image output, as shown in the following formula:

[0090] Continuous wavelet transform formula:

[0091]

[0092] Where a is the translation parameter, b is the scaling parameter, t is the time parameter in the frequency domain, and x(t) is the function for transforming the image to the frequency domain;

[0093] The ψ(t) used is the Halm wavelet function, and the corresponding scaling function is:

[0094]

[0095] Formula for inverse continuous wavelet transform:

[0096]

[0097] ψ1(t) is the dual function of ψ(t) and the dual function must satisfy the following:

[0098]

[0099] Decomposing output1 and output2 yields two sets of wavelet decomposition coefficients as follows:

[0100]

[0101] When j=1, the original image is input on the right, and the high-frequency and low-frequency components are obtained:

[0102]

[0103] The high-frequency components are:

[0104]

[0105] A, B, and C represent the high-frequency components in the horizontal, vertical, and diagonal directions, respectively.

[0106] The low-frequency components are:

[0107]

[0108] This yields high-frequency and low-frequency components. The high-frequency components are weighted and averaged, while the low-frequency components are fused by taking the largest modulus. Finally, the enhanced output is obtained through inverse transformation.

[0109] In this embodiment, Figure 1 The flowchart shown is a blood vessel image enhancement algorithm of the present invention. The input blood vessel image is processed until an enhanced image with richer content and easier for the human eye to perceive is obtained. The detailed operation is as follows:

[0110] Reference Figure 1 Image of oral blood vessels acquired using a laboratory endoscope;

[0111] First, the input image is separated into three components: RGB.

[0112] The G component is enhanced using a multi-scale Gaussian filter, and the steps are shown in the following formula:

[0113] f(x,y)=b(x,y)+d(x,y)

[0114] Where b(x,y) is the base layer, d(x,y) is the detail layer, and f(x,y) is the input graph.

[0115] A two-dimensional Gaussian kernel function g(x,y) with mean 0 and variance σ is:

[0116]

[0117] How the grassroots receive it:

[0118] b(x,y)=f(x,y)⊙g(x,y)

[0119] ⊙ represents convolution;

[0120] The above formula is used to perform multi-scale augmentation calculations by substituting the values ​​w1 = 0.5, w2 = 0.5, and w3 = 0.25 for the three scales.

[0121] b1(x,y)=f(x,y)⊙g1(x,y); b2(x,y)=f(x,y)⊙g2(x,y); b3(x,y)=f(x,y)⊙g3(x,y)

[0122] d1(x,y)=f(x,y)-b1(x,y); d2(x,y)=b1(x,y)-b2(x,y); d3(x,y)=b2(x,y)-b3(x,y)

[0123] D(x,y)=(1-w1*sgn(d1(x,y)))*d1(x,y)+w2*d2(x,y)+w3*d3(x,y)

[0124] Where bi(x,y) is the i-th base layer, di(x,y) is the i-th detail layer, and D(x,y) is the final enhanced detail layer; (x,y) refers to the pixel point at (x,y) coordinates, i = 1, 2, 3.

[0125] The enhanced detail layer D(x,y) is then added to the original image to obtain a new enhanced image output1;

[0126] Next, the G component of the original image is enhanced by bilateral filtering. The G component of the original image and the bilateral filter are convolved to obtain the base layer. Then, the base layer is subtracted from the original image to obtain the detail layer. Finally, the detail layer is multiplied by a gain coefficient to obtain the enhanced detail layer, which is then merged with the base layer to obtain the enhanced image output2. The process is shown in the following formula:

[0127] b(x,y)=f(x,y)⊙Q(x,y,s,l)

[0128] d(x,y)=f(x,y)-b(x,y)

[0129] z(x,y)=b(x,y)+k*d(x,y)

[0130] Where Q(x,y,s,l) is the convolution kernel of the bilateral filter, k is the gain coefficient, b(x,y) is the base layer, d(x,y) is the detail layer, and z(x,y) is the final enhancement result; s and l are the detail parameters of the bilateral filter, which control the weight of the range filter and the spatial domain sensitivity of the edge-preserving property, respectively.

[0131] Next, output1 and output2 are converted to HIS space, and their I components are decomposed into 4-level wavelet decomposition to obtain 1 low-frequency and 12 high-frequency components. The I components are then fused according to the rule of taking the maximum modulus value of the low frequency and taking the weighted average value of the high frequency to obtain the merged I component. This merged I component is then combined with the original HS, and finally converted from HIS space to RGB space to obtain the final enhanced image output. The main process is shown in the following formula:

[0132] Continuous wavelet transform formula:

[0133]

[0134] Where a is the translation parameter, b is the scaling parameter, t is the time parameter in the frequency domain, and x(t) is the function for transforming the image to the frequency domain;

[0135] The ψ(t) used is the Halm wavelet function, and the corresponding scaling function is:

[0136]

[0137] Formula for inverse continuous wavelet transform:

[0138]

[0139] ψ1(t) is the dual function of ψ(t) and the dual function must satisfy the following:

[0140]

[0141] Based on the above formula, output1 and output2 are decomposed into two sets of wavelet decomposition coefficients as follows (when j=1, the original image is input on the right):

[0142]

[0143] High-frequency and low-frequency components are obtained:

[0144]

[0145] The high-frequency components are:

[0146]

[0147] A, B, and C represent the high-frequency components in the horizontal, vertical, and diagonal directions, respectively.

[0148] The low-frequency components are:

[0149]

[0150] This yields high-frequency and low-frequency components. A weighted average is used for the high-frequency components, and the method of maximizing the magnitude is used for the low-frequency components to fuse the high and low frequencies. Finally, an inverse transform is performed to obtain the enhanced output.

[0151] Figure 2 and Figure 3 The images of oral blood vessels from two different volunteers, collected under a laboratory endoscope as an embodiment of the present invention (the enhanced images are shown below), demonstrate that the enhanced blood vessel information is richer and more visible to the human eye, and some subtle blood vessel details are better displayed.

[0152] Tables 1 to 3 show that the algorithm of this invention outperforms the mainstream endoscopic enhancement technologies (FICE, I-SCAN, NBI) in three evaluation metrics (entropy value, psnr (peak signal-to-noise ratio), and AG (average gradient)). This demonstrates that the algorithm of this invention is superior to the three widely used proprietary technologies in terms of entropy value, psnr (peak signal-to-noise ratio), and AG (average gradient). This further illustrates that vascular-related information is better presented not only by visual observation but also by objective evaluation metrics.

[0153] Table 1 compares the algorithm of this invention with that of FICE.

[0154]

[0155] Table 2 compares the algorithm of this invention with I-SCAN.

[0156]

[0157] Table 3 compares the algorithm of this invention with NBI.

[0158]

[0159] The embodiments described in this specification are merely examples of implementations of the inventive concept and are for illustrative purposes only. The scope of protection of this invention should not be considered limited to the specific forms described in these embodiments; rather, it extends to equivalent technical means conceived by those skilled in the art based on the inventive concept.

Claims

1. A method for enhancing blood vessel images based on multi-scale Gaussian filtering, characterized in that, The method includes the following steps: Step 1: Input endoscopic vascular images; Step 2: Separate the input blood vessel image into three channel planes: RGB. Step 3: Enhance the G component using multi-scale Gaussian filtering, as shown in the following formula: f(x,y)=b(x,y)+d(x,y) Where b(x,y) is the base layer, d(x,y) is the detail layer, and f(x,y) is the input graph. A two-dimensional Gaussian kernel function g(x,y) with mean 0 and variance σ is: How the grassroots receive it: b(x,y)=f(x,y)⊙g(x,y) ⊙ represents convolution; Substituting the above formula with the empirically derived values ​​of w1 = 0.5, w2 = 0.5, and w3 = 0.25 for the three scales yields the best results for multi-scale augmentation calculations. b1(x,y)=f(x,y)⊙g1(x,y); b2(x,y)=f(x,y)⊙g2(x,y); b3(x,y)=f(x,y)⊙g3(x,y) d1(x,y)=f(x,y)-b1(x,y); d2(x,y)=b1(x,y)-b2(x,y); d3(x,y)=b2(x,y)-b3(x,y) D(x,y)=(1-w1*sgn(d1(x,y)))*d1(x,y)+w2*d2(x,y)+w3*d3(x,y) Where bi(x,y) is the i-th base layer, di(x,y) is the i-th detail layer, and D(x,y) is the final enhanced detail layer; (x,y) refers to the coordinates of the pixel, i = 1, 2, 3; The enhanced detail layer D(x,y) is then added to the original image to obtain a new enhanced image output1; Step 4: Enhance the original image using bilateral filtering on the G component. Convolve the original image's G component and the bilateral filter to obtain the base layer. Subtract the base layer from the original image to obtain the detail layer. Finally, multiply the detail layer by a gain coefficient to obtain the enhanced detail layer. Combine the enhanced detail layer with the base layer to obtain the enhanced image output2, as shown in the following formula: b(x,y)=f(x,y)⊙Q(x,y,s,l) d(x,y)=f(x,y)-b(x,y) z(x,y)=b(x,y)+k*d(x,y) Where Q(x,y,s,l) is the convolution kernel of the bilateral filter, k is the gain coefficient, b(x,y) is the base layer, d(x,y) is the detail layer, and z(x,y) is the final enhancement result; s and l are the detail parameters of the bilateral filter, which control the weight of the range filter and the spatial domain sensitivity of the edge-preserving property, respectively. Step 5: Convert output1 and output2 to HIS space, perform 4-level wavelet decomposition on their I components to obtain 1 low-frequency and 12 high-frequency components, and fuse them according to the rule of taking the maximum modulus value of the low frequency and taking the weighted average value of the high frequency to obtain the merged I component, then merge it with the original HS, and finally convert it from HIS space to RGB space to obtain the final enhanced image output. Continuous wavelet transform formula: Where a is the translation parameter, b is the scaling parameter, t is the time parameter in the frequency domain, and x(t) is the function for transforming the image to the frequency domain; The ψ(t) used is the Halm wavelet function, and the corresponding scaling function is: Formula for inverse continuous wavelet transform: ψ1(t) is the dual function of ψ(t) and the dual function must satisfy the following: Decomposing output1 and output2 yields two sets of wavelet decomposition coefficients as follows: When j=1, the original image is input on the right, and the high-frequency and low-frequency components are obtained: The high-frequency components are: A, B, and C represent the high-frequency components in the horizontal, vertical, and diagonal directions, respectively. Low-frequency components are This yields high-frequency and low-frequency components. The high-frequency components are weighted and averaged, while the low-frequency components are fused by taking the largest modulus. Finally, the enhanced output is obtained through inverse transformation.

Images