# Inverse Problems 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. 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.

## Authors

• 3 publications
• 2 publications
• 15 publications
• 1 publication
08/20/2018

### Inverse Problems in Asteroseismology

Asteroseismology allows us to probe the internal structure of stars thro...
04/13/2018

### Understanding Soft-Tissue Behavior for Application to Microlaparoscopic Surface Scan

This paper presents an approach for understanding the soft tissue behavi...
07/27/2020

### A concept of a measuring system for probe kinesthetic parameters identification during echocardiography examination

Echocardiography is the most commonly used imaging technique in clinical...
12/09/2017

### A cable-driven parallel manipulator with force sensing capabilities for high-accuracy tissue endomicroscopy

This paper introduces a new surgical end-effector probe, which allows to...
02/04/2017

### Fast and easy blind deblurring using an inverse filter and PROBE

PROBE (Progressive Removal of Blur Residual) is a recursive framework fo...
01/19/2019

### Endoscopic vs. volumetric OCT imaging of mastoid bone structure for pose estimation in minimally invasive cochlear implant surgery

Purpose: The facial recess is a delicate structure that must be protecte...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 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 [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 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 [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 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 [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 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,

[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

 Cc,+(R3;R)={u∈C(R3;R)∣supp(u)≠∅is compact,u≥0}

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 ,

 A(t,x)=C3+R(t)(x−C3+T(t)). (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

 (T,R)↦J[T,R](t,x1,x2)=∫∞−∞u(C3+R(t)(x−C3+T(t)))dx3. (2.3)

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

 Fn[f](k)=(2π)−n2∫Rnf(x)e−i⟨k,x⟩dx

for the

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

###### Definition

The mapping

 P:R3 →R2 x ↦(x1x2)

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 projection-slice theorem, see for example [Bra56] and [Nat01, p.11], that is used quite heavily for orientation reconstruction in Cryo-EM 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:

 F2[J[T,R]]=√2πF3[u](R(t)PTk)ei⟨R(t)PTk,C3⟩ei⟨k,P(T(t)−C3)⟩. (2.4)

Here denotes the two-dimensional Fourier transform with respect to the coordinates .

###### Proof:

We denote by

 ^J:=F2[J[T,R]] (2.5)

the Fourier transform of the attenuation mapping of exposed to a rigid transformation. By definition of and it holds that

 ^J =12π∫R2J[T,R](t,x1,x2)e−ik1x1−ik2x2d(x1,x2) =12π∫R3u(C3+R(t)(x−C3+T(t)))e−i⟨PTk,x⟩dx.

Substituting , we obtain the asserted relation

 ^J =12π∫R3u(y)e−i⟨PTk,R(t)T(y−C3)−(T(t)−C3)⟩dy

This means that the Fourier transform of the two-dimensional 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 three-dimensional 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:

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

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

3. 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

 R(t)TR′(t)x=ω(t)×xfor allx∈R3. (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

 ω(t)=(α(t)v(t)ω3(t))=⎛⎜⎝α(t)cos(φ(t))α(t)sin(φ(t))ω3(t)⎞⎟⎠. (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,

 J[T,R]=J[T+ρe3,R] for every T∈C(R;R),ρ∈C(R;R) and R∈C(R;SO(3)).

###### Proof:

Substituting , we find that

 J[T+ρe3,R](t,x1,x2) =∫∞−∞u(C3+R(t)(x+ρ(t)e3−C3+T(t)))dx3 =∫∞−∞u(C3+R(t)((x1x2~x3)−C3+T(t)))d~x3=J[T,R](t,x1,x2).

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 two-dimensional center of the attenuation mapping ,

 C2:=1∫R2J[T,R](t,x)dx∫R2(x1x2)J[T,R](t,x)dx. (3.3)

Then,

 P(C3−T(t))=C2for everyT∈C(R;R3),R∈C(R;SO(3)),t∈R. (3.4)

###### Proof:

To calculate the center of , we first remark that by substituting it follows that

 (3.5)

Moreover, we find

 ∫R2(x1x2)J[T,R](t,x1,x2)d(x1,x2)=∫R3Pxu(C3+R(t)(x−C3+T(t)))dx.

Substituting again , we get that

 ∫R2(x1x2) J[T,R](t,x1,x2)d(x1,x2)=∫R3P(R(t)T(y−C3)+C3−T(t))u(y)dy =PR(t)T(∫R3yu(y)dy−C3∫R3u(y)dy)+P(C3−T(t))∫R3u(y)dy

Since the first term vanishes according to the definition of , we get with Equation 3.5 the identity

 ∫R2(x1x2)J[T,R](t,x1,x2)d(x1,x2)=P(C3−T(t))∫R2J[T,R](t,x1,x2)d(x1,x2),

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

 ~J:R×R2 →R, (3.6) (t,k) ↦F2[J[T,R]](t,k)ei⟨k,C2⟩.

###### 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

 ~J(t,k)=√2πF3[u](R(t)PTk)ei⟨R(t)PTk,C3⟩,

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

 ~J(s,λt−sP(e3×(R(s)TR(t)e3)))=~J(t,λs−tP(e3×(R(t)TR(s)e3))) (3.7)

for all and with .

###### Proof:

We remark that we have for arbitrary the relation

 (R(t)e3)×(R(s)e3)=R(t)(e3×(R(t)TR(s)e3))=R(t)PTP(e3×(R(t)TR(s)e3)).

Choosing

 k=λs−tP(e3×(R(t)TR(s)e3)) for s≠t

with some arbitrary , we find from Equation 2.4 that

 ~J(t,λs−tP(e3×(R(t)TR(s)e3)))=√2πF3[u](λs−t(R(t)e3)×(R(s)e3))eiλs−t⟨(R(t)e3)×(R(s)e3),C3⟩. (3.8)

Since by the definition of the cross product

 (R(t)e3)×(R(s)e3)=−(R(s)e3)×(R(t)e3),

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 v and the height ω3

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:

 ∂t~J(t,μv(t))=ω3(t)D1[~J](t,μv(t))[[μv⊥(t)]]. (3.9)

###### Proof:

We insert in Equation 3.7 the first order Taylor polynomials

 1t−sP(e3×(R(s)TR(t)e3)) =a0(s)+a1(s)(t−s)+o(t−s)and 1s−tP(e3×(R(t)TR(s)e3)) =b0(s)+b1(s)(t−s)+o(t−s)

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

 D1[~J](s,λa0(s))[[λa1(s)]]=∂t~J(s,λb0(s))+D1[~J](s,λb0(s))[[λb1(s)]].

Inserting the expressions from Equation A.2 for and , , we end up with the identity

 ∂t~J(s,λα(s)v(s))=D1[~J](s,λα(s)v(s))[[λω3(s)α(s)v⊥(s)]] for all λ,s∈R.

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

 μ↦∂t~J(t,μv(t))D1[~J](t,μv(t))[[μv⊥(t)]] (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 non-uniqueness issue as in Cryo-EM, 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

 ˇR(t)=ΣR(t)Σ,

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).

#### 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

 A0(μ)+A02(μ)α(t)2+A1(μ)μα′(t)α(t)=0for allμ∈R, (3.11)

where

 A0(μ) =14μ(ω3+σ)[μ2ω3(ω3−σ)D3[~J](s,μv)[[v⊥,v⊥,v⊥]] +2μD2[~J](s,μv)[[v⊥,ω3σv−ω′3v⊥]]+2D1[~J](s,μv)[[ω23v⊥+ω′3v]] −μ(3ω3−σ)∂tD2[~J](s,μv)[[v⊥,v⊥]]+2∂ttD1[~J](s,μv)[[v⊥]]], A02(μ) =12μ(ω3+σ)D1[~J](s,μv)[[v⊥]], A1(μ) =12(ω3+σ)[μω3D2[~J](s,μv)[[v⊥,v⊥]]−ω3D1[~J](s,μv)[[v]]−∂tD1[~J](s,μv)[[v⊥]]].

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 ω3=−σ)

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

 R(t)=⎛⎜⎝cos(ω3t)−sin(ω3t)0−sin(ω3t)cos(αt)−cos(ω3t)cos(αt)sin(αt)sin(ω3t)sin(αt)cos(ω3t)sin(αt)cos(αt)⎞⎟⎠.

We calculate

 P(e3×(R(s)TR(t)e3))=sin(α(t−s))(cos(ω3s)−sin(ω3s))

and circle back to Equation 3.7, which then reads

 ~J(s,λt−ssin(α(t−s))(cos(ω3s)−sin(ω3s)))=~J(t,λs−tsin(α(t−s))(−cos(ω3t)sin(ω3t)))

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

 u(x)=3∏i=1|x−Pi|2e−14|Dx|2 (4.1)

and define the motion by the choice

 a(t)=⎛⎜⎝cos(6t)cos(12t)cos(6t)sin(12t)sin(t)⎞⎟⎠ and ω(t)=⎛⎜⎝α(t)cos(φ(t))α(t)sin(φ(t))ω3(t)⎞⎟⎠

with

 α(t)=1+10t2,φ(t)=πt+π3, and ω3(t)=12+√12+5t.

We discretise the transformed function by evaluating it at the points

 (j1,j2,j3)δx,j∈{−512,…,511}2×{−256,…,255} with δx=0.05

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

 ℓδt,ℓ∈{0,…,999} for δt=0.0005.

Some attenuation projection images are depicted in Figure 4.

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

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

2. 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

.

3. 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 one-dimensional, 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.

4. 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 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 [Arr99, ArrScho09].

### 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

 f:R×R2 →C, (t,k) ↦f(t,k)

and

 g:R →R2. t ↦(g1(t)g2(t))

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

 h:R →R. t ↦h(t):=f(t,g(t))

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

 Di[f](t,κ)[[v1,v2,⋯,vi]].
• The time derivatives of the composed function up to order three then can be expressed with the simplifying notation , as follows:

 h′(t)= D1[f][[g′]](t)+∂tf(t,g(t)) (A.1) h′′(t)= D2[f][[g′,g′]](t)+∂tD1[f][[g′]](t)+D1[f][[g′′]](t)+∂ttf(t,g(t)) h′′′(t)= D3[f][[g′,g′,g′]](t)+∂tD2[f][[g′,g′]](t)+∂ttD1[f][[g′]](t)+ 3(D2[f][[g′,g′′]](t)+∂tD1[f][[g′′]](t))+D1[f][[g′′′]](t)+∂tttf(t,g(t)).

Note that for the tensor notation simplifies to the usual matrix vector multiplications:

 D1[f][[g′]](t) =∇k[f](t,g(t))⋅g′(t), D2[f][[g′,g′]](t) =g′(t)T∇2k[f](t,g(t))g′(t).

### Higher Order Expansion of Angular Velocities

###### Lemma (Taylor Polynomials)

Let , and be the angular velocity corresponding to . Then, we have

 1t−sP(e3×R(s)TR(t)e3) =3∑j=0aj(s)(t−s)j+~a(s,t), 1s−tP(e3×R(t)TR