Mahmoud Soua
email: mahmoud.soua@esiee.fr

Rostom Kachouri
email: rostom.kachouri@esiee.fr

Mohamed Akil
email: mohamed.akil@esiee.fr

GPU Parallel implementation of the new Hybrid Binarization based-Kmeans method (HBK)

Keywords: OCR, Scanned documents, Binarization, HBK, Kmeans, CUDA, GPU

GPU parallel implementation of the new hybrid binarization based on Kmeans method (HBK)

and converts them to a machine readable text file. The OCR quality is measured according to the character recognition rate [START_REF] Kae | Improving state-of-the-art OCR through highprecision document-specific modeling[END_REF]. Generally, the acceptance rate varies according to the business operation or recognition task. For a common A4 document, scanned at 300 dpi, error rates are often in the range of 1% to 10% of all characters in page [START_REF] Collins-Thompson | Improved String Matching Under Noisy Channel Conditions[END_REF]. Actually, both accuracy and high speed of OCR are required in the real time applications [START_REF] Eikvil | OCR-Optical Character Recognition[END_REF]. For example, in the banking check reading [START_REF] Eikvil | OCR-Optical Character Recognition[END_REF], the acceptable error rate is close to zero and the OCR should be very fast. It may have a very high throughput of up to 150.000 documents per hour [START_REF] Eikvil | OCR-Optical Character Recognition[END_REF]. Also, in the mail sorting [START_REF] Gaceb | A new mixed binarizatin method used in real time application of automatic business document and postal mail sorting[END_REF], the acceptable rate is close to zero. The recognition rate, is about 30.000 letters per hour [START_REF] Eikvil | OCR-Optical Character Recognition[END_REF]. In our work, we focus on the paperless [START_REF] Ashari | FPGA Implementation of realtime Adaptive Image Thresholding[END_REF][START_REF] Fong | Document Imaging: A step toward a paperless office[END_REF] real time application, in which, a high speed scanner Epson Perfection V100 Photo scans over one A4 300 dpi document per 23 second. This application is depicted in Figure 1. Both binarization and recognition are required in the OCR processing [START_REF] Kumar | MAPS: Midline analysis and propagation of segmentation[END_REF]. Binarization plays a very important role in the OCR toolchain. In this operation, documents are converted into bi-level representation distinguishing the foreground from the background. We consider that foreground objects represent printed text, a legend, or a drawing while the complementary objects correspond to the background. Indeed, a high accuracy in extracting untouchable and unfragmented characters is very important to perform efficient recognition.

Actually, we are interested in two points: enhancing the binarization quality to improve the recognition rate and speedup the execution time to increase the number of documents processed per second.

On the first point, we proposed recently an adaptive new method called Hybrid Binarization based on Kmeans (HBK) [START_REF] Soua | A new hibrid binarization method based on Kmeans[END_REF], basically designed for scanned documents. Firstly, the HBK method performs the Kmeans clustering algorithm on local areas seeking more binarization robustness. Then, it gathers local features in a global manner to improve the local binarization quality. Basically, the HBK method outperforms state of the art binarization methods such as SauvolaMS kx [START_REF] Lazzara | Efficient Multiscale Sauvolas Binarization[END_REF] and Niblack [START_REF] Niblack | An Introduction to Digital Image Processing[END_REF]. It scores a character recognition rate of 91% when performing on 125 scanned magazine documents [START_REF] Soua | A new hibrid binarization method based on Kmeans[END_REF].

On the second point, we consider our paperless application time constraint. Given that the binarization step should take, generally, 2% from the total OCR execution time [START_REF] Srinivasan | Performance characterization and acceleration of Optical Character Recognition on handheld[END_REF], the processing time required to binarize one document is 460 ms in our real time system. According to our experiments, the HBK method requires at least 1.9 seconds to process one A4 300 dpi document. Thus, it does not meet the real time constraint of our paperless application. Actually, the running time of Kmeans in the HBK method is time consuming [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF]. Therefore, parallelizing it is a promising approach to overcome the challenge of the huge computational requirement [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] and accelerate HBK. In this context, different implementations of Kmeans on the Graphic Processor Unit (GPU) have been proposed [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF][START_REF] Fang | Parallel Data Mining on Graphics Processors[END_REF][START_REF] Hong-Tao | KMeans on Commodity GPUs with CUDA[END_REF][START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF]18]. Those papers mainly address Kmeans as general purpose clustering algorithm. According to our knowledge, the first popular Kmeans GPU implementation is the UV Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF]. It achieves a speedup of 10x to 40x as compared to four-threaded Minebench [START_REF] Pisharath | Nu-MineBench 2.0[END_REF] running on a dual-core, hyper-threaded CPU. The GPU Miner [START_REF] Fang | Parallel Data Mining on Graphics Processors[END_REF] improves 5x UV Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF] by decreasing the transfer memory latency. HP Kmeans [START_REF] Wu | Clustering billions of data points using GPUs[END_REF] claims another speedup of 2x to 4x as compared with UV Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF] and 20x to 70x speedup as compared with GPU Miner [START_REF] Fang | Parallel Data Mining on Graphics Processors[END_REF]. Y.Li implementation [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] is 4x to 8x faster than UV Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF] when performing on an NVIDIA GTX 280 GPU.

In this paper, we propose an efficient parallel implementation of the HBK binarization method. Our approach is established through two levels of parallelism. The first level is the coarse-grained parallelism when parallelizing the Kmeans algorithm between thread blocks. The second level is the fine-grained parallelism, in which, we parallelize the distance computation and the labeling Kmeans stages between threads. Our HBK implementation ensures an efficient memory management and limits the CPU-GPU memory transfer overheads.

The rest of the paper is organized as follow. In Sect.2 we explain the HBK binarization method and the related quality and real time evaluations. Sect. 3 provides a brief presentation of the GPU architecture. Next, in Sect. 4, we review the related works that parallelized Kmeans on GPU and we describe our proposed HBK GPU implementation. Then, experimental results are discussed in Sect. 5. Finally, conclusion is drawn in Sect.6.

Hybrid Binarization based on Kmeans (HBK)

The HBK method combines local and global approaches to ensure robust binarization. Indeed, it is well adapted for noisy documents thanks to the local approach. Moreover, the global process improves the binarization quality when merging the local binarization features. Figure 2 shows the HBK method and more details are described in Algorithm 1.

To binarize documents, two initial centroids are employed. We note by C G = {C G0 , C G1 } the global centroid set. They are initialized with black and white values. Thus, C G0 = (0, 0, 0) T and C G1 = (255, 255, 255) T design respectively the background and the foreground clusters. The input document is divided into N b blocks. The size of the blocks can be small, medium or large depending on the size of characters [START_REF] Soua | A new hibrid binarization method based on Kmeans[END_REF]. Moreover, each block includes two local centroids given by the C L set:

C L = {(C L0,0 , C L0,1 ), (C L1,0 , C L1,1 ), ..., C Lbl,i , ..., (C LN b-1,0 , C LN b-1,1 )} with C Lbl,i = (C R Lbl,i , C G Lbl,i , C B Lbl,i ) T , i ∈ [0, 1], bl ∈ [0, N b -1]
. Furthermore, each block includes two local accumulators and counters given respectively by the sets Accu L and Count L :

Accu L = {(Accu L0,0 , Accu L0,1 ), (Accu L1,0 , Accu L1,1 ), ..., Accu Lbl,i , ..., (Accu LN b-1,0 , Accu LN b-1,1 )}, Count L = {(Count L0,0 , Count L0,1 ), (Count L1,0 , Count L1,1 ), ..., Count Lbl,i , ..., (Count LN b-1,0 , Count LN b-1,1 )} with Accu Lbl,i = (Accu R Lbl,i , Accu G Lbl,i , Accu B Lbl,i ) T , where, i ∈ [0, 1] and bl ∈ [0, N b -1].
In the beginning, all local centroids are initialized with C G0 and C G1 . Then, Kmeans(bl, CG) 5:

for i ∈ [0, 1] do 6:

AccuGi ← AccuGi + Accu Lbl,i 7:

CountGi ← CountGi + Count Lbl,i 8:

end for 9:

end for // New centroid computation 10:

CGi ← Accu Gi Count Gi , i ∈ [0, 1] 11: end while the Kmeans clustering algorithm [START_REF] Lloyd | Least square quantization in PCM[END_REF] is applied across all blocks as stated in lines 3 and 4. When it converges, the local accumulators Accu L are gathered into global ones defined by

Accu G = {Accu G0 , Accu G1 }, Accu Gi = (Accu R Gi , Accu G Gi , Accu B Gi ) T with i ∈ [0, 1] (line 6
). This operation is based on one addition per each color component. Similarly, local counters are accumulated into the global ones given by Count G = {Count G0 , Count G1 } (line 7). Finally, the global centroids C G are computed according to Accu G and Count G division (line 10). The algorithm stops when the global convergence is reached, else it reiterates while local centroids C L are reset with C G0 and C G1 . Following, we show the robust binarization quality of the HBK method through an OCR evaluation.

