Chourouk Ghodhbane
email: chourouk.ghodhbane@b-com.com

Patrick Savelli
email: patrick.savelli@b-com.com

Cédric Gueguen
email: cedric.gueguen@irisa.fr

Load-Efficiency-Balance Cell Selection Policy for IAB Networks

Keywords: 5G, IAB, wireless connection, cell selection, capacity

The 3rd Generation Partnership Project (3GPP) has proposed recently the Integrated Access and Backhaul (IAB) to simplify the deployment of new base stations, where the connection between the User Equipment (UE) and the Core Network (CN) is provided by a multi-hop 5G wireless connection between the IAB nodes. The main contribution of this paper is to propose a cell selection policy specific to IAB, that aims to guarantee the lowest possible total transmission bandwidth cost (including the backhaul cost which is a significant factor in an IAB context) while respecting the base stations capacities, in order to boost IAB network capacity. In this context, a comparative performance study is performed between our proposed policy and state-ofthe-art policies. System level simulation results show that our solution provides the best system capacity via a trade-off between spectral efficiency and load balancing.

I. INTRODUCTION

Fifth-generation wireless networks present a great challenge in terms of massive number of users with very high data demands, which requires dense deployment of base stations, especially when deployed using Frequency Range 2 (FR2). However, base station densification requires expensive wired infrastructure. To cope with this issue, the Integrated Access and Backhaul (IAB) architecture and its radio protocols have recently been standardized by the 3rd Generation Partnership Project (3GPP) in Release 16 [START_REF]NR and NG-RAN Overall Description[END_REF]. IAB provides a valuable alternative for the densification of cellular networks that avoids the higher cost of a wired backhaul network deployment [START_REF]Study on integrated access and backhaul[END_REF]. In addition to the traditional 5G-based macro cell generally connected with fiber to the Core Network (CN), the IAB architecture allows the operator to easily deploy micro Base Stations (BSs) wirelessly connected to the macro cell. In other words, two types of nodes are defined in an IAB network: the IAB donor, which is a macro cell connected to the CN using wired backhaul, and the IAB node, which is a micro cell using wireless backhaul to its parent node or to the donor. Wireless backhaul links share the same frequency band as access links, which is called inband backhaul. The reuse of the same frequency in IAB allows the operator to avoid paying for a new frequency. The IAB topology offers the possibility for the User Equipement (UE) to connect to the core network either directly through the IAB donor, or indirectly by multi-hop relaying connections through IAB nodes. The direct connection is represented via a direct link between the UE and the IAB donor. The indirect connection is represented via an access link between the UE and a IAB-node and one or more backhaul links between a IAB-node and the IAB donor depending on the number of hops. Even though we can assume better propagation conditions on the backhauling links with respect to the access link, the extra resource cost of transmitting backhaul data in IAB is significant. Consequently, if we go through an IAB-node, the throughput depends on two links, moreover, the IAB donor-IAB node link is shared between several users and the load on this link has an influence. Thus, several strategies for cell selection can be employed, but until now there is no standardized policy. In this paper, we propose a new cell selection policy specific to IAB networks based on a information broadcast between BSs and UEs, that aims to guarantee the lowest possible total transmission bandwidth cost while respecting the base stations capacities. The rest of the paper is organized as follows. In section II, the state of the art is presented. Section III describes our proposal. Sections IV and V include the description of our system model and performance comparison results of our policy with existing policies, respectively. Finally, Section VI concludes the paper and discusses possible extensions of this work.

II. STATE OF ART

The Reference Signal Received Power (RSRP) based cell selection policy [START_REF]User Equipment (UE) procedures in Idle mode and RRC Inactive state[END_REF] is classicaly used in 5G networks. With this policy, the UE compares the strength of signals from neighboring BSs and selects the cell corresponding to the maximal RSRP. The implementation of this policy in an IAB network can result in unbalanced network topology and lower performance [START_REF] Polese | Distributed path selection strategies for integrated access and backhaul at mmWaves[END_REF]. As shown in the example of Fig. 1, if user 4 selects his parent based on RSRP, he will associate with micro 2, which offers the best RSRP value -110 dBm but is also the furthest from the CN in terms of number of hops, which requires more radio resources and increase the latency. The implementation of RSRP-based policy is straightforward since users can measure the quality of the link using the synchronization signals transmitted by the base stations, but has limited performance in an IAB network. For this reason, additional criteria must also be taken into account during the selection of cells to improve performance. In [START_REF] Polese | Distributed path selection strategies for integrated access and backhaul at mmWaves[END_REF], authors are looking for path selection policies that decreases the number of hops without affecting the bottleneck SNR. For this reason, several policies have been studied in this paper:

1) Highest-quality-first (HQF) policy: This policy is based on the calculation of the Signal-to-noise ratio (SNR) offered on link to each parent candidate at each hop in order to select the next Generation Node B (gNB) that offers the highest SNR. 2) Wired-first (WF) policy: This policy is designed to minimize the number of hops required to reach a wired gNB. At each hop, if there is a available link to a wired gNB from the current IAB node with an SNR higher than a threshold Γ th , then the wired gNB is selected, else the HQF policy is applied. 3) Position-aware (PA) policy: This policy is based on the use of IAB nodes and gNB position information. First, the IAB node splits into two half-planes the neighboring region, identified by the line, which is perpendicular to the one that passes through the IAB and the wired gNB positions. Then, for the selection only the parents candidate, which are in the half-plane with also the wired gNB are considered. Finally, the one with the highest SNR is selected. 4) Maximum-local-rate (MLR) policy: This policy is based on the calculation of the achievable rate. At each hop, the user selects the parent candidate with the highest achievable Shannon rate.

In [START_REF] Farhan | An approach for optimal base station selection in 5g hetnets for smart factories[END_REF] authors propose two basic policies that are not specifily dedicated to a backhaul context but can be used as reference policies like the RSRP-based policy :

1) Maximum Signal-to-interference-plus-noise ratio (SINR) association strategy: This policy is based on the calculation of the SINR. Each user is associated to the parent candidate for which it has the highest preliminary SINR. 2) Maximum receive power association strategy: This policy is based on the calculation of the receive power. Each user is associated to the parent candidate for which it has the highest receive power.

The policies proposed in [START_REF] Polese | Distributed path selection strategies for integrated access and backhaul at mmWaves[END_REF] and [START_REF] Farhan | An approach for optimal base station selection in 5g hetnets for smart factories[END_REF], do not take into account the network state in terms of capacity of the base stations, which can introduce an overload on stations and therefore packet losses while there are stations that have the capacity to support these packets. The successful exploitation of the network capacity requires taking into account the network state at the moment of cell selection.

In [START_REF] Ranjan | Backhaul-aware cell selection policies in 5G IAB networks[END_REF], authors propose several cell selection strategies. For the operation of the proposed strategies, the authors propose that the following information will be broadcast by the parent candidate nodes to assist UEs in the cell selection process:

• SNR Γ im offered on the link between the parent candidate i and the user m, refers to the link quality. • The number of hops between the parent candidate node i and the CN h i , refers to the latency. • The number of associated UEs to the parent candidate node i n i , refers to the load.

Achievable Rate of Path (ARP) policies: Each UE m has a parent node BS i, it encounter a set P i of radio links along its path to the CN. A radio link between a child node j (including the UE m) and its potential parent node k, the link normalized capacity is determined by C j = log 2 (1+Γ kj ). The calculating of the normalized capacity rate of each radio link along the path allows to the UE m to identify the achievable rate of the path R i . The UE will associate with the station that offers him the best R i :

i * = argmax i R i (1)
Several ways for the calculation of R i have been proposed:

• ARP using Minimum rate (AM) policy:

R i = 1 1 + h i min j∈pi {C j } (2)
• ARP using Harmonic mean of rates (AH) policy:

R i = 1 1 + h i HM j∈pi {C j } (3)
• ARP using Scaled Minimum rate (ASM) policy (extension of the AM policy):

R i = 1 1 + h i min j∈pi { C j n k } (4)
• ARP using Scaled Harmonic mean of rate (ASH) policy (extension of the AH policy):

R i = 1 1 + h i HM j∈pi { C j n k } (5)
In [START_REF] Ranjan | Backhaul-aware cell selection policies in 5G IAB networks[END_REF], the authors show that taking into account the network state information at the cell selection moment introduces a significant improvement on the network performances. However, the author's consideration of gNB load in terms of the number of associations has the potential to degrade the quality of the decision making. First, the number of associations does not reflect the real load of a parent candidate, i.e. two candidates with the same number of associations do not necessarily have the same load. Second, if the UE has two parent candidates, the first with a very good link capacity and a large number of associations (knowing that it still has the capacity to serve new user) and the second with a low link capacity but it is empty. In this case we can choose to associate with the second candidate, hence the policies ASH and ASM do not always guarantee the best quality of link when it is feasible. Third, there is no guarantee that the selected parent is not saturated.

III. LOAD-EFFICIENCY-BALANCE POLICY

Load-Efficiency-Balance (LEB) is our proposed cell selection policy for IAB networks. LEB aims to achieve a tradeoff between spectral efficiency and load balancing in the IAB context defined in the 3GPP standards. This policy allows each new user joining an IAB network to determine its suitable parent by applying the algorithm summarized by the flow chart presented on Fig. 2.

Fig. 2 -LEB algorithm flow chart

Using LEB, each new user joining the network will first calculate the cost of transmitting a bit through each parent candidate i. This transmission cost is determined by the following formula (expressed in terms of resource block (RB)):

Cost i = j∈Pi 1 W rb × SE j (6)
where P i refers to the set of radio links (access and backhaul links) along the path of the parent node i to the CN, W rb is bandwidth of a RB (expressed in Hz/RB) and SE j is the spectral efficiency of the link j (expressed in Bit/s/Hz) determined by Shannon capacity:

log 2 (1 + SN R).
SNR is a measure that compares the level of a desired signal to the level of background noise:

SN R = P r N 0 × W × N F (7)
where P r is the received signal power (expressed in W), W is the bandwidth (expressed in Hz), N 0 is the noise spectral density and N F is the receiver noise figure.

SINR is a measure of the signal quality: the strength of the wanted signal compared to the unwanted interference and noise.

SIN R

= P r N 0 × W × N F + j∈N P r j ( 8
)
where N is the set of UEs share one or more resources with the current UE. RX P W R represent the received signal power (expressed in dBm) and can be expressed as (see in [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA) radio frequency (rf) system scenarios[END_REF]-section 4.5):

RX P W R = T X P W R -max(pathloss-G T X -G RX , M CL) (9)
where T X P W R is the transmitted signal power, G T X is the transmitter antenna gain = G ant + A(theta) where G ant is the antenna gain and A(theta) is the antenna gain per sector defined in [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA) radio frequency (rf) system scenarios[END_REF]-section 4.2, G RX is the receiver antenna gain and M CL refers to the minimum coupling loss, which is the parameter describing the minimum loss in signal between BS and UE or between two UEs in the worst case. M CL is defined as the minimum distance loss including antenna gains measured between antenna connectors. The "capacity condition" in the algorithm, consists in checking whether or not the parent candidate has the capacity to support a new UE m:

∀j ∈ P i , B j (1 -BO j ) > R m × T f T s × W rb × SE j ( 10
)
where B j is the total number of resource blocks available in a 10-ms frame at the parent node of link j (expressed in RB), BO j is the bandwidth occupancy rate i.e the ratio of allocated RBs to the total available RBs, R m represents the desired rate that has to be transmitted for the user m (expressed in Bit/s), T f is frame duration (expressed in s) and T s is slot duration (expressed in s), considering the allocation of one RB per slot to the same UE. For the operation of the LEB policy, we propose a diffusion of information on the load from BSs to UEs, expressed by a value BO j between 0 and 1. This value will limit the load on the BSs, by preventing overuse of the bandwidth.

IV. SYSTEM MODEL

In this section, we provide the necessary details on the system model used for performance evaluation.

A. Simulation scenario

As a test scenario, a heterogeneous scenario (dense urban) extension as shown in Fig. 3 proposed by 3GPP [START_REF]Study on integrated access and backhaul[END_REF] is adopted. In our scenario an IAB donor is surrounded by 3 hexagonal sectors where each sector contains 3 IAB nodes. The IAB nodes have the IAB donor as a parent. The IAB donor is equipped with 3 antennas, one in each sector. Table I summarizes the main parameters used in this scenario. The simulation results discussed in the next section are obtained by 1000 independent runs, where each execution is running for 100 iterations to reduce interference's impact. Table II summarizes the main system level parameters derived from [START_REF]Study on integrated access and backhaul[END_REF] [9] [10] [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA) radio frequency (rf) system scenarios[END_REF]. To model the path loss, we used the Uma-NLOS model [START_REF]Study on channel model for frequencies from 0.5 to 100 GHz[END_REF] on the direct and backhaul links and the Umi-NLOS model [START_REF] Series | Guidelines for evaluation of radio interface technologies for IMT-Advanced[END_REF] on the access link.

B. Half-duplex implementation

In IAB networks, a base station cannot receive and transmit at the same time due to the half-duplex constraint, to avoid intra-site interference as direct, access and backhaul links share the same spectrum. In the standards there are several techniques to respect this constraint. In our simulations, we considered the Time-division multiplexing (TDM) technique that allows a strict time-domain separation between transmission and reception phases. With this technique, we will have two phases for the downlink and two phases for the uplink as shown in Fig. 4. To classify the macro-UE traffic between phases 1 and 2 we use SNR as creteria. The macro-UE communications with the best SNR are assigned to phase 2 in order to ensure that UEs far from micros are managed on phase 2 to minimize interference during this phase. To evaluate the performance of the proposed policy, we also implemented two policies from the literature presented in the section II: the HQF policy as a simple benchmark policy and the ASH policy since in [START_REF] Ranjan | Backhaul-aware cell selection policies in 5G IAB networks[END_REF] the authors showed that, among their proposed policies, ASH has the best performance. This section represents the main performance evaluation results of the comparison between the LEB and these two policies. Fig. 5 shows the achievable rate as a function of the number of users, considering the packet losses due to interference and capacity. The results show that, at low loads below 160 users, there is a slight difference between the three policies in terms of the achievable rate due to the impact of interference. This impact increases as the number of users grows. Past 160 users, the LEB policy outperforms the HQF and ASH and achieves the best system rate. LEB is more resistant to packet loss due to capacity. This is due to the fact that with LEB, the user chooses the best parent candidate only if the "capacity condition" is validated, if not, another worse parent is chosen instead of losing the packets. Unlike with HQF and ASH, where only the best parent calculated by the policy is directly chosen. Therefore, they assign users to congested cells at high load, which does not allow to take advantage of the residual capacities of the surrounding cells. Fig. 6a indicates the necessary consumption of resource block to reach the total transmission of data that the system aims to transmit as a function of the number of users. On the other hand, it needs less additional resources than HQF and ASH. Indeed, after 220 users, some cells are full and LEB makes users assignments to second choice with lower spectral efficiency but with available resources. Fig. 7a and Fig. 7b present the number of associations of users to parents (macro or micro node) and the total cost of bit transmission given by (6) (on the direct link for a macro parent and on the access plus backhaul links for a micro parent) with these associations as a function of the distance between parent and user, respectively. The results show that, at 160 users, when all three policies provide the same total number of associations (no packet loss due to capacity) LEB always chooses the associations with the lowest cost compared to HQF and ASH, even if it means longer distance in some cases. For example, by choosing to associate with the macro, which is a bit far rather than with the micro, since the cost of the direct link is lower than the cost of the access link plus the cost of the bakchaul link. HQF chooses the micro, which is closer and thus offers a better SNR without taking into account the cost of the backhaul. ASH chooses the micro, Fig. 8a and Fig. 8b illustrate the percentage of requested macro resources during phases 1 and 2 as a function of sector number at 280 users (when the system is losing data due to capacity with all policies), respectively. The results show that, first, LEB requires much more resources for downlink communications on the direct link between phases 1 and 2 than HQF and ASH. This means that LEB did more direct associations with the macro than HQF and ASH. In some cases for example, due to the backhauling cost, the cost of associating with the macro, that is a bit further away and more loaded than micros is lower than the cost of associating with a empty micro, that is closer and offers a better link I). With LEB on sector 1 the most loaded when the macro is saturated, the rest of the users of this sector choose to associate with micros on sectors 2 and 3 as long as the macro has free resources for the backhaul on sectors 2 and 3 unlike HQF and ASH.

VI. CONCLUSION

In this work, a new cell selection strategy, LEB, for IAB networks was proposed. LEB aims to guarantee the lowest possible total transmission cost (including the backhaul cost, which is a significant factor in an IAB context) while respecting the base stations capacities. In order to validate the performance of LEB, a comparative study with two state-ofthe-art policies using different metrics was performed. The results show that, on the one hand, LEB guarantees the best total transmission cost in terms of consumed resources. On the other hand, at high network load, it offers the best system capacity in terms of served users thanks to a good management of the resource usage. To conclude, LEB ensures a trade-off between spectral efficiency and load balancing, which gives the advantage over other policies to achieve our main goal,

Fig. 1 -

 1 Fig. 1 -An example of RSRP-based policy operation

Fig. 3 -

 3 Fig. 3 -3GPP heterogeneous scenario [8] TABLE I Scenario parameters Number of macro BS 1 Number of micro BS 3/sector Inter-Macro Distance 500 Min distance between Micro-UE [m] 10 Min distance between Macro-UE [m] 35 Min distance between Macro-Micro and Micros [m] 40 Number of UEs Several values Percentage of the number of UEs/sector [50,20,30] Position of UEs Random Message size [Mbit] 2 Frame size [ms] 10

Fig. 4 -

 4 Fig. 4 -Downlink and uplink communication phases

Fig. 5 -

 5 Fig. 5 -Achievable Rate as a function of number of users

Fig. 6 -

 6 Fig. 6 -RBs consumption as a function of number of users

Fig. 7 -

 7 Fig. 7 -The functioning of policies at 160 users quality. Second, the results show that with LEB the macro is almost saturated on the 3 sectors. When with HQF and ASH, there are still resources available on sectors 2 and 3 (in our scenario sectors 2 and 3 are less loaded than sector 1 see tableI). With LEB on sector 1 the most loaded when the macro is saturated, the rest of the users of this sector choose to associate with micros on sectors 2 and 3 as long as the macro has free resources for the backhaul on sectors 2 and 3 unlike HQF and ASH.

(a) phase 1 (b) phase 2 Fig. 8 -

 128 Fig. 8 -Percentage of requested macro resources as a function of sector number at 280 users boosting IAB network capacity. As a continuation of this work, we envisage enhancing LEB by optimizing its load balancing feature and, possibly, by managing energy at low load.


Baptiste Oudon

Chanez Belaidi

Alice Prosperi

Jessica Salaün

come

Xflow modeling for investigation of fluid structure interaction of artificial reef: application to burial effect.

Xflow modeling for investigation of fluid structure interaction of artificial reef: application to burial effect.

Baptiste Oudon a ; Chanez Belaidi b ; Alice Prosperi a ; Jessica Salaun c ; Arnaud Coutu a a GéoLab, Institut Polytechnique UniLaSalle, Beauvais, France b Institut Galilée, Villetaneuse; France c Atlantique Landes Récifs, Capbreton, France

The rheology data was obtained from literature (Zhou et al., 2014) and defined the behavior of a saturated sand as a viscous fluid by the power Law. The viscosity model chosen was Herschel-Buckley, which is an extension of the power law (Kang and Kim., 2016), because this viscosity model defines the behavior of a non-Newtonian fluid at threshold. Two types of multiphase modeling were performed, the first (Figure 3) having the same boundary conditions as the bowl simulations. The other (figure 4) having the boundary conditions of a strong current. The modeling figure 3 have shown similar results to the bowl Figure 1, the same transport process was observed. The geometry of the reef was mobile during the first second of simulation, the reef buried of 0.05 m and then stabilized. The simulation (Figure 4) shows a high turbulent kinetic energy which is explained by the high vorticity. The vorticity velocity correlated with the turbulent kinetic energy, induces the transport of a lot of sediment. A basin is created at the depth of 0.12m to 0.54s. Sediment transport with high turbulent energy seems to be consistent with the results of (Mazzuoli et al., 2020).

The purpose of this study is to determine the main cause of the burial of artificial reefs. For this purpose, two types of fluid structure seabed interaction modeling (FSSI) have been performed. A modeling of FSSI of the seabed geomorphology modeled in 3D for the years 2019, 2020 and 2021 was done. The purpose of this work was to determine the reasons of the changing geomorphology through computational fluid dynamics (CFD).

Introduction

The threshold of sediment transport in the ocean environment was evaluated by the shield parameter (Li et al., 2020), and this parameter was compared with the turbulence intensity by (Mazzuoli et al., 2020). The turbulence intensity at ground was displayed in percentage and was calculated using the equation (1) (Dassault Systèmes., 2021):

𝐼 = 2 3
𝐾 𝑈 Where U is the mean velocity and k is the turbulence kinetic energy.

The equation of the turbulent intensity shew that the maximum turbulence intensity, in correlation with the mean velocity, indicated an increase of turbulent kinetic energy as the bowl was increasing. K went from 2.1 J for the 40 % intensity zones to 11 J for the maximum turbulence intensity. (Mazzuoli et al., 2020) have shown that a high turbulent kinetic energy correlated with a very low velocity (< 0.05 m. s -1 ) creates a low turbulence, which creates slow sediment transport. Thus, sand grains roll and slide on the ground. A bigger sediment transport could be created for the same power of turbulent kinetic energy but for a higher speed.

Simulation of the interaction fluid structure on the bow Multiphase Modeling

Figure 1 :Figure 3

 13 Figure 1: Evolution of the turbulence intensity during the year 2019 to 2021 Figure 2: Evolution of the mean velocity during the year 2019 to 2021

Figure 4 :

 4 Figure 4: Interaction turbulence with seabed

The second modeling was a validation of the prediction tool for FSSI modeling, considering the behavior of the surface part of the saturated sand as a pseudo viscous fluid(Zhou et al., 2014)


Johan Aubray

Stéphane Florence Nicol

Stéphane Puechmorel

Geodesic regression on SE(3) and application to estimation of positions

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

