E Koorambas
email: elias.koor@gmail.com

The new interpretation of arithmetic operation symbols

Keywords: MSC number(s):00A30, 03A05, 20B30 Arithmetic Operation Symbols, Permutation Group

We introduce the permutation group of arithmetic operations symbols by getting the permutations of all the common arithmetic operations symbols, with keeping the brackets out of ordering. We find 6 ways of doing the arithmetic operations. Therefore the output of any mathematical formulas depends on which one element of the arithmetical permutation group we work on. We find invariants by the reordering of the arithmetic operation x+y, xy. Working with the irreducible representation of the permutation arithmetic symbols group we define new arithmetic structures called arithmetic particles symbols.

It is certainly both natural and useful to associate positive numbers with length, areas, and volumes. This is an ideal that goes back at least to the ancient Greeks [START_REF] Kline | Mathematical Thought from Ancient to Modern Times[END_REF], and we use it whenever we draw a graph, a histogram or almost any kind of chart. Following Howie.J.M . [START_REF] Howie | Real Analysis Springer -Verag[END_REF] , to take a simple example, if we have a triangle ABC, right angles at A, and with AB and AC of length 1, then by the celebrated theorem of Pythagoras the length BC is a number z (1.1) But how do we find this result? In any set of numbers we follow a standard algorithm by doing the calculation which children know at the first classes of school namely as ordering of operations [START_REF] Hazewinkel | Encyclopedia of Mathematics[END_REF] . 1.0 Axiom: Let two or more numbers x,y,………,  where A is any set of numbers, and let a mathematical formula K made by those numbers, and by the normal operations +, -, / and the brackets ( ) To find K we follow the ordering: A1. The brackets are going first of all the others A2. Second we are doing the squares A3. Third we are doing the multiplications and the divisions A4. Finally we are doing the sums That is the way which gives the result 2 z  in equation (1.1). More formally this equation is the inner product in the set of real numbers R where (1.1) defines the positive length and thus the space which is the Euclidian one for this case. . Make this into the group [START_REF] Weiss | Group and Geometry[END_REF] for some , ab .There are six elements in 3 S . They are:

(1)( 2 )(3) E  , ( 13
)(12) A  , ( 12
)(13) B  , (1)( 2 3) C  , ( 2 )(1 3) D  , (1 2 )(3) F 
This formulation is very important as every element {1,2,3} is denoted by the operation, hence the permutation group 3 S making 5 new ways of doing the arithmetic operation gives 5 new arithmetic's (13)(12), (12) (13),(1)(23), (2)(13),(12)(3) symbols. 2.1. Proposition. The mathematical outcomes from any specific formula K known so far, just belongs to one element of permutation 3 S arithmetic operation symbols group, (1) (2) (3) element, coming from the concept of 3 S arithmetic and respect the 1.0 Axiom. Proof: here we mean that we respect the specific interpretation of arithmetic operations symbols given by the (1) (2) (3) element 3 S arithmetic operation symbols group.

An example of the effect on a mathematical formula by the reordering of the arithmetic operations symbols is the following: let the function

2 ( ) 1, f z z z   
(2.1) And let that respects the arithmetic 3 S group, in the case of the (123) element, the well known as ordinary. This function is said to be conformal in all it's range except the zero [START_REF] Churchill | Complex Variables and Applications[END_REF] . Let the (321) element that also belongs to the arithmetic 3 S group, the function () fz, is given by

( ) ( 1), f z z z z    (2.
2) We see that the same function changes under the reordering of arithmetic operation symbols, hence () fz is conformal respect in all the range except the point 1 / 2  . With this example we see that a property of a function () fz which is the range of where is conformal or not conformal, changes under the action of permutation of arithmetic operation. The possible partitions of the The partition  

3 corresponds only to one Young tableau namely

1 2 3
The symmetrizer corresponding to the tableau is The antisymmetrizer corresponding to the tableau is

(12) (13) (23) (123) (132) A ritm etic E P            (3.1) 3
(12) (13) (23) (123) (132) A ritm etic E Q            (3.2
             (4.1) Its dimension is 3 1 (1 , ) ( 1)( 2) 3 6 n N n n n n        (4.2)
The symmetric tensor corresponding

  1 (12) (13) (23) (123) (132)        (4.3) Its dimension is 1 ( ) ( 1)( 2) 6 N n n n n    (4.4)
For n = 3, we have the following 10 components

111 1 1 2 1 1 3 1 2 2 1 2 3 1 3 3 222 2 2 3 2 3 3 333
The mix symmetric tensor corresponding to arithmetic anion , thus (5.3) is also corresponding sequence in the set X spaces ensures that the above theorem exists.

13 12 12 (1 )(1 ) 3
PP         (4.5) 12 13 13 (1 )(1 ) 2 PP         (4.6) Its dimension is found from the product         11 1 111 21    (4.7) This implies   2 1 ( 2,1 , ) ( 1) 3

Definition:

We defined invariant from the reordering point of view the relations which are not affected from the reordering of arithmetic operations symbols, like relations with only one arithmetic operations symbols , ,...... x y x y  (5.4) By looking the usually arithmetic operations symbols from above, as elements of an ordering set making this set to the symmetric group 3

S and keeping the brackets out of ordering. By 3 S arithmetical group we find the 6 different ways of doing the arithmetical operations yield to different results of mathematical formulas. As we work with the irreducible representation of 3 S arithmetic symbols group we get the following case: the symmetric arithmetical operations symbols (arithmetic bosons) and the antisymmetric of arithmetic operations symbols (arithmetic fermions) and the two mix arithmetic symbols (arithmetic anions).

3 A

 3

  .1 Definition: we define the algebraically structure referring of the partition  

) 3 . 2

 32 Definition: we define the algebraically structure referring to the partition   represent by the mix symmetrizer operators (3.6), (3.4) which is the result of combination of doing the arithmetic operations. The possible Young diagrams are The antisymmetric tensor corresponding 1 (12) (13) (23) (123) (132)

3 )

 3 Let a sequence, or series defined in the space 123 , there are series in the space 321 ? If exists a sequence or series, exists always? The answer is yes for both statements. The reordering map f transfers from sequence or series to sequence or series. Let work with example. Let the series can exist in the set 123


André Körbes
email: korbes@dca.fee.unicamp.br

Giovani Bernardes Vitor
email: giovani@fem.unicamp.br

Roberto De Alencar Lotufo

Vaqueiro Ferreira

Advances on Watershed Processing on GPU Architecture

This paper presents an overview on the advances of watershed processing algorithms executed on GPU architecture. The programming model, memory hierarchy and restrictions are discussed, and its inluence on image processing algorithms detailed. The recently proposed algorithms of watershed transform for GPU computation are examined and briely described. Its implementations are analyzed in depth and evaluations are made to compare them both on the GPU, against a CPU version and on two diferent GPU cards.

Introduction

The watershed transform is a well known image segmentation tool, used to compute connected components, generally out of the gradient of the input image. Several algorithms exist to process the transform sequentially [START_REF] Körbes | Analysis of the watershed algorithms based on the breadthirst and depth-irst exploring methods[END_REF], [START_REF] Roerdink | The watershed transform: deinitions, algorithms and parallelization strategies[END_REF], and since the development of the irst fast transforms, efort is put on research of parallel algorithms and architectures. This has led to the recent development of techniques used on clusters [START_REF] Galilée | Parallel asynchronous watershed algorithm-architecture[END_REF] and FPGA [START_REF] Trieu | Real-time image segmentation based on a parallel and pipelined watershed algorithm[END_REF].

The rise of GPU cards programmable for generic purposes, with low cost and high computational power, and with a trend of improvement of the hardware without severe changes to the architecture and programming model, has established a new platform of interest. These features sets the direction of advances on the watershed algorithms for improving performance, using the concepts and deinitions established on literature, exploring the architecture for the procedures that best it, and showing some new and revisited approaches [START_REF] Kaufmann | A cellular automaton for ultra-fast watershed transform on GPU[END_REF], [START_REF] Wagner | Cell reconstruction on stream computing architectures[END_REF], [START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF].

This paper works on studying this architecture and the advances made on the watershed processing and the algorithms that use it. The usage of such models inluences how the problems must be treated, algorithms designed and programs implemented. The watershed transform is implemented on two forms: one with four steps, with diferent levels of division of tasks, ranging from the neighborhood to global processing; and one that uses only a neighborhood operation, switching values until the optimal solution is found and stabilizes. These different levels of decomposition are analyzed, and the implementation strategies described.

This paper is organized as follows: Sec. 2 discusses the inluence of the GPU architecture on the development of image processing algorithms; Sec. 3 reviews the literature of parallel watershed algorithms and discusses the two algorithms evaluated in this work; Sec. 4 discusses in details the implementation of both algorithms; Sec. 5 discusses the measurements made, and exposes the performance obtained on two GPU cards and a CPU; and lastly, Sec. 6 presents the conclusions drawn from this work.

Architecture Inluence

Given many threads executing, each one having a small computing power, with some shared resources and high bandwidth for communication, the GPU architecture is well suited for parallel algorithms of ine granularity, where each thread will process one or a few pixels of the output image. Regarding the communication between threads, as they share a common memory, this area is used when any kind of information needs to be transferred between them. Also, whenever using GPU memory to communicate with the host CPU, there is a considerable overhead on the copy operation that may degrade the speedup obtained.

