Frédéric Monge
email: frederic.monge@univ-rennes1.fr

Ph.D Dzhoshkun I Shakir

Ph.D Florence Le Jeune

M.D. • Xavier Morandi

M.D. • Nassir Navab

Ph.D • Pierre Jannin

Acquisition Models in Intraoperative Positron Surface Imaging

Keywords: freehand beta probe, functional imaging, intraoperative imaging, neuronavigation, radiation detection physics, reconstruction, validation assessment

Purpose: Intraoperative imaging aims at identifying residual tumor during surgery. Positron Surface Imaging (PSI) is one of the solutions to help surgeons in a better detection of resection margins of brain tumor, leading to an improved patient outcome. This system relies on a tracked freehand beta probe, using 18 F based radiotracer. Some acquisition models have been proposed in the literature in order to enhance image quality, but no comparative validation study has been performed for PSI. Methods: In this study, we investigated the performance of different acquisition models by considering validation criteria and normalized metrics. We proposed a reference-based validation framework to perform the comparative study between acquisition models and a basic method. We estimated the performance of several acquisition models in the light of four validation criteria: efficiency, computational speed, spatial accuracy, and tumor contrast. Results: Selected acquisition models outperformed the basic method, albeit with the real-time aspect compromised. One acquisition model yielded the best performance among all according to the validation criteria: efficiency (1-Spe: 0.1, Se: 0.94), spatial accuracy (max Dice: 0.77) and tumor contrast (max T/B: 5.2). We also found out that above a minimum threshold value of the sampling rate, the reconstruction quality does not vary significantly.

Introduction

According to the World Health Organization (WHO), gliomas are the most frequent primitive brain tumor, affecting 5 persons out of a 100,000 every year. The main treatment option is surgical resection with a careful consideration of the brain eloquent areas while achieving maximal resection, resulting in better outcome for the patient [START_REF] Ius | Low-grade glioma surgery in eloquent areas: volumetric analysis of extent of resection and its impact on overall survival. a single-institution experience in 190 patients: clinical article[END_REF]. Survival prognosis from low grade glioma (LGG) is better when the maximum tumoral tissue has been removed, slowing down evolution into high grade glioma (HGG), which is much more aggressive. The extent of resection [START_REF] Mcgirt | Extent of surgical resection is independently associated with survival in patients with hemispheric infiltrating low-grade gliomas[END_REF], correlated with the success of surgery, is assessed over postoperative MRI of the patient, to eventually detect recurrence. In the last decade, several technologies have been developed to improve the completion of resection by allowing detection of residual tumor during surgery mainly for HGG. Intra-operative MRI yields good results towards maximum resection [START_REF] Hata | Intraoperative tumor segmentation and volume measurement in mriguided glioma surgery for tumor resection rate control[END_REF], better than conventional neuronavigation-guided glioblastoma (HGG) surgery [START_REF] Kubben | Intraoperative mri-guided resection of glioblastoma multiforme: a systematic review[END_REF]. Indeed, many limitations exist by using magnetic field devices in the operating room (OR) such as cost issues. 5-aminolevulinic acid (5-ALA) fluorescence image guided surgery is a relevant alternative to improve surgical resection outcome by giving a real-time visual histology and tumor boundaries [START_REF] Behbahaninia | Intraoperative fluorescent imaging of intracranial tumors: A review[END_REF]. Good results have been achieved in HGG [START_REF] Roessler | Intraoperative tissue fluorescence using 5-aminolevolinic acid (5-ala) is more sensitive than contrast mri or amino acid positron emission tomography (18f-fet pet) in glioblastoma surgery[END_REF]. Recently, intraoperative ultrasound (iUS) systems, performing LGG delineation, have been used to help the surgeon [START_REF] Selbekk | Comparison of contrast in brightness mode and strain ultrasonography of glial brain tumours[END_REF][START_REF] Šteňo | Navigated threedimensional intraoperative ultrasound-guided awake resection of low-grade glioma partially infiltrating optic radiation[END_REF] during surgery. Poor quality of images makes the interpretation subjective without training, especially for small or superficial residuals [START_REF] Gerganov | Twodimensional high-end ultrasound imaging compared to intraoperative mri during resection of low-grade gliomas[END_REF]. In oncology, Positron Emission Tomography (PET) is a gold standard for tumor detection and grading. Considering 18 F-FET (fluoro-ethyl-tyrosine) radiotracer, uptakes in brain tissues is better for detecting LGG by comparison with usual 18 F-FDG (fluoro-desoxy-glucose) [START_REF] Pauleit | Comparison of 18f-fet and 18f-fdg pet in brain tumors[END_REF]. The definition of tumor margins is more accurate with FET than FDG, according to the enhanced tumor contrast. Despite higher sensitivity and specificity in glioma diagnosis compared to MRI or CT scan [5], 18 F-PET scan suffers from a low spatial resolution (between 4 and 7 mm) and is rarely used in the intraoperative context. Moreover, since the gamma radiation (photons) has higher penetration ability than beta particles (positrons), if a gamma sensor is used, radiations emitted from surface or deeper tissues cannot be distinguished. In our case, beta particles is preferable because of its lower penetrating ability. Therefore, all the detected positrons are certainly emitted from the surface tissue (not deeper than 2 mm below). On the contrary, freehand nuclear beta probes, which are pen-sized devices that can detect positrons emitted by radiolabeled tissues, are more suitable in the OR for the detection of residual tumor over the resection cavity. In addition, a surgery performed with such device Fig. 1: (a) Beta probe and its tracking target. will be proceeded according to a normalized protocol considering radioprotection equipment for the surgical team, in agreement with the Nuclear Safety Agency. In this context, such guided surgery will be performed for LGG only, representing less than 10 tumors per year. Freehand Positron Surface Imaging systems [START_REF] Navab | Navigated nuclear probes for intra-operative functional imaging[END_REF] have been introduced to generate functional images using a beta probe combined with a tracking system, as a neuronavigation workstation. Such a system is used to estimate the activity distribution of the radiotracer on 3D iso-surface images (mesh) to check resection margins. However, to perform good quality images from highly undersampled data, radiation detection physics must be modeled. In previous work, acquisition models have been proposed to address this problem, generating reconstruction images with improved spatial localization and tumor contrast. One of the previous studies provided a comparison between different acquisition models [START_REF] Hartl | Detection models for freehand spect reconstruction[END_REF], pointing to the dramatic variations of reconstructions in SPECT (Single Photon Emission Computed Tomography). In Positron Surface Imaging (PSI), this variation between acquisition (or detection) models can lead to bias in the detection of residual tumor during surgery. In this study, we propose a validation framework for comparing the performance of acquisition models. Section 2 presents the methodology of freehand PSI including activity reconstruction using detection models. We also introduce our validation framework and the acquisition models studied in this paper. Section 3 presents the results. Finally, section 4 discusses the implications of our results and concludes with an outlook on future work.

Materials and methods

Freehand PSI System

The proposed system consists of a freehand tracked beta probe [START_REF] Daghighian | Intraoperative beta probe: a device for detecting tissue labeled with positron or electron emitting isotopes during surgery[END_REF]. After a scan of a radiolabeled surface of interest (SOI) with our system (see figure 2), each accumulated probe reading is synchronized with its mean corresponding pose in the reference coordinate system, giving one measurement during a time window (6D + 1D data). A sampling of the SOI is performed by generating a surface mesh with isotropic suxels 1 , based on the probe positions and orientations. Thus, Fig. 2: A Positron Surface Imaging system: an optical tracking system (blue arrow), a beta probe (red arrow) and its workstation (green arrow). measurements are used to determine an estimate of the surface activity distribution of the SOI (i.e. retrieve the activity value in each suxel of the surface: Positron Surface Imaging).

Surface Activity Reconstruction

Radiation detection models (or acquisition models) that allow for iterative reconstruction have been introduced and shown to substantially improve image quality in terms of spatial resolution and tumor-to-background (T/B) signal contrast [START_REF] Shakir | Acquisition model for iterative reconstruction of navigated beta-probe surface images[END_REF][START_REF] Shakir | Evaluation of an ad hoc model of detection physics for navigated beta-probe surface imaging[END_REF]. By taking conventions introduced in [START_REF] Hartl | Detection models for freehand spect reconstruction[END_REF], the sampled SOI suxels x = (x j ) and the measurement vector m = (m i ) are considered here. Each stored measurement m i can be written as a linear combination of the (unknown) activity values of the suxels x j within its intersection with the FOV of the probe as follows:

m i = N j=1 a ij .x j (1)
where coefficients a ij are computed through an acquisition model, corresponding to the probability that the sensor at pose i has detected positron emissions from suxel x j . The unknown reconstruction vector x (i.e. the surface activity reconstructed, composed of N suxels) may be found by solving the system m = A.x, with A = (a ij ) the system matrix (size M * N with M the length of the measurement vector m). Solutions of this system may be approximated using the Maximum Likelihood Expectation Maximization (MLEM) algorithm introduced by [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF], to compute the pseudo inverse of the system matrix A.

