Efficient algebraic tomographic reconstruction of three-dimensional objects

Abstract

The invention relates to a method (100) for determining the volume (3b) of portions (3) of a voxel (2) in a sample (1), which is detected by a beam (4) penetrating the sample (1), as a subvoxel, wherein the beam (4) emanates from a starting point and has a conical geometry, the method having the steps: • a collection (5) of closed formulae (5a) for the calculation of volumes of geometric shapes is provided (110); • intersection points and / or intersection lines (4a) of the beam (4) with the boundaries of the voxel (2) are determined (120); • geometric shapes (3a) of the subvoxels (3) are determined (130) on the basis of these intersection points and / or intersection lines (4a); • for these shapes (3a), it is checked (140) whether the collection (5) contains a closed formula (5a) for calculating the volume (3b) of the respective shape (3a); • if the collection (5) contains a closed formula (5a) for calculating the volume (3b) of the respective shape (3a), this formula (5a) is used to determine (150) the sought volume (3b); • if the collection (5) does not contain such a closed formula (5a), the sought volume (3b) is determined (160) by determining the volumes (6b) of parts (6) of the voxel (2) not covered by the beam (4) and subtracting them from the total volume of the voxel (2).

Description

[0001] Efficient algebraic tomographic reconstruction of three-dimensional objects

[0002] The invention relates to the tomographic reconstruction of three-dimensional objects from a plurality of one- or two-dimensional projections obtained by radiating the object.

[0003] State of the art

[0004] In medicine, materials testing, and access control to sensitive security areas, samples—people, animals, workpieces, or pieces of luggage—are examined using tomography. For this purpose, penetrating radiation is directed through the sample from many different perspectives onto a one- or two-dimensional detector, creating a projection of the sample onto the detector. From a multitude of such projections, the three-dimensional internal structure of the sample is reconstructed.

[0005] The method of filtered back-propagation (FBP) is primarily used for this reconstruction. This method exploits the fact that the internal structure of the sample, on the one hand, and the set of all projections, on the other, are described by the Radon transform, and this transformation is bijective. Thus, there is an inverse Radon transform that leads from the projections back to the internal structure of the sample. The name "filtered back-projection" derives from the fact that the inverse Radon transform involves a convolution integral of the projections with a high-pass filter.

[0006] Alternatively, an algebraic reconstruction technique, ART, can be used. Here, the intensity measured by a detector element is considered the algebraic sum of the contributions of volume elements (voxels) of the sample arranged in a grid, which the radiation has penetrated on its path from the radiation source to the detector. Under the simplifying assumption that the exponential attenuation of the radiation as it passes through the sample can be linearized on the short path through the respective voxel, a linear equation of the form

[0007] Herein is p i the intensity value of the i-th projection. N is the total number of voxels. μ j is the attenuation coefficient of the j-th voxel, which contains the information about the material with which this voxel is filled. w ij - is the

[0008] Weighting coefficient of the i-th projection with respect to the j-th voxel, i.e. the volume fraction of the j-th voxel that is captured by the i-th beam. This results in a linear system of equations that is determined by the attenuation coefficients μ j can be resolved.

[0009] There is a total of M = N p · N d Projections, where N p is the number of all used orientations (e.g. angle steps) and N d is the number of all available detector elements. If M > N, the system of equations is overdetermined. It can theoretically be solved directly, provided that the product W T W of the weight matrix W with its transpose is not singular. In practice, however, iterative reconstruction techniques are typically used.

[0010] The use of ART has several advantages over FBP. For example, ART in medicine requires a single set of projections, which results in a lower total dose for the patient, and is less disrupted by metal implants. Other deviations from idealized assumptions, such as irregular scanner geometries, missing data, finite dimensions of detectors, and detectors with varying sensitivity, can also be better taken into account. This comes at the cost of the very large weight matrix W, also known as the system matrix. It is typically so large that it cannot be stored completely on the hardware platform used. Thus, whenever a weight w ij from this matrix W is needed, this weight w ij These calculations are therefore performed very frequently and are a determining factor in the total computing time required for reconstruction.

[0011] Task and solution

[0012] It is therefore the object of the present invention to provide a method with which the weights w ij for the ART can be calculated more quickly and thus the evaluation of projections using ART can be accelerated overall.

[0013] This object is achieved according to the invention by a method according to the main claim. Further advantageous embodiments emerge from the dependent claims.

[0014] Disclosure of the invention

[0015] The invention provides a method for determining the volume of portions of a voxel in a sample that are captured by a beam penetrating the sample. Portions of the voxel that are captured by the beam are also referred to below as subvoxels. The method assumes that the beam emanates from a starting point and has a conical geometry. Within the framework of the method, a collection of closed formulas for calculating volumes of geometric shapes is provided. These formulas can, in particular, refer to specific types of shapes and be parameterized with reference points or other characteristics of these shapes. For example, the formula for the volume of a polyhedron can be parameterized with coordinates of the vertices of this polyhedron. The formula for the volume of a cone, on the other hand, can be parameterized with the circle radius of the base and the height, for example.

