Inverse Problems of Trapped Objects

TL;DR

This paper introduces a mathematical method for estimating the movement of transparent particles manipulated by optical or acoustical tweezers, drawing parallels to Cryo-EM imaging techniques, to improve contact-free microscopic imaging.

Contribution

It presents a novel mathematical framework for motion estimation of particles in optical and acoustical trapping, assuming fixed internal structure, inspired by Cryo-EM orientation estimation.

Findings

Provides a new mathematical formulation for particle motion estimation.

Shows similarities between trapping particle estimation and Cryo-EM orientation problems.

Lays groundwork for improved non-invasive microscopic imaging techniques.

Abstract

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. The mathematical formulation of the motion estimation shows some similarities to Cryo-EM single particle imaging, where the recording orientations of the probe need to be estimated.

Equations

$$C_{c,+}(\text{mathds}{R}^3;\text{mathds}R)=\{u\in C(\text{mathds}{R}^3;\text{mathds}R)\mid \mathrm{supp}(u)\ne \varnothing \text{is compact},u\ge 0\}$$

$${\mathcal{C}}_3:=\frac{1}{{\int }_{\text{mathds}{R}^3}u(x)\mathrm{d}x}\underset{\text{mathds}{R}^3}{\int }xu(x)\mathrm{d}x\in \text{mathds}{R}^3\text{ with }x=\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right).$$

$$A(t,x)={\mathcal{C}}_3+R(t)(x-{\mathcal{C}}_3+T(t)).$$

$$(T,R)\mapsto\mathcal{J}[T,R](t,x_{1},x_{2})=\int_{-\infty}^{\infty}u\big{(}\mathcal{C}_{3}+R(t)(x-\mathcal{C}_{3}+T(t))\big{)}\,\mathrm{d}x_{3}.$$

$${\mathcal{F}}_n[f](k)=(2\pi {)}^{-\frac{n}{2}}{\int }_{\text{mathds}{R}^n}f(x){\mathrm{e}}^{-\mathrm{i}⟨k,x⟩}\mathrm{d}x$$

$$P:\text{mathds}{R}^3$$

$$x$$

$${\mathcal{F}}_2[\mathcal{J}[T,R]]=\sqrt{2\pi }{\mathcal{F}}_3[u](R(t)P^{T}k){\mathrm{e}}^{\mathrm{i}⟨R(t)P^{T}k,{\mathcal{C}}_3⟩}{\mathrm{e}}^{\mathrm{i}⟨k,P(T(t)-{\mathcal{C}}_3)⟩}.$$

$$\hat{\mathcal{J}}:={\mathcal{F}}_2[\mathcal{J}[T,R]]$$

$$\hat{\mathcal{J}}$$

$$\displaystyle=\frac{1}{2\pi}\int_{\mathds{R}^{3}}u\big{(}\mathcal{C}_{3}+R(t)(x-\mathcal{C}_{3}+T(t))\big{)}\mathrm{e}^{-\mathrm{i}\left<P^{T}k,x\right>}\,\mathrm{d}x.$$

$$\hat{\mathcal{J}}$$

$$=\sqrt{2\pi }{\mathcal{F}}_3[u](R(t)P^{T}k){\mathrm{e}}^{\mathrm{i}⟨R(t)P^{T}k,{\mathcal{C}}_3⟩}{\mathrm{e}}^{\mathrm{i}⟨k,P(T(t)-{\mathcal{C}}_3)⟩}.$$

$$R(t{)}^T{R}^{\prime }(t)x=\omega (t)\times x\text{for all}x\in \text{mathds}{R}^3.$$

$$\omega (t)=\left(\begin{array}{c}\alpha (t)v(t)\\ {\omega }_3(t)\end{array}\right)=\left(\begin{array}{c}\alpha (t)\cos (\phi (t))\\ \alpha (t)\sin (\phi (t))\\ {\omega }_3(t)\end{array}\right).$$

$$\mathcal{J}[T,R]=\mathcal{J}[T+\rho e_3,R]\text{ for every }T\in C(\text{mathds}R;\text{mathds}R),\rho \in C(\text{mathds}R;\text{mathds}R)\text{ and }R\in C(\text{mathds}R;SO(3)).$$

