Emilie Alberola

Julien Chevallier

Benoît Chèze

The EU ETS: CO 2 prices drivers during the learning experience (2005-2007)

Keywords: EU ETS, Cap-and-Trade Program, Climate Change Policy, CO 2 Price, Energy Markets, Weather Influences, Institutional Influences, Energy Policy

This chapter identifies the main price drivers of European Union Allowances (EUAs), valid for compliance under the European Union Emissions Trading Scheme (EU ETS) created in 2005 to regulate CO 2 emissions of more than 10,000 high carbon-intensive installations across Member States. Based on key design features of the EU ETS, this chapter develops carbon pricing strategies based on allowances supply and demand, institutional decisions, and the influence of other energy markets and weather conditions. Finally, we discuss the likely effects on economic growth on CO 2 emissions and carbon prices as a by product. The discussions developed in this chapter focus on Phase I (2005-2007) of the EU ETS, which may described as the "pilot" period for the future development of this environmental market scheme.

Introduction

The primary objective of the European Union Emissions Trading Scheme (EU ETS) is to regulate CO 2 emissions from the most energy-intensive industrial sectors. Installations have to surrender as many allocated allowances as their previous year's emissions. By allowing the emergence of a price -called a European Union Allowance (EUA) -for each emitted ton of CO 2 , the EU ETS aims at conveying appropriate price signals to industrial operators. They can select a combination of capital investments, operating practices, and emissions releases to minimise the sum of abatements costs and allowance expenses. The main characteristic of cap-and-trade systems consists in leaving to each firm the ability to decide its abatement strategies depending on its own risk preferences and perceptions about regulatory uncertainties, as well as technical opportunities (i.e. in terms of fuel-switching).

The EU ETS needed to be functional on January 1, 2005. Its creation required the deployment of new governance processes between the European Commission, national State authorities, and operators in the 10,000 (and over) installations covered by the Directive 2003/87/CE. The development of technical tools was also necessary to deal with the delivery of allowances on industrials' accounts, and to guarantee the monitoring of transactions. This fast-track operational implementation was carried out with some delays (e.g., the start-up of the Polish registry) but without major difficulties. In 2005, the European Union carbon finance sector was born. Within the first three years of its infancy, active players in the market multiplied, organised exchanges were established, and trading activity quickly developed. The volume of over-the-counter trades in emission allowances and on the marketplaces has been constantly growing since then.

This chapter aims at analyzing the price fundamentals of CO 2 allowances during 2005-2007. The remainder of the chapter is organized as follows. Section 2 reviews the key design issues of the EU ETS. Section 3 analyzes the CO 2 price drivers by detailing the emissions shortfall factor, supply and demand mechanisms, the influence of other energy markets and weather events, as well as the likely effects of economic growth. Section 4 concludes.

Key design issues on the EU ETS

This section briefly reviews the main design features of the European carbon market. The EU ETS has been created by the Directive 2003/87/CE. Across its 27 Member States, the EU ETS covers large plants from CO 2 -intensive emitting industrial sectors with a rated thermal input exceeding 20MWh. One allowance exchanged on the EU ETS corresponds to one ton of CO 2 released in the atmosphere. 2.2 billion allowances per year have been distributed during Phase I (2005Phase I ( -2007)). 2.08 billion allowances per year will be distributed during Phase II (2008)(2009)(2010)(2011)(2012). With a value of around 20€ per allowance, the launch of the EU ETS thus corresponds to a net creation of wealth of around 40 billion €. In January 2008, the European Commission extended the scope of the EU trading system to other sectors such as aviation and petro-chemicals by 2013, and confirmed its functioning for a third Phase until 2020. In the next section, we detail the price fundamentals of CO 2 allowances.

EUA Price drivers

In this section, we define first the emissions shortfall factor, second we detail the allowances supply and demand mechanisms, third we detail the influences of policy decisions, other energy markets and weather events. Fourth, we investigate the likely effects of economic growth on the price development of carbon prices.

The "emissions shortfall" factor

Whereas on energy markets the central pricing issue concerns the state of stocks, on the ETS the main question is the expected "emissions shortfall" factor during each compliance year. The emissions shortfall factor, defined as the difference between verified emissions and allocated allowances within a given compliance year, depends on the emissions abatements required by the cap, which are unknown but estimable based on reliable recent data.

Allowances supply and demand mechanisms

As for other commodity markets, the EUA price is driven by the balance between supply and demand, and by other factors related to market structure and institutional policies.

However, EUAs exhibit strong features of being a non-standard commodity: installations do not need to physically hold allowances to produce, but only to match them with verified emissions for their yearly compliance report to the European Commission.

Supply of EUAs

On the supply side, the first period of the market was quite limited, as no other carbon credits than Phase I EU allowances could be used for compliance: neither the use of Kyoto credits -even if transactions started in 2003 -nor Phase II EUAs -through borrowing -were allowed. Thus, the number of allowances available on the market was strictly equal to the number of allowances initially allocated to installations and to new entrants4 .

Demand of EUAs

While allowance supply is politically fixed by each Member States through National Allocation Plans (NAPs), allowance demand is function of expected CO 2 emissions. In the short and long run, the level of emissions depends on a large number of factors.

Short term abatement decisions, and thus the demand for allowances, are mainly driven by unexpected fluctuations in energy demand, energy prices and weather conditions (temperatures, rainfall and wind speed). CO 2 emissions in the EU-ETS are linked one-to-one to the use of fossil fuels (such as oil, gas, coal). In turn, the demand for fossil fuels depends on their absolute and relative prices. The marginal fuel switching costs from high carbonintensive sources of energy to low carbon-intensive sources for power and heat generation are the most important measures for carbon abatement in the short run, since power and heat operators are major actors on the EU ETS5 .

Other factors influencing the demand for allowances are unexpected fluctuations of demand for energy in the short run. As outlined by [START_REF] Alberola | Price Drivers and Structural Breaks in European Carbon Prices 2005-2007[END_REF], CO 2 prices are affected by unexpected weather variations like temperatures, rainfall and wind. Cold winters increase the need for heating by electricity or fuels, whereas warm summers lead to higher electricity demand for cooling and lower rate of utilization of nuclear power plants due to reduced cooling by rivers. Rainfall, wind speeds and sun shine hours affect the share of power generated by carbon-free heat generation from hydropower, wind and solar energy. As a whole, weather is widely acknowledged to have played an important role in explaining CO 2 prices.

In the long run, the demand of allowances can be affected by economic growth and marginal abatement costs. Abstract from very long term investments like carbon capture and storage (CCS), CO 2 emissions cannot be reduced by end-of -pipe technologies. Long-term marginal abatement costs are determined by investment decisions in low carbon-intensive energy utilities and in energy efficiency measures. Overall, with respect to high investment costs and uncertainties about the future of the EU ETS in 2005, firms faced mainly short term abatement costs during the first trading period when deciding on their allowance demand.

Institutional factors, other energy markets and weather influences

One of the most difficult factors to capture the influence on EUA demand is "market sentiment". This expression refers to factors such as uncertainty about future energy prices and policy decisions. These factors are indeed especially important for investors' expectations and the formation of their risk strategies.

There exists an emerging literature analysing the empirical relationships between EU carbon prices and its determinants during 2005-2007 that we propose to detail below.

Evidence on the impacts of the European Commission's voice

Evidence of the relationships between energy prices and carbon prices

Energy prices are the most important short term drivers of the demand for EUAs, due to the ability of power generators to switch between their fuel inputs [START_REF] Mansanet-Bataller | CO 2 Prices, Energy and Weather[END_REF], [START_REF] Delarue | Price determination of ETS allowances through the switching level of coal and gas in the power sector[END_REF], [START_REF] Alberola | Price Drivers and Structural Breaks in European Carbon Prices 2005-2007[END_REF], Bunn and Fezzi (2008), [START_REF] Ellerman | A Top-down and Bottom-up Look at Emissions Abatement in Germany in Response to the EU ETS[END_REF]). As shown by [START_REF] Alberola | Price Drivers and Structural Breaks in European Carbon Prices 2005-2007[END_REF], energy prices forecast errors have basically driven EUA price changes during 2005-2007, but their influence changed depending on institutional events. Brent, natural gas, coal, clean dark and spark spreads, as well as switch prices 8 all impact significantly carbon spot prices, but their influence vary following the two structural breaks detected in the time series, and the corresponding sub-periods under consideration.

