Céline Fouard

Grégoire Malandain

Systematized calculation of optimal coefficients of 3-D chamfer norms

Keywords:

Chamfer distances are widely used in image analysis, and many ways have been investigated to compute optimal chamfer mask coefficients. Unfortunately, these methods are not systematized: they have to be conducted manually for every mask size or image anisotropy. Since image acquisition (e.g. medical imaging) can lead to anisotropic discrete grids with unpredictable anisotropy value, automated calculation of chamfer mask coefficients becomes mandatory for efficient distance map computation. This article presents a systematized calculation of these coefficients based on the automatic construction of chamfer masks of any size associated with a triangulation that allows to derive analytically the relative error with respect to the Euclidean distance, in any 3-D anisotropic lattice.

Introduction

Distance transformations (DTs) are widely used in image analysis since they allow to recover morphometric features of a binary shape. Among other applications, they can be applied to skeleton computation [START_REF] Pudney | Distance-ordered homotopic thinning: A skeletonization algorithm for 3d digital images[END_REF], Voronoï diagram construction, or shape-based interpolation [START_REF] Herman | Shape-based interpolation[END_REF]. Distance transformation transforms a binary image into a grey level image where the value of each foreground pixel corresponds to its shortest distance to the background. Brute-force computation of DT is not compatible with expected image analysis requirements, so DTs are usually computed by propagation. Exact Euclidean maps can be computed through Euclidean Distance Transformations (EDT). Several EDT have been proposed, using morphological operators [START_REF] Shih | A mathematical morphology approach to euclidean distance transformation[END_REF][START_REF] Huang | A euclidean distance transform using grayscale morphology decomposition[END_REF], filters [START_REF] Saito | New algorithms for euclidean distance transformation of an n-dimensional digitized picture with applications[END_REF], several path on rows and columns [START_REF] Hirata | A unified linear-time algorithm for computing distance maps[END_REF], or propagating vectors [START_REF] Danielsson | Euclidean distance mapping[END_REF][START_REF] Ragnemalm | The euclidean distance transform in arbitrary dimensions[END_REF], but lead to time and/or memory consuming algorithms. A good trade-off between precision and computational cost for DT is achieved by chamfer maps that have been made popular by Borgefors [START_REF] Borgefors | Distance transformations in digital images[END_REF]. These maps are computed through two raster-scan on the image that propagate the distance values by the way of chamfer masks. The coefficients of the mask are (proportional) estimation of short-range distances: the larger the chamfer mask is, the closest to the Euclidean map the chamfer map will be. The calculation of optimal coefficients can be done by minimizing either an absolute error [START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF] or a relative one [START_REF] Verwer | Local distances for distance transformations in two and three dimensions[END_REF]. It has first been done for 2-D 3 × 3 masks [START_REF] Borgefors | Distance transformations in arbitrary dimensions[END_REF] in isotropic lattices, then extended to larger masks [START_REF] Borgefors | Distance transformations in digital images[END_REF][START_REF] Verwer | Local distances for distance transformations in two and three dimensions[END_REF] and to higher dimensions [START_REF] Borgefors | On digital distance transforms in three dimensions[END_REF]. Anisotropic lattices have also been considered [START_REF] Coquin | Discrete distance operator on rectangular grids[END_REF][START_REF] Mangin | Chamfer distances in anisotropic 3D images[END_REF][START_REF] Sintorn | Weighted distance transfoms for images using elongated voxel grids[END_REF]. However, those calculations remain tedious and are not systematized: thus they have to be conducted manually for every mask size or anisotropy value.

Our motivation is the computation of DT in 3-D medical images: they are usually acquired on anisotropic lattices (slice thickness is usually larger than the pixel size) and this anisotropy may vary from one acquisition to the other. The efficient computation of chamfer maps requires then the calculation of the chamfer mask's coefficient to be automated. calculation of these coefficients for any mask size and any anisotropy value. In addition to classical error criteria, we also consider norm constraints [START_REF] Remy | Optimizing 3d chamfer masks with norm constraints[END_REF] that guarantee predictable results. Our approach is based on the automatic construction of chamfer masks of any size associated with a triangulation that allows to derive analytically the relative error with respect to the Euclidean distance.

In the following, we first recall some basic definitions. Then we describe error estimation and norm constraints. Some results (coefficients of isotropic 7 3 and anisotropic 3 3 masks) are given before we conclude.

Definitions and notations

We recall here some notations and definitions. We consider the discrete space E = Z 3 . An image I is an application defined on E.

A discrete distance is an application d : E × E -→ N that verifies for all p, q, r ∈ E the following 4 properties:

d(p) ≥ 0,
d(p, q) = 0 ⇐⇒ p = q, d(p, q) = d(q, p), and d(p, q) ≤ d(p, r) + d(r, q).

Given a discrete distance d, the application n :

E -→ N is a discrete norm on Z if and only if ∀p ∈ E, n(λp) = |λ|n(p) ∀λ ∈ Z.
Let us consider a binary image I with foreground X and background X. The distance map D X is an application defined on E such that D X (p) = inf q∈X d(p, q). Distance maps can be approximated by chamfer maps, that can be computed with a two-passes (so-called forward and backward passes) algorithm [START_REF] Rosenfeld | Sequential operations in digital picture processing[END_REF]. To do so, we need to define the chamfer mask which is a set M C = {(v i , ω i ), 1 ≤ i ≤ m} of weighted vectors representing authorized displacements. It is centered in O, symmetrical with respect to its center and contains at least a base of E.

Given a chamfer mask M C and two points p, q ∈ E, we define a path P pq from A to B as a sequence of vectors v i ∈ M C so that -→ pq is expressed as a linear combination of vectors: -→ pq = n i .v i with n i ∈ N. The cost W of a path P pq is defined as W (P pq ) = n i .ω i . The chamfer distance between two points p, q ∈ E as the minimal possible cost, i.e. d C (p, q) = min Ppq W (P pq ).

Computing optimal coefficients for chamfer norms

Calculating the optimal weights for a given chamfer mask is usually achieved by minimizing the error (either absolute or relative) between the chamfer's distance and the Euclidean one. Thanks to symmetry considerations, we can only consider the mask generator M g C , i.e. the part of the anisotropic chamfer mask M C that is included in the first eighth of the space, 1 8 Z 3 , delimited by the half-lines (O, x), (O, y) and (O, z). Moreover, vectors are also chosen so that ∀i, j with i = j, ∃n ∈ N such that v i = nv j .

Estimating the error between a chamfer distance and the Euclidean one is quite awkward when dealing with large masks. This difficulty can be reduced if we are able to triangulate the mask generator M g C into regular cones. A continuous cone, defined by a triplet of vectors and denoted by v i , v j , v k , represents the region of R 3 delimited by the vectors v i , v j and v k , i.e.

v i , v j , v k {M ∈ R 3 | --→ OM = λ i .v i + λ j .v j + λ k .v k , λ i , λ j , λ k ∈ R + } A discrete cone, v i , v j , v k , is the set of points of Z 3 that are included in the continuous cone v i , v j , v k : v i , v j , v k {M ∈ Z 3 |M ∈ v i , v j , v k }. A regular cone is a discrete cone that verifies ∆ i,j,k = ±1 where ∆ i,j,k is the determinant of the matrix |v i v j v k | (first column is vector v i , etc
). Regular cones have the interesting properties that any point of the cone can be expressed in the basis of the 3 vectors defining the cone, i.e.

