Oral cavity CBCT multi-material beam hardening correction method

Abstract

The invention provides an oral cavity CBCT (cone beam computed tomography) multi-material beam hardening correction method, which is particularly suitable for a multi-material scene and can be used for removing dark bands and stripe artifacts around teeth, bones and metal implants to obtain a reconstructed high-precision image. The method comprises the following steps: firstly, removing scattering artifacts in CBCT imaging data to obtain projection data after scattering correction, then determining all materials included in the data, areas where the materials are located and segmentation volumes, and simulating to obtain multispectral simulation multicolor projection and monochromatic projection simulation data; and calculating a difference value between the fitting multicolor projection simulation data and the optimized reference projection data to obtain a non-linear correction item between a CT measurement value and a true value, and then superposing the non-linear correction item to the projection data after scatter correction to recover the FDK measurement value to the true value so as to complete beam hardening correction of the original CBCT imaging data.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of image recognition, and particularly relates to a multi-material beam hardening correction method for oral CBCT. BACKGROUND

[0002] The X-ray in the oral CBCT system is usually a multi-color light mixed energy beam. Beam hardening (BH) is a core physical defect of the oral CBCT (cone beam CT) imaging system, and its essence is derived from the multi-color energy spectrum characteristics of the X-ray source and the nonlinear relationship of the material attenuation coefficient. According to the Lambert-Beer law, when the multi-energy X-ray penetrates the heterogeneous object, the low-energy photons (<50keV) are preferentially absorbed by the high atomic number material, resulting in the phenomenon of "hardening" of the penetrating energy spectrum--the average energy of the energy spectrum is shifted from the initial average energy to the high energy region (the typical shift amplitude can reach 15-30keV). This energy spectrum evolution process makes the linear attenuation coefficient present significant energy dependence, which directly violates the linear assumption of the single-color projection in the FDK (Feldkamp-Davis-Kress) reconstruction algorithm, and then induces two characteristic artifacts:

[0003] (1) Cupping Artifact:

[0004] In the homogeneous density material, the hardening effect causes the systematic underestimation of the projection data in the central region, and the reconstructed image shows the false low density (CT value deviation can reach 200-400HU) in the central region, which seriously affects the quantitative analysis of the bone trabecular structure.

[0005] (2) Streaking Artifact:

[0006] In the high-density material interface (such as the enamel-dentin junction, titanium implant-bone tissue interface), the hardening effect and the photon starvation effect are coupled to produce radial dark band artifacts. Clinical data shows that this kind of artifact seriously reduces the detection rate of periapical lesions, and in the multi-material coexistence area (such as the orthodontic bracket-tooth-bone complex), the artifact superposition effect will cause a diagnostic blind area.

[0007] To improve the imaging clarity, the technical staff has researched many beam hardening correction techniques, such as: hardware filtering method, by adding K-edge filter (commonly used 0.5-3mm aluminum / copper composite filter) at the outlet end of the X-ray tube, 10-40keV low-energy photons can be cut off; but experiments show that even if a 2mm copper filter is used, only the stripe artifacts of the tooth-bone interface can be slightly improved, and the loss of photon flux leads to a decrease in signal-to-noise ratio (SNR), forcing the dose to be increased to 1.8-2.5 times of the conventional scan to maintain the image quality. Polynomial fitting, the second-order polynomial correction based on the water phantom calibration can significantly reduce the cup-shaped artifacts of the uniform phantom, but there is a serious mismatch in the multi-material scene; when applied to the tooth-bearing body, a 6-8mm over-correction ring (CT value abnormally increased by 300-500HU) appears on the edge of the dentin, which is caused by the fact that the polynomial model cannot analyze the transition characteristics of mu(E) with atomic number Z. Iterative reconstruction (statistical model) method, the iterative algorithm based on maximum likelihood expectation maximization (MLEM) needs 200-300 iterations to converge, and the reconstruction time of a single scan is as long as 45-60 minutes (still 8-12 minutes after GPU acceleration); more seriously, the existing medical iterative algorithm relies on the preset base material decomposition model (usually water-bone-iodine three groups), and the dental alloy (such as CoCr Z=27) existing in the oral environment is out of the model coverage range, resulting in a large number of hardening artifacts remaining. Post-correction method, this technology uses a closed-loop correction framework of threshold segmentation-projection simulation, but the metal artifacts in the oral CBCT itself will distort the segmentation boundary (the edge positioning error is 1.2-2.3mm); Monte Carlo simulation shows that when the segmentation error exceeds 0.5mm, the inaccuracy of scattered photon correction will introduce new ring-shaped artifacts, forming a positive feedback of error. Dual-energy method, although dual-energy CT can theoretically eliminate all orders of BH effects by solving the high and low energy spectra (80kVp / 140kVp) simultaneously, there are multiple obstacles in clinical implementation: spatial misalignment caused by patient movement between two scans often occurs in child patients; the heat load of the ball tube exceeds the heat dissipation limit of the oral CBCT, which will trigger the overheat protection; far beyond the capacity of basic clinics. In summary, the existing beam hardening correction techniques have deficiencies in application in multi-material scenes, resulting in inaccurate CBCT imaging in some areas, affecting the accuracy of subsequent reading and checking. SUMMARY

[0008] In order to solve the problem of the existing beam hardening correction techniques in the application in multi-material scenes, the present application provides a kind of oral CBCT multi-material beam hardening correction method, which is especially suitable for multi-material scenes, can remove the dark band, stripe artifact around teeth, bones and metal implants, and obtain high-precision images after reconstruction.