[0016] Intersection points and / or intersection lines of the beam with the boundaries of the voxel are determined. Based on these intersection points and / or intersection lines, geometric shapes of the subvoxels are determined. A check is then carried out to determine whether the collection contains a closed formula for calculating the volume of the respective shape.

[0017] If the collection contains a closed formula for calculating the volume of the respective shape, this formula is used to determine the desired volume.

[0018] If, however, the collection does not contain such a closed-form formula, the desired volume is determined by determining the volumes of those parts of the voxel not covered by the beam and subtracting them from the total volume of the voxel. The idea behind this is that the volume of many seemingly complex shapes can be expressed as the difference between the volumes of two simple shapes. For example, the volume of a cylinder and the volume of a pyramid can both be easily calculated using closed-form formulas, and the volume of a cylinder missing a piece in the shape of the pyramid is the difference between the volume of the cylinder and the volume of the pyramid. However, it would not be so easy to assemble the cylinder with the missing pyramid additively from shapes whose volumes can be calculated using closed-form formulas.It was discovered that, when implemented on a computer, calculating the volume of a geometric shape using a closed-form formula specific to that shape is orders of magnitude faster than calculating that volume using a universal algorithm that is not limited to that shape. The difference in speed is so great that even calculating several parts of the voxel not covered by the beam is still significantly faster than calculating directly with the universal algorithm.

[0019] This is not an approximate calculation, as is often used to save computing time. It still achieves the same exact result as was achieved with previous methods over a significantly longer time.

[0020] A significant overall saving in computing time is achieved even if, in an individual case, a portion of a voxel captured by the beam cannot be sufficiently expressed additively and / or subtractively using shapes whose volume can be calculated using an available closed-form formula. If the calculation using closed-form formulas proposed here is successful for most calculations of volumes captured by the beam during a tomographic reconstruction, this saves so much computing time overall that it is not significant if this approach fails in a small proportion of cases, and in these cases, known but significantly slower algorithms must be used.The speed with which a tomographic reconstruction can be completed as a whole does not depend on whether the method proposed here is successful and faster in every individual case. Rather, it is crucial that the entire to-do list of all calculations of acquired volumes that occur during this reconstruction is processed more quickly overall.

[0021] This can be illustrated using a simple analogy. Imagine a quality control process for mass-produced workpieces where the workpieces were previously placed in an X-ray machine and the X-ray images were evaluated. A new method is now available in which the workpiece is stimulated to acoustic vibrations by a hammer blow, and the quality of the workpiece is determined from the spectrum of these vibrations. This only takes a few seconds, whereas the detour via the X-ray machine previously took several minutes. Most workpieces can be reliably classified using the new method, which saves a huge amount of time. For some workpieces, however, the analysis of the acoustic vibrations does not allow a clear conclusion to be drawn about the quality, so the X-ray machine must be used again. In these cases, the previous use of the new method only cost unnecessary additional time.However, this is hardly significant given the cumulative time savings achieved in the vast majority of cases.

[0022] In a particularly advantageous embodiment, the original voxel is divided into two parts. For each of these parts, detected subvoxels are determined. The total of the subvoxels determined in this way are combined to form the desired subvoxels of the original voxel. By subdividing, detected subvoxels can in particular be divided into two parts for which it is easier to find one or more closed formulas for calculating the volume. For example, if a detected subvoxel is a cuboid with another shape missing in its base and in its lid, this detected subvoxel can be divided down the middle into two cuboids by subdividing the original voxel, with one cuboid missing another shape in the top and the other cuboid missing another shape in the bottom.For these two parts of the subvoxel, the volume calculation using closed formulas, here: the difference between the volume of the cuboid and the volume of the other shape missing within it, is easier to find than for the subvoxel as a whole, i.e. for the cuboid with two missing shapes.

[0023] In another particularly advantageous embodiment, at least one portion of the voxel not covered by the beam is divided into at least two parts. At least one part is selected such that the collection provides a closed formula for calculating its volume. The volume of this at least one part is then determined using this closed formula. In this way, at least a first approximation for the largest part of the desired volume can be obtained very quickly. The remaining missing volume can then be calculated, for example, using a conventional algorithm. The runtime of such a conventional algorithm can depend on the size of the volume to be determined. Therefore, if the largest part of the volume has already been calculated using a closed formula, the remaining, more difficult but smaller part can be calculated using conventional methods more quickly than the volume of the detected subvoxel as a whole.

[0024] This principle can also be continued iteratively. Therefore, in another advantageous embodiment, a part determined during the decomposition, whose volume cannot be determined using a closed formula from the collection, is again divided into two subparts. At least one subpart is chosen such that the collection provides a closed formula for calculating its volume. The volume of this at least one subpart is then determined using this closed formula.

[0025] In another particularly advantageous embodiment, the collection contains at least one closed formula for calculating the volume of a polyhedron. These are the shapes that most frequently occur when rectilinearly propagating beams intersect with cuboid-shaped voxels. The formulas for calculating the volume of polyhedra are also relatively simple and can therefore be processed quickly on a computer. A polyhedron can, for example, be a pyramid or be created by extruding a triangle or trapezoid.

[0026] In a particularly advantageous embodiment, the intersection points and / or intersection lines of the beam with the boundaries of the voxel are determined using

