Separation of superimposed images with subpixel shift

Jailin, C.; Poncelet, M.; Roux, S.

doi:10.1107/S1600577517015892

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 25| Part 1| January 2018| Pages 272-281

https://doi.org/10.1107/S1600577517015892

Separation of superimposed images with subpixel shift

Clément Jailin,^a ^* Martin Poncelet ^a and Stéphane Roux ^a

^aLMT (ENS Cachan/CNRS/Université Paris-Saclay), 61 avenue du Président Wilson, F-94235 Cachan, France
^*Correspondence e-mail: clement.jailin@lmt.ens-cachan.fr

Edited by V. Favre-Nicolin, CEA and Université Joseph Fourier, France (Received 10 July 2017; accepted 1 November 2017)

The problem of the separation of superimposed images is considered in the particular case of a steady background and a foreground that is composed of different patterns separated in space, each with a compact support. Each pattern of the foreground may move in time independently. A single pair of these superimposed images is assumed to be available, and the displacement amplitude is typically smaller than the pixel size. Further, assuming that the background is smoothly varying in space, an original algorithm is proposed. To illustrate the performance of the method, a real test case of X-ray tomographic radiographs with moving patterns due to dust particles or surface scratches of optical elements along the beam is considered. Finally an automatic and simple treatment is proposed to erase the effects of such features.

Keywords: superimposed images; image mixture; image separation; spatial shift; subpixel motion.

1. Introduction