[0009] The technical solution of the present application is as follows: a kind of oral CBCT multi-material beam hardening correction method, characterized in that it comprises the following steps:

[0010] S1: scatter artifact preprocessing on CBCT raw scan data;

[0011] Specifically comprising the following operations:

[0012] Obtain the two-dimensional projection data of the CBCT raw three-dimensional image, denoted as original two-dimensional projection x u ;

[0013] Use multi-GPU accelerated Monte Carlo (MC) photon transport model (FPM) to simulate the spatial distribution of scattered photons in the original image, and obtain the scatter component;

[0014] Subtract the scatter component from the original two-dimensional projection data by iterative subtraction, and output the scatter-corrected pure projection data, denoted as: scatter-corrected projection data x sc ;

[0015] S2: blind estimation of object material;

[0016] Specifically comprising the following operations:

[0017] Filter back projection (FDK) reconstruction is performed on the scatter-corrected projection data to obtain an initial three-dimensional volume; adaptive gray threshold segmentation method is used to separate different material regions in the volume;

[0018] Call the pre-constructed oral material attenuation coefficient lookup table, and iteratively match the material types of each region by the minimum mean square error (MMSE) criterion, label the material types of each region, and output the material-labeled segmented volume;

[0019] S3: multi-spectral simulation of multi-color projection;

[0020] Specifically comprising the following operations:

[0021] Based on the segmented volume corresponding to each material region, use multi-GPU parallel architecture to synchronously simulate multi-color projection data;

[0022] Use least square estimation (LSE) algorithm to combine all the multi-color projections into optimal fitting projections, denoted as: multi-color projection simulation data;

[0023] S4: single-color projection generation and optimization;

[0024] Specifically comprising the following operations:

[0025] Based on the segmented volume corresponding to each material region, simulate multiple single-color projections of different energy levels, denoted as: single-color projection simulation data;

[0026] The mean square error (MSE) of the monochromatic projection simulation data and the scatter-corrected projection data is calculated, and the monochromatic projection data with the highest similarity to the original data is selected as the optimized reference projection data from all the monochromatic projection simulation data.

[0027] S5: original data beam hardening correction;

[0028] Specifically comprising the following steps:

[0029] The difference value of fitting the multicolor projection simulation data and the optimized reference projection data is calculated to obtain a non-linear correction term;

[0030] After superimposing the non-linear correction term on the scatter-corrected projection data, a three-dimensional image is reconstructed by FDK to obtain a corrected CBCT image.

[0031] It is further characterized in that:

[0032] In step S1, the following steps are included in detail:

[0033] S11: obtaining original two-dimensional projection x u of oral CBCT scanning;

[0034] S12: modeling the scatter photon distribution;

[0035] By loading a GPU-accelerated Monte Carlo photon transport model (FPM), the original two-dimensional data x u is back-projected and reconstructed to generate an initial attenuation coefficient volume, simulating the spatial distribution of scatter photons in the original image; to accelerate the scatter photon simulation process, GPU parallelization is used to accelerate the output of the scatter photon distribution; the initial value estimation of the simulated scatter projection distribution x scatter is:

[0036]

[0037] Wherein, x sc (0) is the initial scatter correction projection, α is the initial step size, and x u is the original two-dimensional projection data of the CBCT scanning image;

[0038] S13: subtracting the scatter component from the original two-dimensional projection data by iterative subtraction to output the scatter-corrected pure projection data, denoted as: scatter-corrected projection data x sc ;

[0039] The iterative correction formula is:

[0040]

[0041] Wherein, k is the iteration number; βk is adaptive step, decreasing each iteration; FPM scatter is a GPU-MC scatter estimation subroutine;

[0042] In step S2, the following steps are included in detail:

[0043] S21: scatter-corrected projection data x sc is filtered back projection FDK reconstruction, to obtain the initial three-dimensional volume V sc ;

[0044] S22: using adaptive gray threshold segmentation method, the different material regions in the initial three-dimensional volume V sc are segmented, and the specific method is:

[0045] By using the multi-threshold segmentation method, the threshold T1 and T2 are obtained by Ostu multi-threshold method to automatically segment the reconstruction data into three regions;

[0046]

[0047] In the formula, V soft represents the soft tissue region, V bone represents the bone and tooth region, and V metal is the metal implant region; a represents the initial three-dimensional volume attenuation coefficient value;

[0048] S23: morphological operation closing operation is performed on each segmented region to eliminate segmentation noise;

[0049] S24: reading the pre-constructed oral material attenuation coefficient lookup table;

[0050] The oral material attenuation coefficient lookup table records the density and attenuation coefficient of different types of materials corresponding to different energies;

[0051] S25: based on the oral material attenuation coefficient lookup table, the corresponding material type of each segmented region V i is matched by the minimum mean square error MMSE criterion; it specifically includes the following steps:

[0052] S251: for each segmented region V i , the average projection value is calculated;

[0053]

[0054] In the formula, represents each region V i, mean() is the mean function; FPLS is the fast linear integration model, which directly projects the segmented single-material volume along the ray's incident direction to calculate the sum of the attenuation values of the ray passing through voxels;

[0055] S252: After thickness compensation, the average projection value of the segmented volume V is searched in the oral material attenuation coefficient lookup table to find the minimum Euclidean distance of the closest material, and the material type of each segmented volume V i is determined by the minimum mean square error (MMSE) criterion.