As many image processing algorithms, the watershed transform processes a neighborhood around every pixel. Also, as the blocks of threads divide the image on sub-images that are loaded on shared memory, these must have a border with the pixels from the adjacent blocks. This concept is illustrated on [START_REF] Podlozhnyuk | Image Convolution with CUDA[END_REF] as the apron for each block. The apron is constituted such that blocks are processed with overlapped data, with each block responsible for processing an area of the image, but loading data from adjacent blocks.

The modern GPU architecture evolved to a point where there are several forms of running programs on it, with the most noticeable frameworks and languages being CUDA and OpenCL, which are very similar to C and C++. Both models provide access to the same memory model with similar instructions. The programs developed for this work used CUDA, and the source codes are available on the Internet at http://www.adessowiki.org.

Parallel Watershed Algorithms

The watershed transform is a very data intensive task, that even with quasilinear algorithms is time-consuming. Since its introduction by Vincent and Soille [START_REF] Vincent | Watersheds in digital spaces: An eicient algorithm based on immersion simulations[END_REF] efort is put on researching faster algorithms as well as parallel approaches. An extensive survey on parallel approaches is given by Roerdink and Meijster [START_REF] Roerdink | The watershed transform: deinitions, algorithms and parallelization strategies[END_REF]. Since then, faster algorithms have been created, as well as new parallel approaches [START_REF] Trieu | Real-time image segmentation based on a parallel and pipelined watershed algorithm[END_REF], [START_REF] Galilée | Parallel asynchronous watershed algorithm-architecture[END_REF]. On the past few years, other parallel watershed algorithms have been introduced, speciically for GPU architectures [START_REF] Kaufmann | A cellular automaton for ultra-fast watershed transform on GPU[END_REF], [START_REF] Wagner | Cell reconstruction on stream computing architectures[END_REF], [START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF], using diferent strategies, designed for the restrictions of the platform. This paper analyzes two algorithms: one inspired by the drop of water paradigm and depth-search approaches, based on [START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF], named DW; and one based on cellular automata to process a shortest-path forest with sum cost function [START_REF] Kaufmann | A cellular automaton for ultra-fast watershed transform on GPU[END_REF], named CA. Both algorithms are tuned to better suit the architecture of modern GPUs.

The DW algorithm uses diferent steps of processing, which allows faster processing of operations with diferent levels of data independence. These steps are briely: (i) identiication of the neighbor with the lowest gray level for each pixel (this corresponds to inding the downstream for each pixel, see [START_REF] Meyer | Topographic distance and watershed lines[END_REF], [START_REF] Roerdink | The watershed transform: deinitions, algorithms and parallelization strategies[END_REF]); (ii) propagation of the downstream of the plateau border pixels to inner pixels according to the geodesic distance; (iii) labeling of regional minima by union-ind; and (iv) labeling of non-minima pixels by the path created to every minima. These steps are organized and designed to maximize performance on a GPU architecture.

Step (i) processes each pixel independently, scanning the neighborhood for the lowest gray level; step (ii) is the step less suited for parallelism, caused by the need to uniformly calculate the geodesic distances, creating synchronization barriers on the threads; step (iii) consists of labeling the regional minima pixels using the union-ind strategy for merging connected components; step (iv) consists on compressing paths from each pixel to the regional minima associated with the path.

The CA algorithm is based on the algorithm of Ford-Bellmann to calculate a shortest-path forest, using a single relaxing procedure, performed until convergence of the solution. Therefore, the proposed implementation consists of a single step, executed until stabilization of the solution. Nevertheless, for the CA algorithm produce a correct watersheds transform, its cost function must be adapted to consider the weight of each edge as the topographic distance [START_REF] Meyer | Topographic distance and watershed lines[END_REF]. Also, as this function does not manage plateaus, either the image must be preprocessed by lower completion, or the cost function must be further modiied to consider the lexicographic cost as a secondary component [START_REF] Körbes | On watershed transform: Plateau treatment and inluence of the diferent deinitions in real applications[END_REF]. For greater adherence with the original proposal, the lower completion is used prior to watershed execution.

Implementation Details

This section discuss the implementation of watershed algorithms on the SPMD and SIMD models, considering the access to the memory layers of the GPU, that constitutes the problem of communication between blocks. It is also presented the methodology used for development, highlighting the border treatment, the replication of data, and the data low between memories.

The GPU memory model is hierarchically divided, and programs must consider these diferent levels and its characteristics to take advantage of its features, ultimately to reach the best performance. Using the CUDA framework, the code that is executed on the GPU does not have access to the CPU RAM memory, being an area of access exclusive of the host code. Therefore, for any processing on the GPU, the data must be transferred between host (CPU) memory and device (GPU) memory. The GPU global memory is the area writeable by the CPU, and it also provides special access modes of texture and constant, which are cached and only readable by the device code. There are three other memory areas available to the device code: the registers, the local and the shared memory. These memories have very limited size, but are located close to each stream multiprocessor of the GPU, and provide high speed of access. Registers and local memory are bound to the scope of each thread, whereas the shared memory, with scope on the block, is used for cooperation and communication between threads, with also a high speed access.

Based on the communication relations between diferent types of memory, the host and device code that composes the watershed processing steps are modeled to use as most as possible the shared memory area. Initially the image is loaded by the host on a texture memory, and the neighborhood relation data on the constant memory. As the blocks process parts of the image, it is loaded to the shared memory. From this, the results are computed on registers, and inally written on the global memory. This low is used on all kernels of the algorithms developed, except those that do not use a neighborhood relation, where the global memory is used directly.

With the decomposition of the image on sub-images arises the problem of management of the borders of each block, as each pixel demands data from the neighbors that may be contained on another block, to complete its computation. The strategy used is to load an overlapping area, called the apron, as described on [START_REF] Vitor | Rastreamento de alvo móvel em mono-visão aplicado no sistema de navegação autônoma utilizando GPU[END_REF], [START_REF] Podlozhnyuk | Image Convolution with CUDA[END_REF], according to the neighborhood relation used. The borders values are either the values of pixels of the adjacent block, or in the case of image borders, a constant and predeined mask value is used.

The overlapping scheme, with borders loaded for every block by the device code, is presented on Fig. 1. This scheme is divided on two phases: the irst one loads the data to the shared memory, while the second processes the shared memory and stores the results back on global memory. The loading is performed on blocks of BLOCK_T ILE width and height, which corresponds to the size of blocks of threads used on the CUDA device code invoking. The processing is performed on fewer threads, on a block of REAL_T ILE width and height, which corresponds to the dimensions of the block of data that is actually processed, without overlapping, ignoring border pixels.

Algorithm DW

The DW algorithm is based on the four steps described on Sec. 3 and extensively detailed on [START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF]. To implement this algorithm, six programs of device code (kernels) were created:

-donwstreamCalc: calculates the downstream for each pixel and then uses the data already loaded on memory to propagate it to plateaus internal pixels, until stability inside the block. This kernel is executed once and is associated to the irst and second steps of the algorithm.

-plateauPropagate: continues the propagation of the downstream to plateaus internal pixels inside the block, respecting the geodesic distance to the borders. As the distance propagation may require block communication, this kernel is invoked until stabilization and is associated with the second step. -mergeRegions: creates connected components with the minima pixels, using the index as the reference for the group. The regional minima may spread across diferent blocks, so this kernel is invoked until stability, paired with the pathCompression kernel. This kernel is associated with the third step of the algorithm.

-pathCompression: updates the reference of each pixel directly to the connected component root. This avoids unnecessary memory accesses when labeling the image and merges connected components that spread across blocks. This kernel is associated with the third and fourth steps of the algorithm. On the third step this is invoked until stabilization of the solution, and only once for the fourth step.

-labelPixel: updates the label of each pixel, according to the reference of its group. For regional minima processing, this connects minima that spread on more than a single block, unifying the labels. This kernel is associated with the third step of the algorithm.

-indexAdjust: updates the references matrix to allow the pathCompression kernel to work on all pixels, rather than only on the regional minima, as a preparation for the fourth step.

The minima and paths labeling problem is addressed in many ways on both sequential and parallel algorithms, being a usual bottleneck, where common strategies have not obtained success on the GPU [START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF]. For that matter, these steps on the current approach are performed using the labeling algorithm of [START_REF] Hawick | Parallel graph component labelling with GPUs and CUDA[END_REF] , based on a reference list for path compressing and representative propagation. This algorithm is implemented through the kernels mergeRegions, pathCompression and labelPixel presented before.

To better understand the bottlenecks of this algorithm, each kernel and API call is timed, and a proile of the amount of time spent on each step is produced, shown on Fig. 2. The irst kernel executed, downstreamCalc, takes around 60% of the total time of execution. This is expected, as this step calculates both the downstream and solves small plateaus that are contained inside a block. The next kernel invokes, for plateau propagation, take around 16% of the total time, with 8 calls, each with a smaller computation time, until stabilization. The regional minima and catchment basin labeling are solved with a small percentage of the total time, indicating an algorithm highly adapted for the task in parallel. The irst kernel may look as a severe bottleneck, compared to the second. However, it was observed that this kernel solves many small plateaus, reducing the amount of work to be done on the plateau propagation kernel. Clearly the division of the algorithm on steps described on Sec. 3 does not match exactly the kernels described on this section, as these are designed to minimize the amount of time consumed on memory copies and maximize the processing on each pixel on a single call. However, those steps divide the problem logically, leading to other solutions, such as the sequential algorithm that is used on the comparisons of the next section. This implementation is also based on the drop of water principle and is divided on the aforementioned steps, and is not parallelized on CPU threads, thus running on a single core. The solutions provided by both algorithms are equivalent and valid.

