Oral cavity curved surface tomography method and device based on GPU hardware acceleration

Abstract

The invention relates to an oral cavity curved surface tomography method based on GPU hardware acceleration. The method comprises the following steps: S1, scanning along a pre-designed motion track to obtain strip-shaped projection data; s2, obtaining a plurality of groups of splicing coefficients based on the splicing coefficients and the offset; s3, calculating the size of the curved surface fault image according to the splicing coefficient; s4, determining the number of thread blocks and thread grids in the GPU according to the image size; s5, obtaining coordinates of the to-be-interpolated point based on the number of the strip projection data and a value obtained by accumulating the splicing coefficient; s6, all the points to be interpolated are distributed to each GPU core for interpolation accumulation calculation; and S7, calculating the number of times of accumulating the projection data by each interpolation point, and carrying out standardization processing to obtain a curved surface sectional image. According to the method, each coordinate point of the panoramic image serves as an acceleration object and is placed in a GPU core for parallel calculation, and the speed is higher than that of a traditional moving superposition splicing method; by standardizing the position of each pixel, the problem that spliced stripes and the brightness of each region are inconsistent is solved, and images with consistent gray levels are obtained.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image processing, and in particular to an oral curved surface tomography method and device based on GPU hardware acceleration. Background Art

[0002] In an oral curved surface tomography system, in order to obtain dental images with high resolution and contrast, a pre-set basic stitching coefficient is generally used. Multiple sets of stitching coefficients are obtained by offsetting the basic stitching coefficients. Using a moving superposition method, multiple sets of oral curved surface tomography images are obtained respectively. The multiple sets of curved surface tomography images obtained can be standardized to obtain curved surface tomography images with different depth information. In this process, when multiple sets of images need to be stitched together, a long processing time is required. Summary of the invention

[0003] The technical problem to be solved by the present invention is to design an oral curved surface tomography method and device, based on GPU hardware acceleration, taking each coordinate point of the panoramic image as an acceleration object, placing it in the GPU core for parallel calculation, and solving the existing technical problems.

[0004] In order to solve the above technical problems, the present invention provides an oral curved surface tomography method based on GPU hardware acceleration, which specifically includes the following steps:

[0005] Step S1: Scan along a pre-designed motion trajectory to obtain strip projection data;

[0006] Step S2: obtaining multiple groups of splicing coefficients based on given splicing coefficients and offsets;

[0007] Step S3: Calculate the size of the curved tomographic image according to the stitching coefficient;

[0008] Step S4: Determine the number of thread blocks and thread grids in the GPU according to the image size of the curved surface, and calculate the coordinates of each pixel (iwidth w ,iheight w );

[0009] Step S5: obtaining the coordinates of the point to be interpolated based on the number of strip projection data and the value obtained by accumulating the splicing coefficient;

[0010] Step S6: Allocate all the points to be interpolated to each GPU core for interpolation accumulation calculation, and calculate the coordinates (iwidth w ,iheight w ) grayscale value pano w ;

[0011] Step S7: Calculate the number of times addNum the projection data is accumulated at each interpolation point w , perform the following normalization on the curved tomographic image: panoNormalize w =pano w / addNum w , and obtain the spliced ​​curved tomographic image.

[0012] Furthermore, in step S1, after the projection data is collected, dark field correction and bright field gain correction need to be performed on the projection data. The correction method is as follows:

[0013] Step S11: Collect n air bright field projection images, and the kth bright field image is recorded as Bright k ;

[0014] Step S12: Collect n dark-field projection images, and the kth dark-field image is recorded as Dark k ;

[0015] Step S13: Calculate the average value of n air bright field projection images to obtain an average bright field image:

[0016]

[0017] Step S14: Calculate the average value of n dark field projection images to obtain an average dark field image:

[0018]

[0019] Step S15: Assume that the size of the projected image is M rows*N columns, and the mean value of the air after dark field correction is:

[0020]

[0021] Step S16: Calculate the gain correction coefficient:

[0022] Gain i,j =MeanBrightValue / (MeanBright i,j -MeanDark i,j );

[0023] Step S17: For the collected V s Projection image, the corrected projection image is:

[0024]

[0025] Furthermore, step S2 specifically includes the following steps:

