Motion reconstruction for optical tomography of trapped objects

Optical and acoustical trapping has been established as a tool for holding and moving microscopic particles suspended in a liquid in a contact-free and non-invasive manner. Opposed to standard microscopic imaging, where the probe is fixated, this technique allows imaging in a more natural environment. This paper provides a method for estimating the movement of a transparent particle which is maneuvered by tweezers (assuming that the inner structure of the probe is not subject to local movements) by making use of the assumption of a smooth movement in time. The mathematical formulation of the motion estimation leads to an infinitesimal version of the common line technique used in cryogenic electron microscopy single particle imaging to estimate the orientations of the particles in the probe.


Introduction
In the past decades optical trapping has established itself as a versatile, and easy-to-implement tool for holding and moving microscopic (µm-sized) particles suspended in a liquid in a contact-free and non-invasive manner [8]. 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 cell-clusters, entire micro-organisms, or so-called 'organoids', mini-organs grown in the petri-dish) other trapping modalities, such as acoustic forces, have to be added as auxiliary trapping forces [14].
In this paper, we study the problem of estimating the movement of trapped particles, as a preprocessing step for 3D tomographic imaging. Conveniently, the required tomographic projection images from different view-points 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 where the light propagation can be sufficiently well modeled by straight rays which only undergo spatially varying attenuation; as in the model of x-ray tomography. In this case the optical image resembles a projection image. This is, for instance, fulfilled in low numer ical aperture imaging of biological samples with sufficient amplitude contrast and with large structures compared to the scale of the optical wavelength. For high-resolution imaging with a high numerical aperture objective, or for samples with small structures diffracting the light beams significantly, deviations from geometrical optics are to be expected. Then, optical projection tomography based on the Fourier slice theorem will have to be replaced by optical diffraction tomography based on the Fourier diffraction theorem, as derived in [10].
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 smooth, but complicated motion described by an affine transformation with time-dependent-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 three-dimensional x-ray transform (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, see, for example, [7].
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 single-particle electron microscopy (cryo-EM), see [5,6]. An infinitesimal version of this common line 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 from an infinitesimally short observation time, see section 3.2. As a numerical verification of the algorithm, we provide in section 4 reconstructions of simulated data.

Mathematical model
We consider a bounded object in R 3 , characterized by an attenuation coefficient u ∈ C c,+ (R 3 ; R). Definition 2.1. We denote by C c,+ (R 3 ; R) = {u ∈ C(R 3 ; R) | supp(u) = ∅ is compact, u 0} the space of admissible attenuation coefficients. For u ∈ C c,+ (R 3 ; R) we define its center In this paper we make the assumption that the object does not undergo a deformation during its motion, that is: The attenuation coefficient u is exposed to a rigid motion consisting of a rotation R ∈ C(R; SO(3)) and a translation T ∈ C(R; R 3 ), During such an induced motion, the object is illuminated from the direction e 3 (denoting by e 1 , e 2 , e 3 the standard basis vectors in R 3 ) with a uniform intensity, see figure 1.
In the experiments, the induced motion of the object turned out to be sufficiently smooth on every relevant time scale. As a first approximation for the problem, we will assume to be in a regime where the light propagation can be sufficiently well described by geometrical optics so that we can consider the light to move along straight lines suffering only from attenuation. Then, the measurements can be modeled in the form of attenuation projection mappings. Definition 2.3. The attenuation mapping J maps a rigid motion (described by a translation t → T(t) and a rotation t → R(t)) applied to an object of interest (described by the attenuation function u : R 3 → R) onto the attenuation projection image. That is the attenuation mapping of u is given by (2.3) We will do the reconstruction of the motion in Fourier space, using the convention for the n-dimensional Fourier transform. And it is convenient to use the orthogonal projection on the x 1 x 2 -plane and its adjoint P T : R 2 → R 3 , given by P T k = (k, 0) T , to shorten the notation. In Fourier space, the attenuation mapping J satisfies equation (2.4), below, which is an identity similar to the Fourier projection-slice theorem, see for example [3] and [11, p 11], that is used quite heavily for orientation reconstruction in cryo-EM and x-ray tomography. It also shows that if the motion (T, R) is known, then the Fourier transform of J directly allows us to recover the Fourier transform of the attenuation coefficient u for the corresponding directions, thus providing us a way to fully reconstruct the attenuation coefficient.

Lemma 2.4.
Let u ∈ C c,+ (R 3 ; R) and J [R, T] be the attenuation mapping of a rigid body motion (R, T), which the object of interest u gets exposed to. Then, the following identity holds: (2.4) Here, F 2 denotes the two-dimensional Fourier transform with respect to the coordinates (x 1 , x 2 ).

Proof. We denote bŷ
the Fourier transform of the attenuation mapping of u exposed to a rigid transformation. By the definition of J , see equation (2.3), we then have that which, upon inserting the definition of F 3 [u], is exactly equation (2.4). □ This means that the Fourier transform of the two-dimensional attenuation mapping J of a rigid motion (T, R) of an object of interest u coincides, up to a factor depending on the motion A, with the values of the three-dimensional Fourier transform of u evaluated on the central slice spanned by the first two columns of R(t).

Motion estimation
In this section we present a method for the reconstruction of the horizontal components of the time dependent translation t → T(t) and the rotation t → R(t) from the collected data J [T, R], as defined in equation (2.3). For this purpose we (i) show in section 3.1 how to determine from the centers of the attenuation mappings J [T, R] the two components T 1 and T 2 of the translation; (ii) define accordingly in equation (3.7) the reduced attenuation mapping J which does not depend on the translation T anymore; and (iii) use in section 3.2 the first and third order terms of our formulation of the common line property, see lemma 3.6, to determine the rotation R.
Instead of recovering the rotation matrix function t → R(t) directly, we rather retrieve the corre sponding angular velocity function t → ω(t), which is defined as follows. (3)). Then the angular velocity ω ∈ C(R; R 3 ) corresponding to R is defined via the relation where R (t) denotes the derivative with respect to the time t.
Conveniently, we represent ω in cylindrical coordinates with cylindrical axis e 3 and denote by α ∈ C(R; R + ) the cylindrical radius of ω, and by v ∈ C(R; S 1 ) the cylindrical components of ω, and by ϕ ∈ C(R; R) the cylindrical angle of ω, such that The reconstruction algorithm thus follows the steps shown in figure 2.

