1 Introduction
In tomographic imaging, the aim is to characterize the threedimensional structure of an object from Xray projections. In applications like medical CT, projections are gathered from all directions. This allows for a relatively straightforward reconstruction of the object using socalled filtered backprojection 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 Xray projection on which to base a full threedimensional reconstruction. We refer to this problem as singleshot Xray tomography. This problem is highly relevant for many applications, including industrial quality control and highthroughput imaging [13, 17, 16, 5]. A typical setup consists of a fixed Xray source and detector that collects single Xray projections of the objects of interest. While we specifically target reconstruction from singleangle Xray projections, the techniques we develop are also relevant for limitedangle tomography with applications including highresolution dynamic imaging[21, 25] and cryoelectron microscopy [23, 19].
(a) True  (b) Projection  (c) FBP 
(d) TV  (e) SSC  (f) CoShaRP 
The singleshot Xray tomography problem is extremely underdetermined, making it an illposed inverse problem. A single conebeam projection of a volume containing voxels consists of measurements. Hence, the measurements are undersampled by a factor of . To reduce the illposedness, it is a general practice to incorporate prior information via regularization. However, classic regularization methods fail on singleshot Xray 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 apriori along with their number of repetitions, the image estimation problem can be recast as an estimation of rototranslation parameters of these shapes. We term the process of estimating the shape parameters (i.e., their rototranslation 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 rototranslation parameters (i.e., shape sensing) is a nonlinear problem. This, in turn, makes the inversion process computationally intractable. To avoid such nonlinearity, 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 singleshot Xray 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 shapebased 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 shapedictionary 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 nanoparticles using electron tomography [22]. This method uses a simple norm constraint to recover spherical nanoparticles from their tomographic projections. However, their tomographic projections have a parallelbeam 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 singleshot tomographic shape sensing. However, we demonstrate the failure of SSC in singleshot conebeam tomography in Fig. 1.
1.2 Contributions and Outline
To the best of our knowledge, the singleshot Xray tomography has never been studied, and no reconstruction method exist till date to recover back an image successfully from a singleshot. 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 nonoverlapping shapes from a singleshot. The convex program is novel in the sense that the simplextype constraint enables sharp recovery results from extremely underdetermined singleshot tomographic projections. Although the exact recovery problem is NPhard, our proposed convex program CoShaRP stably recovers the shapes. Moreover, we propose a primaldual 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:

[topsep=0pt,itemsep=1ex,partopsep=1ex,parsep=1ex]

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 nonhomogeneous as well as nonconvex shapes?

How sensitive is CoShaRP to the measurement noise?
We discuss the singleshot 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.
1.3 Notation
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 codomain of , respectively. We denote the convex conjugate of the function by . denotes the proximal of function evaluated at point (for definition, please refer to [9]). We represent the optimal solution to the optimization problem using overline (e.g., ).2 Singleshot Xray Tomography
The acquisition geometry of singleshot Xray 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 conebeam transform of an image function is its integral along a line in the direction passing through . It is mathematically given by
In a singleshot setup, we have a source located at . It sends multiple Xrays 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 underdetermined 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 nonzero element. Determining the image from the measurements is an illposed inverse problem. In general, to resolve this illposedness, 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 rototranslations 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 rototranslations 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 rototranslation parameters of the shapes. However, the image is a nonlinear function of these parameters. Hence, the recovery of these parameters becomes a computationally intractable problem due to the nonconvex structure of the cost function.
To mitigate the nonlinearity associated with the rototranslation parameters, we create a shape dictionary that consists of rototranslations of the shapes. Let the dictionary
where covers possible rototranslations of the shapes. Hence, the target image can be represented as a linear combination of the elements of this dictionary,
with 
where is a coefficient vector, with . Hence, the shape recovery problem is to find a highdimensional 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, NPhard in general [14].
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 leastsquares formulation for the data misfit. Hence, the resulting convex program, which we refer to as the Convex Shape Recovery Program (CoShaRP), reads
(1) 
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 highdimensional vector
closest to the hyperplane
in a Euclidean sense that lies on the intersection of a hyperplane and the hyperbox . In Fig. 2, we show the geometry for a shapesensing 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 postprocessing is required to retrieve the target image. For more, refer to Section 4.3.The CoShaRP consists of constraints that are defined by Ksimplex. The Ksimplex, defined as
is a generalized version simplex (simplex has ). Ksimplex 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 lowdimensional polytope embedded in a highdimensional space.
(a)  (b) 



Remark.
CoShaRP differs significantly from the Sparse Shape Composition (SSC) [1, 22]. SSC formulates the shapesensing problem as
(2) 
The norm ball is bigger in size than the Ksimplex constraint set used in CoShaRP. In particular, the Ksimplex constraint represents the strict ball, i.e., , in the nonnegative region. Since the solution lies on the corners of the Ksimplex constraint, the recovery with CoShaRP is stronger than that of SSC (see, e.g., Fig. 1).
(a) 1simplex  (b) 2simplex  (c) 3simplex 



4 Optimization
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 nonsmooth 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 PrimalDual Algorithm
For simplicity, we express CoShaRP in the following form:
minimize  
and is the indicator function of the set . To solve this optimization problem, we use a primaldual splitting algorithm [24, 8]. The iterates of this primaldual algorithm takes 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 [7]. Moreover, the proximal of both functions are easy to compute.
Remark.
The proposed primaldual 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.
Theorem 1.
The proximal operator of function
where is a known vector, is given by
(3) 
Proof.
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
(4) 
where is a projection onto the box and is a solution of the equation .
Proof.
The orthogonal projection of is the unique solution of
(5) 
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
(6)  
(7) 
Using the expression for the Lagrangian, the relation (6) can be equivalently written as . The feasibility condition (7) then takes a form . ∎
Remark.
The projection onto the box is simple. It is done componentwise as . However, equation (4) consists of finding a root of the nonincreasing function Since is a nonincreasing function for any i, is a nonincreasing 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 nonbinary 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 nonzero 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).
For all the experiments, we run Algorithm 1 with , with . Moreover, we set and . Once we obtain the vector , we run Algorithm 2 to form the image.
5.1 Resolution Analysis
In this experiment, we estimate the required minimum width of the shape sensed by a singleshot. 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 singlepixel 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 semiaxes 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 

