Dark-field and phase-contrast tomography reconstruction algorithms

Abstract

A signal processing system (SPS) and associated methods for tomographic reconstruction of phase contrast or dark field images. An improved weighting model is used to eliminate sensitivity gradient variations to reduce artifacts in the reconstructed images. The system (SPS) does not rely on an iterative reconstruction algorithm. Instead, faster reconstruction algorithms such as filtered backprojection can be used herein.

Description

[Technical field]

[0001] The present invention relates to a system for tomographic reconstruction of phase contrast or dark field images, an imaging apparatus, associated methods, computer program elements, and computer readable media. [Background technology]

[0002] One of the most useful tools in medicine is medical imaging, which allows for the non-invasive acquisition of medical images of a patient's internal structures, organs, and tissues, which can aid in diagnosis and inform treatment strategies.

[0003] One such type of medical imaging is X-ray-based imaging, in which the region of interest ("ROI") to be imaged is exposed to X-rays, detected by a detector, and then converted into an image that is then displayed on a monitor to aid clinicians in treatment or diagnosis, for example.

[0004] X-ray imaging involves radiography, where projection images are acquired and displayed. While this is useful for some applications, imaging based solely on projections has limitations due to occlusion by intervening structures. For example, the image representation of a lung nodule may be lost due to the presence of a rib in the projection path.

[0005] To overcome the limitations of imaging based solely on projections, tomographic imaging has been proposed. In tomographic imaging, cross-sectional images of a region of interest are acquired. In this type of imaging modality, multiple projection images are acquired from different directions around the ROI, and then the projection images are computationally processed to obtain the cross-sectional images. For this purpose, reconstruction algorithms are used.

[0006] Traditional X-ray imaging focuses on attenuation as the main contrast mechanism. However, this has limitations, for example when acquiring images of soft tissues. To overcome the limitations of imaging based on attenuation alone, new X-ray modalities are being investigated, including phase contrast and / or dark field imaging. In both cases, contrast mechanisms other than attenuation are exploited. When X-rays interact with matter, they are not only attenuated, but also affected by small-angle scattering and refraction. The contrast mechanisms of dark field and phase contrast imaging focus on these two types of contrast mechanisms. That is, in dark field imaging, the image contrast is given by the strength with which tissue causes small-angle scattering, while in phase contrast imaging, the contrast is given by the amount by which the X-ray wave is refracted. Attenuation is a measure of the attenuation coefficient, while the quantities reconstructed in dark field or phase contrast tomography imaging are the diffusion coefficient and the refractive index.

[0007] Dark field, or small-angle scattering, is distinguished from other larger-scale scattering effects, in which the scattering experienced by an X-ray photon is small enough that it can still be recorded by a pixel that intersects with the ray that the original incoming X-ray photon traveled before the (small-angle) scattering event.

[0008] Tomographic X-ray, i.e., phase contrast imaging, can suffer from image artifacts. Applicant's International Patent Publication WO2013 / 171657 (hereinafter "WO'657") proposes taking into account magnification effects to reduce the artifacts. Summary of the Invention [Problem to be solved by the invention]

[0009] However, despite these efforts, it has been found that artifacts, especially of the cupping type, sometimes persist.

[0010] Therefore, there is a need for improved phase contrast and dark field tomography. [Means for solving the problem]

[0011] The object of the invention is achieved by the subject matter of the independent claims, further embodiments are incorporated in the dependent claims. It is noted that the aspects of the invention described below apply equally to the imaging apparatus, the associated methods, computer program elements and computer readable media.

[0012] According to a first aspect of the invention, there is provided a system for tomographic reconstruction of phase contrast or dark field images, the system comprising: receiving phase contrast (PC) or dark field (DF) projection data based on measurements along rays through an image domain, the projection data being acquired in a scanning operation by an imaging device; performing data processing operations based on the projection data, the data processing operations including a reconstruction operation for reconstructing a PC or DF image in the image domain; and the data processing operations include at least a weighting based on the average sensitivity of the measurements along each ray.

[0013] In an embodiment, the per-ray weights above depend on the fan angle of the respective ray.

[0014] In an embodiment, the average sensitivity is further related to measurements along a beam complementary to the beam.

[0015] In an embodiment, the weights represent the average of measurements along a ray halfway between two positions of the imaging device's radiation source on the ray assumed in the scanning operation.

[0016] In an embodiment, the imaging device includes at least one imaging facilitator component operable to facilitate converting radiation intensity detectable at a detector of the X-ray imaging device into a DF or PC signal as DF or PC projection data, and the image domain is located between the detector and the at least one imaging facilitator component.

[0017] In an embodiment, the imaging facilitator component is i) an interference grating, or ii) a coded aperture structure.

[0018] In an embodiment, the system implements a filtered backprojection type tomographic reconstruction algorithm.

[0019] In an embodiment, the system is switchable between two modes, one for reconstructing DF images and another for reconstructing PC images using the same backprojection operation, which allows for a streamlined implementation that is, for example, less memory consuming and more efficient to maintain.

[0020] In an embodiment, the system further implements an iterative further tomographic reconstruction algorithm to process the reconstructed image as initial data to reconstruct a second DF or PC image, thereby improving the performance of the iterative reconstruction algorithm.

[0021] In an embodiment, the system is configured such that in the backprojection operation, the contributions of (some or all) rays from radiation source positions that are 360° apart from each other are normalized to 1 / 2.

[0022] In another aspect, there is provided an imaging apparatus comprising a system according to any one of the above claims and an imaging device.

[0023] In another aspect, there is provided an apparatus in which the imaging device is of the cone beam type. Other imaging geometries of the parallel or fan beam type are not excluded herein, and helical geometries are also envisaged.

[0024] In another aspect, a method for tomographic reconstruction of a phase contrast or dark field image is provided, the method comprising: receiving phase contrast (PC) or dark field (DF) projection data based on measurements along rays through an image domain, the projection data being acquired in a scanning operation by an imaging device; performing data processing operations based on the projection data, the data processing operations including a reconstruction operation for reconstructing a PC or DF image in the image domain; and the data processing operations include at least a weighting based on the average sensitivity of the measurements along each ray.