X-ray computed tomography is an imaging technique that provides a 3D image of a specimen from a set of 2D radiographs acquired during its rotation between an X-ray beam and a detector. The reduction in beam intensity at a specific detector site is related to the integral of the local attenuation coefficient along the ray (the line parallel to the beam) passing through that site (i.e. the Beer–Lambert law). Hence, the radiographs have to be normalized by the incident beam intensity, called `flat-fields', i.e. radiographs that would be obtained without a specimen. One practical source of difficulty (Weitkamp et al., 2011 ) is that the beam intensity displays temporal variations and spatial inhomogeneities, especially at synchrotron facilities (Flot et al., 2010 ) [as shown in Fig. 1(a) for a raw radiograph]. Scratches or dust particles along the beam pathway (optics, monochromator, scintillator, mirror, camera) may become visible on the radiographs. They also appear in the flat-fields acquired before and after the experiment [Fig. 1(b)]. Should those defects remain ideally still with respect to the acquisition then they would have no consequence in the reconstructions. However, it is observed that some localized patterns move with respect to the background [Fig. 1(c)]. Although their displacements are small, they are very salient and hence may cause artifacts. Moreover, even if in some cases these moving patterns on the radiographs may not be detected after the reconstruction algorithm, they would still introduce bias for qualitative analyses based on the radiographs (Leclerc et al., 2015 ). Such a real experimental case where such a differential pattern motion takes place will be discussed below, and will be used in §4 to validate our procedure.

Figure 1
(a) Raw radiograph of a cast iron sample with rectangular cross section where the previous inhomogeneities are clearly observed and should be erased by a flat-field normalization provided this background remains temporally steady. (b) Flat-field, signature of the beam inhomogeneity. (c) A flat-field difference in the rectangle drawn in (b) showing a moving particle. The color scale for this residual is expressed in proportion to the flat-field amplitude to provide an order of magnitude.

The separation of image superimpositions with motion into separate layers has received much attention in the past years. Many papers deal with the segmentation and recognition of moving patterns (for video surveillance, pattern recognition, etc.) (Kameda & Minoh, 1996 ; Lipton et al., 1998 ; Rogers et al., 2007 ; Ng & Delp, 2010 ). They consist of segmenting a moving image with respect to a fixed background.

The proposed procedures often make use of the image difference to separate the different layers. Because moving objects are not transparent and mask the scene, large displacements are often considered so that, once the front object shape has been delineated, the foreground and background layers may be separated with binary masks.

Transparent image mixtures with motion are more challenging. Such cases occur, for example, for moving shadows, X-ray images, superimposed semi-transparent objects, partial reflection onto a transparent surface, etc. (Irani & Peleg, 1993 ; Bergen et al., 1990 ; Auvray et al., 2006 ). With transparent images, Be'Ery & Yeredor (2008 ) and Gai et al. (2008 , 2012 ) proposed a separation method based on the inter-correlation of the image mixture (image gradients may be considered to enhance the signal). Remarkable success has been reported in terms of separation; however, it is to be noted that such results require specific conditions for the algorithm to be applicable: images have to be uncorrelated, gradients have to be sparse and with a similar weight, the displacement has to be larger than the correlation length in order to properly separate correlation peaks. In particular, this technique cannot address cases where the displacement is subpixel. Other works aim to separate the background and foreground for a known displacement field (Toro et al., 2003 ; Marcia et al., 2008a ,b ) and wavelet decomposition for the two latter papers. Again, such methods are not suited to dealing with very small displacements. Thus, in spite of the variety of the above-cited problems and algorithms, none appears suited to very small displacements such as those encountered in practice for flat-field corrections in computed tomography.

Ramirez-Manzanares et al. (2007 , 2010 , 2011 ) proposed a method to estimate the displacement field of superimposed transparent images subjected to a small motion. This procedure is based on a variational model for integrating local motion estimations to obtain a multi-valued velocity field. Similar methods were proposed by Stuke et al. (2004 ). However, these procedures assume that some tens of images are available, in between which the velocity is steady. When reduced to a single pair of images, the problem becomes ill-posed and additional assumptions are to be formulated and exploited.

The goal of the present study is to erase moving particles on fixed X-ray transparent images by developing a robust method to separate a small set of superimposed images (i.e. from two different flat-fields) composed of an extended fixed background (with low and high spatial frequencies) and a localized moving foreground (spatial medium frequencies) with the specificity of handling very small unknown displacement amplitudes (subpixel displacement). After having defined our notations in §2, the method to extract the displacement, the background and the foreground with a spatial frequency separation is presented in §3. Then, in §4, the procedure is tested on two X-ray scans. Finally, the correction of other images of the same set and with the extracted results is proposed in §5.

2. Statement of the problem and notations

2.1. Images

A digital image is a collection of gray-level values $[\hat{f}_{i}]$ for each pixel #!i whose centre is at position $[{\bf{r}}_{i}]$ of integer coordinates. Although discrete, it may be seen as the sampling at integer positions $[{\bf{r}}_{i}]$ of an underlying function $[f({\bf{r}})]$ defined for arbitrary real coordinates $[{\bf{r}}]$ = #!(x,y) in a domain Ω.

Because only $[\hat{f}]$ is known and not #!f, an interpolation scheme is proposed to compute an approximation of #!f at any arbitrary point from the knowledge of $[\hat{f}]$ (or #!f at pixel coordinates). Different interpolants with different degrees of regularity can be proposed to extend $[\hat{f}]$ to the continuum. The interpolation scheme relies on a kernel $[h({\bf{r}})]$ defined in the continuum, valued 1 at the origin and 0 at all other integer coordinate points, so that

$[f({\bf{r}})=\sum\limits_{{i\,\in\,{{\cal I}}}} h({\bf{r}}-{\bf{r}}_{i})\,\hat{f}_{i}, \eqno(1)]$

where the sum runs over the set $[{\cal I}]$ of all pixels in the image, i.e. pixel centres, and hence for all pixels $[{\bf{r}}_{i}]$ , $[f({\bf{r}}_{i})]$ = $[\hat{f}_{i}]$ .

Registration of two images captured after a translation of small amplitude can provide estimates of the displacement amplitude with an uncertainty well below $[10^{{-2}}]$ pixel size [down to $[10^{{-3}}]$ pixel size under favourable conditions, as widely documented in the literature (see, for example, Schreier et al., 2009 )].

Although this may appear as paradoxical, this property rests on the ability to propose a subpixel interpolation scheme for images that is more accurate than what can be detected by the noise level of the image. It is noteworthy to add that this property is of statistical nature, and results from average properties over large zones of interest. As their size is reduced, displacement uncertainties increase significantly.

Because #!f can be deduced from $[\hat{f}]$ and vice versa, no difference is made in the following between the discrete image and its extension to the continuum. In a similar spirit, a #!C¹ interpolation scheme (or higher) allows a gradient to be defined, that will be used freely in the following.

2.2. Definition of the problem

An image $[f_{0}({\bf{r}})]$ with Cartesian coordinates $[{\bf{r}}]$ is assumed to be the superimposition of a background $[\varphi({\bf{r}})]$ and a localized pattern $[\psi({\bf{r}})]$ having a compact support (e.g. a dust particle or a surface scratch of a transparent object encountered along the beam),

$[f_{0}({\bf{r}})=\varphi({\bf{r}})+\psi({\bf{r}}).\eqno(2)]$

In addition to #!f₀, a second image is available where the background remains steady, but the localized pattern is translated by some unknown displacement both in orientation and magnitude. The latter will be assumed to be small in the following,

$[f_{1}({\bf{r}})=\varphi({\bf{r}})+\psi({\bf{r}}+{\bf{u}}_{1}).\eqno(3)]$

Moreover, over the entire image, $[{\bf{u}}_{1}]$ is assumed to be uniform (i.e. independent of $[{\bf{r}}]$ ). The two images #!f_n, for #!n = 0, 1, are known but φ and ψ are not. Similarly, the displacement $[{\bf{u}}_{1}]$ is not known; one can conventionally choose $[{\bf{u}}_{0}]$ = 0.

The fact that the moving object is a localized pattern with compact support means that $[\psi({\bf{r}})]$ = 0 away from a region that can be easily circumscribed. Let us note that one may not only encounter X-ray absorption but also phase contrast effects that redirect the beam locally. Hence one may observe an increase of intensity, not only an attenuation, so that a positivity constraint on ψ is not considered. It is also necessary to mention some further property of the background φ. The latter displays both long-wavelength modulation and high-frequency noise. However, it is assumed to be statistically stationary. In other words, if ψ is subtracted from, say, #!f_n, the prior localization of ψ should not be visible in φ. This is to be contrasted with #!f_n where the presence of a specific pattern is manifest.

3. Method

3.1. Ill-posed problem

Without additional assumptions on φ and ψ, the problem is ill-posed, in the sense that its solution is not unique, as shown below.

Let us note that the differences between two frames,

$[\eqalign{ d_{1}({\bf{r}}) & \equiv f_{1}({\bf{r}})-f_{0}({\bf{r}}) \cr& = \psi({\bf{r}}+{\bf{u}}_{1})-\psi({\bf{r}}),} \eqno(4)]$

only depends on ψ and no longer on φ. To simplify the notations, for difference, #!d, and displacement vector, $[{\bf{u}}]$ , the subscript #!n or #!1 is dropped.

The pattern ψ appears to be determined from its finite difference if $[{\bf{u}}]$ is known. Assuming that ψ is null over a strip of width $[{\bf{u}}]$ , integration is straightforward and leads to a particular solution to equation (4). However, any periodic function, of period $[{\bf{u}}]$ , can be added to ψ. Once ψ is determined, then φ is obtained by a mere difference [see equation (2)]. Therefore, in addition to a constant offset as mentioned by Szeliski et al. (2000 ), any periodic function of period $[{\bf{u}}]$ can be added to a particular solution and still fulfills exactly equation (3) for any #!f₀ and #!f₁. Thus, the problem is ill-posed with a large degeneracy.

To make the problem well-posed, some additional constraints have to be prescribed. In our case, the property that ψ is a localized pattern with a compact support provides a simple way to determine the unknown periodic function in order to minimize the power of ψ outside a domain where it is assumed to be non-zero. Thus for any $[{\bf{u}}]$ the degeneracy is reduced to a single solution that achieves the best score in matching the known difference #!d. The remaining question is whether one can determine independently the displacement $[{\bf{u}}]$ .

3.2. Determination of the displacement orientation

Let us consider the Radon transform of #!d, $[\rho\equiv R[d]]$ , given by the line sum of #!d along parallel lines at an angle θ, or

$[\rho(s,\theta)=\int\limits_{{{\cal C}}} d\left(s{\bf{e}}_{{\theta+\pi/2}}+t{\bf{e}}_{{\theta}}\right)\,{\rm d}t, \eqno(5)]$

where $[{\bf{e}}_{\theta}]$ is a unit vector forming an angle θ with the #!x axis, and $[{\cal C}]$ denotes a circular disk of radius #!a containing the domain where #!d is non-zero.

Let us call α the polar angle of the displacement vector, $[{\bf{u}}]$ = $[u{\bf{e}}_{\alpha}]$ , then

$[\eqalignno{ \rho(s,\alpha) & = \textstyle\int\limits_{{{\cal C}}} d(s{\bf{e}}_{{\alpha+\pi/2}}+t{\bf{e}}_{{\alpha}})\,{\rm d}t \cr& = \textstyle\int\limits_{{{\cal C}}} \psi\left[s{\bf{e}}_{{\alpha+\pi/2}}+(t+u){\bf{e}}_{{\alpha}}\right]\,{\rm d}t \cr& \quad -\textstyle\int\limits _{{{\cal C}}} \psi(s{\bf{e}}_{{\alpha+\pi/2}}+t{\bf{e}}_{{\alpha}})\,{\rm d}t \cr& = R[\psi](s,\alpha)-R[\psi](s,\alpha) \cr& = 0 &(6) }]$

for all values of #!s. Thus the #!L₂ norm of $[\rho(.,\theta)]$ , $[l(\theta)]$ ,

$[l(\theta)^{2}\equiv\textstyle\int\limits_{{-r}}^{r}\rho(s,\theta)^{2}\,{\rm d}s,\eqno(7)]$

is a positive function that should reach its minimum (ideally 0) in the direction of the motion $[\theta]$ = $[\alpha]$ . This property is interesting as it allows the determination of the displacement orientation without any further assumption than the compact support of ψ. In order not to bias this criterion along preferred directions, it is natural to clip the integration domain for the Radon transform and the #!L₂ norm to a circular domain centered on the pattern.

3.3. Ill-posedness

Let us note that even if the displacement direction can be inferred, the magnitude itself remains an issue. The problem is still ill-posed as can be easily observed in the limit of a small displacement (a limit that is relevant for our application). Indeed, in such a case, the difference between images can be written

$[\eqalignno{ d({\bf{r}}) & = \psi({\bf{r}}+{\bf{u}})-\psi({\bf{r}}) \cr & \approx |{\bf{u}}|{\boldnabla}\psi({\bf{r}})\cdot{\bf{e}}_{\alpha}. &(8)}]$

Therefore, knowing #!d and $[{\bf{e}}_{\alpha}]$ gives access to the product of the displacement magnitude by the component of the gradient along the displacement direction, but with this sole argument it does not allow $[|{\bf{u}}|]$ to be isolated from $[{\boldnabla}\psi\cdot{\bf{e}}_{\alpha}]$ . Hence the problem remains ill-posed, and the criterion on the compact support of ψ does not help.

3.4. First step

Let us first propose the partial reconstruction of ψ, from the above observations. Because, at this stage, it is not possible to split displacement magnitude and pattern gradient, we propose to compute the pattern for a chosen value of the displacement magnitude, here one pixel. This value is by now conventional and its determination will be discussed later. In order not to introduce any confusion in the latter quantity and ψ, and because it is related to first-order integration of #!d, we call this integral #!D. It solves the following equation,

$[d({\bf{r}})={\boldnabla}D({\bf{r}})\cdot{\bf{e}}_{\alpha}.\eqno(9)]$

Later, when the displacement magnitude is known, ψ will be estimated as

$[\psi({\bf{r}})={{D({\bf{r}})} \over {|{\bf{u}}|}}.\eqno(10)]$

3.5. Integration

Although the problem is now well-posed, there are several ways to implement the integration of #!D numerically. Let us recall that the displacement $[{\bf{u}}]$ may be subpixel, or, even if its magnitude is one pixel as chosen above, the arbitrary angle α requires that a subpixel interpolation scheme be available.

Thus, $[D({\bf{r}})]$ is chosen under the following form,

$[D({\bf{r}})=\textstyle\sum\limits_{{i\,\in\,{\cal I}}} b_{i}\,h\left({\bf{r}}-{\bf{r}}_{i}\right),\eqno(11)]$

where the index #!i runs over the set $[{\cal I}]$ of all pixels in the image, $[{\bf{r}}_{i}]$ designates the coordinate of pixel #!i, #!b_i is the unknown amplitude and #!h is the elementary shape function relative to a pixel centered at the origin. In the following, the interpolation scheme is inspired from finite-element, with a bilinear interpolation. In this case, the shape function suited to the square lattice of pixels, at any arbitrary real point of coordinates $[{\bf{r}}]$ = #!(r_x,r_y), is

$[h({\bf{r}})=\big(1-|r_{x}|\big)\left(1-|r_{y}|\right)\eqno(12)]$

if $[|r_{x}|]$ < 1 and $[|r_{y}|]$ < 1, and otherwise $[h({\bf{r}})]$ = 0. It is classically referred to as Q4P1 (4-noded quadrilateral, polynomial of order 1).

The determination of #!D is to be performed from the minimization of

$[\chi _{1}^{2}(\{ b\}) = \int\limits_{{{\cal C}}} \left\{ \textstyle\sum\limits_{i} b_{i} \left[ h({\bf{r}}-{\bf{r}}_{i}+{\bf{e}}_{\alpha})-h({\bf{r}}-{\bf{r}}_{i})\right]-d({\bf{r}})\right\}^{2}\,{\rm d}{\bf{r}}.\eqno(13)]$

If some additional information is available concerning ψ, one may choose, instead of a pixel representation, #!h, a basis that is tailored to the expectation. The interest of introducing such a form is that, because of noise and subpixel interpolation, the line sum $[\rho(s,\alpha)]$ is not strictly 0 over the domain (disk) of integration. In this case, the above minimization allows the distribution of the additional weight on both sides of the pattern, whereas a direct integration would lead to a dissymmetric #!D.

3.6. Second step

The criterion to find $[|{\bf{u}}|]$ is to assume that if the object is removed from the background then no ghost mark (neither positive nor negative) should appear on the `computed' background. However, one difficulty is that this background itself is unknown. As discussed above, for any $[{\bf{u}}]$ , one can compute $[\psi_{{{\bf{u}}}}({\bf{r}})]$ (where the subscript recalls that this pattern estimate depends on $[{\bf{u}}]$ that remains to be determined). In turn, the background is $[\varphi_{{{\bf{u}}}}({\bf{r}})]$ = $[f_{0}({\bf{r}})-\psi_{{{\bf{u}}}}({\bf{r}})]$ .

