Spiral CT reconstruction method and device, computer equipment, storage medium and computer program product

Abstract

The invention relates to a spiral CT reconstruction method and device, computer equipment, a storage medium and a computer program product. The method comprises the following steps: filtering projection data to obtain filtered data; performing back projection processing on the filtering data to obtain a back projection image, and extracting bone tissue data from the back projection image; performing orthographic projection on the bone tissue data to obtain bone projection data; taking a difference value between the projection data and the bone projection data as bone tissue removed data; and obtaining a CT reconstruction image according to the bone tissue removed data and the back projection image. By adopting the method, the CT reconstruction precision and efficiency can be improved.

Description

Technical Field

[0001] The present application relates to the field of medical image processing technology, and in particular to a spiral CT reconstruction method, apparatus, computer equipment, storage medium and computer program product. Background Art

[0002] In the current spiral CT field, the reconstruction algorithm mainly adopts the WFBP (Weighted Filtered BackProjection) algorithm, which is improved on the basis of the traditional Filtered BackProjection (FBP) algorithm.

[0003] FBP reconstructs images by back-projecting projection data onto the image plane, but usually produces artifacts and noise. The WFBP algorithm reduces noise and artifacts by applying weighted filtering, and the weighting can be adjusted according to different angles or positions to adapt to image characteristics, thereby improving image quality.

[0004] However, the WFBP algorithm often converts the spiral scan into an axial scan during the orthographic projection process. This method requires the original data to be orthographically projected once, resulting in a large amount of CT reconstruction calculations and low reconstruction efficiency. Summary of the invention

[0005] Based on this, it is necessary to provide a spiral CT reconstruction method, apparatus, computer equipment, computer-readable storage medium and computer program product that can improve the efficiency and accuracy of CT reconstruction in order to address the above technical problems.

[0006] In a first aspect, the present application provides a spiral CT reconstruction method, comprising:

[0007] Performing filtering on the projection data to obtain filtered data;

[0008] Performing back-projection processing on the filtered data to obtain a back-projection image, and extracting bone tissue data from the back-projection image;

[0009] Performing orthographic projection on the bone tissue data to obtain bone projection data;

[0010] Taking the difference between the projection data and the bone projection data as bone removal tissue data;

[0011] A CT reconstructed image is obtained according to the bone removal tissue data and the back-projection image.

[0012] In one embodiment, extracting bone tissue data from the back-projection image includes:

[0013] Extracting bone tissue pixel points in the back-projection image, wherein the CT values ​​corresponding to the bone tissue pixel points are within the range of bone tissue CT values;

[0014] For any bone tissue pixel point, determine the bone tissue segment corresponding to the bone tissue pixel point according to the CT value of the bone tissue pixel point, and determine the normalization coefficient of the bone tissue pixel point according to the upper limit CT value and the lower limit CT value of the bone tissue segment, the normalization coefficients corresponding to the upper limit CT value and the lower limit CT value, and the CT value of the bone tissue pixel point, wherein the bone tissue CT value range is composed of a plurality of the bone tissue segments;

[0015] Bone tissue data is obtained according to the normalization coefficient and CT value of each bone tissue pixel point.

[0016] In one embodiment, the step of performing orthographic projection on the bone tissue data to obtain bone projection data includes:

[0017] According to each pixel point in the axial image corresponding to the ray passing through the bone tissue data and the ray length weight corresponding to each pixel point, the CT value of each pixel point is weightedly superimposed to obtain the bone projection data;

[0018] The ray length weight is negatively correlated with the length of the ray in the voxel.

[0019] In one embodiment, obtaining a CT reconstructed image according to the bone removal tissue data and the back-projection image includes:

[0020] The ratio of the bone-removed tissue data to the water reference is used as the equivalent length of each projection pixel point in the bone-removed tissue data relative to water;

[0021] Performing weighted correction on the bone projection data using the artifact correction coefficient corresponding to the equivalent length to obtain corrected bone projection data;

[0022] After the corrected bone projection data is reconstructed, the reconstruction result is superimposed on the back-projection image to obtain a CT reconstructed image.

[0023] In one embodiment, filtering the projection data to obtain filtered data includes:

[0024] The projection data are rearranged using a pre-established initial virtual coordinate system, wherein the initial virtual coordinate system takes the detector irradiated by the central radiation source as the origin, the spiral line as the horizontal axis, and the horizontal displacement bed direction as the vertical axis;

[0025] The rearranged projection data is filtered according to the windowing filter kernel to obtain filtered data.

[0026] In one embodiment, the rearranging the projection data using a pre-established initial virtual coordinate system includes:

[0027] Rearranging the projection data of the fan beam into the projection data of the parallel beam by rearranging the direction angle of the projection data;

[0028] According to the up-sampled data amount and the interval of the direction angle rearrangement, the projection data of the parallel beam is subjected to three-fold Hermite interpolation to obtain a three-dimensional surface, and the three-dimensional surface is subjected to coordinate interpolation processing perpendicular to the initial virtual coordinate system to obtain rearranged projection data.

[0029] In one embodiment, the coordinate interpolation processing of the three-dimensional surface perpendicular to the initial virtual coordinate system to obtain rearranged projection data includes:

[0030] Determine the grid points in the three-dimensional surface and calculate the rotation matrix of each of the grid points;

[0031] The rotation matrix is ​​used to transform the coordinates of each of the grid points from the initial virtual coordinate system to a target virtual coordinate system, wherein the horizontal axis of the target virtual coordinate system is horizontal;

[0032] Interpolation processing is performed on each of the grid points in the target virtual coordinate system to obtain rearranged projection data.

[0033] In one of the embodiments, the interpolation processing is performed on each of the grid points in the target virtual coordinate system to obtain rearranged projection data, including:

[0034] Performing interpolation processing on each of the grid points in the target virtual coordinate system to obtain an interpolation result;

[0035] The interpolation result is filtered according to the initial virtual coordinate system to obtain rearranged projection data.

[0036] In one embodiment, the windowing filter includes an all-pass frequency and a cutoff frequency, and the windowing filter kernel is:

[0037]

[0038] Wherein, the w represents the frequency component of the signal, the wl represents the all-pass frequency, and the wh represents the cut-off frequency.

[0039] In one embodiment, before filtering the projection data to obtain filtered data, the method further includes:

[0040] performing slope correction on the projection data;

[0041] The step of superimposing the reconstruction result with the back-projection image to obtain a CT reconstructed image includes:

[0042] The reconstruction result is superimposed with the back-projection image to obtain a superimposed result;

[0043] The stacking result is intercept-corrected, and the correction result is slice-thickness processed to obtain a CT reconstructed image.

[0044] In one embodiment, performing back-projection processing on the filtered data to obtain a back-projection image includes:

[0045] Determine the extension range according to the extension collimation width, the pitch and the reconstruction interval, and shift the starting position of the horizontal motor to the left and right based on the extension range to obtain a virtual horizontal motor position;

[0046] Interpolation processing is performed based on the virtual horizontal motor position and the filtered data to obtain extended projection data, and weighted back projection processing is performed based on the extended projection data and the filtered data to obtain a back projection image.

[0047] In one embodiment, performing weighted back projection processing based on the extended projection data and the filtered data to obtain a back projection image includes:

[0048] Performing weighted back-projection processing based on the extended projection data and the filtered data to obtain an initial back-projection image;

[0049] For any pixel point in the initial back-projection image, calculate the distance between the coordinate value of the pixel point on the depth axis and the position of the horizontal motor corresponding to the central ray passing through the pixel point, and determine that the pixel point is valid when the distance is less than a preset threshold value;

[0050] According to the Q value weight curve, the rearranged data corresponding to the effective pixel points are processed by weighted summation to obtain the back-projection image.

[0051] In a second aspect, the present application also provides a CT reconstruction device, comprising:

[0052] A filtering module, used for filtering the projection data to obtain filtered data;

[0053] A back-projection module, used for performing back-projection processing on the filtered data to obtain a back-projection image, and extracting bone tissue data from the back-projection image;

[0054] An orthographic projection module, used for performing orthographic projection on the bone tissue data to obtain bone projection data;

[0055] A processing module, used for taking the difference between the projection data and the bone projection data as bone removal tissue data;

[0056] A reconstruction module is used to obtain a CT reconstructed image based on the bone removal tissue data and the back-projection image.

[0057] In one embodiment, the back projection module is specifically used for:

[0058] Extracting bone tissue pixel points in the back-projection image, wherein the CT values ​​corresponding to the bone tissue pixel points are within the range of bone tissue CT values;

[0059] For any bone tissue pixel point, determine the bone tissue segment corresponding to the bone tissue pixel point according to the CT value of the bone tissue pixel point, and determine the normalization coefficient of the bone tissue pixel point according to the upper limit CT value and the lower limit CT value of the bone tissue segment, the normalization coefficients corresponding to the upper limit CT value and the lower limit CT value, and the CT value of the bone tissue pixel point, wherein the bone tissue CT value range is composed of a plurality of the bone tissue segments;

[0060] Bone tissue data is obtained according to the normalization coefficient and CT value of each bone tissue pixel point.

[0061] In one embodiment, the orthographic projection module is specifically used for:

[0062] According to each pixel point in the axial image corresponding to the ray passing through the bone tissue data and the ray length weight corresponding to each pixel point, the CT value of each pixel point is weightedly superimposed to obtain the bone projection data;

[0063] The ray length weight is negatively correlated with the length of the ray in the voxel.

[0064] In one embodiment, the reconstruction module is specifically used to:

[0065] The ratio of the bone-removed tissue data to the water reference is used as the equivalent length of each projection pixel point in the bone-removed tissue data relative to water;

[0066] Performing weighted correction on the bone projection data using the artifact correction coefficient corresponding to the equivalent length to obtain corrected bone projection data;

[0067] After the corrected bone projection data is reconstructed, the reconstruction result is superimposed on the back-projection image to obtain a CT reconstructed image.

[0068] In one embodiment, the filtering module is specifically used for:

[0069] The projection data are rearranged using a pre-established initial virtual coordinate system, wherein the initial virtual coordinate system takes the detector irradiated by the central radiation source as the origin, the spiral line as the horizontal axis, and the horizontal displacement bed direction as the vertical axis;

[0070] The rearranged projection data is filtered according to the windowing filter kernel to obtain filtered data.

[0071] In one embodiment, the filtering module is further used for:

[0072] Rearranging the projection data of the fan beam into the projection data of the parallel beam by rearranging the direction angle of the projection data;

[0073] According to the up-sampled data amount and the interval of the direction angle rearrangement, the projection data of the parallel beam is subjected to three-fold Hermite interpolation to obtain a three-dimensional surface, and the three-dimensional surface is subjected to coordinate interpolation processing perpendicular to the initial virtual coordinate system to obtain rearranged projection data.

[0074] In one embodiment, the filtering module is further used for:

[0075] Determine the grid points in the three-dimensional surface and calculate the rotation matrix of each of the grid points;

[0076] The rotation matrix is ​​used to transform the coordinates of each of the grid points from the initial virtual coordinate system to a target virtual coordinate system, wherein the horizontal axis of the target virtual coordinate system is horizontal;

[0077] Interpolation processing is performed on each of the grid points in the target virtual coordinate system to obtain rearranged projection data.

[0078] In one embodiment, the filtering module is further used for:

[0079] Performing interpolation processing on each of the grid points in the target virtual coordinate system to obtain an interpolation result;

[0080] The interpolation result is filtered according to the initial virtual coordinate system to obtain rearranged projection data.

[0081] In one embodiment, the windowing filter includes an all-pass frequency and a cutoff frequency, and the windowing filter kernel is:

[0082]

[0083] Wherein, the w represents the frequency component of the signal, the wl represents the all-pass frequency, and the wh represents the cut-off frequency.

[0084] In one embodiment, the device further comprises:

[0085] A correction module, used for performing slope correction on the projection data;

[0086] The reconstruction module is further specifically used for:

[0087] The reconstruction result is superimposed with the back-projection image to obtain a superimposed result;

[0088] The stacking result is intercept-corrected, and the correction result is slice-thickness processed to obtain a CT reconstructed image.

[0089] In one embodiment, the back projection module is specifically used for:

[0090] Determine the extension range according to the extension collimation width, the pitch and the reconstruction interval, and shift the starting position of the horizontal motor to the left and right based on the extension range to obtain a virtual horizontal motor position;

