Sofien Ben Sayadia
email: sofien.bensayadia@esiee.fr

Yaroub Elloumi

Rostom Kachouri

Mohamed Hedi Bedoui

Computational Efficiency Optimization of Optic Disc Detection in Fundus Image

Keywords: fundus image; Optic Disc localization, Parallel algorithms, Graphics Processing unit (GPU), CUDA, Real-time GPU implementation

 proposes a performance optic disk detection approach based on blood vessel tracking and optic disk contrast. However, the method is characterized by a higher execution time which is about 10 s in STARE DB images. M oreover, the execution time increases proportionally with the fundus image resolution. This computational performance is a limiting factor to employ the method in diagnosis systems of ophthalmic diseases. This paper aims to optimize the method processing in order to enhance the computational performance. The first contribution consists of optimizing repet itive steps with the aim of reducing times. Thereafter, all processing steps are implemented in GPU architectures. The experimental results indicat e that each one of the contributions insures enhancing computational performance with speedup equal to 1.7 and 2.5, respectively. The implementation with combined contributions leads to a speedup equal to 8.6 which leads to an execution time about 1 second.

Introduction

The Optic Disk (OD) is a main retinal anatomic structure in fundus image. The OD detection is a critical step in many diagnostic systems for ophthalmic diseases. In the case of diabetic retinopathy (DR), the OD has the same color and contrast than Hard Exudates (HEs). Therefore, several works [START_REF] Kittipolwisaeng | Improved fuzzy C-means clustering in the process of exudates detection using mathematical morphology[END_REF][START_REF] Worapankusakunniran | HardExudates Segmentation based on Learned Initial Seeds and Iterative Graph Cut[END_REF][START_REF] Pavleprentašić | Detection of exudates in fundus photographs using deep neural networks and anatomical landmark detection fusion[END_REF][START_REF] Ittal | A generalized method for the segmentation of exudates from pathological retinal fundus images[END_REF][START_REF] Amin | A M ethod for the Detection and Classification of Diabetic Retinopathy Using Structural Predictors of Bright Lesions[END_REF] proceed to detect and remove OD before segment HEs. Moreover, the works described in [START_REF] Sudeshnasilkar | Detection of neovascularization in retinal images using mutual information maximization[END_REF]18,[START_REF] Yu | M achine Learning Based Automatic Neovascularization Detection on Optic Disc Region[END_REF] to detect the neovascularized blood vessels in the OD in order to deduce the proliferative DR. In the case of the Glaucoma, several works aim to extract the OD and the Cup Disc where the size ratio indicate the presence of the disease [START_REF] Soorya | An automated and robust image processing algorithm for glaucomadiagnosis from fundus images using novel blood vessel tracking and bendpoint detection[END_REF][START_REF] Yin1 | Automated Segmentation of Optic Disc and Optic Cup in Fundus Images for Glaucoma Diagnosis[END_REF][START_REF] Cheng | Superpixel Classification Based Optic Disc and Optic Cup Segment ation for Glaucoma Screening[END_REF]] and the pathology level [24]. Other works detect the Glaucoma by figuring out the position of the cup disc with respect to the OD. The Age-related Macular Degeneration (AMD) is always detected through the drusens segmentation which is located in the macula and has the same contrast than the OD. Therefore, the works described in [START_REF] Khai | Wilson :Automatic fovea location in retinal images using anatomical priors and vessel density[END_REF] aim to detect all shape having a higher contrast, and then eliminate the OD in order to detect the drusens. Moreover, several approaches are proposed in the objective of extracting the macula [START_REF] Geetharamani | M acula Segmentation and Fovea Localization employing Image Processing and Heuristic based Clustering for Automated Retinal Screening[END_REF][START_REF] Prakash | An effective fovea detection and automatic assessment of diabetic maculopathy[END_REF], the fovea [START_REF] Kamble | Localization of Optic Disc and Fovea in Retinal Images using Intensity Based Line Scanning Analysis[END_REF] and the blood vessels [START_REF] Sil Kar | Detection of neovascularization in retinal images using multivariatem-M ediods based classifier[END_REF], which OD location is an important processing step.

A significant number of OD localization methods are proposed in the literature. Recent OD localization approaches offer both high and close detection performance.

