WO2012080125A1

WO2012080125A1 - A method and a system for image integration using constrained optimization for phase contrast imaging with an arrangement of gratings

Info

Publication number: WO2012080125A1
Application number: PCT/EP2011/072332
Authority: WO
Inventors: Marco Stampanoni; Thomas THÜRING
Original assignee: Paul Scherrer Institut
Priority date: 2010-12-13
Filing date: 2011-12-09
Publication date: 2012-06-21
Also published as: EP2652708B1; CA2821145A1; JP2014506352A; AU2011344365A1; CN103460251A; US20130279659A1; EP2652708A1; JP5818909B2

Abstract

The present invention yields high-quality, artifact free phase contrast images from an object using an arrangement of gratings. The new method suppresses the need of direct image integration and significantly improves the quality of phase contrast images. In comparison with existing techniques, our approach does not need additional alignment work nor increased exposure time. On the other hand it delivers excellent, direct interpretable information about the phase projection within a radiographic experiment. Due to its general applicability and its simplicity in usage, the suggested invention is expected to become a standard method for a variety of 2D imaging applications using gratings arrangements in particular on medical scanners (for instance mammography), inspection at industrial production lines, non-destructive testing, and homeland security.

Description

A method and a system for image integration using constrained optimization for phase contrast imaging with an arrangement of gratings

The present invention relates to a method and a system for the recovery of an integrated image from a differential image by using constrained optimization.

It is well known that, differently from conventional visible light optics, the refractive index in X-ray optics is very close to and smaller than unity. In first approximation, for small and negligible anisotropy in the medium, the index of refraction characterizing the optical properties of a tissue can be expressed - including X-ray absorption - with its complex form: η = \—δ—ίβ where δ is the decrement of the real part of the refractive index, characterizing the phase shifting property, while the imaginary part β describes the absorption property of the sample. In conventional absorption-based radiography, the X-ray phase shift information is usually not directly utilized for image reconstruction. However, at photon energies greater than 10 keV and for light materials (made up of low-Z elements), the phase shift term plays a more prominent role than the attenuation term because δ is typically three orders of magnitude larger than β . As a consequence, phase-contrast modalities can generate significantly greater image contrast compared to conventional, absorption-based imaging.

Furthermore, far from absorption edges, δ is inversely proportional to the square of the X-ray energy whilst β decreases as the fourth power of energy. A significant consequence of this mechanism is that phase signals can be obtained with much lower dose deposition than absorption, a very important issue when radiation damage has to be taken into account such as in biological samples or in living systems .

Several approaches have been developed in order to record the phase signal. They can be classified as interferometric methods (with crystals), phase propagation methods, techniques based on an analyzer crystal, or on x-ray gratings. The described invention is in particular in context with the latter technique.

Grating based x-ray imaging setups essentially detect the deflections of x-rays in the object. Such deflections can be either caused by refraction on phase shift gradients in the object resulting in differential phase contrast (DPC) or by scattering on in-homogeneities in the sample resulting in the so-called dark-field image (DFI) contrast. The DPC image signal can be used to obtain phase contrast (PC) images by image processing routines. As shown in Figure 1, set-ups with two gratings (Gl and G2 ) or three gratings (GO, Gl, and G2 ) can be applied to record the deflection of the x-rays. In the case of a two-grating set-up, the source needs to fulfill certain requirements regarding its spatial coherence, while in a three grating setup no spatial coherence is required. A GO grating is required, when the source size is bigger than p2*l/d, where p2 is the period of G2, 1 is the distance between the source and Gl, and d is the distance between Gl and G2. Therefore, the three grating set-up is suited for use with incoherent x-ray sources, in particular with x-ray tubes.