[0027] • the relative position of the starting point to the voxel,

[0028] • a direction of propagation of the beam in one or more angles starting from the starting point and

[0029] • an aperture angle of the conical geometry of the beam bundle. These are relatively simple geometric calculations that can also be expressed in closed-form formulas.

[0030] In another particularly advantageous embodiment, at least one volume of a geometric shape calculated using a closed-loop formula is stored in a memory. When the volume of a similar geometric shape needs to be calculated, this volume is reused. In this way, the recalculation of a volume using a complicated formula, the evaluation of which is "expensive" in terms of computing time, can be replaced by simple memory access.

[0031] As previously explained, the ultimate goal of determining the acquired subvoxels is to accelerate the tomographic reconstruction of the sample's internal structure. Therefore, in another particularly advantageous embodiment, volumes of the portions of all voxels captured by a beam penetrating the sample are determined. An equation is then established that relates a sum of voxel attenuation coefficients weighted by these volumes to the intensity of the beam measured by a detector after penetrating the sample. This forms the basis for faster tomographic reconstruction using algebraic reconstruction techniques (ART).

[0032] For this purpose, a large number of equations for beams penetrating the sample via different paths can be combined into a system of equations. The system of equations can then be solved for the attenuation coefficients of the voxels. As previously explained, these attenuation coefficients provide information about the material contained in each voxel of the sample. In particular, the reconstruction reveals a contrast between different materials, which can be visualized as an image of the sample's internal structure.

[0033] In a particularly advantageous embodiment, at least some of the equations are set up for beams that penetrate the sample on paths that are rotated and / or translationally offset from one another. In this way, it can be ensured that the sets of voxels that contribute to the individual equations at least partially overlap, so that a statement can be made for each individual voxel from the interaction of several equations. If, however, for example, a small subset of the equations only contains contributions from a few voxels and these voxels do not contribute to the remaining equations of the system of equations, only this small subset of the equations can be used to determine the desired attenuation coefficients for these voxels. The entire large remainder of the system of equations cannot help here. Ideally, the beams penetrate the sample in a spiral manner, i.e.Both mutually rotated and offset paths. This gives rise to the term "trauma spiral" for a tomographic whole-body examination, which is used primarily in emergency situations, particularly in the form of an X-ray CT scan, due to the high radiation exposure. In another particularly advantageous embodiment, the system of equations is solved using an iterative numerical solution method. This is different from, for example, the direct solution of the matrix equation.

[0034] W · μ = p according to the vector p of the attenuation coefficients according to μ = (W T W)- 1 · W T · p with the Moore-Penrose pseudoinverse, the process of iterative solving can be better controlled and, if necessary, continuously monitored, so that, for example, artifacts that can arise from the amplification of the noise always contained in real measurement data can be better suppressed.

[0035] In another particularly advantageous embodiment, a detector with a flat or annular surface sensitive to the beam is selected. Calculating which portions of a sample voxel are captured by the portion of a beam with a conical geometry that later reaches the respective detector element leads to particularly simple closed formulas for these detector geometries.

[0036] In another particularly advantageous embodiment, the sample is exposed to the beam(s) from a radiation source. The radiation source can, in particular, be an X-ray source, for example. If, for example, electrons are accelerated towards a target in an X-ray tube using high voltage, X-rays are emitted in the conical geometry assumed in the method proposed here. The intensity of the beam after penetrating the sample is registered with at least one detector. This completes the chain of events from the acquisition of the raw data on the sample to the determination of the subvoxels of each voxel covered by the radiation up to the tomographic reconstruction. The method can, in particular, be computer-implemented and thus embodied in the form of software.The invention therefore also relates to a computer program containing machine-readable instructions which, when executed on one or more computers and / or compute instances, cause the computer(s) and / or compute instances to execute the method described above. In this context, control units of X-ray machines and / or tomographs, as well as other embedded systems capable of executing machine-readable instructions, are also considered computers. Compute instances are virtual machines, containers, or other execution environments via which computing power for the execution of machine-readable instructions is provided in a public or private cloud.

[0037] The invention also relates to a machine-readable data carrier and / or a downloadable product containing the computer program. A downloadable product is a digital good that contains the computer program and can be offered for sale, in particular, in an online shop for immediate delivery via download. Furthermore, one or more computers and / or computing instances can be equipped with the computer program, the machine-readable data carrier, and / or the downloadable product.

[0038] Special description section

[0039] The subject matter of the invention is explained below with reference to figures, without limiting the subject matter of the invention. It is shown:

[0040] Figure 1: Embodiment of the method 100 for determining the volume 3b of portions 3 of a voxel 2 captured by a beam 4; Figure 2: Illustration of the use of the collection 5 of closed formulas 5a within the framework of the method 100 using simple examples;

[0041] Figure 3: Schematic subdivision of a pixel by the rays g1, g2, g3 and g4 through the edges V1, V2, V3 and V4 respectively;

[0042] Figure 4: different cases of the contribution of a subvoxel 3 to a detector element;

[0043] Figure 5: Determination of the volume of subvoxels contributing to a detector element by appropriate segmentation;