Those methods proceed to detect the OD based on their retinal characteristics such as the brightness, the contrast and circular shape. The algorithm described in [START_REF] Pourreza | Computationally efficient optic nerve head detection in retinal fundus images[END_REF] e mploys the Radon Transform (RT) to localize OD based on its brightness and roundness. The RT is applied with several angles to each sub -window in the objective of detecting OD circular shape. Hashim et al. [START_REF] Hashim | Optic Disc Boundary Detection from Digital Fundus Images[END_REF] apply a binary mask on the intensity channel to exclude the background pixels. Then, morphological operators and contrast enhancement techniques (Gamma transformations) are used in conjunction with the difference of the Gaussian filter (DOG) to obtain the OD border. In the work of Giraddi et al. [START_REF] Giraddi | Optic Disc DetectionUsing Geometric Properties and GVF snake[END_REF], a thresholding is employed to eliminate false positive based on the OD brightness and roundness shape. Then, the OD segmentation is performed using the vector field gradient (GVF snake). These methods present high success rates in normal images. However, the detection provides inaccurate result s due to the presence lesions having the same brightness or size than the OD.

Others OD detection methods are based on vessel tracking. Foracchia et al. [START_REF] Foracchia | Detection of optic disc in retinal images by means of a geometrical model of vessel structure[END_REF] used a parametric geometric model (parabolic path) to describe the typical direction of the vessel structure. In the work proposed by Zhang et al. [START_REF] Dongbozhang | Novel Accurate and Fast Optic Disc Detection in Retinal Images With Vessel Distribution and Directional Characteristics[END_REF], the density, compactness and uniformity of blood vessels are formulated to find the OD coordinates. Then, the matched filter is applied in various dimensions in order to provide candidate location.

Those approaches tend to be the most Robust in OD appearance change. However, they can provide a wrong OD detection if the vascular network is partially extracted. Some others OD detection methods employ the OD characteristics, and exploit the location and orientation of vessels. For example, Youssif et al. [START_REF] Xiong | An approach to locate optic disc in retinal images with pathological changes[END_REF] uses the directional pattern of retinal blood vessels for the OD detection. Their method involves normalizing contrast and luminosity. Xiong and Li [START_REF] Youssif | Optic disc detection from normalized digital fundus images by means of a vessels' direction matched filter[END_REF] have proposed a method for locating the OD center by extracting a variety of features including vertical and horizontal vessel intensity and the size of the bright object. Soares et al. [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF] proposes an algorithm based on the cumulative sums of successive subdivisions and the vessel enhancement. The next step consists at following vessel convergence to locate the OD. These approaches tend to be the most effective and reliable, even in incomplete appearance and change of OD, and in incomplete construction of the vascular stru cture. The work proposed in [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF] achieves an OD localization accuracy of 99:15%. This performance is provided using eight public datasets including the STARE and DRIVE ones. However, the method is characterized by a higher execution time which is about 10 s in STARE DB images where the resolution is equal to (700 * 605). Moreover, the execution time increases proportionally with the fundus image resolution. As an e x-ample, the current retinographs TRC-NW 7 SF [START_REF] Zanlonghi | Un comparatif de rétinographes non mydriatiques[END_REF] provides fundus images with res olution equal to (3008 x 2000) which is 14 times greater than a STARE dataset image, and hence a similar rise on execution time. This computational performance is a limiting factor to employ the method in diagnosis systems of ophthalmic diseases. This paper proposes optmize and parallelize the processing method proposed in [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF] with the aim of reducing the time of execution. The article is organized as follows. In Section II, we describe the approach of OD detection. In section III, we analyze processing times and complexity in terms of approximate nu mber of operations. Then, we describe our proposed contributions for the acceleration of the OD detection alg orithm. The evaluation of proposed contributions using different retinal image dat abases is done in Section V, followed by Discussion and conclusio n in last section.

Optic Disc Detection method [7]

This section presents a description of the method proposed by Soares et al [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF] for detecting OD in the fundus image. The main idea is based on identifying the conce ntration and the convergence of the main vessels in order to detect the OD location. This method is composed by successive processing blocks as indicated in the flowchart in Fig. 1 and which are described to the following sessions. The first processing block entitled "preporcessing" starts by resizing the green component of the image to a resolution of (900 * 900) using a bicubic interpolation [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF], in the case where the resolution is greater than (900 * 900). Then, the background is separated from the fundus image using a binary mask applied to the red channel where the threshold is equal to 28 whose result is saved in the BW image. The next step eliminates the noise by applying the "Gaussian Blur" filter with a rectangular structure of 13 × 13 pixels and σ = 4, the result is saved in image I.

