Fan-beam CT (Computed Tomography) rapid algebraic reconstruction method and system based on area score

Abstract

The invention belongs to the technical field of medical CT (Computed Tomography) detection, and discloses a fan-beam CT rapid algebraic reconstruction method based on an area score, which is characterized by comprising the following steps: under a planar fan-beam model, taking an area as a weight factor, designing a scheme capable of rapidly and accurately calculating the weight factor, and taking an ART algorithm as a reconstruction method. And reconstructing to obtain a high-quality reconstructed image. The method is a reconstruction method with the area as a weight factor under an ART algorithm framework, the problem of insufficient precision caused by an inaccurate model is avoided, secondly, due to the fact that a calculation mode of recursively calculating the intersection area is adopted, whether intersection exists or not does not need to be determined through pixel-by-pixel traversal, and the calculation efficiency is improved. Area information of all intersected pixels can be obtained line by line (column) only by knowing an entrance point and an exit point where a beam intersects with a reconstruction region, theoretically, the reconstruction speed of the method is obviously superior to that of a traversal method, and due to the fact that no approximate method is adopted in the calculation process, a final image has the reconstruction quality the same as that of the traversal method.

Description

Technical Field

[0001] The present invention belongs to, but is not limited to, the technical field of medical CT detection, and particularly relates to a fast algebraic reconstruction method and system for fan-beam CT based on area integration. Background Art

[0002] CT (computerized tomography), that is, computerized tomography technology, different internal structures and materials of an object will have different effects on the X-rays passing through. For example, in medicine, different human internal tissues have different absorption coefficients for X-rays. Therefore, according to the attenuation of X-rays after passing through an object, the absorption coefficients of the structures on the corresponding paths can be constructed, and the internal structure information of the object can be obtained based on the position information and the magnitude of the absorption coefficients. This is extremely helpful for obtaining the internal structure of a closed object, so it has important applications in medical disease diagnosis, industrial non-destructive testing, and cultural relic protection and reconstruction.

[0003] The ART (algebraic reconstruction) algorithm is based on such a principle that the reconstruction area is discretized into an n*n square grid, and it is stipulated that the row and column numbers start from 0. The grid at the 0th row and 0th column is recorded as number 0, then the number at the nth row and nth column is n 2 -1. Taking j as the number, the value f j of each grid is the average value of the jth pixel. According to the physical process of the beam passing through the reconstruction area, it can be discretized into algebraic equations by integration, and the relationship between the projection data of the beam passing through the reconstruction area and the pixel values of the corresponding area can be represented by an algebraic equation system:

[0004]

[0005] where N = n*n, p i represents the projection value of the ith ray, ω ij represents the contribution of the jth pixel to the projection of the ith ray. There are M rays in total. Therefore, this formula represents the contribution of all pixels in the measurement reconstruction area to a certain ray, and the projection value of the ray passing through the reconstruction area can be obtained. From this, a rectangular matrix can be obtained, and by solving the equation system, the pixel values of the reconstruction area can be obtained. However, due to the large size of the matrix and the existence of contradictory equations in the equation system, it is impossible to solve it by the conventional method of finding the inverse matrix. Usually, the following iterative formula is adopted:

[0006]

[0007] Among them, \(i\) represents the iteration number, and \(\lambda\) is the relaxation factor, whose value is generally between 0 and 2. From the construction of the iterative formula, it can be seen that the \(w\) system matrix is repeatedly used during the iteration process. Therefore, the calculation and access of the system matrix are the main challenges in the iteration process. Some researchers use a simple idea, that is, pre-calculate the system matrix and store it in the corresponding storage device, and then access the storage device to obtain the specific value when calculating. However, this method has two main defects: First, the system matrix is generally very large, and without using compression storage technology, ordinary computers cannot meet the storage requirements; Second, frequent reading and writing of the hard disk will limit the reconstruction speed due to the hard disk read and write speed. Therefore, a more efficient method is to calculate the required system matrix during the line-by-line correction process. The calculation of the system matrix \(w\) will directly affect the operation efficiency of the ART algorithm.

[0008] The selection of the weighted factor calculation model is an important factor affecting the reconstruction result and the consumption of iteration time. The weighted factors can generally be divided into three categories: one is the pixel model, the second is the line integral model, and the third is the area integral model. The pixel model defines the pixel value at the pixel center. Only when the ray passes through the pixel center will the pixel be counted, and its weight value is defined as 1, otherwise it is defined as 0; The line integral model first assumes that the ray width is 0 and defines the pixel weight as the length of the ray passing through the pixel; The area integral model assumes that the beam has a width, and its contribution to a pixel it passes through is directly related to its area in the pixel. Obviously, the area integral model is more in line with the physical reality. In terms of the complexity of calculation, the pixel model is the simplest, and the area integral model is the most complex because it involves the calculation of the area of irregular polygons and the filling inside the beam. In terms of reconstruction quality, there are a large number of artifacts in the pixel model, and the line integral model has a good balance between speed and quality, so it has been studied more deeply, but the area integral model is the most in line with the actual situation.