A spatial frequency separation is performed in the proposed approach. The medium frequencies (5–20 pixels) are considered as moving patterns whereas the low frequencies are related to the background and high frequency (1 pixel) to noise. As a way to estimate the long-wavelength modulation of the background over a region Ω, it is proposed to perform a least-squares fit weighted by a function $[w({\bf{r}})]$ that is null over the expected support of ψ and non-zero in its surrounding. A set of slow modulation functions, $[g_{i}({\bf{r}})]$ , with i = 1,…, N_f, is introduced so that

$[{\bf{a}}=\mathop{}\!{\rm argmin}\left\{\,\int\limits_{\Omega}w({\bf{r}})\left[\varphi _{{{\bf{u}}}}({\bf{r}})-\sum _{{i\,=\,1}}^{{N_{f}}}a_{i}g_{i}({\bf{r}})\right]^{2}\,{\rm d}{\bf{r}}\right\}.\eqno(14)]$

The least-squares solution is obtained from

$[a_{{i}}=M^{\,{-1}}_{{ij}}\,s_{{j}},\eqno(15)]$

where

$[M_{{ij}}=\textstyle\int\limits_{\Omega}w({\bf{r}})\,g_{i}({\bf{r}})\,g_{j}({\bf{r}})\,{\rm d}{\bf{r}}\eqno(16)]$