Vessel Enhancement

The main objective of this block consists at extracting the main vessels from the ret ina vascular network based on their contrast and thickness. The first step, called "bloodvessel extraction", aims to reconstruct the thicker vesselsthat cross the OD.Thus, the Laplace and gradient filter are applied separately to the Iimage. Next, thevascular network is constructed by calculating the difference between the square absolute value of the gradient image and the Laplacian image.The result is saved in the image X (M * N) as shown in fig. 2 (b).

The provided images are characterized by a higher noisy which avoid distinguish between vascular structures and non-vascular ones. Therefore, the Hessian matrix [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF] is employed in order to and then removing the thin vessels, whose treatments are respectively entitled «edgedetection» and «Local Shape Curvature» («LSC»). The second step, called "edge detection", consists of describing an edge detector by using the second-order derivatives. The first-order derivatives are calculated by applying a "Sobel" filter three times on the image I, using three matrices with deriv ative orders respectively in the x direction, the y direction, and the x and y directions. Then, each "Sobel" filter is followed by a Gaussian filter for calculating the second -order derivatives. The results are recorded respectively in the Y, Z and T images. The third step, called "Local Shape Curvature (LSC)" leads to distinguish vascular structures from non-vascular structures by minimizing the impact of false vessels and avoiding lesions such as micro-aneurysms and exudates. In fact, the background pixels are characterized by a small magnitude of the derivatives (the eigenvalues) relative to the values present in the sets of the vessel pixels. Therefore, the spatial derivative matrix H is constructed for each point p of index i, j such that i∈ [0, N-1], j∈ [0, M -1], as indicated in equation [START_REF] Ben Sayadia | Computational Efficiency of Optic Disk Detection on Fundus Image: A survey[END_REF].

𝐻 = ( 𝑌(𝑖, 𝑗) 𝑇(𝑖, 𝑗) 𝑇(𝑖, 𝑗 ) 𝑍 (𝑖, 𝑗) ) (1)
Then, the eigenvalues λ1 and λ2, are calculated by solving the characteristic equation of the Hessian. These values are essential to distinguish, respectively, the minimum and maximum principal curvatures. The author indicates that the minimal eigenvalue λ1 represents a low contrast corresponding to the pixels belonging to the regions of the blood drops or the microaneurysms.However, the maximum Eigenvalue λ2 represents the points of interest corresponding to the pixels belonging to the main vessels that cross the DO [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF].

For this, at each point p, the minimum eigenvalue is removed from the image of the vascular network X as presented in equation ( 2). The result is saved in RES images (Fig. 3 (b)).

𝑅𝐸𝑆 = 𝑋 -𝑚𝑖𝑛 (𝑋, 𝜆1) (2)
Then, to maximize the impact of the main vessels of the vascular network, each pixel of the RES image is multiplied by the corresponding maximum eigenvalue (Fig. 3 (c)).

Vessel Orientation Detection

This processing blockconsists in extracting the orientations of the main vessels. The retinal vessels merge from the DO vertically and horizontally where directions are between 45 ° and 135 °, respectively. Thus, the vessel structures are extracted in each directionamong α = {0 °; 45 °; 90 °; 135 °}. Since the vessel segments may have differentorientations, vessel in each direction α must be evaluated with different angles. Each angles φα for each direction is presented in equation (3).

𝜑0°= {0°, 15°, 30°, 45°, 135°, 150°, 165°};

𝜑45°= {30°, 40°, 50°, 60°}; 𝜑90°= {45°, 60°, 75°, 90°, 105°, 120°, 135°}; 𝜑135°= {120°, 130°, 140°, 150°};

Since the author aims to enhance linear structures, a logical choice for a structuring element is a "line" with a variable length and a variable angle covering both the short and long vessels. The linear structure (M) is performed for all the lengths and the deviation angles φα, where the line lengths designated by l={5; 10; 15; 20;25} pixels. Thus, the first step in this processing block consists at applying Opening in order to extract the vessels having a shape similar to the elements of linear structure M. The application of this operator leadsv to conserve the vessels corresponding to the element structure (Mki)in each image OKi, where k is the length number and i is the angle number.The k * i = 80 images have the same resolution (M * N) of the image input RES. The second step, each image provided by Opening are compared to of the oneswith deviation φα, by fixing the length l. Then the vessel structure of each orientation {0°; 45°; 90°; 135°} is defined as the sum of the maxima obtained for each value of l. The approach of this step is illustrated by equation ( 4).