[0056]

[0057] where j * represents the closest material, and the argmin function returns the argument value when the given objective function takes the minimum value.

[0058] The MMSE criterion is based on the L2 norm, denoted as ||.||2.

[0059] j represents different materials in the oral material attenuation coefficient lookup table, μ j (E ref ) represents the attenuation coefficient of material j at the reference energy E ref ; t i is the average penetration thickness of the segmented volume V i ; E ref is the reference energy, which has the same value as the scanning energy corresponding to the original scanning data.

[0060] S26: After determining the corresponding material of each segmented volume V i , the material type of each segmented volume is labeled, and the segmented volume V seg of each region after material labeling is output.

[0061] Before step S26 is implemented, the corresponding material type of each segmented volume V i needs to be corrected and confirmed, which includes the following steps.

[0062] a1: In the oral material attenuation coefficient lookup table, the density ρ * and linear attenuation coefficient μ est of the j est corresponding material are confirmed, and the density-based Monte Carlo projection x FPM is calculated.

[0063] x FPM = FPM(ρ est );

[0064] where ρest FPM is a Monte Carlo photon transport model, and

[0065] a2: calculate linear projection x based on attenuation coefficient FPLS :

[0066] x FPLS = FPLS(μ est );

[0067] where FPLS is a fast linear integral model, and μ est is a linear attenuation coefficient;

[0068] a3: construct a target optimization function:

[0069]

[0070] a4: minimize the residual of the target function, and output the optimized material label volume V seg ;

[0071] In step S3, the following steps are included in detail:

[0072] S31: based on the segmentation volume V seg corresponding to each material region, use the Monte Carlo photon transport model FPM to simulate the projection of polychromatic X-rays passing through the object, generate 4 sets of polychromatic projections under GPU acceleration, and simulate different filter combinations by SpekCalc software; wherein the 4 energy spectrums used in the polychromatic projection are centered on the energy corresponding to the CBCT original scanning data, and the values are taken in the energy range of 100-200keV;

[0073] S32: perform least squares estimation LSE fitting, first establish a projection fitting model;

[0074]

[0075] where x p,i is the i-th set of simulated projection, i is the serial number of the polychromatic projection, 1≤i≤4; c i is the weight coefficient of the i-th set of projection;

[0076] S33: to solve the optimal weight c i , by constructing the L2 norm of the minimum real projection x μ and the simulated projection x p :

[0077]

[0078] Solve by matrix:

[0079]

[0080] wherein, is the simulated projection autocorrelation matrix, is the simulated projection cross-correlation vector with the real projection;

[0081] S34: After obtaining the optimal weight through the matrix inversion fast calculation, the multi-color projection simulation data x is obtained p ;

[0082] In step S4, the following steps are specifically included:

[0083] S41: Energy discretization and initial simulation;

[0084] Discretize the X-ray energy spectrum into K energy intervals; perform density conversion according to the closest material j * corresponding linear attenuation coefficient μ, and convert it into mass density ρ;

[0085] Use the GPU accelerated Monte Carlo model FPM to simulate the projection of monochromatic X-rays through the object, use the adaptive GPU Monte Carlo model FPM input requirements to simulate monochromatic projection; for each energy bin E k , use FPM to generate monochromatic projection x m (E k );

[0086] The calculation formula is as follows:

[0087]

[0088] wherein, I mono,0 (E k ) = Nη(E k ), represents the intensity before penetrating the object, N is the number of incident photons, and η(E k ) represents the energy response of the detector;

[0089] is the intensity after penetrating the object, wherein L is the path length of the ray through the material;

[0090] S42: For each simulated monochromatic projection x m (E k ), calculate the mean square error MSE with the real projection x u ;

[0091]

[0092] wherein, N is the total number of projection pixels;

[0093] S43: Select the energy E opt with the minimum MSE as the monochromatic projection reference, denoted as: optimal reference projection data;

[0094]

[0095] where argmin function returns the argument value at which the given objective function takes the minimum value;

[0096] In step S5, the following is included:

[0097] S51: record the optimization reference projection data as: x m (E opt ), and the polychromatic projection simulation data x p Calculate the nonlinear correction term;

[0098] Φ CT = x m (E opt )-x p ;

[0099] where Ф CT represents the compensation of the attenuation nonlinearity caused by polychromatic spectrum;

[0100] S52: Perform projection data beam hardening correction: superimpose Ф CT to the scatter corrected data x sc , to obtain the final corrected projection data x c :

[0101] x c = Φ CT +x sc ;

[0102] S53: Use the filtered back projection FDK reconstruction algorithm to reconstruct the corrected projection data x c into a three-dimensional image V c , to the corrected CBCT image.