and

$[s_{{i}}=\textstyle\int\limits_{\Omega}w({\bf{r}})\,g_{i}({\bf{r}})\,\varphi_{{{\bf{u}}}}({\bf{r}})\,{\rm d}{\bf{r}}.\eqno(17)]$

It is to be noted that this expression can be rewritten as

$[\eqalign{ a_{i} & = \textstyle\int\limits_{\Omega}\gamma_{i}({\bf{r}})\,\varphi_{{{\bf{u}}}}({\bf{r}})\,{\rm d}{\bf{r}}, \cr \gamma_{i}({\bf{r}}) & = w({\bf{r}})\,M^{\,{-1}}_{{ij}}\,g_{j}({\bf{r}}),} \eqno(18)]$

and hence the remainder is

$[\eqalignno{ \widetilde{\varphi}_{{{\bf{u}}}}({\bf{r}}) & = \varphi_{{{\bf{u}}}}({\bf{r}})-\left[\,\textstyle\int\limits_{\Omega}\gamma_{i}({\bf{r}}^{{\prime}}) \,\varphi_{{{\bf{u}}}}({\bf{r}}^{{\prime}})\,{\rm d}{\bf{r}}^{{\prime}}\right]\,g_{i}({\bf{r}}) \cr& = \textstyle\int\limits_{\Omega}\left[\delta({\bf{r}}-{\bf{r}}^{{\prime}})-\gamma _{i}({\bf{r}}^{{\prime}})\,g_{i}({\bf{r}})\right]\,\varphi _{{{\bf{u}}}}({\bf{r}}^{{\prime}})\,{\rm d}{\bf{r}}^{{\prime}}, &(19)}]$

where the above writing illustrates that $[\widetilde{\varphi}_{{{\bf{u}}}}]$ is related to $[\varphi_{{{\bf{u}}}}]$ through a linear operator (projector),

$[\widetilde{\varphi}_{{{\bf{u}}}}({\bf{r}})= \textstyle\int\limits_{\Omega} P({\bf{r}},{\bf{r}}^{{\prime}})\, \varphi_{{{\bf{u}}}}({\bf{r}}^{{\prime}})\,{\rm d}{\bf{r}}^{{\prime}}, \eqno(20)]$

where

$[P({\bf{r}},{\bf{r}}^{{\prime}})\equiv\delta({\bf{r}}-{\bf{r}}^{{\prime}})-\gamma _{i}({\bf{r}}^{{\prime}})\,g_{i}({\bf{r}}). \eqno(21)]$

