Corentin Sanchez

Philippe Xu

Alexandre Armand

Philippe Bonnifait

Spatial Sampling and Integrity in Lane Grid Maps

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

A key capability of autonomous vehicles (AV) is to take the best decision in every situation it encounters. In order to achieve this task, the decision-making (DM) module in the architecture of the AV needs to have a good understanding of the surrounding of the ego vehicle. Operating in urban areas remains challenging for autonomous driving systems. Thus, the information provided to the DM module must not be misleading so that proper decisions can be made. For example, characterizing a driving space as free when it is actually occupied may lead to an hazardous decision. Our work stands between the DM module and its raw input data provided by the perception and localization modules. Indeed, as data is provided by sensors, information is extracted by the perception module but still needs to be understood in order to take the best decision. This paper aims at providing an augmented information representation that takes into account information brought by the perception as well as information stored in a multi-layer map. The purpose is to get a better understanding of the situation that is supplied to the DM module. Such a representation can be generated within our This work has been carried out within SIVALab, a shared laboratory between Renault and Heudiasyc (UTC/CNRS), and financed by the CNRS.

1 Université de technologie de Compiègne, CNRS, Heudiasyc, CS60319 -60203 Compiègne Cedex, France 2 Renault S.A.S, Guyancourt, France In several architectures, the DM system can be seen as a tripartite module [START_REF] Ulbrich | Towards a Functional System Architecture for Automated Vehicles[END_REF]. It works under three abstraction levels, which correspond to different levels of decision. At the lowest level (operational level), there is the local planning system that handles local and short term trajectories [START_REF] Hundelshausen | Driving with tentacles: Integral structures for sensing and motion[END_REF]. At the intermediate level (tactical level), the maneuver planner takes decisions on tasks that have to be achieved at a lane level [START_REF] Orzechowski | Decision-making for automated vehicles using a hierarchical behavior-based arbitration scheme[END_REF]. Finally, at the most abstracted layer (strategic level) the global planner plans long-term decisions such as the itinerary to follow to complete the mission. In order to plan trajectories at a low level, two types of approaches can be highlighted. The first one is to use objects lists and to navigate through the obstacles [START_REF] Althoff | Model-Based Probabilistic Collision Detection in Autonomous Driving[END_REF], [START_REF] Papp | World modeling for cooperative intelligent vehicles[END_REF]. In this case there is no explicit information on the absence of objects, and occlusions have to be handled [START_REF] Sun | Behavior Planning of Autonomous Cars with Social Perception[END_REF], [START_REF] Schorner | Predictive Trajectory Planning in Situations with Hidden Road Users Using Partially Observable Markov Decision Processes[END_REF]. In addition, there is a need to get spatial information. Another approach is to use occupancy grids [START_REF] Mouhagir | Using evidential occupancy grid for vehicle trajectory planning under uncertainty with tentacles[END_REF] or free space polygons [START_REF] Schreier | Compact Representation of Dynamic Driving Environments for ADAS by Parametric Free Space and Dynamic Object Maps[END_REF] to plan the trajectories. However there is a limitation if the vehicle operates in a complex environments with multiple interactions. The use of maps is necessary to add context, and thus to facilitate situation understanding [START_REF] Ulbrich | Defining and Substantiating the Terms Scene, Situation, and Scenario for Automated Driving[END_REF]. Using maps implies to have reliable localization of the ego vehicle, in order to project perception information onto the map. Therefore, the integrity of the localization needs to be properly considered [START_REF] Hage | Localization Integrity for Intelligent Vehicles Through Fault Detection and Position Error Characterization[END_REF]. At the intermediate level, the maneuver planner enables to take into account this context as it is shown in the "environment model" mentioned in [START_REF] Orzechowski | Decision-making for automated vehicles using a hierarchical behavior-based arbitration scheme[END_REF].

In this work, we propose a lane level information representation at tactical level, that handles occlusions and enables to improve confidence on information. In addition, it reduces the risk of providing misleading information to the DM module. An HD map is used to infer the lanes of interest as presented in [START_REF] Sanchez | Lane level context and hidden space characterization for autonomous driving[END_REF], and to infer their occupancy through an along-track sampling as illustrated in figure 1. The sampling step will be the key stone to handle the integrity of the generated information. We propose an evaluation method which combines both spatial information and distance metrics. These are our main contributions. Related work is discussed in Section II. Section III introduces the Lane Grid Map and its principles. The evaluation method is explained in section IV and is followed by the integrity bought by this representation in section IV-B. Finally, results are presented in section V.

II. RELATED WORK

In our world model representation, as previously introduced in [START_REF] Sanchez | Lane level context and hidden space characterization for autonomous driving[END_REF], in a similar correspondence with the DM module, the notion of abstraction levels is early introduced with [START_REF] Michon | A Critical View of Driver Behavior Models: What Do We Know, What Should We Do?[END_REF] and reused for decision-making [START_REF] Ulbrich | Towards a Functional System Architecture for Automated Vehicles[END_REF]. There are three levels. From the lowest to the highest we have respectively the operational level, the tactical level and the strategic level. Each level operates in its own representation area from local with the lowest level to global with the highest.

Information representation used as inputs to these levels of the decision-making is therefore a key point. First, we need to dissociate discrete representation like with objects with spatial representations. In [START_REF] Papp | World modeling for cooperative intelligent vehicles[END_REF] or in [START_REF] Lundquist | Estimation of the Free Space in Front of a Moving Vehicle[END_REF], objects and features are used to represent the vehicle surroundings. Maps also contain at each level of abstraction different information: topological, semantic and geometric information as in [START_REF] Poggenhans | Lanelet2: A high-definition map framework for the future of automated driving[END_REF]. However, unlike objects, spatial information enables to characterize surfaces and to give information on missing information. Among this kind of representation we could cite grids: probabilistic grids [START_REF] Chen | Dynamic Environment Modeling with Gridmap: A Multiple-Object Tracking Application[END_REF], evidential grids [START_REF] Moras | Credibilist occupancy grids for vehicle perception in dynamic environments[END_REF], interval maps [START_REF] Weiherer | An interval based representation of occupancy information for driver assistance systems[END_REF] or even more than 2D representations as with [START_REF] Danescu | A Particle-Based Solution for Modeling and Tracking Dynamic Digital Elevation Maps[END_REF]. There is also a more compact representation with polygon approaches as with the parametric free space presented in [START_REF] Schreier | Compact Representation of Dynamic Driving Environments for ADAS by Parametric Free Space and Dynamic Object Maps[END_REF]. All these representations are built thanks to the perception system.

For decision-making evaluation, the metric often used is distances based on the situation encountered by the vehicle. It can be expressed as Time-To-Collision (TTC) or Enhanced-Time-To-Collision (ETTC) [START_REF] Minderhoud | Extended time-to-collision measures for road traffic safety assessment[END_REF] or even Time-To-Entrance (TTE) [START_REF] Noh | Decision-Making Framework for Autonomous Driving at Road Intersections: Safeguarding Against Collision, Overly Conservative Behavior, and Violation Vehicles[END_REF]. Even more complex evaluation models are used for safety evaluation such as the Responsibility-Sensitive Safety (RSS) model [START_REF] Shalev-Shwartz | On a Formal Model of Safe and Scalable Self-driving Cars[END_REF]. We therefore want to combine the advantages of both representations, i.e using distances and spatial occupancy, to supply non-misleading information at the maneuver level.

In order to manage confidence of information and risk to provide misleading information, the notion of integrity can be cited. This criterion has been developed and formalized in the aviation field. It has been explored for urban environments as presented in [START_REF] Zhu | GNSS Position Integrity in Urban Environments: A Review of Literature[END_REF]. This notion of integrity has also been applied for localization systems for instance [START_REF] Hage | Localization Integrity for Intelligent Vehicles Through Fault Detection and Position Error Characterization[END_REF]. Given an integrity risk, the system should not underestimate the uncertainty we can have on input values.