[0044] Figure 6: Integration of the weighting coefficients for the respective detector elements D1 to D6 for a single voxel 2;

[0045] Figure 7: Consideration of blind spots B1, B2 between detector elements D1, D2 and D3 during the evaluation.

[0046] Figure 1 is a schematic flow diagram of an exemplary embodiment of the method 100 for determining the volume 3b of portions 3 of a voxel 2 in a sample 1, which is captured by a beam 4 penetrating the sample 1, as a subvoxel. The beam 4 emanates from a starting point and has a conical geometry.

[0047] In step 110, a collection of 5 closed-form formulas 5a for calculating the volumes of geometric shapes is provided. According to block 111, this collection 5 may contain at least one closed-form formula 5a for calculating the volume of a polyhedron.

[0048] In step 120, intersection points and / or intersection lines 4a of the beam 4 with the boundaries of the voxel 2 are determined.

[0049] According to block 121, the intersection points and / or intersection lines 4a of the beam 4 with the boundaries of the voxel 2 can be determined using

[0050] • the relative position of the starting point to voxel 2,

[0051] • a direction of propagation of the beam 4 given at one or more angles from the starting point and

[0052] • an opening angle of the conical geometry of the beam bundle 4 can be determined.

[0053] In step 130, geometric shapes 3a of the subvoxels 3 are determined based on these intersection points and / or intersection lines 4a.

[0054] For this purpose, the original voxel 2 can be divided into two parts, for example, according to block 131. Subvoxels 3' detected for each of these parts can then be determined according to block 132. The subvoxels 3' thus determined can be combined according to block 133 to form the desired subvoxels 3 of the original voxel 2.

[0055] In step 140, a check is performed for the previously determined shapes 3a to determine whether the collection 5 contains a closed formula 5a for calculating the volume 3b of the respective shape 3a. If this is the case (truth value 1), the desired volume 3b is determined using this formula 5a in step 150. If, however, no closed formula 5a exists for the shape 3a (truth value 0), the desired volume 3b is determined in step 160 by determining the volumes 6b of portions 6 of the voxel 2 not captured by the beam 4 and subtracting them from the total volume of the voxel 2.

[0056] According to block 161, at least one portion 6 of the voxel not captured by the beam 4 can be divided into at least two parts. At least one part is selected such that the collection 5 provides a closed formula 5a for calculating its volume 6b. The volume 6b of this at least one part can then be determined using this closed formula 5a according to block 162.

[0057] This can be further iteratively nested by subdividing a part determined during the decomposition, whose volume 6b cannot be determined using any closed formula 5a from the collection 5, into two subparts according to block 161a. At least one subpart is selected such that the collection 5 provides a closed formula 5a for calculating its volume. According to block 161b, the volume of this at least one subpart can then be determined using this closed formula 5a.

[0058] In step 170a, at least one volume of a geometric shape calculated using a closed-form formula 5a can be stored in a memory. It can then be reused in step 170b whenever the volume of a similar geometric shape needs to be calculated.

[0059] In the example shown in Figure 1, in step 180, volumes 3b of the portions 3 of all voxels 2 that are detected 4 by a beam penetrating the sample 1 are determined. Then, in step 190, at least one equation 7 is established that relates a sum of attenuation coefficients 2a of voxels 2, weighted with these volumes 3b, to an intensity 4b of the beam 4 measured by a detector after penetrating the sample 1. According to block 191, a plurality of equations 7 for beams 4 that penetrate the sample 1 via different paths can be combined into a system of equations. This system of equations is solved for the attenuation coefficients 2a of the voxels 2 according to block 192.

[0060] In this case, in particular, for example according to 191a, at least part of the equations 7 can be set up for beams 4 which penetrate the sample 1 on paths which are rotated and / or translationally offset relative to one another.

[0061] According to block 193, the system of equations can be solved using an iterative numerical solution method.

[0062] According to block 194, a detector with a flat or annular surface sensitive to the beam can be selected.

[0063] In the example shown in Figure 1, in step 105, sample 1 is exposed to the beam(s) from a radiation source. In step 106, the intensity 4b of the beam 4 is registered with at least one detector after penetrating sample 1. This step can also be performed, for example, after the volumes 3b have been determined and before the system of equations 7 has been solved. This means that the irradiation of sample 1 can follow a precise plan for which the essential preparatory work for evaluation has already been performed.

[0064] Figure 2 uses two simple (freely chosen) examples to illustrate how the collection of 5 closed-form formulas 5a can contribute to the faster calculation of the volumes 3b of subvoxels 3. Intersection points and / or intersection lines 4a of the beam 4 with the voxel 2 are shown in each case. In the example shown in Figure 2a, the portion 3 of the voxel 2 that is captured by the beam 4 is located in the upper left corner of the front side of the voxel 2 facing the observer. The geometric shape 3a of this portion 3 is a triangular cylinder, i.e., a triangle extruded into the plane of the drawing. For calculating the volume 3b of this geometric shape 3a, the collection of 5 contains a suitable closed-form formula 5a.