[0103] This application provides a multi-material beam hardening correction method for oral CBCT. The fundamental reason for the formation of beam hardening artifacts is that the multicolor energy spectrum characteristics of the X-ray source and the material attenuation coefficient are nonlinear. When the FDK method reconstructs CBCT images, it uses approximately one million independent detector measurements and treats the data of each detector as monochromatic light based on a linear assumption for reconstruction. This leads to an error between the CT measurement values ​​in the reconstructed image and the true attenuation coefficient of the scanned object. These measurement errors are reflected in the reconstructed image as artifacts. Therefore, this method first removes scattering artifacts from CBCT imaging data to obtain scatter-corrected projection data. Then, using the scatter-corrected projection data, a pre-constructed dedicated attenuation coefficient lookup table is used to determine all materials included in the data, their regions, and segmentation volumes. Based on parameters such as material type, attenuation coefficient, and segmentation volume, multispectral simulated multicolor projection and monochromatic projection simulation data are obtained. Since multispectral simulated multicolor projection and the original CBCT imaging data have similar nonlinear characteristics, monochromatic projection data is used as the optimized reference projection data. The difference between the fitted multicolor projection simulation data and the optimized reference projection data is calculated, which yields the error value between the CT measurement value and the true value. These error values ​​are used as nonlinear correction terms for the original CBCT imaging data. After superimposing the nonlinear correction terms onto the scatter-corrected projection data, the FDK measurement value can be restored to the true value. This can simultaneously suppress cup-shaped artifacts around the metal crown and bright streaks in bone tissue caused by secondary scattering hardening, thus completing the beam hardening correction of the original CBCT imaging data. Attached Figure Description

[0104] Figure 1 This is the overall flowchart of this method;

[0105] Figure 2 A flowchart for blind estimation of object materials;

[0106] Figure 3 Example of a lookup table structure for attenuation coefficients of dental materials. Detailed Implementation

[0107] like Figure 1 As shown, this invention includes a multi-material beam hardening correction method for oral CBCT, used to remove dark bands and streaks around teeth, bones, and metal implants, providing doctors with high-precision images for diagnosis. It includes the following steps.

[0108] S1: Preprocess the raw CBCT scan data to remove scattering artifacts. This includes the following steps:

[0109] The two-dimensional projection data of the original three-dimensional image of the oral cavity CBCT is obtained and denoted as the original two-dimensional projection.

[0110] The spatial distribution of scattered photons in the original image is simulated by using a multi-GPU accelerated Monte Carlo photon transport model FPM (Forward Projection Model), to obtain a scattered component;

[0111] The scattered component is subtracted from the original two-dimensional projection data by iterative subtraction, and the scatter-corrected pure projection data is output, denoted as: scatter-corrected projection data. Specifically, a fast iterative scatter correction algorithm is used to ensure efficient removal of scatter artifacts.

[0112] In step S1, the following steps are included in detail:

[0113] S11: Obtain the original two-dimensional projection x of the oral CBCT scan u ;

[0114] S12: Model the distribution of scattered photons;

[0115] The GPU-accelerated Monte Carlo photon transport model FPM is loaded, and the model physical configuration parameters are: 120KV X-ray source tungsten target spectrum (generated by SpekCalc software), 2.5mm thick aluminum filter, and detector non-uniform response gain correction table.

[0116] According to the original two-dimensional data x u , the initial attenuation coefficient volume is generated by back-projection reconstruction, simulating the spatial distribution of scattered photons in the original image; to accelerate the scattered photon simulation process, GPU parallelization is used to accelerate the output of the scattered photon distribution; the simulated scattered projection distribution x scatter is obtained.

[0117]

[0118] Wherein, x sc (0) is the initial scatter correction projection, α is the initial step size, x u is the original two-dimensional projection data of the CBCT scan image;

[0119] S13: The scattered component is subtracted from the original two-dimensional projection data by iterative subtraction, and the scatter-corrected pure projection data is output, denoted as: scatter-corrected projection data x sc ;

[0120] The iterative correction formula is:

[0121]

[0122] Wherein, k is the iteration number, the iteration number in this embodiment is three, and k takes values of 0, 1 and 2; β kFor adaptive step size, each iteration is reduced; FPM scatter For GPU-MC scatter estimation subroutines.

[0123] S2: blind estimation of object material; as Figure 2 shown, specifically includes the following operations.

[0124] The filtered back-projection FDK reconstruction is performed on the scatter-corrected projection data to obtain an initial three-dimensional volume; an adaptive gray threshold segmentation method is used to separate different material regions in the volume; the adaptive gray threshold segmentation method can be based on various algorithms in the prior art (based on neighborhood mean calculation, adaptive threshold method, local dynamic adjustment), and in the embodiment, the Otsu algorithm is used; the oral cavity is a material region included in the scan image, such as metal implants, dentin, bone tissue, etc.

[0125] A pre-constructed oral material attenuation coefficient lookup table is called, each region material type is labeled by iteratively matching the material type of each region through the minimum mean square error MMSE criterion, projection residual optimization is performed, and a segmented volume after material labeling is output.

[0126] The method pre-constructs an oral material attenuation coefficient lookup table in the system, as Figure 3 shown, the oral material attenuation coefficient lookup table records the density and attenuation coefficient of different types of materials under different energies, and the data is obtained according to actual measurement; the "energy" column in the table corresponds to the energy value that can be used in the scan, and the materials recorded in the table include not only the oral cavity inherent materials such as teeth and bones, but also dental special materials such as titanium alloy and hydroxyapatite; the oral material attenuation coefficient lookup table is stored in the system and is maintained by a dedicated person, and once a new material is applied in the oral cavity, the related data is supplemented and updated to the table to ensure that basic data support can be provided for calculation. Note Figure 3 that the data in the table are all pseudo-data, which are only used to show the structure of the table. The method iteratively matches the material type of each region through the minimum mean square error (MMSE) criterion; and outputs a segmented volume after material labeling. The blind material estimation is quickly and accurately realized through the lookup table and the minimization of the projection residual.

[0127] The method breaks through the limitation that the traditional correction method must pre-set material parameters through a dynamic expansion type oral material blind estimation mechanism, pre-constructs an oral special attenuation coefficient lookup table, and iteratively matches the material type in combination with the minimum projection residual criterion. When the material estimation error is ≤20%, the correction effect can still be maintained, the problem that the traditional method fails for unknown materials is solved, and adaptive correction without material prior knowledge is realized.

