Graph cutting and its application to biological data

2.1 Maximal Flow and Minimal Cut

Maximal flow is the flow with maximal value, [1, 6, 12]. Maximal flow is obtained step by step saturating darts and augmenting paths in the network, so that the precise conditions are satisfied. Let A be a subset V(G) and also s ∈ A, t ∈ V(G) − A. The set < A, V(G) − A > denotes all darts emanating from the vertex v ∈ A to w ∈ V(G)−A. The value f (< A, V(G) − A >) presents the sum of flows on darts from the set < A, V(G)−A >. Capacity c(< A, V(G)−A >) denotes the sum of capacity of all darts from the set (< A, V(G)−A > ). For each size of flow it holds:

The set of darts < A, V(G) − A > is the cut defined by the set A. It means that the cut separates source from the sink and c(< A, V(G) − A >) is the capacity of the cut. Ford-Fulkerson theorem [6] says that the size of maximal flow from the source s to the sink t is equal the sum of capacity of the edges of minimal cut [2, 4]. The reason is that maximal flow saturates all edges belonging to the set < A, V(G) − A >, which belong to minimal cut:

The sequence of darts is called an augmenting path P, if it creates a path and we do not consider orientation of darts. The augmenting path can contain direct as well as reverse darts. We denote the dart as direct, if its orientation is the same like the orientation of the augmenting path. To every dart of this augmenting path there are assigned weights, which are called "reserves" r(u, v). The value of the "reserve" depends on the orientation of the dart:

We will choose from all "reserves" the smallest one and call it the reserve of the augmenting path r(P) and if its value is positive, then we call this value "increasing augmenting path". The increasing s − t augmenting path is a path starting in the vertex s and ending in the vertex t. If there is no increasing s−t augmenting path, then there is a maximal flow in the network. If there is still some flow, in the next step we increase the flow in the network and for the same value of the network there holds:

For finding maximal flow in the network and its minimal cut there are usually used two types of algorithm based on two approaches. One is a family of augmenting paths algorithms [4, 6, 13], specially Ford - Fulkerson algorithm [6] another one is Push relabel algorithms [7, 11, 13].

2.2 Ford-Fulkersonov algorithm

It is the basic algorithm for finding maximal flow in the network. It was described in 1956 by L. R. Ford Jr. and D. R. Fulkerson [6]. Searching of maximal flow in the network is by iterations, finding for augmenting s − t paths and increasing of the flow W(f ). Augmenting paths are possible be defined and searched by different approaches. Ford - Fulkerson algorithm is general, it is not direct way how to find them.

If there is a solution, that we reach t, it means we found maximal flow in the network. All vertices, which were not discovered and reached in this last iteration, will create a set A. Then the number of darts 〈A, V(G) − A〉 will create maximal flow with the value c(〈A, V(G) − A〉) = W(f ). We discuss Edmond’s Karp algorithm [1, 12] and its modification in the next chapter.

2.3 Edmonds-Karp algorithm

Each dart (v_iv_j) has assigned two value, the flow f and the capacity c, as in Chapter 2.1. The flow f is an actual amount passing the edge and the capacity c is a maximum value of the flow f that can pass the edge.

In Figure 1 we present the "flow chart diagram" for Edmonds-Karp algorithm as we implemented it in our program. In the queue Q we add vertices we have found in breadth-first search. We also remove vertex from Q, when no darts satisfy certain conditions (2.6). Specifically, the queue Q contains indexes of vertices so it is an array of integers. Each vertex v has its own label l(v). This label shows the index of a vertex, which we came to v from. The label for source vertex s contains value l(s) = 0, because we start from this vertex and s is known at the beginning of the breadth-first search. We set labels for other vertices as v_i as l(v_i) = −1. The queue Q contains just s at position Q[0] at the beginning.

Figure 1

The flow chart diagram: The flow chart diagram for the Edmonds-Karp algorithm. The Breadth-first search algorithm is marked as green, the back-tracing of augmenting s-t path is marked pink.

Then we pick the vertex v_i at position Q[0] and check two types of dart containing the vertex v_i:

If we find an edge satisfying these conditions, then we add vertex v_j at the end of Q and set l(v_j) = v_i. If v_j = t , we found sink t and it means we can increase flow on an augmenting path.

