1. Introduction
In the following, the goal is to determine the paths of particles in space over a period of moments in time from images taken by a a fixed number 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 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 moments in time, a solution of the discrete tomography task at time is explicitly available (and is then considered the correct solution regardless whether it is uniquely determined by its X-rays). For , two (affine) lines in 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 -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 -hard even for ; 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 -hard already for 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 -hard for 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 ); 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 -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 -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 -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 -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 and denote the set of reals, rational numbers, integers, and natural numbers, respectively. For , let and . With
we denote the all-ones vector. Further, let
, , and setThe point sets model the sets of particles (in a static
environment). Let us point out that, mathematically, the restriction to
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
With we denote the set of all lattice lines, i.e., the set of -dimensional linear subspaces of that are spanned by vectors from . (Note that as coincides with the set of lines spanned by rational vectors, there is no further restriction here.) The lines in 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 we set Then, for and the (discrete -dimensional) X-ray of parallel to is the function
defined by
for each where, as usual, denotes the cardinality of and
is the characteristic function of
.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 are called tomographically equivalent with respect to if
for all
Given different lines , the X-ray data is given in terms of functions
with finite support , represented by appropriately chosen data structures; see [29].
Suppose, we are given an instance of measurements of an otherwise unknown set Then, of course, , where is the usual -norm. From the data functions we can infer that must be contained in the grid
of candidate points. Note that, in general, even when the given instance has a unique solution the grid can be a proper superset of
2.2. Dynamic Discrete Tomography
In tomographic point tracking, we consider consecutive moments in time. Hence, in particular, we are interested in the sets of particles at each moment in time . Since, in general, these sets are accessible only through their X-rays taken from a given finite set of directions, they represent uncoupled instances of discrete tomography.
Of course, we also need to track the particles over time. Let denote the (abstract) set of particles. Then, for each , we are actually interested in a one-to-one mapping that identifies the points of with the particles. Hence, the particle tracks are given by 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 already the static reconstruction problem of discrete tomography is -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 . In real-world particle tracking this corresponds to images taken from two cameras which may be regarded as two planes in orthogonal to two different directions , , respectively.
For a set satisfying the X-ray constraints with respect to the two given directions and can be efficiently determined for each ; see [27, 46]. Also, it can be checked efficiently, whether each 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 the correct set we still have to identify the paths of all individual particles over time. This turns the 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 denote the set of all -tuples of potential particle tracks over moments of time for fixed point sets . 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, is then still a lower bound for the number of different potential solutions. Now note that , and even if the order of the individual paths is irrelevant, this number reduces only to which is still exponential in and .
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’ of all explicitly. Even if we assume that all sets are given and such costs can be computed as a simple function of the costs of the individual paths, e.g., as
the number of different costs of a particle path (which still need to be available to determine the cost of a solution) reduces only to
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 its cost 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 will be called an objective function oracle. In the special case that provides the values and 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 given different directions and based on an objective function oracle , can now be formulated as follows.
TomTrac.
-
X-ray data functions for with .
-
Decide whether, for each , there exists a set such that for all . If so, find particle tracks of minimal cost for among all couplings of all tomographic solutions .
Let us add a remark concerning the assumption, that the -norms of the X-ray data functions in TomTrac 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 particles are recorded at each moment in time
In the following we will distinguish the cases that for none, some or all , the correct solution 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 reduces to the following problem.
Trac.
-
, sets with .
-
Find particle tracks of minimal cost for among all couplings of the sets .
Note that each instance of Trac can be viewed as a -dimensional assignment problem. For a comprehensive survey on assignment problems see [16] and the references quoted therein.
As pointed out before, for , the positionally determined case is the generic situation even for as two lines parallel to and do only meet in points of . However, for , any two non-parallel lines in the plane meet in a point, whence forming a candidate grid of size . (Note that, generally, no fixed number of X-ray images suffices for unique determination.) Also, even in with there are particle constellations for which the candidate grids do not consist of just 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 ) is a pair with denoting a map where denotes the displacement vector for a particle located at position at time In particular, if denotes the position of the -th particle at time then for and Now, the tomographic tasks for the time step consists of determining the positions of the particles at time such that the solution at time determined by the displacement field matches the tomographic data for We say that and are -compatible if both sets satisfy the corresponding tomographic constraints and
For the following problem we suppose that a for any a displacement field is known. While for our main -hardness result, the displacement field is given explicitly, all that our model of computation actually requires is that for any and either a point is provided such that 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 is as follows.
TomDisplaceTrac.
-
X-ray data functions for with ; displacement field for .
-
Find sets such that for all and such that and are -compatible for all or decide that no such sets exist.
2.4. Main results
Our main results are as follows. First, we consider TomTrac for Markov-type models, i.e., for objective function oracles 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 for the positionally determined case, called Trac can be solved in polynomial time (Theorem 1).
For the tomographic case, we show that even if all instances are restricted to , and the solution is given explicitly, TomTrac is -hard (Theorem 2). (Also the corresponding uniqueness problem is -hard and the counting problem is -hard.) This result is complemented by Theorem 3, which shows that TomDisplaceTrac is -hard for and certain conditions imposed on the specified displacement field.
Theorem 4 shows that the problem Trac is -hard, even if all instances are restricted to a fixed and is a path value oracle. Also in the case of straight line movement and fixed it is -complete problem to decide whether a solution of Trac 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 -hard (Theorem 8). On the other hand, we show in Theorem 9 that it is in 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 at moments in time. As before, will denote the corresponding grid at .
While the grids are subsets of , we need to distinguish a point from a point when even if, physically, both points occupy exactly the same position in . So, formally, must be regarded as a point . In the following, however, we will not always ‘verbally’ distinguish between the interpretations and if there is no risk of confusion.
In the positionally determined case where the correct sets are known, we can, of course, directly work with Hence, whenever we speak of the positionally determined case, it is in the following tacitly assumed that
3.1. A Markov-type integer programming model
In this section we will consider the problem TomTrac for objective function oracles 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 as for . Also, we set , and refer to such pairs as tracking edges. The objective function will then depend only on weights associated with the tracking edges, which measure the cost of assigning a point to a point . Explicitly expressed, we have
Setting for , we will refer to as the weight function realizing the Markov-type objective function oracle . Here, of course, the full list of all weights is part of the input of our problem. Note that the number of entries in this list is bounded from above by 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 and Trac
Note that the weights 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 --variables associated with the grid points and the tracking edges, respectively. More precisely, for and the variable corresponds to while the variable is associated with
The tomographic variables will be used to describe the tomographic constraints; signifies that the grid point is actually present in the computed solution The tracking variable has value if, and only if, the particle that, at time is located at moves to For a compact notation, we set write for the corresponding coefficient matrix, and encode the X-ray information in the ‘right-hand side vector’
With this notation we can in principle solve the problem as an Ilp.
Algorithm 1 (Tomographic Tracking-Ilp).
Solve the following integer linear program:
(3.1) |
Note that the coefficient matrices 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 - 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 . (2) From each point in a considered solution, i.e., when , exactly one ‘leaving’ tracking edge is selected. Similarly, when , there must be exactly one ‘entering’ tracking edge.
Note that by adding the two constraints
for and we derive the set of constraints
(3.2) |
which express that, if a particle enters 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. (Recall that now is a Markov-type objective function oracle which is realized by a rational weight function .)
Theorem 1.
Trac 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
(3.3) |
Note that, here, in the positionally determined case, we have . Hence we need to solve the independent minimum weight perfect bipartite matching tasks
(3.4) |
In (3.4) the condition has been replaced by the non-negativity constraints . Since, by the other constraints, , 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 constraints and binary variables, which, in the positionally determined case, reduces to constraints and 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 up to and 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 as the coefficient matrix can involve 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 , 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., , and is known while the set of particle positions for is only accessible through its two X-rays and .
Theorem 2.
Even if all instances are restricted to the case , where the solution is given explicitly, TomTrac is -hard. Also the corresponding uniqueness problem is -complete and the counting problem is -complete.
Proof.
We use a reduction from the following problem, where are three different lattice lines.
Consistency.
-
Data functions
-
Does there exist a set such that for all ?
This problem and its uniqueness and counting versions have been shown in [29] to be -complete or -complete, respectively. In fact, the proof of [29] reveals, that Consistency remains -complete even if all instances are restricted to those, where, for two directions, say , the nonzero X-rays are all and for each of these X-ray lines the grid contains exactly two candidate points. So, let be the data functions of such an instance, and let .
Let consist of points on a line parallel to . Hence the support of is a single line, and the corresponding value is while the support of consists of lines, and the corresponding values are . Note that is uniquely determined by its X-rays and .
Now, we consider the given (restricted) instance of Consistency as ‘living’ at time , regard the first two data functions as the input at 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
So, let be the grid coming from the data function and Further, let be the lines parallel to that meet the candidate grid of and let , , denote the corresponding two grid points on , . Finally, for and we define the objective function vector by setting
Now consider a minimum weight bipartite matching of cardinality on the complete bipartite graph with the partition of the vertex set and edge costs for any edge In this matching let denote the endpoints of the tracking edges that are not in Then, of course, each point is assigned to exactly one point of If the objective function value of the given matching is then must be assigned to or for each Hence
On the other hand, if satisfies for then, in particular, for each , the set contains exactly one of the points or . Hence there is a matching of cardinality with objective function value .
Thus, the given instance is feasible if, and only if, there exist a set satisfying the X-ray constraints with respect to and that allows a perfect matching between and of weight 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 -hard. Also the corresponding uniqueness problem is -complete while the counting problem is -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 , the -hardness of the problem trivially means that, unless , 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 Simply set for and 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 from below.
Corollary 2.
Unless , there is no polynomial-time algorithm that solves the following problem up to any polynomial-space multiplicative constant in polynomial-time. Given , let be the grid associated with and let encode the edge weights of the complete bipartite graph with bipartition Find a cardinality matching of and a set with that is minimal with respect to for all such sets
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 . For instance, if we choose the points in such a way that they have equal Euclidean distance to and 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.
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 are unknown, but the particle displacement field is given for all 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 -hard. If, on the other hand, and the displacement field are known, then 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 -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 in the context of compressed motion sensing. Our next theorem shows, however, that, unless 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 we call a displacement field proper (for and ) if is an affine transformation with and Note that, in particular, this notion of properness encompasses several different types of circular displacements; for instance, might represent a -rotation in the -plane that is followed by a factor 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 is -hard, even if all instances are restricted to the case . Also the corresponding uniqueness problem is -hard and the counting problem is -hard.
Proof.
We construct a parsimonious transformation from Consistency for certain lattice lines . This problem is the analogue to Consistency and involves four given data functions The task is to decide whether there exist a set such that for all
Again, this problem and its uniqueness and counting versions have been shown in [29] to be -complete or -complete, respectively for any four different different lattice lines .
Now, let with , and set and . Since the given displacement field is proper, are four different lattice lines. Hence, the hardness results of [29] do apply to Consistency. So, suppose, we are given an instance of this problem by means of the X-ray functions , and let be the corresponding candidate grid.
By associating a binary variable with each point of we can formulate the task as deciding feasibility of the ILP
where and encode the point-line incidences along and respectively, while and encode the data functions and respectively.
Setting this problem is, of course, equivalent to the ILP
(3.5) |
Now, suppose we associate with each variable the point This does not change the ILP, but only its geometric interpretation. The problem can now be viewed as to involve only X-rays taken from the directions
Consider now the ILP (3.1) for and the following settings. Let denote the grid defined by the data functions (along ), and let where denotes the grid defined by the data functions (along ). Let the points of each grid be labeled for and set
Comments
There are no comments yet.