[0128] In step S2, the following steps are included in detail:

[0129] S21: projection data x after scatter correction sc FDK reconstruction is performed to obtain an initial three-dimensional volume V sc ;

[0130] S22: using an adaptive gray threshold segmentation method, the initial three-dimensional volume V sc is segmented into different material regions, and the specific method is as follows:

[0131] By using a multi-threshold segmentation method, the threshold T1 and T2 are obtained by Ostu multi-threshold method to automatically segment the reconstruction data into three regions;

[0132]

[0133] In the formula, V soft represents a soft tissue region, V bone represents a bone and tooth region, and V metal is a metal implant region; a represents the attenuation coefficient value in the initial three-dimensional volume.

[0134] S23: morphological operation closing operation is performed on each segmented region to eliminate segmentation noise;

[0135] S24: reading a pre-constructed oral material attenuation coefficient lookup table;

[0136] S25: based on the oral material attenuation coefficient lookup table, each segmented region V i is iteratively matched with the corresponding material type based on the minimum mean square error (MMSE) criterion; the specific steps include the following:

[0137] S251: for each segmented region V i , the average projection value is calculated;

[0138]

[0139] In the formula, represents the average projection value of each region V i , and mean() is the average value function; FPLS is a fast linear integration model.

[0140] Wherein, FPLS (Forward Projection based on Line-Simple method) is a forward projection model based on line-simple method, and the specific projection method is as follows: for the segmented single material volume, the forward projection is directly performed along the incident direction of the ray to calculate the sum of the attenuation values of the rays passing through the voxels.

[0141] S252: after thickness compensation, the average projection value The minimum Euclidean distance to the closest material, by the Minimum Mean Square Error, MMSE, criterion, for each segmented region V i The closest material is found by the following method:

[0142]

[0143] Where j * represents the closest material, and the argmin function returns the argument value at which the given function takes its minimum value;

[0144] The MMSE criterion is based on the L2 norm, denoted as ||.||2;

[0145] j represents different materials in the lookup table of oral material attenuation coefficients, μ j (E ref ) represents the attenuation coefficient of material j at the reference energy E ref ; t i is the average penetration thickness of the segmented region V i ; and E ref is the reference energy, which has the same value as the scanning energy corresponding to the original scan data.

[0146] When the matching error is less than or equal to 10%, it is confirmed that the corresponding material is found.

[0147] Before step S26 is implemented, it is also necessary to correct and confirm the corresponding material type for each segmented region V i . Through projection residual optimization, it is estimated whether the material is matched according to the double-model projection verification.

[0148] Specifically, the following steps are included:

[0149] a1: In the lookup table of oral material attenuation coefficients, it is confirmed that the material density ρ * and the linear attenuation coefficient μ est corresponding to j est are calculated, and the density-based Monte Carlo projection x FPM is calculated;

[0150] x FPM = FPM(ρ est );

[0151] Where ρ est is the material density, and FPM is the Monte Carlo photon transport model;

[0152] a2: The linear projection based on the attenuation coefficient x FPLS is calculated:

[0153] x FPLS = FPLS(μ est );

[0154] where FPLS is a fast linear integration model, μ est is a linear attenuation coefficient;

[0155] a3: build the target optimization function:

[0156]

[0157] a4: minimize the target function residual, output the optimized material label volume V seg .

[0158] Blind estimation of materials is achieved through a dynamic lookup table + thickness compensation matching, and various material implants such as titanium / zirconia are identified under the condition of no prior knowledge.

[0159] where the FPM model is based on the initial guessed material density ρ est (through the oral material attenuation coefficient lookup table, which can be obtained), simulates the primary photon projection x FPM (only contains primary rays, does not contain scattering); the FPLS model: directly performs forward projection on the segmented single material volume, calculates the sum of attenuation values of the rays passing through the voxels x FPLS ; iteratively adjust the material estimation, build the target function x FPLS The projection does not depend on the attenuation table and the spectrum, and is regarded as a reference standard for correcting the deviation of the simulation. The FPM model depends on the material attenuation table and the spectrum model, but these data may have errors (such as inaccurate attenuation table, spectrum simulation deviation). Iterative optimization takes the FPLS model as a reference (does not depend on the attenuation table and the spectrum), forces the FPM projection to approach the real scan data, and compensates for the model defects. FPM

[0160] S26: determine the material of each segmented region V i After determining the corresponding material, the material type of each segmented region is labeled, and the segmented volume V seg of each region after material labeling is output.

[0161] S3: multispectral simulation of polychromatic projection; specifically including the following operations.

[0162] Based on the segmented volume corresponding to each material region, the multispectral simulation of polychromatic projection data is simulated synchronously using a multi-GPU parallel architecture, and a least squares estimation (LSE) algorithm is used for fitting. All polychromatic projections are combined into an optimal fitting projection after being weighted, denoted as: polychromatic projection simulation data. The LSE fitting method is used to significantly reduce the simulation projection residual, significantly overcome the system mismatch problem caused by inaccurate energy spectrum modeling, and ensure high matching between the fitted projection and the measured data.

[0163] ​In a specific implementation, the spectrum configuration is parallel to the GPU task, 4 sets of polychromatic projections of different energy spectra are generated under GPU acceleration using the Monte Carlo photon transport model (FPM), and different filter combinations (e.g., 1.5 mm Cu, 2.5 mm Al, etc.) are simulated by SpekCalc software to cover the energy range of 100-200 keV.