When we look at label l(t) of t we get vertex v_j. Edge (v_j t) will be included in the augmenting path P. Then we look at label l(v_j) and so on. In the end we find vertex with label l(v_j) = 0 what means that v_j is our source s. Such way we find all edges creating the augmenting path P. On these edges we have to change flow. Each edge has own reserve r(v_iv_j) - amount of the flow possible to be added on this edge. The value r(P) represents the reserve of the augmenting path P. It is defined as: r ( P ) = min ( v i v j ) ∈ P r ( v i v j ) .

Then the flow in the network changes in the following way:

After checking all edges meeting mentioned conditions we can remove vertex v_i from Q and repeat these steps for vertex situated on the next place in Q. When Q no longer contains a vertex, then we cannot find the augmenting path and therefore there is now actually maximum flow in the network.

Edmonds-Karp algorithm finds the shortest augmenting paths in the network. Our implementation uses exactly one augmenting path in each iteration. Therefore, the computantional complexity of these types is used for solving maximum flow-minimum cut problem is O(VE²), where V is the number of vertices in the network and E is the number of edges.

For optimization we used one modification, which increases the speed of the algorithm. We saturate s − t augmenting paths containing exactly two darts (sv_i) and (v_i t) before starting the Edmonds-Karp algorithm. We compare capacities of these darts and set the flow f on them as: f (sv_i) = f (v_i t) = min(c(sv), c(vt)). Therefore, one of these edges will be saturated. We apply this on every vertex v_i. Consequently, there will be augmenting paths of length, at least three in the network. Then the Edmonds-Karp algorithm finds maximum flow.

3.1 Transformation of real data into network

The data comes from Comenius University. They are actually examining a production of iron in cells of living organisms. The theory states that in a certain part of the brain with Purkyne cells more iron cells appear than in another part. Our task is to confirm this statement using our segmentation techniques. Therefore we acquired grayscale and coloured images obtained by microscope. These images contain brain cells of rabbits after irradiation with the GSM the signal.

We created a program satisfying the conditions for this specific problem and working with these specific sorts of images. We use segmentation, histogram counting and percentual ratio of segmented area to confirm this statement.

An important remark is that the obtained segmentation by the Edmonds-Karp algorithm is not enough for obtaining requested results. The mentioned algorithm gives us just very basic results. But we needed to create other new algorithms and new mathematical tools for obtaining "good" segmentation, as well as detection in which region more cells appear. For this we needed to develop other algorithms and counting histograms as well as vector analysis.

First, we need to transform a biomedical image into the network. We show the transformation on a grayscale image containing pixels p with intensities I_p in the range of 0 to 255. This approach was first presented by Boykov in [3, 13].

Single pixels of the image become vertices and boundaries of pixels become edges in the network. Edges represent relationships between neighboring pixels. We consider four pixels as neighbors to pixel p. In that way we receive an unordered graph. Each network contains two special vertices source s and sink t. We add these two in our graph and connect them through arrows with each vertex in the graph. The orientation of these darts is su and ut and they are called t−links. An orientation of remaining edges has to satisfy the following condition. When it is possible to get from vertex u to v, then there has to be an opportunity to get from v to vertex u. Therefore we replace the original edge with two darts with opposite orientation and the same capacity. These darts are called n − links. The capacities of darts are based on the difference of intensities of neighboring pixels.

At the beginning we have to select some source pixels. We want to segment an image into two regions - an object and background region. We choose pixels best describing the object and these pixel then will belong to the object set O. Pixels that best describe background will belong to the background set B. At least one pixel must be chosen to each set, but no pixel can belong to both sets. For us, the pixels in the first set will represent iron in cells. To source pixels there belong certain vertices in the network - source vertices.

The minimum s − t cut algorithm divides our network into two sets of vertices. The vertices v assigned to set O remain reachable from the source s, so there exists the s −v path while there does not exists v − t path, because the reserve (vt) = 0. Similarly for the vertices v in set B there exists v − t path and does not exists s − v path.

We set capacities within a simple way. Let have an image with pixel intensities I_p. Then find maximum value M and minimum value L of intensity in the image. The difference between maximum and minimum value is D = M − L. The average value of source pixels in set O is I s ¯ and in set B it is I t ¯ . The capacities of t − links will be defined the following way:

These t − link capacities are called regional property and represent a penalty for incorrect assignment of pixel to set O or B. The coefficient λ sets the weight of this regional property and affects a result of segmentation. A smaller value leads to objects with smoother boundaries.

Next we will ensure that selected source pixels remain in the set we selected at the beginning. For the t − links pixels satisfies a special property: the vertex v ∈ O has to be always reachable from the source s and sink t cannot be reached from v. Therefore the capacities for v must satisfy a special property:

For v ∈ B it is opposite:

The last thing to solve is a capacity of n − links. It can be easily defined as:

The n −link capacities are called "boundary property" and represents penalty for discontinuity between neighboring pixels. We get large value pixels that have similar color and smaller value for pixels of different colors. These two properties use Boykov and Jolly [3] for determining the cost function, which minimization is equivalent to final segmentation.

3.2 Finding the minimum s–t cut

The theorem of Ford and Fulkerson states that a maximum flow from s to t saturates a set of edges corresponding to a minimum s − t cut. Moreover, the cost of the minimum cut is equal to the value of the maximum flow. Thus, the minimum cut and the maximum flow are dual problems. By finding the maximum flow we receive the minimum cut. The Ford–Fulkerson algorithm is a basic algorithm that solves the maximum flow problem. It belongs to the group of augmenting path algorithms. The augmenting path is a path from s to t that has available capacity on all its edges. The main idea is that we start with zero flow and as long as an augmenting path exists, we push the maximum possible flow along it – we saturate it. By saturating each augmenting path the total flow in the graph is increased. Once no an augmenting path exists, the maximum flow is reached and the algorithm terminates.

3.3 Marking procedure

This procedure performs a search for an augmenting path. Let f (e) denote current flow in the edge e, c(e) is the capacity of e.

1. Mark source node s, all other nodes are not marked.

2. (exists e = (p, q) ∈ E, p is marked, q is not marked, f (e) < c(e)) mark→ q

3. (exists e = (p, q) ∈ E, pi is not marked, q is marked, −c(e) < f (e))→mark p

4. If t is marked, an augmenting path exists. STOP.

   If no node has been marked since last loop iteration, no augmenting path exists. STOP.

5. GOTO 2.

The steps 2 and 3 can be implemented in various ways. The Edmonds–Karp algorithm uses breadth-first search, which results in finding the shortest augmenting path in each step. For the reconstruction of the augmenting path we need to record which node caused the marking and in which direction; step 2 “marks forwards”, step 3 “marks backwards”. Then, the augmenting path can be reconstructed by tracing back from t to s.

3.4 Modifications of the algorithm

Here we present two major changes in the implementation of the Edmonds–Karp algorithm. The first one is easy to perform, significantly reduces running time with graphs based on real data and does not change the result compared to non-modified algorithm. Second modification further improves the speed of the algorithm. While the resulting segmentation is no longer the same, we show that it corresponds better to the expectations in our application.

3.5 Improvement of the marking procedure

The original marking procedure stops right after the sink has been marked. Then we are able to reconstruct only one augmenting path. However, if we allow multiple marks of the sink, more than one augmenting path can be found at once after performing one marking procedure. As a breadth-first search is used, all these paths found in one step have the same length. Because the marking procedure is a very time consuming operation, this modification reduces the total count of marking procedures before the algorithm terminates. Note that now not all of these paths found in one step will actually increase the total flow in the graph. The reason is that multiple paths can share some of their edges. There is a chance that by saturating one path we also saturate others. In the worst case only one of these paths will increase the total flow while others will not; this modification would not bring any advantage in this case. But tests with graphs constructed in our application show that the implementation modified in this way completes the segmentation process about ten times faster than the original algorithm while the results of both implementations are identical.

3.6 Further modifications

The original algorithm used for image segmentation starts by marking the source. In the next step it marks all nodes connected to the source by a non-saturated edge, which often means vast majority of the set P. So the algorithm spends a lot of time in those parts of the image that are far from expected result. Following alternation of the algorithm significantly reduces this problem. First, we define a new function that will help to describe the alternated algorithm. It is clear that for all p ∈ P a path only consisting of two t-links (s, p), (p, t) is the augmenting path if neither of these links is saturated. The saturation of all these paths will always be the first step of the algorithm and leads to the saturation of at least one t-link in each path. At most one t-link to each node p preserves free capacity. For all p ∈ P we define a new function R(p) : P −→ R as follows:

This function fully describes the state of all t − links after previously introduced first step of the algorithm. A positive value means that the t-link from the source has some remaining capacity R(p), while the t-link to the sink is fully saturated. A negative value indicates that the t-link from the source is saturated, while the t-link to the sink preserves capacity |R(p)|.R(p) = 0 means that both t-links got saturated. We evaluate this function during the construction of the graph. Now we propose to start the marking procedure from all object seeds instead of the source s. The procedure will end by marking any p ∈ P so that R(p) < 0. Because the starting node is connected to the source by a non-saturable edge (with infinite capacity) and the node at the end is connected to the sink by a nonsaturated edge, we have found one augmenting path from the source to the sink. However, this path never includes t-links from the source with a finite capacity. These links would be meaningless. Therefore we also modify the procedure reconstructing the augmenting paths. Instead of one path we reconstruct more of them. Each of them starts in the source, follows a different t-link, then they unite one after another and finally they reach the sink together using common t-link. They form a system of paths that can increase the total flow at most by the value −R(p), where p was the last node marked at the end of the marking procedure. While following the path back we bring flow.

3.7 Reconstruction of the system of paths

Results of the algorithm modified in this way will differ from the unchanged algorithm. In general, we perform the marking in a smaller area, thus resultant object segment will be smaller, especially in the case of the result consisting of several disjoint areas. The resultant segmentation of the original algorithm can contain more disjoint segments even if there was only one object seed. Our modified algorithm requires at least one object seed in every disjoint segment, otherwise marking would neither start in this segment nor reaches it. This requirement is not a limitation in our application, our goal is to segment only one object.

Note: CPU after all this optimization is reduced to O(V²E) and for some dense graphs it can be reduced to O(V³).

3.8 The program functionality and implementation

Basic outline of functions needed to implement the consist of an image data loading, creating a network and its representation, finding a maximum flow and finally creating a segmented image based on reserves remaining on darts in the network after the last step of the algorithm. But just to see the resulting image is not enough for us and we need to get numerical results to confirm the theory about gathering the iron cells in the region close to the Purkyne cells. We suggest our own methods to calculate amount of segmented pixels in the image. The first method better visualizes concentration of cells and gives us a grayscale result on square net of the image with brighter squares representing higher concentration. Next method uses histogram of image describing numbers of pixels segmented on a user defined line and lines parallel to it. The line is named as the direction vector. Consequently percentage of segmented pixels in a user defined region is counted based on this histogram. Last method implemented is our interpretation of region-growing method to determine the amount of iron cells in the image, because one cell can be formed from more pixels. Detailed description of them is listed in following rows. These methods and our program was designed for the specific purpose of counting iron cells and determining the region with more cells.

We implemented the maximum flow-minimum cut method for solving segmentation problem in Microsoft Visual Studio environment in C++ language. We find minimum cut with the Edmonds-Karp algorithm described above. Our program has a graphic user interface with more control elements and two main windows as shown in Figure 2. The user loads a gray scale or color image in the left window. Then the user can choose the source pixel or pixels of the object and the background, a direction vector and a threshold value. One pixel in each set will be sufficient for a maximum flow-minimum cut algorithm, but for reaching more accurate results we recommend to choose more source pixels. The direction vector determines the direction of counting of segmented pixels. Threshold value sets the upper threshold in our case, because we have segment images with bright background and dark objects and bright pixels have higher values of intensity. The only parameter in our program is λ-value, which has an influence on the result of segmentation. The last thing to set is a length of the side of the square. There the image will be divided into squares colored according to the number of the segmented pixels as one of the results. All these settings are situated in the upper part of the program highlighted in blue colour.

Figure 2

Program user interface: The part highlighted in the blue colour contains initialization settings and tools. The part highlighted in the green colour serves first for clicking on the object and the background pixels or for choosing the direction vector. Later it displays image with segmentation and regions selected in the histogram. In the orange rectangle histogram will be displayed. Numerical results will be shown in the purple rectangle. And finally the region with brown boundary defines the histogram processing part and after segmentation displays the percentage of segmented pixels in the selected region.