[0009] The light source beam model is also divided into two types: one is the parallel beam, and the other is the fan beam. The difference lies in the different ways of obtaining projection data. The parallel beam is realized by the translation and rotation of the point ray source and the point detector. The ray utilization rate is low, and the acquisition speed is slow. Also, because it involves two operations of translation and rotation, the frequent change of the motion mode limits the scanning speed and is only used in the first-generation CT. In view of the shortcomings of the parallel beam model, the fan beam model is proposed. The core of this model lies in two points: single X-ray and multiple detectors, and the continuous rotation of the detection rack to avoid mechanical motion changes. It is the current mainstream method.

[0010] After the reconstruction is completed, a quantitative evaluation of the reconstruction process is carried out. The reconstruction speed is obtained according to the average running time of one or more iterations. For the reconstruction quality, we use two indicators for evaluation, the Normalized Root Mean Square (NRMS) and the Normalized Mean Absolute (NMA). The calculation expressions of NRMS and NMA are as follows:

[0011]

[0012] where \(t\) u，v , \(r\) u，v represent the pixel values of the \(u\)-th row and \(j\)-th column in the original image and the reconstructed image respectively. refers to the average value of all pixel values in the source image. NRMS can more sensitively reflect some small error situations that many points have, while NMS reflects the situation where relatively large errors occur at some points. The smaller these two judgment bases are, the smaller the difference between the reconstructed image and the original image, and the better the reconstruction quality.

[0013] According to the above-mentioned iterative formula of ART, various improved or fast algorithms have been proposed. Mesquita et al. can obtain images with higher reconstruction quality in fewer iterations by improving the selection of the relaxation factor; Wang et al. change the selection order of the projection sequence in the reconstruction process and adopt the method combining WDS (Weighted Distance Scheme) and orthogonality to obtain the weighted distance orthogonality method, thus achieving a faster convergence rate and better reconstruction quality; Yang proposed a method for quickly calculating pixel indices and penetration lengths recursively under the framework of the distance model, eliminating the sorting process in the Siddon algorithm and obtaining higher execution efficiency; in terms of the calculation of the area weight factor, the commonly used calculation scheme is to determine the endpoints of the intersecting pixels, then use the Sutherland-Hodgman clipping algorithm to calculate the endpoints of the intersecting polygon and sort the endpoints, and finally use the polygon area calculation formula in the form of endpoints to obtain the area of the polygon formed by the beam intersecting the pixel. Sungsoo established a lookup table for the area obtained by the Sutherland-Hodgman algorithm, replacing a large number of calculation processes with a small amount of memory operations and obtaining the accurate area by interpolation. Further, a regression model was established based on the lookup table to approximately calculate the intersecting area, improving the reconstruction speed with less quality loss; Qiao proposed a scheme that does not require calculating the system matrix under the framework of the TV-minimization algorithm. Based on the idea that the rotation of the detector and the rotation of the object to be detected are equivalent, only the rotated image needs to be iterated each time, and the rotation scheme is carried out by interpolation. Different interpolation methods will significantly affect the quality and speed of the final reconstructed image; I.K. Hong et al. made full use of various symmetries to reduce the data to be stored, reducing the amount of computational tasks to 1 / 8 of the original.

[0014] In the algorithms proposed above, the research on the improvement of the relaxation factor and the selection of the projection order cannot improve the running efficiency of the program in one iteration. Therefore, the problem of time consumption of the area model cannot be fundamentally solved. The S-H algorithm used to calculate the area has a high time complexity. Multiple calls will make the iteration time of the ART algorithm too long to meet the requirements of the reconstruction speed. The methods of rotating the image and establishing the lookup table or regression formula both adopt the interpolation method. Although this approximate method can simplify the calculation and improve the efficiency, the final reconstruction quality depends on the accuracy of the approximation. In scenarios that require higher accuracy, more complex interpolation methods are needed, which leads to an increase in the computational amount offsetting the speed advantage brought by the approximation. The method using symmetry needs to store the information of the symmetric reference system matrix in the memory. As the image expands, the memory capacity is not enough to support the simultaneous storage of a large amount of data, resulting in performance degradation. Secondly, in the symmetric model, it is assumed that the entire CT system is strictly symmetric. However, in actual situations, due to mechanical errors, the positions of the detectors are not strictly symmetric, and the reconstructed image obtained by symmetry will have a gap with the real image. Summary of the Invention

[0015] Aiming at the problems existing in the prior art, the present invention provides a fast algebraic reconstruction method for fan-beam CT based on area integration.

[0016] The present invention is implemented as follows. A fast algebraic reconstruction method for fan-beam CT based on area integration, characterized in that the method includes, under the planar fan-beam model, taking the area as a weight factor, designing a scheme that can quickly and accurately calculate the weight factor, and using the algebraic (ART) algorithm as the reconstruction method to obtain a high-quality reconstruction image through the intersection of the beam and the pixel.