Reference-based assessment of reconstruction

In this section, we present our validation scheme for reconstructions generated from Positron Surface Imaging System. The need for a dedicated assessment scheme in this novel imaging approach drives us to develop a specific framework, according to the methodology proposed in [START_REF] Jannin | Assessment of image-guided interventions[END_REF]. Our focus is to propose a comparative study allowing the performance characterization of acquisition models and the identification of parameters which influence reconstruction outputs, for intraoperative use in neurosurgery. We focused on assessment levels 1 and 2, dedicated to technical parameters and clinical reliability setting respectively. Biological phantom datasets acquired in simulated OR conditions with a minimal setup are considered, where the clinical realism is close to a real surgery with very few controlled phantom parameters. Based on this methodology, assessed features concern the reconstruction efficiency, computation speed and tumoral characterization accuracy. Our validation protocol is based on [START_REF] Jannin | Model for defining and reporting reference-based validation protocols in medical image processing[END_REF], where a scheme for medical image processing validation procedure is described via a structured framework. It relies on the comparison between the results of a method and a reference, assumed to be close or equal to the correct solution. This scheme has been successfully used for designing and reporting medical image processing validation studies, such as in [START_REF] Maier-Hein | Comparative validation of single-shot optical techniques for laparoscopic 3-d surface reconstruction[END_REF]. The main components of this reference-based validation scheme are presented in figure 3.

The validation objective of our study was to evaluate detection performances and clinical relevance (V C ) of reconstructed images generated with acquisition model F M , on in-vitro data sets in the context of neurosurgery, for the detection of resection margins (C). Following sections provide an instantiation of the framework, where each component is detailed in figure 4.

D I M P Ii M D I ref F Mi F ref R Mi E Mi R ref E ref F NRref R Nref E Nref F NRMi R NMi E NMi F C1 F C2 O D1i O D2i
Fig. 4: Adapted reference-based validation process according to our evaluation context.

Validation data sets D I

Our validation datasets consist of acquisitions with an in-vitro phantom as described in [START_REF] Shakir | Acquisition model for iterative reconstruction of navigated beta-probe surface images[END_REF].

Bioluminescence images D ref I

Cancer cells were supplied with luciferine (0.1 ml per dish) to perform imaging with a photographic camera of each Petri dish. We define D ref I the obtained 2D images considered as the reference to which reconstructions RM i (see below) will be compared for each configuration (shown in figure 5, lower row). Petri dish 1 (figure 5d) : two high-contrast tumors. Petri dish 2 and 3 (figures 5e and 5f) : three tumors within regular contrast.

Raw positron surface data D M I

Petri dish surfaces were scanned with a 3axis step motor setup in 1 mm steps. For each position of the step motor, probe readings were accumulated for 3 seconds. Due to the acquisition frequency of the beta probe (≈ 10 Hz), each position is associated with radioactive counts in the time window (i.e. ≈ 30 probe readings during each 3 second acquisition). A long acquisition time leads to good statistical data consistency in laboratory conditions, but this might not be feasible in intraoperative conditions. In order to observe the impact of lower acquisition times on data quality, raw data were undersampled by a selected factor (see input parameters section 2.3.5). Below, D M I is used to represent each configuration of multi-sampled raw positron surface data.

(a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l)

Surface activity reconstruction method F M i

To perform reconstruction, raw data D M I were sampled into a 2D rectangular grid (50 x 50, 1 mm isotropic pixels), for a fixed scanning depth value. Also, the measurement vector m was created from D M I . Here, objective of reconstruction is to define each pixel value of the image, within the rectangular 2D grid. When the basic approach activity projection (proj) is considered, pixel values are defined using a nearest neighbor interpolation over measurement vector m for each position in the grid. By contrast, when acquisition models (presented below) are used to estimate matrix A, reconstruction is obtained through equation x = A -1 .m, with x the reconstruction vector, m the measurement vector and A -1 the pseudo inverse of the system matrix. The MLEM algorithm with a selected number of iterations (see section 2.3.5) was used to compute its pseudo inverse. RM i represents the image from the reconstruction vector x estimated with method F M i .

Solid angle model (SA)

As baseline to define each contribution of suxel x j to measurement m i , the solid angle model (SA) introduced by [START_REF] Oezguer | Phenomenological models for intraoperative positron emission surface imaging using handheld probes[END_REF] considers the geometrical attenuation based on the relative position and orientation of the probe to the source. An implementation of the SA model, used in [START_REF] Shakir | Acquisition model for iterative reconstruction of navigated beta-probe surface images[END_REF] and defined in equation 2, computes each attenuation coefficient a ij , using the following parameters: d ij the distance between the detector and the suxel x j for measurement m i (probe to source distance), r the detector radius and α ij the tilt angle between the detector axis and the detector to source axis (see figure 6).

a ij = f (d ij , α ij ) = 1 2 .   1 - 1 (r 2 /d 2 ij ) + 1   . cos α ij (2)
Look-up-table model (LUT) In the previous analytical model, each coefficient of the system matrix A is directly computed from the measurement vector m i . An alternative solution consists in selecting precomputed a ij coefficients stored in a look-up table (LUT). As proposed in [START_REF] Shakir | Two new ad-hoc models of detection physics and their evaluation for navigated beta probe surface imaging[END_REF], attenuation coefficients a ij are defined by equation 3 for respective LUT positions in the FOV.

a ij = L(i, j) (3)
Here, L refers to the LUT stored in memory and the position in the FOV is represented by the pair (i, j), shown in figure 7a. Implicitly, a mapping function is used to relate the source-to-detector distance d ij and the tilt angle α ij to the LUT coefficients. This model is based on coefficients stored in memory. Partition model (PM) Although the LUT model yields good results in defining structures better than the SA model when comparing with a groundtruth [START_REF] Shakir | Two new ad-hoc models of detection physics and their evaluation for navigated beta probe surface imaging[END_REF], looking up the a ij in the LUT for each (i, j) pairs involves additional computational time for defining the whole system matrix A; possibly using nine neighbor averaging interpolation. A partition model (PM) is proposed to overcome this constraint, based on the division of the FOV into sectors where particular homogeneous coefficients a k are assigned for each partition k [START_REF] Monge | Partition-based acquisition model for speed up navigated beta-probe surface imaging[END_REF]. PM model is analytically defined in equation 4 and example of a partitioning scheme is shown in figure 7b. Here, a k factors are defined through a k-means algorithm performed over the LUT, where the number of cluster k corresponds to the wanted number of partition k. Thereby, values of the LUT defined within each cluster k are averaged to obtain a k attenuation coefficients. In the same way as the LUT model, a k coefficients are assigned according to respective (i, j) positions in the FOV.

a ij = a k source within partition k, 0 source not within any partition.

(4)

Error ÊM i of reconstructions

The reconstruction process involves estimating the pseudo-inverse of the matrix A (computed via acquisition model F M i ) by the MLEM iterative method. Thus, reconstruction error ÊM i , correlated with iteration number, was defined as the norm of residual given ÊM = ‖A.ĉ -g‖with ĉ the final estimated reconstruction vector, A the coefficient matrix and g the probe reading vector.

Reference method F ref and associated error Êref

In oncology, presence of residual tumor is generally observed with the radiotracer uptake on functional PET images, combined with a CT scan or an MRI when available (high performance) for better localization using co-registered images (CT in this case is primarily needed for attenuation correction). However, this modality requires long acquisition time (30 min. for a PET scan). In this study, we used bioluminescence imaging as a reference for determining the tumor location in Petri dishes [START_REF] O'neill | Bioluminescent imaging: a critical tool in pre-clinical oncology research[END_REF], without affecting initial cell growth behavior [START_REF] Tiffen | Luciferase expression and bioluminescence does not affect tumor cell growth in vitro or in vivo[END_REF]. This method gives high resolution functional images, where accuracy is assumed better than reconstruction with methods F M i . Reference images were obtained by cropping bioluminescence images D ref I following a square (size 180 x 180 pixels), fitted and tangential to the Petri dish edges. The symbol RRref is used to represent reference bioluminescence images for each configuration. In addition, the associated error Êref is related to the loss of background signal in bioluminescence image (absence of peri-tumoral contrast, see figure 5), which cannot be quantified.

Input parameters P I

In this section, parameters that can possibly influence the performance of reconstruction are elaborated. According to the model of reference-based validation [START_REF] Jannin | Model for defining and reporting reference-based validation protocols in medical image processing[END_REF], parameters related to data acquisition and the reconstruction method are considered.

Sampling raw data D M I

The variation of this data-related parameter allows the observation of data quality over different simulated acquisition times. The range [0.1; 3] seconds has been selected for simulating scans between a quick scan (0.1 sec per position) and a more detailed one (3 sec per position), corresponding to 3% and 100% of probe readings considered per position, respectively. In figure 5, row 3 and row 4 show reconstructed images generated with the projection method using 3 sec. and 1 sec. acquisitions, respectively.

MLEM iteration number