Algorithm CA

The CA is based on the algorithm described on Sec. 3 and detailed on the work of Kaufmann and Piche [START_REF] Kaufmann | A cellular automaton for ultra-fast watershed transform on GPU[END_REF]. This algorithm, inspired on the algorithm of Ford-Bellmann for calculating shortest-path forests with the sum path cost function, iterates performing changes on the solution whenever a path of less cost is found, until stabilization. To apply this concept on the watershed transform, the weight of each edge of the inner graph must be associated with the topographical distance, as described by Meyer [START_REF] Meyer | Topographic distance and watershed lines[END_REF]. To address the issue of the cost on each edge depending on the minimal value of the neighborhood, the implementation is based on two kernels:

initialization: initializes the cost for each pixel, depending on whether they are a regional minima or not. Also, scans the neighborhood to store the minimal value, necessary to compute the cost on each edge, during the processing stage. iterate: process each pixel looking for paths of less cost. This kernel process each block until stabilization of the sub-image solution, and also is invoked as many times as necessary for stabilization of the whole image solution.

To use this algorithm, it is necessary to provide markers, either deined by the user or the regional minima. However, as the kernel that propagates the solution is invoked as many times as necessary, the greatest path on the inner graph will bind the execution time. Thus, the performance is proportional to the ratio of the area of catchment basins by the area of regional minima. To implement this algorithm eiciently, the greatest improvement is the preprocessing of the cost, which leads to several less memory accesses.

Performance Measurements

This section presents the experiments performed with the watershed algorithms discussed on this paper, divided on three subsections: irst, the algorithms are measured and compared one against the other on the GPU; second, algorithm DW is compared against a sequential implementation also based on the drop of water principle; third the algorithm DW is timed on two diferent GPU cards, to evaluate the impact of the evolution of hardware.

The hardware used is: an NVIDIA GTX 295 GPU card with 240 cores running CUDA 2.30; an NVIDA GTX 470 GPU card with 448 cores running CUDA 3.20; and a CPU AMD Phenom II X3 CPU of 2.6Ghz clock and 7.5Mb of cache with 4 GB of RAM. The measurements of GPU algorithm execution also consider memory transfer from CPU to GPU. For the CA algorithm, the computation of regional minima is discarded. This computation is not measured as it is out of scope of this work and would implicate on an additional processing that may not be useful if it is supplied markers obtained from other methods, such as user input. The algorithms were run on the images lena, cameraman, peppers and baboon. These images have diferent proiles, with varying number of regional minima and extension of plateaus. The measurements were averaged, and the standard deviation presented indicates the variation on each image size. The images were resized to 64x64, 128x128, 256x256, 512x512, 1024x1024 and 2048x2048.

GPU comparison

This experiment intends to compare the algorithms DW and CA executing on the GTX 295 card. Tab. 1 shows the average times obtained for each image size, and the standard deviation (STD) for both algorithms. Clearly, the algorithm DW produces results much faster and with smaller variation of time. The large standard deviation observed for CA is a direct consequence of its single-step propagation that will take as long to execute as the longer path on the image. The data of Tab. 1 is also shown as a chart, on Fig. 3, where the diference of times measured is better visualized. This chart presents the measured times in milliseconds (ms) by the total number of pixels on the image. The results of this experiment may also be compared with the results of the algorithm of Wagner and Godehardt [START_REF] Wagner | Cell reconstruction on stream computing architectures[END_REF], which was executed on a GTX280, a card very similar to the GTX295. This comparison is possible only for the 1024x1024 image size, which has nearly the same number of pixels of the cube of dimension 100 reported by the authors, with an execution time of 550ms.

GPU and CPU comparison

This experiment compares the algorithm DW running on the GTX 470 card against the CPU version. The measurements obtained on the executions are shown on Tab. 2. The speedup calculated is the ratio of the CPU by the GPU time. It is seen that for images smaller than 256x256, the speedup is less than two, and given that the absolute times measured are very small, the usage of a GPU algorithm may not prove useful. However, for images larger than that, the speedup is increased up to 6.5, and the absolute times make this acceleration even more useful, which can be seen on the chart of Fig.

GPU evolution comparison

This experiment intends to show how the evolution of GPU cards impact the execution of the same code. On Tab. 3 it is shown the execution times of the algorithm DW on both cards. The speedup is calculated as the ratio of the measurement on the GTX 295 by the timing of GTX 470. The average of the speedup of the images is 1.780. This might be associated with the increase of the number of cores from 240 to 448, which gives a ratio of 1.86. However, it must be noted that there are several other improvements on the architecture of GTX 470, such as enhancing cache capabilities and the operation of atomic functions. This experiment focuses on how the algorithms DW and CA perform on the same images, iltered to have less regional minima. To evaluate this, the four images were used on the size of 512x512 and iltered to have diferent number of regional minima, ranging from 562 to 24848. However, as the images have diferent proiles, only the image Lena is used to show the efect of iltering, as its results were the most afected. The times measured are presented on the chart of Fig. 5. This chart shows the measured times by the number of regional minima. This chart conirms the irst experiment, that algorithm DW is faster that algorithm CA, and shows that the algorithm DW is more stable, with most times bound between 10 and 20 ms, with a standard deviation of 4.4 ms, while algorithm CA is bound between 20 and 50 ms, with a standard deviation of 9.3 ms. The efect of severe variation of running time observed on algorithm CA is not observed for every image tested. Nevertheless, such variation was not observed for algorithm DW for any of the images.

Conclusion

This paper presented an analysis of modern watershed algorithms designed for GPU architectures with some considerations on the implementation and design of such algorithms. The issues of implementing a parallel watershed algorithm were discussed, and two implementations using the CUDA framework were detailed.

On the performance measurements, three scenarios were used, comparing algorithms on the GPU, the GPU against the CPU and GPU cards of diferent conigurations. The comparison on the GPU showed that having steps specialized for each task, and less dependence on stabilization produces faster results. The comparison with the CPU showed that studies of GPU algorithms -i.e. algorithms that use the SPMD model -may lead to reasonable speedups. The comparison of evolution of cards showed that using the number of cores of a board is a reasonable measure of normalization of times, and that the evolution of hardware may further reduce the execution times and restrictions that currently exist. Also, it has been observed that, because of the levels of speedup achieved, especially when considering the evolution of the GPUs, the cost of implementation of these algorithms is rewarded. In fact, the usage of CUDA and/or OpenCL technologies, with comprehension of the general architecture, enables the development of programs that are not strictly bound to the hardware, and may run on several cards with the beneits of speedup. As a consequence, future works are seen on the investigation of more algorithms focused on the SPMD model, with implementations not dependent on special hardware coniguration and also suitable for 3D volumes.

Fig. 1 .

 1 Fig. 1. Overlapping of data loaded on shared memory to process border pixels correctly.It is used extra threads inside each block to load data, which are discarded and do not perform the computation.

Fig. 2 .

 2 Fig. 2. Percentage of the total execution time taken by each kernel and API call of the DW algorithm.

Fig. 3 .

 3 Fig. 3. Chart comparison between algorithms DW and CA running on GTX295.

Fig. 4 .

 4 Fig. 4. Chart comparison of times between GPU and CPU implementations.

Fig. 5 .

 5 Fig. 5. Chart comparison of measured times by the number of regional minima on image Lena of size 512x512.

Table 1 .

 1 Comparison of times of algorithm DW and CA running on GTX295

	Image size Alg. DW (ms) STD Alg. CA (ms) STD
	64x64	1.065	0.064	1.099	0.094
	128x128	1.224	0.126	1.871	0.320
	256x256	2.431	0.624	4.872	0.946
	512x512	5.801	1.035	19.290	5.224
	1024x1024	19.918	3.866	109.532	33.655
	2048x2048	88.262	23.322	657.115	253.307

Table 2 .

 2 4. Comparison of times between GPU GTX470 and CPU implementations

	Image size GPU (ms) STD CPU (ms) STD Speedup
	64x64	0.646	0.043	0.311	0.003	0.481
	128x128	0.776	0.068	1.237	0.059	1.594
	256x256	1.213	0.156	4.838	0.337	3.988
	512x512	3.365	0.285	18.744	1.414	5.571
	1024x1024	10.889	0.873	72.087	3.510	6.620
	2048x2048	46.487	8.200	307.067	14.208	6.605

Table 3 .

 3 Comparison of times between two diferent GPU cards

		Image size GTX 295 (ms) STD GTX 470 (ms) STD Speedup
		64x64	1.065	0.064	0.646	0.043 1.649
		128x128	1.224	0.126	0.776	0.068 1.578
		256x256	2.431	0.624	1.213	0.156 2.004
		512x512	5.801	1.035	3.365	0.285 1.724
		1024x1024	19.918	3.866	10.889	0.873 1.829
		2048x2048	88.262	23.322	46.487	8.200 1.899
	5.4	Algorithm scalability comparison


Giovani Bernardes Vitor
email: giovani@fem.unicamp.br

André Körbes
email: korbes@dca.fee.unicamp.br

Roberto De Alencar Lotufo