$$\mathcal{J}[T+\rho e_3,R](t,x_1,x_2)$$

$$={\int }_{-\infty }^{\infty }u\left({\mathcal{C}}_3+R(t)\left(\left(\begin{array}{c}{x}_1\\ x_2\\ {\tilde{x}}_3\end{array}\right)-{\mathcal{C}}_3+T(t)\right)\right)\mathrm{d}{\tilde{x}}_3=\mathcal{J}[T,R](t,x_1,x_2).$$

$${\mathcal{C}}_2:=\frac{1}{{\int }_{\text{mathds}{R}^2}\mathcal{J}[T,R](t,x)\mathrm{d}x}\underset{\text{mathds}{R}^2}{\int }\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\mathcal{J}[T,R](t,x)\mathrm{d}x.$$

$$P({\mathcal{C}}_3-T(t))={\mathcal{C}}_2\text{for every}T\in C(\text{mathds}R;\text{mathds}{R}^3),R\in C(\text{mathds}R;SO(3)),t\in \text{mathds}R.$$

$$\int_{\mathds{R}^{2}}\mathcal{J}[T,R](t,x_{1},x_{2})\,\mathrm{d}(x_{1},x_{2})=\int_{\mathds{R}^{3}}u\big{(}\mathcal{C}_{3}+R(t)(x-\mathcal{C}_{3}+T(t))\big{)}\,\mathrm{d}x=\int_{\mathds{R}^{3}}u(y)\,\mathrm{d}y.$$

$$\int_{\mathds{R}^{2}}\begin{pmatrix}x_{1}\\ x_{2}\end{pmatrix}\mathcal{J}[T,R](t,x_{1},x_{2})\,\mathrm{d}(x_{1},x_{2})=\int_{\mathds{R}^{3}}Px\,u\big{(}\mathcal{C}_{3}+R(t)(x-\mathcal{C}_{3}+T(t))\big{)}\,\mathrm{d}x.$$

$${\int }_{\text{mathds}{R}^2}\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$$

$$=PR(t{)}^T\left({\int }_{\text{mathds}{R}^3}yu(y)\mathrm{d}y-{\mathcal{C}}_3{\int }_{\text{mathds}{R}^3}u(y)\mathrm{d}y\right)+P({\mathcal{C}}_3-T(t)){\int }_{\text{mathds}{R}^3}u(y)\mathrm{d}y$$

$${\int }_{\text{mathds}{R}^2}\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)\mathcal{J}[T,R](t,x_1,x_2)\mathrm{d}(x_1,x_2)=P({\mathcal{C}}_3-T(t)){\int }_{\text{mathds}{R}^2}\mathcal{J}[T,R](t,x_1,x_2)\mathrm{d}(x_1,x_2),$$

$$\tilde{\mathcal{J}}:\text{mathds}R\times \text{mathds}{R}^2$$

$$(t,k)$$

$$\tilde{\mathcal{J}}(t,k)=\sqrt{2\pi }{\mathcal{F}}_3[u](R(t)P^{T}k){\mathrm{e}}^{\mathrm{i}⟨R(t)P^{T}k,{\mathcal{C}}_3⟩},$$

$$\tilde{\mathcal{J}}\left(s,\frac{\lambda }{t-s}P(e_3\times (R(s{)}^{T}R(t)e_3))\right)=\tilde{\mathcal{J}}\left(t,\frac{\lambda }{s-t}P(e_3\times (R(t{)}^{T}R(s)e_3))\right)$$

$$(R(t)e_3)\times (R(s)e_3)=R(t)(e_3\times (R(t{)}^{T}R(s)e_3))=R(t)P^{T}P(e_3\times (R(t{)}^{T}R(s)e_3)).$$

$$k=\frac{\lambda }{s-t}P(e_3\times (R(t{)}^{T}R(s)e_3))\text{ for }s\ne t$$