In this paper we want to enhance spatial information by using a low level spatial information representation but also give a distance metric and confidence over the information supplied to the decision-making module. We reuse the notion of integrity to guarantee a functional operating domain for the AV through the sampling method of the Lane Grid Map. We show that oversampling can lead to integrity loss.

III. LANE GRID MAP

In a world model system, the three levels of abstraction provide a hierarchy of levels with several types of information representation. Each level has a specific type of task to achieve. As this paper focuses on the spatial representation layer, it aims at describing and evaluating a tactical level representation. This representation shown in figure 1 is called a Lane Grid Map (LGM). This section aims at describing our approach on an intermediate level representation that contains information with a high integrity level for handling non-misleading situation representation in order to enhance the DM process. Indeed, in order to be able to get on the road, an absence of vehicle is not enough. It is also necessary that the space is free and that this information has a high confidence value. We show that this representation enables to handle and guarantee a confidence level of integrity of the information taking into account localization uncertainty thanks to a sampling step parameter.

A. Principles

In order to build the LGM, we need to define the areas that are of interest for the AV. At the strategic level we are able to construct the Areas Of Interest (AOI) tree of the AV under the Interaction Graph named representation (right part of figure 2). Each node of the Interaction Graph is a portion of the AOI that describes in two ways the most important part of the vehicle surroundings. These AOI are built by the primary order lanes and secondary order lanes. Primary lanes are road lanes that are in direct interaction with the ego vehicle path. Secondary lanes correspond to lanes that are in direct interaction with the primary lanes but not in direct interaction with the ego vehicle path. These primary and secondary order represent road lanes extracted from the topological layer of an HD map where the ego vehicle must focus its attention in order to take best decisions and should not provide misleading information. These concepts are detailed in [START_REF] Sanchez | Lane level context and hidden space characterization for autonomous driving[END_REF].

At the tactical level, we use the Interaction Graph nodes and the geometric information of the HD Map to build a lane level spatial representation (center part of figure 2). Thus, the LGM depends on the context in which the AV drives and is a discrete representation of the AOI made of cells.

B. Spatial sampling process

The cells are contiguous and have a sampling step parameter. In our case, at the tactical level, the width of cells corresponds to the road lane width. The polyline of the center of the road lane is provided by the HD map. We select a sampling step parameters to split this polyline in several parts. Each segment is then extrapolated on its both sides to build single cells. Through this sampling step parameter, we are able to refine more or less the level of detail of the is primary of

Operational Level

Tactical Level Strategic Level information. Bigger are the cells,i.e the less number of cells we have, the less information we have about the situation.

Occupancy Grid Lane Grid Map Interaction Graph

The LGM aims at augmenting the operational spatial representation (left part of figure 2) to provide non-misleading information. It gives information of spatial occupancy but also with a notion of along-track distances.

C. Characterization process

Within the true world, the cells of the LGM are either Free (F ) or Occupied (O). However, in practice, some of the cells may be hidden because of occlusions and should be handled carefully. Therefore, in the characterization process of the cells of LGM, the Hidden (H) class is considered explicitly along with F and O. Note that some cells are not observed not because they are hidden but because they are out of the field of view of the sensors. These cells will be ignored during the evaluation process.

The classification of the cells into one of the three aforementioned classes can be done via the use of probabilistic occupancy grids [START_REF] Chen | Dynamic Environment Modeling with Gridmap: A Multiple-Object Tracking Application[END_REF] or more robustly with evidential grids [START_REF] Moras | Credibilist occupancy grids for vehicle perception in dynamic environments[END_REF]. Another simpler way to achieve this characterization is using a free space polygon [START_REF] Schreier | Compact Representation of Dynamic Driving Environments for ADAS by Parametric Free Space and Dynamic Object Maps[END_REF]. This type of approach has been implemented in this work.

D. Efficiency and scalability advantages

In classical occupancy grids (OG), all the space surrounding the vehicle is characterized regardless of their relevance. On the contrary, the LGM takes into account the situation in which the AV drives, i.e. in a given road layout, therefore it only stores information within the relevant lanes (AOI) provided by the Interaction Graph at the strategic level. As a consequence, the LGM encodes the information of the space surrounding the vehicle in a more compact and efficient way. Each cell of the LGM has two values: the width c w , which corresponds to a lane level width and a sampling step length c l . At the opposite, the OG is a spatial sampling over the entire space given square cells of width s w . Let us consider an area of 50 m over 50 m (A = 2500m 2 ) that needs to be encoded in both representations with the same cells size.

The OG has a number of cells C OG = A visibility distance d LGM = A cw×c l ×c l = A cw m. For example:

c w = 1m c l = s w = 1m ⇒ d LGM = 50 × d OG with C LGM = C OG (1)
The increase of the visibility distance is easily explained by the fact that cells are located in a focused area for the autonomous driving. Computationally we could raise the number of lanes by decreasing the visibility distance that would be dispatched on each lane.

It should be noticed that the LGM representation is easily scalable. Indeed, we can vary the sampling step and easily aggregate cells with low difficulty. On top of that, the fact that cells fit properly the road layout enables to have a good fidelity information representation over distances. Thus, this representation has the advantage to supply a coupled information. According to the situation encountered, it supplies occupancy information with distances.

In order to compare, if we look at OG used in the literature, these phenomena are two limitations. As cells are independent from the geometrical layer of the map, the delimitation of a number of cells to take into account for a portion of the road would not be easily defined. Distances over a grid are also not easy to fit as it does not reason in the coordinate system of the road lane.

IV. PERFORMANCE METRICS

This section presents the evaluation method used to demonstrate the effectiveness of the Lane Grid Map. This methodology introduces the integrity criterion. The metrics presented aim at providing two types of information in a single representation. At maneuver level for the DM, the LGM supplies both distances and occupancy for the task to be fulfilled.

A. Evaluation metrics 1) Ground truth generation: In order to validate the LGM model and efficiency, we want to compare the characterization process using the ground truth position of several vehicles. Given the precise position of several vehicles on the road and their dimensions, the occupancy of these vehicles are approximated as rectangles. All the cells overlapped by these rectangles are considered as occupied on the ground truth LGM. Then, the rest of the cells are considered as being free. The cells of the ground truth LGM are either Free (F ) or Occupied (O). The top row in figure 3 shows the characterization of the LGM by the ground truth.

2) Observation generation: At the same time, we apply the characterization process with the on-board perception system. Using polygons approach, we have the two types of information: the convex bounding polygons of the perceived objects and a free-space polygon. The perception system is

Observation Ground Truth

f o h F F : f → n 1 F : o → n 2 F : h → n 3 O O : f → n 4 O : o → n 5 O : h → n 6
described in section V. The bottom row in figure 3 shows the characterization by the ego vehicle of the cells of the LGM into one of the three states: F , O and H.

Over the ego vehicle trajectory, we record at each time step t ∈ [0, T ] the LGM (t). At each time step we have the LGM characterized by the AV perception and also by the ground truth of the objects. Combining the different configurations, for each cell of the LGM, there are six possible indicators N i obtained from the combination of both characterization as defined in table I. The correspondence is shown as an example on the indicator row of figure 3. As a cell depicts the occupancy of a fixed step distance, we can express these results not only in terms of cells but also in terms of distance and more precisely with a percentage. All the cells c k of each LGM are associated with a distance n k j ∈ {n 1 , . . . , n 6 }. The indicators can be written as follows:

N i = t∈[0,T ] c k ∈LGM (t) δ ij * n k j with δ ij = 1 if i = j 0 if i = j (2)
Two indicators are particularly relevant for the integrity analysis. N 2 corresponds to miss-classified free cells considered as occupied. It is an overcautious indicator. N 4 is the most important indicator as it corresponds to the missclassified occupied cells considered as free. It is a highly miss-leading indicator that should be as low as possible. Indeed, if we give information of free space around the vehicle whereas it is occupied it can leads to a hazardous decision-making. At the opposite, if occupied information is given instead of a free one, it remains a cautious situation for the DM. In section V we analyze the rates of these two indicators. The False Negative Rate (F N R) is the rate of N 4 over all occupied cells whereas the False Positive Rate (F P R) is the rate of N 2 over all free cells.