As detailed above, weather conditions have an impact on EUA prices by influencing energy demand. During the winter 2006, colder temperatures than the decennial average have 7 Phase I spot prices decreased towards zero until the end of Phase I due to banking restrictions implemented between 2007 and 2008 (Alberola and Chevallier (2009)). 8 The dark spread is the theoretical profit that a coal-fired power plant makes from selling a unit of electricity having purchased the fuel required to produce that unit of electricity. The spark spread refers to the equivalent for natural gas-fired power plants. The equilibrium between these clean spreads represents the carbon price above which it becomes profitable for an electric power producer to switch from coal to natural gas, and below which it is beneficial to switch from natural gas to coal. As long as the market carbon price is below this switching price, coal plants are more profitable than gas plants -even after taking carbon costs into account. This switching price is most sensitive to changes in natural gas prices than to coal prices changes [START_REF] Kanen | Carbon Trading and Pricing[END_REF]. These three profitability indicators are used to determine the preferred fuel in power generation. had the expected statistical positive impact on EUA price changes. During the summer 2006 and the winter 2007, hotter temperatures than their decennial averages have affected negatively carbon price changes. [START_REF] Alberola | Price Drivers and Structural Breaks in European Carbon Prices 2005-2007[END_REF]).

The puzzle of economic activity effects on CO 2 emissions and EUA prices

Economic activity is maybe the most obvious and the less understood driver of CO 2 price changes. Economic growth leads to an increase in energy demand, and to higher industrial production in general. The first empirical analysis devoted to this topic may be related to Alberola et al. (2009a).

Conclusion

This chapter features an analysis of the CO 2 price determinants in the EU ETS during the "pilot" Phase (2005)(2006)(2007). Following a brief review of key design issues in the EU ETS, we have detailed the main carbon price drivers based on classic supply and demand analysis, but also with respect to the salient characteristics of the European cap-and-trade program (especially in terms of aggregated level of verified emissions and compliance).

Finally, we have emphasized that the main price drivers are linked to institutional decisions, other energy market prices (oil, gas, coal, electricity, and more generally fuelswitching in the power sector), and unanticipated temperatures changes. The last puzzle that remains to be solved concerning Phase I EUA price drivers lies in the determination of the influence of economic activity, through changes in CO 2 emissions levels, as initiated by recent empirical analyses.

  Institutional decisions on the overall stringency of the cap have an impact on the carbon price formation through initial allocation. Besides, any decision or announcement from regulators (being Member States or the European Commission) may induce changes in market players' behaviour. Official communications by the European Commission were essential to reach a better information flow on installations' net short/long positions 6 during Phase I. The disclosure of compliance results may be seen as the cornerstone of changes in market participants' expectation and strategies. In April 2006, when Member States and the European Commission disclosed the 2005 verified emissions for all EU ETS installations, the gap between initial allocation to industrials and business-as-usual emission forecasts appeared positive, thereby causing a sharp fall in carbon prices of all maturities (more than 50% in four days). Political choices can also impact the functioning of the market: in October 2006, announcements by the European Commission to validate more strictly NAPs during Phase II reinforced the general feeling that the market could be short during 2008-2012. From this date, Phase II prices increased and stabilized over 20€/ton of CO 2 . Alberola et al. (2008) identify statistically two structural endogenous breaks in the time-series of spot prices during 2005-2007: the first compliance period break goes from April 24, 2006 to June 12, 2006; and the second break due to the EU Commission stricter Phase II validation occurred from October 26, 2006 onwards 7 .

Figure 1 :

 1 Figure 1: Industrial production in EU 27 during 2005-2007

Figure

  Figure 2: Industrial productions in EU ETS sectors

Mission Climat Caisse des Dépôts, Paris. Emilie.alberola@caissedesdepots.fr

Corresponding author. Imperial College London (Grantham Institute for Climate Change), j.chevallier@imperial.ac.uk. South Kensington Campus, London SW7 2AZ, UK.

University Paris 10 (EconomiX-CNRS), benoit.cheze@u-paris10.fr

Note the equilibrium market price may also be affected by flexibility mechanisms, such as banking and borrowing provisions and the inflow of project mechanism credits, but these features were not applicable during Phase I.

The main sector whose share is the largest in terms of emissions covered by the EU ETS in Phase I is the power sector, which accounts for over 50% of emissions capped by the scheme. Their emissions abatement costs are assumed to be the lowest compared to others sectors, notably through fuel switching from coal to gas.

An installation is defined as short (long) when it records a deficit (surplus) of allowances allocated with respect to actual emissions[START_REF] Ellerman | Over-allocation or abatement? A preliminary analysis of the EU ETS based on the 2005-06 emissions data[END_REF]).

grew in all MS, from 2.24% in the UK to 12.95% in Italy. 9 The combustion sector represents the largest share of installations in the EU ETS with 70\% of the EU allocation. The combustion sector was defined in a different way by each MS and contains too many subactivities. Based on the CITL data and the Classification NACE Rev.1 C-F, [START_REF] Trotignon | Technical Memorandum on Analysis of the EU ETS Using the Community Independent Transaction Log[END_REF] and [START_REF] Trotignon | European CO 2 Market and the CITL: The Trial Period Under Scrutiny[END_REF] classify between large electricity production plants, district heating facilities (cogeneration when details were available) and other installations.


Nolann Lainé
email: nolann.laine@creatis.insa-lyon.fr

Herve Liebgott
email: herve.liebgott@creatis.insa-lyon.fr

Guillaume Zahnd
email: guillaume.zahnd@helmholtz-muenchen.de

Maciej Orkisz
email: maciej.orkisz@creatis.insa-lyon.fr

Hervé Liebgott

Carotid artery wall segmentation in ultrasound image sequences using a deep convolutional neural network

Keywords: segmentation, ultrasound images, deep learning

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

Introduction

Atherosclerosis, the "silent killer" that can progress to acute events such as stroke or heart attack, begins to develop in the intima-media complex (IMC), the innermost structure of the artery wall. IMC thickening is widely considered as a marker of the atherosclerosis onset. Its screening is usually performed in ultrasound (US) images of the carotid artery. In the clinical routine, the measurement of the intima-media thickness (IMT) is performed manually, which may lead to substantial inter-and intra-observer variability. Numerous methods, compared in several surveys [START_REF] Molinari | A state of the art review on intima-media thickness (IMT) measurement and wall segmentation techniques for carotid ultrasound[END_REF][START_REF] Loizou | A review of ultrasound common carotid artery image and video segmentation techniques[END_REF][START_REF] Meiburger | Automated localization and segmentation techniques for B-mode ultrasound images: A review[END_REF], have been developed to automate the measurement of the carotid artery IMT. Most of them infer the IMT as the distance between the lumen-intima (LI) and media-adventitia (MA) contours extracted within a region of interest (ROI) cropped manually or automatically by seeking the window where the characteristic dual-line pattern is discernible (Fig. 1). Contour extraction is performed using various "conventional" algorithms, such as snakes or dynamic programming. Virtually all assume that the contours are smooth curves passing through points with strong intensity gradients expected between artery lumen and intima, and between media and adventitia layers. Lately, deep-learning-based (DL) methods have also been proposed for this task and obtained promising results, but the comparisons [START_REF] Biswas | A review on joint carotid intima-media thickness and plaque area measurement in ultrasound for cardiovascular/stroke risk monitoring: Artificial intelligence framework[END_REF] were performed on different datasets. A recent study [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): An open multicenter analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF] compared five state-ofthe-art conventional methods on an open access dataset CUBS1 [START_REF] Meiburger | DATASET for "Carotid Ultrasound Boundary Study (CUBS): an open multi-center analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF] containing 2176 US images from two centers, and concluded that their accuracy is comparable to skilled human experts, while presenting substantially smaller variability.

