Method for enhancing quality of low-dose CBCT (cone beam computed tomography) image

Abstract

The invention provides a low-dose CBCT (cone beam computed tomography) image quality enhancement method, which more comprehensively considers the noise type in an image, can effectively suppress noise and can keep image details. According to the method, a three-dimensional low-dose CBCT image is firstly sliced to generate a single 2D image, and then noise estimation is carried out on the single 2D image by using a composite noise model; in the noise reduction stage, after removal of a spatial domain is achieved for additive noise, multiplicative noise is converted into additive noise and then separated into a high-frequency component and a low-frequency component, and a filtering method is designed in a targeted mode for precise noise reduction; after the original data structure of the image after noise reduction is recovered, edge recovery is carried out on edge pixels lost in noise reduction, the data after edge recovery is reconstructed into a three-dimensional image, and consistency adjustment is carried out on the three-dimensional image data based on the brightness between layers, so that the quality enhancement of the image is realized.

Description

Technical Field

[0001] The present invention relates to the technical field of image processing, and in particular to a method for enhancing low-dose CBCT image quality. Background Art

[0002] In dentistry, particularly in orthodontics and implantology, high-resolution oral cavity images are often acquired using cone-beam computed tomography (CBCT) technology for treatment planning. High-resolution CBCT images typically require a higher radiation dose, which also increases the risk of cancer. While low-dose CBCT can ensure patient safety, it also results in higher noise in the acquired images, compromising image quality. To ensure usable images, image denoising and quality enhancement are essential.

[0003] Existing low-dose CBCT image quality enhancement methods are primarily divided into traditional image processing methods and deep learning approaches. Most traditional methods assume that noise in low-dose images is additive Gaussian noise and then perform denoising on the additive noise. However, low-dose CBCT images often have more complex noise components. For example, the noise may be a composite type of speckle noise, which includes not only additive noise but also multiplicative noise. Assuming only one type of noise will affect the image enhancement effect. Summary of the Invention

[0004] To address the problem of poor low-dose CBCT image enhancement in the prior art, the present invention provides a method for low-dose CBCT image quality enhancement that more comprehensively considers the types of noise in the image, effectively suppressing noise while preserving image details.

[0005] The technical solution of the present invention is as follows: a method for enhancing low-dose CBCT image quality, characterized in that it comprises the following steps:

[0006] S1: Acquire CBCT noise images;

[0007] Low-dose head model projection data were collected, and the corresponding noisy CBCT three-dimensional images were reconstructed using the FDK algorithm.

[0008] S2: extracting a 2D axial slice based on the CBCT 3D image, which is recorded as the image to be processed;

[0009] S3: For each image to be processed, analyze the noise distribution type and construct a composite noise distribution model:

[0010] y=n multi *x+n add ;

[0011] Among them, x represents the clean image, y represents the noisy image, and n multi represents multiplicative noise, n add represents additive noise;

[0012] S4: Based on the noise distribution model corresponding to the image to be processed, after removing the additive noise in the image to be processed, the multiplicative noise in the image is converted into additive noise to obtain: a preliminarily processed image;

[0013] S5: Decomposing the preliminarily processed image into a high-frequency component and a low-frequency component using discrete wavelet transform, which are respectively recorded as: a high-frequency component image and a low-frequency component image;

[0014] S6: De-noising the low-frequency component image by using a low-frequency adaptive Wiener filtering method to obtain the de-noised low-frequency feature I ld , after replacing the corresponding pixel points, the denoised low-frequency component data is obtained;

[0015] Low-frequency characteristics after noise reduction I ld The calculation method is:

[0016]

[0017] Among them, I ld Represents the low-frequency characteristics after noise reduction, μ l Represents the pixel mean value of the local area weighted by w in the low-frequency component; σ l 2 represents the variance of pixels in the local area; ν lf 2 is a hyperparameter, indicating the noise variance; max() is the maximum value calculation; I n Indicates the pixel value of each pixel in the traversal low-frequency component;

[0018] ν lf 2 The calculation method is:

[0019]

[0020] Where N represents the local area window, σ lf 2 Represents the set of all local pixel variances in the low-frequency component; mean() is the mean calculation;

[0021] S7: performing noise reduction processing on the high-frequency component image by improving the high-frequency adaptive wavelet threshold denoising method, specifically comprising the following steps:

[0022] a1: Combined with noise estimation to generate adaptive soft threshold T hi ;