[0017] Furthermore, the method specifically includes:

[0018] A beam is composed of two homologous rays. The intersection of the beam and the reconstruction area can be regarded as the intersection of the two rays and the area between the two rays and the reconstruction area. Therefore, the problem of the intersection of the beam and the reconstruction area can be transformed into the problem of the intersection of the straight line and the reconstruction area and the filling of the area inside the two rays. The internal filling can be performed according to the boundary conditions after the beam edge is determined. Here, according to the slope of the ray, the intersection is divided into three cases: the slope |k|≥1, |k|<1, and two special cases, namely horizontal and vertical. The horizontal and vertical cases are relatively simple, and the index and weight of the row or column can be determined only according to the positions of the exit point and the entry point. When k is the other two cases, the operations are consistent. Therefore, only one of the cases needs to be considered. Here, taking |k|<1 as an example, at this time, the ray steps faster along the x direction, and at this time, row-by-row processing is selected. Conversely, column-by-column processing is performed along the y direction.

[0019] A line segment AB, where A and B are both endpoints of pixel grids, and the side length of each pixel square is denoted as 1. Now, it is required to find the area of the polygon formed by AB and each pixel it passes through. Let the slope of the line be k. Notice that the areas of their intersections are all triangles or right trapezoids. Therefore, one method is to find the equation of the line where AB lies, intersect it with the lines x = a, where a = 0, 1, 2, 3, 4, 5 respectively, obtain all the intersection points, and then use the triangle or trapezoid area formula to find the intersection area of each pixel respectively. However, this method requires obtaining all the intersection points to get the corresponding length information and using the area formula, which consumes a large amount of resources. Consider the method of translation. If the line is translated one unit length to the right, notice that each indexed area has the same difference from the previous area, which is the size of a parallelogram. From the slope relationship, the area of this quadrilateral is a constant M = k. For example, if the area to be calculated now is S n , and the previous area is S n-1 , then there is always:

[0020] S n = S n-1 + M

[0021] Among them, the area of the first triangle can be directly obtained from the slope, and the areas of all subsequent intersection regions can be obtained from this recurrence relation. Considering that the exit point may not be at the endpoint of the pixel grid, the calculation of the last intersecting pixel of the reference triangle can be carried out in a way of decreasing area, that is, the total area is gradually subtracted by the intersection areas of each region. The above is the calculation method of the reference triangle. Here, mainly a simple one-step addition is used to replace the steps of obtaining intersection points and using the area formula. Extend the line segment in the direction of point B. Except for the change of pixel index, the second calculation can be regarded as the translation process of the reference triangle obtained in the first calculation. Under continuous translation, the loop length of each row is Then process each row one by one until the end point of the ray in the reconstruction region is reached, and enter the next line segment to be processed. This means that only the intersection situation of the first reference triangle needs to be obtained, and the intersection situations of triangles in all subsequent rows can be obtained by translation means. According to this method, the unilateral area of the pixels where all line segments intersect with the line can be obtained faster;

[0022] Considering the general situation of translation, it is divided into two cases according to whether point A is at the endpoint of the pixel block: When the starting point A is at the endpoint of the pixel block, the result is the same as the calculation of the reference triangle at this time, and only the new pixel number needs to be obtained; When point A is not at the endpoint of the pixel block, as shown in Figure 3 , this situation can be regarded as being translated a distance m from the region formed by parallel lines starting from the origin as the reference. The number of pixels passed through by the ray and the areas formed in each pixel have changed due to translation. Except for the first and the last pixels, the other pixels all differ from those before translation by a constant area M 1= m * |k|. The area of the first pixel can be calculated jointly by the translation distance m and the slope k. Whether in the calculation of the reference triangle or the translated triangle, the area of the last pixel can be obtained by decreasing the area of the total area in that row. Since only a translation motion is performed, the total area of each row is a constant value. The increase in the number of pixels is only determined by the translation distance m relative to the reference triangle and the position of point B. Using the ceiling function ceil(), if ceil(B.x) - B.x > m, then this translation does not change the number of pixels passed through. Otherwise, the number of pixels passed through increases by one relative to the reference state.

[0023] Therefore, if we look at the intersection of the ray within the reconstruction area row by row, it can be found that the triangles formed by the ray and the upper and lower boundaries of each row are congruent. The intersection situation of each row can be regarded as a reference triangle starting from the origin being translated. Therefore, only by recursively obtaining the information of the reference triangle and performing translation can we obtain the unilateral area information of the intersection of the ray and the reconstruction area. There will be pixels within the beam that do not intersect with either of the two boundary rays. At this time, it can be regarded as an internal filling process when the boundary lines are known. After that, it can be quickly obtained which pixels in the reconstruction area the beam passes through and the intersection area value.

[0024] Combined with the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solution to be protected by the present invention are as follows:

[0025] First, the effect of the present invention is verified using the Sheep-Logan head simulation model. Corresponding code is written using Visual Studio 2022 on the windows11 platform and compared and verified with the following two different methods.