[0065] Conversely, in the example shown in Figure 2b, a central region 3 of voxel 2, bounded by two uncaptured triangular cylinders 6, 6a, is captured by the beam 4. Collection 5 does not contain a suitable closed-form formula 5a for the shape 3a of this central region 3. Collection 5, however, contains a suitable closed-form formula 5a for calculating the volumes 6b of the triangular cylinders 6a. By subtracting these volumes 6b from the total volume of voxel 2, the desired volume 3a of the central region 3 can be determined.

[0066] The following provides an exemplary mathematical derivation of closed formulas for calculating the volumes of geometric shapes that most frequently occur as subvoxels captured by the beam bundle or as parts of a voxel not captured by the beam bundle.

[0067] This derivation assumes that the intensity value p ithe i-th projection, the attenuation coefficient μ j of the j-th voxel and the weighting coefficient w ij the i-th projection with respect to the j-th voxel via (1) are linked to each other. Each weighting coefficient is the ratio of the sub-voxel volume S covered by the beam j of voxel j to the total volume V j this voxel j:

[0068] Derivation for a flat panel detector

[0069] The Cartesian system used has its origin at the center of the three-dimensional voxel grid. As is common in practically implemented reconstruction methods, the X-ray source S is assumed to be a point with the coordinate s = (5, 0, 0), where S is the distance from the grid center. In the case of a flat-panel detector, the square detector units are defined by the sides d y and d zin the y- and z-direction, respectively. The pixel response of the detector is assumed to be uniform. The distance of the detector from the origin of the coordinate system on the opposite side of the source S in the x-axis is given as D. Then, a vector d can be assigned to a point at the edge of a detector element as d = (-D, m y d y ,m z d z ) T be assigned, where m y and m z are integers. A trajectory from S to the detector point d can be described by the vector d s be displayed.

[0070] When the X-ray source and detector are rotated by the angle φ counterclockwise around the z-axis (as in computer tomography of the human body), which corresponds to a clockwise rotation of the sample in a fixed tube-detector system, the s be written as follows where λ is a real number. Equation (3) can also be written as The vector d s and a vector from S to the center of the detector thus include the

[0071] Angle β m a.

[0072] From this we can conclude that for any point (x,y,z) along a ray path with angle β the following relations are valid:

[0073] Equations (11) to (14) apply to any beam path with angle β.

[0074] Two-dimensional case: XY plane

[0075]

[0076] The angle β i of the beam path from the source S to the i-th intersection point of the

[0077] Detector elements can be calculated from:

[0078] The integers m yi can be derived from the following relationships: where [x] integeris an integer immediately less than x. It can be assumed that For this purpose, the fourfold symmetry of the rectangular or square pixels with respect to the 360° circumference is used: In the second and fourth cases, sides a and b must be swapped.

[0079] To calculate the pixel contributions to each detector element, which is generated by a

[0080] beam can be hit, the relationships between the angles β i and β gs can be determined. The resulting area elements can be calculated directly or by combining two or three area elements positively or negatively. The formulas for two area elements are sufficient: a triangle and a trapezoid. The triangle is bounded by a corner of the pixel and a ray to a detector intersection point. The trapezoid is bounded by two trajectories that intersect two opposite pixel sides.

[0081] Triangle in the XY plane

[0082] Figure 3 shows, by way of example, the schematic subdivision of a pixel by the rays g1, g2, g3, and g4 through the edges V1, V2, V3, and V4, respectively. (a) Rays through Vi and within g2 bound a triangle (b) The rays between g2 and g4 define a trapezoid. The lines b1 to b 5 represent trajectories to the edges of the detector elements, which delimit subvoxels contributing to the respective detector elements.

[0083] In the following derivation it is assumed that the triangle is defined by the

[0084] Corner V1of the pixel and the intersection points of the i-th ray path with the enclosing sides of this corner, which are called P x and P y can be assigned (see Figure 3a). Using equations (11) and (12), the corresponding vectors can be written as follows:

[0085] Using equations (22) and (20), the two sides of the resulting right-angled triangle can be calculated as follows: The resulting area A tri can therefore be described as follows:

[0086] Since the area of ​​the rectangular pixel is ab, the relative area factor is:

[0087] Trapezoid in the XY plane

[0088] When the i-th and the (il)-th ray paths intersect two opposite sides of the pixel, a trapezoid is formed. The intersection points can be denoted as P i , P i+1 , Q i and Q i+1 (see Figure 3b). Since the height of the trapezoid is a, the area

[0089] Using equations (21) and (22), equation (32) can be written as follows:

[0090] Area factors

[0091] When determining the weighting coefficients, the relationships between the angles β i and β g s are taken into account. For this purpose, in a first step, the angles β g s are sorted by their magnitudes. Two cases can be distinguished: β g > β 2g and β g < β 2g This can be achieved by introducing the angles β G = min (β g β 2g ) and β 2G = max (βg , β 2g ) can be achieved. Analogously, β 3G = min(βs g , β 4g ) and β 4G = max (β 3 g , β4 g ). Note that for square pixels The lines spanned by the angles divide the pixel into three parts: two triangles, the enclose a trapezoid. The rays with angles β i can be assigned to these three subcells according to their magnitudes, so that for the first triangle and for the second triangle