F N R = N4 N4+N5+N6 , F P R = N2 N1+N2+N3 (3)

B. Integrity with localization and perception uncertainties

For autonomous navigation tasks, the spatial occupancy characterization must not provide misleading information. The goal is to provide information with a high level of confidence with the notion of integrity. For a given integrity risk level, the characterization error represented by the F N R must be below the set risk level. We will show that the main contribution brought by the LGM is its ability to keep the integrity of the information thanks to its capacity to control the F N R via the sampling step which has an impact on the overall LGM characterization performance.

In addition to that, we need to evaluate the integrity with the impact of the localization errors. As HD maps

Free Occupied Hidden

Step 1

Step i

Step i with standardization Fig. 4: Aggregation and standardization of the LGM from initial step to another.

are used to create the LGM in a global frame coordinates, perceived information in the AV local frame is transformed into the global frame thanks to the localization position. As a consequence, when the uncertainty of the localization is taken into account, the characterization process has to provide, whatever the situation is, a confidence level of the information stored. In order to evaluate the robustness of our model, we added some Gaussian noises with different standard deviations on the pose of our ground truth. Thus, we can study the impact of the localization uncertainty on our system.

At first, we consider that there is a noise on the localization but we do not model it in the perception outputs. We want to show that despite uncertainties, we are still able to manage the characterization process through the sampling step.

Usually localization estimation also quantifies its uncertainty, typically with a covariance matrix. This uncertainty can be propagated to the polygons by enlarging them. As localization errors on the AV position are introduced, the uncertainty is propagated on the object hulls when they are transformed from the local coordinates of the car to the world coordinates. We want to show that if we take a maximum allowed threshold of error for the misleading indicator, i.e. N 4 value, we are able to remain under this threshold thanks to the cells sampling step value. Thus, despite the uncertainty of the perception and the localization, we are able to guarantee a functional domain for the vehicle.

In order to vary the sampling step, we applied a methodology to replay the scenario for each step (figure 4). To proceed, we took the results for a minimum step (step 1 of figure 4), which is typically used for the ground truth, and combined them to increase the step. For instance, with a sampling step of 0.2 m we combined information of two cells of 0.1 m by applying an aggregation rule that gives an identical result as if we had run the process in real time (step i of figure 4). The aggregation rule is the same as in table II. Then, to be able to compare the ground truth with the initial step, the aggregated cells are split back to the standard minimal sampling step of the LGM used for the ground truth as shown in figure 4.

V. EXPERIMENTS

A. Setup

In order to validate our method, we analyze experiments carried out with three vehicles in real road conditions. The experiments were conducted with two experimental Renault ZOE vehicles and one Renault Master of the Heudiasyc laboratory. An HD map previously constructed in the city of Compiègne has been used to lead the experiments. A NovAtel SPAN-CPT IMU with post-processed GNSS PPK corrections was used to have a centimeter-level accurate localization of each vehicle. A Velodyne VLP32-C LiDAR was used for the environment perception only in the ego vehicle. This sensor has a 360 • field of view with a theoretical range of 100 meters and runs at 10 hz. It enables to build several types of spatial occupancy representation.

1) Occupancy Grid (OG): With this representation, there are in the literature two methods to encode the space state: probabilistic grids [START_REF] Chen | Dynamic Environment Modeling with Gridmap: A Multiple-Object Tracking Application[END_REF] or evidential grids [START_REF] Mouhagir | Using evidential occupancy grid for vehicle trajectory planning under uncertainty with tentacles[END_REF]. As illustrated in figure 2, the intersection between each cell of the LGM and the OG is computed such that one cell of the LGM encompasses N cells of the OG. When merging two cells information, the O state is more conservative than the H state which is more cautious than the F state. Therefore, we consider the following aggregation rule, table II, when characterizing cells of the LGM. The resulting cell is occupied if at least one cell is occupied and is free only if both cells are free. In other cases, the resulting cell is hidden.

2) Polygons: Instead of a grid, an approach based on polygons can be used. One can compare it for instance to a parametric free space [START_REF] Schreier | Compact Representation of Dynamic Driving Environments for ADAS by Parametric Free Space and Dynamic Object Maps[END_REF]. The free space polygon F p considers that the space inside is free (illustrated with a purple contour in figure 5). The occupancy of objects are given under a list of polygons O p considered as being occupied. Thus, we can compute the intersection between each cell Cell i of the LGM and the different polygons. The results of these intersections provide the state S(•) of the cells:

if (O p ∩ Cell i ) = ∅ then S(Cell i ) = O ; else if (F p ∩ Cell i ) = Cell i then S(Cell i ) = F ; (4) else S(Cell i ) = H
In our approach, the polygon method has been implemented for its practicality. A simple geometric ground fitting based on Principal Component Analysis algorithm and a clustering algorithm was used to measure the free space and object surrounding the ego vehicle [START_REF] Zermas | Fast segmentation of 3D point clouds: A paradigm on LiDAR data for autonomous vehicle applications[END_REF]. A clustering has been used to extract objects clusters. Finally, the geometric layer of the HD map has been used to filter only objects on the road (figure 5). The ego-vehicle and paths of interest are extracted recursively from the HD map as they can be pre-processed thanks to the topological layer. The AOI are represented by polygons using the center line with its width. The Robot Operating System (ROS) middle-ware was used for the implementation and the Boost Geometry library for the geometric operations between areas of interest and polygons.

Cell i+1

Cell i F O H F F O H O O O O H H O H
A single dataset with sensors raw data has been recorded on an open road along a trajectory of 1.5 km in an urban area. The ego vehicle was following the two others all along the trajectory. That is why, during this scenario, only the clusters corresponding to the two perceived vehicles were used to compute the results.

When replaying data in real time, the LGM has been evaluated over 50 m in front of the ego vehicle as the two others were driving just in front of it (figure 5) and generated at a rate of 10hz. Firstly, the LGM has been recorded when replaying the dataset with a minimum sampling step of 0.1 m. In a second time, data were replayed in real time in order to simulate the same scenario but applying a Gaussian noise on the AV pose. It permits to add an error on the position with the same sensors and perception setup. We report six replays with standard deviation varying from 0 to 0.5 m. For each replay, there are two scenarios. The first one, without uncertainty propagation, the clusters size does not take into account the noise applied on the localization. The second one, with uncertainty propagation, a convex hull of each polygon object is augmented taking into account the localization uncertainty. We have used the method presented in [START_REF] Bernardi | High Integrity Lanelevel Occupancy Estimation of Road Obstacles Through LiDAR and HD Map Data Fusion[END_REF]. Each vertex of the polygon of an object generates several other vertices with a transformation that takes three times the standard deviation on the pose in each direction. Given all the resulting points, a new object convex hull is obtained and used in the characterization process.

In addition to that, given the fact that for each scenario we added simulations with a variation of the sampling step from 0.1 m to 5 m with a step of 0.1 m, it gives the equivalent of 6 × 2 × 50 = 600 simulations to study the influence of the sampling step, the noise and the uncertainty propagation on the perception polygons.

Figure 6 depicts values that are computed at the end of one scenario, replayed in real time and for a sampling step p s ∈ [0.1, 5] m. It shows cumulative sum of each indicator N i for several sampling steps of the LGM cells.

B. Results

Figure 7 depicts results obtained by the replay of data without propagating the uncertainty on the polygons of the detected vehicles. On the left of the figure, the F N R of the misleading occupation is plotted in function of the sampling step of the LGM. On the right, the F P R of the over cautious occupation is also plotted in function of the sampling step. Each curve shows the evolution of the rate value depending on the standard deviation of the noise added on the pose.