In that study, the best results were achieved by a method based on dynamic programming [START_REF] Zahnd | Simultaneous extraction of carotid artery intima-media interfaces in ultrasound images: assessment of wall thickness temporal variation during the cardiac cycle[END_REF]. A subsequent study [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): Technical considerations on an open multi-center analysis of computerized measurement systems for intima-media thickness measurement on common carotid artery longitudinal B-mode ultrasound scans[END_REF] confronted the same methods with two DL approaches on another dataset CUBS2 [START_REF] Meiburger | DATASET for "Carotid Ultrasound Boundary Study (CUBS): Technical considerations on an open multicenter analysis of computerized measurement systems for intima-media thickness measurement on common carotid artery longitudinal B-mode ultrasound scans[END_REF] containing 500 US images from five different centers. This comparison has demonstrated the sensitivity of several methods to the arterial wall morphology, namely a decreased accuracy in curved or inclined arteries. One of the DL-based methods, previously unpublished, outperformed the others in almost all comparisons except for processing time. While the selected network architecture was the well-known U-net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], the original contribution of the method resides in a specific patch-based strategy devised to cope with inclined and curved arteries. In this article, we describe for the first time the details of this method, and then report main results in terms of accuracy and robustness on both open-access datasets CUBS1 [START_REF] Meiburger | DATASET for "Carotid Ultrasound Boundary Study (CUBS): an open multi-center analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF] and CUBS2 [START_REF] Meiburger | DATASET for "Carotid Ultrasound Boundary Study (CUBS): Technical considerations on an open multicenter analysis of computerized measurement systems for intima-media thickness measurement on common carotid artery longitudinal B-mode ultrasound scans[END_REF], compared to available expert annotations and the other methods. We also assess how much the processing time can be improved by changing the parameter setting, without degrading the accuracy.

Datasets

In this work, for training and evaluation purposes, we have used both previously mentioned publicly-available datasets CUBS13 and CUBS24 , i.e., a total of 2676 static US images acquired with different clinical equipment and including a subset (n = 100) of simulated images. These images have pixel sizes ranging from 29 to 267 µm (mean size 64 µm in CUBS1 and 60 µm in CUBS2) and come along with LI and MA contours independently traced by two or more experts, one of them having traced them twice. The experts independently selected a ROI, where they considered the LI and MA interfaces as sufficiently perceptible, and traced control points within it. To obtain smooth contours, piecewise cubic Hermite interpolating polynomial (PCHIP ) was applied using MATLAB, Version 2020b (The Math Works, Inc.). The same annotations as in the previous studies [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): An open multicenter analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF][START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): Technical considerations on an open multi-center analysis of computerized measurement systems for intima-media thickness measurement on common carotid artery longitudinal B-mode ultrasound scans[END_REF] were used as reference. These will be referred to as A1. According to A1, the mean IMT was of 725 µm in CUBS1 and 857 µm in CUBS2. The same expert's second set of annotations (A1') was used to calculate the intra-observer variability, while the inter-observer variability was calculated between A1 and the second expert's annotations (A2).

Segmentation method

Similarly to many existing methods, e.g. [START_REF] Zahnd | Simultaneous extraction of carotid artery intima-media interfaces in ultrasound images: assessment of wall thickness temporal variation during the cardiac cycle[END_REF], the first step is a very simple user interaction (two mouse clicks) defining the exploitable ROI, where the IMC is perceptible. The remainder is fully automatic and split into two main steps: localization of the far wall and actual IMC segmentation (Fig. 2). The proposed solution is based on U-net architecture [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF], which has widely demonstrated its ability to produce accurate results in medical image segmentation with limited annotated data available for supervised training. We have kept standard components, such as ReLU activation function, and experimentally selected the following settings: number of layers 5, initial number of convolution kernels 32, and size of the kernel 3 × 3 (with batch normalization and bias). We implemented dilated convolutions on the bottleneck, to increase the receptive field [START_REF] Meshram | Deep learning for carotid plaque segmentation using a dilated U-net architecture[END_REF]. The outputs of the filters with different dilation factors {1, 2, 3, 4, 5, 6} are concatenated into one single tensor. The U-net operates in fixed-size (128pixel width W , 512-pixel height H) overlapping patches distributed within the ROI. The two first layers contract only the height, so as to achieve a square shape (128 × 128) of the feature space starting from the third layer. The horizontal overlap between patches is equal to W -∆x, where ∆x is the horizontal shift. A post-processing combines the predictions made within the patches to extract smooth contours over the entire ROI regardless of its width. The core of the method consists of two steps: approximately localizing the far wall (Fig. 2c), and precisely segmenting the IMC around this location (Fig. 2d).

Detection of the far wall

Like in many state-of-the-art methods [START_REF] Zahnd | A fully-automatic method to segment the carotid artery layers in ultrasound imaging: Application to quantify the compression-decompression pattern of the intima-media complex during the cardiac cycle[END_REF][START_REF] Wang | Fully automatic measurement of intimamedia thickness in ultrasound images of the common carotid artery based on improved Otsu's method and adaptive wind driven optimization[END_REF][START_REF] Menchón-Lara | Early-stage atherosclerosis detection using deep learning over carotid ultrasound images[END_REF][START_REF] Qian | Segmentation of the common carotid intima-media complex in ultrasound images using 2-D continuous max-flow and stacked sparse autoencoder[END_REF], the initial localization of the far wall is performed to focus the subsequent actual segmentation. In this step, the algorithm attempts to separate the ROI in two regions: above and below the median axis of the IMC, respectively. Here, the patches are of full image height, and the corresponding U-net will be referred to as Θ F W . We first describe how the images were pre-processed and how Θ F W was trained, then we specify the post-processing chosen to obtain the curve approximately localizing the far wall on the entire ROI width from patch-wise predictions inferred using Θ F W .

Pre-processing and training: All images of the database were resampled to a constant height of 512 pixels; as the native height of all images in the database is around 600 pixels, the distortion thus introduced was minimal (regions containing alphanumerical information and/or ECG curve overlaid onto the actual image were not clipped). For training data, the median axis of the IMC was defined as the line halfway between LI and MA annotations, interpolated across the entire width of the ROI, and a reference mask (M ROI ) was generated by setting all pixels below the median axis to one and the others to zero. Then the ROI and M ROI were identically cut into patches; a 100-pixel overlap (∆x = 28) between patches aimed at data augmentation. Thus obtained patches with their associated masks (Fig. 3) were fed into the training process, which used the ADAM optimizer and the sum of the binary cross-entropy and the Dice loss as the loss function. The latter was experimentally chosen to minimize the boundary distance and maximize the overlay with respect to the reference masks. Training was performed with on-the-fly data augmentation: horizontal and vertical flip applied on 50% of the images, and affine transformation with random rotation [-2, 2] degrees, shearing [-2, 2] degrees, and translation (vertical [-20, 20] pixels, horizontal [-5, 5] pixels). The batch size was of 32, and the number of epochs 300 with early stopping switched on, so that the process stopped if the loss value on the validation subset did not improve during 50 epochs.

Inference and post-processing: Prior to inference, each image is resampled as described above, and then the corresponding ROI is cut into overlapping 128 × 512-pixel patches. Next, all patches are segmented using Θ F W . Knowing the location and the size of each patch, two maps are created:

-prediction map: contains, for each pixel, the sum of values predicted by Θ F W in all patches. -overlay map: contains, for each pixel, the number of overlapping patches it belonged to.

Dividing the prediction map by the overlay map provides, for each pixel, an average value in the range [0, 1], which is then binarized by using a threshold of 0.5, to obtain the segmentation map. The latter is cleaned by retaining the largest connected component. The median axis we seek is the upper boundary of thus segmented region. Eventually, this boundary is smoothed by using a third-order polynomial regression.

Segmentation of the IMC