[0092] Using the above formulas, the area factors f i (i = 1, 2, ... 6) in Figure 3c as follows: with f3= f 31 + f 32 and f5= f 51 + f 52This example clearly demonstrates the general algorithm for generating the two-dimensional weighting factors. The matrix of the third group can be derived from the transposed matrix of the first two equations. In three dimensions, this structural algorithm remains in the xy plane.

[0093] Three-dimensional case

[0094] Figure 4 shows different cases of the contribution of a subvoxel 3 to a detector element. Rays to the edges of the detector element are shown in dashed lines, (a) subvoxel VPQW with the coefficient vector C m = (1, 0, 0, 0). (b) Subvoxel V1V2W2Q2Q1W1with the coefficient vector C m = (1, 1, 0, 0). The volume is the difference between the volumes V1PQ1W1 and V2PW2Q2. (c) Subvoxel with the coefficient vector C m = (1, 1, 1, 1). (d)

[0095] Subvoxel with the coefficient vector C m = (l, 1, 1).

[0096] For the following analysis, an arbitrarily chosen orthorhombic voxel V with the side lengths a, b and c can have the edges are integers. Then it is necessary to determine the z-position of the detector elements involved in addition to the y-component. From equation (7) we can derive where z max and z min can be chosen as the upper or lower limit of the voxel element. In the case of V 10 and V 01 , where β = β g , taking into account equation (21), the following relationships apply:

[0097] First, we consider the case where only the intersection point V = ((n x + 1)a, y, z) T of the i-th beam has a contribution from the z-component. The resulting volume Vol pyrcan be calculated by determining the points where the beam paths to the upper boundary of the detector meet its surface (see Figure 4a). The vectors to these points P and Q can be written as

[0098] The volume Vol pyr of the resulting pyramidal element can be calculated as follows

[0099] Using equations (41) - (43), equation (44) can be written as follows: and tan If only the intersection point of the i-th ray has a contribution from the z-component, the resulting volume Vol pyr be calculated as follows If the intersection points of the i-th ray path and V 2 of the (i+1)-th ray path have contributions due to the z-component (see Figure 4b), the resulting volume can be considered a composite of the individual pyramid volumes of both trajectories. It can be determined by subtracting the smaller volume from the larger volume.

[0100] Therefore, the volume Vol can be calculated as follows:

[0101] In summary, the relative volume (in relation to the orthorhombic voxel volume Vol voxel = abc) Vol rel express as follows: where the components of the vector are either 1 or 0 are when an intersection occurs or not. For simplicity, only the variables that change within the voxel are listed in the Vol function. In Figures 4a and 4b, (1, 0, 0, 0) T or (1, 0, 0, 1) T . The The resulting volume elements can be calculated directly or by the negative combination of two volume elements:

[0102] In a next step, the second fundamental volume element, which in the two-dimensional case corresponds to the triangle, is further analyzed. If all vertices contribute to the volume Vol (see Figure 4d), one obtains: Equations (58) and (61) correspond to equations (28) and (29), which describe the

[0103] represent the sides of the two-dimensional triangle. After some algebraic work and using equation (65), we get:

[0104] In general, the relative volume (relative to the orthorhombic voxel volume Vol voxel = abc) Vol rel be expressed as follows: where either 1 or 0 if an overlap takes place or not. The resulting volume elements can be calculated directly or by the negative combination of two volume elements:

[0105] Figure 5 illustrates the developed algorithm.

[0106] To calculate the relative volumes of the individual voxel elements, the cell is first divided into two halves by a plane E that intersects the front and back planes. Then, the volumes from above and below are calculated. For the upper and lower elements E1 and E2, the vectors C m (1, 1, 0, 1) T or (0, 0, 1, 0) T To obtain the volume elements of the next cells, they are separated from the subsequent cells using the vectors C m = (1, 1, 1, 1) T subtracted that reach the upper and lower surfaces. Finally, the volumes of the divided middle cell must be added.

[0107] Ring-shaped static detector In the case of a closed, curved, i.e. ring-shaped detector as in CT of the human body, a vector d can be assigned to a point at the edge of a detector element as follows

[0108] (70) where m y and m z are integers, d z is the length of the detector element in the z-direction and Ω is the aperture angle of a detector element. If only the X-ray source is rotated counterclockwise around the z-axis by the angle φ, the s be written as follows

[0109] (71) where A is a real number. Equation (71) can also be written as

[0110] Therefore, for any point (x,y,z) along a ray with angle β, the following relations hold: Here, equations (78) to (80) apply to any ray path. Let z = 0. For the following analysis, it is convenient to use the angles and β 4g the four edges where a and b are the sides of the pixel.

[0111] From equation (78) the following relationships can be derived:

[0112] From equations (81) and (82) the following expression can be derived:

[0113] Combining equations (82) and (81) with equation (86), we obtain:

[0114] The angle β i of the beam from the source S to the i-th intersection point of the detector elements can be written as follows:

[0115] The following relationships can be derived from equation (89):

[0116] Using equation (90), the following relationship can be derived from equation (74): The integers m yi can be derived from the following relationships:

[0117] Two-dimensional case: XY plane

[0118] Triangle in XY plane