[0023]

[0024] Where, T hi represents the soft threshold of the i-th high-frequency component, v hi 2 Represents the noise variance of the i-th high-frequency component, i represents the sequence number of the high-frequency component separated by wavelet transform, σ i 2 represents the variance of the i-th high-frequency component;

[0025]

[0026] Where, median(|H i |) represents the middle pixel value in the i-th high-frequency component, σ hi 2 Represents the set of all local pixel variances in the i-th high-frequency component, N is the local area window, and M represents the size of the i-th high-frequency component;

[0027] a2: Transform the wavelet threshold function through soft thresholding to obtain the wavelet coefficient w after soft thresholding λ,hi ;

[0028]

[0029] Where sign() is the sign function, w hi represents the wavelet coefficient of the i-th high-frequency component; w λ,hi Represents the wavelet coefficients after threshold function processing;

[0030] a3: wavelet coefficient w after soft threshold processing λ,hi Perform signal reconstruction to generate denoised high-frequency component data;

[0031] S8: performing an inverse transformation on the denoised low-frequency component data and the denoised high-frequency component data to restore the previous data structure, thereby obtaining a denoised image to be processed;

[0032] S9: performing edge restoration processing on the denoised image to be processed to obtain an edge-restored image to be processed;

[0033] The edge restoration process includes: image smoothing and image sharpening performed in sequence;

[0034] S10: stitching all the 2D edge-restored images to be processed into 3D data again according to the order of the previous slices to obtain a stitched CBCT 3D image; performing 3D data equalization processing on each 2D image layer of the stitched CBCT 3D image to obtain a final quality-enhanced image;

[0035] The three-dimensional data balancing process includes: brightness consistency adjustment and balancing the image brightness between layers.

[0036] It is further characterized by:

[0037] In step S4, the additive noise in the image to be processed is removed by using an adaptive median filtering method;

[0038] In step S4, the multiplicative noise is converted into additive noise using a logarithmic transformation;

[0039] In step S5, the discrete wavelet transform is used to decompose the signal into high-frequency components and low-frequency components, thereby obtaining one low-frequency component and three high-frequency components.

[0040] In step S6, the local area window N is: N=3×3; the corresponding local mean calculation weight w is calculated as follows:

[0041]

[0042] In step S8, the inverse transformation of the denoised low-frequency component data and the denoised high-frequency component data includes: sequentially performing an inverse wavelet transform and an inverse logarithmic transform;

[0043] In the edge restoration process, the image smoothing process uses the first step of the basic estimation of the BM3D algorithm to perform image noise reduction, and uses CUDA for parallel acceleration; the image sharpening uses the USM algorithm to enhance the edge information that was smoothed in the previous noise reduction process;

[0044] The three-dimensional data equalization processing specifically includes the following calculations:

[0045] I c =min(max(r·I′ filter ,min(I′ noise )),max(I′ noise ))+min(I noise );

[0046] Where, I c Indicates the image after brightness adjustment, I filter is the image after noise reduction filtering, I noise is the source noise image;

[0047] I' filter and I' noise is the median value and is calculated as follows:

[0048] I filter =I filter -min(I noise ;

[0049] I′noise =I noise -min(I noise ;

[0050] r represents the weight coefficient of each layer of image, and the calculation method is as follows:

[0051]

[0052] The present application provides a method for enhancing the quality of low-dose CBCT images, which first evaluates the noise distribution of the image to be processed based on a composite noise distribution model to ensure that the noise composition of the image to be processed can be more faithfully represented, thereby ensuring that the subsequent noise reduction processing can be more targeted and effectively improve the noise reduction effect; for the noise of the image to be processed, the additive noise is first removed, and then the multiplicative noise is converted into additive noise, and then the converted additive noise is decomposed into high-frequency components and low-frequency components through discrete wavelet transform. The low-frequency component mainly contains the main visual information, and the high-frequency component mainly contains the edge details and concentrated noise information. Classifying the high-frequency component and the low-frequency component separately for noise reduction can take into account both the noise reduction effect and the edge preservation. edge details requirements; in this method, a low-frequency adaptive Wiener filtering method is proposed for low-frequency components, and the noise level is calculated more accurately by using local statistical methods. It is superior to the traditional Wiener filtering method in denoising low-frequency components; for high-frequency components, this method proposes an improved high-frequency adaptive wavelet threshold denoising method, which calculates the thresholds of different high-frequency components through noise estimation to achieve the function of adaptive denoising. Compared with the traditional wavelet threshold denoising method, it can effectively suppress noise and retain image details; through the combination of the low-frequency adaptive Wiener filtering method and the improved high-frequency adaptive wavelet threshold denoising method, the denoising effect of multiplicative noise is effectively improved, thereby improving the quality enhancement effect of low-dose CBCT images. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 A flow chart of the method for low-dose CBCT image quality enhancement in this application;

[0054] Figure 2 Schematic diagram of noise processing process;

[0055] Figure 3 A flow chart of the median filtering method used in additive noise processing;

[0056] Figure 4 Schematic diagram of edge recovery processing flow. DETAILED DESCRIPTION

[0057] like Figure 1 As shown, the present application includes a method for low-dose CBCT image quality enhancement, which includes three stages: a noise reduction stage, an edge recovery stage, and a three-dimensional data equalization stage.

[0058] The noise reduction stage involves acquiring CBCT noise images, estimating noise models, and performing a noise processing module. Acquiring CBCT noise images aims to generate a single 2D image; estimating noise models assesses the noise type; and the noise processing module uses a series of targeted high- and low-frequency filtering methods to improve noise reduction performance. Specifically, the following detailed steps are included.

[0059] S1: Acquire CBCT noise images;

[0060] Low-dose head model projection data are acquired, and the corresponding noisy CBCT three-dimensional image is reconstructed using the FDK algorithm to provide a single 2D image for subsequent processing.

[0061] S2: Extracting 2D axial slices based on CBCT 3D images, denoted as "images to be processed." A specific method for extracting 2D axial slices from CBCT 3D images is implemented based on existing technology.

[0062] S3: For each image to be processed, analyze the noise distribution type and build a composite noise distribution model:

[0063] y=n multi *x+n add ;

[0064] Among them, x represents the clean image, y represents the noisy image, and n multi represents multiplicative noise, n add represents additive noise.

[0065] This method evaluates the noise type of the processed image based on a composite noise model. The evaluation includes both additive and multiplicative noise, resulting in a more accurate representation of the actual noise distribution. The specific evaluation method is based on existing technologies and selects an appropriate filter method for evaluation by analyzing the noise distribution of low-dose CBCT images. In this embodiment, the noise distribution model of low-dose CBCT images is simulated by analyzing the histogram distribution of a specified region. The local region for noise evaluation can be a relatively stable area, such as soft tissue or teeth; the region size can be freely selected based on the size of the teeth.

[0066] Additive noise is directly superimposed on the original signal. Additive noise exists independently of the signal and has nothing to do with the signal; multiplicative noise is the multiplication of noise and the original signal. It changes with the signal amplitude and is directly related to the signal strength. The two types of noise are associated with the signal in different ways and act in different ways, and the denoising methods are also different. In this method, additive noise and multiplicative noise are removed in two steps: first, the additive noise is removed in the spatial domain by using the adaptive median filtering method, and then the multiplicative noise is converted into additive noise. The converted additive noise is then removed by using two improved wavelet domain filtering methods. This method combines the filtering methods of the spatial domain and the wavelet domain, and combines the adaptive denoising filtering method based on noise estimation, the edge restoration and the brightness consistency adjustment method to achieve the goal of removing noise while retaining image details. As Figure 2 As shown, the noise processing process performed in the noise processing module includes three stages: noise simplification processing, high and low frequency adaptive filtering and image inverse transformation. The noise simplification processing includes additive noise processing and multiplicative noise conversion.

[0067] S4: Based on the noise distribution model corresponding to the image to be processed, the additive noise in the image to be processed is removed by using the adaptive median filtering method.

[0068] The adaptive median filtering method can be implemented based on existing technology. Figure 3 As shown in Figure 2, the additive noise processing uses the adaptive median filtering method to process the additive noise in the speckle noise. First, initialize the sliding window size S xy and the maximum window size S max ; Then determine the median pixel value z med Is it at the minimum pixel value z min and the maximum pixel value z max If the value range is between 0 and 1, the pixel value z at the center position (x, y) will be determined. xy Is it at the minimum pixel value z min and the maximum pixel value z max In the range of values ​​between , the value is z xy Otherwise, use the intermediate value z med Replace the center pixel value z xy ; if z med If it is not within the range, increase S xy Size, S xy Greater than S max When directly using the intermediate value z med Replace the center pixel value z xy Otherwise, repeat the above process. Additive noise processing can effectively deal with granular salt and pepper noise in images.

[0069] After filtering, the CBCT image still contains multiplicative noise, which needs to be converted into additive noise so that it can be effectively removed using a filter designed specifically for additive noise. For the image to be processed after the additive noise has been removed, a logarithmic transformation is used to convert the multiplicative noise in the image into additive noise, resulting in the following image after preliminary processing. The process of converting multiplicative noise into additive noise using a logarithmic transformation is specifically expressed as:

[0070] log(n multi x)=log(n multi )+log(x);

[0071] Where n multi x represents the image to be processed after removing the additive noise. After conversion, it is represented as the clean CBCT image log(x). There is new additive noise log(n multi ).

[0072] The high- and low-frequency adaptive filtering designed by this method includes frequency separation, low-frequency adaptive Wiener filtering and high-frequency adaptive wavelet threshold denoising.

[0073] S5: After preliminary processing, the image is decomposed into high-frequency components and low-frequency components using discrete wavelet transform, which are respectively recorded as: high-frequency component image and low-frequency component image.

[0074] This method first removes additive noise to achieve preliminary noise reduction in the spatial domain, primarily to simplify the noise model. Because multiplicative noise is directly related to signal strength, it cannot be directly denoised in the spatial domain. Instead, it must be separated into high-frequency and low-frequency data, each of which undergoes targeted noise reduction. Therefore, the low-frequency and high-frequency signals are separated in the wavelet domain for signals containing multiplicative noise. Low-frequency signals primarily contain the main visual information, while high-frequency signals primarily contain edge details and concentrated noise. Separating these signals allows for both effective noise reduction and the preservation of edge details.

[0075] When performing wavelet transforms, the higher the number of stages, the greater the computational effort and complexity. The selection should be based on actual needs. In this embodiment, a first-order discrete wavelet transform is used to decompose the image into high-frequency and low-frequency components, resulting in one low-frequency component and three high-frequency components. This ensures that noise can be effectively removed with limited computational effort. The three high-frequency components represent the horizontal detail component, the vertical detail component, and the diagonal detail component, respectively.

[0076] The low-frequency component denoising in this method adopts an improved adaptive Wiener filter, uses a new noise variance estimation and weighted local mean calculation to estimate the noise level more accurately.

[0077] S6: De-noise the low-frequency component image using the low-frequency adaptive Wiener filtering method to obtain the de-noised low-frequency feature Ild , after replacing the corresponding pixel points, the denoised low-frequency component data is obtained;

[0078] Low-frequency characteristics after noise reduction I ld The calculation method is:

[0079]

[0080] Among them, I ld Represents the low-frequency characteristics after noise reduction, μ l Represents the pixel mean value of the local area weighted by w in the low-frequency component; σ l 2 represents the variance of pixels in the local area; ν lf 2 is a hyperparameter, indicating the noise variance; max() is the maximum value calculation; I n Indicates the pixel value of each pixel in the low-frequency component.

[0081] ν lf 2 The calculation method is:

[0082]

[0083] Where N represents the local area window, σ lf 2 Represents the set of all local pixel variances in the low-frequency component; mean() is the mean calculation.

[0084] In this embodiment, the local area window N is: N=3×3; the corresponding calculation method of the weight w for local mean calculation is:

[0085]

[0086] This method proposes an improved low-frequency adaptive Wiener filtering method, which divides the low-frequency component image into local area windows N. Denoise each pixel in each local area window, and calculate the pixel mean μ after weighting w in the local area of ​​the low-frequency component based on the local area window N. l , to ensure that the calculated low-frequency feature after noise reduction I ld More accurate. This method uses local statistical methods to more accurately calculate the noise level and is superior to traditional Wiener filtering methods in removing low-frequency components.

[0087] This method adopts an improved adaptive wavelet threshold denoising method for high-frequency component denoising. By combining noise estimation to generate an adaptive soft threshold to transform the wavelet threshold function, it can effectively suppress noise while retaining detail information.

[0088] S7: By improving the high-frequency adaptive wavelet threshold denoising method, the high-frequency component image is subjected to noise reduction processing, which specifically includes the following steps.

[0089] a1: Combined with noise estimation to generate adaptive soft threshold T hi ;

[0090]

[0091] Where, T hi represents the soft threshold of the i-th high-frequency component, v hi 2 represents the noise variance of the i-th high-frequency component, σ i 2 represents the variance of the i-th high-frequency component; i represents the sequence number of the high-frequency component obtained by wavelet transformation separation. In this embodiment, there are three high-frequency components, and the value of i is i=0, 1, 2.

[0092]

[0093] Where, median(|H i |) represents the middle pixel value in the i-th high-frequency component, σ hi 2 represents the set of all local area pixel variances in the i-th high-frequency component; N is the local area window, and in this embodiment, N=5×5; M represents the size of the i-th high-frequency component.

[0094] a2: Transform the wavelet threshold function through soft thresholding to obtain the wavelet coefficient w after soft thresholding λ,hi ;

[0095]

[0096] Where sign() is the sign function, w hi represents the wavelet coefficient of the i-th high-frequency component; w λ,hi Represents the wavelet coefficients after threshold function processing.

[0097] a3: Wavelet coefficient w after soft threshold processing for each high-frequency component λ,hi Perform signal reconstruction to generate denoised high-frequency component data.

[0098] In multiplicative noise, the noise variance of high-frequency components will affect the distribution of noise coefficients. In this method, the soft threshold of high-frequency components is calculated by the noise variance of high-frequency components, and then the wavelet coefficients w are transformed by the soft threshold. λ,hi, ensuring that the calculation results are more targeted. The improved high-frequency adaptive wavelet threshold denoising method proposed in this method calculates the thresholds of different high-frequency components through noise estimation to achieve the function of adaptive denoising. Compared with the traditional wavelet threshold denoising method, it can effectively suppress noise while preserving image details.

[0099] S8: Performing an inverse transform on the denoised low-frequency component data and the denoised high-frequency component data to restore the previous data structure and obtain the denoised image to be processed. Specifically, for the previous signal separation process and noise conversion process, the inverse transforms implemented include: sequentially performing an inverse wavelet transform and an inverse logarithmic transform.

[0100] During the denoising process, the edge information of the image may be damaged, so it is necessary to perform an edge restoration operation on the image to be processed to restore the structure.

[0101] S9: performing edge restoration processing on the denoised image to be processed to obtain an edge-restored image to be processed.

[0102] Edge restoration involves sequential image smoothing and sharpening. Smoothing uses the first step of the BM3D (Block-Matching and 3D Filtering) algorithm for image noise reduction and smoothing, using CUDA for parallel acceleration. Image sharpening uses the Unsaturated Mask (USM) algorithm to enhance edge information smoothed during the noise reduction process and improve edge detail. This process includes the following steps.

[0103] The specific steps of the BM3D algorithm in the first step of edge recovery processing are as follows Figure 4 As shown, the specific steps include:

[0104] b1: For the denoised image to be processed, select a batch of k×k reference blocks; search in the m×m area around the reference blocks to find several blocks with the smallest difference, denoted as arbitrary blocks, and integrate these arbitrary blocks and the reference blocks into a three-dimensional matrix G(P);

[0105]

[0106] Where P is the reference block, Q is any block found in the search area, d(P,Q) represents the Euclidean distance between the two blocks, τ represents the difference threshold, and G(P) represents the three-dimensional matrix obtained by integrating similar blocks.

[0107] In this embodiment, in order to better utilize CUDA computing resources, the CUDA thread warp method is used to improve the operating efficiency of similar block grouping. First, k is fixed to 8. Considering the operating speed and noise level, m can be selected between 8 and 32, and the number of blocks with the smallest difference can be selected between 8 and 16. At the same time, the search step size is increased, and the optional step size is between 3 and 6.

[0108] b2: Next, collaborative filtering is performed, which can be performed based on several three-dimensional matrices.

[0109] A two-dimensional discrete cosine transform is performed on each two-dimensional block in the three-dimensional matrix G(P); in this embodiment, a two-dimensional discrete cosine transform (DCT transform) is used; a one-dimensional Hadamard transform is then performed on the third dimension of the three-dimensional matrix; in this embodiment, a one-dimensional Hadamard transform (WHT transform) is used; a hard threshold process is performed on the transformed three-dimensional matrix, and coefficients less than the threshold are set to 0; and then a one-dimensional inverse transform and a two-dimensional inverse transform are performed on the third dimension to obtain the processed image block.

[0110] The process of step b2 is expressed as:

[0111]

[0112] Among them, T 3Dhard It is uniformly expressed as a two-dimensional transformation and a one-dimensional transformation; Q(P) is the three-dimensional matrix after collaborative filtering processing;

[0113] γ represents the hard thresholding operation, which is as follows:

[0114]

[0115] Among them, σ represents the standard deviation of noise, which represents the noise intensity; λ 3D Represents a hyperparameter, which can be selected as 2.7;

[0116] b3: The similar block aggregation operation fuses the two-dimensional denoised image blocks processed by collaborative filtering to their original positions. The grayscale value of each pixel can be obtained by weighted averaging the blocks at each corresponding position. The aggregation weight used in the weighted averaging calculation can be calculated by reciprocally calculating the number of non-zero values ​​obtained in the hard thresholding operation to generate the weight coefficient of each similar block.

[0117] Image sharpening uses the Unsharpen Mask (USM) algorithm to enhance edge information that may have been smoothed during the previous noise reduction process. The specific process can be implemented based on existing technology. The calculation formula of this algorithm is as follows:

[0118]

[0119] Where I represents the input source image, which corresponds to the image after noise reduction and image smoothing output in step b3. G Represents the image after Gaussian filtering, I s Represents the image after USM sharpening; w represents the sharpening factor, which ranges from 0.1 to 0.9.

[0120] Since adaptive noise reduction methods applied to a single 2D CBCT image can lead to brightness inconsistencies between layers, this method corrects these inconsistencies by adjusting the brightness of each layer. The 3D data equalization phase of processing the image after edge restoration involves brightness consistency adjustment to equalize image brightness between layers. Specifically, the following steps are included.

[0121] S10: After all 2D edges are restored, the images to be processed are stitched into 3D data again in the order of the previous slices to obtain a stitched CBCT 3D image; and 3D data equalization processing is performed on each layer of the 2D image of the stitched CBCT 3D image to obtain a final quality-enhanced image.

[0122] The three-dimensional data equalization processing includes calculations:

[0123] I c =min(max(r·I′ filter ,min(I′ noise )),max(I′ noise ))+min(I noise );

[0124] Where, I c Indicates the image after brightness adjustment; I filter is the image after noise reduction filtering, corresponding to the image to be processed after edge restoration output in step S9; noise is the source noise image, corresponding to the 2D image extracted based on the CBCT 3D image slices in step S2;

[0125] I' filter and I' noise is the median value and is calculated as follows:

[0126] I′ filter =I filter -min(I noise );

[0127] I′ noise =I noise -min(I noise );

[0128] r represents the weight coefficient of each layer of image, and the calculation method is as follows:

[0129]

[0130] After using the technical solution of the present invention, for three-dimensional low-dose CBCT images, slices are first generated into a single 2D image, and then a composite noise model is used to accurately estimate the noise of the generated single 2D image to ensure the accuracy of the subsequent noise reduction process; in the noise reduction stage, the additive noise is removed in the spatial domain, and after the multiplicative noise is converted into additive noise, it is separated into high-frequency components and low-frequency components, and targeted filtering methods are designed for accurate noise reduction to ensure that both the noise reduction effect and the requirement of maintaining edge details are taken into account; after the original data structure of the denoised image is restored, the edge pixels lost in the noise reduction are restored, the edge-restored data is reconstructed into a three-dimensional image, and the three-dimensional image data is adjusted for consistency based on the brightness between layers to ensure that a high-quality quality-enhanced image is finally obtained.

Claims

1. A method for enhancing low-dose CBCT image quality, characterized in that: It includes the following steps: S1: Acquire CBCT noise images; Low-dose head model projection data were collected, and the corresponding noisy CBCT three-dimensional images were reconstructed using the FDK algorithm. S2: extracting a 2D axial slice based on the CBCT 3D image, which is recorded as the image to be processed; S3: For each image to be processed, analyze the noise distribution type and construct a composite noise distribution model: y=n multi *x+n add 4 Among them, x represents the clean image, y represents the noisy image, and n multi represents multiplicative noise, n add represents additive noise; S4: Based on the noise distribution model corresponding to the image to be processed, after removing the additive noise in the image to be processed, the multiplicative noise in the image is converted into additive noise to obtain: a preliminarily processed image; S5: Decomposing the preliminarily processed image into a high-frequency component and a low-frequency component using discrete wavelet transform, which are respectively recorded as: a high-frequency component image and a low-frequency component image; S6: De-noising the low-frequency component image by using a low-frequency adaptive Wiener filtering method to obtain the de-noised low-frequency feature I ld , after replacing the corresponding pixel points, the denoised low-frequency component data is obtained; Low-frequency characteristics after noise reduction I ld The calculation method is: Among them, I ld Represents the low-frequency characteristics after noise reduction, μ l Represents the pixel mean value of the local area weighted by w in the low-frequency component; σ l 2 represents the variance of pixels in the local area; ν lf 2 is a hyperparameter, indicating the noise variance; max() is the maximum value calculation; I n Indicates the pixel value of each pixel in the traversal low-frequency component; ν lf 2 The calculation method is: Where N represents the local area window, σ lf 2 Represents the set of all local pixel variances in the low-frequency component; mean() is the mean calculation; S7: performing noise reduction processing on the high-frequency component image by improving the high-frequency adaptive wavelet threshold denoising method, specifically comprising the following steps: a1: Combined with noise estimation to generate adaptive soft threshold T hi ; Where, T hi represents the soft threshold of the i-th high-frequency component, v hi 2 Represents the noise variance of the i-th high-frequency component, i represents the sequence number of the high-frequency component separated by wavelet transform, σ i 2 represents the variance of the i-th high-frequency component; Where, median(|H i |) represents the middle pixel value in the i-th high-frequency component, σ hi 2 Represents the set of all local pixel variances in the i-th high-frequency component, N is the local area window, and M represents the size of the i-th high-frequency component; a2: Transform the wavelet threshold function through soft thresholding to obtain the wavelet coefficient w after soft thresholding λ,hi ; Where sign() is the sign function, w hi represents the wavelet coefficient of the i-th high-frequency component; w λ,hi Represents the wavelet coefficients after threshold function processing; a3: wavelet coefficient w after soft threshold processing λ,hi Perform signal reconstruction to generate denoised high-frequency component data; S8: performing an inverse transformation on the denoised low-frequency component data and the denoised high-frequency component data to restore the previous data structure, thereby obtaining a denoised image to be processed; S9: performing edge restoration processing on the denoised image to be processed to obtain an edge-restored image to be processed; The edge restoration process includes: image smoothing and image sharpening performed in sequence; S10: stitching all the 2D edge-restored images to be processed into 3D data again according to the order of the previous slices to obtain a stitched CBCT 3D image; performing 3D data equalization processing on each 2D image layer of the stitched CBCT 3D image to obtain a final quality-enhanced image; The three-dimensional data balancing process includes: brightness consistency adjustment and balancing the image brightness between layers.

2. The method for enhancing low-dose CBCT image quality according to claim 1, characterized in that: In step S4, an adaptive median filtering method is used to remove additive noise in the image to be processed.

3. The method for enhancing low-dose CBCT image quality according to claim 1, characterized in that: In step S4, the multiplicative noise is converted into additive noise using a logarithmic transformation.

4. The method for enhancing low-dose CBCT image quality according to claim 1, characterized in that: In step S5, the discrete wavelet transform is used to decompose the signal into high-frequency components and low-frequency components, thereby obtaining one low-frequency component and three high-frequency components.

5. The method for enhancing low-dose CBCT image quality according to claim 1, characterized in that: In step S6, the local area window N is: N=3×3; the corresponding local mean calculation weight w is calculated as follows:

6. The method for enhancing low-dose CBCT image quality according to claim 1, characterized in that: In step S8, the inverse transformation performed on the denoised low-frequency component data and the denoised high-frequency component data includes: sequentially performing an inverse wavelet transformation and an inverse logarithmic transformation.

7. The method for enhancing low-dose CBCT image quality according to claim 1, characterized in that: In the edge restoration process, the image smoothing process uses the first step basic estimation of the BM3D algorithm to perform image noise reduction, and uses CUDA for parallel acceleration; the image sharpening uses the USM algorithm to enhance the edge information that was smoothed in the previous noise reduction process.

8. The method for enhancing low-dose CBCT image quality according to claim 1, characterized in that: The three-dimensional data equalization processing specifically includes the following calculations: IN c =min(max(r·I′ filter ,mini' noise )),maxi' noise ))+min(I noise ); Where, I c Indicates the image after brightness adjustment, I filter is the image after noise reduction filtering, I noise is the source noise image; I' filter and I' noise is the median value and is calculated as follows: IN' filter =I filter -mini noise ); IN' noise =I noise -mini noise ); r represents the weight coefficient of each layer of image, and the calculation method is as follows:

Images