Geodesic regression on SE(3) and application to estimation of positions

Johan Aubray 1 Florence Nicol 1 Stéphane Puechmorel 1 1 Ecole Nationale de l'Aviation Civile, Toulouse, France, e-mail: firstname.surname@enac.fr

Our application

Estimating the true position of a mobile object (aircraft, drone, ...) from measured (therefore noisy) raw positions Two cases :

From past positions: prediction : frequently solved with linear / extended / unscented Kalman filters, or particle filters [START_REF] Bar-Shalom | Estimation with applications to tracking and navigation[END_REF] From past and future positions: smoothing This may be difficult with a highly maneuvering target [START_REF] Pilte | Dynamic management of tracking ressources for hyper-manoeuvring targets[END_REF]

Our framework

Instead of considering the trajectory in the state space (R 3 ), the movement of the mobile object is considered as a sequence of displacements: rotations and translations in R 3 This set is called the Special Euclidean group SE(3) SE( 3) is a Lie group, which is both: a differentiable manifold a group in which both the product and inverse maps are smooth

Underlying Lie group: SE(3)

An element of SE(3) can be represented by a 4-by-4 matrix

R T 0 1 , R ∈ SO(3), T ∈ R 3 ∶
An element of its Lie algebra se(3) can be represented by a 4-by-4 matrix:

Ω U 0 0 , Ω skew-symmetric, U ∈ R 3

Our method

A minimizing geodesic between two points g 0 , g 1 with respect to ⟪⋅, ⋅⟫ g is a curve γ∶ [0, 1] → G, such that γ(0) = g 0 , γ(1) = g 1 minimizing the functional:

L(γ) = 1 2 1 0 ⟪ γ(t), γ(t)⟫ g dt (1)
Left invariance of the metric implies it is enough to seek after a solution of the reduced problem:

L(γ) = 1 2 1 0 ⟨ ξ(t), ξ(t)⟩ dt ( 2
)
where ξ is a curve in g. ξ satisfies the so-called Euler-Arnold equation [START_REF] Marsden | Introduction to Mechanics and Symmetry: A Basic Exposition of Classical Mechanical Systems[END_REF]:

ξ(t) = ad † ξ(t) ξ(t) (3)
An extended notion of parametric linear regression on Riemannian manifolds, called intrinsic geodesic regression, has been independently developed in [START_REF] Fletcher | Geodesic Regression on Riemannian Manifolds[END_REF] and [START_REF] Niethammer | Geodesic regression on image time series. In Medical image computing and computer-assisted intervention : MICCAI[END_REF] Let y 1 , . . . , y N be a set of points representing the trajectory of the mobile on a Riemannian manifold (M, g): the idea is to find a geodesic curve γ(t) on the manifold that best fits those points in a "linear" manner at known times t 1 , . . . , t N The geodesic regression model is expressed as follows in [START_REF] Fletcher | Geodesic Regression on Riemannian Manifolds[END_REF] Y = Exp Exp(p, tv), ε

where γ(t) = Exp(p, tv) is the geodesic curve given by the initial conditions γ(0) = p and γ(0) = v, and ε is a Gaussian random variable taking values in the tangent space at γ(t).

In [START_REF] Hinkle | Polynomial regression on riemannian manifolds[END_REF], this geodesic regression model is generalized to a polynomial one where the curve γ(t) is a Riemannian polynomial of order k, as defined by

∇ γ(t) ∇ γ(t) . . . ∇ γ(t) γ(t) = ∇ k γ(t) γ(t) = 0,
where ∇ denotes the Levi -Civita connection The estimate is found by minimizing a least squares criterion based on a Riemannian distance between the model γ(t) and the data

E(γ) = 1 N N j=1 d(γ(t j ), y j ) 2
under the constraint that γ shall be a geodesic One can prove, by taking the differential of E with respect to the initial conditions that it is equivalent to cancelling out the Riemannian gradient

Figure 1 .

 1 Figure 1. A straight line Figure 2. An helix


Chinnaraji Annamalai
email: anna@iitkgp.ac.in

Upper Limits for Velocity, Momentum, and Energy of Motion

Keywords: classical mechanics, electromagnetic energy, light waves, speed of light

The velocity of a moving object with its non-zero mass or a moving particle with its non-zero mass is less than the speed of light because light has no mass. So, the speed of light is the upper limit for all moving objects in the Universe. In this article, the upper limits for momentum and energy of motion are derived for scientific problem-solving.

Introduction

The Sun is the main source of natural light on the Earth. Light waves are generally electromagnetic and do not require a medium for their propagation. Thus, the light can travel in a vacuum. Quanta of light are the smallest discrete packets of electromagnetic energy. A quantum of light that has no mass is called a photon. The relativistic motion [START_REF] Annamalai | A Novel Mass-Energy Equation from Lorentz Factor and Energy of Motion[END_REF][START_REF] Annamalai | Energy-Momentum Equivalence[END_REF][START_REF] Annamalai | Momentum-Velocity Equivalence[END_REF][START_REF] Annamalai | Energy-Work Equivalence[END_REF][START_REF] Annamalai | Lorentz Factor and Time Dilation on the Special Theory of Relativity[END_REF] of a particle is concerned with the motion of the particle whose velocity approaches the speed of light.

Velocity-Momentum Relation

In classical mechanics, a variable-mass object is an object whose mass varies with time. A rocket, a moving object that burns completely at regular time intervals, and an object in motion that loses its mass at every unit of time are examples of variable-mass. In this section, the author introduces definitions of momentum-velocity equivalence [START_REF] Annamalai | Energy-Work Equivalence[END_REF][START_REF] Annamalai | Lorentz Factor and Time Dilation on the Special Theory of Relativity[END_REF][START_REF] Annamalai | A New Mass-Energy Equivalence from Lorentz Factor and Energy of Motion[END_REF][START_REF] Annamalai | Momentum-Velocity Equivalence[END_REF][START_REF] Annamalai | Computation of Mass-Energy Equation from Lorentz Factor and Kinetic Energy[END_REF][START_REF] Annamalai | Mass-Energy Equivalence: Light Energy[END_REF].

Definition for constant-mass:

The momentum (P) of a moving object with its constant-mass (m) is directly proportional to its velocity (v).

𝑃 𝛼 𝑣 ⇒ 𝑃 = 𝑚𝑣, where 𝑚 is constant.

Here, 𝑑𝑃 𝑑𝑡 = 𝐹 = 𝑑(𝑚𝑣) 𝑑𝑡 = 𝑚 𝑑𝑣 𝑑𝑡 Definition for variable-mass: Suppose a moving object loses its mass (m) at every unit of time and increases its velocity (v) proportionately with its mass (m). Then, the momentum (P) of the moving object is directly proportional to its velocity (v).

𝑃 𝛼 𝑣 ⇒ 𝑃 = 𝑚𝑣 = 𝑥 𝑖 𝑣 𝑖 at the unit 𝑡 𝑖 of time for 𝑖 = 1, 2, 3, 4, ⋯ , where 0 < 𝑥 < 𝑚.

Then, 𝑑𝑃 𝑑𝑡 = 𝐹 = 𝑣 𝑑𝑥 𝑑𝑡 + 𝑥 𝑑𝑣 𝑑𝑡 = 𝑑(𝑚𝑣) 𝑑𝑡

The above equation can be written as follows for our understanding:

𝑑𝑃 𝑑𝑡 = 𝐹 = 𝑣 𝑑𝑚 𝑑𝑡 + 𝑚 𝑑𝑣 𝑑𝑡 = 𝑑(𝑚𝑣) 𝑑𝑡 = 𝑚 𝑑𝑣 𝑑𝑡 (1)

Kinetic Energy

Kinetic energy is a fundamental concept in classical mechanics that quantifies the mechanical energy performed by an object due to its motion.

The equation of kinetic energy is given below:

𝐸 = 1 2 𝑚𝑣 2 , ( 2
)
where m is the mass and v is the velocity.

By differentiating the equation of kinetic energy with respect to time t, we obtain

𝑑𝐸 𝑑𝑡 = 1 2 (𝑚2𝑣 𝑑𝑣 𝑑𝑡 )
By simplifying the above differential equation, we receive 𝑑𝐸 𝑑𝑡 = 𝑚𝑣 𝑑𝑣 𝑑𝑡 where the mass (m) of an object in motion is constant. From the equation ( 1), we conclude that 𝑑𝐸 = 𝑚𝑣𝑑𝑣 = 𝑣𝑑(𝑚𝑣).

(3)

Derivation of the Lorentz Factor

The Lorentz factor is the factor by which time, length, and mass change for an object moving at a speed close to the speed of light. By the following experiment, the author of this article derives the equation of the Lorentz factor [START_REF] Annamalai | Lorentz Factor and Time Dilation on the Special Theory of Relativity[END_REF][START_REF] Annamalai | Mass-Energy Equivalence: Light Energy[END_REF].

Figure 1. depicts the relativistic speed of the metal tube (AD). Let's consider that the metal tube shown in the Figure 1 is at a stationary position and the light takes time t to travel from the point A to the point D. The length of the metal tube is measured as follows: AB = BC = h = ct, where c is the speed of light and t is the time taken to move from the point A to the point D.

Suppose the metal tube with mass m and velocity v moves from the point A to the point B horizontally and the light travels from the point A along with the moving metal tube. When the light reaches the point C at time 𝑡′, the metal tube moves to point B at time 𝑡′. 𝐴𝐵 = 𝑣𝑡 ′ and 𝐴𝐶 = 𝑐𝑡 ′ According to the definition of the Pythagoras theorem, we can derive the Lorentz factor as follows:

𝐴𝐶 2 = 𝐴𝐵 2 + 𝐵𝐶 2 (4)
By substituting 𝐴𝐵 = 𝑣𝑡 ′ and 𝐴𝐶 = 𝑐𝑡 ′ in the equation ( 4), we obtain

(𝑐𝑡′) 2 = (𝑣𝑡′) 2 + (𝑐𝑡) 2 𝑡′ 2 (1 - 𝑉 2
𝐶 2 ) = 𝑡 2 By simplifying the above equation, we get

𝑡 ′ = 𝑡 √ 1 - 𝑣 2 𝑐 2 𝑡′ 𝑡 = 1 √ 1 - 𝑣 2 𝑐 2 𝛾 = 1 √ 1 - 𝑣 2 𝑐 2 (5)
Hence, 𝛾 is the Lorentz factor.

Mass-Energy Relation

The equation of relativistic mass is m=𝑚 0 𝛾, where the Lorentz factor 𝑖𝑠 𝛾 = 1

√ 1 - 𝑣 2 𝑐 2 Now, m = 𝑚 0 √ 1 - 𝑣 2 𝑐 2 (6)
Squaring both sides of the equation ( 6):

𝑚 2 = 𝑚 0 2 1 - 𝑣 2 𝑐 2 ⟹ 𝑚 2 (𝑐 2 -𝑣 2 ) = 𝑚 0 2 𝑐 2
By simplifying the above equation, we obtain

𝑚 2 𝑐 2 -𝑚 2 𝑣 2 = 𝑚 0 2 𝑐 2 (7)
By differentiating the equation ( 7) with respect to time, we get

