optimal multileaf collimator eld segmentationkalinowski/papers/tk_diss.pdf · 2013. 3. 6. ·...

Optimal multileaf collimator

field segmentation

Dissertation

zur Erlangung des akademischen Grades

doctor rerum naturalium (Dr. rer. nat.)

der Mathemathisch–Naturwissenschaftlichen Fakultat

der Universitat Rostock

vorgelegt von

Thomas Kalinowski, geb. am 27.11. 1980 in Pasewalk

aus Rostock

Rostock, 1. September 2004

Dekan: Prof. Dr. Udo Kragl

Gutachter: Prof. Dr. K. Engel (Rostock)

Gutachter: Prof. Dr. H. Hamacher (Kaiserslautern)

Gutachter: Prof. Dr. R. Ahlswede (Bielefeld)

Tag der offentlichen Verteidigung: 28.01.2005

ACKNOWLEDGEMENT

It is a great pleasure for me to thank Prof. Konrad Engel for his reliableadvice concerning all questions arising in connection with my dissertation,and especially for his patience with several versions of this work. I also thankmy family and friends for always supporting me. Special thanks to MarkusStoy for many useful LATEX hints.

CONTENTS

1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2. TNMU–optimal field segmentation . . . . . . . . . . . . . . . . . . 152.1 A Linear Programming formulation . . . . . . . . . . . . . . . 152.2 The lower bound . . . . . . . . . . . . . . . . . . . . . . . . . 172.3 The algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.4 Test results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3. A heuristic for the reduction of the number of segments . . . . . . . 493.1 The algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.2 Test results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4. Exact minimization of the number of segments . . . . . . . . . . . . 594.1 Single row intensity maps . . . . . . . . . . . . . . . . . . . . 614.2 Multiple row intensity maps . . . . . . . . . . . . . . . . . . . 724.3 Test results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

5. Further Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 855.1 Using the MLC in two directions . . . . . . . . . . . . . . . . 85

5.1.1 A lower bound . . . . . . . . . . . . . . . . . . . . . . 855.1.2 Heuristic results . . . . . . . . . . . . . . . . . . . . . . 87

5.2 Two orthogonal MLCs . . . . . . . . . . . . . . . . . . . . . . 92

1. INTRODUCTION

In cancer treatment high energetic radiation is used to destroy the tumor. Toachieve this goal the irradiation process must be planned in such a way thatthe tumor (target volume) receives a sufficiently high dose while the organsclose to it (organs at risk) are not damaged. In clinical practice the radiationis delivered by a linear accelerator and the beam head is part of a gantry thatcan be rotated about the treatment couch (see Figure 1.1). In order to design

Fig. 1.1: A treatment couch with gantry and beam head.

a treatment plan a set of gantry angles and corresponding beam intensitieshave to be determined. See [13, 21, 22] for reviews of different approaches tothis problem. Because of technological developments recently intensity mod-ulated radiotherapy (IMRT) became an important method to improve thequality of the treatment. Here instead of one rectangular homogeneous fielda number of differently shaped fields are superimposed, and so an intensitymodulated field is produced. The shapes of the fields are determined by amultileaf collimator (MLC). A multileaf collimator consists of two opposingbanks of metal leaves (see Figure 1.2). The leaves can be moved independentof each other and each leaf configuration corresponds to a radiation field of aparticular shape that is realized by inserting the MLC between the radiationsource and the patient. Now the realization of an intensity modulated field

10 1. Introduction

Fig. 1.2: The leaf pairs of a multileaf collimator.

using an MLC can be modeled as follows. The beam head is discretized intobixels and the modulated field is described by giving the desired intensityfor each bixel. So the intensity map can be considered as an m× n−matrix,where each row corresponds to a leaf pair. There are two essentially differentways to generate intensity modulated fields with multileaf collimators: in thestatic mode (step–and–shoot) [6–8, 10, 11, 16, 20, 23, 28] the beam is switchedoff when the leaves are moving while in the dynamic mode [7, 9, 17, 19, 24–27]the beam is switched on during the whole treatment and the modulation isachieved by varying the speed of the leaf motion. Two important criteriafor the quality of a treatment plan are the total irradiation time and thetotal treatment time. The total irradiation time should be small since thereis always a small amount of radiation transmitted through the leaves, andif the used model ignores this leaf transmission the error increases with thetotal irradiation time. A small total treatment time is desirable for efficiencyreasons. In the dynamic mode the two criteria coincide. So here the problemis to determine a velocity function for each leaf such that the given intensityis realized in the shortest possible time. In the static mode the whole treat-ment consists of the irradiation and the intervals in between when the leavesare moved. Thus we have two parameters which influence the total treat-ment time: the irradiation time and the number of homogeneous fields thatare needed. How these parameters have to be weighted depends on the usedtechnology: the longer the time intervals between the different fields are, themore the reduction of the number of fields becomes important. The lengthsof these time intervals are influenced by the leaf velocity and by the so calledverification and record overhead, which is the time necessary to check thecorrect positions of the leaves. In a more realistic model one should also take

11

the shapes of the fields into account, because clearly the necessary leaf traveltime between two fields depends on the shapes of these fields [4, 23]. Thedynamic mode has the advantage of a smaller total treatment time, but thestatic mode involves no leaf movement with radiation on and so the verifi-cation of the correct realization of the treatment plan is easier which makesthe method less sensitive to malfunctions of the technology.

There are additional machine–dependent restrictions which have to beconsidered when determining the leaf positions:

Interleaf collision constraint (ICC): In some widely used MLCs it is forbid-den that opposite leaves of adjacent rows overlap, because otherwisethese leaves collide. So leaf positions as illustrated in Figure 1.3 (wherethe shading indicates the area that is covered by the leaves) are notallowed.

Fig. 1.3: Leaf position that is excluded by the ICC.

Tongue and groove constraint: In order to reduce leakage radiation betweenadjacent leaves the commercially available MLCs use a tongue–and–groove (or similar) design with the effect that there is a small overlapof the regions that are covered by adjacent leaves. This is illustratedin Figure 1.4 which shows a cut through the leaf bank perpendicularto the direction of leaf motion.

Radiation

Fig. 1.4: The principle of the tongue–and–groove design.

12 1. Introduction

Consider two bixels x and y that are adjacent along a column andtwo homogeneous fields, where in the first field x is irradiated and y iscovered and in the second field y is irradiated and x is covered. Thenin the composition of these fields along the border of x and y there is anarrow strip (the overlap of the regions that are covered by the leavesin the rows of x and y, respectively) that receives no radiation at all.Figure 1.5 illustrates this for the intensity map ( 2 3

3 4 ).

3 MU

1 MU

1 MU

2 MU

1 MU

1 MU

a) The overlap of bixels (1, 1) and(2, 1) receives no radiation becauseof the tongue and groove effect.

b) The overlaps of bixels that are ad-jacent along a column receive thesmaller one of the relevant doses.

Fig. 1.5: Two different realizations of the same intensity map.

To avoid this effect one may require that two bixels that are adjacentalong a column are irradiated simultaneously for the time the lower ofthe two doses is delivered. Then the border region receives this lowerdose. If this is the case for all the relevant pairs of adjacent bixels thetreatment plan is said to satisfy the tongue and groove constraint.

In this work we consider the problem of constructing segmentations, tak-ing into account the interleaf collision constraint, while we neglect the tongueand groove constraint. Describing the possible leaf positions of the MLCby certain (0, 1)–matrices, called segments, this amounts to the search fora realization of the given intensity matrix A with a small total number ofmonitor units (TNMU) and a small number of segments (NS). In generalit is not possible to minimize both parameters simultaneously [15]. For thecase of an MLC without ICC there are several segmentation algorithms [6–8, 11, 20, 23, 28], some of them providing the optimal TNMU but a large NS,others reducing the NS heuristically at the price of an increased TNMU.The exact minimization of the NS is NP–hard already for intensity matricesconsisting of just one row [2, 15]. In principle both, TNMU and NS, can be

13

optimized by integer programming and this method can be adapted to addi-tional restrictions like ICC [18]. But clearly this is applicable only for smallproblem sizes. Another approach is the reformulation of the segmentationproblem in a network flow setting. In [4] this is done for MLC–segmentationwith ICC, yielding a network flow algorithm for the TNMU–optimal seg-mentation. In [2] this approach is developed further and a heuristic for thereduction of the NS is added. The method of [23] yielding TNMU–optimalsegmentations without ICC is modified in [16] such that it takes into accountthe ICC (and an even more general condition). In [10] there is presented anefficient segmentation algorithm yielding the optimal TNMU and a very smallNS for the segmentation problem without ICC. See [15] for a survey and acomparison of the different segmentation algorithms. Engel’s algorithm [10]for the segmentation without ICC is derived from an explicit formula for thesmallest possible TNMU. Evaluating this formula is equivalent to solving alongest path problem in a properly constructed layered digraph. Theorem 1,one of the main results of this thesis, is a generalization of that constructionsuch that the longest path problem in the new digraph corresponds to theevaluation of the minimal TNMU for a segmentation with ICC. Chapter 2 isdevoted to the proof of this theorem which consists of two parts: using a du-ality argument we establish that the maximal weight of a path in the digraphto be described below is a lower bound for the number of monitor units ina segmentation, and then we describe a method to achieve this bound. InChapter 3 we develop a heuristic method to construct a segmentation withminimal number of monitor units and a small number of segments. In Chap-ters 4 and 5 we discuss some results for MLC segmentation without ICC.In Chapter 4 we show that a segmentation with in first instance minimalTNMU and in second instance minimal NS can be determined in polynomialtime if the entries of the intensity matrix are bounded by a constant. Finally,in Chapter 5 we suggest two more flexible ways of using the MLC, thus in-creasing the number of possible shapes for the homogeneous fields. In 5.1 westudy the situation when it is allowed to rotate the MLC about 90◦ betweenthe delivery of two segments and in 5.2 we consider a pair of two MLCs withorthogonal directions of leaf motion.

14 1. Introduction

2. AN ALGORITHM FOR TNMU–OPTIMAL FIELDSEGMENTATION WITH INTERLEAF COLLISION

CONSTRAINT

2.1 A Linear Programming formulation

Throughout m and n are positive integers and for positive integers k we usethe notation [k] := {1, 2, . . . , k}. Let A = (ai,j) 1≤i≤m

1≤j≤nbe an m×n-matrix with

nonnegative integer entries. In addition we put ai,0 = ai,n+1 = 0 (i ∈ [m]). Asegment is a matrix that corresponds to a position of an MLC with interleafcollision constraint. This is made precise in the following definition.

Definition 1. A segment is an m×n-matrix S = (si,j), such that there existintegers li, ri (i ∈ [m]) with the following properties:

li ≤ ri + 1 (i ∈ [m]), (2.1)

si,j =

{

1 if li ≤ j ≤ ri

0 otherwise(i ∈ [m], j ∈ [n]), (2.2)

ICC: li ≤ ri+1 + 1, ri ≥ li+1 − 1 (i ∈ [m − 1]). (2.3)

The interpretation is that li−1 and ri +1 are the positions of the i–th leftand right leaf, respectively. A segmentation of A is a representation of A asa sum of segments, that is A =

∑ki=1 uiSi with segments Si (i = 1, 2, . . . , k)

and positive integers ui (i = 1, 2, . . . , k). The TNMU of this segmentation is∑k

i=1 ui and our goal is to find a segmentation of A with minimal TNMU.

Example 1. A segmentation with 10 MU for a benchmark matrix from [18]is

(

4 5 0 1 4 52 4 1 3 1 42 3 2 1 2 45 3 3 2 5 3

)

= 3

(

0 0 0 0 1 10 0 0 0 0 10 0 0 0 0 10 0 0 0 1 1

)

+ 3

(

1 1 0 0 0 00 1 0 0 0 00 0 0 0 0 01 0 0 0 0 0

)

+ 1

(

0 0 0 0 0 10 0 0 1 1 10 0 0 0 1 00 1 1 1 1 0

)

+ 1

(

0 0 0 1 1 10 0 0 1 0 00 1 1 0 0 01 1 1 0 0 0

)

+ 1

(

0 1 0 0 0 01 0 0 0 0 01 1 0 0 0 01 1 1 1 1 0

)

+ 1

(

1 1 0 0 0 01 1 1 1 0 01 1 1 1 1 10 0 0 0 0 0

)

.

(2.4)

16 2. TNMU–optimal field segmentation

By F we denote the family of subsets of V := [m] × [n] that correspondto segments, precisely

F = {T ⊆ V : There exists a segment S with (i, j) ∈ T ⇐⇒ si,j = 1}.

With a segmentation A =∑k

i=1 uiSi we can associate a function f : F → IN :for 1 ≤ i ≤ k we put f(T ) = ui for the T ⊆ V corresponding to the segmentSi, and for the remaining T we put f(T ) = 0. Now the LP-relaxation of thesegmentation problem is:

(P )

minimize∑

T∈F

f(T ) subject to

f(T ) ≥ 0 ∀T ∈ F ,

∑

T∈F :(i,j)∈T

f(T ) = ai,j ∀(i, j) ∈ V.

The dual variables (one variable for each (i, j) ∈ V ) can be considered as afunction g : V → IR and in this formulation the dual program is

(D)

maximize∑

(i,j)∈V

ai,jg(i, j) subject to

∑

(i,j)∈T

g(i, j) ≤ 1 ∀T ∈ F .

To solve the segmentation problem we proceed in two steps: first we constructa feasible solution for the program (D) which yields a lower bound for theTNMU, and in the second step we construct a sequence of segments that

realizes this lower bound. We define a directed acyclic graph−→G = (V ∪

{0, 1}, E). For E we take all possible arcs of the forms (0, (i, 1)) and ((i, n), 1),as well as all the arcs between the j−th and the (j + 1)−th column (j =1, 2, . . . , n − 1), precisely E = E1 ∪ E2 ∪ E3, where

E1 = {(0, (i, 1)) : i ∈ [m]},

E2 = {((i, n), 1) : i ∈ [m]},

E3 = {((i, j), (i′, j + 1)) : i, i′ ∈ [m], j ∈ [n − 1]}.

The next step is the definition of a weight function δ (depending on A) onE. δ reflects the structure of certain dual feasible solutions g (to be definedin the next section), in such a way that the objective value of the program

2.2. The lower bound 17

(D) for a solution g equals the weight of a certain (0, 1)–path in−→G . We put

di,j = ai,j − ai,j−1 (i ∈ [m], j ∈ [n + 1]) and

δ(0, (i, 1)) = ai,1 (i ∈ [m]),

δ((i, n), 1) = 0 (i ∈ [m]),

δ((i, j), (i, j + 1)) = max{0, di,j+1} (i ∈ [m], j ∈ [n − 1]),

δ((i, j), (i′, j + 1)) = max{0, di′,j+1} −i′−1∑

k=i

ak,j

(i, i′ ∈ [m], i < i′, j ∈ [n − 1]),

δ((i, j), (i′, j + 1)) = max{0, di′,j+1} −i∑

k=i′+1

ak,j

(i, i′ ∈ [m], i > i′, j ∈ [n − 1]).

Since the considered matrix will be clear from the context we omit it in thenotation for the weight function. For a path P = (v0, v1, . . . , vl) in

−→G its

weight is δ(P ) =∑l

i=1 δ(vi−1, vi). Now we are prepared to formulate themain result of this chapter.

Theorem 1. The minimal TNMU of a segmentation of a nonnegative matrix

A equals the maximal weight of a (0, 1)−path in−→G .

Note that the minimal TNMU in a segmentation without ICC can beinterpreted analogously. In this case the minimal TNMU equals (see [10])

max1≤i≤m

n∑

j=1

max{0, di,j},

that is the maximal weight of a (0, 1)−path in the graph that is obtained

from−→G by deleting all the arcs ((i, j), (i′, j + 1)) with i 6= i′. For notational

convenience we put

c(A) = max{δ(P ) : P is a (0, 1) − path in−→G},

so that in order to prove the theorem we have to show that c(A) is a lowerbound for the TNMU of a segmentation of A and that this bound is sharp.

2.2 The lower bound

In this section we show how the (0, 1)−paths in−→G correspond to certain

feasible solutions for the program (D) and from this we derive the lower


bound part of Theorem 1. A (0, 1)−path P is uniquely determined by theindices of the columns in which P changes the row and the indices of therows in which P runs between the row changes. So let x1, x2, . . . , xk−1 with0 < x1 < x2 < · · · < xk−1 < n denote the indices of the columns where P

changes the row, i.e.

(i, xt), (i′, xt + 1) ∈ P with i 6= i′ (t ∈ [k − 1]),

and let i∗t be the row index with (i∗t , xt) ∈ P (t = 1, 2, . . . , k − 1) and i∗k theindex with (i∗k, n) ∈ P . Finally, we put x0 = 0, and xk = n + 1. Thus

P = (0, (i∗1, 1), (i∗1, 2), . . . , (i∗1, x1), (i∗2, x1 + 1), . . . , (i∗2, x2), . . . , . . . , (i

∗k, n), 1),

and P is uniquely determined by its parameters (i∗1, x1), . . . , (i∗k, xk). Now we

define g : V → {1,−1, 0} by

g(i, j) =

1 if xt−1 < j < xt − 1, i = i∗t , di,j ≥ 0, di,j+1 < 0,1 if xt−1 < j = xt − 1, i = i∗t , di,j ≥ 0,−1 if xt−1 < j < xt − 1, i = i∗t , di,j < 0, di,j+1 ≥ 0,−1 if j = xt, i∗t ≤ i < i∗t+1 or i∗t+1 < i ≤ i∗t ,

−1 if j = xt, i = i∗t+1, di,j+1 ≥ 0,0 otherwise.

(2.5)

For i ∈ [m] we put J(i) := {j ∈ [n] : g(i, j) 6= 0}. Fix some i ∈ [m] anddenote the elements of J(i) by j1, j2, . . . , jp such that

j1 < j2 < · · · < jp.

Then the following observations follow immediately from (2.5).

1. If i = i∗1 and k > 1 then x1 ∈ J(i) and the sequence

g(i, j1), g(i, j2), . . . , g(i, x1)

is an alternating (1,−1)−sequence ending with −1.

2. If i = i∗k, k > 1 and J(i) ∩ {xk−1, xk−1 + 1, . . . , n} 6= ∅ then for

q = min{τ : jτ ≥ xk−1},

g(i, jq), g(i, jq+1), . . . , g(i, jp) is an alternating (1,−1)−sequence start-ing with −1 and ending with 1.

3. If i = i∗1 and k = 1 then the sequence g(i, j1), g(i, j2), . . . , g(i, jp) isempty or an alternating (1,−1)−sequence starting and ending with 1.


· //· //

3

��***

****

****

*** · //· 2 //· 2 //·

��· //· //· //· //· //·

''·

3;;wwwwwwwwwww

1 33

2++

##

·· //· //· //· //· //·

77

· //· 6 //· //· 1 //· //·

DD

(gi,j) =

1 −1 0 0 0 00 −1 0 0 0 00 −1 0 0 0 00 −1 1 −1 0 1

Fig. 2.1: A path P and the corresponding g.

4. If i = i∗t for 2 ≤ t ≤ k − 1 then xt ∈ J(i) and for

q = min{τ : jτ ≥ xt−1},

g(i, jq), g(i, jq+1), . . . , g(i, xt) is an alternating (1,−1)−sequence start-ing and ending with −1.

5. If j ∈ J(i) and (i, j) does not correspond to a term in one of thesequences described in the first 4 cases then j = xt for some t ∈ [k− 1]with i 6= i∗t and i 6= i∗t+1 and g(i, j) = −1.

Example 2. Fig. 2.1 shows a path P of weight 7 with respect to the matrix

A =

(

3 0 0 0 2 41 1 1 2 3 32 2 2 1 1 10 0 6 0 1 1

)

and the corresponding function g. The dotted lines are some more arcslabeled with their weight, where the unlabeled arcs have weight 0.

In order to prove that for every (0, 1)−path P the corresponding function g

is a feasible solution for the program (D), we have to show that, for everyT ∈ F ,

∑

(i,j)∈T

g(i, j) ≤ 1.

Lemma 1. Let P be a (0, 1)−path with parameters (i∗1, x1), . . . , (i∗k, xk), and

let g be defined according to (2.5). In addition let 1 ≤ l ≤ r + 1 ≤ n + 1.Then, for every i ∈ [m],

r∑

j=l

g(i, j) ≤ 1,

and equality implies xt−1 < l ≤ r < xt for some t ∈ [k] with i∗t = i.


Proof. We choose an arbitrary i ∈ [m]. As above we denote the elements ofJ(i) = {j ∈ [n] : g(i, j) 6= 0} by j1, . . . , jp such that

j1 < j2 < · · · < jp.

As a consequence of the observations before Example 2 we obtain, for 1 ≤q ≤ p − 1,

g(i, jq) = 1 ⇒ g(i, jq+1) = −1.

Now the first part of the lemma follows from

r∑

j=l

g(i, j) =

q′∑

τ=q

g(i, jτ ) for some q, q′ ∈ [p].

Supposer∑

j=l

g(i, j) =

q′∑

τ=q

g(i, jτ ) = 1.

By construction the sequence g(i, jq), g(i, jq+1), . . . , g(i, jq′) has to be an al-ternating (1,−1)−sequence starting and ending with 1. This implies

xt−1 < l < xt and xt′−1 < r < xt′

for some t, t′ ∈ [k] with i = i∗t = i∗t′ . Assume t 6= t′ and put

t′′ = min{σ > t : i∗σ = i} and q′′ = min{τ : jτ ≥ xt′′−1}.

Thenjq < xt < jq′′ < jq′, g(i, xt) = g(i, jq′′) = −1,

and g(i, j) ≤ 0 for all j with xt < j < jq′′ . So g(i, jq), g(i, jq+1), . . . , g(i, jq′)contains two consecutive (−1)−terms, which is a contradiction. Hence t = t′

and the second part of the lemma follows. �

The next lemma gives a condition that must hold if the sum of the g(i, j)over a row of a segment vanishes. (By a row of a segment we mean the partof the row that is left open by the MLC in the corresponding leaf position.)


let g be defined according to (2.5). Assume i ∈ [m], 1 ≤ l ≤ r + 1 ≤ n + 1and

r∑

j=l

g(i, j) = 0.

Suppose in addition that for some t ∈ [k − 1] one of the following conditionsholds


-8 -7 -6 -5 -4 -3 -20

1

2

3

4� � � �

� �

i∗t

i

i∗t+1

li ri xt

a) i∗t < i < i∗t+1 and l ≤ xt

2 3 4 5 6 7 80

1

2

3

4� �

� � �

�

i∗t−1

i∗t = i

i∗t+1

li=xt−1 ri xt

b) t ≥ 2, i∗t = i < i∗t+1 and l ≤ xt−1.

Fig. 2.2: Illustration of Lemma 2.

1. i∗t < i < i∗t+1 and l ≤ xt

2. t ≥ 2, i∗t = i < i∗t+1 and l ≤ xt−1

3. i∗t > i > i∗t+1 and l ≤ xt

4. t ≥ 2, i∗t = i > i∗t+1 and l ≤ xt−1

Then r < xt.

Proof. We consider only the first two cases that are illustrated in Fig. 2.2.The other two are treated analogously. Assume r ≥ xt. In order to derivea contradiction we use the following observation several times. If P leavesrow i in (i, j) then g(i, j) = −1, and if P enters row i′ in (i′, j ′), j ′ > 1, theneither g(i′, j ′ − 1) = −1 or the first nonvanishing g(i′, j ′′) we meet on thesubpath starting with (i′, j ′) equals −1. We put

J = {j : l ≤ j ≤ r, g(i, j) 6= 0},

and denote the elements of J by j1, j2, . . . , jp (j1 < j2 < · · · < jp). Inparticular jq = xt for some q ∈ [p].

Case 1: g(i, j1) = −1.

By assumption g(i, j1), . . . , g(i, jp) is an alternating (1,−1)−sequencestarting with −1 and ending with 1. From g(i, xt) = −1 follows q < p