The above-described approximation of the far-wall median axis is used to initialize the actual segmentation of the IMC. This segmentation process presents many similarities with concepts explained in Section 3.1: overlapping patches of 128 × 512 pixels, an overlay map, a prediction map, a similar post-processing except that two contours are extracted (the LI and MA interfaces), as well as the same optimizer, loss function, and U-net architecture. The dilated U-net trained here will be referred to as Θ IM C . Hereafter, we emphasize the specific choices made for this step.

Pre-processing and training:

The segmentation task has to be as accurate as possible, hence the algorithm works at a sub-pixel resolution. To this purpose, the vertical pixel size of the images was homogenized to 5 µm using a linear interpolation. According to this physical size, the patch height H = 512 pixels roughly corresponds to 2.6 mm, which aims to encompass the IMC, knowing that the average IMC thickness is about 0.8 mm. For training, the ground truth was then deduced from thus interpolated images (Fig. 4): each pixel located between the annotated LI and MA interfaces was set to one, and the others to zero. Unlike the far wall detection, the patches were extracted along the median axis and, in addition to the 100-pixel horizontal overlap (i.e., ∆x = W -100 = 28-pixel shift), at each abscissa x i , the mean ordinate y i of the median axis was computed on the patch width W = 128, and three patches were extracted, respectively centered at y i and y i ± ∆y (∆y = 128). This data augmentation attempted to cope with possibly inaccurate far-wall approximation as well as with tilted arteries.

Inference and post-processing: During inference, the patches are extracted along the far wall approximation resulting from the first step (Section 3.1). At each abscissa x i three or more patches are captured at different ordinates, depending on the tilt of the median axis, with the goal to cover all the expected extent of the IMC. The predictions made by Θ IM C in all patches are combined into a prediction map, and then the segmentation map is derived thereof, as described above. Finally, the LI and MA interfaces are respectively defined as the upper and lower boundaries of the region thus segmented.

Results

The evaluation was carried out using 5-fold cross-validation, so as to assess each network on data not seen during its training. In each fold, the database (combined data from the open-access datasets CUBS1 [START_REF] Meiburger | DATASET for "Carotid Ultrasound Boundary Study (CUBS): an open multi-center analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF] and CUBS2 [START_REF] Meiburger | DATASET for "Carotid Ultrasound Boundary Study (CUBS): Technical considerations on an open multicenter analysis of computerized measurement systems for intima-media thickness measurement on common carotid artery longitudinal B-mode ultrasound scans[END_REF]) was split into training (60%), validation (20%), and testing (20%) subsets. Thus, five pairs of networks Θ F W and Θ IM C were trained and tested independently, and the results reported here are the merging of the test sets of these five pairs, thus evaluating the method on the entire database.

In the proposed cascade approach, a failure of the first step (far wall detection) will trigger a failure of the second step (IMC segmentation). To conduct Table 1. Errors of segmentation results (mean ± standard deviation): mean absolute difference (MAD) for thickness quantification (IMT) and Hausdorff distance (HD) for contour locations (LI and MA). The comparisons are reported for the entire database and separately for CUBS1 and CUBS2 datasets. For each dataset, the following differences are reported: the proposed method against reference annotations, the inter-and intra-observer variabilities. a fair evaluation of both steps, we first quantified the success rate of the first step alone, and then we quantified the accuracy of the second step by manually enforcing valid initial conditions when needed.

Robustness of the far wall detection: After visual inspection, 36 out of 2676 predicted median axes (1.3% of the database) were considered as failures, i.e., curves unusable to initialize the IMC segmentation step. Hence, the success rate was of 98.7% and, in the 36 images with failures, the median axis was manually redrawn using a home-made graphical interface.

Accuracy of the IMC segmentation:

The segmentation error was quantified in two ways. To enable a fair comparison of our method with five state-of-theart methods evaluated on the CUBS1 and CUBS2 datasets in the previously mentioned recent studies [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): An open multicenter analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF][START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): Technical considerations on an open multi-center analysis of computerized measurement systems for intima-media thickness measurement on common carotid artery longitudinal B-mode ultrasound scans[END_REF], we asked the first authors of the latter (see Acknowledgements) to apply onto the LI and MA contours extracted by our method exactly the same evaluation protocol and metrics as used in their publications. Hence, the IMT values and all errors were calculated on a region restricted to a common support where all participating methods succeeded in extracting the contours. The IMT was measured as the polyline distance between LI and MA contours, and its error was calculated as the mean absolute difference (MAD) between the method output and the reference values calculated from the expert's (A1) annotations. The (worst case) errors for LI and MA were separately estimated by calculating the Hausdorff distances between the respective contours extracted by the method and the annotations performed by A1. These results are summarized in Table 1.

The reported values were obtained with the same implementation and parameter settings as in the previous study restricted to the CUBS2 dataset [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): Technical considerations on an open multi-center analysis of computerized measurement systems for intima-media thickness measurement on common carotid artery longitudinal B-mode ultrasound scans[END_REF]. In that study, the processing time was the only drawback of our method: 0.92 sec for far-wall localization and 1.77 sec for IMC segmentation, i.e., a total of 2.69 sec per image, and was achieved with a very large horizontal overlap of 124 pixels (i.e., shift ∆x = 4). Therefore, without modifying the implementation, we varied ∆x ∈ {8, 16, 32, 64, 96} to check how much the horizontal overlap -and thereby the number of processed patches -can be decreased without degrading the results, and how much the processing time can thus be improved. The results obtained on the entire database are summarized in Table 2. Here the IMT and the errors in contour extraction were assessed using in-house tools, by measuring the column-wise distances instead of the polyline and Hausdorff distances. Also, the comparison was not restricted to the common support for all methods, and used the full extent of the extracted contours. In addition, the failures of the automatic far-wall localization step were not manually corrected. For all these reasons, the results for ∆x = 4 differ from those in Table 1, but they can be consistently compared with those obtained for the remaining shifts.

Discussion

We developed and assessed a deep-learning method to extract the contours of the intima-media complex in longitudinal B-mode ultrasound images of the carotid artery. The method first approximately localizes the far wall, and then segments the anatomical interfaces of interest. The proposed approach allows for variable width of the region of interest without rescaling the images. Robustness of the far-wall localization step is a prerequisite for overall correct segmentation. This step was successful in 98.7% of the images, therefore attesting of the method robustness. The actual segmentation step achieved good accuracy, with errors smaller than the inter-observer and intra-observer variability. The results reported in Section 4 were obtained on combined data from the open-access datasets CUBS1 [START_REF] Meiburger | DATASET for "Carotid Ultrasound Boundary Study (CUBS): an open multi-center analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF] and CUBS2 [START_REF] Meiburger | DATASET for "Carotid Ultrasound Boundary Study (CUBS): Technical considerations on an open multicenter analysis of computerized measurement systems for intima-media thickness measurement on common carotid artery longitudinal B-mode ultrasound scans[END_REF] and can be compared with the results of five state-of-the-art methods evaluated on the same datasets in recent studies [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): An open multicenter analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF][START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): Technical considerations on an open multi-center analysis of computerized measurement systems for intima-media thickness measurement on common carotid artery longitudinal B-mode ultrasound scans[END_REF]. On both datasets, our method outperformed the remaining ones. In particular, on the CUBS1 dataset, the compared methods produced IMT errors (MAD) ranging from 114 ± 117 to 255 ± 230 µm, while for our method the MAD was of 99 ± 89 µm. Similarly, on the CUBS2 dataset, the other methods produced IMT errors ranging from 139 ± 118 to 224 ± 178 µm, while for our method the MAD was of 106 ± 89 µm.

On average, the MAD produced by the proposed method represented 1.6 pixels (which corresponds to 13.7% of the mean IMT) in CUBS1, and 1.8 pixels (12.4%) in CUBS2, respectively. Although these inaccuracies remain relatively large compared to the target measurement (IMT), they are respectively twice and 1.5 times smaller than the inter-observer and intra-observer variabilities.