v i , v j , v k = {M ∈ Z 3 | --→ OM = λ i .v i + λ j .v j + λ k .v k , λ i , λ j , λ k ∈ N}
only holds for regular cones [START_REF] Hardy | An Introduction to the Theory of Numbers[END_REF]. Having a mask generator that can be triangulated into regular cones allows us to reduce the calculation of the error into independent calculations in each regular cone. In the following, we will only deal with such mask generators. To ensure that they can be triangulated into regular cones, we build them with the Farey triangulation [START_REF] Remy | Optimizing 3d chamfer masks with norm constraints[END_REF]. This technique allows us to recursively (and automatically) built large mask generators M g C with their associated regular triangulation

T g C = { v i , v j , v k } (see appendix A).

Error definition and calculation

We have chosen to minimize the relative error between the chamfer distance and the Euclidean one, computed on planes x = cste, or y = cste, or z = cste. Let us consider a point P = (x, y, z). According to above section, we know that its chamfer distance to the origin O, denoted d C (P ), is a linear combination of the weights of the 3 vectors defining the regular cone it belongs to. We have then d

C (P ) = a.ω i + b.ω j + c.ω k with --→ OP = a.v i + b.v j + c.v k .
Solving the latter expression yields (recall that ∆ i,j,k = ±1 for regular cones)

a = 1 ∆ i,j,k x x j x k y y j y k z z j z k , b = 1 ∆ i,j,k x i x x k y i y y k z i z z k , and c = 1 ∆ i,j,k x i x j x y i y j y z i z j z and allows us to obtain d C (x, y, z) = α.x + β.y + γ.z with α = (y j z k -y k z j ).ω i + (y k z i -y i z k ).ω j + (y i z j -y j z i ).ω k β = (z j x k -z k x j ).ω i + (z k x i -z i x k ).ω j + (z i x j -z j x i ).ω k γ = (x j y k -x k y j ).ω i + (x k y i -x i y k ).ω j + (x i y j -x j y i ).ω k
Chamfer distances are usually computed with integer weights, and have to be scaled with a real factor, ε (typically the displacement associated with the smallest voxel size), to be compared to the Euclidean distance d E (x, y, z) =

d 2 x x 2 + d 2 y y 2 + d 2 z z 2
where d x , d y and d z denote the voxel size in the x, y, and z direction. The relative error to minimize is then defined by

E relative (x, y, z) = d C /ε -d E d E = 1 ε α.x + β.y + γ.z d 2 x x 2 + d 2 y y 2 + d 2 z z 2 -1. (1)
Depending on the orientation of the cone, this error has to be minimized on either the plane x = M , or y = M , or z = M . Without loss of generality, we will only go into details for the case x = M ,M = 0. The error has then to be estimated on the triangle defined as the intersection between the cone v i , v j , v k and the plan x = M . The vertices of this triangle are the points

V l = (M, M y l x l , M y l x l ) = (M, y ′ l , x ′ l ) for l = i, j, k.
In plane x = M , the relative error can be rewritten

E x (y ′ , z ′ ) = 1 ε α + β.y ′ + γ.z ′ d 2 x + d 2 y y ′ 2 + d 2 z z ′ 2 -1 with y ′ = y M and z ′ = z M . ( 2
)
E x is continuous on a closed and bounded interval (the triangle V i V j V k ), it is then bounded and reaches its bounds. Its extrema can be located either inside the triangle, or on the edges of the triangle, or at the vertices of the triangle. Let us consider the three cases.

1. The extremum is inside the triangle. By derivating equation 2, it comes that the extremum will be located at (y ′ max , z ′ max ) = . If this point is inside the triangle V i V j V k , it yields an extreme value

E x (V i V j V k ) = 1 ε α 2 d 2 x + β 2 d 2 y + γ 2 d 2 z -1.
2. The extremum is on an edge. There are three edges, but we will only present the calculation for V i V j . In this case, a point M belonging to the edge can be represented by M = aV i + (1 -a)V j yielding the relative error along the edge

E x (a) = 1 ε (β.Y +γ.Z)a+(α+β.yj +γ.zj ) √ (d 2 y Y 2 +d 2 z Z 2 )a 2 +2(d 2 y yj Y +d 2 z zj Z)a+d 2 x +d 2 y y 2 j +d 2 z z 2 j -1 with Y = y i -y j Z = z i -z j .
(3) After derivation, it can be shown that the extreme value is reached for

a max = - (β(yj Z-zj Y )+αZ)zj d 2 z +(γ(zj Y -yj Z)+αY )yj d 2 y -(βY +γZ)d 2 x (β(yj Z-zj Y )+αZ)Zd 2 z +(γ(zj Y -yj Z)+αY )Y d 2 y .
If 0 ≤ a max ≤ 1, the extreme value E x (V i V j ) is given by E x (a max ) whose form is not simple enough to be displayed here.

3. The extremum is reached on one of the triangle's vertices, V l where the relative error value is given by

E x (V l ) = 1 ε ω l ||v l || R -1 with ||v l || R = x 2 l d 2 x + y 2 l d 2 y + z 2 l d 2 z
Thus we are now able to compute both the minimum and maximum relative errors, τ min and τ max (please recall that they depend on ε), for a mask generator M g C by

τ min = min vi,vj ,vk ∈T g C min E u (V i V j V k ), min l,m E u (V l V m ), min l E u (V l ) τ max = max vi,vj ,vk ∈T g C max E u (V i V j V k ), max l,m E u (V l V m ), max l E u (V l )
u denoting the plane (x, y, or z) where the error is minimized; E u (V i V j V k ) and E u (V l V m ) being estimated only if the corresponding extremum lies in the correct interval.

The global relative error is defined by τ (ε) = max(τ min (ε), τ max (ε)). According that τ min (ε) < 0 and τ max (ε) > 0, we can make them equal in absolute value by changing the value of ε into ε opt = ε τmin+τmax 2 + 1 [START_REF] Thiel | Les distances de chanfrein en analyse d'images : fondements et applications[END_REF]. We obtain the optimal relative error τ opt by τ opt = -τ min (ε opt ) = τ max (ε opt ). It can be shown that the chamfer distance d C induced by any chamfer mask M C is a discrete distance [START_REF] Remy | Normes de chanfrein et axe médian dans le volume discret[END_REF]. However, a distance that is not a norm is not invariant by homothety and this may not be desirable (for instance when comparing skeletons of the same object at different scales). Therefore, we introduce additional criteria to ensure that the computed weights will define a discrete norm.

Norm constraints

A distance is a norm if and only if its ball is convex, symmetric, and homogeneous. For chamfer masks, symmetry is achieved by construction, homogeneity is due to the regular triangulation (also obtained by construction) while convexity can be assessed on the equivalent rational ball of the chamfer mask [START_REF] Remy | Optimizing 3d chamfer masks with norm constraints[END_REF].

Given a chamfer mask

M C = {(v i ∈ Z 3 , ω i ∈ N)} its equivalent rational mask is defined by M ′ C = {(v i /ω i ∈ Q 3 , 1)}.
The polyhedron defined by this equivalent rational mask is the equivalent rational ball (see figure 1).

To check the convexity of the ball, we have to check whether the ball is convex at each of its edges [START_REF] Remy | Optimizing 3d chamfer masks with norm constraints[END_REF]: each edge must be "turned to the outside" of the ball. It turns out that we only have to check a local convexity criterion (LCC) at each edge of the equivalent rational ball. Given 2 faces (P, Q, S) and (Q, R, S) of a triangulation sharing edge (Q, S), the LCC can be expressed as

LCC(P, Q, R, S) = 1 ω P .ω Q .ω R .ω S . x Q x R x S x P y Q y R y S y P z Q z R z S z P ω Q ω R ω S ω P ≥ 0. ( 4
)

Automatic calculation of chamfer mask coefficients

The computation of optimal coefficients for a mask of size (2n + 1) 3 is done in three steps: generation of the Farey triangulation, generation of the norm constraints, and iterative computation of the optimal sets of weights.