According to [START_REF] Hartl | Detection models for freehand spect reconstruction[END_REF], reconstruction performed with large number of iteration (e.g. 100) increases artifact formation, in SPECT imaging systems. The variation of the number iteration in the range [1; 50] allows us to observe the influence of this method-related parameter on reconstruction.

Number of clusters in PM Finally, in the partition model PM proposed to reduced computational power by dividing the LUT FOV, according to a k-means algorithm. The variation of the number of clusters k in the range [2; 6], allows us to observe influence of this method-related parameter on the computation time in comparison with other reconstruction methods. In the following, we use the symbol PM-k (PM-2, PM-3, ..., PM-6).

Normalization function of reconstruction F N RM i

In order to compare reconstructions with their respective reference RRef , RM i images are resized to fit reference image spatially (180 x 180 pixels) and amplitudes are normalized within the range [0; 1]. Thus RNM i represents RM i images normalized through F N RM i function and ÊNM i represents visual blurring errors due to the interpolation method used in the resizing step.

Normalization function of reference F N Rref

Due to the different modalities used for generating the validation datasets D ref I and D M I , a scale shift exists between bioluminescence reference images RRef and reconstructions RM i , respectively acquired with a camera and the PSI system. For each Petri dish configuration, the reconstructed image RM i generated through the projection method F M i , parametrized with the minimum sampling rate (3 sec per position, no undersampling) is considered. The visual scaling shift was evaluated by overlapping reference Rref with its respective considered reconstruction. Experimentally, a scale of 1.1 is observed for reconstructed image by comparison with bioluminescence images, in each case (see figure 8, upper row). Then, a 2D affine spatial transformation (with 1.1 scaling parameter) was applied to cropped bioluminescence images (figure 8 lower row). A binarization step was performed by thresholding images with a determined level. Heuristically, a 0.3 level value gives the global size-position (see figure 9). Finally, remaining Petri dish edges were removed with a morphological opening operation. RNref represents segmented bioluminescence images in the SOI coordinate system obtained with both level values, and ÊNref represents the qualitative error corresponding to the possible mismatch between Petri dish edges (shown in figure 8, lower row).

Validation criterion and comparison function F C

From normalized outputs RNref and RNM i , we estimated the following validation criterion: reconstruction efficiency, computation speed, spatial accuracy and tumor characteristics (contrast and size). Efficiency For each reconstructed image, compared to reference image, we calculated the number of false negative (FN), false positive (FP), true positive (TP) and true negative (TN) pixels. Then, sensitivity-specificity curves were computed using these metrics. Sensitivity has been defined as the true positive rate (T P/(T P + F N )), which represents the ability of our system to detect tumor when a tumor is present in the reference. Specificity has been defined as the true negative rate (T N/(T N + F P )), which represents the ability of our system to detect a tumor when a tumor is not present in the reference. Efficiency of reconstruction O D = F C ( RNref , RNM i ) can be quantified with the optimal point observed on sensitivity-specificity curves (i.e. optimal point obtained by maximizing sensitivity and specificity).

Computation speed

The computation time of each reconstruction performed with acquisition model was determined. The speed performance is a relevant criterion for an efficient intraoperative use.

Spatial accuracy

Correlation of reconstructions RNM i to respective reference RNref were estimated. Due to the nature of RNM i and RNref outputs, the Dice coefficient (DSC) was used as comparison function F C , shown in equation 5.

DSC(A, B) = 2 | A ∩ B | | A | + | B | (5)
Tumor contrast Similarly with PET scan images, the functional information over reconstructions is provided by the tumor to background ratio (T/B). From normalized reconstructed images RNM i , mean T/B was computed over normalized reconstruction RNM i using normalized reference RNref as binary mask. Based on the boolean values in the mask definition, each pixel with a 'true' value was considered part of the tumor region of interest (tROI) and each 'false' value part of the background region of interest (bROI). Each ROI was averaged by its number of pixels and thus T/B ratio was computed, defined as result O D .

Quality index O QI and hypothesis testing F H

From comparison functions F C defined above according to the validation criteria, a quality index can be deduced for characterizing the overall performance of acquisition models. In each sensitivity-specificity curve, Area Under the Curve (AUC) was computed for each acquisition model using the trapezoidal numerical integration method (F QI ), representing the quality index O QI . In addition, no statistical evaluation F H was performed in this scheme, due to insufficient amount of data.

Visual comparisons have been done between reconstruction and related reference.

Results

Two independent studies were performed using comparison function F C along two input parameters P I . The first comparison was made by varying the sampling rate while fixing the number of MLEM iteration at 25, an empirical value used in previous work on PSI systems [START_REF] Shakir | Evaluation of an ad hoc model of detection physics for navigated beta-probe surface imaging[END_REF][START_REF] Shakir | Two new ad-hoc models of detection physics and their evaluation for navigated beta probe surface imaging[END_REF]. The second one was made by varying configuration. Visually, performance of acquisition models can be observed over sensitivity-specificity curves derived from reconstructions. In figure 12, mean ROC (Receiver Operating Characteristic) curves over Petri dish configurations were defined by considering optimal sampling and MLEM input parameters P I (MLEM = 25 iterations and sampling = 3 seconds). In ROC representation, the best curve is found by minimizing the 1-specificity term while maximizing the sensitivity.

Values of optimal points are presented in table 1. Here, two acquisition models outperformed the projection method: SA and LUT models.

In figures 13a and 13b, AUC values are shown for ROC curves varying the sampling rate and MLEM iterations respectively. We observed that the projection method results are nearly constant along the sampling rate (AUC around 0.74) and invariant along the number of MLEM iterations. From 1-sec sampling onwards, each acquisition model tends to 0.77 mean AUC value, except for method PM-2 (mean AUC = 0.67). Moreover, mean AUC values tend to decrease linearly along MLEM iterations, with a steeper decrease when using the PM-2 model: from a 0.86 starting AUC value for one MLEM iteration down to 0.79 around for 50 MLEM iterations. Globally, both AUC curves emphasize three efficient and stable acquisition models by comparison with the projection method: SA, LUT and PM-5 models. Considering computation time required by the projection method (3 ms independently of sampling rate), the reconstruction process was hardly time consuming. This computation time of reconstruction methods is linearly dependent on the sampling rate of data, as shown in figure 14. The LUT was the most time consuming method with a maximum of 31 minutes observed. Other methods reached a maximum around 28 minutes of computation required to perform reconstructions. The quantification of spatial accuracy of each reconstruction method is given by the Dice score metric. In figure 15a, the projection method yields the most stable Dice score, around 0.6, overtaking the PM-2 model. Starting with the 0.5 sec sampling value, all other reconstruction models outperformed the projection method, where higher Dice values are obtained with the LUT model. The mean Dice score variation along MLEM number of iterations is shown in figure 15b. Results show a slight variation of the behavior of the Dice value for each reconstruction method starting with 5 MLEM iterations. Two methods are emphasized here: LUT and SA models.

Finally, variations of the tumor contrast along the sampling rate and number of MLEM iterations are shown in figure 16a and 16b respectively. The projection method yields constant mean T/B contrast rate over varying sampling rates or number of MLEM iterations. For reconstruction methods, T/B values outperformed the projection method, increasing slightly with the sampling, starting at 0.5 (figure 16a). On the contrary, mean T/B value decreases with increasing MLEM number of iterations (figure 16b), but remains better than the projection method. LUT and SA models yields the highest T/B values for a sampling rate of 3 combined with 5 MLEM iterations.

Discussion and conclusion

Discussion

We have compared detection models by assessing respective reconstructed images from a chosen dataset, relying on a validation framework. We have outlined benefits and drawbacks of each acquisition model according to selected input parameters P I (i.e. section 2.3.5) and validation criteria (i.e. section 2.3.8). Images reconstructed using acquisition models yield global detection improvement along chosen comparison functions and quality index (i.e. section 2.3.9) with respect to the basic projection method. According to the clinical applicability involved here in the selection of an adequate acquisition model, the trade-off between reconstruction quality and reconstruction time has been outlined and influences this selection. The quality of reconstruction has impact on decision making of the surgeon on the detection of residual tumor. Reconstruction time is important for respect of the surgical workflow. Additional time for PSI acquisition needs be as short as possible in order to avoid potential complications or infections [START_REF] Daley | How slow is too slow? correlation of operative time to complications: an analysis from the tennessee surgical quality collaborative[END_REF]. Our previous work [START_REF] Shakir | Acquisition model for iterative reconstruction of navigated beta-probe surface images[END_REF], involved long acquisition times. Although long acquisitions provide statistically good data, it is not realistic to expect such data in clinical conditions. In addition, acquisition time and scan quality depend on the operator (i.e. neurosurgeon). In this work, we simulated the acquisition time variations based on the sampling rate input parameter. Results show that this variation slightly affects data quality. However, transient states can be observed in figures 13a, 15a and 16a in the range of [0.1; 0.7] seconds before reaching a stable value. This, suggests a threshold value of 0.7 second per location that should be considered as a measure of the data acquisition speed. A sampling rate above this value does not improve the data quality significantly, but still increases acquisition times. Empirically, 25 iterations for the MLEM reconstruction algorithm were fixed to perform reconstructions with freehand PSI systems [START_REF] Shakir | Evaluation of an ad hoc model of detection physics for navigated beta-probe surface imaging[END_REF][START_REF] Shakir | Two new ad-hoc models of detection physics and their evaluation for navigated beta probe surface imaging[END_REF], close to values used in freehand SPECT systems [START_REF] Matthies | Mini gamma cameras for intra-operative nuclear tomographic reconstruction[END_REF][START_REF] Wendler | Towards intra-operative 3d nuclear imaging: reconstruction of 3d radioactive distributions using tracked gamma probes[END_REF] (i.e. 20 iterations). We investigated this value by varying it as an input parameter and looking at the respectively generated reconstructions. Increasing the MLEM number of iterations has slightly digressive effects on validation criteria. Peaks were observed at five iterations in figures 13b, 15b and 16b. Based on our experiments, we propose to use this value in future experiments for optimal validation criteria values. Intraoperative use involves high performance in terms of computational speed. One factor that affects it is the sampling rate, which is directly correlated with the amount of data acquired per location. In figure 14, the linearity between the computational speed and sampling rate can be observed. Another factor that affects computational speed is the choice of the acquisition model. In this context, the partition model was introduced as a faster LUT model. In this work, we have proposed using a range of k clusters in the validation of the partition model ( where similar reconstruction results as the off-the-shelf LUT have been expected). In figure 14, each generated PM-k model falls under the LUT experimental curve, i.e. is better in terms of computational speed. Results given by each PM-k were slightly faster than other methods, especially with the PM-2 method.