$$\tilde{\mathcal{J}}\left(t,\frac{\lambda }{s-t}P(e_3\times (R(t{)}^{T}R(s)e_3))\right)=\sqrt{2\pi }{\mathcal{F}}_3[u]\left(\frac{\lambda }{s-t}(R(t)e_3)\times (R(s)e_3)\right){\mathrm{e}}^{\mathrm{i}\frac{\lambda }{s-t}⟨(R(t)e_3)\times (R(s)e_3),{\mathcal{C}}_3⟩}.$$

$$(R(t)e_3)\times (R(s)e_3)=-(R(s)e_3)\times (R(t)e_3),$$

Full text

1Faculty of Mathematics

University of Vienna

Oskar-Morgenstern-Platz 1

A-1090 Vienna, Austria

2Division for Biomedical Physics

Medical University of Innsbruck

Müllerstraße 44

A-6020 Innsbruck, Austria

3Johann Radon Institute for Computational

   and Applied Mathematics (RICAM)

Altenbergerstraße 69

A-4040 Linz, Austria

1. Introduction

In the past decades optical trapping has established itself as a versatile, and easy-to-implement tool for holding and moving microscopic ($\mu$m-sized) particles suspended in a liquid in a contact-free and non-invasive manner [7]. 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 [10].

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 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 for which the light propagation can be sufficiently well modelled by straight rays which only undergo spatially varying attenuation, as in X-ray tomography [9]. 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 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, see, for example, [6], (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 [5, 4]. 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 $\mathds{R}^{3}$, characterised by an attenuation coefficient $u\in C_{c,+}(\mathds{R}^{3};\mathds{R})$.

Definition \thedefinition.

In this paper

[TABLE]

denotes the space of admissible attenuation coefficients. For $u\in C_{c,+}(\mathds{R}^{3};\mathds{R})$ we define its center

[TABLE]

In this paper we make the assumption that the object does not undergo a deformation during its motion, that is:

Assumption \theassumption.

The attenuation coefficient $u$ is exposed to a rigid motion consisting of a rotation $R\in C(\mathds{R};SO(3))$ and a translation $T\in C(\mathds{R};\mathds{R}^{3})$,

[TABLE]

During such an induced motion, the object is illuminated from the $e_{3}$-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 \thedefinition.

The attenuation mapping $\mathcal{J}$ maps a rigid body motion (described by translation $t\mapsto T(t)$ and rotation $t\mapsto R(t)$) applied to an object of interest (described by the attenuation function $u:\mathds{R}^{3}\to\mathds{R}$) onto the attenuation projection image. That is the attenuation mapping of $u$ at time $t$ is given by

[TABLE]

We will do the reconstruction of the motion in Fourier space, using the convention

[TABLE]

for the $n$-dimensional Fourier transform. For deriving the adequate formulas we need to introduce some notation first:

Definition \thedefinition.

The mapping

[TABLE]

denotes the orthogonal projection onto the $x_{1}x_{2}-$plane. Moreover, the adjoint of $P$, $P^{T}:\mathds{R}^{2}\to\mathds{R}^{3}$ is given by $P^{T}k=(\begin{smallmatrix}k\\ 0\end{smallmatrix})$.

With this notation the attenuation mapping $\mathcal{J}$ satisfies Equation 2.4, below, which is an identity similar to the Fourier projection-slice theorem, see for example [3] and [9, p.11], that is used quite heavily for orientation reconstruction in Cryo-EM and $X$-Ray tomography.

Lemma 2.1**.**

Let $u\in C_{c,+}(\mathds{R}^{3};\mathds{R})$ and $\mathcal{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:

[TABLE]

*Here $\mathcal{F}_{2}$ denotes the two-dimensional Fourier transform with respect to the coordinates $(x_{1},x_{2})$. *

Proof*:*

We denote by

[TABLE]

the Fourier transform of the attenuation mapping of $u$ exposed to a rigid transformation. By definition of $\mathcal{\hat{J}}$ and $\mathcal{J}$ it holds that

[TABLE]

Substituting $y=\mathcal{C}_{3}+R(t)(x-\mathcal{C}_{3}+T(t))$, we obtain the asserted relation

[TABLE]

$\square$

This means that the Fourier transform of the two-dimensional attenuation mapping $\mathcal{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)$.

3. Motion Estimation

In this section we present a method for reconstruction of the time dependent translation $t\mapsto T(t)$ and the rotation parameters of $t\mapsto R(t)$, respectively, from collected data of $\mathcal{\hat{J}}$, as defined in Equation 2.5. For this purpose we proceed iteratively:

1. (i)

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

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

In Section 3.2 we use Lemma 3.3 in order to find the remaining rotational part.

Instead of recovering the rotation matrix function $t\mapsto R(t)$ directly, we rather retrieve the corresponding angular velocity function $t\mapsto\omega(t)$, which is defined as follows:

Definition \thedefinition.

Let $R\in C^{1}(\mathds{R};SO(3))$. Then the angular velocity $\omega\in C(\mathds{R};\mathds{R}^{3})$ corresponding to $R$ is defined via

[TABLE]

Moreover, we represent $\omega$ in cylindrical coordinates with cylindrical axis $e_{3}$ and denote by $\alpha\in C(\mathds{R};\mathds{R}^{+})$ the cylindrical radius of $\omega$, and by $v\in C(\mathds{R};\mathbb{S}^{1})$ the cylindrical components of $\omega$ and $\varphi\in C(\mathds{R};\mathds{R})$ the cylindrical angle of $\omega$, such that

[TABLE]

Moreover, we set $v^{\perp}(t)=(-v_{2}(t),v_{1}(t))^{T}$.

The general workflow is depicted in Figure 2.

3.1. Reconstruction of the Translation

Since the attenuation operator integrates along the $e_{3}$ 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 $t\mapsto T(t)$ from the attenuation projection images:

Lemma 3.1**.**

Let $u\in C_{c,+}(\mathds{R}^{3};\mathds{R})$ and $\mathcal{J}$ be the attenuation mapping of a rigid motion of $u$. Then,

[TABLE]

Proof*:*

Substituting $\tilde{x}_{3}=x_{3}+\rho(t)$, we find that

[TABLE]

$\square$

According to Lemma 3.1 it is impossible to reconstruct the third component of the translation $T$ of object of interest (that is the attenuation function), and thus we aim to recover only $P(T)$. It turns out that the orthogonal projection of the translated center of $u$ is equal to the center of the corresponding attenuation mapping $\mathcal{C}_{2}$, which can be directly computed from the projection data.

Proposition \theproposition.

Let $u\in C_{c,+}(\mathds{R}^{3};\mathds{R})$ and $\mathcal{J}$ be the attenuation mapping of a rigid motion an object of interest $u$ is exposed to. Moreover, let $\mathcal{C}_{3}$ be the center of $u$ defined in Equation 2.1 and $\mathcal{C}_{2}$ the two-dimensional center of the attenuation mapping $\mathcal{J}$,

[TABLE]

Then,

[TABLE]

Proof*:*

To calculate the center of $\mathcal{J}[T,R](t)$, we first remark that by substituting $y=\mathcal{C}_{3}+R(t)(x-\mathcal{C}_{3}+T(t))$ it follows that

[TABLE]

Moreover, we find

[TABLE]

Substituting again $y=\mathcal{C}_{3}+R(t)(x-\mathcal{C}_{3}+T(t))$, we get that

[TABLE]

Since the first term vanishes according to the definition of $\mathcal{C}_{3}$, we get with Equation 3.5 the identity

[TABLE]

which is just Equation 3.4. $\square$

Remark*:*

*From Equation 3.4 we can reconstruct the translation part of $T$ in the $x_{1}x_{2}$-plane. According to Lemma 3.1 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 $t\mapsto v(t)$ ($t\mapsto\varphi(t)$, respectively), $t\mapsto\alpha(t)$ and $t\mapsto\omega_{3}(t)$ corresponding to the rotation function $t\mapsto R(t)$ from the reduced attenuation mapping $\tilde{\mathcal{J}}$, which is the Fourier transform of $\mathcal{J}$ after translation correction with Equation 3.4:

Definition \thedefinition.

We define the reduced attenuation map corresponding to $u$ as

[TABLE]

Lemma 3.2**.**

*The reduced attenuation mapping is independent of the translation $t\mapsto T(t)$ and only dependent on the rotation $t\mapsto R(t)$. *

Proof*:*

Multiplying Equation 2.4 with $\mathrm{e}^{\mathrm{i}\left<k,\mathcal{C}_{2}\right>}$ and taking into account the relation between $\mathcal{C}_{2}$ and $\mathcal{C}_{3}$ from Equation 3.4 we get

[TABLE]

where the right hand side is independent of the translation $T$. $\square$

The reduced attenuation mapping possesses the following symmetry property.

Lemma 3.3**.**

Let $u\in C_{c,+}(\mathds{R}^{3};\mathds{R})$ and let $\tilde{\mathcal{J}}$ be the reduced attenuation mapping. Then, for arbitrary $R\in C(\mathds{R};SO(3))$ the following identity holds

[TABLE]

*for all $\lambda\in\mathds{R}$ and $s,t\in\mathds{R}$ with $s\neq t$. *

Proof*:*

We remark that we have for arbitrary $s,t\in\mathds{R}$ the relation

[TABLE]

Choosing

[TABLE]

with some arbitrary $\lambda\in\mathds{R}$, we find from Equation 2.4 that

[TABLE]

Since by the definition of the cross product

[TABLE]

the right hand side of Equation 3.8 is invariant with respect to interchanging $s$ and $t$, and we get Equation 3.7. $\square$

The above result is the starting point for the recovery of the rotation matrix function $t\mapsto R(t)$ from measurement of the attenuation mapping $\mathcal{J}$ 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 $v$ and the height $\omega_{3}$ of the angular velocity $\omega$ (see Section 3.2).

Reconstruction of the cylindrical component $v$ and the height $\omega_{3}$

We show below that by differentiation of Equation 3.7 with respect to $t$ we can recover the cylindrical components $v$ and $\omega_{3}$ of the angular velocity $\omega$ 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 \theproposition.

Let $u\in C_{c,+}(\mathds{R}^{3};\mathds{R})$ and $\tilde{\mathcal{J}}$ as defined in Equation 3.6. Moreover, let $R\in C^{2}(\mathds{R};SO(3))$ and $\omega\in C^{1}(\mathds{R};\mathds{R}^{3})$ the associated angular velocity as defined in Equation 3.1. Then, for all $t\in\mathds{R}$ satisfying $\alpha(t)\neq 0$ and all $\mu\in\mathds{R}$ the following relation holds:

[TABLE]

Proof*:*

We insert in Equation 3.7 the first order Taylor polynomials

[TABLE]

with the coefficients $a_{j}$ and $b_{j}$, $j=0,1$, calculated in Lemma A.1. Then, by differentiating Equation 3.7 with respect to $t$ at the position $t=s$, we find for every $\lambda\in\mathds{R}$ the relation

[TABLE]

Inserting the expressions from Equation A.2 for $a_{0}=\alpha v=b_{0}$ and $2a_{1}=\omega_{3}\alpha v^{\bot}+(\alpha v)^{\prime}$, $2b_{1}=-\omega_{3}\alpha v^{\bot}+(\alpha v)^{\prime}$, we end up with the identity

[TABLE]

If we choose for an $s\in\mathds{R}$ with $\alpha(s)\neq 0$ the parameter $\lambda=\frac{\mu}{\alpha(s)}$, this gives us Equation 3.9. $\square$

It is possible to recover for every time $t\in\mathds{R}$ the cylindrical components $v(t)$ and the cylindrical height $\omega_{3}(t)$ using Section 3.2. Since Equation 3.9 holds for every $\mu\in\mathds{R}$ it is sufficient look for a vector $v(t)\in\mathbb{S}^{1}$ such that the function

[TABLE]

is constant for every $\mu\in\mathds{R}$. This is in principle possible as we show in Section 4. The value of this constant function will then be $\omega_{3}(t)$. Of course we have to assume that there does not exist a vector $v\in\mathbb{S}^{1}$ such that $\partial_{t}[\tilde{\mathcal{J}}](t,\mu v)=0$ and ${\bf D}^{1}[\tilde{\mathcal{J}}](t,\mu v)[\![\mu v^{\perp}]\!]=0$ for all $\mu\in\mathds{R}$. The question remains open under what conditions it is possible to determine $v$, and consequently $\omega_{3}$, uniquely? It turns out that we have to face a similar non-uniqueness issue as in Cryo-EM, where the reconstruction of the rotations is only possible up to an orthogonal transformation [8]. This follows from the fact that the object of interest, in our case described by the attenuation coefficient $u$, is unknown as well and a reflection of $u$ in the $x_{1}x_{2}$-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 $v(t)$ and $\check{v}(t)=-v(t)$. Denoting by $\omega(t)$ and $\check{\omega}(t)$ the corresponding angular velocities and by $R(t)$ and $\check{R}(t)$ the solutions to Equation 3.1 with initial conditions $R(0)=\check{R}(0)=\mathds{1}$, we get the following connection

[TABLE]

where $\Sigma=\mathrm{diag}(1,1,-1)$. When the object is moved with respect to $R$ and the reflected object with respect to $\check{R}$, we get exactly the same attenuation projection data (see Figure 3).

Reconstruction of the cylindrical radius $\alpha$

So far we have explained how to recover the cylindrical parameters $t\mapsto v(t)$ (or equivalently $t\mapsto\varphi(t)$) and $t\mapsto\omega_{3}(t)$. For the full reconstruction of the motion will still need to estimate the cylindrical radius $t\mapsto\alpha(t)$, which can be done by using Section 3.2 below. Note that we go directly to the third derivative of Equation 3.7, as the second derivative provides no new information.

Proposition \theproposition.

Let $u\in C_{c,+}(\mathds{R}^{3};\mathds{R})$, $\tilde{\mathcal{J}}$ be the reduced attenuation mapping of a rigid motion of $u$ in Fourier space. Let further $R\in C^{4}(\mathds{R};SO(3))$, $t\in\mathds{R}$ and $\omega\in C^{3}(\mathds{R};\mathds{R}^{3})$ be the angular velocity corresponding to $R$ and let $\sigma(t)=\varphi^{\prime}(t)$. Then, for all $t\in\mathds{R}$ such that $\alpha(t)\neq 0$ and $\sigma(t)\neq-\omega_{3}(t)$, we have

[TABLE]

where

[TABLE]

*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. $\square$

It is possible to recover for every time $t\in\mathds{R}$ the parameter $\alpha(t)$ using Section 3.2. Since Equation 3.11 holds for every $\mu\in\mathds{R}$, we can consider it as an overdetermined linear system for $\alpha^{2}(t)$ and $\frac{\alpha^{\prime}(t)}{\alpha(t)}$, see Section 4.

Example** (The case $\omega_{3}=-\sigma$).**

Apparently, the coefficients $A_{0},A_{02}$ and $A_{1}$ in Section 3.2 vanish, if $\omega_{3}=-\sigma$, in which case Equation 3.11 becomes trivial. Let us consider this case in more detail. For given cylindrical components $t\mapsto\omega(t)$, $t\mapsto v(t)$ ($t\mapsto\varphi(t)$, respectively), defined in Equation 3.2 we get in this special situation $\sigma(t)=\varphi^{\prime}(t)=-\omega_{3}(t)$. Now, we calculate the rotation $R$, defined in Equation 3.1 for given $t\mapsto\omega(t)$. To simplify the calculations, we assume that $\alpha^{\prime}=0$ and $\omega_{3}^{\prime}=0$. Solving Equation 3.1 with the initial condition $R(0)=\mathds{1}$ then gives

[TABLE]

We calculate

[TABLE]

and circle back to Equation 3.7, which then reads

[TABLE]

and holds for all $\lambda\in\mathds{R}$. Since the term $\sin(\alpha(t-s))$ can be absorbed into the variable $\lambda$, this equation contains no information about $\alpha$. Therefore, it is not possible to recover this special motion completely with our technique.

4. Numerics

We consider for the points $P_{1}=(1,\frac{1}{2},-1)$, $P_{2}=(-\frac{1}{2},1,1)$, $P_{3}=(0,-1,\frac{1}{2})$ and the diagonal matrix $D=\mathrm{diag}(\sqrt{2},1,1)$ as an example the attenuation coefficient

[TABLE]

and define the motion $T$ by the choice

[TABLE]

with

[TABLE]

We discretise the transformed function by evaluating it at the points

[TABLE]

and calculate the attenuation projection along the third direction at every time step

[TABLE]

Some attenuation projection images are depicted in Figure 4.

For the reconstruction, we proceed for each time step as follows:

1. (i)

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

2. (ii)

To reconstruct the third component $\omega_{3}(t)$ as function of the yet unknown angular value, we interpret Equation 3.9 as a least square minimisation problem for the function

[TABLE]

where we remark that $\tilde{\mathcal{J}}$ 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 $\tilde{\omega}_{3}(\varphi(t))$. 3. (iii)

To obtain the angular function $\varphi(t)$, we minimise the function

[TABLE]

on $[0,\pi)$ with some tiny $\varepsilon>0$. Since the function is quite complicated and the minimisation problem is only one-dimensional, see Figure 6, we chose to minimise by brute force, that is, we simply evaluate it for sufficiently many values $\phi$ and pick the one with smallest value. The minimiser gives us $\varphi(t)$ and thus $\omega_{3}(t)=\tilde{\omega}_{3}(\varphi(t))$, see Figure 7.

4. (iv)

To obtain $\alpha$, we consider Equation 3.11 as overdetermined linear system (one equation for each value $\mu\in\{j\delta_{x}\mid j\in\{-512,\ldots,511\}\}$) for $\alpha^{2}(t)$ and $\frac{\alpha^{\prime}(t)}{\alpha(t)}$, where the coefficients can be explicitly calculated with the values of $\varphi$ and $\omega_{3}$ obtained so far. Solving the corresponding normal equations, gives us the values $\alpha(t)$, 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 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 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 single-particle Cryo-EM, however, here we can avoid the time-consuming step of labelling of attenuation projection images, because the labels can be identified with a time-stamp of the recordings over time. Nevertheless, one could envision that our method could serve as a corrector method, when images in Cryo-EM 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 phase-samples in a microscope may significantly differ from attenuation projections, as explained in [1, 2].

Acknowledgements

The authors are supported by the Austrian Science Fund (FWF), with SFB F68, project F6804-N36 (Coupled Physics Imaging), project F6806-N36 (Inverse Problems in Imaging of Trapped Particles), and project F6807-N36 (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

[TABLE]

and

[TABLE]

In this paper we use derivatives of up to order $3$ of the composed function

[TABLE]

These derivatives will be expressed via a tensor notation:

• •

We denote by ${\bf D}^{i}[f](t,\kappa):\underbrace{\mathds{R}^{2}\times\mathds{R}^{2}\cdots\mathds{R}^{2}}_{i\text{ times}}\to\mathds{C}$ the derivative of order $i=0,1,2,\ldots$ of the function $k\in\mathds{R}^{2}\mapsto f(t,k)$ with respect to the variable $k$ for fixed $t$ at a point $\kappa\in\mathds{R}^{2}$. We write the evaluation of the tensor ${\bf D}^{i}[f](t,\kappa)$ in the form

[TABLE]

• •

The time derivatives of the composed function $t\to h(t)$ up to order three then can be expressed with the simplifying notation ${\bf D}^{i}[f]={\bf D}^{i}[f](t,g(t))$, as follows:

[TABLE]

Note that for $i=1,2$ the tensor notation simplifies to the usual matrix vector multiplications:

[TABLE]

Higher Order Expansion of Angular Velocities

Lemma A.1** (Taylor Polynomials).**

Let $R\in C^{4}(\mathds{R};SO(3))$, $t\in\mathds{R}$ and $\omega\in C^{3}(\mathds{R};\mathds{R}^{3})$ be the angular velocity corresponding to $R$. Then, we have

[TABLE]

with remainder terms $a,b$ fulfilling $\lim_{t\to s}\frac{\tilde{a}(s,t)}{(t-s)^{3}}=\lim_{t\to s}\frac{\tilde{b}(s,t)}{(t-s)^{3}}=0$ and with the coefficients

[TABLE]

Proof*:*

We abbreviate the matrix-valued function $R^{\mathrm{T}}R^{\prime}$ by $B$ and calculate from the identity $R^{\prime}=RB$ further derivatives by differentiation:

[TABLE]

With this, we can calculate the coefficients

[TABLE]

where we can use the antisymmetry $B=-B^{\mathrm{T}}$ for the calculatuon of the coefficients $b_{j}$.

Plugging in Equation 3.1 to write $Bx=\omega\times x$ and applying the identity $f\times(g\times h)=\left<f,h\right>g-\left<f,g\right>h$ for $f,g,h\in\mathds{R}^{3}$, we find

[TABLE]

Inserting $P(\omega)=\alpha v$ and using Taylor’s theorem, we arrive at the stated Taylor polynomials. $\square$

Proof* (of Section 3.2):*

We use again in Equation 3.7 the Taylor polynomials

[TABLE]

with the coefficients $a_{j}$ and $b_{j}$, $j=0,1$, calculated in Lemma A.1. Then, by taking three derivatives with respect to $t$ and evaluating at $t=s$ (omitting the argument $s$ in the coefficients $a$ and $b$) we find with the notation of Equation A.1 that

[TABLE]

We calculate some of the terms separately:

• •

The third derivatives of $\tilde{\mathcal{J}}$ yield

[TABLE]

• •

For the second derivatives of $\tilde{\mathcal{J}}$, we find

[TABLE]

• •

For the first derivatives of $\tilde{\mathcal{J}}$, we get with

[TABLE]

that

[TABLE]

• •

On the other side, we find for the mixed derivatives of $\tilde{\mathcal{J}}$

[TABLE]

The terms with the second order derivatives in $\alpha$ are given by

[TABLE]

for $\alpha\neq 0$, since, if we take the derivative of Equation 3.9 with respect to $\lambda$, we find that

[TABLE]

Taking the time derivative of Equation 3.9, we arrive at

[TABLE]

Another derivative with respect to time gives us

[TABLE]

Subtracting the identity (by our choice of the functions $v$ and $\omega_{3}$) of Equation A.4 from Equation A.3, we find, by putting everything together, the condition

[TABLE]

We sort the terms with respect to the derivative $\alpha^{\prime}$ and write the equation in the form

[TABLE]

• •

The coefficient $A_{12}$ of $(\alpha^{\prime})^{2}$ is given by

[TABLE]

We take the derivative of Equation 3.9 with respect to $\lambda$, giving us

[TABLE]

and taking another derivative with respect to $\lambda$, we find the identity

[TABLE]

Comparing this with Equation A.6, we see that for $\alpha\neq 0$:

[TABLE]

• •

The coefficient $A_{1}$ of $\alpha^{\prime}$ is

[TABLE]

If we compare this with the time derivative of the derivative of Equation 3.9 with respect to $\mu$ with $\lambda=\frac{\mu}{\alpha}$, that is, the time derivative of

[TABLE]

which is

[TABLE]

then we see that some of the coefficients cancel each other and we are left with

[TABLE]

We rewrite it in the form

[TABLE]

Remarking that $|v|^{2}=1$ and thus $\left<v,v^{\prime}\right>=0$, we have a scalar function $\sigma$ such that

[TABLE]

With this, $A_{1}$ becomes (using $(v^{\perp})^{\prime}=(v^{\prime})^{\perp}=-\sigma v$)

[TABLE]

Thus, $A_{1}(\mu)$ is zero if either $\sigma=-\omega_{3}$, that is $v^{\prime}=-\omega_{3}v^{\perp}$, in which case the bracket can be written as

[TABLE]

or if the bracket vanishes, that is,

[TABLE]

Taking the time derivative of Equation 3.9 with $\lambda=\frac{\mu}{\alpha}$, we get

[TABLE]

• •

The coefficient $A_{02}$ of $\alpha^{2}$ consists of the terms

[TABLE]

With Equation A.7 this becomes

[TABLE]

which again vanishes for $\sigma=-\omega_{3}$.

• •

Finally, the coefficient $A_{0}$ is given by

[TABLE]

Writing again $v^{\prime}=\sigma v^{\perp}$, this becomes

[TABLE]

Subtracting Equation A.8, we get rid of one term:

[TABLE]

$\square$