The processing time was below 1 sec for the far-wall localization step and below 2 sec for the actual IMC segmentation step, when using a very large horizontal overlap and thereby a large number of patches. It can be seen that, on average, the errors were stable when decreasing the overlap, particularly up to the shift ∆x = 16, but the resulting reduction of the processing time was limited. Of note, when segmenting a sequence of images, the far-wall localization step needs to be executed only once, at the beginning of the sequence, while the actual IMC segmentation needs to be repeated for each image. In the current implementation, the processing time of this step remains greater than 0.5 sec even with a strongly reduced overlap. Further reduction of this time would require an optimization of the implementation, which has not yet been attempted.

The largest errors occurred in the presence of atherosclerotic plaques. As the work presented here was oriented towards asymptomatic plaque-free subjects, images with plaques were not expected. Nevertheless, we anticipate that results might be improved by increasing the number of such images in the database and re-training the networks.

In conclusion, with a 98.7% success rate and errors smaller than the interand intra-observer variability, the proposed method likely is robust and accurate enough to be used in clinical practice, as well as to study the periodic compression-decompression of the arterial wall, which is another potential biomarker of atherosclerosis.

Code and examples are available at https://github.com/nl3769/caroSegDeep.

Compliance with ethical standards information

The data from human subjects used in this work were obtained and treated in line with the principles of the Declaration of Helsinki. Approval was granted by the Ethics Committees of the institutions involved in creating the multicentric database, from which these data were accessed.

Fig. 1 .

 1 Fig. 1. Longitudinal US images of the common carotid artery, examples. The curves represent the lumen-intima (red) and media-adventitia (green) boundaries of the far wall interpolated from expert's annotations (dots). The annotations were restricted to the exploitable region (ROI) where the characteristic double-line pattern of the IMC is discernible. Yellow arrows point examples of inconsistent annotations.

Fig. 2 .

 2 Fig. 2. Flowchart of the proposed method: a) Input image. b) User delimitation of the left and right borders of the ROI. c) Far wall detection: patches are extracted through a sliding window with overlap within ROI borders; post-processing of predicted overlapping masks leads to extraction of the median axis (green). d) IMC segmentation: overlapping patches are picked along the median axis and post-processing of the predicted masks leads to extraction of the LI (red) and MA (yellow) interfaces.

Fig. 3 .

 3 Fig. 3. Data used during the far wall detection training phase: patch cut from the ROI (left) and the corresponding patch from MROI . Red and cyan dots represent the annotations for LI and MA interfaces, respectively. The green curve is the median axis calculated from the interpolated annotations.

Fig. 4 .

 4 Fig. 4. Data used during the training phase for IMC segmentation. Patches and their associated masks located at: a) (xi, yi -128), b) (xi, yi), and c) (xi, yi + 128). Red and blue dots represent the corresponding annotations for LI and MA interfaces, respectively. The contours were obtained by interpolating the annotations. Of note: the actual size of the patches is 128 × 512 pixels, but they have been vertically compressed for display purposes.

Table 2 .

 2 Estimation errors (mean of absolute differences ± their standard deviation) with respect to reference annotations A1, and processing time (FW stands for farwall localization and IMC for the actual segmentation), as a function of the horizontal shift ∆x. The horizontal overlap between patches is equal to 128 -∆x. The reported processing time was achieved on the following hardware: CPU Intel(R) Core(TM) i7-6700, 32 GB RAM,3.40 GHz; GPU NVIDIA GeForce GTX 1070, 8 GB RAM. std (µm) 158 ± 88 158 ± 88 158 ± 88 159 ± 88 159 ± 89 160 ± 87 LI MAD ± std (µm) 116 ± 68 116 ± 68 116 ± 68 113 ± 65 114 ± 66 118 ± 71 MA MAD ± std (µm) 104 ± 62 104 ± 62 104 ± 63 107 ± 68 107 ± 68 105 ± 62

	Shift ∆x (pixels)	4	8	16	32	64	96
	IMT MAD ± FW time (sec)	0.92	0.73	0.64	0.59	0.58	0.56
	IMC time (sec)	1.77	1.14	0.80	0.65	0.58	0.55
	total time (sec)	2.69	1.87	1.44	1.24	1.16	1.11

https://data.mendeley.com/datasets/fpv535fss7/1

https://data.mendeley.com/datasets/m7ndn58sv6/1

Acknowledgments

The authors thank K. M. Meiburger and F. Marzola for their help in calculating the comparisons with reference annotations on CUBS1 and CUBS2 datasets. This work was partly supported, via NL's doctoral grant, by the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).


Nolann Lainé

Guillaume Zahnd

Olivier Bernard

Maciej Orkisz

Herve Liebgott

Gan

Th Maciej Orkisz

Hervé Liebgott

Generation of realistic simulated B-mode image texture with a GAN

Keywords: Ultrasound simulation, carotid artery, deep learning, generative adversarial network, domain adaptation

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

I. INTRODUCTION

Cardiovascular diseases-particularly atherosclerosis, known as the "silent killer"-are considered the leading cause of death worldwide [START_REF] Kaptoge | World Health Organization cardiovascular disease risk charts: revised models to estimate risk in 21 global regions[END_REF]. We seek to detect early manifestations of the disease by extracting biomarkers from the common carotid artery (CCA), considered the sentinel of atherosclerosis [START_REF] Yousefi Rizi | Carotid wall longitudinal motion in ultrasound imaging: An expert consensus review[END_REF].

In particular, we aim to quantify the displacement field of the CCA far wall in ultrasound (US) image sequences and to deduce thereof its kinematic behavior through a cardiac cycle, in order to predict if a patient is ill. Nevertheless, obtaining reliable references to validate motion tracking algorithms remains difficult, as ground truth information is unavailable. Moreover, newly developed algorithms based on fully-supervised deep learning generally perform better than traditional methods in other applications [START_REF] Jiang | Learning to estimate hidden motions with global motion aggregation[END_REF], [START_REF] Teed | Raft: Recurrent all-pairs field transforms for optical flow[END_REF], but the paucity of said ground-truth is an obstacle to the generation of an adequate training set. Therefore, we propose a simulation pipeline that aims to produce realistic in silico sequences, along with associated deformation parameters that can serve as fully-determined references.

Most methods to simulate US images of arteries [START_REF] Makūnaitė | Simulation of ultrasound rf signals backscattered from a 3D model of pulsating artery surrounded by tissue[END_REF], [START_REF] Meiburger | Validation of the carotid intima-media thickness variability: Can manual segmentations be trusted as ground truth?[END_REF] consist of physics-based ultrasound simulators such as Field II [START_REF] Jensen | Field: A program for simulating ultrasound systems[END_REF], [START_REF] Jensen | Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers[END_REF] or SIMUS [START_REF] Garcia | Make the most of MUST, an open-source Matlab UltraSound Toolbox[END_REF]. These methods take as input in silico phantoms, which are typically generated using specific statistical distributions of reflection coefficients. Alternatively, in silico phantoms have also been simulated from in vivo images [START_REF] Alessandrini | A pipeline for the generation of realistic 3D synthetic echocardiographic sequences: Methodology and open-access database[END_REF], [START_REF] Sun | A pipeline for the generation of synthetic cardiac color Doppler[END_REF]. The latter solution is attractive to generate realistic samples because it directly uses a genuine in vivo image as the source, but its full potential remains locked by the incomplete knowledge of the modeled tissues and of the image-acquisition system parameters. In this work, we present a simulation pipeline devised to counteract these limitations via the inclusion of a post-processing operation based on a generative adversarial network (GAN, [START_REF] Mirza | Conditional generative adversarial nets[END_REF]).

II. INITIAL SIMULATION PHASE

A classical simulation pipeline takes an in vivo clinical static image of the CCA as input and the output is a temporal sequence of synthetic images. We first construct an in silico phantom (part II-A), which is fed into a realistic US simulator (in our case, Field II) along with probe and acquisition parameters (part II-B). These steps are detailed below.

A. Phantom construction