According to our experiments, one acquisition model outperformed others: the LUT model. This model offers the best result in terms of efficiency (1-Spe: 0.1; Se: 0.94), spatial accuracy (max Dice: 0.77) and tumor contrast (max T/B: 5.2), but requires the highest computational speed among all models (max computation time: 31 min). However, considering the stability of data quality over the sampling time threshold (i.e. 0.7 second per location), the LUT model seems appropriate for intraoperative use. Moreover, the current LUT model developed in [START_REF] Shakir | Two new ad-hoc models of detection physics and their evaluation for navigated beta probe surface imaging[END_REF] has been generated from planar acquisition. In real cases, the SOI would represent the resection cavity as a 3D surface. A realistic model would consider the inclination of the probe to the surface axis, the tilt angle. In our validation framework, qualitative errors have been identified in the output of the reference method F Ref and the normalization functions F N RM i and F N Rref . Such errors inducing possible bias in the results of our comparison functions F C should be taken into account. According to [START_REF] Jannin | Model for defining and reporting reference-based validation protocols in medical image processing[END_REF], bias can be studied with respect to two criteria: accuracy and metric realism. Firstly, ÊRef has been identified as the loss of background signal in bioluminescence images. No difference can be observed between inner and outer contrast of the Petri Dish (background in blue, figure 5, row 2) by comparison with reconstructed images (outer background in red and inner background in blue, figure 5, row 3 and 4). Thus, bioluminescence images suffer from a low realism due to radiotracer uptake in both tumoral and healthy tissue. In terms of accuracy, this can induce a few errors in the ROC curve creation for the efficiency validation criteria, without impacting the global behavior of curves. Secondly, ÊNM i corresponds to the blurring effect of the interpolation performed on reconstructions. The smaller the size of the reconstructed images, the more the blurring effect increases, leading to a high bias in comparison functions. A possible impact of this is the overestimation of residual tumor size. And thirdly, ÊNRef gives the overall idea of the possible misalignment after the normalization method of references. In a real case, the imaging devices used for the different modalities that give the actual acquisition and the groundtruth reference, do not have the same coordinate system. Here, since edges of Petri dish are visible in both images (reconstruction and bioluminescence), it is possible to perform circle fitting to estimate the difference between the radius of the actual and reference images, which could subsequently be used as an estimate of the misalignment error ÊNRef . While a Petri dish radius larger in reference than in reconstruction will lead to an underestimation of the real tumor size, a larger radius in reconstruction than in reference will lead to an overestimation. We have introduced several validation criteria for selecting an optimal acquisition model. They are strongly dependent to the nature of the reference dataset. The Dice score metric is adapted for validating the global localization of tumor spots with a reference extracted for an imaging modality of a different nature than the one used (i.e. bioluminescence imaging). In contrast, for a reference from an emission-based imaging modality (PET, SPECT, etc..), the normalized cross correlation metric (NCC) would give more information in checking resection margins. In addition, the tumor size over reconstruction (in cm 2 ) would provide information on the under/over estimation of the real tumor size in reference.

Conclusion

For the first time, we have compared acquisition models, dedicated to Positron Surface Imaging based on a dedicated reference-based validation scheme. Associated to our input parameters, new validation criteria can be integrated, considering the nature of reference datasets. Combined with AUC values, quality criteria such as Youden index can also be used for the characterization of ROC curves in terms of detection. In future work, the validation protocol will be extended to 3D realistic reference datasets (e.g. realistic neurosurgical phantom), as a necessary validation step involving realistic conditions for clinical use in neurosurgery.

Fig. 3 :

 3 Fig. 3: Main steps of the reference-based validation protocol for medical image processing (figure reproduced from [10]): from a defined validation objective related to a clinical context C, results of the method to be assessed F M applied on validation data sets D I with parameters of the method P I are compared with a reference through a comparison function F C . Then, the final result is obtained by testing the output of F C against the validation objective.

Fig. 5 :

 5 Fig. 5: Row 1: Brightfield images of the three Petri dish configurations, where right oriented red arrows show BCC grown locations (tumor) and left oriented green arrow show grown HFF locations (spot or homogeneous background), (a) Petri dish 1, (b) Petri dish 2 and (c) Petri dish 3. Row 2: Reference bioluminescence images D ref I of (d) Petri dish 1, (e) Petri dish 2 and (f) Petri dish 3. Row 3: reconstructed images with projection method for 3 second sampling rate (no undersampling). Raw 4: reconstructed images with projection method for 1 second sampling rate.

Fig. 6 :

 6 Fig. 6: Plot of the solid angle model function a ij = f (d, α ij ) with d the probe to source distance and α ij the tilt angle.

Fig. 7 :

 7 Fig. 7: (a) Representation of the LUT proposed in [24], where a ij are represented by a color-map, according to their position in the FOV. Here, distance from center of the source (X axis) and depth from the source to the detector (Y axis) are defined through d ij and α ij , (b) Partition model generated from k = 4 clusters with a p coefficients represented along distance from center of the source (X axis) and depth from the source to the detector (Y axis): a 1 in dark red, a 2 in green, a 3 in blue and a 4 in dark blue.

Fig. 8 :

 8 Fig. 8: Upper row: Difference of cropped bioluminescence and raw images, (a) Petri dish 1, (b) Petri dish 2 and (c) Petri dish 3. Lower row: bioluminescence images in the SOI coordinate system Rref (d) Petri dish 1, (e) Petri dish 2 and (f) Petri dish 3.

Fig. 9 :

 9 Fig. 9: Normalized reference, giving global size and position of tumors. (a) Petri dish 1, (b) Petri dish 2, (c) Petri dish 3.

Fig. 10 :

 10 Fig. 10: Reconstructed images varying MLEM iteration number with SA model for sampling time fixed at 3 seconds, in Petri dish 2. (a) 10 iterations, (b) 30 iterations, (c) 50 iterations.

Fig. 11 :Fig. 12 :

 1112 Fig. 11: Blurring effect of normalization function F N RM i , example with projection method, optimal input parameters P I on Petri dish 2 (a) Reconstructed image RM i (50 x 50 pixels), (b) Normalized reconstructed image RNM i (180 x 180 pixels).

Fig. 13 :

 13 Fig. 13: Mean AUC values along input parameters P I : (a) Varying sampling rate and MLEM = 25 iterations, (b) Varying MLEM iterations, sampling = 3 seconds.

Fig. 14 :

 14 Fig. 14: Mean computation time observed for reconstruction methods varying sampling rate.

Fig. 15 :

 15 Fig. 15: Mean Dice values along input parameters P I : (a) Varying sampling rate and MLEM = 25 iterations, (b) Varying MLEM iterations, sampling = 3 seconds.

Fig. 16 :

 16 Fig. 16: Mean T/B values along input parameters P I : (a) Varying sampling rate and MLEM = 25 iterations, (b) Varying MLEM iterations, sampling = 3 seconds.

Table 1 :

 1 Optimal points extracted from mean ROC curves, computed for the different Petri dish configurations.

	PM-k 3

1

Specificity.

2

Sensitivity.

3

Partition model witk k ∈ [2; 6] clusters selected.

A suxel is defined as the smallest element of the SOI.

Acknowledgments

The authors would like to thank CAMPUS FRANCE (PROCOPE Program), Rennes Métropole and Le Ministère des Affaires étrangères et du Développement international for French-Bavaria student exchanges and La Ligue Contre le Cancer for financial support.

Compliance with Ethical Standards