𝛼 = ∑ max φ α 𝑶 𝒍 φ α 5 𝑙=1 (4)
The third step leads to extract separately the horizontal and vertical coordinates of the OD position (px, py). For this purpose, two images IH and IV are created which correspond respectively to the horizontal axis and to the vertical axis.The first image IH contains the structure of the vertical vessels, which is determined by subtracting the orientations {0 °; 90}. Similarly, the second (IV) contains the structure of the vertical vessels, which is determined by adding the orientations {45 °; 90 °; 135 °}.

Vessel concentration («VC»)

This processing block consists at detecting the converging points and finding the concentration zone of the vessels network. Thus, the images IH and IV are subdivided successively d times along the vertical and horizontal direction, respectively. The maximu m number of divisions dmax, is calculated such as formulated in equation [START_REF] Amin | A M ethod for the Detection and Classification of Diabetic Retinopathy Using Structural Predictors of Bright Lesions[END_REF].

d 𝑚𝑎𝑥 = round ( max ( N,M ) µ ) (5)
Where µ=70 or 45 respectively for IH and IV images. At each subdivision d, the regions ri (i = 1; ...; d+1), disjoint vertical, are created on the image IH (each resolution (N, M / d + 1)).Then, a vertical average of each region moy_r i is calculated based on the image IH and BW (binary mask image), as indicated in equation ( 6). Then, a vertical division image (dv) is created, whose pixel values are equal to the average value of the corresponding region.Finally, the images extracted at each subdivision are added together.Similarly, this is extended to image IV to create the horizontal division image (dh), as indicated in equation [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF].

𝑚𝑜𝑦_𝑟 𝑖 = ∑ ∑ 𝐼𝐻 𝑖 * 𝑀/𝑑 𝑖 𝑁 1 ∑ ∑ 𝐵𝑊 𝑖 * 𝑀/𝑑 𝑖 𝑁 1 (6) 𝑚𝑜𝑦_𝑟 𝑖 = ∑ ∑ 𝐼𝑉 𝑀 1 𝑖 * 𝑁/𝑑 𝑖 ∑ ∑ 𝐵𝑊 𝑀 1 𝑖 * 𝑁/𝑑 𝑖 (7)

Post-processing

In this processing block, the OD region is determined based on the highest vessel concentration index.Therefore, a horizontal projection on the vertical division image (dv) is performed to identify the pxposition. Similarly, a vertical projection on the horizontal division image (dh) is performed to identify thepyone. The point selected as point [px; py], is marked in the retinal image as the location of OD.

In some retinal images, (px; py) may be slightly outside the OD region. To ove rcome these situations, the maximum point of vessel convergenceand the point of maximum intensity, designated respectively by (cx; cy) and (bx; by), are calculated within these regions. The final OD position is given by the average of the three points p(x,y),c(x,y) and b(x,y) (Fig. 4).

Fig. 4. Final OD localisation

Complexity analysis & parallelism principles

The objective of this section is to study the computing performance of the OD dete ction method proposed by Soares [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF]. Therefore, the computational complexitiesare determined for each processing blocks and for each step of the vessel enhancement and the vessel orientation detection blocks, in terms of the input image resolution M × N. similarly, the execution times are provided by implementing the method in C ++ &OpenCV and run using STARE database images whose resolution is (605 * 700). The computational complexity in terms of approximate number of operations and execution time values are indicated in Table 1.