[0026] One is a relatively basic algorithm, roughly as follows. When calculating the situation of each beam passing through the reconstruction area, all pixels within the reconstruction area are traversed, and each pixel is judged whether it intersects with the beam. If it intersects, the Sutherland-Hodgman clipping algorithm is used to obtain all the intersection points of the triangular beam and the pixel block boundary, and the polygon area formula is used to calculate the area of the intersection region of the two figures. Each ray is corrected in this way, which is called the traversal clipping method. The other is the method for the line integral model, which is used to verify the advantage of this area integral method in the quality of the reconstructed image. Here, the traditional Siddon algorithm is adopted.

[0027] In the traversal clipping algorithm, the calculation of each ray requires traversing all the pixels in the reconstruction area to determine whether it intersects with the beam. However, a single beam intersects only a very limited number of pixels in the reconstruction area, which obviously causes waste. Moreover, the method of obtaining the intersecting polygon using the Sutherland-Hodgman clipping algorithm requires finding all the endpoints of the polygon, resulting in a long calculation time. For the Siddon algorithm used in the line integral model, in terms of model selection, the beam is only regarded as a ray with zero width, ignoring the true physical characteristics of the beam, reducing the number of weight factors required for iteration and the calculation amount of a single weight factor. This causes the inaccuracy of the model, resulting in more artifacts in the finally reconstructed image when the number of iterations is small. The newly proposed recursive rule does not have the above problems. First, this method is a reconstruction method with area integral as the weight factor under the ART algorithm framework, avoiding the problem of insufficient accuracy caused by inaccurate models. Second, due to the use of a recursive calculation method for the intersecting area, there is no need to traverse each pixel to determine whether there is an intersection. Only by knowing the entry point and exit point where the beam intersects the reconstruction area can the area information of all intersecting pixels be obtained row by row (column). Theoretically, the reconstruction speed of this method will be significantly better than the traversal method, and since no approximation method is used in the calculation process, the final image will have the same reconstruction quality as the traversal method. The results are as follows:

[0028] Figure 1 They are respectively the reconstruction results obtained by iterating three times in three ways: the line integral Siddon method, the area integral traversal clipping method, and the area integral recursive method. It can be seen that compared with the latter two area integral methods, there are still relatively serious artifacts in the reconstructed image of the line integral Siddon method after three iterations, and the reconstruction effects of the latter two methods are basically the same. From this, it can be shown that: the new area integral recursive method is superior to the reconstruction method with line integral as the weight factor in terms of reconstruction quality.

[0029] In terms of speed, when using the same computing platform (Intel Core i7 8750H, 8G of memory) and code structure, the time consumption of each method is shown in Table 1. The new recursive method is significantly better than the traversal clipping method and also better than the Siddon method of the line integral model, greatly improving the computing power of the area integral model and having a certain speed competitiveness.

[0030] Table 1

[0031]

[0032]

[0033] Second, the expected benefits and commercial value after the transformation of the technical solution of the present invention are as follows: This algorithm significantly improves the speed and quality of the reconstructed image generation, making the internal structure of the object clearer, and the dynamic changes of the internal structure can also be obtained through rapid reconstruction, which is of great significance for medical diagnosis. This fast algorithm creates the possibility for the imaging of dynamic scanning of the internal structure of the object and is expected to become an important technology in the process of the next generation of CT image reconstruction.

[0034] The technical solution of the present invention overcomes the technical prejudice: In recent years, due to the wide application of machine learning and deep learning technologies, the field of CT image reconstruction has increasingly focused on the above two technologies. However, deep learning and machine learning technologies also have certain deficiencies, such as the lack of a strict mathematical proof process and the reconstruction results being greatly affected by the training data set. Currently, most of the actually applied algorithms are FBP algorithms. However, this algorithm has high requirements for the acquired data and is greatly affected by noise, resulting in a hidden risk of high radiation dose. The algebraic reconstruction algorithm does not have the above defects, but has always been considered a method with extremely high computational consumption. Especially, the area integral model has great defects in terms of time consumption and cannot be applied in actual detection. However, the rapid algorithm, acceleration technology, and iterative upgrade of computer hardware have greatly shortened the time consumption of this method, making it have practical application value, proving the effectiveness of the iterative algorithm in actual application and having the value of further development. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 It is the reconstruction result obtained by projecting and iterating three times in the angular order by the three methods of line integral Siddon method, area integral traversal and clipping method, and area integral recurrence method provided by the embodiment of the present invention;

[0036] Figure 2 It is a schematic diagram of the area of the polygon formed by AB and each pixel it passes through provided by the embodiment of the present invention;

[0037] Figure 3 It is a schematic diagram showing that point A is not at the end point of the pixel block provided by the embodiment of the present invention;

[0038] Figure 4 It is a statistical chart of the running time of the embodiment of the present invention. The running time takes the running time of the core function of different methods in a single iteration and the total running time of a single iteration as the comparison standard, and the unit is seconds;

