Dynamic Discrete Tomography

1. Introduction

In the following, the goal is to determine the paths $P_{1}, \dots, P_{n}$ of $n$ particles in space over a period of $t \in N$ moments in time from images taken by a a fixed number $m$ of cameras. Each spot in each camera image is the projection of a particle perpendicular to the plane of the corresponding camera. The number of particles that project on the same spot can be detected from the brightness of the image. Hence, in effect, for each moment in time we have the information how many particles lie on each line perpendicular to the planes of the cameras. We will refer to this information as the X-ray images of the set of particles in $m$ directions; details will be given in Section 2.

This problem lies at the core of dynamic discrete tomography, and we will refer to it as tomographic point tracking or tomographic particle tracking. As it turns out, the problem comprises two different but coupled basic underlying tasks, the reconstruction of a finite set of points from few of their X-ray images (discrete tomography) and the identification of the points over time (tracking). For surveys on various aspects of discrete tomography see [28, 31, 32, 34, 35]. Tracking is known in the literature as data association or object tracking and, for more specific applications, as multi-target tracking, multisensor surveillance or, in physics, as particle tracking; see [14, 21, 41, 42, 43, 44, 45, 50] and, e.g., [1, 16, 18] for surveys.

Tomographic particle tracking was already considered in [24, 26, 38, 40, 54] and [8]. In fact the linear programming

based heuristic introduced in

[8] was successfully applied to determine 3D-slip velocities of a gliding arc discharge in [56]. Another linear programming approach, different from [8, 56], was recently proposed in [24].

In the present paper we study the problem from a mathematical and algorithmic point of view with a special focus on the interplay between discrete tomography and tracking. Therefore we will distinguish the cases that for none, some or all of the $t$ moments $τ_{1}, \dots, τ_{t}$ in time, a solution of the discrete tomography task at time $τ_{1}, \dots, τ_{t}$ is explicitly available (and is then considered the correct solution regardless whether it is uniquely determined by its X-rays). For $d \geq 3$ , two (affine) lines in $R^{d}$ in general position are disjoint. Hence, generically, X-ray lines meet only in points of the underlying set. Therefore even the latter situation is of considerable practical relevance. We will refer to it as the positionally determined case while, in the more general situation, we will speak of the (partially) or (totally) tomographic case of point tracking.

In the following, we discuss several models for dynamic discrete tomography, give various algorithms and complementing $N P$ -hardness results. In particular, we show for the positionally determined case that the tracking problem can be solved in polynomial time if it exhibits a certain Markov-type property (which, effectively, allows only dependencies between any two consecutive time steps); see Theorem 1. The partially tomographic case, however, is $N P$ -hard even for $t = 2$ ; see Theorem 2. Complementing this result, we consider the tomographic tracking problem for two directions where a so-called displacement field is assumed to be given (the displacement field uniquely determines the particle’s next position). Again, this problem turns out to be $N P$ -hard already for $t = 2$ and certain realistic classes of displacement fields; see Theorem 3.

We then turn to the rolling horizon approach introduced in [8] which proceeds successively from step to step. After giving a short account of this method in Section 3.4 we will show that, while being quite successful in practice, it is not guaranteed to always yield the correct solution; see Example 1. Then, in Section 3.5 we study the issue of of how to incorporate additional prior knowledge of particle history into the models. Under rather general assumptions we show that already in the positionally determined case the tracking problems becomes $N P$ -hard for $t > 2;$ see Theorem 4. In particular, we show that even if the particles are known to move along straight lines, this prior knowledge cannot efficiently be exploited algorithmically (unless $P = N P$ ); see Theorem 5 and Corollary 3 and 4.

We proceed by introducing three algorithms that can be viewed as rather general paradigms of heuristics that involve prior knowledge about the movement of the particles and which can be used in the tomographic case.

We then discuss combinatorial models. In these models the positions of the particles in the next time step are assumed to be known approximatively in the sense that the candidate positions are confined to certain windows, which are finite subsets of positions. Again, under rather general conditions, we show $N P$ -hardness of the respective tomographic tracking tasks; see Theorem 8. However, we also identify polynomial-time solvable special cases of practical relevance; see Theorem 9 and 10.

The paper is organized as follows: Section 2 will introduce the relevant notation, provide some basic background, discuss various modeling issues including the question of how to rigorously incorporate a notion of ‘physically most likely’ solution, and state our main results.

Then main focus of the present paper lies on optimization models. They are studied in quite some detail in Section 3. In particular, we consider Markov-type integer programming and rolling horizon models and discuss the issue of utilizing the particle history from a structural algorithmic point of view. Also, we determine the computational complexity of the problems. We prove various $N P$ -hardness results, both, for the partially tomographic and the totally tomographic case, in the latter, even when the displacement field is known. Then we give various heuristic algorithms (or, to be more precise, classes of algorithmic paradigms) which are designed with a view towards balancing the running time and the achieved quality of the solutions. In the prevailing presence of $N P$ -hardness such a balancing is mandatory for a successful use in practical applications (of realistic size).

Section 4 deals with a different way of encoding a priori information about the movement of the particles. The studied combinatorial models utilize such information in terms of bounds on the number of particles in certain windows. We provide $N P$ -hardness results but also identify polynomially solvable classes of instances. The paper concludes with some final remarks in Section 5.

2. Notation, basics and main results

2.1. Discrete Tomography

Let $R,$ $Q,$ $Z,$ and $N$ denote the set of reals, rational numbers, integers, and natural numbers, respectively. For $q \in N$ , let $[q] = {1, \dots, q}$ and $[q]_{0} = [q] \cup {0}$ . With $11$

we denote the all-ones vector. Further, let

d, m \in N

d \geq 2

, and set

F^{d} = {F : F \subseteq Q^{d} \land F is finite} .

The point sets $F \in F^{d}$ model the sets of particles (in a static environment). Let us point out that, mathematically, the restriction to $Q^{d}$ is merely a tribute to the model of computation that we are going to apply. In fact, we will use the