[0025] In another aspect, a method is provided for calculating weights for use in the weighting. The weights are calculated for each ray angle. For a given ray angle, the weight is calculated as the average of measurements along the ray halfway between two positions of the imaging device's radiation source on the ray assumed in the scanning operation.

[0026] In another aspect, there is provided a computer program element which, when executed by at least one processing unit, causes the processing unit to perform the method.

[0027] In another aspect, at least one computer readable medium having stored thereon a program element is provided, the weights being dependent on the ray angle.

[0028] The projection data are preferably combined in a combination operation to make the projection data more isotropic for improved processing by a filtered backprojection type reconstruction algorithm. Weighting is then applied to the projection data thus combined. In the combination operation, the measurements of each set of complementary rays are combined, such as by calculating an average (such as an arithmetic mean) of the measurements of the set of complementary rays. The weighting corrects for a sensitivity gradient along the rays of the measurements. This sensitivity is caused by the presence of the imaging facilitator structure in the X-ray beam. The sensitivity measures the different signal responses along a given ray.

[0029] The system allows for excellent reconstruction results with few or no artifacts, since it takes into account a certain kind of anisotropy in the sensitivity gradient that was not addressed before. Combining the projection data for the rays by calculating the average or arithmetic mean removes the sensitivity gradient along each ray, but the proposed system goes further than this. By only combining the projection data, the sensitivity still varies with the ray angle. The proposed system corrects this ray angle-dependent angular variation by using weighting, so that the modeling assumptions on which filtered backprojection type reconstruction algorithms are based are better met, which reduces artifacts. The system does not rely on iterative reconstruction algorithms, and faster reconstruction algorithms such as filtered backprojection can be used here.

[0030] The proposed weights can be used for axial and helical reconstructions.

[0031] The same proposed weights can be used for phase contrast and dark field imaging. Thus, the same weighting can be applied before a suitable filtered back projection (FBP) algorithm, the weighting preferably being applied over, for example, a 360° scan.

[0032] In dark-field imaging, the newly proposed weights allow the use of standard FBP as used in attenuated contrast imaging. For differential phase-contrast imaging, some basic methods are also available that differ only slightly from the corresponding attenuated FBP versions.

[0033] Unlike previous iterative approaches, the proposed system allows the reconstruction stage to reconstruct only phase contrast or only dark field images, as required, without the need for one or the other. This is because an improved weighting model allows the use of filtered backprojection type algorithms for phase contrast or dark field. In contrast, iterative reconstruction approaches (e.g., so-called intensity-based iterative reconstruction (IBSIR)) rely on at least two-channel reconstruction, where both phase contrast and dark field images must be reconstructed in parallel with each other, which increases the computational load. However, in the system proposed herein, depending on the desired configuration or output, the system can be adjusted by a user interface to reconstruct either phase contrast or dark field images. However, this does not exclude embodiments in which both dark field and phase contrast images are still reconstructed, where the converter provides two projection data sets, one sinogram for the phase contrast and one sinogram for the dark field signal. Even if only one reconstruction is desired, the conversion stage can still provide both sets of sinograms. This is because some phase contrast / dark field recovery algorithms still essentially operate on two channels each. However, since the transformation step operates in the projection domain, no backprojection operation is required, which is computationally saving compared to the two-channel iterative reconstruction algorithm.

[0034] The proposed weighting can be used in particular for interferometers with inverse geometry, which are less expensive to manufacture than interferometers with direct geometry, although the use of direct geometry interferometers is not excluded herein.

[0035] The system is operable to apply the weighting in the projection domain or during reconstruction. Specifically, in some embodiments, the weighting is applied to the projection data in a data processing operation prior to a reconstruction operation to obtain weighted projection data, which is then reconstructed in a reconstruction operation. Alternatively, the weighting is applied during and / or in the data reconstruction operation. [Brief description of the drawings]

[0036] Exemplary embodiments of the invention are described with reference to the following drawings, which are not to scale unless specifically indicated, in which:

[0037] [Figure 1] FIG. 1 shows a schematic block diagram of a tomographic imaging instrument for dark field and / or phase contrast imaging. [Diagram 2] FIG. 2 shows two different interferometer geometries. [Diagram 3] FIG. 3 shows a schematic diagram of data acquisition based on the line integral model. [Figure 4] FIG. 4 shows a schematic diagram of a scan trajectory with a divergent imaging geometry. [Diagram 5] FIG. 5 illustrates the calculation of sensitivity correction weights as contemplated in embodiments herein. [Figure 6] Figure 6 shows a schematic block diagram of a signal processing system for tomographic phase contrast and / or dark field imaging. [Figure 7] FIG. 7 shows details of the axial cone beam imaging geometry. [Figure 8] FIG. 8 shows details of the helical tomographic imaging geometry. [Figure 9] FIG. 9 shows a flow chart of a signal processing method for tomographic phase contrast and / or dark field imaging. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0038] Reference is first made to the block diagram of Figure 1, which shows a tomographic imaging apparatus AR. The apparatus AR includes an X-ray based tomographic imaging device IA (hereafter simply referred to as "imager") and a signal processing system SPS for processing data generated by the imager IA.

[0039] Broadly speaking, a tomographic imager IA generates projection data that is processed by a signal processing system into a cross-sectional dark-field or phase-contrast image of the object OB.

[0040] The cross-sectional images are then stored in a database DB or other memory and / or processed by a visualizer VS to visualize the images on a display device DD, which will be examined in more detail below.

[0041] The signal processing system SPS implements improved reconstruction techniques that reduce artifacts in tomographic dark field and / or phase contrast images. The objects imaged herein include in particular human or animal patients or parts thereof to obtain cross-sectional images of regions of interest such as the patient's anatomy. However, the following principles are not necessarily limited to medical applications and the imaging of non-living objects such as non-destructive testing is also envisaged herein. However, reference will be made herein primarily to the medical field and where appropriate the "object OB" will be referred to as a patient.

[0042] The imaging device includes an X-ray source XS capable of generating an X-ray beam XB. The beam includes photons propagating along rays that constitute the beam. In operation, the X-ray beam traverses an examination region ER in which a patient OB is present. The X-ray beam interacts with tissue of the patient OB and is modified. The modified radiation is then detected as projection data λ at a detector D. The detector D is positioned opposite the source XS, and between the two an examination region (within which the object OB is located) is defined.