[0091] Interpolation processing is performed based on the virtual horizontal motor position and the filtered data to obtain extended projection data, and weighted back projection processing is performed based on the extended projection data and the filtered data to obtain a back projection image.

[0092] In one embodiment, the back projection module is further used for:

[0093] Performing weighted back-projection processing based on the extended projection data and the filtered data to obtain an initial back-projection image;

[0094] For any pixel point in the initial back-projection image, calculate the distance between the coordinate value of the pixel point on the depth axis and the position of the horizontal motor corresponding to the central ray passing through the pixel point, and determine that the pixel point is valid when the distance is less than a preset threshold value;

[0095] According to the Q value weight curve, the rearranged data corresponding to the effective pixel points are processed by weighted summation to obtain the back-projection image.

[0096] In a third aspect, the present application further provides a computer device, comprising a memory and a processor, wherein the memory stores a computer program, and the processor implements any of the above spiral CT reconstruction methods when executing the computer program.

[0097] In a fourth aspect, the present application further provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program implements any of the above spiral CT reconstruction methods when executed by a processor.

[0098] In a fifth aspect, the present application also provides a computer program product, including a computer program, which implements any of the above spiral CT reconstruction methods when executed by a processor.

[0099] The above-mentioned spiral CT reconstruction method, device, computer equipment, storage medium and computer program product, the spiral CT reconstruction method provided in the embodiment of the present application, can filter the projection data to obtain filtered data, and perform back projection processing on the filtered data, after obtaining the back projection image, extract the bone tissue data from the back projection image. Perform forward projection on the bone tissue data to obtain bone projection data, and use the difference between the projection data and the bone projection data as bone-removed tissue data, and then obtain a CT reconstructed image based on the bone-removed tissue data and the back projection image. Using the spiral CT reconstruction method, device, computer equipment, storage medium and computer program product provided in the embodiment of the present application, in the post-reconstruction process, there is no need to perform a forward projection on the original data, and the reconstructed projection data can be directly used, which can effectively reduce the amount of calculation and improve the efficiency and accuracy of CT reconstruction. BRIEF DESCRIPTION OF THE DRAWINGS

[0100] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the related technologies, the drawings required for use in the embodiments or the related technical descriptions are briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0101] Figure 1 is a schematic flow chart of a spiral CT reconstruction method in one embodiment;

[0102] Figure 2 is a schematic diagram of a positive projection path in one embodiment;

[0103] Figure 3 is a flow chart of step 104 in one embodiment;

[0104] Figure 4a The relationship between the ray and the axial image pixels in the orthographic projection is shown;

[0105] Figure 4b The relationship between the ray and the voxel diagonal within a voxel of the target reconstructed object is shown;

[0106] Figure 5 is a flow chart of step 108 in one embodiment;

[0107] Figure 6 is a flow chart of step 102 in one embodiment;

[0108] Figure 7 is a schematic diagram of a spiral CT reconstruction method in one embodiment;

[0109] Figure 8is a flow chart of step 602 in one embodiment;

[0110] Fig. 9 is a schematic diagram of a spiral CT reconstruction method in one embodiment;

[0111] Fig.10 is a schematic diagram of channel spacing in one embodiment;

[0112] Fig.11a is a schematic diagram of a three-dimensional curved surface in one embodiment;

[0113] Fig.11b is a schematic diagram of a UV direction plane in one embodiment;

[0114] Fig.12 is a schematic diagram of a spiral CT reconstruction method in one embodiment;

[0115] Fig.13 is a schematic diagram of step 702 in one embodiment;

[0116] Fig.14 is a schematic diagram of step 1206 in one embodiment;

[0117] Fig.15 is a schematic diagram of a spiral CT reconstruction method in one embodiment;

[0118] Fig.16 is a flow chart of step 506 in one embodiment;

[0119] Fig.17a A schematic diagram of a spiral CT reconstruction method in one embodiment;

[0120] Fig.17b is another schematic diagram of a spiral CT reconstruction method in one embodiment;

[0121] Fig.18 is a flow chart of step 104 in one embodiment;

[0122] Fig.19 is a schematic diagram of a virtual horizontal motor position in one embodiment;

[0123] Fig. 20 is a flow chart of step 1704 in one embodiment;

[0124] Fig.21 is a schematic diagram of a Q value weight curve in one embodiment;

[0125] Fig.22a A schematic diagram of a spiral CT reconstruction method in one embodiment;

[0126] Figure 22bis another schematic diagram of a spiral CT reconstruction method in one embodiment;

[0127] Fig.23 is a structural block diagram of a spiral CT reconstruction device in one embodiment;

[0128] Fig.24 is a schematic diagram of a spiral CT reconstruction method in one embodiment;

[0129] Fig.25 is a schematic diagram of a spiral CT reconstruction method in one embodiment;

[0130] Fig.26 FIG. 4 is a diagram showing the internal structure of a computer device in one embodiment. DETAILED DESCRIPTION

[0131] In order to make the purpose, technical solution and advantages of the present application more clearly understood, the present application is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0132] In one embodiment, Figure 1 As shown, a spiral CT reconstruction method is provided. This embodiment uses the method applied to a terminal as an example for illustration. It can be understood that the method can also be applied to a server, and can also be applied to a system including a terminal and a server, and is implemented through the interaction between the terminal and the server. In this embodiment, the method includes the following steps 102 to 110, wherein:

[0133] Step 102: filter the projection data to obtain filtered data.

[0134] In the embodiment of the present application, during the CT reconstruction process, the projection data may be filtered to remove noise and improve the CT reconstruction effect. For example, the projection data may be filtered using filtering algorithms such as Laplace filtering and Hamming window filtering.

[0135] Exemplarily, projection data can be obtained from a CT detection instrument. These projection data are the intensity information of X-rays after penetrating an object at multiple angles. First, a filter is selected, for example, a Hamming window filter is selected. The Hamming window filter can smooth the data in the frequency domain and reduce artifacts. Then the projection data can be Fourier transformed to convert the projection data to the frequency domain. The Fourier transform can decompose the signal into different frequency components to facilitate filtering. The aforementioned selected filter is applied in the frequency domain, including: multiplying the Fourier transform result of the projection data with the frequency response of the filter, thereby suppressing unnecessary frequency components (such as noise and artifacts), selectively retaining low-frequency components, and removing high-frequency components. Finally, the filtered frequency domain data is inversely Fourier transformed and converted back to the time domain (or space domain), so that the filtered data (filtered projection data) can be obtained. Exemplarily, after obtaining the filtered data, the filtered data can be further processed, such as normalization, smoothing, etc., to improve image quality and usability.

[0136] Through the above filtering process, not only the noise and artifacts are removed, but also the contrast and edge features of the image can be improved, providing better input data for subsequent CT image reconstruction.

[0137] Step 104, back-projecting the filtered data to obtain a back-projected image, and extracting bone tissue data from the back-projected image.

[0138] Backprojection is the process of backprojecting each projection data into the image space according to its imaging angle. The specific operation is to backproject each projection data into a two-dimensional image along the angle at which it was acquired. Exemplarily, for the filtered data at any angle, each pixel value of the filtered data can be accumulated along the angle to gradually construct a backprojection image. The filtered data at all angles are accumulated during the backprojection process to form a backprojection image, which usually contains mixed information of various tissues. In fact, weighted backprojection methods (such as normalizing the filter values) can also be used during the accumulation process to improve the quality of backprojection.

[0139] After obtaining the back-projection data, clear bone tissue data can be extracted from the back-projection image. For example, the bone tissue data can be extracted from the back-projection image through the CT value range corresponding to the bone tissue. The bone tissue data extraction process will be introduced in the following embodiments of the present application, and no further description will be given in the embodiments of the present application.

[0140] Step 106, performing orthographic projection on the bone tissue data to obtain bone projection data.

[0141] In the embodiment of the present application, when the bone tissue data is orthographically projected, the path of the orthographic projection can be regarded as a ray following the scanning path at a conical angle (that is, the whole follows a spiral line). Figure 2 As shown, Figure 2 The figure shows the forward projection path of a frame of image. The scanning path follows a spiral line as a whole, which simulates the process of imaging an image in three-dimensional space. This scanning method can ensure that each frame of the original image has corresponding projection data, thus eliminating the rearrangement step in the reconstruction process. Eliminating the rearrangement step can significantly increase the speed of data processing, making the imaging process more efficient.

[0142] Moreover, the projection line can also be downsampled during the orthographic projection process, which can reduce the amount of data processing, for example, by reducing the resolution of the image to half (or one quarter) of the original image, which can reduce the amount of calculation and improve the efficiency of the orthographic projection while ensuring the imaging quality. When downsampling to extract the original image value, in the embodiment of the present application, the x-coordinate can also be simply moved one position to the left, so that the corresponding projection value can be quickly obtained, and downsampling can be achieved without complicated processing.

[0143] In an exemplary embodiment, in step 106, performing orthographic projection on the bone tissue data to obtain bone projection data may include the following steps:

[0144] According to each pixel point in the axial image corresponding to the data of the ray passing through the bone tissue and the ray length weight corresponding to each pixel point, the CT value of each pixel point is weightedly superimposed to obtain the bone projection data; wherein the ray length weight is negatively correlated with the length of the ray in the voxel.

[0145] In the embodiment of the present application, the ray length weight corresponding to the ray can be determined according to the length of the ray in the voxel, wherein the ray length weight is inversely proportional to the length of the ray in the voxel, that is, the longer the ray in the voxel, the smaller the influence of the point on the axial plane (the orthogonal weight is the largest). Determine the pixel points in the axial image corresponding to the bone tissue data through which the ray passes, and perform weighted superposition of the pixel value corresponding to each pixel point and the ray length weight, so as to obtain the orthogonal projection value of each ray in the orthogonal projection, and then obtain the bone projection data. Exemplary, Figure 4a The relationship between the ray and the axial image pixels in the orthographic projection is shown; Figure 4b The relationship between the ray and the voxel diagonal in a voxel of the target reconstructed object is shown.

[0146] Step 108: taking the difference between the projection data and the bone projection data as bone removal tissue data.

[0147] In an embodiment of the present application, after obtaining the bone projection data, the difference between the projection data and the bone projection data can be determined. This difference represents the signal of other tissues (such as soft tissue or liquid) other than bone tissue at the pixel position. The difference is used as bone-removed tissue data, and the bone-removed tissue data is an image of tissues other than bone tissue.

[0148] Step 110, obtaining a CT reconstructed image based on the bone removal tissue data and the back-projection image.

[0149] In the embodiment of the present application, the bone projection data can be corrected by removing the bone tissue data, and the corrected bone projection data can be used for reconstruction to obtain a new bone reconstruction image. After the new bone reconstruction image is superimposed with the back-projection image, a CT reconstruction image can be obtained.

[0150] The spiral CT reconstruction method provided in the embodiment of the present application can filter the projection data to obtain the filtered data, and perform back-projection processing on the filtered data. After obtaining the back-projection image, the bone tissue data can be extracted from the back-projection image. The bone tissue data is forward projected to obtain the bone projection data, and the difference between the projection data and the bone projection data is used as the bone-removed tissue data, and then the CT reconstructed image is obtained according to the bone-removed tissue data and the back-projection image. Using the spiral CT reconstruction method provided in the embodiment of the present application, in the post-reconstruction bone sclerosis correction process, there is no need to perform a forward projection on the original data, and the reconstructed projection data can be directly used, which can reduce the amount of calculation and improve the efficiency of CT reconstruction.

[0151] In an exemplary embodiment, referring to Figure 3 As shown, in step 104, extracting bone tissue data from the back-projection image may include the following steps 302 to 306, wherein:

[0152] Step 302, extracting bone tissue pixels in the back-projection image, wherein the CT values ​​corresponding to the bone tissue pixels are within the range of bone tissue CT values;

[0153] Step 304, for any bone tissue pixel point, determine the bone tissue segment corresponding to the bone tissue pixel point according to the CT value of the bone tissue pixel point, and determine the normalization coefficient of the bone tissue pixel point according to the upper and lower CT values ​​of the bone tissue segment, the normalization coefficients corresponding to the upper and lower CT values, and the CT value of the bone tissue pixel point, wherein the range of the bone tissue CT value is composed of multiple bone tissue segments;