The above linear regression is expected to capture the slow modulation of the background and hence, for the appropriate value of the displacement magnitude $[|{\bf{u}}|]$ , $[\widetilde{\varphi}_{{{\bf{u}}}}({\bf{r}})]$ should only consist of white noise. It is therefore proposed to estimate $[|{\bf{u}}|]$ from the minimization of $[\chi_{2}^{2}(|{\bf{u}}|)]$ = $[(1/2)\|\widetilde{\varphi}_{{{\bf{u}}}}\|^{2}]$ .

Using

$[\eqalignno{ \varphi _{{{\bf{u}}}}({\bf{r}}) & = f_{0}({\bf{r}})-\psi_{{{\bf{u}}}}({\bf{r}}) \cr& = f_{0}({\bf{r}})-{{D({\bf{r}})} \over {|{{\bf{u}}}|}}, &(22)}]$

we can write the stationarity condition as

$[{{\partial\chi_{2}^{2}(|{\bf{u}}|)} \over {\partial|{\bf{u}}|}} = {{\partial\chi_{2}^{2}(|{\bf{u}}|)} \over {\partial\varphi _{{{\bf{u}}}}}}\,{{\partial\varphi_{{{\bf{u}}}}} \over {\partial|{\bf{u}}|}} = 0. \eqno(23)]$

Because the piece-wise linear interpolation of the discrete values of pixels is used, the integration over the domain Ω can be written as a discrete sum. Hence the previous equation can be expressed in matrix notation,

$[\left(f_{0}-{{D} \over {|{\bf{u}}|}}\right)^{\!\!\top} PP^{\top}{{D} \over {|{\bf{u}}|^{2}}}=0, \eqno(24)]$

with $[|{\bf{u}}|]$ a scalar homogeneous displacement, #!f₀ and #!D vectors composed of the $[N_{\rm{pix}}]$ pixels of the integration domain, the superscript #!^T denotes transposition and where the projection matrix #!P is of size $[[N_{\rm{pix}}\times N_{\rm{pix}}]]$ . Thus, finally, we arrive at the expression of the displacement modulus,

$[|{\bf{u}}|={{D^{\top}P^{\top}PD} \over {f_{0}^{\top}P^{\top}PD}}.\eqno(25)]$

The pattern shape can finally be obtained using equation (10), and the background from the difference.

In the case of a large displacement $[{\bf{u}}]$ typically larger than the correlation length or the spot length for a compact image, the linear approximation of equation (8) is meaningless and thus $[d({\bf{r}})]$ cannot be approximated by the derivative of ψ. The displacement can be decomposed into a first-guess displacement $[{\bf{u}}_{0}]$ (e.g. obtained from a previous computation in an iterative scheme or a maximum cross-correlation) and a correction $[\delta{\bf{u}}]$ , hence $[{\bf{u}}]$ = $[{\bf{u}}_{0}]$ + $[\delta{\bf{u}}]$ . The above procedure can be extended to estimate the displacement correction $[\delta{\bf{u}}]$ . In this case, the problem becomes affine in $[\delta{\bf{u}}]$ (rather than linear as discussed above) as a result of

$[d({\bf{r}})\simeq\psi({\bf{r}}+{\bf{u}}_{0})-\psi({\bf{r}})+\nabla\psi\left({\bf{r}}+{\bf{u}}_{0}\right)\delta{\bf{u}}.\eqno(26)]$

4. Case study

4.1. Test case presentation

To illustrate the performance of the proposed method, a real tomographic acquisition performed at the ESRF synchrotron, beamline ID19, is chosen. A nodular graphite cast iron is imaged at a voxel size of 5.1 µm. A complete scan corresponds to 600 radiographs, and flat-fields are recorded before the scan and after every 100 radiographs, resulting in $[N_{\rm f}]$ = 7 flat-fields. The radiographs are denoted as $[I({\bf{r}},t)]$ and the flat-fields as $[f_{i}({\bf{r}})]$ , i = 1,…, N_f.

The noise is supposed to be white and Gaussian with a standard deviation of 0.24% as assessed from regions in the images where no mobile patterns are present and thus image difference is assumed to be essentially due to noise.