Janito Vaqueiro Ferreira
email: janito@fem.unicamp.br

Analysis of a step-based watershed algorithm using CUDA

This paper proposes and develops a parallel algorithm for the watershed transform, with application on graphics hardware. The existing proposals are discussed and its aspects briefly analysed. The algorithm is proposed as a procedure of four steps, where each step performs a task using different approaches, inspired by the existing techniques that best adhere to the problem addressed. The algorithm is implemented using the CUDA libraries and its performance is measured on the GPU and compared to a sequential algorithm running on the CPU, achieving an average speedup of two times the execution time of the sequential approach. This work improves on previous results of hybrid approaches and parallel algorithms with many steps of synchronisation and iterations between CPU and GPU.

Introduction

The identification of objects on images needs in most cases a pre-processing step, with algorithms based on segmentation by discontinuity or the opposite, by similarity. Segmentation itself is not a trivial task, being among the hardest ones in image processing [START_REF] Gonzalez | Digital Image Processing[END_REF]. Some of its inherent problems are illumination variation through image sequences, dynamic change of background where the camera and/or target move, as well as changes on the topology of the target or its partial or full occlusion, resulting on a higher complexity of the scene and therefore a more difficult problem.

The watershed transform is a useful tool that is a part of some segmentation processes, for tasks of region labelling. Since its introduction by [START_REF] Digabel | Iterative algorithms[END_REF] constant efforts have been made to improve its performance, in sequential architectures such as the first fast algorithm by [START_REF] Vincent | Watersheds in digital spaces: An efficient algorithm based on immersion simulations[END_REF] until the recent ones by [START_REF] Osma-Ruiz | An improved watershed algorithm based on efficient computation of shortest paths[END_REF] and Cousty et al. (2009), as well as in parallel, such as the framework developed by [START_REF] Moga | Implementation of a distributed watershed algorithm[END_REF] until the algorithm of [START_REF] Galilée | Parallel asynchronous watershed algorithm-architecture[END_REF]. These efforts, combined with the recent dissemination of out of the box SIMD (Single Instruction Multiple Data) parallel architectures through GPUs (Graphics Process Unit) and the languages and libraries that ease its development, suggest a new approach on research for accelerating the watershed implementations.

This paper aims for the development of a watershed transform algorithm suited for GPU processing, developed with NVIDIA's CUDA (Compute Unified Device Architecture) platform, in four major steps. A fully parallel SIMD algorithm is proposed, taking advantage of the massive parallelism using specialised algorithms for each necessary step. The proposed algorithm is tested on a set of images and compared to the faster sequential algorithm that the authors built, based on the works by [START_REF] Lin | Comparison between immersion-based and toboggan-based watershed image segmentation[END_REF] and [START_REF] Osma-Ruiz | An improved watershed algorithm based on efficient computation of shortest paths[END_REF], that is very similar in function to the algorithm proposed.

The paper is organised as follows: Sec. 2 revises the state of the art on the watershed transform algorithms both on sequential and parallel fields. Sec. 3 goes through the implementation and architecture of GPGPU (General-Purpose computation on Graphics Processing Units) processing with the CUDA libraries. Sec. 4 exposes the algorithm developed, explaining its operation. Sec. 5 presents the timing results obtained, comparing the sequential and parallel algorithms. Lastly, Sec. 6 depicts conclusions obtained with this work.

Watershed Transform

Taking different approaches, the watershed transform is defined diversely on the literature, each of which producing a different set of solutions [START_REF] Vincent | Watersheds in digital spaces: An efficient algorithm based on immersion simulations[END_REF][START_REF] Meyer | Topographic distance and watershed lines[END_REF][START_REF] Lotufo | The ordered queue and the optimality of the watershed approaches[END_REF][START_REF] Bieniek | A connected component approach to the watershed segmentation[END_REF]Cousty et al., 2009). The definitions are based on regional or global elements, such as influence zones and shortest-path forests with max or sum of weights of edges, or on local elements, such as the steepest descent paths, where the neighbour's information is used to create a path to the corresponding minimum. These definitions are thoroughly revised on the literature [START_REF] Audigier | Watershed by image foresting transform, tie-zone, and theoretical relationships with other watershed definitions[END_REF][START_REF] Roerdink | The watershed transform: definitions, algorithms and parallelization strategies[END_REF]Cousty et al., 2009). In this work the local condition definition is used, henceforth named LC-WT (Local Condition Watershed Transform).

The LC-WT definition was introduced by [START_REF] Bieniek | A parallel watershed algorithm[END_REF] with the purpose of achieving speedups on the parallel watershed transform, once it mimics the behaviour of a drop of water on a surface, requiring less steps of global processing, and thus requiring less communication and synchronisation. Next, we discuss the algorithms that implement this definition on their sequential and parallel versions.

Sequential Watershed Algorithms

The recent faster sequential watershed transform algorithms are the result of the evolution of the arrowing technique for watershed of [START_REF] Bieniek | An efficient watershed algorithm based on connected components[END_REF] and the union-find of [START_REF] Meijster | A disjoint set algorithm for the watershed transform[END_REF]. The arrowing is the algorithmic representation of the drop of water, where for every pixel of an image, an arrow is drawn from the current to the next pixel creating a path that ultimately leads to a regional minimum. The arrow indicates the direction that a drop of water would slide, considering the image as a surface. For the LC-WT definition, every pixel that does not belong to regional minima will have one and only one arrow. The union-find technique for watershed transform is based on the algorithm for disjoint sets of [START_REF] Tarjan | Data structures and network algorithms[END_REF]. Given that the regions of the output image form a partition and those are disjoint by definition, the union-find algorithm processes paths to identify the roots -or representativesfor every pixel.

Several algorithms based on this preliminary works have been proposed, using variations of the previous procedures, achieving considerable speedups without loss of precision [START_REF] Sun | A fast watershed algorithm based on chain code and its application in image segmentation[END_REF][START_REF] Lin | Comparison between immersion-based and toboggan-based watershed image segmentation[END_REF][START_REF] Osma-Ruiz | An improved watershed algorithm based on efficient computation of shortest paths[END_REF]Cousty et al., 2009). These algorithms are all based on evaluating the neighbourhood to identify the direction of sliding of a drop of water on a surface until it reaches a minimum, and label the regions where the drops falls into the same minimum. However, these procedures are performed using different types of data structures and approaches to the intrinsic model, such as processing the watershed on a graph built upon the image but where pixel values are addressed to edges and cuts are made on these in order to separate regions (Cousty et al., 2009).

The work by Cousty et al. (2009) introduced one of the fastest sequential watershed algorithms due to its linear complexity. However, its constraints and style of switching between breadth-first and depth-first propagation compose a hard problem for a parallel version. [START_REF] Osma-Ruiz | An improved watershed algorithm based on efficient computation of shortest paths[END_REF] proposed an improved algorithm on the sense of pixel visitation by optimising on previous works. Nevertheless, these improvements required a massive use of queues for synchronisation of pixel visitation. This feature, interesting for a sequential algorithm, in order to minimise memory access, demands several strategies of flow control on a parallel implementation, translated into synchronisation steps. The algorithms of [START_REF] Sun | A fast watershed algorithm based on chain code and its application in image segmentation[END_REF], [START_REF] Lin | Comparison between immersion-based and toboggan-based watershed image segmentation[END_REF], [START_REF] Bieniek | An efficient watershed algorithm based on connected components[END_REF] and [START_REF] Meijster | A disjoint set algorithm for the watershed transform[END_REF] are very similar in the sense of problem approach. Their difference is that Lin's and Meijster and Roerdink's procedure solution is unique, where Sun, Yang and Ren's and Bieniek and Moga's solution depend on scan order. These algorithms are divided in clear steps, where an analysis to identify independence of data to process simultaneously works better. Also, as these steps might be seen as a pipeline, different algorithms might be applied when the used approach does not fit on a parallel architecture.

Even though different, these implementations fall on common problems, such as regional minima detection, plateau resolution, paths construction and pixel labelling. These problems are the same for both sequential and parallel watershed algorithms, once they are parts of the main problem to be solved, i.e. the application of the definition. There are several ways to address these problems, as noted before on (Körbes & Lotufo, 2009), with inner consequences to the other steps of the algorithm. The plateau problem in particular is discussed in depth by [START_REF] Roerdink | The watershed transform: definitions, algorithms and parallelization strategies[END_REF], with a proposal of pre-processing the image. However, the algorithm applied for that matter works with the same principles as the others referenced on this section, thus creating an unnecessary overhead on the whole process.

Parallel Watershed Algorithms

The watershed transform is a useful tool for segmenting images on computer vision processes, and since its introduction studies on parallelism are performed. Initially, the focus was on typical image processing systems, where the segmentation represented a time consuming operation [START_REF] Meijster | A Proposal for the Implementation of a Parallel Watershed Algorithm -CAIP[END_REF][START_REF] Bieniek | A parallel watershed algorithm[END_REF]. [START_REF] Roerdink | The watershed transform: definitions, algorithms and parallelization strategies[END_REF] extensively surveyed the literature on this subject. More recently, the evolution of sequential algorithms along with hardware minimised the performance issue. However, for fast applications (e.g. surveillance and navigation) the sequential algorithms are not fast enough, and the focus of works on parallel watershed algorithms have changed to this area using specialised architecture such as FPGAs and clusters [START_REF] Trieu | Real-time image segmentation based on a parallel and pipelined watershed algorithm[END_REF][START_REF] Galilée | Parallel asynchronous watershed algorithm-architecture[END_REF].