[0026] Step S21: Calculate the reference splicing coefficient: Calculate Nviews-1 sets of splicing coefficients according to the number of acquired frames Nviews, and add 0 to obtain Nviews splicing coefficients, where x0=0, x s+1 =x, then we can get a set of benchmark splicing coefficients: X base ={x s};

[0027] Step S22: obtaining multiple groups of splicing coefficients based on the reference splicing coefficients: assuming that the variation of the splicing coefficients is δ, in order to obtain W groups of splicing coefficients, W is an odd number, then the number of half groups is halfW=(W-1) / 2, and the splicing coefficients are:

[0028] X w =(w-halfW)*δ*X base +X base , where w = 0, 1, 2,…, W-1.

[0029] Further, in step S3, according to the calculated splicing coefficient X w ={x w,s}, where w = 0, 1, 2, ..., W-1, s = 0, 1, 2, ..., Nviews-1, and the size of the multi-layer panoramic image PanoWidth is calculated w *PanoHeight w ,in:

[0030] N is the column width of the strip projection image;

[0031] PanoHeight w =M, where M is the row height of the strip projection image.

[0032] Furthermore, in step S4, CUDA is used to accelerate parallel computing:

[0033] Set thread block Threads w =(BLOCKDIM_COL, BLOCKDIM_ROW), where BLOCKDIM_COL is the number of thread blocks allocated in the column direction of the surface fault panorama, and BLOCKDIM_ROW is the number of thread blocks allocated in the row direction of the surface fault panorama.

[0034] The number of grids is calculated based on the number of thread blocks:

[0035] Grids w =(Ceil(PanoWidth w / BLOCKDIM_COL),Ceil(PanoHeight w / BLOCKDIM_ROW)), where Ceil represents the round-up function, PanoWidth w PanoHeight is the width of the stitched image size. w The height of the image after stitching.

[0036] Furthermore, the stitching image coordinates (iwidth w ,iheight w ) is calculated based on thread blocks and thread grids, specifically:

[0037] iwidth w =threadIdx.x+blockIdx.x*bolckDim.x;

[0038] iheight w =threadIdx.y+blockIdx.y*bolckDim.y;

[0039] Among them, threadIdx.x is the index of the current thread in the thread block in the column direction, blockIdx.x is the index of the current thread block in the column direction, and blockDim.x is the number of thread blocks in the column direction;

[0040] threadIdx.y is the index of the current thread in the thread block in the row direction, blockIdx.y is the index of the current thread block in the row direction, and blockDim.y is the number of thread blocks in the row direction.

[0041] Furthermore, step S5 specifically includes:

[0042] Step S51: Calculate the cumulative sequence cumSumX of all splicing coefficients w :

[0043] Let cumSumX w = {cumSum_x w,s}, Where w=0,1,2,…,W-1,s=0,1,2,…,Nviews-1;

[0044] Step S52: Calculate the difference between the horizontal coordinate of the pixel and the accumulated sequence of the splicing coefficients:

[0045] ε w,s =iwidth w -cumSum_x w,s ;

[0046] Step S53: Determine the range of the difference and determine the coordinates to be interpolated in the projection. If 0≤ε w,s ≤N-1, then εw,s For the projected image V s Overlay to the stitched image coordinates (iwidth w ,iheight w ) is the horizontal coordinate value at .

[0047] Furthermore, in step S6, the interpolation weight of the horizontal coordinate linear interpolation is:

[0048] weight Ceil =ε w,s -floor(ε w,s ), weight Ceil is the interpolation weight of the coordinate grayscale after rounding up;

[0049] weight Floor =1-weight Ceil , weight Floor is the interpolation weight of the coordinate grayscale after rounding down;

[0050] The stitching surface tomographic image is in coordinate (iwidth w ,iheight w ) is:

[0051] Among them, pano w (iwidth w ,iheight w ) is the coordinate in the curved tomographic image obtained after the wth group of splicing coefficients are spliced ​​(iwidth w ,iheight w )’s grayscale value.

[0052] Further, in step S7, when 0≤ε w,s ≤N-1, addNum w (iwidth w ,iheight w )=addNum w (iwidth w ,iheight w )+1 until all the projection data are cycled through.

[0053] The present invention also provides an oral curved surface tomography device based on GPU hardware acceleration, comprising:

[0054] at least one processor; and

[0055] at least one memory in communication with the processor;

[0056] The memory stores instructions executable by the processor, and the instructions are executed by the processor to enable the device to perform the aforementioned oral curved surface tomography imaging method based on GPU hardware acceleration.