[0164] Among them, the polychromatic projection simulates the X-ray penetration effect of 4 different filtered spectra (such as aluminum and copper filters), and the energy spectrum is generated by SpekCalc software; the least square estimation (LSE) algorithm is used to combine the 4 sets of polychromatic projections into the optimal fitting projection to overcome the error caused by inaccurate energy spectrum and attenuation coefficient in traditional simulation. The specific implementation simulated by SpekCalc software is, for example, 100 keV without filter, 125 keV with 2.5 mm Al, 150 keV with 1.5 mm Cu, 175 keV with 1.5 mm Cu and 2.5 mm Al.

[0165] In step S3, the following steps are included in detail:

[0166] S31: Based on the segmentation volume V seg , using the Monte Carlo photon transport model FPM, simulate the projection of polychromatic X-rays passing through the object, generate 4 sets of polychromatic projections of different energy spectra under GPU acceleration, and simulate different filter combinations by SpekCalc software; among them, the 4 energy spectra used in polychromatic projection are centered on the energy corresponding to the CBCT original scan data, and the values are taken in the energy range of 100-200 keV;

[0167] S32: Perform least square estimation (LSE) fitting, first establish a projection fitting model;

[0168]

[0169] Among them, x p,i is the i-th set of simulated projections, i is the serial number of the polychromatic projection, 1≤i≤4; c i is the weight coefficient of the i-th set of projections;

[0170] S33: To solve the optimal weight c i , the L2 norm of the real projection x μ and the simulated projection x p is minimized:

[0171]

[0172] Solve by matrix, specifically, take the derivative of the objective function and set the gradient to zero to get the coefficient solution:

[0173]

[0174] In the formula, is the simulated projection autocorrelation matrix, is the simulated projection and real projection cross-correlation vector;

[0175] S34: After obtaining the optimal weight through the fast calculation of matrix inversion, the multi-color projection simulation data x is obtained p ;

[0176] S4: Single-color projection generation and optimization, mainly including energy discretization and initial simulation, optimal energy selection; Specifically including the following operations:

[0177] Based on the segmentation volume corresponding to each material region, simulate multiple single-color projections of different energy levels, denoted as x (Ei) ; Single-color projection simulation data; Among them, different energy levels of single-color projection, for example: 80keV, 100keV, 120keV, 150keV.

[0178] By calculating the mean square error MSE of the single-color projection simulation data and the scatter-corrected projection data, select the single-color projection data with the highest similarity among all single-color projection simulation data as the optimal reference projection data, and preferentially match the high metal contrast energy level.

[0179] In step S4, the following steps are specifically included:

[0180] S41: Energy discretization and initial simulation;

[0181] Energy binning, discretize the X-ray energy spectrum (such as 100-200keV) into K energy intervals (for example, K=10, step 10keV); Perform density conversion according to the linear attenuation coefficient μ of the closest material j * Corresponding to the material j

[0182] Use the GPU-accelerated Monte Carlo model FPM to simulate the projection of single-color X-rays passing through the object, and use the input requirements of the GPU-accelerated Monte Carlo model FPM to simulate the projection of single-color X-rays passing through the object; For each energy bin E k , use FPM to generate the single-color projection x m (E k ) corresponding to each energy;

[0183] The calculation formula is as follows:

[0184]

[0185] Among them, I mono,0 (E k ) = Nη(E k), represents the intensity before penetrating the object, N is the number of incident photons, and η(E k ) represents the energy response of the detector;

[0186] is the intensity after penetrating the object, where L is the path length of the ray through the material.

[0187] Next, the optimal energy selection (MMSE principle) optimizes the single-energy projection.

[0188] S42: For each simulated monochromatic projection x m (E k ), the mean square error MSE with the true projection x u is calculated;

[0189]

[0190] where N is the total number of projection pixels;

[0191] S43: Select the energy E opt with the smallest MSE as the reference for the monochromatic projection, denoted as: optimized reference projection data;

[0192]

[0193] where the argmin function returns the argument value when the given objective function takes the minimum value.

[0194] E opt is usually located near the mean value of the energy spectrum, but it dynamically changes with the thickness of the material. This method selects monochromatic projections simulated at different energies, and iteratively selects the most accurate single-energy simulated projection. Compared with the fixed energy method, the MMSE principle optimizes the energy to reduce the projection error by 30%.

[0195] S5: Beam hardening correction of the original data; specifically including the following steps:

[0196] Nonlinear correction term generation: calculate the difference between the fitted polychromatic projection simulation data and the optimized reference projection data to obtain the nonlinear correction term;

[0197] Projection data beam hardening correction: after superimposing the nonlinear correction term on the scatter-corrected projection data, a three-dimensional image is reconstructed by FDK to obtain the corrected CBCT image.

[0198] This method constructs a joint framework for scatter preprocessing and hardening correction: the scatter component is accurately removed based on the GPU Monte Carlo model; the metal crown peripheral cup-shaped artifacts and bone tissue bright stripes caused by secondary hardening of scatter are simultaneously suppressed through the nonlinear correction term; after the composite artifacts are eliminated, the bone trabecular structure in the apical region is clearly defined, providing artifact-free images for subsequent lesion diagnosis.

[0199] In step S5, the following is included:

[0200] S51: The optimized reference projection data E opt is obtained in step S4 is used to calculate the nonlinear correction term m (E opt ) and the polychromatic projection simulation data x p obtained in step S3.

[0201] Φ CT = x m (E opt ) - x p ;

[0202] wherein Ф CT represents the compensation for the attenuation nonlinearity caused by the polychromatic spectrum, and the difference between the monochromatic and polychromatic projections is directly reflected by the correction term.

[0203] S52: Perform projection data beam hardening correction: superimpose Ф CT on the scatter-corrected data x sc to obtain the final corrected projection data x c :

[0204] x c = Φ CT + x sc ;

[0205] S53: Use the filtered back-projection FDK reconstruction algorithm to reconstruct the corrected projection data x c into a three-dimensional image V c , which is a corrected CBCT image, and the beam hardening artifacts in the reconstructed image have been suppressed.

[0206] After using the technical solution of the present application, the dependence of the conventional statistical model on the regularization parameter / spectrum weight is abandoned: the original projection data is input to automatically output the corrected DICOM image; the weighted coefficient is automatically solved by the LSE algorithm; the material estimation, projection fitting, and energy optimization are all self-determined, which is suitable for complex clinical scenarios.

Claims

1. A multi-material beam hardening correction method for oral CBCT, characterized in that, It includes the following steps: S1: Preprocessing the raw CBCT scan data for scattering artifacts; Specifically, the following operations are included: Obtain the two-dimensional projection data of the original three-dimensional image of the oral cavity CBCT, denoted as the original two-dimensional projection x. u ; Using multi-GPU acceleration based on the Monte Carlo (MC) photon transport model (FPM), the spatial distribution of scattered photons in the original image is simulated to obtain the scattering components. The scattered component is subtracted from the original two-dimensional projection data using iterative subtraction, and the purified projection data after scattering correction is output, denoted as: scattering-corrected projection data x. sc ; S2: Blind estimation of material for reconstructed objects; Specifically, the following operations are included: The initial three-dimensional volume is obtained by filtering back projection FDK reconstruction of the scattering-corrected projection data; the adaptive gray-scale threshold segmentation method is used to separate the different material regions in the volume. The pre-built lookup table of attenuation coefficients for oral materials is invoked, and the material type of each region is iteratively matched using the minimum mean square error (MMSE) criterion. The material type of each region is labeled, and the segmented volume after material labeling is output. S3: Multispectral simulation of multicolor projection; Specifically, the following operations are included: Based on the segmented volume corresponding to each material region, multi-color projection data is synchronously simulated using a multi-GPU parallel architecture. The least squares estimation (LSE) algorithm is used to weight and combine all the multicolor projections into the best-fit projection, denoted as: multicolor projection simulation data; S4: Monochrome projection generation and optimization; Specifically, the following operations are included: Based on the segmented volume corresponding to each material region, multiple monochromatic projections of different energy levels are simulated, denoted as: Monochromatic projection simulation data; By calculating the mean square error (MSE) between the monochrome projection simulation data and the scattering-corrected projection data, the monochrome projection data with the highest similarity to the original data is selected as the optimized baseline projection data from all the monochrome projection simulation data. S5: Raw data beam hardening correction; Specifically, the following steps are included: The difference between the fitted multicolor projection simulation data and the optimized reference projection data is calculated to obtain the nonlinear correction term; After superimposing the nonlinear correction term onto the scattering-corrected projection data, the three-dimensional image is reconstructed using FDK to obtain the corrected CBCT image.

2. The method for multi-material beam hardening correction in oral CBCT according to claim 1, characterized in that: Step S1 includes the following steps in detail: S11: Obtain the original two-dimensional projection x of the oral CBCT scan u ; S12: Model the distribution of scattered photons; By loading the GPU-accelerated Monte Carlo photon transport model FPM, based on the original two-dimensional data x u Back-projection reconstruction generates an initial attenuation coefficient volume to simulate the spatial distribution of scattered photons in the original image. To accelerate the scattered photon simulation process, GPU parallelization is used to speed up the output of the scattered photon distribution. The resulting simulated scattered projection distribution x scatter The initial value estimate is: Where, x sc (0) The initial scattering correction projection is given by α, where α is the initial step size and x is the initial step size. u The raw two-dimensional projection data of the CBCT scan image; S13: Subtract the scattering component from the original two-dimensional projection data through iterative subtraction, and output the scatter-corrected pure projection data, denoted as: scatter-corrected projection data x. sc ; The iterative correction formula is: Where k is the number of iterations; β k For adaptive step size, it decreases with each iteration; FPM scatter This is a GPU-MC scattering estimation operator.