Building the Farey triangulation

The recursive automated construction of the Farey triangulation of order n is described in appendix A. This triangulation T g C corresponds to isotropic chamfer mask generator M g C . When dealing with anisotropic lattice, one has to add extra vectors to the mask generator and extra cones to the triangulation. This is achieved by symmetry considerations. For instance, for a 3 3 mask, if the voxel size d z along z is different from the ones along x and y, d x and d y , we have to consider in the mask generator, in addition to the vectors {(1,0,0), (1,1,0), (1,1,1)}, the two extra vectors {(0,0,1), (1,0,1)} that corresponds to weights induced by the anisotropy. These extra vectors belongs to the two extra cones, (1,0,0), (1,0,1), (1,1,1) and (0,0,1), (1,0,1), [START_REF] Pudney | Distance-ordered homotopic thinning: A skeletonization algorithm for 3d digital images[END_REF][START_REF] Pudney | Distance-ordered homotopic thinning: A skeletonization algorithm for 3d digital images[END_REF][START_REF] Pudney | Distance-ordered homotopic thinning: A skeletonization algorithm for 3d digital images[END_REF] , that are to be considered for the error computation and the local convexity constraints.

Generating convexity criteria

The triangulation T g C has been built as described above. It allows us to generate all the local convexity constraints (equation 4) that are to be verified. They have to be generated for every edge inside the mask generator, but also for the edges that are at the border of the mask generator. For the latter, the fourth point (see figure 2) is derived from symmetry considerations.

Please notice that each of the generated LCC depends on 4 weights ω i .

Finding the optimal coefficients

This is the tough part. We have to identify the m-tuples (ω 1 . . . ω m ) of weights corresponding to the chamfer mask generator M g C = {v i , 1 ≤ i ≤ m} to find the optimal ones that yield optimal error. These sets of optimal coefficients are searched by a brute-force method. However, we try to reduce this computationally expensive search by throwing away m-tuples (ω 1 . . . ω m ) as soon as part of them do not satisfy the local convexity constraints (as sketched by below recursive algorithm1 ).

1: procedure Test( n ) 2: if some LCCs can be verified with (ω i , . . . ω n ) then 3:

test these LCCs and return if one of them is not verified 4: if n equals to m then {All ω i are set.}

5:

Compute the error τ opt 6:

if this τ opt is smaller than the previous one then (ω i , . . . , ω m ) is an optimal set of coefficients 8:

return 9: for ω n+1 from ω 1 ||v i || ∞ to ω 1 ||v i || 1 do {Iteratively set a value to ω n+1 .} 10:
Test( n + 1 ) 11: {Main Program} 12: for ω 1 from 1 to some user provided value do 13:

Test( 1 )
ω 1 , the coefficient corresponding to the direction of smallest voxel size, varies from 1 to some maximal value provided by the user, while the other coefficients are searched in the interval [ω

1 ||v i || ∞ , ω 1 ||v i || 1 ].
Error computation is only performed on coefficients sets that verify all the local convexity constraints. As a result, this algorithm gives all the optimal m-tuples in lexicographical order.

Results

Table 1 presents optimal sets of weights of a 3 × 3 × 3 chamfer mask for an anisotropic grid with dx = 1, dy = 1.5, dz = 3.0. The points belonging to this mask are: aX(1, 0, 0), aY (0, 1, 0), aZ(0, 0, 1), bY Z(0, 1, 1), bXZ(1, 0, 1), bXY (1, 1, 0), and c(1, 1, 1). The time needed to compute these sets is 958 ms.

Table 2 presents optimal sets of weights the associated maximum relative error for 7 × 7 × 7 isotropic chamfer masks. The points belonging to this mask are: a(1, 0, 0), b(1, 1, 0), c(1, 1, 1), d(2, 1, 0), e(2, 1, 1), f (2, 2, 1), g(3, 1, 0), h(3, 1, 1), i(3, 2, 0), j(3, 2, 1), k(3, 2, 2), l(3, 3, 1), m(3, 3, 2). The computational times needed to examine all the m-tuples with ω 1 less or equal to 5, 7, 10, and 14 are respectively of 2 min, 25 min, 6 h 37 mn, and 102 h.

Conclusion

We have proposed an automated approach to compute optimal chamfer norm coefficients for mask of any size and for lattice of any anisotropy. It is based on the Farey triangulation that permits us to recursively build large masks while ensuring a regular triangulation of the chamfer mask generators. It allows us to automatically compute the error of any mask, thanks to analytical expressions of errors we can derive on regular cones. In addition, the coefficients we calculate verify norm constraints, thus yields scale invariant chamfer maps.

A Recursive Farey triangulation construction

A.1 Farey Set Points

A Farey set F n of order n is a set of all the irreducible points y x , z x in Q ∩ [0, 1] whose denominator does not exceed n. It is built only with visible points (this means the greatest common divisor of (x, y, z) is 1).

A Farey set of order n correspond to the vectors of the generator of a 3-D chamfer mask of size (2n + 1) 3 . For example, the ordered (lexicographical order) Farey set of order 1

F 1 = {( 0 1 , 0 1 ), ( 1 1 , 0 1 ), ( 1 1 , 1
1 )} correspond to the set of vectors {(1, 0, 0), (1, 1, 0), (1, 1, 1)} which is the generator of an isotropic chamfer mask of size 3 3 . Other vectors that are involved in an anisotropic chamfer mask are deduced from the previous ones by symmetries.

The Farey set of order n + 1, F n+1 , can be built from Fn by

F n+1 = F n ∪ y x , z x + y ′ x ′ , z ′ x ′ with x + x ′ ≤

A.2 Recursive construction of Farey Set triangulations

The triangulation T 1 associated to F 1 is composed of a single cone (1,0,0), (1,1,0), (1,1,1) , or equivalently a Farey triangle ( 0 1 , 0 1 ), ( 11 , 0 1 ), ( 11 , 1 1 ) , that is regular. To build T n+1 from T n , we first put all the Farey triangle in a list L. We now examine successively the triangle in L, and will try to build new triangles by splitting the existing one into two triangles.

Let us consider the triangle A, B, C of L. We try to add a new vertex along its largest edge 2 , say AB. Such a vertex belongs to F n+1 if and only if x a + x b ≤ n + 1. If the latter is not true, the triangle is put again in the list but will no more be considered. If The construction of T n+1 stops when no more triangles, whose vertices are Farey points of order n + 1, can be inserted into L.

Figure 3 shows the different steps of the construction of T 2 from T 1 . T 3 and T 4 are displayed in Figure 4. 2 We consider that large discrepancies between the chamfer distance and the Euclidean one are more likely to occur along the largest edges.

Fig. 1 .

 1 Fig. 1. Equivalent rational ball of a 3D 5 × 5 × 5 isotropic chamfer mask.

Fig. 2 .

 2 Fig. 2. Notations for equation 4.

  a b c d e f g h i j k l m εopt τopt(%)

  x a + x b ≤ n + 1 is true, let us denote C ′ = A +B the new Farey point: the two triangles A, C, C ′ and B, C, C ′ are put into the list L. It can also recursively be shown that those two triangles are regular.

Fig. 3 .

 3 Fig. 3. Construction of T2 from T1.

Fig. 4 .

 4 Fig. 4. T3 and T4.

Table 1 .

 1 3 × 3 × 3 chamfer mask coefficients for anisotropic grid.

	7	13 4.592 14.81

Table 2 .

 2 7 × 7 × 7 chamfer mask coefficients.

	6 6	1.207 20.71

Java code is available from http://www-sop.inria.fr/epidaure/personnel/Celine.Fouard/.