HBK OCR-based Evaluation

The Optical Character Recognition (OCR) is a process by which text characters contained in an image can be recognized. The binarization quality has a direct influence on the OCR result. Commonly, the computer uses an OCR Engine to recognize characters in an image. In our evaluation, we employ the well-known Tesseract 3.02 engine [START_REF] Smith | An Overview of the Tesseract OCR Engine[END_REF]. The OCR is performed directly on binarized scanned documents. For this we use the LRDE-DBD1 dataset [23]. It is composed of 125 scanned and nonscanned documents extracted from "Le Nouvel Observateur 2 " magazine. These documents have A4 format and 300-dpi resolution. According to HBK, the block size of 32x32 guarantees a robust binarization quality [START_REF] Soua | A new hibrid binarization method based on Kmeans[END_REF]. Table 1 shows the OCR recognition rate including the HBK and several binarization methods.

The OCR results are close: The recent SauvolaMs kx [START_REF] Lazzara | Efficient Multiscale Sauvolas Binarization[END_REF] gives an acceptable OCR rate by scoring 89% of

Methods

Accuracy (%)

HBK [START_REF] Soua | A new hibrid binarization method based on Kmeans[END_REF] 91 SauvolaMS kx [START_REF] Lazzara | Efficient Multiscale Sauvolas Binarization[END_REF] 89 Lelore [START_REF] Lelore | Super-resolved binarization of text based on the FAIR algorithm[END_REF] 85 Niblack [START_REF] Niblack | An Introduction to Digital Image Processing[END_REF] 80 TMMS [START_REF] Fabrizio | Text segmentation in natural scenes using toggle-mapping[END_REF] 73

OCR accuracy. It ensures a robust text binarization however in some document areas, artifacts may appear leading to character miss-recognition. The HBK method gives the best OCR rate scoring 91% of accuracy thanks to its robust binarization results. Indeed, the local approach gives robust binarization quality and the global approach clears the artifacts generated by the local one. The main issue encountered by outperformed methods is related to large objects and the fact that parameters do not fit to the contents. The HBK approach succeeds in improving the SauvolaMS kx algorithm because text with colored background is correctly retrieved. Indeed, according to our experiments, it ensures an F-Measure [START_REF] Chen | On the statistical properties of the Fmeasure, Quality Software[END_REF] of 98% while SauvolaMS kx reaches only 95% [START_REF] Lazzara | Efficient Multiscale Sauvolas Binarization[END_REF], with F-Measure is a metric that evaluates the pixel-based accuracy. The TMMS [START_REF] Fabrizio | Text segmentation in natural scenes using toggle-mapping[END_REF] method gives degraded results despite the efficient binarization of the text edges. This is due to the uneven illumination introduced by the scanned documents.

In this work, we aim to perform a binarization that can be integrated in a real time OCR application. In the next subsection, we check the HBK ability to satisfy the real time constraint of our paperless application.

HBK Real Time Evaluation

The experiments were conducted on an Intel Core i3 CPU performing at 3.07 GHz with a memory of 4 GB. Evaluations were performed using an average of 5 subsequent executions. In our target paperless application, the Epson Perfection V100 Photo scanner digitizes one A4 300 dpi document per 23 seconds. Actually, when the scanner digitizes a document, the OCR is performed simultaneously on the previous scanned document. Hence, the OCR time constraint is equal to the scan time (23 seconds). The HBK method binarizes one document per 1.9 seconds. Unfortunately, this performance does not meet to real time constraint, because it exceeds 2% [START_REF] Srinivasan | Performance characterization and acceleration of Optical Character Recognition on handheld[END_REF] from the total OCR execution time required. Following, we profile the HBK method to overview which processed parts are the most time consuming. Three parts were profiled: Initialization, Local phase and Global phase. The initialization includes memory allocations and data initialization. The local phase process includes the running time of Kmeans in all blocks of the document. The global phase process includes the check of the global , the Graphic Processing Unit (GPU) develops continuously and provides a promising platform for parallelizing K-means [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF]. Therefore, the GPU is an adequate alternative to accelerate Kmeans, and thus HBK. In the next section, we present the GPU architecture.

GPU Architecture

GPU is a plenty powerful many-core processor that support parallel data processing and high precision computation. It includes a set of Single Instruction Multiple Data (SIMD) Streaming Multiprocessors (SM). The SMs perform simultaneously the same operation on multiple data. Each SM is composed of several processors (cores). GPUs include a texture, constant, global, local, shared and register memories as shown in Figure 3 [START_REF]CUDA C best practices guide (version 4.0)[END_REF]. GPU grid, blocks and threads are explained a bit later. GPUs have dedicated texture hardware optimizing the 2D spatial locality. In addition, GPUs include a constant memory. It is a fast memory readable by all GPU cores. We note that, constant and texture memories are read-only which constitutes a major drawback if we need to write data. The main memory in GPUs is the global one. It can be accessed from all device cores. Depending on the GPU model, it can accommodate one Giga Bytes (GB) or more which let process large amount of pixels. The global memory is optimized to simultaneously handle a lot of memory requests providing a considerable large bandwidth. The access to this memory is cached which makes it very fast [START_REF]NVidia CUDA C Programming Guide[END_REF]. Besides, local memory is an abstraction of global one. It is accessed only by thread blocks. Moreover, it resides off chip, which makes it as expensive to access. The compiler makes use of local memory when it determines that there is not enough register space to hold variables. In the other side, the GPU includes shared memory which is 100x faster than global one. It can be used as a user-managed cache L1 to reduce the number of slow global memory accesses. This kind of memory is common to a group of cores and can be concurrently accessed by several ones. In such case, cores may collaborate in a given task, or simply store temporally data and reduce then register pressure. Eventually, the fastest memory access in GPU is provided by registers. Unfortunately, it has small memory storage. In the case of Fermi architecture, there are about two Mega Bytes (MB), while current CPUs only have few Kilo Bytes (KB).

To exploit the GPU platform, it is possible to use multiple technologies. The most popular ones are the Open Computing Language (OpenCL) [START_REF] Khronos | OpenCL: The open standard for parallel programming of heterogeneous systems[END_REF] and the Compute Unified Device Architecture (CUDA) [START_REF]Parallel programming and computing platform, CUDA[END_REF]. In one side, OpenCL is an open standard offering cross-platform capabilities. The most important drawback of OpenCL is the high complexity of the Application Programming Interface (API) and the execution model [START_REF] Tompson | An introduction to the programming Model[END_REF]. In the other side, CUDA is an API that allows high C and C++ programming level on NVIDIA GPUs. In addition, the CUDA environment is more user friendly and able to port the programming code much more consistently. As shown in Figure 3, CUDA notes the SMs and the cores respectively by blocks and threads. The container of all CUDA blocks is called grid. According to the CUDA environment, a virtual architecture is considered defining an execution model as described in Figure 4 [START_REF]NVidia CUDA C Programming Guide[END_REF].

The execution model is then well suitable for parallelizing treatments in image processing because blocks group thousands of independent threads ready to process a huge number of pixels in parallel. In the next section, we explain our HBK GPU parallelization.

HBK Parallelization on GPU

The previous analysis, in section 2.2 shows that the HBK execution time is bounded with the Kmeans clustering one. For this, we aim to accelerate Kmeans. Following, we introduce the related Kmeans GPU acceleration methods. Then, we detail our proposed HBK implementation strategy based on our Kmeans GPU version.

Related Works

Following, we give an overview about existing Kmeans implementations on the GPU. The first fast and very popular Kmeans GPU implementation is the UV Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF]. In this work, centroids and input data are stored respectively in the constant and texture memories. Each thread processes a single data point and the new cluster centroids are computed in the CPU device. In the UV-Kmeans implementation, the accumulation operations, are done on the CPU. However it increases the memory latency when transferring data between CPU and GPU. In addition, faster memories such as registers and shared memory are not used. W. Fang [START_REF] Fang | Parallel Data Mining on Graphics Processors[END_REF] accelerates 5x UV-Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF] by minimizing memory latency transfer between CPU and GPU. Indeed, he uses a bitmap grid to determine, inside GPU, the number of data points belonging to each cluster. This technique avoids large amount of atomic additions. However, the grid bitmap uses a lot of memory when the centroid number is high. Moreover, the memory management is not optimal. Indeed, registers and shared memory, are not employed.