and by construction of g the contradiction

g(i, jq) = g(i, jq+1) = −1.

Case 2: g(i, j1) = 1.


This implies l < xt′ for some t′ < t with i∗t′ = i, and consequentlyjq′ = xt′ for some q′ ∈ [p − 1], q′ < q. Thus

g(i, jq′) = g(i, jq′+1) = −1,

and g(i, j1), g(i, j2), . . . , g(i, jp) contains two consecutive (−1)−terms.By assumption this implies g(i, jp) = 1, hence p > q, and by construc-tion of g,

g(i, jq) = g(i, jq+1) = −1.

But now g(i, j1), g(i, j2), . . . , g(i, jp) contains two pairs of consecutive(−1)−terms (if q′+1 < q) or three consecutive (−1)−terms (if q ′ +1 =q). Again this yields a contradiction. �

The following lemma is the crucial step in the proof of the feasibility ofg. We show that the ICC implies that in any segment, between two rows inwhich the values of g add up to 1 there is a row in which this sum is at most−1.


let g be defined according to (2.5). Suppose S is a segment described byl1, l2, . . . , lm,r1, r2, . . . , rm and there are row indices i0, i1 (1 ≤ i0 < i1 ≤ m)such that

ri0∑

j=li0

g(i0, j) = 1 and

ri1∑

j=li1

g(i1, j) = 1.

Then there exists a row index i with i0 < i < i1 andri∑

j=li

g(i, j) ≤ −1.

Proof. W.l.o.g. we may assume that there is no row i with i0 < i < i1 and

ri∑

j=li

g(i, j) = 1.

Suppose that for all i with i0 < i < i1,ri∑

j=li

g(i, j) = 0. By Lemma 1 there

are t, t′ ∈ [k] such that

xt−1 < li0 ≤ ri0 < xt, i∗t = i0 and

xt′−1 < li1 ≤ ri1 < xt′ , i∗t′ = i1.

W.l.o.g. we may assume t < t′. Now let i0 = z0 < z1 < · · · < zp = i1be an increasing sequence of row indices such that there is a corresponding


-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 70

1

2

3

4

5

6

7� �

� �

� �

� �

� �

� �

� �

z0 = i0

z1

z2

z3 = i1

ri0 xt0 xt1 xt2 xt3 xt′−1 li1

Fig. 2.3: Situation in the proof of Lemma 3 with p = 3.

sequence t = t0 < t1 < · · · < tp ≤ t′ with i∗tq = zq (0 ≤ q ≤ p) and in additionfor 0 ≤ q ≤ p − 1 there is no τ with tq < τ < tq+1 and zq < i∗τ ≤ zq+1.These sequences always exist, are uniquely determined, and can be obtainedrecursively as follows. We put

t0 = t and z0 = i0,

and for q ≥ 1, if zq−1 < i1,

tq = min{τ : zq−1 < i∗τ ≤ i1} and zq = i∗tq .

So for some q we obtain zq = i1, and then we put p = q (see Fig. 2.3).

Claim 1: For 0 ≤ q ≤ p − 1,

rzq< xtq ⇒ ri < xtq+1−1 for all i with zq ≤ i < zq+1.

Claim 2: For 0 ≤ q ≤ p − 1 we have rzq< xtq .

Proof of Claim 1. We proceed by induction on i. By assumption

rzq< xtq ≤ xtq+1−1.

So let zq < i < zq+1 and assume ri−1 < xtq+1−1. The ICC implies

li ≤ ri−1 + 1 ≤ xtq+1−1,

and by Lemma 2 we obtain ri < xtq+1−1.


Proof of Claim 2. Here we use induction on q. Clearly,

rz0 = ri0 < xt = xt0 .

So let q > 0 and assume by induction rzq−1 < xtq−1 . Then by Claim 1,

rzq−1 < xtq−1.

Thus by the ICC

lzq≤ rzq−1 + 1 ≤ xtq−1,

and hence, again by Lemma 2, rzq< xtq .

Combining Claims 1 and 2 we obtain

ri1−1 < xtp−1 ≤ xt′−1 < li1 ,

thus ri1−1 < li1 − 1 in contradiction to the ICC. �


let g be defined according to (2.5). Then g is feasible for (D).

Proof. Let T ∈ F be arbitrary and let S be the corresponding segment withparameters li, ri (i ∈ [m]). Then

∑

(i,j)∈T

g(i, j) =

m∑

i=1

ri∑

j=li

g(i, j).

By Lemma 1, for all i ∈ [m],ri∑

j=li

g(i, j) ≤ 1, and by Lemma 3 between two

rows i and i′′ with i < i′′ and

ri∑

j=li

g(i, j) =

ri′′∑

j=li′′

g(i′′, j) = 1

there is always a row i′ with i < i′ < i′′ andri′∑

j=li′

g(i′, j) ≤ −1. Consequently,

∑

(i,j)∈T

g(i, j) ≤ 1,

that is the feasibility of g. �



let g be defined according to (2.5). Then

∑

(i,j)∈V

g(i, j)aij =

k∑

t=1

xt−1∑

j=xt−1+1

max{0, di∗t ,j}−k−1∑

t=1

i∗t+1−1∑

i=i∗t

ai,xt+

i∗t∑

i=i∗t+1+1

ai,xt

.

For brevity of notation, here and for the rest of this work we use theconvention that an empty sum is zero, i.e.

∑s

i=r zi = 0 if s < r.

Proof. Immediately from (2.5) it follows that

xt−1∑

j=xt−1

g(i∗t , j)ai∗t ,j =xt−1∑

j=xt−1+1

max{0, di∗t ,j} (t = 1, 2, . . . , k).

The remaining nonzero g(i, j) correspond to the row changes of P , and wehave to add for t = 1, 2, . . . , k − 1,

i∗t+1−1∑

i=i∗t

g(i, xt)ai,xt= −

i∗t+1−1∑

i=i∗t

ai,xtif i∗t < i∗t+1 and

i∗t∑

i=i∗t+1+1

g(i, xt)ai,xt= −

i∗t∑

i=i∗t+1+1

ai,xtif i∗t > i∗t+1. �

For the weight of P to be equal to the value of the program (D) forthe corresponding g we need an additional restriction on P . We call the(0, 1)−path P with parameters (i∗1, x1), . . . , (i

∗k, xk) feasible (with respect to

A) if di∗t ,xt< 0 for t = 1, 2, . . . , k − 1, which in particular implies that the

last arcs of the horizontal parts of P have weight 0.

Lemma 6. Let P be a feasible (0, 1)−path and let g be defined according to(2.5). Then

∑

(i,j)∈V

g(i, j)ai,j = δ(P ).

Proof. Let P be given by parameters (i∗1, x1), . . . , (i∗k, xk). For t ∈ [k] we

denote by Pt the subpath from (i∗t , xt−1 + 1) to (i∗t , xt). Thus

δ(P ) =

k∑

t=1

δ(Pt) + δ(0, (i∗1, 1)) +

k−1∑

t=1

δ((i∗t , xt), (i∗t+1, xt + 1)).


From the feasibility of P follows that the last arc of Pt has weight 0 for allt ∈ [k], and we obtain

δ(Pt) =xt−1∑

j=xt−1+2

max{0, di∗t ,j}.

In addition, δ(0, (i∗1, 1)) = ai∗1 ,1 = max{0, di∗1,1}, and for t ∈ [k − 1],

δ((i∗t , xt), (i∗t+1, xt + 1)) = max{0, di∗t+1,xt+1} −

i∗t+1−1∑

i=i∗t

ai,xt−

i∗t∑

i=i∗t+1+1

ai,xt.

Thus

δ(P ) =

k∑

t=1

xt−1∑

j=xt−1+1

max{0, di∗t ,j} −k−1∑

t=1

i∗t+1−1∑

i=i∗t

ai,xt+

i∗t∑

i=i∗t+1+1

ai,xt

,

and the claim follows by Lemma 5. �

Lemma 7. There exists a feasible (0, 1)−path P with δ(P ) = c(A).

Proof. For any (0, 1)−path P with parameters (i∗1, x1), . . . , (i∗k, xk) denote by

R(P ) ⊆ [k − 1] the subset of indices that destroy the feasibility of P , i.e.

R(P ) = {t ∈ [k − 1] : di∗t ,xt≥ 0}.

Then

λ(P ) =∑

t∈R(P )

∣

∣i∗t − i∗t+1

∣

∣

measures how far P is from being feasible. In particular, λ(P ) = 0 is a neces-sary and sufficient condition for the feasibility of P . Let P0 be a (0, 1)−pathwith parameters (i∗1, x1), . . . , (i

∗k, xk) and weight δ(P0) = c(A). If λ(P0) = 0

then P0 is feasible and there is nothing to do. So we assume that for r ≥ 1we have a (0, 1)−path Pr−1 with parameters

(i∗1, x1), . . . , (i∗k, xk),

δ(Pr−1) = c(A) and λ(Pr−1) > 0. From this we construct a (0, 1)−path Pr

with δ(Pr) = c(A) and λ(Pr) ≤ λ(Pr−1) − 1. This will prove the lemma,since after finitely many steps we obtain a path P with δ(P ) = c(A) andλ(P ) = 0. Let t be the smallest element of R(Pr−1).


Case 1: di∗t ,j ≥ 0 for xt−1 < j < xt.

We define Pr as follows.

1. If i∗t < i∗t+1 −1 and i∗t−1 6= i∗t +1 the parameters of Pr are (see Fig.2.4 and 2.6)

(i∗1, x1), . . . , (i∗t−1, xt−1), (i

∗t + 1, xt), (i

∗t+1, xt+1), . . . , (i

∗k, xk).

2. If i∗t > i∗t+1 +1 and i∗t−1 6= i∗t −1 the parameters of Pr are (see Fig.2.9 and 2.11)

(i∗1, x1), . . . , (i∗t−1, xt−1), (i

∗t − 1, xt), (i

∗t+1, xt+1), . . . , (i

∗k, xk).

3. If i∗t−1−1 = i∗t < i∗t+1−1 or i∗t−1 +1 = i∗t > i∗t+1 +1 the parametersof Pr are (see Fig. 2.7 and 2.12)

(i∗1, x1), . . . , (i∗t−2, xt−2), (i

∗t−1, xt), (i

∗t+1, xt+1), . . . , (i

∗k, xk).

4. If i∗t + 1 = i∗t+1 6= i∗t−1 or i∗t − 1 = i∗t+1 6= i∗t−1 the parameters of Pr

are (see Fig. 2.5 and 2.10)

(i∗1, x1), . . . , (i∗t−1, xt−1), (i

∗t+1, xt+1), . . . , (i

∗k, xk).

5. If i∗t + 1 = i∗t+1 = i∗t−1 or i∗t − 1 = i∗t+1 = i∗t−1 the parameters of Pr

are (see Fig. 2.8 and 2.13)

(i∗1, x1), . . . , (i∗t−2, xt−2), (i

∗t+1, xt+1), . . . , (i

∗k, xk).

Case 2: di∗t ,j < 0 for some j with xt−1 < j < xt.

We put

x := max{j ≤ xt : di∗t ,j < 0, di∗t ,j+1 ≥ 0},

and define Pr as follows.

1. If i∗t < i∗t+1 − 1 the parameters of Pr are (see Fig. 2.14)

(i∗1, x1), . . . , (i∗t−1, xt−1), (i

∗t , x), (i∗t + 1, xt), (i

∗t+1, xt+1), . . . , (i

∗k, xk).

2. If i∗t > i∗t+1 + 1 the parameters of Pr are (see Fig. 2.16)

(i∗1, x1), . . . , (i∗t−1, xt−1), (i

∗t , x), (i∗t − 1, xt), (i

∗t+1, xt+1), . . . , (i

∗k, xk).


3. If i∗t = i∗t+1 − 1 or i∗t = i∗t+1 + 1 the parameters of Pr are (see Fig.2.15 and 2.17)

(i∗1, x1), . . . , (i∗t−1, xt−1), (i

∗t , x), (i∗t+1, xt+1), . . . , (i

∗k, xk).

We have to show that δ(Pr) = c(A) and λ(Pr) ≤ λ(Pr−1) − 1. The lastassertion follows from the fact that either

R(Pr) = R(Pr−1) or R(Pr) = R(Pr−1) \ {t},

and consequently,

λ(Pr) = λ(Pr−1) − 1 or λ(Pr) = λ(Pr−1) −∣

∣i∗t − i∗t+1

∣

∣ .

Now we check that in any case δ(Pr) ≥ δ(Pr−1) and hence δ(Pr) = c(A). In

the following let the vertices of−→G be denoted as in the corresponding figures.

In addition for two vertices X and Y on a path P we denote by DP (X, Y )the weight of the (X, Y )−subpath of P . Then in any case,

δ(Pr) = δ(Pr−1) − DPr−1(U, A) − DPr−1(A, B) − DPr−1(B, V )

+ DPr(U, A′) + DPr

(A′, B′) + DPr(B′, V ). (2.6)

Cases 1.1(a), 1.4(a): (Fig. 2.4, 2.5)

xt−1 xt

//· U

��???

?

·A

//·B

//·

��///

//// i∗t

·V

//· //· i∗t+1

Pr−1

xt−1 xt

//· U

��///

////

·A′

//·B′

//·��?

??? i∗t + 1

·V

//· //· i∗t+1

Pr

Fig. 2.4: Transition from Pr−1 to Pr in Case 1.1.a) i∗t−1 < i∗t .

Using di∗t ,j ≥ 0 for xt−1 < j ≤ xt we obtain

DPr−1(A, B) = ai∗t ,xt− ai∗t ,xt−1+1,

DPr(B′, V ) = DPr−1(B, V ) + ai∗t ,xt

,

DPr(U, A′) = DPr−1(U, A) −

(

ai∗t ,xt−1+1 − ai∗t ,xt−1

)

− ai∗t ,xt−1 + max{0, di∗t +1,xt−1+1}

= DPr−1(U, A) − ai∗t ,xt−1+1 + max{0, di∗t +1,xt−1+1}.


xt−1 xt

//· U

��///

////

·A

//·B

//·��?

??? i∗t

·V

//· //· i∗t+1

Pr−1

xt−1 xt

//· U

��***

****

****

·A′

//·B′

//· //·V

//· //· i∗t+1

Pr


Substituting into (2.6) yields

δ(Pr) = δ(Pr−1) − DPr−1(U, A) −(

ai∗t ,xt− ai∗t ,xt−1+1

)

− DPr−1(B, V )

+(

DPr−1(U, A) − ai∗t ,xt−1+1 + max{0, di∗t +1,xt−1+1})

+ DPr(A′, B′) +

(

DPr−1(B, V ) + ai∗t ,xt

)

,

that is

δ(Pr) = δ(Pr−1) + max{0, di∗t +1,xt−1+1} + DPr(A′, B′)

≥ δ(Pr−1).

Cases 1.1(b), 1.3(a), 1.5(a): (Fig. 2.6, 2.7, 2.8)

xt−1 xt

·A //· //·B

��***

****

**** i∗t

//· U

GG��

·V

//· //· i∗t+1

Pr−1

xt−1 xt

·A′

//· //·B′

��///

//// i∗t + 1

//·U

??��

·V

//· //· i∗t+1

Pr

Fig. 2.6: Transition from Pr−1 to Pr in Case 1.1.b) i∗t−1 > i∗t .

Again,

DPr−1(A, B) = ai∗t ,xt− ai∗t ,xt−1+1,


.

But in these cases

DPr(U, A′) = DPr−1(U, A) −

(


)

+ ai∗t +1,xt−1 + max{0, di∗t +1,xt−1+1}.


xt−1 xt

·A //· //·B

��///

//// i∗t

//·U

??��

·V

//· //· i∗t+1

Pr−1

xt−1 xt

//·U

//·A′

//· //·B′

��???

? i∗t + 1

·V

//· //· i∗t+1

Pr

Fig. 2.7: Transition from Pr−1 to Pr in Case 1.3.a) i∗t−1 > i∗t .

xt−1 xt

·A //· //·B

��???

? i∗t//·

U

??�� ·V

//· //· i∗t+1

Pr−1

xt−1 xt

//·U

//·A′

//· //·B′

//·V

//· //· i∗t+1

Pr


And substituting into (2.6) yields

δ(Pr) = δ(Pr−1) − DPr−1(U, A) −(

ai∗t ,xt− ai∗t ,xt−1+1

)

− DPr−1(B, V )

+[

DPr−1(U, A) −(


)

+ ai∗t +1,xt−1+

max{0, di∗t +1,xt−1+1}]

+ DPr(A′, B′) +

(


)

,

that is

δ(Pr) = δ(Pr−1) + ai∗t ,xt−1 + ai∗t +1,xt−1 + max{0, di∗t +1,xt−1+1} + DPr(A′, B′)

≥ δ(Pr−1).

Cases 1.2(a), 1.4(b): (Fig. 2.9, 2.10)

xt−1 xt ·V //· //· i∗t+1

·A //· B //·

GG�� i∗t//·

U

??��

Pr−1

xt−1 xt ·V //· //· i∗t+1

·A′

//· B′//·

??�� i∗t − 1

//· U

GG��

Pr


The computation is the same as in Case 1.1(a) but in the formula forDPr

(U, A′) we have to replace di∗t +1,xt−1+1 by di∗t−1,xt−1+1.


xt−1 xt

·V //· //· i∗t+1

·A

//·B//·

??�� i∗t//·

U

??��

Pr−1

xt−1 xt

·A′

//·B′//· //·V //· //· i∗t+1

//· U

GG��

Pr


xt−1 xt ·V //· //· i∗t+1

//· U

��///

////

·A

//·B//·

JJ�� i∗t

Pr−1

xt−1 xt ·V //· //· i∗t+1

//· U

��???

?

·A′

//·B′//·

GG�� i∗t − 1

Pr

Fig. 2.11: Transition from Pr−1 to Pr in Case 1.2.b) i∗t−1 < i∗t .

xt−1 xt ·V //· //· i∗i+1

//· U

��???

?

·A

//·B//·

GG�� i∗t

Pr−1

xt−1 xt ·V //· //· i∗i+1

//·U //·A′

//·B′//·

??�� i∗t

Pr


xt−1 xt

//· U

��???

? ·V //· //· i∗i+1

·A

//·B//·

??�� i∗t

Pr−1

xt−1 xt

//·U //·A′

//·B′//· //·V //· //· i∗t+1

Pr



Cases 1.2(b), 1.3(b), 1.5(b): (Fig. 2.11, 2.12, 2.13)

The computation is the same as in Case 1.1(b) but in the formula forDPr

(U, A′) we have to replace di∗t +1,xt−1+1 by di∗t−1,xt−1+1.

Cases 2.1, 2.3(a): (Fig. 2.14, 2.15)

//·��?

??? x xt

· //·U

//·A //·

��///

////B i∗t

·V

//· i∗t+1

Pr−1

//·��?

??? x xt

· //·U ��?

??? i∗t

·A′

//·B′

��???

? i∗t + 1

·V

//· i∗t+1

Pr

Fig. 2.14: Transition from Pr−1 to Pr in Case 2.1.

//·��?

??? x xt

· //·U

//·A //·��?

???B i∗t

·V

//· i∗t+1

Pr−1

//·��?

??? x xt

· //·U ��?

??? i∗t

·A′

//·B′

//·V

//· i∗t+1

Pr


Using di∗t ,j ≥ 0 for x < j < xt, we obtain

DPr−1(U, A) = ai∗t ,x+1 − ai∗t ,x,

DPr−1(A, B) = ai∗t ,xt− ai∗t ,x+1,


,

DPr(U, A′) = max{0, di∗t +1,x+1} − ai∗t ,x

= max{0, di∗t +1,x+1} + DPr−1(U, A) − ai∗t ,x+1,

and so with (2.6)

δ(Pr) = δ(Pr−1) − DPr−1(U, A) −(

ai∗t ,xt− ai∗t ,x+1

)

− DPr−1(B, V )

+(

max{0, di∗t +1,x+1} + DPr−1(U, A) − ai∗t ,x+1

)

+ DPr(A′, B′) +

(


)

,

that is

δ(Pr) = δ(Pr−1) + max{0, di∗t +1,x+1} + DPr(A′, B′)

≥ δ(Pr−1).

2.3. The algorithm 33

x xt ·V //· i∗t+1

· //·U

//·A

//·

GG��B

i∗t//·

??��

Pr−1

x xt ·V //· i∗t+1

·A′

//·

??��B′

i∗t − 1

· //·

??��U

i∗t//·

??��

Pr

Fig. 2.16: Transition from Pr−1 to Pr in Case 2.2.

x xt

·V //· i∗t+1

· //·U

//·A

//·

??��B

i∗t//·

??��

Pr−1

x xt

·A′

//·B′

//·V //· i∗t+1

· //·

??��B

i∗t//·

??��

Pr


Cases 2.2, 2.3(b): (Fig. 2.16, 2.17) The computation is the same as in Case2.1 but in the formula for DPr

(U, A′) we have to replace di∗t +1,x+1 bydi∗t−1,x+1. �

From Lemmas 4, 6 and 7 we deduce by duality that c(A) is a lower boundfor the sum of the coefficients of a segmentation of A and thus we have alreadyproved the first half of the theorem.

2.3 The algorithm

In this section we assume c(A) > 0 and construct a segment S such that A−S

is still nonnegative and c(A − S) ≤ c(A) − 1. Iterating this construction weobtain a sequence of c(A) segments whose sum is A. For (i, j) ∈ V wedenote by α1(i, j) the maximal weight of a (0, (i, j))−path, by α2(i, j) themaximal weight of an ((i, j), 1)−path and by α(i, j) the maximal weight ofa (0, 1)−path through (i, j), that is

α1(i, j) = max{δ(P ) : P (0, (i, j)) − path in−→G},

α2(i, j) = max{δ(P ) : P ((i, j), 1) − path in−→G},

α(i, j) = α1(i, j) + α2(i, j).

Now we define two subsets V1, V2 ⊆ V . In V1 we collect the pairs (i, j)that determine local maxima or right ends of plateaus in the sequences


ai,1, ai,2, . . . , ai,n (i = 1, 2, . . . , m), precisely

V1 = {(i, j) ∈ V : di,j ≥ 0, di,j+1 < 0}.

The second subset V2 is defined to be the set of pairs (i, j) ∈ V1 with thefollowing properties

1. There exists a (0, 1)−path P of weight c(A) through (i, j).

2. The sequence ai,1, . . . , ai,j is nondecreasing, i.e. ai,1 ≤ · · · ≤ ai,j.

3. The horizontal (0, (i, j))−path is a (0, (i, j))−path of maximal weight.

In other words,

V2 = {(i, j) ∈ V1 : α(i, j) = c(A) and α1(i, j) = ai,j}.

Observe that for (i, j) ∈ V1, δ((i, j), (i, j + 1)) = 0 and thus, for j ′′ > j,

δ((0, (i, 1), (i, 2), . . . , (i, j ′′))) =

j′′∑

j′=1

max{0, di,j′}

≥

j∑

j′=1

di,j′ +

j′′∑

j′=j+2

di,j′ = ai,j + (ai,j′′ − ai,j+1)

> ai,j′′,

and hence α1(i, j′′) > ai,j′′. In particular, for any fixed row i there is at

most one column index j with (i, j) ∈ V2. In order to see that c(A) > 0implies V2 6= ∅ consider a feasible (0, 1)−path P with δ(P ) = c(A). If P

is a horizontal path without any row change then δ(P ) > 0 implies that P

contains an element of V1. Otherwise let ((i, j), (i′, j + 1)) be the first rowchange of P . Then by the feasibility of P , di,j < 0 and thus the subpath0, (i, 1), . . . , (i, j) contains an element of V1. In both cases the first vertex onP which is in V1 is in V2 as well. Note that the α1(i, j) ((i, j) ∈ [m] × [n])can be determined as follows.

for i = 1 to m do α1(i, 1) := ai,1

for j = 2 to n do

for i = 1 to m do α1(i, j) := max1≤i′≤m

α1(i′, j − 1) + δ((i′, j − 1), (i, j))