[0043] The detector D comprises a radiation sensitive surface made up of detector sensitive pixel elements which convert radiation intensities into measured values. The layout of the detector pixels which make up the radiation sensitive surface can be either a 1D layout or, preferably, a 2D (matrix) layout. A Data Acquisition Unit (DAQ, not shown) comprises an A / D conversion stage which converts the readings into numerical values ​​to obtain a projection image in digital form.

[0044] The projection images are then transferred to a signal processing system SPS for processing as described in more detail below. The tomographic imager IA, also referred to herein simply as CT, performs a scanning operation during which projection data Λ is collected. In a scanning operation, a focal point FS of an x-ray source XS is moved through the object to obtain a series of projection images along different directions.

[0045] All types of CT imagers with various imaging geometries are contemplated herein. Figure 1 shows, as an example, a third generation CT imager in which an X-ray source including a focal spot traces an imaging trajectory σ. The trajectory σ is either an arch or a closed line (such as a complete circle) as required. The angular coverage is preferably at least 180°, but also includes a full scan in 360°. The rotation axis Z of the scan trajectory σ passes through the region of interest to be imaged, preferably at the isocenter IS. The imaging / rotation axis extends perpendicular to the plane of the drawing in Figure 1.

[0046] The exemplary imaging geometry shown in FIG. 1 is of the diverging beam type, such as a fan or cone beam. However, other imaging geometries, such as a parallel beam, are also contemplated herein, as are CT imagers of first, second or fourth generation. In the third generation type, the X-ray source and a detector located on the opposite side of the X-ray source across the examination region in which the object is present are mounted on a rotating gantry RG. The rotating gantry is located on a fixed (i.e. non-rotatable) gantry SG. The rotating gantry RG rotates about an axis of rotation, and with it the X-ray source and some types of scanners, such as third generation. As mentioned above, during the rotation, projection images are acquired from various angular orientations. However, since in fourth generation scanners a fixed detector ring is located around the examination region, such mechanical rotation of the detector is not necessarily required herein. Furthermore, since in fifth generation scanners the scan trajectory σ is achieved by magnetic deflection of the electron beam to trace the trajectory σ, the mechanical movement of the X-ray source shown in FIG. 1 and all known generations is also not required herein. The electron beam is generated by an electron gun in the X-ray source XS.

[0047] The examination region ER is formed as a bore of the gantry SG, RG. The imager therefore generally has an annular ("donut" shaped) appearance, as shown in the front view of FIG. 1. The examination region ER is the space between the X-ray source and the detector in which the object OB to be imaged resides. The object OB is supported in the examination region during imaging by a support surface SS. The scan trajectory σ is not necessarily limited to one plane per revolution, since helical trajectories are also envisaged herein. In a helical imaging geometry, there is a relative lateral movement between the source and the object induced along the imaging axis Z. This is done, for example, by advancing the patient support SS carrying the patient OB through the bore of the imager while the focal spot traces its (for example circular) trajectory. Alternatively, the stationary gantry SG may be moved laterally. The movement of the rotating gantry (if present) and / or the support surface SS is achieved by an actuator (for example an electric motor, etc.). Control of imaging operations, such as setting imaging parameters (tube current / voltage, etc.) and controlling the actuators, is performed from a user-operable console CS.

[0048] The tomographic imager IA is configured for dark-field and / or phase-contrast imaging. This configuration can be realized by arranging an imaging facilitator component IFC in the examination region as well as an object OB to be imaged, as shown diagrammatically in FIG. 1. In some embodiments (though not necessarily in all embodiments), the imaging facilitator component IFC (referred to herein for brevity as "imaging facilitator IFC") comprises an interferometer. The interferometer comprises at least one, two, or sometimes three gratings G0-G2. The principles of grating-based phase-contrast or dark-field imaging have been described above and will be elaborated below.

[0049] The detector system D used in most imaging systems is only responsive to the intensity of the impinging X-rays. In some detector setups, the phase contrast and the dark field contrast are not originally recorded per se but rely on a computational processing of the data that can be obtained by the operation of the imaging facilitator IFC. In particular, the function of the imaging facilitator IFC is to convert the detected intensity into the dark field projection data λ ε and / or phase contrast projection data λ δ Therefore, the dark field projection data λ can be easily converted into data that can be converted into the data. ε and / or phase contrast projection data λ δ can capture the small angle scattering and refraction contributions of the radiation wavefront, respectively.

number

[0050] Broadly speaking, an imaging facilitator IFC, such as one or more gratings G0-G2, is used as follows: incoming X-rays received from an X-ray source XS interact with the patient's tissue and the imaging facilitator IFC, resulting in a fringe pattern being detected at the detector. The detector may be of a 1D type, such as a column of detector pixels, but is preferably of a 2D type and thus comprises a number of detector pixels arranged in a matrix layout of rows and columns.

[0051] The fringe pattern is known to encode the dark-field and / or phase contrast contribution. For example, the amplitude or loss of visibility of the fringe pattern is known to encode the dark-field signal, and the phase shift of the fringe pattern is known to be related to the phase contrast contribution. Reference data is used, which includes a reference visibility and a reference phase. The recorded fringe patterns are computationally processed by a phase contrast or dark-field recovery algorithm to generate dark-field or phase contrast projection data. That is, for any given angular orientation α0 of the acquisition beam, a respective dark-field and / or phase contrast projection data corresponds. The respective projection data of all angular orientations are calculated by a set of all dark-field projection data λ ε and / or phase contrast projection data λ σ , which is also referred to herein as a dark field or phase contrast sinogram.