[0119] In the following derivation, it is assumed that the triangle is bounded by the vertex Vi of the pixel and the intersection points of the i-th ray with the enclosing sides of this vertex, which are denoted as P x and P y can be assigned (see Figure 3a). Using equations (79) and (80), the corresponding vectors can be written as follows:

[0120] Using equations (86) and (88), the two sides of the resulting right-angled triangle can be calculated as follows:

[0121] The resulting area A trican therefore be described as follows: Since the area of ​​the rectangular pixel is ab, the relative area factor is:

[0122] Trapezoid in the XY plane

[0123] When the i-th and (i+1 )-th rays intersect two opposite sides of the pixel, a trapezoid is formed. The intersection points can be denoted as P i , P i+1 , Q i and Q i+1 (see Figure 3b). Since the height of the trapezoid is a, the area A tra be expressed as follows Then the relative area factor is:

[0124] Using equations (87) and (88), equation (102) can be written as follows:

[0125] Three-dimensional case

[0126] Using equation (92), equation (75) can be rewritten as follows: where z max and z mincan be chosen as the upper or lower limit of the voxel element. In the case of V 11 and V 01 , where β = β g For example, the following relationships apply, taking into account equation (87): Using equation (105), equations (76) and (77) can be rewritten as follows:

[0127] Further algebraic work leads to the following relationships (equations (111) and (117)), which are analogous to equations (54) and (67) for the flat detector:

[0128] Figure 6 shows the integration of the weighting coefficients for the respective detector elements D1 to D6 for a single voxel 2 during rotation (a) and rotational and translational motion (b). The X-ray beams to the corners of the detector elements D1 to D6 are shown in dashed lines.

[0129] This means that with the algorithm presented here, it is possible to incorporate the data acquired during the rotation and longitudinal translation of the beam 4 relative to sample 1 into the tomographic reconstruction practically seamlessly and without interpolation. Figure 7 shows how blind spots B1, B2 between detector elements D1, D2, and D3 can be taken into account during the evaluation: The calculation of the subvoxels 3, which are acquired by the beam 4 and contribute to the intensity on the detector elements D1, D2, and D3, respectively, can then be carried out based only on the active parts of these detector elements D1, D2, and D3. Therefore, if portions of the beam that pass through certain areas B1 and B2 of voxel 2 do not land in the active area of ​​a detector element D1, D2, or D3 and are therefore not registered by any of the detector elements D1, D2, and D3, these areas B1 and B2 can be excluded from the outset as "blind spots."

[0130] Matlab code

[0131] The following is a Matlab code with which an embodiment of the invention was realized:

[0132]

[0133] Conclusion

[0134] An efficient and easy-to-implement algorithm for cone-beam CT has been developed to decompose the voxels into their respective sub-voxels contributing to the corresponding detector elements. Explicit mathematical expressions have been derived to calculate the exact volumes. Two mathematical functions, i.e., equations (54) & (67) for flat-panel detectors and (111) & (117) for ring-shaped detectors, are sufficient to determine the weighting coefficients of all contributing subvoxels. Data acquisition during rotation and longitudinal translation, as in spiral CT techniques, can be incorporated without interpolation. By properly choosing the voxel parameters, i.e., the angles in the xy-plane and the z-axis, the algorithm can be easily extended to account for system-inherent conditions, such as blind spots between the detector elements and collimators in CT scanners for the human body.

[0135] List of reference symbols

[0136] 1 sample

[0137] 2 voxels

[0138] 2a Weighting factor of voxel 2

[0139] 3, 3' captured portion, subvoxel, of a voxel 2

[0140] 3 a geometric shape of subvoxel 3

[0141] 3b Volume of subvoxel 3

[0142] 4 beam bundles

[0143] 4a Intersection points / intersection lines of beam 4 with voxel 2

[0144] 4b intensity of beam 4 recorded with detector

[0145] 5 Collection of closed formulas 5a

[0146] 5a closed formulas for volumes 3b of shapes 3a

[0147] 6 areas of voxel 2 not covered by beam 4

[0148] 6a Forms of areas 6

[0149] 6b Volumes of areas 6

[0150] 7 equations

[0151] 100 Procedure for determining the volumes 3b of shares 3a

[0152] 110 Providing collection 5 of formulas 5a

[0153] 111 Choice of a collection 5 with polyhedra

[0154] 120 Determining intersection points and / or intersection lines 4a

[0155] 121 specific methods for calculating intersection pointsZ-lines 4a

[0156] 130 Determining geometric shapes 3a from subvoxels 3

[0157] 131 Dividing the original voxel 2 into two parts

[0158] 132 Determine detected subvoxels 3' for the two parts of voxel 2

[0159] 133 Union of subvoxel 3' to subvoxels 3 of voxel 2

[0160] 140 Check if Collection 5 contains Formula 5a for Form 3a

[0161] 150 Calculating the volume 3b using formula 5a

[0162] 160 Subtractive determination of volume 3b using proportions 6 161 Dividing an unrecorded area 6 into two parts 161a Dividing a part into two subparts without formula 5a 161b Calculating the volume of a subpart using formula 5a

[0163] 162 Calculating the volume 6b of a part using formula 5a 170a Storing calculated volumes in memory