Figure 7a shows a decreasing level of the F N R when the sampling step raises. Indeed, the cells classified as occupied get bigger and, therefore, they cover a larger space that encompasses nearby cells that were incorrectly classified as free because of the localization error. As the standard deviation of the noise increases, the F N R decreases but at a higher level. This shows that it is possible to decrease the F N R under a given integrity risk by increasing the sampling step. For example, if we set a maximum integrity risk of 0.35% for the F N R, the integrity requirement can be reached by setting the sampling step to the value depicted by the red dots in figure 7. To compute these red dots, the different curves are interpolated in a logarithmic scale (dashed curves) and intersected with the required integrity risk. For a localization noise with a standard deviation of 0.5 m, we need a sampling step of around 3.5 m while for a standard deviation smaller than 0.2, a sampling step lower than 1 m is enough. The disadvantage of subsampling is that navigation quantities (such as distances to other vehicles, for example) may be too discrete to properly perform the current task. In other words, this graph shows that if we are not able to model the localization uncertainty, a correct sampling allows to keep the information safe as defined by the functional domain. For a given integrity risk, higher the uncertainty, higher the sampling step. The question can also be addressed the other way around: oversampling the LGM brings more misleading information since the curves increase as we get closer to zero.

On figure 7b, one can observe that the noise does not have influence on the F P R. As the aggregation rule is conservative on the occupied state, the F P R increases as the size of each occupied cell gets bigger. The superposition of each curve can be explained by the fact that, observing an occupied polygon, whatever the noise is, the size of the occupied polygon is not modified which does not change the impact on the cells.

Figure 8a shows the same curves as the previous figure with the difference that the uncertainty is now propagated onto the polygons which are enlarged. As a consequence, they cover a bigger space to ensure covering the correct occupied state. All the curves have the same shape with a much lower F N R. When the uncertainty is properly handled, we are therefore able to keep integrity at a given level for a shorter sampling step.

As it can be observed on figure 8b, the consequence of making the polygons bigger is a higher variation on the F P R depending on the standard deviation of the noise. It is explained by the fact that the perception system is much more conservative as occupied area gets bigger as localization uncertainty increases. When the noise model is estimated and propagated, subsampling leads to a loss of efficiency as the F P R grows rapidly for no variation of the F N R. When uncertainty is taken into account and propagated to the polygons, it is possible to decrease the F N R while the F P R increases with the noise.

For a given integrity risk, with uncertainty propagation, it is therefore possible to lower the sampling step compared to no propagation. It brings us closer to a finer representation. If localization uncertainty is not well defined, in order to keep a target integrity at 0.35% for the FNR, the sampling step parameter should be equal or below 3.6 m. However, if uncertainty is well bounded, the sampling step parameter can be lowered to 1m. Thus, we are able to refine information stored by the LGM. As a consequence, if the localization uncertainty is well defined, a lower sampling step is suitable in order not to supply misleading information. A low integrity risk level can be chosen. However, if the estimation error is not well bounded, we need to use a higher sampling step to keep the same integrity risk level.

VI. CONCLUSION

In this paper, an intermediate lane level information representation has been deepened with Lane Grid Maps. This representation aims at providing spatial information at a lane level to the maneuver planner of a decision-making module. We have demonstrated the capability of such a representation to help providing non-misleading information given an integrity risk. A LGM provides a spatial information representation in the road lane frame. This grid is sampled and each cell supplies spatial state information among: free, occupied or hidden. Based on its sampling step parameters, a confidence value of the information has been evaluated and compared to the ground truth of vehicles with real experiments. We demonstrated the need to find a proper sampling step for keeping information integrity under a given target. Indeed, oversampling can lead to increase misleading information, in addition to requesting heavier calculations to be done in real time. A proper sampling step is thus necessary depending on uncertainty management. The methodology we have presented in this article allows to answer this question.

In future work, situation predictions based on this representation will also be explored.

Fig. 1 :

 1 Fig. 1: Lane Grid Map representation. Areas Of Interest of the ego vehicle in blue are sampled. Cells are characterized by different states: Free (green), Occupied (red), Hidden (gray), out field of view (black). .

Fig. 2 :

 2 Fig. 2: Spatial information representation at each level of abstraction. The cells of the LGM are represented by the squared nodes.

s 2 wFig. 3 :

 23 Fig.3: Correspondence between evaluation parameters and the LGM of both the observation and the GT.

Fig. 5 :

 5 Fig.5: Lane Grid Map displayed in RViz with three vehicles. The ego vehicle (blue) is following the two other ones (gray and white). Each cell is characterized as one of the following states: free (green), occupied (red), hidden (white). The free space polygon is shown in purple and objects polygons in green.The global path is displayed at the bottom right.

Fig. 6 :

 6 Fig. 6: Cumulative percentages of the six N i indicators

Fig. 7 :

 7 Fig. 7: FNR and FPR of the LGM without uncertainty propagation on the polygons of the detected vehicles. Random errors with different standard deviations are added to the localization estimates. The x-axis corresponds to the spatial sampling of the LGM.

Fig. 8 :

 8 Fig.8: FNR and FPR of the LGM with uncertainty propagation on objects hulls as described in section IV-B.

TABLE I :

 I Parameters obtained with the combination of the ground truth and the observation.

TABLE II :

 II Aggregation rule of two OG cells. The information is condensed in the LGM cell.


Charly Girot
email: charly.girot@universite-paris-saclay.fr

Andreas Volk

Christine Walczak

Nathalie Lassau

Stéphanie Pitre-Champagnat

New Method for Quantification of Intratumoral Heterogeneity: A feasibility study on K trans maps from Preclinical DCE-MRI

Keywords: Tumor biomarker, Cluster Analysis, Magnetic Resonance Imaging, Angiogenesis, Treatment outcome 1. INTRODUCTION

Object: To develop new imaging biomarkers of therapeutic efficacy through the quantification of intratumoral microvascular heterogeneity.

Materials and Methods: The described method was a combination of non-supervised clustering and extraction of intratumoral complexity features (ICF): number of non-connected objects, volume fraction. It was applied to a set of 3D DCE-MRI K trans maps acquired previously on tumor bearing mice prior to and on day 4 of anti-angiogenic treatment. Evolutions of ICF were compared to conventional summary statistics (CSS) and to heterogeneity related whole tumor texture features (TF) on treated (n=9) and control (n=6) mice.

Results: Computed optimal number of clusters per tumor was 4. Several intratumoral features extracted from the clusters were able to monitor a therapy effect. Whereas no feature significantly changed for the control group, 6 features significantly changed for the treated group (4 ICF, 2 CSS). Among these, 5 also significantly differentiated the two groups. (3 ICF, 2 CSS). TF failed in demonstrating differences within and between the two groups.

Discussion: ICF are potential imaging biomarkers for anti-angiogenic therapy assessment. The presented method is expected to have advantages with respect to texture analysis based methods regarding interpretability of results and setup of standardized image analysis protocols.

Tumor heterogeneity has been identified as an important factor in need of evaluation for clinical diagnosis, prognosis, prediction of treatment efficacy and relapse of the disease. It is well known that tumor heterogeneity can be found at different levels: at population level, between patients with the same tumor type not responding equally to the treatment, which has triggered personalized therapeutic strategies over the last years; at patient level, between different lesions, or between primary tumor and metastasis that may have different characteristics [START_REF] Marusyk | Tumor heterogeneity: Causes and consequences[END_REF]; at tumor level, due to the presence of cell sub-populations with different genotypes and phenotypes that may differ in sensitivity to treatment [START_REF] Melo | Cancer heterogeneity-a multifaceted view[END_REF], as well as due to microvascular heterogeneity [START_REF] Rose | Quantifying spatial heterogeneity in dynamic contrast-enhanced MRI parameter maps[END_REF], this latter aspect being addressed in the present article.