[0057] Beneficial effects of the present invention:

[0058] The oral curved surface tomography imaging method of the present invention is based on GPU hardware acceleration, takes each coordinate point of the panoramic image as an acceleration object, and places it in the GPU core for parallel calculation, which is faster than the traditional mobile superposition stitching method. The present invention solves the problem of inconsistent brightness of the splicing stripes and each area by standardizing the position of each pixel, and can obtain an image with consistent grayscale, solving the technical problem that the mobile superposition stitching method causes large grayscale differences at different positions when the splicing coefficients are large during stitching. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] The specific implementation of the present invention will be further explained below in conjunction with the accompanying drawings.

[0060] Figure 1 It is a flow chart of the oral curved surface tomography method based on GPU hardware acceleration of the present invention;

[0061] Figure 2 A schematic diagram for collecting projection data;

[0062] Figure 3 A schematic diagram of basic splicing coefficient calculation;

[0063] Figure 4 It is a schematic diagram of the movement of a certain point A in the dental arch line on the detector;

[0064] Figure 5 This is a schematic diagram comparing the strip projection data (intercepted part) before and after correction;

[0065] FIG6( a ) is a schematic diagram of a set of rotation speeds in Example 1;

[0066] FIG6( b ) is a schematic diagram of a set of rotation speeds in Example 1;

[0067] FIG6(c) is a schematic diagram of a set of reference splicing coefficients calculated based on the rotation speed in FIG6(a) and the translation speed in FIG6(b) in Example 1;

[0068] Figure 7 Schematic diagram of the comparison of stitched images before and after standardization. DETAILED DESCRIPTION

[0069] Example 1

[0070] Combination Figure 1The oral curved surface tomography method based on GPU hardware acceleration of this embodiment includes the following steps:

[0071] Step S1: Scan along a pre-designed motion trajectory to obtain strip projection data, and number the projection data according to the order in which they are obtained.

[0072] Specifically, Figure 2 A schematic diagram of collecting projection data is given in the figure. The X-Ray light source and detector are fixed on a rotating arm. The distance between them remains unchanged. The center of rotation can be moved along the panoramic axis. The entire motion trajectory is a combination of rotation and linear motion.

[0073] In this embodiment, the sampling frame rate is set to 200 Hz and the exposure time is 14 seconds. Then, 2800 frames of projection data are collected in the whole process. Let Nviews = 2800. s Indicates the projection of the sth frame, 0≤s <Nviews。

[0074] In this embodiment, preferably, in step S1, after the projection data is collected, dark field correction and bright field gain correction need to be performed on the projection data, and the correction method is as follows:

[0075] Step S11: Collect n air bright field projection images, and the kth bright field image is recorded as Bright k ;

[0076] Step S12: Collect n dark-field projection images, and the kth dark-field image is recorded as Dark k ;

[0077] Step S13: Calculate the average value of n air bright field projection images to obtain an average bright field image:

[0078]

[0079] Step S14: Calculate the average value of n dark field projection images to obtain an average dark field image:

[0080]

[0081] Step S15: Assume that the size of the projected image is M rows*N columns, and the mean value of the air after dark field correction is:

[0082]

[0083] Step S16: Calculate the gain correction coefficient:

[0084] Gain i,j =MeanBrightValue / (MeanBright i,j -MeanDarki,j );

[0085] Step S17: For the collected V s Projection image, the corrected projection image is:

[0086]

[0087] In this way, the processed projection data is obtained. After dark field correction and gain correction, the corrected projection data is obtained as follows Figure 5 As shown, Figure 5 (a) is the strip projection data before correction, and (b) is the strip projection data after correction. The comparison before and after correction shows that the constant lines in the strip image are corrected and the image contrast is higher.

[0088] Step S2: Based on the given stitching coefficients and offsets, multiple groups of stitching coefficients are obtained, and each stitching coefficient in each group corresponds to a strip projection sequence number.

[0089] Specifically, it is necessary to calculate the benchmark splicing coefficient first, such as Figure 3 As shown in the figure, assuming that the ray source S1 moves a small angle θ, the ray source is at position S2 after the movement, the rotation center moves from O1 to O2, and the intersection of the ray and the dental arch moves from point A to point B. Then the projection position of point A in the detector after the movement becomes point E, and the stitching coefficient is Figure 2 The x shown in is the length of line segment DE. According to the similarity principle of triangles, we can get:

[0090] AM = AC sinθ;

[0091] S2D = source_to_detector;

[0092] S2M=S2SO2+O2M=source_to_rotation+O2M;

[0093] AC=d+c, where d is a known quantity, namely the distance from the center of rotation to the dental arch line.

[0094] O2M=f, h+f=d, then f=dh;

[0095] O1N=esinα, e is the distance between two adjacent rotation centers, which is a known quantity;

[0096]

[0097] a=ecosα;

[0098]

[0099] h=CB-CM=(d+a+b)-(d+c)cosθ.

[0100] Finally, the splicing coefficient can be calculated as:

[0101]

[0102] where θ, e and α can be calculated based on the rotation speed R vot and the speed P of the panoramic axis vot It is calculated, where α is the angle between the vertical line from the light source to the detector and the panoramic axis Y axis. In this embodiment, the initial angle of this angle is -112 degrees. The angle between each sampling position and the Y axis can be calculated according to the rotation speed.

[0103] In this embodiment, preferably, step S2 specifically includes the following steps:

[0104] Step S21: Calculate the reference splicing coefficient: Calculate Nviews-1 groups of splicing coefficients according to the above calculation method and the number of frames collected Nviews, and add 0 to obtain Nviews splicing coefficients, where x0=0, x s+1 =x, then we can get a set of benchmark splicing coefficients: X base ={x s}; Combined with Figure 6(a) and Figure 6(b), according to a given set of rotation speeds and translation speeds, the horizontal axis in Figure 6(a) is the number of sampling points, the speed of every 10 points is the same, and the vertical axis is the rotation speed; the horizontal axis in Figure 6(b) is the number of sampling points, the speed of every 10 points is the same, and the vertical axis is the translation speed; a set of reference stitching coefficients can be calculated through the above speeds, as shown in Figure 6(c), where the horizontal axis in Figure 6(c) represents the number of scanned sampling points, and the vertical axis represents the number of pixels translated between each frame.

[0105] Step S22: obtaining multiple groups of splicing coefficients based on the reference splicing coefficients: assuming that the variation of the splicing coefficients is δ, in order to obtain W groups of splicing coefficients, W is an odd number, then the number of half groups is halfW=(W-1) / 2, and the splicing coefficients are:

[0106] X w =(w-halfW)*δ*X base +X base , where w = 0, 1, 2,…, W-1.

[0107] Step S3: Calculate the size of the curved tomographic image according to the stitching coefficient.

[0108] In this embodiment, preferably, in step S3, according to the calculated splicing coefficient X w ={x w,s}, where w = 0, 1, 2, ..., W-1, s = 0, 1, 2, ..., Nviews-1, because stitching is only done horizontally and not vertically, it is only necessary to calculate the width of the image after horizontal stitching. According to the principle of moving superposition, N is the column width of the strip projection image; PanoHeight w =M, M is the row height of the strip projection image. Calculate the size of the multi-layer panoramic image PanoWidth w *PanoHeight w .

[0109] Step S4: Determine the number of thread blocks and thread grids in the GPU according to the image size of the curved surface, and calculate the coordinates of each pixel (iwidth w ,iheight w ).

[0110] Preferably, this embodiment uses CUDA for parallel computing acceleration. When CUDA is executed, each Kernel function in the Host program is executed on the GPU according to the concept of thread grid. One Kernel function corresponds to one thread grid. The tasks in each thread grid are certain. When these tasks need to be executed, each thread grid divides the tasks into parts of thread blocks. There are several threads in the middle of the thread block, which are then distributed to the threads to complete.

[0111] Set thread block Threads w =(BLOCKDIM_COL, BLOCKDIM_ROW), where BLOCKDIM_COL is the number of thread blocks allocated in the column direction of the surface fault panorama, and BLOCKDIM_ROW is the number of thread blocks allocated in the row direction of the surface fault panorama. For example, these two values ​​can be set to 8.

[0112] The number of grids is calculated based on the number of thread blocks:

[0113] Grids w =(Ceil(PanoWidth w / BLOCKDIM_COL),Ceil(PanoHeight w / BLOCKDIM_ROW)), where Ceil represents the round-up function, PanoWidth w PanoHeight is the width of the stitched image size. w The height of the image after stitching.