Recently, works [START_REF] Hong-Tao | KMeans on Commodity GPUs with CUDA[END_REF][START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] perform more Kmeans steps on the GPU. For example in [START_REF] Hong-Tao | KMeans on Commodity GPUs with CUDA[END_REF], authors use the CPU only for centroid initialization and the control flow. Also, in [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF], several reduction steps are done on GPU before performing the centroid update stage on the CPU. This method, improves the works [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF] and [START_REF] Fang | Parallel Data Mining on Graphics Processors[END_REF] by optimizing the memory management. Actually, the input data is stored in registers that allows very fast memory access. Y.Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] designs an algorithm that exploits GPU on-chip fast memories to significantly decrease the data access latency.

Some works perform the centroid updating on the CPU [START_REF] Hong-Tao | KMeans on Commodity GPUs with CUDA[END_REF]18]. This usage is due to the employment of GPU architecture that are unable to achieve atomic floating point operations. For example, in the HP Kmeans [START_REF] Wu | Clustering billions of data points using GPUs[END_REF] method, the centroid updating step is performed on the CPU, but some accelerations are achieved like asynchronous transfer and cuda streams. In [18], researchers have implemented the Kmeans algorithm in the context of image segmentation. They parallelize only the labeling phase. However, the updating of centroids is performed on the CPU. Indeed, they employ a small centroid number compared to pixel number [18]. Hence the data transfer time between CPU and GPU is masked by the huge distance computation task.

In our work, we perform an efficient HBK GPU implementation based on a powerful Kmeans algorithm. Actually, HBK processes documents according to small blocks. Then, we apply Kmeans independently on each block. Following section explains our HBK GPU implementation.

HBK GPU Implementation

In our implementation, we perform the coarse and the fine grained levels of parallelism: In the coarse-grained level, we parallelize the document block treatment between blocks. In the fine-grained level, we parallelize the pixel processing inside blocks [START_REF]NVidia CUDA C Programming Guide[END_REF]. Figure 5 shows this concept.

Actually, in our HBK implementation, we parallelize Kmeans between the GPU blocks. In the other hand, we parallelize the distance computation and the labeling stages between threads (Figure 6). For this aim, we design two CUDA kernels: The Local and the Global Kernels. They handle respectively the local and the global HBK phases (Figure 2). All of the centroids C G , accumulators Accu G and counters Count G are allocated in the global memory. After the CPU fetches image into the global memory, it launches the Local Kernel on the GPU device. Thus, the Kmeans algorithm is performed on each block, in parallel.

An efficient memory management is ensured by employing fast GPU memories. We limit the global memory utilization by dispatching the data to shared memory and registers. Indeed, the local centroids C L , accumulators Accu L and counters Count L of each block In the Kmeans initialization step, each two local centroids C L are initialized with the RGB values (0, 0, 0) T and (255, 255, 255) T . We note that, the distance computation between centroids and pixels is performed independently between pixels. Therefore, we process it in parallel between w threads. Then, the complexity decreases from O(N pxk) to O(N pxk/w) with N p is the total number of pixels in the scanned document and k is the centroid number in each block. We should carefully map threads and pixels depends on their positions in the image. Moreover the thread number should match the total number of pixels to ensure a reliable binarization.

Similarly, the labeling step of Kmeans, where each pixel is assigned to its cluster, is parallelized on w threads.

In our HBK implementation, we decrease the CPU-GPU communication overheads by performing a fully Kmeans in the GPU. Indeed, the centroid updating is performed on the GPU device. One master thread in each block perform this operation. Actually, the scanned document contains high number of blocks because their sizes are relatively small compared to the image one. This way, the total number of master threads in the whole GPU remains high. Then, GPU resources are efficiently exploited. In addition, updating centroids on the GPU is possible thanks to the availability of integer and floatingpoint atomic addition on recent GPU architectures like Fermi and Kepler [START_REF]NVIDIAs Next Generation CUDA TM Compute Architecture: Kepler TM GK[END_REF].

During the labeling stage of the Kmeans algorithm, a reduction is performed on each block. Thus, local features Accu L and Count L are computed before being gathered into the global ones: Accu G and Count G , as seen in Figure 7. Indeed, two steps are made. Firstly, local features are gathered using the shared atomic addition within the threads blocks. Secondly, all local features are added to the global ones using the global atomic additions. This technique decreases the global atomic operation overhead. Indeed atomic operations on global memory can be very costly, as they need to serialize a potentially large number of threads in the kernel. To reduce this overhead, we usually addition operations first to variables declared in the shared memory of each thread block [START_REF] Yi | Efficient parallel CKY parsing using GPUs[END_REF]. Every block process stops when centroids of the current iteration are equal to the previous iteration ones. This means that the convergence criterion of Kmeans is reached. At that time, the CPU launches the Global Kernel. We do not need to transfer the accumulators Accu G and the counters Count G to the CPU because we keep them on the global memory. Based on these features, the Global Kernel computes centroids C G . Similarly to the local convergence, global centroids are submitted to the convergence test. If the centroids of the current global iteration are equal to those of the previous one, the global convergence is reached. In this case, a flag is set to 1, otherwise it is set to 0. Subsequently, the CPU checks whether the computed flag is equal to 1 value to decide ending the program as shown in Figure 8. The transfer of the flag variable maintains a low CPU-GPU overhead which enhances the HBK performances. In the next section, we evaluate our HBK GPU implementation.

Experimental Results

In previous work [START_REF] Soua | A new hibrid binarization method based on Kmeans[END_REF], we have demonstrated that the HBK method is effective in binarizing scanned documents. In this section, we evaluate the execution time of HBK. After presenting the employed materials, we discuss the efficiency of our GPU-based Kmeans implementation. Then, we compare the HBK processing time between GPU and CPU to show the GPU acceleration. For this, several experiments were performed by varying parameters such as; block sizes, amount of noise (Scanned, NScanned documents) and document content (text areas, colored background). Finally, we evaluate the HBK binarization quality on CPU and GPU to show our implementation reliability.

Materials

In this evaluation, both scanned and non-scanned LRDE-DBD dataset [23] are employed. In the rest of the paper, we note the non-scanned documents by "NScanned". Our experiments were conducted on PC with an Intel Core i3 CPU performing at 3.07 GHz. In addition, we employ an NVIDIA GeForce GTX 660 GPU. This one has 5 SIMD multiprocessors. Each one contains 192 processors and performs at 980 MHz. In addition, every multiprocessor has 64 KB of configurable RAM including 16 KB of shared memory and 48 KB of L1 cache. The global memory has 2 GB with a peak bandwidth of 144.7 GB/s. Moreover, CUDA version 5.0 is used. In order to visualize clearly the speedup effect, the HBK execution time is computed after input and output data transfer between CPU and GPU. Before we evaluate our HBK GPU implementation, we study our fully GPU-based Kmeans. We name this implementation Kmeans SD because it handles small datasets. Then, we adapt our Kmeans parallelization strategy to handle large datasets.

GPU-based Kmeans Time Evaluation

The Kmeans acceleration has been widely discussed in the literature [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF]. Table 3 shows a comparative study between several GPU based Kmeans implementations such as; Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF], HP Kmeans [START_REF] Wu | Clustering billions of data points using GPUs[END_REF], UV Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF] and GPU Miner [START_REF] Fang | Parallel Data Mining on Graphics Processors[END_REF]. The minimum and maximum execution time of every method are represented according to [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF]. Large Data Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] 0,22 -2,26 0,34 -1,15 HP Kmeans [START_REF] Wu | Clustering billions of data points using GPUs[END_REF] 1,45 -9,03 -UV Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF] 2,84 -34,54 1,86 -8,67 GPU Miner [START_REF] Fang | Parallel Data Mining on Graphics Processors[END_REF] 61,39 -474,83 4,26 -40,6

We can see from this study that the Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] outperforms all of the evaluated methods in both small and large data sizes. We propose in the following section to compare our Kmeans GPU based implementation to the Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] method which is the best implementation. Evaluation was conducted within different data sizes. We call our Kmeans implementation, respectively, Kmeans SD and Kmeans LD according to the use of small or large datasets.

Kmeans Time Evaluation on Small Dataset (Kmeans SD)

HBK relies on parallelizing Kmeans on small datasets. For this data size, Kmeans SD fully processes on the GPU. However, Kmeans Li processes the centroid updating on the CPU. In this subsection, we note Kmeans SD and Kmeans Li respectively by Kmeans-based GPU update and Kmeans-based CPU update. We compared the two Kmeans implementations on 8x8, 16x16 and 32x32 cropped documents. According to Figure 9, our Kmeansbased GPU update provides the lowest execution time in all cropped documents. Actually, the CPU-GPU memory transfer cost is high on small data compared to the computation gained benefits. Thus, it is better to keep the centroid update stage of Kmeans on GPU. Indeed, our Kmeans SD succeeds in avoiding a valuable CPU-GPU data transfer time.

In addition, we observe that the smaller processed data size is, the lower Kmeans execution time is. Thus, 32x32 data size provides the highest execution time, however, 8x8 data size provides the lowest one. In the following, we adapt our Kmeans parallelization strategy to handle large dataset and show its efficiency.