Jessica Goya Outi

Raphael Calmon

Fanny Orlhac

Cathy Philippe

Nathalie Boddaert

Stephanie Puget

Irène Buvat

Vincent Frouin

Jacques Grill

Frédérique Frouin

Jessica Goya-Outi

Stéphanie Puget

Can Structural MRI Radiomics Predict DIPG Histone H3 Mutation and Patient Overall Survival at Diagnosis Time?

Keywords: Structural MRI, Radiomics, Image Standardization, Machine Learning, Rare Cancer

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Radiomics [START_REF] Lambin | Radiomics: Extracting more information from medical images using advanced feature analysis[END_REF] is currently widely investigated in oncology. It aims at extracting multiple quantitative imaging features to identify tumor phenotypes with some predictive values. In this study, we investigate the contributions of radiomics to the diagnosis and prognosis of patients with diffuse intrinsic pontine glioma (DIPG). DIPG is a rare inoperable lethal pediatric cancer frequently associated with histone H3 mutations (H3.1K27M or H3.3K27M). These mutations are currently identified following biopsy and are associated with patient response to therapy [START_REF] Castel | Histone H3F3A and HIST1H3B K27M mutations define two subgroups of diffuse intrinsic pontine gliomas with different prognosis and phenotypes[END_REF].

In this context, we analyzed the ability of radiomic models to distinguish H3 mutation types non-invasively and to predict patient overall survival (OS). The ultimate goal will be to define whether this could avoid biopsy, or replace it when it is not feasible, and guide patient care from diagnosis time. For these prediction tasks, two methods for computing imaging features inside a spherical region of interest included in the tumor were tested, a stringent feature selection procedure was proposed and radiomic signatures were built using different machine learning methods.

CLINICAL AND IMAGE DATA

This retrospective monocentric study included 38 patients with DIPG, scanned at diagnosis with four structural MRI modalities: T1-weighted (T1w), T2-weighted (T2w), T1-weighted post-contrast injection (T1wc) and T2-weighted FLAIR (FLAIR) with the same scanning parameters. Among these patients, 22 were H3.3K27M mutated, 9 were H3.1K27M mutated, 4 were H3.3/1 wild-type and 3 had unknown H3 mutation status. One of the 38 patients had a censored OS (last follow-up 390 days after diagnosis). There was no significant difference in sex between mutation types, H3.3K27M patients were older at diagnosis (Wilcoxon p=0.0009) and survived less time (Wilcoxon p=0.03) than H3.1K27M patients (Table 1). These data were similarly to [START_REF] Castel | Histone H3F3A and HIST1H3B K27M mutations define two subgroups of diffuse intrinsic pontine gliomas with different prognosis and phenotypes[END_REF] and thus showed that the studied sample is representative. Only the 31 patients with H3 mutation were included in the mutation discrimination model while all the 38 patients were included in the OS prediction modeling.

IMAGE FEATURE EXTRACTION

Images were pre-processed by a dedicated pipeline [START_REF] Goya-Outi | Computation of reliable textural indices from multimodal brain MRI: suggestions based on a study of patients with diffuse intrinsic pontine glioma[END_REF] including intensity standardization according to the hybrid white stripe approach, resampling to isotropic voxels (1 mm 3 ) and multi-modal images registration to each T2w scan. A total of 79 indices including first-order and texture features were computed for each MRI modality using PyRadiomics [START_REF] Van Griethuysen | Computational radiomics system to decode the radiographic phenotype[END_REF]. For each patient, a large spherical region was drawn inside the tumor (globT ) on the T2w scans. The same globT was used in T1w, T1wc and FLAIR scans. Two sets of radiomic features were extracted using either global (F globT ) or local (F locT ) approaches (Fig. 1):

• F globT : A total of 316 (79 textural indices x 4 MRI modalities) imaging features were computed within globT ROI;

• F locT : To eliminate the influence of the volume of the ROI on textural values, a small sphere (locT , 5-mm radius) was used to scan every globT , with 1 voxel step, and textural indices were computed within locT ROI, for each locT position. Each globT voxel (v) was characterized by the feature values computed in locT centered on v. The means of each index were used to define the imaging feature set F locT .

For the robustness analysis, further explained in the next section, the same feature extraction procedures were used in two additional scenarii, corresponding to two global regions derived from globT : a) globT dilated, globT d and b) globT eroded, globTe, with a 5-mm radius sphere.

LEARNING MODELS

Three feature sets were used as input in the development of predictive models: Imaging feature set (316 features); Clinical feature set composed of age at diagnosis, sex and the globT volume (approximation of tumor volume); and a Combined set of clinical and imaging features (316 + 3 features) was also investigated.

Since the number of patients is too small for effective training/test sample splitting, Leave-One-Out Cross-Validation (LOOCV) was applied in feature selection and machine learning steps for model performances estimation in both H3 mutation prediction and OS prediction. It is important to note that the feature selection is performed inside each LOOCV fold, preventing selection bias [START_REF] Ambroise | Selection bias in gene extraction on the basis of microarray gene-expression data[END_REF]. For each LOOCV fold, all 316 imaging and 3 clinical features were standardized by mean subtraction and unit variance scaling in the training set and the same normalization parameters were then used to normalize the validation set. Models using F globT and F locT were estimated separately.

H3 mutation prediction

In order to prevent over-fitting and make results interpretable, a small number of features should be selected, given the limited number of patients. Aiming to select robust, informative and non redundant features, a three steps selection procedure was applied to the imaging features.

Step 1: Features were selected according to their robustness to the spherical ROI delineations. Using the three definitions of tumor region (globT , globTe and globT d ), the absolute agreement intraclass correlation coefficient (ICC) of each feature was computed. Only features with ICC > 0.9 were kept. Step 2: Features presenting an individual Area Under the Receiver Operating Characteristic curve (AUC) < 0.75 were excluded. This threshold was a compromise between keeping features that could combine well with others and excluding those that could degrade the model. Step 3: To reduce redundancy, hierarchical clustering was performed, keeping the minimum absolute Spearmans correlation coefficient (|r|) between cluster members greater than 0.85. The feature with the greatest AUC of each cluster was finally selected. For the clinical feature set, only steps 2 and 3 were applied. If none of the features was selected, an univariate model was built with the feature presenting the greatest AUC.

The minority class H3.1 was resampled using regular Synthetic minority Over-sampling Technique (SMOTE) [START_REF] Chawla | Smote: Synthetic minority over-sampling technique[END_REF] in the training set. Support Vector Machine (SVM), K-Nearest Neighbors (KNN) and Random Forest (RF) were applied. Briefly, linear kernel was used in SVM and its penalty parameter C was set to 1. The parameter K was set to 3 in KNN. The number of trees was set to 100 in RF.

OS prediction

For the OS regression, the same feature selection steps were applied. However, in Step 2 and in Step 3, AUC < 0.75 criterion was replaced by C-index < 0.55. The C-index is the equivalent of the AUC for regression problems.

Cox proportional hazards (COXPH) combined with inner-LOOCV Least Absolute Shrinkage and Selection Operator (LASSO) and with inner-LOOCV Ridge regression as well as Random Survival Forest (RSF) were used to estimate patients risks. For each LOOCV fold, the median estimated risk value in the training set was used as a threshold to classify the validation patient into low-risk or high-risk groups. The number of trees was set to 1000 and minimal node size to 5 in RSF.

STATISTICAL ANALYSIS

H3 mutation prediction

Performances were estimated by the LOOCV F1-weighted score, grouping the validation patient results of every fold. The LOOCV was repeated 30 times to estimate the standard deviation of the prediction scores due to randomness introduced by SMOTE and RF.

OS prediction