The post processing complexity is modeled in terms of w which corresponds to the OD diameter. Based on [START_REF] Ben Sayadia | Computational Efficiency of Optic Disk Detection on Fundus Image: A survey[END_REF], the w value can be substituted by 1/7 M, which implies a whole complexity equal to 245*M*N+3.5*M. The implementation leads to an exec ution time equal to 10.25s. This approach is applied to image where the maximal resolution is 900 * 900 pixels. However, actual ophthalmologic devices provide fundus images an important higher resolution, such as the ones described in [START_REF] Zanlonghi | Un comparatif de rétinographes non mydriatiques[END_REF] where the resolution 4 to 7.4 times greater. Based on the whole computational complexity in terms of approximate number of operations, such rising on fundus image resolution implies a similar increase on the execution time. Consequently, the computational performance is a limiting factor to employ the OD method [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF]. Therefore, parallelizing the OD detection proc essing is primordial in order to reduce execution time. Based on the method description in Section II, each step employs the result provided by the previous processing. Consequently, a parallelism on the processing step level is inadequate. All processing steps can be performed in (n × MN) instructions where n is an integer value n > 0. Each step corresponds to an iterative processing with a higher iteration number. Therefore, a parallelism strategy is able to be applied for each step separately. The SIMD principle is the adequate principle of parallelism where the implementation is to run on GPU architectures .

Processing optimization

Our first contribution aims to optimize the processing in order to reduce the execution time. The "LSC" processing consists at defining the H matrix of spatial derivatives for each pixel, as indicated in the equation (1) in section III.A. The Hessian matrix processing requires defining λ parameter by resolving the determinant of the matrix, for each pixel as indicated in equation [START_REF] Yu | M achine Learning Based Automatic Neovascularization Detection on Optic Disc Region[END_REF].

Det(|H -λA| = 0) ( 8
)
Where A is the identity matrix. Thereafter, the formulations of λ1 and λ2 parameters are to be computed. The resolution of the matrix H and hence the formulation of λ1 and λ2 equations are done for each image pixel I whose number is M.N. However, those tasks are performed with the same size (2 * 2) of the H matrix. Therefore, λ1 and λ2 will have the same formulation whatever the pixel is , where their equations are indicated respectively in ( 9) and [START_REF] Yin1 | Automated Segmentation of Optic Disc and Optic Cup in Fundus Images for Glaucoma Diagnosis[END_REF]

Thereby, we proceed to determine the equations ( 8) and ( 9) only once in order to optimize the "LSC" processing time. The Figure 5 illustrates the flowcharts of the "LSC" step respectively before and after optimization, where the formulation steps of λ1 and λ2 modeled with a green background, are removed from the iterative loop. Equations λ1 and λ2 are provided as constant in the implementation. In this way, The "LSC" is performed in O(22MN) times instead of O(123xMN) times, where the input image size is (M*N). The proposed Vesselness measure [START_REF] Aslani | A new supervised retinal vessel segmentation method based on robust hybrid features[END_REF], Hessian multiscale features [START_REF] Thangaraj | Retinal vessel segmentation using neural network[END_REF] and Feature Extraction [START_REF] Zhuc | Retinal vessel segmentation in colour fundus images using Extreme Learning M achine[END_REF] also based in the Hessian matrix eigenvalues . It can be observed that the optimization principle of "LSC" can be in applied order to reduce the exec ution time of these treatments.

Substract matrices

Compute λ1 and λ2

Select min values between λ1 and λ2 delete min eigenvalue multiply max eigenvalue

For the vertical subdivision in the "VC" processing block, the column pixel sums Vv and Bv are performed respectively from the image IH and BW, as indicated in (Eq.6). Thereafter, the successive subdivision is applied. Those steps are repeated 10 times. However, V v and B v computing leads to the same results whatever the iteration is. Similarly, the V h and V bh in the horizontal subdivision are performed 20 times to result identical values.

Therefore, we proceed to optimize the processing by computing once the Vh, Vv, Bv and Bvt vectors. Theirs processing are moved outside the loops, as mod eled with a red background in Fig 6 . In this way, the treatment «VC» is performed in O(MN+10.M+ 20.N) times.

GPU implementation

In this section, we aim to parallelize the implementation of the pipeline image processing proposed by [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF] on GPU architecture. In fact, the Opencv / GPU library proposes a predefined set of image processing functions that are run in GPU architecture, proposing a higher computational performance. Therefore, our parallel implementation principle consists of implementing the steps using directly OpenCV/GPU library if they have corresponding functions. These processing steps are joined in the GPU architectures to avoid communication time b etween host and device. The intermediate data is directly integrated into the memory of the GPUs architecture. In the opposite case, the processing is implemented using CUDA kernel as described in the following section. The implementation of all steps is modeled in fig. 7 where OpenCV/GPU functions and kernels are modeled respectively by yellow and red rectangles.

Vv=Sum cols IH (∑ N 1 )
Bv=Sum cols BW (∑ N 1)