Reconstruction of the translation
Since the attenuation operator integrates along the direction e 3 , we have an invariance of the attenuation projection data with respect to translations of the object in this direction. Proof. Substituting It is therefore impossible to reconstruct the third component of the translation T of the object of interest (that is the attenuation function), and thus we aim to recover only the two components T 1 and T 2 . Now, the orthogonal projection of the center of the translated and rotated object u is equal to the center of the corresponding attenuation mapping, which can be directly computed from the projection data.

Proposition 3.3.
Let u ∈ C c,+ (R 3 ; R) and J be the attenuation mapping of a rigid motion an object of interest u is exposed to. Moreover, let C 3 ∈ R 3 be the center of u defined in equation (2.1) and C 2 : R → R 2 be the two-dimensional center of the attenuation mapping J [T, R] for some translation T ∈ C(R; R 3 ) and some rotation R ∈ C(R; SO(3)), that is, Then, we have Proof. To calculate the center of J [T, R], we first remark that the denominator in equation (3.3) is simply the total attenuation of the object, which can be seen explicitly by inserting the definition of J , see equation (2.3), and by substituting y (3.5) In the same way, we find Since the first term vanishes according to the definition of C 3 , see equation (2.1), we get with equation (3.5) the identity which is just equation (3.4). □

From equation (3.4), we immediately get that
which determines the components T 1 and T 2 up to a global translation P(T(0)) (which is due to the fact that the placement of the object u is not fixed in any way).