[0052] Apply phase or dark-field retrieval to each projection direction to obtain a complete set of dark-field or phase-contrast projection data, λ d , λ pSome such recovery algorithms require a dedicated phase stepping operation at the beginning, in which the grating (at least part of the general imaging facilitator) is subsequently displaced (with respect to the beam XB). In CT applications, such explicit phase stepping is not excluded herein, but is less preferred, since it slows down image acquisition and is cumbersome to implement mechanically. Instead of performing explicit phase stepping, it is preferable to perform a phase step of adjacent angular directions α centered on the current direction α0. j,j=-N,…-1 , α k,k=1…N The signals detected at adjacent pixels along are processed and used in place of explicit phase stepping to capture a good approximation of the respective fringe pattern for each projection direction α0. The phase and visibility losses of the respective fringe patterns are then extracted as described above by any of the appropriate recovery algorithms. Such algorithms are sometimes called phase retrieval algorithms, but this is a misnomer since the dark-field signal is also recovered in such operations as well. Fourier-based techniques can be used to process the fringe patterns and to separate and extract therefrom the dark-field and / or phase contrast contributions. Some such retrieval algorithms recover both the dark-field and phase contrast contributions, while others recover only one contribution (dark-field or phase contrast) as required. Suitable phase / dark field retrieval algorithms are described, for example, by F Pfeiffer et al. in "X-ray dark-field and phase-contrast imaging using a grating interferometer," JOURNAL OF APPLIED PHYSICS, Vol. 105, 102006 (2009) or by I Zanette et al. in "Tri-modal low-dose X-ray tomography," PNAS, PNAS, Vol. 109(26), pp. 10199-10204 (2012).

[0053] Reference is now made to the schematic diagrams of Figures 2A and 2B, which show in more detail an imaging facilitator structure IFC of the grating-interferometer type as envisaged in the embodiments herein.

[0054] In particular, a Talbot-Lau interferometer can be used in different geometries, either in a reciprocal geometry as shown in FIG. 2A or in a direct geometry as shown in FIG. 2B. The interferometer comprises two gratings G1 (also called phase grating) and G2 (also called analyzer grating). Preferably, G2 is an absorption grating and G1 is a phase grating or an absorption grating. G2 may be integrated into the detector D.

[0055] An interferometric setup for tomographic phase-contrast and dark-field imaging is described by M von Teuffenbach et al. in “Grating-based phase-contrast and dark-field computed tomography: a single-shot method,” Sci Rep7, 7476, (2017).

[0056] In the direct geometry (FIG. 2B), the two gratings G1, G2 are arranged downstream of an examination area (in particular of the object OB), so that during acquisition of the projection data the two gratings are located between the object and the detector D. The examination area in this direct geometry is between the X-ray source and the grating pack formed by the two gratings G1 and G2.

[0057] In the inverse geometry (FIG. 2A), the examination region / imaged object OB is between the phase grating G1 and the analyzer grating G2.

[0058] As shown in Figures 2A and 2B, the pitches of at least the G0 and G2 gratings are different for the two different geometries.

[0059] If the X-rays provided by the source XS are incoherent, there is an additional absorption grating, called the source grating G0, placed between the focal point FS of the XR source and the object OB to increase the coherence of the emitted radiation. If the X-rays generated by the source G0 are inherently coherent, the source grating G0 is not required.

[0060] The distances between G0 and G1 and between G1 and G2 are specifically adjusted according to the Talbot-Lau setup described elsewhere. The distances between G0 and G1 and between G1 and G2 need to be fine-tuned to fit the Talbot distance requirement of the appropriate order, which is a function of the "pitch" of the respective gratings (i.e., the spatial period in the grating scale) and the average energy of the X-rays generated by the X-ray source XS.

[0061] Gratings can be fabricated by photolithography, or by cutting a silicon wafer to define a periodic pattern of trenches. In the case of absorption gratings, the gaps between the trenches are filled with lead, gold, or other high-Z material.

[0062] The distances L0, L1, and L2 shown in Figures 2A and 2B represent the distances of the source grating, phase grating, and analyzer grating, respectively. Preferably, the distances L0, L1, and L2 are measured from the focal point, as defined herein. Other interferometers, not necessarily Talbot-Lau type, are contemplated herein, and the present disclosure is not limited to Talbot-Lau type interferometers or interferometers in general, as previously discussed.

[0063] Phase contrast and / or dark field tomography imaging has traditionally been found to suffer from various imaging artifacts, such as cupping type artifacts. In Applicant's WO'657, it has been found that some of these artifacts can be addressed by a magnification effect. It has been found that the imaging signal is magnified depending on the position that each image point in the image domain has relative to a grid, or more generally, an imaging facilitator structure IFC.

[0064] Another effect distinct from magnification is observed: a certain non-uniform measurement sensitivity s. While magnification is a purely geometric effect that spreads the image signal across adjacent pixels, sensitivity s is a physical effect, not a geometric one. The (measurement) sensitivity s describes the different signal strengths caused by a virtual test object moving in the image domain along an imaging ray L. Some image points have higher signal strengths than others along the same ray. Sensitivity s is a function of the distance along the ray L to the imaging facilitator IFC of the image point. The mere presence of the imaging facilitator IFC in the ray appears to induce a sensitivity gradient along the ray.

[0065] In more detail, a tomographic measurement collection operation is shown diagrammatically in Figure 3. The detector D includes detector pixels that detect respective measurements m. A measurement m can be modeled as a line integral of a quantity of interest q along a ray L that intersects the detector at a given pixel where the measurement is made:

number

number

[0066] Equation (1b) is the forward model of the refractive index depletion σ, and equation (1a) is the forward model of the linear diffusion coefficient ε. The operator ∂ x refers to the differential coefficient taken perpendicular to the grating trench. L is the direction

number

number

[0067] A particular pixel of the detector D records a line integral of measurements collected along a light ray L. The light ray extends from a focal point FS at a given position on the trajectory σ, through an image object OB in the image domain, to the respective detector pixel. Conceptually, the image domain ER in which the object OB resides (and thus the object OB itself) can be thought of as being made up of image points (also called voxels).

[0068] It is known that measurements collected at various image points in the image domain along the ray L will yield different contributions with different intensities for the measurement m detected at the detector D. More specifically, for two different image points on a given ray L, an image point l proximal to the source XS, p and the distal image point l d Consider the following: d , l p A test object, such as a diffuser or refractor, placed at each of the 3D optics generates a signal with a different respective intensity. Thus, the measurement sensitivity s, i.e., the signal intensity, varies along the ray with the distance to the imaging facilitator IFC.