[0114] In this embodiment, preferably, the stitching image coordinates (iwidth w,iheight w ) is calculated based on thread blocks and thread grids, specifically:

[0115] iwidth w =threadIdx.x+blockIdx.x*bolckDim.x;

[0116] iheight w =threadIdx.y+blockIdx.y*bolckDim.y;

[0117] Among them, threadIdx.x is the index of the current thread in the thread block in the column direction, blockIdx.x is the index of the current thread block in the column direction, and blockDim.x is the number of thread blocks in the column direction;

[0118] threadIdx.y is the index of the current thread in the thread block in the row direction, blockIdx.y is the index of the current thread block in the row direction, and blockDim.y is the number of thread blocks in the row direction.

[0119] Step S5: based on the number of strip projection data and the value obtained by accumulating the splicing coefficient, the coordinates of the point to be interpolated are obtained. The coordinates of the moving superposition interpolation point are calculated based on the splicing coefficient.

[0120] In this embodiment, preferably, step S5 specifically includes:

[0121] Step S51: Calculate the cumulative sequence cumSumX of all splicing coefficients w :

[0122] Let cumSumX w = {cumSum_x w,s}, Where w=0,1,2,…,W-1,s=0,1,2,…,Nviews-1;

[0123] Step S52: Calculate the difference between the horizontal coordinate of the pixel and the accumulated sequence of the splicing coefficients:

[0124] ε w,s =iwidth w -cumSum_x w,s ;

[0125] Step S53: Determine the range of the difference and determine the coordinates to be interpolated in the projection. If 0≤ε w,s ≤N-1, then ε w,s For the projected image V s Overlay to the stitched image coordinates (iwidth w ,iheightw ), wherein N represents the size of the bar graph, that is, the number of columns N of the size of the projection image in step S15.

[0126] Step S6: Allocate all the points to be interpolated to each GPU core for interpolation accumulation calculation, and calculate the coordinates (iwidth w ,iheight w ) grayscale value pano w .

[0127] In this embodiment, preferably, in step S6, the interpolation weight of the horizontal coordinate linear interpolation is:

[0128] weight Ceil =ε w,s -floor(ε w,s ), weight Ceil is the interpolation weight of the coordinate grayscale after rounding up;

[0129] weight Floor =1-weight Ceil , weight Floor is the interpolation weight of the coordinate grayscale after rounding down;

[0130] The stitching surface tomographic image is in coordinate (iwidth w ,iheight w ) is:

[0131] Among them, pano w (iwidth w ,iheight w ) is the coordinate in the curved tomographic image obtained after the wth group of splicing coefficients are spliced ​​(iwidth w ,iheight w )’s grayscale value.

[0132] Step S7: Calculate the number of times addNum the projection data is accumulated at each interpolation point w ,In order to eliminate the problem of large image grayscale differences caused by different stacking times, the curved tomographic image is standardized as follows: panoNormalize w =pano w / addNum w , and obtain the spliced ​​curved tomographic image.

[0133] In this embodiment, preferably, in step S7, when 0≤ε w,s ≤N-1, addNum w(iwidth w ,iheight w )=addNum w (iwidth w ,iheight w )+1, until all the projection data are cycled through, wherein N represents the size of the bar graph, that is, the number of columns N of the size of the projection image in step S15.

[0134] Specifically, Figure 4 As shown, suppose there is a point A on the dental arch curve (this point corresponds to a pixel point in the curved tomographic image). With the rotation of the light source and the detector, and the movement of the panoramic axis, point A moves at the detector position along with the rotation angle and the movement of the rotation center of the panoramic axis. Since the size of the detector is limited, point A will not be displayed in all projection images. Since the rotation speed and the translation speed of the rotation center are non-uniform, the grayscale of different points on the dental arch curve (different pixel points in the curved tomographic image) after superposition will be different.

[0135] Assume that the coordinates of point A in the curved tomographic image are (iwidth w ,iheight w ), you can use the column-wise coordinate iwidth w Calculate ε w,s =iwidth w -cumSum_x w,s , judge ε w,s Is it 0≤ε? w,s ≤N-1, if the above conditions are met, point A is in the Vth s There is a projection in the projection, set addNum w (iwidth w ,iheight w ) is the number of images of the stored point A in the projection images. For example, for point A, there are m projection images of point A, then addNum w (iwidth w ,iheight w )=m, let the gray level of point A after superimposing m projection images be G A , then the standardized value NormalG A =G A / m.