Heterogeneity of tumor microvasculature is promoted by the neoangiogenic process which initializes and supports the creation of new blood vessels with abnormal structure and function, possibly leading to local tumor-cell hypoxia which in turn supports malignancy [START_REF] Jain | Antiangiogenesis Strategies Revisited: From Starving Tumors to Alleviating Hypoxia[END_REF]. Anti-angiogenic treatments modify blood vessel morphology, network, and function, therefore impacting tumor vascular distribution, which is related to treatment outcome. In a previous study, our team demonstrated that the evaluation of tumor vasculature with dynamic contrast-enhanced ultrasonography (DCE-US) provided a reliable imaging biomarker of anti-angiogenic treatments efficacy [START_REF] Lassau | Validation of Dynamic Contrast-Enhanced Ultrasound in Predicting Outcomes of Antiangiogenic Therapy for Solid Tumors: The French Multicenter Support for Innovative and Expensive Techniques Study[END_REF].

Non-invasive imaging is indeed a technique of choice to evaluate macroscopic aspects of intratumoral heterogeneity [START_REF] Connor | Cancer heterogeneity and imaging[END_REF]. It enables monitoring of the entire tumor and may guide or even replace invasive biopsies when using validated imaging biomarkers. Today, one of the radiological challenges is to identify and qualify new potential biomarkers based on vascular heterogeneity. These new biomarkers would reinforce the reliability of those currently used in clinical practice, mainly anatomic, and possibly provide pertinent predictive information earlier after therapy onset [START_REF] O'connor | Imaging biomarker roadmap for cancer studies[END_REF].

Two basic approaches are currently used to quantify intratumoral heterogeneity as an imaging biomarker. These two approaches consist in the extraction of features, on the one hand, from the whole tumor volume and, on the other hand, from various tumor sub-regions.

When extracted from the whole tumor volume, features describing image heterogeneity can be classified into two categories: Ia) The spatial gray-level methods (SGLM), which quantify spatial voxel distribution patterns based on intensity values, and Ib) Fractal features and geometry properties, which provide information on the ability of a shape to spread and occupy a defined surface or volume and inform on its complexity. For example, as a tool to depict microvasculature heterogeneity, category Ia) was used to compute texture features (TF) on dynamic contrast-enhanced MRI (DCE-MRI) parametric maps by Thibault et al., who showed these features to be potential biomarkers for early prediction in breast cancer [START_REF] Thibault | DCE-MRI Texture Features for Early Prediction of Breast Cancer Therapy Response[END_REF], and by Alic et al. for monitoring treatment response in limb sarcomas [START_REF] Alic | Heterogeneity in DCE-MRI parametric maps: a biomarker for treatment response?[END_REF]. Category Ib) was used by Rose et al. who employed the Rényi fractal-dimensions and geometric properties to describe the complexity of DCE-MRI parametric maps. They compared these features with conventional summary statistics (CSS) and demonstrated they were providing new additional information to discriminate between low and high-grade gliomas [START_REF] Rose | Quantifying spatial heterogeneity in dynamic contrast-enhanced MRI parameter maps[END_REF].

Concerning the sub-regions approach, tumor are subdivided, in a supervised or non-supervised manner, into multiple regions of interest (ROI) also called clusters. Features extracted from these clusters can be classified into two categories: IIa) CSS features for histogram description and IIb) features describing image organization, like the clusters' voxel count or volume fraction. As examples for category IIa), Pitre-Champagnat et al. matched histological measures of the mean vessel density within multiple ROI to the mean area under the curve extracted from DCE-US derived parametric maps [START_REF] Pitre-Champagnat | Dynamic Contrast-Enhanced Ultrasound Parametric Maps to Evaluate Intratumoral Vascularization[END_REF], and Longo et al. evaluated CSS from the pharmacokinetic model parameters K trans and vp in multiple clusters to assess tumor response to therapy [START_REF] Longo | Cluster analysis of quantitative parametric maps from DCE-MRI: application in evaluating heterogeneity of tumor response to antiangiogenic treatment[END_REF]. As an example for category IIb), Andersen et al. associated the volume fraction of clusters from K trans and ve maps with local tumor control [START_REF] Andersen | Pharmacokinetic analysis and k-means clustering of DCEMR images for radiotherapy outcome prediction of advanced cervical cancers[END_REF].

Here we present a new method for quantification of tumor microvasculature heterogeneity. This method is a combination of the previous approaches and is based on the analysis of tumor sub-regions with a specific focus on the quantification of the spatial organization complexity of these sub-regions, quantified by features we called intratumoral complexity features (ICF). Our method was applied to DCE-MRI derived 3D K trans parametric maps acquired on experimental tumors in mice during longitudinal follow-up with or without anti-angiogenic therapy. These parametric maps depict tissue perfusion and capillary permeability, and their heterogeneity is therefore related to microvascular heterogeneity. The dataset, considered as a test dataset for our method, had been previously acquired during a separate methodological study [START_REF] El | Anti-angiogenic therapy follow-up in a mouse tumor model by a novel 3D radial multi-gradient echo DCE MRI technique with individual AIF measurement[END_REF]. More specifically, the aim of the present study was threefold: I) to develop a new method for heterogeneity quantification based on ICF providing potential imaging biomarkers of therapeutic efficacy, II) to investigate if this method, applied to a set of DCE-MRI 3D K trans parametric maps, can differentiate the treated from the control group, III) to compare ICF with CSS as well as with SGLM based heterogeneity features extracted from the whole tumor.

MATERIALS AND METHODS

Our method is based on the combination of two successive image processing techniques, clustering and ICF extraction, and its schematic overview is given in Fig. 1. First (2.1), the clustering process was applied to parametric maps to highlight voxel clusters with similar K trans values, thus corresponding to tissues with similar perfusion and capillary permeability. Second (2.2), features were extracted from each cluster after binarization to quantify its complexity and thereby the spatial heterogeneity of corresponding tumor sub-regions. These features were independent of the voxels' K trans values. The method was retrospectively applied to data originating from a previous study [START_REF] El | Anti-angiogenic therapy follow-up in a mouse tumor model by a novel 3D radial multi-gradient echo DCE MRI technique with individual AIF measurement[END_REF]. These data correspond to K trans parametric maps of manually segmented entire tumor volumes, hereafter termed ROI. Main data characteristics and the MRI acquisition technique (described previously [START_REF] El | Anti-angiogenic therapy follow-up in a mouse tumor model by a novel 3D radial multi-gradient echo DCE MRI technique with individual AIF measurement[END_REF][START_REF] Vautier | Radial multigradient-echo DCE-MRI for 3D K trans mapping with individual arterial input function measurement in mouse tumor models: 3D Radial DCE-MRI in Mice with Individual AIF[END_REF]) are summarized in section 2.3.

Image Clustering

The K-means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] was used for the clustering process. K-means randomly selects k voxels as cluster centers, then assigns each voxel of the ROI to the cluster with lowest quadratic distance between voxel and cluster center values. After assigning all voxels within the ROI to a cluster, each cluster center is redefined as the mean value of all voxels belonging to the cluster. The process is repeated until convergence, which is reached when the center of each cluster does not change from one iteration to the next. Prior to the clustering process, the number of clusters k, which is the hyper-parameter for K-means, was optimized with the procedure mentioned in [START_REF] Andersen | Pharmacokinetic analysis and k-means clustering of DCEMR images for radiotherapy outcome prediction of advanced cervical cancers[END_REF] and described in [START_REF] Rokach | Clustering Methods[END_REF]. The optimization compares the mean intra-cluster distance (MICD) and the inter-cluster minimum distance (ICMD) after the clustering process. MICD and ICMD respectively decreases and increases when increasing the number of clusters. The optimal number of clusters is the one associated to the minimum value of the sum of MICD and ICMD normalized to their respective maximum value. This optimization process allows to select the number of clusters with the best compromise between a low intra-cluster variance and a large inter-cluster distance. It was run 30 times for each number of clusters from k=2 to 8 on every K trans map, and the final number of clusters was the one with the highest appearance rate among all maps. Finally, each parametric map was submitted 30 times to the clustering process, to avoid falling into local minima, and the resulting cluster centers were those with the maximum appearance rate among the 30 computations. The clustering process was performed in MATLAB (The MathWorks Inc., Natick, MA, USA).