The flat-field correction procedure with stationary intensity correction has been presented in detail by Jailin et al. (2017 ). The principle can be summarized as follows: taking into account the multiplicative nature of the corrections, the logarithm of the radiographs $[G({\bf{r}},t)]$ = $[\log[I({\bf{r}},t)]]$ , i.e. hereafter called the `projections', should first be computed and the logarithm of the beam intensity should be subtracted off. The standard Beer–Lambert law relates the sum of the absorption coefficient along the ray hitting the detector at position $[{\bf{r}}]$ and time #!t to $[s({\bf{r}},t)]$ = $[G({\bf{r}},t)-F({\bf{r}},t)]$ where $[F({\bf{r}},t)]$ = $[\log[\,f({\bf{r}},t)]]$ and $[f({\bf{r}},t)]$ is the flat-field at time #!t. If the beam were steady, $[F({\bf{r}},t)]$ should be equal to $[F_{i}({\bf{r}})]$ = $[\log[\,f_{i}({\bf{r}})]]$ , for all #!i. However, in real life, $[F({\bf{r}},t)]$ varies in time as the X-ray beam is not steady. Since it is not possible to measure simultaneously #!I and #!F, one has to estimate $[F({\bf{r}},t)]$ from $[F_{i}({\bf{r}})]$ .

Because the edges of the radiographs [ $[\boldrho\in\Omega _{l}\cup\Omega _{r}]$ , defined by the two rectangles shown in Fig. 2(a)] are never masked by the scanned specimen they provide information about the intensity variation at all instants #!t. The two sub-images from the projections clipped to these edge regions at time #!t can be approximated as a linear combination of the corresponding sub-images extracted from the (logarithm of the) $[N_{\rm f}]$ flat-fields $[F_{i}(\boldrho)]$ using least-squares [some extra fields can be included in this database of all flat-fields (hereafter called the library), as will be the case in §5 with the addition of a flat-field gradient],

$[G({\boldrho},t)=F({\boldrho},t)=\textstyle\sum\limits_{{i\,=\,1}}^{{N_{\rm f}}} \beta_{i}(t)\,F_{i}({\boldrho}).\eqno(27)]$

Because the flat-fields are known over the entire detector, the above decomposition allows the extension of the logarithm of the raw beam intensity $[F({\bf{r}},t)]$ to any $[{\bf{r}}]$ .

Figure 2
(a) Flat-field #!F₁ with the edge areas defined by the two rectangles. (b) Flat-field #!F₂. (c) Differences between the two weighted flat-fields where the five largest moving spots are shown in their bounding box.

This procedure provides an accurate correction in most of the detector area. Nevertheless, a few features remain visible because the radiographs are polluted by moving dust particles or surface scratches of optical elements along the beam.

The flat-fields acquired before and after the experiment and those used for the flat-field library #!F are shown in Figs. 2(a) and 2(b). These images are composed of a low-frequency background (i.e. mostly a vertical gradient), a high-frequency noise and a few medium-frequency bright spots. The difference, Fig. 2(c), shows:

(i) A stationary background composed of low frequencies and a few spots with negative values. The intensities of these spots are different from that of the background and are thus very clearly visible in the difference. However, because of the (logarithm of the) flat-field sampling, $[N_{\rm f}]$ = 7 linear combinations can capture those spots and the current procedure of allowing any combination of fields allows us to account for them.

(ii) Other patterns composed of positive and negative values (labeled from 1 to 5). These patterns are moving spots and remain visible in the corrected radiographs. They can be automatically selected by a mere thresholding procedure.

The following study will successively extract the shape of each of these spots, focusing on areas around each moving pattern (Fig. 2c) with two images weighted to have a zero-mean value on these regions. The difference of two flat-fields zoomed in the first region is shown Fig. 3. A diverging color bar is used for residual maps to highlight positive and negative patterns. It can be seen that the background disappears and a positive and negative pattern due to the moving spot in a vertical direction remains. The previous procedure can be applied to extract the spot shape of these areas.

Figure 3
(a) Bright spot of the first flat-field #!F₁. (b) Bright spot of the second flat-field #!F₂. (c) Differences of the two images with a zero-mean value. The circle shows the boundary of the disk used for the Radon transform.

4.2. Separation of the two images

The proposed procedure has been first applied on spot 1 [see Fig. 2(c)] and the results for the other spots follow. The extracted area 1 is shown in Fig. 3 for the initial (a) and the final (b) flat-fields.

The first step is to obtain the direction of the displacement by using the Radon transform. The difference image, #!d, is clipped to a disk $[{\cal C}]$ [i.e. for the surface lying outside the circle shown in Fig. 3(c) #!d is set equal to 0], the Radon transform $[\rho=R[d]]$ is computed and the mean quadratic intensity for every angle, $[l(\theta)^{2}]$ , is computed as in equation (7). The results are displayed in Fig. 4. The minimum obtained in Fig. 4(b) gives the direction of the displacement: $[\theta]$ ≃ 80° from the horizontal x-axis.

Figure 4
(a) Radon transform (i.e. projection at every angle). (b) Mean quadratic amplitude of the Radon transform, $[l(\theta)^{2}]$ . The minimum, reached for $[\theta]$ ≃ 80°, indicates the displacement direction.

Because the background is essentially composed of low frequencies, it can be extracted with the projector #!P using low-order polynomials (up to second order) $[g_{i}({\bf{r}})]$ . This projection has to be weighted by $[w({\bf{r}})]$ to not be affected by the bright spot shape. The indicator function of the complement to the disk used for the Radon transform was chosen for the studied spot. For elongated marks such as the one labeled `4', a rectangle was selected.

The computation of $[D({\bf{r}})]$ composed of #!26×32 pixels is regularized by reverting to a square-shaped (Q4) mesh composed of #!10×15 nodes. Equation (25) allows the estimation of the displacement amplitude as $[|{\bf{u}}|]$ = 0.54 pix = 2.72 µm. The foreground $[\psi_{{{\bf{u}}}}]$ [equation (10)] and finally the background $[\varphi_{{{\bf{u}}}}]$ are evaluated. The $[D({\bf{r}})]$ results linked to the shape of the pattern may be sensitive to the high-frequency noise. In our case, the spot is composed of medium frequencies so $[D({\bf{r}})]$ can be regularized by a finite-element mesh composed of a small number of degrees of freedom (the choice of a mesh based on elements larger than the pixel size leads to a lower sensitivity to noise).

The procedure can be summarized in the following four steps:

(i) Radon transform for the measurement of the motion direction [equation (7)] to evaluate the direction of motion.

(ii) Computation of #!D and #!P.

(iii) Measurement of $[|{\bf{u}}|]$ [equation (25)].

(iv) Foreground and background $[\psi_{{{\bf{u}}}}]$ and $[\varphi_{{{\bf{u}}}}]$ extraction [equation (10)].

The different images, #!F, $[\psi_{{{\bf{u}}}}]$ and $[\varphi_{{{\bf{u}}}}]$ , are shown in Fig. 5.

Figure 5
Analysis of the pattern shown in Fig. 3