Each patient was classified into low and high risk group using a model computed on the training set that they were not part of. All the patients classified as low-risk in any of the folds of the LOOCV were grouped together and a single Kaplan-Meier curve was computed for that low-risk group. The same was done for the high-risk group.

A permutation based strategy was used to assess the statistical significance since it is not possible to use the cross-validated survival curves directly because the observations are not independent. The statistical significance of the LOOCV log-rank statistic was then obtained from the permutation distribution of the LOOCV log-rank statistic. It was estimated by 500 repetitions of the whole learning model process with random permutations of patients' OS time [START_REF] Rm Simon | Using cross-validation to evaluate predictive accuracy of survival risk classifiers based on high-dimensional data[END_REF].

RESULTS

H3 mutation prediction

In average 10±1 local imaging features and 7±1 global imaging features were selected in each fold after the 3 steps selection. Features selected in each fold tended to be the same, 9 local features and 6 global features were selected in at least 20 out of the 31 folds. Concerning the clinical features, age at diagnosis was selected in every fold, volume in 7 out of 31 folds and sex in none of the folds. The F1-weighted scores for all three machine learning methods using local and global imaging features are shown in Fig. 2(a). Each point correspond to one of the 30 trials. For global extraction, F globT , F1-weighted scores for SVM were slightly better than for the KNN and RF methods. For local extraction, F locT , F1-weighted scores for KNN were the best followed by RF and SVM. Furthermore, performances of models with local and global features are comparable. Overall, all machine learning models performed well for H3 mutation classification based on imaging features (F1-weighted score around 0.74). For the sake of comparison, the F1-weighted score for a model classifying all the samples as the most frequent one (H3.3) is 0.59.

In Fig. 2(b), the performances of SVM models using imaging features are compared to SVM models using clinical and combined (clinical + imaging) features. For global extraction, the imaging model is better than the clinical one and for both local and global extraction, the combined models are better than the imaging model alone. The SVM model using combined features was better with F globT than with F locT . The same trend (mean (standard deviation)) was observed for the F1-weighted scores values of KNN and RF models using combined features: KNN local, 0.76 (0.04), KNN global 0.81 (0.02); RF local 0.74 (0.04), RF global 0.79 (0.04).

Based on the performances estimation with the LOOCV, the final proposed model using all 31 patients was computed by SVM using global extracted imaging features combined with the clinical features. The resulting selected features were: 1) Age at diagnosis, 2) T1wc GLSZM GrayLevelVariance, 3) T1w first-order Median, 4) T1wc first-order 10Percentile, 5) FLAIR GLCM Homogeneity1, 6) FLAIR first-order Median and 7) FLAIR first-order 90Percentile.

OS prediction

After the 3 steps feature selection 15±2 local imaging features and 7±1 global imaging features were selected in each fold. Features selected in each fold tended to be the same, 10 local features and 6 global features were selected in more than 26 out of the 38 folds. Concerning the clinical features, sex and volume were selected in 22 out of 38 folds while age at diagnosis only in 11 folds. In 5 folds, none of the clinical features were selected leading to a combined model identical to the imaging model. The cross-validated Kaplan-Meier curves (using the validation patient of each trial) for all the three regression methods, clinical, imaging (global and local), and clinical + imaging (global and local) feature dataset are shown in Fig. 3. COXPH Ridge regression and RSF invert the low and high risk groups in validation for all four sets of features including imaging features, while COXPH LASSO presents this behavior for the two sets of features including F globT . COXPH LASSO regression presented adequate survival curves for the two sets of features using F locT . For the sake of comparison, training Kaplan-Meier curves for the model learned were computed using all 38 patients. All the RSF and most of the ridge regressions showed significant log-rank pvalues, suggesting over-fitting. Features selected in the final model using F locT imaging set presenting non-zero weight after LASSO (a) Imaging feature set models (results using F globT are in red, F locT in green) (b) SVM models using clinical (in blue), imaging and combined feature set (results using F globT are in red, F locT in green) reduction (the method with the most satisfactory LOOCV result) were: 1) FLAIR GLSZM LowIntensitySmallAreaEmphasis and 2) T2w first-order Maximum.

DISCUSSION

Recently, radiomics has been promising to reveal patients prognosis or phenotypes profiles non-invasively. In this study, we investigated if radiomic features from diagnostic structural MRI could reveal the H3 mutation associated with each DIPG tumor and predict patient OS. Results showed that good H3 mutation classification scores were found by three different machine learning methods when using imaging features previously selected by our proposed selection pipeline. Moreover, they performed better than a clinical model and F1-weighted scores increased when combining imaging and clinical features. Three first-order and textural imaging features issued from FLAIR, two from T1wc and one from T1w were selected in the final proposed model along with age at diagnosis. Future tests in an independent data set should be performed to confirm these scores. Local feature extraction was motivated by the presumption that features reflecting texture information could be formerly hidden by volume bias [START_REF] Goya-Outi | Computation of reliable textural indices from multimodal brain MRI: suggestions based on a study of patients with diffuse intrinsic pontine glioma[END_REF]. However, even though more local features were selected after our feature selection pipeline, models using global features presented equivalent or better F1-weighted scores for H3 mutation classification. Therefore, for this classification task, the interest of local extraction is not demonstrated. As the computation burden is high for the local extraction, we could suggest removing these local features for our particular application. But this recommendation cannot be generalized, and the potential impact of local features deserves to be tested for some other predictive tasks.

Considering OS prediction, none of the 15 different modeling configurations could separate significantly low and high risk patient groups at diagnosis time. Nevertheless, the problem that we attempted to address is complex since different mutation populations were mixed up and only non-invasive diagnostic features (i.e. ignoring data derived from biopsy) were used. Actually, DIPG longterm survivals are defined as patients having survived longer than two years after diagnosis [START_REF] Lm Hoffman | Clinical, radiologic, pathologic, and molecular characteristics of long-term survivors of diffuse intrinsic pontine glioma (DIPG): A collaborative report from the international and european society for pediatric oncology dipg registries[END_REF]. Since too few patients (five in our dataset) had long-term survival times, we could not use this classification. As several methods, especially RSF, presented good training performances but poor validation performances, we investigated some possible over-fitting. We tested a modified version of RSF [START_REF] Chen | Pathway hunting by random survival forests[END_REF] with minimal depth pathway hunting, which was more drastic than our three-step feature selection. As shown in Fig. 4, this RSF procedure provided better validation results. This study was undertaken with a low number of patients, which makes the use of radiomics and machine learning methods very challenging. This low number could not be increased at the present time for the three following reasons. First, DIPG is a rare disease. Second, the use of retrospective homogeneous datasets, including the four MRI modalities, has reduced the number of admissible patients. Third, as newly diagnosed patients are involved in on-going therapeutic trials, their data are not currently available. However, this first discovery study was necessary to define the potential interest of the different structural images. For instance, the interest of FLAIR and post-contrast T1 weighted images to understand the different types of H3 mutation was first illustrated by our results. This interest remains to be confirmed using data of an on-going clinical trial, when they will be available.

CONCLUSION

The combination of radiomic features, including first-order and textural indices derived from structural MRI at diagnosis, and of clinical data were found to be predictive of two types of histone H3 mutation in patients with DIPG. However, when using these features, it was not possible to validate a model discriminating two subgroups of patients with statistically significant differences in their overall survival.

Fig. 1 .

 1 Fig. 1. Exemplification of local (F locT ) and global (F globT ) feature extraction approaches.

Fig. 2 .

 2 Fig. 2. Box-and-whisker plots of the F1-weighted validation scores. Values for each of the 30 LOOCV trials are displayed. Mean (red square) and standard deviation values are displayed. The black line represents the F1-weighted score, 0.59, corresponding to all patients classified as H3.3.

