CC BY
56
Вычислительные технологии
Том 9, № 5, 2004
ON THE APPLIED MATHEMATICS OF DISCRETE TOMOGRAPHY
O. Ya§ar, C. Diner, A. Dogan Institute of Applied Mathematics, Mathematics Department,
METU, Ankara, Turkey e-mail: [email protected]
G.-W. Weber Institute of Applied Mathematics, METU, Ankara, Turkey F. Ozbudak, A. Tiefenbach Institute of Applied Mathematics, Mathematics Department, METU, Ankara, Turkey
В прикладной математике обратные проблемы реконструирования в задачах геологии, геодезии, медицины, экономики и технологии в последнее время дополнились задачами томографии. Это привело к развитию компьютерной томографии. Обычная компьютерная томография имеет дело с непрерывной пространственной переменной с ограниченным числом проекций. В дискретной томографии делается попытка реконструкции дискретного множества из проекций. В статье приводится обзор предметной области и наших разработок. Практическое применение других разделов математики также рассмотрено. Определены проблемы непротиворечивости, единственности и реконструкции. Приведены последние достижения, обозначены нерешенные проблемы и темы для будущих исследований.
Introduction
The starting point of computerized tomography (CT) might be the need to construct the density distribution within the human body by means of X-ray projections. Let us consider the problem of locating a tumor. We often need an estimate of the location on the basis of noninvasively available data to plan the treatment or an operation. In our case, the available information consists of the projections. The algorithm of CT reconstructs the volume data at a resolution limited by the number of projections. However, it is possible to reconstruct the data at much higher resolutions, if it may take a limited number of values from a discrete set. The accuracy of localization depends on the resolution, while larger number of projections cost higher doses of ionizing radiation. Thus Discrete Tomography, abbreviated as DT, advances CT whenever it is applicable. In DT, we try to solve such problems in an ideal, at least approximate way and develop algorithms.
© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.
Fig. 1. A picture illustrating brain tomography taken from Nevzat Genger [13, 23].
Taking the tomography of a 3 dimensional object, e.g., human brain, means observing the 2 dimensional slides, whereas in two dimensions the projections are obtained by rays, e.g. X-rays. In brain tomography, the system developed by Nevzat Gencer from Brain Research Laboratory and Biomedical Group in METU, uses magnetic excitation to induce currents inside the body and measure the resulting magnetic fields [13, 23] (Fig. 1). They investigate both the biology of the brain and the technical devices, and look for improving the methods used from numerics and optimization. A good reference book on inverse problems is given by Aster et al. [1].
We assume that we are given a domain which can be discrete or continuous and a function f with a discrete range. Our aim is to reconstruct f from weighted sums which are the projections of f in the chosen horizontal and vertical directions in the image (Fig. 2). This is an inverse problem in which we want to reconstruct a lattice set from its X-rays or projections [5]. Batenburg [2] presents an improvement of a reconstruction algorithm in order to minimize the time complexity where the optimized version is 50 times faster than the existing approach [16].
In this paper, we are mainly concerned with the mathematical theory of Discrete Tomography and about how we can make it useful in applications by means of computer programs or algorithms.
Let us say a few words about the history of DT. Although the first consistency result was given by H.J. Ryser in 1957 [31], the name Discrete Tomography was given at 1994 by Larry Shepp. The question of uniquely determining a planar convex object was proposed by P.C. Hammer in 1961 [18]. DT strongly developed and valuable scientific exchange took place at a series of workshops, three of which were held in Germany (1994), Hungary (1997) and France (1999). The forthcoming meeting will be held in New York in 20051. For further information on
1 http:// www .dig.cs.gc.cuny.edu/workshops / workshop. html
the theory, algorithms and applications of DT please refer to the book of Herman and Kuba [19].
When the subject to be reconstructed is assumed to be continuous, real analysis is used to develop the theory. As the name suggests, Discrete Tomography has its theory based on discrete mathematics [11, 14]. In addition, it deals with many other fields of mathematics, namely combinatorics, functional analysis, geometry, coding theory and graph theory.
In recent years, various applications of DT is reported. It has been applied to diverse areas such as medical sciences, image processing, electron microscopy, scheduling, statistical data security, game theory and material sciences. For instance considered as a first result on medical applications, Reiber et al. [29] reconstructed the right coronary artery from two cineangiograms. For various medical applications of DT such as reconstruction from sparse radiographic data, enhancement of tomographic images, reconstruction of human organs, e.g., blood vessels, please refer to the survey paper of Kuba et al. [25]. We will provide a survey of some important applications in Section 5.
After having introduced the problem here, in the rest of the paper we will give the notation, basic problems of DT and their complexity results. Consistency problem will be presented in detail. Uniqueness and additivity will be examined. We will finish the paper with some open problems and applications.
1. Preliminaries
Let us first give the basic notation and definitions which we will use in this paper. We will refer to the notation used by Shepp et al. [7].
Here Z stands for the set of integers, and N0 stands for the set of natural numbers including
0. Lattice sets are discrete sets F C Zd which are finite subsets of integers (Fig. 3), d > 2.
Lattice directions are nonzero vectors of Zd over the field Q, rationals. A finite sequence of distinct lattice directions will be denoted by D, hence,
D = (vi,...,vq), q > 2. (1)
A lattice line l is parallel to a vector vk £ D and furthermore l fl Zd is nonempty. For a visualization of lattice lines see li and l2 in Fig. 2.
The set of all lattice lines that are parallel to vk £ D is denoted by Lk and E will be the class of finite sets in Zd.
The collection of the set of lattice lines determined by D is given by
L =(Li,...,Lq), q > 2. (2)
Fig. 3. A lattice is given at the left and at the right hand side is a very simple lattice set.
1
1
2
X 2
* X 3 9 X 4 I 1
X 6 W
V. t 1 J
2 2
Fig. 4. An instance illustrating the problem as a linear programming problem. Projection on two lattice lines in the directions (1,0) and (0,1).
Hence, a lattice line l will correspond to an X-ray and it will be chosen from the set Lk if l is parallel to vk.
Our reconstruction problem can be redefined and considered as a feasibility problem in two dimensions. In this case, any two distinct lattice directions can be considered, but since a linear transformation is enough to transform any two dimensional lattice to a lattice with directions (1, 0) and (0,1), we will consider that latter case. The aim is to find the binary vector which satisfies a matrix equation, which may be given as (cf. also [1])
Px = b, (3)
where P £ {0,1}MxN, b £ . Here, if the smallest rectangular box containing the finite set
to be reconstructed has dimensions n1 x n2, then M = n1 + n2 and N = n1 ■ n2. Hence, M is the number of lattice lines, in this case parallel to the directions (1, 0) and (0,1) on which there is at least one element from out discrete set and N is the total number of points to be reconstructed.
We could as well extend our feasibility problem by maximizing f (x), the sum of the Xj’s subject to equation (3). A solution of this optimization problem implies feasibility, but it additionally aims at a maximal atom density. Of course, other objective functions f (x) are possible as well, and they can be selected depending on the application from science, technology, ecology, social science or medicine. For the optimization problem we mentioned, there are polynomial time interior point methods [33].
Example 1. Consider the lattice set given in Fig. 4. It is contained in an 3 x 2 rectangle, hence M = 5 and N = 6.
For this instance we have the following system of equations
Xi =1
x4 =1
X5 + xe = 2
xi + x5 =2
x4 + x6 = 2.
Here, P, X and b are
( 1 0 0 0 0 0 \
P
0 0 0 0 0 0
0 0 11
1 0 0 0 1 0 \ 0 0 0 0 0 1 /
x =
( Xi \
X2 X3 x4 X5 V X6 /
1 1 2 2
V2/
Projection of a lattice set in direction vk is pp^ : L(k) ^ N0 such that
pFk)(i) = |F n i| = £ f (X)
(4)
x€l
where f is the characteristic function of F.
Two lattice sets F and F' are said to be tomographically equivalent with respect to the directions D if the following equality is satisfied
(k) (k)
Pf =Pf > ,
(5)
Now, let us state the three main problems that discrete tomography is concerned with.
Consistency (E, L)
Given: For k = 1,... , q a function : L(k) ^ N0 with finite support.
Question: Does there exist an F 6 E such that p^ = pk for k = 1,... , q?
Uniqueness (E, L)
Given: An F 6 E.
Question: Does there exist a different F' 6 E such that F and F' are tomographically equivalent with respect to the directions of D?
If a set is non-unique, then it means that it can not be distinguished from any other set in Zd by its projections.
Reconstruction (E, L)
Given: For k = 1,... , q a function p^/1 : L(k) ^ N0 with finite support.
Task: Construct a finite set F 6 E such that p^ = pk for k = 1,... , q.
It should be clear from the definitions that if the reconstruction problem is solvable, then the consistency problem is solved as well.
Complexity results for lattice lines are given below [10, 22]. In the following table, PTA stands for polynomial time algorithm.
One of the results by Gritzmann and Gardner [36] states that any discrete set of Z2 is unique with respect to D if the cardinality of D exceeds 6.
The complexity of the problem is not in general developed yet for r-dimensional X-rays when r is greater than 1. Only a few results are known. For instance, when r = 2, n = 3 and L consists of the 3 coordinate planes, the problems remain open.
If we only have two lattice directions, then due to Ryser’s theorem, by checking whether a known sub-matrix is contained in a binary matrix or not, we can say whether it is unique or not. On the other hand, it is known that one can not check uniqueness with the same procedure if there are more than two directions.
b
Consistency Uniqueness Reconstruction
q = 2 ,n = 2 PTA PTA PTA
q > 3 NP Complete NP Complete NP Hard
2. Two Projections
As we learned, the consistency problem deals with the problem, given a projection function, P, whether we can find a discrete set F which satisfies P. The solution depends on the number of direction vectors and the dimension of our space. Here we are going to explain and illustrate the results given in [19].
In this section, we restrict the problem to two dimensional space and two-direction vectors,
i.e., the projection function gives the values of the lines in two directions, namely, in x-axis direction and in y-axis direction. In this case, we can model the problem by linear algebra. So, we are actually trying to find a binary m x n matrix F, which corresponds to our discrete set, whenever a projection function and two vectors R and S representing the value of row sums and column sums, are given. The 1’s in the binary matrix denote the existence of a tissue in that entry and 0’s denote nonexistence.
But the given projection function should satisfy some conditions in order to be meaningful; this property is called compatibility:
(R, S): pair of vectors is compatible if
1. R £ , S £ N£.
2. r% < n for 1 < i < m, Sj < m for 1 < j < n.
3 r. — Vn s.
°. i=i 1 % Z^j=i sj.
These conditions are actually straightforward to deduce. Since the matrix F is binary, row
and column sums should be natural numbers as a result of first condition. The second condition
implies that any row of F cannot exceed n since there are n columns; similar bound is requested for column sums. The third condition is also trivial because it does not matter how one add the numbers in the matrix.
Let (R, S) be given compatible vectors. Using R — (r1,...,rm) we define the modified column sum S — (Sr1,... , Sn) as follows:
Let A be the m x n matrix whose i-th row consists of r% 1’s, followed by n — r% 0’s for each 1 < i < m. Then S is the column sums of A.
Using S — (s1,... , sn) we define S' — (si,... , sn) as follows: S' is nothing but the permutation
of S such that
si > 4... > sn.
Example 2. Let (R, S) be given as R = (2, 3, 4,1)
4 3 2 1 0
2 1 1 0 0 0 \
3 1 1 1 0 0
4 1 1 1 1 0
1 1 0 0 0 0
This is A since in each row all 1’s are shifted to
(4, 2,1, 2,1); then,
A.
e very left of the matrix. Hence, S — (4, 3, 2,1, 0) and S' — (4, 2, 2,1,1).
By U(R, S) we denote the class of binary matrices whose vector of row sums is equal to R and whose vector of column sums is equal to S.
Now, we can state the following theorem that gives a solution for the consistency problem.
Theorem 1. Let R = (ri,..., rm) and S class U(R, S) is non-empty if and only if
(s1,... , sn) be a pair of compatible vectors. The
j > ^ Sj, for 2 < l < n. j=i
(6)
Example 3. Using (R, S), given in the previous example, we see that
s5
s5 + S4
S5 + S4 + s3
s5 + s4 + s3 + s'
2 —
s5,
s 5 + s 4 j
s 5 + s 4 + s 3,
s 5 + s 4 + s 3 + s 2-
Hence, by theorem 1, U(R, S) is non-empty.
The advantage of this theorem is that it solves the consistency problem. Now, we give an algorithm solving the reconstruction problem which also proves the sufficiency of theorem 1. This means that we can find a matrix which satisfies U(R, S) if the projection vectors, R and
S, satisfy the conditions above. Our algorithm works as follows:
Algorithm 1.
Input Step 1 Step 2 Step 3
a compatible pair of vectors (R, S) satisfying (6); construct S' from S by permutation n; let B = A and k = n; while (k > 1),
{
m
while (sk >^bifc),
¿=1
{
Step 4 Output
let jo = maxi<i<m{j < k | by = 1, bij+i = ••• = bifc = 0};
let row i0 be where such a j0 was found;
set bi0j0 = 0 and biok = 1 (i.e., shift the 1 to the right)
};
reduce k by 1 };
construct the matrix A from B by permutation n-1 of the columns; matrix A.
We note that the complexity of this algorithm is O(n(m + logn)) [19]. Furthermore, it does not give all the solutions of the problem, instead it gives only one solution which satisfies (6). Another remark is that, the determination of the precise number of matrices in U(R, S) is an open problem, for which only lower bounds are found [20].
In addition to this theorem, there is an alternative solution approach for the consistency problem. We should give a definition first:
n
Definition 1. For all index sets I C {1,... , m} and J C {1,... , n} , let
i(/,J )=| I|| J|+ £ r, - £ sj.
jeJ
In fact, t(I, J) = CT0(A[I, J]) + a1(A[J, J]), where a0(A[I, J]) is the number of 0’s in the augmented matrix A[I, J] and a1(A[J, J]) denotes the number of 1’s in the matrix (A[J, J]) if such a matrix exists. Using this function, Ford and Fulkerson gave another necessary and sufficient condition for the consistency problem:
Theorem 2. Let R = (r1,... , rm) and S = (s1,... , sn) be a pair of compatible vectors. The class U(R, S) is non-empty if and only if
t(I, J) > 0, for all I C {1,..., m} and J C {1,..., n}.
The difference between these theorems is that the first one can be checked in a polynomial time, whereas the second one can be computed in an exponential time algorithm. Because it is enough to consider n — 1 inequalities in the first theorem, whereas in the other one 2m+n inequalities should be checked. The second theorem is not only an alternative solution of the existence problem, but with the help of the function defined in the theorem, we can also solve the existence problem of a restricted class of U(R, S), which will be introduced in the next definition.
Now, we will think about what can be done if we have further information about the class U(R, S). With further information we mean the case of not only knowing the row and column sums, but also some elements of the matrix are prescribed to be 0 or 1. Now, let us consider the general problem.
Definition 2. Let P and Q be binary matrices of size m x n. We say that Q covers P if
p,j < q,j for i = 1,... , m, j = 1,... , n. We denote this relation as P < Q. We define
UQ(R,S) := {A | P < A < Q,A e U(R,S)}.
Clearly, UQ(R, S) = U(R, S) if P = (0)mxn =: 0 and Q = (1)mxn. By suitable matrices we are able to prescribe 0 or 1 to any position. Now, without loss of generality, we can restrict ourselves to the classes UQ(R, S) =: (R, S).
Theorem 3. Let R = (r1,... , rm) and S = (s1,... , sn) be a pair of compatible vectors. The
class UQ(R, S) is nonempty if and only if
EE q,j > max \ — 5] sj ,Ssj — S r,|I C {1,... , m}, J C {1,... , n}
ie/ j€J ’ [ jJ j&J
for all I C {1,... , m} and J C {1,... , n}.
As we can see easily, this theorem is the prescribed version of theorem 2. Hence, this theorem solves the consistency problem for the class UQ(R, S).
Definition 3. Let A G UQ (R,S). A switching chain is a finite sequence (ii,ji), (*i,j2),
(i2, j2), (i2, j3), • • • , (ik, jk), (ik, j\) of free positions of the matrix A such that
ailjl = ai2j2 = ''' = aikjk =
— 1 — aiij2 = 1 — ai2j3 = ' ' ' = 1 — aik-ljk =
= 1 - aik jl
(k > 2). The corresponding switching (operation) is defined as changing the 0’s and the 1’s at all positions in the switching chain.
Theorem 4. A binary matrix with prescribed values is nonunique if and only if it has a switching chain.
Now, we will give another condition for uniqueness with the help of the first algorithm. Instead of the inequality in theorem 1, we will examine the situation of equality. Suppose A G U(R, S) and
n n
sj = ^ Sj for 1 < l < n. j=1 j=1
In this case it is obvious that = S. Hence, we do not have to change anything in the third step of Algorithm 1. Therefore, S' has no switching component since a matrix in the form of A has all of its 1’s in the leftmost position of it. Then, what we get at the end of the algorithm from the inverse of the permutation n cannot have a switching component. Hence, it is unique. So, we have showed that the above equality implies uniqueness. We also have shown that, a matrix in the form of A is a unique matrix.
We will have one more observation for the end of the uniqueness discussion.
Lemma 5. If A is a unique binary matrix, then
aij = 1 Sj > |{k | rk > ri}|.
By the help of this lemma, we will give our last condition for uniqueness of binary matrices. However, we should first define another concept of ’additivity’.
Definition 4. A = (aij) is additive if there are vectors X = (x1,... ,xm) G Rm and Y = (y1,... ,yn) G Rn such that, for i = 1,... ,m and j = 1,... ,n,
aij = 1 Xi + yj > 0.
The following theorem gives our last condition for binary matrices being unique.
Theorem 6. A binary matrix is unique if and only if it is additive.
If a binary matrix is unique, taking yj = Sj and xi = |{k | rk > ri}| implies that it
is additive by the previous lemma. Also, since the sum of any two real numbers cannot be
negative and at the same time non-negative, no binary matrix can have a switching component if it is additive. This yields the above theorem.
Summarizing what we did up to here in this section of uniqueness, gives us the following theorem.
Theorem 7. Let A G U(R, S). The following conditions are equivalent:
(1) A is unique with respect to R and S;
(2) A has no switching component;
nn
(3) sj = S s, for 2 < ^ <n;
j= j=
(4) A is additive.
Now, for more detailed observations, we will classify the elements of a binary matrix into three sets.
Definition 5. Let U(R, S) be a nonempty class. The position (i,j) is variant if there exist matrices A,B G U(R, S) such that aij = 1 — bij. A position (i,j) is an invariant 0 or an
invariant 1 if aij = 0 or aij = 1 for all A G U(R, S), respectively.
Also, we will use non-increasing permutations of the row and column sums, without loss of generality.
Definition 6. A class U'(R', S') is called normalized if the elements of the vectors R' = (r1/,... , rm') and S' = (s1',... , sn') are ordered as nonincreasing. In a normalized class, we define the (m+1) x (n+1) structure matrix T = (tkl) by
tki = min{t(1, J) | |11 = k, |J| = l} =
= k ■ l + Y^ ri' — Y1 sj'
i>k j<l
for all k = 0,1,... , m and l = 0,1,... , n.
By theorem 2, we see that U'(R', S') is non-empty if and only if T has no non-negative
elements.
Now, we state two lemmas by which we will show that invariant sets are unions of rectangles. Lemma 8. Let U'(R',S') be a normalized class.
(i) If there is a matrix A in the class U'(R', S') such that aij = 0 and aij> = 1 for some 1 < i < m and 1 < j < j' < n, then both (i, j) and (i, j') are variant positions of a switching component in A.
(ii) If there is a matrix A in the class U'(R', S') such that aij = 0 and ai/j = 1 for some 1 < i < i' < m and 1 < j < n, then both (i, j) and (i', j) are variant positions of a switching component in A.
What this lemma says is that, invariant 1’s come before variants and variants come before invariant 0’s from left to right and from top to down.
Lemma 9. Let U'(R',S') be a normalized class, let ¿1 and i2 be two rows which contain variant positions such that 1 < i1 < i2 < m and let [j^, jY] and [j2, j2'] be the corresponding integer ranges of the (consecutive) variant positions. Then, either
ji > j2
or
both ji = j'2 and ji' = jV.
Joining these two lemmas and some more observations, we conclude the following theorem.
Theorem 10. The variant set of a normalized class U' = U'(R',S') = 0, can always be written as
V(U') = up_i/, X Jq
(p=0 if there are no variant elements), where
, • • •, iq'}, 1 < ¿1 < ¿1' < ¿2 < ¿2' < ■ ■ ■ < < ip' < m,
Jq {jq, • • • , jq }, 1 < jp < jp < jp-1 < jp-1 < ■ ■ ■ < j1 < j1 < n.
By this theorem, we can deduce that invariant sets are unions of discrete rectangles between invariant 1’s and invariant 0’s. Then, we can also conclude that invariant 1’s and invariant 0’s are unions of discrete rectangles as well.
Here, we will define another concept called compatibility, which is close to additivity and will be useful to decide whether a position is variant or not.
Definition 7. Let A = (aj) be an mxn binary matrix. A pair of vectors X = (x1, • • • ,xm) £ Rm and Y = (y1,...,yn) £ Rn is said to be compatible with A if for i = 1, •••,■№, and j = ^.• •,n
/ > 0, if aij = 1,
Xi + < 0. if a„ =0.
Theorem 11. Let R = (r1,... , rm) and S = (s1, • • • , sn) be vectors of nonnegative integers
such that there is a binary matrix A £ U(R, S). Then, for i = 1,... , m and j = 1,... , n, (i, j)
is not a variant position if and only if there exists a pair of vectors X = (x1,... ,xm) and Y = (y1,... , n) which is compatible with A such that
Xi + Vj = 0.
Therefore, if we have a pair of vectors which are compatible with the matrix, it is easy to find the structure of the class U'(R', S'). Accordingly, by reversing the permutations of the column and row sums, we can obtain the structure of the class U(R, S).
To finalize this section, we give the algorithm to reconstruct the unique binary matrix in U(R, S) from its row and column sums. Since we know that it is unique, just determining the rows with maximum 1’s and putting enough 1’s to the columns with highest sums is sufficient.
As it is seen, the complexity of this algorithm is much better than the Algorithm 1 as it is
expected.
Algorithm 2.
Input : a compatible pair of vectors, (R, S), satisfying
nn
s'j = sj, for 2 < ^ <n; j= j=^
Step 1. A = 0 (zero matrix);
Step 2. find i1, i2,..., im such that
ril > ri2 > ■ ■ ■ > rim ;
Step 3. for j = 1 to n,
for k =1 to Sj,
aifcj = 1;
Output : matrix A.
Example 4. Let R = (4, 3, 2, 0) and S = (2, 3, 3,1) be given. There is a unique matrix in the class U(R, S) because S = S' hence
sj = Sj, for 2 < l < n.
The matrix is
jj j=< j=i
4 3 2 0
2 /1 1 0 0 \
3 1 1 1 0
3 1 1 1 0
1 1 0 0 0
3. More Projections in Higher Dimensions
In this section, we will first give the definition of the terms used to develop the problem for more than two projections [7]. These definitions will often be generalizations of the ones given in Section 2. We will explain their correspondence where it is necessary. Here, instead of lattice lines our projections will be taken from a set H = {h1,... , hm}, called discrete Radon base. It consists of linear subspaces of Rn with dimensions 1 to n — 1. These subspaces should satisfy the following properties:
1. Their intersection should only consist of the origin:
m
Pi hk={0}. k=1
2. Each one of them does not contain any other:
j = k ^ hj ^ hk •
3. Their intersection with the integers should contain at least two points:
n
n
|hj n Zn| > 2.
The family of translates of hk will be denoted by Hk. These translates have nonempty intersections with Zn. A translate h is
h = hk + a where a £ Rn
such that hk + a fl Zn is nonempty.
Finally, we can define the family H as
m
H = y Hk •
k=1
We will define here the projection function by a transformation. The discrete Radon transform is the mapping vS : H ^ N0 such that for all h £ H we have
vS(h) = |S if h|.
Consider a lattice set S and a discrete Radon base H. We say that S is a set of uniqueness with respect to S, or S is H-unique if vS = vT for some T C Zn then S = T.
Suppose that there exists a mapping g : H ^ R such that for all x £ Zn we have
X £ S ^ (h)g(h) > °. (7)
Then, S is additive with respect to H, or S is H-additive.
Now, let us see that the definition of additivity in two dimensions (R2) and in more dimensions (Rn) are the same. First, let us recall the first definition.
Recall. An m x n binary matrix A = (aij) is additive if there are vectors X = (x1,... , xm) £ Rm and Y = (y1, • • • , yn) £ Rn such that for i = 1,... , m and j = 1,... , n, aij = 1 if and only if xi + Vj > 0.
In the first definition, we mention matrices, while in the second we mention lattice sets. However, we can consider the 1’s of the matrix as our lattice set. Hence, reconstructing a matrix A is in fact nothing but reconstructing the lattice set S, which corresponds to the 1’s of the matrix A, in the R2 plane.
Now, assume that an m x n binary matrix A = (aij) is additive. Therefore, there are vectors X = (x1,... ,xm) £ Rm and Y = (y1,... ,yn) £ Rn such that, for i = 1,..., m and j = 1,... , n, aij = 1 if and only if xi + Vj > 0. Also, since A exists, the corresponding lattice set S C Z2, which consists of the 1’s of A, exists as well. Now, we should find a function g : H ^ R such that for all x £ Zn we have
x £ S ^ (h)g(h) > 0,
heH
where H has two elements h1 and h2, namely h1 being the x-axis and h2 being the y-axis. These are our two directions, and we refer to H = {hk + a : hk + a If Z2 = 0} for k =1, 2. We will construct the function g in the following way.
For the vectors X = (x1,... ,xm) and Y = (y1,..., yn), let xi = g(h1i) and yj = g(h2j) for i = 1,..., m and j = 1,..., n, where h1i = h1 + ai and h2j = h2 + bj for some ai, bj £ R2 and i = 1,... , m, j = 1, • • • , n. For other h £ H let g(h) = 0. Since there are only two directions, there are only two lattice lines that pass through any point in Z2. Therefore, for the point
aij £ S only g(h1i) and g(h2j) pass through aij. Hence, only vaij(h1i) = 1 = vaij(h2j). For any other lattice line, say h', (h') = 0. Hence,
X] (h)g(h) = g(h1i) + g(h2j) = xi + yj > 0,
heH
since we know that xi + yj > 0.
Also, by a similar argumentation, we can also construct the vectors X and Y from the function g : H ^ R.
Hence, the two definitions of additivity are in fact the same if we are working in Z2 and with two lattice directions which are the coordinate axes.
Note that additivity is a stronger concept than uniqueness. In other words, if a set is additive then it is unique, but the converse is not true, which will be explained later on.
In the remainder of this section, we will make use of the mappings of Zn either into the set {0,1} or to the interval [0,1]. Namely,
E : Zn ^ {0,1},
F : Zn ^ [0,1],
where for a mapping f in the union of E and F, f must take positive values for a finite number of x’s. Hence, the set
{x £ Zn : f (x) > 0} is finite. For a given S and H, we define
Es,h = {f £ E : f (h f Zn) = v*(h) Vh £ H},
Fs,h = {f £ F : f (h f Zn)= v*(h) Vh £ H}.
First of all, let us observe that a discrete set S is H-unique if and only if ESH = {xS}, where is the characteristic function of S. If S is not H-unique, then there exists a discrete set T distinct from S with the same discrete radon transform. Therefore, ES,H will contain the characteristic functions of both T and S. Conversely, if ESH does not consist of the characteristic function of S only, there exists a set adjoint from S with the same discrete transform. Hence, S is not H-unique.
Consider Fig. 5 in which we assume that (0, 0, 0) and (1,1,1) define the set S and (0, 0,1)
and (1,1, 0) make up the set T. Observe that when H consists of the three planes perpendicular
to the three axis, S and T will have the same projections with respect to the elements of H. This substructure, which we want to avoid, is called a bad H-rectangle.
To generalize the notion of bad rectangles, let K-bad H-configuration for S, K > 2, be a pair of lists x1,... , xK of K distinct points in S and y1,... , yK of K distinct points in Zn\S such that
K K
(h) = X (h) Vh £ H-
k=1 k=1
A weakly K-bad H-configuration is defined in the same way except for the condition that the elements would be distinct, hence they can appear with multiplicity in the sum above. A bad H-configuration for S is any K-bad H-configuration for S and the same argument holds for a weakly bad H-configuration. We are then faced with the connections of additivity and uniqueness. It is known that additivity implies uniqueness but the converse is valid only when m =2.
Theorem 12. The following conditions are equivalent:
1. S is H-unique.
2. S has no bad H-configuration.
3. es,h = {Xs}.
Theorem 13. The following conditions are equivalent:
1. S is H-additive.
2. S has no weakly bad H-configuration.
3. Fs,h = {Xs}•
The following theorem states that uniqueness and additivity are equivalent when m = 2. This property is validated when m is increased.
Theorem 14. Suppose H = {h1,h2} is a discrete Radon base and S is a nonempty finite subset of Zn. Then, the following conditions are equivalent:
1. S is H-unique.
2. S is H-additive.
3. Fs,H = {Xs}.
4. S has no weakly bad H-configuration.
5. S has no bad H-configuration.
6. S has no bad H-rectangle.
This theorem is extremely important in reducing the time complexity. As it states that, when m = 2, instead of checking all subsets of S to see if they are bad H-configurations for S, one only needs to look at bad H-rectangles. This is not valid when m is greater than two, which is illustrated with an example by Fishburn et al. [7]. This instance lies in 3 dimensions, H consists of the three planes perpendicular to an axis of R3 and the set S is
S = {111,222,333,113,121,122,123,131,133,221, 223, 231, 233, 323}.
Then, it is pointed out that S has no bad H-rectangle but it is not H-unique.
Furthermore, Fishburn and Shepp have shown that for the case of a plane when H consists of three lines through the origin, one can construct a set of 11 points which are unique but not additive.
4. Applications
Although the birth of Discrete Tomography is not long a time ago, it has been used in many applications. Mostly, DT is applied for electron microscopy-techniques, for quality control of industrial products in industrial imaging, for reconstruction of the shape of heart chambers from orthogonal biplane cardiac angiograms in medical imaging and for quality control in VLSI chip design. These applications promise to result in significant improvements in their fields.
The study of Discrete Tomography with medical applications is emerging as an important new research area. Significant applications have been found in medical research, for instance in radiosurgical treatment planning, virtual endoscopy, the reconstruction of the shape of hearth chambers from orthogonal biplane cardiac angiograms [17, 26, 27]. To recover cross-section images from a number of projections, the object is illuminated by X-ray beams from a number of directions. The transmitted rays convey information about the density distribution inside the body. The problem is reconstructing the best approximation of the real cross-section. In order to model the 3 dimensional shape of the left or right heart chambers from the density distributions of orthogonal biplane ventriculograms, Onnasch and Prause [28] introduce a reconstruction method. In [28], techniques of image acquisition and restoration are also presented.
One can make use of coding theory to represent the problem under consideration applied to Very Large Scale Integration (VLSI) systems chip design (Weber [35]). If we represent any existing atom cluster inside of a grid containing our discrete points to be reconstructed and interpret it as a word over the alphabet in the {0,1} field, linear coding theory may help us in error detection and correction in the reconstruction process [21]. To this end, one looks for Hamming codes or asks for the extent of cyclicity. To detect roughness in layer structures on chips wavelets can be used [15, 35]. There are also other approaches coming from optimal experimental design. Bertram et al. [3] use methods from experimental design, which were introduced by Gaffke and Heiligers [8, 9], in the optimization of crystal structures (cf. [35]).
In industrial imaging, to obtain shape and dimensional information of industrial parts, CT has been used as an important and powerful tool [6], for instance in reconstructing radioactive materials. Pointing out that most of the objects are made up of one material, one can use DT by representing the related material with 1 and air with 0. In [4], the authors point out how DT may act as a valuable application in industrial imaging and manufacturing. This may be applied to a wide range of materials [32].
The revival of interest in DT problems in the past decade was motivated by new developments in electron-microscopy [24, 34]. Being able to scan objects at the atomic scale, is not just important for research purposes. Quoting from [30] let us explain the usage of DT in this field. It has now become possible to count the number of atoms lying on each line in a crystal lattice, measuring along the lattices principal directions. When attempting to reconstruct the crystal structure from these partial data, continuous tomographic techniques cannot be used, because the unknown image is binary instead of real-valued: A lattice cell either contains one or no atom, but it certainly will not contain half an atom. Introducing this extra constraint transforms the problem from an analytic problem into a combinatorial problem. Because of the additional constraints on the unknown image, one may hope that it is possible to reduce the number of measurement directions that is necessary for reconstruction. In the particular application of electron-microscopy, this is even a necessity, since measuring the crystal lattice destroys some of its structure. The new theory helps to improve existing electron-microscopy techniques.
Conclusion
Discrete Tomography is a promising field of mathematics which is developing rapidly. As we have seen in the previous section, it has various applications such as in medicine, brain tomography, DNA-microarray chips [12]. Nevertheless, its use requires appropriate modeling of the problem and efficient implementation of the algorithms. Hence, further research in this field suggests significant improvements.
In this paper, we have reflected the main problems that Discrete Tomography is concerned with and their state-of-the-art. The existence problem in two projections and two dimensions is closer explained and algorithms for solving it are given, as we mentioned in Sections 2 and 3. On the other hand, if the number of projections increases, the problem gets more complex. The theory needed and the concepts of additivity and uniqueness are given in Section 4.
Finally, let us mention some open problems. It has been pointed out by Fishburn and Shepp [7] that for a discrete Radon base with at least two elements it is difficult to construct sets S which are unique but not additive. It stimulated the idea that the number of non-similar sets that are unique but not additive tends to zero as the cardinality of S increases.
Another open problem considers a grid with a known number of objects of the same size to be put on each row or column. The aim is to cover the whole area with the objects. The problem is solvable in polynomial time if there is only one object which was considered in detail in Section 3. The computational complexity of the two object case is open. For three or more objects the problem is NP-complete. It is also known that it is NP hard for six or more objects.
After examining the theory of the problem, our future work will concentrate on applications of Discrete Tomography, especially on medical applications. In addition, by developing practical algorithms we will work on improving and extending the existing applications further.
Acknowledgement: We would like to express our gratitude to Institute of Applied Mathematics of Middle East Technical University, especially, to Prof. Dr. Aydin Aytuna and Prof. Dr. Bulent Karasozen, for encouragement and support of our Working Group of Inverse Problems (IAM, METU), and likewise to Prof. Dr. Yurii Shokin from Journal of Computational Technologies. For many fruitful discussions we thank Dr. Hakan Oktem, Prof. Dr. Nevzat Gencer and Dr. Yesim Serinagaoglu.
References
[1] Aster A., Borchers B., Thürber C. Parameter Estimation and Inverse Problems. N. Y.: Acad. Press, to appear in 2004.
[2] BatenbürG K.J. Analysis and optimization of an algorithm for discrete tomography. Preprint Submitted to Elsevier Science. 2003. 12 p.
[3] Bertram A., Bohlke T., Gaffke N. et al. On the generation of discrete isotropic orientation distribution for linear elastic crystals // J. Elast. 2000. P. 233-248.
[4] Browne J.A., Koshy M. Discrete Tomography Foundations, Algorithms and Applications // Assisted Engineering and Manufacturing. Birkhauser, 1999. P. 345-360.
[5] Dalinger E. Aktuelle Fragestellungen der Diskreten Tomographie und ihre komplexitatstheoretische Behandlung. Diploma Thesis, Darmstadt Univ. of Technology, Department of Mathematics, 2003.
[6] DüSAüSSOY N.J., Cao Q., Yancey R.N., Stanley J.H. Image Proc. for ct-assisted reverse engineering and part characterization // Proc. 1995 IEEE Intern. Conf. on Image Proc. 1995. Vol. 3. P. 33-36.
[7] Fishbürn P.C., Shepp L.A. Discrete tomography foundations, algorithms and applications // Sets of Uniqueness and Additivity in Integer Lattices. Birkhauser, 1999. P. 35-57.
[8] Gaffke N., Heiligers B. Algorithms for optimal design with applications to multiple polynomial regression // Metrika. 1995. Vol. 42. P. 173-190.
[9] Gaffke N., Heiligers B. Approximative designs for polynomial regression: invariance, admissibility, and optimality // Handbook of Statistics / S. Ghosh, C.R. Rao. Elsevier Science. 1996. Vol. 13. P. 1149-1199.
[10] Gardner R.J., Gritzmann P. Discrete tomography foundations, algorithms and applications // Uniqueness and Complexity in Discrete Tomography, Birkhauser, 1999. P. 85-111.
[11] Gardner R.J., Gritzmann P. Reconstructing polyatomic structures from discrete X-rays: NP-completeness proof for three atoms, with Marek Chrobak // Theor. Computer Sci. 2001. Vol. 259. P. 81-98.
[12] Gebert J., Latsch M., Pickl S.W., Weber G.W., Wünschiers R. Towards DNA chip analysis and discrete tomography by reverse problems and optimization // Proc. of 4th MATHMOD Vienna - Fourth Intern. Symp. on Math. Modelling. 2003. P. 1298-1302.
[13] Gencer N.G., Tek M.N. Electrical conductivity imaging via contact-less measurements // IEEE Trans. on Medical Imaging. 1999. Vol. 18, N 7. P. 617-627.
[14] Gritzmann P., Gardner R.J., Prangenberg D. On the computational complexity of reconstructing lattice sets from their X-rays // Discrete Math. 1999. N 202. P. 45-71.
[15] Haase M. An exploration of self-similar structures with wavelets // Similarity Methods / B. Kroplin, St. Rudolph, St. Bruecker (Eds). Institut fur Statik und Dynamik der Luft und Raumfahrtkonstruktionen. Univ. of Stuttgart. 2000. P. 183-192.
[16] Hajdü L., Tijdeman R. An algorithm for discrete tomography // Lin. Alg. Appl. 2001. Vol. 339. P. 147-169.
[17] Hall P. A model for learning human vascular anatomy // DIMACS Series in Discrete Mathematics and Theoretical Computer Science. Title: Discrete Math. Problems with Medical Appl. 2000. Vol. 55. P. 11-27.
[18] Hammer P.O. Problem 2 // Proc. Symp. Pure Math. Vol. VII: Convexity. 1963. Vol. 7. P. 498-499.
[19] Herman G.T., Kuba A. Discrete Tomography Foundations, Algorithms and Applications. Birkhauser, 1999.
[20] Honghui W. Structure and cardinality of the class u(r, s) of (0,1)-matrices // J. Math. Research and Exposition. 1984. Vol. 4. P. 87-93.
[21] Ihringer Th. Diskrete Mathematik. Teubner Verlag, 1999.
[22] Irwing R.W., Jerrum M.R. Three-dimensional statistical data security problems // SIAM J. Comput. 1994. Vol. 23 P. 170-184.
[23] KARBEYAZ B.U., Gencer N.G. Electrical conductivity imaging via contactless measurements: An experimental study // IEEE Trans. on Medical Imaging. 2003. Vol. 22, N 5. P. 627-635.
[24] Kisielowski C., Schwander P., Baumann F.H., Siebt M. An approach to quantitative highresolution transmission electron microscopy of crystaline materials // Ultramicroscopy. 1995. Vol. 58. P. 417-432.
[25] Kuba A., Herman G.T., Matej S., Todd-Pokropek A. Medical applications of discrete tomography. Preprint Univ. of Szeged, Hungary, 2003.
[26] Lee Y.-J., Mangasarian O.L., Wolberg W.H. Breast cancer survival and chemotherapy: A support vector machine analysis // DIMACS Series in Discrete Math. and Theor. Computer Sci. Title: Discrete Mathematical Problems With Medical Appl. 2000. Vol. 55. P. 1-10.
[27] Nystrom I., Smedby O. Analysis of magnetic resonance angiography images using skeletonization and distance transforms // Ibid. 2000. P. 75-90.
[28] Onnasch D.G.W., Prause G.P.M. Discrete tomography foundations, algorithms and applications // Heart Chamber Reconstruction from Biplane Angiography. Birkhauser, 1999. P. 385-401.
[29] Reiber J.H.C., Gerbrands J.J., Troost G.J. et al. 3d reconstruction of coronararterial segments from two projections // Digital Imaging in Cardiovascular Radiology. 1983. P. 151-163.
[30] Riele H., Batenburg J. Research project: Discrete Tomography. Centrum voor Wiskunde en Informatica Annual Report, Mathetical Institute in Leiden. Netherlands, 2002.
[31] Ryser H.J. Combinatorial properties of matrices of zeros and ones // Can. J. Mathematics. 1957. Vol. 9. P. 371-377.
[32] Salzberg P., Figueroa R.R. Discrete tomography foundations, algorithms and applications // Tomography on the 3D-Torus and Cyrstals. Birkhauser, 1999. P. 385-401.
[33] Schwander P., Fishburn L., Shepp P., Vanderbei. The discrete radon transform and its approximate inversion via, linear programming // Discrete Appl. Math. 1997. Vol. 75, N 1. P. 39-61.
[34] Schwander P., Kisielowski C., Baumann F.H. et al. Mapping projected potential, interfacial roughness, and composition in general crystaline solids by quantitative transmission electron microscopy // Phys. Review Lett. 1993. Vol. 71. P. 4150-4153.
[35] Weber G.-W. On the topology of generalized semi-infinite optimization // J. of Convex Analysis. 2002. Vol. 9, N 2. P. 665-691.
[36] Yung Kong T., Herman G.T. Discrete Tomography Foundations, Algorithms and Applications // Tomographic Equivalence and Switching Operations. Birkhauser, 1999. P. 59-83.
Received for publication April 23, 2004