The α2(i, j) can be determined analogously by running through the matrixfrom right to left. Obviously, this gives a method to determine c(A) and theset V2 in time O(m2n). We denote the elements of V2 by

(i1, j1), (i2, j2), . . . , (it, jt),


such that i1 < i2 < · · · < it. A segment S (given by the parametersl1, l2, . . . , lm, r1, r2, . . . , rm) is constructed according to the following strat-egy. In row ik (k ∈ [t]) we choose the open part maximal under the conditionthat the right boundary is jk, i.e. we put

rik = jk and lik = max{j ≤ jk : aik ,j = 0} + 1.

In the remaining rows we choose the open part minimal (in a sense to bemade precise below) under the condition that the final result is a segment.The rows i < i1 and i > it remain closed. If lik > rik+1

+ 1 we choose theopen part in row ik + 1 maximal with rik+1 = lik − 1. If necessary we repeatthis step in the following rows, until finally li ≤ rik+1

+ 1 for some i withik < i < ik+1. If lik+1

> rik + 1 we proceed analogously, starting in rowik+1 − 1. For the details of the construction see Algorithm 1.

Example 3. Let

A =

0 0 0 0 0 0 0 5 90 0 0 0 1 1 2 4 20 0 2 2 3 3 3 2 11 1 2 2 1 1 1 1 11 3 4 2 2 2 4 4 72 2 2 2 1 2 2 3 30 2 2 7 2 2 2 1 1

.

Then c(A) = 9, V2 = {(1, 9), (5, 3), (7, 4)} and the algorithm yields the seg-ment

S =

0 0 0 0 0 0 0 1 10 0 0 0 1 1 1 0 00 0 1 1 0 0 0 0 00 0 0 0 0 0 0 0 01 1 1 0 0 0 0 0 00 0 0 0 0 0 0 0 00 1 1 1 0 0 0 0 0

,

where the bold 1’s correspond to the elements of V2. For the resulting matrix

A − S =

0 0 0 0 0 0 0 4 80 0 0 0 0 0 1 4 20 0 1 1 3 3 3 2 11 1 2 2 1 1 1 1 10 2 3 2 2 2 4 4 72 2 2 2 1 2 2 3 30 1 1 6 2 2 2 1 1

.

we have c(A − S) = 8.

Note that in Algorithm 1 each row i is considered exactly once and thatli and ri are defined either depending on some already defined values or bysearching for a zero entry running through the row from right to left. So thetotal time complexity of Algorithm 1 is at most O(mn). This is dominatedby the construction of V2, hence the complexity of the construction of S isO(m2n). To prove the correctness of the algorithm we need an alternative

description of paths in−→G that yields some insight into the relation between

the constructed segment S and the path weights. For this let−→H be a directed


Algorithm 1 Segment

Input: A = (ai,j) 1≤i≤m1≤j≤n

and V2 = {(i1, j1), . . . , (it, jt)}

Output: li, ri (i = 1, . . . , m)for (i, j) ∈ V2 do

li := max{j ′ ≤ j : ai,j′ = 0}+1ri := j

end for5: for i = 1 to i1 − 1 do

li := li1; ri := li − 1end forfor i = it + 1 to m do

li := lit; ri := li − 110: end for

for k = 1 to t − 1 doif jk > jk+1 then

i := ikwhile i < ik+1 and li > rik+1

+ 1 do15: i := i + 1

ri := li−1 − 1li := max{j ≤ ri : ai,j = 0} + 1

end whilefor i′ = i + 1 to ik+1 − 1 do

20: ri′ := rik+1; li′ := ri′ + 1

end forelse

i := ik+1

while i > ik and li > rik + 1 do25: i := i − 1

ri := li+1 − 1li := max{j ≤ ri : ai,j = 0} + 1

end whilefor i′ = ik + 1 to i − 1 do

30: ri′ := rik ; li′ := ri′ + 1end for

end ifend forreturn li, ri (i = 1, . . . , m)


graph with vertex set V ∪ {0, 1}. As the arc set of H we take E0 = E(1)0 ∪

E(2)0 ∪ E

(3)0 ∪ E

(4)0 , where

E(1)0 = {(0, (i, 1)) : i ∈ [m]} ∪ {((i, n), 1) : i ∈ [m]},

E(2)0 = {((i, j), (i, j + 1) : i ∈ [m], j ∈ [n − 1])},

E(3)0 = {((i, j), (i + 1, j)) : i ∈ [m − 1], j ∈ [n − 1]},

E(4)0 = {((i, j), (i − 1, j)) : 2 ≤ i ≤ m, j ∈ [n − 1]}.

Let the weight function δ0 on E0 be defined by

δ0(0, (i, 1)) = ai,1 (i ∈ [m]),

δ0((i, n), 1) = 0 (i ∈ [m]),

δ0((i, j), (i, j + 1)) = max{0, di,j+1} (i ∈ [m], j ∈ [n − 1]),

δ0((i, j), (i + 1, j)) = −ai,j (i ∈ [m − 1], j ∈ [n − 1]),

δ0((i, j), (i − 1, j)) = −ai,j (2 ≤ i ≤ m, j ∈ [n − 1]).

Example 4. Figure 2.18 shows−→H corresponding to the matrix

A =

(

4 5 0 1 4 52 4 1 3 1 42 3 2 1 2 45 3 3 2 5 3

)

.

· 1 //

−4��

· 0 //

−5

��

· 1 //

0

��

· 3 //

−1

��

· 1 //

−4

��

·

0

1��444

4444

4444

4

· 2 //

−2

��

−2YY

· 0 //

−4

��

−4YY

· 2 //

−1

��

−1YY

· 0 //

−3

��

−3YY

· 3 //

−1

��

−1

YY

· 0

''OOOOOOO

·

4

0

;;wwwwwwwwwwwwwwwww2

33gggggggggggggg2

++WWWWWWWWWWWWWW

5

##GGGGG

GGGGGG

GGGGG

G ·

· 1 //

−2

��

−2YY

· 0 //

−3

��

−3YY

· 0 //

−2

��

−2YY

· 1 //

−1

��

−1YY

· 2 //

−2

��

−2

YY

·0

77ooooooo

· 0 //

−5YY

· 0 //

−3YY

· 0 //

−3YY

· 3 //

−2YY

· 0 //

−5

YY

·

0

DD

Fig. 2.18: The digraph corresponding to the matrix A in Example 4.

It is easy to see that there is a bijection between the paths in−→G and the

paths in−→H with the additional restriction that the last arc is in E

(1)0 ∪ E

(2)0 .

In addition this bijection preserves the weight, that is for a path P in−→G and

the corresponding path Q in−→H we have

δ(P ) = δ0(Q).


In particular, there is a weight–preserving bijection between the (0, 1)−paths

in−→G and

−→H . The advantage of Q compared to P is that possibly existing

”long, skew” arcs in P are replaced by a sequence of vertical arcs and onehorizontal arc, and the weights of these arcs are easier to control. Analogousto α, α1 and α2 we define for (i, j) ∈ V ,

β1(i, j) = max{δ0(Q) : Q (0, (i, j)) − path in−→H},

β2(i, j) = max{δ0(Q) : Q ((i, j), 1) − path in−→H},

β(i, j) = β1(i, j) + β2(i, j).

We need some information about the connection between the weights in−→G

and−→H . Obviously β2(i, j) = α2(i, j) for all (i, j) ∈ V . The next lemma is an

analogous result about α1 and β1 for the vertices on (0, 1)−paths of maximalweight.

Lemma 8. For all (i, j) ∈ V with α(i, j) = c(A), we have β1(i, j) = α1(i, j).

Proof. β1(i, j) ≥ α1(i, j) is trivial, since for every (0, (i, j))−path in−→G there

is the corresponding (0, (i, j))−path in−→H of the same weight. Let Q1 and

Q2 be a (0, (i, j))−path in−→H and an ((i, j), 1)−path in

−→H , respectively, with

δ0(Q1) = β1(i, j) and δ0(Q2) = β2(i, j) = α2(i, j).

By concatenating Q1 and Q2 we obtain a (0, 1)−path Q in−→H with

δ0(Q) = β1(i, j) + α2(i, j).

Since the last arc of Q is in E(1)0 , this implies the existence of a (0, 1)−path

P in−→G with δ(P ) = β1(i, j) + α2(i, j). So

β1(i, j) + α2(i, j) ≤ c(A) = α1(i, j) + α2(i, j),

and thus β1(i, j) ≤ α1(i, j). �

Lemma 9. Let (i, j), (k, l) ∈ V , i > k and put p = i − k.

1. If j < l and there are column indices j ′1, j′2, . . . , j

′p such that

j ≤ j ′1 ≤ j ′2 ≤ · · · ≤ j ′p < l and

ai−q,j′q= 0 for q = 1, 2, . . . , p,

then there exists a ((i, j), (k, l))−path P in−→G with

δ(P ) ≥ ak,l − ai,j.


j ′1 = j ′2 j ′3 j ′4...

.... . . ∗ ∗ ∗ ∗ 0 ak,l . . .

∗ ∗ 0 ∗ ∗ ∗

∗ 0 ∗ ∗ ∗ ∗

∗ 0 ∗ ∗ ∗ ∗. . . ai,j ∗ ∗ ∗ ∗ ∗ . . .

......

· //· //·(k,l)

//·

· //·

OO

·

OO

· //·

OO

· (i,j)

OO

Fig. 2.19: A part of a matrix A and the corresponding path Q.

2. If j > l and there are column indices j ′1, j′2, . . . , j

′p such that

l ≤ j ′1 ≤ j ′2 ≤ · · · ≤ j ′p < j and

ak+q,j′q= 0 for q = 1, 2, . . . , p,

then there exists a ((k, l), (i, j))−path P in−→G with

δ(P ) ≥ ai,j − ak,l.

Proof. We consider only the first case that is illustrated in Fig. 2.19. Thesecond one is treated analogously. First we construct a ((i, j), (k, l))−path

Q in−→H . We take ((i, j), (i − 1, j)) with weight −ai,j as the first arc and

complete this arc to a ((i, j), (k, l))−path Q in such a way, that row changesoccur only along the arcs

((i − q, j ′q), (i − q − 1, j ′q)) (1 ≤ q ≤ p − 1).

This is possible by our assumption on the j ′q. Thus the vertical arcs of Q,except for the first one, have weight 0 and since the horizontal arcs havenonnegative weight in any case we conclude that the ((i, j), (k, j ′p))−subpathof Q has weight at least −ai,j. Finally the weight of the path

(k, j ′p), (k, j ′p + 1), . . . , (k, l)

is at least ak,l and from l > j ′p follows that the last arc of Q is in E(2)0 and

thus there exists a ((i, j), (k, l))−path P in−→G with

δ(P ) = δ(Q) ≥ ak,l − ai,j.

�


Lemma 10. Algorithm 1 yields a segment S.

Proof. Suppose the algorithm does not yield a segment. This is possible onlyif for some k ∈ [t − 1] the condition of the while–loop in line 14 (resp. line24) holds for all i ∈ {ik, ik + 1, . . . , ik+1 − 1} (resp. for all i ∈ {ik + 1, ik +2, . . . , ik+1}). If jk = jk+1 then

lik ≤ rik+1and lik+1

≤ rik .

So we may assume jk 6= jk+1. Let jk > jk+1. (The case jk < jk+1 can betreated analogously.) We put p = ik+1 − ik and

j ′q = lik+1−q − 1 (q = 1, 2, . . . , p).

The assumption that the while–condition is fulfilled for all ik+1 − q (q =1, 2, . . . , p) implies

rik+1+ 1 ≤ j ′1 ≤ j ′2 ≤ · · · ≤ j ′p < jk and

aik+1−q,j′q= 0 (q = 1, 2, . . . , p).

Thus by Lemma 9 there is a ((ik+1, jk+1 + 1), (ik, jk))−path P0 in−→G of

weight at least aik,jk− aik+1,jk+1+1. Using (ik+1, jk+1) ∈ V2, i.e. aik+1,jk+1

>

aik+1,jk+1+1, this yields

δ0(P0) > aik ,jk− aik+1,jk+1

.

Now we concatenate the path 0, (ik+1, 1), (ik+1, 2), . . . , (ik+1, jk+1 + 1) withP0 to obtain a (0, (ik, jk))−path of weight at least

aik+1,jk+1+ δ(P0) > aik ,jk

,

in contradiction to (ik, jk) ∈ V2. �

Let S = (si,j) be the result of Algorithm 1. By construction, si,j = 1implies ai,j ≥ 1 and so the entries of A − S are nonnegative. We put

a′i,j = ai,j − si,j (i ∈ [m], j ∈ [n]),

a′i,0 = ai,n+1 = 0 (i ∈ [m]),

d′i,j = a′

i,j − a′i,j−1 (i ∈ [m], j ∈ [n]).


By δ′ and δ′0 we denote the weight functions on−→G and

−→H , respectively, which

correspond to A′ = (a′i,j). For (i, j) ∈ V we put

α′1(i, j) = max{δ′(P ) : P (0, (i, j)) − path in

−→G},

α′2(i, j) = max{δ′(P ) : P ((i, j), 1) − path in

−→G},

α′(i, j) = α′1(i, j) + α′

2(i, j),

β ′1(i, j) = max{δ′(Q) : Q (0, (i, j)) − path in

−→H},

β ′2(i, j) = max{δ′(Q) : Q ((i, j), 1) − path in

−→H},

β ′(i, j) = β ′1(i, j) + β ′

2(i, j).

By T we denote the subset of V which corresponds to the segment S, that is

T = {(i, j) ∈ V : si,j = 1}.

The next lemma asserts that for (i, j) ∈ T the sequence ai,1, . . . , ai,j is non-decreasing and the horizontal path from 0 to (i, j) has maximal weight with

respect to A in both of−→G and

−→H .

Lemma 11. For (i, j) ∈ T we have

β1(i, j) = α1(i, j) = ai,j and α(i, j) = c(A).

Proof. Let (i, j) ∈ T . Clearly, β1(i, j) ≥ α1(i, j) ≥ ai,j. Assume P0 is a

(0, (i, j))−path in−→G with δ(P0) > ai,j. Recall that V2 = {(i1, j1), . . . , (it, jt)}

with i1 < i2 < · · · < it. We claim that for some k ∈ [t] there is an

((i, j), (ik, jk))−path P1 in−→G of weight at least aik ,jk

− ai,j. To see thiswe distinguish three types of vertices in T :

1. i = ik and j ≤ jk for some k ∈ [t]:

The path (ik, j), (ik, j + 1), . . . , (ik, jk) has weight aik,jk− aik,j.

2. ik < i < ik+1 for some k ∈ [t − 1] with jk > jk+1:

By construction of S there are column indices j ′1, j′2, . . . , j

′p, where p =

i − ik, such that

j ≤ j ′1 ≤ j ′2 ≤ · · · ≤ j ′p < jk and

ai−q,j′q= 0 (q = 1, 2, . . . , p).

Thus the claim follows by Lemma 9.


3. ik−1 < i < ik for some k ∈ {2, 3, . . . , t} with jk−1 < jk:

By construction of S there are column indices j ′1, j′2, . . . , j

′p, where p =

ik − i, such that

j ≤ j ′1 ≤ j ′2 ≤ · · · ≤ j ′p < jk and

ai+q,j′q= 0 (q = 1, 2, . . . , p).

Thus the claim follows by Lemma 9.

But now we can concatenate P0 and P1 to obtain a (0, (ik, jk))−path P in−→G

withδ(P ) = δ(P0) + δ(P1) > ai,j + (aik,jk

− ai,j) = aik,jk

in contradiction to (ik, jk) ∈ V2. This proves α1(i, j) = ai,j. In addition,concatenating the paths (0, (i, 1), (i, 2), . . . , (i, j)), P1 and a ((ik, jk), 1)−pathof maximal weight yields α(i, j) = c(A) and thus also β1(i, j) = α1(i, j) byLemma 8. �

Now we want to prove that for (i, j) ∈ T the horizontal (0, (i, j))−path isstill maximal with respect to A′. We need the following necessary conditionfor β1(i, j) > ai,j.

Lemma 12. Suppose β1(i, j) > ai,j and Q is a (0, (i, j))−path in−→H with

δ0(Q) = β1(i, j). Then there exists a vertex (i′, j ′) ∈ V1 such that either

• j ′ = 1 and ((i′, 1), (i′, 2)) is an arc of Q or

• 1 < j ′ < n and ((i′, j ′ − 1), (i′, j ′)), ((i′, j ′), (i′, j ′ + 1) are arcs of Q.

If in addition β(i, j) = c(A) then we can choose (i′, j ′) even in V2.

Proof. Let Q be a (0, (i, j))−path with δ0(Q) = β1(i, j) and assume thereis no such vertex in V1. We show δ0(Q) = ai,j which gives the desiredcontradiction. Clearly, δ0(Q) ≥ ai,j. The first arc of Q is of the form (0, (i′, 1))and has weight ai′,1. So we may assume that Q has more than one arc andproceed by induction on the number of arcs of an initial subpath of Q.

Case 1: The last arc of Q is in E(3)0 ∪ E

(4)0 .

W.l.o.g. the last arc is ((i − 1, j), (i, j)) with weight −ai−1,j . Since byinduction δ0(Q \ {(i, j)}) = ai−1,j, we obtain δ0(Q) = 0 ≤ ai,j.

Case 2: The last arc of Q is in E(2)0 , and the second last arc is in E

(3)0 ∪E

(4)0 .

W.l.o.g. the last two arcs of Q are ((i− 1, j − 1), (i, j − 1)) and ((i, j −1), (i, j)). By induction the weight of the (0, (i− 1, j − 1))−subpath ofQ is ai−1,j−1. Thus the weight of the (0, (i, j − 1))−subpath is 0 andby maximality of Q follows ai,j−1 = 0, hence δ0(Q) = ai,j.


Case 3: The last two arcs of Q are in E(1)0 ∪ E

(2)0 .

By induction the (0, (i, j − 1))−subpath of Q has weight ai,j−1. Bymaximality of Q this implies di,j′ ≥ 0 for all j ′, 1 ≤ j ′ ≤ j − 1. Nowdi,j ≥ 0, since otherwise (i, j − 1) is a vertex in V1 that fulfills theconditions of the lemma. Thus δ0(Q) = ai,j.

Now suppose β(i, j) = c(A). Then we can complete Q to a (0, 1)−path Q′ of

weight c(A). Let P be the corresponding (0, 1)−path in−→G , and let (i′, j ′) ∈

V1 be the first vertex on Q that has the claimed properties. Then (i′, j ′) ∈ P

and the (0, (i′, j ′))−subpath of P has weight ai′,j′, that is (i′, j ′) ∈ V2. �

Lemma 13. For (i, j) ∈ T we have β ′1(i, j) = α′

1(i, j) = a′ij.

Proof. Again trivially,

β ′1(i, j) ≥ α′

1(i, j) ≥ a′i,j.

Let (i, j) ∈ T and assume β ′1(i, j) > a′

i,j. In particular, j > 1 since obviously

β ′1(i, 1) = a′

i,1 for all i ∈ [m]. There is a (0, (i, j))−path Q in−→H with

δ′0(Q) = β ′1(i, j) > a′

i,j.

W.l.o.g. we may assume that (i, j) is the first counterexample to the lemmaon Q, i.e.

β ′1(i0, j0) = a′

i0,j0for all (i0, j0) ∈ (Q \ {(i, j)}) ∩ T.

Case 1: (Q \ {(i, j)}) ∩ T = ∅.

Let e be the last arc of Q. Then δ0(e1) = δ′0(e1) for all arcs e1 6= e ofQ.

Case 1.1: e ∈ E(2)0 .

Then δ0(e) = δ′0(e)+1, hence δ0(Q) = δ′0(Q)+1, and consequently(using Lemma 11),

β ′1(i, j) = δ0(Q) − 1 ≤ β1(i, j) − 1 = a′

i,j.

Case 1.2: e ∈ E(3)0 ∪ E

(4)0 .

W.l.o.g. e = ((i − 1, j), (i, j)) and δ0(e) = δ′0(e) = −ai−1,j , andthus

δ0(Q) = δ′0(Q) = β ′1(i, j).


Assume δ0(Q) = β1(i, j) = ai,j. Then δ0(Q) > 0, and thus

δ0(Q \ {(i, j)}) > ai−1,j .

By Lemma 11, β(i, j) = α(i, j) = c(A) and consequently byLemma 12, Q \ {(i, j)} contains a vertex (i0, j0) ∈ V2 ⊆ T . Thisis a contradiction and we conclude

β ′1(i, j) = δ0(Q) < β1(i, j) = ai,j,

and thus β ′1(i, j) = a′

i,j.

Case 2: (Q \ {(i, j)}) ∩ T 6= ∅.

Let (i0, j0) be the last vertex on Q \ {(i, j)} that is in T and denoteby Q1 and Q2 the (0, (i0, j0))−subpath and the ((i0, j0), (i, j))−subpathof Q, respectively. By assumption δ′0(Q1) = a′

i0,j0, so w.l.o.g. we may

assume Q1 = (0, (i0, 1), (i0, 2), . . . , (i0, j0)), and then

δ0(Q1) = β1(i0, j0) = ai0,j0 = δ′0(Q1) + 1. (2.7)

We denote the arcs of Q2 by e1, e2, . . . , ep. For p = 1 we obtain

δ′0(Q) = δ′0(Q1) − a′i0,j0

= 0 if e1 ∈ E(3)0 ∪ E

(4)0 and

δ′0(Q) = δ′0(Q1) + max{0, d′i,j} if e1 ∈ E

(2)0 .

Since e1 ∈ E(2)0 implies (i, j), (i, j − 1) ∈ T and thus d′

ij = di,j ≥ 0(Lemma 11), we obtain δ′0(Q) ≤ a′

i,j and consequently β ′1(i, j) = a′

i,j.So let p > 1. Then

δ′0(ei) = δ0(ei) (2 ≤ i ≤ p − 1), (2.8)

δ′0(e1) =

δ0(e1) + 1 if e1 ∈ E(3)0 ∪ E

(4)0 ,

δ0(e1) + 1 if e1 ∈ E(2)0 and di0,j0+1 ≥ 0,

δ0(e1) if e1 ∈ E(2)0 and di0,j0+1 < 0 and

(2.9)

δ′0(ep) =

{

δ0(ep) if ep ∈ E(3)0 ∪ E

(4)0 ,

δ0(ep) − 1 if ep ∈ E(2)0 ,

(2.10)

and in particular,

δ′0(Q2) ≤ δ0(Q2) + 1.


·(i0,j0)//·

e1��·

AA�� · //·(i1,j1)

//·ep

(i,j)��·

a) e1 ∈ E(3)0 ∪ E

(4)0

·(i0,j0)

//·e1 //·

��·

AA�� ·(i1,j1)//· //·

ep

(i,j)��·

b) e1 ∈ E(2)0 and di0,j0+1 ≥ 0

Fig. 2.20: Paths Q as in Case 2.2. of Lemma 13.

Case 2.1: δ′0(Q2) ≤ δ0(Q2).

δ′0(Q) = δ′0(Q1) + δ′0(Q2) ≤ a′i0,j0

+ δ0(Q2) < δ0(Q) implies

β ′1(i, j) < β1(i, j) = ai,j,

and thus β ′1(i, j) = a′

i,j.

Case 2.2: δ′0(Q2) = δ0(Q2) + 1.

In this case (2.7) and (2.8)–(2.10) imply

δ′0(Q) = δ0(Q) and ep ∈ E(3)0 ∪ E

(4)0 ,

w.l.o.g. ep = ((i − 1, j), (i, j)) with weight −ai−1,j . Assumeδ0(Q) = β1(i, j). Then δ0(Q) > 0 and thus

β1(i − 1, j) > ai−1,j .

By Lemma 11, β(i, j) = α(i, j) = c(A), and by Lemma 12 there isa vertex (i1, j1) ∈ V2 such that Q contains the arc ((i1, j1), (i1, j1 +1)). From (2.8)–(2.10) it follows that δ′0(Q2) = δ0(Q2) + 1 ispossible only if