Apply vertical subdivision 10

Vv=Sum cols IH (∑ N 1 )

Bv=Sum cols BW (∑ N 1)

Apply vertical subdivision 10

Vh=Sum Row IV (∑ M 1 )

Bvt=Sum Row BW (∑ M 1)

Apply horizontal subdivision 20

Vh=Sum Row IV (∑ M 1 )

Bvt=Sum Row BW (∑ M 1)
Apply horizontal subdivision 20

"LSC" processing kernel

"LSC"processing consists at computing the RES[i,j] in terms of the pixels X[i,j], Y[i,j], Z[i,j] and T[i,j] where Res in the output image, X, Y, Z and T in the intputs images. The "LSC" treatment can be processed independently for each pixel (N*M values), giving significant computational effort. The result is also that all pixels can be processed simultaneously by independent computing threads. we have previously uploaded the four images X,Y,Z and T, from the host memory to the device global memory and a CUDA kernel is launched to make the LSC treatment.

Vessel concentration

Preprocessing

Firstly we designed an algorithm similar a ''pixel by pixel'' approach where each thread will do the computations concerning one pixel and add the resultat to Image RES. With this approach, The images are divided into several sub-images xi, yi, zi and ti and each is processed into a thread block, as indicated by the yellow cell in the image of Fig. 8 (a). Each thread block (N threads) computes the LSC for a particular sub-image RESi of the result image, where i correspndante to the number of thread block is determined as indicated in equation [START_REF] Kittipolwisaeng | Improved fuzzy C-means clustering in the process of exudates detection using mathematical morphology[END_REF].

NBblock = round(NBpixels/(NBSM * NBgpu )) (11)
Thereafter, we proceed to parallelize the processing of the resulting image pixels. Therefore, each thread in the thread block provides single-pixel "LSC" processing of the resulting image. The pixels of the same indices of images X, Y, Z and T represent the parameters of single-thread Input, as indicated by cells in the input sub-images of Fig. 8 (b). In the second kernel «horiz_subdivision», Vh and Bvt are the input vectors for the subdivision application along the horizontal direction. Similarly, the same parallelism principle is performed for the 20 subdivisions. The main difference consists at modifying the line 6 by TR = Rows / d. The times of parallel processing between threads is decreased from O (20 N) to O (N/4).

6

Experimental Results

Experiment principles

As mentioned earlier, three contributions are used to reduce the execution time. After introducing the soft and hard environment employees described in session 6.2, we conducted three experiments to evaluate the execution time and the speed up of each contribution. The first experiment evaluates the impact of the algorithmic optimization on the OD detection execution time. The second experiment examines the GPU parallel implementation. Thus, we compare the perfo rmance of the Opencv/CPU functions and OpenCV/GPU ones. Hence, the impact of the kernel implementation is studied. Finally, we quantify the rising of the execution time of the whole impleme n-tation after the contributions. To insure a credible experimentation, the implementations are applied using 10 images of the STARE database, chosen randomly.

6.2

Hardware &software resources All implementations are tested on an I7 architecture having a processor frequency equal to 3.67 GHz with 8 GB of main memory and Windows 8.1 running. The parallel version is implemented using the CUDA v8.0 programming environment on NVIDIA Geforce GTX 980. This architecture belongs to the Maxwell family. It co ntains 16 streaming multiprocessors (SM). Each containing 128 processors (GPU) t hat operates at 1216 MHz.

OpenCV (Open Source Computer Vision) is an open source library originally d eveloped by Intel, which provides functions for creating real time applications of co mputer vision and image processing. This library is written in C and C++ and can be run in environments such as Linux, Windows and Mac OS X. Initially, the impleme ntation of the sequential OD detection algorithm was performed by combining the OpenCVversion 3.2 and C ++ programming languages.

In the parallel version, the operations are performed on the GPU using the appropriate OpenCV/GPU extension. OpenCV/GPU are open source libraries that provide an interface for video input, display and programming on GPU using a bunch of highlevel implementations of various image processing and computer vision algorithms [START_REF] Gpucv | A GPU-Accelerated Framework for Image Processing and Computer Vision YannickAllusse 23[END_REF]. The processing time with OpenCV/cuda up to 18 times faster than native OpenCV function [23].

6.3