[0154] Step 306, obtaining bone tissue data according to the normalization coefficient and CT value of each bone tissue pixel point.

[0155] In the embodiment of the present application, a bone tissue CT value range is preset, and the bone tissue CT value range can be composed of multiple bone tissue segments, for example: the bone tissue CT value range is hu1-hu2, including 3 bone tissue segments, namely: the first segment hu1-hua (such as 850-1150), which is a CT value between brain tissue and bone tissue, the second segment hua-hub (such as 1150-2075), which is a CT value range of bone tissue, and the third segment hub-hu2 (such as 2075-3000), which is a CT value between bone tissue and high-density body. The upper limit CT value and the lower limit CT value of each segment have a preset normalization coefficient, for example: hu1 corresponds to s1 (such as 0), hua corresponds to s2 (such as 0.4), hub corresponds to s3 (such as 0.7), and hu2 corresponds to s4 (such as 1).

[0156] For each pixel in the back-projection image, the pixel whose corresponding CT value is within the range of bone tissue CT values ​​can be regarded as a bone tissue pixel. For example, when the CT value hui of pixel i is between hu1 and hu2, the pixel i can be determined as a bone tissue pixel.

[0157] After obtaining each bone tissue pixel point, the CT value of each bone tissue pixel point can be obtained, and the bone tissue segment corresponding to each bone tissue pixel point can be determined based on the CT value of each bone tissue pixel point, and the normalization coefficient of each bone tissue pixel point can be determined based on the upper and lower CT values ​​of the bone tissue segment corresponding to each bone tissue pixel point, the normalization coefficient corresponding to the upper and lower CT values, and the CT value of the bone tissue pixel point. Exemplarily, the determination process of the normalization coefficient can refer to the following formula (I):

[0158]

[0159] Among them, s represents the normalization coefficient corresponding to the bone tissue pixel, s i is the normalization coefficient corresponding to the lower limit CT value of the bone tissue segment where the bone tissue pixel is located, s j is the normalization coefficient corresponding to the upper limit CT value of the bone tissue segment where the bone tissue pixel point is located, hu is the CT value of the bone tissue image point, and hu i is the lower limit CT value of the bone tissue segment where the bone tissue image point is located, hu j It is the upper limit CT value of the bone tissue segment where the bone tissue image point is located.

[0160] After obtaining the normalization coefficient corresponding to the bone tissue pixel point, the CT value of each bone tissue pixel point is multiplied by the normalization coefficient corresponding to each bone tissue pixel point, which is the bone tissue data extracted by segmentation.

[0161] By adopting the spiral CT reconstruction method provided in the embodiment of the present application and setting multiple bone tissue segments, the accuracy of bone tissue data extraction can be improved, and then bone sclerosis correction can be performed based on the extracted bone tissue data, which can further improve the CT reconstruction accuracy.

[0162] In an exemplary embodiment, referring to Figure 5 As shown, in step 108, obtaining a CT reconstructed image based on the bone removal tissue data and the back-projection image may include the following steps 502 to 506, wherein:

[0163] Step 502, taking the ratio of the bone-removed tissue data to the water reference as the equivalent length of each projection pixel point in the bone-removed tissue data relative to water;

[0164] Step 504, performing weighted correction on the bone projection data using the artifact correction coefficient corresponding to the equivalent length to obtain corrected bone projection data;

[0165] Step 506: After reconstructing the corrected bone projection data, the reconstruction result is superimposed on the back-projection image to obtain a CT reconstructed image.

[0166] In the embodiment of the present application, the bone-removed tissue data is first divided by the water reference, and the ratio is used as the equivalent length of each projection pixel in the bone-removed tissue data relative to water. Then, based on the equivalent length, the corresponding artifact correction coefficient is searched in a pre-set coefficient table, and the artifact correction coefficient is used to perform polynomial weighted correction on the above-mentioned bone projection data to obtain the corrected bone projection data. The correction process can be shown in the following formula (II):

[0167]

[0168] Where prj represents the bone projection data before correction, prj′ represents the bone projection data after correction, N is the number of polynomials, α i is the artifact correction coefficient.

[0169] After obtaining the corrected bone projection data, reconstruction can be performed based on the corrected bone projection data to obtain a reconstruction result. The specific reconstruction process is subject to the selected reconstruction algorithm. The reconstruction result is superimposed with the back-projection image to obtain a CT reconstructed image.

[0170] Since the scanning path of the bone projection data in the embodiment of the present application is a rearranged path, there is no need to perform rearrangement, so the efficiency of bone sclerosis correction can be greatly improved.

[0171] In an exemplary embodiment, referring to Figure 6 As shown, in step 102, filtering the projection data to obtain filtered data may include the following steps 602 to 604, wherein:

[0172] Step 602, rearrange the projection data using a pre-established initial virtual coordinate system, wherein the initial virtual coordinate system takes the detector irradiated by the central radiation source as the origin, the spiral line as the horizontal axis, and the horizontal displacement bed direction as the vertical axis;

[0173] Step 604 , filtering the rearranged projection data according to the windowing filter core to obtain filtered data.

[0174] In the embodiment of the present application, the detector position irradiated by the central radiation source in the radiation source irradiation range (cone beam) can be firstly used as the origin of the UV coordinate system (initial virtual coordinate system), and then the horizontal axis (U axis) of the coordinate system is defined as the direction along the spiral line, which is usually a rotation trajectory around the scanning center and is used to represent the movement of the detector during the scanning process; the vertical axis (V axis) of the coordinate system is defined as the direction parallel to the horizontal displacement bed direction, which is used to represent the direction in which the patient moves along the bed during the scanning process, with reference to Figure 7 shown.

[0175] Exemplarily, a two-dimensional UV coordinate system is established with the center point of the detector (the point corresponding to the source is the origin. It should be noted that if the detector is not precisely placed at the center, then this center point will not coincide with the point corresponding to the source). According to the preset sampling size, △U (the amount of change on the U axis) and △V (the amount of change on the V axis) in the coordinate system are calculated. Among them, the U axis is called a virtual coordinate channel, and its characteristic is that it is not a straight line, but a curved shape, which means that in this coordinate system, the value of the U axis is not arranged according to a conventional linear relationship.

[0176] In each frame, the position of the detector irradiated by the central ray of the source irradiation range (cone beam) is taken as the origin, the spiral line is taken as the U axis, and the same direction as the horizontal displacement bed is taken as the V axis. Then the UV wrapped area is an irregular image, see Figure 7 In the midpoint matrix, the values ​​are taken at equal intervals according to the wrapping distance in the V direction, and at equal curve distances in the spiral direction in the U direction.

[0177] The original projection data is obtained from the CT detector. These data are usually recorded at different angles and positions. Each projection data point usually includes information such as the projection angle, the detector position, and the corresponding ray intensity. The obtained projection data is rearranged using a pre-established UV coordinate system. For example, the transformation parameters of the projection data can be calculated, and the position of each projection data point in the new coordinate system can be calculated based on the definition of the initial virtual coordinate system and the transformation parameters of the projection data. Each projection data point is transformed from the original coordinate system to the new virtual UV coordinate system. This process often needs to consider information such as the angle, position, and spiral trajectory of the detector.

[0178] In an exemplary embodiment, referring to Figure 8 As shown, in step 602, the projection data is rearranged using a pre-established initial virtual coordinate system, which may include steps 802 to 804, wherein:

[0179] Step 802, rearranging the projection data of the fan beam into the projection data of the parallel beam by rearranging the direction angle of the projection data;

[0180] Step 804, according to the up-sampled data amount and the interval of the direction angle rearrangement, the projection data of the parallel beam is subjected to three-dimensional Hermite interpolation to obtain a three-dimensional surface, and the three-dimensional surface is subjected to coordinate interpolation processing perpendicular to the initial virtual coordinate system to obtain rearranged projection data.

[0181] In the embodiment of the present application, the projection data is first rearranged in azimuth to rearrange the projection data of the fan beam into a parallel beam: the non-central rays of the fan beam are rearranged into a parallel beam pattern so that these rays are parallel to the central rays under the azimuth. Exemplarily, the azimuth rearrangement may include azimuth adjustment and geometric transformation, wherein the azimuth adjustment includes adjusting the fan beam ray data under different azimuths to be consistent with the target parallel beam direction; the geometric transformation includes geometrically transforming the data to recalculate and adjust each ray in the fan beam to the form of a parallel beam, and the process may include rotation and scaling to match the spacing of the parallel beam.

[0182] For example, refer to Fig. 9 As shown, the direction angle is first rearranged to rearrange the fan-shaped beam into a parallel beam. Fig. 9 The rays parallel to the central ray at the middle direction angle (set as the starting angle 0 degrees) are rearranged to the non-central ray at this direction angle, that is, the non-central ray is parallel to the central ray, and the interval after rearrangement is smaller than the original interval, because the non-central position is at any position, such as the oblique ray hitting the detector at the coordinate (iu, iv), and is rotated through the rotation center to the ray parallel to the central ray at the direction angle IV ( Figure 8 The example in the figure is rotated 20.35 degrees) and is closer to the center: AI <DJ。

[0183] It should be noted that the channel spacing is not equal after the azimuth angle is rearranged. The channel spacing is related to the installation rules of the detector. For example, assuming that the number of channels is 384, refer to Fig.10 As shown, the intervals between channels after the 384 channel interval direction angles are rearranged (the horizontal axis is the channel interval, the vertical axis is the interval distance in mm, from left to right are the interval between the 1st channel and the 0th channel, the interval between the 2nd channel and the 1st channel, ...).

[0184] After the azimuth angles are rearranged, the intervals between the rearranged channels (hereinafter referred to as the azimuth angle rearrangement intervals in the present application embodiment) can be obtained. According to the upsampled data volume and the azimuth angle rearrangement intervals, the projection data of the parallel beams are interpolated three times with Hermite (Hermite interpolation) to obtain a three-dimensional surface. Fig.11a Since the rearranged single-view detector data is not a standard two-dimensional matrix, but a curved matrix related to the spiral track, the cubic interpolation used in the embodiment of the present application has a better interpolation effect than the general linear interpolation.

[0185] Since the interpolation results in a three-dimensional surface, it is similar to a parallelogram in two dimensions. That is, after the Hermite interpolation in the previous step, the V axis has the same spacing upward, but because the U axis is not horizontal, directly using the data on the three-dimensional surface during back projection will result in inappropriate data for interpolation, especially when the pitch is low. For example, refer to Fig.12 As shown (it should be noted that Fig.12 What is actually displayed in the figure is a curved surface rather than a plane). The points used for bilinear interpolation of point G1 in the virtual UV coordinate system are G, M, L, and H. It can be seen from the figure that point H is obviously inappropriate, and the most appropriate data points should be G, M, L, and R.

[0186] Therefore, after performing the three-dimensional Hermite interpolation, the three-dimensional surface obtained after the interpolation can be subjected to a coordinate interpolation process perpendicular to the initial virtual coordinate system. The purpose of the interpolation process is to make the U axis horizontal.

[0187] In an exemplary embodiment, referring to Fig.13 As shown, in step 802, coordinate interpolation processing of the three-dimensional surface perpendicular to the initial virtual coordinate system is performed to obtain rearranged projection data, which may include steps 1302 to 1306, wherein:

[0188] Step 1302, determining the grid points in the three-dimensional surface and calculating the rotation matrix of each grid point;

[0189] Step 1304, using a rotation matrix to transform the coordinates of each grid point from the initial virtual coordinate system to the target virtual coordinate system, where the horizontal axis of the target virtual coordinate system is horizontal;

[0190] Step 1306: interpolate each grid point in the target virtual coordinate system to obtain rearranged projection data.

[0191] In the embodiment of the present application, after obtaining the three-dimensional surface, the UV coordinates need to be interpolated to ensure that the U axis is horizontal. Exemplarily, evenly distributed grid points are selected on the three-dimensional surface according to the needs. Usually, they can be distributed according to the geometric features of the object to ensure that the entire area within the effective range is covered. Exemplarily, the parameters u and v can be set in three-dimensional space to form a UV grid, and these grid points will be interpolated in the subsequent process.

[0192] The target direction of the U axis can be determined, and the rotation matrix R can be constructed by calculating the rotation angle θ based on the initial UV coordinates and the target direction of the grid point. For example:

[0193] Furthermore, this rotation matrix is ​​applied to the UV coordinates of each grid point to adjust the direction of the U axis, for example: (u,,v,,1) = R(u,v,1), where (u,v) is the original coordinate in the standard initial virtual coordinate system, and (u,,v,) represents the new coordinate after the original coordinate is converted to the target virtual coordinate system. In this way, the coordinates of each grid point can be transformed from the virtual coordinate system to the virtual coordinate system horizontal to the U axis.

[0194] After the coordinate transformation, each grid point can be interpolated in the target virtual coordinate system. For example, bilinear interpolation or bicubic interpolation can be used to obtain smoother projection data.

[0195] In an exemplary embodiment, referring to Fig.14 As shown, in step 1306, interpolation processing is performed on each grid point in the target virtual coordinate system to obtain rearranged projection data, which may include the following steps 1402 to 1404, wherein:

[0196] Step 1402, performing interpolation processing on each grid point in the target virtual coordinate system to obtain an interpolation result;

[0197] Step 1404 , filtering the interpolation result according to the initial virtual coordinate system to obtain rearranged projection data.

[0198] In the embodiment of the present application, each grid point can be interpolated in the target virtual coordinate system to obtain an interpolation result. The interpolation method can include methods such as bilinear interpolation or bicubic interpolation. The interpolation method is not specifically limited in the embodiment of the present application. After obtaining the interpolation result, the data within the boundary range of the original coordinate system can be saved as valid data, and the peripheral data is non-valid data. Exemplarily, the boundary range of the three-dimensional surface in the initial virtual coordinate system can be determined, and the data that exceeds the boundary range in the interpolation result can be marked as invalid data to avoid affecting the result during the reconstruction process. Only the data within the boundary range is extracted as valid data to form the final rearranged projection data for subsequent CT reconstruction.

[0199] Since the initial virtual coordinate system provided in the embodiment of the present application has a spiral line as the horizontal axis and a horizontal displacement bed direction as the vertical axis, which is consistent with the scanning path, curve fitting interpolation is performed according to the U-axis of the spiral line during rearrangement, and grid node interpolation is performed according to the UV coordinate system that approximates the horizontal U-axis after rearrangement, which can greatly improve the accuracy of back projection.

[0200] In an exemplary embodiment, the windowing filter includes an all-pass frequency and a cutoff frequency, and the windowing filter kernel Π(w) is shown in the following formula (III):

[0201]

[0202] Among them, w represents the frequency component of the signal, wl represents the all-pass frequency, and wh represents the cutoff frequency.

[0203] In one example, the window function can be a sinc window, which has clearer details at the same cutoff frequency than the hanning window, and the hanning is too smooth. Take two limits, low frequency (all-pass frequency) and high frequency (cutoff frequency). When the frequency component of the projection data is less than or equal to the all-pass frequency, the window value of the corresponding windowing filter kernel can be determined to be 1; or, when the frequency component of the projection data is greater than the cutoff frequency, the window value of the corresponding windowing filter kernel can be determined to be 0; or, when the frequency component of the projection data is greater than the all-pass frequency and less than the cutoff frequency, the window value of the corresponding windowing filter kernel is

[0204] In one example, referring to Fig.15 As shown, the positions of the full-pass frequency and the cutoff frequency are shown, the horizontal axis is the sampling point (positive and negative normalization, that is, the center is the center of the image sampling width), and the vertical axis is the window function value. Fig.15 Where WL is the full-pass frequency (for example, it can take values ​​of 1, 1 / 16, 1 / 32), WH is the cutoff frequency (for example, it can take values ​​of 1, 2 / 3, 1 / 2), WL<=WH; the horizontal axis is the sampling width of the convolution (the value of 0 represents the sampling center), and the vertical axis is the normalized window width value.

[0205] In one example, according to the requirements of reconstructing the image, the full-pass frequency and the cutoff frequency can take different values ​​and be divided into multiple groups, for example 7 groups, corresponding to different sharpening and smoothing effects. One group can be selected from the multiple groups to calculate the windowing filter kernel according to the needs, and the projection data can be subjected to windowing filtering processing.

[0206] In an exemplary embodiment, before filtering the projection data to obtain the filtered data, the method may further include the steps of: performing slope correction on the projection data;

[0207] In this embodiment, refer to Fig.16 As shown, in step 506, the reconstruction result is superimposed with the back-projection image to obtain a CT reconstructed image, which may include steps 1602 to 1604, wherein:

[0208] Step 1602, superimposing the reconstruction result and the back-projection image to obtain a superimposed result;

[0209] Step 1604: perform intercept correction on the superposition result and perform layer thickness processing on the correction result to obtain a CT reconstructed image.

[0210] In the related art, slope correction and intercept correction are performed in the post-reconstruction process to stretch the dynamic range of the reconstruction result so that the air reaches -1000 and the water reaches 0, that is, linear correction, and the slope and intercept are obtained to correct each pixel point. In the embodiment of the present application, the projection data is slope corrected before back projection. Exemplarily, for the projection data, the projection data can be corrected based on the air average value. The difference between the air average value under the initial angle period of the projection data and the air average value under the initial angle period of the air can be considered during the correction process. The specific process refers to the following formula (IV):

[0211]

[0212] Where Δ=I air -J air , J air Refers to the value of the initial angle of the projection data in one cycle (for example, one cycle is 360 degrees, and the projection data is at 30 degrees). air It refers to the value of the initial angle of the air in one cycle (for example, one cycle is 360 degrees, and the air data is at 70 degrees). Δ represents the deviation of the angle. j It represents the projection data before correction, air j+Δ It's prj j The corresponding air data, where j represents the current angle position, is the projection data takes the data at the angle position j, then the air takes the data at the position j+Δ, prj j' represents the projection data after correction.

[0213] After air correction, the air-corrected projection data can be further water-hardened corrected. The projection data after water-hardening correction can be slope-corrected, and the reconstruction result can be intercept-corrected in the post-reconstruction process. Since the projection value is a superposition of pixel points, the slope part of the CT value is corrected before reconstruction, which can be scaled equally, effectively expanding the dynamic range of the projection data and improving data accuracy and reconstruction accuracy.

[0214] In the post-reconstruction process, the reconstruction result and the back-projection image can be superimposed to obtain a superimposed result, and the superimposed result can be intercept-corrected to obtain a corresponding corrected result. The corrected result can be further processed by layer thickness to obtain a CT reconstructed image.

[0215] The process of layer thickness processing includes: layer thickness superposition calculation, obtaining effective layer thickness according to reconstruction interval, display layer spacing and preset layer thickness, layer thickness superposition calculation, and reference Fig.17a and Fig.17b As shown, schematic diagrams are shown without layer thickness processing and with layer thickness processing respectively.

[0216] In the present application, refer to Fig.18 As shown, in step 104, back-projection processing is performed on the filtered data to obtain a back-projection image, which may include the following steps 1802 to 1804, wherein:

[0217] Step 1802, determining an extension range according to the extension collimation width, the pitch and the reconstruction interval, and shifting the starting position of the horizontal motor to the left and right based on the extension range to obtain a virtual horizontal motor position;

[0218] Step 1804 , interpolation processing is performed based on the virtual horizontal motor position and the filtered data to obtain extended projection data, and weighted back projection processing is performed based on the extended projection data and the filtered data to obtain a back projection image.

[0219] In the embodiment of the present application, the horizontal motor is used to drive the displacement bed, so that the patient and the CT device move relative to each other in the horizontal direction. During the spiral scanning, the CT device performs a rotational motion, and driven by the horizontal motor, the patient moves relative to the CT device along the axis of the rotational motion, i.e., the horizontal direction, to realize the spiral scanning process.

[0220] During the CT reconstruction process, the projection data can also be extended in the Z direction, including setting the virtual horizontal motor position. For example, the extension range can be determined according to the extension collimation width, pitch and reconstruction interval, and then the horizontal motor starting position is shifted left and right respectively based on the extension range to obtain the virtual horizontal motor position. Fig.19As shown, the virtual horizontal motor position (that is, the position of the spiral sweep extending upward and downward) is shown. The process of determining the extension range can refer to the following formula (V):

[0221]

[0222] Among them, HarfOffset represents the extension range, ecol represents the extension collimation width, inc represents the reconstruction interval, and p represents the pitch.

[0223] In the embodiment of the present application, interpolation can be performed based on the virtual motor horizontal motor position and projection data to obtain extended projection data. For the extended projection data, CT weighted back projection is performed using the original projection data and the extended projection data, thereby increasing the number of reconstruction frames and reducing artifacts of incomplete projection.

[0224] Exemplarily, each pixel point and each rearrangement perspective of the target coordinate system XY layer (each axial diagram) is traversed to determine the corresponding starting and ending positions of the horizontal motor under the rearrangement perspective. The starting and ending positions are shifted left and right according to the previously determined extension range, and the upper and lower boundaries of the effective Z layer under the rearrangement perspective can be determined. The pixel points on XY are converted to the initial virtual coordinate system to determine the UV coordinate value, and the channel position is determined according to the U coordinate value, that is, the position of the horizontal channel for collecting rearrangement data (only the valid horizontal channel is calculated); according to the V-direction position of the ray, the V coordinate value determined by converting the pixel points on XY to the initial virtual coordinate system and the V position (i.e., Z position) of the desired image layer, the row number is determined, that is, the position of the longitudinal channel for collecting rearrangement data (only the valid longitudinal channel is calculated). According to the position of the valid longitudinal channel for collecting rearrangement data, the weight W is calculated according to the following formula (six): Q (q):

[0225]

[0226] Among them, W Q (q) is the Q value weight, q is the normalized corresponding value of the layer direction (row number), and Q is a threshold (within 1). For example, Fig.21 A schematic diagram of the Q value weight curve is shown. Assuming Q is 0.6, Fig.21 It can be seen that W Q (q) weight is the vertical axis, q is the horizontal axis, q is the normalized value according to the row direction, that is, the V direction. For example, if there are 32 rows in the V direction, then the center is the 16.5th row, and the corresponding point is the horizontal 0 point in Figure 22, which is the position of q=0. Then the 1st row and the 32nd row are the positions of q=-1 and q-1.

[0227] Iterate the rearranged data corresponding to the corresponding XYZ value and multiply it by the corresponding weight W Q(q). Perform weighted averaging on the weights in the same direction, i.e., perform interpolation reconstruction using 180-degree interpolation. In the embodiment of the present application, an extended image can be obtained by extension, and the weight of the position of the extended image can be adjusted (setting virtual horizontal motor parameters, and extending the weight of the effective weight W value), so as to increase the proportion of rays and reduce artifacts of incomplete projection. Finally, the superimposed axial image is multiplied by the angle weight, which is the back-projection image obtained after back-projection.

[0228] In an exemplary embodiment, referring to Fig. 20 As shown, in step 1804, weighted back projection processing is performed based on the extended projection data and the filtered data to obtain a back projection image, including the following steps 2002 and 2006, wherein:

[0229] Step 2002, performing weighted back-projection processing based on the extended projection data and the filtered data to obtain an initial back-projection image;

[0230] Step 2004, for any pixel point in the initial back-projection image, calculate the distance between the coordinate value of the pixel point on the vertical axis and the position of the horizontal motor corresponding to the central ray passing through the pixel point, and when the distance is less than a preset threshold value, determine that the pixel point is valid;

[0231] Step 2006, performing weighted sum processing on the rearranged data corresponding to the valid pixel points according to the Q value weight curve to obtain a back-projection image.

[0232] In an embodiment of the present application, a weighted back projection process is performed based on the extended projection data and the filtered data, and then an initial back projection image can be obtained. During the reconstruction process, effective reconstruction points can be screened from the initial back projection image. Exemplarily, the distance between the V coordinate value (coordinate value on the vertical axis) of each coordinate point (reconstruction point) and the position of the horizontal motor of the center ray of the ray passing through it can be calculated. If this distance is greater than or equal to the preset threshold value, the pixel point is considered invalid, otherwise the pixel point is considered valid.