Reconstruction of the rotation
In a second step we aim at recovering the cylindrical corresponding to the rotation function t → R(t) from the reduced attenuation mapping J , which is the Fourier transform of J after translation correction with equation (3.4): Definition 3.4. We define the reduced attenuation mapping corresponding to an admissible attenuation coefficient u under a motion (T, R) as Lemma 3.

The reduced attenuation mapping of an admissible attenuation coefficient u under a motion (T, R) is independent of the translation t → T(t).
Proof. Multiplying equation (2.4) with e i k,C2(t) and taking into account the relation be- for the reduced attenuation mapping, which is seen to be independent of the translation T. □ The reduced attenuation mapping possesses the following symmetry property, which is the common line property known from methods in cryo-EM, see, for example, [5,6]. Lemma 3.6. Let u ∈ C c,+ (R 3 ; R) and R ∈ C(R; SO (3)). Then, the reduced attenuation mapping J of u under the rotation R fulfills the relatioñ Proof. We remark that we have for arbitrary s, t ∈ R the relation Choosing for s = t the vector with some arbitrary λ ∈ R in equation (3.8), we find that (3.10) Now, the right hand side of equation (3.10) is because of the skew-symmetry of the cross product invariant with respect to interchanging s and t, and we therefore get equation (3.9). □ Remark 1. Since lim t→s P(e 3 × (R(s) T R(t)e 3 )) = 0, we chose the scaling λ t−s for the free scalar parameter in equation (3.9) so that the argument λ t−s P(e 3 × (R(s) T R(t)e 3 )) converges in the later considered limit t → s in general (that is, for α(s) = 0) to a finite value.
Based on this identity, the relative orientation of the object at different times can be determined, which is also extensively used in cryo-EM, see, for example, [4,13,15]. In the following, we want to show a way to explore the smooth time dependence to calculate from this relation locally at every time the change in movement (that is, calculating the angular velocity ω) by taking time derivatives of equation (3.9), which is an information not available in cryo-EM.

Reconstruction of the cylindrical and third components of the angular velocity.
We show below that by differentiation of equation (3.9) with respect to t we can recover the cylindrical components v and the third component ω 3 of the angular velocity ω as defined in equations (3.2) and (3.1), respectively. In this section we use the tensor derivative notation as summarized in the appendix (see section A.1 in the appendix).

Proposition 3.7.
Let u ∈ C c,+ (R 3 ; R) be an admissible attenuation coefficient and R ∈ C 2 (R; SO(3)) be a rotation with associated angular velocity ω ∈ C 1 (R; R 3 ) defined by equation (3.1). Let further J be the reduced attenuation map corresponding to u under the rotation R as defined in equation (3.7).
Then, we have for all t ∈ R satisfying α(t) = 0 and all µ ∈ R the relation Proof. We insert in equation (3.9) the first-order Taylor polynomials If we choose for an s ∈ R with α(s) = 0 the parameter λ = µ α(s) , this gives us equation (3.11). □ Thus, the cylindrical components v(t) and the third component ω 3 (t) satisfy for every time t ∈ R the equation (3.11) and we can therefore attempt to reconstruct these parameters by looking for a vector ṽ ∈ S 1 such that the function is finite and constant on the set of all values µ ∈ R for which numerator and denominator do not both vanish. Ideally, there would exist a unique vector ṽ with this property, which would then necessarily give us the cylindrical components v(t) =ṽ of the angular velocity and the constant value of the function in equation (3.13) would equal to the third component: Now, the uniqueness can at best only be up to a choice of sign, since with v(t) also −v(t) fulfills the property of proposition 3.7 (with the same value ω 3 (t)). In fact, this non-uniqueness is already present in equation (3.9), since it is not only invariant under global rotations (3) (reflecting that the orientation of the attenuation coefficient can be arbitrarily chosen), but also under the orientation changing coordinate transform figure 3 for an illustration, and for which the angular velocity ω corresponding to the transformed rotation Ř is, because ofŘ , ω 3 (t)) T . This nonuniqueness is also observed in cryo-EM, where the reconstruction of the rotations is only possible up to a global orthogonal transformation [9].
In general, one may, of course, also have other solutions; in the extreme case of a spherically symmetric attenuation coefficient, for example, the attenuation projection data contains no information on the orientation of the object and every rotation R will fulfill equation (3.9). However, for a sufficiently asymmetric object, we can expect a unique solution ṽ up to sign so that Ψṽ is constant, which is also indicated by our numerical tests, see section 4.