[0136] Based on the above steps, W stitched curved tomographic images can be quickly obtained.

[0137] Based on the projection data and reference stitching coefficients obtained in this embodiment, 21 groups of stitching calculations are compared, and the speed comparison between the traditional CPU serial calculation and the GPU parallel calculation of the present invention is as follows:

[0138]

[0139] As can be seen from the table, the GPU-based acceleration algorithm of the present invention processes a single-layer image 15 times faster than the CPU serial method; and processes 21-layer images 7 times faster than the CPU serial method.

[0140] In addition, the present invention can obtain an image with consistent grayscale by standardizing the position of each pixel. Figure 7 As shown, Figure 7 (a) is the unstandardized stitched image, and (b) is the standardized image. It can be seen that before standardization, there are obvious stitching stripes and light and dark boundaries. The standardized image can eliminate the problems of stitching stripes and inconsistent brightness of each area.

[0141] Example 2

[0142] The oral curved surface tomography device based on GPU hardware acceleration of this embodiment includes:

[0143] at least one processor; and

[0144] at least one memory in communication with the processor;

[0145] The memory stores instructions that can be executed by the processor, and the instructions are executed by the processor to enable the device to perform the oral curved surface tomography imaging method based on GPU hardware acceleration in Example 1.

[0146] Many specific details are described in the above description to facilitate a full understanding of the present invention. However, the above description is only a preferred embodiment of the present invention. The present invention can be implemented in many other ways different from those described herein, so the present invention is not limited to the specific implementation disclosed above. At the same time, any person familiar with the art can make many possible changes and modifications to the technical solution of the present invention using the methods and technical contents disclosed above without departing from the scope of the technical solution of the present invention, or modify it into an equivalent embodiment of equivalent changes. Any simple modification, equivalent change and modification made to the above embodiments according to the technical essence of the present invention without departing from the content of the technical solution of the present invention still falls within the scope of protection of the technical solution of the present invention.

Claims

1. An oral curved surface tomography method based on GPU hardware acceleration, characterized in that: The steps include: Step S1: Scan along a pre-designed motion trajectory to obtain strip projection data; Step S2: obtaining multiple groups of splicing coefficients based on given splicing coefficients and offsets; Step S3: Calculate the size of the curved tomographic image according to the stitching coefficient; Step S4: Determine the number of thread blocks and thread grids in the GPU according to the image size of the curved surface, and calculate the coordinates of each pixel (iwidth w ,iheight w ); Step S5: obtaining the coordinates of the point to be interpolated based on the number of strip projection data and the value obtained by accumulating the splicing coefficient; Step S6: Allocate all the points to be interpolated to each GPU core for interpolation accumulation calculation, and calculate the coordinates (iwidth w ,iheight w ) grayscale value pano w ; Step S7: Calculate the number of times addNum the projection data is accumulated at each interpolation point w , perform the following normalization on the curved tomographic image: panoNormalize w =pano w / addNum w , and obtain the spliced ​​curved tomographic image.

2. The method for oral curved surface tomography based on GPU hardware acceleration according to claim 1, characterized in that: In step S1, after the projection data is collected, dark field correction and bright field gain correction need to be performed on the projection data. The correction method is as follows: Step S11: Collect n air bright field projection images, and the kth bright field image is recorded as Bright k ; Step S12: Collect n dark-field projection images, and the kth dark-field image is recorded as Dark k ; Step S13: Calculate the average value of n air bright field projection images to obtain an average bright field image: Step S14: Calculate the average value of n dark field projection images to obtain an average dark field image: Step S15: Assume that the size of the projected image is M rows*N columns, and the mean value of the air after dark field correction is: Step S16: Calculate the gain correction coefficient: Gain i,j =MeanBrightValue / (MeanBright i,j -MeanDark i,j ); Step S17: For the collected V s Projection image, the corrected projection image is:

3. The method for oral curved surface tomography based on GPU hardware acceleration according to claim 1, characterized in that: Step S2 specifically includes the following steps: Step S21: Calculate the reference splicing coefficient: Calculate Nviews-1 sets of splicing coefficients according to the number of acquired frames Nviews, and add 0 to obtain Nviews splicing coefficients, where x0=0, x s+1 =x, then we can get a set of benchmark splicing coefficients: X base ={x s }; Step S22: obtaining multiple groups of splicing coefficients based on the reference splicing coefficients: assuming that the variation of the splicing coefficients is δ, in order to obtain W groups of splicing coefficients, W is an odd number, then the number of half groups is halfW=(W-1) / 2, and the splicing coefficients are: X w = (w - halfW) * δ * X base + X base , where w = 0, 1, 2, …, W - 1.