Funding: This study was funded by Campus France (FR), Rennes Métropole (FR), Le Ministère des Affaires Étrangères (FR) and la Ligue Contre le Cancer (FR).

Conflict of Interest: The authors declare that they have no conflict of interest. Ethical approval: This article does not contain any studies with human participants or animals performed by any of the authors. Informed consent: This articles does not contain patient data.


Pierre Matri
email: pmatri@fi.upm.es

Yevhen Alforov
email: alforov@dkrz.de

Alvaro Brandon
email: abrandon@fi.upm.es

María S Pérez
email: mperez@fi.upm.es

Alexandru Costan
email: alexandru.costan@irisa.fr

Gabriel Antoniu
email: gabriel.antoniu@inria.fr

Michael Kuhn

Philip Carns
email: carns@mcs.anl.gov

Thomas Ludwig
email: ludwig@dkrz.de

Álvaro Brandon

Mission Possible: Unify HPC and Big Data Stacks Towards Application-Defined Blobs at the Storage Layer

published or not. The documents may come

Introduction

HPC and Big Data platforms are carving new data storage models. This is made necessary by the everincreasing scale of the computation and of the datasets ingested and produced by large-scale applications. The success of key-value stores [START_REF]Apache Cassandra[END_REF][START_REF]Project Voldemort[END_REF] or block storage systems [START_REF]Openstack Swift[END_REF][START_REF] Weil | Ceph: A scalable, high-performance distributed file system[END_REF] on Clouds, and the advent of burst buffers [START_REF] Liu | On the role of burst buffers in leadershipclass storage systems[END_REF][START_REF] Romanus | Challenges and considerations for utilizing burst buffers in high-performance computing[END_REF] or advanced I/O libraries [START_REF] Liu | Hello ADIOS: the challenges and lessons of developing leadership class I/O frameworks[END_REF][START_REF] Bartz | A best practice analysis of HDF5 and NetCDF-4 using lustre[END_REF] for HPC clearly highlight this need.

At the heart of these different methods is the move from legacy POSIX-compliant storage systems towards simple storage paradigms designed especially for one purpose, trading versatility for performance. Indeed, POSIX-IO imposes functionality such as hierarchical namespaces or file permissions. While these features are often provided for convenience, they are in practice rarely needed by modern applications and can significantly hinder the storage performance. Indeed, the libraries and frameworks commonly used to access the storage on HPC [START_REF] Forum | Mpi: A message-passing interface standard[END_REF] and Big Data platforms [START_REF] Zaharia | Apache spark: A unified engine for big data processing[END_REF][START_REF] Shvachko | The Hadoop distributed file system[END_REF] provide relaxed semantics (i.e., the set of rules and guarantees provided by the system regarding the behavior of its storage operations) compared to those of the underlying file system.

Yet, deploying new storage models on HPC platforms used to be hard or simply impossible. Indeed, parallel file systems such as Lustre or GPFS on HPC have been the cornerstone of HPC storage for decades and are likely to remain so in the next few years. This is largely explained by the high level of versatility and support for legacy applications, which is without comparison with that of purpose-built storage systems. In contrast, this is easy on cloud computing platforms such as [START_REF]Microsoft Azure[END_REF]13], which enable users to deploy and configure exactly the storage system they need on compute nodes.