To separate the conventional attenuation contrast (AC) from the DPC and DFI contrast, a phase-stepping approach is applied. One of the gratings is displaced transversely to the incident beam whilst acquiring multiple images. The inten- sity signal at each pixel in the detector plane oscillates as a function of the displacement. The average value of the oscillation represents the attenuation contrast (AC) . The phase of the oscillation can be directly linked to the first derivative of the wave-front phase profile and thus to the DPC signal. The amplitude of the oscillation depends on the scattering of x-rays in the object and thus yields the DFI signal . For the (two or three) gratings, several approaches have been proposed and applied. The grating GO (if required) is the one closest to the source. It usually consists of a transmission grating of absorbing lines with the period pO . It can be replaced by a source that emits radiation only from lines with the same period. The grating Gl is placed further downstream of the source. It consists of lines with a period pi. The grating G2 is the one most downstream of the setup. It usually consists of a transmission grating of absorbing lines with the period p2. It can be alternatively replaced by a detector system that has a grating-like sensitivity with the same period.

Two regimes of setups can be distinguished: in the so called "near field regime" and the "Talbot regime". In the "near field regime", the grating period p, grating distances d and the x-ray wavelength λ are chosen such, that diffraction effects are negligible. In this case, all gratings need to consist of absorbing lines. In the "Talbot regime", diffraction on the grating structures is significant. A sharp dis- tinction between the two regimes is not easily given, as the exact criterion depends on the duty cycle of the grating structure, and whether the gratings are absorbing or phase shifting. E.g., for a grating with absorbing lines and a duty cycle of 0.5, the condition for the "near field regime" is d > ρ²/2λ. Here, Gl should consist of grating lines that are either absorbing or, preferentially, phase shifting. Several amounts of phase shift are possible, preferentially π/2 or multiples thereof. The grating periods must be matched to the relative distances between the gratings. In case of setups in the "Talbot regime" the Talbot effect needs to be taken into account to obtain good contrast. The formulae for the grating periods and distances are known in the prior art. The sample is mostly placed between GO of Gl (or upstream of Gl in case of a two-grating set-up) , however it can be advantageous to place it between Gl and G2.

The presented invention is relevant in all of the abovemen- tioned cases, i.e. in the two- and three-grating case, in the case of the "nearfield regime" and the "Talbot regime", and for the sample placed upstream or downstream of Gl .

In addition, the invention presented here also works in com- bination with scanning-based systems, for parallel and quasi parallel geometries using planar gratings or for compact fan-beam or cone-beam geometries using cylindrically or spherically curved gratings. An imaging experiment with the abovementioned system is sensitive to X-ray refraction. The relation between the refractive angle a (relative to the optical axis) and the differential phase shift measured with the system is given by:

where d is the inter-grating distance and 2 is the period of G2. The x-axis corresponds to the phase stepping direction. Furthermore, the relationship between (x) and the phase profile φ(χ) of the wave after the sample is given by [12] :

Substituting Eq. (2) into Eq. (1) leads

Therefore, Αφ is proportional to the first derivative of φ(χ) . According to this equation, the reconstruction of φ(χ) requires an integration of the DPC measurement in x direction,

In general, the integration of a noisy signal accumulates the noise errors and the noise variance. As a result, horizontal stripe artifacts with increasing amplitudes in the direction of integration appear. Figure 2 shows the integration of a noisy DPC image, generated from a modified Shepp- Logan phantom, where in Figure 2c, the stripe artifact are clearly visible. It is therefore an objective of the present invention to provide a system and the method to improve the quality of the integrated phase contrast image.

This objective is achieved according to the present inven- tion by the features of the independent claims 1 and 19.

Preferred embodiments of the present invention are given in the dependent claims 2 to 18 and 20 to 36 resp. The inventive system and the inventive method for the recovery of an integrated image from a differential image by using constrained optimization are comprising:

a) means for identifying an image that minimizes a cost function \\Tf\\ subject to a constraint

||^(£>_χ/-φ)||^ , which is written as

minimize

subject to : ≤ε

b) wherein f is an image vector, T is a transform operator matrix, D_x is a differentiation operator matrix, W a diagnoal matrix containing the reciprocal noise standard deviations, φ is the measured image vector and ε a boundary for the noise power. ... indicates the p-norm of a vector, where p is a positive number, preferably an integer in the interval [0,2] .