𝑐 2 2𝑚 𝑑𝑚 𝑑𝑡 -[2𝑚𝑣
Substituting the equation ( 8) in the differential equation ( 2) of kinetic energy:

𝑑𝐸 = 𝑐𝑑(𝑐𝑚) = 𝑣𝑑(𝑚𝑣), ( 9
)
where 𝑐 is a numeric constant and also one of the values of 𝑣.

Therefore, 𝑣 = 𝑐.

(10)

By substituting the equation (10) in the equation ( 1), we conclude that

𝐸 = 1 2 𝑚𝑐 2 .
(11)

Energy-Work Relation

The equations of displacement (𝑠), initial velocity (𝑢), final velocity (v), acceleration (𝑎), force (𝐹), work (W), and energy (E) of an moving object with mass (𝑚) and time (𝑡) are shown below:

𝑎 = 𝑣 -𝑢 𝑡 ⇒ 𝑣 = 𝑢 + 𝑎𝑡 ⇒ 𝑑𝑠 𝑑𝑡 = 𝑢 + 𝑎𝑡 ⇒ 𝑑𝑠 = 𝑢𝑑𝑡 + 𝑎𝑡𝑑𝑡 ⇒ 𝑠 = 𝑢𝑡 + 𝑎 𝑡 2 𝑡 𝑣 2 = (𝑢 + 𝑎𝑡) 2 ⇒ 𝑣 2 = 𝑢 2 + 2𝑎 (𝑢𝑡 + 𝑎 𝑡 2 𝑡
) ⇒ 𝑣 2 = 𝑢 2 + 2𝑎𝑠 and 𝐹 = 𝑚𝑎.

Let us derive the energy-work equivalence [6-10] from the velocity equation. 𝑣 2 = 𝑢 2 + 2𝑎𝑠, (12) where u is the initial velocity, v is the final velocity, a is the acceleration, and s is the displacement of a moving object. Now, let us consider the velocity of a moving particle approaching the speed of light. Then, 𝑢 = 0 and 𝑣 = 𝑐. (13) By substituting the equation (13) in the equation (12), we obtain

𝑐 2 = 0 + 2𝑎𝑠 ⇒ 𝑐 2 = 2𝑎𝑠 ⇒ 𝑎𝑠 = 1 2 𝑐 2 ⇒ 𝑚𝑎𝑠 = 1 2 𝑚𝑐 2 ⇒ 𝐹 × 𝑠 = 𝑊 = 1 2 𝑚𝑐 2 , ( 14
)
where W denotes the work .

From the equation ( 11) and ( 14), we conclude that

𝐸 = 1 2 𝑚𝑐 2 , ( 15
)
The equation (15) denotes the mass-energy equivalence.

Upper Limits

The relativistic motion is concerned with an object whose velocity approaches the speed of light. The relativistic momentum (𝑃) of a particle with its non-zero mass (𝑚) cannot be equal to 𝑚𝑐. 𝑖. 𝑒. 𝑃 ≠ 𝑚𝑐 because 𝑐 is the speed of light with zero mass. Therefore, 𝑃 = 𝑚𝑐 𝑣 , ( 𝑐 𝑣 < 𝑐 ), is the upper limit for momentum.

We know that the speeds of subatomic particles with mass are less than the speed of light. For example, neutrino speed is not faster than the speed of light because it has mass.

Since 𝑃 ≠ 𝑚𝑐, the relativistic energy (𝐸) ≠ 𝑚𝑐 2 .

If 𝑃 = 𝑚𝑐 𝑣 , then 𝐸 = 1 2 𝑚𝑐 𝑣 2 < 1 2 𝑚𝑐 2 .

Hence, 𝐸 = 1 2 𝑚𝑐 2 is the upper limit for the energy of motion.

(17)

Conclusion

In this article, the upper limits for momentum and the energy of motion have been proved mathematically for scientific problem-solving in the relativistic mechanics. The relativistic mechanics is concerned with the motion of bodies whose velocity approaches the speed of light. This idea can enable the researchers for further involvements in the scientific research and development.

https://doi.org/10.33774/coe-2023-1bvs4 ORCID: https://orcid.org/0000-0002-0992-2584 Content not peer-reviewed by Cambridge University Press. License: CC BY 4.0


Cyril Meyer

Étienne Baudrier

Patrick Schultz

Benoit Naegel

Benoît Naegel

Combining max-tree and CNN for segmentation of cellular FIB-SEM images

Keywords: Mathematical Morphology, Connected Operators, Maxtree, Segmentation, Deep Learning, Convolutional Neural Network, Electron Microscopy, FIB-SEM

come

Introduction

The max-tree structure [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF][START_REF] Najman | Building the Component Tree in Quasi-Linear Time[END_REF] allows to perform efficiently attribute filtering [START_REF] Jones | Connected Filtering and Segmentation Using Component Trees[END_REF] and has been involved in many image processing applications. The resulting operators are called connected [START_REF] Salembier | Connected operators[END_REF] since they do not create new contours nor modify ⋆ IdEx Doctoral contract, Université de Strasbourg their position. In [START_REF] Cabrera | Segmentation of Axillary and Supraclavicular Tumoral Lymph Nodes in PET/CT: A Hybrid CNN/Component-Tree Approach[END_REF], Farfan et al. have suggested that max-tree attributes could be used to feed a deep convolutional neural network (CNN) in order to improve the results of detection and segmentation tasks. Following and generalizing this approach, the aim of this paper is to provide a reproducible framework enabling to perform various experiments involving max-trees and CNN in a context of semantic segmentation of cellular FIB-SEM images.

Our contributions are twofold: in a first approach, input images are preprocessed using various attribute filters [START_REF] Meyer | Automatic Multi Class Organelle Segmentation For Cellular Fib-Sem Images[END_REF] and then concatenated as additional inputs of a CNN. In a second approach, maps of attributes are computed from the max-tree and then added in a CNN, following Farfan et al. approach.

Finally, our work aims to be handy. For this purpose, all the methods we propose can be used on a high-end workstation and do not require large GPU/TPU clusters. Also, our source code, datasets (original images, annotations) and documentation are publicly available, allowing everybody to reproduce the results, but also to reuse the code for their own needs.

State of the art

To address the segmentation of cellular electron microscopy images, the state-ofthe-art methods are currently based on CNN [START_REF] Heinrich | Whole-cell organelle segmentation in volume electron microscopy[END_REF][START_REF] Liu | Automatic Reconstruction of Mitochondria and Endoplasmic Reticulum in Electron Microscopy Volumes by Deep Learning[END_REF][START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF][START_REF] Çiçek | 3D U-Net: Learning Dense Volumetric Segmentation from Sparse Annotation[END_REF][START_REF] Oztel | Mitochondria segmentation in electron microscopy volumes using deep convolutional neural network[END_REF][START_REF] Xiao | Automatic Mitochondria Segmentation for EM Data Using a 3D Supervised Convolutional Network[END_REF][START_REF] Meyer | Automatic Multi Class Organelle Segmentation For Cellular Fib-Sem Images[END_REF] and the U-Net architecture remains mainly used. However, despite the good accuracy that can be obtained using these methods, the resulting segmentations can still suffer from various imperfections. In particular, thin and elongated objects such as endoplasmic reticulum can be disconnected and some parts may be distorted [START_REF] Meyer | Automatic Multi Class Organelle Segmentation For Cellular Fib-Sem Images[END_REF]. These effects may be the result of the context window that is fixed and too narrow in the first layer of the CNN, preventing to capture sufficient global information.

To overcome this, Farfan et al. [START_REF] Cabrera | Segmentation of Axillary and Supraclavicular Tumoral Lymph Nodes in PET/CT: A Hybrid CNN/Component-Tree Approach[END_REF] have proposed to enrich a CNN with attributes computed from the max-tree, enabling to capture at a pixel level, an information that may be non-local.

In the sequel of this paper, we will explore various strategies in order to incorporate max-tree attributes into CNN with the aim of potentially improving segmentation results.

Methods

Max-tree

Let I : E → V be a discrete, scalar (i.e. grayscale) image, with E ⊆ Z n and V ⊆ Z. A cut of I at level v is defined as: X v (I) = {p ∈ E|I(p) ≥ v}. Let C[X]
be the set of connected components of X. Let Ψ be the set of all the connected components of the cuts of I:

Ψ (I) = v∈V {C[X v (I)]}
The relation ⊆ is a partial order on Ψ . The transitive reduction of the relation ⊆ on Ψ induces a graph called the Hasse diagram of (Ψ, ⊆). This graph is a tree, the root of which is E. The rooted tree T = (Ψ, L, E) is called the max-tree of I, with Ψ, L, E being respectively the set of nodes, the set of edges and the root of T . The parent node of N , denoted P ar(N ), is the unique node such that: (P ar(N ), N ) ∈ L and N ⊆ P ar(N ) for N ̸ = E. The branch associated to a node is the set of its ancestors and is defined for a node N ∈ Ψ by: Br

(N ) = {X ∈ Ψ | X ⊇ N }.
In this work, the computation of the max-tree is based on the recursive implementation of Salembier [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF], and node attributes are computed during the construction of the tree. In the rest of this paper, we will focus on the following attributes which have been proposed in the literature:

-The height H is the minimum gray level of the connected component [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF].

-The area A is the number of pixels in the connected component [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF].

-The contour CT represents the number of pixels that have both a neighbor inside and outside the component [START_REF] Salembier | Morphological operators for image and video compression[END_REF]. -The contrast C is the difference between the max and minimum in the connected component [START_REF] Salembier | Antiextensive connected operators for image and sequence processing[END_REF]. -The complexity CP L represents the contour length |CT | divided by the area [START_REF] Salembier | Morphological operators for image and video compression[END_REF]. -The compacity (sometimes named compactness or circularity) CP A is the area divided by the square of the contour length |CT | 2 [START_REF] Salembier | Connected operators[END_REF]. -The volume V is the sum of the difference between the pixels values in the node and the node height [START_REF] Najman | Building the Component Tree in Quasi-Linear Time[END_REF]. -The mean gradient border M GB represents the mean value of the gradient magnitude for the contour pixels [START_REF] Cabrera | Segmentation of Axillary and Supraclavicular Tumoral Lymph Nodes in PET/CT: A Hybrid CNN/Component-Tree Approach[END_REF].

The tree attributes can be merged in order to compute an image, by associating to each pixel an attribute value computed from its corresponding nodes [START_REF] Cabrera | Segmentation of Axillary and Supraclavicular Tumoral Lymph Nodes in PET/CT: A Hybrid CNN/Component-Tree Approach[END_REF]. Each pixel p belongs to several nodes: the connected component N including p in the level-set X I(p) (I) and all the nodes belonging to its branch Br(N ). To associate a unique value to each pixel, different policies can be implemented: for example, by keeping the maximum, the minimum or the mean value of the attributes of the branch nodes [START_REF] Cabrera | Segmentation of Axillary and Supraclavicular Tumoral Lymph Nodes in PET/CT: A Hybrid CNN/Component-Tree Approach[END_REF].

In this work, we propose the following strategy. For each pixel p, the set of nodes belonging to the branch of p is retrieved, and only a subset of nodes having an attribute value in a certain range (given as a parameter) is kept. From this set, the node N best optimizing a certain stability criterion is kept. Finally, the value of p in the resulting image is set to the attribute value of N best . The resulting image is normalized in the range V = 0, 255 .

The criterion used to retrieve the optimal node is based on the concept of Maximally Stable Extremal Regions proposed by Matas et al. [START_REF] Matas | Robust wide-baseline stereo from maximally stable extremal regions[END_REF]. The idea is to retrieve the most stable regions based on the area variation between successive nodes since these regions represent salient objects of the image. For each node N ∈ Ψ , with N ̸ = E (i.e. different from the root), we define two stability attributes as follows:

∇ A (N ) = | A(P ar(N )) -A(N ) | |H(P ar(N )) -H(N )| • 1 A(N ) ∆ A (N ) = |∇ A (P ar(N )) -∇ A (N )|
where P ar(N ) defines the parent node of N .

Segmentation

We base our segmentation method on fully convolutional networks for semantic segmentation. We use attribute filtered images or max-tree attribute maps as additional input of our network. We feed the network with the original and the preprocessed images, concatenating them in the color (spectral) channels.

For model architectures, we use a 2D U-Net [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF], which is a reference for biomedical image segmentation. A 3D U-Net [START_REF] Çiçek | 3D U-Net: Learning Dense Volumetric Segmentation from Sparse Annotation[END_REF] could also be used, but the results are not necessarily better [START_REF] Srikrishna | Comparison of Two-Dimensional-and Three-Dimensional-Based U-Net Architectures for Brain Tissue Classification in One-Dimensional Brain CT[END_REF][START_REF] Nemoto | Efficacy evaluation of 2D, 3D U-Net semantic segmentation and atlasbased segmentation of normal lungs excluding the trachea and main bronchi[END_REF][START_REF] Zettler | Comparison of 2D vs. 3D Unet Organ Segmentation in abdominal 3D CT images[END_REF] and the computational cost of training the model is much higher. In a preliminary experiment, we have compared 2D and 3D models, using an equal number of parameters and same input size and the 2D U-Net perform equal, if not better as the 3D one.

Each block of the network is composed of convolutions with ReLU activation, followed by batch normalization [START_REF] Ioffe | Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift[END_REF] and a residual connection [START_REF] He | Deep Residual Learning for Image Recognition[END_REF], see Figure 1. We use a 50% dropout entering the deepest block to avoid over fitting. We always use padded convolution to maintain the spatial dimension of the output. The model starts with 64 filters at the first level, for a total of 32.4 millions parameters. We train our models by minimizing the Dice loss. For the binary segmentation task, we use classical Dice loss and for multi-class segmentation problems, we use a weighted mean of the loss for each class, with the same weight (W = 0.5) for our two classes. We note X the ground truth, Y the prediction, W the weight list and C the classe list. The ε term is used for stability when (X + Y ) = 0 and is set to 10 -4 .

L Dice (X, Y ) = 1 - 2 • X • Y + ε (X + Y ) + ε L DiceM ean (X, Y, W ) = 1 |C| c∈C W c • L Dice (X c , Y c )
The model is trained for 128 epochs, each epoch is composed of 512 batches, a batch is composed of 8 patches and a patch is a 256 × 256 × C subpart of the image, with C the number of channels. We use random 90 • rotation, horizontal and vertical flips for data augmentation on the patches. We trained our model using Adam [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF] optimizer with the following parameters: α = 0.001, β 1 = 0.9, β 2 = 0.999, ε = 1e-07.

Evaluation Metrics

To evaluate our models, we binarize the model prediction with a threshold at 0.5. To predict a slice, we use the whole slice to avoid the negative border effects of padding.

To evaluate our results, we use the F1-Score which is a region-based metric, and the average symmetric surface distance (ASSD) which is a boundary-based metrics. In fact, the F1-Score is equivalent to the Dice Score. We note respectively T P , F P and F N the cardinal of the sets of true positives, false positives and false negatives. X is the ground truth, Y the binary prediction and ∂X the boundary of X.

F1-Score = 2 × T P 2 × T P + F P + F N ASSD(X, Y ) = x∈∂X d(x, Y ) + y∈∂Y d(y, X) |∂X| + |∂Y | with d(x, A) = min y∈A ∥x -y∥ 2 .

Experiments

In this section, we test the improvement of the segmentation thanks to the addition of filtered images in the input. We compare the original image input with the enriched version. We also compare in the same time a multi-class segmentation model and two binary segmentations models. Finally, we repeat each of these configurations 11 times. In total, 363 models have been trained for this experiment, each training lasts about 5 hours with an NVIDIA GeForce RTX 2080 Ti. First, we define the filters we use as our experiment variable. Attributes are selected for their potential help for the segmentation, we list the selected filters in Table 1 (attribute maps strategy) and 2 (connected operators strategy), and renamed them for sake of simplicity.

Table 1. List of used attribute maps filters and their names. Before processing the image with a max-tree, we apply a low pass filter (9 × 9 mean filter).

Name Attribute Criterion Limit Contrast∆

A C ∆A 0 ≤ A ≤ 1024 2 Complexity ∆ A CP L ∆A 0 ≤ A ≤ 1024 2 Compacity ∆ A CP A ∆A CP A ≤ 50 Volume∆ A V ∆A 0 ≤ A ≤ 1024 2 MGB M GB M GB 0 ≤ A ≤ 1024 2

Data

We perform our experimentation on a stack of 80 slices from a 3D FIB-SEM image. Each slice has a size of 1536×1408. The image represents a HeLa cell and has an (x×y ×z) resolution of 5 nm × 5 nm × 20 nm. A ground truth is available on the stack for two kinds of organelles (i.e. cell subunits): mitochondria and endoplasmic reticulum. A default background class is affected to non-assigned pixel. An example slice with label is available in Figures 3 and4 in Section A.

Results

Figures 2a and2b depict the F1-score on the two classes of segmentation as box plots. Detailed mean and variation score are available in Table 3 and4 Baseline segmentation results Mitochondria are well segmented with a median F1-score up to 95% in multi-class segmentation and 94% in binary segmentation, which let a little possible improvement. Median F1-scores for reticulum up to 72% and 71% in multi-class and binary segmentation respectively. This thin organelle is indeed more difficult to segment.

Additional-input segmentation results

On the mitochondria, the additional inputs improve the results for 11 cases out of 20 and 17 out of 20 for the reticulum. The gain on the reticulum is very interesting since it is the more difficult to segment for the baseline setup. The following additional inputs improve the result in all the four tests (binary and multi-class segmentations on mitochondria and reticulum): Contrast ∆ A , Complexity ∆ A , Contrast β. Moreover, the whiskers show that only Contrast ∆ A and Contrast β have a good stability in this experiment. These two inputs are therefore good candidates to be additional inputs to improve segmentation results. On the contrary, Compacity ∆ A and MGB do not yield to any improvement.

Reproducibility

In this section, we present the required steps to reproduce the results we presented.

We follow the ACM definition of reproducibility: "The measurement can be obtained with stated precision by a different team using the same measurement procedure, the same measuring system, under the same operating conditions, in the same or a different location on multiple trials. For computational experiments, this means that an independent group can obtain the same result using the author's own artifacts." [1] For this purpose, our codes and datasets are publicly available. The project is split between three subprojects. First, the experimentation detailed documentation, training scripts, evaluation scripts, preprocessing script, results logs [START_REF] Meyer | CTAISegmentationCNN. GitHub repository[END_REF]. Second, the max-tree related functions and preprocessing script [START_REF] Meyer | ComponentTreeAttributeImage. GitHub repository[END_REF]. Finally, the dataset with images and annotations [START_REF] Meyer | DGMM2022-RRPR-MEYER-DATA[END_REF] The following information are also available on the documentation repository with more details.

Requirements A system with Ubuntu 18.04.6 (or compatible) with g++ and git installed. Python 3.6.9 with an environment including TensorFlow 2.6.2, NumPy, SciPy, scikit-image and MedPy.

Image preprocessing

-Prepare data for extraction with low pass filter. python 01_mean_filter.py -Extract attribute image from pre-processed images.

./build_bin_requirements.sh ./02_attribute_image.sh -Crop the image to the annotated area and construct a tiff stack. python 03_crop_roi.py

Network training and evaluation For the following commands, $ID is a unique identifier for the train, $INPUT is the folder containing the dataset, $OUTPUT is the folder containing the trained models and evaluation metrics, $DATASET select the dataset to use (in our case, binary or multi-class dataset), $SETUP select the experiment to run. An automation bash script is available on the repository to run all the 33 setups once.

-

Conclusion

In this paper, we have presented a detailed experimental setup to evaluate the use of additional input in a CNN based segmentation task. The additional inputs are attribute maps obtained from a max-tree representation of the image. The evaluation is made on segmentation tasks in the context of 3D electronic microscopy. If most of the additional inputs improve one segmentation task, two of them -namely Contrast ∆ A and Contrast β-improve all the tested segmentation tasks in terms of median F1-score and stability. Further than the segmentation results, the setup inspired from [START_REF] Cabrera | Segmentation of Axillary and Supraclavicular Tumoral Lymph Nodes in PET/CT: A Hybrid CNN/Component-Tree Approach[END_REF] has been entirely implemented in C++ and Python and is proposed in open access to make it reproducible. As a perspective of this work, the feature extraction method based on maxtree attributes presented in this paper could be used for other applications. For example, it would be interesting to compute the max-tree directly inside the model and to use the attributes images as a nonlinear filter. Also, the attribute maps could be used as feature maps for more simple and explainable classifier as random forests or even decision tree. Besides, the ∆ A attributes we defined could be used in an interactive segmentation setup, where a single pixel will allow to select an interesting node, selecting an object connected component over a background. Finally, we proposed here a max-tree based method, but an extension to the tree of shapes [START_REF] Monasse | Fast computation of a contrast-invariant image representation[END_REF] could be interesting and add more information to the image. images and V. Mallouh for providing the annotations. We acknowledge the use of resources of the French Infrastructure for Integrated Structural Biology FRISBI ANR-10-INBS-05 and of Instruct-ERIC. We acknowledge that this work is supported by an IdEx doctoral contract, Université de Strasbourg.

A.2 Example preprocessing visualization

Fig. 1 .

 1 Fig.1. U-Net architecture and its backbone block. The number in the box corresponds to the number of filters at this level for the convolutions blocks. The network is fed with the original and the preprocessed images, concatenating them in the color channels.

2 .

 2 Figures 5 to 14 depict the slice for each applied filter. We divide the stack into 3 sets: training (first 40 slices), validation (next 20 slices) and test (last 20 slices). The training set is used to train the network, the validation set to select the best model during the training and the test set to provide evaluation metrics.

  Figures 2a and 2b depict the F1-score on the two classes of segmentation as box plots. Detailed mean and variation score are available in Table 3 and 4 in Section A.1.

Fig. 2 .

 2 Fig. 2. The box shows the quartiles of the results, while the whiskers extend to show the rest of the distribution.

Fig. 3 . 4 .Fig. 7 .Fig. 9 .Fig. 13 .

 347913 Fig. 3. Original image Fig. 4. Label, mitochondria in gray, endoplasmic reticulum in white

Table 2 .

 2 List of used connected filters and their names. The "inverse" column is checked when filter is applied on the inverted image.

	Name	Attribute	Range	Inverse
	Contrast α	C	[10, 150]
	Contrast β	C	[10, 150]	✓
	Area α	A	[2 × 32 2 , 1024 2 ]
	Area β	A	[64 2 , +∞]
	Area γ	A	[16 2 , 512 2 ]	✓

  Train the networks python train.py $ID $INPUT $OUTPUT $DATASET $SETUP -Evaluate the networks python eval.py $ID $INPUT $OUTPUT $DATASET $SETUP BEST Result analysis and figures reproduction Since the output of each model evaluation is a comma separated values (CSV) file, the analysis of the results can be done using various tools. We propose to use a Jupyter notebook with Pandas and Seaborn, merging the CSV files into a single dataframe. An example analysis.ipynb notebook is provided on the GitHub, which we use to produce our result figures and tables.

Acknowledgements We acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements d'Avenir) and the CPER Alsacalcul/Big Data. We thank D. Spehner from the Institut de Génétique et de Biologie Moléculaire et Cellulaire for providing the

A Appendix

A.1 Results

In the following tables, the mean and deviation are computed over the 8 best model out of 11, selected using the F1-score on the validation set.


Charles Voivret
email: charles.voivret@univ-reunion.fr

Dimitri Bigot
email: dimitri.bigot@univ-reunion.fr

Garry Rivière
email: garry.riviere@univ-reunion.fr

A Method to Compute Shadow Geometry in Open Building Information Modeling Authoring Tools: Automation of Solar Regulation Checking

Keywords: shading factor, building information modeling, solar protection, shadow geometry, regulation checking, automation

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

Introduction

The incorporation of renewable energy sources is essential for attaining enduring sustainability [START_REF] Ang | A comprehensive study of renewable energy sources: Classifications, challenges and suggestions[END_REF]. Limitless solar energy is necessary to mitigate the construction sector's greenhouse gas (GHG) emissions [START_REF] Li | Solar energy utilisation: Current status and roll-out potential[END_REF]. There are two methods for converting radiative energy. The first is thermal transfer, where the energy is used for heating purposes or to generate hot water by passing it through the building envelope. The other method involves the absorption of solar radiation by solar panels, which then convert it into electric energy using photovoltaic (PV) systems [START_REF] Hayat | Solar energy-A look into power generation, challenges, and a solar-powered future[END_REF]. Optimizing both conversion systems is crucial to achieving optimal performance throughout the seasons, particularly in tropical climates [START_REF] Vassiliades | Building integration of active solar energy systems: A review of geometrical and architectural characteristics[END_REF][START_REF] Rabani | Achieving zero-energy building performance with thermal and visual comfort enhancement through optimization of fenestration, envelope, shading device, and energy supply system[END_REF][START_REF] Hwang | Identifying relative importance of solar design determinants on office building façade for cooling loads and thermal comfort in hot-humid climates[END_REF]. It is essential to consider the unique characteristics of each building, like shadows. Considering shadows is crucial to guarantee the dependable functioning of systems. While building designers may prioritize assessing the overall performance of a structure, such as energy consumption per square meter, they can also utilize more minor scale performance indicators, such as the shading factor, to evaluate specific components. This factor quantifies the proportion of energy that reaches the glass surface of a window, taking into account any potential barrier caused by the surrounding environment, compared to the energy that comes to the surface without any obstruction. Subsequently, it can be employed, specifically in regions with tropical climates, to assess the efficacy of shade devices. Similarly to the regulations followed in Indonesia and India [START_REF] Bhatia | An approach to calculate the equivalent solar heat gain coefficient of glass windows with fixed and dynamic shading in tropical climates[END_REF], the French regulation for tropical climate (RTAADOM) presents a parametric model and tables of Buildings 2023, 13, 3120 2 of 14 acceptable values for determining the shading factor of some building components such as walls or windows. These values are contingent upon the location of the building and the orientation of the windows [START_REF] Grosdemouge | Passive design in tropical climates: Key strategies implemented in a French certified sustainable neighbourhood[END_REF][START_REF] Garde | Implementation and experimental survey of passive design specifications used in new low-cost housing under tropical climates[END_REF][START_REF] Garde | Elaboration of thermal standards for french tropical islands.Presentation of the PERENE Project[END_REF][START_REF] Garde | Integrated building design in tropical climates: Lessons learned from the ENERPOS net zero energy building[END_REF].

Figure 1 shows the framework of this paper. The left part of the figure illustrates this point of view and highlights each stage's important input and output. The right part shows different implementation alternatives regarding the computation stages. The left part of Figure 1 shows that regulation based on the shading factor relies on three stages:

(1) evaluation of shadow geometric features such as area, [START_REF] Li | Solar energy utilisation: Current status and roll-out potential[END_REF] integration in a shading factor model to compute the value of shading factor, and, finally, (3) comparison to a reference value to check a regulation. This paper focuses on this framework of regulation checking thanks to a shading factor. While the second stage is nowadays straightforward, thanks to accessible numerical tools such as a spreadsheet or a scripting language, the computation of 3D shadow geometry features (first stage) remains a challenging task. The RTAADOM offers a specialized tool in a secure spreadsheet that combines all steps necessary to calculate shading factor values. The system utilizes a basic parametric geometric model of solar protection, requiring the user to manually map their building design onto the parameter space of the geometric model. Figure 2 shows both the geometric configurations modeled in the RTAADOM and a snapshot of the user interface. Figure 2a shows all the parameters used to parametrize the overhang protection. If vertical fins are added, Figure 2b shows all the supplementary parameters to parametrize them. The official and protected spreadsheet that implements this model can be found on the French government website [START_REF]Textes réglementaires et fiches d'application de la RTAADOM[END_REF] and a snapshot of the interface is shown in Figure 2c. While the user directly inputs most parameters to describe the solar protection geometry, the parameters ψ and φ correspond to the sun's position in the computation and are not exposed to the user. Real designs are generally more complex than the geometric model and the projection represents a critical step in the process. To overcome this point, in other words, to compute the shading factor with complex geometry, two alternative approaches are iden-tified. The first one is to use research tools that implement various approaches over the first two stages of the process. While very flexible and capable thanks to cutting-edge algorithms, by design, its purpose is not oriented through checking and relies on tailored inputs. Designers can use Building Energy Modeling (BEM) approaches and established and reliable tools like EnergyPlus. EnergyPlus incorporates a shading factor as an intermediate outcome in the overall performance assessment. Proficiency in software development is necessary to make any adjustments to this implementation based on many modeling assumptions. Therefore, it is essential to ensure that the shading factor used in BEM matches the one used in the software. Nevertheless, this alignment needs to be more consistently apparent, as exemplified by the instance of RTAADOM. Integrating Building Information Modeling (BIM) into building performance modeling is a significant step in sustainable design. This collaboration enhances the ability to use parametric design, allowing architects and engineers to quickly evaluate multiple potential solutions, leading to more knowledgeable and effective decision-making processes [START_REF] Casini | Chapter 5-Building performance simulation tools[END_REF]. In addition, the seamless integration of BIM with building rating systems can transform the green certification process by automating several components and enhancing a building's sustainability credentials. The BIM approach, known for extracting material quantities from models automatically, makes green building assessments easier and faster by expediting sustainability evaluations [START_REF] Seghier | A BIM-Based Method to Automate Material and Resources Assessment for the Green Building Index (GBI) Criteria[END_REF][START_REF] Marzouk | Green building system integration into project delivery utilising BIM[END_REF]. Furthermore, integrating BIM into evaluations of environmentally friendly construction dramatically improves project completion efficiency and data organization while supporting sustainability goals. BIM also has a crucial function in optimizing the efficiency of regulatory compliance and certification procedures. An example is the emphasized advantages of integrating Building Information Modeling (BIM) with energy simulation as part of the LEED certification procedure [START_REF] Ryu | A study on the LEED energy simulation process using BIM[END_REF]. Nevertheless, exporting the BIM model amplifies the likelihood of flaws inside the file, which could compromise the outcomes' precision. Another method involves directly processing BIM model data within the modeling software to avoid any data loss, although this necessitates using proprietary software. To tackle more extensive issues related to accessibility, some authors have proposed using tools that rely on the IFC format. The utilization of IFC files simplifies the implementation of external tools, especially in cloud environments. An instance of a BIM platform that operates in the cloud and includes a calculation tool to verify the Envelope Thermal Transfer Value (ETTV) calculation showcases the potential of this technique to improve accessibility and scalability while guaranteeing the precision of essential calculations [START_REF] Liu | Envelope thermal performance analysis based on building information model (BIM) cloud platform-Proposed green mark collaboration environment[END_REF].

Nowadays, Building Information Modeling (BIM) workflow and associated authoring tools have become standard practice in designing buildings. Beyond software, BIM is a way to structure and efficiently share the building data's complexity. Over the last two decades, the Industry Foundation Classes (IFC) data scheme has become the standard for open BIM. The Industry Foundation Classes (IFC) file format is a crucial Building Information Modeling (BIM) standard. IFC, developed by BuildingSMART, aims to enhance interoperability within the architectural, engineering, and construction (AEC) sector. The open data model, compliant with the Standard [18] , facilitates the transfer and dissemination of information often utilized in BIM workflows, surpassing the constraints imposed by proprietary software. An IFC file is a complete data container containing detailed metadata and geometric information about building elements. The file format encompasses specific information such as the geometric depiction, material characteristics, and spatial connections of building elements. It facilitates interoperability in AEC but true interoperability between BIM software using IFC has not been achieved. Nonetheless, based on open or closed data schema, BIM authoring tools should be the original source of information to compute the shading factor. In other words, the inputs of the three previously presented approaches (Dedicated Tool, Research Tool, and BEM) should be extracted and/or converted from BIM to a suitable data format. Due to their specificity, manual extraction prevails for the first two. Concerning BEM, automating the generation of the BEM model from BIM data is an active topic of research [START_REF] Ciccozzi | BIM to BEM for Building Energy Analysis: A Review of Interoperability Strategies[END_REF] but remains challenging due to the complexity of both data models. Even if automation works, the benefits of BIM workflow could be limited by any data flow without automatic feedback from BEM to BIM. To overcome that challenge, a solution is to embrace a total BIM approach where the needed shading factor is computed using a BIM authoring tool with native data. This is the primary goal of this paper.

The primary obstacle lies in calculating the shadow geometry, specifically within a BIM authoring tool. Although they may be seen, they lack practicality and must be recalculated as simple geometric entities. This research introduces a novel technique dubbed "solid clipping" that utilizes common elements of the Computer Aided Design (CAD) kernel found in BIM authoring tools. The method is easily comprehensible and uncomplicated to execute with either a visual or textual API. This paper demonstrates implementing a technique in openBIM using IFCOpenshell and Python. The method is then employed to calculate the RTAADOM shading factor in a BIM model. The results of the geometric reference model are discussed.

This paper is structured as follows. Section 2.1 provides information regarding the calculation of the shading factor, encompassing both Stage 1 and Stage 2. Section 2.2 entails an exposition of the solid clipping technique. Section 3 covers applying solid clipping to calculate the RTAADOM shading factor in BIM.

Methodology

This section provides a brief state of the art about the computation of shading factor. It also introduces the shading factor used in the RTAADOM that has to be checked on the BIM model. A new method to compute the BIM authoring tool's shadow geometry easily is introduced.

An Overview of Shading Factor Computation

Shadow Geometry Computation

Computation of shadow is a mature field in video games, virtual reality, and 3D modeling. Over the last few decades, numerous techniques have been created to, most of the time, render better and faster [START_REF] Akenine-Möller | Real-Time Rendering[END_REF], in other words, to optimize software and hardware to achieve realistic views for human eyes through pixels. The goal here is different. It is to compute the geometry of the direct shadow of a single light source for quantitative analysis needed in design and engineering. Unfortunately, standard BIM authoring tools or viewers compute and display such shadows but do not expose them as actionable objects for quantitative analysis. Intuitively, lighting simulation software would be a perfect candidate to compute shadow. Still, most of the effort is oriented to rendering a realistic view and predicting lighting conditions inside buildings [START_REF] Tzempelikos | The impact of shading design and control on building cooling and lighting demand[END_REF][START_REF] Montiel-Santiago | Sustainability and Energy Efficiency: BIM 6D. Study of the BIM Methodology Applied to Hospital Buildings. Value of Interior Lighting and Daylight in Energy Simulation[END_REF][START_REF] Shikder | Evaluation of four artificial lighting simulation tools with virtual building reference[END_REF][START_REF] Natephra | Integrating building information modeling and virtual reality development engines for building indoor lighting design[END_REF]. Moreover, direct shadow geometries are intermediate and unreachable computation results as a step of a building energy analysis workflow. In the same way, direct shadows are an intermediate result of Building Environmental Modeling (BEM) based on zone model [START_REF] Shin | Thermal zoning for building HVAC design and energy simulation: A literature review[END_REF]. Beyond the expertise needed to extract such results, these modeling approaches need to convert and augment BIM data to a suitable dedicated model [START_REF] Gao | Building information modelling based building energy modelling: A review[END_REF][START_REF] Sun | A review of the-state-of-the-art in data-driven approaches for building energy prediction[END_REF][START_REF] Farzaneh | Review of using Building Information Modeling for building energy modeling during the design process[END_REF]. Even if progress has been made in this task, the cost and expertise needed to build such a model to extract shadow are disproportionate.

In the context of solar radiation modeling, shadow computation methods can be classified into four categories: trigonometric (TgM), polygon clipping (PgC), pixel counting (PxC), and ray tracing (RT). Trigonometric models use parametric geometry of obstructions to analytically project them on the surface of interest [START_REF] Cascone | Calculation procedure of the shading factor under complex boundary conditions[END_REF]. While computationally efficient, they are limited in terms of geometry [START_REF] Elmalky | Computational procedure of solar irradiation: A new approach for high performance façades with experimental validation[END_REF]. The French regulation RTAADOM proposes such a model, implemented in a spreadsheet, to compute the shading factor for a parametric model of solar protection [START_REF]Outil de Calcul du Coefficient Cm[END_REF]. Polygon clipping (PgC) computes shadow as the intersection of polygons on the surface of interest [START_REF] Melo | Using a shading matrix to estimate the shading factor and the irradiation in a three-dimensional model of a receiving surface in an urban environment[END_REF]. The polygons are the projection of the faces of obstructing objects. This method is much more versatile but shows limitations with non-convex or holed shapes in standard implementations [START_REF] Rocha | Experimental validation and comparison of direct solar shading calculations within building energy simulation tools: Polygon clipping and pixel counting techniques[END_REF]. Despite these limitations, PgC is primarily implemented in the BEM tools. Pixel counting computes the number of lighted (oppositely shaded) pixels representing the surface of interest from the sun's viewpoint. While this method provides an estimation due to sampling through pixels, it is versatile and precise [START_REF] Rocha | Experimental validation and comparison of direct solar shading calculations within building energy simulation tools: Polygon clipping and pixel counting techniques[END_REF][START_REF] Rocha | Domus method for predicting sunlit areas on interior surfaces[END_REF] thanks to modern graphic hardware that can efficiently deal with high-resolution images. In order to achieve optimal results, it is necessary to apply this technique using a visual programming language like OpenGL, which consequently demands specialized expertise. Ray tracing involves the casting and propagating of numerous basic rays that eventually undergo reflection [START_REF] Wang | A novel and efficient method for calculating beam shadows on exterior surfaces of buildings in dense urban contexts[END_REF][START_REF] Robledo | From video games to solar energy: 3D shading simulation for PV using GPU[END_REF][START_REF] Erdélyi | Three-dimensional SOlar RAdiation Model (SORAM) and its application to 3-D urban planning[END_REF]. In the forward version, the rays are cast from the sun and those that hit the surface of interest can be used to evaluate the shadow area. In the backward version, the rays are cast from the surface of interest and the counts of the number of rays intercepted by the obstruction can be used to estimate shadows. As a sampling method based on randomness, precision should assessed [START_REF] Arias-Rosales | Shadow modeling in urban environments for solar harvesting devices with freely defined positions and orientations[END_REF]. One significant advantage is that the multireflection of rays can be modeled.

Shading Factor

For a given surface, the shading factor describes the ratio between incoming irradiance with obstruction I s and irradiance without I. The shading factor comprises two models: a sky model to describe the different irradiance components caught by a surface of interest and a shadow model to compute the shaded area. The sky model generally decomposes the total irradiance into direct, diffuse, and reflected components, and projects it on a tilted surface. The shadow model describes how each irradiance component is impacted by obstructions in front of the surface of interest.

Sky model decomposes the total irradiance I on a tilted surface in the direct I b and the diffuse I d ones for a given location and time on earth. Numerous models have been developed [START_REF] Antonanzas-Torres | Clear sky solar irradiance models: A review of seventy models[END_REF], varying in terms of physical and mathematical complexity to offer better precision in predictions [START_REF] Loutzenhiser | Empirical validation of models to compute solar irradiance on inclined surfaces for building energy simulation[END_REF]. Models mainly differ in how they decompose the diffuse component I d into different contributions: isotropic from the sky dome, circumsolar, horizon brightening, and reflected from the ground. This work will only use a basic model with isotropic diffuse radiation for both sky and ground [START_REF] Liu | The interrelationship and characteristic distribution of direct, diffuse and total solar radiation[END_REF][START_REF] Duffie | Solar Engineering of Thermal Processes[END_REF]. Such a model allows us to consider essential physical traits of the problem without much mathematical complexity. Moreover, the thermal regulation that has to be checked automatically is based on this model. The total irradiance I is then written as I = I b + I d + I r with:

I b = I h,b R b (1)
I d = I h,d 1 + cos β 2 (2)
I r = I h ρ 1 -cos β 2 (3)
where I h,b is the direct-normal component of solar irradiance on the horizontal surface, R b is a geometric factor, I h,d is the global diffuse horizontal solar irradiance, β is the tilt angle of the surface from the horizon, I h is the global horizontal solar irradiance, and ρ is the soil albedo.

To compute a shading factor, including diffused irradiance, one needs to consider the effect of shadow on diffuse component I d . One common approach is to consider the nondirect source of radiance as a sum of direct sources from different directions with associated direct shadows. In other words, the shaded diffuse component is computed via the integration of discretized shaded sources [START_REF] Elmalky | Computational procedure of solar irradiation: A new approach for high performance façades with experimental validation[END_REF][START_REF] Melo | Using a shading matrix to estimate the shading factor and the irradiation in a three-dimensional model of a receiving surface in an urban environment[END_REF][START_REF] Loutzenhiser | Empirical validation of models to compute solar irradiance on inclined surfaces for building energy simulation[END_REF][START_REF] Maestre | Influence of selected solar positions for shading device calculations in building energy performance simulations[END_REF]. Following [START_REF] Cascone | Calculation procedure of the shading factor under complex boundary conditions[END_REF], the instantaneous shading factor S f is defined for one solar position by the ratio between the irradiance with shadow and without:

S = I s I = F b I b + F d I d + F r I r I b + I d + I r (4)
where F b , F d , and F r are the geometric shading factors for direct, diffuse, and reflected irradiance components, respectively. For direct irradiance, F b is given by the ratio of the sunlit area and its total area: F b = A s /A. For diffuse component and a given hemisphere discretization, F d is given by:

F d = ∑ i ∑ j F b,ij .R ij . cos θ ij .δΩ ij ∑ i ∑ j R ij . cos θ ij .δΩ ij (5)
where i and j denote the discretized height and latitude coordinates of the considered portion of sky hemisphere, R ij is the radiance passing through the portion of hemisphere calculated according a sky model, θ ij is the angle between the sky element and the normal of the surface, and δΩ ij is the solid angle of the sky portion seen from the surface. A similar approach can compute a geometric shading ratio for the reflected component depending on the use case and the modeling hypothesis about ray reflection. While the instantaneous shading factor can be interesting to assess predictions minutely, the average overtime period is often used in application, S f t . Alternatively, the RTAADOM uses an integral form to define the shading factor:

S R = ∑ t I s t ∑ t I t (6)

A New Method to Compute Shadow Geometry: Solid Clipping

BIM authoring tools generally implement a Computer Aided Design (CAD) kernel suitable to model the shapes of building components using surface and solid modeling. The proposed method is based on essential solid modeling operations such as extrusion and intersection to compute shadow, inspired by the polygon clipping approach. It is an oversimplified version of the shadow volume technique [START_REF] Mccool | Shadow volume reconstruction from depth maps[END_REF], a well-known method capable of rendering complex scenes in terms of geometry and daylighting. Indeed, our goal is to compute shadow on one face F of the model for a unique direct light source.

As it can be seen in Figure 3, where the method is applied to cast shadows on a window, the shadow computation method consists of four steps:

1.

Identification of the building walls or elements that can cast a shadow regarding sun exposition, as in Figure 3a; 2.

Extrusion of a three-dimensional (3D) sun path in the direction of the wall where shadows have to be evaluated, as in Figure 3b; 3.

Calculation of 3D intersections of the sun path with the building to identify the parts that can generate shadows on the considered walls, as in Figure 3c; 4.

Extrusion of previous 3D intersections, reduced to intersection surfaces, up to the considered walls to trace and calculate the shadows, as in Figure 3d. As it can be seen in Figure 3, where the method is applied to cast shadows on a 241 window, the shadow computation method consists of four steps : Identification of the building walls or elements that can cast a shadow regarding sun 243 exposition, as in Figure 3a;

244 2.
Extrusion of a three-dimensional (3D) sun path in the direction of the wall where 245 shadows have to be evaluated, as in Figure 3b;

246 3.
Calculation of 3D intersections of the sun path with the building to identify the parts 247 that can generate shadows on the considered walls, as in Figure 3c; As it can be seen in Figure 3, where the method is applied to cas window, the shadow computation method consists of four steps :

1.

Identification of the building walls or elements that can cast a shado exposition, as in Figure 3a; 2.

Extrusion of a three-dimensional (3D) sun path in the direction of shadows have to be evaluated, as in Figure 3b; 3.

Calculation of 3D intersections of the sun path with the building to i that can generate shadows on the considered walls, as in Figure 3c;

(dark yellow), intersection with F gives F shadow as a patch of surface (dark blue).

This method has inspired the algorithm described in Algorithm 1.

Line 1

The -→ sun points downward, i.e., from the sun through the face of interest F. F is extruded to the infinite along --→ sun to create the "solid of light" Light. It represents the light flux without obstruction, as shown in Figure 3a.

Line 2

The Boolean intersection between the solid Light and all the other solid of the BIM model (Building) is computed and stored in CatchedLight. Shown in black in Figure 3c.

Line 5 to 8 Each face of CatchedLight with normal oriented through the sun (red in Figure 3d) is extruded back to F (dark yellow in Figure 3d) to build PartialShadowVolume. The face resulting from the intersection between PartialShadowVolume and F is computed and stored in PartialShadowSur f ace.

Line 10

The Boolean union of all the PartialShadowSur f ace (dark blue in Figure 3d) defines the face shadow F shadow .

The order of solid operation has been optimized to propose a robust algorithm. Indeed, one can think of uniting the extrusion to perform one unique intersection, thus saving time. It has been found that the intersection of a unique volume resulting from the union of multiple possibly overlapping extrusions is more error-prone than the union of multiple intersections on the face F. In other words, for the sake of robustness, it is better to deal with the most straightforward atomic geometry. The advantages and limitations of this method regarding application results will be discussed at the end of the paper.

Application: Checking the RTAADOM on BIM Models

The solid clipping method has been applied to compute the RTAADOM shading factor. Thanks to the IFCOpenShell library, a Python tool has been developed to implement the solid clipping method to compute the RTAADOM shading factor S bim R .

Validation of RTAADOM Shading Factor Computed Using BIM

In this section, a comparison of shading factor values given by the RTAADOM official Excel sheet S re f R with those computed using our automatic analysis of an IFC file S bim R has been made. Varying window configuration regarding protection geometry and orientation, S R values were quantitatively compared. The first model comprises four walls and 10 identical windows with solar protection hosted on the same linear wall. An overhang of different depths is associated with each window, as shown in Figure 4a. For each window, the overhang depth d increases from 0.1 m to 1.0 m by 0.1m steps. This BIM model, and thus these 10 protection configurations, are strictly geometrically equivalent to those managed by the parametric model of the RTAADOM (see Figure 2). The window's orientation can be controlled by setting the angle between the project north and the true north. The comparison of the S R values for a window oriented to the four cardinal orientations can be achieved.

Figure 4b shows S R values as a function of d o for four orientations. Table 1 shows the corresponding raw data A second model was created by adding vertical fins to each overhang, as shown in Figure 5a. Table 2 shows the corresponding raw data. Figure 5b shows an even better agreement between S bim R and S re f R when vertical fins are added. Table 3 shows the Mean Absolute Error (MAE) for each orientation and protection configuration. The MAE is below five percent except for the overhang facing south. Merging all the data, MAE is equal to 0.038. The raw data can be found in Tables 1 and2. The discrepancy can be attributed to the modeling of the diffuse part. While the RTAADOM gives an explicit formula, our implementation requires an integration process sensitive to many parameters. Nonetheless, the results are very close and such small errors should be balanced by the automation gain, both in time and quality.

Conclusions and Perspectives

This paper introduces a new method to compute the exact shadow geometry necessary to evaluate shading factors. This method, called solid clipping, uses simple computational geometry tools implemented in BIM authoring tools and actionable through textual or visual API. Implementing the shading factor with a reference tool associated with French regulation has been validated. It has also been qualitatively shown that the design is unverifiable with the regulation and can be easily checked with our implementation. From an engineering perspective, two main benefits can be identified. Firstly, proper automation makes the cost of regulation checking plummet. Thus, checking a regulation becomes a quality assessment tool throughout the whole life of the project not just at the end. Secondly, implementing solid clipping using the BIM API allows us to deal with the geometry of shadows as complex as those of the building elements that generate them. In other words, if a solar protection design is created with a BIM authoring tool, the generated shadow can be leveraged to compute a shading factor. Combining these benefits will lead to a better cost-quality-delay balance in delivering sustainable buildings designed in BIM.

Extending Regulation Using BIM

Having assessed the precision and thus validated the implementation of shading factors, it is possible to extrapolate the computation to a configuration not handled by the regulation. So, S bim R was computed for protection with triangle-shaped fins instead of rectangular ones. Intuitively, this is an intermediate solution in terms of performance between the previously studied configurations: overhang with and without rectangular fins. Figure 6b confirms this by showing S bim R for the three studied protection configurations. As expected, the triangle-shaped fins generate intermediate values of S bim R . Beyond such a design as is currently implemented in the tropical French territory, the original and creative design can be assessed usingour method. For example, Figure 7 shows three examples of the solar protection of different architectural complexity for illustration. While inefficient, they highlight the capability of the solid clipping method to deal with complex geometries.

Research Directions

From a research perspective, this work paves the way in multiple directions. For example, advanced components such as phase-changing material or thermochromic glasses are exposed to non-uniform irradiance on their surface. Having easy access to realistic shadow geometries thanks to solid clipping and BIM modeling will improve models and thus their deployment for energy saving. More broadly, such exact and potentially complex shadows can be used as an input table in BEM tools to assess the impact of advanced solar protection. Smart solar protection could also use precise knowledge of shadow geometry. Indeed, having an internal reference in addition to sensors could help smart systems make better decisions and detect faults.

Figure 1 .

 1 Figure 1. Checking of a regulation based on shading factor computation is split into three stages: (1) shadow geometry computation, (2) shading factor computation, and (3) checking regulation.

Figure 2 .

 2 Figure 2. RTAADOM geometric model: (a) configuration and parameters to describe overhang, (b) configuration and parameters to describe fins, and (c) snapshot of the translated protected spreadsheet implementing the model with overhang in orange, vertical left fin in green and vertical right fin in yellow.

Figure 3 .

 3 Figure 3. Successive stages of computation for a simple 2D model: (a) The model of interest with a window and a few wall elements that can cast a shadow. The vector

Figure 3 .

 3 Successive stages of computation for a simple 2D model : (a) The model of interest with a window and few wall elements that can cast a shadow. The vector -→ sun is plotted. (b) Extrusion (c) Intersection with volumeIntersection in black (d) Face selection (red), extrusion of faces along -→ sun (dark yellow), intersection with F gives F shadow as a patch of surface (dark blue)

248Figure 3 .

 3 is plotted. (b) Extrusion. (c) Intersection with volumeIntersection in black. (d) Face selection (red), extrusion of faces along Version December 11, 2023 submitted to Buildings of rendering complex scenes in terms of geometry and daylighting. Indee compute shadow on one face F of the model and for a unique direct light Successive stages of computation for a simple 2D model : (a) The mode window and few wall elements that can cast a shadow. The vector -→ sun is plotted Intersection with volumeIntersection in black (d) Face selection (red), extrusion (dark yellow), intersection with F gives F shadow as a patch of surface (dark blue)

Algorithm 1 :then 6 PartialShadowVolume 7 8 List 11

 167811 Algorithm to compute shadow Input : A flat face F Output : A flat face F shadow 1 Light = Extrude(F, --→ sun, ∞) 2 CatchedLight = Intersect(Light, Building) 3 List = [ ] 4 foreach f ace in CatchedLight. f aces() do 5 if --→ n f ace • -→ sun ≤ 0 = Extrude( f ace, -→ sun, ∞) PartialShadowSur f ace = Intersect(PartialShadowVolume, f ace) ShadowSur f ace = Union(List) 12 return ShadowSur f ace Let us detail line by line Algorithm 1 along with a 2D illustration in Figure 3. The extension to 3D is straightforward.

  . The solid line represents values computed from the IFC file with solid clipping S bim R and the dashed lines correspond to the values given by the regulation Excel tool S re f R . As expected, S R decreases as d o increases. Except for the south orientation, there is an excellent agreement between S bim R and

Figure 4 .

 4 Figure 4. (a) Snapshot of a window with an overhang of depth d. (b) Shading factor S R as a function of the overhang depth d for the four cardinal orientations (color) and the two computation methods (line style).

Figure 5 .

 5 Figure 5. (a) Snapshot of a window with an overhang and rectangular vertical fins of depth d. (b) Shading factor S R as a function of the overhang and vertical fins depth d for the four cardinal orientations (color) and two computation methods (line style).

Figure 6 .

 6 Figure 6. (a) Snapshot of a window with an overhang and triangular vertical fins of depth d. (b) Shading factor S R as a function of the overhang and triangular vertical fins depth d for the north orientation is represented by the dotted curve (labeled ht). The overhang with (labeled hr) and without (labeled h) vertical rectangular fins is added for comparison.

Figure 7 .

 7 Figure 7. Snapshot of shadow computed by solid clipping for original solar protection design: (left) curved boundary parametrized by a spline, (right) rectangular shape with holes

Table 1 .

 1 Shading factor values computed with the spreadsheet (left) and with our approach (right) for solar protection with overhang.

			Reference					BIM
	d	North South	East	West	d	North South	East	West
	0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0	0.91 0.80 0.72 0.67 0.64 0.61 0.60 0.60 0.59 0.59	0.94 0.90 0.88 0.86 0.84 0.82 0.81 0.80 0.79 0.78	0.95 0.88 0.81 0.76 0.70 0.66 0.63 0.60 0.58 0.56	0.96 0.89 0.83 0.78 0.73 0.69 0.66 0.63 0.61 0.59	0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0	0.92 0.82 0.75 0.70 0.67 0.66 0.65 0.65 0.64 0.64	0.96 0.89 0.83 0.78 0.73 0.74 0.69 0.67 0.66 0.64	0.96 0.94 0.92 0.91 0.90 0.90 0.89 0.89 0.88 0.88	0.96 0.91 0.85 0.80 0.76 0.77 0.72 0.70 0.68 0.66

Table 2 .

 2 Shading

			Reference					BIM
	d	North South	East	West	d	North South	East	West
	0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0	0.66 0.54 0.47 0.42 0.38 0.37 0.35 0.34 0.34 0.33	0.71 0.64 0.59 0.56 0.52 0.50 0.48 0.47 0.45 0.44	0.80 0.70 0.63 0.56 0.51 0.47 0.44 0.41 0.38 0.37	0.81 0.73 0.66 0.60 0.55 0.51 0.47 0.44 0.42 0.40	0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0	0.70 0.56 0.46 0.40 0.36 0.33 0.32 0.31 0.30 0.29	0.84 0.74 0.65 0.58 0.53 0.48 0.45 0.42 0.39 0.37	0.79 0.71 0.64 0.59 0.56 0.53 0.50 0.48 0.47 0.45	0.87 0.78 0.70 0.64 0.58 0.53 0.49 0.46 0.43 0.40

factor values computed with the spreadsheet (left) and with our approach (right) for solar protection with overhang and vertical fins.

Table 3 .

 3 Mean absolute error between computed value of S

	re f R and S bim R .

Funding: This research was funded by the French Deposits and Consignments office grant number Fund: PIA3-MCENT

Data Availability Statement:

The data presented in this study are available in the article.

Conflicts of Interest:

The authors declare no conflict of interest.