3.2.2.
Reconstruction of the cylindrical radius of the angular velocity. So far we have explained how to recover the cylindrical parameters t → v(t) (or equivalently t → ϕ(t)) and t → ω 3 (t). For the full reconstruction of the motion, we will still need to estimate the cylindrical radius t → α(t), which can be done by using propositions 3.8 below. Note that we go directly to the third derivative of equation (3.9), since the second derivative provides no new information. Proposition 3.8. Let u ∈ C c,+ (R 3 ; R), J be the reduced attenuation mapping of a rigid motion of u in Fourier space. Let further R ∈ C 4 (R; SO(3)), t ∈ R and ω ∈ C 3 (R; R 3 ) be the angular velocity corresponding to R and let σ(t) = ϕ (t). Then, for all t ∈ R such that α(t) = 0 and σ(t) = −ω 3 (t), we have where Here, we omitted for the sake of readability the arguments of the functions v, v ⊥ , σ, and ω 3 and their derivatives: they are always to be evaluated at the time t.
Proof. The proof is provided in appendix. □ It is possible to recover for every time t ∈ R the parameter α(t) using propositions 3.8. Since equation (3.14) holds for every µ ∈ R, we can consider it as an overdetermined linear system for α 2 (t) and α (t) α(t) , see section 4. Example. Apparently, the coefficients A 0 , A 02 , and A 1 in propositions 3.8 vanish if ω 3 = −σ, in which case equation (3.14) becomes trivial and thus a reconstruction of α from this equation impossible.
We want to study this in the simple case of an angular velocity ω for which the cylindrical radius α and the third component ω 3 are constant functions. Then, the critical condition In this case, the rotation R can be explicitly calculated by solving equation (3.1) and we find with the normalization R(0) = 1 that We therefore have which means that the common line property, equation (3.9), upon which the whole reconstruction is based, reduces tõ for all λ ∈ R. Since the term sin(α(t − s)) can be absorbed into the variable λ, this equation contains no information about α. Therefore, it is not possible to recover this special motion completely from the common line property alone.