[0039] Figure 5 It is a statistical chart for comparing the quality of the reconstructed images of the embodiment of the present invention. The above-mentioned NRMS and NMA values are used as criteria, and the values after projecting and iterating six times in the angular order are adopted here. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0040] In order to make the objectives, technical solutions and advantages of the present invention more clear and understandable, the present invention will be further described in detail below in conjunction with embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0041] Application Example 1: Medical Tomography (CT) Image Reconstruction

[0042] The fast algebraic reconstruction method for fan-beam CT based on area integration of the present invention can be widely applied to medical CT imaging. In practical applications, traditional reconstruction algorithms are often limited by high computational complexity and low reconstruction efficiency when dealing with tomographic images with complex tissue structures. This method significantly improves the efficiency and accuracy of CT image reconstruction by using area as a weighting factor and a recurrence relation to quickly calculate the intersection area between the beam and the pixel. For example, when scanning the thoracic or abdominal region, this method can quickly generate high-resolution images, accurately reconstruct complex anatomical structures, and help doctors detect lesions earlier and make diagnoses.

[0043] Application Example 2: Industrial Non-Destructive Testing

[0044] This method can be used in non-destructive testing (NDT) systems in the industrial field, especially for the detection of internal defects in mechanical components. For example, when detecting whether there are cracks or pores inside an aero-engine blade or a complex metal part, an industrial CT device can use the method of the present invention to quickly reconstruct the scanned data. Traditional algebraic reconstruction methods are slow when dealing with a large number of detection tasks, while the present invention significantly reduces the calculation time by efficiently calculating the intersection area between the beam and the pixel. Specifically, this method can quickly generate high-precision internal structure images, enabling inspectors to identify material defects in a timely manner, improving the detection efficiency and reducing production risks.

[0045] An embodiment of the present invention provides a fast algebraic reconstruction method for fan-beam CT based on area integration, characterized in that the method includes, under a planar fan-beam model, using area as a weighting factor, designing a scheme that can quickly and accurately calculate the weighting factor, and using the ART algorithm as the reconstruction method to reconstruct a high-quality reconstructed image.

[0046] The method specifically includes:

[0047] Recursive establishment of the reference triangle and calculation of the area of the translated triangle

[0048] A beam is composed of two homologous rays. The intersection of the beam and the reconstruction area can be regarded as the intersection of the two rays and the area between the two rays and the reconstruction area. Therefore, the problem of the intersection of the beam and the reconstruction area can be transformed into the problem of the intersection of a straight line and the reconstruction area and the filling of the area within the two rays. The internal filling can be performed according to the boundary conditions after the beam edge is determined. Here, according to the slope of the ray, the intersection is divided into three cases: the slope |k| >= 1, |k| < 1, and two special cases, namely horizontal and vertical. The horizontal and vertical cases are relatively simple, and the index and weight of the row or column can be determined only according to the positions of the exit point and the entry point. When k is in the other two cases, the operations are consistent. Therefore, only one of the cases needs to be considered. Here, |k| < 1 is taken as an example. At this time, the ray steps faster in the x direction, and the rows are processed one by one. Conversely, the columns are processed one by one along the y direction.

[0049] As Figure 2 shown, there is a line segment AB, where both A and B are the endpoints of pixel grids. The side length of each pixel square is denoted as 1. Now, it is required to obtain the area of the polygon formed by AB and each pixel it passes through. Let the slope of the straight line be k. Note that the areas of their intersections are all triangles or right trapezoids. Therefore, one method is to obtain the equation of the straight line where AB is located, intersect with the straight lines x = a, a = 0, 1, 2, 3, 4, 5 respectively, obtain all the intersection points, and then use the triangle or trapezoid area formula to obtain the intersection area of each pixel respectively. However, this method requires obtaining all the intersection points to obtain the corresponding length information and using the area formula, which consumes a large amount. Consider the method of translation. If the straight line is translated one unit length to the right, note that the area of each index area has the same difference from the previous area, and they are all the size of a parallelogram. From the slope relationship, the area of this quadrilateral is a constant M = k. For example, if the area to be calculated now is S n , and the previous area is S n-1 , then there always exists:

[0050] S n = S n-1 + M

[0051] where the area of the first triangle can be directly obtained from the slope, and the areas of all subsequent intersection regions can be obtained from this recurrence relation. Considering that the exit point may not be at the endpoint of the pixel grid, the calculation of the last intersection pixel of the reference triangle can be performed in the way of decreasing area, that is, the total area is gradually subtracted by the intersection areas of each region. The above is the calculation method of the reference triangle. Here, mainly a simple one-step addition is used to replace the steps of obtaining the intersection points and using the area formula. The line segment is extended in the direction of point B. Except that the pixel index changes, the second calculation can be regarded as the translation process of the reference triangle obtained in the first calculation. Under continuous translation, the loop length of each row is Then process line by line until reaching the end point of the ray in the reconstruction area, and move on to the next line segment to be processed. This means that only by obtaining the intersection situation of the first reference triangle can the intersection situations of the triangles in all subsequent lines be obtained through translation. According to this method, the unilateral area of the pixels where all line segments intersect with the straight line can be obtained more quickly;