[0164] 170b Reusing calculated volumes

[0165] 180 Determining the volumes 3 for all parts 3 of voxels 2

[0166] 190 Setting up at least one equation 7

[0167] 191 Combining many equations 7 to form a system of equations 191a Setting up equations 7 for rotated / offset bundles 4

[0168] 192 Solving the system of equations for attenuation coefficients 2a

[0169] 193 Solving with iterative numerical solution methods

[0170] 193 Choosing a flat or ring-shaped detector

[0171] B1, B2 blind spots D1-D6 Detector elements

[0172] E1, E2 Volume element

Claims

Claims 1. Method (100) for determining the volume (3b) of portions (3) of a voxel (2) in a sample (1) which is detected by a beam (4) penetrating the sample (1) as a subvoxel, wherein the beam (4) emanates from a starting point and has a conical geometry, comprising the steps: • a collection (5) of closed formulas (5a) for calculating volumes of geometric shapes is provided (110); • intersection points and / or intersection lines (4a) of the beam (4) with the boundaries of the voxel (2) are determined (120); • based on these intersection points and / or intersection lines (4a), geometric shapes (3a) of the subvoxels (3) are determined (130); • for these shapes (3a) it is checked (140) whether the collection (5) contains a closed formula (5a) for calculating the volume (3b) of the respective shape (3a); • if the collection (5) contains a closed formula (5a) for calculating the volume (3b) of the respective shape (3a), the desired volume (3b) is determined using this formula (5a) (150); • if the collection (5) does not contain such a closed formula (5a), the desired volume (3b) is determined (160) by determining the volumes (6b) of parts (6) of the voxel (2) not covered by the beam (4) and subtracting them from the total volume of the voxel (2).

2. Method (100) according to claim 1, wherein the original voxel (2) is divided into two parts (131), for each of these parts respectively detected subvoxels (3') are determined (132) and the subvoxels (3') thus determined are combined to form the desired subvoxels (3) of the original voxel (2) (133).

3. Method (100) according to one of claims 1 to 2, wherein • at least one portion (6) of the voxel not captured by the beam (4) is divided (161) into at least two parts, wherein at least one part is selected such that the collection (5) provides a closed formula (5a) for calculating its volume (6b); and • the volume (6b) of said at least one part is determined using said closed formula (5a) (162).

4. The method (100) according to claim 3, wherein • a part determined during the decomposition, the volume (6b) of which cannot be determined using a closed formula (5a) from the collection (5), is again divided into two sub-parts (161a), wherein at least one sub-part is selected such that the collection (5) provides a closed formula (5a) for calculating its volume; and • the volume of said at least one sub-part is determined using this closed formula (5a) (161b).

5. Method (100) according to one of claims 1 to 4, wherein the collection (5) contains at least one closed formula (5a) for calculating the volume of a polyhedron (111).

6. Method according to one of claims 1 to 5, wherein the intersection points and / or intersection lines (4a) of the beam (4) with the boundaries of the voxel (2) are determined based on the relative arrangement of the starting point to the voxel (2), • a direction of propagation of the beam (4) given at one or more angles from the starting point and • an opening angle of the conical geometry of the beam bundle (4) can be determined (121).

7. Method (100) according to one of claims 1 to 6, wherein at least one volume of a geometric shape calculated using a closed formula (5a) is stored in a memory (170a) and, when the calculation of the volume of a same geometric shape is pending, is reused (170b).

8. Method (100) according to one of claims 1 to 7, wherein • Volumes (3b) of the portions (3) of all voxels (2) that are captured (4) by a beam penetrating the sample (1) are determined (180) and • at least one equation (7) is established (190) which relates a sum of attenuation coefficients (2a) of voxels (2) weighted with these volumes (3b) to an intensity (4b) of the beam (4) measured with a detector after penetrating the sample (1).

9. The method (100) according to claim 8, wherein • a multitude of equations (7) for beams (4) which penetrate the sample (1) in different ways are combined into a system of equations (191) and • the system of equations is solved for the attenuation coefficients (2a) of the voxels (2) (192).

10. The method (100) according to claim 9, wherein at least part of the equations (7) is set up (191a) for beams (4) which penetrate the sample (1) on paths which are rotated and / or translationally offset relative to one another.

11. The method (100) according to any one of claims 9 to 10, wherein the system of equations is solved using an iterative numerical solution method (193).

12. Method (100) according to one of claims 8 to 11, wherein a detector with a flat or annular surface sensitive to the beam is selected (194).

13. Method (100) according to one of claims 1 to 12, further comprising: applying (105) the sample (1) with the beam(s) (4) from a radiation source and registering (106) the intensity (4b) of the beam (4) after penetrating the sample (1) with at least one detector.

14. A computer program containing machine-readable instructions which, when executed on one or more computers and / or compute instances, cause the computer(s) and / or compute instances to carry out the method (100) according to any one of claims 1 to 12.

15. Machine-readable data carrier and / or download product with the computer program according to claim 14.

16. One or more computers and / or compute instances with the computer program according to claim 14, and / or with the machine-readable data carrier and / or download product according to claim 15.

Images