[0069] The spatial dependence along the ray can be modeled as a ramp function of the distance between G1 and G2 or between G0 and G1, depending on whether the inverse geometry (FIG. 2A) or the direct geometry (FIG. 2B) is used. The dependence of the sensitivity along the ray is different for the direct and the inverse geometry. The ramp function of the signal intensity or sensitivity s and its dependence on the distance from the imaging facilitator structure is described by Donath et al. in "Inverse geometry for grating-based x-ray phase-contrast imaging", JOURNAL OF APPLIED PHYSICS 106, 054703, (2009)

[0070] The sensitivity s is given by, for the direct geometry,

number

number

[0071] Herein, we propose to use an improved weighting in the reconstruction of cross-sectional dark-field or phase-contrast images that takes into account better accuracy the sensitivity along the ray experienced by the virtual sample test object as it moves along the ray L. In addition, the sensitivity varies with the ray angle, and this variation is also taken into account in the improved weighting. For this reason, we propose a high-fidelity weighting model that allows to reduce artifacts in the reconstructed phase-contrast and / or dark-field images caused at least in part by the sensitivity gradient. As mentioned above, this measured sensitivity gradient along the ray is different from and in addition to any magnification effect that may occur due to, for example, a diverging beam geometry. In particular, the sensitivity effect does not usually mean that the signal is spread over several pixels, as is the case for example with magnification. Since this measured sensitivity is a per detector pixel phenomenon, there is a corresponding sensitivity effect in the fringe pattern detected during acquisition at the detector D.

[0072] Before detailing the improved weighting model, reference is first made to FIG. 4, which illustrates the concept of complementary measurements and associated anisotropy in tomographic sinograms.

[0073] Specifically, FIG. 4 provides a view along the imaging axis Z of a CT imaging geometry with a diverging beam shape, such as a fan angle or cone beam type. The system fan angle represents the spread of the beam used for each exposure. Two complementary positions of the focal point on the scan trajectory σ are shown at angular positions β and β'. A scan trajectory σ that covers a sufficient angular range, such as 180°, gives rise to the concept of complementary measurements m, m' along different rays. Measurements (and their rays) are collected along two different rays emanating from different positions of the focal point on the trajectory σ, and are complementary if the two positions are defined by two intersections of a line intersecting the trajectory. For the two exemplary positions β, β' in FIG. 4, the two complementary rays

number

number

[0074] Conceptually, the complementary rays are divided into equivalence classes, and each equivalence class is assigned a respective combined measurement from the original sinogram (phase contrast or dark field sinogram). Thus, a new sinogram (typically λ ) is generated, which we call the combined sinogram. ε ', λ σ For example, as shown in Figure 4, we can obtain the complementary ray

number

number

[0075] Therefore, for divergent beam shapes and similar shapes used herein, the projection sinograms are combined (e.g., averaged). It is known that the combined (e.g., averaged) value depends on the respective ray angle or fan angle (distinct from the system fan angle), which is the direction defined by the position of the focal point and the inclination of the individual rays emanating from the focal point position. In a fan beam setup, the complementary rays of each setup are divided into two rays.

number

number

number

number

[0076] The ray angle α is shown in FIG. 5. The ray angle is measured relative to the central ray R of the fan beam. The central ray R extends from the source position on the trajectory to the rotation axis Z. FIG. 5 also shows a view along the imaging / rotation axis Z mentioned above through the isocenter IS, together with the scan trajectory σ. As mentioned above, a combination operation such as averaging eliminates the sensitivity variation along the ray. However, there still remains a generally undesirable sensitivity variation with the ray angle α. To correct this variation, a weight c(α) of each of the combined complementary measurements for a given angle α is calculated. For example, as shown in FIG. 5, the sensitivity correction weight c(α) is calculated as the average of the complementary and direct measurements. Specifically, as also shown in FIG. 5, the weight c(α) is calculated midway along a given ray L. This is calculated by using the cosine of the ray angle α. The weight c(α) proposed here is therefore a purely geometric quantity based on the imaging geometry at hand.

[0077] In particular, the arithmetic mean / average or other suitable combination of measurements of the direct and complementary rays appears to be insensitive to the sensitivity gradient since the sum of the sensitivities along the direct and complementary rays is a constant. The value of this constant depends on the ray angle α. The value of the constant is calculated, for example, by calculating the value of s(l) halfway between the actual and complementary source positions (see Figure 5): c(α)=s(L(α))=s(R cos α) (2)

[0078] Thus, for the inverse geometry, the weight c(α) is calculated as follows:

number

[0079] For direct geometry, the weight c(α) is calculated as follows:

number

[0080] The calculated weights are different and depend on the ray angle α. Thus, preferably, various such weights can be calculated, one weight for each α. The weights are applied to the combined (e.g., averaged) sinogram. For example, in an embodiment, the measurements calculated for the combined dark field or phase contrast sinogram are divided by a weight c(α) and c(α) is applied pointwise to the combined measurements for each corresponding angle α.

[0081] That is, the combined phase contrast and dark field data is ε ', λ σ ', the improved weighting model can be constructed as follows:

number

number

[0082] We now refer to Fig. 6 which shows a block diagram of the new proposed signal processing system SPS. The signal processing section SPS acts as a transformation stage which converts the projection images recorded in the projection domain at the detector into image values ​​in the image domain in the examination region to calculate the cross-sectional images. The imaging domain, i.e. the examination region, can be considered as consisting of a virtual grid of voxels as mentioned above. In the reconstruction operation, an image value is calculated for each of these voxels / grid locations.

[0083] The signal processing system includes a pre-processing stage PS operating in the projection domain and a data processing stage DP performing a transformation from the projection domain to the image domain. In this stage, an analytical direct reconstruction algorithm such as filtered back projection FBP is preferably used. FBP type reconstruction algorithms and the like can be implemented directly by using Fourier-based methods, as opposed to algebraic or iterative reconstruction methods that use a discretized version of the system matrix. The system matrix is ​​an explicit description of the contribution of the image domain voxel I to the ray measurement m. Also, iterative reconstruction algorithms, as the name suggests, perform iterative processing, whereas analytical direct reconstruction algorithms such as filtered back projection can directly calculate cross-sectional images without iterations. Iterative reconstruction algorithms are not excluded here and can be used in combination with analytical direct reconstruction such as FBP, as will be described in more detail below.