binary Turing machine model

where each element is encoded in binary and hence its size is the number of binary digits needed for this representation; see e.g. [30]. (The binary size of a positive integer

q

is therefore essentially

log (q)

.) Note that this restriction is also in accordance with standard practical measurements which do not allow infinite precision.

With $S^{d}$ we denote the set of all lattice lines, i.e., the set of $1$ -dimensional linear subspaces of $R^{d}$ that are spanned by vectors from $Z^{d}$ . (Note that as $S^{d}$ coincides with the set of lines spanned by rational vectors, there is no further restriction here.) The lines in $S^{d}$ are (in a slight abuse of language) often referred to as directions. In fact, in an experimental environment, these lines are the ‘viewing directions’ under which the cameras ‘see’ the particles i.e., they are perpendicular to the image planes of the cameras.

For $S \in S^{d},$ we set $A (S) = {v + S : v \in Q^{d}} .$ Then, for $F \in F^{d}$ and $S \in S^{d},$ the (discrete $1$ -dimensional) X-ray of $F$ parallel to $S$ is the function

X_{S} F : A (S) \to N_{0} = N \cup {0}

defined by

X_{S} F (T) = | F \cap T | = \sum x \in T χ_{F} (x)

for each $T \in A (S),$ where, as usual, $| F \cap T |$ denotes the cardinality of $F \cap T$ and $χ_{F}$

is the characteristic function of

F