In order to calculate the amplitudes of the reflection coefficients of the scatteres in the in silico phantom from an in vivo B-mode image I, we first inverse the log-compression and introduce variability through Gaussian distribution (see [START_REF] Sun | A pipeline for the generation of synthetic cardiac color Doppler[END_REF]). Reflection coefficients RC m , m = 1 . . . M thus calculated are then assigned to M points at locations (x m , y m ) drawn normally at random:

RC m = I (x m , y m ) max(I) 1 γ • N (0, 1) . ( 1
)
Here γ is the gamma-compression constant (set to 0.3), and N (0, 1) is the normal distribution. The scatterer map is constructed with a density of 10 per square wavelength, so M 300000 for one image in our case.

B. Modeling of the system

When using clinical images, it is difficult to know all the image-acquisition and reconstruction parameters, particularly for images from open-access multi-center databases such as those used in our work: CUBS 1 [10] 1 and CUBS 2 [9] 2 . Therefore, we had to design our probe model and specify empirical simulation parameters. We have selected a linear probe L12-3v from MUST [START_REF] Garcia | Make the most of MUST, an open-source Matlab UltraSound Toolbox[END_REF], consistent with our application, on which we just changed the number of elements from 128 to 192 (Fig. 2). The acquisition parameters are summarized in Table I. The in vivo images were acquired with scanlinebased beamforming. The sonographer therefore chose a depth of focus located around the IMC of the far wall. For each acquisition, the depth varies (in meters) and we should adjust it for each simulation. To ensure the focus at the IMC, we decided to use a synthetic aperture, thus all depths are at focus in both transmission and reception. This choice may lead to a slight lack of realism away from the IMC.

C. Considerations about the current modeling limitations

We noticed four main limitations with the conventionally-used modeling strategy:

1) The compression parameters used to produce the envelope image are unknown. Applying an inverse gamma compression to return to the gray level distribution of the envelope image (Eq. 1) is likely inaccurate. Consider the upper part as the in vivo system. For ease of explanation, the medium is represented as a single point scatterer. The resulting B-mode image is therefore the PSF. The simulation (bottom) does not take a single point scatterer as above (because we don't know where it really is and what its reflection coefficient is), but a cloud of scatterers with coefficients calculated from the diffuse spot representing the PSF in the B-mode image, which amounts to adding the double convolution of the in vivo system to the double convolution simulated.

2) The inverse reconstruction operator used to obtain the envelope image (namely, the delay-and-sum algorithm) only achieves imperfect fidelity when recovering the initial medium. Therefore, the resulting gray level distribution of the pixels is expected to not exactly correspond quantitatively to the underlying distribution of reflection coefficients.

3) The actual spatial coordinates and the spatial density of the scatterers in the in vivo medium are unknown. Populating the scatterer distribution via drawing at random from the normal distribution from an empirical range is an approximation of the model. 4) The in vivo images are acquired by physical systems that systematically introduce a "double spread" effect: the point spread function (PSF) of the system is defined by the transmission and the reception properties of the system (P SF := pulse * IR trans * IR trans ). Therefore, the reflection coefficients extracted from in vivo images do not correspond to true values. When subsequently applying an ultrasound simulator onto such scatterer distribution, the PSF of the simulator actuates, for a second time, the response of the physical system, which results in exaggeratedly smoothing the simulated signals (Fig. 3). Note that we also do not know with precision the exact excitation signal and impulse response of the transducer used to acquire in vivo images.

The appropriate analytical solution to construct in silico phantoms would be to apply a mathematically-rigorous modelbased inverse reconstruction operator, or at minima, to utilize a proper deconvolution algorithm to extract the scatterers distribution from the in vivo image, however these approaches are still exploratory in the field of ultrasound. An alternative solution may be the use of an artificial intelligence technique named "domain adaptation". The key idea is to train a neural network for in silico to in vivo image translation. To this end, we suggest to use a Generative Adversarial Network (GAN) [START_REF] Mirza | Conditional generative adversarial nets[END_REF], which is well suited for domain adaptation.

Recently, GANs were used to generate synthetic ultrasound images from masks [START_REF] Cronin | Using deep learning to generate synthetic B-mode musculoskeletal ultrasound images[END_REF], or to mimic the image texture with a physics-based simulator [START_REF] Tom | Simulating patho-realistic ultrasound images using deep generative networks with adversarial learning[END_REF]. Nevertheless, to the best of our knowledge, no GAN-based approach was proposed to generate ultrasound image sequences.

III. GAN-BASED POST PROCESSING PHASE A. Dataset

We used a subset (n = 769) of in vivo images from the CUBS 1 [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): An open multicenter analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF] and CUBS 2 [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): Technical considerations on an open multi-center analysis of computerized measurement systems for intima-media thickness measurement[END_REF] databases. As the images also contain text and other graphical inlays, we clipped them out and cropped each sample down to only a region of interest (ROI). Then, the scatterer maps were constructed for each image and simulated using the pipeline presented in Figure ,1, bypassing step 5 (post-processing). Thus, the collection of in vivo images (R) is matched with a corresponding collection of simulated in silico images (S). It is noteworthy that no postprocessing is applied during the simulation. The GAN takes as input S and learns how to predict R. These pairs of images (S, R) made up the data set, which was split into training (70 %), validation (20 %), and test (10 %) subsets.

B. Network & training

A GAN consists of two different networks during the training process: a generator and a discriminator. The generator performs the image-domain adaptation while the discriminator helps the generator by assessing if the generator output resembles an image from the target domain (here in-vivo image). Learning therefore relies on two different loss functions, one for the discriminator (Eq. 2) and one for the generator (Eq. 3).

L Discr = - 1 2 (log (1 -ŷ R) + log (ŷ R )) , (2)
L Gen = -λ discr log (ŷ R) + λ pixel ||R -R|| 2 , (3)
where R is the in vivo image, R is the GAN output, ŷR and ŷ R are the probabilities returned by the discriminator that R and R respectively are real images, while λ discr and λ pixel are weighting factors set to 10 -3 and 1, respectively.

During training, we used a batch size of 6. Concerning the pre-processing, the ROIs were normalized between [0, 255] and resized to 256 × 512 pixels. We evaluated three different architectures for the generator (namely, U-net, dilated U-net, and ResNet-based super-resolution GAN).

IV. RESULTS

Peak signal-to-noise ratio (PSNR, Eq. 4) was used to measure the similarity between images I 1 and I 2 (S vs. R or S vs. R):

P SN R(I 1 , I 2 ) (dB) = 10 • log 10 255 2 ||I 1 -I 2 || 2 , (4)
where 255 is the maximal gray level of the pixels. Quantitative results are reported in Table II, while Figure 4 qualitatively illustrates example results obtained with the best architecture. V. DISCUSSION

The networks, and particularly the dilated U-net, help to provide in silico images that are more similar to the source clinical images compared to the intermediate simulation result (Fig. 4). In the initial simulated images, we observed that the proximal wall often lacks luminance, which likely means that the simulation fails to reproduce the presence-or lack thereof-of time gain compensation. Conversely, the GAN succeeds in correcting this axial profile, which explains the gain in dB. Nevertheless, current GAN results do not seem to improve the image texture, as the GAN output appears smoother than both simulated and in vivo images. Moreover, when applying the GAN to an in silico sequence, in which the scatterers obeyed a controlled displacement field, we observed a jitter in the movie appearance.

There are several avenues to explore to improve our method. First, the proposed GAN is not scanner-specific and has learned to average the features from different scanners. Conditional GAN architectures could be a solution to take scanner information into account. Second, temporal information was not explored. A cycle-GAN architecture fed with a 2D + t volume could learn the spatio-temporal filtering and temporal coherence and likely eliminate jitter effects. Last but not least, we noticed a potential issue in the pre-processing step:

The ROIs extracted from the in vivo images have different initial sizes, but are reshaped to a fixed size. Therefore, the appearance (speckle size and radial-to-axial proportions) of the texture is altered, and the network is forced to operate under

Figure 1 .

 1 Figure 1. Simulation pipeline with a GAN-based post-processing (step 5).

Figure 2 .

 2 Figure 2. Linear probe used in our simulation.