[0052] Considering the general situation, it is divided into two cases according to whether point A is at the end point of the pixel block: when the starting point A is at the end point of the pixel block, the calculation result is the same as that of the reference triangle, and only the new pixel number needs to be obtained; when point A is not at the end point of the pixel block, as Figure 3 shown, this situation can be regarded as being translated by a distance m with the area formed by the parallel straight lines starting from the origin as the reference. The number of pixels passed through by the ray and the area formed in each pixel change due to translation. Except for the first and the last pixels, the other pixels differ from those before translation by a constant area M 1 = m * |k|, and the area of the first pixel can be calculated jointly by the translation distance m and the slope k Whether in the calculation of the reference triangle or the translated triangle, the area of the last pixel can be obtained by subtracting the area of the last pixel from the total area in this row. Since only a translation movement is made, the total area of each row is a constant value, and the increase in the number of pixels is only determined by the translation distance m and the position of point B. Using the ceiling function ceil(), if ceil(B.x) - B.x > m, then this translation does not change the number of pixels passed through, otherwise, the number of pixels passed through increases by one compared to the reference state;

[0053] If looking at the intersection situation of the ray in the reconstruction area in terms of rows, it can be found that the triangles formed by the ray and the upper and lower boundaries of each row are congruent. The intersection situation of each row can be regarded as being translated from a reference triangle with the starting point at the origin. Therefore, only by obtaining the information of the reference triangle through recursion and performing translation can the unilateral area information of the intersection of the ray and the reconstruction area be obtained; there will be pixels in the beam that do not intersect with either of the two boundary rays. At this time, it can be regarded as an internal filling process under the condition of known boundary straight lines. After that, it can quickly obtain the intersection area values of the pixels passed through by the beam in the reconstruction area.

[0054] Example 1: Reconstructing medical CT images using area integration

[0055] For a certain medical scenario, it is necessary to perform a CT scan on a patient's chest and use a fan-beam CT system for fast algebraic reconstruction in order to obtain high-precision tomographic images for doctors to diagnose lung diseases.

[0056] 1. Establishment of the beam model

[0057] Using a planar fan-beam model, the beam emitted by the X-ray source is defined as the region between two homologous rays.

[0058] Take the intersection area of each beam with the pixels in the reconstruction region as the weight factor.

[0059] 2. Weight factor calculation

[0060] Classify according to the ray slope, and set the cases where the slope |k| > 1, |k| < 1, and horizontal and vertical rays.

[0061] Taking |k| < 1 as an example, process row by row in the x direction, and use the area recurrence method to quickly calculate the intersection area of each pixel passed by the ray:

[0062] First, calculate the area information of the reference triangle under this k value.

[0063] Then, calculate the area of the triangle at the initial intersection of the ray and the pixel, calculate the area of the constant parallelogram through the length of the straight line translated to the right from the reference, and quickly obtain the intersection area and position number of the subsequent pixels through the difference in area.

[0064] 3. Tail adjustment

[0065] When the starting point of the ray is not at the end point of the pixel block, use the total area decreasing method to quickly calculate the area of the end pixel.

[0066] 4. Reconstructed image

[0067] Input the calculated weight factor into the ART algorithm to perform algebraic reconstruction on the scanned data and obtain a clear chest CT image.

[0068] Example 2: Efficient image reconstruction in industrial inspection

[0069] In industrial non-destructive testing, it is necessary to perform fan-beam CT scanning on a large mechanical part to detect internal cracks or defects. Due to the large size and complex structure of the part, a fast algebraic reconstruction method is required to improve the detection efficiency.

[0070] 1. Scanning model establishment

[0071] Establish a fan-beam CT scanning model for the mechanical part, and take the intersection area of the fan-beam ray and the detection region as the weight factor.

[0072] According to the size of the part, divide it into multiple pixel units, and set the side length of each pixel to 1.

[0073] 2. Fast weight factor calculation

[0074] For the case where the slope |k| > 1, process column by column in the y direction, recursively obtain the reference triangle, and gradually calculate the intersection area of the ray and the pixel.

[0075] Using the translation recurrence method, quickly obtain the intersection area and position number of adjacent pixels, avoiding the steps of calculating intersection points and area formulas one by one.

[0076] When the starting point of the ray is not at the non-pixel end point, the total area decreasing method is used to quickly correct the intersection area.

[0077] 3. Interior filling optimization

[0078] For pixels where the beam boundaries do not directly intersect, use the interior filling algorithm to supplement the full authority factor based on the beam boundary conditions and recurrence relations.

[0079] 4. Image reconstruction and analysis

[0080] Input the weight factor and scan data into the ART algorithm for algebraic reconstruction to generate high-resolution CT slice images.

[0081] Analyze the reconstructed image to locate the crack position and size inside the mechanical part.