The initial work of [START_REF] Meijster | A Proposal for the Implementation of a Parallel Watershed Algorithm -CAIP[END_REF] is based on the work of [START_REF] Vincent | Watersheds in digital spaces: An efficient algorithm based on immersion simulations[END_REF], in the sense of definition and problem approach. The algorithm proposed relies on a graph transformation, performing the watershed on three steps: convert the image into a graph, where each vertex is a connected component of same grey level -a plateau; perform the watershed on the graph; convert the graph into the output image. The most important task, the watershed itself, is done by performing a breadth-first search from the minima, on an iterative flooding. The speedup is a consequence of the reduction of nodes on the graph and of the independence of those. [START_REF] Bieniek | A parallel watershed algorithm[END_REF] proposed another parallel algorithm with a modified definition, which is the LC-WT. In fact, this proposal settles a framework for parallelisation, as it approaches the problem on how it should divide the image, generate unique labels for each regional minimum and merge regions, depending on a sequential algorithm executed on each region. The steps of the algorithm are: split the image on regions; for each region find the regional minima and generate labels; set temporary labels to pixels of the borders of regions; run a sequential watershed algorithm on each region; merge the regions processing the pixels with temporary labels. [START_REF] Galilée | Parallel asynchronous watershed algorithm-architecture[END_REF] introduces a new algorithm for parallel watershed transform, stating to be the first that does not require minima detection as a first step prior to definition of catchment basins. However, a first sequential step for attributing a distinct label for each pixel is necessary, which is in fact its address in raster scan. The algorithm itself is defined as a single procedure with state management and message passing, for status updating of the pixels. This way, pixels that cannot be processed only with local information wait for neighbourhood data messages to arrive in order to decide which label is taken.

Following the line of work on fast applications, [START_REF] Trieu | Real-time image segmentation based on a parallel and pipelined watershed algorithm[END_REF] developed an algorithm for FPGA architectures based on the sequential algorithm of [START_REF] Sun | A fast watershed algorithm based on chain code and its application in image segmentation[END_REF]. The use of queues and stacks for synchronisation on this algorithm is substituted by a process of stabilisation of geodesic distance values and labelling; reading the memory always sequentially, on raster or anti-raster scan order. The plateau resolution is done by explicitly calculating the geodesic distances. Minima are labelled independently on a pixel basis, with merging of these sub-regions deciding for the minor label. Previous works on parallel watershed algorithms for GPU architectures exists [START_REF] Kauffmann | A cellular automaton for ultra-fast watershed transform on gpu[END_REF]Körbes et al., 2009;[START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF] and emphasise the problem that stabilisation steps introduce when not appropriately designed. A hybrid approach also is considered, using the fastest processing unit for each of the steps of the algorithm. However, entirely parallel proposals are possible to improve the slower steps using adequate algorithms, especially when considering the labelling of minima and regions.

Another approach to the problem is to take it as a shortest path graph problem. For that matter, there are several sequential algorithms that are suited for parallelisation. [START_REF] Kauffmann | A cellular automaton for ultra-fast watershed transform on gpu[END_REF] developed a modification of the Ford-Bellman algorithm for calculation of such paths and propagation of labels. This algorithm works as a single kernel, executed on each pixel of the image, thus massively parallel, until stabilisation, with synchronisation after every iteration. On the next section the GPU processing is discussed on the context of its implications on parallel programming and on the algorithm proposed on this paper.

GPU Processing

The performance gain in the capacity of GPU devices in the last years is no longer restricted to graphics processing, as it becomes an excellent alternative for general purposes, known as GPGPU, where the CPU's known speed limits are broken with the insertion of this many-core programming paradigm. Also, the range of applications of parallel nature, the gain in speed compared to conventional computers, the simplicity and practicality in use and acquisition of this technology has drawn great attention from the scientific community.

As an example, the GeForce 8800, a card of the NVIDIA's G8 series of graphics cards has 128 units of calculation (called multiprocessors), distributed on 8 vector processors. It is an architecture similar to that found in clusters of CPUs, but confined to a single hardware device, with a style of programming called SIMD (Single Instruction Multiple Data), or, data level parallelism. For this card, execution times of up to 100 times faster than the CPU time in classical programs such as a multiplication of matrices have been obtained [START_REF] Cooperman | Gpu programmingsyllabus[END_REF]. A study by NVIDIA presented a chart comparing the computational power of an Intel CPU, measured in peak GFlops, to the NVIDIA graphics cards, where the most modern GPU architecture delivers performance up to 6 times higher than the CPUs (NVIDIA, 2009). However, this chart must be evaluated carefully, as it is problem dependent, once the parallelism level is dependent on how data are managed. Also, for some problems there might not be any speedup, if its nature is synchronous or the algorithm is not suited for the architecture. That is, before implementation on GPU, the algorithm must be analysed and rearranged to suit on a data parallel hardware [START_REF] Kirk | Programming Massively Parallel Processors: A Hands-on Approach, chapter Parallel Programming and Computational Thinking[END_REF].

To explore the potential of the GPU, a different paradigm of programming must be used, called programming flow. Data are packaged in streams and the arithmetic calculations are kernels operating on them. Excerpts of programming that have enormous arithmetic rates can be shared in order to use the most of the GPU. [START_REF] Baggio | Gpu based image segmentation livewire algorithm implementation[END_REF] characterises the structure of algorithms as follows: (1) the parallel sections of the program are identified and implemented with a kernel, which is a share of the GPU to process arrays of data in parallel on different processors; (2) the organisation of these data should follow a hierarchy for the best arrangement of processing cores.

A note to consider is that the data within the block must perform the same process, as SIMD architecture. This simple and necessary routine for GPU programming does not always fit into some problems because they cannot be arranged in parallel, thus these may not benefit of the acceleration of GPUs. On the other hand when working together with CPU and GPU, the gains are significant. On this line, NVIDIA developed an architecture called CUDA that enables the development involving sequential and parallel programming, where some sections are sequential and others parallel, depending on the problem.

CUDA

The CUDA architecture is a language binding to the C/C++ language for general purpose parallel implementation. CUDA consists of a run-time library extended from C. Its main abstraction is based on the hierarchy of thread groups, memory sharing and synchronisation. Key elements of CUDA are: common C/C++ source code with different compiler forks for CPUs and GPUs; function libraries that simplify programming; and a hardware-abstraction mechanism that hides the details of the GPU architecture from programmers [START_REF] Halfhill | Parallel Processing with CUDA: Nvidia's high-performace computing platform uses massive multithreading[END_REF]. As a complex topic out of the scope of this paper, the reader is referred to the work of [START_REF] Nickolls | Scalable parallel programming with cuda[END_REF] for a detailed view on CUDA programming and modelling.

Algorithm Proposal

In this section we propose an algorithm for fully parallel SIMD implementation. Firstly, the motivations and inspirations are presented. Next, the notation used is explained, and the algorithm depicted and explained. Lastly, implementation details are pointed out as well as performance considerations.

As presented on Sec. 2, there are several approaches for watershed algorithms. Even though the evolution of techniques achieved a great speedup in comparison with the first fast transforms, those are still not enough for applications such as automatic navigation. Thus, the parallel approach seems obvious to achieve an even greater speedup. Along with the recent development on GPU parallel processing, presented on Sec. 3, with further optimisation for image processing, a new field for the watershed transform is seen, taking advantage of this massive many-core architecture. Also, experiments with hybrid approaches already shown a speedup [START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF].

This algorithm is inspired in both sequential and parallel previous algorithms presented on Sec. 2. Among the sequential ones, the most influential is by [START_REF] Bieniek | An efficient watershed algorithm based on connected components[END_REF], which inspired on storing the pixel addresses for finding the steepest descent paths on a single step. The strategy of [START_REF] Galilée | Parallel asynchronous watershed algorithm-architecture[END_REF] of waiting for neighbourhood data to split non-minima plateaus is used. On this case, since there is no message system to exchange data between concurrent threads on the GPU, this communications is handled through synchronisation steps. Iterating until the plateaus are resolved, this algorithm requires memory copies to and from the GPU until it is stable. Other approaches are available for solving the plateau problem, however they all tend to use some kind of ordered propagation, even when transforming the image in the lower completion process. Some algorithms, such as the watershed cut of Cousty et al. (2009) choose to not split plateaus evenly, propagating over these almost randomly.

Another relevant step of any watershed algorithm is the minima and paths labelling. This problem is itself addressed in many ways on both sequential and parallel algorithms, being a usual bottleneck, where common strategies fail on the GPU as noted before [START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF]. For that matter, these steps on the current approach are performed using the labelling algorithm of K.A. Hawick et al. (2009), based on a reference list for path compressing and representative propagation. This algorithm itself is inspired on the union-find approach for path compressing [START_REF] Tarjan | Data structures and network algorithms[END_REF], though applied in parallel.