[0233] Among them, V coordinate value = frame number * inc + startPos, wherein the frame number indicates which layer the reconstruction point is in, that is, the frame number, inc indicates the layer interval, that is, the distance between two frames, startPos indicates the distance for starting reconstruction (that is, the starting position of the horizontal motor movement), indicating the distance from the 0th image to the starting position of the horizontal motor; the position of the horizontal motor of the center ray of the ray passing through the reconstruction point is the position value corresponding to the frame number of the center ray read from the horizontal motor parameters, representing the distance from the frame to the initial horizontal motor position; the difference between the V coordinate value of the reconstruction point and the position of the horizontal motor of the center ray of the ray passing through it is the distance between the V coordinate value of the reconstruction point and the position of the horizontal motor of the center ray of the ray passing through it.

[0234] In fact, refer to Fig.22a As shown in the figure, the reconstructed target field of view is a sphere for each frame ray (it may also be an elliptical sphere in one example). If the reconstructed point is inside the sphere, it means that the reconstructed point is valid and can be used for reconstruction. If the reconstructed point is outside the sphere, it can be considered that the reconstructed point is too far from the source and is an edge ray. Due to errors and other reasons, the accuracy of the reconstructed point is poor. In this case, the reconstructed point can be considered invalid and will not be included in the superposition calculation (including not considering it when calculating the weight). Figure 22b As shown in the figure, the same reconstruction point of the two balls is a valid point for one ray and enters the reconstruction, but is an invalid point for the other ray and does not enter the reconstruction. This is equivalent to an effective three-dimensional simulated collimator, which can effectively improve the cone angle artifact and windmill artifact of the image.

[0235] In order to enable those skilled in the art to better understand the embodiments of the present disclosure, the embodiments of the present application are described below by means of specific examples. In this example, the spiral CT reconstruction method includes four stages: pre-reconstruction, coarse reconstruction, fine reconstruction and post-reconstruction. Specifically:

[0236] 1. In the pre-reconstruction stage, it includes calculation of the number of reconstruction layers, calculation of the starting position of the horizontal motor according to the positioning frame selection position, allocation of video memory, advance reading of data from various parameter files, establishment of a virtual coordinate system, customization of missing parameters, texture binding of weights, and obtaining of windowed filter kernels.

[0237] Among them, pre-reconstruction can be placed before data acquisition and does not take up acquisition and reconstruction time. Because the storage space must be opened up in advance (opening up video memory takes time), the story first needs to know the number of layers to be reconstructed. Among them, the reconstruction interval and reconstruction layer thickness are data used in the reconstruction process, and they are different in the coarse reconstruction and fine reconstruction processes. The number of reconstruction layers can include the number of fine reconstruction layers and the number of coarse reconstruction layers, where the number of fine reconstruction layers can be calculated based on the reconstruction parameters, including: according to the scan length scanlenth, collimation angle colangle, collimation width colwidth, extended collimation width ecol, pitch pitch, reconstruction interval inc, and display layer spacing thk, the number of fine reconstruction layers Num_bp is calculated. Among them, the reconstruction interval of coarse reconstruction is pre-designed, the reconstruction interval of fine reconstruction is selected by the user interface during shooting (or the user selects which protocol, the protocol defaults), the collimation angle, collimation width and extended collimation width are set before the image is acquired, and the collimation angle and collimation width are determined by the machine mechanics and are related to the mechanical design. The extended collimation width is set according to empirical values.

[0238] The process of determining the number of fine reconstruction layers is shown in the following formula (VII):

[0239]

[0240] Since there is a certain proportional relationship between the number of fine reconstruction layers and the number of coarse reconstruction layers, the number of coarse reconstruction layers can be further determined according to the proportional relationship. The process of determining the number of coarse reconstruction layers num_rough is shown in the following formula (VIII):

[0241] num_rough=(num_bp*inc-inc+5) / 5Formula (8)

[0242] Wherein, 5 is the reconstruction interval of 5 mm for the coarse reconstruction, and may be replaced by other values ​​when other reconstruction intervals are used for the coarse reconstruction.

[0243] The number of reconstructed images after the fine reconstruction layer thickness processing can be further determined based on the fine reconstruction layer number and the reconstruction interval, as well as the display layer spacing, as shown in formula (IX):

[0244]

[0245] Calculating the starting position of the horizontal motor movement according to the selected position of the positioning piece includes: the user selects the shooting area through the positioning piece, and the starting point of the scan needs to be moved to the starting point selected by the user during spiral scanning, and the scanning distance is the distance specified by the user in the positioning piece. The initial horizontal motor position startcouphpos, the distance scoutdis from the initial position selected in the positioning piece to the initial position of the positioning piece, the collimation angle colangle, the collimation width colwidth, the extended collimation width ecol, and the pitch pitch can be obtained, and the starting position startpos of the horizontal motor movement is calculated, wherein the initial horizontal motor position is the mechanical position of the displacement bed, that is, the zero point position of the displacement bed, and the starting position of the horizontal motor movement is the position of the horizontal motor when scanning starts from the starting position selected by the positioning piece, and the starting position of the horizontal motor movement is calculated using formula (10) or formula (11):

[0246]

[0247] In the scheme corresponding to formula (10), the moving direction of the predefined horizontal motor is opposite to the V-axis of the initial virtual coordinate system, and the opposite concept is reflected by the negative value; in the scheme corresponding to formula (11), the moving direction of the predefined horizontal motor is consistent with the V-axis of the initial virtual coordinate system.

[0248] Among them, the user selects the shooting area through the positioning film, and the starting point of the scan needs to be moved to the starting point selected by the user during spiral scanning. Therefore, the position can be selected according to the positioning film (the positioning film usually shows the geometric layout of the entire scanning area to help operators and doctors determine the specific position and range of the scan. The starting position of the scan is pre-defined on the positioning film, which is called the initial position of the positioning film. The operator will use the mouse or other tools to select another area on the positioning film. This selected area is usually to determine the region of interest (ROI) for further CT scanning. The selected position can help ensure that the scanner acquires data at the correct position and angle. The selected position also includes a starting position for the scan, which is called the initial position selected in the positioning film. For example: when scanning the chest, the positioning film may show the outline of the entire chest, and the radiologist will select the lung area so that the CT scanner can focus on acquiring data from the lungs), determine the starting position of the horizontal motor movement, and move the displacement bed according to the starting position of the horizontal motor movement.

[0249] In the process of developing video memory, the video memory is developed and initialized in advance to avoid developing and releasing video memory during the reconstruction process, which greatly improves the reconstruction time. Specifically, the opened video memory space includes space for storing the following data: tomographic images after layer thickness superposition reconstruction; tomographic images after back projection; total projection data images; detector geometric position parameter index table (X, Y, Z of each detector point, in a virtual coordinate system); actual horizontal motor correction parameter index table (the position of the patient support, that is, the bed, in each frame of the scan, in mm value in the Z direction); uniform horizontal motor correction parameter index table (used for coarse reconstruction, average calculated based on pitch, collimation width, etc.); air data image (storing air images after dark field correction); angle data of each frame of rays relative to the central ray (in the UV coordinate system); coordinate value distance data of each frame of rays in the V direction (in the UV coordinate system); directional angle interpolation data (coordinate value distance data of each ray in the V direction after directional angle rearrangement); rearranged image; filter kernel data; filtered image; water uniformity correction parameter index table; water CT value correction parameter index table; bone sclerosis correction parameter index table; detector odd-even frame correction ratio parameter index table; segmented projection data images, etc.

[0250] The data of each parameter file is read in advance, including air correction data, detector position data, inconsistency parameters, water correction parameters, bone sclerosis parameters, etc. The parameter data is read into the GPU in advance to avoid reading the hard disk and switching from CPU to GPU during the reconstruction process, which greatly increases the reconstruction time.

[0251] During the establishment of the virtual coordinate system, the center point of the detector (the point corresponding to the source, such as the detector offset placement, non-center point) is used as the origin to establish a two-dimensional UV coordinate system. △U and △V are calculated according to the preset sampling size, and the virtual coordinate U channel size is calculated. The U axis of the virtual coordinate system is not a linear axis, but a curved axis. Under each frame of the image, the position of the detector irradiated by the central ray of the source irradiation range (cone beam) is used as the origin, the spiral line is the U axis, and the V axis is in the same direction as the horizontal displacement bed. The UV wrapped area is an irregular image, see Figure 7 In the midpoint matrix, the values ​​are taken at equal intervals according to the wrapping distance in the V direction, and at equal curve distances in the spiral direction in the U direction.

[0252] Among them, like oral scan registration and fusion, CT also needs coordinate system transformation, but the coordinate system of CT is simpler, there are only two, one is the target coordinate system, which can be regarded as the world coordinate system, that is, the x, y, z of the reconstructed image, and the other is the acquired geometric coordinate system, which is a frame of image, that is, the coordinate system of the detector at a fixed position. The length and width directions u and v of the detector are the coordinate systems x and y, and the direction from the vertical detector to the focus of the source is z. This coordinate system can be a virtual coordinate system. Back projection includes the transformation from the virtual coordinate system to the target coordinate system. The back projection process is the superposition of the target position in each virtual coordinate system, so there is this transformation. One purpose of rearrangement is to facilitate back projection interpolation. Another rearrangement process is also an interpolation process, that is, for example, the original 384 points are rearranged into 512 points, which is equivalent to making the data dense, which is beneficial to image quality.