shapes 
Density

Nonhomogeneous

SNR 
SNR  
Rotation

Nonconvex


SNR  
5.3 Nonhomogeneous and Nonconvex Shapes
The top figure in Fig. 5(c) provides the reconstruction from CoShaRP for nonhomogeneous 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 nonconvex 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 abovementioned nonconvex shape. From these two experiments, we conclude that the CoShaRP can recover the nonhomogeneous and nonconvex 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 abovementioned 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.
6 Conclusions
We introduced a singleshot tomographic shape sensing problem that aims to recover shapes from a single conebeam projection. To solve this problem, we develop a convex program CoShaRP. CoShaRP is novel in the sense that the simplextype constraint enables sharp recovery results from extremely underdetermined singleshot tomographic projections. Moreover, we propose a primaldual 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 nonhomogeneous as well as nonconvex shapes, (iv) CoShaRP tolerates only a moderate amount of measurement noise.
The limitations of CoShaRP are as follows: (i) The rototranslations 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.
Acknowledgments
The authors thank Nick Luiken for stimulating discussions. This work was supported by the Dutch Research Council (grant no. OCENW.XS.039).
References
 [1] Alireza Aghasi and Justin Romberg. Sparse shape reconstruction. SIAM Journal on Imaging Sciences, 6(4):2075–2108, January 2013.
 [2] Alireza Aghasi and Justin Romberg. Convex cardinal shape composition. SIAM Journal on Imaging Sciences, 8(4):2887–2950, January 2015.
 [3] 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.
 [4] Amir Beck. FirstOrder Methods in Optimization. Society for Industrial and Applied Mathematics, October 2017.
 [5] Tekin Bicer, Doğa Gürsoy, Vincent De Andrade, Rajkumar Kettimuthu, William Scullin, Francesco De Carlo, and Ian T. Foster. Trace: a highthroughput tomographic reconstruction engine for largescale datasets. Advanced Structural and Chemical Imaging, 3(1), January 2017.
 [6] Stephen Boyd, Stephen P Boyd, and Lieven Vandenberghe. Convex optimization. Cambridge university press, 2004.

