1 Introduction
In the past decades optical trapping has established itself as a versatile, and easytoimplement tool for holding and moving microscopic (msized) particles suspended in a liquid in a contactfree and noninvasive manner [JonMarVol15]. The optical forces needed to balance the weight scale with the particle volume. Thus, for increasing particle size heating of the sample due to absorption of the trapping light leads to a problem, except for very transparent specimens. To hold heavier samples (for example cellclusters, entire microorganisms, or socalled “organoids”, miniorgans grown in the petridish) other trapping modalities, such as acoustic forces, have to be added as auxiliary trapping forces [ThaSteMeiHilBer11].
In this paper, we study the problem of estimation the movement of trapped particles, as a preprocessing step for 3D tomographic imaging. Conveniently, the required tomographic projection images from different viewpoints can be taken by rotating the trapped specimen directly by means of the acoustic or optical trapping forces. Here we will assume to be in a setting for which the light propagation can be sufficiently well modelled by straight rays which only undergo spatially varying attenuation, as in Xray tomography [Nat01]. In this case the optical image resembles a projection image. This is, for instance, fulfilled in low numerical aperture imaging of biological samples with sufficient amplitude contrast.
In contrast to classical computerized tomography, however, the motion of the particle in the trap is irregular: Since the particle is not completely immobilized in the trap and since the locally acting forces are typically inhomogeneous during an entire revolution, it undergoes a complicated motion described by an affine transformation with timedependent  and unknown  parameters, the momentary translation and angular velocity vectors. The focus of this work is to provide, as a first step of the reconstruction, a recovery of the motion of the particle without knowledge of the inner structure. As soon as the motion is determined, the attenuation projection data can be properly aligned to give the threedimensional Xraytransform, see, for example,
[Hel99], (at least for those directions which were attained during the motion) and from this, a standard reconstruction formula could be applied to recover the inner structure of the object.A main ingredient for our reconstruction of the motion is the common line method known for the determination of the relative orientations of different projection images as it is used, for example, in cryogenic electron microscopy [Hee87, Gon88a]. An infinitesimal version of this technique, allows us (after correcting for a potential translation by tracking the center of the images, as explained in Section 3.1) to find at every time step the local angular velocity of the motion, see Section 3.2. Moreover, we provide a numerical test by presenting reconstructions of simulated data (in Section 4).
2 Mathematical model
We consider a bounded object in , characterised by an attenuation coefficient .
Definition
In this paper
denotes the space of admissible attenuation coefficients. For we define its center
(2.1) 
In this paper we make the assumption that the object does not undergo a deformation during its motion, that is:
Assumption
The attenuation coefficient is exposed to a rigid motion consisting of a rotation and a translation ,
(2.2) 
During such an induced motion, the object is illuminated from the direction with a uniform intensity, see Figure 1.
Assuming in a rough approximation that the light moves along straight lines and only suffers from attenuation, we consider measurements in the form of attenuation projection mappings.
Definition
The attenuation mapping maps a rigid body motion (described by translation and rotation ) applied to an object of interest (described by the attenuation function ) onto the attenuation projection image. That is the attenuation mapping of at time is given by
(2.3) 
We will do the reconstruction of the motion in Fourier space, using the convention
for the
dimensional Fourier transform. For deriving the adequate formulas we need to introduce some notation first:
Definition
The mapping
denotes the orthogonal projection onto the plane. Moreover, the adjoint of , is given by .
With this notation the attenuation mapping satisfies Equation 2.4, below, which is an identity similar to the Fourier projectionslice theorem, see for example [Bra56] and [Nat01, p.11], that is used quite heavily for orientation reconstruction in CryoEM and Ray tomography.
Lemma
Let and be the attenuation mapping of a rigid body motion , which the object of interest gets exposed to. Then, the following identity holds:
(2.4) 
Here denotes the twodimensional Fourier transform with respect to the coordinates .
Proof:
We denote by
(2.5) 
the Fourier transform of the attenuation mapping of exposed to a rigid transformation. By definition of and it holds that
Substituting , we obtain the asserted relation
This means that the Fourier transform of the twodimensional attenuation mapping of a rigid motion of an object of interest coincides, up to a factor depending on the motion , with the values of the threedimensional Fourier transform of evaluated on the central slice spanned by the first two columns of .
3 Motion Estimation
In this section we present a method for reconstruction of the time dependent translation and the rotation parameters of , respectively, from collected data of , as defined in Equation 2.5. For this purpose we proceed iteratively:

In a first step, explained in Section 3.1, we exploit the centers of the attenuation mappings in order to find the translations.

We then define in Equation 3.6 reduced attenuation mappings which do not depend on the translation anymore.