A preferred embodiment of the present invention is achieved when the differential image is obtained from an arrangement for x-rays, in particular hard x-rays, for obtaining quanti- tative x-ray images from a sample comprising:

a. an X-ray source;

b. three or at least two gratings dubbed GO, Gl and G2 or Gl and G2 resp.;

c. a position-sensitive detector with spatially modu- lated detection sensitivity having a number of individual pixels;

d. means for recording the images of the detector;

e. means for evaluating the intensities for each pixel in a series of images in order to identify the char- acteristic of the object for each individual pixel as an absorption dominated pixel and/or a differen- tial phase contrast dominated pixel and/or an x-ray scattering dominated pixel;

f. wherein the series of images is collected by continuously or stepwise rotating from 0 to n or 2n ei- ther the sample or the arrangement and the source relative to the sample.

Preferably, the system and the method can be operated either in the so-called "near field regime" or in the "Talbot- regime " .

A preferred embodiment for the grating Gl provides for Gl as a line grating being either an absorption grating or a phase grating, wherein the phase grating is a low absorption grating but generating a considerable X-ray phase shift, the latter preferably of n/2 or multiples thereof.

Accordingly, the grating G2 can be realized as a line grating having a high X-ray absorption contrast with its period being the same as that of the self image of Gl; G2 being placed in front of the detector with its lines parallel to those of Gl .

The system and the method allow for a certain freedom of operation where an operation can be chosen to be either in parallel-beam, quasi parallel, fan-beam or cone-beam mode and GO, Gl and G2 have a corresponding planar, cylindrical or spherical shape resp. Accordingly, the operation can be chosen to be either in fullfield mode with two dimensional gratings or in scanning mode with one dimensional gratings.

For the setup of the system and the method, both types of operation can be defined as follows: a) for near-field- regime operation, the distance between the gratings is chosen freely within the regime, and b) for the Talbot-regime is chosen according to

where «

= 1,3,5 and η =

L +D„,_sph p_Y

2 if the phase shift of G_l is (21 - 1) π, p₂

L 2 where / = 1,2,3 , D„ is an odd fractional Talbot distance when the parallel X-ray beam is used, while D_{n sph} is that when the fan or cone X-ray beam is used, L is the distance between the source and the Gl .

In order to benefit from the full advantages of the differential phase contrast approach, the system and the method are executed in a way that the phase stepping is performed by mechanical shift of one grating GO, Gl or G2 with respect to the others.

With respect to the grating structure, the grating structure may be advantageously manufactured by planar technology.

In order to separate the differential phase information in the image, such as a medical image from a CT, MRI or ultrasonic facility, the differential phase information may be obtained according to the European Patent application EP 10167569.2 which is herewith incorporated by reference. In a further preferred embodiment of the present invention, the phase relation between Gl and G2 can correspond exactly to the value for which the intensity curve can be expanded by a first order Taylor series and the differential phase information can be obtained preferably according to the In- ternational Patent application PCT/EP2010/051291 (WO

2010/089319) which is herewith incorporated by reference. During the solution of the mathematical cost function, the operator D_x may be designed as a differentiation operator of any order in the phase stepping direction of the gratings. Accordingly, the transform operator T may be chosen to be a differentiation operator of any order in the direction perpendicular to the stepping direction of the gratings. Further, the weighting operator W may be chosen to be a diagonal weighting matrix containing the inverse standard deviation l/cr_DPC of the DPC image in each pixel.

Preferably, the constrained optimization problem can be solved by recasting it to a second order cone program

(SOCP) . Alternatively, the constrained optimization problem may be solved by casting it to an unconstrained form possi- bly according to: minimize +∑λ,

Alternatively, the unconstrained optimization problem may be solved by means of a Gradient descent or a (non-linear) Conjugate Gradient algorithm.

Thereby, the new invention addresses the problem of stripe artifacts upon direct integration of noisy DPC images. The fundamental idea is to suppress the variations in the vertical image direction by solving a constrained optimization problem. While maintaining consistency with the measured data, the phase image is retrieved by minimizing a cost function. Applied to the case of noisy DPC measurements, this forces the integration to generate lower variations in the image and therefore improve image quality.