Features extraction and selection

To describe for each cluster the complexity of spatial organization, heterogeneity related ICF were computed. ICF were selected based on the definition of spatial heterogeneity of a categorical map by Li and Reynolds, stating that heterogeneity can be defined as the complexity of the patches (i.e. objects) composing the categorical map, defined by their configuration and composition [START_REF] Li | On Definition and Quantification of Heterogeneity[END_REF]. Each of the k clusters was converted into a binary image (categorical map with two levels) and two ICF were computed to describe its complexity: the number of objects, computed with a 6-connected component analysis (corresponding to the voxel's faces), and the volume fraction of the cluster with respect to the volume of the ROI. The number of objects describes the spatial heterogeneity of a given cluster.

A high number of objects corresponds to a cluster split into several parts within the ROI, whereas a low number of objects indicates that the cluster is spatially gathered, thus representing low spatial heterogeneity. The volume fraction defines the contribution of the cluster's related heterogeneity to the overall degree of tumor heterogeneity. Four cases can be distinguished and are summarized in Table . 1. All in all, a total of k clusters × 2 features was extracted using a MATLAB algorithm. These ICF were compared to TF, i.e. second order statistical features, evaluated over the whole 3D tumor ROI. Amongst available features, our attention was focused on those being specifically related to spatial heterogeneity. The following 6 TF were extracted from the gray-level co-occurrence matrix (GLCM): Homogeneity, Energy, Contrast, Correlation, Entropy, Dissimilarity. LIFEx open access software was used for this analysis, and feature extraction process and mathematical definitions are described in more detail in the software documentation (http://www.lifex.soft.org , Orsay, France, [START_REF] Nioche | LIFEx: A Freeware for Radiomic Feature Calculation in Multimodality Imaging to Accelerate Advances in the Characterization of Tumor Heterogeneity[END_REF]). Feature extraction was performed with the following computational parameter settings: number of gray levels 256, absolute intensity rescaling with a maximum value of K trans = 0.8 min -1 corresponding to the overall maximum K trans value amongst all available parametric maps.

Moreover, CSS were used to describe the distributions of the K trans values both for 3D tumor ROI and within each cluster. Mean and median values, as well as the volume, were computed over the 3D tumor ROI and over the clusters. Mean and median highlights global differences in perfusion and permeability distribution patterns. The whole tumor volume was used as a common biomarker for treatment.

Finally, correlation between all individual features was assessed. They were considered to be correlated for Pearson correlation coefficient ≥ 0.95. In each group of correlated features only one was kept for further analysis.

Dataset properties

Investigated data originated from a previous study [START_REF] El | Anti-angiogenic therapy follow-up in a mouse tumor model by a novel 3D radial multi-gradient echo DCE MRI technique with individual AIF measurement[END_REF]. Main properties of the dataset and the data acquisition protocol are summarized hereafter. Data had been acquired on 15 male swiss nude mice carrying a human colorectal tumor xenograft (TC302, Institut Curie, France) subcutaneously implanted at the abdominal level. Animals had been treated for 10 days by single oral administration of Sunitinib® (50 mg/kg, n=9) or vehicle only solution (n=6). Mice underwent DCE-MRI immediately prior to (day 1) and at day 4 after therapy initiation. DCE-MRI experiments had been carried out on a 4.7 T small animal MRI system (Bruker, Biospin, Germany) using the previously described radial 3D DCE-MRI technique RETIA [START_REF] Vautier | Radial multigradient-echo DCE-MRI for 3D K trans mapping with individual arterial input function measurement in mouse tumor models: 3D Radial DCE-MRI in Mice with Individual AIF[END_REF], with a home built quadrature birdcage probe (D=35 mm, l=50 mm). Animals had been anesthetized with isoflurane (AErrane®, Baxter, France) and the body temperature had been stabilized. For localization purpose, pre-contrast high-resolution images had been acquired with the same 3D FOV as for radial DCE-MRI covering the whole tumor (respiratory triggered multi-slice spin echo sequence TR~1s, TE = 10.7 ms, matrix = 128 × 128, slice thickness = 0.47 mm, 64 slices). An experimental medium molecular weight contrast agent P846 (3.5 kDa, Guerbet, France) had been injected at a dose of 0.089 mmol/kg (injection rate 600 µL/min, injected volume 200 µL). During DCE-MRI acquisition, 15 3D datasets (tumor) and 960 2D datasets (Arterial Input Function (AIF) measured on a single slice positioned on the heart) had been acquired continuously and simultaneously including one 3D pre-contrast image (64 2D pre contrast images). Main 3D parameters were TE1 = 0.9 ms, ∆TE = 1.11 ms, 10 echoes, matrix = 64×64×64, FOV = 30 mm × 30 mm × 30 mm, 64 readout points, 4096 projections, temporal resolution = 2 min. Main 2D parameters were TE1 = 1.9 ms, ∆TE = 1.11 ms, 10 echoes, matrix = 64×64, FOV = 30 mm × 30 mm, slice thickness = 2 mm, 64 readout points, 64 projections, temporal resolution = 2 s, TR = 31 ms. Images had been reconstructed using a C++ standard regridding algorithm [START_REF] Vautier | Radial multigradient-echo DCE-MRI for 3D K trans mapping with individual arterial input function measurement in mouse tumor models: 3D Radial DCE-MRI in Mice with Individual AIF[END_REF][START_REF] Jackson | Selection of a convolution function for Fourier inversion using gridding (computerised tomography application)[END_REF]. Sliding window reconstruction had been performed for 3D data providing a virtual temporal resolution of 30 seconds. CA concentrations in voxels (tumor) or AIF-ROI (heart) had been estimated from R2* corrected R1(t) using r1 = 15 s -1 mM -1 . AIF time constants had been measured by fitting a biexponential decay to the concentration time-curve. K trans values had been computed voxel wise using the Tofts-Kermode pharmacokinetic model [START_REF] Tofts | Measurement of the blood-brain barrier permeability and leakage space using dynamic MR imaging. 1. Fundamental concepts[END_REF]. The process of reconstructing K trans maps produced not-a-number (NaN) voxels. NaN voxels were generated when fitting results did not fall in a physical range or in the case of fit failures due to poor signal to noise ratio and/or lack of R1 increase [START_REF] Vautier | Radial multigradient-echo DCE-MRI for 3D K trans mapping with individual arterial input function measurement in mouse tumor models: 3D Radial DCE-MRI in Mice with Individual AIF[END_REF]. The mean fraction of NaN voxels was 23.8 ± 13.2 %. They were excluded from further analysis. The entire tumor ROI were manually segmented on the high-resolution spin echo images and then transferred to the K trans maps. Previous CSS analyses highlighted a significant decrease of median K trans between day 1 and day 4 in the treated group [START_REF] El | Anti-angiogenic therapy follow-up in a mouse tumor model by a novel 3D radial multi-gradient echo DCE MRI technique with individual AIF measurement[END_REF].

Statistical Analysis

Selected features were analyzed individually to characterize their behavior regarding three comparisons: (i) evolution within the control group between the examinations at day 1 and day 4; (ii) evolution within the treated group between the examinations pre-and post-treatment; (iii) evolution differences between the two groups. For (i) and (ii), the feature values pre-(day 1) and post-treatment (day 4) were analyzed using a Wilcoxon paired test. For (iii), feature differences were computed between post-and pre-treatment, then analyzed using a Wilcoxon unpaired test. Resulting p-values ≤ 0.05 were considered demonstrating statistically significant differences. A correction for multiple comparisons was applied to all selected features using the Benjamini-Hochberg procedure. The false discovery rate (FDR) was set to 0.05. Analyses were conducted with R software (R Core Team, 2019).

RESULTS

Clustering

The optimization process resulted in an optimal number of k = 4 clusters per parametric map. Table 2 represents the normalized sums of normalized MICD and ICMD and shows that 4 clusters were the best compromise between low intra-cluster variance and high inter-cluster distance. Clusters were numbered #1 to #4 in ascending order of mean K trans with mean centroids over the whole dataset of 0.8, 2.2, 4.4 and 9.9 (10 -2 min -1 ). Fig. 2 shows the clustering result for a single tumor slice. It illustrates how the clusters were used, either as binary images to compute ICF, or as masks for selecting sub-populations of voxels and extract CSS.

Visual inspection of heterogeneity evolution

Visual inspection of 2D maps allowed qualitative evaluation of heterogeneity evolution in the images. As illustrated in Fig. 3 on representative slices, tumor K trans heterogeneity evolves over the 4 days follow-up period. While heterogeneity appears to decrease in the treated mouse, with an evolution from clusters with intermediate (#2 and #3) towards low K trans values (#1), this is not the case for the control mouse, for which cluster organization complexity increases. Our method intended to quantify these spatial heterogeneity modifications within the clusters through ICF extraction.

Feature extraction and selection

A total of 29 features were computed from each parametric map: 8 ICF were computed with our method after the clustering process; 15 CSS were extracted from the clusters (4 clusters × 3 features) and from the whole tumor ROI (3 features); 6 GLCM were extracted from the whole tumor ROI.

Feature correlations were evaluated to define a limited number of independent features. None of the ICF were correlated with each other. Means and medians were correlated for each cluster and for the ROI, therefore, only the medians were kept, except the median of cluster #1 which was correlated to the ROI median. Amongst GLCM features, contrast and entropy were respectively correlated to dissimilarity and homogeneity and therefore not selected.

As a result, a total of 21 features were kept for further analysis: 8 ICF (number of objects and volume fraction of each cluster), 4 GLCM features (homogeneity, energy, correlation, and dissimilarity), and 9 CSS (median and volume of each cluster and of the ROI -except median of cluster #1). We checked that these 21 features were not correlated with each other.

Statistical Analysis

Wilcoxon tests were performed individually on the 21 features and results with p-value ≤ 0.1 are summarized in Tables 3 to 5. Tables 3 and4 report mean and standard error values respectively for the control group at day 1 and day 4 and for the treated group for pre-and post-treatment examinations. Statistical significance of the difference between day 1 and day 4 values as well as the adjusted p-values (q-values) after Benjamini-Hochberg correction are also reported. Table 5 reports the delta features (post -pre-treatment) for control and treated group and statistical analysis between these groups. In Tables 3 to 5, features that remained significant after correction with FDR=0.05, are labelled with "*". Feature values were similar for the control and the treated group before treatment initiation with no significant differences. Also, within the control group, no feature demonstrated significant differences between day 1 and day 4 (Table 3).

Amongst the ICF evaluated at cluster-level, after correction for multiple comparisons, 4 demonstrated significant differences between pre-and post-treatment evaluations within the treated group (Table 4): number of objects in cluster #3 (p=0.0039) and volume fractions in clusters #1 (p=0.0039), #2 (p=0.0039), #3 (p=0.0117). 3 demonstrated significant differences in their delta-values between the two groups (Table 5): delta-number of objects in cluster #3 (p=0.0016), and delta-volume fractions in clusters #1 (p=0.0048), #2 (p=0.0008). Amongst the 4 GLCM features, after correction for multiple comparisons none presented significant differences within the treated group, neither between both groups. Amongst the 9 CSS, after correction for multiple comparisons, 3 demonstrated significant differences within the treated group (Table 4): volumes in clusters #2 (p=0.0039), #3 (p=0.0039) and median in ROI (p=0.0078). 3 demonstrated significant differences between the two groups (Table 5): delta volumes in clusters #2 (p=0.0004), #3 (p=0.012) and delta volume in ROI (p=0.0028). Features from Table 5 that were statistically discriminant after correction between both groups are represented as boxplot in Fig. 4.

DISCUSSION

In oncology, a growing interest for targeted therapies and particularly for anti-angiogenic drugs led to a growing number of studies assessing tumor microvasculature and vascular heterogeneity due to the ability of these therapies to normalize and reorganize microvascular structure. Nonetheless, in clinical practice, tumors are still assessed using morphological criteria without considering the modification of heterogeneity patterns which may arise earlier during anti-angiogenic treatments. Some studies tried to exploit tumor heterogeneity using texture analysis to assess therapeutic response [START_REF] Fave | Delta-radiomics features for the prediction of patient outcomes in non-small cell lung cancer[END_REF][START_REF] Liu | Use of texture analysis based on contrast-enhanced MRI to predict treatment response to chemoradiotherapy in nasopharyngeal carcinoma: MRI in Treatment Response Prediction[END_REF][START_REF] Panzeri | Prediction of Chemoresistance in Women Undergoing Neo-Adjuvant Chemotherapy for Locally Advanced Breast Cancer: Volumetric Analysis of First-Order Textural Features Extracted from Multiparametric MRI[END_REF], but none of these approaches succeeded yet to reach use in clinical practice, mainly because of a lack of standardization and guidelines [START_REF] Lambin | Radiomics: the bridge between medical imaging and personalized medicine[END_REF].

In this study, we developed a new method for assessing intratumoral heterogeneity, and we described a first application to the longitudinal follow-up of microvascular heterogeneity on an experimental tumor model. The method associates a clustering process and the extraction of features, we called ICF, to identify and quantify different levels of complexity of intratumoral heterogeneity. The development of this method was guided by the following strategy: First, processed images should be parametric maps. This allows to directly interpret heterogeneity of the parametric maps in terms of physiological heterogeneity. Second, a non-supervised clustering process was applied to these parametric maps. This procedure allows to subdivide the tumor into several functional sub-regions characterized by the measured physiological parameter levels. It has indeed been suggested that, in some clinical scenarios, measurements performed in tumor sub-regions may provide more useful biomarkers than measurements averaged over the whole tumor [START_REF] O'connor | Imaging Intratumor Heterogeneity: Role in Therapy Response, Resistance, and Clinical Outcome[END_REF]. Third, the features used to quantify heterogeneity should provide a more intuitive spatial representation (qualitatively cross-checkable by visual inspection of 2D images) than second order statistical analysis approaches in texture analysis.

According to this strategy, the method was applied to a set of 3D DCE-MRI K trans parametric maps, previously acquired in mice during a longitudinal study. Heterogeneity related features extracted from DCE-MRI parametric maps are of interest, as they should monitor treatment effects by assessing microvasculature modifications. Information from these features is expected to be particularly useful in the absence of morphological changes. In this study, heterogeneity was assessed by ICF. They were compared to GLCM features describing heterogeneity over the whole tumor ROI and to classical non-spatial CSS features.

The ability of ICF, GLCM, and CSS features to monitor therapeutic efficacy was evaluated by the following comparisons: i) day 1 vs day 4 within the control group, ii) pre-(day 1) vs post-treatment (day 4) within the treated group, and iii) evolution differences for the control vs the treated group. It was checked that all features had similar values in the treated and control groups at day 1 confirming similar initial status of both groups.