[0084] First, reference is made to the pre-processing section PS. An imaging facilitator IFC operates to receive intensity-based projection data from original sinograms including fringe patterns and pass them to a conversion stage CV. The conversion stage CV processes the fringe patterns and converts the intensity projection images into dark-field and / or dark phase projection data (phase-contrast or dark-field sinograms). As explained above, a phase / dark-field recovery algorithm is used at this stage. The conversion stage is optional if the data acquisition setup allows for the acquisition of phase-contrast or dark-field projection data originally.

[0085] The dark field or phase contrast projection data is then processed by a combination component CB to establish at least partial data isotropy for the complementary rays, such that there is constant sensitivity along each ray. The combination component CB can combine the projection data for each fan angle, for example, by calculating an average or other arithmetic mean for the complementary measurements. For example, if m represents a direct measurement of the angle α and m' represents a complementary measurement, the combined projection data for this angle is calculated as (m+m') / 2. Other combination aspects may alternatively be contemplated. The combined (e.g., averaged) projection data / sinogram of the dark field or phase contrast signal is then passed to a corrector DP.

[0086] The sensitivity gradient corrector CP processes the combined projection data by applying correction weights c(α) depending on the ray angle α. For example, the combined projection data is divided point-wise by respective weights c(α) depending on the ray angle α. The weights are calculated by the weight calculation unit WCU based on the imaging geometry according to any one of the equations (2)-(4). The corrector CP corrects the variation of the respective constants with the ray angle α, so that perfect isotropy is achieved.

[0087] The combined projection data thus weighted are then passed to a data processing stage DP via an input IN. A reconstructor RECON of the data processing stage DP preferably implements an analytical tomographic reconstruction algorithm, such as of the FBP type. FBP type algorithms use filters to combine the weighted projection data.

number

[0088] Some application scenarios are as follows: For example, in the type of axial reconstruction, preferably full (360°) scan data λ is acquired, which is combined by a combiner CB to obtain λ′. Then, a corrector CP divides the combined data by a weight c(α) according to equations (5a, b).

[0089] It should be noted that the order of operation of the corrector CP and the combiner CB is interchangeable. In particular, it is envisaged to apply the corrector CP first and then the combiner CB, or vice versa. While the order of operations is not actually important for 2D reconstruction, in cone beam geometry it is preferred to apply the corrector CP first and the combiner CB afterwards, since the selection of the appropriate complementary rays to combine depends on the particular given voxel to be reconstructed. Furthermore, it is preferred that the combination operation is applied during reconstruction, rather than as a pre-processing before the reconstruction operation.

[0090] In cone-beam geometry, FDK reconstruction can be applied to dark-field CT. The FDK algorithm is an FBP-type reconstruction algorithm and is described by LA Feldkamp et al. in "Practical cone-beam algorithm," Journal of the Optical Society of America A, Vol. 1, No. 6, pp. 612-619 (1984)). Modified FDK-type reconstruction algorithms can also be used for phase-contrast CT, as described, for example, by Z Qi and G Chen in "Direct Fan-Beam Reconstruction Algorithm via Filtered Backprojection for Differential Phase-Contrast Computed Tomography," Hindawi Publishing Corporation, X-Ray Optics and Instrumentation, Vol. 2008, Paper ID 835172, p. 8. The FDK algorithm can be used for planar or other detector geometries by using appropriate filtering.

[0091] In Fig. 7, the operation of the combiner CB for the cone beam axial trajectory and the corresponding reconstruction is shown. In particular, Fig. 7 shows a ray R originating from a source position S to be backprojected. Two exemplary voxels v1 and v2 are shown. There are multiple rays originating from complementary source positions C, all in the same plane containing the ray R. The plane is parallel to the rotation axis Z. For the reconstruction operation, for each voxel, from the multiple rays, a ray is selected that passes through the voxel of interest. For example, the complementary ray of ray R to voxel v1 is ray CR1, and for voxel v2, ray CR2 is complementary to R. That is, in a cone beam, the complementarity depends on the voxel position on a given ray. Different voxels on the same ray define different complementary rays. Therefore, in a cone beam or other similar geometry, it is preferable to combine the data during reconstruction, rather than as a pre-processing step before reconstruction.

[0092] The operation of the combiner CB for the helical trajectory and the corresponding reconstruction is shown in Fig. 8. As in Fig. 7, two source positions D1, S are shown along with complementary positions C1 and C2 of complementary rays CR1, CR2 relative to ray R.

[0093] A common aspect of the preferred operation of the combiner in both axial and helical cone beam CT is that a so-called 360° normalization is applied. This concept was introduced in another context for high-resolution helical CT in "High-Resolution Images of Cone Beam Collimated CT Scans" by Shechter et al., IEEE Transactions on Nuclear Science 52(1), vol. 247 (2005). The basic concept is to apply a special normalization during backprojection. For dose efficiency reasons, one principle of cone beam FBP is that the reconstruction of a particular voxel V needs to consider all projections representing rays passing through that voxel V, and the manner in which rays contribute to a given voxel V is preferably normalized to 1. This normalization takes into account imaging events where a given voxel V may have been "seen" more frequently from some projection directions than from other directions.

[0094] In the general case, when performing a backprojection of a ray R from a source position S, there is a complementary ray CR to be considered. i (i=1, …, N c ) to consider. i (i=1, ..., ND). In such a case, the combiner unit CB calculates the normalized weights w0 of the ray R and the normalized weights w i (i=1, . . . , N D ), as well as the complementary ray normalization weights v i (i=1, …, N C ) satisfies the following condition:

number

[0095] In the 2D case mentioned above, there is only one complementary ray (N C =1) and there are no further direct redundant rays. That is, N D = 0. The normalization condition (6) above therefore reduces to the operation of the combiner CB with w0 = v1 = 0.5. Thus, in this embodiment, there is only an averaging of two rays, as described in the previous embodiment.