Numerical simulation
Finally, we want to give a proof of principle that the results in section 3 allow us under some conditions to reconstruct the motion correctly. To this end, we choose an attenuation coefficient without obvious symmetries to prevent multiple solutions of equation (3.11) for the components v and ω 3 of the angular velocity ω = (αv 1 , αv 2 , ω 3 ) T . We thus consider for the points P 1 = (1, 1 2 , −1) T , P 2 = (− 1 2 , 1, 1) T , P 3 = (0, −1, 1 2 ) T and the diagonal matrix D = diag( √ 2, 1, 1) as an example the attenuation coefficient and define the motion (T, R) by the arbitrary choice where we pick for the rotation angles the functions ψ 1 (t) = t 2 , ψ 2 (t) = t + 3t 3 , and ψ 3 (t) = sin(3t). Here, we denoted by R 2 (φ) and R 3 (φ) the rotations by the angle φ ∈ R around the standard basis vectors e 2 and e 3 , respectively: for all x ∈ R 3 , we can write the derivative of the rotation R in the form for every x ∈ R 3 . This gives us according to equation (3.1) the angular velocity which we can write in cylindrical coordinates to obtain the parameters α, ϕ, and ω 3 as defined in equation (3.2); see figure 7 for a plot of the graphs of these functions. We discretize the transformed function x → u(C 3 + R(t)(x − C 3 + T(t))), where C 3 denotes the center of u, by evaluating it at the points Following the reconstruction formulas presented in propositions 3.3, 3.7, and 3.8, we will attempt to reconstruct the first two components T 1 and T 2 of the translation T and the rotation R, parametrized by the functions α, ϕ, and ω 3 , from the data J [T, R]. In view of the instability of this reconstruction algorithm due to the differentiations of the reduced attenuation mapping J and of the recovered versions of the functions ϕ and ω 3 with respect to time, see propositions 3.7 and 3.8, we consider here as a basic test of the practicability of the reconstruction only the ideal case of data without noise and with a sufficiently fine discretization so that we avoid the problems caused by numerical differentiation. Under more realistic settings, appropriate regularization methods will have to be used to mitigate this ill-posedness, but this is out of the scope of this article. For the reconstruction of T and ω (from which we can then recover the rotation R by solving the differential equation in equation (3.1) up to the arbitrary initial rotation R(0)), we proceed for each time step t = δ t as follows.
We therefore define ω 3 (t, φ) for each φ ∈ [0, π) as the minimum point of the least squares functional with the (for noiseless data arbitrary) cutoff J 0 = 50, which is given explicitly bỹ ; so that we have The values of J and its derivatives, which appear in this formula, can be nicely interpolated (with the help of the Whittaker-Shannon interpolation formula, for example) as they are Fourier transforms of very rapidly decreasing functions. (iii) To recover the cylindrical angle ϕ(t), we remark that the minimum of the functional W 3 from equation (4.5), which is for every choice of φ attained at the point ω 3 (t, φ), should theoretically become zero at the correct value ϕ(t). We therefore reconstruct ϕ(t) as the minimum point of the functional Since Φ t is in general non-convex, but only a one-dimensional function, we choose to simply minimize it on the sufficiently fine grid {mδ φ | m ∈ {0, . . . , 2047}}, δ φ = π 2048 , by evaluating it at all points (and by additionally refining the grid in the vicinity of the minimum point to improve the accuracy), see figure 6 for a sketch of the graph of Φ t at the time step 50δ t . From the so reconstructed value of ϕ(t), we then calculate the value of ω 3 (t) via equation (4.6). (iv) To obtain α, we consider equation (3.14) as an overdetermined linear system (one where the coefficients A 0 , A 02 , and A 1 can be explicitly calculated with the values of ϕ and ω 3 obtained so far. This means, we reconstruct the vector (α 2 (t), α (t) α(t) ), and thus α(t) > 0, as the minimum point of the functional with the matrix A(t) ∈ C (2J0+1)×2 and the vector b(t) ∈ C 2J0+1 given by for j ∈ {−J 0 , . . . , J 0 }. Explicitly, this minimum point is calculated as a solution of the normal equation The results of this reconstruction algorithm (applied to the motion specified in equations (4.2) and (4.3)) for the parameters α, ϕ, and ω 3 of the angular velocity ω are shown in figures 7 and 8. They agree up to the expected error caused by the discretization well with the exact values calculated from equation (4.4).

Comparison to ASPIRE
There exists a large number of software tools or software applications that are used in the field of single-particle cryo-EM, which also take on the task of orientation estimation. One of them is the ASPIRE package [12], which provides common-line based algorithms for orientation assignment. Here, we want to draw a comparison between (the MATLAB ® version of) ASPIRE and our algorithm using the simulated data defined via equations (4.1), (4.2) and (4.3). Since we are not familiar with all the details of ASPIRE, some of the obtained results might still be improvable.
The orientation estimation done by ASPIRE, which is based on [15], takes as an input a stack of projection images and gives back the corresponding rotation matrices. We have calculated the translations with our algorithm and used the reduced attenuation projection images independent of the translations as an input for ASPIRE, so that only the rotations are left to be determined. As it is designed to solve problems from cryo-EM, the orientation estimation algorithm of ASPIRE expects projection images from different viewing angles (class averages) that are not very similar to each other. As the temporal step size in our generated data is very small, we only used every 105th image for the reconstruction. This seemed to yield the best results. Figure 6. The evaluated points for the graph of the function Φ t on the initial grid at the time step t = 50δ t are shown on the left. On the right, we plot the function on the refined grid around the minimum point, which is estimated to be at 1.555844 and is our numerical approximation for the exact value ϕ(t) ≈ 1.555847. In figure 9 the reconstructed rotation matrices are depicted via their Euler angle representations. Finally, in figure 7 the corresponding angular velocity parameters are shown next to the exact as well as our reconstructed values. We used a Gaussian filter for computing the derivatives of the rotation matrices from ASPIRE, which explains the deviations of the reconstructed values close to the boundaries of the time interval. The fact that our algorithm takes, in contrast to ASPIRE, the temporal coherence of the data into account is clearly an advantage. However, being built for problems from cryo-EM, the ASPIRE package is designed to deal with a lot of noise, while our algorithm is based on numerical differentiation and is therefore very sensitive to faulty data. An algorithm that combines the best of both worlds would be ideal for data from optical microscopy of trapped specimens. We plan to introduce appropriate regularization methods to our algorithm, which should improve performance with noisy data while still utilizing the time component.

Conclusions and outlook
Tomographic reconstruction of a trapped object has the intrinsic problem that the object's rotation is typically not as smooth and well-defined (and therefore known) as applied during tomographic data acquisition in other settings. Instead, the object undergoes a more irregular motion described as time-dependent translations and rotations.
In the present work we deliver a first step into the direction of tomographic reconstruction of optically 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 single-particle cryo-EM. However, in our trapping application, we have the additional information that there is a continuous relation between the time and the orientation of the object. In the future, the proposed motion estimation will be tested on video data acquired from biological samples held in optical 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 phase-samples in a microscope may significantly differ from attenuation projections, as explained in [1,2]. Note that for i = 1, 2 the tensor notation simplifies to the usual matrix vector multiplications: , .

A.2. Higher order expansion of angular velocities
As we would like to expand the common line property in lemma 3.6 in a Taylor polynomial, we first calculate the relevant Taylor coefficients of the arguments of the function J therein. (3)), t ∈ R and ω ∈ C 3 (R; R 3 ) be the angular velocity corresponding to R. As in equation (3.2), we write ω = (αv, ω 3 ) T with α ∈ C 3 (R; R + ), v ∈ C 3 (R; S 1 ), and ω 3 ∈ C 3 (R; R). Then, we have a 0 = αv, Proof. According to Taylor's theorem, the coefficients a j and b j , j = 0, 1, 2, 3, are given by To calculate the derivatives of R, we abbreviate the matrix-valued function R T R by B and calculate from the identity R = RB the higher derivatives Plugging in equation (3.1) to write Bx = ω × x and applying the identity (A.6) • For the first derivatives of J , we get with • And for the mixed derivatives of J , we have Looking only at the terms including a second derivative of α, we see that they cancel out:  We sort the terms with respect to the term α and write the equation in the form A 0 (s, λα) + α 2 A 02 (s, λα) + A 1 (s, λα)λα + A 12 (s, λα)(λα ) 2 = 0. (A.12) • The coefficient A 12 of (α ) 2 is then given by (A.13) By taking another derivative with respect to λ of equation (A.10), we find the identity Comparing this with equation (A.13), we see that • The coefficient A 1 of α is seen to be then we see that some of the coefficients cancel each other and we are left with We rewrite this in the form Remarking that |v| 2 = 1 and thus v, v = 0, we have a scalar function σ such that v = σv ⊥ . With this, A 1 becomes (using (v ⊥ ) = (v ) ⊥ = −σv) • The coefficient A 02 of α 2 consists of the terms