Algorithmic optimization evaluation In this session, the algorithmic optimizations described in session V are evaluated. Each processing is coded before and after optimization, and run for all image set. The execution times are illustrated in Fig. 10.(a) and Fig. 10.(b) respectively for "LSC" and "VC". In this optimization phase an increased speedup, compared to the normal implementation phase, can be determined. Then, the whole method is implemented before and after algorithmic optimization, where execution times are illustrated in Fig. 11. This optimization allowed achieving averages speed up of 1.7. Fig. 11. Time needed to perform the whole OD detection algorithm without "LSC" & "VC"optimized code Vs. T ime with "LSC" and "VC" optimized code.T he speedup of the OD detection algorithm depending on the "LSC" and "VC" optimized code. The speedup is always given as relative to a CPU runtime.

GPU implementati on evaluation

Evaluati on of Implementing OpenCV functions on GPU architecture.

The OpenCV function set is implemented by default on the CPU and there after implemented in GPU-architecture. The execution times of the Opencv function set, on CPU and GPU are indicated in Fig. 12.(a). Similarly, the whole execution time is shown in Fig. 12.(b) where values prove that running OpenCV functions on GPU architecture leads to a speedup equal to 2.

Evaluati on of GPU kernels

To implement the "LSC" processing on GPU, we proceed to split the input image into 16 sub-images in order to distribute then to 16 blocks, with respect to the SMs number. Thus, each sub-image with (57 * 900) resolution is processed on a separate block. In fact, each SM is composed of Id.x *Id.y =128 GPUs. Thereby, each thread of the same block generates (1 * 900) pixels of the result image. For the distribution of images X, Y, Z and T between the threads, the pseudo code "*App_Hessien" will be called by the host with the following instructions : *App_Hessien<<< 16, 128>>(INPUT:*X,*Y,*Z,*T, cols, rows, OUTPUT:*RES);

The LSC processing is run on both CPU and GPU architectures where execution times are shown in Fig. 13 (a). To implement the "VC" on GPU architecture, we proceed to perform consecutively tow kernels of subdivision image. For the kernel «verical subdivision », the maximum number of subdivisions performed on the vector Vh and Vv is less than the number of GPUs in a single SM of the GTX980 graphics card. Thereby, the processing of each subdivision is assigned to a single SM. In such a way, each thread consists of applying a single vertical subdivision. Similarly, in the second kernel «horiz_subdivision », the number of subdivisions performed horizo ntally does not exceed the number of GPUs in single SM. Thereby; each thread consists of applying a single horizontal subdivision. The blood vessel concentration execution times in both CPU and GPU architectures are illustrated in Fig. 13 (b). The implementation of both "LSC" and "VC"processing on GPU architecture allows to enhance considerably the execution times, where values is illustrated in Fig. 13 (c). The execution time improvement leads to a speedup is equal to 2.5.

Conclusion

The paper objective consists at optimizing the method processing, proposed in [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF], in order to enhance the computational performance. First, we proceeded to optimize processing by shifting steps outs ide loops to reduce the time. Afterwards, all processing are implemented in GPU architectures. Fig. 13. Time needed to perform "LSC" sequential optimized code Vs "LSC" parallel optimized (b) Time needed to perform "VC" sequential optimized code Vs "VC" parallel optimized code.(c) Time needed to perform whole OD detection optimized algorithm with Opencv-cuda and with : the sequential "LSC" & "VC" code Vs. The parallel "LSC" and "VC" code. The speedup is given as relative to a GTX 980 runtime. Based on experimental results, the contributions insure enhancing computational performance where the speedups are respectively equal to 1.7 and 2.5. Then, all c ontributions provide a significant computational efficiency enhancement where the speedup is equal to 8.6, as indicated in Fig. 14. Consequently, the average execution time is reduced from 10.2 to 1.1 seconds. The robustness of the technique proposed in [START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF] is guaranteed by evaluating the met hod in eight publicly-available datasets. Experiments revealed an OD localization accuracy of 99.15%, the proposed model reduces the average computation time 8.6 times when compared to the sequential method

Fig. 1 .

 1 Fig.1. M ethod proposed by Soares et al[START_REF] Soares | Optic Disc Localization in Retinal Images based on Cumulative Sum Fields[END_REF]

Fig. 2 .

 2 Fig. 2. Extract blood vessels networks.(a) Retinal imageI(M *N). (b)Image X(M *N).