Application-defined storage on HPC comes as a solution. It leverages local storage on compute nodes available in a growing number of leadership-class supercomputers [14,[START_REF]Theta[END_REF] to allow scientists to deploy transient data services alongside the application. Such services offer the application exactly the semantics and fine-grained tuning it needs. Multiple examples of such services exist in the literature [START_REF] Zheng | Deltafs: Exascale file systems scale better without dedicated servers[END_REF][START_REF] Carns | Enabling NVM for dataintensive scientific services[END_REF]. This integration of the application with the storage it needs greatly eases its containerization and hence application portability across platforms.

We argue that deploying storage alongside the application additionally obviates the need for the aforementioned features of distributed file systems by removing the multi-user constraint. Blob-, or Object-based storage [START_REF] Weil | Rados: A scalable, reliable storage service for petabyte-scale storage clusters[END_REF][START_REF] Nicolae | Blobseer: Nextgeneration data management for large scale infrastructures[END_REF][START_REF] Matri | Týr: Blob storage meets built-in transactions[END_REF] has been demonstrated to provide an alternative to file-based storage on HPC and Big Data platforms. The reason is twofold. First, the flat namespace and simple semantics they provide enables performance improvements that are simply inaccessible to distributed file systems. Second, the data model they provide is close enough to that of file systems so most applications could use it with little to no modification.

In this article we leverage application-defined storage to assess the applicability and benefits of object-based storage for both HPC and Big Data platforms. Our contributions can be summarized as follows:

• After briefly describing our goals and reviewing related work (Section 2), we propose blobs as candidates for addressing the storage needs of HPC and Big Data (Section 3).

• We leverage a representative set of HPC and Big Data applications to prove that the vast majority of the I/O calls performed can be covered by stateof-the-art blob storage systems (Section 6).

• We describe the modifications necessary in the storage stack for HPC and Big Data applications to run atop blob storage (Section 7).

• We use an experimental testbed and a leadershipclass supercomputer (Section 5) to evaluate the performance benefits and trade-offs running these applications atop the same blob storage systems rather than traditional file-based storage (Section 8). We highlight a completion time improvement of up to 25% with blobs.

We finally conclude with future work that enhances our proposal (Section 9).

Related work and motivation

In this section we start by reviewing the state of the art regarding relaxing POSIX semantics on both HPC and Big Data applications, convergence between both these worlds and application-defined storage for HPC.

HPC: Relaxing POSIX-IO API and semantics

Increasingly large amounts of data are generated by HPC applications as the result of simulations and largescale experiments. Thus, storage systems need to provide concurrent access to the data for large numbers of tasks and processes. Such parallel storage operations rely on the usage of a parallel file system (PFS) implementing the POSIX-IO interface as the storage layer. Typical examples of such file systems used on most HPC platforms are Lustre [START_REF]The Lustre File System[END_REF] and OrangeFS [START_REF] Moore | OrangeFS: Advancing PVFS[END_REF].

Beyond the POSIX-IO interface lies the POSIX-IO semantics that a fully compliant file systems must implement. This standard has advantages regarding portability, but its inflexibility can cause considerable performance degradation [START_REF] Kimpe | Storage models: Past, present, and future[END_REF]. For example, this standard requires that changes made to a shared file must be visible immediately by all processes. Because an application has no way of telling the file system that POSIX-IO semantics are unnecessary or unwanted, it cannot avoid this performance penalty. For many file systems, these performance issues are noticeable even for small numbers of client processes and straightforward I/O patterns [START_REF] Cope | Bridging HPC and grid file I/O with IOFSL[END_REF][START_REF] Huang | UNIO: A unified I/O system framework for hybrid scientific workflow[END_REF]. The issues also affect higher levels of the I/O stack because an underlying POSIX-compliant file system effectively forces POSIX-IO semantics upon all other layers . For instance, this applies to the common HPC I/O stack with Lustre [START_REF] Kuhn | A Semantics-Aware I/O Interface for High Performance Computing[END_REF].

Big data: from file systems to object storage

Big data applications require a storage model that follows a write-once, read-many model. This requirement drove the design of many distributed file systems by sacrificing some of the POSIX-IO operations in order to gain data throughput. GFS [START_REF] Ghemawat | The Google file system[END_REF] implements only a set POSIX-compliant operations needed by data-intensive applications, namely, create, delete, open, close, read, and write. HDFS [START_REF] Shvachko | The Hadoop distributed file system[END_REF] is based on GFS and is designed to work in commodity hardware. It implements some additional POSIX-IO requirements such as directory operations and file permissions, but it discards some others such as concurrent reads and writes. Ceph [START_REF] Weil | Ceph: A scalable, high-performance distributed file system[END_REF] follows the same trend, discarding some POSIX-IO semantics and implementing only those that allow a distributed file system to work with most applications. GlusterFS eliminates the metadata server and claims to be fully POSIX compliant. However, work has shown that this compliance can impact throughput [START_REF] Mikami | Using the gfarm file system as a POSIX compatible storage platform for Hadoop MapReduce applications[END_REF].

We can clearly see a trend where the file system POSIX-IO API or semantics such as providing a hierarchical namespace, file permissions, or strict file access parallelism are unnecessary. Thus, they can be traded for performance and adaptability for Big Data.

Storage convergence between HPC and Big Data

During the past decade many research projects and workshops were dedicated to the opportunities and possibilities of running large scale scientific applications by using cloud computing technologies [START_REF]BDEC -Big Data and Extreme-Scale Computing[END_REF][START_REF] Zhang | Scientific computing meets big data technology: An astronomy use case[END_REF]. In general, many efforts there were made to investigate the performance of HPC applications (mostly form life sciences [START_REF] Duran-Limon | Efficient execution of the WRF model and other HPC applications in the cloud[END_REF]) on clouds (with and without virtualization [START_REF] Jaikar | Cloud computing: Read before use[END_REF]) highlighting cost efficiency or trade-offs [START_REF] Gupta | Evaluation of HPC applications on cloud[END_REF].

Several research efforts focus on building optimized or customized distributed-computing platforms that meet the requirements of HPC applications and scientific simulations [START_REF] Ledyayev | High performance computing in a cloud using OpenStack[END_REF]. Many of those are based on big data frameworks such as Spark or Hadoop / MapReduce. In contrast, Pan et al. [START_REF] Pan | Integrating high performance file systems in a cloud computing environment[END_REF] propose to port parallel file systems to cloud environments in order to support a wide range of applications expecting POSIX-IO on cloud applications. However, the application use cases considered in that work are rarely data-intensive. In the same way other researchers also aimed to provide the features of PFS in the cloud storage. For instance, Y. Abe and G. Gibson in [START_REF] Abe | pwalrus: Towards better integration of par-allel file systems into cloud storage[END_REF] presented a storage model which gives data access to a user through the storage service layer (S3 interface) and directly through a PFS.

Could blobs be the enabling factor?

Although the set of tools and techniques used for HPC and big data environment differ, many objectives are similar. The most important is probably to provide the highest-possible data access performance and parallelism. As such, the storage stack for HPC and Big Data looks similar. Indeed, the related work showed that a common trend for both HPC and big data is to relax many of the concurrent file access semantics, trading such strong guarantees for increased performance. Nevertheless, some differences remain. Specifically, while the big data community increasingly drops POSIX-IO altogether, the HPC community tends to provide this relaxed set of semantics behind the same API. Although this choice increases backwards compatibility with legacy applications, it also has significant performance impact. Very few HPC applications actually rely on strong POSIX-IO semantics. For instance, the MPI-IO standard does not only relaxes many of these semantics but also drops many of its operations altogether. For example, it does not expose the file hierarchy or permissions to the end user. Therefore, applications leveraging these libraries do not need these features to be provided.

Consequently, we ask ourselves whether the underlying storage technologies from both worlds could be unified, leading to specific software stacks and libraries running atop a low-level, low-opinionated storage paradigm, such as the one provided by blob storage systems (Figure 1). Indeed, blob-based storage systems such as Týr [START_REF] Matri | Týr: Blob storage meets built-in transactions[END_REF] or Rados [START_REF] Weil | Rados: A scalable, reliable storage service for petabyte-scale storage clusters[END_REF] could provide a strong alternative to file-based storage on both sides. These systems typically have a much more limited set of primitive compared with file systems:

• Blob Access: object read, object size,

• Blob Manipulation: object write, truncate,

• Blob Administration: create object, delete object,

• Namespace Access: scan all objects.

These operations are similar to those permitted by the POSIX-IO API on a single file. Therefore, most file operations performed on a file system can be mapped directly to the corresponding primitives of blob storage systems. In that model we classify file open and unlink as file operations.

In contrast, directory-level operations do not have their blob counterpart, because of the flat nature of the blob namespace. Should applications need them, such operations can be emulated using the scan operation. Obviously, this emulation is far from optimized. Yet, since we expect these calls to be vastly outnumbered by blob-level operations, this performance drop is likely to be compensated by the gains permitted by using a flat namespace and simpler semantics. While some of these features (e.g., permissions) are usually not used at the application level, they may be necessary at the system management level, for multiuser systems as provided on HPC environments. In the context of blob-based storage systems, we advocate that such access control can be provided using multiple keyspaces (e.g., per-application keyspaces), in which these features can be configured at the keyspace level.

Legacy application, which could rely on a fully compliant POSIX-IO interface, could leverage a POSIX-IO interface implemented atop such blob storage. This is proven possible by the Ceph file system, a file-system interface to Rados.

A representative set of applications

The challenges posed by convergence between HPC and Big Data applications have raised many discussions in the community [START_REF]BDEC -Big Data and Extreme-Scale Computing[END_REF]. One of the emerging ideas from these discussions is that the applications cannot be considered separately from the underlying software stack; the fuel for convergence could enable a wide variety of HPC and Big Data applications to leverage converging services and underlying infrastructure. Although designing a full converging stack is not the objective of this paper, we focus on storage as one of the critical milestones that could make such convergence possible. Specifically, we choose I/O-intensive applications, which could benefit most from such converged storage.

Accordingly, we base our experiments on a number of I/O-intensive applications extracted from the literature [START_REF] Fox | Towards a comprehensive set of big data benchmarks[END_REF][START_REF] Wang | BigDataBench: A big data benchmark suite from Internet services[END_REF] that cover the diversity of I/O workloads commonly encountered on both HPC (Section 4.1) and Big Data platforms (Section 4.2).

HPC Applications

The HPC applications we use are based on MPI. They all leverage either large input of output datasets associated with large-scale computation atop centralized storage usually provided by a distributed, POSIX-IOcompliant file system such as Lustre [START_REF]The Lustre File System[END_REF].

mpiBLAST [START_REF]mpiBLAST: Open-Source Parallel BLAST[END_REF] is a parallel MPI implementation of NCBI BLAST [START_REF] Altschul | Basic local alignment search tool[END_REF]. It is a read-intensive biomolecular tool searching for regions of similarity between biological sequences. It compares nucleotide or protein sequences to sequence databases and calculates the statistical significance.

MOM (Modular Ocean Model) [START_REF]Ocean circulation models[END_REF] is a read-intensive three-dimensional ocean circulation model targeted at understanding the ocean climate system.

ECOHAM5 (ECOlogical model, HAMburg, version 5) [42] is a write-intensive three-dimensional biogeochemical ecosystem model with the focus on the North Sea [START_REF] Lorkowski | Interannual variability of carbon fluxes in the North Sea from 1970 to 2006competing effects of abiotic and biotic drivers on the gasexchange of CO 2[END_REF][START_REF] Große | Looking beyond stratification: A model-based analysis of the biological drivers of oxygen depletion in the North Sea[END_REF]. It modelizes the pelagic and benthic cycles of carbon, nitrogen, phosphorus, silicon and oxygen on the northwest European continental shelf.

Ray Tracing is a balanced read-write workload extracted from the BigDataBench [START_REF] Wang | BigDataBench: A big data benchmark suite from Internet services[END_REF], itself derived from [START_REF] Stone | An efficient library for parallel ray tracing and animation[END_REF]. It generates images by tracing the path of light and simulating the effects of its encounters with virtual objects.

Big Data Applications

As the leading open-source Big Data processing and analytics framework, Apache Spark [START_REF] Zaharia | Resilient distributed datasets: A fault-tolerant abstraction for in-memory cluster computing[END_REF] appears as an ideal candidate for this research. Chosen applications are extracted from SparkBench [START_REF] Li | Sparkbench: a comprehensive benchmarking suite for in memory data analytic platform Spark[END_REF], a benchmarking suite for Spark. It comprises a representative set of workloads belonging to four application types: machine learning, graph processing, streaming, and SQL queries.

Sort is a widely used benchmark that reads input data and sorts it based on a given key. It is I/O-intensive since all the data read will be processed and written back to the file system. For example, it can be used to sort a series of readings from sensors by date.

Grep is a filtering workload that searches in the input data for lines containing a given word and saves these lines into HDFS. In contrast to Sort, the size of the input and output will not be equal, and some data will be filtered out.

Decision Tree is a machine learning workload that reads a dataset containing rows with a series of features and a class they belong to. This dataset is then split into a training and a test set. The workload creates a predictive model with the training set that is able to predict the class of the elements in the test set. These predictions are written back to disk.

Connected Component is an algorithm that finds the subgraphs in a graph in which any two vertices are connected by a path but are not connected to any other node on the supergraph. This can be seen as a way of finding clusters of nodes. For example, in a social network it can be used to find communities of users. The implementation in the benchmark will read a dataset and write the labels of each component back to disk.

Tokenizer is a Spark application we developed that reads a text file, tokenizes each line into words, and for each line calculates the NGrams=2 (i.e., contiguous sequences of 2 words from each line). These Ngrams are saved in a text file. This is a common preprocessing step in topic modeling for documents where word patterns have to be extracted as an input to a machine learning model. This application shows a write-intensive workload.

Experimental configuration

We run experiments using the Grid'5000 [START_REF] Balouek | Adding virtualization capabilities to the Grid'5000 testbed[END_REF] experimental testbed, which spans 11 sites in France and Luxembourg. In this paper we use 32 nodes of the parapluie cluster in Rennes. Each node embeds 2 x 12-core 1.7 Ghz 6164 HE, 48 GB of RAM, and 250 GB HDD. We use Gigabit Ethernet connectivity (MTU = 1500 B) for Big Data applications, and 4 x 20G DDR Infiniband for HPC applications. HPC applications ran atop Lustre 2.9.0 and MPICH 3.2, using multiple ratios of storageto-compute nodes. Big Data application ran atop Spark 2.1.0, Hadoop / HDFS 2.7.3 and Ceph Kraken.

All storage systems are configured with similar parameters to allow for a fair comparison. Specifically, each system is configured with a replication factor of 1 (no replication), and using a stripe size of 64 MB.

In addition, we prove that the results obtained with HPC applications are replicable to a high-end supercomputer. To do so, we performed extra experiments on the Theta supercomputer [START_REF]Theta[END_REF] hosted at the Argonne Leadership Computing Facility (ALCF). Theta is a lastgeneration 9.65-petaflop Cray XC40 system. It is composed of 3,624 nodes, each containing a 64 core Intel Knights Landing processor with 16 GB of highbandwidth in-package memory (MCDRAM), an additional 192 GB of DDR4 RAM, and a 128 GB local SSD.

Analyzing the distribution of storage calls

In this section we demonstrate that the actual I/O calls made by both HPC and big data applications are not incompatible with the set of features provided by stateof-the-art blob storage systems. Our intuition is that read and write calls are vastly predominant in the workloads of those applications and that other features of distributed file systems such as directory listings are rarely used, if at all.

Tracing HPC applications

Figure 2 summarizes the relative count of storage calls performed by our set of HPC applications. The most important observation for all four applications is the predominance of reads and writes. Except for ECO-HAM, no application performed any other call to the storage system that reads or writes files, thus confirming our first intuition. This was expected because the MPI-IO standard does not permit any other operation.

The few storage calls other than read and write (mainly extended attributes reads and directory listings) are due to the run script necessary to prepare the run and collect results after it finishes. These steps can be performed offline from the I/O-heavy MPI part of the application. This results in only reads and writes being performed (EH / MPI).

We conclude that the only operations performed by our set of HPC applications, namely, file I/O, can be mapped to blob I/O on a blob storage system. Consequently, these applications appear to be suited to run unmodified atop blob storage. Similar to what we observed with HPC applications, the storage calls are vastly dominated by reads and writes to files. In contrast with HPC, however, all applications also cause Spark to perform a handful of directory operations (86 in total across all our applications). These directory operations are not related to the data processing because input / output files are accessed directly by using read and write calls. Analyzing these directory operations, we notice that they are related solely to (a) creating the directories necessary to maintain the logs of the application execution, (b) listing the input files before each application runs if the input data is set as a directory, and (c) maintaining the .sparkStaging directory. This directory is internally used by Spark to share information related to the application between nodes and is filled during the application submission. It contains application files such as the Spark jar or the application jar, as well as distributed cache files [START_REF]Running Spark on YARN[END_REF].

Tracing Big Data applications

We analyze in detail the directory operations performed by big data applications. Table 2 shows the breakdown of all such directory operations across all applications by storage call. We note that only the input data directories are listed, meaning that Spark accesses directly all the other files it needs with their path. Consequently, a flat namespace such as the one provided by blob storage systems could probably be used.

Big Data: Hierarchical to flat namespace

We previously observed that Spark storage calls included several very rare directory operations. Our observations concluded that such calls are only performed for separating data files from temporary files and are not strictly necessary for the application itself. Yet, one of our objectives is to prove that Spark applications can run unmodified atop blob storage. We prove that Spark can run atop a flat namespace providing only object-level operations, as provided by blob storage systems; such a proof would confirm our previous assertions.

Since file operations dominate directory operations, we choose to optimize the former at the expense of the latter. Thus, for any given hierarchical file path we generate a predictable flat path. Subsequent operations to that path are translated to a file operation on the rewritten path. In HDFS we achieve this by storing all files at the root. That is, we store on path /foo bar a file that would normally be stored on /foo/bar in a hierarchical namespace. Listing or deleting a directory is implemented by scanning and filtering the whole set of files in the system, selecting only the matching files that would be contained in that folder. Although such scan operations are costly, they are infrequent. Consequently the performance gains by flattening the file system (i.e., not managing permissions) should outweigh the cost of such operations. Table 3 summarizes the rewrite rules we apply on incoming storage calls from Spark.

We intercept and rewrite storage calls by modifying Hadoop / HDFS file storage interface. We run all Big Data applications 100 times and compare the average completion time of each benchmark suite on HDFS us- ing the original hierarchical namespace or simulating a flat namespace. We plot the results in Figure 4. Flattening the namespace on HDFS does not result in any significant task completion time variation despite the higher completion time of directory operations, plotted in Figure 5. Indeed, these calls are diluted in the vastly superior amount of file read and write calls. We demonstrated that running our set of Big Data applications over a simulated flat namespace not only is possible but also does not cause any significant performance variation. Consequently, these applications also appear to be suited to run atop a blob storage system, thus further enhancing the performance of the application by leveraging the reduced complexity of managing a flat namespace.

Replacing file-based by blob-based storage

In this section we demonstrate the potential of blobbased storage to suit the storage needs of both HPC and Big Data applications. We deploy each application listed in Section 4 atop state-of-the-art blob storage systems. We detail these systems in Section 8.1. We prove that the performance of these applications running atop converged blob-based storage matches or exceeds that of the same applications running atop Lustre for HPC (Section 8.2) and HDFS for Big Data (Section 8.4).

Overview of the blob storage systems

We run our applications atop two state-of-the-art blob storage systems: Týr [START_REF] Matri | Týr: Blob storage meets built-in transactions[END_REF] and Rados [START_REF] Weil | Rados: A scalable, reliable storage service for petabyte-scale storage clusters[END_REF]. We are using only the basic blob storage functionality they provide, and do not make use of any advanced features they may support. Although their high-level design has similarities, these two systems have different strengths and weaknesses resulting from design decisions made for each to support specific use cases.

Týr is a large-scale blob storage system designed around the same design principles as the Dynamo key-value database [START_REF] Decandia | Dynamo: Amazon's highly available key-value store[END_REF]. It is targeted at high access parallelism using multiversion concurrency control (MVCC) associated with built-in multiobject transactions. Týr offers fine-grained random write access to data, as well as single-hop reads (i.e., accessing the storage server without prior communication with any metadata server).

Rados is the storage layer for Ceph FS [START_REF] Weil | Ceph: A scalable, high-performance distributed file system[END_REF]. Rados has the ability to scale to thousands of hardware devices by using management software that runs on each of the individual storage nodes. The software provides features such as thin provisioning, snapshots and replication.

We assess the performance impact of replacing filebased with blob-based storage by observing three metrics. The job completion time is the total execution time of the application, from submission to completion. Read bandwidth and write bandwidth respectively represent the average data transferred per unit of time for read and write requests. We collect all these metrics on the compute nodes by instrumenting the adapter, and we aggregate the results.

Replacing Lustre with blob-based storage on HPC

In this section, we demonstrate how blob-based storage system can be used to support HPC applications while matching or exceeding Lustre I/O performance by replacing the latter with Týr and Rados.

We experiment using three storage-to-compute node configurations in order to ensure that our results are independent of the cluster configuration. We run the same experiments respectively with 28 compute / 4 storage nodes, 24 / 8 and 20 / 12. We average the results of 100 experiment runs.

On each node, we deploy a small interceptor to redirect POSIX storage calls to the blob storage system. It is based on FUSE [START_REF] Ishiguro | Optimizing local file accesses for fuse-based distributed storage[END_REF], which is supported on most Linux kernels today. In that configuration, this interceptor acts as the HPC adapter as presented in Figure 1. This adapter translates file operations to blob operations according to Table 4. Directory operations are not supported as we showed previously that they are unnecessary for HPC applications. The APIs of the blob storage systems we consider allow for a direct mapping between file-based and blob-based storage operations. We partially implement the stat function. Specifically, the file size is mapped to the blob size, the permissions are set to 777, the block size and allocated block size are set to 512 bytes, and the inode number is set to the hash of the blob key. The remaining information is set to 0. Our implementation does not support symbolic or hard links, which are not needed by our applications.

In Figure 6 we plot the average aggregate read and write bandwidth for all applications while varying the compute-to-storage node ratio. We note that for our configuration the 24 compute node / 8 storage node setup results in the higher bandwidth for all storage systems. Hence, the following experiments are performed with that configuration. This ratio is much lower than on common HPC platforms (3:1 vs. ∼ 70:1 at ORNL, for instance [START_REF] Oral | Best practices and lessons learned from deploying and operating large-scale data-centric parallel file systems[END_REF]) mainly because the jobs we run are significantly more data-intensive than compute-intensive. We note from these results that blob storage systems constantly outperform Lustre in all configurations for both reads and writes. We will detail these results in the following experiments. For read-intensive applications such as BLAST and MOM, this performance increase allows blob storage systems with 4 storage nodes to achieve a bandwidth comparable to Lustre's with 8 storage nodes.

In Figure 7 we plot the average read bandwidth for each of our HPC applications with Lustre file-based storage and Týr or Rados blob-based storage. We note an average 14% reduction of the total read time when using blob-based storage compared to Lustre. This is because of the optimized write path of the two blob storage systems considered. Indeed, both enable clients to locate and access any piece of data directly without prior communication with any dedicated metadata node. Although both blob storage systems behave similarly with respect to read performance, Rados shows a higher read performance due to its lower read consistency. We plot in Figure 8 the average write bandwidth obtained in similar conditions. These results show that blob storage systems outperform Lustre for BLAST, MOM and ECOHAM. Interestingly, for Ray Tracing, Rados performace drops behind that of other systems. This is the result of the relatively small writes (8 KB on average) performed by this application compared to the other systems, causing high lock contention on concurrent writes that hinders the throughput of the system.

In Figure 9 we plot the average application completion time improvement. The I/O performance gains are here diluted in compute operations. As expected considering the previous results, read-intensive applications exhibit the greatest decrease. BLAST and MOM show a completion time reduction of nearly 8% with both blob storage systems. In contrast, write-intensive applications such as ray tracing show a lower 3% completion time decrease with Týr or Rados as the underlying storage when compared with Lustre.

Replicating results on a high-end supercomputer

In this section, we seek to prove that the experiments obtained in Section 8.2 are reproducible on a high-end supercomputer. To do so, we leverage the Theta supercomputer hosted at Argonne National Laboratory, and run the same experiments as in the aforementioned section. Although arguably not an easy task due to the strong limitations of the platform and the lack of superuser rights, blob storage systems such as Týr and Rados are deployable on such platform.

We deploy the applications as described in the previous section, using 32 nodes (totaling 2048 cores), using 24 nodes for computation and 8 nodes for storage, and measure the completion time for each application. Experiments with Lustre were using the file system available to the computer, totaling 170 storage nodes and shared across all users.

We plot the results in Figure 10. We show the performance improvement to be significantly higher than on our testbed. The reason for this significant perfor- mance increase is to be found in the technical characteristics of Theta, which offers significantly more RAM than SSD space. As such, for the most part, the storage systems deployed on these nodes behave as in-memory storage systems. We acknowledge that the setup of this platform is particular in this regard and by consequence that the results are not representative of those of another platform with different setup. Yet, we advocate that the results reach the goal of demonstrating that deploying blob storage systems on high-end HPC platforms is possible without requiring any application modification. We observe a higher variance in the results compared to Grid'5000, which we attribute to the shared nature of the centralized storage.

In Figure 11 we scale all four applications on up to 128 nodes of Theta, or 8192 cores. Because of allocation limitations, each experiment was performed only 5 times. We notice a near-linear decrease of computation time as the size of the cluster increases. With applications applications such as ECOHAM or Ray Tracing, the performance improvement is slightly lower than with purely read-intensive applications due to the higher cost of write operations compared to read operations.

Running Spark applications atop blob storage

In this section we run the same set of experiments for the set of Big Data applications. We demonstrate that Týr and Rados significantly outperform HDFS for all applications. In order to provide an additional baseline of the performance of file systems, we also run these applications atop CephFS [START_REF] Weil | Ceph: A scalable, high-performance distributed file system[END_REF], itself based on Rados.

We use the same configuration as in Section 7, running computation alongside storage on 32 nodes. We integrate the storage adapter for blob storage directly inside HDFS. The Hadoop installation has been modified to redirect storage calls to blob storage systems. The translation between POSIX-like calls and flatnamespace blob operations is done by using the translation rules defined in Table 3. We implement the CRUSH algorithm to provide Spark with the physical data location on Rados and CephFS. We use the client API capabilities to provide that information with Týr.

In Figure 12 we plot the average read bandwidth achieved for each of the Spark benchmarks. We notice a striking read bandwidth improvement when using Týr and Rados over HDFS (the read bandwidth is increased by 28% and 32% respectively). Same as in HPC, this is due mainly to the direct read feature of the two storage systems that, unlike HDFS, enable the applications to bypass any centralized metadata service and access the storage servers directly. CephFS performance highlights the cost of file-based storage by showing a degraded performance compared to Rados. As with HDFS, this is mostly due to the additional communication required with dedicated metadata servers in the critical path for read requests, made necessary by the file hierarchy management.

Figure 13 shows the average write bandwidth for these applications. Similar to what was observed with HPC, we note a constant write bandwidth improvement with blob storage over HDFS and CephFS. We also note the same pattern we observed with HPC. Specifically, Rados outperforms Týr on read-intensive applications, whereas Týr enables higher throughput on writeintensive applications. This is visible with the Tokenizer application, where lock contention due to lack of multiversion concurrency control in Rados causes significant performance loss on concurrent write access.

In Figure 14 we plot the relative improvement in the total application completion time, diluted in computation. Running Big Data applications atop blobs improves application cpmpletion time, up to 22% compared to HDFS and 7% compared to CephFS. For Big Data, the highest gains are obtained with read-intensive applications such as Grep and Decision Tree. In comparison, write-intensive applications such as Tokenizer also benefit from improved performance, although relatively smaller due to the globally greater complexity of the write protocols for each storage system.

Conclusion

In this paper we argue that blob storage is a strong candidate for replacing traditional storage for both HPC and Big Data. Its simple data model is enough to map directly file operations to blob operations. Based on the previous observation that simple file reads and writes constitute the vast majority of the storage calls made by both HPC and Big Data applications, we factually prove that this convergence is possible by mapping both HPC and Big Data applications to blob storage. This does not require any modification in the application thanks to a thin adapter layer between the application and the persistent storage. We leverage 4 real-life HPC applications as well as 5 Big Data benchmarks to prove on an experimental testbed that not only such convergence is possible, but that it also significantly improve performance by up to 25% with read-intensive applications. We confirm on the Theta supercomputer that this setup is reasonable and applicable to a high-end supercomputer with near-linear scalability up to 8,192 cores.

In future work we will experiment both Týr and Rados at the same time on a real supercomputer as well as on one of the leading cloud computing platforms, both with a larger set of applications and frameworks. We will compare our approach with an extensive set of competitor storage systems and platform configurations.

Figure 1 :

 1 Figure 1: Converged side-to-side storage stack for HPC and Big Data.

Figure 3 Figure 2 :Figure 3 :

 323 Figure3shows the relative count of storage calls performed by our set of big data applications to HDFS.

Figure 4 :Figure 5 :

 45 Figure 4: Average Big Data task completion time with and without flat namespace simulation, with 95% confidence intervals.

Figure 7 :

 7 Figure 7: Comparison of read throughput for each HPC application with Lustre, Týr and Rados, with 95% confidence intervals.

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Comparison of write throughput for each HPC application with Lustre, Týr, and Rados, with 95% confidence.

Figure 11 :

 11 Figure 11: Average performance improvement at scale relative to 32 nodes setup for HPC applications using blob-based storage on Theta.

Figure 14 :

 14 Figure 14: Average performance improvement relative to HDFS for Big Data applications using blob-based storage.

Table 1 :

 1 Application summary

	Platform	Application	Usage	Total reads	Total writes	R / W ratio	Profile
		mpiBLAST (BLAST)	Protein docking	27.7 GB	12.8 MB	2.1 × 10 3	Read-intensive
	HPC / MPI	MOM ECOHAM (EH)	Oceanic model Sediment propagation	19.5 GB 0.4 GB	3.2 GB 9.7 GB	6.01 4.2 × 10 -2	Read-intensive Write-intensive
		Ray Tracing (RT)	Video processing	67.4 GB	71.2 GB	0.94	Balanced
		Sort	Text Processing	5.8 GB	5.8 GB	1.00	Balanced
		Connected Component (CC)	Graph Processing	13.1 GB	71.2 MB	0.18	Read-intensive
	Cloud / Spark	Grep	Text Processing	55.8 GB	863.8 MB	64.52	Read-intensive
		Decision Tree (DT)	Machine Learning	59.1 GB	4.7 GB	12.58	Read-intensive
		Tokenizer	Text Processing	55.8 GB	235.7 GB	0.24	Write-intensive

Table 2 :

 2 Spark directory operation breakdown

	Operation	Action	Count
	mkdir	Create directory	43
	rmdir	Remove directory	43
	opendir (Input data directory)	Open / List directory	5
	opendir (Other directories)	Open / List directory	0

Table 3 :

 3 Big Data storage call translation rules

	Original operation	Rewritten operation
	create(/foo/bar)	create(/foo bar)
	open(/foo/bar)	open(/foo bar)
	read(fd )	read(bd )
	write(fd )	write(bd )
	mkdir(/foo)	Dropped operation
	opendir(/foo)	scan(/), return all files matching /foo *
	rmdir(/foo)	scan(/), remove all files matching /foo *

Table 4 :

 4 HPC storage call translation rules

	POSIX Call	Translated Call
	create(/foo/bar)	create(/foo bar)
	open(/foo/bar)	open(/foo bar)
	read(fd )	read(bd )
	write(fd )	write(bd )
	mkdir(/foo)	Unsupported operation
	opendir(/foo)	Unsupported operation
	rmdir(/foo)	Unsupported operation

Acknowledgments

This work is part of the "BigStorage: Storagebased Convergence between HPC and Cloud to handle Big Data" project, H2020-MSCA-ITN-2014-642963, funded by the European Commission within the Marie Skłodowska-Curie Actions framework. This work was supported by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research under Contract DE-AC02-06CH11357. It is also supported by the ANR OverFlow project, ANR-15-CE25-0003. Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER, and several universities and organizations.