Fig. 1 shows an example of the steps applied iteratively for connected component labelling. On (a) it is shown the greyscale image with only its regional minima values. Fig. 1 (b) shows the index positions, used for calculation of the reference list, and shown on (c) as the masked data to be processed for the regional minima. Fig. 1 (d) shows the first step of the first iteration, where every pixel looks on its neighbourhood (inside the mask) for the smallest reference, and store as its own reference. Next, those indexes are propagated, compressing the path from every pixel to the minor reference on the same connected component. As seen on (e), a single iteration does not merge the whole component, thus, a second one is necessary, processing the same steps to set the same reference for every pixel, as shown on (f). Algorithm 1 presents our proposal for a parallel watershed transform. On the course of it, the statement for <condition> in parallel do denotes a kernel execution on GPU architecture, indicating that each iteration can be processed independently, in a SIMD way. The algorithm is divided in four major steps: (1) find the lowest neighbour of each pixel (direct path of steepest descent); (2) find the nearest border of internal pixels of plateaus, propagating uniformly from the borders; (3) minima labelling by minor neighbour address and (4) pixel labelling by path compression. The input image is called I, with a domain D and the output labelled image is called lab, which is also used for working. As noted before, a similar algorithm has been proposed recently, though with modified strategies for labelling (Körbes et al., 2009), that had shown worst results than a sequential approach [START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF].

Algorithm 1: Parallel watershed transform (GPU-WS)

// First Step 1: PLATEAU = 0 2: for p in D in parallel do 3:

if exists q in N(p): I(q) < I(p) and I(q) = min{I(u), u in N(p)} then 4: lab(p) = -q 5: else 6: lab(p) = PLATEAU 7:

end if 8: end for // Second step 9: while lab is not stable do 10: Tlab = lab 11:

for p in D: lab(p) = PLATEAU in parallel do 12:

if exists q in N(p): lab(q) < 0 and I(q) = I(p) then 13:

Tlab(p) = -q 14: end if 15:

end Using this algorithm, an image is taken as a grid divided in 16x16 pixels blocks, resulting in up to 256 threads. The image loaded on the CPU is copied to the texture memory of the GPU, where it is then processed via kernel access. The kernels were developed in three ways. The first is for step 1, which is based on texture memory access with block division as mentioned above, outputting its results to global memory. For the further steps, the second kernel demands a special approach, once pixel evaluation is highly guided by its neighbours. In effect, each block must include a border, containing neighbour pixels of the block, for correct processing. Fig. 2 shows the modelling used for each image block. Considering coalescent data access, the block dimension is kept on 16x16, though with a processing range of 14x14 centred on it. At that point, each block processes 196 threads, with other 60 inactive threads, where these data are processed by adjacent blocks. Once the data is loaded from the global memory to shared memory, the results are stored on output global memory, ensuring that border results do not overlap [START_REF] Podlozhnyuk | Image Convolution with CUDA[END_REF].

The third step was developed using two kernels, where the first propagates the minimal index for the neighbours on the connected component, and the second transforms the image into the inner reference list. This list is processed for equivalence of indexes until stabilisation, that is, the minor label is propagated through the connected path obtained on the first kernel. These two kernels iterate until every minima is merged with its pixels, with these iterations controlled on the CPU. Fig. 1 (d)-(f) shows these iterations of each kernel. Along with these specifications, the algorithm is dependent on stabilisation steps, which are verified on the CPU. As the experimental results presented on Sec. 5 show, these steps contribute heavily for increasing the performance of Algorithm 1, compared to a previous work [START_REF] Vitor | Fast image segmentation by watershed transform on graphical hardware[END_REF]. The last step use one kernel to solve the equivalence list, propagating the label of the regional minima through the connected path between the pixels obtained on the first step. Fig. 3 shows each of the steps applied on a sample image. On this image, there are two regional minima, with levels 0 and 10. On the first step, every pixel locates its lower neighbour, as is indicated by the arrows. The plateau's inner pixels are then resolved on step 2, with the created paths indicated by four new arrows. On step 3, the regional minima are labelled, and finally on step 4 the labels of the minima are propagated to every other pixel.

Experimental Results

The algorithm was tested on several images on sizes ranging from 64x64 to 2048x2048, and had its performance measured on each step. The data presented following is the average of the experiments, scaled to milliseconds (ms). The sequential algorithm implemented is a modified version of [START_REF] Lin | Comparison between immersion-based and toboggan-based watershed image segmentation[END_REF], which operates also on four steps very similar in purpose to those of Algorithm 1, but without the watershed label, as the authors originally proposed, for matters of equivalence of definitions and better performance. Table 1 shows the performance for the sequential algorithm, processed on the CPU. These results were obtained on a computer with an AMD Phenom II X3 (2.6Ghz 7.5Mb Cache) with 4 GB RAM and a GeForce GTX295 with 1792Mb. Only one GPU core was used for the GPU processing. 1 shows that for the sequential algorithm, the first two steps are severe bottlenecks, greatly degrading the performance, especially on larger images. Moving to the parallel implementation on the GPU, developed according to Algorithm 1, Table 2 shows that the bottleneck on step 1 is completely mitigated, once it is processed on single kernel execution, whereas for step 2 the performance depends on image size, though showing small improvements. Also the memory copies to and from the GPU must be considered, as it is a time consuming task that does not exist for the sequential algorithm. This is easily seen as the sum of each step of the algorithm on Table 2 does not match the total, where this difference corresponds to memory allocation and copying on the GPU between global and texture memories and also between CPU and GPU memories. With the data collected on these tables, a comparison may be done, checking the effective speedup in time performance. This evaluation is seen on Table 3, where the algorithms are put aside and the speedup ratio between GPU and CPU is calculated for the average total time for each image size tested. The analysis of Table 3 must highlight some points of GPU parallel processing. One is the fact that for small images, such as 64x64 and even 128x128 pixels, the memory operations produce an overhead that eliminates or greatly minimises advantages of a parallel algorithm. Once image sizes increase, this effect is reduced and the speedup is more sensible. For average sizes of images -typically those produced by video cameras -the parallel approach shows an average speedup of 2.7, meaning a reduction of time of more than half of the sequential algorithm. These results are better than the previous obtained with a hybrid algorithm (Vitor et al., 2009). Nevertheless, the algorithm might be further improved, especially on step 2, where a different approach for the plateau resolution could increase the speedup.

Noting that a considerable speedup has been achieved, it is also necessary to compare this work with the one by [START_REF] Kauffmann | A cellular automaton for ultra-fast watershed transform on gpu[END_REF]. For that matter, the total number of pixels is used as the base for comparison. However, as the authors were interested in 3D images, with orders of tens of millions of pixels, there is data available only for the size of 512x512 and 1024x1024 pixels. There must be used some caution when comparing these data. Given that the GPUs are in constant evolution, with increasing performance, this must be taken into consideration. [START_REF] Kauffmann | A cellular automaton for ultra-fast watershed transform on gpu[END_REF] used a Radeon X1950 Pro with 512MB of memory whilst in this paper was used a GTX295 with 1792MB of memory. In order to make a reasonable comparison, and since both GPUs have the same internal clock of 595 MHz, the number of cores is used as a basis. The Radeon X1950 Pro have 48 pixel shaders (AMD, 2010), and the GTX295 have 240 CUDA cores (NVIDIA, 2010). Thus, the measure obtained on the proposed algorithm is multiplied by this ratio, simulating its operation on 48 cores. The measures of both algorithms are presented on Table 4, where GPU-WS is the current proposal, and CA-WS is the data of [START_REF] Kauffmann | A cellular automaton for ultra-fast watershed transform on gpu[END_REF].

GPU-WS (ms)

CA-WS (ms) 512x512 12.62 x 5 = 63.1 600 1024x1024 50.85 x 5 = 254.25 1177

Table 4: Performance data of algorithm 1 and the algorithm of [START_REF] Kauffmann | A cellular automaton for ultra-fast watershed transform on gpu[END_REF] Once all the algorithms compared here yield very similar results, due to it's conceptually equivalence regarding watershed transform definition, there is no need to show the produced images.

Conclusion

This paper presents a watershed transform algorithm with implementation on the GPU using the CUDA run-time libraries. With this implementation, a speedup of more than two times was obtained in comparison with a CPU algorithm. A significant part of this improvement on previous works is due to the labelling algorithm utilised, reducing greatly the iterations necessary for the task completion. However, more effort is necessary as there is room to improve on the step that splits plateaus, which still demands a stabilisation that degrades performance when the image has large areas of such type. A comparison with another watershed algorithm based on the GPU was made, and even though based on different platforms, shows a reasonable improvement. We believe that the spreading of such many-core architectures provides the environment for great improvements on the image segmentation fields, using complex and global algorithms, such as the watershed one.

As future works, we intend to continue on the work on parallel algorithms specialised to many-core architectures, such as morphological reconstruction and watershed from markers.

Figure 1 :

 1 Figure 1: Steps for connected component labelling based on reference list

Figure 2 :

 2 Figure 2: Border behaviour to include neighbours outside of the block

Figure 3 :

 3 Figure 3: Representation of the processing performed by each of the steps of the proposed algorithm

Table 1 :

 1 Sequential algorithm processed on the CPU Table

		64x64	128x128	256x256	512x512	1024x1024	2048x2048
	Step 1	0.18	0.72	3.16	11.94	45.80	171.07
	Step 2	0.27	1.10	4.36	17.84	72.69	306.14
	Step 3	0.03	0.18	0.46	1.68	4.46	17.07
	Step 4	0.05	0.15	0.93	3.56	15.58	70.06
	Total	0.55	2.17	8.92	35.04	138.54	564.36