Figure 3 .

 3 Figure3. Expansion of the PSF function. Consider the upper part as the in vivo system. For ease of explanation, the medium is represented as a single point scatterer. The resulting B-mode image is therefore the PSF. The simulation (bottom) does not take a single point scatterer as above (because we don't know where it really is and what its reflection coefficient is), but a cloud of scatterers with coefficients calculated from the diffuse spot representing the PSF in the B-mode image, which amounts to adding the double convolution of the in vivo system to the double convolution simulated.

Figure 4 .

 4 Figure 4. Example images. Top: in vivo US common carotid images (R). Middle: corresponding in silico images (S), simulated with Field II. Bottom: corresponding in silico images ( R), post-processed by the GAN.

Table II PSNR

 II ON EACH SUBSET FOR THREE DIFFERENT ARCHITECTURES. THIRD COLUMN: PSNR BETWEEN THE in silico IMAGE WITHOUT POST-PROCESSING. FOURTH COLUMN: PSNR BETWEEN IN VIVO IMAGE AND THE IN SILICO IMAGE POST-PROCESSED BY THE GAN.different receptive fields (in meter). Since we aim to mimic the texture, the spatial information captured by the network is essential, and a pre-processing step to transform the database into isotropic and homogeneous data could be beneficial.VI. ACKNOWLEDGMENTSThis work was partly supported, via NL's doctoral grant, by the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). The authors have no relevant financial or nonfinancial interests to disclose.

	Model	Subset	PSNR original vs. simulated (dB) PSNR original vs. GAN output (dB)
		training	15.5 ± 3.5	24.7 ± 1.6
	U-net	validation	15.6 ± 3.4	22.2 ± 2.1
		test	15.5 ± 3.1	22.3 ± 1.8
		training	15.5 ± 3.5	25.7 ± 1.8
	Dilated U-net	validation	15.6 ± 3.4	23.4 ± 2.1
		test	15.5 ± 3.1	23.3 ± 2.0
		training	15.5 ± 3.5	17.1 ± 2.9
	Super-resolution GAN (ResNet-based) validation	15.6 ± 3.4	16.4 ± 3.6
		test	15.5 ± 3.1	16.2 ± 3.2

VII. COMPLIANCE WITH ETHICAL STANDARDS

INFORMATION

The data from human subjects used in this work were obtained and treated in line with the principles of the Declaration of Helsinki. Approval was granted by the Ethics Committees of the institutions involved in creating the multicentric database, from which these data were accessed.


Lainé

Guillaume Zahnd

Hervé Liebgott

Th Maciej Orkisz

Segmenting the carotid-artery wall in ultrasound image sequences with a dual-resolution U-net 1 st Nolann

Keywords: Atherosclerosis, carotid artery, segmentation, ultrasound, deep learning

Thickening of intima-media complex in the common carotid artery is a biomarker of atherosclerosis. To automatically measure this thickness, we propose a region-based segmentation method, involving a supervised deep-learning approach based on the dilated U-net architecture, named caroSegDeep. It was trained and evaluated using 5-fold cross-validation on two open-access databases containing a total of 2676 annotated images. Compared with the methods already evaluated on these databases, caroSegDeep established a new benchmark and achieved a mean absolute error twice smaller than the interobserver variability.

I. INTRODUCTION

World Health Organization reports the cardiovascular diseases, particularly atherosclerosis, as the leading cause of death worldwide [START_REF] Kaptoge | World health organization cardiovascular disease risk charts: revised models to estimate risk in 21 global regions[END_REF]. Their prevention [START_REF] Mcgill | Preventing heart disease in the 21st century: implications of the pathobiological determinants of atherosclerosis in youth (PDAY) study[END_REF] requires screening by means of a non-ionising and inexpensive imaging modality. Ultrasound (US) imaging has these characteristics and is routinely used to explore the common carotid artery (CCA), which is often considered as the sentinel of atherosclerosis [START_REF] Yousefi Rizi | Carotid wall longitudinal motion in ultrasound imaging: An expert consensus review[END_REF]. An early sign of this disease onset is arterial wall thickening. To measure the thickness of interest, the contours of the intimamedia complex (IMC), namely, lumen-intima (LI) and mediaadventitia (MA) interfaces, need to be identified (Fig. 1). This task is most frequently addressed by contour-based approaches [START_REF] Delsanto | Characterization of a completely user-independent algorithm for carotid artery segmentation in 2-D ultrasound images[END_REF], [START_REF] Loizou | Snakes based segmentation of the common carotid artery intima media[END_REF], [START_REF] Raj | Automated measurement of compression-decompression in arterial diameter and wall thickness by image-free ultrasound[END_REF], [START_REF] Zahnd | A fully-automatic method to segment the carotid artery layers in ultrasound imaging: Application to quantify the compressiondecompression pattern of the intima-media complex during the cardiac cycle[END_REF], which exploit intensity changes caused by the echoes at acoustic interfaces. Region-based segmentation methods based on texture, adaptive thresholding [START_REF] Nagaraj | Segmentation of intima media complex from carotid ultrasound images using wind driven optimization technique[END_REF], [START_REF] Wang | Fully automatic measurement of intima-media thickness in ultrasound images of the common carotid artery based on improved Otsu's method and adaptive wind driven optimization[END_REF] or clustering have also been proposed. Recently, deep- learning (DL) has been successfully used in vascular USimage segmentation to enhance the structures of interest prior to the actual delineation by more conventional contour-based methods [START_REF] Menchón-Lara | Early-stage atherosclerosis detection using deep learning over carotid ultrasound images[END_REF], [START_REF] Qian | Segmentation of the common carotid intima-media complex in ultrasound images using 2-D continuous maxflow and stacked sparse auto-encoder[END_REF], [START_REF] Shin | Automating carotid intima-media thickness video interpretation with convolutional neural networks[END_REF]. The drawback of these approaches is the necessity to combine a learnable pre-processing operation with an analytic segmentation task.

II. METHOD

The proposed solution builds on the U-net [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] architecture, with dilated convolutions on the bottleneck to increase the receptive field [START_REF] Meshram | Deep learning for carotid plaque segmentation using a dilated U-net architecture[END_REF]. It works in two steps and uses two identical networks with different weights (Fig. 2). First, the user's two mouse clicks define the region of interest (ROI) to be segmented. Next the far wall is detected. Eventually the IMC segmentation is done in the vicinity of the detected location. The code is available at https://github.com/nl3769/caroSegDeep.

A. Detection of the far wall

Similarly to many state-of-the-art methods [START_REF] Menchón-Lara | Early-stage atherosclerosis detection using deep learning over carotid ultrasound images[END_REF], [START_REF] Qian | Segmentation of the common carotid intima-media complex in ultrasound images using 2-D continuous maxflow and stacked sparse auto-encoder[END_REF], [START_REF] Wang | Fully automatic measurement of intima-media thickness in ultrasound images of the common carotid artery based on improved Otsu's method and adaptive wind driven optimization[END_REF], [START_REF] Zahnd | A fully-automatic method to segment the carotid artery layers in ultrasound imaging: Application to quantify the compressiondecompression pattern of the intima-media complex during the cardiac cycle[END_REF], the far wall is first detected to drive the subsequent actual IMC segmentation. The corresponding U-net, referred to as Θ F W , uses patches of full image height and 128-pixel width. The next paragraph describes the data preparation for the training phase of Θ F W , while the subsequent one specifies how the patch-wise predictions inferred using Θ F W are post-processed to obtain the curve approximately localizing the far wall on the entire ROI.

Pre-processing and training: Each image was resampled to a 512-pixel height. For training, the IMC median axis was defined as the line halfway between LI and MA annotations, extending across the entire width of the ROI. Pixels below the median axis were set to 1 and the others to 0, to generate a reference mask (M ROI ). Image ROI and M ROI were identically cut into patches with a 100-pixel overlap aiming at data augmentation. These gray-level and binary patches (Fig. 3) were the inputs of the training process.

Inference and post-processing: The same resampling and cutting of the ROI into overlapping 128 × 512-pixel patches is performed to prepare each image for inference. Then, each patch is segmented using Θ F W . Two maps are created:

• overlay map: containing, for each pixel, the number of overlapping patches covering this pixel, • prediction map: containing, for each pixel, the sum of values predicted by Θ F W within these patches.

Subsequently, an average value in the range [0, 1] is calculated for each pixel by dividing the prediction map by the overlay map. The averaged prediction map is binarized using a threshold of 0.5 and then cleaned by retaining the largest connected component. The median axis we seek is the upper boundary of the region thus obtained.

B. Segmentation of the IMC

The actual segmentation of the IMC uses the same network architecture (the dilated U-net used here is referred to as Θ IM C ) and several concepts from Section II-A: overlapping 128 × 512-pixel patches, overlay and prediction maps, as well as a similar post-processing except that two contours are extracted (the LI and MA interfaces). The following paragraphs focus on the specific choices made for this step.

Pre-processing and training: Here, the height of the patches just needs to encompass the IMC, which is 0.8 -millimeter thick, on average. Therefore, each image is vertically resampled to 5 µm/pixel, so that the 512-pixel patch height roughly corresponds to 2.6 mm. The ground truth for training was deduced from the images thus interpolated (Fig. 4): pixels located between the annotated LI and MA interfaces were set to 1, and the others to 0. The patches were picked along the median axis: at each abscissa x i , three patches were extracted, centered at y i and y i ± 128, where y i is the mean ordinate of the median axis on the patch width. This choice aimed at data augmentation coping with possibly inaccurate far-wall approximation, as well as with tilted arteries.

Inference and post-processing: The patches are also extracted along the far-wall approximation, which results here from the first step (Section II-A), and the number of vertically overlapping patches extracted at each abscissa x i may be three or more, depending on the tilt of the median axis. After combining the predictions made by Θ IM C in all patches into a prediction map, the segmentation map is derived thereof, as described before. The LI and MA interfaces are eventually deduced as the respective upper and lower boundaries of the region thus segmented.

III. RESULTS

To assess each network on data not seen during its training, a 5-fold cross-validation was performed: the database was divided into five subsets of equal size, and five pairs of Θ F W and Θ IM C networks were trained and tested independently using in turn three subsets for training (60%), one subset for validation (20%) and one subset for testing (20%). The results obtained on the test subsets of these five pairs were merged to assess the performance of the method on the entire database.

The second step (actual IMC segmentation) of the proposed cascade approach assumes that the median axis resulting from the first step (far wall detection) is approximately correct: if the far wall is wrongly detected, the ICM segmentation is likely to fail. Therefore, we first quantified the success rate of the first step alone, then incorrect median axes were manually redrawn using a home-made graphical interface before assessing the accuracy of the second step. Robustness of the far wall detection: Images where the distance between the predicted and reference median axes was greater than 0.8 mm were selected for visual inspection, which confirmed that 36 of them (1.3% of the database) were failures, i.e. curves unusable to initialize the IMC segmentation step. Hence, the success rate was of 98.7%. In the 36 images with failures, the median axis was interactively reinitialized, in preparation for the subsequent segmentation step. Accuracy of the IMC segmentation: The IMC segmentation inaccuracies were calculated with respect to the annotations performed by the expert A1. We assessed the median absolute difference (MAD) on IMT, which is the target measure, as well as the Hausdorff distance on the LI and MA contours. These errors were compared with the inter-and intra-observer variabilities, A2 vs. A1 and A1' vs. A1, respectively. Table I summarizes the segmentation errors separately calculated for images from the CUBS 1 and CUBS 2 databases, which allows a comparison with other methods evaluated on these databases in previous studies [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): Technical considerations on an open multi-center analysis of computerized measurement systems for intimamedia thickness measurement[END_REF], [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): An open multicenter analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF]. Thus, the previous benchmark on the CUBS 1 database was a MAD value on IMT equal to 114 ± 117 µm, while on the CUBS 2 database the best conventional method achieved 139 ± 119 µm and another Unet-based method obtained 178 ± 120 µm.

IV. DISCUSSION

We proposed and fairly evaluated on open-access databases an almost-automatic (two user mouse-clicks) deep-learning method devised to extract the contours of the intima-media complex in longitudinal B-mode ultrasound images of the carotid artery. The proposed approach is patch-based, which allows for segmenting variable-width ROIs of without resizing the images.

The far-wall localization step was successful in all but 1.3% of the images. This robustness is a prerequisite for overall correct segmentation. The latter achieved errors smaller than the interobserver variability, both in terms of estimated thickness and contour location. Compared to the methods already evaluated on the same data [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): Technical considerations on an open multi-center analysis of computerized measurement systems for intimamedia thickness measurement[END_REF], [START_REF] Meiburger | Carotid ultrasound boundary study (CUBS): An open multicenter analysis of computerized intima-media thickness measurement systems and their clinical impact[END_REF], caroSegDeep established a new benchmark. As the proposed method is based on supervised learning, it has the potential to increase its performance by using larger and more diverse data for training.

In particular, the largest errors occurred in the presence of calcified plaques, which were present in few images. This work was oriented towards asymptomatic plaque-free subjects, images with plaques were not expected. Nevertheless, we expect that results might be improved by re-training the networks on a database enriched with such images. This avenue deserves investigation. Another direction might exploit multiple annotations available for each image to increase the robustness by learning how to account for such uncertainties.

VI. COMPLIANCE WITH ETHICAL STANDARDS

Fig. 1 .

 1 Fig. 1. Example longitudinal B-mode US image of a CCA. a) Far wall encompassed by a rectangle. b) Enlarged region detailing the IMC with its interfaces segmented by caroSegDeep: LI (red) and MA (blue).