[0253] Three-dimensional geometric position display reference after spiral CT rearrangement Fig.11a , the detector UV direction plane display reference after spiral CT rearrangement Fig.11b As shown, taking the plane display as an example, the U axis of the coordinate system is the horizontal axis (but it is actually a curved surface, the curved surface where the spiral line is located, or the spiral line is attached to the outer surface of the cylinder, and the U axis is the curve axis of the outer surface of the cylinder (the curve axis is perpendicular to the central axis of the cylinder), and the V axis of the coordinate system is the vertical axis (parallel to the central axis of the cylinder). The advantage of this is that the interpolation change is accurate, that is, the points rearranged to the detector are not straight lines with completely equal intervals, but points on an oblique surface. This position interpolation is more in line with reality than the point interpolation position change of the approximate square dot matrix.

[0254] Customization of missing parameters, including: customizing horizontal motor parameters according to pitch, position parameters (that is, the mm value of the horizontal motor encoder for each frame of image, how many values ​​there are for how many frames of image are taken), etc. It is used for coarse reconstruction, that is, coarse reconstruction takes approximate horizontal motor parameters.

[0255] During the texture binding process of the weight, the specific calculation formula of the weight is shown in formula (VI).

[0256] In the process of obtaining the window filter kernel, the window function is a sinc window. Compared with the hanning window with the same cutoff frequency, the details are clearer, and the hanning is too smooth. Take the low frequency and high frequency limits. When the window value is less than the low frequency, it is 1, that is, the limit is fully open, which is called the full-pass frequency. When it is greater than the high frequency, the window value is 0, that is, the high-frequency limit is the cutoff frequency. According to the requirements of reconstructing the image, the full-pass frequency and the cutoff frequency take different values, which are divided into 7 groups, corresponding to different sharpening and smoothing effects.

[0257] The positions of the full-pass frequency and the cutoff frequency, the horizontal axis is the sampling point (positive and negative normalization, that is, the center is the center of the image sampling width), and the vertical axis is the window function value. As shown in Figure 4, WL is the full-pass frequency (values ​​1, 1 / 16, 1 / 32), WH is the cutoff frequency (values ​​1, 2 / 3, 1 / 2), WL<=WH. The horizontal axis is the sampling width of the convolution (the value of 0 represents the sampling center), and the vertical axis is the normalized window width value. The formula for the windowed filter kernel is shown in formula (three). In this way, the filter kernel is prepared in the preprocessing, that is, before acquisition. The texture binding above is also done before acquisition, and does not take up acquisition and reconstruction time.

[0258] 2. In the coarse reconstruction stage, in order to display the reconstruction results in real time and facilitate the doctor to judge the acquisition situation, coarse reconstruction is performed. The reconstruction interval, display layer spacing and layer thickness of coarse reconstruction are fixed to 5mm. It includes: logarithmic calculation, water hardening correction, CT value correction, rearrangement, FFT, fast back projection and HU value calculation.

[0259] Coarse reconstruction physically divides the original acquisition according to the acquisition time to ensure that the reconstruction time is less than the acquisition time (acquisition and reconstruction are two independent threads). The horizontal motor speed (bed stepping speed) is calculated according to the average of the pitch and other parameters set by the user, without using the real-time speed (this way, there is no need to read the horizontal motor file, and the reconstruction speed is faster).

[0260] When calculating the logarithm, the initial angle of the air correction parameter shooting and the initial angle of the collected data are considered, and the offset value is added to make the angles correspond. The specific process is shown in formula (IV), and this example will not be repeated here.

[0261] And considering that the ray is not always stable after it is turned on, for example, it shows a slow attenuation trend. Although the amplitude is not large and meets the requirements of source stability, it still has a certain impact on the image brightness level correction. Therefore, the influence factor Δ is added, that is, the difference between the air average value under the prj initial angle period and the air average value under the air initial angle period.

[0262] In the water hardening correction process, polynomial fitting correction is adopted, and the coefficients of the polynomial are determined according to the water models taken multiple times.

[0263] In the CT value correction process, the advantage of placing the slope correction before the back projection is that the dynamic range of prj can be enlarged and the data accuracy can be improved (when the same float type). Generally, the CT value correction is placed after the reconstruction, and the dynamic range of the reconstruction result is stretched so that the air reaches -1000 and the water reaches 0, that is, linear correction, and the slope and intercept are obtained to correct each pixel point. In the disclosed embodiment, it is divided into two steps. Before reconstruction, the slope part k of the CT value is corrected first. Because the projection value is the superposition of pixel points, it can be equally scaled, and the implementation of the slope method effectively expands the dynamic range of the data and improves the reconstruction accuracy. After the reconstruction is completed, the intercept part b=1000 is processed. The correction of the CT value is all linear stretching.

[0264] The CT value correction formula is as follows:

[0265] In the rearrangement process, first: Direction angle rearrangement: fan beam is rearranged into parallel beam. With this direction angle IV (refer to Fig. 9 , the direction angle is directly above, set as the starting angle 0 degrees) and the rays parallel to the central ray are rearranged to the non-central ray of the direction angle, that is, the non-central ray is parallel to the central ray, and the interval after rearrangement is smaller than the original interval, because the non-central position is at any position, such as hitting the coordinate (iu, iv) on the detector obliquely. The ray is rotated through the rotation center to the ray parallel to the central ray under the direction angle IV ( Fig. 9 The example shown is rotated 20.35 degrees) closer to the center: AI <DJ。

[0266] After the direction angles are rearranged, the channel spacing is not equal. The 384 channel spacings are as follows: Fig.10 As shown, it is related to the installation rules of the detector. The intervals between channels after the azimuth angle is rearranged (the horizontal axis is the channel interval, the vertical axis is the interval distance in mm, from left to right, it is the interval between the 1st channel and the 0th channel, the interval between the 2nd channel and the 1st channel, ...).

[0267] Secondly: According to the up-sampled data volume and the azimuth rearrangement interval, perform cubic Hermite interpolation (because the rearranged single-view detector data is not a standard two-dimensional matrix, but a curved matrix related to the spiral track, the cubic spline interpolation effect is better than general linear interpolation). After the azimuth rearrangement, the detector direction is a three-dimensional surface, which is similar to a parallelogram in two dimensions. That is, after the Hermite interpolation in the previous step, the V-direction interval is the same, but because the U-direction is not horizontal, directly using the data during back-projection will result in inappropriate data for interpolation. This is especially obvious when the pitch is low. Fig.12Approximate diagram (because it is not a plane but a curved surface). The point used for bilinear interpolation of point G1 in this coordinate system is GMLH. Point H is obviously not suitable, and GMLR is the best. So here we add a step to interpolate the vertical UV coordinates to make the U axis horizontal. That is, set the horizontal grid according to the sampling density. And save the data within the boundary range of the original coordinate system as valid data, and the peripheral data as invalid data.

[0268] In the FFT process, first, 2-power zero-padding extension is performed; second, FFT transformation is performed; third, frequency domain filtering is performed (using a pre-reconstruction filter kernel, refer to the aforementioned embodiment); then, FFT inverse transformation is performed; and finally, multiplication by the coefficient △ is performed.

[0269] In the assignment stage, the results of the projection processing are stored in the overall pre-back-projection data, so that these processing do not need to be considered during fine reconstruction, and only back-projection is required.

[0270] In the fast back-projection process, since the purpose of coarse reconstruction is only to display the approximate shooting position, fast back-projection is used in this example, and the image accuracy requirement is not high. Fast back-projection is achieved by the following means: back-projection with a thicker reconstruction interval, 360-degree interpolation can be used instead of 180-degree interpolation, so that the weight calculation is reduced by half, the position of the horizontal motor corresponding to each frame can be replaced by a uniform average value calculated according to the pitch, that is, the intervals are equal, and the horizontal interval can be directly multiplied according to the frame number, and there is no layer thickness superposition calculation, and the display result is directly displayed (the layer thickness is equal to the reconstruction interval).

[0271] During the HU value correction process, refer to the aforementioned CT value correction formula and perform a subtraction process of 1000, or perform a window position addition process of 1000 during display.

[0272] 3. In the fine reconstruction stage, the fine reconstruction display layer spacing and layer thickness are re-assigned according to the preset values ​​in the user interface. (It is set to 5mm during the coarse reconstruction. If the layer thickness in the user agreement is relatively small and the calculation amount is relatively large, the fine reconstruction is based on the user settings.) The specific process of fine reconstruction includes:

[0273] Read the horizontal motor parameters, and the physical position of the horizontal displacement bed corresponding to each frame of the acquired image. Coarse reconstruction is an approximate horizontal motor parameter, that is, based on the average pitch interval, but in fact the horizontal displacement of each frame is not equal. Fine reconstruction is based on the actual position and has higher accuracy.

[0274] Reconstruct the extension of the Z-direction field of view, expand the reconstruction range, and improve the dose utilization rate. The extended half-width range is calculated using formula (V). For example, the position of the upper and lower extension of the spiral scan is referenced Fig.19As shown. By adjusting the weight of the position where the extended image is located (setting the virtual horizontal motor parameters, and extending the effective weight W value), the proportion of rays can be increased and the artifacts of incomplete projection can be reduced. The weight that needs to be calculated for extension is to consider the extended horizontal motor, such as the w weight value at the -10 position, and add the denominator weighted weight when superimposing the projection.

[0275] According to the pitch and collimation width, the horizontal corresponding position of the largest reconstruction point to the center line of the superposition ray is determined. This process is equivalent to the reconstructed target field of view being a sphere for each frame of ray (of course it may be an elliptical sphere). That is, the reconstruction point is valid if it is inside the sphere, and the reconstruction point is invalid if it is outside the sphere, and does not enter the superposition calculation (of course it is not considered when calculating the weight), which is equivalent to an effective three-dimensional simulated collimator. This method can effectively improve the cone angle artifact and windmill artifact of the image. The main reason for this is that the reconstruction point is too far away from the source and the edge ray itself is of low quality due to errors, so it does not enter the superposition calculation. Figure 22b For the same reconstruction point on the two spheres, one ray is a valid point and enters the reconstruction, while the other ray is a non-valid point and does not enter the reconstruction.

[0276] In the back-projection process, each pixel of the XY plane (each axial diagram) of the target coordinate system and each acquisition rearrangement perspective is traversed to determine the corresponding starting and ending positions of the horizontal motor under the rearrangement perspective, and the upper and lower boundaries of the effective Z plane under the acquisition rearrangement perspective are determined according to the semi-extension range. The pixel points on XY are converted to the virtual coordinate system to determine the UV coordinate value, and the channel position is determined according to the U coordinate value, that is, the position of the horizontal channel of the acquisition rearrangement data (only the effective horizontal channel is calculated); according to the V-direction position of the ray, the V coordinate value determined by converting the pixel points on XY to the virtual coordinate system and the V position (i.e., Z position) of the desired image plane, the row number is determined, that is, the position of the longitudinal channel of the acquisition rearrangement data (only the effective longitudinal channel is calculated). And according to the position of the effective longitudinal channel of the acquisition rearrangement data, the Q value weight W is calculated according to the following formula (six): Q (q).

[0277] Iterate the rearranged data corresponding to the corresponding XYZ values ​​and multiply them by the corresponding weights. The weights in the same direction are weighted averaged, that is, 180-degree interpolation is performed (compared to the 360-degree interpolation of rough reconstruction, 180-degree interpolation has more Fig.24 Finally, the superimposed axial image is multiplied by the angle weight to obtain the back-projected image.

[0278] In the post-reconstruction stage, the bone sclerosis correction, CT value correction and layer thickness processing processes are included.

[0279] Bone sclerosis correction is actually the optimization of high-density body artifacts. It is generally divided into two methods, one is pre-processing and the other is post-processing. Pre-processing generally requires modeling based on a preset model. Since the bone density and morphological steps of each patient are different, there will be certain deviations and the application is limited. Post-processing is also roughly divided into two types, one is the orthographic projection reconstruction method of the remaining tissues after removing the bone tissue, and the other is the orthographic projection reconstruction method of the bone tissue after extracting the bone tissue. In the present disclosed embodiment, first, the bone tissue segmentation adopts segmented correction segmentation, and secondly, the orthographic projection is not a projection according to the rotation of the axial scan, but a spiral scan projection method. The steps are as follows:

[0280] First: segment the bone tissue (first of all, it should be noted that the CT value -1000 has not been corrected in this step, so these values ​​are 1000 more than the normal CT value). The bone image is segmented. Since the CT value of the bone image varies in each axial layer, each part, and each different person, it is difficult to accurately segment it with the same threshold. Therefore, a two-step method is adopted. First, the bone tissue part is extracted, that is, the part with CT values ​​between hu1 and hu2 (such as 850 to 3000). This part tries to include all situations and has a wide range. The CT values ​​of other parts are set to 0. That is, all bone parts are included in hu1 to hu2. Then the CT values ​​between hu1 and hu2 are segmented. We divide it into three segments here, or it can be divided into multiple segments. The first segment hu1~hua (such as 850~1150) is the CT value between brain tissue and bone tissue, the second segment hua~hub (such as 1150~2075) is the CT value range of bone tissue, and the third segment hub~hu2 (such as 2075~3000) is the CT value between bone tissue and high-density body. By presetting the normalization coefficient of each node, for example, hu1 corresponds to s1 (such as 0), hua corresponds to s2 (such as 0.4), hub corresponds to s3 (such as 0.7), and hu2 corresponds to s4 (such as 1). Then, the coefficient s is calculated for each point of the segmented bone tissue part of hu1~hu2, si is the low-end coefficient value of the segment where the point is located, sj is the high-end coefficient value of the segment where the point is located, hu is the CT value of the point, hui is the low-end CT value of the segment where the point is located, and huj is the high-end CT value of the segment where the point is located. The CT value of the bone tissue part is multiplied by this coefficient to obtain the final segmented and extracted bone tissue.

[0281] Secondly, the bone tissue is projected orthographically. First, the projection path is as follows: Figure 2As shown, it can be regarded as a parallel beam of rays with a cone angle. Its scanning path follows a spiral line as a whole. The advantage of this is that a pair of original images have a projection of this position (equivalent to the result of the projection image after rearrangement), which can be directly matched. Second, this saves the step of rearrangement in the reconstruction and increases the speed. Third, the positive projection of this position is downsampled to take the projection line according to half of the rearrangement (the sampling can be reduced according to the demand for the amount of calculation), which can reduce the amount of calculation, and there is no need to process the projection of the original image. It only needs to shift the x-coordinate to the left when taking the value of the original image.

[0282] Iterate according to the axial image pixels and ray length weights, and calculate the total value of each ray in the orthographic projection. The ray length weight is inversely proportional to the ray length, that is, the longer the ray in the voxel, the smaller the influence of the point on the axial plane (the orthographic weight is the largest). That is, the weight is AG / IJ (such as Figure 4b , where AG is the most diagonal length of the voxel and IJ is the path length of one of the rays passing through the voxel).

[0283] Again, bone sclerosis parameter correction. Weighted correction is performed on the bone projection image. First, the result of the above positive projection is subtracted from the original image projection (of course, there is a left shift, that is, division by 2). The result is the projection of the tissue without bone tissue. This projection is divided by the water reference to obtain the equivalent length of each projection pixel relative to water. The equivalent length is used to determine the artifact correction coefficient, and this value is used to correct the result of the above positive projection.

[0284] Finally: the corrected bone projection image is reconstructed (there is no need to rearrange the previous steps in the reconstruction step, because the scanning path of the projection image is always the rearranged path) to obtain a new bone reconstruction image. The original image and the new bone reconstruction image are superimposed to obtain the final corrected image.

[0285] During the CT value correction process, a 1000 offset correction is used, that is, the values ​​of all pixels are subtracted by 1000.

[0286] In the layer thickness processing process, layer thickness superposition calculation is performed. The effective layer thickness is obtained according to the reconstruction interval, display layer spacing and preset layer thickness, and layer thickness superposition calculation is performed, such as Fig.17a and Fig.17b The data is stored in the specified memory (the space previously allocated in the pre-reconstruction).

[0287] Repeat the above steps until all the collected data are processed. The time for one acquisition is set to 1 second or 2 seconds according to the user, and the time for one rough reconstruction is no more than 0.2s (ordinary computer configuration, high-end computer configuration is much less than 0.2a). The reconstruction time is much shorter than the acquisition time, so the axial reconstruction image of the current shooting position can be displayed almost in real time, thereby achieving real-time reconstruction display, which is convenient for doctors to understand the shooting situation in time. If there is a problem with the shooting position, it can be cancelled in time without confirmation after the shooting is completed.

[0288] The overall reconstruction process includes coarse reconstruction, fine reconstruction and post-reconstruction, where coarse reconstruction, fine reconstruction and post-reconstruction are displayed simultaneously, including the following steps:

[0289] After the ray collection button is pressed, the first thing is delayed exposure. For the delayed exposure interface, refer to Fig.25 As shown, the left picture is the picture of the positioning film, and the right picture is the picture of the spiral scan. At this time, the right picture is completely black.

[0290] Complete pre-reconstruction before acquisition. This can be done during the delay or before multi-thread synchronization. This time is less than 4 seconds. If the delay time is too short, it is recommended to put it in front to ensure that this time does not delay the subsequent acquisition and reconstruction. In this example, it is actually placed after the device preparation button is clicked. After the device is prepared, some doctors have time to check, such as making sure that the patient has not moved, etc. The device is ready to pop up the delayed exposure interface above, and the delay starts when you click to shoot.

[0291] After the positioning film is taken, the spiral scanning parameters are set and the device can be prepared. Pre-reconstruction is placed in this process, that is, pre-reconstruction is completed before acquisition, which does not take up acquisition and reconstruction time and the user does not feel it.

[0292] After scanning a circle, the right picture starts to display the results, that is, the acquisition is performed while the coarse reconstruction is performed. In this example, it is set to coarse reconstruction after one circle. According to the acquisition time setting, if the acquisition time of one circle is too long, it can be shortened accordingly (half circle, etc.). In this example, it takes 1-2 seconds to acquire one circle, and only 0.1+ seconds to coarsely reconstruct one circle, so it is completely sufficient.

[0293] How many frames are collected in one circle depends on the computer configuration. This example sets less than 2000 frames (1440 frames in one circle). If the computer performance is good, you can set more. If the computer performance is low, you can set less to ensure real-time effect. It can be seen as segmenting the image to achieve the purpose of simultaneous acquisition and reconstruction. The callback function can be used to inform the system software of how many images have been reconstructed. The callback function in the coarse reconstruction can return the total number of axial images of the coarse reconstruction and the number of axial images that have been reconstructed in real time. Because the coarse reconstruction also needs to consider the acquisition time (in fact, the reconstruction time can be ignored relative to the acquisition), that is, the number of reconstructed axial images needs to be displayed according to the acquisition time of 2s. For example, if 2 new images are reconstructed, then 1 image is displayed in 1s. The coarse reconstruction time interval is set to -1, because the main time is not the reconstruction time, but the acquisition time (the coarse reconstruction truncation reconstruction time is less than the acquisition time).

[0294] In the post-reconstruction process, the post-reconstruction module is bone sclerosis, CT value intercept correction, slice thickness calculation, etc. The time of post-reconstruction is synchronously placed in the callback of fine reconstruction, that is, after the time interval of fine reconstruction is displayed, the post-reconstruction has been completely reconstructed. Therefore, it is recommended to refresh the displayed data after the fine reconstruction progress bar is fully completed.

[0295] Using the spiral CT reconstruction method provided by the embodiment of the present disclosure, in order to achieve real-time reconstruction (reconstruction while acquiring), spiral CT currently generally adopts methods such as using a faster processor (high-end reconstruction computer), uploading reconstruction to a cloud server, or a new reconstruction algorithm. In this example, the reconstruction process is segmented, low-precision coarse reconstruction is used for display during acquisition, and high-precision fine reconstruction is used to save the reconstruction results. This ensures the needs of real-time acquisition and reconstruction, and the final accuracy. In addition, a pre-reconstruction stage is added, and part of the reconstruction content is completed before acquisition. The spiral CT reconstruction method provided in this example can achieve real-time reconstruction even with a low-end computer, and does not affect the doctor's real-time supervision of the shooting process. Since sudden situations such as patient movement may occur during the shooting process, it is necessary to detect in real time whether the scanning position is normal.

[0296] Before the rough reconstruction, judgment is made, and the rough reconstruction step can be selected: according to the different parameters selected by the operator, it is judged whether the rough reconstruction is needed. When the reconstruction calculation amount is small, it directly enters the fine reconstruction. In addition, the rough reconstruction includes correction, that is, scanning, acquisition, correction, reconstruction, and display are all completed in real time.

[0297] Setting of UV virtual coordinate system: The whole coordinate system is an approximate rectangular surface, and the U axis is the curve axis. The coordinate system can be considered as two, one is the former UV virtual coordinate system, which performs curve fitting interpolation according to the U axis (not the final U axis, the U axis of the former UV virtual coordinate system, that is, the axis of the spiral line) when rearranged, not linear interpolation. The other is the final UV virtual coordinate system, and the points after reshooting perform grid node interpolation according to the final UV coordinate system, that is, the oblique points are made straight, which improves the accuracy of back projection.

[0298] The fast back-projection of coarse reconstruction and the back-projection of fine reconstruction have greatly improved the speed. And the coarse reconstruction algorithm can be adjusted according to the acquisition time. For example, if the acquisition time is 1-2 seconds, and the acquisition time is much longer than the time for one circle of coarse reconstruction, the coarse reconstruction can be divided into more than one circle, and the coarse reconstruction can be completed within the acquisition time. If the acquisition time is subsequently increased, such as the fastest spiral CT in the industry is 0.28s, the first six steps of coarse reconstruction can be divided into shorter frames, such as half a circle, or even one frame, and only one full circle is performed during back-projection. It can also ensure that coarse reconstruction and acquisition are carried out at the same time, and do not take up time separately.

[0299] The spiral CT reconstruction provided in this example effectively reduces cone beam artifacts and also improves windmill artifacts. Collimation is usually added by hardware, but here software simulation is used to increase collimation to achieve the same effect of removing cone beam artifacts.

[0300] It should be understood that, although the various steps in the flowcharts involved in the above-mentioned embodiments are displayed in sequence according to the indication of the arrows, these steps are not necessarily executed in sequence according to the order indicated by the arrows. Unless there is a clear explanation in this article, the execution of these steps does not have a strict order restriction, and these steps can be executed in other orders. Moreover, at least a part of the steps in the flowcharts involved in the above-mentioned embodiments can include multiple steps or multiple stages, and these steps or stages are not necessarily executed at the same time, but can be executed at different times, and the execution order of these steps or stages is not necessarily to be carried out in sequence, but can be executed in turn or alternately with other steps or at least a part of the steps or stages in other steps.

[0301] Based on the same inventive concept, the embodiment of the present application also provides a spiral CT reconstruction device for implementing the spiral CT reconstruction method involved above. The implementation scheme for solving the problem provided by the device is similar to the implementation scheme recorded in the above method, so the specific limitations in one or more spiral CT reconstruction device embodiments provided below can refer to the limitations on the spiral CT reconstruction method above, and will not be repeated here.

[0302] In an exemplary embodiment, Fig.23As shown, a spiral CT reconstruction device is provided, including: a filtering module 2302, a back projection module 2304, a forward projection module 2306, a processing module 2308 and a reconstruction module 2310, wherein:

[0303] A filtering module 2302 is used to filter the projection data to obtain filtered data;

[0304] A back-projection module 2304 is used to perform back-projection processing on the filtered data to obtain a back-projection image, and to extract bone tissue data from the back-projection image;

[0305] An orthographic projection module 2306 is used to perform orthographic projection on the bone tissue data to obtain bone projection data;

[0306] A processing module 2308, configured to use the difference between the projection data and the bone projection data as bone removal tissue data;

[0307] The reconstruction module 2310 is used to obtain a CT reconstructed image based on the bone removal tissue data and the back-projection image.

[0308] The above-mentioned spiral CT reconstruction device can filter the projection data to obtain filtered data, and perform back-projection processing on the filtered data to obtain the back-projection image, and then extract the bone tissue data from the back-projection image. The bone tissue data is forward projected to obtain bone projection data, and the difference between the projection data and the bone projection data is used as bone-removed tissue data, and then a CT reconstructed image is obtained based on the bone-removed tissue data and the back-projection image. When the spiral CT reconstruction device provided in the embodiment of the present application is used, in the post-reconstruction process, there is no need to perform a forward projection on the original data, and the reconstructed projection data can be directly used, which can effectively reduce the amount of calculation and improve the efficiency and accuracy of CT reconstruction.

[0309] In one embodiment, the back projection module 2304 is specifically used to:

[0310] Extracting bone tissue pixel points in the back-projection image, wherein the CT values ​​corresponding to the bone tissue pixel points are within the range of bone tissue CT values;

[0311] For any bone tissue pixel point, determine the bone tissue segment corresponding to the bone tissue pixel point according to the CT value of the bone tissue pixel point, and determine the normalization coefficient of the bone tissue pixel point according to the upper limit CT value and the lower limit CT value of the bone tissue segment, the normalization coefficients corresponding to the upper limit CT value and the lower limit CT value, and the CT value of the bone tissue pixel point, wherein the bone tissue CT value range is composed of a plurality of the bone tissue segments;

[0312] Bone tissue data is obtained according to the normalization coefficient and CT value of each bone tissue pixel point.

[0313] In one embodiment, the orthographic projection module 2306 is specifically used to:

[0314] According to each pixel point in the axial image corresponding to the ray passing through the bone tissue data and the ray length weight corresponding to each pixel point, the CT value of each pixel point is weightedly superimposed to obtain the bone projection data;

[0315] The ray length weight is negatively correlated with the length of the ray in the voxel.

[0316] In one embodiment, the reconstruction module 2310 is specifically used to:

[0317] The ratio of the bone-removed tissue data to the water reference is used as the equivalent length of each projection pixel point in the bone-removed tissue data relative to water;

[0318] Performing weighted correction on the bone projection data using the artifact correction coefficient corresponding to the equivalent length to obtain corrected bone projection data;

[0319] After the corrected bone projection data is reconstructed, the reconstruction result is superimposed on the back-projection image to obtain a CT reconstructed image.

[0320] In one embodiment, the filtering module 2302 is specifically used for:

[0321] The projection data are rearranged using a pre-established initial virtual coordinate system, wherein the initial virtual coordinate system takes the detector irradiated by the central radiation source as the origin, the spiral line as the horizontal axis, and the horizontal displacement bed direction as the vertical axis;

[0322] The rearranged projection data is filtered according to the windowing filter kernel to obtain filtered data.

[0323] In one embodiment, the filtering module 2302 is further configured to:

[0324] Rearranging the projection data of the fan beam into the projection data of the parallel beam by rearranging the direction angle of the projection data;

[0325] According to the up-sampled data amount and the interval of the direction angle rearrangement, the projection data of the parallel beam is subjected to three-fold Hermite interpolation to obtain a three-dimensional surface, and the three-dimensional surface is subjected to coordinate interpolation processing perpendicular to the initial virtual coordinate system to obtain rearranged projection data.

[0326] In one embodiment, the filtering module 2302 is further configured to:

[0327] Determine the grid points in the three-dimensional surface and calculate the rotation matrix of each of the grid points;

[0328] The rotation matrix is ​​used to transform the coordinates of each of the grid points from the initial virtual coordinate system to a target virtual coordinate system, wherein the horizontal axis of the target virtual coordinate system is horizontal;

[0329] Interpolation processing is performed on each of the grid points in the target virtual coordinate system to obtain rearranged projection data.

[0330] In one embodiment, the filtering module 2302 is further configured to:

[0331] Performing interpolation processing on each of the grid points in the target virtual coordinate system to obtain an interpolation result;

[0332] The interpolation result is filtered according to the initial virtual coordinate system to obtain rearranged projection data.

[0333] In one embodiment, the windowing filter includes an all-pass frequency and a cutoff frequency, and the windowing filter kernel is:

[0334]

[0335] Wherein, the w represents the frequency component of the signal, the wl represents the all-pass frequency, and the wh represents the cut-off frequency.

[0336] In one embodiment, the device further comprises:

[0337] A correction module, used for performing slope correction on the projection data;

[0338] The reconstruction module 2310 is further specifically used for:

[0339] The reconstruction result is superimposed with the back-projection image to obtain a superimposed result;

[0340] The stacking result is intercept-corrected, and the correction result is slice-thickness processed to obtain a CT reconstructed image.

[0341] In one embodiment, the back projection module 2304 is specifically used to:

[0342] Determine the extension range according to the extension collimation width, the pitch and the reconstruction interval, and shift the starting position of the horizontal motor to the left and right based on the extension range to obtain a virtual horizontal motor position;

[0343] Interpolation processing is performed based on the virtual horizontal motor position and the filtered data to obtain extended projection data, and weighted back projection processing is performed based on the extended projection data and the filtered data to obtain a back projection image.

[0344] In one embodiment, the back projection module 2304 is further configured to:

[0345] Performing weighted back-projection processing based on the extended projection data and the filtered data to obtain an initial back-projection image;

[0346] For any pixel point in the initial back-projection image, calculate the distance between the coordinate value of the pixel point on the depth axis and the position of the horizontal motor corresponding to the central ray passing through the pixel point, and determine that the pixel point is valid when the distance is less than a preset threshold value;

[0347] According to the Q value weight curve, the rearranged data corresponding to the effective pixel points are processed by weighted summation to obtain the back-projection image.

[0348] Each module in the above spiral CT reconstruction device can be implemented in whole or in part by software, hardware or a combination thereof. Each module can be embedded in or independent of a processor in a computer device in the form of hardware, or can be stored in a memory in a computer device in the form of software, so that the processor can call and execute the operations corresponding to each module.

[0349] In an exemplary embodiment, a computer device is provided. The computer device may be a terminal, and its internal structure diagram may be as shown in FIG. Fig.26 As shown. The computer device includes a processor, a memory, an input / output interface, a communication interface, a display unit and an input device. The processor, the memory and the input / output interface are connected via a system bus, and the communication interface, the display unit and the input device are connected to the system bus via the input / output interface. The processor of the computer device is used to provide computing and control capabilities. The memory of the computer device includes a non-volatile storage medium and an internal memory. The non-volatile storage medium stores an operating system and a computer program. The internal memory provides an environment for the operation of the operating system and the computer program in the non-volatile storage medium. The input / output interface of the computer device is used to exchange information between the processor and an external device. The communication interface of the computer device is used to communicate with an external terminal in a wired or wireless manner, and the wireless manner can be implemented through WIFI, a mobile cellular network, NFC (near field communication) or other technologies. When the computer program is executed by the processor, a spiral CT reconstruction method is implemented. The display unit of the computer device is used to form a visually visible picture, which can be a display screen, a projection device or a virtual reality imaging device. The display screen can be a liquid crystal display screen or an electronic ink display screen, and the input device of the computer device can be a touch layer covering the display screen, or a button, trackball or touchpad set on the computer device shell, or an external keyboard, touchpad or mouse.

[0350] Those skilled in the art will understand that Fig.26 The structure shown in the figure is only a block diagram of a part of the structure related to the solution of the present application, and does not constitute a limitation on the computer device to which the solution of the present application is applied. The specific computer device may include more or fewer components than those shown in the figure, or combine certain components, or have a different arrangement of components.

[0351] In one embodiment, a computer device is further provided, including a memory and a processor, wherein a computer program is stored in the memory, and the processor implements the steps in the above method embodiments when executing the computer program.

[0352] In one embodiment, a computer-readable storage medium is provided, on which a computer program is stored. When the computer program is executed by a processor, the steps in the above-mentioned method embodiments are implemented.

[0353] In one embodiment, a computer program product is provided, including a computer program, which implements the steps in the above method embodiments when executed by a processor.

[0354] It should be noted that the user information (including but not limited to user device information, user personal information, etc.) and data (including but not limited to data used for analysis, stored data, displayed data, etc.) involved in this application are all information and data authorized by the user or fully authorized by all parties, and the collection, use and processing of relevant data must comply with relevant regulations.

[0355] Those skilled in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by instructing the relevant hardware through a computer program, and the computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it can include the processes of the embodiments of the above-mentioned methods. Among them, any reference to the memory, database or other medium used in the embodiments provided in the present application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetoresistive random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory can include random access memory (RAM) or external cache memory, etc. As an illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM). The database involved in each embodiment provided in this application may include at least one of a relational database and a non-relational database. Non-relational databases may include distributed databases based on blockchains, etc., but are not limited to this. The processor involved in each embodiment provided in this application may be a general-purpose processor, a central processing unit, a graphics processor, a digital signal processor, a programmable logic device, a data processing logic device based on quantum computing, etc., but are not limited to this.

[0356] The technical features of the above embodiments may be combined arbitrarily. To make the description concise, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0357] The above-described embodiments only express several implementation methods of the present application, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the present application. It should be pointed out that, for a person of ordinary skill in the art, several variations and improvements can be made without departing from the concept of the present application, and these all belong to the protection scope of the present application. Therefore, the protection scope of the present application shall be subject to the attached claims.