Kmeans Time Evaluation on Large Dataset (Kmeans LD)

Commonly, the Kmeans algorithm is processed on large datasets like A4 documents. In this context, Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] provides one of the most efficient Kmeans implementations. In this work, the data is divided into blocks of 256 data points. Moreover, centroids are stored in the constant memory. In the beginning of the algorithm, the data is dispatched once from the global memory to registers. Each thread computes the distance between each data point and all the centroids. Thus, it determines the closest center to the corresponding point. For the new centroid computation, a divide and conquer strategy is employed. Indeed, the data is divided into n groups, where n is initialized with n M . With n is the data size belonging to the corresponding centroid, and M is a multiple of stream multiprocessor number in the GPU. Each group is reduced to get a temporary centroid. n is then updated, by dividing the temporary centroids into n groups and reduce them iteratively on the GPU until n is smaller than M , indicating the GPU has no advantage than the CPU for further computing. In the end, the final centroids are computed on the CPU.

We adapt our Kmeans parallelization strategy to handle large datasets. Then, we compare this Kmeans LD implementation to the Kmeans Li one. Actually, we use registers to store data points. Registers provide fast memory access during the process execution which will accelerate Kmeans. In another hand, we use shared memory to store the centroid values. In the beginning, each thread determines its closest centroid. Then, centers are computed directly from the first pixel distribution on the GPU blocks and store them in the shared memory. Finally, we reduce these centroids on the GPU until we get the final ones.

Y.Li evaluated his implementation using a dataset of 4 million points on a GTX 280 GPU architecture and an i5 CPU. We have implemented Kmeans Li on our available hardware given in Section 5.1. In addition, we used the LRDE-DBD dataset composed with images of 8 million pixels which is larger than the dataset employed by Y.Li. Table 4 shows the time comparison between Kmeans LD and Kmeans Li on the LRDE-DBD dataset. Kmeans LD 425 Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] 600

Kmeans LD is 1.45x faster than the Kmeans Li. Actually, in Kmeans LD we do not redistribute the pixels on the GPU after distance computation. Therefore, we save processing time and ensure a coalescent memory read. However in the Y.Li Kmeans implementation, after distance computation the data is redistributed on the GPU to reduce the temporary centroids which is time consuming. Indeed, the GPU is setup with a new block number and a new pixel indexes are computed and dispatched to threads. Hence, pixel indexes become disorderly and cause a non coalescent memory read. In addition, in Kmeans LD all centroids are computed in parallel inside blocks unlike Kmeans Li that do rearrangement subsequently for each centroid.

Table 5 shows a comparison between our Kmeans GPU-based implementation and several state of the art methods in which we describe their strengths and areas of improvement. Actually, our method outperforms Kmeans Si [18] because we parallelize efficiently the centroid updating stage on the GPU. The same strategy is done by UV Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF] but the main limitation encountered, as in GPU Miner [START_REF] Fang | Parallel Data Mining on Graphics Processors[END_REF] and Kmeans Ta [START_REF] Hong-Tao | KMeans on Commodity GPUs with CUDA[END_REF], is the lack of fast memory use. However, our method employs fast memories like registers and shared memory. Moreover, we ensure a Kmeans implementation that is more efficient than Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] according to two sides. Firstly, we reduce the transfer time between CPU and GPU when processing all Kmeans steps on the GPU for small datasets. Secondly, on the large datasets, we ensure more coalescent global memory access than Processing very large data Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] Decreases memory latency by using fast registers

Costly data transfer between CPU and GPU GPU Miner [START_REF] Fang | Parallel Data Mining on Graphics Processors[END_REF] Avoids large atomic operations Poor use of fast GPU memories UV Kmeans [START_REF] Che | A Performance Study of General-Purpose Applications on Graphics Processors Using CUDA[END_REF] Uses the cache mechanism for high speed memory access

Fast memories such as registers and shared are not used Kmeans Ta [START_REF] Hong-Tao | KMeans on Commodity GPUs with CUDA[END_REF] Parallelizes all complex Kmeans steps on the GPU including the centroid updating stage

The use of global memory to permanently store centroids and data increases the memory latency

Kmeans Si [18]

Offloads the Kmeans labeling stage to the GPU Do not parallelize the Kmeans centroid updating stage

Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] because we do not rearrange data pixels. Thus, pixel indexes remain ordered in each block.

In addition, all centroids are processed in parallel unlike Kmeans Li [START_REF] Li | Speeding up k-Means Algorithm by GPUs[END_REF] that processes and reduces one centroid at a time. In fact, our Kmeans implementation handles efficiently small and large datasets. Indeed, it decreases the memory transfer latency by using a carefully designed workload between CPU and GPU. Moreover, it ensures a coalescent and fast memory access thanks to the effective memory management. Following, we evaluate our HBK GPU implementation which includes the studied Kmeans proposal.

GPU-based HBK Time Evaluation

In the first part of this section, we evaluate the HBK GPU implantation to show the speedup effect when compared to HBK CPU-based one. For this, we study our HBK GPU implementation to determine the suitable block size that satisfies our paperless real time constraint. Moreover, the time performance of our GPU implementation is evaluated on different document parameters (noise amount and content). In the second part of this section, we compare the HBK GPU real-time score with that of different binarization GPU-based methods in the literature.

GPU-based HBK Speedup

Following, we evaluate our proposed HBK GPU implementation performances. In the beginning, we vary the block size parameters to see its effect on the GPU speedup. In this context, Table 6 shows the HBK execution time average across 125 Scanned documents.

We observe that the larger block size is, the less is the HBK processing time on CPU and the greater it is on GPU. Moreover, the speedup of the HBK GPU implementation decreases when increasing the size of the blocks. Actually, the execution time on the GPU is faster when the number of pixels in the image blocks is close to the physical number of threads in the GPU blocks. In this case, each pixel is processed directly by a thread. Otherwise, the pixel follows a queue until a resource (thread) gets free. Moreover, the number of divergence branches is lower using smaller data. In another hand, by scoring 425 ms and 388 ms, the HBK is able to be integrated in the OCR paperless application that require to binarize one document per 460 ms. For these reasons, we consider that the block size of 16x16 pixels is the most suitable for our paperless application because it ensures both fast and acceptable binarization quality [START_REF] Soua | A new hibrid binarization method based on Kmeans[END_REF].

Next, we discuss the HBK time performance on GPU when varying the amount of noise. Figure 10 visualizes the HBK processing time on the CPU and the GPU architectures, represented respectively by purple color and dotted red one. This evaluation is performed across 125 Scanned and NScanned LRDE-DBD documents. We note that, the noise is greater on the Scanned documents (Figure 10.b) than on the NScanned ones (Figure 10.a) because of the degradation in the scan process. The execution time varies greatly between the two document types on CPU. However, the process remains nearly constant on GPU. Indeed, the noise complexity induced in the Scanned documents requires more computation tasks which are efficiently parallelized on the GPU. To make a relevant comparison, we consider the HBK average time on the LRDE-DBD 125 documents as shown in Table 7. The GPU version is 6x faster than the CPU one. In the other hand, the HBK GPU implementation is faster on NScanned documents, reaching 347 ms compared to 425 ms on Scanned ones. Indeed, the difference between execution time on Scanned and NScanned documents is due to the process complexity that affects the overall time. In fact, unlike the CPU implementation, our HBK parallelization on GPU is slightly sensitive to the noise introduced in documents after the scanning process.

Following, we evaluate the HBK processing time while varying the document content (number of text areas, colored background). In, Figure 11, the document1 contains a few text areas. The document2 includes more text areas and document3 contains many text areas and a colored background. The CPU time increases across the three document types. The more a document includes text and background the more processing time is. It reaches 4000 ms on document3 compared to 1100 ms in document1. The GPU execution time still nearly constant around 425 ms over the three document types. Moreover, the GPU speedup increases with the increase of the document complexity. The GPU reaches a speedup of 7,5x when binarizing the document3 compared to 3,5x when binarizing the document1. Indeed, the more information documents include, the more devices requires computation tasks. The GPU scores higher speedup because it is efficient for parallelizing high amount of computation. Indeed, the GPU implementation of HBK is able to establish a binarization time lower than our paperless application real time constraint of 460 ms, whatever are the noise or contents in the documents. Following, to establish the conceptual merit of our HBK implementation, we compare its real time with a set of binarization methods on GPU.

Real-time evaluation of HBK and different binarization methods

Several binarization methods where studied in our previous work [START_REF] Soua | A new hibrid binarization method based on Kmeans[END_REF]. A few of them have been implemented on GPU in the literature. In this section, we compare our GPU-based HBK to the GPU implementations of Niblack [START_REF] Niblack | An Introduction to Digital Image Processing[END_REF], Sauvola [START_REF] Sauvola | Adaptive Document Binarization, 4th Int. Conf. On Document Analysis and Recognition[END_REF] and Sauvola M skx [START_REF] Lazzara | Efficient Multiscale Sauvolas Binarization[END_REF] algorithms. We give a brief description of each algorithm and its GPU implementation.