Preferred example of the present invention are described hereinafter in more detail. Hereafter, images, such as medical images or images from structure analysis or material testing and the like, are represented as vectors. An image vector f(i) is obtained from the column-wise extraction of pixel values in an image

where n_x x n_Y is the image size. The size of an image vector is 1 x n with n = n_x-n_y. Image transformations are represented by operators (matrices) which can be applied to an image vector. Operator matrices are of the size m x n, where m is the size of the transformed vector. In most cases, m = n holds .

The new method is based on a standard linear regression

model, where the measurement of the DPC image is given by φ = D_x ·φ +w .

D_x is the measurement operator modelling the differential measurement of φ and w is a random vector modelling the

noise in a DPC image. D_x can be implemented with a finite difference transform of φ in x-direction.

Regarding the noise model, the noise variance in a pixel a DPC image is given by

DPC (7)

NV where N is the number of photons and V is the mean fringe visibility at this pixel. N and V can be calculated by using the AC and DFI images of the measurement. For the phase retrieval from a DPC measurement, a constrained optimization problem is defined: minimize Ife/IL

p

subject to : <ε

In this optimization problem, f is an image vector, T is a transform operator, φ is the measured DPC image vector and is a boundary for the noise power. W is a diagonal weightin matrix for taking into account the noise model, with w±_ri = l/o_DPCfi. ... indicates a p-norm of a vector, which is de- fined as

The problem expressed in (8) seeks for the minimal p-norm of the vector Tf for any f which meets the data consistency constraint |^ _χ/-φ)|^ <ε .

In the minimization term of (8), the image can be transformed into any linear transform domain represented by the matrix T (e.g. Fourier transform, wavelet, finite differences, etc.) . This makes the utilization of data constraints extremely flexible. In the case of DPC measurements, a finite differences transform

T = D_yl (10) is preferable, because the integrated image is distorted by high intensity variations (horizontal stripes) in the direction perpendicular to phase stepping.

There is a variety of possibilities to solve problem (8) . In the convex case (p ≥ 1 ) , it can be recasted and solved as a second order cone program. An alternative is to recast (8) to an unconstrained form using a Lagrange multiplier, minimize ψ(ϋ_χ + , (11)

or in general, with an arbitrary amount of data constraints, minimize ψ _χί (12)

This problem formulation is also known as regularization, which has mainly been applied for the inversion of ill-posed problems. The Lagrangian multiplier (or regularization parameter) Ai controls the weighting of the regularization

term compared to the data conistency term ||^ _χ/-φ)||^ .

Problem (12) has the powerful property to allow any number of regularization terms, which is particularly useful if more a-priori knowledge about the object is available.

The choice of the norm parameter p in the regularization

term depends on the used transform operator T. A typical choice is p=2 , since in this case, problems (11)/ (12) are linear and an explicit solution exists. For the case T=D_y, the /₂-norm can lead to blurring and therefore reduce image resolution. On the other hand, the minimization of the i - norm (p=l) is well known to preserve edges in the image and is thus expected to yield a better image resolution than the

/2-norm [1] . An fi-norm minimization leads to a non-linear optimization problem where no explicit solution exists.

Here, an iterative algorithm is used to solve the non-linear optimization problem. It is based on a non-linear Conjugate Gradients (NLCG) method, which is characterized by a fast convergence for the inversion of large-scale linear systems

[2] .

References [1] L. Rudin, S. Osher, and E. Fatemi, "Nonlinear total variation based noise removal algorithms," Phys . D Nonlinear Phenom. 60, 2597268 (1992) [2] J. Nocedal and S. Wright, Numerical optimization

(Springer verlag, 1999)

Claims

1. A method for the recovery of an integrated image from a differential image by using constrained optimization;

comprising the steps of:

identifying an image that minimizes a cost function \\Tf\\_f subject to the constraint |^(¾./^"-φ)|^ , which is written as nunimize ϊ\_ρ

p

subject to : < ε

where f is an image vector, T is a transform operator matrix, D_x is a differentiation operator matrix, W a diag- noal matrix containing the reciprocal noise standard deviations, φ is the measured image vector and ε a boundary for the noise power. ... indicates the p-norm of a vec- tor, where p is a positive number, preferably an integer in the interval [0,2] .

2. A method according to claim 1, where the differential

data is obtained from an arrangement for x-rays, in par- ticular hard x-rays, for obtaining quantitative x-ray images from a sample comprising:

a. an X-ray source;

b. three or at least two gratings dubbed GO, Gl and G2 or Gl and G2.

c. a position-sensitive detector with spatially modulated detection sensitivity having a number of individual pixels;

d. means for recording the images of the detector;

e. means for evaluating the intensities for each pixel in a series of images in order to identify the characteristic of the object for each individual pixel as an absorption dominated pixel and/or a differential phase contrast dominated pixel and/or an x-ray scattering dominated pixel;

f. wherein the series of images is collected by continuously or stepwise rotating from 0 to n or 2n either the sample or the arrangement and the source relative to the sample.

The method according to claim 1 or 2, operated either in the so-called "near field regime" or in the "Talbot- regime " .

The method according to any of the preceding claims 1 to

3, wherein Gl is line grating (Gl) either an absorption grating or a phase grating, that is a low absorption grating but generating a considerable X-ray phase shift, the latter preferably of n/2 or multiples thereof.

The method according to any of the preceding claims 1 to

4, wherein G2 is a line grating having a high X-ray absorption contrast with its period being the same as that of the self image of Gl. G2 is placed in front of the detector with its lines parallel to those of Gl .

The method according to any of the preceding claims 1 to

5, wherein an operation is chosen to be either in parallel-beam, quasi parallel, fan-beam or cone-beam mode and GO, Gl and G2 have a corresponding planar, cylindrical or spherical shape resp.

The method according to any of the preceding claims 1 to

6, wherein an operation is chosen to be either in full- field mode with two dimensional gratings or in scanning mode with one dimensional gratings.

The method according to any of the preceding claims 1 to 7, wherein for near-field-regime operation, the distance between the gratings is chosen freely within the regime, and for the Talbot-regime is chosen according to

w

where / = 1,2,3 , D_n is an odd fractional Talbot distance when the parallel X-ray beam is used, while D_nsph is that when the fan or cone X-ray beam is used, L is the distance between the source and the Gl .

The method according to any of the preceding claims 1 8, wherein phase stepping is performed by mechanical shift of one grating GO, Gl or G2 with respect to the others .

10. The method according to any of the preceding claims 1 to 9, wherein the grating structure is manufactured by planar technology.

11. The method according to any of the preceding claims 1 to 10, wherein the differential phase information is obtained according to the European Patent application EP 10167569.2.

12. The method according to any of the preceding claims 1 to 11, wherein the phase relation between Gl and G2 corresponds exactly to the value for which the intensity curve can be expanded by a first order Taylor series and the differential phase information is obtained preferably according to International Patent application

PCT/EP2010/051291 (WO 2010/089319) .

13. The method according to any of the preceding claims 1 to 12, wherein the operator D_x is a differentiation operator of any order in the phase stepping direction of the gratings .

14. The method according to any of the preceding claims 1 to 13, wherein the transform operator T is a differentiation operator of any order in the direction perpendicular to the stepping direction of the gratings.

15. The method according to any of the preceding claims 1 to 14, wherein the weighting operator W is a diagonal weighting matrix containing the inverse standard devia- tion l/c¾pc of the DPC image in each pixel.

16. The method according to any of the preceding claims 1 to 15, wherein the constrained optimization problem is solved by recasting it to a second order cone program (SOCP) .

17. The method according to any of the preceding claims 1 to 15, wherein the constrained optimization problem is solved by casting it to an unconstrained form possibly according to: minimize

18. The method according to any of the preceding claims 1 to 17, wherein the unconstrained optimization problem is solved by means of a Gradient descent or a (non-linear) Conjugate Gradient algorithm or others.

19. A system for the recovery of an integrated image from a differential image by using constrained optimization; comprising

||^(£>_χ/-φ)||^ , which is written as

minimize 7/|

subject to : ≤ε

b) wherein f is an image vector, T is a transform operator matrix, D_x is a differentiation operator matrix, W a diagnoal matrix containing the reciprocal noise standard deviations, φ is the measured image vector and ε a boundary for the noise power.