Table 2 :

 2 Parallel algorithm processed on the GPU with time measurements in ms

		64x64	128x128	256x256	512x512	1024x1024	2048x2048
	Step 1	0.05	0.08	0.19	0.60	2.31	9.18
	Step 2	0.22	0.52	1.53	6.00	26.61	186.80
	Step 3	0.24	0.67	1.52	4.02	18.95	79.81
	Step 4	0.02	0.04	0.06	0.16	0.69	3.75
	Total	1.52	1.94	3.54	12.62	50.85	284.96

Table 3 :

 3 Comparison between CPU and GPU algorithms using the speedup as the ratio of CPU by GPU measured time.

		CPU (ms)	GPU (ms)	Speedup
	64x64	0.55	1.52	0.36
	128x128	2.17	1.94	1.11
	256x256	8.92	3.54	2.52
	512x512	35.04	12.62	2.78
	1024x1024	138.54	50.85	2.72
	2048x2048	564.36	284.96	1.98


1

COURTS

In Medina, the Prophet Muḥammad acted as an arbiter within the early community (ummah). His role was sanctioned by the Qurʾān, which prescribed that he judge according to the divine revelation (5:48-49). After his death in 632, the first caliphs took over his role and dispensed justice in Medina. The date when the first judge (qāḍī) was delegated is still the subject of controversy. Reports that the Prophet himself sent judges to Yemen are rejected by modern scholarship. According to Islamic sources, the second caliph, ʿUmar (r. 634-644), appointed qāḍīs in conquered territories soon after the beginning of the conquest. Although the status of these first judges and of their legal rulings is unclear, the first Muslim armies certainly needed men appointed by an authority to arbitrate the inevitable conflicts among the soldiers. However, because of the contradictory sources, some scholars dismiss these reports and defer the appearance of the first judiciary system until the reign of .

The Sharīʿah Court

One judge alone dispensed justice in every large Muslim city, except in a capital such as Baghdad, which from the late eighth century was divided into three districts and from the eleventh century on into four districts. In the aftermath of the conquest, only cities with Friday mosques had a qāḍī. As an increasing number of conquered peoples converted to Islam, especially from the ninth century on, jurisdictions were created in secondary towns, where the judge was usually a deputy (nāʾib) of the qāḍī of a larger city. Judgeships were held exclusively by adult Muslim males of unsullied reputation.

The main task of the qāḍī was to dispense justice between litigants. During the classical period, he sat in a public place, usually the mosque (some qāḍīs preferred holding court in their home), where a plaintiff could come and lodge a complaint against an adversary. If the defendant refused to appear, the judge could summon him and even have him forcibly brought. If the defendant did not confess his guilt, the burden was on the plaintiff to prove his claim by providing at least two honorable witnesses (a process known as bayyinah, or testimonial proof). Alternatively, Medinan doctrine would accept one witness if the plaintiff took an oath. If there were no witnesses (or if the witnesses were not reliable), the judge would ask the defendant to swear that he was innocent. Ḥanafīs regarded the defendant's refusal to take an oath as a proof of his culpability. However, according to Mālikīs and Shāfiʿīs, the defendant's decline to swear was no proof and in this situation the judge had to refer the oath to the plaintiff. Some circumstantial evidence could also be taken into consideration in favor of the defendant, on the word of an expert (Johansen, 2002, "Signs as Evidence," p. 174). On the basis of this legal evidence, the judge then issued a written and binding decision, which could only be repealed if there were a serious breach of law. The judge was assisted in his task by one or more clerks and by a body of professional witnesses-these first appeared in Egypt at the turn of the ninth century, and thereafter this custom seems to have spread to other provinces. With the help of other assistants, the qāḍī also performed "administrative" tasks, such as the supervision of endowed foundations (waqfs) and the control of the property of orphaned, disabled, or bankrupt people.

Relationships with the Rulers

From the very beginning of Islam, justice was the prerogative of the ruler. He delegated his judicial authority to qāḍīs, who were appointed for an indeterminate length of time. At any moment, the ruler could dismiss them and appoint other judges in their place. The qāḍī received a salary and was therefore considered a civil servant.

Under the Umayyads (661-750), who ruled from Syria, provincial governors were usually responsible for the appointment of their qāḍīs. Although some Umayyad caliphs, such as ʿUmar II (r. 717-720), Yazīd II (r. 720-724), and Hishām b. ʿAbd al-Malik (r. 724-743), attempted to overrule the authority of some governors by appointing provincial qāḍīs themselves, centralization of the judgeship under the authority of the caliph was only completed under the ʿAbbāsids. The second ʿAbbāsid caliph, al-Manṣūr (r. 754-775), began to appoint judges in the principal Iraqi cities, and later in the provinces. This reform was part of a larger political agenda that was meant to reinforce central power, the main symbolic expression of which was the foundation of the round city of Madīnat al-Salām (Baghdad) in 762. Provincial governors did not easily accept the undermining of their power, however, and attempts were made to regain control of the judiciary, especially when the central power was at its most vulnerable (as, e.g., during and after the civil war between the caliphs al-Amīn and his brother al-Maʾmūn, in 811-813). On the whole, ʿAbbāsid control of the qāḍīs became widely accepted except in territories where the caliph was not recognized, as in al-Andalus (Muslim Spain, where a dynasty of Umayyad governors came to power in 756). The clearest representation of the centralization of the judiciary was the creation of the office of chief judge (qāḍī al-quḍāt, literally, judge of judges) in Baghdad around 790, under the caliph Hārūn al-Rashīd (r. 786-809). Although the chief judge rarely appointed other qāḍīs himself, he was a close advisor of the caliph, and he had a decisive influence on the selection of judges. In al-Andalus, the judge of Córdoba, who had held the title of qāḍī al-jamāʿah (judge of the community) since the beginning of the Umayyad rule, tended to be regarded as the western counterpart of the qāḍī al-quḍāt.

From the late ninth century onward, when provinces began to claim political and financial autonomy from the ʿAbbāsid caliphate, asserting centralized authority over local judgeships became a major issue. With the political decline of the ʿAbbāsid caliphate in the early tenth century, the vizier (chief of the central administration) played an increasingly important role in the appointment of judges, but after 945, the Būyid amirs soon extended their authority over the judiciary. They were replaced in 1055 by a Seljuk sultan who claimed the same prerogative over the judgeship, even though the vizier Niẓām al-Mulk (d. 1092) acknowledged that the qāḍī was still theoretically a deputy of the caliph.

The selection of qāḍīs according to their schools of law (sing. madhhab) was also a political issue. The ʿAbbāsids selected their first qāḍīs from among the Hejazī jurists, perhaps because promoting Medinan law (which relied on ʿamal, the uninterrupted practice of the Medinan community since the death of the Prophet) would reinforce their legitimacy. From the reign of the caliph al-Mahdī (r. 775-785) onward, the ʿAbbāsids gravitated toward the followers of the Iraqi scholar Abū Ḥanīfah (d. 767), the eponym of the Ḥanafī school of law, whom they appointed as qāḍīs in the principal regions of the empire. This policy met resistance by urban communities who were still attached to their local legal traditions. During the ninth century, qāḍīs were increasingly selected for their theological beliefs, and the Ḥanafīs lost their hegemony after 871. Thereafter, rulers appointed judges mainly from the Ḥanafī, Mālikī, and Shāfiʿī legal schools (except for the Shiite Fāṭimids, r. 909-1171, who preferred Ismāʿīlī qāḍīs), with each dynasty favoring one or another of these schools for the chief judgeship. In 1265 the legal pluralism at play was formalized under the Mamlūks (1250-1517) when Sultan Baybars appointed four chief judges, one for each Sunnī school (Shāfiʿī, Ḥanafī, Mālikī, and Ḥanbalī), first in Cairo and then in Damascus. This reform, which from a legal point of view was quite notable, has also been interpreted as a political act undertaken to weaken the authority of a single chief judge who had been showing annoying signs of independence.

Tension and Judicial Independence

Tension between the ruler and the judiciary existed since the Umayyad period. The qāḍī was a mere deputy of the ruler and depended on him for the execution of some of his judgments. For example, judges had no authority over the prisons in which they incarcerated debtors, and local governors or chiefs of police could release prisoners without the judge's consent. Rulers could easily overrule or interfere with the judiciary, whereas the qāḍīs, who were hired from among the scholars who were the authors of the Islamic jurisprudence, considered themselves the main judicial authority. The centralization of legal appointments under the caliph al-Manṣūr allowed judges to be free of the authority of local governors, but they were beholden to the caliph, who also interfered in the judicial process.

Under the ʿAbbāsids, qāḍīs became a major political tool. The dynasty had seized power by denouncing the impiety and the injustice of the Umayyads. In order to maintain their legitimacy, at a time of intense legal development that led to the birth of the classical schools of law, the first ʿAbbāsid caliphs surrounded themselves with scholars and used qāḍīs to embellish the official facade of the state. Qāḍīs were used as advisors, witnesses, and emissaries. Caliphs also used them to conduct political trials, or to issue fatwās allowing the execution of political enemies. The political use of qāḍīs reached its peak during the miḥnah ("examination"), instituted at the end of the caliph al-Maʾmūn's reign and lasting until the beginning of that of . In an effort to restore his authority, which was being challenged by the traditionalist movement, the caliph sought to impose his control over the judicial system and through it impose his control over the whole Muslim community. He determined that only qāḍīs who adhered to the official theological dogma that the Qurʾān was a part of God's creation were to remain in office, and they had to check that all witnesses (who usually belonged to the social and religious elite) professed the same dogma.