Niblack [START_REF] Niblack | An Introduction to Digital Image Processing[END_REF] is a local thresholding algorithm. It computes a pixel-wise threshold by sliding a rectangular window over the gray level image. The threshold value for each pixel is decided by local mean and local standard deviation over a specific window size around each pixel. Thus, the local threshold T (x, y) for pixel (x, y) is calculated by formula:

T (x, y) = m(x, y) + k.s(x, y) (1)
Where m(x, y) and s(x, y) are the local mean and the local standard deviation of the pixels within the local window region. The value of k controls the amount of is that it always identifies the text regions correctly as foreground but on the other hand tends to produce a large amount of binarization noise in non-text regions [START_REF] Khurshid | Comparison of Niblack inspired Binarization methods for ancient documents[END_REF]. Sauvola [START_REF] Sauvola | Adaptive Document Binarization, 4th Int. Conf. On Document Analysis and Recognition[END_REF] proposed an algorithm similar to Niblack [START_REF] Niblack | An Introduction to Digital Image Processing[END_REF]. It made some assumptions based on the distribution of grey values associated with foreground and background pixels. The Sauvola threshold is computed as:

T (x, y) = m(x, y) + [1 + k.( s(x, y) R -1)] (2)
Where m(x, y), s(x, y) and k are the same parameters given in Niblack. R is the dynamic range of standard deviation. B.Singh parallelized Sauvola on GPU [START_REF] Singh | Parallel Implementation of Souvola's Binarization Approach on GPU[END_REF]. This GPU implementation achieved an average speed-up of 20.8x compared to the serial program. Both Niblack and Sauvola GPU versions employ the same parallelization strategy. Indeed, the input image is stored as texture in the device memory. After that, block and grid sizes were computed according to the dimensions of the input images. A single thread calculates the threshold for a single pixel in the output image. In our work, we have implemented Niblack and Sauvola on our GTX 660 GPU device.

Sauvola M skx [START_REF] Lazzara | Efficient Multiscale Sauvolas Binarization[END_REF] was proposed to improve Sauvola to handle the multi-scale text. It is composed of four steps. First, the image is sampled on different scales. Then, each image is binarized according to the conventional Sauvola method. Next, a threshold is set for each object according to its scale belonging. A final image including thresholds is produced and then binarized. According to our knowledge, Sauvola M Skx was not implemented on the GPU yet. Thus, we have implemented our own version. Indeed, we parallelized the threshold computation on each scale and we kept the other steps performing on the CPU.

Generally, the adjustment of parameter k and w of these three methods requires prior knowledge about the set of documents. For our evaluation, we employ the parameters recommended by [START_REF] Lazzara | Efficient Multiscale Sauvolas Binarization[END_REF]. Indeed, we adjust w = 51 and k = -0.2 for the Niblack method. For Sauvola and Sauvola M Skx we use w = 51 and k = 0.34. Table 8 shows the real time comparison between HBK, Niblack, Sauvola, and SauvolaM skx methods.

The evaluation result shows that HBK gives the best execution time. Indeed, it is 3.9x faster than Sauvola M skx and 1.4x faster than Niblack and Sauvola. Actually, our HBK GPU-based implementation is able to binarize pixels rapidly because all thread needed vari- ables are stored in the fast on-chip registers and shared memories. In the other hand, both Niblack and Sauvola need to compute a threshold for each pixel before binarizing the data. For this, they require fetching a great number of neighbor values from the texture memory to compute the mean and the standard deviation for each pixel. Otherwise, HBK employs faster memories compared to Sauvola M skx which uses further huge texture accesses to compute pixel thresholds on several scales. Therefore, HBK outperforms Sauvola M skx on the GPU. Niblack and Sauvola times are close, this is because the workflow is similar. However, Sauvola is slower because of its more complex threshold formula. Sauvola M Skx gives the slowest execution time. Indeed, the multi-scale processing implies a cost on computation time. Following, we evaluate the binarization quality of our HBK GPU-based and CPU-based implementations to show the reliability of our implementation.

GPU-based HBK Quality Evaluation

Generally, the GPU architecture has some memory constraints and model execution rules to make faster implementations. According to Figure 12, the HBK GPU version does not show any loss in the binarization quality compared to the CPU one. To prove this visual result, we compute the pixel based accuracy of the HBK on both CPU and GPU. The result gives the same F-Measure, for both devices. Indeed, our GPU program does not encounter any memory limitation or conflict that produces different binarization results.

Conclusion

In the paperless applications, the OCR systems are used to recognize text in the digital documents. In our work, we focused on the paperless real time application, in which, a scanner processes one A4 300 dpi paper per 23 seconds. The binarization is a very important component in the OCR tool chain. In this context, recently, we have proposed the HBK binarization method. It offers a high OCR accuracy compared to well known binarization methods. However, HBK does not respect our paperless real time constraint. It processes one document per 1.9 seconds, exceeding the binarization real time constraint of 460 ms. For this, we have proposed a parallel implementation of the HBK method on the Graphic Processing Unit (GPU). Our proposed CUDA implementation is based on combining fine-grained and coarse-grained parallel strategies and uses efficient memory management. According to our experiments, the GPU implementation of HBK performs one document per 425 ms which satisfy the paperless real time application constraint while offering a high OCR accuracy using 16x16 block size. Indeed, we have effectively avoided the communication overhead between the CPU and the GPU thanks to the full processing of Kmeans algorithm on GPU. In addition, our HBK GPU implementation is slightly sensitive to noise and document contents. The performed GPUbased comparisons of HBK to recent and well known binarization methods show that our proposed strategy implementation has a better performance and runs on GPU much faster than the compared methods.

Fig. 1

 1 Fig. 1 Paperless application toolchain based on OCR system

Fig. 2

 2 Fig.2The HBK binarization method

Fig. 3

 3 Fig. 3 GPU Architecture in CUDA environment

Fig. 4 GPU

 4 Fig. 4 GPU Execution model in CUDA environment

Fig. 5

 5 Fig. 5 Fine-grain and coarse grain parallelism levels in the HBK GPU implementation

Fig. 6

 6 Fig. 6 Kmeans implementation strategies on GPU. a. The GPU thread blocks (multiprocessors) b. The Kmeans clustering algorithm

Fig. 7

 7 Fig. 7 Optimization of accumulator and counter variables computation

Fig. 8

 8 Fig.8The HBK implementation workflow

Fig. 9

 9 Fig. 9 Time evaluation of GPU Kmeans implementation with CPU and GPU based centroid update when processing different cropped documents

Fig. 10

 10 Fig. 10 Time evaluation of the HBK binarization on GPU and CPU. a NScanned documents. b Scanned documents

Fig. 11

 11 Fig. 11 Time evaluation of the HBK method on GPU and CPU when varying content. a Original scanned documents. b Processing time on GPU and CPU. c The speedup of GPU against CPU

Fig. 12

 12 Fig. 12 GPU (a) and CPU (b) based binarization quality comparison of the HBK method on cropped Scanned document

Table 1

 1 OCR accuracy evaluation of HBK and well-known binarization methods.

Table 2

 2 Profiling of the HBK stages

	Processes	Time (%)
	Initialization	0.001
	Local phase (Kmeans)	96.593
	Global phase	0.013
	convergence and the centroid updating. According to
	table 2, the overall HBK time execution is bounded with
	the Kmeans one.Indeed, reducing the Kmeans processing
	time will naturally decrease the HBK one. According to
	[18]

Table 3

 3 Kmeans Li comparison with several state of the art Kmeans GPU implementations

	Methods	Time (s)
		Small Data

Table 4

 4 Time evaluation of Kmeans LD and Kmeans Li on 125 (2516x3272) LRDE-DBD images

	Methods	Time (ms)

Table 5

 5 Strengths and areas of improvement comparison between different GPU-based kmeans implementations.

	Methods	Strengths	Areas for improvement
	Our Kmeans implemen-	Avoids costly CPU-GPU transfer time and
	tation	ensures a coalescent memory access

Table 6

 6 Time evaluation of HBK method on CPU and GPU with different block sizes

	Method Blocs (pix) Device Time (ms) Speedup
		8x8	CPU GPU	2716 388	-7x
	HBK	16x16	CPU GPU	2546 425	-6x
		32x32	CPU GPU	2411 1050	-2,29x

Table 7

 7 Time evaluation of the HBK method on CPU and GPU with different noise complexity images

	Method	Images	Device Time (ms) Speedup
	HBK	Scanned NScanned	CPU GPU CPU GPU	2546 425 1930 347	-6x -5.5x

Table 8

 8 GPU Real time evaluation of HBK and three binarization methods on the LRDE-DBD Dataset

	Methods	Time (ms)
	HBK	425
	Niblack	594
	Sauvola	595
	Sauvola M Skx	1660