[0096] The complete cone-beam FBP reconstruction of a voxel V is formulated as follows: V=BP CB (F(CP(c(α),λ))) (7) where (CP(c(α),λ) describes the operation of the corrector CP on the data with the value c(α), F(·) is the usual filter operation of the FBP algorithm, and BP CB (·) is the backprojection operation that includes the combiner CB by ensuring that the normalized weights of the rays satisfy the normalization condition (6).

[0097] In both helical and axial CT, to facilitate normalization of the complementary rays, it is desirable to acquire data such that each voxel to be reconstructed is within a cone of at least 360°. Thus, the projection data is preferably low pitch (<1) in helical CT, and acquired using at least a full scan to ensure that all voxels are within a cone of at least 360°.

[0098] As mentioned above, the above reconstruction operations by the processor DP are applied to either the phase contrast projection data or the dark field projection data as required, whichever modality is more useful to the clinical user at any given time. In the proposed method, there is no need to reconstruct both the phase contrast and dark field image projection data. Therefore, a user interface UI is provided that allows the user to select whether to reconstruct the phase contrast image or the dark field image, or both, as required and where latency is not a major concern.

[0099] Although the forward models may differ slightly for the phase-contrast and dark-field image channels, the exact same backprojection operation can be used for either the dark-field or phase-contrast reconstruction algorithms, as a testament to the proposed computational efficiency of the proposed system.

[0100] The above principles can be applied not only to parallel beam geometries but also to diverging beam geometries such as fan beams and cone beams. The fan beam geometry is transformed to achieve a parallel beam geometry by relying on the projection data (also known as a rebinning operation). In this case, the above-mentioned pre-processing chain PS in the projection domain can include a rebinner RB.

[0101] Unlike previous approaches that focus only on magnification correction, such as WO'657, the proposed system can be used in both inverse and direct geometry interferometric setups. Advantageously, in the inverse geometry, the largest grating G2 at detector D can be manufactured with a larger pitch, as shown by the larger gap in Figure 2A, which results in lower manufacturing costs.

[0102] The proposed system recognizes the distinct nature of the measurement sensitivity as opposed to the magnification effect: while the measurement sensitivity induced by the presence of the imaging facilitator IFC is similar to the magnification effect in a direct geometry far from the source grid G0, the two effects are significantly different near the source grid and even opposite in a reverse geometry. The proposed system takes both effects into account.

[0103] As mentioned above, the proposed reconstruction using the improved weighted forward model can be used in combination with an iterative reconstruction algorithm. In this embodiment, the data processor stage includes two reconstruction modules in series, one module RECON for FBP or other analysis algorithms, and a second reconstruction module RECON' downstream from it, which implements a suitable iterative reconstruction scheme. In this combination embodiment, cross-sectional images reconstructed by the first reconstructor, e.g. via FBP, are used as initial values ​​to the second iterative reconstructor RECON'. This is because it has been observed that certain iterative reconstruction algorithms get trapped in local minima or otherwise take a long time to converge to a final result. The selection of initial images calculated with FBP can improve performance. These images are then refined during the iterations, allowing faster convergence. The noise behavior can be improved by using an iterative reconstruction algorithm.

[0104] Reference is now made to the flow chart of Fig. 9. The above system can be realised using the method steps described below, although it will be understood that the following method steps can be understood on their own and are not necessarily associated with the above system. The proposed system implements an improved signal processing method, in particular for use in the reconstruction of tomographic phase contrast and / or dark field images.

[0105] In step S910, projection data acquired in a scanning operation by a tomographic imaging device for dark field and / or phase contrast imaging is acquired. Preferably, the device comprises one or more image facilitator structures, such as gratings, arranged in the examination region / imaging domain.

[0106] In step S920, the acquired intensity-based projection data is converted into phase contrast or dark field signal projection data, although step S920 is optional since some imaging devices can inherently provide phase contrast or dark field contrast projection images.

[0107] In step S930, the dark field or phase contrast data is combined as described above for measurements along the complementary beams to achieve at least partial isotropy. In particular, the variation in sensitivity along the beams is eliminated or at least reduced. For example, the measurements for each complementary beam are averaged or otherwise combined. In this manner, combined projection data / sinograms are obtained.

[0108] In step S940, measurement sensitivity weights c(α) are calculated depending on the ray angle α of the combined projection data. Preferably, a weight is calculated for each α within the angle range covered by the scanning operation, but this may not be required in all embodiments. Therefore, to increase efficiency, in some application scenarios it is sufficient to calculate the weights only for a subset of the angle range. The weights may be pre-calculated before the scanning operation. The weights are calculated according to any of the above equations (2), (3), (4), etc.

[0109] The method may further include some or all of the following data processing steps: In particular, in step S950, weights are applied to the combined projection data obtained in step S930. The order of steps S940 and S950 may be reversed, especially in a helical scan setup. The weight application step S950 also eliminates the variation of sensitivity with ray angle, thereby achieving full isotropy.

[0110] In step S960, the weighted combined projection data is processed by a reconstruction algorithm, such as filtered backprojection or other analytical direct reconstruction algorithm, resulting in cross-sectional images for each cross-sectional plane or in full 3D.

[0111] In step S970, the reconstructed cross-sectional images are output.

[0112] In optional step S980, the output cross-sectional phase contrast or dark-field tomographic image is provided as initial data to an iterative reconstruction algorithm to calculate a new improved version of the tomographic cross-sectional image.

[0113] It will be appreciated that the step of applying weights S950 may occur prior to reconstruction in S960. Alternatively, depending on the linearity of the backprojection operation, the application of weights in S950 may occur together with S960 during the reconstruction operation. The weights may be applied in the filtering step or in the backprojection step.

[0114] Preferably, a combination step is applied during reconstruction, such as in a helical scan setup, as shown in Figures 7-8.

[0115] The components of the signal processing system SPS may be implemented as one or more software modules executed on one or more general purpose processing units PU, such as a workstation associated to an imager IA, or on a server computer associated to a group of imagers. The system SPS may also be integrated into the console CS of the imager.

[0116] Alternatively, some or all of the components of the SPS may be located in hardware, such as a suitably programmed microcontroller or microprocessor, such as an FPGA (Field Programmable Gate Array), or as a hardwired IC chip that is an application specific integrated circuit (ASIC) integrated into the imaging system IA. In further embodiments, the system SPS may be realized both partly in software and partly in hardware.

[0117] The various components of the SPS may be implemented on a single data processing unit PU, or alternatively, several or several components are implemented on different processing units PU, possibly located remotely in a distributed architecture and connectable in a suitable communication network, such as in a cloud configuration or a client-server setup.

[0118] One or more features described herein may be configured or implemented as or with circuitry encoded in a computer-readable medium, and / or combinations thereof, including discrete and / or integrated circuits, systems-on-a-chip (SOC), machines, computer systems, processors and memories, computer programs.

[0119] In another exemplary embodiment of the invention, a computer program or a computer program element is provided, characterized in that it is adapted to carry out, on a suitable system, the method steps of the method according to one of the above embodiments.

[0120] Thus, the computer program element may be stored in a computing unit which may be part of an embodiment of the present invention. This computing unit may execute or induce the execution of the steps of the above-mentioned method. Furthermore, it may operate the above-mentioned device components. The computing unit may operate automatically and / or execute user instructions. The computer program may be loaded into the working memory of a data processor. The data processor is thus ready to execute the method of the present invention.

[0121] This exemplary embodiment of the present invention covers both computer programs that use the present invention from the beginning and computer programs that convert existing programs, via updates, into programs that use the present invention.

[0122] Moreover, the computer program element may provide all the steps required to carry out the procedures of the exemplary embodiments of the methods described above.

[0123] According to a further exemplary embodiment of the present invention, a computer readable medium, such as a CD-ROM, is presented having stored thereon computer program elements, as described in the previous section.

[0124] The computer program may be stored / distributed on a suitable medium (in particular, but not necessarily a non-transitory medium), such as an optical storage medium or a solid-state medium, supplied together with or as part of other hardware, but may also be distributed in other forms, such as via the Internet or other wired or wireless communication systems.

[0125] However, the computer program may also be presented via a network such as the World Wide Web and can be downloaded from such a network into the working memory of a data processor. According to a further exemplary embodiment of the invention, a medium making available a computer program element for downloading is provided, the computer program element being configured to perform a method according to one of the aforementioned embodiments of the invention.

[0126] It should be noted that the embodiments of the present invention are described with reference to different subject matters. In particular, some embodiments are described with reference to method type claims, while other embodiments are described with reference to device type claims. However, a person skilled in the art can infer from the above and following description that, unless otherwise specified, any combination of features belonging to one type of subject matter, as well as any combination of features related to different subject matters, are considered to be disclosed in the present application. However, all features can be combined if they provide a synergistic effect that is more than a mere collection of features.

[0127] While the invention has been illustrated and described in detail in the drawings and the foregoing description, such illustration and description are to be considered illustrative or exemplary and not restrictive. The invention is not limited to the disclosed embodiments. Other variations of the disclosed embodiments can be understood and effected by those skilled in the art in practicing the claimed invention, from a study of the drawings, the disclosure, and the dependent claims.

[0128] In the claims, the word "comprising" does not exclude other elements or steps, and the singular elements do not exclude a plurality. A single processor or other unit may fulfill the functions of several items recited in the claims. The mere fact that certain measures are recited in mutually different dependent claims does not indicate that a combination of these measures cannot be used to advantage. Any reference signs in the claims should not be interpreted as limiting the scope. Such reference signs may be composed of numbers, letters or any alphanumeric combination.

Claims

1. A system for tomographic reconstruction of phase contrast or dark field images, the system comprising: Receiving projection data of phase contrast or dark field based on measurements along rays passing through the image domain, the projection data being acquired in a scan operation by an X-ray imaging device, receiving the projection data; Performing a data processing operation including a reconstruction operation for reconstructing a phase contrast image or a dark field image in the image domain based on the projection data; And performing The data processing operation includes at least weighting based on the average sensitivity of the measurements along each ray; The weight for each ray depends on the fan angle of the respective ray; System.

2. The system according to claim 1, wherein the average sensitivity is further related to measurements along rays complementary to the ray.

3. The system according to claim 1, wherein the weight represents the average of the measurements along the ray between two positions of the radiation source of the X-ray imaging device on the ray, and the two positions are assumable by the radiation source in the scan operation.

4. The X-ray imaging device includes at least one imaging facilitator component operable to convert the radiation intensity detectable at the detector of the X-ray imaging device into a dark field signal or a phase contrast signal as the projection data of the dark field or phase contrast, and the image domain is located between the detector and the at least one imaging facilitator component. The system according to claim 1.

5. The system according to claim 4, wherein the imaging facilitator component is i) an interference grating or ii) an encoded aperture structure.