Religious scholars in general and qāḍīs in particular argued increasingly against such manipulation of the judiciary. They disseminated prophetic reports (ḥadīths) insisting on the individual responsibility of the judge, who would be held responsible for his actions in the hereafter. They emphasized that the judge's decision relied on scriptural sources-including the consensus (ijmāʿ) of jurists-or on his own assessment (ijtihād), and that left no place for caliphal intervention. During the ninth century, Ḥanafī jurists began to question the link between the caliph and his qāḍīs. According to the jurist al-Khaṣṣāf (d. 874), the qāḍī could issue a judgment against the caliph, despite his being a deputy of the latter. A century later, al-Jaṣṣāṣ (d. 981) asserted that the caliph did not appoint a judge in his own name but that the judge acted as a deputy for the whole community. The qāḍī served the Muslims (and not the caliph) on the basis of the law formulated by the legal scholars (fuqahāʾ). According to this theory, a qāḍī could be appointed by rebels (Khārijites), and even by the community itself in the absence of a ruler. Al-Jaṣṣāṣ's theory was later adopted by other scholars, such as the Mālikī al-Bāqillānī (d. 1013). Although the caliphate recognized a greater independence of the qāḍīs after the end of the miḥnah and the victory of Sunnism under al-Mutawakkil (r. 847-861), Muslim rulers never ceased in their attempts to manipulate the judiciary.

Shiite Discourse on the Judiciary

Imāmī Shiites, who believe that the only true rulers (imams, in Shiite discourse) were the fourth caliph, ʿAlī ibn Abī Ṭālib (r. 656-661), and eleven of his descendants, regarded the Umayyads and the ʿAbbāsids as illegitimate rulers. According to al-Kulaynī (d. 940), the imam Jaʿfar al-Ṣādiq (d. 765) forbade people to appeal to official qāḍīs. Instead, Shiite followers were advised to seek the judgment of a Shiite scholar, who would be considered a deputy of the true imam. The Imāmī doctrine evolved under the Būyid dynasty (945-1055), which was of Shiite persuasion. According to al-Ṭūsī (d. 1067), an Imāmī Shiite appointed qāḍī by an illegitimate ruler could dispense justice as long as he remained faithful to his own doctrine and regarded himself as a deputy of the true imam. It was the duty of an Imāmī jurist to accept his appointment, lest an unjust (Sunnī) scholar be appointed qāḍī in his place.

Parallel Judicial Institutions

The administration of law was exercised in forums other than the Sharīʿah court, which allowed rulers to assert their judicial authority. Himself a source of judicial power, the ruler could preside over the litigations and give binding decisions, or he could appoint deputies who unlike qāḍīs were not bound by the rules of the Sharīʿah (even if these deputies could be Sharīʿah court judges at the same time). Often called "secular" justice because it took place outside the Sharīʿah court, it manifested itself in different ways.

The Maẓālim Court

The maẓālim (literally, "wrongs") court symbolized the discretionary authority vested in the ruler. It had a broad jurisdiction-it received petitions against officials and abuse of power, and could occasionally serve as court of appeal against the qāḍīs' judgments. The judge (theoretically the ruler, but more often his deputy) was usually assisted by jurists and qāḍīs, who guaranteed the legitimacy of the decision, as well as by clerks and witnesses. Procedure was not restricted by legal doctrine, however. Whereas the qāḍī was bound by the accusatory procedures formulated by the religious law and could not violate the evidentiary rules, the maẓālim judge had procedural powers that went beyond those-he could use the inquisitorial method, order investigations, not wait for a case to be brought, and rely on noncanonical types of proof. The ruler could thus base his decision on principles of equity, or follow the dictates of the state rather than Sharīʿah law as such.

Following the Sassanian model, maẓālim court sessions were probably organized on a regular basis during the early ʿAbbāsid period. Caliphs al-Mahdī and al-Hādī (r. 785-786) presided in person over a court. Later, the maẓālim were mostly supervised by the vizier, until the Shiite Būyids handed them to the Imāmī naqīb al-ashrāf (who ran the affairs of ʿAlī's descendants). From the eleventh century on, the maẓālim court was characterized by increasing bureaucratization and placed under the responsibility of different officials (sultans, governors, military officers). No specific locale was reserved for these public sessions until Nūr al-Dīn Zankī (r. 1146-1174) established a "house of justice" (dār al-ʿadl) in Damascus. A similar structure was built in Cairo under the Mamlūks.

From the early ʿAbbāsid period on, the maẓālim competed with the qāḍīs' jurisdiction, even though qāḍīs were occasionally appointed to head this court. ʿAbbāsid caliphs used the maẓālim forum to reinstate their authority when confronted with a qāḍī's excessive autonomy or noncompliance with official ideology, as during the miḥnah. Through this institution, provincial governors seeking autonomy from the caliphate could also reinforce their power. In the late ninth century, when local dynasties first appeared in the east, provincial qāḍīs were still appointed by the caliph. Governors who sought autonomy had to regain control of the judiciary, but they could not easily dismiss a qāḍī dispensing justice in the name of the caliph. Rulers such as Ibn Ṭūlūn (r. 868-884), his son Khumārawayh (r. 884-896), and later the Ikhshīds (935-969) in Egypt therefore regularly used the maẓālim court as a way of bypassing the ordinary judicial system. Moreover, during the Mamlūk period, the maẓālim tended to encroach on matters regulated by the Sharīʿah, which fell under the jurisdiction of the qāḍī, such as cases involving pious foundations.

The proliferation of centers of power after the eleventh century increased the importance of justice dispensed by rulers. Other courts, such as the yarghu of the Īl-Khān Mongols (1256-1335), which dealt specifically with disputes among Mongols as well as state affairs, were assimilated into the maẓālim when the dynasty converted to Islam.

The Police

The shurṭah (also called maʿūnah from the ninth century onward), was an elite unit of the army that served as a security and police force. It was established under the first caliphs or the early Umayyads to protect cities and villages against riots and banditry and to ensure the night watch. Its chief, the ṣāḥib al-shurṭah (also called shiḥnah after the end of the ninth century), was appointed by the caliph or the governor, and he was often the head of the ruler's personal bodyguard. In early Islamic Egypt, the ṣāḥib al-shurṭah was the second highest official of the province after the governor, and he acted as his deputy when the latter was away. The shurṭah could be called upon to maintain order at the qāḍī's court or to carry out his rulings.

The shurṭah was not only a military unit with control and security functions. The ṣāḥib al-shurṭah also presided over a court where he dispensed criminal justice. Like the maẓālim, his court was not bound by the rules of the Sharīʿah, and the ṣāḥib al-shurṭah operated much more freely than the qāḍī. He could conduct investigations and rule on the grounds of physical evidence. His discretionary powers allowed him to inflict punishments harsher than those prescribed by Islamic law. The ṣāḥib al-shurṭah probably judged criminals arrested and brought before him by his troops. He prosecuted on his own authority, without the necessity of a plaintiff.

Finally, the ṣāḥib al-sūq or muḥtasib, chief of the "market police," also heard complaints and dispensed justice. His principle task was the supervision of moral behavior in public. He intervened on his own accord and passed punitive judgment based on custom (ʿurf) against offenders.

Summary

In early Islam, the qāḍī was no more than a legal official under the ruler's supreme judicial power. Between the eighth and the tenth centuries, as Islamic law developed into a specific field governed by legal scholars, the qāḍīs were increasingly identified with a religious jurisdiction that necessarily had to escape from under the authority of the ruler. Therefore, two sets of judicial institutions developed, which came to complement and sometimes compete with each other. For a ruler who needed to govern according to the public interest, to ensure security beyond the prescriptions of the Sharīʿah, or simply to serve state interests, it was necessary to rely on institutions that were not bound by the strict prescriptions of Islamic law and could be monitored more easily. During the Mamlūk period, authors such as Ibn Taymīyah (d. 1328) and Ibn Qayyim al-Jawzīyah (d. 1350) attempted to repair the rupture between the two types of institutions by developing the concept of siyāsah sharʿīyah (governance in accordance with the sacred law). According to their theory, many practices of the extra-Sharīʿah institutions conformed to the spirit of Islamic law and could be justified on this ground; Sharīʿah legitimacy thus extended to actual states. Their theory would strongly influence the Ottoman legal system.

The question of the boundaries between the various courts is still unresolved. Most scholars today think the qāḍī had little jurisdiction in penal law, although the Sharīʿah holds that the qāḍī has authority in criminal cases. As Émile Tyan remarked (Tyan, 1938-1943, II, p. 411-2), procedural law limited the involvement of the qāḍī in criminal casesboth a claimant and a defendant were required, and the judge could neither act on his own nor conduct inquiries. Most penal law cases were probably heard by the security forces (especially the police). Only on rare occasions, such as when a crime was witnessed and an accusation could be leveled, did the judge's competence in criminal matters become a reality. As Christian Lange has argued (Lange, 2008, p. 48), the situation may be described as a network of overlapping jurisdictions, subject to variations according to era, place, and specific cases.

[

  See also Ḥisbah; Justice; and Sharīʿah.] Bibliography Bligh-Abramski, Irit. "The Judiciary (Qāḍīs) as a Governmental-Administrative Tool in Early Islam." Journal of the Economic and Social History of the Orient 35 (1992): 40-71.