Copyright(c)

EPITA and Development Laboratory (LRDE) with permission from Le Nouvel Observateur. LRDE-DBD is available online on the web site: http://www.lrde.epita.fr/cgi-bin/twiki/view/Olena/ DatasetDBD 2 Le Nouvel Observateur. Issue 2402, November 18-24, 2010 and available on the website: http://tempsreel.nouvelobs.com


T Bonald

M-A Tran

Flow vs. Time Sampling for Throughput Performance Evaluation

Keywords:

We introduce two throughput metrics referred to as flow-and time-sampled throughputs. The former gives the throughput statistics of an arbitrary flow while the latter weights these throughput statistics by the flow durations. Under fair sharing assumptions, the latter is shown to coincide with the steady-state instantaneous throughput weighted by the number of flows, which provides a useful means to measure and estimate it. We give some generic properties satisfied by both metrics and illustrate their difference on a few examples.

Introduction

Throughput is a key performance indicator for Internet access. Users perceive the quality of their Internet access through the time necessary to download and upload various files. Providers need throughput measurements to control quality of service and to take appropriate decisions about equipment upgrades. They also need simple analytical results that relate user throughput to the degree of network congestion in order to determine target traffic loads and to assess the efficiency of various control schemes like routing, scheduling and buffer management.

Although throughput is a very common performance metric, it is actually not immediately obvious how to measure and estimate it. One of the difficulties comes from the random nature of traffic. Throughput is most often the result of the statistical multiplexing of data flows that share the same network resources. These flows arrive at random and leave the network at the end of the corresponding document transfer. They adapt elastically to congestion, typically under the control of TCP [START_REF] Stevens | TCP/IP Illustrated[END_REF]. Thus throughput is inherently random, determined by the complex interaction of a dynamic set of elastic flows through the congestion control algorithms.

Another difficulty is related to the well-known observer's paradox. The simplest way to measure the throughput achieved by users during a file download is to actually download this file, which changes the network state. Similarly, calculating user throughput by means of a traffic model requires the study of the transient behaviour of the system, given the presence of the considered file transfer. In both cases, it would be much more convenient to derive throughput estimates from the steady state of the network, either measured through passive monitoring tools or calculated using an appropriate traffic model.

In this paper, we introduce two throughput metrics referred to as flow-and time-sampled throughputs. The former gives the throughput statistics of an arbitrary flow. It turns out to be hard to evaluate analytically even for the simplest traffic models. This is why it is most often approximated by the steady-state instantaneous throughput [START_REF] Berger | Dimensioning bandwidth for elastic traffic in highspeed data networks[END_REF][START_REF] Kherani | Stochastic models for throughput analysis of randomly arriving elastic flows in the Internet[END_REF][START_REF] Litjens | Throughputs in processor sharing models for integrated stream and elastic traffic[END_REF]. We shall see that this approximation is in fact far from accurate. The latter weights the throughput statistics by the flow durations. A key result of the paper is that, under fair sharing assumptions, this throughput metric coincides with the steady-state instantaneous throughput weighted by the number of flows. This yields explicit, insensitive expressions for a large class of traffic models [START_REF] Ben Fredj | Statistical bandwidth sharing: a study of congestion at flow level[END_REF][START_REF] Bonald | A queueing analysis of maxmin fairness, proportional fairness and balanced fairness[END_REF][START_REF] Delcoigne | Modeling integration of streaming and data traffic[END_REF][START_REF] Heyman | A new method for analysing feedback-based protocols with applications to engineering Web traffic over the Internet[END_REF][START_REF] Litjens | Elastic calls in an integrated services network: the greater the call size variability the better the QoS[END_REF].

We highlight the fact that both throughput metrics are random variables that represent the whole throughput distribution. This is particularly relevant for adaptive streaming flows that typically require a minimum throughput. Though we mainly focus on elastic traffic, the results may be used as a conservative approximation for a more realistic scenario where both elastic and adaptive streaming flows are multiplexed [START_REF] Bonald | On performance bounds for the integration of elastic and adaptive streaming traffic[END_REF]. It is also worth noting that many elastic flows in fact correspond to streamed audio and video files with long playback delays (a few seconds, say). It is clear that the mean throughput is not a good metric for assessing the quality of such flows. Asymptotic results on flow duration and throughput, see e.g. [START_REF] Bekker | Performance of TCP-friendly streaming sessions in the presence of heavy-tailed elastic flows[END_REF][START_REF] Borst | Tail asymptotics for discriminatory processor-sharing queues with heavy-tailed service requirements[END_REF][START_REF] Guillemin | Tail asymptotics for processor-sharing queues[END_REF], are also not sufficient for this purpose. What is really needed is either the throughput distribution during the flow, as considered in the present paper, or some bounds or approximations on the flow duration conditionally to the flow size, as described in [START_REF] Cheung | Insensitive bounds for the moments of the sojourn time distribution in the M/G/1 processor-sharing queue[END_REF] for instance.

We describe the traffic model in the next section. Sections 3 and 4 are devoted to the flow-and time-sampled throughputs, respectively. These throughput metrics are compared with the steady-state instantaneous throughput in Section 5. The results are illustrated by a few examples in Section 6. Section 7 concludes the paper.

Traffic model

Consider a tagged flow of potentially infinite size starting at time 0. We view this flow as a fluid stream and denote by ϕ(t) its instantaneous throughput at time t. This stochastic process is typically not stationary due to the impact of the tagged flow on the system state, which may be viewed as a manifestation of the observer's paradox. We assume it reaches a steady state characterized by some random variable ϕ ∞ equal to the weak limit of ϕ(t) when t tends to infinity. We also assume that the process ϕ(t) takes positive values in some discrete, bounded set X , is regular, right-continuous with left-hand limits, and ergodic in the sense that:

∀x ∈ X , lim T →∞ 1 T T 0 I {ϕ(t)=x} dt = Pr(ϕ ∞ = x) a.s. (1)
We denote by ϕ 0 = ϕ(0) the throughput of a virtual flow of null size.

If the tagged flow is of size s, its duration D(s) is the random variable defined by the equality:

s = D(s) 0 ϕ(t)dt.
We have the following simple result.

Proposition 1 The throughputs of flows of infinitely small and infinitely large sizes are respectively given by:

lim s→0 s D(s) = ϕ 0 a.s. (2)
and

lim s→∞ s D(s) = E[ϕ ∞ ] a.s. (3)
Proof. Let θ be the minimum time t > 0 such that ϕ(t) ̸ = ϕ(0). We have ϕ(t) = ϕ 0 for all t < θ so that

D(s) = s ϕ 0 if s < ϕ 0 θ.
We deduce (2) from the fact that ϕ 0 θ > 0 a.s. Now, since X is bounded, we have:

lim s→∞ D(s) = ∞ a.s. ( 4
)
so that (3) follows by ergodicity. ✷

In the following, we denote by S a random flow size and by D = D(S) the corresponding flow duration:

S = D 0 ϕ(t)dt.
3 Flow-sampled throughput

Definition

Assume the tagged flow has random size S. For all x ∈ X , consider the mean fraction of time the tagged flow has throughput x:

E ⎡ ⎣ 1 D D 0 I {ϕ(t)=x} dt ⎤ ⎦ . (5)
This defines a distribution on X . We refer to the corresponding random variable X flow as the flow-sampled throughput.

The definition extends to any given flow size s > 0. The corresponding random variable X flow (s) satisfies for all x ∈ X :

Pr(X flow (s) = x) = E ⎡ ⎢ ⎣ 1 D(s) D(s) 0 I {ϕ(t)=x} dt ⎤ ⎥ ⎦ . (6)

Properties

The mean flow-sampled throughput, denoted by γ flow , is given by:

γ flow = E[X flow ], = x∈X xE ⎡ ⎣ 1 D D 0 I {ϕ(t)=x} dt ⎤ ⎦ , = E ⎡ ⎣ 1 D D 0 x∈X xI {ϕ(t)=x} dt ⎤ ⎦ , = E ⎡ ⎣ 1 D D 0 ϕ(t)dt ⎤ ⎦ , = E S D . (7)
Similarly, we get the mean flow-sampled throughput of a flow of size s:

γ flow (s) = E[X flow (s)] = E s D(s) . (8)
We have the following result for the limiting cases of flows of infinitely small and infinitely large sizes:

Proposition 2 The flow-sampled throughput X flow (s) tends in distribution to ϕ 0 when s tends to 0 and to ϕ ∞ when s tends to ∞.

Proof. Recall that:

Pr(X flow (s) = x) = E ⎡ ⎢ ⎣ 1 D(s) D(s) 0 I {ϕ(t)=x} dt ⎤ ⎥ ⎦ .
As in the proof of Proposition 1, we denote by θ be the minimum time t > 0 such that ϕ(t) ̸ = ϕ(0). We have:

1 D(s) D(s) 0 I {ϕ(t)=x} dt = I {ϕ 0 =x} if s < ϕ 0 θ.
Since ϕ 0 θ > 0 a.s., we deduce:

lim s→0 1 D(s) D(s) 0 I {ϕ(t)=x} dt = I {ϕ 0 =x} a.s. and lim s→0 Pr(X flow (s) = x) = Pr(ϕ 0 = x).
The other limit follows by ergodicity, using (1) and ( 4):

lim s→∞ 1 D(s) D(s) 0 I {ϕ(t)=x} dt = Pr(ϕ ∞ = x) a.s. ✷ 4 Time-sampled throughput

Definition

Assume the tagged flow has random size S. For all x ∈ X , consider the ratio of the mean duration at throughput x to the mean flow duration:

E ⎡ ⎣ D 0 I {ϕ(t)=x} dt ⎤ ⎦ E[D] . (9)
This defines a distribution on X provided the mean flow duration is finite. We refer to the corresponding random variable X time as the time-sampled throughput.

The definition extends to any given flow size s > 0. The corresponding random variable X time (s) satisfies for all x ∈ X :

Pr(X time (s) = x) = E ⎡ ⎢ ⎣ D(s) 0 I {ϕ(t)=x} dt ⎤ ⎥ ⎦ E[D(s)] . ( 10
)
To illustrate the difference with the other throughput metric, consider K flows of size s, respective durations d 1 , . . . , d K and fractions of time at throughput x f 1 , . . . , f K . For the flow-sampled throughput metric, the probability that a flow of size s has throughput x is simply estimated as the empirical mean:

1 K K k=1 f k .
For the time-sampled throughput metric, it is estimated as the empirical mean weighted by the flow durations:

K k=1 d k f k K k=1 d k .
This corresponds to the fraction of time flows have throughput x. Note that, provided the K samples are i.i.d., the previous two expressions tend to ( 6) and [START_REF] Fayolle | Sharing a processor among many classes[END_REF], respectively, as K grows to infinity.

Properties

The mean time-sampled throughput, denoted by γ time , is given by:

γ time = E[X time ], = x∈X x 1 E[D] E ⎡ ⎣ D 0 I {ϕ(t)=x} dt ⎤ ⎦ , = 1 E[D] E ⎡ ⎣ D 0 x∈X xI {ϕ(t)=x} dt ⎤ ⎦ , = 1 E[D] E ⎡ ⎣ D 0 ϕ(t)dt ⎤ ⎦ , = E[S] E[D] . (11)
Similarly, we get the mean time-sampled throughput of a flow of size s:

γ time (s) = E[X time (s)] = s E[D(s)] . (12)
Again, to illustrate the difference with the other throughput metric, consider K flows of size s and respective durations d 1 , . . . , d K . The mean flow-sampled throughput is simply estimated as the empirical mean:

1 K K k=1 s d k ,
while the mean time-sampled throughput is estimated as the empirical mean weighted by the flow durations:

K k=1 d k s d k K k=1 d k = s 1 K K k=1 d k .
Provided the K samples are i.i.d., these expressions tend to ( 8) and ( 12), respectively, as K tends to infinity. They correspond to the arithmetic and harmonic empirical means, respectively.

We have the following result:

Proposition 3 For all flow sizes, the mean time-sampled throughput is less than the mean flow-sampled throughput:

∀s > 0, γ time (s) ≤ γ flow (s).
Proof. The proof follows from the convex inequality:

E[D(s)]E 1 D(s) ≥ 1.

✷

For the limiting cases of flows of infinitely small and infinitely large sizes, we get:

Proposition 4 The time-sampled throughput X time (s) tends in distribution to ϕ ∞ when s tends to ∞. Moreover,

lim s→0 Pr(X time (s) = x) = Pr(ϕ 0 = x) xE[1/ϕ 0 ] .
Proof. Recall that:

Pr(X time (s) = x) = E ⎡ ⎢ ⎣ D(s) 0 I {ϕ(t)=x} dt ⎤ ⎥ ⎦ E[D(s)]
.

By ergodicity, it follows from (1), ( 3) and (4) that:

lim s→∞ 1 s D(s) 0 I {ϕ(t)=x} dt = Pr(ϕ ∞ = x) E[ϕ ∞ ] a.s. and lim s→∞ E[D(s)] s = 1 E[ϕ ∞ ] , so that lim s→∞ 1 E[D(s)] D(s) 0 I {ϕ(t)=x} dt = Pr(ϕ ∞ = x) a.s.
For the other limit, we use again the fact that ϕ(t) = ϕ 0 for all t < θ so that 1 D(s)

D(s) 0 I {ϕ(t)=x} dt = I {ϕ 0 =x} if s < ϕ 0 θ.
Using (2), we deduce:

lim s→0 1 s D(s) 0 I {ϕ(t)=x} dt = 1 x I {ϕ 0 =x} a.s. and lim s→0 D(s) s = 1 ϕ 0 a.s. so that lim s→0 E ⎡ ⎢ ⎣ D(s) 0 I {ϕ(t)=x} dt ⎤ ⎥ ⎦ E[D(s)] = Pr(ϕ 0 = x) xE[1/ϕ 0 ] . ✷ 5

Steady-state instantaneous throughput

The stochastic process ϕ(t) involves the transient behaviour of the system and is therefore very hard to characterize even for the simplest traffic models.

It is tempting to derive flow-and time-sampled throughput statistics from the steady-state behaviour of the system, as observed in practice by means of passive monitoring.

We here consider a stream of tagged flows arriving at rate λ. As above, we denote by S their random size and by D their random duration. Let n be the number of active tagged flows in steady state. We assume a perfectly fair sharing among these flows and denote by ϕ their common instantaneous throughput, defined when n ≥ 1.

Conditional distribution

The steady-state instantaneous throughput, denoted by Y , is commonly used as an approximation for the flow-sampled throughput [START_REF] Kherani | Stochastic models for throughput analysis of randomly arriving elastic flows in the Internet[END_REF][START_REF] Litjens | Throughputs in processor sharing models for integrated stream and elastic traffic[END_REF]. Since it is defined if and only if n ≥ 1, its distribution is the conditional distribution:

∀x ∈ X , Pr(Y = x) = Pr(ϕ = x|n ≥ 1).
This throughput metric turns out to be biased: it depends on the proportion of tagged flows in the overall traffic stream. As shown in the example of §6.2 below, it is a loose approximation for the flow-sampled throughput when this proportion is equal to 1 and it coincides with the time-sampled throughput when this proportion tends to 0, which is in fact a simple consequence of the following more general result.

Weighted distribution

It turns out that the time-sampled throughput has the same distribution as the steady-state instantaneous throughput ϕ weighted by the number n of tagged flows. This weighting reflects the fact that congestion events impact the throughput metric in proportion to the number of active flows.

Proposition 5

The time-sampled throughput satisfies:

∀x ∈ X , Pr(X time = x) = E[nI {ϕ=x} ] E[n] .
Proof. The result is a simple consequence of Little's law [START_REF] Baccelli | Elements of Queueing Theory[END_REF]. We have:

E[n] = λE[D].
Moreover, considering those time periods where ϕ = x, we get

E[n|ϕ = x] = λ Pr(ϕ = x) E ⎡ ⎣ D 0 I {ϕ(t)=x} dt ⎤ ⎦ , so that: E[nI {ϕ=x} ] = λE ⎡ ⎣ D 0 I {ϕ(t)=x} dt ⎤ ⎦ . (13)

✷

In view of Proposition 5, we have:

γ time = E[X time ] = E[nϕ] E[n] . (14)
Using ( 13), we get:

E[nϕ] = E n x∈X xI {ϕ=x} , = x∈X xλE ⎡ ⎣ D 0 I {ϕ(t)=x} dt ⎤ ⎦ , = λE ⎡ ⎣ D 0 ϕ(t)dt ⎤ ⎦ , = λE[S],
which corresponds to the traffic intensity generated by tagged flows (in bit/s). We deduce the following expression for the mean time-sampled throughput:

γ time = λE[S] E[n] . ( 15
)
Note that this expression in fact directly follows from ( 11) and Little's law.

Examples

We now apply previous results to a few key examples. In the following, we simply refer to flows of infinitely small and infinitely large sizes as small and large flows, respectively.

Fair sharing

Consider a unit capacity link. Flows arrive as a Poisson process of intensity λ and have i.i.d. sizes characterized by the random variable S. Sharing is assumed to be perfectly fair so that the instantaneous throughput ϕ of each flow is equal to 1/n in the presence of n flows, with n ≥ 1. In particular, the throughput takes values in the discrete set X = {1, 1/2, 1/3, . . .}.