4. The method for oral curved surface tomography based on GPU hardware acceleration according to claim 3, characterized in that: In step S3, according to the calculated splicing coefficient X w ={x w,s }, where w = 0, 1, 2, ..., W-1, s = 0, 1, 2, ..., Nviews-1, and the size of the multi-layer panoramic image PanoWidth is calculated w *PanoHeight w ,in: N is the column width of the strip projection image; PanoHeight w =M, where M is the row height of the strip projection image.

5. The method for oral curved surface tomography based on GPU hardware acceleration according to claim 1, characterized in that: In step S4, CUDA is used to accelerate parallel computing: Set thread block Threads w =(BLOCKDIM_COL, BLOCKDIM_ROW), where BLOCKDIM_COL is the number of thread blocks allocated in the column direction of the surface fault panorama, and BLOCKDIM_ROW is the number of thread blocks allocated in the row direction of the surface fault panorama. The number of grids is calculated based on the number of thread blocks: Grids w =(Ceil(PanoWidth w / BLOCKDIM_COL),Ceil(PanoHeight w / BLOCKDIM_ROW)), where Ceil represents the round-up function, PanoWidth w PanoHeight is the width of the stitched image size. w The height of the image after stitching.

6. The method for oral curved surface tomography based on GPU hardware acceleration according to claim 5, characterized in that: Stitching image coordinates (iwidth w ,iheight w ) is calculated based on thread blocks and thread grids, specifically: iwidth w =threadIdx.x+blockIdx.x*bolckDim.x; iheight w =threadIdx.y+blockIdx.y*bolckDim.y; Among them, threadIdx.x is the index of the current thread in the thread block in the column direction, blockIdx.x is the index of the current thread block in the column direction, and blockDim.x is the number of thread blocks in the column direction; threadIdx.y is the index of the current thread in the thread block in the row direction, blockIdx.y is the index of the current thread block in the row direction, and blockDim.y is the number of thread blocks in the row direction.

7. The method for oral curved surface tomography based on GPU hardware acceleration according to claim 6, characterized in that: Step S5 specifically includes: Step S51: Calculate the cumulative sequence cumSumX of all splicing coefficients w : Let cumSumX w = {cumSum_x w,s }, Where w=0,1,2,…,W-1,s=0,1,2,…,Nviews-1; Step S52: Calculate the difference between the horizontal coordinate of the pixel and the accumulated sequence of the splicing coefficients: ε w,s =iwidth w -cumSum_x w,s ; Step S53: Determine the range of the difference and determine the coordinates to be interpolated in the projection. If 0≤ε w,s ≤N-1, then ε w,s For the projected image V s Overlay to the stitched image coordinates (iwidth w ,iheight w ) is the horizontal coordinate value at .

8. The method for oral curved surface tomography based on GPU hardware acceleration according to claim 7, characterized in that: In step S6, the interpolation weight of the horizontal coordinate linear interpolation is: weight Ceil =ε w,s -floor(ε w,s ), weight Ceil is the interpolation weight of the coordinate grayscale after rounding up; weight Floor =1-weight Ceil , weight Floor is the interpolation weight of the coordinate grayscale after rounding down; The stitching surface tomographic image is in coordinate (iwidth w ,iheight w ) is: Among them, pano w (iwidth w ,iheight w ) is the coordinate in the curved tomographic image obtained after the wth group of splicing coefficients are spliced ​​(iwidth w ,iheight w )’s grayscale value.

9. The method for oral curved surface tomography based on GPU hardware acceleration according to claim 7, characterized in that: In step S7, when 0≤ε w,s ≤N-1, addNum w (iwidth w ,iheight w )=addNum w (iwidth w ,iheight w )+1 until all the projection data are cycled through.

10. An oral curved surface tomography device based on GPU hardware acceleration, characterized in that: include: at least one processor; as well as at least one memory in communication with the processor; The memory stores instructions executable by the processor, and the instructions are executed by the processor so that the device executes the oral curved surface tomography imaging method based on GPU hardware acceleration as described in any one of claims 1 to 9.

Images