In tomographic imaging, the aim is to characterize the three-dimensional structure of an object from X-ray projections. In applications like medical CT, projections are gathered from all directions. This allows for a relatively straight-forward reconstruction of the object using so-called filtered back-projection methods (FBP)[10, 15, 26]. When only a limited angular sampling is available, more advanced iterative reconstruction techniques that use prior information about the structure of the object have been developed [11, 18, 20, 3, 12]. In this paper, we consider an extreme case where we acquire only a single X-ray projection on which to base a full three-dimensional reconstruction. We refer to this problem as single-shot X-ray tomography. This problem is highly relevant for many applications, including industrial quality control and high-throughput imaging [13, 17, 16, 5]. A typical setup consists of a fixed X-ray source and detector that collects single X-ray projections of the objects of interest. While we specifically target reconstruction from single-angle X-ray projections, the techniques we develop are also relevant for limited-angle tomography with applications including high-resolution dynamic imaging[21, 25] and cryo-electron microscopy [23, 19].
|(a) True||(b) Projection||(c) FBP|
|(d) TV||(e) SSC||(f) CoShaRP|
The single-shot X-ray tomography problem is extremely under-determined, making it an ill-posed inverse problem. A single cone-beam projection of a volume containing voxels consists of measurements. Hence, the measurements are undersampled by a factor of . To reduce the ill-posedness, it is a general practice to incorporate prior information via regularization. However, classic regularization methods fail on single-shot X-ray tomography, as demonstrated in Fig. 1.
Since it is evident that a strong prior is necessary to recover the target image from a single projection, we consider the class of objects that are composed of a limited number of known building blocks. This is a reasonable assumption when imaging materials that are made up of basic structures, for example, a 3D structure comprising of interlocking bricks or a protein structure consisting of repeated amino acid groups. Hence, if such shapes are known a-priori along with their number of repetitions, the image estimation problem can be recast as an estimation of roto-translation parameters of these shapes. We term the process of estimating the shape parameters (i.e., their roto-translation in the image) from the linear measurements of the image that constitutes them as shape sensing. Unlike the image estimation problem – which has a simple linear structure – the estimation of roto-translation parameters (i.e., shape sensing) is a non-linear problem. This, in turn, makes the inversion process computationally intractable. To avoid such non-linearity, we use a shape dictionary approach that expresses the target image as a linear combination of shapes from the available dictionary. Due to the linear structure, we show that it is possible to recover the shapes from a single-shot X-ray projection (see Fig. 1) by solving a convex problem. We regard the recovery of shapes from their tomographic projections as tomographic shape sensing.
1.1 Related Work
The shape sensing problem has been studied in the context of shape-based characterization, object tracking and optical character recognition. Inspired by the compressed sensing, a recent approach called Sparse Shape Composition (SSC) imposes an -norm constraint on the shape-dictionary coefficients [1, 2]. SSC has the advantage that it can form new shapes from the intersection and union of basic shapes. However, the main drawback comes from large feasible solution space inherent to the -norm constraint in high dimension (see Remark 3.2). This large feasible space may lead to an incorrect solution. A simplistic version of SSC performs the 3D characterization of nano-particles using electron tomography . This method uses a simple -norm constraint to recover spherical nano-particles from their tomographic projections. However, their tomographic projections have a parallel-beam geometry and require measurements from more than one projection angle. Although SSC has been extended to tomographic inverse problem, it has never been tested for single-shot tomographic shape sensing. However, we demonstrate the failure of SSC in single-shot cone-beam tomography in Fig. 1.
1.2 Contributions and Outline
To the best of our knowledge, the single-shot X-ray tomography has never been studied, and no reconstruction method exist till date to recover back an image successfully from a single-shot. We introduce the tomographic shape sensing problem that assumes the prior knowledge about the shapes in the image. The principal contribution of this paper is the development of the convex program CoShaRP to reconstruct images composed of non-overlapping shapes from a single-shot. The convex program is novel in the sense that the simplex-type constraint enables sharp recovery results from extremely under-determined single-shot tomographic projections. Although the exact recovery problem is NP-hard, our proposed convex program CoShaRP stably recovers the shapes. Moreover, we propose a primal-dual algorithm to find the optimal solution of CoShaRP. The algorithm does not rely on the inversion of large matrices and uses a simple proximal operator to project onto a simplex constraint. Using numerical experiments, we answer the following questions:
What is the minimum resolution of the shape that CoShaRP can sense?
Is CoShaRP robust to the number of repetitions and the possible rotations of the shapes in the target image?
Can CoShaRP recover non-homogeneous as well as non-convex shapes?
How sensitive is CoShaRP to the measurement noise?
We discuss the single-shot tomographic inverse problem in Section 2. Section 3 discusses the tomographic shape sensing problem and introduces a convex program CoShaRP. We describe an efficient iterative scheme to find an optimal solution to CoShaRP in Section 4. We illustrate the numerical experiments in Section 5 and conclude the paper in Section 6.
Throughout this paper, boldface letters (
) denote vectors in. The identity and zero elements are denoted by and respectively. The Euclidean inner product is denoted by for with a corresponding norm . However, for all other norms, we use subscripts (e.g., , ). To represent the matrices, we use uppercase letters (e.g., ). The elements of a matrix are denoted by . All the functions are represented as , where and are the domain and co-domain of , respectively. We denote the convex conjugate of the function by . denotes the proximal of function evaluated at point (for definition, please refer to ). We represent the optimal solution to the optimization problem using overline (e.g., ).
2 Single-shot X-ray Tomography
The acquisition geometry of single-shot X-ray tomography consists of one source and an array of regularly spaced detectors. Let be a directional vector, and be any position vector, where is the dimension of the scene. The cone-beam transform of an image function is its integral along a line in the direction passing through . It is mathematically given by
In a single-shot setup, we have a source located at . It sends multiple X-rays through the object (compactly supported on ) in a cone with a vertex at and spanning angles in a set that determines the geometry of cone. Let these angles be , , then the measurement is
where denotes the value of in the voxel and is the contribution of the voxel to the ray. The measurements can now be expressed as a linear system of equations
The above linear system of equations is extremely under-determined since the number of measurements is far smaller than the number of unknowns . We do assume that each voxel is intersected by at least one ray, so that each column of the matrix has at least one non-zero element. Determining the image from the measurements is an ill-posed inverse problem. In general, to resolve this ill-posedness, regularization needs to be added in the inversion procedure to incorporate the prior knowledge about the target image. However, conventional regularization techniques are not sufficient to resolve the true image, as was illustrated in Fig. 1.
3 Convex Shape Recovery
In this section, we discuss the tomographic shape sensing problem. Our formulation hinges on the formation of a dictionary that consists of possible roto-translations of the known shapes and the representation of target image as a convex combination of dictionary elements.
3.1 Image model and Dictionary
Let the functions , , denote the compactly supported shape functions. The image is now assumed to be composed of roto-translations of these shapes
where and are the angle and the shift of copy of shape respectively, and is a rotation matrix that depends on the angle . The total number of shapes in an image are . Hence, from the knowledge of the shapes, the image estimation translates to finding the roto-translation parameters of the shapes. However, the image is a non-linear function of these parameters. Hence, the recovery of these parameters becomes a computationally intractable problem due to the non-convex structure of the cost function.
To mitigate the non-linearity associated with the roto-translation parameters, we create a shape dictionary that consists of roto-translations of the shapes. Let the dictionary
where covers possible roto-translations of the shapes. Hence, the target image can be represented as a linear combination of the elements of this dictionary,
where is a coefficient vector, with . Hence, the shape recovery problem is to find a high-dimensional binary vector from its linear measurements
Here, contains the projections of the individual dictionary elements, sampled at the appropriate points. The binary constraints on make the recovery problem an integer program, and hence, NP-hard in general .
3.2 Convex Shape Recovery Program (CoShaRP)
The binary constraints on the coefficients can be relaxed using the bounds constraints. Moreover, a Gaussian assumption on the noise leads to a least-squares formulation for the data misfit. Hence, the resulting convex program, which we refer to as the Convex Shape Recovery Program (CoShaRP), reads
Here, the inequality between vectors is imposed elementwise. Note that we have used the Euclidean norm instead of its square to measure the misfit.
The geometric interpretation of CoShaRP is as follows: We are trying to find a high-dimensional vector
closest to the hyperplanein a Euclidean sense that lies on the intersection of a hyperplane and the hyperbox . In Fig. 2, we show the geometry for a shape-sensing problem with two possible shapes (). In Fig. 2(a), the hyperplane corresponding to tomographic measurements intersects the hyperplane corresponds to equality constraints () outside the hyperbox. Hence, the solution to CoShaRP in this case is binary. However, a binary solution can not always be guaranteed as these hyperplanes may intersect inside the hyperbox (cf. Fig. 2(b)). In such cases, further post-processing is required to retrieve the target image. For more, refer to Section 4.3.
The CoShaRP consists of constraints that are defined by K-simplex. The K-simplex, defined as
is a generalized version simplex (simplex has ). K-simplex represents a polytope in -dimension with vertices. Moreover, these polytopes are regular, i.e., they posses highest level of symmetry. We plot -simplex in three dimension in Fig. 3. These simplices are equilateral triangle except for . However, it is important to note that the number of shapes in the target image will be much smaller than the number of dictionary elements (i.e., ). Hence, we will frequently encounter feasible regions to be extremely low-dimensional polytope embedded in a high-dimensional space.
The -norm ball is bigger in size than the K-simplex constraint set used in CoShaRP. In particular, the K-simplex constraint represents the strict ball, i.e., , in the non-negative region. Since the solution lies on the corners of the K-simplex constraint, the recovery with CoShaRP is stronger than that of SSC (see, e.g., Fig. 1).
|(a) 1-simplex||(b) 2-simplex||(c) 3-simplex|
We discuss a fast iterative scheme to find an approximate solution of CoShaRP numerically. The iterative scheme is based on splitting strategy that separates the non-smooth part from the smooth. We also introduce a thresholding method to recover the image from the coefficient vector, in case the solution is not binary.
4.1 Primal-Dual Algorithm
For simplicity, we express CoShaRP in the following form:
for , where , with , are parameters that controls the speed of convergence. The main characteristic of this algorithm is that we avoid an inversion of a large matrix which often occurs in other splitting methods such as alternating direction method of multipliers . Moreover, the proximal of both functions are easy to compute.
The proposed primal-dual algorithm does not require user to store a large dictionary matrix as well as tomography matrix . Hence, the algorithm can utilize the functional forms of the dictionary as well as tomography operator since it only requires the forward and the adjoint operation with the operator.
4.2 Proximal operators
The conjugate function of is
We refer to [6, Example 3.26] for the derivation of the conjugate function. The conjugate function is linear inside the Euclidean norm ball of size 1 and outside. Hence, the conjugate function is convex. Its proximal operator is given in the following theorem.
The proximal operator of function
where is a known vector, is given by
The proximal operator for function reads
The Euclidean norm constraints enforces two cases: (i) The optimal point without the constraints is . This optimal solution holds when . (ii) When , the optimal solution lies on the surface of the Euclidean norm ball with size 1. Moreover, the optimal solution is in the direction of . Hence, the proximal point is
This concludes the proof. ∎
We use the following theorem to compute the proximal of , adapted from [4, Theorem 6.27].
Theorem 2 (projection onto the intersection of a hyperplane and a box).
Let be a set. The proximal operator of an indicator function to the set, , is given by
where is a projection onto the box and is a solution of the equation .
The orthogonal projection of is the unique solution of
A Lagrangian of this minimization problem reads
where is a Lagrange multiplier. It follows from the strong duality that is an optimal solution of problem (5) if and only if there exists a dual variable for which
The projection onto the box is simple. It is done component-wise as . However, equation (4) consists of finding a root of the non-increasing function Since is a non-increasing function for any i, is a non-increasing function. Its root can be found using the Newton procedure, where derivative is
4.3 Image Formation
The convex program CoShaRP does not always lead to a binary solution (refer to Fig. 2). Moreover, if the optimization procedure is terminated early, we may not have a binary solution. Hence, an accurate image formation process is essential to retrieve the target image from the non-binary shape coefficient vector resulting from CoShaRP. We propose the image formation procedure based on sorting of the coefficients. We first sort the coefficients in descending order, and selectively form the image consistent with the measurements. Algorithm 2 enumerates the steps in the image formation process. Here, is a natural basis vector with non-zero element located at position. To make sure the shapes do not overlap, we also add necessary conditions (see step 4 in Algorithm 2).
5 Numerical Experiments
In this section, we try to answer questions regarding resolution, sparsity, rotations and the performance under noise using 2D numerical experiments. For all the experiments, we have images of size discretized on pixels, and the tomography matrix has at least measurements. The typical tomography setup is shown in Fig. 4. In the CoShaRP performance plots, we generate different realizations of the target image with given constraints (for examples, size, rotations and repetitions of the shape), and the success rate is measured from the average over all instances. We say an instance is successful if the recovered image is close to the target image (in Euclidean norm).
5.1 Resolution Analysis
In this experiment, we estimate the required minimum width of the shape sensed by a single-shot. For simplicity, we consider circular disc of constant intensity with size varying from to pixels. Fig. 5(a) shows the performance of CoShaRP against varying sizes of the disc. As the number of measurements (i.e., detector pixels) is increased, the success rate increases implying that the recovery of even single-pixel shapes is possible with CoShaRP.
5.2 Invariance with respect to Density and Rotation
We first look at the success of CoShaRP with multiple repetitions of the shape. We consider a circular disc with size pixels and the number of repetitions in the image from to . The top figure in Fig. 5(b) shows the performance of CoShaRP with density. This shows that the CoShaRP is insensitive to the number of repetitions of the shapes in the image. Next, we take an ellipsoidal disc with semi-axes and . This image is rotated for angles making sure that each angle represents a different shape on the pixels. The bottom figure in Fig. 5(b) demonstrates the performance of CoShaRP with the number of possible rotations. This implies that the CoShaRP is insensitive to the number of possible rotations of the shapes.
|(a) Resolution||(b) Invariance||(c) Shapes||(d) Noise|
5.3 Non-homogeneous and Non-convex Shapes
The top figure in Fig. 5(c) provides the reconstruction from CoShaRP for non-homogeneous shape. We consider a circular shape with four different intensities varying radially. The true image consists of 5 repetitions of this shape. The CoShaRP recovers these 5 copies successfully as shown from the shape coefficients (given below the figure). For non-convex shape, we consider ellipsoidal shell with outer axes and , and inner axes and . The bottom figure in Fig. 5(c) gives the reconstruction with CoShaRP for a true image with 6 repetitions of the above-mentioned non-convex shape. From these two experiments, we conclude that the CoShaRP can recover the non-homogeneous and non-convex shapes.
5.4 Measurement Noise
We consider the true image shown in Fig. 1 and consider three noisy scenarios where the Gaussian noise of strength , and is added to 1024 measurements. In Fig. 5(d), we plot the measurements (on the left) and the reconstructed image (on the right) for the above-mentioned noise values. In the measurements plots, the true noisy data is mentioned as ‘true’, while the forward projected data from the reconstructed image is denoted by ‘rec’. The plots below them show the difference between the two. From Fig. 5(d), it is evident that the CoShaRP is stable till noise, while fails for extremely noisy measurements.
We introduced a single-shot tomographic shape sensing problem that aims to recover shapes from a single cone-beam projection. To solve this problem, we develop a convex program CoShaRP. CoShaRP is novel in the sense that the simplex-type constraint enables sharp recovery results from extremely under-determined single-shot tomographic projections. Moreover, we propose a primal-dual algorithm to find an approximately optimal solution to CoShaRP quickly. The numerical results demonstrate that (i) the resolution limit to sense the shape depends on the number of measurements, (ii) CoShaRP is insensitive to the number of repetitions of the shape and the number of possible rotations of the shape, (iii) CoShaRP can sense the non-homogeneous as well as non-convex shapes, (iv) CoShaRP tolerates only a moderate amount of measurement noise.
The limitations of CoShaRP are as follows: (i) The roto-translations of the shapes must be included in the dictionary for the exact recovery of the target image. This inclusion requirement makes CoShaRP a computationally expensive approach due to the large dictionary size. (ii) CoShaRP also requires the correct knowledge of shapes and their intensity. If the shape is not known accurately, the CoShaRP may fail. (iii) CoShaRP relies on the knowledge of the number of shapes in the target image. CoShaRP relies on the knowledge of the number of shapes in the target image. Its estimation may be a costly procedure since one has to solve CoShaRP several times.
The authors thank Nick Luiken for stimulating discussions. This work was supported by the Dutch Research Council (grant no. OCENW.XS.039).
-  Alireza Aghasi and Justin Romberg. Sparse shape reconstruction. SIAM Journal on Imaging Sciences, 6(4):2075–2108, January 2013.
-  Alireza Aghasi and Justin Romberg. Convex cardinal shape composition. SIAM Journal on Imaging Sciences, 8(4):2887–2950, January 2015.
-  K. J. Batenburg and J. Sijbers. DART: A practical reconstruction algorithm for discrete tomography. IEEE Transactions on Image Processing, 20(9):2542–2553, September 2011.
-  Amir Beck. First-Order Methods in Optimization. Society for Industrial and Applied Mathematics, October 2017.
-  Tekin Bicer, Doğa Gürsoy, Vincent De Andrade, Rajkumar Kettimuthu, William Scullin, Francesco De Carlo, and Ian T. Foster. Trace: a high-throughput tomographic reconstruction engine for large-scale datasets. Advanced Structural and Chemical Imaging, 3(1), January 2017.
-  Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.
Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato, and Jonathan Eckstein.
Distributed optimization and statistical learning via the alternating
direction method of multipliers.
Foundations and Trends® in Machine Learning, 3(1):1–122, 2010.
-  Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120–145, December 2010.
-  Patrick L. Combettes and Jean-Christophe Pesquet. Proximal splitting methods in signal processing. In Springer Optimization and Its Applications, pages 185–212. Springer New York, 2011.
-  M. Defrise and R. Clack. A cone-beam reconstruction algorithm using shift-variant filtering and cone-beam backprojection. IEEE Transactions on Medical Imaging, 13(1):186–195, March 1994.
-  A.H. Delaney and Y. Bresler. Globally convergent edge-preserving regularized reconstruction: an application to limited-angle tomography. IEEE Transactions on Image Processing, 7(2):204–221, 1998.
-  Jürgen Frikel. Sparse regularization in limited angle tomography. Applied and Computational Harmonic Analysis, 34(1):117–141, January 2013.
-  Max F. Hantke, Dirk Hasse, Filipe R. N. C. Maia, Tomas Ekeberg, Katja John, Martin Svenda, N. Duane Loh, Andrew V. Martin, Nicusor Timneanu, Daniel S. D. Larsson, Gijs van der Schot, Gunilla H. Carlsson, Margareta Ingelman, Jakob Andreasson, Daniel Westphal, Mengning Liang, Francesco Stellato, Daniel P. DePonte, Robert Hartmann, Nils Kimmel, Richard A. Kirian, M. Marvin Seibert, Kerstin Mühlig, Sebastian Schorb, Ken Ferguson, Christoph Bostedt, Sebastian Carron, John D. Bozek, Daniel Rolles, Artem Rudenko, Sascha Epp, Henry N. Chapman, Anton Barty, Janos Hajdu, and Inger Andersson. High-throughput imaging of heterogeneous cell organelles with an x-ray laser. Nature Photonics, 8(12):943–949, November 2014.
-  Richard M. Karp. Reducibility among combinatorial problems. In Complexity of Computer Computations, pages 85–103. Springer US, 1972.
-  Hiroyuki Kudo, Frédéric Noo, and Michel Defrise. Cone-beam filtered-backprojection algorithm for truncated helical data. Physics in Medicine and Biology, 43(10):2885–2909, October 1998.
-  Nouamane Laanait, Wittawat Saenrang, Hua Zhou, Chang-Beom Eom, and Zhan Zhang. Dynamic x-ray diffraction imaging of the ferroelectric response in bismuth ferrite. Advanced Structural and Chemical Imaging, 3(1), March 2017.
-  Daniël M. Pelt and Vincent De Andrade. Improved tomographic reconstruction of large-scale real-world data by filter optimization. Advanced Structural and Chemical Imaging, 2(1), December 2016.
-  M Persson, D Bone, and H Elmqvist. Total variation norm for three-dimensional iterative reconstruction in limited view angle tomography. Physics in Medicine and Biology, 46(3):853–866, February 2001.
-  Sjors H.W. Scheres. RELION: Implementation of a bayesian approach to cryo-EM structure determination. Journal of Structural Biology, 180(3):519–530, December 2012.
-  Emil Y Sidky and Xiaochuan Pan. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization. Physics in Medicine and Biology, 53(17):4777–4807, August 2008.
-  J C H Spence, U Weierstall, and H N Chapman. X-ray lasers for structural and dynamic biology. Reports on Progress in Physics, 75(10):102601, September 2012.
-  Daniele Zanaga, Folkert Bleichrodt, Thomas Altantzis, Naomi Winckelmans, Willem Jan Palenstijn, Jan Sijbers, Bart de Nijs, Marijn A. van Huis, Ana Sánchez-Iglesias, Luis M. Liz-Marzán, Alfons van Blaaderen, K. Joost Batenburg, Sara Bals, and Gustaaf Van Tendeloo. Quantitative 3d analysis of huge nanoparticle assemblies. Nanoscale, 8(1):292–299, 2016.
-  X. Zhang, E. Settembre, C. Xu, P. R. Dormitzer, R. Bellamy, S. C. Harrison, and N. Grigorieff. Near-atomic resolution using electron cryomicroscopy and single-particle reconstruction. Proceedings of the National Academy of Sciences, 105(6):1867–1872, January 2008.
-  Xiaoqun Zhang, Martin Burger, and Stanley Osher. A unified primal-dual algorithm framework based on bregman iteration. Journal of Scientific Computing, 46(1):20–46, August 2010.
-  Wenjuan Zhu, Wenbo Ma, Yirong Su, Zeng Chen, Xinya Chen, Yaoguang Ma, Lizhong Bai, Wenge Xiao, Tianyu Liu, Haiming Zhu, Xiaofeng Liu, Huafeng Liu, Xu Liu, and Yang (Michael) Yang. Low-dose real-time x-ray imaging with nontoxic double perovskite scintillators. Light: Science & Applications, 9(1), June 2020.
-  Yu Zou and Xiaochuan Pan. Image reconstruction on PI-lines by use of filtered backprojection in helical cone-beam CT. Physics in Medicine and Biology, 49(12):2717–2731, June 2004.