e1 ∈ E(3)0 ∪ E

(4)0 or (e1 ∈ E

(2)0 and di0,j0+1 ≥ 0).

Hence, using di0,j′ ≥ 0 for 1 ≤ j ′ ≤ j0, (i1, j1) 6∈ Q1 and we obtainthe contradiction

(i1, j1) ∈ (Q \ {(i, j), (i0, j0)}) ∩ V2.

Thus δ′0(Q) = δ0(Q) < β1(i, j) = ai,j, and so β ′1(i, j) = a′

i,j. �

Now we are prepared for the final step.

Lemma 14. c(A′) ≤ c(A) − 1.


Proof. Let Q be a (0, 1)−path in−→H with δ′0(Q) = c(A′) and let (i0, j0) be

the last vertex on Q that is in T . We denote the (0, (i0, j0))−subpath andthe ((i0, j0), 1)−subpath of Q by Q1 and Q2, respectively. By Lemmas 11and 13,

β1(i0, j0) = ai0,j0 = a′i0,j0

+ 1 = β ′1(i0, j0) + 1,

and w.l.o.g. we may assume Q1 = (0, (i0, 1), (i0, 2), . . . , (i0, j0)). For the firstarc e0 of Q2 we have δ0(e0) = δ′0(e0) or δ0(e0) = δ′0(e0) − 1, and for all arcse 6= e0 of Q2, δ0(e) = δ′0(e).

Case 1: δ0(e0) = δ′0(e0).

δ0(Q) = δ0(Q1) + δ0(Q2) = δ′0(Q1) + 1 + δ′0(Q2)

= δ′0(Q) + 1 = c(A′) + 1,

and thus c(A) ≥ c(A′) + 1.

Case 2: δ0(e0) = δ′0(e0) − 1.

By the same argument as in the first case we only get

δ0(Q) = c(A′).

Assume δ0(Q2) = α2(i0, j0). From

α(i0, j0) = c(A) and α1(i0, j0) = ai0,j0

we deduce δ0(Q) = c(A). Now consider two cases:

1. If Q has a vertex (i, j) with β1(i, j) > ai,j, then by Lemma 12, Q

contains an arc ((i1, j1), (i1, j1 + 1)) with (i1, j1) ∈ V2.

2. If β1(i, j) = ai,j for every (i, j) ∈ Q, let (i1, j1) be the second lastvertex of Q, i.e. j1 = n and ((i1, j1), 1) is the last arc of Q. Notethat β1(i1, j1) = ai1,j1 implies (i1, j1) ∈ V2.

From δ0(e0) = δ′0(e0) − 1 follows that either

e0 ∈ E(3)0 ∪ E

(4)0 or (e0 ∈ E

(2)0 and di0,j0+1 ≥ 0).

Hence, using di0,j′ ≥ 0 for 1 ≤ j ′ ≤ j0, (i1, j1) 6∈ Q1 and we obtain thecontradiction

(Q2 \ {(i0, j0)}) ∩ V2 6= ∅.

Consequently, δ0(Q2) < α2(i0, j0) and there exists an ((i0, j0), 1)−pathQ∗

2 with δ0(Q∗2) > δ0(Q2). By concatenating Q1 and Q∗

2 we obtain a(0, 1)−path Q∗ with δ0(Q

∗) > c(A′), and thus

c(A) ≥ c(A′) + 1.

�

2.4. Test results 47

Now we collect the lemmas to prove Theorem 1.

Proof of Theorem 1. That the maximal weight of a path is a lower bound forthe TNMU is an immediate consequence of Lemmas 4, 6 and 7 and duality.The existence of a segmentation with

∑k

i=1 ui = c(A) is proved by inductionon c(A). If c(A) = 0 then A = 0 and there is nothing to do. For c(A) > 0 weapply Algorithm 1 to construct a segment S with c(A − S) ≤ c(A) − 1. Byinduction there are segments S2, S3, . . . , Sk and positive integers u2, u3, . . . , uk

such that

A − S =

k∑

i=2

uiSi and

k∑

i=2

ui = c(A − S) ≤ c(A) − 1,

and thus with S1 = S and u1 = 1,

A =k∑

i=1

uiSi andk∑

i=1

ui = c(A − S) + 1 ≤ c(A).

�

As observed after Example 3 each segment can be determined in timeO(m2n), so the time needed for the whole segmentation is bounded by

c(A)O(m2n). On any 0 − 1−path in−→H the number of arcs with positive

weight is bounded by n, because only the horizontal arcs can have positiveweight. And for each of these arcs the weight is bounded by

L = max{ai,j : i = 1, . . . , m; j = 1, . . . , n}.

So c(A) ≤ nL and the complexity of the whole segmentation algorithm isO(m2n2L).

2.4 Test results

Table 2.1 shows some test results of our algorithm in comparison with otheralgorithms. Each row shows the average TNMU for a 15 × 15–matrix withrandomly chosen entries from {0, . . . , L}. The columns labeled ’Xia–Verhey’,’Bortfeld’ and ’Galvin’ contain the results for the algorithms of Xia andVerhey [28], Bortfeld et al [6] and Galvin et al [11], respectively. The numbersin these columns are taken from Xia and Verhey [28]. The last column showsthe average TNMU obtained by Engel’s algorithm [10], which is TNMU–optimal for the segmentation problem without ICC. To obtain the resultsof the column labeled ’new’ we implemented Algorithm 1 in C++. For a


L new Xia–Verhey

Bortfeld Galvin Engel (with-out ICC)

3 15.4 19.5 17.7 19.7 14.04 19.5 29.6 22.8 40.5 17.95 23.6 30.9 27.9 40.1 21.76 27.6 46.8 32.8 44.2 25.67 31.7 45.6 37.9 67.1 29.48 35.7 63.4 42.8 72.3 33.29 39.8 67.1 47.8 72.3 37.010 43.8 68.6 52.6 76.5 40.911 47.7 68.6 57.6 81.4 44.712 51.8 101.1 62.4 106.8 48.513 55.7 100.6 67.3 101.1 52.314 59.8 100.0 72.2 112.7 56.215 63.8 98.0 77.1 116.0 59.816 67.7 124.9 82.0 154.5 63.3

Tab. 2.1: Average TNMU for random 15 × 15–matrices with maximal entry L.

matrix A the segment S was determined and subtracted from A, and thiswas iterated until the zero matrix was reached. For each L this was donefor 10000 random matrices A and the average TNMU was determined. Ona 1.3 GHz PC the computation for the whole column (i.e. the segmentationof 140000 matrices) took 206 seconds.

3. A HEURISTIC FOR THE REDUCTION OF THENUMBER OF SEGMENTS

In this chapter we present a greedy–heuristic that can be used to find asegmentation with minimal TNMU and a small NS. For brevity of notation

we slightly modify the digraph−→H from Chapter 2. We add vertices (0, i)

and (n + 1, i) ( for i ∈ [m]) and replace every arc (0, (i, 1)) by the two arcs(0, (i, 0)), ((i, 0), (i, 1)) and every arc ((i, n), 1) by the two arcs ((i, n), (i, n +1)), ((i, n + 1), 1). The weights of the new arcs are determined by

δ0(0, (i, 0)) = δ0((i, n), (i, n + 1)) = δ0((i, n + 1), 1) = 0,

δ0((i, 0), (i, 1)) = ai,1

for all i ∈ [m]. The resulting digraph is called−→H again. Figure 3.1 illustrates

the modification for m = n = 4. By the results of Chapter 2, we may assume

· //

��

· //

��

· //

��

·0

1��999

9999

9999

· //

��

[[

· //

��

[[

· //

��

[[

· 0

((RRRRRRR

·

a1,1

0

<<zzzzzzzzzzzzz

a2,1 33hhhhhhhhhh a3,1

++VVVVVVVVVV

a4,1

""DDDD

DDDD

DDDD

D ·· //

��

[[

· //

��

[[

· //

��

[[

·0

66lllllll

· //

[[

· //

[[

· //

[[

·

0

BB��

·a1,1 //· //

��

· //

��

· //

��

· 0 //·0

1��;;;

;;;;

;;;

·a2,1 //· //

��

]]

· //

��

]]

· //

��

]]

· 0 //· 0

))SSSSSSS

·

0

0

AA��0 55kkkkkkk 0

))SSSSSSS

0��;

;;;;

;;;;

; ··

a3,1 //· //

��

]]

· //

��

]]

· //

��

]]

· 0 //·0

55kkkkkkk

·a4,1 //· //

]]

· //

]]

· //

]]

· 0 //·

0

AA��

Fig. 3.1: The old and the new digraph.

that we have already determined the minimal TNMU which equals

c(A) = max{δ0(P ) : P is a (0, 1) − path in−→H},

and for every (i, j) ∈ [m] × [n] the values

β1(i, j) = max{δ0(P ) : P is a (0, (i, j)) − path in−→H}, (3.1)

β2(i, j) = max{δ0(P ) : P is a ((i, j), 1) − path in−→H}. (3.2)

50 3. A heuristic for the reduction of the number of segments

3.1 The algorithm

Adopting the terminology of [10] we call the pair (u, S) of a positive integeru and a segment S an admissible segmentation pair if

A′ = A − uS is nonnegative and

c(A′) = c(A) − u.

The essential step of our algorithm is to determine the maximal coefficientu with the property that there exists a segment S, such that (u, S) is anadmissible segmentation pair. Iterating this step with A′ = A−uS we clearlyobtain a segmentation of A with c(A) monitor units. In order to derive anupper bound for the coefficient u in an admissible segmentation pair (u, S),we identify, according to [4], the set of segments with the set of paths fromD to D′ in the layered digraph Γ = (W, F ), constructed as follows. Thevertices in the i−th layer correspond to the possible leaf positions in row i

(1 ≤ i ≤ m) and two additional vertices D and D′ are added:

W = {(i, l, r) : i = 1, . . . , m, l = 1, . . . , n + 1, r = l − 1, . . . , n} ∪ {D, D′}.

Between two vertices (i, l, r) and (i+1, l′, r′) there is an arc if the correspond-ing leaf positions are consistent with the ICC, i.e. if l′ ≤ r +1 and r′ ≥ l− 1.In addition, the arc set F contains all arcs from D to the first layer and fromthe last layer m to D′, so

F = F+(D) ∪ F−(D′) ∪m−1⋃

i=1

F+(i), where

F+(D) = {(D, (1, l, r)) : (1, l, r) ∈ W},

F−(D) = {((m, l, r), D′) : (m, l, r) ∈ W},

F+(i) = {((i, l, r), (i + 1, l′, r′)) : l′ ≤ r + 1, r′ ≥ l − 1}.

There is a bijection between the possible leaf positions and the paths fromD to D′ in Γ. This is illustrated in Fig. 3.2 which shows the paths in Γ form = 4, n = 2, corresponding to the segments

(

1 00 11 11 0

)

(straight lines) and

(

0 11 11 00 1

)

(dotted lines).

Assume, for every triple (i, l, r), 1 ≤ i ≤ m, 1 ≤ l ≤ r + 1 ≤ n + 1, wehave already determined some upper bound u0(i, l, r) for the coefficient u in


/. -,() *+D

vvlllllllllllll