6. The system according to claim 1, wherein a data processor implements a filtered backprojection type tomographic reconstruction algorithm.

7. The system according to claim 6, which is switchable between two modes, one mode for reconstructing a dark field image and another mode for reconstructing a phase contrast image, using the same backprojection operation.

8. The system according to claim 1, further implementing an iterative tomographic reconstruction algorithm to process the reconstructed image as initial data for reconstructing a second dark field image or a second phase contrast image.

9. The system according to claim 6, wherein in the back-projection operation, the contribution of rays from radiation source positions 360° apart from each other is normalized to 1 / 2.

10. The system according to any one of claims 1 to 9, said X-ray imaging device An imaging device comprising: said X-ray imaging device is preferably a cone beam type, imaging device.

11. A method for tomographic reconstruction of a phase contrast or dark field image, comprising: Receiving projection data of phase contrast or dark field based on measurements along rays passing through an image domain, wherein the projection data is acquired in a scan operation by an X-ray imaging device, the step of receiving the projection data; Performing a data processing operation including a reconstruction operation for reconstructing a phase contrast image or a dark field image in the image domain based on the projection data and including the data processing operation includes at least weighting based on the sensitivity of the average of the measurements along each ray, the weight for each ray depends on the fan angle of each ray, method.

12. A method for supporting tomographic reconstruction of a phase contrast or dark field image, comprising: Calculating weights for use in weighting in a reconstruction algorithm, the weights are calculated as the average of the measurements for each ray angle, between two positions of the radiation source of the imaging device on the ray, which are assumable in a scan operation.

13. A computer program that, when executed by at least one processing unit, causes the at least one processing unit to perform the method according to claim 11 or 12.

14. At least one computer-readable medium storing the computer program according to claim 13.

Images