with (a) original image #!F₁, (b) reconstructed spot $[\psi_{{{\bf{u}}}}({\bf{r}}_{i})]$ = $[D({\bf{r}}_{i})/|{\bf{u}}|]$ , and (c) corresponding background $[\varphi_{{{\bf{u}}}}]$ for $[|{\bf{u}}|]$ = 0.54 pixel with a #!10×15 Q4 mesh regularization.

Results for the four other proposed spots of Fig. 2, #!F₁, $[F_{1}-F_{2}]$ , $[\psi_{{{\bf{u}}}}]$ and $[\varphi_{{{\bf{u}}}}]$ , are shown, respectively, in Figs. 6(a), 6(b), 6(c) and 6(d). It can be seen that, for the fourth pattern, a `ghost mark' still remains on the background but the general shape of the pattern is well captured.

Figure 6
Results for the four other spots. The four rows correspond, respectively, to spots 2 to 5. Column (a): initial image #!F. Column (b): difference image with positive and negative pattern due to motion. Column (c): reconstructed spot $[\psi_{{{\bf{u}}}}]$ . Column (d): reconstructed background $[\varphi_{{{\bf{u}}}}]$ .

5. Application to automatic flat-field correction

The previous section showed that it is possible to segment each pattern, estimate its shape and intensity as well as its displacement from one image to another one. Yet, this treatment requires a number of manual operations, from the detection of each feature to its analysis. Although feasible, it is a lengthy treatment that appears as affordable only for exceptionally important data. Thus a key question is whether one can benefit from such an analysis to remove artefacts from tomographic projections at a much lower cost. It turns out that the answer is positive, at least when the displacement range remains small.

The key observation relies on the fact that the absolute shape of the pattern is not strictly needed, nor is the displacement amplitude. Only an estimate of the pattern gradient is necessary, and, for the component of the gradient along the translation direction, this is part of flat-field differences. However, because the motion is different for the patterns and the background, it is necessary to limit those differences to the immediate surrounding of the pattern. In fact, the creation of such a binary mask $[\mu({\bf{r}})]$ , valued 1 around each spot and 0 elsewhere, is easily performed from standard image analysis techniques as the patterns are isolated salient objects in each flat-field image. Then masked flat-field differences $[[f_{i}({\bf{r}})-f_{0}({\bf{r}})]\,\mu({\bf{r}})]$ directly contain (i.e. without further treatment) an estimate of the pattern gradient along the direction of motion. Such a field is shown in Fig. 7. Hence, in the same way as in §4.1, by minimizing the quadratic difference of any projection with a linear combination of flat-fields and, in addition, of the masked differences it should be possible to account for both the background correction and spurious mobile features. This procedure is in line with the proposed scheme to estimate the raw beam intensity, and turns out to be a very simple extension. In other words, the flat-field gradient has to be added to the flat-field library $[F_{i}({\bf{r}})]$ .

Figure 7
Additional flat-field gradient component along the motion direction masked over the five spots. The rectangles represent the left and right edges never masked by the scanned sample.

Note, however, that this procedure relies on the fact that only the motions sampled by the different flat-fields are generated. When the number of flat-fields is limited to a small number, one may not achieve enough freedom. However, the partition of the mask into different masks (one for each pattern) is an easy way to allow for different motions at different places. Yet, this forces the direction of motion to be aligned with the one that is sampled. To provide even more freedom, the mask can be applied to each of the two components of a flat-field gradient.

In order to validate the above procedure, it is now applied to the actual test case. The edges of the projections are composed of only two spots: (2) and (5) (see Fig. 2). The residual in these regions will obviously be low, thus it cannot be considered as a criterion. The criterion to judge whether the procedure is accurate is the comparison of another spot not included in the minimization process: the fourth, which has the highest intensity value. The gradient field composed of the different extracted pattern gradients is shown Fig. 7. The result of the flat-field correction around the fourth spot is shown in Fig. 8 for the standard correction [proposed by Jailin et al. (2017)] (Fig. 8a) (i.e. with a library composed of two flat-fields) and the presently proposed correction with the gradient field (Fig. 8b).

Figure 8
Residual field of the flat-field correction of the intermediate flat-field $[F_{\rm m}]$ , (a) with the stationary method and (b) with the additional gradient field. Note that the same color bar is used in both images.

Because the standard procedure is weighted by the global intensity variations, it does not totally take into account the moving patterns. With the actual procedure, the displacement amplitude is obtained with the two moving spots on the edges and then used to correct the central area (and hence the other spots including the above-shown fourth one). The low and uniform residual with the proposed method shows that the separated pattern shape is well estimated. The obtained displacement amplitude $[|{\bf{u}}|]$ is 0.42 pixels and, because the previous analysis has been carried out, one can estimate the motion of spot 1 to about 1.15 µm.

Although the initial motivation was based on the reduction of artefacts in reconstructed images, the present procedure, albeit successful, turns out not to display an appreciable benefit. The main reason for this comes from the fact that a spurious static feature in the radiographs (if away from the projection of the rotation axis) is smeared over half or a complete ring in the reconstructed volume and becomes difficult to distinguish. One may note that this tolerance to local bias is one of the reasons for the tremendous success of tomographic reconstruction. Even when algebraic reconstruction techniques are used, the improvement due to erasing of the mobile patterns is difficult to evidence. Although the difference in the reconstructed fields with or without mobile pattern removal clearly shows the impact of the preprocessing, finding objective measurements to decide on the best reconstruction is subtle. Attempts made with entropy evaluation showed at most a 1% difference, a level than can hardly be considered as meaningful.

However, it has recently been proposed to track the time evolution of a deformable 3D object from a prior full tomographic reconstruction and the inspection of only a few of its projections at later times (Leclerc et al., 2015; Taillandier-Thomas et al., 2016a ,b ). Such algorithms are much less resilient to inaccurate or unfaithful projections, as just a few of them are used to extract all the required information about sample motion. For such demanding applications, the proposed treatment is expected to be much more beneficial.