Statistical analysis showed that several intratumoral features extracted from the clusters were able to monitor a therapy effect. Whereas no feature significantly changed for the control group (comparison (i)), 6 features significantly changed for the treated group (4 ICF, 2 CSS -comparison (ii)). Among these, 5 also significantly differentiated the two groups. (3 ICF, 2 CSS -comparison (iii)).

On the contrary, heterogeneity related GLCM features extracted from the whole tumor ROI didn't change significantly for either comparison. Therefore, the presented method seems to be more sensitive than classical texture analysis for monitoring therapy effects under our experimental conditions. Median K trans measured over the whole tumor ROI decreased for the treated group (comparison (ii)), as previously described [START_REF] El | Anti-angiogenic therapy follow-up in a mouse tumor model by a novel 3D radial multi-gradient echo DCE MRI technique with individual AIF measurement[END_REF]. This result is compatible with feature evolution at the cluster level: an increase of the volume fraction in cluster #1, representing the lowest K trans values, associated to a decrease of the volume fractions in clusters #2 and #3, representing intermediate K trans values. These findings indicate a shift of tumor microvasculature characteristics to lower perfusion and capillary permeability, as may be expected under antiangiogenic therapy.

It has to be noted that, even though the tumor volume tends to decrease for the treated group at day 4 of treatment, this change is not significant in our study. Whole tumor volume is generally used as a reference biomarker in oncology. However, in clinical practice during anti-angiogenic treatment, morphological evaluation of tumors based on the RECIST criterion has been shown to significantly underestimate tumor response to treatment [START_REF] Smith | Assessing Tumor Response and Detecting Recurrence in Metastatic Renal Cell Carcinoma on Targeted Therapy: Importance of Size and Attenuation on Contrast-Enhanced CT[END_REF]. For example, some studies showed that no tumor size modification was observed during the course of treatment for respectively 50 to 75% of patients [START_REF] Motzer | Nivolumab plus Ipilimumab versus Sunitinib in Advanced Renal-Cell Carcinoma[END_REF][START_REF] Luo | Early evaluation of sunitinib for the treatment of advanced gastroenteropancreatic neuroendocrine neoplasms via CT imaging: RECIST 1.1 or Choi Criteria?[END_REF]. This is mainly because vascular modifications often occur before any size reduction of the tumor. Our preclinical results are in line with these clinical findings.