The corresponding queueing system is the processor-sharing queue. Provided the traffic intensity ρ = λE[S] is less than 1, the queue is stable and has the stationary distribution:

π(n) = ρ n (1 -ρ),
independently of the flow size distribution beyond the mean [START_REF] Serfozo | Introduction to Stochastic Networks[END_REF]. In the presence of a permanent flow, the stationary distribution of the number n of other flows becomes:

π ′ (n) = (n + 1)ρ n (1 -ρ) 2 .
The result is again insensitive to the flow size distribution beyond the mean. We deduce the following expressions for the throughput of a virtual flow of null size and for the throughput of a permanent flow of infinite size, where x = 1/m, m ≥ 1, denotes any element of X :

Pr(ϕ 0 = x) = π(m -1) = ρ m-1 (1 -ρ), Pr(ϕ ∞ = x) = π ′ (m -1) = mρ m-1 (1 -ρ) 2 .
In particular, we have:

E[ϕ 0 ] = - ln(1 -ρ) ρ (1 -ρ) and E[ϕ ∞ ] = E[1/ϕ 0 ] -1 = 1 -ρ. ( 16
)

Mean throughput

It turns out that the mean flow-sampled throughput has no simple expression, even in the simplest case of an exponential flow size distribution. We only have explicit expressions in the limiting cases of small flows and large flows. In view of Proposition 2 and ( 16),

lim s→0 γ flow (s) = - ln(1 -ρ) ρ (1 -ρ) and lim s→∞ γ flow (s) = 1 -ρ. (17)
The mean time-sampled throughput, on the other hand, follows from [START_REF] Litjens | Elastic calls in an integrated services network: the greater the call size variability the better the QoS[END_REF]:

γ time = ρ E[n] = 1 -ρ. (18)
In view of Proposition 4 and ( 16), this is also the mean time-sampled throughput for small flows and large flows. The expression is in fact valid for all flow sizes. This follows from [START_REF] Heyman | A new method for analysing feedback-based protocols with applications to engineering Web traffic over the Internet[END_REF] and the fact that the mean sojourn time of a customer in a processor-sharing queue is proportional to its service requirement [START_REF] Kleinrock | Computer Applications[END_REF].

Figure 1 below compares both throughput metrics for an exponential flow size distribution of unit mean with ρ = 0.8. The mean flow-sampled throughput and the mean time-sampled throughput are estimated using 10,000,000 simulation points (only 10,000 of them are shown on the graph). We observe that, whereas the time-sampled throughput is equal to 0.2 independently of the flow size, as expected, the mean flow-sampled throughput decreases from around 0.4 to 0.2, as predicted by [START_REF] Serfozo | Introduction to Stochastic Networks[END_REF], when the flow size grows from 0 to ∞.

Throughput distribution

There is again no explicit expression for the flow-sampled throughput distribution except in the limiting cases of small flows and large flows, cf. Proposition 2. For the time-sampled throughput distribution, it follows from Proposition 5 that for all x = 1/m, m ≥ 1:

Pr(X time = x) = mπ(m) E[n] = mρ m-1 (1 -ρ) 2 .
This is also the time-sampled throughput distribution for small flows and large flows, cf. Proposition 4. Again, this expression is in fact valid for all flow sizes, as explained in §6.2 below.

The throughput distribution of small flows is shown in Figure 2 for the same scenario as that of Figure 1. Note that the difference between both throughput metrics is significant. In particular, the probability that the throughput is maximum is equal to 1 -ρ (here 0.2) for the flow-sampled throughput and only (1 -ρ) 2 (here 0.04) for the time-sampled throughput.

Sampling bias

We now illustrate the sampling bias associated with the steady-state instantaneous throughput Y defined in §5.1. We consider the model of §6.1 with a proportion p of flows tagged at random. We refer to those flows as class-1 flows and to the other flows as class-2 flows. Denoting by ρ 1 = pρ and ρ 2 = (1 -p)ρ the corresponding traffic intensities, the stationary distribution of the number of flows of each class is given by:

π(n 1 , n 2 ) = n 1 + n 2 n 1 ρ n 1 1 ρ n 2 2 (1 -ρ).
We deduce the mean steady-state instantaneous throughput as evaluated through class 1:

γ 1 = E[ϕ|n 1 ≥ 1] = n 1 ≥1,n 2 π(n 1 , n 2 ) n 1 + n 2 n 1 ≥1,n 2 π(n 1 , n 2 ) = - ln 1 - ρ 1 1 -ρ 2 ρ 1 1 -ρ 2 (1 -ρ).
In view of ( 17) and [START_REF] Stevens | TCP/IP Illustrated[END_REF], this coincides with the mean flow-sampled throughput of small flows for p = 1 and with the mean time-sampled throughput for p → 0.

The bias vanishes when the probability measure is weighted by the number of flows. In view of ( 14), the mean time-sampled class-1 throughput is given by:

γ time 1 = E[n 1 ϕ] E[n 1 ] = n 1 ≥1,n 2 n 1 π(n 1 , n 2 ) n 1 + n 2 n 1 ≥1,n 2 n 1 π(n 1 , n 2 ) = 1 -ρ.
Similarly, for all x = 1/m, m ≥ 1:

Pr(X time 1 = x) = E[n 1 I {ϕ=x} ] E[n 1 ] = n 1 +n 2 =m n 1 π(n 1 , n 2 ) n 1 ≥1,n 2 n 1 π(n 1 , n 2 ) = mρ m-1 (1 -ρ) 2 .
This expression is independent of the choice of class 1. Tagging all flows whose size is in the interval (s -δ, s + δ), for all s > 0 and δ > 0, we conclude by letting δ tend to 0 that the time-sampled throughput of a flow of size s is independent of s.

Rate limits

In many practical cases, flows do not have full access to the link capacity but are imposed a rate limit. In the simple case of a common rate limit c, the instantaneous throughput ϕ of each flow is equal to min(c, 1/n) in the presence of n flows, with n ≥ 1. The stationary distribution becomes:

π(n) = ρ nc π(n -1) if nc ≤ 1, π(n) = ρπ(n -1) otherwise.
In the presence of a permanent flow, we obtain:

π ′ (n) = ρ nc π ′ (n -1) if (n + 1)c ≤ 1, π ′ (n) = n + 1 n ρπ ′ (n -1) otherwise.
We deduce the distributions of ϕ 0 and ϕ ∞ as above. We verify that the timesampled throughput distribution is again the same for all flow sizes. Figures 3 and4 below are the analogues of Figures 1 and2 for a rate limit c = 0.1. Though the difference between both throughput metrics is much lower than without rate limit, it is still significant for short flows. For instance, the probability that short flows get the maximum throughput is equal to 0.59 under flow sampling and 0.49 under time sampling. As a comparison, the probability that the steady-state instantaneous throughput Y is maximum is equal to 0.67.

These results extend to more general models with different access rates and several resources shared according to balanced fairness [START_REF] Bonald | A queueing analysis of maxmin fairness, proportional fairness and balanced fairness[END_REF]. The underlying queuing system is a Kelly-Whittle network. The same properties then hold: the time-sampled throughput distribution is insensitive to the flow size distribution beyond the mean and independent of the considered flow size.

Unfair sharing

Finally, consider a unit capacity link under unfair sharing. Specifically, we consider two flow classes such that the throughput of a class-1 flow is equal to w 1 /w 2 times that of a concurrent class-2 flow, for some weights w 1 , w 2 . The corresponding model is the discriminatory processor-sharing queue [START_REF] Fayolle | Sharing a processor among many classes[END_REF].

There is no simple expression for the steady state distribution, even for an exponential flow size distribution. Figures 5 and6 below are the analogues of Figures 1 and2 for a weight ratio w 1 /w 2 = 10. Note that both the flow-sampled throughput and the time-sampled throughput now depend on the considered flow size. The difference between both throughput metrics is significant.

Conclusion

We have introduced two metrics for throughput performance evaluation, we refer to as flow-and time-sampled throughputs. The former gives the throughput statistics of an arbitrary flow while the latter weights these throughput statistics by the flow durations. We have seen that, under fair sharing assumptions, the time-sampled throughput may be simply interpreted as the steady-state instantaneous throughput weighted by the number of flows. The weighting reflects the fact that congestion events impact the throughput metric in proportion to the number of active flows. In the absence of weighting, the steady-state instantaneous throughput metric is biased and provides only a loose approximation of the flow-sampled throughput.

Fig. 1 .

 1 Fig. 1. Mean flow-and time-sampled throughputs (load ρ = 0.8).

Fig. 2 .

 2 Fig. 2. Flow-and time-sampled throughput distributions of small flows (load ρ = 0.8).

Fig. 3 .Fig. 4 .

 34 Fig. 3. Mean flow-and time-sampled throughputs (load ρ = 0.8, rate limit c = 0.1).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Mean flow-and time-sampled throughputs of low priority flows (load ρ = 0.8).