3. The method for multi-material beam hardening correction in oral CBCT according to claim 1, characterized in that: In step S2, details Includes the following steps: S21: For the scattering-corrected projection data x sc Perform filtered back projection FDK reconstruction to obtain the initial 3D volume V. sc ; S22: An adaptive grayscale threshold segmentation method is used for the initial three-dimensional volume V. sc The different material regions in the image are divided using the following method: By employing a multi-threshold segmentation method, thresholds T1 and T2 are obtained using the Ostu multi-threshold method to automatically divide the reconstructed data into three regions. In the formula, V soft V represents a soft tissue region. bone V represents the skeletal and dental region. metal The region is the metal implantation area; 'a' represents the attenuation coefficient value within the initial three-dimensional volume. S23: Perform morphological closing operations on each segmented region to eliminate segmentation noise; S24: Read the pre-constructed lookup table of attenuation coefficients for dental materials; The dental material attenuation coefficient lookup table records the density and attenuation coefficient of different types of materials at different energies; S25: Based on the oral material attenuation coefficient lookup table, iteratively match each segmented region V using the minimum mean square error (MMSE) criterion. i Corresponding material types; Specifically, the following steps are included: S251: For each segmented region V i Calculate the average projection value; In the formula, Indicates each region V i The average projection value, mean() is the average value function; FPLS is a fast linear integral model that directly projects the segmented single material volume along the incident direction of the ray and calculates the sum of the attenuation values ​​of the ray passing through the voxel; S252: After thickness compensation, find the average projection value in the oral cavity material attenuation coefficient lookup table. The minimum Euclidean distance of the closest materials, obtained through the minimum mean square error (MMSE) criterion, is used for each segmented region V. i The method for finding the closest material is as follows: Where, j * The argmin function returns the value of the independent variable when the given objective function reaches its minimum, representing the closest material. The MMSE criterion is based on the L2 norm and is expressed as ||.||2. j represents different materials that exist in the dental material attenuation coefficient lookup table, μ j (E ref () indicates that material j at the reference energy E ref The attenuation coefficient used below; t i To divide the region V i Average penetration thickness; E ref The reference energy is the same as the scan energy value corresponding to the original scan data; S26: Determine each segmented region V i After specifying the corresponding material, the material type of each segmented region is labeled, and the segmented volume V of each region after material labeling is output. seg .

4. The method for multi-material beam hardening correction in oral CBCT according to claim 3, characterized in that: Before step S26 is performed, it is also necessary to process each segmented region V. i The corresponding material type needs to be corrected and confirmed, specifically including the following steps; a1: In the dental material attenuation coefficient lookup table, confirm according to j * The corresponding material density ρ est and linear attenuation coefficient μ est Calculate the density-based Monte Carlo projection x FPM ; x FPM =FPM(ρ est ); Where, ρ est Where FPM is the material density, and FPM is the Monte Carlo photon transport model. a2: Calculate the linear projection x based on the attenuation coefficient FPLS : x FPLS =FPLS(μ est ); Wherein, FPLS is the fast linear integral model, μ est The linear attenuation coefficient is used. a3: Construct the objective optimization function: a4: Minimize the residual of the objective function, and output the optimized material label volume V. seg .

5. The method for multi-material beam hardening correction in oral CBCT according to claim 3, characterized in that: Step S3 includes the following steps in detail: S31: Based on the segmented volume V corresponding to each material region seg Using a multi-GPU parallel architecture and the Monte Carlo photon transport model (FPM), the projection of multicolor X-rays through an object was simulated. Four sets of multicolor projections with different energy spectra were generated under GPU acceleration, and different filter combinations were simulated using SpekCalc software. The four energy spectra used for the multicolor projection were centered on the energy corresponding to the original CBCT scan data and were taken in the energy range of 100-200 keV. S32: Perform least squares estimation LSE fitting, first establish a projection fitting model; Where, x p,i Let c be the i-th simulated projection, where i is the index of the multicolor projection, 1≤i≤4; i Let be the weighting coefficients for the i-th projection; S33: To solve for the optimal weight c i By constructing a model that minimizes the true projection x μ With analog projection x p L2 norm: Solving by matrix: In the formula, To simulate the projected autocorrelation matrix, The cross-correlation vector between the simulated projection and the real projection; S34: After quickly calculating the optimal weights through matrix inversion, the multicolor projection simulation data x is obtained. p .

6. The method for multi-material beam hardening correction in oral CBCT according to claim 5, characterized in that: Step S4 specifically includes the following steps: S41: Energy Discretization and Initial Simulation; The X-ray energy spectrum is discretized into K energy ranges; density conversion is performed, based on the closest material j. * The corresponding linear decay coefficient μ is converted into mass density ρ; A GPU-accelerated Monte Carlo model FPM is used to simulate the projection of monochromatic X-rays through an object. The input requirements for the GPU-accelerated Monte Carlo model FPM are used to simulate monochromatic projection; for each energy box E... k Using FPM to generate a monochrome projection x m (E k ); The calculation formula is as follows: Among them, I mono,0 (E k )=Nη(E k ), representing the intensity before penetrating the object, where N is the number of incident photons, η(E k () indicates the detector's energy response; The intensity of the ray after penetrating the object is L, where L is the length of the path the ray travels through the material. S42: For each analog monochrome projection x m (E k ), calculate and the true projection x u The mean square error (MSE); Where N is the total number of projected pixels; S43: Select the energy E with the minimum MSE. opt As a monochrome projection reference, it is denoted as: optimized reference projection data; The argmin function returns the value of the argument when the given objective function reaches its minimum.

7. The method for multi-material beam hardening correction in oral CBCT according to claim 6, characterized in that: Step S5 includes the following: S51: The optimized reference projection data is denoted as: x m (E opt ), and the multicolor projection simulation data x p Calculate the nonlinear correction term; F CT =x m (E opt )-x p ; Among them, Ф CT This indicates the attenuation nonlinearity caused by compensating for the multicolor energy spectrum; S52: Perform beam hardening correction on the projection data: adjust Ф CT Superimposed on the scattering-corrected data x sc The final corrected projection data x is obtained. c : x c =Φ CT +x sc ; S53: Use the filtered back projection FDK reconstruction algorithm to reconstruct the corrected projection data x c Reconstructed into a 3D image V c , to the corrected CBCT image.

Images