As mentioned above, concerning heterogeneity description, in our study only ICF were demonstrating significant differences for comparisons (ii) and (iii) after BH correction. More precisely, 2 GLCM features with p < 0.05 (Energy, comparison (ii); Correlation, comparison (iii)) were discarded after BH correction. In our study BH correction was applied with FDR = 0.05 in order to discard false positives in the case of multi-testing. Nevertheless, the use of a correction for multi-testing is still a topic of discussion as some authors advise not to use any correction when a study is exploratory with unplanned comparisons being considered as potential ones for further investigation [START_REF] Armstrong | When to use the Bonferroni correction[END_REF]. Further studies will be necessary to compare (or to combine) ICF and GLCM features as potential imaging biomarkers related to tumor microvascular heterogeneity.

Evolution of ICF can be relatively easily interpreted in terms of microvascular properties. Indeed, for comparison (ii) within the treated group, observed ICF evolutions can be considered indicative of K trans map homogenization: Volume fraction for cluster #1, which had the smallest and also constant number of objects amongst clusters #1 to #3, increased by 37%, and the number of objects in cluster #3 was reduced by 43 %. Volume fraction of cluster #2, the most heterogeneous one, decreased by 38% (100-100*22.5/35.9) with a stable number of objects. Thus, by considering clusters #1 and #3, a total of 76.1% of the analyzed tumor volume can be considered more homogenous in terms of perfusion and capillary permeability with lower K trans at day 4 post-treatment, which, again, may be expected under anti-angiogenic therapy. Furthermore, it is noteworthy that, for the control group, while no ICF demonstrated a significant evolution between day 1 and day 4, numbers of objects in all clusters tended to increase with volume fractions staying nearly identical. This observation is compatible with an overall increase of microvascular heterogeneity, which may be expected in growing tumors with neo-angiogenic activity.

Results also indicated that none of the significant features were found within the cluster #4 with highest K trans values, which may be considered to be "hotspots". This is in contrast with previous studies which focused on such "hotspots" considering they were bearing more information than values averaged over the whole tumor. For example, a greater range of K trans changes was identified by Hayes et al. in these "hotspots" when compared to the tumor ROI [START_REF] Hayes | Assessing changes in tumour vascular function using dynamic contrast-enhanced magnetic resonance imaging[END_REF]. The difference with respect to our results may be explained by the definition of the "hotspots". While they were defined as the 75 th percentile of the tumor histogram by Hayes et al., our non-supervised clustering approach resulted in "hotspot" clusters with few voxels and volume fractions around only 2%.

From a technical point of view, the main characteristics and advantages of the presented approach appear to be the following:

Our approach quantifies heterogeneity by the description of tumor sub-regions complexity, related, in our study, to K trans map heterogeneity. This is in contrast with previous studies who assessed either the complexity of the tumor as a whole [START_REF] Rose | Quantifying spatial heterogeneity in dynamic contrast-enhanced MRI parameter maps[END_REF], or only the histogram distribution in tumor sub-regions without considering spatial distribution patterns [START_REF] Longo | Cluster analysis of quantitative parametric maps from DCE-MRI: application in evaluating heterogeneity of tumor response to antiangiogenic treatment[END_REF].

Moreover, a potential specific advantage comes from the use of ICF extracted from binary images. As they do not depend on absolute voxel values after the non-supervised clustering process, they should be less sensitive to acquisition related variations or bias in parameter quantitation (e.g. due to poor AIF quality or pre-contrast T1 mapping in the case of DCE-MRI). This should be of advantage for the setup of multicentric studies where measurement results need to be independent of the MRI scanner and protocol. Furthermore, in this context, it is noteworthy that the only user defined computational parameter for the proposed method is the "connectivity" parameter set during computation of object numbers within the clusters. We checked, that, when varying connectivity from 6 (all voxel faces), as used in the study, to 18 (6 faces and 12 edges), results still discriminated both groups of mice with the same significant features. However, this was no longer the case for the maximal 3D connectivity of 26 (6 faces, 12 edges and 8 vertices). Indeed, in this latter case, all voxels are merged into few objects thus no longer capturing image heterogeneity. Consequently, this connectivity should be avoided using our method.

On the contrary, in the case of texture analysis, TF depend on several user-defined computational parameters and settings: number of bins for discretization, size of bins, absolute (user defined min, max values) or relative (ROI defined min, max values) intensity rescaling. These computational parameters and settings may also be handled differently by different software, making standardization a more complicated task [START_REF] Traverso | Repeatability and Reproducibility of Radiomic Features: A Systematic Review[END_REF].

Consequently, when analyzing data acquired with different MRI systems and protocols, our approach may be expected to be more robust in comparison to texture analysis. This aspect will be addressed in future studies.

In conclusion, in this preclinical feasibility study, we presented a new method for image heterogeneity analysis, and we showed that, when applied to DCE-MRI derived K trans parametric maps, this method may provide potential imaging biomarkers for anti-angiogenic therapy assessment in agreement with the consensus on Imaging Biomarkers [START_REF] O'connor | Imaging biomarker roadmap for cancer studies[END_REF]. The method is expected to have advantages with respect to texture analysis based methods regarding interpretability of results and setup of standardized image analysis protocols. Further studies applying this method to different parametric maps across multiple imaging modalities will be needed to evaluate its usefulness in a clinical setting for treatment monitoring or, more generally, tissue characterization.

Number of Objects

LowTable 1 :Table 2 :Table 3 : 4 :Table 5 :Fig. 1 :Fig. 2 :Fig. 3 :Figure 4 :

 123451234 Fig. 1: Schematic overview of the method proposed in this study to investigate intratumoral heterogeneity from parametric maps. The method is a combination of clustering and extraction of complexity features (ICF). This method was compared to heterogeneity features from the gray-level co-occurrence matrix (GLCM) and to conventional summary statistics (CSS)

Acknowledgements

For MRI data acquisition the tumor model had been provided by Institut Curie (Preclinical Investigation Laboratory, LIP); Sunitinib® had been provided by Pfizer; MRI contrast agent had been provided by Guerbet.

Authors' Contributions

Charly Girot contributed to study conception and design, analysis and interpretation of data, drafting of manuscript and critical revisions. Andreas Volk contributed to data acquisition, analysis and interpretation of data and critical revision. Stéphanie Pitre-Champagnat contributed to analysis and interpretation of data and critical revision. Christine Walczak contributed to acquisition and analysis of data. Nathalie Lassau contributed to critical revision.

Compliance with Ethical Standards

Funding: the authors did not receive support from any organization for the submitted work.

Conflict of Interest: all authors certify that they have no affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in the manuscript.

Ethical approval: Care and use of animals had been in compliance with national and institutional regulations and guidelines