II...I indicates the p-norm of a vector, where p is a positive number, preferably an integer in the interval [0,2] .

20. The system according to claim 19, where the differential image is obtained from an arrangement for x-rays, in particular hard x-rays, for obtaining quantitative x-ray images from a sample comprising:

a. an X-ray source;

b. three or at least two gratings dubbed GO, Gl and G2 or Gl and G2 resp.;

d. means for recording the images of the detector;

e. means for evaluating the intensities for each pixel in a series of images in order to identify the characteristic of the object for each individual pixel as an absorption dominated pixel and/or a differen- tial phase contrast dominated pixel and/or an x-ray scattering dominated pixel; f. wherein the series of images is collected by continuously or stepwise rotating from 0 to n or 2n either the sample or the arrangement and the source relative to the sample.

21. The system according to claim 19 or 20, being operated either in the so-called "near field regime" or in the "Talbot-regime " .

22. The system according to any of the preceding claims 19 to 21, wherein Gl is a line grating (Gl) being either an absorption grating or a phase grating, wherein the phase grating is a low absorption grating but generating a considerable X-ray phase shift, the latter preferably of n/2 or multiples thereof.

23. The system according to any of the preceding claims 19 to 22, wherein G2 is a line grating having a high X-ray absorption contrast with its period being the same as that of the self image of Gl; G2 being placed in front of the detector with its lines parallel to those of Gl .

24. The system according to any of the preceding claims 19 to 23, wherein an operation is chosen to be either in parallel-beam, quasi parallel, fan-beam or cone-beam mode and GO, Gl and G2 have a corresponding planar, cylindrical or spherical shape resp.

25. The system according to any of the preceding claims 19 to 24, wherein an operation is chosen to be either in fullfield mode with two dimensional gratings or in scanning mode with one dimensional gratings.

26. The system according to any of the preceding claims 19 to 25, wherein for near-field-regime operation, the dis- tance between the gratings is chosen freely within the regime, and for the Talbot-regime is chosen according to

w

27. The system according to any of the preceding claims 19 to 26, wherein phase stepping is performed by mechanical shift of one grating GO, Gl or G2 with respect to the others .

The system according to any of the preceding claims 19 to 27, wherein the grating structure is manufactured by planar technology.

29. The system according to any of the preceding claims 19 to 28, wherein the differential phase information is obtained according to the European Patent application EP 10167569.2.

30. The system according to any of the preceding claims 19 to 29, wherein the phase relation between Gl and G2 corresponds exactly to the value for which the intensity curve can be expanded by a first order Taylor series and the differential phase information is obtained preferably according to International Patent application

PCT/EP2010/051291 (WO 2010/089319) .

31. The system according to any of the preceding claims 19 to 30, wherein the operator D_x is a differentiation operator of any order in the phase stepping direction of the gratings .

32. The system according to any of the preceding claims 19 to 31, wherein the transform operator T is a differentiation operator of any order in the direction perpendicular to the stepping direction of the gratings.

33. The system according to any of the preceding claims 19 to 32, wherein the weighting operator W is a diagonal weighting matrix containing the inverse standard deviation l/c¾pc of the DPC image in each pixel.

34. The system according to any of the preceding claims 19 to 33, wherein the constrained optimization problem is solved by recasting it to a second order cone program (SOCP) .

35. The system according to any of the preceding claims 19 to 33, wherein the constrained optimization problem is solved by casting it to an unconstrained form possibly according to: minimize

36. The system according to any of the preceding claims 19 to 35, wherein the unconstrained optimization problem is solved by means of a Gradient descent or a (non-linear)

Conjugate Gradient algorithm.