[0082] As Figure 1 shown, it is the comparison result of the simulation experiment of the 256*256-sized Sheep-Logan head model. Among them, the Siddon algorithm is a classic method of algebraic reconstruction of the line integral model. This method uses the length of the line passing through the pixel as the weight factor. Compared with the area integral, the overall calculation amount is less and the running speed is faster. However, in the projection reconstruction of sequential angles, this method will produce artifacts, affecting the final imaging result. The area integral model effectively avoids this problem, such as the traversal clipping method. However, due to the higher computational complexity of area calculation, the time consumption is too long. This recurrence method can quickly calculate the pixel numbers and intersection areas passed by the beam, and does not lose the accuracy on the reconstructed image. It is also superior to the two comparison methods in terms of reconstruction speed.

[0083] Specifically, in terms of running speed, as Figure 4 shown, the running time of the core function and the total running time of this method in a single iteration are both less than those of the Siddon method; in terms of running quality, as Figure 5 shown, using the NMA and NRMS values obtained after six iterations of the two methods as the evaluation criteria, the smaller these two values are, the better the reconstruction quality and the closer to the original image. It can be seen that the two criteria of the new method are both better than the traditional Siddon algorithm. This shows that this new method has certain advantages in both the speed and quality of the reconstructed image.

[0084] The operation of the system provided by the embodiments of the present invention starts from the data acquisition module. This module acquires projection data through a fan-beam CT device and converts it into a matrix format suitable for calculation. The acquired data contains the scanning results at different projection angles, providing complete information on the beam passing through the reconstruction area. In this stage, preprocessing of the acquired data is also performed, including denoising and normalization operations, to improve the accuracy of subsequent calculations.

[0085] The weight factor calculation module quickly calculates the intersection area between the beam and the pixel according to the planar fan-beam model by the area recurrence method. First, the slope classification unit divides the rays into four cases according to the slope: |k|>1, |k|<1, horizontal, and vertical, and processes the intersection relationships in different cases. For the rays with |k|<1, the area recurrence unit quickly obtains the intersection area of the reference triangle and its subsequent regions by the translation method, avoiding calculating the intersection points pixel by pixel. The endpoint adjustment unit corrects the boundary conditions according to the pixel where the starting point or the ending point of the ray is located to ensure the calculation accuracy of the weight factor.

[0086] The reconstruction algorithm module processes the projection data by combining the algebraic reconstruction technique (ART algorithm). The beam filling unit fills the pixel regions that do not directly intersect within the beam boundary through geometric constraint conditions, thereby complementing the missing weight factor information. The pixel weight optimization unit adjusts the pixel weights of the projection data according to the weight factors, optimizes the gray value distribution of each pixel in the reconstructed image, and improves the detail resolution and overall quality of the image.

[0087] The high-quality reconstructed image is displayed in real time through the display module for users to analyze and diagnose. The display module supports multiple display modes, including grayscale images, pseudo-color enhanced images, and magnified display of specific regions, to meet different application requirements. In addition, users can also adjust the display parameters (such as contrast and brightness) of the image through the interactive interface of the module to further improve the visualization effect. The entire process realizes the efficient processing of fan-beam CT scan data and the rapid generation of high-quality images.

[0088] It should be noted that the embodiments of the present invention can be implemented by hardware, software, or a combination of software and hardware. The hardware part can be implemented using dedicated logic; the software part can be stored in a memory and executed by an appropriate instruction execution system, such as a microprocessor or dedicated designed hardware. Those of ordinary skill in the art can understand that the above devices and methods can be implemented using computer-executable instructions and / or included in processor control code, for example, such code is provided on a carrier medium such as a disk, CD, or DVD-ROM, a programmable memory such as a read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. The devices and modules of the present invention can be implemented by hardware circuits of programmable hardware devices such as very large scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, etc., or field programmable gate arrays, programmable logic devices, etc., can also be implemented by software executed by various types of processors, or can be implemented by a combination of the above hardware circuits and software, such as firmware.

[0089] As described above, the above are only specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention, any modifications, equivalent replacements, and improvements made within the spirit and principle of the present invention shall all be covered by the protection scope of the present invention.

Claims

1. A fast algebraic reconstruction method for fan-beam CT based on surface integral, characterized in that: Under the planar fan beam model, the area is used as a weight factor to calculate the intersection area of ​​the beam and the reconstruction area; the following steps are included: a) The intersection situations are divided into three types according to the slope of the ray: when the absolute value of the slope is greater than or equal to 1, the intersection area is calculated column by column; when the absolute value of the slope is less than 1, the intersection area is calculated row by row; for horizontal or vertical rays, the intersection area is determined directly according to the pixel start and end positions; b) For row-by-row or column-by-column calculations, obtain the one-sided area by following these steps: The intersection area of ​​the initial pixel is calculated by the starting position and slope of the ray entering the pixel, and the area is determined by the geometric relationship of the triangle or trapezoid formed by the ray passing through the pixel; For the intersection area of ​​subsequent pixels, it is recursively calculated by adding a fixed increment, which is the area of ​​the parallelogram formed by the ray and the pixel boundary; The area of ​​the last pixel is obtained by subtracting the previously calculated area from the total intersection area of ​​the row or column; c) After the calculation of the beam boundary is completed, the area surrounded by the boundary ray is filled with calculations to obtain the complete area distribution of the beam passing through the reconstruction area.