Fig. 3 .

 3 Fig. 3. Cross-validated Kaplan-Meier curves for the three overall survival models and five different combinations of features. Highrisk in red, low-risk in blue.

Fig. 4 .

 4 Fig. 4. Cross-validated Kaplan-Meier curves for the pathway hunting RSF procedure. High-risk in red, low-risk in blue.

Table 1 .

 1 Characteristics of 38 DIPG patients included in this study. Age at diagnosis is given in years, overall survival (OS) in days.

		H3.1	H3.3	WT/unknown
	Patients	9	22	4/3
	Age	5.0±.3	8.6±3.3	6.7±3.5
	Boys/girls 4/5	8/14	5/2
	OS	531±281 328±170 367±221


Nathali A Franchina Vergel

Post ‡, ⊥ L Christiaan

Davide Sciacca

Maxime Berthe

François Vaurette

Yannick Lambert

Dmitri Yarekha

David Troadec

Christophe Coinon

Guillaume Fleury

Gilles Patriarche

Tao Xu

Ludovic Desplanque

Xavier Wallart

Daniel Vanmaekelbergh

Christophe Delerue

Bruno Grandidier

Engineering a robust flat band in III-V semiconductor heterostructures

Keywords: Two-dimensional lattice, III-V semiconductor, quantum well, band engineering, flat band, disorder, block copolymer lithography, scanning tunneling spectroscopy, tight binding calculations

Electron states in semiconductor materials can be modified by quantum confinement. Adding to semiconductor heterostructures the concept of lateral geometry offers the possibility to further tailor the electronic band structure, with the creation of unique flat bands.

Using block copolymer lithography, we describe the design, fabrication and characterisation of multi-orbital bands in a honeycomb In0.53Ga0.47As/InP heterostructure quantum well with a lattice constant of 21 nm. Thanks to an optimized surface quality, scanning tunnelling spectroscopy reveals the existence of a strong resonance localized between the lattice sites, signature of a porbital flat band. Together with theoretical computations, the impact of the nanopatterning imperfections on the band structure is examined. We show that the flat band is protected against the lateral and vertical disorder, making this industry-standard system particularly attractive for the study of exotic phases of matter.

The propagation of waves interacting with a periodic potential landscape in a two-dimensional (2D) system leads to the formation of energy bands for wavelengths on the order of the lattice constant. For some particular periodicities, non-dispersive and macroscopically degenerated flat bands emerge (FBs). Well-known systems with a flat band include the Dice, Lieb's, Kagome lattices [START_REF] Leykam | Artificial flat band systems: from lattice models to experiments[END_REF] and lattices formed from a Moiré pattern such as twisted double-layer graphene. [START_REF] Bistritzer | Moiré bands in twisted double-layer graphene[END_REF] A special case is the honeycomb lattice: flat bands exist due to the incompatibility of the in-plane px,py orbitals and the trigonal lattice, [START_REF] Wu | Flat bands and Wigner crystallization in the honeycomb optical lattice[END_REF] but they only appear in the absence of s-p hybridisation. All these lattices are theoretically described with tight-binding Hamiltonians and flat band systems have experimentally been realized by designing lattices with micropillars for polaritons, [START_REF] Jacqmin | Direct observation of Dirac cones and a flatband in a honeycomb lattice for polaritons[END_REF] optical waveguides for photons, [START_REF] Mukherjee | Observation of a localized flat-band state in a photonic Lieb lattice[END_REF] optical lattices for ultracold atoms, [START_REF] Taie | Coherent driving and freezing of bosonic matter wave in an optical Lieb lattice[END_REF] and array of vacancies or molecules for electrons [START_REF] Drost | Topological states in engineered atomic lattices[END_REF][START_REF] Slot | Experimental realization and characterization of an electronic Lieb lattice[END_REF][START_REF] Gardenier | Orbital Flat Band and Dirac Cone in the Electronic Honeycomb Lattice Article[END_REF] and magic-angle twisted bilayer graphene for electrons. [START_REF] Cao | Correlated insulator behaviour at half-filling in magic-angle graphene superlattices[END_REF] Transposing this approach to semiconductor heterostructures, in which the Fermi level, the carrier density and the lattice constant can be tuned, is particularly attractive. This would open the gate to novel electronic phases and band structure topology, even more so when geometric lattice concepts could be combined with spin-orbit coupling, and charge-and spin-interactions.

For standard vertical semiconductor heterostructures grown by the semiconductor industry, a suitable potential landscape would consist of repulsive potential barriers arranged in a lattice inside a quantum well (QW), so that the 2D electrons are periodically scattered by the array of potential barriers. [START_REF] Gibertini | Engineering artificial graphene in a two-dimensional electron gas[END_REF][START_REF] Wang | Observation of Dirac bands in artificial graphene in small-period nanopatterned GaAs quantum wells[END_REF] Empty pores perforated in the QW, acting as potential barriers, would ensure the highest energy potential and dielectric contrasts, satisfying the need for a strong lateral quantum confinement. A further increase of the confinement could be obtained by the selection of semiconductor materials with a small effective mass, such as indium-containing III-V compounds.