Since we are doing a binary segmentation, a result will be two sets of pixels - object and background pixels. After finishing the algorithm, more information about the segmentation and results in form of images will be displayed, see Figure 3. The most important information is displayed under labels Area, Percent and Cells in the right part of the program highlighted in the purple colour in the Figure 2. A total number of segmented pixels from the whole image is displayed in a row labeled as Area, how many percents it is from all the pixels in the image is displayed in a row labeled as Percent and in a row labeled as Cells we can find number of independent segmented objects. When two neighbouring pixels (squares with common edge, not just the vertex) are marked in process of segmentation as object pixels, they will belong to the same object. Then two regions that touch just in diagonal pixels are considered as two objects.

Figure 3

Displaying results in our program: The part in the green rectangle shows segmented iron cells and user selected region. Histogram with two columns selected is displayed in the orange rectangle. Two red columns represent two green lines in the image in the left window. Numerical results are in the purple rectangle. The percentage of segmented pixels in the region between the two green lines is written in the brown rectangle.

An image in the left window highlighted in the green colour will be redrawn and will show highlighted objects founded in the process of segmentation. If the user has selected a direction vector, the histogram will be displayed in the right window highlighted in the orange colour. This histogram represents absolute numbers of segmented pixels on lines parallel to the direction vector. It means the program passes the whole image line after line in the selected direction so that each pixel belongs to exactly one line and count absolute numbers of segmented pixels on each line. These numbers are shown on the y-axis of the histogram, on the x-axis there are indexes of lines. Then the histogram processing part of program highlighted in brown under the image of histogram becomes active. The user can select two columns in histogram that represent number of segmented pixels on two lines in the resulting image. After that the program returns number of segmented pixels in the region between these two lines in the image compared to all segmented pixels in the image. This result is shown in a percentage form and can acquires values from range 0 to 1. At the same time the lines forming boundaries of a selected region are displayed in image in the left window (highlighted in the green colour). When the user is satisfied with the results he/she can save these images or he/she can select other lines in the histogram or a different direction vector.

The results, saved automatically, belongs to the black-white image showing white objects on a black background, while boundaries of these objects can be also be highlighted in this image, see Figure 7. The next automatically saved result shows highlighted objects in an original image. The special image is created as one of the results. This grayscale image of the size the original image shows square net with a square edge of a size chosen by the user. Each square has an intensity ranging 0 to 255 according to the number of segmented pixels in the original image in region falling to this square. This number is normalized what means that in a region of the original image in the region falling to this square. This number is normalized which means that in a region of the original image belonging to the white square there occurs the most segmented pixel. The darker a square is the least segmented pixels to appear in this region. The result in this form is situated in the Figure 8.

For representing our image as a network (transformation was mentioned in the section 3.1) we use structures and classes. We have classes named Vertex, Arrow (a dart), Seed and a structure named Histogram. All of the vertices and darts in the network are grouped in arrays, so every vertex and dart has a unique index. The class Vertex consists of information such as indexes of incident darts and at the same these darts are divided into arrays of direct and reverse darts. This step increases memory requirements, but also increases the speed of searching of the augmenting path. We also need an index of previous vertex during the reconstructiong of the augmenting path, so this information will be stored too. The class Arrow (darts) contains two indexes - vertex where a dart starts and ends, a capacity and a reserve of a dart. The program contains two objects of the class Seed. The one is an array of source vertices of an object and the second is an array of source vertices of the background.

The structure Histogram will contain information about starting and ending point of a direction vector and every line on which histogram is counted. Also the number of segmented pixels on every line is saved.

After the network is constructed, that means all objects of classes Vertex and Arrow (a dart) are initialized, source pixels are set, capacities of darts are set, maximum flow is counted. When working with color images, capacities are set three times for each color channel separately and maximum flow is also counted three times. After some experiments with selecting the final segmentation, we select the second option below, because it is better for our type of images containing a background of mainly pink, red and white color. The first one is to create a final segmentation from pixels that are segmented at least on two color channels. The second one is to select just segmentation done on the red channel. There are other types of color image processing, too. But we coded just this option for our purpose.

For increasing the speed of the Edmonds-Karp algorithm we use optimization mentioned in Section 2.3. After the algorithm terminates, we get final the segmentation 4

Figure 4

Results: Final segmentation and corresponding histogram