Claims

1. A spiral CT reconstruction method, characterized in that: The method comprises: Performing filtering on the projection data to obtain filtered data; Performing back-projection processing on the filtered data to obtain a back-projection image, and extracting bone tissue data from the back-projection image; Performing orthographic projection on the bone tissue data to obtain bone projection data; Taking the difference between the projection data and the bone projection data as bone removal tissue data; A CT reconstructed image is obtained according to the bone removal tissue data and the back-projection image.

2. The method according to claim 1, characterized in that The step of extracting bone tissue data from the back-projection image comprises: Extracting bone tissue pixel points in the back-projection image, wherein the CT values ​​corresponding to the bone tissue pixel points are within the range of bone tissue CT values; For any bone tissue pixel point, determine the bone tissue segment corresponding to the bone tissue pixel point according to the CT value of the bone tissue pixel point, and determine the normalization coefficient of the bone tissue pixel point according to the upper limit CT value and the lower limit CT value of the bone tissue segment, the normalization coefficients corresponding to the upper limit CT value and the lower limit CT value, and the CT value of the bone tissue pixel point, wherein the bone tissue CT value range is composed of a plurality of the bone tissue segments; Bone tissue data is obtained according to the normalization coefficient and CT value of each bone tissue pixel point.

3. The method according to claim 1 or 2, characterized in that: The step of performing orthographic projection on the bone tissue data to obtain bone projection data includes: According to each pixel point in the axial image corresponding to the ray passing through the bone tissue data and the ray length weight corresponding to each pixel point, the CT value of each pixel point is weightedly superimposed to obtain the bone projection data; The ray length weight is negatively correlated with the length of the ray in the voxel.