Fig. 3 .

 3 Fig. 3. Extract LSC:(a) Retinal image; (b) thin vessels removed;(c) M aximize the impact of principal vessels.

ForFig. 5 .

 5 Fig. 5. Optimization of treatment «LSC» (a) Computingλ1 and λ2 inside loop (b) computing λ1 and λ2 outside loop.

Fig. 6 .

 6 Fig. 6. Optimization of treatment «VC» (a) Computing Vv, Bv, Vh and Bvt inside loop (b) computing Vv, Bv, Vh and Bvt outside loop.

Fig. 7 .

 7 Fig. 7. Diagram block of the proposed GPU implementation.

Fig. 8 .

 8 Fig. 8. Fragmentation of images: (a) Ordering on blocks of threads. (b) Ordering on threads

Fig. 10 .

 10 Fig. 10. Execution time: (a) Without "LSC" optimized code Vs With"LSC" optimized code. (b) Without "VC" optimized code VsWith"VC" optimized code.

Fig. 12 .

 12 Fig. 12. (a)Time needed to perform all OpenCVfuction Vs. Time needed to perform all OpenCV GPU fuction. (b) Time needed to perform the whole optimized OD detection algorithm: With OpenCV Vs. With OpenCV-GPU. The speedup is given as relative to a GTX 980 runtime.

Fig. 14 .

 14 Fig. 14. whole implementation of OD detection: first version (c++/OpenCv version without optimization) vs. last version (parallel/Opencv-Cuda with Optimization).

Table 1 .

 1 OD detection treatments profiling.

	Complexity in terms of approximate num-ber of operations	Execution time (seconds)

  ).

	Algorithm 2 : Vertical_subdivision KERNEL
	_global_void Vertical_subdivision (IN :Float* Vv,
	Float* Bv, ,int cols; OUT :Float *dv)
	tx ← thread x position
	bx ← block x position
	Idx←tx+ bx* block x dimensions
	d← bx+4
	TR ←cols/d;
	for i = (Idx *TR) to ((Idx + 1)*TR) do
	If(i<cols)
	S1 ←S1+Vv[q] %Somme of region in Vv vector
	S2 ←S2+ Bv[q] %Somme of region in Bv vector
	end if
	end for
	S3 ← S1 / S2
	for i = (tx *TR) to ((tx + 1)*TR) do
	if (i<cols)
	dv[i + Idx*cols] ← S3
	end if
	end for
	end KERNEL

«VC» processing kernel

The implementation of "VC" consists of creating subdivision images along the vert ical and horizontal direction. The creations of the subdivision images are performed respectively according to the input images BW and IH or IV.

Based on the VC processing optimization in Section III, we note that the sums of the pixels of the columns of the images IH and BW «Compute ∑ N 1 (IH) & ∑ N 1 (BW) », are performed in parallel using the "cuda :: reduce ()" function of the OpenCV/CUDA library. The results are recorded respectively in the vectors Vv and Bv. The reduce() function can be used to compute horizontal and vertical sums of an image.Thus, the sum of the pixels of the lines of the IV and BW images «Com-pute∑ M 1 (IV) & ∑ M 1 (BW) », is performed using the same function.The results are recorded respectively in the vectors Vh and Bvt.

Therefore, the parallelism of the "VC" processing consists in using two consecutive kernels «Vertical_subdivision » and « horizontal_subdivision».

Based on section II.C (Eq (5)), the maximum division number dmax is computed based on the largest dimension of the input image. Thus the images are proportionally scaled in a way that the largest dimension is 900 pixels. Consequently, in the "verti-cal_subdivision" processing allows only 10 subdivisions. Similarly, in the "horizon-tal_subdivision" processing allows 20 subdivisions.

The implementation of the kernels on GPU involves the transfer of four vectors Vh, Vv, Bv and Bvt into the global memory of the GPU. Vv and Bv are the input vectors for the «Vertical_subdivision » kernel. Firstly, the parallelism principle leads to perform the 10 subdivisions, which are applied along the vertical direction in parallel and each is processed into a thread block. The number of thread blocks is equal to the number of subdivision. Each thread block applies a single subdivision. Thereafter, at each vertical subdivision, we proceed to parallelize the processing of the creation of rp (p = 0; ...; d) disjunct regions. Each thread in the thread block is responsible for creating a single region by dividing the sum of vector region Vv by the sum of region Bv. The size of each region is calculated by dividing