Fig. 2 .

 2 Fig. 2. Outline of the proposed method. a) Input image. b) Interactive delimitation of the left and right borders of the ROI (two mouse clicks). c) Far wall detection: median axis (magenta) deduced from overlapping masks predicted within full image-height patches. d) IMC segmentation: LI (cyan) and MA (green) interfaces deduced from predicted overlapping masks within high-resolution patches distributed along the median axis.

1 Fig. 3 .Fig. 4 .

 134 Fig. 3. Schematic representation of image patches and corresponding masks used during the training phase of the far-wall detection network.

  partly supported, via NL's doctoral grant, by the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). The authors have no relevant financial or nonfinancial interests to disclose.

TABLE I SEGMENTATION

 I INACCURACIES (MEAN ± STD) FOR CUBS 1 AND CUBS 2 DATASETS: MEAN ABSOLUTE DIFFERENCE (MAD) FOR THICKNESS QUANTIFICATION (IMT) AND HAUSDORFF DISTANCE (HD) FOR CONTOUR LOCATIONS (LI AND MA). THE ERRORS OF THE PROPOSED CAROSEGDEEP METHOD WITH RESPECT TO REFERENCE ANNOTATIONS ARE COMPARED TO THE INTER-AND INTRA-OBSERVER VARIABILITIES.

	Measure	Compare	Dataset	CUBS 1	CUBS 2
		caroSegDeep vs. A1	99 ± 89	106 ± 89
	IMT: MAD (µm)	A2 vs. A1	206 ± 168	192 ± 166
		A1' vs. A1	144 ± 123	160 ± 140
		caroSegDeep vs. A1	320 ± 193	305 ± 197
	LI: HD (µm)	A2 vs. A1	380 ± 207	327 ± 138
		A1' vs. A1	357 ± 204	352 ± 140
		caroSegDeep vs. A1	287 ± 153	289 ± 147
	MA: HD (µm)	A2 vs. A1	351 ± 161	338 ± 184
		A1' vs. A1	319 ± 155	346 ± 185

http://dx.doi.org/10.17632/fpv535fss7.1

http://dx.doi.org/10.17632/m7ndn58sv6.1

INFORMATION

The data from human subjects used in this work were obtained and treated in line with the principles of the Declaration of Helsinki. Approval was granted by the Ethics Committees of the institutions involved in creating the multicentric database, from which these data were accessed.