In Section 3.2 we use Lemma in order to find the remaining rotational part.
Instead of recovering the rotation matrix function directly, we rather retrieve the corresponding angular velocity function , which is defined as follows:
Definition
Let . Then the angular velocity corresponding to is defined via
(3.1) 
Moreover, we represent in cylindrical coordinates with cylindrical axis and denote by the cylindrical radius of , and by the cylindrical components of and the cylindrical angle of , such that
(3.2) 
Moreover, we set .
The general workflow is depicted in Figure 2.
3.1 Reconstruction of the Translation
Since the attenuation operator integrates along the direction, we are facing a translational invariance of the attenuation projection data in this direction. This makes it impossible to uniquely reconstruct the third component of the translation from the attenuation projection images:
Lemma
Let and be the attenuation mapping of a rigid motion of . Then,
Proof:
Substituting , we find that
According to Lemma it is impossible to reconstruct the third component of the translation of object of interest (that is the attenuation function), and thus we aim to recover only . It turns out that the orthogonal projection of the translated center of is equal to the center of the corresponding attenuation mapping , which can be directly computed from the projection data.
Proposition
Let and be the attenuation mapping of a rigid motion an object of interest is exposed to. Moreover, let be the center of defined in Equation 2.1 and the twodimensional center of the attenuation mapping ,
(3.3) 
Then,
(3.4) 
Proof:
To calculate the center of , we first remark that by substituting it follows that
(3.5) 
Moreover, we find
Substituting again , we get that
Since the first term vanishes according to the definition of , we get with Equation 3.5 the identity
which is just Equation 3.4.
Remark:
From Equation 3.4 we can reconstruct the translation part of in the plane. According to Lemma this is the most we can hope for.
3.2 Reconstruction of the Rotation
In a second step we aim at recovering the cylindrical parameters (, respectively), and corresponding to the rotation function from the reduced attenuation mapping , which is the Fourier transform of after translation correction with Equation 3.4:
Definition
We define the reduced attenuation map corresponding to as
(3.6)  
Lemma
The reduced attenuation mapping is independent of the translation and only dependent on the rotation .
Proof:
Multiplying Equation 2.4 with and taking into account the relation between and from Equation 3.4 we get
where the right hand side is independent of the translation .
The reduced attenuation mapping possesses the following symmetry property.
Lemma
Let and let be the reduced attenuation mapping. Then, for arbitrary the following identity holds
(3.7) 
for all and with .
Proof:
We remark that we have for arbitrary the relation
Choosing
with some arbitrary , we find from Equation 2.4 that
(3.8) 
Since by the definition of the cross product
the right hand side of Equation 3.8 is invariant with respect to interchanging and , and we get Equation 3.7.
The above result is the starting point for the recovery of the rotation matrix function from measurement of the attenuation mapping defined in Equation 2.3. The basic idea is to differentiate Equation 3.7 with respect to time up to order three. We see below that from the first order derivative of the equation it is possible to find the cylindrical component and the height of the angular velocity (see Proposition).
Reconstruction of the cylindrical component and the height
We show below that by differentiation of Equation 3.7 with respect to we can recover the cylindrical components and of the angular velocity as defined in Equation 3.2 and Equation 3.1
, respectively. In this section we use the tensor derivative notation as summarized in
Appendix A.First we derive the derivative of Equation 3.7:
Proposition
Let and as defined in Equation 3.6. Moreover, let and the associated angular velocity as defined in Equation 3.1. Then, for all satisfying and all the following relation holds:
(3.9) 
Proof:
We insert in Equation 3.7 the first order Taylor polynomials
with the coefficients and , , calculated in Lemma. Then, by differentiating Equation 3.7 with respect to at the position , we find for every the relation
Inserting the expressions from Equation A.2 for and , , we end up with the identity
If we choose for an with the parameter , this gives us Equation 3.9.
It is possible to recover for every time the cylindrical components and the cylindrical height using Proposition. Since Equation 3.9 holds for every it is sufficient look for a vector such that the function
(3.10) 
is constant for every . This is in principle possible as we show in Section 4. The value of this constant function will then be . Of course we have to assume that there does not exist a vector such that and for all . The question remains open under what conditions it is possible to determine , and consequently , uniquely? It turns out that we have to face a similar nonuniqueness issue as in CryoEM, where the reconstruction of the rotations is only possible up to an orthogonal transformation [Lam08]. This follows from the fact that the object of interest, in our case described by the attenuation coefficient , is unknown as well and a reflection of in the plane through the origin leads to the same attenuation projection data. In fact, this ambiguity can directly be observed from Equation 3.9, which is obviously solved by and . Denoting by and the corresponding angular velocities and by and the solutions to Equation 3.1 with initial conditions , we get the following connection
where . When the object is moved with respect to and the reflected object with respect to , we get exactly the same attenuation projection data (see Figure 3).
(B) and (D): Reflected object rotated with respect to at the same time steps.
The attenuation projection images are the same.
Reconstruction of the cylindrical radius
So far we have explained how to recover the cylindrical parameters (or equivalently ) and . For the full reconstruction of the motion will still need to estimate the cylindrical radius , which can be done by using Proposition below. Note that we go directly to the third derivative of Equation 3.7, as the second derivative provides no new information.
Proposition
Let , be the reduced attenuation mapping of a rigid motion of in Fourier space. Let further , and be the angular velocity corresponding to and let . Then, for all such that and , we have
(3.11) 
where
Note that we omitted the dependence on time of the coefficients for the sake of readability.
Proof:
The proof is provided in the Appendix A.
It is possible to recover for every time the parameter using Proposition. Since Equation 3.11 holds for every , we can consider it as an overdetermined linear system for and , see Section 4.
Example (The case )
Apparently, the coefficients and in Proposition vanish, if , in which case Equation 3.11 becomes trivial. Let us consider this case in more detail. For given cylindrical components , (, respectively), defined in Equation 3.2 we get in this special situation . Now, we calculate the rotation , defined in Equation 3.1 for given . To simplify the calculations, we assume that and . Solving Equation 3.1 with the initial condition then gives
We calculate
and circle back to Equation 3.7, which then reads
and holds for all . Since the term can be absorbed into the variable , this equation contains no information about . Therefore, it is not possible to recover this special motion completely with our technique.
4 Numerics
We consider for the points , , and the diagonal matrix as an example the attenuation coefficient
(4.1) 
and define the motion by the choice
with
We discretise the transformed function by evaluating it at the points
and calculate the attenuation projection along the third direction at every time step
Some attenuation projection images are depicted in Figure 4.
For the reconstruction, we proceed for each time step as follows:

We calculate the center of the projections and read off the first two components of the displacement via Equation 3.4, see Figure 5.

To reconstruct the third component as function of the yet unknown angular value, we interpret Equation 3.9 as a least square minimisation problem for the function
where we remark that
and its derivatives can be nicely interpolated, with the help of the Whittaker–Shannon interpolation formula, for example, as they are Fourier transforms of a function which is very rapidly decreasing. We call the minimisation point
. 
To obtain the angular function , we minimise the function
on with some tiny . Since the function is quite complicated and the minimisation problem is only onedimensional, see Figure 6, we chose to minimise by brute force, that is, we simply evaluate it for sufficiently many values and pick the one with smallest value. The minimiser gives us and thus , see Figure 7.
Figure 7: Reconstruction of on the left and on the right (the blue crosses are the reconstructed values, the red curve is the exact function). 
To obtain , we consider Equation 3.11 as overdetermined linear system (one equation for each value ) for and , where the coefficients can be explicitly calculated with the values of and obtained so far. Solving the corresponding normal equations, gives us the values , see Figure 8.
Since we chose a sufficiently high resolution in time and space, the simulations did not require any sort of regularisation (we did not enforce continuity of the reconstructions in any way) and serve more as a test for the correctness of the formulas, see the plot of the reconstruction errors in Figure 9.
Conclusions and Outlook
Tomographic reconstruction of a trapped object has the intrinsic problem that the object’s rotation is typically not as smooth and welldefined (and therefore known) as applied during tomographic data acquisition in other settings. Instead, the object undergoes a more irregular motion described as timedependent translations and rotations.
In the present work we deliver a first step into the direction of tomographic reconstruction of optically and/or acoustically trapped particles. We assume imaging at low numerical aperture of absorptive samples with amplitude contrast, which means that the attenuation projection describes the imaging process reasonably well.
For this case, we demonstrate—by explicit reconstruction—how the motional parameters can be inferred, wherever permitted uniquely. However, not all parameters are uniquely retrievable, for example, the component of the translation in the direction the projection images are taken cannot be seen, or, in case of symmetries of the sample, we encounter fake solutions of our equations for the angular velocity (in the extreme case of a spherically symmetric object, no information on the motion is available). More uniqueness studies are on the way.
The mathematical problem formulation reveals similarities to singleparticle CryoEM, however, here we can avoid the timeconsuming step of labelling of attenuation projection images, because the labels can be identified with a timestamp of the recordings over time. Nevertheless, one could envision that our method could serve as a corrector method, when images in CryoEM are labelled by neighboring projection images, and difference quotients are used to approximate artificial time derivatives. In the future, the proposed motion estimation will be tested on video data acquired from biological samples held in optical and/or acoustic traps, see Figure 10. Moreover, it will be necessary to study corrections or alternative approaches required when going from attenuation projection images to optical images. Image formation of amplitude or phasesamples in a microscope may significantly differ from attenuation projections, as explained in [Arr99, ArrScho09].
Acknowledgements
The authors are supported by the Austrian Science Fund (FWF), with SFB F68, project F6804N36 (Coupled Physics Imaging), project F6806N36 (Inverse Problems in Imaging of Trapped Particles), and project F6807N36 (Tomography with Uncertainties). The authors thank Mia Kvåle Løvmo and Benedikt Pressl for providing the video of the trapped pollen grain.
Appendix A Appendix
The main ingredient of our results for motion estimation in Section 3 are based on calculating higher order derivatives of composed functions. For this we require a tensor notation:
Tensor Derivatives
Let
and
In this paper we use derivatives of up to order of the composed function
These derivatives will be expressed via a tensor notation:

We denote by the derivative of order of the function with respect to the variable for fixed at a point . We write the evaluation of the tensor in the form

The time derivatives of the composed function up to order three then can be expressed with the simplifying notation , as follows:
(A.1) Note that for the tensor notation simplifies to the usual matrix vector multiplications:
Higher Order Expansion of Angular Velocities
Lemma (Taylor Polynomials)
Let , and be the angular velocity corresponding to . Then, we have
Comments
There are no comments yet.