4. The method according to claim 1 or 2, characterized in that: The step of obtaining a CT reconstructed image according to the bone removal tissue data and the back-projection image comprises: The ratio of the bone-removed tissue data to the water reference is used as the equivalent length of each projection pixel point in the bone-removed tissue data relative to water; Performing weighted correction on the bone projection data using the artifact correction coefficient corresponding to the equivalent length to obtain corrected bone projection data; After the corrected bone projection data is reconstructed, the reconstruction result is superimposed on the back-projection image to obtain a CT reconstructed image.

5. The method according to claim 1 or 2, characterized in that: The filtering of the projection data to obtain filtered data includes: The projection data are rearranged using a pre-established initial virtual coordinate system, wherein the initial virtual coordinate system takes the detector irradiated by the central radiation source as the origin, the spiral line as the horizontal axis, and the horizontal displacement bed direction as the vertical axis; The rearranged projection data is filtered according to the windowing filter kernel to obtain filtered data.

6. The method according to claim 5, characterized in that The method of rearranging the projection data by using a pre-established initial virtual coordinate system includes: Rearranging the projection data of the fan beam into the projection data of the parallel beam by rearranging the direction angle of the projection data; According to the up-sampled data amount and the interval of the direction angle rearrangement, the projection data of the parallel beam is subjected to three-fold Hermite interpolation to obtain a three-dimensional surface, and the three-dimensional surface is subjected to coordinate interpolation processing perpendicular to the initial virtual coordinate system to obtain rearranged projection data.

7. The method according to claim 6, characterized in that The performing coordinate interpolation processing perpendicular to the initial virtual coordinate system on the three-dimensional curved surface to obtain rearranged projection data includes: Determine the grid points in the three-dimensional surface and calculate the rotation matrix of each of the grid points; The rotation matrix is ​​used to transform the coordinates of each of the grid points from the initial virtual coordinate system to a target virtual coordinate system, wherein the horizontal axis of the target virtual coordinate system is horizontal; Interpolation processing is performed on each of the grid points in the target virtual coordinate system to obtain rearranged projection data.

8. The method according to claim 7, characterized in that The interpolation processing is performed on each of the grid points in the target virtual coordinate system to obtain rearranged projection data, including: Performing interpolation processing on each of the grid points in the target virtual coordinate system to obtain an interpolation result; The interpolation result is filtered according to the initial virtual coordinate system to obtain rearranged projection data.

9. The method according to claim 5, characterized in that The windowing filter includes an all-pass frequency and a cutoff frequency, and the windowing filter kernel is: Wherein, the w represents the frequency component of the signal, the wl represents the all-pass frequency, and the wh represents the cut-off frequency.

10. The method according to claim 4, characterized in that Before filtering the projection data to obtain filtered data, the method further includes: performing slope correction on the projection data; The step of superimposing the reconstruction result with the back-projection image to obtain a CT reconstructed image includes: Superimposing the reconstruction result with the back-projection image to obtain a superimposed result; The stacking result is intercept-corrected, and the correction result is slice-thickness processed to obtain a CT reconstructed image.

11. The method according to claim 1, characterized in that: The back-projection processing is performed on the filtered data to obtain a back-projection image, including: Determine the extension range according to the extension collimation width, the pitch and the reconstruction interval, and shift the starting position of the horizontal motor to the left and right based on the extension range to obtain a virtual horizontal motor position; Interpolation processing is performed based on the virtual horizontal motor position and the filtered data to obtain extended projection data, and weighted back projection processing is performed based on the extended projection data and the filtered data to obtain a back projection image.

12. The method according to claim 11, characterized in that The step of performing weighted back projection processing based on the extended projection data and the filtered data to obtain a back projection image comprises: Performing weighted back-projection processing based on the extended projection data and the filtered data to obtain an initial back-projection image; For any pixel point in the initial back-projection image, calculate the distance between the coordinate value of the pixel point on the depth axis and the position of the horizontal motor corresponding to the central ray passing through the pixel point, and determine that the pixel point is valid when the distance is less than a preset threshold value; According to the Q value weight curve, the rearranged data corresponding to the effective pixel points are processed by weighted summation to obtain the back-projection image.

13. A spiral CT reconstruction device, characterized in that: The device comprises: A filtering module, used for filtering the projection data to obtain filtered data; A back-projection module, used for performing back-projection processing on the filtered data to obtain a back-projection image, and extracting bone tissue data from the back-projection image; An orthographic projection module, used for performing orthographic projection on the bone tissue data to obtain bone projection data; A processing module, used for taking the difference between the projection data and the bone projection data as bone removal tissue data; A reconstruction module is used to obtain a CT reconstructed image based on the bone removal tissue data and the back-projection image.

14. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the steps of the method according to any one of claims 1 to 12 are implemented.

15. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 12 are implemented.

16. A computer program product comprising a computer program, characterized in that When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 12 are implemented.

Images