2. The method according to claim 1, characterized in that The intersection area of ​​the beam with the reconstruction region is calculated as follows: a) The initial area of ​​the reference pixel is determined by the ray slope and the starting point position, where the area is equal to the area of ​​the triangle or trapezoid formed within the pixel; b) The subsequent intersection area is calculated by recursion, where the area of ​​the current pixel is equal to the area of ​​the previous pixel plus a fixed increment, which is equal to the area of ​​the parallelogram corresponding to the slope of the ray; c) When the starting point of the ray is not on the pixel boundary, the area of ​​the pixel where the starting point is located is calculated by the relative position of the starting point and the pixel boundary, and the subsequent area is adjusted by the translation formula; d) After the intersection area of ​​the beam and the pixel boundary is calculated, the area of ​​the inner region is determined by the boundary conditions and filled to obtain a complete reconstruction result.

3. The method for fast algebraic reconstruction of fan-beam CT based on surface integral according to claim 1, characterized in that: The method specifically includes: The beam is decomposed into two homologous rays. The classification method based on the slope of the line is used to divide the ray slope |k| into four cases: greater than or equal to 1, less than 1, and horizontal and vertical. The intersection area of ​​the ray and the pixel is calculated step by step. In the horizontal and vertical case, the weight is directly calculated according to the position of the ray entry and exit points. In other cases, the weight is calculated step by step in the x direction row by row or in the y direction column by column.

4. The method for fast algebraic reconstruction of fan-beam CT based on surface integral according to claim 3, characterized in that: By translating the straight line equation of the ray to the right by a unit length, the area recursion relationship is used to quickly calculate the intersection area of ​​each pixel in the reference triangle. In each recursive process, the area change is equal to the fixed size area of ​​the parallelogram determined by the slope k, avoiding the calculation of all intersection points and the repeated use of area formulas.

5. The method for fast algebraic reconstruction of fan-beam CT based on surface integral according to claim 4, characterized in that: For the case where the starting point is at the end point of the pixel block, it is calculated in the following way: When the starting point is at the pixel endpoint, this is consistent with the intersection of the reference triangle. It is only necessary to change the pixel number and use the total area decreasing method to obtain the area of ​​the pixel at the end point. When the starting point is not at the pixel endpoint, the area of ​​the first pixel is calculated by combining the translation distance m and the slope k relative to the reference triangle, and the number of pixels and the area are adjusted recursively according to m.

6. The method for fast algebraic reconstruction of fan-beam CT based on surface integral according to claim 5, characterized in that: The intersection of rays in each row is regarded as a translation of the reference triangle starting at the origin. For the internal area calculation of the intersecting pixels in each row, a row-by-row recursive method is used to quickly obtain the intersection area of ​​each pixel, avoiding direct intersection calculation of each pixel.

7. The method for fast algebraic reconstruction of fan-beam CT based on surface integral according to any one of claims 1 to 6, characterized in that: For pixels inside the beam that do not intersect with the boundary rays, an internal filling method based on boundary conditions is used to determine their pixel numbers, and algorithm optimization is used to quickly calculate the intersection area of ​​all pixels passed by the beam, thereby improving the accuracy and speed of reconstructed images.

8. A fast algebraic reconstruction system for fan-beam CT based on surface integral, characterized in that: The system includes: A data acquisition module, used to obtain projection data of fan-beam CT scanning; The weight factor calculation module uses the area recursion method to calculate the intersection area between the beam and the pixels in the reconstruction area based on the planar fan beam model; The reconstruction algorithm module reconstructs the image of the projection data based on the algebraic reconstruction technology (ART algorithm) combined with the weight factor; The display module is used to display the reconstructed high-quality image.

9. The surface integral-based fan-beam CT fast algebraic reconstruction system according to claim 1, characterized in that: The weight factor calculation module further comprises: The slope classification unit is used to classify the rays into four cases: |k|>1, |k|<1, horizontal and vertical according to the ray slope |k|; An area recursive unit is used to calculate the area change of adjacent pixels by translation method to avoid pixel-by-pixel intersection calculation; The endpoint adjustment unit is used to correct the final weight factor according to the pixel position of the starting point or end point of the ray.

10. The surface integral-based fan-beam CT fast algebraic reconstruction system according to claim 1, characterized in that: The reconstruction algorithm module further comprises: A beam filling unit, for internally filling pixel areas that do not intersect with the boundary based on beam boundary conditions; The pixel weight optimization unit is used to adjust the pixel weight of the reconstructed image according to the weight factor and the projection data to improve the image reconstruction quality and accuracy.

Images