A tight binding prediction for an ultrathin InAs film nanopatterned with different geometric lattices has shown that resolving distinct energy-momentum dispersions is feasible for lateral geometries with lattice constants smaller than 20 nm. [START_REF] Tadjine | From lattice Hamiltonians to tunable band structures by lithographic design[END_REF] So far, these dimensions were out of reach for topdown technologies used to process III-V QWs, [START_REF] Wang | Observation of Dirac bands in artificial graphene in small-period nanopatterned GaAs quantum wells[END_REF] the smallest achievable periodicity with e-beam lithography for the nanoperforation of a triangular antidot lattice being 40 nm. 14 As a model III-V semiconductor heterostructure, we consider here a 10 nm-thick In0.53Ga0.47As QW grown by molecular beam epitxy (MBE) on an InP (001) substrate. As shown in Figure 1, band theory predicts that the step-like 2D density of states (DOS) of the QW in the conduction band gets featured upon its perforation with tiny pores arranged in a triangular lattice with a periodicity D of 36 nm. This arises from the interplay between band folding in the Brillouin zone of the superlattice and the gap opening at the edges or at the centre of the Brillouin zone due to scattering of the electronic waves by the periodic potential of the pores (Figure 1(b)). As the pore size increases, these features get more pronounced and separate in energy due to a stronger lateral confinement (Figure S3). At a pore diameter d of 20 nm, the DOS consists of well separated mini-bands (Figure 1(c)). This triangular lattice of pores is equivalent to a honeycomb lattice with a lattice constant a=D/3 ~ 21 nm. In this honeycomb lattice, we can conceptually define lattice sites at the barycenter of three adjacent pores. Each such lattice site can be considered as a quantum dot, characterized by s and px,y-orbitals in the lowest part of the spectrum, the pz orbitals being only involved in the high-energy states coming from the second subband of the QW. If s-p hybridization can be avoided, the electronic coupling between the lattice sites leads to two dispersive s-bands forming a Dirac cone, and a set of four px,y-orbital bands. Among the latter bands, the two dispersive px,y-bands form Dirac cones (DC) at the K points, with a flat-band beneath the DC, and one above. The hexagonal symmetry is not compatible with the orthogonal in-plane px,y-orbitals, resulting in destructive interferences and the creation of FBs. [START_REF] Wu | y-orbital counterpart of graphene: Cold atoms in the honeycomb optical lattice[END_REF] Due to the porbital character, the FBs should give rise to a strong peak in the local density of states (LDOS) at positions between the lattice sites (which we will call bridge sites) and zero-density of states on the lattice sites themselves (Figure 1(a) and Figure S5). Crucial to the existence of the lower porbital FB is the energy separation with the s-band, which takes place for ratios d/D higher than 0.5 (Figures S4 andS6). At this stage, a major question is to what extent the disorder in the lattice will cloud the existence of such a FB. Despite numerous theoretical predictions of the influence of disorder on specific energy-momentum dispersions obtained with translationally invariant Hamiltonians, [START_REF] Bodyfelt | Flatbands under correlated perturbations[END_REF][START_REF] Soriano | Three-dimensional models of topological insulators: Engineering of dirac cones and robustness of the spin texture[END_REF][START_REF] Bilitewski | Disordered flat bands on the kagome lattice[END_REF] the actual fate of the FB is still unknown, as a gradual increase of disorder can lead first to the destruction of specific momentum-dispersions and then to their reappearance. [START_REF] Schubert | Fate of topological-insulator surface states under strong disorder[END_REF] Aiming at engineering well-separated mini-bands with pure orbital character and band gaps larger than the thermal fluctuations at room temperature (25 meV), our approach takes advantage of the wafer-scale and low-cost self-assembled block copolymer (BCP) lithography. [START_REF] Suh | Sub-10-nm patterning via directed self-assembly of block copolymer films with a vapour-phase deposited topcoat[END_REF] This technique enables the creation of a honeycomb lattice in an In0.53Ga0.47As/InP heterostructure QW with a lattice constant of 21 nm. Moreover, the selective removal of the (poly)methylmethacrylate (PMMA) domains in the BCP leads to a poly(styrene) (PS) honeycomb mask with an intended mean pore size of 20 nm, yielding a ratio d/D of 0.56 for a periodicity of 36 nm. While such a ratio is suitable for the creation of a FB, it also renders the s-DC bandwidth to 7 meV at the M points of the Brillouin zone (Figure S4). Direct access to the local density of states (LDOS) in the lattice is obtained with scanning tunnelling spectroscopy (STS), which has beautifully revealed the existence of FBs in 2D atomic crystals [START_REF] Xie | Spectroscopic signatures of many-body correlations in magic-angle twisted bilayer graphene[END_REF][START_REF] Zhang | Flat bands in twisted bilayer transition metal dichalcogenides[END_REF][START_REF] Mao | Evidence of Flat Bands and Correlated States in Buckled Graphene Superlattices[END_REF] or atomic and molecular lattices on metal surfaces [START_REF] Drost | Topological states in engineered atomic lattices[END_REF][START_REF] Slot | Experimental realization and characterization of an electronic Lieb lattice[END_REF][START_REF] Gardenier | Orbital Flat Band and Dirac Cone in the Electronic Honeycomb Lattice Article[END_REF] .

The STS results that we present on the In0.53Ga0.47As honeycomb lattice show not only the formation of well-separated dispersive s and p bands, but also the existence of a p-orbital FB.

Comparison with tight binding calculations, which are particularly suitable to take into account the effect of a random disorder on the lattice sites, reveals that the lattice imperfections perturb the bands, but does not erase them. Especially, the flat band is omnipresent on all bridge sites of the honeycomb lattice, even for imperfections producing local potential fluctuations of ± 5 meV, an energy perturbation comparable to the designed width of the dispersive bands and gaps.

The top-down patterning methodology to produce the honeycomb lattice is illustrated in Figure 1(d) and, with more details, in Figure S10. The plasma-enhanced chemical vapour deposition of a 16 nm-thick SiO2 layer on the QW, prior to the spin coating of the BCP thin film, warrants a successful pattern transfer to the In0.53Ga0.47As QW. As seen in Figure 2(a), after the removal of this protective layer, the QW shows periodic pores, arranged into the micro-domain structure of the BCP mask (Figure S11). The generated antidot lattice has a periodicity of D = 36.0 ± 2.8 nm and a pore size of d = 20.0 ± 1.4 nm (Figure 2(b,c)). As the deposition of the SiO2 layer leads to the oxidation of the QW surface, which is further altered by the etching of the SiO2 mask at the end of the whole workflow, a cleaning of the surface with atomic hydrogen and As4 flux was performed in the MBE chamber. 24 It was followed by a thin amorphous As capping, 25 which is essential to recover a clean and well-ordered In0.53Ga0.47As (001) surface, 26 as demonstrated in the STM image of Figure 2d. The surface structure is similar to the one observed on the surface of the pristine QW (Figure S12), consistent with an As-rich (2x4) surface. Although the STM tip is not sharp enough to reach the bottom of the pore, cross-sectional scanning transmission electron microscopy (STEM) of the heterostructure indicates the full perforation of the QW (Figure 2(f)). In order to resolve the DOS in the conduction band, the QW was doped with Be at a concentration higher than 10 19 cm -3 . The large hole density in the QW ensures the fast recombination of the tunnelling electrons, avoiding in this way any charging of the QW. Above all, the dense hole gas minimizes the potential fluctuations caused by the random distribution of charged impurities and defects. 26 Efficient screening is preserved, even after nanoperforation of the QW, as the pinning of the Fermi level at the surface remains close to the top of the valence band with deviations smaller than 85 meV for remote regions of the sample (Figure S12). A mean position of the Fermi level of 160 ± 11 meV above the valence band edge indicates a low density of defects, which is remarkable in view of the technological processes in the nanopatterning.

Spatially-resolved tunnelling spectra were first acquired on positions between two lattice sites (Figure 3(a)). Right at the bridge site (B), the differential conductance spectrum reveals a strong peak centred at 788 ± 3 meV, with two tiny shoulders pointed by arrows at smaller bias (770 meV and 754 meV) and another one at higher bias (798 meV) (Figure 3(b)). The peak gradually vanishes away from the bridge site, and is absent at the lattice sites L1, L2. On the honeycomb lattice site L2, a doublet appears at lower energy, and this is also the case for the position next to the lattice sites L2 and L1. At position L1 and further away from L1, we do not detect any signal below 800 meV.

As shown in Figure 3(c), we observe a height difference of two bilayers between lattice sites L1 and L2, very probably the reason that the peak doublet cannot be observed at position L1 and away from L1 with the chosen STS feedback parameters (Figure S13).

A careful analysis of the surface corrugation (Figures 2(a), 3(a) and Figure S14) reveals that the etching of the QW is usually more pronounced at the bridge sites, making the centre of mass of the square wave function of the confined states closer to the surface and allowing a stronger coupling of these states with the tip states. Focussing the STS measurements in the vicinity of bridge sites for many other regions of the sample, some of them being separated by hundreds of micrometers from each other, we consistently measured a strong peak at the bridge sites gradually vanishing and disappearing 2-3 nm away from the bridge site (Figure 3(d) and Figure S13). In line with the theoretical predictions of Figure 1(a) and Figure S5, we assign this peak to the lowest p-orbital FB (notice the energy shift as the experimental reference is given by the position of the Fermi level). Then, the lower-energy peak doublet at (and close to) the lattice sites should correspond to the s-orbital Dirac cone; the peak maxima reflecting the maximum LDOS at the M points. We observe a strong variation in the shape and peak separation of this doublet. For instance, the energy difference between the two maxima is often substantially larger than the 7 meV expected for an ideal honeycomb lattice.

The structural analysis of the nanoperforated QW reveals obvious variation in the pore geometry and in the thickness t of individual lattice sites. The conical shape of the pores does not significantly affect the band structure. [START_REF] Post | Triangular nanoperforation and band engineering of InGaAs quantum wells: a lithographic route toward Dirac cones in III-V semiconductors[END_REF] But variation in the thickness and deviations from the ideal trigonal geometry of pores alter the vertical and lateral confinement of each lattice site with direct impact on the s and p-orbital on-site energies and the coupling between two adjacent sites.