((/. -,() *+110 /. -,() *+111

,,XXXXXXXXXXXXXXXXXXXXXXXXXXX /. -,() *+112 /. -,() *+121 /. -,() *+122

tt

/. -,() *+132

/. -,() *+210 /. -,() *+211 /. -,() *+212

yy

/. -,() *+221 /. -,() *+222

ttiiiiiiiiiiiiiiiii/. -,() *+232

/. -,() *+310 /. -,() *+311

,,

/. -,() *+312

yyssssss

ss/. -,() *+321 /. -,() *+322 /. -,() *+332

/. -,() *+410 /. -,() *+411

((RRRRRRRRRRRRR/. -,() *+412 /. -,() *+421 /. -,() *+422

vv

/. -,() *+432

/. -,() *+D′

Fig. 3.2: The vertices of Γ for m = 4, n = 2 and two (D,D ′)–paths.

an admissible segmentation pair (u, S), where S is a segment with li = l andri = r. In other words, u ≤ u0(i, li, ri) for all i if (u, S) is an admissiblesegmentation pair and li, ri (i = 1, . . . , m) are the parameters of S. We put

u = max{u : There is a path D, (1, l1, r1), . . . , (m, lm, rm), D′

in Γ with u0(i, li, ri) ≥ u for i = 1, . . . , m}.

Clearly, u is an upper bound for the coefficient u in an admissible segmen-tation pair (u, S). Now we describe an algorithm which constructs an ad-missible segmentation pair (u, S) with maximal u. Fix u and assume wehave already determined the first i − 1 rows of a segment. If it is possi-ble to complete these i − 1 rows to obtain a segment S such that (u, S) isan admissible segmentation pair, then procedure Complete Segment(i) (Al-gorithm 2) determines li, . . . , lm and ri, . . . , rm realizing such a completion.

Here MaxWeight(i) denotes the maximal weight of a path in−→H that has all

its vertices in the first i rows, where the arc weights are determined accord-ing to A − uS for any segment S with parameters l1, . . . , li and r1, . . . , ri inthe first i rows. Clearly MaxWeight(i) depends only on the values that arealready determined and the condition in line 4 is necessary for the possibilityto continue with the candidates li, ri and obtain an admissible segmentationpair. Now the pair (u, S) is constructed by procedure Construct Segment

(Algorithm 3). Clearly, the efficiency of the backtracking depends very muchon the quality of the bounds u0(i, l, r). We give some bounds that turnedout to be quite good in numerical experiments. Trivially, in an admissiblesegmentation pair (u, S) we have, for all i,

u ≤ v1(i, li, ri) := min{ai,j : li ≤ j ≤ ri}.


Algorithm 2 Complete Segment(i)

for (li, ri) with 1 ≤ li ≤ ri−1 + 1,max{li, li−1} − 1 ≤ ri ≤ n andu0(i, li, ri) ≥ u do

if MaxWeight(i) ≤ c(A) − u thenif i < m thenComplete Segment(i + 1)

elsefinished:=true

end ifend if

end for

Algorithm 3 Construct Segment

u := u

finished:=falsel0 := 1, r0 := n + 1while not finished doComplete Segment(1)if not finished then

u := u − 1end if

end while


Fix an admissible segmentation pair (u, S), denote by δ′0 the weight function

on−→H corresponding to A′ = A − uS and let

β ′1(i, j) = max{δ′0(P ) : P is a (0, (i, j)) − path in

−→H}, (3.3)

β ′2(i, j) = max{δ′0(P ) : P is a ((i, j), 1) − path in

−→H}. (3.4)

The upper bounds below are based on the following simple observations.

1. The only arcs e with δ′0(e) < δ0(e) are of the form e = ((i, li − 1), (i, li))(1 ≤ i ≤ m), and for these arcs δ′0(e) ≥ δ0(e) − u.

2. For arcs of the form e = ((i, j), (i ± 1, j)) with li ≤ j ≤ ri we haveδ′0(e) = δ0(e) + u.

3. If j < lk for some k ∈ [m] then, on every (0, (k, j))−path P , the numberof arcs of the form ((i, li − 1), (i, li)) is equal to or less than the numberof arcs of the form ((i, j), (i ± 1, j)) with li ≤ j ≤ ri.

4. If j ≥ lk for some k ∈ [m] then, on every ((k, j), 1)−path P , the numberof arcs of the form ((i, li − 1), (i, li)) is equal to or less than the numberof arcs of the form ((i, j), (i ± 1, j)) with li ≤ j ≤ ri.

The third observation is valid since for a fixed (0, (k, j))−path P there is aninjective mapping from the set of arcs of the form ((i, li − 1), (i, li)) on P

to the set of arcs of the form ((i, j), (i ± 1, j)) with li ≤ j ≤ ri on P : anarc ((i, li − 1), (i, li)) is mapped to the arc ((i′ ± 1, j ′), (i′, j ′)) where (i′, j ′)is the first vertex on the ((i, li), (k, j))−subpath of P which is covered bysome left leaf. Here the ICC assures that li′±1 ≤ j ′ ≤ ri′±1. Similarly, inthe fourth observation the arcs of the form ((i, li − 1), (i, li)) can be mappedinjectively to the arcs of the form ((i, j), (i ± 1, j)) with li ≤ j ≤ ri bymapping ((i, li − 1), (i, li)) to ((i′ ± 1, j ′), (i′, j ′)) where (i′ ± 1, j ′) is the lastvertex on the ((k, j), (i, li))−subpath of P which is not covered by some leftleaf. This is illustrated in Figure 3.3. It follows, for 1 ≤ i ≤ m,

β ′1(i, j) ≥ β1(i, j) for j < li,

β ′2(i, j) ≥ β2(i, j) for j ≥ li.

Lemma 15. Let (u, S) be an admissible segmentation pair with li = l andri = r. Then u ≤ v2(i, l, r) where

v2(i, l, l − 1) = c(A) − β1(i, l − 1) − max{0, di,l} − β2(i, l),


-5 -4 -3 -2 -1 0 1 2 3 4 50

1

2

3

4

5

6

�

� � � �

�

� �

� �

� �

�

�

� � � �

� �

�

Fig. 3.3: Illustration of observations 3 and 4.

and if r ≥ l then v2(i, l, r) = min{γ1, γ2, γ3, γ4}, where

γ1 = c(A) − β1(i, l − 1) − β2(i, l),

γ2 = c(A) − β1(i, l − 1) −r∑

j=l+1

max{0, di,j} − β2(i, r + 1),

γ3 = c(A) − β1(i, l − 1) − di,l −r∑

j=l+1

max{0, di,j} − di,r+1 − β2(i, r + 1),

γ4 =1

2

(

c(A) − β1(i, l − 1) −r∑

j=l+1

max{0, di,j} − di,r+1 − β2(i, r + 1)

)

.

Proof. Let P be the the concatenation of the paths P1, P2 and P3, whereP1 is a (0, (i, l − 1))−path with δ0(P1) = β1(i, l − 1), P2 is the path ((i, l −1), (i, l), . . . , (i, r+1)), and P3 is an ((i, r+1), 1)−path with δ0(P2) = β2(i, r+1).

Case 1: r = l − 1. Using the above observations, we obtain

c(A) − u = c(A′) ≥ δ′0(P ) ≥ β1(i, l − 1) + max{0, di,l} + β2(i, l),

and thus u ≤ c(A) − β1(i, l − 1) − max{0, di,l} − β2(i, l).

Case 2: r ≥ l. Now

δ′0(P ) = δ′0(P1) + max{0, di,l − u} +r∑

j=l+1

max{0, di,j}

+ max{0, di,r+1 + u} + δ′0(P2),


and thus

β1(i, l − 1) + max{0, di,l − u} +r∑

j=l+1

max{0, di,j}

+ max{0, di,r+1 + u} + β2(i, r + 1) ≤ c(A) − u,

or

u + max{0, di,l − u} + max{0, di,r+1 + u} ≤ c(A) − β1(i, l − 1)

−r∑

j=l+1

max{0, di,j} − β2(i, r + 1),

which implies u ≤ γi (i = 2, 3, 4). To see u ≤ γ1, consider the pathQ that is the concatenation of P1, the arc ((i, l − 1), (i, l)) and an((i, l), 1)−path P4 with δ0(P4) = β2(i, l). Then

δ′0(Q) ≥ β1(i, l − 1) + β2(i, l),

and thus u ≤ γ1. �

Lemma 16. Suppose (u, S) is an admissible segmentation pair, fix some i,2 ≤ i ≤ m − 1, and put

λ1 = maxli≤t≤ri

{β1(i − 1, t) − ai−1,t − ai,t + β2(i + 1, t)},

λ2 = maxli≤t≤ri

{β1(i + 1, t) − ai+1,t − ai,t + β2(i − 1, t)}.

Thenu ≤ v3(i, li, ri) := c(A) − min{λ1, λ2}.

Proof. By symmetry, w.l.o.g. λ1 ≤ λ2. Assume u > c(A) − λ1, and let t

be an index where the maximum in the definition of λ1 is attained. LetP be the concatenation of the three paths P1, P2 and P3, where P1 is an(0, (i− 1, t))−path with δ0(P1) = β1(i− 1, t), P2 = ((i− 1, t), (i, t), (i + 1, t))and P3 is an ((i + 1, t), 1)–path with δ0(P3) = β2(i + 1, t). Then

δ′0(P ) ≤ c(A′) = c(A) − u < λ1 = δ0(P ).

By the above observations, we have δ′0(P1) ≥ δ0(P1)− u, δ′0(P3) ≥ δ0(P3)− u

and

δ′0(P2) =

{

δ0(P2) + 2u if li−1 ≤ t ≤ ri−1,

δ0(P2) + u otherwise.


So δ′0(P ) < δ0(P ) implies

δ′0(P1) < δ0(P1),

δ′0(P2) = δ0(P2) + u,

δ′0(P3) < δ0(P3).

And from this follows

li−1 ≤ t and li+1 > t.

Now denote by t′ the index where the maximum in the definition of λ2 isattained. Since u > c − λ1 ≥ c − λ2, by the same argument as above weobtain

li+1 ≤ t′ and li−1 > t′.

But this is a contradiction to li+1 > t if t′ ≤ t and to li−1 ≤ t if t′ > t. �

Thus we may put

u0(i, l, r) = min{vk(i, l, r) : k = 1, 2, 3}, (3.5)

and obtain the following result.

Theorem 2. If the u0(i, l, r) are determined according to (3.5) the algorithmConstruct Segment yields an admissible segmentation pair (u, S) such thatu′ ≤ u for any admissible segmentation pair (u′, S ′).

Example 5. For the benchmark matrix from [18] our algorithm yields thesegmentation (2.4) from Example 1.

3.2 Test results

To test our algorithm we computed segmentations for 15× 15–matrices withrandom entries from {0, 1, . . . , L} for 3 ≤ L ≤ 16. Table 3.1 shows theresults. The numbers in the columns TNMU (new) and NS (new) are theaverage total number of monitor units and the average number of segments,where we have averaged over 10000 matrices with randomly chosen entriesfrom {0, . . . , L} (uniformly distributed). The remaining columns show thecorresponding results from [28]: the columns labeled X–V, B, G contain theresults for the algorithms of Xia and Verhey [28], Bortfeld et al [6] and Galvinet al [11], respectively. On an 1.3 GHz–PC the computation of the two newentries in a row of the table, i.e. the segmentation of 10000 matrices, tookapproximately 1 hour. But it should be mentioned that the algorithm isfast for the vast majority of the matrices, while there are some very rareexceptions. We also tested the algorithm on 13 clinical sample matrices(10 × 10–matrices with entries between 0 and 10). The results are shown inTable 3.2.


L TNMU TNMU TNMU TNMU NS NS NS NS(new) (X–V) (B) (G) (new) (X–V) (B) (G)

3 15.4 19.5 17.7 19.7 12.6 13.3 17.7 13.44 19.5 29.6 22.8 40.5 14.5 18.6 22.8 20.45 23.6 30.9 27.9 40.1 16.0 19.0 27.9 20.46 27.6 46.8 32.8 44.2 17.2 20.3 32.8 21.57 31.7 45.6 37.9 67.1 18.2 20.0 37.9 27.18 35.7 63.4 42.8 72.3 19.1 24.3 42.8 28.29 39.8 67.1 47.8 72.3 19.9 24.3 47.8 28.310 43.8 68.6 52.6 76.5 20.7 25.7 52.6 28.911 47.7 68.6 57.6 81.4 21.3 25.7 57.6 30.912 51.8 101.1 62.4 106.8 21.9 27.0 62.4 34.813 55.7 100.6 67.3 101.1 22.5 26.9 67.3 35.514 59.8 100.0 72.2 112.7 23.0 26.9 72.2 35.615 63.8 98.0 77.1 116.0 23.5 26.7 77.1 35.916 67.7 124.9 82.0 154.5 24.0 30.0 82.0 41.7

Tab. 3.1: Average TNMU and NS for random 15 × 15–matrices with maximalentry L.

no. MU NS CPU–time1 18 10 0.05 s2 16 8 0.06 s3 16 8 0.05 s4 20 10 0.10 s5 19 11 0.05 s6 18 9 0.05 s7 17 9 0.11 s8 23 12 0.22 s9 24 11 0.17 s10 22 10 0.22 s11 30 15 0.17 s12 23 13 0.22 s13 22 11 0.22 s

Tab. 3.2: Test results for clinical matrices.

4. EXACT MINIMIZATION OF THE NUMBER OF

SEGMENTS FOR COLLIMATORS WITHOUT INTERLEAFCOLLISION CONSTRAINT

The problem of minimizing the number of segments is NP–complete in thestrong sense even for single row matrices. The NP–hardness was shown in [2]by reduction of 2–Partition [12]. Woeginger gave an unpublished proof ofthe NP–hardness in the strong sense by a reduction of 3–Partition [12]. In[15] the NS–minimization for one row has been reduced to the bipartite caseof Minimum Edge–Cost Flow [12]. For special instances of Minimum

Edge–Cost Flow there is a reduction in the reverse direction and thisyields a new point of view on Woegingers argument which is presented below.The following special case of Minimum Edge–Cost Flow has been shownto be strongly NP–complete in [3] by a reduction of 3–Partition.

Instance: A complete bipartite graph G = (U ∪ V, E) with |U | = 3|V | and afunction w : U → IN \ {0}.

Question: Is there a flow function f : E → IN such that

∀x ∈ U∑

y∈V

f(xy) = w(x), (4.1)

∀y ∈ V∑

x∈U

f(xy) = 3w, where w =1

|U |

∑

x∈U

w(x), (4.2)

|{xy ∈ E : f(xy) > 0}| ≤ |U |. (4.3)

This problem can be reduced to the NS–minimization problem as follows. Weput q = |V |, n = 4q, denote the elements of U by u1, . . . , u3q, the elementsof V by v1, . . . , vq, and define the row vector a = ( a1 ... an ) as follows.

ai =

i∑

j=1

w(uj) for 1 ≤ i ≤ 3q,

ai = 3(n − i)w for 3q + 1 ≤ i ≤ n.

60 4. Exact minimization of the number of segments

Theorem 3 (Woeginger). There is a segmentation of a with 3q segmentsiff there is a function f : E → IN satisfying (4.1)–(4.3).

Proof. “⇒”: Suppose there is a segmentation

a =

3q∑

j=1

cjs(j) (4.4)

where the segments are described by

s(j)i =

{

1 if lj ≤ i ≤ rj

0 otherwise.(j ∈ [3q], i ∈ [n]).

By Lemma 1 from [15] we may assume that lj ≤ 3q ≤ rj for all j ∈ [3q].Moreover, ai > ai−1 for all i ∈ [3q] (where a0 = 0) implies that for eachi ∈ [3q] there is some j ∈ [3q] with lj = i, hence we may assume li = i

for i ∈ [3q]. Let

f(uivri−3q+1) = ci (i ∈ [3q])

and f(xy) = 0 for all the remaining edges xy. Observe that an = 0, sori < n for all i and 1 ≤ ri − 3q + 1 ≤ q. Clearly, (4.3) is satisfied. Nowfix i ∈ [3q]. From (4.4) and the fact that j = i is the only index withlj = i we obtain that

w(ui) = ai − ai−1 = ci = f(uivri) =

∑

y∈V

f(uiy),

so (4.1) is satisfied. Now fix i, 3q + 1 ≤ i ≤ n. From (4.4) we obtain

3w = ai−1 − ai =∑

j∈[3q]:rj=i−1

cj

=∑

j∈[3q]:rj=i−1

f(ujvi−3q) =

3q∑

j=1

f(ujvi−3q),

thus (4.2) is satisfied.

“⇐”: Suppose there is a function f satisfying (4.1)–(4.3). By (4.1) and (4.3),for each j ∈ [3q] there is exactly one k(j) ∈ [q] with f(ujvk(j)) > 0. Forj ∈ [3q], put

cj = f(ujvk(j)), lj = j, rj = 3q + k(j) − 1,

4.1. Single row intensity maps 61

and define segments s(j) by

s(j)i =

{


0 otherwise.(j ∈ [3q], i ∈ [n]).

This yields a segmentation of a: for i ≤ 3q we have s(j)i = 1 iff lj ≤ i,

and so

3q∑

j=1

cjs(j)i =

i∑

j=1

cj =

i∑

j=1

∑

v∈V

f(uj, v) =

i∑

j=1

w(uj) = ai,

and for i > 3q we have s(j)i = 1 iff rj ≥ i, so

3q∑

j=1

cjs(j)i =

n−1∑

t=i

∑

j∈[3q]:rj=t

cj =

n−1∑

t=i

∑

j∈[3q]:k(j)=t−3q+1

cj

=n−1∑

t=i

3q∑

j=1

f(ujvt−3q+1) = (n − i)w = ai.

�

This shows that the NS–minimization is NP–hard. But the reductionessentially depends on the fact that the entries can become arbitrary large.In this chapter we show that the NS–minimization problem can be solved intime polynomial in the matrix dimensions m and n if the maximal entry L ofthe intensity matrix is bounded. This seems to be a reasonable assumption inpractice: for instance Xia and Verhey [28] report that they obtained matriceswith 7 nonzero intensity levels when they applied a preliminary version ofthe CORVUS inverse treatment planning system (NOMOS corporation) toa very complex head and neck tumor case. The algorithm proposed here isan application of the dynamic programming principle (see [5]).

4.1 Single row intensity maps

First we give an exact formulation of the problem L–One Row–Min MU–

Min NS:

Instance: A vector a = ( a1 a2 ... an ) of integers with 0 ≤ ai ≤ L (i = 1, . . . , n).

Problem: Find a segmentation with in first instance minimal TNMU and insecond instance minimal NS!


As before, we put a0 = an+1 = 0. Let

P = {i ∈ [n] : ai ≥ ai−1 and ai > ai+1},

Q = {i ∈ [n] : ai < ai−1 and ai ≤ ai+1}.

Clearly, |P | = |Q|+ 1 if an 6= 0 and |P | = |Q| if an = 0. If an 6= 0 denote theelements of P and Q by p1, . . . , pt and q1, . . . , qt−1 such that

p1 < q1 < p2 < q2 < · · · < qt−1 < pt,

and put q0 = 0 and qt = n + 1. If an = 0 denote the elements of P and Q byp1, . . . , pt and q1, . . . , qt such that

p1 < q1 < p2 < q2 < · · · < qt−1 < pt < qt.

From the results of [10] it follows that in a TNMU–optimal segmentation

a =

k∑

j=1

cjs(j)

every segment is of the form

s(j)i =

{

1 for lj ≤ i ≤ rj,

0 otherwise,

with qτ−1 < lj ≤ pτ and pτ ′ ≤ rj < qτ ′ for some τ, τ ′ ∈ [t]. Since theorder of the segments is not relevant, we may order them in such a way thatr1 ≤ · · · ≤ rk. For τ ∈ [t − 1], let k0(τ) be the unique index with rj < qτ forj ≤ k0(τ) and rj ≥ qτ for j > k0(τ), and put

a(τ) = a −

k0(τ)∑

j=1

cjs(j).

Also put k0(0) = 0, k0(t) = k, a(0) = a and a(t) = 0. For j > k0(τ), fromrj ≥ qτ it follows that for i ≤ qτ ,

s(j)i = 1 ⇐⇒ lj ≤ i.

In particular, for i = 1, . . . , qτ − 1 and j = k0(τ) + 1, . . . , k,

s(j)i = 1 =⇒ s

(j)i+1 = 1. (4.5)


For 0 ≤ τ ≤ t − 1, we have

a(τ) =

k∑

j=k0(τ)+1

cjs(j),

hence (4.5) implies that

a(τ)1 ≤ a

(τ)2 ≤ · · · ≤ a(τ)

qτ,

and the multisets

Uτ = {a(τ)i − a

(τ)i−1 : 1 ≤ i ≤ qτ , a

(τ)i 6= a

(τ)i−1}, (4.6)

Vτ = {a(τ)i − a

(τ)i−1 : qτ < i ≤ pτ+1, a

(τ)i 6= a

(τ)i−1}, (4.7)

Wτ = {a(τ)i − a

(τ)i+1 : pτ+1 ≤ i < qτ+1, a

(τ)i 6= a

(τ)i+1} (4.8)

are partitions of aqτ, apτ+1 − aqτ

and apτ+1 − aqτ+1, respectively. Observe that

a(τ)i = ai for i ≥ qτ , hence Vτ and Wτ depend only on a, while Uτ depends

also on the pairs(s(1), c1), . . . , (s

(k0(τ)), ck0(τ)).

Considering the sequence (Uτ , Vτ , Wτ ) (τ = 0, . . . , t), where we add Ut =Vt = Wt = ∅, we will derive a method to construct the desired segmentation.

Definition 2. For integers u, v and w with 0 ≤ u ≤ v ≤ L and 0 ≤ w < v,a (u, v, w)−peak is a triple (U, V, W ) of unordered partitions of u, v − u andv − w, i.e. a triple of multisets of positive integers with

∑

x∈U

x = u,∑

x∈V

x = v − u,∑

x∈W

x = v − w.

In addition, the triple (∅, ∅, ∅) is called (0, 0, 0)−peak.

Thus for τ = 0, . . . , t, (Uτ , Vτ , Wτ ) is an (aqτ, apτ+1, aqτ+1)−peak (where

apt+1 = aqt+1 = 0), and for τ ≤ t − 1, the choice of the pairs

(s(k0(τ)+1), ck0(τ)+1), . . . , (s(k0(τ+1)), ck0(τ+1))

can be considered as the choice of a way to go from the peak (Uτ , Vτ , Wτ ) tothe peak (Uτ+1, Vτ+1, Wτ+1). We claim that the number of segments neededfor this step does not depend on the particular a(τ), but only on the multisetsUτ ∪ Vτ , Wτ and Uτ+1. To prove this we associate with a (u, v, w)–peak(U, V, W ) a vector b = ( b1 ... bβ ) as follows. Put α = |U | + |V |, β = α +


|W |, denote the elements of U ∪ V by d1, . . . , dα and the elements of W bydα+1, . . . , dβ, such that

d1 ≥ d2 ≥ · · · ≥ dα and dα+1 ≥ dα+2 ≥ · · · ≥ dβ.

So, for U = Uτ , V = Vτ and W = Wτ the di (i = 1, . . . , β) are the absolutevalues of the nonzero differences of consecutive entries of the initial part( a

(τ)1 ... a

(τ)qτ+1 ) of a(τ). Now b is defined by

bi =

i∑

j=1

dj for 1 ≤ i ≤ α,

v −i∑

j=α+1

dj for α + 1 ≤ i ≤ β.

In addition, let b0 = 0.

Example 6. The associated vector for any peak with U ∪ V = {4, 2, 1, 1}and W = {2, 2, 1} is b = ( 4 6 7 8 6 4 3 ).

Lemma 17. Fix some τ , 0 ≤ τ ≤ t − 1, and let b = ( b1 ... bβ ) be the vectorassociated with the (aqτ

, apτ+1, aqτ+1)–peak (Uτ , Vτ , Wτ ), defined according to(4.6)–(4.8), where α = |Uτ ∪Vτ | and β = α+ |Wτ |. Also let U ′ be a partitionof aqτ+1, and let c1, . . . , cρ be positive integers with

ρ∑

j=1

cj = apτ+1 − aqτ+1. (4.9)

Then the following statements are equivalent.

1. There exist integers lj, rj with 1 ≤ lj ≤ pτ+1 ≤ rj < qτ+1 (j = 1, . . . , ρ),

such that for a′ = a(τ) −ρ∑

j=1

cjs(j), where

s(j)i =

{


0 otherwise(j = 1, . . . , ρ; i = 1, . . . , n)

we have

(a) 0 ≤ a′1 ≤ a′

2 ≤ · · · ≤ a′qτ+1

(b) {a′i − a′

i−1 : 1 ≤ i ≤ qτ+1, a′i 6= a′

i−1} = U ′ (where a′0 = 0).


2. There exist integers l′j, r′j with 1 ≤ l′j ≤ r′j ≤ β − 1 for j = 1, . . . , ρ,

such that for b′ = b −∑ρ

j=1 cjf(j), where

f(j)i =

{

1 if l′j ≤ i ≤ r′j0 otherwise

(j = 1, . . . , ρ; i = 1, . . . , β)

we have

(a) b′1 ≤ b′2 ≤ · · · ≤ b′β = bβ

(b) {b′i − b′i−1 : 1 ≤ i ≤ β, b′i 6= b′i−1} = U ′ (where b′0 = 0).

Proof. Let

R1 = {i : 1 ≤ i ≤ pτ+1, a(τ)i 6= a

(τ)i−1},

R2 = {i : pτ+1 ≤ i < qτ+1, a(τ)i 6= a

(τ)i+1}.

Clearly,

Uτ ∪ Vτ = {a(τ)i − a

(τ)i−1 : i ∈ R1} and

Wτ = {a(τ)i − a

(τ)i+1 : i ∈ R2}.

But by construction of b we also have

Uτ ∪ Vτ = {bi − bi−1 : 1 ≤ i ≤ α} and

Wτ = {bi − bi+1 : α ≤ i ≤ β − 1}.

Together this implies that there are bijections

ϕ1 : R1 → {1, . . . , α}, ϕ2 : R2 → {α, . . . , β − 1},

such that

a(τ)i − a

(τ)i−1 = bϕ1(i) − bϕ1(i)−1 for i ∈ R1 and

a(τ)i − a

(τ)i+1 = bϕ2(i) − bϕ2(i)+1 for i ∈ R2.

It is an easy consequence of the results of [10], that from the assumption(4.9) it follows that for lj, rj (j = 1, . . . , ρ) as in the first statement, wehave lj ∈ R1 and rj ∈ R2 for all j and for l′j, r′j (j = 1, . . . , ρ) as in thesecond statement we have l′j ≤ α and r′j ≥ α for all j. Suppose that lj, rj

(j = 1, . . . , ρ) satisfy the conditions of the first statement. The difference ofthe entries number i and i − 1 changes only when lj = i or rj = i − 1 forsome j. Thus, if i 6∈ R1 and i − 1 6∈ R2 we have

a′i − a′

i−1 = a(τ)i − a

(τ)i−1 = 0.


Hence, for i = 1, . . . , qτ+1,

a′i − a′

i−1 6= 0 =⇒ i ∈ R1 or i − 1 ∈ R2.

Put

C1(i) = {j ∈ [ρ] : lj = i} for i ∈ R1,

C2(i) = {j ∈ [ρ] : rj = i} for i ∈ R2.

Then

a′i − a′

i−1 = a(τ)i − a

(τ)i−1 −

∑

j∈C1(i)

cj for i ∈ R1

a′i − a′

i+1 = a(τ)i − a

(τ)i+1 −

∑

j∈C2(i)

cj for i ∈ R2.

By condition (a) of the first statement we have a′i − a′

i+1 ≤ 0 for i =0, . . . , qτ+1 − 1. For i ∈ R2 this yields

∑

j∈C2(i)

cj ≥ a(τ)i − a

(τ)i+1,

and together with

∑

i∈R2

∑

j∈C2(i)

cj =

ρ∑

j=1

cj = apτ+1 − aqτ+1 =∑

i∈R2

(

a(τ)i − a

(τ)i+1

)

we obtain for i ∈ R2,∑

j∈C2(i)

cj = a(τ)i − a

(τ)i+1.

and thus a′i − a′

i+1 = 0 for i ∈ R2. So the only nonzero differences a′i − a′

i−1

come from indices i ∈ R1. Now put l′j = ϕ1(lj) and r′j = ϕ2(rj) (j = 1, . . . , ρ)and let b′ be defined as in the second statement. Then l′j = ϕ1(i) iff j ∈ C1(i)and r′j = ϕ2(i) iff j ∈ C2(i), hence for i ∈ R1 we have

b′ϕ1(i) − b′ϕ1(i)−1 = bϕ1(i) − bϕ1(i)−1 −∑

j : l′j=ϕ1(i)

cj

= bϕ1(i) − bϕ1(i)−1 −∑

j∈C1(i)

cj

= ai − ai−1 −∑

j∈C1(i)

cj

= a′i − a′

i−1,


and for i ∈ R2,

b′ϕ2(i) − b′ϕ2(i)+1 = bϕ2(i) − bϕ2(i)+1 −∑

j : r′j=ϕ2(i)

cj

= bϕ2(i) − bϕ2(i)+1 −∑

j∈C2(i)

cj

= ai − ai+1 −∑

j∈C2(i)

cj

= a′i − a′

i+1 = 0.

So the second statement holds, and since all the arguments are reversible,we have proved that lj, rj (j = 1, . . . , ρ) satisfy the conditions of the firststatement iff l′j = ϕ1(lj), r′j = ϕ2(rj) (j = 1, . . . , ρ) satisfy the conditions ofthe second statement, and this proves the lemma. �

In fact the proof shows even more than just the equivalence of the twostatements: knowing l′j and r′j (j = 1, . . . , ρ) and R1 and R2, we can deter-

mine the lj, rj (j = 1, . . . , ρ) and R′ = {i : 1 ≤ i ≤ qτ+1, a(τ+1)i 6= a

(τ+1)i−1 }

in a number of steps that is bounded by a constant.

Example 7. Suppose a(τ) = ( 2 2 3 7 7 9 8 5 5 12 ) with Uτ = {2, 1}, Vτ ={4, 2}, Wτ = {3, 1}, R1 = {1, 3, 4, 6} and R2 = {6, 7}. The associated vectoris b = ( 4 6 8 9 6 5 ) and bijections as in the proof of Lemma 17 are given by

ϕ1 : 1 7→ 2, 3 7→ 4, 4 7→ 1, 6 7→ 3,

ϕ2 : 6 7→ 5, 7 7→ 4.

Now from

( 4 4 4 5 5 5 ) = ( 4 6 8 9 6 5 ) − ( 0 2 2 2 0 0 ) − ( 0 0 1 1 1 0 ) − ( 0 0 1 1 0 0 ),

where we have

l′1 = 2, r′1 = 4, l′2 = 3, r′2 = 5, l′3 = 3, r′3 = 4,

we obtain

l1 = 1, r1 = 7, l2 = 6, r2 = 6, l3 = 6, r3 = 7,

corresponding to

( 0 0 1 5 5 5 5 5 5 12 ) = ( 2 2 3 7 7 9 8 5 5 12 ) − ( 2 2 2 2 2 2 2 0 0 0 )

− ( 0 0 0 0 0 1 0 0 0 0 ) − ( 0 0 0 0 0 1 1 0 0 0 ).


Lemma 17 motivates the following definitions.

Definition 3. Let b = ( b1 ... bβ ) be the vector associated with some (u, v, w)–peak (U, V, W ) where α = |U ∪V | and β = α+ |W |, and let U ′ be a partitionof w. Let T be the set of positive integers ρ such that there are integersl1, . . . , lρ, r1, . . . , rρ and coefficients c1, . . . , cρ ∈ IN \ {0} such that

1.ρ∑

j=1

cj = v − w,

2. 1 ≤ lj ≤ rj ≤ β − 1 for j = 1, 2, . . . , ρ.

and for b′ = b −ρ∑

j=1

cjf(j), where

f(j)i =

{

1 if lj ≤ i ≤ rj,

0 otherwise,(j = 1, . . . , ρ; i = 1, . . . , β)

we have

3. b′1 ≤ b′2 ≤ · · · ≤ b′β = bβ = w and

4. {b′i − b′i−1 : 1 ≤ i ≤ β, b′i 6= b′i−1} = U ′ (with b′0 = 0).

Then we define

ρ(b, U ′) =

{

min T if T 6= ∅,∞ if T = ∅.

Definition 4. Let (U, V, W ) and (U ′, V ′, W ′) be a (u, v, w)–peak and a(u′, v′, w′)−peak, respectively, where u′ = w. Then we put

δ((U, V, W ), (U ′, V ′, W ′)) = ρ(b, U ′),

where b is the vector associated with (U, V, W ).

In order to model the segmentation process we define a digraph G =(V, E). The vertex set is

V = {(τ, U, Vτ , Wτ ) : 0 ≤ τ ≤ t, U is a partition of aqτ},

where

Vτ = {ai − ai−1 : qτ < i ≤ pτ+1, ai 6= ai−1},

Wτ = {ai − ai+1 : pτ+1 ≤ i < qτ+1, ai 6= ai+1}


·

2

d

g""EE

EEEE

EEEE

EEEEE

EEE

·b

∞

55lllllllllllllll

1

))RRRRRRRRRRRRRRR

∞

��:::

::::

::::

::::

::::

:::

·

1

55llllllllllllllla 1

))RRRRRRRRRRRRRRR · 3e //·

·

∞

BB��

155lllllllllllllll

1c

))RRRRRRRRRRRRRRR

·

4

f

<<yyyyyyyyyyyyyyyyyy

Fig. 4.1: The digraph for the vector a.

for 0 ≤ τ ≤ t. Observe that there is only one vertex with first component 0,namely (0, ∅, V0, W0) corresponding to a(0) = a and there is only one vertexwith first component t, namely (t, ∅, ∅, ∅) corresponding to the zero vector.In general, the vertices with first component τ represent the possibilities for(Uτ , Vτ , Wτ ), and by the observation before Definition 2 for each τ there isonly one choice for Vτ and Wτ , depending only on a. In the arc set E weinclude all arcs of the form

((τ, U, Vτ , Wτ ), (τ + 1, U ′, Vτ+1, Wτ+1))

for τ = 0, . . . , t − 1. Figure 4.1 shows G for a = ( 1 3 2 4 3 4 ), where thevertices are labeled as follows.

a = (0, ∅, {1, 2}, {1}), b = (1, {2}, {2}, {1}), c = (1, {1, 1}, {2}, {1}),

d = (2, {3}, {1}, {4}), e = (2, {2, 1}, {1}, {4}), f = (2, {1, 1, 1}, {1}, {4})

g = (3, ∅, ∅, ∅).

We define the arc weights in G to be the distances of the correspondingpeaks, i.e.

δ((τ, U, Vτ , Wτ ), (τ + 1, U ′, Vτ+1, Wτ+1)) = δ((U, Vτ , Wτ ), (U′, Vτ+1, Wτ+1))

for 0 ≤ τ ≤ t − 1 and all partitions U and U ′ of aqτand aqτ+1, respec-

tively. Observe that in this definition we used the fact that (U, Vτ , Wτ ) and(U ′, Vτ+1, Wτ+1) are an (aqτ

, apτ+1, aqτ+1)–peak and an (aqτ+1, apτ+2, aqτ+2)–peak, respectively. This assures that the condition u′ = w in the definitionof δ is satisfied. For instance, the segment ( 1 1 0 0 0 0 ) corresponds to the arc(a, b), since

a − ( 1 1 0 0 0 0 ) = ( 0 2 2 4 3 4 ),


while ( 0 1 0 0 0 0 ) corresponds to the arc (a, c), since

a − ( 0 1 0 0 0 0 ) = ( 1 2 2 4 3 4 ).

In general, an arc of weight ρ corresponds to a linear combination of ρ seg-ments. Now with a segmentation we can associate a path

(0, ∅, V0, W0), (1, U1, V1, W1), . . . , (t, ∅, ∅, ∅) (4.10)

in G.

Example 8. The segmentation

a = ( 1 1 0 0 0 0 ) + ( 0 1 1 1 0 0 ) + ( 0 1 1 1 1 1 ) + 2( 0 0 0 1 1 1 ) + ( 0 0 0 0 0 1 )

corresponds to the path (a, b, e, g) in Figure 4.1 as follows.

a = ( 1 3 2 4 3 4 )

(a, b) = −( 1 1 0 0 0 0 )

b = =( 0 2 2 4 3 4 )

(b, e) = −( 0 1 1 1 0 0 )

e = =( 0 1 1 3 3 4 )

(e, g) = −( 0 1 1 1 1 1 )

−( 0 0 0 2 2 2 )

−( 0 0 0 0 0 1 )

g = =( 0 0 0 0 0 0 ).

With these definitions the minimal number of segments needed to realizea segmentation corresponding to (4.10) equals the weight of this path.

Lemma 18. In time O(1) we can determine the values ρ(b, U ′) for all vectorsb that are associated with some (u, v, w)–peak and for all partitions U ′ of w.In addition we obtain values cj, l′j, r′j (j = 1, . . . , ρ(b, U ′)) satisfying theconditions of Definition 3.

Proof. The total number of vectors b associated with some (u, v, w)−peakswhen u, v and w run through all the possible values is

L∑

v=1

v−1∑

w=0

PvPv−w

where Pi is the number of partitions of i ∈ IN . Fix one of these vectors b.We consider all the sets S = {(l′j, r

′j, cj) : j = 1, . . . , ρ} (ρ ∈ IN), such


that the vectors f (1), . . . , f (ρ), defined as in Definition 3 and the coefficientsc1, . . . , cρ satisfy the conditions in Definition 3. We claim that there are atmost

vv−w ≤ LL

possibilities for S. Writing∑cj

k=1 f (j) for cjf(j) we can express

∑ρj=1 cjf

(j)

as a sum of∑ρ

j=1 cj = v − w (0, 1)−vectors. In order to satisfy Conditions1 and 3 of Definition 3, for i = α, . . . , β − 1, in exactly bi − bi+1 of these(0, 1)−vectors there must be 0 at position i + 1 and a 1 at position i. So wemay assume that the v−w right leaf positions are fixed. Since for each rightleaf position there are at most v left leaf positions the claim follows. For eachS the resulting partition U ′ of w can be computed in O(1) steps, since ρ isbounded by v −w ≤ L, and β is bounded by 2L. Thus the number of peaksis bounded by a constant, the number of sets S to be checked for each peakis bounded by a constant, for each of these sets the number of steps for thechecking is bounded by a constant, and this completes the proof. �

Lemma 19. In time O(n) we can determine the arc weights δ(e) for alle ∈ E and for each arc e a sequence

(s(1), c1), . . . , (s(δ(e)), cδ(e))

realizing its weight.

Proof. By Lemma 18 we may assume that we know all the ρ(b, U ′). First wedetermine in time O(n) the sets

P = {p1, . . . , pt},

Q = {q0, . . . , qt},

R1,τ = {i : qτ < i ≤ pτ+1, ai 6= ai−1} (τ = 0, . . . , t − 1),

R2,τ = {i : pτ+1 ≤ i < qτ+1, ai 6= ai+1} (τ = 0, . . . , t − 1),

and the partitions Vτ and Wτ (τ = 0, . . . , t). By induction, we assume thatwe have already determined the weights of the arcs up to layer τ for someτ , 0 ≤ τ ≤ t − 1. The number of vertices in layers τ and τ + 1 are boundedby Paqτ

and Paqτ+1, respectively. So the number of arcs is bounded by P2

L.Fix some (τ, Uτ , Vτ , Wτ ) and (τ +1, Uτ+1, Vτ+1, Wτ+1). Also by induction, weassume that we know the set

R1 = {i : 1 ≤ i ≤ pτ+1, a(τ)i 6= a

(τ)i−1}

for some possible a(τ) corresponding to (τ, Uτ , Vτ , Wτ ). Now by Lemma 17(and its proof) we obtain

δ((τ, Uτ , Vτ , Wτ ), (τ + 1, Uτ+1, Vτ+1, Wτ+1))


and a sequence realizing this value in constant time from the correspondingdata for b and U ′ where b is the vector associated with (Uτ , Vτ , Wτ ) andU ′ = Uτ+1. If τ ≤ t − 2 this also yields

R′1 = {i : 1 ≤ i ≤ pτ+2, a

(τ+1)i 6= a

(τ+1)i−1 }

for some possible a(τ+1) corresponding to (τ + 1, Uτ+1, Vτ+1, Wτ+1). So theweights for all arcs between adjacent layers can be determined in time O(1).And since the number of layers t+1 is bounded by n, the lemma is proved. �

Now the search for a segmentation with minimal NS amounts to the searchfor a path of minimal weight in a layered digraph with at most n layers wherethe number of vertices per layer is bounded by the constant PL. This can bedone in time O(n) ([14]). Thus we have proved

Theorem 4. L–One Row–Min MU–Min NS can be solved in time O(n).

4.2 Multiple row intensity maps

In this subsection we generalize the basic idea of the preceding subsection toprove that for bounded L and an MLC without ICC the NS–minimization ispolynomially solvable also for multiple row matrices. The problem L–Min

MU–Min NS is:

Instance: An integer matrix A = (ai,j) 1≤i≤m1≤j≤n

with 0 ≤ ai,j ≤ L (i ∈ [m], j ∈

[n]).

Problem: Find a segmentation of A with in first instance minimal TNMUand in second instance minimal NS!

Assume we have already determined the minimal TNMU c. From a seg-mentation of A we obtain a partition c = c1 + c2 + · · · + ck where ci is thecoefficient of the i–th segment (i = 1, . . . , k). First we consider the problemto check for a given partition if there is a segmentation of A with coefficientsc1, . . . , ck. This problem can be solved by checking the rows of A indepen-dently. For the moment we omit the row index and denote by a = ( a1 ... an )a fixed row of A and we put a0 = an+1 = 0. Compared to the single row casean additional difficulty in the multiple row case arises from the fact that theminimal TNMU that would be sufficient for a segmentation of a might besmaller than c. As a consequence we can not use Lemma 17, where condition(4.9) is essential. Here the order of the elements of the considered partition

4.2. Multiple row intensity maps 73

must be taken into consideration. For instance, for b = ( 2 5 0 ) there is asegmentation with coefficients 4, 1 and 1, namely

b = 4( 0 1 0 ) + ( 1 1 0 ) + ( 1 0 0 ),

while there is no segmentation with these coefficients for b′ = ( 3 5 0 ). Soinstead of peaks we have to consider ordered peaks to be defined below.Also, in order to describe the segmentation, we attach to a peak a multisetX of coefficients, and call the result an extended ordered peak. This is madeprecise in the following definition.

Definition 5. For integers v and w with 0 ≤ w < v ≤ L an extended ordered(v, w)–peak is a pair (b, X) of an integer vector b = ( b1 b2 ... bβ ), such thatthere is an integer α with 1 ≤ α < β and

0 < b1 < b2 < · · · < bα = v,

v = bα > bα+1 > · · · > bβ = w,

and a multiset X of positive integers. In addition, a pair (b, X), where b = ()is the empty tuple and X is a multiset of positive integers is called extendedordered (0, 0)–peak.

Example 9. (( 2 5 7 4 3 ) , {1, 2, 2, 3, 3}) is an extended ordered (7, 3)–peak(with α = 3, β = 5).

Let p1, . . . , pt and q0, . . . , qt be defined as in the preceding section. Thenfor a segmentation

a =k∑

j=1

cjs(j)

we can define k0(τ) and a(τ) (τ = 0, . . . , t) as before. Now for τ = 0, . . . , t,we associate with a(τ) an extended ordered (apτ+1, aqτ+1)–peak (b(τ), Xτ ) asfollows. For τ < t, let

Iτ = {i : 1 ≤ i ≤ pτ+1, a(τ)i 6= a

(τ)i−1},

Jτ = {i : pτ+1 < i ≤ qτ+1, a(τ)i 6= a

(τ)i−1},

denote the elements of Iτ by i1, . . . , iα and the elements of Jτ by iα+1, . . . , iβsuch that i1 < i2 < · · · < iβ, and put

b0 = 0, bl = ail (l = 1, . . . , β).

Let X0 = {c1, . . . , ck} and

Xτ+1 = Xτ \ {ck0(τ)+1, ck0(τ)+2, . . . , ck0(τ+1)} (τ = 0, . . . , t − 1).


Now for τ < t, (b(τ), Xτ) describes the initial part of a(τ) (up to columnqτ+1) together with the coefficients available for the remaining segments.In the final state (τ = t) we have the zero row a(t) = 0 and a multisetXt of coefficients, that are not needed for the considered row. With thezero row we associate the empty tuple b(t) = (), and thus we obtain fromany segmentation a sequence (b(0), X0), (b(1), X1), . . ., (b(t), Xt) of extendedordered peaks.

Example 10. Suppose a = ( 2 4 3 1 6 3 0 6 1 ) is a row in an intensity matrixwith minimal TNMU c = 18, and we are checking the partition c = 5 + 3 +2 + 2 + 2 + 1 + 1 + 1 + 1. Then from the segmentation

( 2 4 3 1 6 3 0 6 1 )

= ( 2 2 2 0 0 0 0 0 0 )

+ ( 0 1 0 0 0 0 0 0 0 )

+ ( 0 0 0 0 3 0 0 0 0 )

+ ( 0 0 0 0 2 2 0 0 0 )

+ ( 0 1 1 1 1 1 0 0 0 )

+ ( 0 0 0 0 0 0 0 5 0 )

+ ( 0 1 0 0 0 0 0 1 1 )

we obtain

τ a(τ) b(τ) Xτ

0 ( 2 4 3 1 6 3 0 6 1 ) ( 2 4 3 1 ) {5,3,2,2,2,1,1,1,1}1 ( 0 1 1 1 6 3 0 6 1 ) ( 1 6 3 0 ) {5,3,2,2,1,1,1}2 ( 0 0 0 0 0 0 0 6 1 ) ( 6 1 0 ) {5,2,1,1}3 ( 0 0 0 0 0 0 0 0 0 ) () {2,1}

That the vectors b(τ) provide enough information to construct the seg-mentation, follows from the simple observation, that w.l.o.g. a sequence ofconsecutive entries of equal value

ai1 = ai1+1 = · · · = ai2

can be considered as one single entry. This is intuitively clear and provedformally in the next lemma.

Lemma 20. Let a =∑k

j=1 cjs(j) be a segmentation with

s(j)i =

{


0 otherwise(j = 1, . . . , k).

There are integers l′j and r′j (j = 1, . . . , k) with the following properties.


1. We have a =k∑

j=1

cjs′(j) where

s′(j)i =

{


(j = 1, . . . , k).

2.

ai = ai−1 =⇒ s′(j)i = s

′(j)i−1 (i = 2, . . . , n; j = 1, . . . , k). (4.11)

Proof. In order to satisfy the last condition, we have to replace the segmentswith s

(j)i 6= s

(j)i−1 but ai = ai−1 for some i. Our strategy is to modify the given

segments as follows. For each sequence of consecutive entries of a of equalvalue we choose one representative, for instance the rightmost one, and adaptthe entries for each segment to the chosen column. This corresponds to thefollowing shifting of the leaves: if the left leaf covers a part of the plateau itis shifted to the right until the whole plateau is open, and if the right leafcovers a part of the plateau it is shifted to the left until the whole plateau iscovered.

First observe that s(j)i can differ from s

(j)i−1 only if i = lj or i − 1 = rj. So

for (4.11) it is sufficient that, for all j, we have

al′j6= al′j−1 and ar′j

6= ar′j+1. (4.12)

Suppose alj = alj−1 for some j. Then i1 < lj ≤ i2 for some i1, i2 with

ai1 = ai1+1 = · · · = ai2 = a and ai1−1, ai2+1 6= a. (4.13)

Since we want to adapt the entries of the segment to the rightmost columni2 we have to shift the left leaf to the left and put l′j = i1. Similarly, ifarj

= arj+1, then i1 ≤ rj < i2 for some i1, i2 with (4.13), and in order toadapt the entries of the segment to column i2, we have to shift the right leafto the left and put r′j = i1 − 1. In summary, for j ∈ [k] we put

l′j =

{

lj if alj 6= a′lj−1

max{i < lj : ai 6= alj} + 1 if alj = a′lj−1

r′j =

{

rj if arj6= a′

rj+1

max{i < rj : ai 6= arj} if arj

= a′rj+1

Then (4.12) is valid for all j, hence (4.11) is satisfied. In order to check the

first condition of the lemma, fix some i ∈ [n]. If s′(j)i = s

(j)i for all j, then

k∑

j=1

cjs′(j)i =

k∑

j=1

cjs(j)i = ai.


So assume s′(j)i 6= s

(j)i for some j. By construction this can be the case only

if ai = ai−1 or ai = ai+1. Now let i1 and i2 be the indices with i1 ≤ i ≤ i2,

ai1 = ai1+1 = · · · = ai = · · · = ai2 and ai1−1, ai2+1 6= ai.

We claim that s′(j)i = s

(j)i2

(j = 1, . . . , k). If s(j)i2

= 0, lj > i2 or rj < i2. Byconstruction, in the first case l′j > i2 and in the second case r′j < i1, so in

both cases s′(j)i = 0. If s

(j)i2

= 1, lj ≤ i2 and rj ≥ i2. By construction, l′j ≤ i1

and r′j ≥ i2, hence s′(j)i = 1 and the claim is proved. From this follows

k∑

j=1

cjs′(j)i =

k∑

j=1

cjs(j)i2

= ai2 = ai,

and since this argument works for any i ∈ [n] the first condition of the lemmais satisfied. �

By Lemma 20 applied to a(τ), w.l.o.g. we may assume that a(τ)lj

6= a(τ)lj−1

and a(τ)rj 6= a

(τ)rj+1 for all j > k0(τ). With this assumption the next lemma,

whose proof is obvious, justifies that we use the b(τ) instead of the a(τ).

Lemma 21. For fixed τ , 0 ≤ τ ≤ t−1, let b(τ) and Xτ be defined as describedabove and let {c1, . . . , cρ} ⊆ Xτ be fixed. If aqτ+1 6= 0 let g = ( g1 ... gγ ) besome vector with

0 < g1 < · · · < gγ = aqτ+1.

Then the following statements are equivalent.

1. There exist integers lj, rj with 1 ≤ lj ≤ rj < qτ+1, a(τ)lj

6= a(τ)lj−1 and

a(τ)rj 6= a

(τ)rj+1 (j = 1, . . . , ρ) such that for a′ = a(τ) −

ρ∑

j=1

cjs(j), where

s(j)i =

{


0 otherwise(j = 1, . . . , ρ; i = 1, . . . , n)

we have

(a) 0 ≤ a′1 ≤ a′

2 ≤ · · · ≤ a′qτ+1

= aqτ+1

(b) If aqτ+1 6= 0 there are exactly γ indices 1 ≤ i1 < · · · < iγ ≤ qτ+1

with a′i∗6= a′

i∗−1 (where a′0 = 0) and we have

(

ai1 ai2 . . . aiγ

)

= g.


2. There exist integers l′j, r′j with 1 ≤ l′j ≤ r′j ≤ β − 1 for j = 1, . . . , ρ,

such that for b′ = b −∑ρ

j=1 cjf(j), where

f(j)i =

{


(j = 1, . . . , ρ; i = 1, . . . , β)

we have

(a) b′1 ≤ b′2 ≤ · · · ≤ b′β = bβ = aqτ+1

(b) If aqτ+1 6= 0 there are exactly γ indices 1 ≤ i1 < · · · < iγ ≤ β withb′i∗ 6= b′i∗−1 (where b′0 = 0) and we have

(

bi1 bi2 . . . biγ

)

= g.

Now for τ = 0, 1, . . . , t − 1 the choice of the pairs(

sk0(τ)+1, ck0(τ)+1

)

, . . . ,(

sk0(τ+1), ck0(τ+1)

)

can be viewed as a way to go from the extended ordered (apτ+1, aqτ+1)–peak

(b(τ), Xτ ) to the extended ordered (apτ+2, aqτ+2)–peak (b(τ+1), Xτ+1) (withapt+1 = aqt+1 = 0).

Definition 6. Let 0 ≤ w < v and let (b, X) be an extended ordered (v, w)–peak, and let v′, w′ be integers with w ≤ v′ ≤ L and 0 ≤ w′ < v′ orv′ = w′ = 0. In addition let X ′ be a submultiset of X and denote theelements of X ′ by x1, . . . , x|X′|. We call an extended ordered (v′, w′)–peak(b′, X \X ′) accessible from (b, X) if there are integers l′1, . . . , l

′|X′|, r′1, . . . , r

′|X′|

such that

1. 1 ≤ l′j ≤ r′j ≤ β − 1 for j = 1, . . . , |X ′| (where b = ( b1 ... bβ )).

and for b′′ = b −|X′|∑

j=1

xjf(j), where

f(j)i =

{

1 if l′j ≤ i ≤ r′j,

0 otherwise,(j = 1, . . . , |X ′|; i = 1, . . . , β)

we have b′′ = 0 if v′ = w′ = 0 and otherwise

2. b′′1 ≤ b′′2 ≤ · · · ≤ b′′β = bβ = w and

3. If i1 < i2 < · · · < iγ′ are the indices with b′′i∗ 6= b′′i∗−1 (where b′′0 = 0),then

b′1 < b′2 < · · · < b′γ′ = w,

and we have(

b′′i1 b′′i2 . . . b′′iγ′

)

=(

b′1 b′2 . . . b′γ′

)

.


The definition can be interpreted as follows. Assume ap1 = v, aq1 = w,

ap2 = v′, aq2 = w′, let b(0) be associated with a(0) as above, and let b′ =( b′1 ... b′

β′ ) be a vector with

0 < b′1 < · · · < b′α′ = v′, v′ = b′α′ > · · · > b′β′ = w′.

Then (b′, X \ X ′) is accessible from (b(0), X) iff we can assign segments s(j)

to the elements of X ′, described by lj, rj (j = 1, . . . , |X ′|) with rj < q1 forall j, such that for

a(1) = a(0) −

|X′|∑

j=1

xjs(j)

we have a(1)1 ≤ a

(1)2 ≤ · · · ≤ a

(1)p2 and the extended ordered (v′, w′)–peak

associated with a(1) is (b′, X \ X ′).

Example 11. Let a = ( 0 2 5 5 7 4 3 3 5 6 8 2 ), X = {5, 3, 2, 2, 2, 1, 1, 1} andX ′ = {3, 1}. The associated extended ordered (7, 3)–peak is (b, X) whereb = ( 2 5 7 4 3 ). Now we want to determine the extended ordered (8, 0)–peaks(b′, X \ X ′) that are accessible from (b, X), where

b′ =(

b′1 . . . b′γ = 3 5 6 8 2)

.

We obtain that (b′, X \X ′) and (b′′, X \X ′) are accessible from (b, X), whereb′ = ( 2 3 5 6 8 2 ) and b′′ = ( 1 3 5 6 8 2 ):

( 2 2 3 3 3 ) = b − ( 0 3 3 0 0 ) − ( 0 0 1 1 0 ),

( 1 1 3 3 3 ) = b − ( 0 3 3 0 0 ) − ( 1 1 1 1 0 ).

This corresponds to the following possible beginnings of a segmentation.

( 0 2 5 5 7 4 3 3 5 6 8 2 )

− ( 0 0 3 3 3 0 0 0 0 0 0 0 )

− ( 0 0 0 0 1 1 0 0 0 0 0 0 )

= ( 0 2 2 2 3 3 3 3 5 6 8 2 )

and

( 0 2 5 5 7 4 3 3 5 6 8 2 )

− ( 0 0 3 3 3 0 0 0 0 0 0 0 )

− ( 0 1 1 1 1 1 0 0 0 0 0 0 )

= ( 0 1 1 1 3 3 3 3 5 6 8 2 ) .


On the other hand one can check that (( 3 5 6 8 2 ) , X \ X ′) is not accessiblefrom (b, X) and this corresponds to the fact that it is not possible to find(l1, r1) and (l2, r2) with r1, r2 < 7 such that after subtracting the correspond-ing segments with coefficients 3 and 1 from a we obtain a row vector a′ witha′

1 = · · · = a′i = 0, a′

i+1 = · · · = a′7 = 3 for some i, 1 ≤ i ≤ 6. Similar

statements can be made for b′ = ( 1 2 3 5 6 8 2 ).

Lemma 22. Let (b, X) be an extended ordered (v, w)−peak. Then the setof all (b′, X \ X ′) that are accessible from (b, X) can be determined in timeO(1).

Proof. Observe that the accessibility does not depend on the whole vec-tor b′ but only on the initial part ( b′1 ... b′

γ′=w ). So in order to determine

the accessible extended ordered peaks it is sufficient to determine the pairs(( b′1 ... b′

γ′ ), X \ X ′) of initial parts and multisets of coefficients. Let b =( b1 ... bβ ) and let α be the unique index with bα = v. We have b1 < · · · < bα

and bα > · · · > bβ. So for 1 ≤ k ≤ v − 1 there are at most two indices i

and i′ with 1 ≤ i, i′ ≤ β − 1 and bi = k, bi′ = k (namely the first one with1 ≤ i ≤ α − 1 and the second one with α + 1 ≤ i′ ≤ β − 1). The only indexi with bi = v is i = α, and so we have

β−1∑

i=1

bi ≤ v + 2

v−1∑

k=1

k ≤ L2.

Hence it is sufficient to consider at most PL2 candidates for X ′, where eachof these has at most L2 elements. Fix one of these X ′. Labeling the elementsof X ′ as in Definition 6, for each xj ∈ X ′ there are at most

(

2L−12

)

choices for

f (j). So the total number of choices for the pairs (f (j), xj) that have to beconsidered is bounded by

[(

2L − 1

2

)]|X′|

≤

[(

2L − 1

2

)]L2

.

For each of these choices the time needed to determine the resulting b′′ isbounded by a constant. Precisely, in order to subtract one of the xjf

(j) wehave to do at most 2L subtractions. So after at most L2 · 2L subtractionswe have determined b′′. Finally, in order to determine the corresponding( b′1 ... b′

γ′ ) according to condition 3 of Definition 6, we have to run throughthe at most 2L entries of b′. This proves the lemma, since the number ofsteps to determine the required data is bounded by


PL2

[(

2L − 1

2

)]L2

(L2 + 1)2L.

�

In order to model the segmentation we construct sets V0, . . . ,Vt of ex-tended ordered peaks. Put V0 = {(b(0), X0)} and suppose we have alreadyconstructed V0, . . . ,Vτ for some τ with 0 ≤ τ < t. Now we put

Vτ+1 = {(b′, X ′) : (b′, X ′) is an (apτ+2, aqτ+2) − peak that

is accessible from some (b, X) ∈ Vτ}.

Here for brevity of notation we put apt+1 = 1 and aqt+1 = 0. The elements of

Vτ represent the possibilities for (b(τ), Xτ). There is a segmentation of therow with coefficients c1, . . . , ck iff Vt 6= ∅. Note that a natural interpretationof this construction is a breadth first search (BFS) in the tree with vertex setV0 ∪ . . .∪Vt starting at (b(0), X0), where two vertices (b, X) and (b′, X ′) areconnected by an edge iff (b, X) ∈ Vτ , (b′, X ′) ∈ Vτ+1 for some τ and (b′, X ′)is accessible from (b, X).

Lemma 23. For given Vτ , the set Vτ+1 can be determined in time O(nL+1).

Proof. According to [10], the sum of the elements of X0, which is the minimalTNMU equals

c = max1≤i≤m

n∑

j=1

max{0, ai,j − ai,j−1} ≤ nL.

Now in any partition c = c1+· · ·+ck where the ci (i ∈ [k]) are the coefficientsof a segmentation of A, we have ci ≤ L for i ∈ [k]. Hence such a partition canbe described by an L–tuple (λ1, . . . , λL) of integers, where λr is the numberof summands equal to r for r ∈ [L]. Then

λr ≤nL

r(r ∈ [L]),

and so there are O(nL) choices for X0. Now the multiset X in

(b, X) ∈ Vτ

is a partition of some c′ with 0 ≤ c′ ≤ c ≤ nL with all summands less thanor equal to L. So there are nL possibilities for c′, and for each of these thereare O(nL) possible partitions. Thus the number of choices for X is bounded


by O(nL+1). The vectors b in the elements of Vτ differ only in the initialpart ( b1 ... bγ ), where bγ = aqτ

. But these initial parts are in bijection to theordered partitions of aqτ

, and of these there are (see for instance [1])

aqτ∑

i=1

(

aqτ− 1

i − 1

)

≤ L

(

L

bL2c

)

.

Since L is bounded by a constant we obtain that |Vτ | is bounded by O(nL+1).By Lemma 22, for each (b, X) ∈ Vτ the set of accessible (b′, X \ X ′) can bedetermined in time bounded by a constant, and this yields the claim. �

Lemma 24. For a fixed partition c = c1 + · · ·+ ck, it can be checked in timeO(nL+2) if there is a segmentation of a with coefficients c1, . . . , ck.

Proof. We only have to check if Vt 6= ∅. Since t ≤ n the claim is an immediateconsequence of Lemma 23. �

Now we can prove

Theorem 5. The problem L–Min MU–Min NS can be solved in timeO(mn2L+2).

Proof. Obviously,

c = max1≤i≤m

n∑

j=1

max{0, ai,j − ai,j−1}

can be determined in time O(mn). As in the proof of Lemma 23 the numberof partitions of c = c1 + · · · + ck that have to be considered is bounded byO(nL). By Lemma 24, for a fixed partition c = c1 + · · ·+ ck it can be checkedin time O(mnL+2) if there is a segmentation of A with coefficients c1, . . . , ck,and this concludes the proof. �

We finish this section with a remark concerning practical aspects of thisresult. Though the time complexity of the NS–minimization is polynomial inm and n the exponent grows linearly with L and also the L–dependent con-stants that were used to estimate the time–complexities of the different stepsof the algorithm, grow rapidly with L. So we expect an efficient algorithmonly for very small L. In the proof of the polynomiality we constructed thewhole sets Vτ (τ = 1, . . . , t), i.e. we performed a BFS as described beforeLemma 23. But in order to decide if there is a segmentation with the con-sidered coefficients we need to know only if Vt is nonempty, and in order toreconstruct a segmentation basically one path from the unique element of V0

to some element of Vt is sufficient. So for practical purposes it is natural touse depth first search (DFS) instead of BFS.


4.3 Test results

We implemented the algorithm described in the preceding section and Tables4.1 and 4.2 show test results for random 10 × 10– and 15 × 15–matrices,respectively. The computations where done on a 2 GHz workstation andwe determined the minimal NS for 1000 randomly generated matrices withmaximal entry L. The entry in column ’max. time’ is the maximal timeneeded for one single matrix, and the entry in column ’total time’ is the timeneeded for all the 1000 matrices. For comparison the tables also containheuristic results that were obtained with a slightly improved version of thealgorithm described in [10]. In order to evaluate the performance of the

exact heuristicL NS max. time total time NS total time3 6.9 1 s 9 s 6.9 0.9 s4 7.6 1 s 13 s 7.8 1.0 s5 8.1 1 s 29 s 8.4 1.1 s6 8.5 21 s 99 s 8.9 1.2 s7 8.8 50 s 5.6 min 9.3 1.2 s8 9.1 66 s 6.2 min 9.7 1.3 s9 9.3 3.4 min 16.1 min 10.0 1.3 s10 9.5 5.6 min 41.3 min 10.3 1.4 s11 9.8 11.0 min 1.3 h 10.6 1.4 s12 9.9 24.0 min 2.0 h 10.9 1.5 s13 10.0 1.4 h 7.0 h 11.1 1.5 s

Tab. 4.1: Average number of segments for random 10×10–matrices with maximalentry L. Each entry is averaged over 1000 matrices.

heuristic we determined the differences between the heuristic values and theexact minimums. Tables 4.3 and 4.4 show the frequencies of the values ofthe differences when 1000 matrices where treated for each value of L. Weconclude that for the considered range of parameters the exact algorithmyields only small improvements in terms of the number of segments, whilethe computational effort is extremely high already for small values of L. Sofor practical purposes the heuristic seems to be a good compromise betweencomputation time and accuracy of the optimization. Finally, we also testedour algorithm with the 13 clinical matrices, and the results are shown inTable 4.5.


exact heuristicL NS max. time total time NS total time3 9.7 1 s 16 s 9.8 4.8 s4 10.7 1 s 31 s 10.9 5.4 s5 11.3 12 s 175 s 11.7 5.8 s6 11.8 54 s 18.6 min 12.4 6.5 s7 12.3 3.1 min 0.8 h 13.0 6.8 s8 12.6 4.5 h 14.7 h 13.5 7.1 s9 12.9 24.1 h 37.9 h 14.0 7.4 s10 13.2 10.0 h 44.7 h 14.5 7.6 s

Tab. 4.2: Average number of segments for random 15×15–matrices with maximalentry L. Each entry is averaged over 1000 matrices.

L 0 1 2 33 969 31 0 04 876 123 1 05 780 218 2 06 663 331 2 07 525 456 19 08 437 516 47 09 335 603 62 010 306 584 104 611 262 615 121 212 168 654 173 513 141 641 213 5

Tab. 4.3: Frequencies of the differencesbetween the heuristic num-ber of segments and the exactminimum for 10×10–matrices.

L 0 1 2 3 43 940 60 0 0 04 809 189 2 0 05 609 379 12 0 06 453 509 37 1 07 327 585 86 2 08 250 594 151 5 09 150 609 230 11 010 85 551 335 28 1

Tab. 4.4: Frequencies of the differencesbetween the heuristic num-ber of segments and the exactminimum for 15×15–matrices.


exact heuristicno. MU NS CPU–time NS CPU–time1 16 7 0.04 s 8 0.01 s2 16 7 0.19 s 7 0.00 s3 20 8 0.39 s 8 0.01 s4 19 7 0.04 s 8 0.00 s5 15 7 0.01 s 7 0.00 s6 17 8 0.70 s 9 0.00 s7 18 7 0.03 s 7 0.00 s8 22 9 1.30 s 9 0.01 s9 26 9 25.77 s 10 0.00 s10 22 8 0.62 s 9 0.00 s11 22 10 7.88 s 10 0.00 s12 23 9 1.96 s 10 0.01 s13 23 9 2.36 s 9 0.01 s

Tab. 4.5: Test results for clinical matrices

5. FURTHER RESULTS

In this chapter we discuss some more results concerning the use of multi-leaf collimators without interleaf collision constraint. In the first section wepropose to change the direction of leaf movement between the different seg-ments and present some numerical results for a heuristic approach to theTNMU–minimization in this setup. A possible further step in this directionis discussed in the second section: one could use two MLC’s, perpendicularto each other, so that there is a leaf pair for each row and a leaf pair for eachcolumn. The numerical tests presented here suggest that this might yieldconsiderable savings in terms of the TNMU.

5.1 Using the MLC in two directions

Here and in the next section we are concerned only with the minimization ofthe TNMU, neglecting the number of segments. So for simplicity we assumethat every segment has coefficient 1. In this section we search for a realizationof the given intensity matrix A in the following setup: the multileaf collimatorhas no interleaf collision constraint and can be rotated about 90◦. So thereare two different types of segments, called horizontal and vertical segments,according to the choice of the direction of leaf motion. In the first subsectionwe derive a backtracking algorithm for the TNMU–minimization problem.Due to the computational complexity this algorithm is applicable only forsmall problem sizes (10 × 10–matrices with entries between 0 and 7), butin any case we obtain a lower bound for the TNMU by interrupting thebacktracking after some time. In the second subsection these lower boundsare compared with heuristic results.

5.1.1 A lower bound

Theorem 6. Let c1, c2 be nonnegative integers and put c = c1 + c2. Thena segmentation of A with c1 horizontal and c2 vertical segments exists iffB := cJ −A (where J is the all–one–matrix of size m× n) can be written asa sum of four nonnegative integer matrices P = (pi,j), Q = (qi,j), R = (ri,j),S = (si,j) with the following properties.

86 5. Further Results

1. pi,j + qi,j ≤ c1, ri,j + si,j ≤ c2 for (i, j) ∈ [m] × [n]

2. pi,j ≥ pi,j+1, qi,j ≤ qi,j+1 for i ∈ [m], j ∈ [n − 1]

3. ri,j ≥ ri+1,j, si,j ≤ si+1,j for i ∈ [m − 1], j ∈ [n]

Proof. ”⇒”: By construction, bi,j is the number of times the bixel (i, j) hasto be covered in a segmentation with c monitor units. Suppose there isgiven a segmentation with c1 horizontal and c2 vertical segments. For(i, j) ∈ [m] × [n], let pi,j, qi,j, ri,j and si,j be the number of segmentsin which bixel (i, j) is covered by the left, the right, the upper and thelower leaf, respectively. This yields the desired decomposition of B.

”⇐”: Suppose B = P +Q+R+S where P , Q, R and S satisfy the conditionsof the theorem. Now we define segments S(1), S(2), . . . , S(c) as follows.For 1 ≤ k ≤ c1, let

σ(k)i,j =

{

0 if k ≤ pi,j,

1 otherwise,τ

(k)i,j =

{

0 if k > c1 − qi,j,

1 otherwise.

For c1 + 1 ≤ k ≤ c, let

σ(k)i,j =

{

0 if k − c1 ≤ ri,j,

1 otherwise,τ

(k)i,j =

{

0 if k − c1 > c2 − si,j,

1 otherwise.

Finally, for all k ∈ [c], put s(k)i,j = σ

(k)i,j τ

(k)i,j . To conclude the proof we

have to check that the matrices S(k) = (s(k)i,j ) are segments and that

their sum is A. Condition 1 in the theorem implies that, for each k andeach (i, j), at most one of the numbers σ

(k)i,j and τ

(k)i,j is 0, corresponding

to the fact that a bixel can not be covered by both leaves of a leaf pairat the same time. Thus the number of indices k with s

(k)i,j = 0 equals

pi,j + qi,j + ri,j + si,j = bi,j, hence

c∑

k=1

s(k)i,j = c − bi,j = ai,j.

We show that S(k) is a horizontal segment for k ≤ c1. For k > c1 oneobtains similarly that S(k) is a vertical segment. If s

(k)i,j = 0, either

σ(k)i,j = 0 or τ

(k)i,j = 0. By construction of the σ

(k)i,j , τ

(k)i,j and by Condition

2 in the theorem this implies σ(k)i,j′ = 0 for all j ′ < j or τ

(k)i,j′ = 0 for all

j ′ > j, and consequently s(k)i,j = 0 for all j ′ < j or for all j ′ > j. �

5.1. Using the MLC in two directions 87

Example 12. For A =

(

1 4 2 51 3 3 21 3 2 56 4 6 0

)

and c = 6 = 4 + 2 we have

B =

(

5 2 4 15 3 3 45 3 4 10 2 0 6

)

=

(

3 2 2 02 1 0 03 2 0 00 0 0 0

)

+

(

0 0 0 01 1 1 20 0 0 00 0 0 4

)

+

(

2 0 2 12 0 2 12 0 2 00 0 0 0

)

+

(

0 0 0 00 1 0 10 1 0 10 2 0 2

)

.

This yields

σ(1) =

(

0 0 0 10 0 1 10 0 1 11 1 1 1

)

, τ (1) =

(

1 1 1 11 1 1 11 1 1 11 1 1 0

)

, S(1) =

(

0 0 0 10 0 1 10 0 1 11 1 1 0

)

,

σ(2) =

(

0 0 0 10 1 1 10 0 1 11 1 1 1

)

, τ (2) =

(

1 1 1 11 1 1 11 1 1 11 1 1 0

)

, S(2) =

(

0 0 0 10 1 1 10 0 1 11 1 1 0

)

,

σ(3) =

(

0 1 1 11 1 1 10 1 1 11 1 1 1

)

, τ (3) =

(

1 1 1 11 1 1 01 1 1 11 1 1 0

)

, S(3) =

(

0 1 1 11 1 1 00 1 1 11 1 1 0

)

,

σ(4) =

(

1 1 1 11 1 1 11 1 1 11 1 1 1

)

, τ (4) =

(

1 1 1 10 0 0 01 1 1 11 1 1 0

)

, S(4) =

(

1 1 1 10 0 0 01 1 1 11 1 1 0

)

,

σ(5) =

(

0 1 0 00 1 0 00 1 0 11 1 1 1

)

, τ (5) =

(

1 1 1 11 1 1 11 1 1 11 0 1 0

)

, S(5) =

(

0 1 0 00 1 0 00 1 0 11 0 1 0

)

,

σ(6) =

(

0 1 0 10 1 0 10 1 0 11 1 1 1

)

, τ (6) =

(

1 1 1 11 0 1 01 0 1 01 0 1 0

)

, S(6) =

(

0 1 0 10 0 0 00 0 0 01 0 1 0

)

.

Based on Theorem 6, Algorithm 4 finds the minimal TNMU c and adecomposition of B = cJ − A corresponding to a segmentation with thisTNMU. The step in line 7 can be realized using backtracking. Of course thismethod is very time–consuming. Our implementation solved the problemfor 10× 10–matrices with random entries between 0 and 7 in a few seconds,but for larger problems this algorithm is not practicable. By interruptingthe backtracking when some fixed time limit is reached without getting aconstruction or a contradiction in line 7, we still obtain a lower bound forthe TNMU.

5.1.2 Heuristic results

For notational convenience we add a 0−th and an (m + 1)−th row witha0,j = am+1,j = 0 (j ∈ [n]), and put

α(i) =

n∑

j=1

max{0, ai,j − ai,j−1} (i ∈ [m]),

β(j) =m∑

i=1

max{0, ai,j − ai−1,j} (j ∈ [n]).


Algorithm 4 Minimal TNMU for 2–directional segmentation

c = max{ai,j : (i, j) ∈ [m] × [n]}finished:=falsewhile not finished do

B := cJ − A

5: c1 := c; c2 := 0while not finished and c1 ≥ 0 do

Construct a decomposition of B as in Theorem 6 or derive a contra-diction from the assumption that such a decomposition existsif Construction successful then

finished:=true10: else

c1 := c1 − 1; c2 := c2 + 1end if

end whileif not finished then

15: c := c + 1end if

end while

Let also

chor(A) = max1≤i≤m

α(i), cvert(A) = max1≤j≤n

β(j).

According to [10] there is a segmentation without vertical segments withchor(A) monitor units, and there is a segmentation without horizontal seg-ments with cvert(A) monitor units. Hence an upper bound for the minimalnumber of monitor units needed to realize A is c(A) = min{chor(A), cvert(A)}.

As in Chapter 2 we construct a segmentation by successively subtractingsegments until the zero matrix is reached. For the choice of the segment S

we suggest a heuristic method. First chor(A) and cvert(A) are determined andthe direction of leaf motion is chosen to be horizontal if chor(A) ≤ cvert(A)and vertical otherwise. We describe the construction of S for the horizontalcase, the vertical case is treated analogously. For the new matrix A′ := A−S,we have two aims. Firstly, we want that

chor(A′) = chor(A) − 1, (5.1)

and secondly cvert(A′) should be small. Let us consider the second condition

first. With the segments we associate (D, D′)–paths in a layered digraph


Γ = (V, E) similar to the digraph Γ in Chapter 3:

V = {D, D′} ∪ {(i, l, r) : i ∈ [m], 1 ≤ l ≤ r + 1 ≤ n + 1,

ai,j > 0 for l ≤ j ≤ r},

and E consists of all possible arcs between row i and row i+1 (i = 1, . . . , m−1), all arcs between D and row 1, and all arcs between row m and D′. Nowwe define a weight function on E and determine the segment S as a pathof maximal weight in Γ. This approach seems to be natural because thevalues a′

i,j − a′i−1,j, and thus cvert(A

′) depend on the combinations of the leafpositions in adjacent rows. Experiments have shown that for the consideredrange of parameters the weight function w described below works quite well.

w(D, (1, l, r)) =r∑

j=l

max{1, β(j) + 5 − cvert(A)} (1 ≤ l ≤ r + 1 ≤ n)

w((m, l, r), D′) = 0 (1 ≤ l ≤ r + 1 ≤ n + 1)

In order to define the weights for the arcs ((i − 1, l, r), (i, l′, r′)) we put, fori = 2, 3, . . . , m and j = 1, 2, . . . , n,

ρ1(i, j) =

{

max{1, β(j) + 5 − cvert(A)} if ai,j ≥ ai−1,j

0 otherwise,

ρ2(i, j) =

{

max{1, β(j) + 5 − cvert(A)} if ai,j > ai−1,j

0 otherwise.

ρk(i, j) (k = 1, 2) measure the influence of a potential segment with si,j 6=si−1,j on ai,j − ai−1,j and thus on β(j). ρ1(i, j) is nonzero if a segment withsi,j = 0 and si−1,j = 1 increases the value of max{0, ai,j − ai−1,j}, and sim-ilarly ρ2(i, j) is nonzero if a segment with si,j = 1 and si−1,j = 0 decreasesmax{0, ai,j − ai−1,j}. The magnitudes of the ρk(i, j) (k = 1, 2) are chosenaccording to the idea that columns j with β(j) close to cvert(A) should beconsidered more important than columns j with a small value of β(j). Now


we define w as follows.

w((i − 1, l, r), (i, l′, r′)) =

−l′−1∑

j=l

ρ1(i, j) −r∑

j=r′+1

ρ1(i, j)

if l ≤ l′ ≤ r′ + 1 ≤ r + 1,

−l′−1∑

j=l

ρ1(i, j) +r′∑

j=r+1

ρ2(i, j)

if l ≤ l′ ≤ r + 1 ≤ r′,

−r∑

j=l

ρ1(i, j) +r′∑

j=l′ρ2(i, j)

if l ≤ r + 1 < l′ ≤ r′ + 1,

l−1∑

j=l′ρ2(i, j) −

r∑

j=r′+1

ρ1(i, j)

if l′ < l ≤ r′ + 1 ≤ r + 1,

l−1∑

j=l′ρ2(i, j) +

r′∑

j=r+1

ρ2(i, j)

if l′ < l ≤ r + 1 ≤ r′,

r′∑

j=l′ρ2(i, j) −

r∑

j=l

ρ1(i, j)

if l′ ≤ r′ + 1 < l ≤ r + 1.

Example 13. Let rows i − 1 and i of matrix A be

( 2 4 1 4 3 4 3 31 5 3 2 7 4 3 2 ) ,

and consider the arc e = ((i−1, 4, 7), (i, 2, 6)), corresponding to the two rows

( 0 0 0 1 1 1 1 00 1 1 1 1 1 0 0 ) .

If we choose a segment with these rows the corresponding part of A′ is

( 2 4 1 3 2 3 2 31 4 2 1 6 3 3 2 )

The relevant columns for w(e) are columns number 2, 3 and 7. Assume

cvert(A) = 23, β(2) = 20, β(3) = 21 and β(7) = 17.

Then

ρ1(i, 2) = 2, ρ1(i, 3) = 3, ρ1(i, 7) = 1

ρ2(i, 2) = 2, ρ2(i, 3) = 3, ρ2(i, 7) = 0


and we obtain

w(e) = 2 + 3 − 1 = 4.

Here the positive terms correspond to

max{0, a′i,2 − a′

i−1,2} = 0 < 1 = max{0, ai,2 − ai−1,2} and

max{0, a′i,3 − a′

i−1,3} = 1 < 2 = max{0, ai,3 − ai−1,3},

while the negative term corresponds to

max{0, a′i,7 − a′

i−1,7} = 1 > 0 = max{0, ai,7 − ai−1,7}.

Finally, we delete all vertices (i, l, r) that lead to segments violating (5.1).Fix some vertex (i, l, r) and put

a′i,j =

{

ai,j − 1 if l ≤ j ≤ r,

ai,j otherwise.

A segment corresponding to a path through (i, l, r) can not satisfy (5.1) ifα′(i) :=

∑n

j=1 max{0, a′i,j −a′

i,j−1} > chor(A)− 1. But the only terms in α′(i)that could be different from the corresponding terms in α(i) are the termsfor j = l and j = r + 1 (if r < n), and for these terms we have, if l ≤ r,

max{0, a′i,l − a′

i,l−1} =

{

max{0, ai,l − ai,l−1} − 1 if ai,l > ai,l−1

max{0, ai,l − ai,l−1} if ai,l ≤ ai,l−1,

max{0, a′i,r+1 − a′

i,r} =

{

max{0, ai,r+1 − ai,r} + 1 if ai,r ≤ ai,r+1

max{0, ai,r+1 − ai,r} if ai,r > ai,r+1,

So

α′(i) =

α(i) − 1 if l ≤ r, ai,l > ai,l−1 and ai,r > ai,r+1,

α(i) if l = r + 1,α(i) if l ≤ r, ai,l > ai,l−1 and ai,r ≤ ai,r+1,

α(i) if l ≤ r, ai,l ≤ ai,l−1 and ai,r > ai,r+1,

α(i) + 1 if l ≤ r, ai,l ≤ ai,l−1 and ai,r ≤ ai,r+1.

Consequently, we have to delete all vertices satisfying one of the followingconditions.

1. α(i) = chor(A) and r = l − 1

2. α(i) = chor(A) and (ai,l ≤ ai,l−1 or ai,r ≤ ai,r+1)

3. α(i) = chor(A) − 1 and (l ≤ r, ai,l ≤ ai,l−1 and ai,r ≤ ai,r+1)


Now we choose a segment S corresponding to a (D, D′)−path of maximalweight. In Table 5.1 the results of the heuristic are compared to the averagelower bound obtained by the backtracking method from Subsection 5.1.1 andto the optimal TNMU when the MLC is used in only one direction. On a 2GHz workstation the computation times for the whole heuristic columns were172 seconds (10×10) and 926 seconds (15×15). For the clinical matrices our

10 × 10 15 × 15L TNMU TNMU Lower TNMU TNMU Lower

(old) (new) bound (old) (new) bound3 9.8 7.2 6.7 14.0 10.0 8.34 12.6 9.0 8.3 17.9 12.4 10.25 15.5 10.8 10.1 21.7 14.8 12.36 18.1 12.6 11.8 25.6 17.2 14.57 20.8 14.3 13.5 29.4 19.6 16.28 23.6 15.9 14.7 33.2 21.9 18.09 26.4 17.7 16.6 37.0 24.2 20.110 29.0 19.5 18.4 40.9 26.6 21.811 31.8 21.3 19.6 44.7 29.0 24.012 34.5 23.0 21.9 48.5 31.4 25.713 36.9 24.6 22.8 52.3 33.7 27.814 39.8 26.3 24.8 56.2 36.1 29.515 42.4 28.2 26.2 59.8 38.3 31.316 45.2 29.8 27.9 63.3 40.7 32.6

Tab. 5.1: Average TNMU for random 10 × 10– and 15 × 15–matrices withmaximal entry L when the MLC is used in both directions (column’TNMU(new)’), compared to the average lower bound and the resultsfor the one–directional usage of the MLC (column ’TNMU(old)’).

implementation of Algorithm 4 yields the optimal solution in a reasonabletime and the results are shown in Table 5.2.

5.2 Two orthogonal MLCs

In this section we suppose that two MLCs without ICC are arranged in sucha way that the leaf pairs of the one are perpendicular to the leaf pairs of theother, and the segments are the (0, 1)–matrices describing any combinationof leaf positions. For instance,

( 0 1 1 0 1 01 0 1 1 1 11 0 1 1 0 01 0 1 1 1 11 0 1 0 0 01 0 0 0 1 1

)

(5.2)

5.2. Two orthogonal MLCs 93

no. MU (new) MU (old) CPU–time1 11 16 0.48 s2 11 16 0.27 s3 12 20 0.54 s4 12 19 50.06 s5 13 15 0.18 s6 12 17 0.19 s7 13 18 0.91 s8 14 22 1.05 s9 16 26 4.04 s10 14 22 13.07 s11 15 22 93.68 s12 16 23 1.59 s13 14 23 201.07 s

Tab. 5.2: Test results for the clinical matrices: the TNMU when the MLC is usedin both directions (’MU (new)’) compared with the TNMU when theMLC is used in only one direction (’MU (old)’), and the time needed todetermine the segmentation according to Algorithm 4

.

is a segment corresponding to the leaf positions shown in Figure 5.1. The LP–formulation of the TNMU–segmentation problem given in Section 2.1 is stillvalid, we just have a larger set F of allowed segments, and correspondinglythere are more inequalities that have to be satisfied by a dual feasible solutiong. Here we suggest a heuristic segmentation method based on dual feasiblesolutions of a particular type. For pairs (i, j) with 2 ≤ i ≤ m − 1 and2 ≤ j ≤ n − 1 we put

γ1(i, j) = max

{

1

3(ai,s + ai,t + au,j + av,j − ai,j) : 1≤s<j<t≤n,

1≤u<i<v≤m

}

. (5.3)

Thenc(A) = max{γ1(i, j) : 2 ≤ i ≤ m − 1, 2 ≤ j ≤ n − 1}

is a lower bound for the TNMU, because for any (i, j) with 2 ≤ i ≤ m − 1,2 ≤ j ≤ n − 1 and any numbers s, t, u and v with 1 ≤ s < j < t ≤ n and1 ≤ u < i < v ≤ m we can define a dual feasible solution g by putting

g(i, s) = g(i, t) = g(u, j) = g(v, j) =1

3, g(i, j) = −

1

3

and g(p, q) = 0 for all other (p, q) ∈ [m] × [n]. This is illustrated in Figure5.2. If we choose s, t, u and v so that the maximum in (5.3) is attained the


Fig. 5.1: Leaf positions corresponding to the segment (5.2).

j

13

u

i13

−13

13

13

v

s t

Fig. 5.2: The structure of the dual feasible solution g, where the empty spaces arefilled with zeros.


objective value of g is γ1(i, j) and to see that g is feasible we just have toobserve that in order to cover the bixel (i, j) we have to cover at least oneof the bixels (i, s), (i, t), (u, j), (v, j). Of course, c(A) does not have to bean integer, and since we are using only integer coefficients even dc(A)e canbe used as a lower bound for the TNMU. As in Chapter 2, our algorithmis based on the general principle of extracting segments from A. Precisely,depending on A we determine a segment S, such that

A′ = A − S

is still nonnegative, and then we iterate this step with A′ until the zeromatrix is reached. The main idea underlying our heuristic approach to theconstruction of S is that we try to decrease the value of c. Observe that

ai,j = 0 =⇒ si,j = 0 (i ∈ [m], j ∈ [n]) (5.4)

is a necessary and sufficient condition for the nonnegativity of A′ = A − S.Now we start with S equal to the all–one matrix of dimension m×n and coversuccessively all the zero–entries of A. The construction of S is described inAlgorithm 5. Lines 2 to 9 ensure that condition (5.4) is satisfied for bixels

Algorithm 5 Segment for two orthogonal MLCs

si,j := 1 (i ∈ [m], j ∈ [n])for j = 1 to n do

if a1,j = 0 then s1,j = 0if am,j = 0 then sm,j = 0

5: end forfor i = 1 to m do

if ai,1 = 0 then si,1 = 0if ai,n = 0 then si,n = 0

end for10: while (5.4) is violated do

Choose an (i, j) with ai,j = 0 and si,j = 1Choose a covering direction from {left, right, up, down}Cover bixel (i, j) with the leaf from the chosen directionLet S be the segment corresponding to the new leaf positions

15: end while

(i, j) with i ∈ {1, m} or j ∈ {1, n}. These bixels can be covered withoutinfluencing other bixels, for instance (i, 1) can be covered by the left leafof row i, so lines 2 to 9 imply no loss of generality. We still have to makeprecise, how we choose the bixel (i, j) in line 11 and the direction from which


we cover it in line 12. For the choice of the bixel we follow the strategy tocover bixels (i, j) with a high value of γ1(i, j) first. Suppose we have alreadychosen the bixel to cover. Now we have to choose the direction from whichwe want to cover it. For instance by covering bixel (i, j) by the left leaf wealso cover all bixels (i, j ′) with 1 ≤ j ′ ≤ j. So it is natural to cover bixel(i, j) from a direction with the property that the maximal entry of a bixelthat is covered although the entry is nonzero, is as small as possible. Let

γ(left)2 (i, j) = max{ai,j′ : 1 ≤ j ′ ≤ j},

γ(right)2 (i, j) = max{ai,j′ : j ≤ j ′ ≤ n},

γ(up)2 (i, j) = max{ai′,j : 1 ≤ i′ ≤ i},

γ(down)2 (i, j) = max{ai′,j : i ≤ i′ ≤ m},

and choose the covering direction of bixel (i, j),

dir(i, j) ∈ {left, right, up, down}

to be the direction with the smallest value of γ(∗)2 (i, j). To decide between

directions with equal value of γ(∗)2 (i, j) we consider the number of bixels (i′, j ′)

violating the condition for the segment (i.e. with ai′,j′ = 0 and si′,j′ = 1)that are covered when we cover bixel (i, j) from the respective direction. Tomake this precise, let

γ(left)3 (i, j) = |{(i, j ′) : 1 ≤ j ′ < j, ai,j′ = 0 and si,j′ = 1}| ,

γ(right)3 (i, j) = |{(i, j ′) : j < j ′ ≤ n, ai,j′ = 0 and si,j′ = 1}| ,

γ(up)3 (i, j) = |{(i′, j) : 1 ≤ i′ < i, ai′,j = 0 and si′,j = 1}| ,

γ(down)3 (i, j) = |{(i, j ′) : i < i′ ≤ m, ai′,j = 0 and si′,j = 1}| ,

and choose for dir(i, j) the direction with in first instance minimal value of

γ(∗)2 (i, j), and in second instance maximal value of γ

(∗)3 (i, j). If there is still

a tie it can be decided randomly. Finally, we put

γ3(i, j) = γ(dir(i,j))3 (i, j).

We choose among the bixels (i, j) with ai,j = 0 and si,j = 1 one with in firstinstance maximal value of γ1(i, j), and in second instance maximal value of


γ3(i, j). Now the chosen bixel (i, j) is covered by the leaf that is given bydir(i, j). Precisely, we put

si,j′ := 0 for

{

1 ≤ j ′ ≤ j if dir(i, j) = left,j ≤ j ′ ≤ n if dir(i, j) = right,

si′,j := 0 for

{

1 ≤ i′ ≤ i if dir(i, j) = up,

i ≤ i′ ≤ m if dir(i, j) = down.

This covering step is repeated until (5.4) is satisfied.

Example 14. Consider the matrix

A =

( 1 4 3 1 5 14 0 3 1 2 12 1 5 4 0 04 4 5 2 3 32 2 2 1 1 15 0 0 0 5 3

)

.

γ1(2, 2) = 5 is the maximal value of γ1(i, j), and (2, 2) can be covered fromthe right with maximal value of a covered bixel equal to 3, while from allother directions we would have to cover a bixel with entry 4. For all theother zero–bixels it is obvious how to cover them without covering a nonzerobixel. So the first segment is

( 1 1 1 1 1 11 0 0 0 0 01 1 1 1 0 01 1 1 1 1 11 1 1 1 1 11 0 0 0 1 1

)

,

and continuing we obtain a segmentation with 5 segments.

( 1 4 3 1 5 14 0 3 1 2 12 1 5 4 0 04 4 5 2 3 32 2 2 1 1 15 0 0 0 5 3

)

−

( 1 1 1 1 1 11 0 0 0 0 01 1 1 1 0 01 1 1 1 1 11 1 1 1 1 11 0 0 0 1 1

)

=

( 0 3 2 0 4 03 0 3 1 2 11 0 4 3 0 03 3 4 1 2 21 1 1 0 0 04 0 0 0 4 2

)

,

( 0 3 2 0 4 03 0 3 1 2 11 0 4 3 0 03 3 4 1 2 21 1 1 0 0 04 0 0 0 4 2

)

−

( 0 1 1 0 1 01 0 1 1 1 11 0 1 1 0 01 0 1 1 1 11 0 1 0 0 01 0 0 0 1 1

)

=

( 0 2 1 0 3 02 0 2 0 1 00 0 3 2 0 02 3 3 0 1 10 1 0 0 0 03 0 0 0 3 1

)

,

( 0 2 1 0 3 02 0 2 0 1 00 0 3 2 0 02 3 3 0 1 10 1 0 0 0 03 0 0 0 3 1

)

−

( 0 1 1 0 1 01 0 0 0 0 00 0 1 1 0 01 1 1 0 1 10 1 0 0 0 01 0 0 0 1 1

)

=

( 0 1 0 0 2 01 0 2 0 1 00 0 2 1 0 01 2 2 0 0 00 0 0 0 0 02 0 0 0 2 0

)

,

( 0 1 0 0 2 01 0 2 0 1 00 0 2 1 0 01 2 2 0 0 00 0 0 0 0 02 0 0 0 2 0

)

−

( 0 1 0 0 1 00 0 1 0 1 00 0 1 1 0 01 1 1 0 0 00 0 0 0 0 01 0 0 0 1 0

)

=

( 0 0 0 0 1 01 0 1 0 0 00 0 1 0 0 00 1 1 0 0 00 0 0 0 0 01 0 0 0 1 0

)

.


Table 5.3 shows the average TNMU obtained by our algorithm, the aver-age lower bound c(A) and the optimal TNMU when only one MLC is usedin one direction. On a 2 GHz workstation the computation for the whole10 × 10–column took 91 seconds and for the 15 × 15–column it took 610seconds. The results for the clinical sample matrices are shown in Table 5.4.

10 × 10 15 × 15L TNMU TNMU Lower TNMU TNMU Lower

(old) (new) bound (old) (new) bound3 9.8 5.4 3.9 14.0 7.4 4.04 12.6 6.7 5.2 17.9 9.0 5.35 15.5 8.0 6.4 21.7 10.5 6.76 18.1 9.2 7.6 25.6 12.0 7.97 20.8 10.4 8.8 29.4 13.6 9.38 23.6 11.7 10.0 33.2 15.1 10.59 26.4 12.9 11.2 37.0 16.6 11.810 29.0 14.1 12.4 40.9 18.1 13.011 31.8 15.4 13.6 44.7 19.6 14.212 34.5 16.6 14.8 48.5 21.1 15.713 36.9 17.9 16.0 52.3 22.5 16.914 39.8 19.1 17.1 56.2 24.1 18.115 42.4 20.3 18.3 59.8 25.5 19.616 45.2 21.5 19.5 63.3 27.0 20.4

Tab. 5.3: Average TNMU for random 10×10– and 15×15–matrices mith maximalentry L when two orthogonal MLCs are used (column ’TNMU(new)’),compared to the average lower bound and the results when one MLC isused in only one direction (column ’TNMU(old)’).


no. MU (1) MU (2) MU (3) Bound CPU–time1 16 11 11 11 0.00 s2 16 11 11 11 0.00 s3 20 12 11 11 0.00 s4 19 12 11 10 0.00 s5 15 13 10 10 0.00 s6 17 12 12 11 0.01 s7 18 13 11 11 0.00 s8 22 14 13 12 0.01 s9 26 16 13 12 0.01 s10 22 14 14 13 0.01 s11 22 15 14 13 0.01 s12 23 16 14 12 0.02 s13 23 14 14 12 0.02 s

Tab. 5.4: Test results for clinical matrices: the TNMU with one MLC in one di-rection (’MU (1)’), with one MLC in two directions (’MU (2)’) and withtwo orthogonal MLCs (’MU (3)’), the lower bounds for the last case(’Bound’) and the times needed to determine the segmentations for twoMLCs with our algorithm.

BIBLIOGRAPHY

[1] M. Aigner. Diskrete Mathematik. Vieweg, Braunschweig/Wiesbaden,4th edition, 2001.

[2] D. Baatar and H.W. Hamacher. New LP model for multileaf collimatorsin radiation therapy. contribution to the conference ORP3, Universityof Kaiserslautern, 2003.

[3] T. Benoist and F. Chauvet. Complexity of some FPP related problems.Technical report, Bouygues’ e–lab, 2001.

[4] N. Boland, H.W. Hamacher, and F. Lenzen. Minimizing beam-on timein cancer radiation treatment using multileaf collimators. NETWORKS,43(4):226–240, 2004.

[5] I.M. Bomze and W. Grossmann. Optimierung – Theorie und Algorith-men. BI–Wissenschaftsverlag, Mannheim, 1993.

[6] T.R. Bortfeld, D.L. Kahler, T.J. Waldron, and A.L. Boyer. X–ray fieldcompensation with multileaf collimators. Int. J. Radiat. Oncol. Biol.Phys., 28:723–730, 1994.

[7] A.L. Boyer and C.Y. Yu. Intensity–modulated radiation therapy withdynamic multileaf collimators. Semin. Radiat. Oncol., 9:48–59, 1999.

[8] J. Dai and Y. Zhu. Minimizing the number of segments in a deliv-ery sequence for intensity–modulated radiation therapy with a multileafcollimator. Med. Phys., 28:2113–2120, 2001.

[9] M.L.P. Dirkx, B.J.M. Heijmen, and J.P.C. van Santvoort. Leaf trajec-tory calculation for dynamic multileaf collimation to realize optimizedfluence profiles. Phys. Med. Biol., 43(8):1171–1184, 1998.

[10] K. Engel. A new algorithm for optimal multileaf collimator field seg-mentation. Preprint 03/5, Fachbereich Mathematik, Uni Rostock, underrevision for Discr. Appl. Math., 2003.

102 Bibliography

[11] J.M. Galvin, X.G. Chen, and R.M. Smith. Combining multileaf fieldsto modulate fluence distributions. Int. J. Radiat. Oncol. Biol. Phys.,27:697–705, 1993.

[12] M.R. Garey and D.S. Johnson. Computers and intractability, a guide tothe theory of NP–completeness. W.H. Freeman, 1979.

[13] A. Holder and B. Salter. A tutorial on radiation oncology and opti-mization. Technical Report 86, Trinity University, San Antonio, Texas,2004.

[14] D. Jungnickel. Graphen, Netzwerke und Algorithmen. BI–Wissenschaftsverlag, Mannheim, 1994.

[15] T. Kalinowski. Realization of intensity modulated radiation fields us-ing multileaf collimators. In R. Ahlswede, with the assistance of L.Baumer, and N. Cai, editors, General Theory of Information Transferand Combinatorics. Shannon Foundation, to be published 2004. Reporton a Research Project at the ZIF (Center of interdisciplinary research)in Bielefeld Oct. 1, 2002 – August 31, 2003.

[16] S. Kamath, S. Sahni, J. Li, J. Palta, and S. Ranka. Leaf sequenc-ing algorithms for segmented multileaf collimation. Phys. Med. Biol.,48(3):307–324, 2003.

[17] S. Kamath, S. Sahni, J. Palta, and S. Ranka. Algorithms for optimalsequencing of dynamic multileaf collimators. Phys. Med. Biol., 49(1):33–54, 2004.

[18] M. Langer, V. Thai, and L. Papiez. Improved leaf sequencing reducessegments of monitor units needed to deliver IMRT using multileaf colli-mators. Med. Phys., 28:2450–2458, 2001.

[19] L. Ma, A.L. Boyer, L. Xing, and C.M. Ma. An optimized leaf settingalgorithm for beam intensity modulation using dynamic multileaf colli-mators. Phys. Med. Biol., 43(6):1629–1643, 1998.

[20] W. Que. Comparison of algorithms for multileaf collimator field seg-mentation. Med. Phys., 26:2390–2396, 1999.

[21] I.I. Rosen, R.G. Lane, S.M. Morill, and J.A. Belli. Treatment planoptimization using linear programming. Med. Phys., 18:141–152, 1991.

Bibliography 103

[22] D. Shepard, M. Ferris, G. Olivera, and T. Mackie. Optimizing the de-livery of radiation therapy to cancer patients. SIAM Review, 41(4):721–744, 1999.

[23] R.A.C. Siochi. Minimizing static intensity modulation delivery timeusing an intensity solid paradigm. Int. J. Radiat. Oncol. Biol. Phys.,43:671–680, 1999.

[24] S.V. Spirou and C.S. Chui. Generation of arbitrary intensity profiles bydynamic jaws or multileaf collimators. Med. Phys., 21:1031–1041, 1994.

[25] R. Svensson, P. Kallman, and A. Brahme. An analytical solution for thedynamic control of multileaf collimators. Phys. Med. Biol., 39(1):37–61,1994.

[26] J.P.C. van Santvoort and B.J.M. Heijmen. Dynamic multileaf collima-tion without ’tongue-and-groove’ underdosage effects. Phys. Med. Biol.,41(10):2091–2105, 1996.

[27] S. Webb, T. Bortfeld, J. Stein, and D. Convery. The effect of stair–step leaf transmission on the ’tongue–and–groove problem’ in dynamicradiotherapy with a multileaf collimator. Phys. Med. Biol., 42(3):595–602, 1997.

[28] P. Xia and L. Verhey. Multileaf collimator leaf–sequencing algorithm forintensity modulated beams with multiple static segments. Med. Phys.,25:1424–1434, 1998.

Hiermit erklare ich, dass ich die vorliegende Arbeit selbstandig und ohnefremde Hilfe verfaßt, andere als die von mir angegebenen Quellen nicht be-nutzt und die den benutzten Werken wortlich oder inhaltlich entnommenenStellen als solche kenntlich gemacht habe.

Thomas Kalinowski

LEBENSLAUF

Thomas Kalinowskigeb. am 27.11.1980 in PasewalkFamilienstand: ledig

Schulbildung

09/86 – 08/92 : Polytechnische Oberschule in Grambow09/92 – 08/94 : Gymnasium Locknitz09/94 – 07/98 : Abitur an der CJD–Jugenddorf–

Christophorusschule Rostock

Studien und akademische Ausbildung

10/99 – 09/03 : Mathematikstudium an der Universitat Rostock(Nebenfach Informatik),Stipendiat der Studienstiftung des deutschen Volkes

09/01 – 04/02 : Studienaufenthalt als ERASMUS–Stipendiatan der University of Glasgow (Schottland)

09/03 : Abschluß des Diplomstudiums,Thema der Diplomarbeit: Realization of intensitymodulated radiation fields using multileaf collimators

Berufliche Tatigkeit

Seit 10/03 : Wissenschaftlicher Mitarbeiter am FachbereichMathematik der Universitat Rostock

Schulerforderung

Seit 09/00 : Mitglied im Mathematik–Olympiaden e.V.09/00 – 07/01 : Leitung einer Mathematik–AG fur

Schuler der Klassen 9–1009/02 – 07/03 : Leitung einer Mathematik–AG fur

Schuler der Klassen 11–12

Thesen zur Dissertationvon Thomas Kalinowski

1. In der Krebsbehandlung werden Mehrlamellenkollimatoren eingesetztum intensitatsmodulierte Bestrahlungsfelder zu erzeugen. Ein Schrittbei der Erstellung eines Behandlungsplans ist die Realisierung einer ge-gebenen Fluenzverteilung mit einer geringen Gesamtbestrahlungsdauerund einer kleinen Anzahl von homogenen Feldern. Die moglichen La-mellenpositionen des Kollimators konnen durch (0, 1)−Matrizen, soge-nannte Segmente, beschrieben werden. Damit lasst sich das Problemals Suche nach einer Segmentierung einer nichtnegativen ganzzahligenm × n–Matrix A formulieren, d.h. nach einer Darstellung

A =

k∑

i=1

uiSi

von A als positive Linearkombination von Segmenten S1, . . . , Sk.

2. In vielen Mehrlamellenkollimatoren verbietet die Lamellenkollisionsbe-dingung (interleaf collision constraint) das Uberlappen von gegenuber-liegenden Lamellen in benachbarten Zeilen. Fur diese Kollimatoren sinddie m × n–Segmente S = (si,j) charakterisiert durch die Existenz vonganzen Zahlen li, ri (i = 1, . . . , m), so dass

li ≤ ri + 1 (i ∈ [m]),

si,j =

{

1 falls li ≤ j ≤ ri

0 sonst(i ∈ [m], j ∈ [n]),

li ≤ ri+1 + 1, ri ≥ li+1 − 1 (i ∈ [m − 1]).

3. Die Minimierung der Gesamtbestrahlungsdauer wird als lineares Opti-mierungsproblem formuliert:

(P )

∑

S∈S

f(S) → min

f(S) ≥ 0 ∀S ∈ S,

∑

S∈S:(i,j)∈S

f(S) = ai,j ∀(i, j) ∈ V.

Hier ist S die Menge der Segmente, wobei ein Segment S mit der Mengeder Paare (i, j) ∈ [m] × [n] mit si,j = 1 identifiziert wird.

4. Mit Hilfe des dualen Problems

(D)

∑

(i,j)∈[m]×[n]

ai,jg(i, j) → max

∑

(i,j)∈S

g(i, j) ≤ 1 ∀S ∈ S

wird der minimale Wert der Zielfunktion im Problem (P ) als maximalesGewicht eines 0 − 1−Weges in folgendem gewichteten Digraphen G =(V, E) charakterisiert:

V = [m] × [n] ∪ {0, 1},

E = E1 ∪ E2 mit

E1 = {(0, (i, 1)) : i ∈ [m]} ∪ {((i, n), 1) : i ∈ [m]},

E2 = {((i, j), (i′, j + 1)) : i, i′ ∈ [m], j ∈ [n − 1]},

mit der Gewichtsfunktion

δ(0, (i, 1)) = ai,1 (i ∈ [m]),

δ((i, n), 1) = 0 (i ∈ [m]),

δ((i, j), (i, j + 1)) = max{0, ai,j+1 − ai,j} (i ∈ [m], j ∈ [n − 1]),

δ((i, j), (i′, j + 1)) = max{0, ai′,j+1 − ai′,j} −i′−1∑

k=i

ak,j

(i, i′ ∈ [m], i < i′, j ∈ [n − 1]),

δ((i, j), (i′, j + 1)) = max{0, ai′,j+1 − ai′,j} −i∑

k=i′+1

ak,j

(i, i′ ∈ [m], i > i′, j ∈ [n − 1]).

5. Aus dem Beweis der Min–Max–Aussage in These 4 wird ein effizienterAlgorithmus zur Bestimmung einer optimalen Losung des Problems(P ) abgeleitet.

6. Sei c(A) die minimale Koeffizientensumme einer Segmentierung von A.Ein Paar (u, S) aus einer naturlichen Zahl u und einem Segment S heißezulassig, falls A′ := A − uS nichtnegativ ist und c(A′) = c(A) − u gilt.Zur Konstruktion einer Segmentierung mit minimaler Koeffizienten-summe und kleiner Segmentanzahl wird die folgende Greedy–Strategieverwendet.

1. Bestimme ein zulassiges Paar (u, S) mit maximalem u.

2. Setze A := A − uS.

3. Falls A 6= 0, gehe zu Schritt 1.

7. Die Minimierung der Segmentanzahl ist bereits fur m = 1 streng NP–schwer. Es wird jedoch bewiesen, dass dieses Problem in Linearzeitgelost werden kann, wenn die Matrixeintrage durch eine Konstante be-schrankt sind.

8. Fur m > 1 wird ein Algorithmus der Komplexitat O(mn2L+2) angege-ben, der fur einen Kollimator ohne Lamellenkollisionsbedingung eineSegmentierung mit minimaler Gesamtbestrahlungsdauer und minima-ler Segmentanzahl bestimmt. Hier ist L = max{ai,j : (i, j) ∈ [m]×[n]}.

9. Eine deutliche Reduzierung der Gesamtbestrahlungsdauer ist moglich,wenn der Kollimator zwischen der Verabreichung der einzelnen Felderum 90◦ gedreht werden darf. Hierzu wird eine Heuristik implementiertund getestet.

10. Eine weitere Reduzierung der notwendigen Gesamtbestrahlungsdauerkann durch die Verwendung von zwei zueinander orthogonal angeord-neten Kollimatoren erreicht werden, wie wiederum mittels Implemen-tierung eines heuristischen Algorithmus gezeigt wird.

optimal multileaf collimator eld segmentationkalinowski/papers/tk_diss.pdf · 2013. 3. 6. ·...

Documents