Let us point out that this definition is the basis of the paramount grid model of discrete tomography (see [28, 31, 32, 34, 35] which we will use throughout the paper. Of course, the term X-ray is meant generically here, i.e., does not necessarily refer to a specific imaging technique. However, the main motivation underlying the present paper is that of tracking physical particles which are ‘visible’ only through the camera image of their projections.

Two sets $F_{1}, F_{2} \in F^{d}$ are called tomographically equivalent with respect to $S_{1}, \dots, S_{m} \in S^{d}$ if

X_{S_{i}} F_{1} = X_{S_{i}} F_{2},

for all $i \in [m] .$

Given $m$ different lines $S_{1}, \dots, S_{m} \in S^{d}$ , the X-ray data is given in terms of functions

f_{i} : A (S_{i}) \to N_{0} (i \in [m]),

with finite support $T_{i} \subseteq A (S_{i})$ , represented by appropriately chosen data structures; see [29].

Suppose, we are given an instance $(f_{1}, \dots, f_{m})$ of measurements of an otherwise unknown set $F \in F^{d} .$ Then, of course, $| F | = ∥ f_{1} ∥_{(1)} = \dots = ∥ f_{m} ∥_{(1)}$ , where $∥ \cdot ∥_{(1)}$ is the usual $1$ -norm. From the data functions we can infer that $F$ must be contained in the grid

G = m ⋂ i = 1 ⋃ T \in T_{i} T

of candidate points. Note that, in general, even when the given instance has a unique solution $F,$ the grid $G$ can be a proper superset of $F .$

2.2. Dynamic Discrete Tomography

In tomographic point tracking, we consider $t$ consecutive moments in time. Hence, in particular, we are interested in the sets $F^{(τ)}$ of particles at each moment in time $τ \in [t]$ . Since, in general, these sets are accessible only through their X-rays taken from a given finite set of directions, they represent $t$ uncoupled instances of discrete tomography.

Of course, we also need to track the particles over time. Let $P = {p_{1}, \dots, p_{n}}$ denote the (abstract) set of $n$ particles. Then, for each $τ \in [t]$ , we are actually interested in a one-to-one mapping $π^{(τ)} : P \to F^{(τ)}$ that identifies the points of $F^{(τ)}$ with the particles. Hence, the particle tracks are given by $P_{i} = (π^{(1)} (p_{i}) \dots, π^{(t)} (p_{i})),$ $i \in [n] .$ We will refer to this identification as coupling.

Note that the coupling between two consecutive moments in time bears the character of a matching in a bipartite graph. In tomographic point tracking, however, the tasks of discrete tomography and matching strongly interact.

As it is well-known that for $m \geq 3$ already the static reconstruction problem of discrete tomography is $N P$ -hard [29] (see also [15]), highly instable [3, 9], and since in practical experiments space for installing cameras is typically strictly limited (see, e.g, [8]) we will mainly focus on the case $m = 2$ . In real-world particle tracking this corresponds to images taken from two cameras which may be regarded as two planes in $R^{3}$ orthogonal to two different directions $S_{1}$ , $S_{2}$ , respectively.

For $m = 2,$ a set $F^{(τ)}$ satisfying the X-ray constraints with respect to the two given directions $S_{1}$ and $S_{2}$ can be efficiently determined for each $τ \in [t]$ ; see [27, 46]. Also, it can be checked efficiently, whether each $F^{(τ)}$ is unique. (For related stability results see [4, 51].) Since a set is only very rarely uniquely determined by just two of its X-rays, the reconstruction will typically rely on additional physical knowledge. (For various uniqueness results the reader is referred to [33] and the literature quoted therein.)

But even if we know for each $τ \in [t]$ the correct set $F^{(τ)},$ we still have to identify the paths of all individual particles over time. This turns the $t$ otherwise uncoupled systems into a single coupled system. Here we will in particular discuss the question of how to utilize additional information via the coupling, which thus might reduce the ambiguity of the uncoupled tomographic tasks.

But how can the goal of determining a ‘physically most likely’ solution be modeled? Let $Π (n, t)$ denote the set of all $n$ -tuples $(P_{1}, \dots, P_{n})$ of potential particle tracks over $t$ moments of time for fixed point sets $F^{(τ)}$ . If the given tomographic particle tracking data is feasible there may, of course, exist many different but tomographically equivalent solutions for each $τ$ . In any case, $| Π (n, t) |$ is then still a lower bound for the number of different potential solutions. Now note that $| Π (n, t) | = (n!)^{t}$ , and even if the order of the individual paths is irrelevant, this number reduces only to $(n!)^{t - 1},$ which is still exponential in $n$ and $t$ .

Of course, this means that any enumerative algorithm will fail in practice even for quite moderate numbers of particles and time steps. But even worse, we cannot encode as input the ‘physical costs’ $c (P_{1}, \dots, P_{n})$ of all $(P_{1}, \dots, P_{n})$ explicitly. Even if we assume that all sets $F^{(τ)}$ are given and such costs can be computed as a simple function of the costs $w (P_{i})$ of the individual paths, e.g., as

n \sum i = 1 w (P_{i}) or max i \in [n] w (P_{i}),

the number of different costs $w (P)$ of a particle path (which still need to be available to determine the cost of a solution $(P_{1}, \dots, P_{n})$ ) reduces only to $n^{t} .$

In the following we will therefore in general refrain from assuming that such values are explicit parts of the input. We will instead assume that an algorithm is available which computes for any solution $(P_{1}, \dots, P_{n})$ its cost $c (P_{1}, \dots, P_{n})$ in time that is polynomial in all the other input data. In some situations such an algorithm will be based on additional specific data which reflect how appropriate certain (local) choices are and which will then be regarded as be part of the input, too. Such an algorithm $O$ will be called an objective function oracle. In the special case that $O$ provides the values $w (P_{i})$ and $c (P_{1}, \dots, P_{n}) = \sum_{i = 1}^{n} w (P_{i}),$ the oracle will be referred to as path value oracle.

For all practical algorithms, such an oracle will of course be specified explicitly. (It will typically be based on weights on grid points or pairs of grid points which augment the input of our problem.)

2.3. Basic algorithmic problems

Algorithmically, our general problem of tomographic particle tracking for $m$ given different directions $S_{1}, \dots, S_{m} \in S^{d},$ and based on an objective function oracle $O$ , can now be formulated as follows.

TomTrac $(O; S_{1}, \dots, S_{m})$ .

$t \in N,$ X-ray data functions $f_{1}^{(τ)}, \dots, f_{m}^{(τ)}$ for $τ \in [t]$ with $∥ f_{1}^{(τ)} ∥_{(1)} = \dots = ∥ f_{m}^{(τ)} ∥_{(1)} = n$ .
Decide whether, for each $τ \in [t]$ , there exists a set $F^{(τ)} \in F^{d}$ such that $X_{S_{i}} F^{(τ)} = f^{(τ)}$ for all $i \in [m]$ . If so, find particle tracks $P_{1}, \dots, P_{n}$ of minimal cost for $O$ among all couplings of all tomographic solutions $F^{(1)}, \dots, F^{(t)}$ .

Let us add a remark concerning the assumption, that the $1$ -norms of the X-ray data functions in TomTrac $(O; S_{1}, \dots, S_{m})$ are known and assumed to coincide. First, note that the latter condition is quite natural if all particles are detected at each moment in time. Also, both conditions can be checked in polynomial time. Since in practice, cameras have a finite field of view it may, however, happen in an experimental setting that several previously detected particles may move out of the field of view of the camera. Also particles may appear (or reappear) during the measurements. In this situation we might at least aim at obtaining partial results by reconstructing partial tracks or by including dummy nodes, which account for invisible particles (for a discussion, see, e.g., [36]). In the following, our prime focus will lie, however, on the case that all $n$ particles are recorded at each moment in time $τ \in [t] .$

In the following we will distinguish the cases that for none, some or all $τ \in [t]$ , the correct solution $F^{(τ)}$ is explicitly given. The former will be referred to as the (partially) or (totally) tomographic case while we speak of the latter as positionally determined. In the positionally determined case the problem TomTrac $(O; S_{1}, \dots, S_{m})$ reduces to the following problem.

Trac $(O; d)$ .

$t \in N$ , sets $F^{(1)}, \dots, F^{(t)} \in F^{d}$ with $| F^{(1)} | = \dots = | F^{(t)} | = n$ .
Find particle tracks $P_{1}, \dots, P_{n}$ of minimal cost for $O$ among all couplings of the sets $F^{(1)}, \dots, F^{(t)}$ .

Note that each instance of Trac $(O; d)$ can be viewed as a $t$ -dimensional assignment problem. For a comprehensive survey on assignment problems see [16] and the references quoted therein.

As pointed out before, for $d \geq 3$ , the positionally determined case is the generic situation even for $m = 2$ as two lines parallel to $S_{1}$ and $S_{2}$ do only meet in points of $F^{(τ)}$ . However, for $d = 2$ , any two non-parallel lines in the plane meet in a point, whence forming a candidate grid $G^{(τ)}$ of size $O (n^{2})$ . (Note that, generally, no fixed number $m$ of X-ray images suffices for unique determination.) Also, even in $R^{d}$ with $d \geq 3,$ there are particle constellations for which the candidate grids $G^{(τ)}$ do not consist of just $n$ points. And, of course, if the resolution of the images is low, such constellations become relevant even in practice.

In some applications, there is strong prior knowledge about the possible couplings available. Most notably is the case where a displacement field is prescribed. Formally a particle displacement field (for the time step $τ \to τ + 1$ ) is a pair $(p, φ^{(τ)} (p) - p)$ with $φ^{(τ)} : Q^{d} \to Q^{d}$ denoting a map where $φ^{(τ)} (p) - p$ denotes the displacement vector for a particle located at position $p$ at time $τ .$ In particular, if $π^{(τ)} (p_{i})$ denotes the position of the $i$ -th particle at time $τ,$ then $π^{(τ + 1)} (p_{i}) = φ^{(τ)} (π^{(τ)} (p_{i})),$ for $i \in [n]$ and $τ \in [t - 1] .$ Now, the tomographic tasks for the time step $τ \to τ + 1$ consists of determining the positions of the particles at time $τ$ such that the solution at time $τ + 1$ determined by the displacement field matches the tomographic data for $τ + 1.$ We say that $F^{(τ)}$ and $F^{(τ + 1)}$ are $φ^{(τ)}$ -compatible if both sets satisfy the corresponding tomographic constraints and $φ^{(τ)} (F^{(τ)}) = F^{(τ + 1)} .$

For the following problem we suppose that a for any $τ \in [t - 1]$ a displacement field $(p, φ^{(τ)} (p) - p),$ $p \in Q^{d},$ is known. While for our main $N P$ -hardness result, the displacement field is given explicitly, all that our model of computation actually requires is that for any $τ \in [t - 1]$ and $p \in G^{(τ)},$ either a point $q \in G^{(τ + 1)}$ is provided such that $φ^{(τ)} (p) = q$ or it is reported that no such point exists. (It is not difficult to model this formally again as some displacement field oracle). With this understanding, the tomographic point tracking problem with given displacement field by the oracle $O$ is as follows.

TomDisplaceTrac $(O; S_{1}, \dots, S_{m})$ .

$t \in N,$ X-ray data functions $f_{1}^{(τ)}, \dots, f_{m}^{(τ)}$ for $τ \in [t]$ with $∥ f_{1}^{(τ)} ∥_{(1)} = \dots = ∥ f_{m}^{(τ)} ∥_{(1)} = n$ ; displacement field $φ^{(τ)}$ for $τ \in [t - 1]$ .
Find sets $F^{(1)}, \dots, F^{(t)} \in F^{d}$ such that $X_{S_{i}} F^{(τ)} = f^{(τ)}$ for all $i \in [m],$ and such that $F^{(τ)}$ and $F^{(τ + 1)},$ are $φ^{(τ)}$ -compatible for all $τ \in [t - 1],$ or decide that no such sets exist.

2.4. Main results

Our main results are as follows. First, we consider TomTrac $(O; S_{1}, \dots, S_{m})$ for Markov-type models, i.e., for objective function oracles $O$ that can be given explicitly as the sum of all costs of assigning points between consecutive moments in time. We show that the corresponding version of TomTrac $(O; S_{1}, \dots, S_{m})$ for the positionally determined case, called Trac $(d),$ can be solved in polynomial time (Theorem 1).

For the tomographic case, we show that even if all instances are restricted to $t = 2$ , and the solution $F^{(1)}$ is given explicitly, TomTrac $(S_{1}, S_{2})$ is $N P$ -hard (Theorem 2). (Also the corresponding uniqueness problem is $N P$ -hard and the counting problem is $# P$ -hard.) This result is complemented by Theorem 3, which shows that TomDisplaceTrac $(O; S_{1}, S_{2})$ is $N P$ -hard for $t = 2$ and certain conditions imposed on the specified displacement field.

Theorem 4 shows that the problem Trac $(O; d)$ is $N P$ -hard, even if all instances are restricted to a fixed $t \geq 3,$ and $O$ is a path value oracle. Also in the case of straight line movement and fixed $d \geq 2,$ $t \geq 3$ it is $N P$ -complete problem to decide whether a solution of Trac $(O; d)$ exists where all particles move along straight lines (Corollary 3).

We conclude Section 3 by introducing and discussing algorithmic paradigms that allow to include prior knowledge of the possible particle tracks.

Section 4 then studies combinatorial models. In particular, we discuss Tomography under Window Constraints, which, under rather mild conditions, turns out to be $N P$ -hard (Theorem 8). On the other hand, we show in Theorem 9 that it is in $P$ if all windows are disjoint horizontal and vertical windows of width 1. Another polynomial-time solvable variant of Tomography under Window Constraints

, which arises in super-resolution imaging applications, is discussed in Theorem

10.

3. Optimization Models

In the following we assume that we have X-ray measurements from two different directions $S_{1}, S_{2} \in S^{d}$ at $t$ moments in time. As before, $G^{(τ)}$ will denote the corresponding grid at $τ \in [t]$ .

While the grids $G^{(τ)}$ are subsets of $R^{d}$ , we need to distinguish a point $g^{(τ_{1})} \in G^{(τ_{1})}$ from a point $g^{(τ_{2})} \in G^{(τ_{2})}$ when $τ_{1} \neq τ_{2}$ even if, physically, both points occupy exactly the same position in $R^{d}$ . So, formally, $g^{(τ)}$ must be regarded as a point $(g, τ) \in R^{d} \times [t]$ . In the following, however, we will not always ‘verbally’ distinguish between the interpretations $G^{(τ)} \subseteq R^{d}$ and $G^{(τ)} \subseteq R^{d} \times [t]$ if there is no risk of confusion.

In the positionally determined case where the correct sets $F^{(τ)}$ are known, we can, of course, directly work with $F^{(τ)} .$ Hence, whenever we speak of the positionally determined case, it is in the following tacitly assumed that $G^{(τ)} = F^{(τ)} .$

3.1. A Markov-type integer programming model

In this section we will consider the problem TomTrac $(O; S_{1}, \dots, S_{m})$ for objective function oracles $O$ that can be given explicitly since their values are just the sums of all costs of assigning points between consecutive moments in time. We call these models Markov-type since the objective function reflects only dependencies that occur between neighboring layers. In Section 3.5 we will discuss more general models.

Let us now present an integer linear programming (Ilp) model for this problem. In order to make the notation as transparent as possible we number the points of each grid $G^{(τ)}$ as $g_{i}^{(τ)}$ for $i \in I^{(τ)} = [| G^{(τ)} |]$ . Also, we set $e_{i, j}^{(τ)} = (g_{i}^{(τ)}, g_{j}^{(τ + 1)})$ , and refer to such pairs as tracking edges. The objective function will then depend only on weights $ω_{i, j}^{(τ)},$ associated with the tracking edges, which measure the cost of assigning a point $g_{i}^{(τ)} \in G^{(τ)}$ to a point $g_{j}^{(τ + 1)} \in G^{(τ + 1)}$ . Explicitly expressed, we have

c (P_{1}, \dots, P_{n}) = n \sum k = 1 \sum τ \in [t] \sum e_{i, j}^{(τ)} \in P_{k} ω_{i, j}^{(τ)} .

Setting $w (e) = ω_{i, j}^{(τ)}$ for $e = e_{i, j}^{(τ)}$ , we will refer to $w$ as the weight function realizing the Markov-type objective function oracle $O$ . Here, of course, the full list of all weights $w (e)$ is part of the input of our problem. Note that the number of entries in this list is bounded from above by $(t - 1) n^{4}$ even in the totally tomographic case. Hence there are at most polynomially many rational numbers to be presented. With this specification we refer to our respective problems as TomTrac $(S_{1}, \dots, S_{m})$ and Trac $(d) .$

Note that the weights $w (e)$ can be used to model additional ‘local’ physical knowledge such as information about the potential ranges of velocities or directions of the particles.

Now, we introduce two sets of $0$ - $1$ -variables associated with the grid points and the tracking edges, respectively. More precisely, for $τ \in [t]$ and $i \in I^{(τ)},$ the variable $ξ_{i}^{(τ)}$ corresponds to $g_{i}^{(τ)}$ while the variable $η_{i, j}^{(τ)}$ is associated with $e_{i, j}^{(τ)} .$

The tomographic variables $ξ_{i}^{(τ)}$ will be used to describe the tomographic constraints; $ξ_{i}^{(τ)} = 1$ signifies that the grid point $g_{i}^{(τ)}$ is actually present in the computed solution $F^{(τ)} .$ The tracking variable $η_{i, j}^{(τ)}$ has value $1$ if, and only if, the particle that, at time $τ$ is located at $g_{i}^{(τ)},$ moves to $g_{j}^{(τ + 1)} .$ For a compact notation, we set $x^{(τ)} = (ξ_{1}^{(τ)}, \dots, ξ_{| G^{(τ)} |}^{(τ)})^{T},$ write $A^{(τ)}$ for the corresponding coefficient matrix, and encode the X-ray information in the ‘right-hand side vector’ $b^{(τ)} .$

With this notation we can in principle solve the problem as an Ilp.

Algorithm 1 (Tomographic Tracking-Ilp).

Solve the following integer linear program:

\omit\span\omit\span\omit\span\omitmin∑τ∈[t−1]∑i∈I(τ)∑j∈I(τ+1)ω(τ)i,jη(τ)i,js.\,t. A(τ)x(τ)=b(τ),(τ∈[t]),ξ(τ)i−∑j∈I(τ+1)η(τ)i,j=0,(i∈I(τ),τ∈[t−1]),ξ(τ+1)j−∑i∈I(τ)η(τ)i,j=0,(j∈I(τ+1),τ∈[t−1]),η(τ)i,j∈{0,1},(i∈I(τ),j∈I(τ+1),τ∈[t−1]),x(τ)∈{0,1}|G(τ)|,(τ∈[t]).

(3.1)

Note that the coefficient matrices $A^{(τ)}$ are totally unimodular (see, e.g., [2, 31]) and the first set of equality constraints contains only the tomographic variables. The second and third set of equalities in (3.1) couple the tomographic and the tracking variables. As the tomographic variables are $0$ - $1$ they guarantee two properties: (1) A tracking edge can only connect grid points that are present in the tomographic solution. The second set of constraints corresponds to the edges ‘leaving’ time $τ$ while the third set corresponds to those ‘entering’ time $τ + 1$ . (2) From each point $g_{i}^{(τ)}$ in a considered solution, i.e., when $ξ_{i}^{(τ)} = 1$ , exactly one ‘leaving’ tracking edge is selected. Similarly, when $ξ_{i}^{(τ + 1)} = 1$ , there must be exactly one ‘entering’ tracking edge.

Note that by adding the two constraints

ξ_{i}^{(τ + 1)} - \sum j \in I^{(τ + 1)} η_{j, l}^{(τ + 1)} = 0 % and - ξ_{i}^{(τ + 1)} + \sum j \in I^{(τ)} η_{j, i}^{(τ)} = 0,

for $i \in I^{(τ)}$ and $τ \in {2, \dots, t - 1},$ we derive the set of constraints

\sum j \in I^{(τ)} η_{j, i}^{(τ)} - \sum j \in I^{(τ + 1)} η_{i, j}^{(τ + 1)} = 0,

(3.2)

which express that, if a particle enters $g_{i}^{(τ)}$ it must exit it again, and vice versa. Hence, the equations (3.2) can be viewed as path or flow constraints, and the corresponding coefficient matrix is again totally unimodular. When we replace appropriate conditions in (3.1) by these flow constraints, we obtain an equivalent but different Ilp-formulation which decomposes into two totally unimodular parts and a reduced set of coupling constraints which, however, still destroy the total unimodularity of the whole system. As it will turn out in the next subsection, this is not merely a nuisance but a severe obstacle for efficient algorithms.

We prefer the formulation (3.1) because it reveals the uncoupled structure in the positionally determined case and leads to the following result for Trac $(d)$ . (Recall that now $O$ is a Markov-type objective function oracle which is realized by a rational weight function $w$ .)

Theorem 1.

Trac $(d)$ decomposes into uncoupled minimum weight perfect bipartite matching problems and can hence be solved in polynomial time.

Proof.

Under the assumptions of this theorem, (3.1) reads

\omit\span\omit\span\omit\span\omitmin∑τ∈[t−1]∑i∈I(τ)∑j∈I(τ+1)ω(τ)i,jη(τ)i,js.\,t. ∑j∈I(τ+1)η(τ)i,j=1,(i∈I(τ),τ∈[t−1]),∑i∈I(τ)η(τ)i,j=1,(j∈I(τ+1),τ∈[t−1]),η(τ)i,j∈{0,1},(i∈I(τ),j∈I(τ+1),τ∈[t−1]).

(3.3)

Note that, here, in the positionally determined case, we have $| I^{(τ)} | = | G^{(τ)} | = | F^{(τ)} | = n$ . Hence we need to solve the $t - 1$ independent minimum weight perfect bipartite matching tasks

\omit\span\omit\span\omit\span\omitmin∑i,j∈[n]ω(τ)i,jη(τ)i,js.\,t. ∑j∈[n]η(τ)i,j=1,(i∈[n]),∑i∈[n]η(τ)i,j=1,(j∈[n]),η(τ)i,j≥0,(i,j∈[n]).

(3.4)

In (3.4) the condition $η_{i, j}^{(τ)} \in {0, 1}$ has been replaced by the non-negativity constraints $η_{i, j}^{(τ)} \geq 0$ . Since, by the other constraints, $η_{i, j}^{(τ)} \leq 1$ , we have just switched to the Lp-relaxation. The feasible points of this Lp-relaxation form a polytope with integral vertices [12] (see also [31]), hence each of these problems can be solved in polynomial time; see, e.g., [47, Sect. 16 and 19]. ∎

Let us close this section by turning again to the ILP (3.1). Clearly, the ILP contains $O (t n^{2})$ constraints and $O (t n^{4})$ binary variables, which, in the positionally determined case, reduces to $O (t n)$ constraints and $O (t n^{2})$ binary variables. As the computation time depends strongly on the structure of the problem instance, it is not clear which computation times will occur for a given practical instance. For the positionally determined case, however, the computational study in [48] shows that random instances with $n$ up to $20, 000$ and $t = 2$ can be solved with state-of-the-art algorithms in reasonable time. Note, however, that without resorting to sparsity and advanced optimization techniques it seems hopeless to solve instances in the tomographic case with $n = t = 30$ as the coefficient matrix can involve $54, 000 \times 23, 517, 000$ entries (amounting to more than 1 Terabyte of storage space). Of course, it is always possible to resort to LP relaxations of (3.1), see [24]. In general, however, the returned solutions will then not be binary.

3.2. On the complexity of the partially tomographic case

We consider the question of when, for $m = 2$ , the tomographic case can be solved efficiently. The following result shows the limitations already for the following quite restricted partially tomographic case. There is only one time step, i.e., $t = 2$ , and $F^{(1)}$ is known while the set $F^{(2)}$ of particle positions for $τ = 2$ is only accessible through its two X-rays $X_{S_{1}} F^{(2)}$ and $X_{S_{2}} F^{(2)}$ .

Theorem 2.

Even if all instances are restricted to the case $t = 2$ , where the solution $F^{(1)}$ is given explicitly, TomTrac $(S_{1}, S_{2})$ is $N P$ -hard. Also the corresponding uniqueness problem is $N P$ -complete and the counting problem is $# P$ -complete.

Proof.

We use a reduction from the following problem, where $S_{1}, S_{2}, S_{3}$ are three different lattice lines.

Consistency $_{F^{d}} (S_{1}, S_{2}, S_{3})$ .

Data functions $f_{1},$ $f_{2},$ $f_{3} .$
Does there exist a set $F \in F^{d}$ such that $X_{S_{i}} F = f_{i}$ for all $i \in [3]$ ?

This problem and its uniqueness and counting versions have been shown in [29] to be $N P$ -complete or $# P$ -complete, respectively. In fact, the proof of [29] reveals, that Consistency $_{F^{d}} (S_{1}, S_{2}, S_{3})$ remains $N P$ -complete even if all instances are restricted to those, where, for two directions, say $S_{2}, S_{3}$ , the nonzero X-rays are all $1$ and for each of these X-ray lines the grid $G$ contains exactly two candidate points. So, let $(f_{1}, f_{2}, f_{3})$ be the data functions of such an instance, and let $n = ∥ f_{1} ∥_{(1)} = ∥ f_{2} ∥_{(1)} = ∥ f_{3} ∥_{(1)}$ .

Let $F^{(1)}$ consist of $n$ points $p_{1}, \dots, p_{n}$ on a line parallel to $S_{1}$ . Hence the support of $X_{S_{1}} F^{(1)}$ is a single line, and the corresponding value is $n,$ while the support of $X_{S_{2}} F^{(1)}$ consists of $n$ lines, and the corresponding values are $1$ . Note that $F^{(1)}$ is uniquely determined by its X-rays $X_{S_{1}} F^{(1)}$ and $X_{S_{2}} F^{(1)}$ .

Now, we consider the given (restricted) instance $(f_{1}, f_{2}, f_{3})$ of Consistency $_{F^{d}} (S_{1}, S_{2}, S_{3})$ as ‘living’ at time $τ = 2$ , regard the first two data functions as the input $(f_{1}^{(2)}, f_{2}^{(2)})$ at $τ = 2$ of our dynamic 2-direction X-ray problem, and introduce an objective function that allows to model the X-ray information in the third direction via a minimum weight bipartite matching of cardinality $n .$

So, let $G^{(2)}$ be the grid coming from the data function $f_{1}^{(2)} = f_{1}$ and $f_{2}^{(2)} = f_{2} .$ Further, let $T_{1}, \dots, T_{n}$ be the lines parallel to $S_{3}$ that meet the candidate grid $G$ of $(f_{1}, f_{2}, f_{3}),$ and let $g_{i, j}$ , $i = 1, 2$ , denote the corresponding two grid points on $T_{j}$ , $j \in [n]$ . Finally, for $j \in [n]$ and $g \in G^{(2)},$ we define the objective function vector $w = (ω_{p_{j}, g})$ by setting

ω_{p_{j}, g} = {\begin{matrix} 0 & for g \in {g_{1, j}, g_{2, j}}; 1 & otherwise. \end{matrix}

Now consider a minimum weight bipartite matching of cardinality $n$ on the complete bipartite graph with the partition $(F^{(1)}, G^{(2)})$ of the vertex set $F^{(1)} \cup G^{(2)}$ and edge costs $ω_{p_{j}, g}$ for any edge ${p_{j}, g} \in F^{(1)} \times G^{(2)} .$ In this matching let $F^{(2)}$ denote the endpoints of the tracking edges that are not in $F^{(1)} .$ Then, of course, each point $p_{j} \in F^{(1)}$ is assigned to exactly one point of $G^{(2)} .$ If the objective function value of the given matching is $0$ then $p_{j}$ must be assigned to $g_{1, j}$ or $g_{2, j}$ for each $j \in [n] .$ Hence $X_{S_{3}} F^{(2)} = f_{3} .$

On the other hand, if $F^{(2)}$ satisfies $X_{S_{i}} F^{(2)} = f_{i},$ for $i \in [3],$ then, in particular, for each $j \in [n]$ , the set $F^{(2)}$ contains exactly one of the points $g_{1, j}$ or $g_{2, j}$ . Hence there is a matching of cardinality $n$ with objective function value $0$ .

Thus, the given instance $(f_{1}, f_{2}, f_{3})$ is feasible if, and only if, there exist a set $F^{(2)}$ satisfying the X-ray constraints with respect to $S_{1}$ and $S_{2}$ that allows a perfect matching between $F^{(1)}$ and $F^{(2)}$ of weight $0.$ Of course, the transformation runs in polynomial-time and is parsimonious. ∎

The following corollary is merely a reformulation of Theorem 2 highlighting its interpretation in terms of matchings.

Corollary 1.

Min weight max cardinality bipartite Matching under totally unimodular constraints for one set of the bipartition is $N P$ -hard. Also the corresponding uniqueness problem is $N P$ -complete while the counting problem is $# P$ -complete.

It is clear that Theorem 2 can also be interpreted as a non-approximability result. In fact, since, in the case of feasibility, the optimal objective function value is $0$ , the $N P$ -hardness of the problem trivially means that, unless $P = N P$ , there is no polynomial-time approximation algorithm of relative accuracy at most $α$ for any factor $α$ . Note, however, that the definition of the objective function in the proof of Theorem 2 can be altered in order to guarantee that the objective function is bounded away from $0.$ Simply set $ω_{p_{j}, g} = 1$ for $j \in [n]$ and $g \in {g_{1, j}, g_{2, j}},$ and let all other values be sufficiently large. This fact can be used to interpret the previous result in the context of approximation complexity even if we assume that the optimum is bounded away from $0$ from below.

Corollary 2.

Unless $P = N P$ , there is no polynomial-time algorithm that solves the following problem up to any polynomial-space multiplicative constant in polynomial-time. Given $F^{(1)}$ , $f_{1}^{(2)},$ $f_{2}^{(2)},$ let $G^{(2)}$ be the grid associated with $f_{1}^{(2)},$ $f_{2}^{(2)},$ and let $w \geq 1$ encode the edge weights of the complete bipartite graph with bipartition $(F^{(1)}, G^{(2)}) .$ Find a cardinality $| F^{(1)} |$ matching of $F^{(1)}$ and a set $F^{(2)} \subseteq G^{(2)}$ with $X_{S_{1}} F^{(2)} = f_{1}^{(2)},$ $X_{S_{2}} F^{(2)} = f_{2}^{(2)},$ that is minimal with respect to $w$ for all such sets $F^{(2)} .$

In view of the comments after Corollary 1, this result may seem somewhat artificial. It does, however, allow for non approximability results that are rather close to the observed practical difficulty. In fact, the objective function in the proof of Theorem 2, and hence for obtaining the result of Corollary 2, can be adapted so as to accommodate different experimental settings. This can be combined with different choices of $F^{(1)}$ . For instance, if we choose the points $p_{j} \in F^{(1)}$ in such a way that they have equal Euclidean distance to $g_{1, j}$ and $g_{2, j}$ but different distance to all other grid points, then the objective function can be interpreted as rewarding conformity with certain known velocities of the particles. Hence the hardness proof does to a certain extent reflect practical difficulties in tracking the particles.

We remark that, in contrast to Corollary 2, there are several approximability results for the positionally determined case and specific types of cost functions, see, e.g., [10, 22, 25, 39].

3.3. On the complexity of the tomographic case when the displacement field is known

Let us now turn to the other extreme, assuming that the $F^{(1)}, \dots, F^{(t)}$ are unknown, but the particle displacement field is given for all $τ \in [t] .$ Note that this is quite different from the setting considered in the proof of Theorem 2. In fact, in the proof of Theorem 2 we make explicit use of the fact that there are two possible next positions for each particle, which renders the problem $N P$ -hard. If, on the other hand, $F^{(1)}$ and the displacement field are known, then $F^{(2)}, \dots, F^{(t)}$ and the particle tracks can be determined in polynomial time. In general, however, we will show that the reconstruction problem with a given displacement field is again $N P$ -hard.

Particle displacement fields are often (approximatively) known in experimentally controlled environments (for instance, when charged particles in a particle accelerator or plasma fusion device move in the presence of experimentally controlled electromagnetic fields; see, e.g., [13, 53]). The particle tracking task reduces in this case to that of reconstructing the particle positions from the measurements. Of course, the general aim is to reduce ambiguities by taking several time steps into account. For instance, in [24] this problem is considered for $t = 2$ in the context of compressed motion sensing. Our next theorem shows, however, that, unless $P = N P,$ there is generally no algorithm that can make efficient use of the additional data.

In order to state the theorem we introduce the following notion. Given $S_{1}, S_{2} \in S^{d},$ we call a displacement field $(p, φ (p) - p),$ $p \in Q^{d},$ proper (for $S_{1}$ and $S_{2}$ ) if $φ$ is an affine transformation with $φ (Q^{d}) \subseteq Q^{d}$ and $φ (S_{i}) - φ (0) \notin {S_{1}, S_{2}},$ $i = 1, 2.$ Note that, in particular, this notion of properness encompasses several different types of circular displacements; for instance, $φ$ might represent a $45^{\circ}$ -rotation in the $(S_{1}, S_{2})$ -plane that is followed by a factor $2 / \sqrt{2}$ dilation.

Theorem 3.

Let $φ$ denote the affine transformation (encoded by a rational matrix and a rational translation vector) of a proper displacement field. Then, the problem TomDisplaceTrac $(O; S_{1}, S_{2})$ is $N P$ -hard, even if all instances are restricted to the case $t = 2$ . Also the corresponding uniqueness problem is $N P$ -hard and the counting problem is $# P$ -hard.

Proof.

We construct a parsimonious transformation from Consistency $_{F^{d}} (S_{1}, S_{2}, S_{3}, S_{4})$ for certain lattice lines $S_{1}, S_{2}, S_{3}, S_{4}$ . This problem is the analogue to Consistency $_{F^{d}} (S_{1}, S_{2}, S_{3})$ and involves four given data functions $f_{1},$ $f_{2},$ $f_{3},$ $f_{4} .$ The task is to decide whether there exist a set $F \in F^{d}$ such that $X_{S_{i}} F = f_{i}$ for all $i \in [4] .$

Again, this problem and its uniqueness and counting versions have been shown in [29] to be $N P$ -complete or $# P$ -complete, respectively for any four different different lattice lines $S_{1}, S_{2}, S_{3}, S_{4}$ .

Now, let $S_{1}, S_{2} \in S^{d}$ with $S_{1} \neq S_{2},$ , and set $S_{3} = φ (S_{1}) - φ (0)$ and $S_{4} = φ (S_{2}) - φ (0)$ . Since the given displacement field is proper, $S_{1}, S_{2}, S_{3}, S_{4}$ are four different lattice lines. Hence, the hardness results of [29] do apply to Consistency $_{F^{d}} (S_{1}, S_{2}, S_{3}, S_{4})$ . So, suppose, we are given an instance of this problem by means of the X-ray functions $(f_{1}, f_{2}, f_{3}, f_{4})$ , and let $G$ be the corresponding candidate grid.

By associating a binary variable with each point $g_{1}, \dots, g_{| G |}$ of $G$ we can formulate the task as deciding feasibility of the ILP

\begin{matrix} A x & = & b A^{'} x & = & b^{'} x & \in & {0, 1}^{| G |}, \end{matrix}

where $A$ and $A^{'}$ encode the point-line incidences along ${S_{1}, S_{2}}$ and ${S_{3}, S_{4}},$ respectively, while $b$ and $b^{'}$ encode the data functions $f_{1}, f_{2},$ and $f_{3}, f_{4},$ respectively.

Setting $A^{(1)} = A,$ $A^{(2)} = A^{'},$ $b^{(1)} = b,$ $b^{(2)} = b^{'},$ this problem is, of course, equivalent to the ILP

\begin{matrix} A^{(1)} x^{(1)} & = & b^{(1)}, A^{(2)} x^{(2)} & = & b^{(2)}, x^{(1)} - x^{(2)} & = & 0, x^{(1)} & \in & {0, 1}^{| G |}, x^{(2)} & \in & {0, 1}^{| G |} . \end{matrix}

(3.5)

Now, suppose we associate with each variable $ξ_{i}^{(2)},$ $i \in [| G |],$ the point $φ^{- 1} (g_{i}) \in Q^{d} .$ This does not change the ILP, but only its geometric interpretation. The problem $A^{(2)} x^{(2)} = b^{(2)}$ can now be viewed as to involve only X-rays taken from the directions

φ^{- 1} (S_{i + 2}) - φ^{- 1} (0) = φ^{- 1} (φ (S_{i}) - φ (0)) - φ^{- 1} (0) = S_{i}, (i \in [2]) .

Consider now the ILP (3.1) for $t = 2$ and the following settings. Let $G^{(1)}$ denote the grid defined by the data functions $f_{1}, f_{2}$ (along $S_{1}, S_{2}$ ), and let $G^{(2)} = φ^{- 1} (G^{'})$ where $G^{'}$ denotes the grid defined by the data functions $f_{3}, f_{4}$ (along $S_{3}, S_{4}$ ). Let the points of each grid $G^{(τ)}$ be labeled $g_{i}^{(τ)}$ for $i \in I^{(τ)} = [| G^{(τ)} |],$ $τ \in [t],$ and set

Dynamic Discrete Tomography

Authors

On the Reconstruction of Static and Dynamic Discrete Structures

A mathematical approach towards THz tomography for non-destructive imaging

Atomic Super-Resolution Tomography

Fast 3D Image Moments

The Modulo Radon Transform: Theory, Algorithms and Applications

Rejection sampling from shape-constrained distributions in sublinear time

Spatiotemporal tomography based on scattered multiangular signals and its application for resolving evolving clouds using moving platforms

1. Introduction

2. Notation, basics and main results

2.1. Discrete Tomography

2.2. Dynamic Discrete Tomography

2.3. Basic algorithmic problems

2.4. Main results

3. Optimization Models

3.1. A Markov-type integer programming model

Algorithm 1 (Tomographic Tracking-Ilp).

Theorem 1.

Proof.

3.2. On the complexity of the partially tomographic case

Theorem 2.

Proof.

Corollary 1.

Corollary 2.

3.3. On the complexity of the tomographic case when the displacement field is known

Theorem 3.

Proof.

Dynamic Discrete Tomography

Authors

Related Research

On the Reconstruction of Static and Dynamic Discrete Structures

A mathematical approach towards THz tomography for non-destructive imaging

Atomic Super-Resolution Tomography

Fast 3D Image Moments

The Modulo Radon Transform: Theory, Algorithms and Applications

Rejection sampling from shape-constrained distributions in sublinear time

Spatiotemporal tomography based on scattered multiangular signals and its application for resolving evolving clouds using moving platforms

This week in AI

1. Introduction

2. Notation, basics and main results

2.1. Discrete Tomography

2.2. Dynamic Discrete Tomography

2.3. Basic algorithmic problems

2.4. Main results

3. Optimization Models

3.1. A Markov-type integer programming model

Algorithm 1 (Tomographic Tracking-Ilp).

Theorem 1.

Proof.

3.2. On the complexity of the partially tomographic case

Theorem 2.

Proof.

Corollary 1.

Corollary 2.

3.3. On the complexity of the tomographic case when the displacement field is known

Theorem 3.

Proof.