6. Conclusion

The separation of a slowly varying background and a localized pattern is addressed in the difficult case of subpixel motion. Although the general problem is ill-posed, using assumptions on the compact support of one image and a smoothly varying background, a methodology is proposed to achieve the partition and to estimate the motion in orientation and magnitude.

The proposed algorithm was tested on a set of flat-fields acquired on a synchrotron beamline. Five different features were analyzed independently, and revealed subpixel translations. Based on this result, a very simple and generic treatment is proposed to correct for such artefacts prior to reconstruction. The proposed methodology reduces to the construction of a masked image difference as an enrichment of a library built to account for inhomogeneous beam intensity.

Funding information

This work was made possible by an ESRF grant for the experiment MA-501 on beamline ID19.

References

Auvray, V., Bouthemy, P. & Lienard, J. (2006). Proceedings of the 2006 International Conference on Image Processing, pp. 1057–1060. IEEE. Google Scholar
Be'ery, E. & Yeredor, A. (2008). Trans. Image Process. 17, 340–353. Google Scholar
Bergen, J. R., Burt, P. J., Hingorani, R. & Peleg, S. (1990). Proceedings of the Third International Conference on Computer Vision, pp. 27–32. IEEE. Google Scholar
Flot, D., Mairs, T., Giraud, T., Guijarro, M., Lesourd, M., Rey, V., van Brussel, D., Morawe, C., Borel, C., Hignette, O., Chavanne, J., Nurizzo, D., McSweeney, S. & Mitchell, E. (2010). J. Synchrotron Rad. 17, 107–118. Web of Science CrossRef CAS IUCr Journals Google Scholar
Gai, K., Shi, Z. & Zhang, C. (2008). Proceedigns of the 2008 IEEE Conference on Computer Vision and Pattern Recognition, pp. 1–8. Google Scholar
Gai, K., Shi, Z. & Zhang, C. (2012). Trans. Pattern Anal. Mach. Intell. 34, 19–32. Google Scholar
Irani, M. & Peleg, S. (1993). J. Vis. Commun. Image Representation, 4, 324–335. CrossRef Google Scholar
Jailin, C., Buffière, J.-Y., Hild, F., Poncelet, M. & Roux, S. (2017). J. Synchrotron Rad. 24, 220–231. Web of Science CrossRef IUCr Journals Google Scholar
Kameda, Y. & Minoh, M. (1996). Proceedings of the 1996 International Conference on Virtual Systems and Multimedia, pp. 135–140. Google Scholar
Leclerc, H., Roux, S. & Hild, F. (2015). Exp. Mech. 55, 275–287. CrossRef Google Scholar
Lipton, A. J., Fujiyoshi, H. & Patil, R. S. (1998). Workshop on Applications of Computer Vision, pp. 8–14. IEEE. Google Scholar
Marcia, R. F., Kim, C., Eldeniz, C., Kim, J., Brady, D. J. & Willett, R. M. (2008a). Opt. Express, 16, 16352–16363. CrossRef PubMed Google Scholar
Marcia, R. F., Kim, C., Kim, J., Brady, D. J. & Willett, R. M. (2008b). Proceedings of the IEEE International Conference on Image Processing, San Diego, CA, USA, pp. 2620–2623. IEEE. Google Scholar
Ng, K. K. & Delp, E. J. (2010). Proc. SPIE, 7543, 75430M. CrossRef Google Scholar
Ramirez-Manzanares, A., Palafox-Gonzalez, A. & Rivera, M. (2010). In Advances in Artificial Intelligence, MICAI 2010, Lecture Notes in Computer Science, Vol. 6437. edited by G. Sidorov, A. Hernández Aguirre & C. A. Reyes García Berlin, Heidelberg: Springer. Google Scholar
Ramírez-Manzanares, A., Rivera, M., Kornprobst, P. & Lauze, F. (2007). International Conference on Scale Space and Variational Methods in Computer Vision, pp. 227–238. Berlin: Springer. Google Scholar
Ramírez-Manzanares, A., Rivera, M., Kornprobst, P. & Lauze, F. (2011). J. Math. Imaging Vis. 40, 285–304. Google Scholar
Rogers, S. S., Waigh, T. A., Zhao, X. & Lu, J. R. (2007). Phys. Biol. 4, 220–227. CrossRef PubMed CAS Google Scholar
Schreier, H., Orteu, J. & Sutton, M. (2009). Image Correlation for Shape, Motion and Deformation Measurements. New York: Springer. Google Scholar
Stuke, I., Aach, T., Barth, E. & Mota, C. (2004). Intl J. Comput. Inf. Sci. 5, 141–152. Google Scholar
Szeliski, R., Avidan, S. & Anandan, P. (2000). Proceedings of the 2000 IEEE Conference on Computer Vision and Pattern Recognition, Vol. 1, pp. 246–253. IEEE. Google Scholar
Taillandier-Thomas, T., Jailin, C., Roux, S. & Hild, F. (2016a). Proc. SPIE, 9896, 98960L. Google Scholar
Taillandier-Thomas, T., Roux, S. & Hild, F. (2016b). Phys. Rev. Lett. 117, 025501. PubMed Google Scholar
Toro, J., Owens, F. & Medina, R. (2003). Pattern Recognit. Lett. 24, 597–605. CrossRef Google Scholar
Weitkamp, T., Haas, D., Wegrzynek, D. & Rack, A. (2011). J. Synchrotron Rad. 18, 617–629. Web of Science CrossRef CAS IUCr Journals Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 25| Part 1| January 2018| Pages 272-281

https://doi.org/10.1107/S1600577517015892