Calculation of the increase of confinement due to a bilayer reduction, as observed in Figure 2(e), yields an energy of 5.6 meV for the s-state. Similarly, calculations of the s on-site energies and of the hopping integral between the s and p orbitals show an increase of the on-site energy by 3.4 meV and a reduction of the hopping integral by 1.2 meV as the ratio d/D increases from 0.55 to 0.60 (Figure S7). These energy variations have a magnitude in the same order as the separation between the s-orbital doublet and the FB peak, and the peak separations observed in the s-orbital doublet. We are hence in a situation in which disorder affects the energetics but does not cloud the intrinsic band features of the lattice. Gaining a more accurate understanding of the influence of disorder was achieved by considering a simple (effective) tight-binding Hamiltonian, where each site of the honeycomb lattice is described by one s and two in-plane px,y-orbitals (see Methods in the Supporting Information). For that purpose, the electronic structure of the ideal honeycomb lattice was first calculated with an effective mass approximation (EMA) and compared with a fully atomistic tightbinding calculation. Good agreement between both calculations was obtained when non-parabolic corrections were included in EMA. The parameters of the simple tight-binding Hamiltonian, defined by two on-site energies ES and EP and four nearest-neighbour hopping parameters were then derived from the fit of the band structure calculated with the EMA model (Figure 4(a)). As the on-site energies are more affected by the structural disorder than the hopping terms (Figure S7), a random energy perturbation was applied on ES and EP of each site to explicitly inform on the effect of disorder (Figure 4(b)). The amplitude of the on-site energy variation was raised by steps of 1 meV. For an amplitude of ± 5 meV, which is in line with the observed structural disorder, the LDOS calculated between two particular lattice sites shows that the disorder has only minute effects on the p-orbital FB (Figure 4(c)). The simulated disorder has more effect on the peak doublet of the s-orbital Dirac cone. It pulls the peaks in the first doublet (M points of the s band) further apart, adding 4 meV to their natural separation. Calculations for tens of different sites (Figure 4(d) and Figure S8) shows a substantial variation in the doublet peak shapes and peak separation for the dispersive bands, whereas the p-FB is preserved. Moreover, the simulated disorder shifts states of the dispersive p-band towards the p-FB. These effects are in line with the experimental observations. We simulate a spatial variation in the energy of the p-FB of 8 meV, while the experimental variation is around 19 meV in Figure 3(c). We attribute this difference between simulations and experiment to small variations of the Fermi level pinning. Indeed, a variation in the position of the Fermi level induces a long-range potential shift, which is not included in the calculations. The resilience of the FB peak to disorder is attributed to the combination of several effects detailed in the SI. The main one is that FB states can be described as localized states on hexagonal plaquettes. [START_REF] Wu | y-orbital counterpart of graphene: Cold atoms in the honeycomb optical lattice[END_REF] It can be shown that the coupling of these states with the other states induced by the disorder is weak so that the broadening of the FB peak comes essentially from the average disorder potential in these states which, given the large number of states available in this band, results in a LDOS peak that remains very strongly marked.

To conclude, using the honeycomb lattice as a prominent example, we presented the basis of a future technological platform for III-V semiconductor heterostructures in which the band structure can be engineered by nanoperforation with lattice constants ~ 20 nm. The intrinsic honeycomb resonances are perturbed by nanoscale disorder, but still remain visible over the entire honeycomb structure. In contrast to the s-orbital Dirac cone, which was not unambiguously resolved in the tunnelling spectra, the p-orbital flat band, arising from the interference of in-plane p-orbitals, is an intrinsic lattice resonance that is strongly resilient against the disorder induced by imperfect nanoperforation. Combining well-known semiconductor heterostructures that allow for electronic gating with band-engineering by nanopatterning is one of the most rational paths toward exotic states and topological protection.

ASSOCIATED CONTENT Supporting Information. The Supporting Information are available free of charge. Whole electronic structure calculations; details of the sample fabrication; additional scanning tunneling spectroscopic measurements.

Figure 1 .

 1 Figure 1. Band structure engineering upon the periodic nanoperforation of a quantum well. (a) Variation of the conduction-band density of states of an In0.53Ga0.47As quantum well grown on a InP (001) substrate, which is nanoperforated with increasing pore sizes. The thickness t of the quantum well, the periodicity D of the lattice and the diameter d of the pores are indicated. The local density of states at the lattice (L) and bridge (B) sites are superimposed to the (total) density of states for a pore diameter of 20 nm. Vertical arrows point to key features of the density of states, such as Dirac cones (DC) and a flat band (FB). The energy reference is the top of the bulk valence band. (b) Band structure of the quantum well in the first Brillouin zone and in the reduced zone scheme obtained for a periodicity of 36 nm in the real space. (c) Band structure of the nanoperforated quantum well. (d) Illustration of the fabrication process using a self-assembled block copolymer (BCP) thin film as a mask to obtain a honeycomb lattice with a lattice constant a = D/3.

Figure 2 .

 2 Figure 2. Structural properties of the nanoperforated In0.53Ga0.47As quantum well. (a) SEM image showing the successful pattern transfer. (b) and (c), Distributions of the pore periodicity and pore diameter measured from SEM images. (d) Atomically-resolved STM image of the nanoperforated In0.53Ga0.47As (001) surface, showing variability in the pore size and height of the lattice sites. (e) Height profile measured along the dash segment of the STM image. The corrugation seen on the atomic terraces arises from the As-rich (2x4) reconstruction, which is highlighted in the graph. (f) High resolution cross-sectional high angle annular dark field STEM image of the nanoperforated heterostucture. Inset: Larger cross-sectional view.

Figure 3 .

 3 Figure 3. Spatially-resolved scanning tunnelling spectroscopy of the honeycomb lattice. (a) topographic and (b) spectroscopic measurements at varying positions across two adjacent lattices via a bridge site, labelled L1, L2 and B respectively. The vertical segment indicates the position of the p-orbital flat band, whereas the arrows point to spectral features of the s and p-dispersive bands. The dotted line spectrum was measured on the pristine QW for comparison. Inset: Magnified view on the spectral region showing the s and p-flat bands at the bridge site. (c) Contour followed by the STM tip scanning at constant height between the acquisition of the spectra. The electronic coupling of the tip states with the wave function Ψs of the antidot lattice depends on the thickness of the lattice sites L1 or L2. (d) Check of the occurrence of the p-orbital FB peak on the bridge site concomitant with its absence on the adjacent lattice site in four regions of the sample. Insets: STM images of the corresponding regions, which show various degrees of patterning disorder.

Figure 4 .

 4 Figure 4. (a) Disorder-free band structure calculated with the effective tight-binding model, where each site of the lattice is described by one s and two px,y-orbitals. (b) Portion of the hexagonal

ACKNOWLEDGMENT We thank D. Stiévenard for fruitful discussions. We gratefully acknowledge Arkema and Digital Surf for providing the materials to prepare the block copolymer mask and routines to process the STS data. Funding: This work was supported by the French National Research Agency (Dirac III-V project ANR-16-CE24-0007-01, Germanene project ANR-17-CE09-0021-03), the EQUIPEX program Excelsior (Grant No. ANR-11-EQPX-0015), the RENATECH network, the H2020 program (ERC Advanced Grant 692691-"FIRST STEP"), the Dutch Research Council (NWO Chemistry -Toppunt "Superficial superstructures") and the Natural Science Foundation of Shanghai (19ZR1419500).

AUTHOR INFORMATION

Corresponding Author *christophe.delerue@iemn.fr; bruno.grandidier@iemn.fr; D.Vanmaekelbergh@uu.nl; xtld@shu.edu.cn.

Author Contributions