[7]
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.  [8] Antonin Chambolle and Thomas Pock. A firstorder primaldual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120–145, December 2010.
 [9] Patrick L. Combettes and JeanChristophe Pesquet. Proximal splitting methods in signal processing. In Springer Optimization and Its Applications, pages 185–212. Springer New York, 2011.
 [10] M. Defrise and R. Clack. A conebeam reconstruction algorithm using shiftvariant filtering and conebeam backprojection. IEEE Transactions on Medical Imaging, 13(1):186–195, March 1994.
 [11] A.H. Delaney and Y. Bresler. Globally convergent edgepreserving regularized reconstruction: an application to limitedangle tomography. IEEE Transactions on Image Processing, 7(2):204–221, 1998.
 [12] Jürgen Frikel. Sparse regularization in limited angle tomography. Applied and Computational Harmonic Analysis, 34(1):117–141, January 2013.
 [13] 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. Highthroughput imaging of heterogeneous cell organelles with an xray laser. Nature Photonics, 8(12):943–949, November 2014.
 [14] Richard M. Karp. Reducibility among combinatorial problems. In Complexity of Computer Computations, pages 85–103. Springer US, 1972.
 [15] Hiroyuki Kudo, Frédéric Noo, and Michel Defrise. Conebeam filteredbackprojection algorithm for truncated helical data. Physics in Medicine and Biology, 43(10):2885–2909, October 1998.
 [16] Nouamane Laanait, Wittawat Saenrang, Hua Zhou, ChangBeom Eom, and Zhan Zhang. Dynamic xray diffraction imaging of the ferroelectric response in bismuth ferrite. Advanced Structural and Chemical Imaging, 3(1), March 2017.
 [17] Daniël M. Pelt and Vincent De Andrade. Improved tomographic reconstruction of largescale realworld data by filter optimization. Advanced Structural and Chemical Imaging, 2(1), December 2016.
 [18] M Persson, D Bone, and H Elmqvist. Total variation norm for threedimensional iterative reconstruction in limited view angle tomography. Physics in Medicine and Biology, 46(3):853–866, February 2001.
 [19] Sjors H.W. Scheres. RELION: Implementation of a bayesian approach to cryoEM structure determination. Journal of Structural Biology, 180(3):519–530, December 2012.
 [20] Emil Y Sidky and Xiaochuan Pan. Image reconstruction in circular conebeam computed tomography by constrained, totalvariation minimization. Physics in Medicine and Biology, 53(17):4777–4807, August 2008.
 [21] J C H Spence, U Weierstall, and H N Chapman. Xray lasers for structural and dynamic biology. Reports on Progress in Physics, 75(10):102601, September 2012.
 [22] Daniele Zanaga, Folkert Bleichrodt, Thomas Altantzis, Naomi Winckelmans, Willem Jan Palenstijn, Jan Sijbers, Bart de Nijs, Marijn A. van Huis, Ana SánchezIglesias, Luis M. LizMarzá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.
 [23] X. Zhang, E. Settembre, C. Xu, P. R. Dormitzer, R. Bellamy, S. C. Harrison, and N. Grigorieff. Nearatomic resolution using electron cryomicroscopy and singleparticle reconstruction. Proceedings of the National Academy of Sciences, 105(6):1867–1872, January 2008.
 [24] Xiaoqun Zhang, Martin Burger, and Stanley Osher. A unified primaldual algorithm framework based on bregman iteration. Journal of Scientific Computing, 46(1):20–46, August 2010.
 [25] 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. Lowdose realtime xray imaging with nontoxic double perovskite scintillators. Light: Science & Applications, 9(1), June 2020.
 [26] Yu Zou and Xiaochuan Pan. Image reconstruction on PIlines by use of filtered backprojection in helical conebeam CT. Physics in Medicine and Biology, 49(12):2717–2731, June 2004.
Comments
There are no comments yet.