Christophe Vernay
email: cvernay@solais.fr

Sébastien Pitaval

Philippe Blanc

Review of satellite-based surface solar irradiation databases for the engineering, the financing and the operating of photovoltaic systems

Keywords: surface solar irradiation, satellite-based database, meteorological ground station, photovoltaic, bankable yield report, monitoring

This paper explores the possibilities provided by satellite-based surface solar irradiation databases, for the attention of the photovoltaic actors, from the engineering and the financing to the operating of photovoltaic systems. First, the problematic of using meteorological ground stations for the determining of solar dataset or time-series is addressed with the example of the French national meteorological network and the focus made on the South-West part of France: the heterogeneity in terms of spatial and temporal representativeness along with the relevance of the delivered measurement are questioned. Then, an innovative synthesis of 17 satellite-based surface solar irradiation databases available so far is presented through the distribution of their corresponding features within three categories: spatial and temporal representativeness, data type in terms of component and format, and finally operating mode for the retrieving of the data in terms of price and accessibility. The results of this review shall help photovoltaic actors making the correlation between the available satellite-based databases and their specific needs.

Introduction

Solar resource assessment is a recurrent problematic faced by all the actors of the solar industry, either photovoltaic (PV) or thermal solar. This paper focuses on PV but can surely be extended to thermal solar technology such as concentrating solar power (CSP).

PV project developers first need to identify the geographical areas where solar resource would match first basic criteria, by focusing on the yearly assessment of global horizontal irradiance (GHI) for classical PV projects or direct normal irradiance (DNI) for concentrated PV (CPV) technology. Once a site has been identified, the same developers along with engineering offices and bank institutes, through their own technical advisors, need to have an estimation of the electrical yield that is likely to occur during the life-time of the project, typically 20 years. For that purpose, simulation software such as PVSyst (www.pvsyst.com) considers the behavior of the whole system, from the modules to the transformers, when facing specific conditions in terms of ambient temperature, wind velocity and incident irradiance, i.e. global tilted irradiance (GTI) for PV and DNI for CPV. Hourly time-step or even less (10-min) is presently the baseline for CPV projects by means of typical meteorological years (TMY) [START_REF] Vignola | Building a bankable solar radiation dataset[END_REF][START_REF] Espinar | D404.1: Report on the product S4: Typical meteorological year[END_REF], i.e. 8,760 values for a typical year. For PV projects, empirical models are commonly used to derive hourly values of GHI from daily values (cf. Aguiar's autoregressive Gaussian model [START_REF] Aguiar | TAG: a time-dependent, autoregressive, Gaussian model for generating synthetic hourly radiation[END_REF]), daily values from monthly values (cf. Aguiar's model [START_REF] Aguiar | A simple procedure for generating sequences of daily radiation values using a library of markov transition matrices[END_REF]) and finally GHI to GTI (Hays's or Perez's models [START_REF] Hay | Calculation of monthly mean solar radiation or horizontal and inclined surfaces[END_REF][START_REF] Perez | A new simplified version of the Perez diffuse irradiance model for tilted surfaces[END_REF]), as widely explained in [START_REF] Ineichen | Global irradiance on tilted and oriented planes: model validations[END_REF]. Therefore, for the pre-feasibility studies of classical PV projects, the monthly GHI is the minimum information that is required; such a minimum data consists in a long-term average 12-month dataset, i.e. one mean value per month. Moreover, during the operating of the solar plant, the same actors along with operators and maintenance companies may need to evaluate the estimation of the local irradiance for a specific period, on a daily or hourly basis, in order to correlate the monitored electrical production with the effective resource and therefore to assess the overall performance of the system. As an illustration, these timeseries may be required to face the outage in the monitoring of the irradiance sensors or to help in the extrapolation of the local measurement from one set of azimuth and tilt to another one. Finally, the emergence of harsh dispute between insurance companies magnifies the need for surface solar irradiance (SSI) database that must be reliable, flexible and up-to-date; this is typically the case when the abnormal sub-performance of a PV plant leads to the triggering of performance guarantee.

The need for assessing SSI has been addressed for several decades by the scientific community which has setup different methods and databases based on different types of observations: physical measurements from meteorological ground stations such as pyranometer and pyrheliometer, reanalysis of numerical weather prediction models and climate models, and finally images from geostationary satellites. Even though many papers have been published regarding most of the databases, their heterogeneity make the synthesis complex for PV actors who need SSI information: representativeness of the database, nature and sources of the data that is made available, and finally conditions for retrieving the data. This paper aims at addressing this problematic by making the link between the supply and the demand through a review of the available satellite-based products in order to eventually help PV actors in their choices.

Meteorological ground stations

Main networks measuring solar radiation

Solar radiation has been measured by national networks for several decades. Part of the collected data has been processed, in the scope of research programs, in order to setupalasfew reliable networks with heterogeneous spatial coverage [START_REF] Budig | Comparison of different sources of meteorological data for Central Asia and Russia[END_REF]: from about 50 ground stations for Baseline surface radiation network (BSRN [START_REF] Ohmura | Baseline Surface Radiation Network (BSRN/WCRP): new precision radiometry for climate research[END_REF]) and International daylight measurement program (IDMP, http://idmp.entpe.fr/) up to more than 1,000 for Global energy balance archive (GEBA [START_REF] Gilgen | Means and Trends of Shortwave Irradiance at the Surface Estimated from Global Energy Balance Archive Data[END_REF]) and World radiation data center (WRDC, http://wrdc.mgo.rssi.ru/). As part of scientific research programs, the added-value for such databases mainly consists in the high-quality of the available data, considering that stringent requirements have been defined:  The quality of the sensors;  The periodic on-site maintenance along with sensor calibration;  The quality-check procedure to ensure the reliability of the measurement [START_REF] Roesch | Assessment of BSRN radiation records for the computation of monthly means[END_REF].

Despite such a quality check, data aggregation (yearly, monthly, daily and even hourly) always needs to be viewed with a critical eye because of the missing data and the shading problematic that may differ from one ground station to another.

Major obstacles prevent PV actors from working with these databases. Firstly, the spatial distribution is likely not appropriate to the need; as an example, in continental France (543,965 km²), there are respectively 2, 2, 40 and 131 ground stations belonging to BSRN, IDMP, GEBA and WRDC networks. The use of such ground-based SSI dataset involves in general the use of spatial and temporal interpolation and extrapolation whose quality drastically decreases with the distance of the site of consideration with the nearest ground stations. Moreover the operating mode for retrieving the data is not straightforward and requires skillful and time-consuming processing which is not compliant with the time and cost constraints of any PV project. Finally, as part of Research programs, these databases are exclusively reserved for non-commercial activities with the exception of GEBA (whose data can be accessed for commercial applications via meteonormwww.meteonorm.com), thus mitigating the added-value for project developers or engineering offices with commercial activities.

Such obstacles may partially disappear when considering the original national meteorological networks which more recently operate commercial activities based on their large number of ground stations measuring solar radiation. The case of the French national network operated by Météo France illustrates that the promising added-value might not be there due to several problems: spatial density of the station network, continuity and persistence of the temporal coverage, shadowing effects and maintenance quality.

Problematic of national networks: the example of Météo France network

Météo France is the French meteorological network which consists in 3,811 meteorological ground stations located in the country, among which only 327 used to or still measuring GHI. It is also to be noted that only 7 of them measured in the past the direct component of the SSI and, presently, only one is still measuring it. The analysis presented below focuses on the South-West quarter of the country (130,132km², 20% of the territory) whose irradiance level is suitable to PV project and where 61 ground stations measure this component (cf. Fig 1). [START_REF] Perez | Comparing satellite remote sensing and ground network measurements for the production of site/time specific irradiance data[END_REF][START_REF] Muselli | Improved procedure for stand-alone photovoltaic systems using Meteosat satellite images[END_REF][START_REF] Zelenka | Effective accuracy of satellite-derived hourly irradiances[END_REF]; nevertheless it is to be noted that in some specific locations, for example without noticeable orography, monthly and yearly temporal aggregations may lead to a typical acceptable distance to the nearest ground station greater than 25 km. The spatial heterogeneity of Météo France network in the South-West of France is therefore questionable when considering the basic engineering of PV projects. In addition to spatial and temporal representativeness, one shall always keep in mind the problematic associated with the delivery of ground based measurements: besides the quality of the pyranometers and the maintenance of the ground station which must be blameless, the location of the sensor within its local environment is of big concern. The impact of any close shading is obviously important and must be considered These figures show that only a quarter of the 61 sites are part of the highest class and that information is missing for less than 20% of the sites, thus reducing the representativeness of this database. This is all the more relevant when considering that such a classification does not take into account the impact of the surrounding relief [START_REF] Leroy | Classification d'un site. Note technique n°35[END_REF].

Finally, the quality of the measurement must be systematically questioned, whatever the database is. What is the protocol used for the quality check? What is the protocol used for data aggregation when values are missing (e.g. hour-to-day)? These questions are addressed by the 4 renowned meteorological networks presented before (BSRN, etc.) with a promising transparency through scientific publications. Such may not be the case for national meteorological networks which must therefore be handled with care.

Satellite-based databases

Description

The main part of the present study consisted in analyzing the available satellite-based databases that may partially or totally fulfill PV actors' needs. Table 1 below presents the overview of 17 databases available so far and which are still in operating condition. This list is surely not exhaustive and does not consider the databases with national coverage, nor the databases performing the reanalysis of weather forecasts such as NCEP / NCAR (National centers for environmental prediction / National center for atmospheric research) or ERA-Interim.

Each database has been studied based on the information provided by different sources: web service and peerto-peer publication as listed in Table 1, other publications and finally either phone call or email exchange for the missing information. The gathered information is structured using three categories in order to address the main concerns of the actors in the PV domain:

"Representativeness" of the database deals with the geographical area covered by the satellite, its spatial and temporal resolution and the temporal coverage; the latter is related to both the duration, in year, and the updating strategy which may reveal the potential obsolescence of the database. "Data type" deals with the SSI data itself, i.e. the product that is delivered to the end-user. It consists in the available component of the solar radiation (GHI, GTI and diffuse -DHIfor PV projects, DNI for CPV), the possibility to retrieve a dataset for a typical year through 12 monthly values and/or to work on time-series for a specific period and time sampling. The possibility to get TMY is not addressed in this study because it involves a complex processing of the raw solar radiation estimation along with the retrieval of other relevant meteorological parameters such as the air temperature or the wind speed. "Operating mode" finally focuses on the process used for retrieving the data, from the availability of a graphical user interface (GUI), web services and/or offline request through email exchange, to the financial conditions for accessing the data. The accuracy of each database has been intentionally omitted because this indicator may vary from one study to another, depending on the considered location (spatial variability) and the time period (inter-annual variability); it is therefore very difficult to analyze and almost impossible to draw general conclusions. However, a large number of publications will help people having a good view of database accuracy, among which [START_REF] Ineichen | Global irradiation: average and typical year, and year to year annual variability[END_REF][START_REF] Ineichen | Five satellite products deriving beam and global irradiance validation on data from 23 ground stations[END_REF].

The detail of the analysis of the 17 databases is available in the appendix.

Synthesis

Table 2 below presents the scoring scheme that has been applied for each category, thus resulting to Fig. 3 where the color gives an indication on how a database stands among the whole databases, from green for the highest score to red for the lowest score.

Conclusion

The first conclusion that can be drawn from the overall analysis of these 17 SSI databases is the absence of free database that would fulfill the whole needs presented before. Therefore, PV actor must consider the commercial products presented in this paper as a powerful mean for answering part or the totally of its needs. However, it can be mentioned that:  For the performance assessment of operating PV plants, 6 out of 17 databases (35%) provide GTI time-series at a time step lower or equal to one day;  For the processing of TMY in the scope of CPV studies, 10 out of 17 databases (59%) provide DNI timeseries at a time step lower or equal to one hour.

The temporal representativeness of any database must be emphasized in order to match the required analysis. Long-term databases allow catching the typical on-site irradiance through dataset or TMY along with the interannual variability [START_REF] Vignola | Building a bankable solar radiation dataset[END_REF][START_REF] Ineichen | Interannual variability and global irradiance evaluation[END_REF]. These are key information for the processing of bankable PV yield reports through, for example, P50 and P90 indicators. The added-value of a long-term database that is no longer refreshed is also questionable considering that a yield report aims at estimating the solar resource for the lifetime of a future PV project, i.e. the next 20 years to come. Finally, the ability to deliver recent data through, for example, a day+1 updating scheme is crucial for the online monitoring of a PV plant for in situ calibration process for which local measurement campaign must be properly monitored and potentially ended as soon as the convergence criteria of the calibration process is reached [START_REF] Vernay | Characterizing measurements campaigns for an innovative calibration approach of the global horizontal irradiation estimated by HelioClim-3[END_REF].

The list of SSI databases presented in this paper may not be exhaustive but it is certainly representative of the available databases so far. When assessing the typical irradiance of a given site, the averaging of several databases is therefore strongly questionable: what is the rationale for averaging the results of several databases which have different spatial or temporal resolutions and temporal coverage? Indeed, the uncertainty of the resulting SSI assessment is directly impacted by the interannual variability in case of databases with different temporal coverage and/or by the intra-pixel spatial variability in case of databases with different spatial resolution. In order to be scientifically consistent, one shall definitely focus on a single database whose characteristics are understood and then refine the satellite-based estimation through a calibration using local measurements [START_REF] Vernay | Characterizing measurements campaigns for an innovative calibration approach of the global horizontal irradiation estimated by HelioClim-3[END_REF].

For all these reasons, it is highly recommended to go deeper in the exhaustive results provided in the appendix in order to better understand the differences between each database and to objectively exploit the most appropriate one depending on the need.

Fig. 1 .

 1 Fig. 1. Location of the 61 ground stations measuring GHI in the South-West quarter of France.

Fig

  Fig 2a belowdepicts in blue the temporal coverage during which GHI measurements are available on an hourly basis; it starts in 1988 for 5 stations while the measurement of daily GHI starts in 1964. The figure highlights the heterogeneity in the duration which is between 2.5 and 25 years, with median equal to 10.5 years, along with, in some rare cases, the potential obsolescence with 10 stations which don't measure GHI anymore; discontinuity of the measurement is also highlighted with 8 stations facing several interruption periods. The spatial representativeness of the network is also addressed in Fig2bwhich presents the histogram of the geodesic distance between one ground station and its nearest neighbor. The median distance is equal to 40 km but it can go up to 166 km; this is typically the case for rural areas which are paradoxically more suitable to groundmounted PV plants. Literature claims that, considering the spatial variability of the SSI compared to the typical satellite-based SSI estimation error, it is better retrieving hourly and daily irradiation from satellite databases rather than interpolating measurements from surface stations 20 to 25 km apart[START_REF] Perez | Comparing satellite remote sensing and ground network measurements for the production of site/time specific irradiance data[END_REF][START_REF] Muselli | Improved procedure for stand-alone photovoltaic systems using Meteosat satellite images[END_REF][START_REF] Zelenka | Effective accuracy of satellite-derived hourly irradiances[END_REF]; nevertheless it is to be noted that in some specific locations, for example without noticeable orography, monthly and yearly temporal aggregations may lead to a typical acceptable distance to the nearest ground station greater than 25 km. The spatial heterogeneity of Météo France network in the South-West of France is therefore questionable when considering the basic engineering of PV projects.

Fig. 2 (

 2 Fig. 2 (a) temporal coverage for hourly GHI measurement; (b) histogram of the geodesic distance.

Fig. 3 .

 3 Fig. 3. Synthesis of the scoring of the 17 SSI databases.

Table 1 .

 1 Overview of the satellite-based databases

	Name	Provider	URL	Coverage	Ref.
	3TIER	3TIER	3tier.com	Global	-
	CM SAF	DWD	cmsaf.eu	Europe, Africa	[16]
	DLR-ISIS	DLR	www.pa.op.dlr.de/ISIS/	Global	[17]
	EnMetSol	University of Oldenburg	energy-meteorology.de	Europe, part of Africa	[18]
	focus solar	focus solar GmbH	focussolar.de	North-America, Europe,	-
				Africa, part of Asia,
				Oceania
	HelioClim-1	MINES ParisTech -Armines	soda-is.com	Europe, Africa, part of	[19]
				South-America
	HelioClim-3	MINES ParisTech -Armines	soda-is.com	Europe, Africa, part of	[20, 21]
				South-America
	IrSOLaV	IrSOLaV SL	irsolav.com	Global	[22]
	Land SAF	European consortium	landsaf.meteo.pt/	Europe, Africa, part of	[23]
				South-America
	meteonorm	METEOTEST	meteonorm.com	Global	[24]
	NASA SSE Release 6.0	NASA	eosweb.larc.nasa.gov/sse/	Global	-
	OSI SAF	European consortium	osi-saf.org	Europe, Africa,	-
				America
	PVGIS CM-SAF	JRC	re.jrc.ec.europa.eu/pvgis/	Europe, Africa	[25]
	Satel-Light	ENTPE	satel-light.com	Europe, North-Africa	-
	SolarAnywhere	Clean Power Research	solaranywhere.com	North-America, Hawaii,	[26]
				Caribbean
	SolarGIS	GeoModel Solar	solargis.info	Global	-
	Solemi	DLR	dlr.de/tt/solemi	Europe, Africa, part of	-
				South America,
				Indian Ocean, Asia

Table 2 .

 2 Details of the scoring scheme

	Category	Sub-category	score=1	score=2	score=3	score=4
			lowest			highest
	Representativeness Spatial coverage	Smaller area	Mono-satellite	Multi-satellite	Global
		Spatial resolution	Lowest to highest resolution using a logarithm scale
		Timeframe	Lowest to highest period using a logarithm scale
		Updating schedule	None	Potentially	Periodically	Day+1
	Data type	Components	GHI only	GHI & DNI	Partially, among	All components
					GHI, DHI, GTI & DNI
		Time step	Monthly	Dayly	Hourly	Sub-hourly
		Aggregation	Dataset only	-	Time-series only	Dataset & time-
						series
	Operating mode	Accessibility	Offline request	FTP access	Partial GUI/web service	GUI / web
			only (Email)	only	otherwise offline	service
					request
		Price	Commercial	-	Free for scientific usage Free

Appendix A. Detailed results

A. [START_REF] Vignola | Building a bankable solar radiation dataset[END_REF]


Razmig Kéchichian

Sébastien Valette

Michel Desvignes

Rémy Prost

Keywords: inverse problems, image data compression, multiresolution algorithms, wavelets, shape modeling and meshes processing

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Shortest-Path Constraints for 3D Multiobject Semiautomatic Segmentation via Clustering and Graph Cut

Razmig Kéchichian, Sébastien Valette * , Michel Desvignes, and Rémy Prost, Member, IEEE Abstract-We derive shortest-path constraints from graph models of structure adjacency relations and introduce them in a joint centroidal Voronoi image clustering and Graph Cut multiobject semiautomatic segmentation framework. The vicinity prior model thus defined is a piecewise-constant model incurring multiple levels of penalization capturing the spatial configuration of structures in multiobject segmentation. Qualitative and quantitative analysis and comparison with a Potts prior-based approach and our previous contribution on synthetic, simulated and real medical images show that the vicinity prior allows for the correct segmentation of distinct structures having identical intensity profiles and improves the precision of segmentation boundary placement while being fairly robust to clustering resolution. The clustering approach we take to simplify images prior to segmentation strikes a good balance between boundary adaptivity and cluster compactness criteria furthermore allowing to control the trade-off. Compared to a direct application of segmentation on voxels, the clustering step improves the overall runtime and memory footprint of the segmentation process up to an order of magnitude without compromising the quality of the result.

Index Terms-Image segmentation, image clustering, Markov random field, spatial prior, Graph Cut.

EDICS categories-TEC-BIP, SMR-SMD, SMR-STM.

I. INTRODUCTION

C OMPUTER-AIDED medical image analysis has tradi- tionally focused on single-organ or pathology applications. Recent technological and algorithmic advances have brought increasing interest in simultaneous analysis and segmentation of multiple anatomical structures for comprehensive diagnosis and preoperative planning [1]- [START_REF] Linguraru | Statistical 4D graphs for multi-organ abdominal segmentation from multiphase CT[END_REF]. Moreover, the creation of full-body patient-specific models for the semantic navigation of anatomy is becoming a popular application [START_REF] Blume | Google body: 3D human anatomy in the browser[END_REF].

Clinical practice, especially whole-body Computed Tomography (CT) scanning, often generates large numbers of highresolution images, which makes tasks of efficient data access, transfer, analysis and visualization challenging, especially in Manuscript received on September 26, 2012; revised on January 21 and April 10, 2013. This work is supported by the Région Rhône-Alpes (France) via the SIMED project of the research cluster ISLE. R. Kéchichian, S. Valette and R. Prost are with Creatis research center: Université de Lyon; CNRS UMR 5220; Inserm U1044; INSA-Lyon; Université Lyon 1, Villeurbanne, France (e-mail: firstname.lastname@creatis.insalyon.fr).

M. Desvignes is with Gipsa-lab: CNRS UMR 5216; Grenoble-INP; Université Joseph Fourier; Université Stendhal, Saint-Martin-d'Hères, France (email: michel.desvignes@gipsa-lab.grenoble-inp.fr).

Digital Objet Identifier today's distributed computing environments which have seen increasing use of handheld terminals for interactive data access and visualization of anatomy.

CT and Magnetic Resonance (MR) images have intrinsic characteristics that render automatic analysis particularly difficult. In theory, they are regarded as piecewise-constant intensity maps over a number of tissue classes. However, contrast between tissue class intensities, which facilitates segmentation, depends on imaging conditions which are seldomly ideal. In practice, the piecewise-constance property is considerably degraded by various noise sources and the partial volume effect (multiple tissue-class occupancy within a voxel) which creates diffuse edges making it difficult to identify true structure boundaries. MR images may, in addition, suffer from spatial distortion of tissue intensity due to main magnetic field inhomogeneity. Regardless of the imaging modality and related artifacts, many anatomically and functionally distinct structures can have similar intensity levels in images and, furthermore, blend into surrounding tissues having intensities close to their own. It is impossible to identify and segment such structures on the basis of intensity information only.

Semiautomatic segmentation methods are usually resorted to in order to alleviate prohibitive time and effort requirements of manual delineation and as an alternative to fully-automatic methods when the latter are not applicable, for example due to the singularity of the studied case or the lack of training sets. Furthermore, automatic methods which rely on image and anatomical features learned beforehand from presegmented data, while effective in many cases, are prone to be bound to the characteristics of the training set and its variabilities, and may therefore have difficulty in accounting for previously unencountered but possible pathological cases and photometric characteristics of images acquired under different conditions.

A wide variety of semiautomatic segmentation approaches relying on different interaction paradigms have been proposed (see references in [START_REF] Rother | GrabCut": interactive foreground extraction using iterated graph cuts[END_REF], [START_REF] Vezhnevets | Grow-Cut" -interactive multi-label N-D image segmentation[END_REF]) where a reasonable amount of user input can initialize and steer the segmentation process refining its result by resolving erroneous decisions made by the algorithm. In particular, when segmentation is formulated as a voxel labeling problem where each segmented object is represented by a label class, user interaction can be carried out by the attribution of "seed" labels to voxels inside targeted structures. This provides clues on what the user intends to segment and can be used, for example, to collect intensity 0000-0000/00$00.00 c 201X IEEE statistics of targeted objects and their relative positions, and to constrain the solution space of the algorithm by prohibiting the attribution of other labels to seed-marked voxels.

A. Graph Cut in segmentation and its complexity

A fundamental problem of image processing and computer vision, image segmentation is the back-bone of many applications, especially in medical imaging. Graph Cut methods, which have been widely applied to single-object segmentation problems [START_REF] Boykov | Graph cuts and efficient N-D image segmentation[END_REF], rely on a maximum-flow binary global optimization scheme capable of finding the optimal solution in mathematical sense [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF]. With few exceptions [START_REF] Ishikawa | Exact optimization for markov random fields with convex priors[END_REF], multiobject generalizations of this approach do not, in general, have such optimality properties due to the NP-hardness of the underlying problem, but some algorithms [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] are nevertheless capable of finding a local solution optimal within a known factor of the global optimum. Simultaneous multiobject segmentation approaches, on the other hand, are superior to their sequential counterparts in that they raise questions neither on the best segmentation sequence to follow nor on how to avoid the propagation of errors on individual segmentations [START_REF] Fouquier | Sequential model-based segmentation and recognition of image structures driven by visual features and spatial relations[END_REF].

Graph Cut algorithms usually employ rectangular grid graphs for data structure. In applications on high-resolution 3D medical images where several tens of millions of voxels are common, such graphs have to be defined on a number of vertices equal to the number of image voxels and a number of edges at least five times as much. With data structure size touching on the limit or exceeding the storage capacity of computer memory, numerical optimization becomes impractical.

Approaches which address the complexity of Graph Cut methods on high-resolution data can be broadly classified into two categories; algorithm parallelization and image simplification. Parallel versions of some maximum flow algorithms have been devised for multiprocessor architectures [START_REF] Delong | A scalable graph-cut algorithm for N-D grids[END_REF] and graphics processing units [START_REF] Vineet | CUDA cuts: Fast graph cuts on the GPU[END_REF]. While these methods attain good speedups, they do not reduce memory footprint since they operate on the entire image. Image simplification approaches address this problem directly observing that the pixel representation is often redundant because objects usually comprise many similar pixels that could be grouped. In [START_REF] Lombaert | A multilevel banded graph cuts method for fast image segmentation[END_REF], Graph Cut segmentation is applied to a low-resolution image propagating the solution to the full resolution where Graph Cut segmentation is applied in a narrow band surrounding the projected object boundary. In order to accelerate user feedback in interactive segmentation, [START_REF] Li | Lazy snapping[END_REF] oversegment the image via the Watershed algorithm. Graph Cut segmentation is then formulated according to the watershed regions graph. Starting from a grid partition, [START_REF] Cigla | Efficient graph-based image segmentation via speeded-up turbo pixels[END_REF] cluster image pixels by an iterative k-means algorithm augmented by color similarity and shape compactness criteria. The cluster graph is then partitioned based on color information yielding a coarse segmentation of the image. Recent conditional random field-based (CRF) segmentation approaches such as [START_REF] Fulkerson | Class segmentation and object localization with superpixel neighborhoods[END_REF], [START_REF] Lucchi | Are spatial and global constraints really necessary for segmentation?[END_REF] rely on image clustering not only to reduce memory overhead but also to collect image features from clusters and their neighborhood.

B. Prior information in Graph Cut segmentation

Most successful segmentation methods usually incorporate some image or domain-specific prior information. Medical image segmentation methods should ideally be able to incorporate a wide range of anatomical and physiological prior information to ensure the consistency of the result with respect to anatomical properties and variabilities thereof. From the standpoint of Graph Cut optimization, prior information can be regarded as a set of constraints defined by a reference model.

Many forms of prior information have been introduced in Graph Cut segmentation for different applications. We broadly classify the most relevant of such approaches in three categories according to the nature of the introduced prior.

a) Location prior, where spatially-varying prior probability maps of structure locations, also called probabilistic atlases, are used to guide segmentation [START_REF] Linguraru | Statistical 4D graphs for multi-organ abdominal segmentation from multiphase CT[END_REF], [START_REF] Song | Adaptive graph cuts with tissue priors for brain MRI segmentation[END_REF]. In simplest cases, these are defined as distance maps [START_REF] Freedman | Interactive graph cut based segmentation with shape priors[END_REF]. Such probabilistic maps need to be learned in advance and most have to be registered to an image before it is segmented. Registration cannot trivially handle topological changes, due to pathologies for example, and is often computationally intensive.

b) Shape prior, where explicit or implicit shape constraints are used to match the segmented object with predefined or learned shapes, such as in [START_REF] Vu | Shape prior segmentation of multiple objects with graph cuts[END_REF] where an affine-invariant shape similarity measure is used in sequential multiobject Graph Cut segmentation, and in [START_REF] Malcolm | Graph cut segmentation with nonlinear shape priors[END_REF] where a statistical model for a shape and its variabilities is learned from a set of training shapes via kernel principal component analysis. Shape-prior segmentation approaches entail shape location initialization and alignment mechanisms, and often require model learning.

c) Spatial configuration prior, the use of which arises in the context of multiobject segmentation, where (usually pairwise) relations define interobject geometric and spatial relationships. It has appeared in scene layout estimation and segmentation where spatial relationships are predefined on a limited number of classes [START_REF] Liu | Order-preserving moves for graph-cut-based optimization[END_REF], and in part-based object detection and segmentation where such relationships are learned from annotated training data [START_REF] Lucchi | Are spatial and global constraints really necessary for segmentation?[END_REF], [START_REF] Winn | The layout consistent random field for recognizing and segmenting partially occluded objects[END_REF]. [START_REF] Delong | Globally optimal segmentation of multiregion objects[END_REF] segment multiple or compound objects encountered in relative positions of containment and exclusion, also introducing constraints on preferred distances between the boundaries of regions. Some useful configurations however cannot be represented by definitions therein, such as the inclusion of an adjacent pair of objects in a third.

There are a number of prior introduction approaches closely related to our work which have appeared outside Graph Cut segmentation literature. [START_REF] Cuadra | Comparison and validation of tissue modelization and statistical classification methods in T1-weighted MR brain images[END_REF] use a simple three-level piecewise-constant prior model for tissue classification in brain MR images. According to a brain tissue model, it encourages the identical classification of pairs of neighboring pixels, allowing different class adjacency patterns which are consistent with the model while penalizing adjacencies which are not. We use a similar prior for Graph Cut multiobject segmentation in our earlier work [START_REF] Jacinto | A web interface for 3d visualization and interactive segmentation of medical images[END_REF]. In [1], pairwise geometric relationships are defined by an anatomical network on easily detectable body landmarks and used in the organ detection phase of a sequential segmentation framework for thoracic and abdominal organs in whole-body CT scans. [START_REF] Okada | Abdominal multi-organ segmentation of CT images based on hierarchical spatial modeling of organ interrelations[END_REF] address the same problem with an automatic two-stage method whereby once relatively stable organs are identified and segmented, organ interrelations are used to proceed with the segmentation of variable organs by fitting previously learned statistical organ atlases and shape models. Organ interrelations are defined according to canonical correlation analysis and represented by a directed graph. Intuitively, these relations represent the degree of influence of an organ's surface upon another. A generic, graph-based scene model representing spatial information explicitly is also used by [START_REF] Fouquier | Sequential model-based segmentation and recognition of image structures driven by visual features and spatial relations[END_REF] in the problem of segmenting multiple subcortical structures in brain MR images. Fuzzy representations of distance, adjacency and spatial directionality relations are used to account for model imprecisions and variabilities. However, the segmentation algorithm is sequential and requires error detecting and segmentation order rearranging strategies.

C. Contribution

We propose a piecewise-constant prior model incurring multiple levels of penalization capturing the spatial configuration of structures in multiobject segmentation. The vicinity prior, as we term it, is defined as shortest-path pairwise constraints on a graph model of interobject adjacency relations.

We introduce the vicinity prior, along with user-supplied constraints, into a joint centroidal Voronoi image clustering and Graph Cut multiobject semiautomatic segmentation framework. The clustering method itself we use to simplify images prior to segmentation provides a good balance between the conflicting goals of boundary adaptivity and cluster compactness furthermore allowing to control the trade-off. Compared to a direct application of segmentation on voxels, the clustering step improves overall runtime and memory footprint of the segmentation process up to an order of magnitude without compromising the quality of segmentation in practice.

Qualitative and quantitative analysis and comparison with a standard Potts prior-based approach and our previous method [START_REF] Jacinto | A web interface for 3d visualization and interactive segmentation of medical images[END_REF] on synthetic, simulated and real medical images confirm the advantages of the vicinity prior particularly in the correct segmentation of distinct structures having similar intensities, the accurate placement of segmented structure boundaries even with coarsely clustered images subject to high levels of noise, and the robustness of segmentation to clustering resolution.

D. Paper organization

Section II outlines the image-adaptive centroidal Voronoi tessellation and compares it with two closely-related algorithms [START_REF] Levinshtein | TurboPixels: fast superpixels using geometric flows[END_REF], [START_REF] Veksler | Superpixels and supervoxels in an energy optimization framework[END_REF]. The statistical segmentation approach is described in Section III. In Section III-A, we present our graphbased vicinity prior model followed by data-cost definition in Section III-B. We briefly describe the numerical optimization method we employ in Section III-C and establish the necessary condition of its applicability. Results of qualitative and quantitative evaluation on synthetic, simulated and real medical images are given in Section IV. In Section IV-A we present evaluation results on a dataset generated from our synthetic phantom (Fig. 3a). Next, in Section IV-B1 we give quantitative evaluation results on the BrainWeb simulated MRI dataset [START_REF] Collins | Design and construction of a realistic digital brain phantom[END_REF] and compare them to brain tissue classification results reported in [START_REF] Gorthi | Comparison of energy minimization methods for 3-D brain tissue classification[END_REF]. Finally, in Section IV-B2 we present a segmentation case study comprising 34 thoracic and abdominal structures in a CT image and conclude the paper in Section V with a brief discussion outlining future venues of research.

II. IMAGE CLUSTERING BY CENTROIDAL VORONOI

TESSELLATION Several image clustering algorithms have been developed by vision and pattern analysis communities for image classification and segmentation. Of special interest is a class of algorithms which produce a dense oversegmentation of compact clusters, often called superpixels, having relatively uniform size and shape, which furthermore adapt to local intensity edges. Two recent examples are the TurboPixels algorithm of [START_REF] Levinshtein | TurboPixels: fast superpixels using geometric flows[END_REF] and the Graph Cut superpixels algorithm of [START_REF] Veksler | Superpixels and supervoxels in an energy optimization framework[END_REF] which achieve a better balance between the conflicting goals of compactness and boundary adherence than some wellknown image partitioning algorithms which produce segments of highly variable shape and size, like Watershed [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF] and Mean-shift [START_REF] Comaniciu | Mean shift: a robust approach toward feature analysis[END_REF] algorithms. A large segment of irregular shape is more likely to span more than one object, especially in the absence of boundary cues with insufficient contrast. On the other hand, a compact regular shape is less likely to cross object boundaries unless they are too wiggly. However, compactness comes at the expense of boundary adherence, it is therefore desirable to be able to control the trade-off.

Before presenting our clustering approach, let us define a grayscale image I as a mapping I :

V → I from a voxel domain V = {v | v = (x, y, z)} to a set of gray levels I ⊂ R.
To simplify notation, we shall write v ∈ I, and denote the gray-level of a voxel by I v .

Given a grayscale image I and n sites c i ∈ I such that 1 ≤ i ≤ n, a Voronoi tessellation partitions the image I into n disjoint clusters C i associated with each site c i :

C i = {v ∈ I | d(v, c i ) < d(v, c j ); 1 ≤ j ≤ n, j = i}. (1)
In a centroidal Voronoi tessellation (CVT), each site c i corresponds to the mass centroid of the associated cluster, and the tessellation minimizes the following energy:

E(v; c i ) = n i=1 v∈Ci ρ(v) v -c i 2 (2)
where ρ(v) is a density function defined below. We extend the cluster geometry compactness property expressed by ( 2) by adding an intensity-space norm as well:

F (v; c i ) = n i=1 v∈Ci ρ(v) v -c i 2 + α I v -I i 2 ( 3
)
were α is a positive scalar and I i is the gray-level of the cluster C i defined as the mean intensity of its voxels. Intuitively, minimizing (3) now corresponds to maximizing cluster compactness in terms of both geometry and gray-level. In (3), we define ρ(v) as a linear function of intensitygradient magnitude at voxel v, ρ(v) = m|∇I v | + b. This encourages the alignment of clusters with intensity edges in the image and the formation of relatively small clusters nearby, thus allowing fine-grained placement of segmentation boundaries. In all results reported in this paper, we set α = m = b = 1. Refer to Fig. 1b for an illustration in 2D.

To minimize (3), we use a variant of the clustering algorithm in [START_REF] Dardenne | Variational tetraedral mesh generation from discrete volume data[END_REF] which approximates a CVT in a computationallyefficient manner, involving only local queries on voxels located on boundaries of pairs of clusters. To assess its runtime performance and output quality, we have compared it with the TurboPixels and Graph Cut superpixels algorithms, both of which have similar cluster compactness criteria, on the Berkeley segmentation dataset [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF]. This dataset comprises 300 generic photos of size 321 × 481 segmented by different human subjects. In quantitative analysis, we evaluate runtime, boundary recall and undersegmentation error on 80 grayscale images from the dataset, choosing, for each such image, one ground-truth segmentation that does not contain small segments of few pixels. Such segments can produce large outliers and disrupt mean undersegmentation error measurements. Boundary recall is defined as the fraction of ground-truth boundary pixels which fall within a small distance from at least one cluster boundary. We use a distance of one pixel since we seek high precision in segmentation boundary placement in medical applications. Given a CVT clustering of an image C = {C 1 , . . . , C n } and a ground-truth partition of it T = {g 1 , . . . , g k }, the undersegmentation error for segment g j is:

segerr j (C, T ) = i:Ci∩gj =∅ |C i | -|g j | |g j | (4)
For a given clustering, we average (4) over all segments of the corresponding ground-truth partition and give the mean error. Figs. 2a, 2c and 2e present results of quantitative evaluation of the CVT clustering algorithm, the TurboPixels algorithm and the Graph Cut superpixels algorithm variant using cluster intensity homogeneity constraints. For a given clustering resolution, we report mean performance measures over all images. We can see from the graphs that the CVT clustering algorithm outperforms both TurboPixels and Graph Cut superpixels algorithms in terms of runtime and produces output of comparable quality to the best performing algorithm of the remaining two. Qualitative comparison of the output of the three clustering algorithms on a Berkeley dataset photo is given in Figs. 2b, 2d and 2f. We note that evaluations were conducted on a 4-core processor running at 2.84 GHz using software implementations provided by respective authors of the three algorithms.

For discussion purposes in subsequent sections, we define the dual graph of a CVT, illustrated in Fig. 1b. In 2D, denote the boundary (the surface, in 3D) of a cluster C i by ∂C i . Given a CVT clustering C where clusters are indexed by S, let G = S, E be an undirected graph on cluster centroids where pairs of clusters having nonzero length (area) common boundary (surface) define the set of edges E = {i, j} | i, j ∈ S, |∂C i ∩ ∂C j | = 0 . Consequently, the neighborhood of a node i ∈ S is defined as

N i = j | j ∈ S, ∃ {i, j} ∈ E .

III. MULTIOBJECT SEGMENTATION AS BAYESIAN

LABELING

We shall formulate segmentation as a labeling problem, defined as the assignment of a label from a set of m labels L = {l 1 , . . . , l m } representing objects to be segmented to each of the variables in a set of n variables corresponding to the clusters of a CVT-clustered image C = {C 1 , . . . , C n } indexed by S. Assume furthermore that each variable i ∈ S is associated with the corresponding node in the dual graph G of the CVT defined at the end of Section II.

An assignment of labels to all variables is called a configuration, and is denoted by ℓ ∈ L, where L = L n is the space of all possible configurations. The assignment of a label to a single variable i ∈ S is denoted by ℓ i . The space of admissible configurations, which corresponds to possible segmentations, can be identical to L or, if user-specified constraints are taken into account, to a subset of L. Let X = {X 1 , . . . , X n } be a set of random variables on S taking its values in L. A configuration ℓ can also be expressed in probabilistic terms as the joint event X 1 = ℓ 1 , . . . , X n = ℓ n , abbreviated as X = ℓ.

The "best" segmentation is defined as the configuration corresponding to the highest mode of the posterior probability distribution of configurations Pr(X = ℓ | O) conditional on some image-derived observation O. According to the Bayes rule, the posterior probability can be stated in terms of a likelihood distribution Pr(O | X = ℓ) and a prior distribution of configurations Pr(X = ℓ). If the latter can be expressed as a first-order Markov random field (MRF) for label configurations with respect to G and image-derived likelihood densities can be defined for each label l ∈ L from user-supplied evidence, then the maximum probability configuration can be equivalently found by minimizing the following energy function [START_REF] Li | Markov random field modeling in image analysis[END_REF]:

E(ℓ) = t i∈S D i (ℓ i ) + i∈S j∈Ni V i,j (ℓ i , ℓ j ) (5)
where t is a hyperparameter called the temperature and N i is the neighborhood of the variable i ∈ S. Pairwise terms of (5) encode prior information on interactions between labels assigned to pairs of neighboring sites encouraging the spatial consistency of labeling with respect to a reference model. Unary terms, also called data terms, are negative log-likelihood functions derived from observed data and measure the cost of assigning a label to respective variables.

A. The vicinity prior

For image segmentation, the pairwise interaction prior in (5) has to encourage piecewise-constant labeling, partitioning the image into several regions where voxels within a region are identically labeled. Furthermore, it has to be discontinuity preserving, in that it should attribute identical labels only to image regions which do not contain intensity edges.

The earliest discontinuity-preserving piecewise-constant prior model is the Ising model in binary labeling where |L| = 2 [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF]. It is defined as V i,j (ℓ i , ℓ j ) = 1δ(ℓ i , ℓ j ), where δ is the Kronecker delta. Its generalization for multiple labels where |L| > 2 is called the Potts model, which has widely been used in image segmentation and continues to enjoy popularity in recent CRF-based approaches [START_REF] Fulkerson | Class segmentation and object localization with superpixel neighborhoods[END_REF], [START_REF] Shotton | TextonBoost for image understanding: Multi-class object recognition and segmentation by jointly modeling texture, layout, and context[END_REF]. We can see that this model incurs a single level of penalization corresponding to the attribution of different labels to pairs of neighboring variables and incurs no penalty for the identical labeling of such pairs. This penalization scheme stems from the very nature of binary labeling where only a single type of boundary exists between regions labeled by either label. We extend the definition of the Potts prior model to account for multiple types of boundaries in multiobject segmentation, incurring multiple levels of penalization. To motivate our approach, let us first give an illustration on a synthetic example. Consider the image in Fig. 3a. It is further described in Section IV-A where we give evaluation results that use it as ground-truth estimate. Notice how the structure labeled "9" is nested in structure "4" which, in turn, is adjacent to structure "6". Due to this adjacency pattern, it makes perfect sense to encourage the assignment of the pair of labels {9, 4} to neighbor variables in a segmentation task while penalizing the attribution of the pair {9, 6}. Hence, ∀i, j ∈ S, V i,j (9, 6) > V i,j [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF][START_REF] Linguraru | Statistical 4D graphs for multi-organ abdominal segmentation from multiphase CT[END_REF]. Furthermore, since the structure "8" is farther still, we can say V i,j (9, 8) > V i,j (9, 6) > V i,j [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF][START_REF] Linguraru | Statistical 4D graphs for multi-organ abdominal segmentation from multiphase CT[END_REF]. By contrast, the simple structure adjacency prior models employed in our earlier work [START_REF] Jacinto | A web interface for 3d visualization and interactive segmentation of medical images[END_REF] and in [START_REF] Cuadra | Comparison and validation of tissue modelization and statistical classification methods in T1-weighted MR brain images[END_REF] penalize nonadjacent pairings {9, 6} and {9, 8} equally implying that they are equally unfit.

Recall that the set of labels L represents the objects to be segmented. Let Given the graph A, we define the pairwise term in (5) as

V i,j ℓ i , ℓ j = ̟ l a , l b , ℓ i = l a , ℓ j = l b . (6)
where ̟ l a , l b is the shortest-path weight from l a to l b in A.

At present, we let relations r ∈ R represent structure adjacency. A structure is considered to be adjacent to another in 2D if they share a boundary (a surface, in 3D). For two such structures a and b, the obvious definition of the adjacency relation's characteristic function is r = 1 when ∃r ∈ R, l a r l b . This results in a connected undirected graph A with unitweight edges between adjacent vertices. Shortest paths from a vertex to all others, conveniently calculated by an application of breadth-first search to the graph A, define the vicinity of the associated structure with respect to all other structures, which is introduced in the energy function of segmentation (5) via pairwise interaction terms [START_REF] Rother | GrabCut": interactive foreground extraction using iterated graph cuts[END_REF]. Clearly, the graph A constitutes a prior model of explicit spatial configuration of objects. As an example, the graph prior model for the image structures in Fig. 3a is given in Fig. 3b.

Note that the graph A is not limited to unit-weight edges. It can encode any other quantity that accumulates linearly along a path that we would want to minimize.

In 2D, we weigh (6) by the length (area, in 3D) of the common boundary (surface) of adjacent clusters |∂C i ∩ ∂C j | so that the sum of pairwise energies in [START_REF] Blume | Google body: 3D human anatomy in the browser[END_REF] for any pair of labels l a , l b ∈ L is equal to the length (area) of the common boundary (surface) between the corresponding pair of objects multiplied by V (l a , l b ). Along with our data energy definition [START_REF] Vezhnevets | Grow-Cut" -interactive multi-label N-D image segmentation[END_REF], this ensures that the energies of identical label configurations on CVT clusterings of different resolution are equal. Due to limited space, we cannot give the formal proof. Nevertheless, we establish the advantage of this weighting scheme in Section IV-A by quantitative comparison with the inverse distance weighting method frequently encountered in MRF-based segmentation, for example in [START_REF] Boykov | Graph cuts and efficient N-D image segmentation[END_REF] and [START_REF] Cuadra | Comparison and validation of tissue modelization and statistical classification methods in T1-weighted MR brain images[END_REF].

B. Likelihood-based data cost

The likelihood probability for a variable i ∈ S given the label assignment ℓ i is given by Pr(O i = I i | X i = ℓ i ), where the observation O constitutes the spatially distinct records of image intensity levels for each variable. Given user-supplied "seeds", 1 we estimate the conditional probability distribution of intensity Pr(I | L = l) for every label l ∈ L as a Gausssmoothed and normalized intensity histogram.

To facilitate the understanding of different cases, we define the data term D i (ℓ i ) of the energy (5) as follows:

D i (ℓ i ) =        0, ∃s k ∈ C i , k = ℓ i (7a) ∞, ∃s k ∈ C i , k = ℓ i (7b) -ln v∈Ci Pr(I v | ℓ i ), ∃s k ∈ C i , ∀k ∈ L (7c)
where s k represents a voxel marked with a seed of label type k ∈ L. Note that ( 7) is a continuous function where (7a) and (7b) correspond to negative log-likelihood values at maximum and zero probability respectively. We represent the infinite cost in (7b) by a large constant several orders of magnitude (10 9 -10 12 ) greater than the value of (7c) corresponding to the smallest estimated likelihood probability, so that, for practical settings of t > 0, (7a) and (7b) amount to imposing hard constraints for acceptable configurations. By comparison, [START_REF] Boykov | Graph cuts and efficient N-D image segmentation[END_REF] avoid multiplying (7b) by t in (5) and define a constant representing its cost such that prior terms do not accidentally impose a labeling decision on corresponding variables. In our formulation, as t greatly decreases, (7b) approaches the order of magnitude of prior terms until the data sum of (5) vanishes at t = 0, which is hardly an interesting setting in practice.

In subsequent refinements of an initial segmentation, we allow the user to introduce additional "corrective" seeds to constrain the label preference of incorrectly segmented image regions. In order not to disrupt initial appearance statistics, we do not take such seeds into account for likelihood estimation.

1 Recall that seeds are label attributions to voxels inside targeted structures.

C. Multilabel optimization via Graph Cut

Having defined the prior model and data-cost terms, we now turn our attention to the numerical optimization of the energy function of segmentation [START_REF] Blume | Google body: 3D human anatomy in the browser[END_REF]. We do not give the details of the optimization algorithm itself, but establish a necessary and sufficient condition of its applicability for optimizing the aforementioned energy function using our prior definition [START_REF] Rother | GrabCut": interactive foreground extraction using iterated graph cuts[END_REF].

In [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], it has been shown that the multilabel minimization of ( 5) is NP-hard even with the simple discontinuity-preserving Potts prior model. In binary labeling with |L| = 2 using this model, the global minimum of ( 5) can be computed in polynomial time by a maximum flow algorithm [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF]. In general, however, and in our multiobject segmentation problem with |L| > 2, efficient approximation algorithms must be resorted to. We use the Expansion Moves multilabel Graph Cut algorithm described in [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], which has been shown to outperform popular multilabel optimization algorithms on benchmarks in terms of both speed and quality of obtained solutions [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields with smoothness-based priors[END_REF].

The Expansion Moves algorithm produces a labeling locally-optimal within a known factor of the global minimum equal to 2 × max V i,j (l a , l b )/min V i,j (l a , l b ), ∀i, j ∈ S, l a , l b ∈ L, a = b. For example, this factor is equal to two for the Potts model. According to [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF], the algorithm can minimize (5) if the pairwise term is submodular, that is, if the following condition holds ∀i, j ∈ S, ∀l a , l b , l c ∈ L:

V i,j (l a , l a ) + V i,j (l b , l c ) ≤ V i,j (l a , l c ) + V i,j (l b , l a ). (8)
The vicinity prior ( 6) is submodular by definition. This can be easily verified by observing that the shortest-path distance between any three vertices on the graph A satisfies the triangle inequality. By comparison, the adjacency prior in [START_REF] Jacinto | A web interface for 3d visualization and interactive segmentation of medical images[END_REF] is essentially limited to definitions where the penalty of nonadjacency has to be twice the penalty of adjacency so that the prior is submodular. In other words, it is impossible to penalize nonadjacencies by some large constant without penalizing adjacencies half as severely.

IV. EVALUATION RESULTS AND APPLICATIONS

This section presents the results of qualitative validation and quantitative analysis of our vicinity prior-based Graph Cut segmentation approach. In all evaluations, comparisons are drawn with the Potts prior-based approach, and in some cases with our earlier adjacency prior [START_REF] Jacinto | A web interface for 3d visualization and interactive segmentation of medical images[END_REF]. We note that our software implementation uses the multilabel energy optimization library developed by the authors of [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF], [START_REF] Boykov | An experimental comparison of mincut/max-flow algorithms for energy minimization in vision[END_REF].

In all evaluations, we measure the quality of segmentation with respect to ground-truth via the well-known Dice similarity metric (DSM), which measures the overlap between the segmentation and the ground-truth. Let S l and T l represent the sets of voxels labeled with l ∈ L in the segmented I and the ground-truth T images respectively, and denote the DSM for label class l ∈ L by dsm l (I, T ). We define an "overall" segmentation quality metric, calculated from mean weighted DSM measures for all labels in L as:

dsm L (I, T ) = l∈L dsm l (I, T ) |T l | |T | . ( 9
) (a) (b) (c) (d) (e)
rel(d) = Pr(D > d) = |I d | |I| . ( 10
)
For the purpose of evaluating the precision of segmentation boundary placement on the medical case study in Section IV-B2, we use an error metric to measure the Euclidean distance of each point on the surface of a segmented structure from the surface of the corresponding ground-truth structure. Let M l S and M l T be respectively the triangle meshes representing isosurfaces of structure volumes labeled by l ∈ L in the segmented I and the ground-truth T image. We obtain such meshes by applying the Marching Cubes algorithm to labeled volumes [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF]. The error at x ∈ M l S is defined as:

err(x) = inf {d(x, y) | y ∈ M l T }. (11)

A. Synthetic images

We have generated a dataset from a synthetic phantom we have created, inspired by the Shepp-Logan phantom [START_REF] Shepp | The fourier reconstruction of a head section[END_REF], for the purposes of quantitative and qualitative evaluation of segmentation with vicinity prior. Due to the small number of structures and their simplistic spatial arrangement, the Shepp-Logan phantom is not well-suited for such evaluation. The layout of the 10 structures our phantom features (Fig. 3a) mimics the spatial configuration of some thoracic and abdominal organs in a coronal cross-section of a CT scan. In fact, the intensity levels of these structures were derived from mean intensity values of corresponding structures in a CT image.

Graph Cut segmentation with either prior was applied to the pixels2 as well as to CVT clusterings of 512 × 512 images generated from the phantom in Fig. 3a by the introduction of additive zero-mean Gaussian noise with standard deviation σ ∈ [0, 0.2], image intensity values being clamped to the range [0, 1]. 40 images were generated for each noise level. Of ranges swept by the hyperparameter t of ( 5), we report results in the representative window [0, 1]. We report clustering resolution as a percentage of the number of image pixels. The set of "seeds" used in likelihood estimation and to impose constraints is given in Fig. 4a. The number of seed-marked pixels for each label is in the range [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF]160] and roughly proportional to the size of the corresponding structure, otherwise constraints of equal support for all labels would create uneven bias of segmentation ease for smaller structures, such as "9" and "10".

Fig. 4 gives a qualitative comparison of best results produced by Potts and vicinity prior-based Graph Cut segmentation applied to the pixels and a 10% clustering of an image generated from the phantom with σ = 0.09, which can be considered a relatively high noise level. The vicinity prior was defined according to the graph model presented in Fig. 3b. Improvements brought by the vicinity prior are obvious on both segmentations, especially on thin structures "2" and "9".

In the remaining of this section as well as in Fig. 4, 5 and 6, any reference to DSM is actually to overall DSM measures for all structures as defined by [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF]. Mean overall DSM values are calculated on segmentations corresponding to best temperature settings for all 40 images generated at a given noise level.

In Fig. 5a we give mean DSM maps for segmentation with either prior applied to pixels. We calculate mean DSM values for segmentations corresponding to a temperature setting on all images generated at a given noise level. Improvements brought by the vicinity prior are remarkable, particularly for noise levels up to σ = 0.1. Note that low mean DSM values corresponding to low temperature settings are due to high regularization by the prior since data terms would not contribute as much at that level. Also recall that the vicinity prior incurs higher regularization than the Potts prior since its highest penalty in this case is four times greater than its Potts counterpart. To observe the difference between the two priors on a finer scale, we give a magnified view of DSM maps on the low temperature range [0.005, 0.05]. We do not give mean DSM maps for segmentation applied to clustered images which display similar performance patterns. ters with respect to increasing noise levels at the temperature setting yielding the best segmentation for each image. Fig. 5c and 6d compare reliability curves of segmentation with these priors on pixels and clusters respectively. Since evaluations with other clustering resolutions display similar behavior, we present the results of those corresponding to the level 10%. We also evaluate the robustness of Potts, adjacency and vicinity prior segmentation to clustering resolution. Due to limited space, results corresponding to a single noise level are given in Fig. 6b. We can see that vicinity prior segmentation is quite robust to clustering resolution. This means that similar segmentations can be obtained using lower clustering resolutions thus improving both runtime and memory footprint.

Overall, on evaluations with respect to varying noise and clustering resolution, we notice that the vicinity prior outperforms the Potts prior, and the adjacency prior to a lesser extent.

Lastly, we present the results of evaluating prior weighting schemes for varying clustering resolution. Results corresponding to clustering level 5% are given in Fig. 6a. Let C i and C j be a pair of clusters sharing a common boundary (a surface in 3D) with centroids c i and c j respectively. We compare the following four weighting schemes: 1) cluster common-boundary length weighting, where pairwise terms (6) are weighted by |∂C i ∩ ∂C j |, 2) inverse centroid-distance weighting, where pairwise terms are weighted by 1/d(c i , c j ), 3) both weighting schemes applied together, and 4) no weighting. Fig. 6a shows that the first weighting scheme is a clear winner.

B. Simulated and real medical images

1) BrainWeb simulated MRI dataset: We have performed a comprehensive quantitative evaluation and comparison on the BrainWeb simulated MRI dataset [START_REF] Collins | Design and construction of a realistic digital brain phantom[END_REF]. In the creation of this dataset, a semirealistic anatomical head and brain phantom was used by an MRI simulator of T1-, T2-and PD-weighted modalities to generate images according to MR acquisition physics with 6 noise levels (0%, 1%, 3%, 5%, 7% and 9%) and 3 intensity inhomogeneity levels (0%, 20% and 40%). For evaluation purposes, we have used 18 images from the T1-weighted modality covering the entire noise and inhomogeneity ranges. Image volumes are of size 217 × 181 × 217 voxels with isotropic 1 mm voxel size. The phantom itself used to generate images comprises 10 structure and tissue classes, of which we use 9 classes3 listed in Table I.

It should be noted that the phantom exhibits certain anatomical inconsistencies. For example, some structures corresponding to tissue classes "muscle" and "skin" are located inside the skull within the volume of tissue class "cerebrospinal fluid" and adjacent to "gray matter" structures. Therefore, we have defined our graph-based anatomical model as a compromise between phantom imprecisions and true anatomical properties.

In Table I, we give means and standard deviations for overall as well as individual-class DSM measures calculated on voxel and 10%-clustered image segmentations corresponding to best temperature settings for all images in the dataset. The vicinity prior is defined according to the graph model in Fig. 8a. It is easy to observe in Table I that Graph Cut segmentation attains better average performance levels with the vicinity prior on the BrainWeb dataset. It also performs slightly better on clusters. Improvements significant at level α = 0.05 are highlighted. Mean runtime and memory footprint figures for voxel and 10%-clustered image segmentations are given in Table III.

We have compared our results on the BrainWeb dataset with those reported in [START_REF] Gorthi | Comparison of energy minimization methods for 3-D brain tissue classification[END_REF] which is the only recent work we are aware of that evaluates the Expansion Moves algorithm on brain tissue classification using the Potts prior model. There are, however, some differences between the segmentation approach in [START_REF] Gorthi | Comparison of energy minimization methods for 3-D brain tissue classification[END_REF] and ours. Firstly, it uses Gaussian likelihood functions the parameters of which are learned during an iterative application of the algorithm. Secondly, it discards all but brain structures i.e. cerebrospinal fluid, gray and white matters, which exhibit good contrast in images in general [44]. Thus, computational complexity is reduced and the segmentation task is made easier. In our evaluation, we have used the same subset of the BrainWeb dataset as [START_REF] Gorthi | Comparison of energy minimization methods for 3-D brain tissue classification[END_REF] corresponding to the 12 T1-weighted images generated with inhomogeneity levels 20% and 40% and all noise levels. In Fig. 7, we compare means and standard deviations for DSM, sensitivity and specificity measures 4 for the 3 brain structures from full-image vicinity prior segmentations with corresponding measures in [START_REF] Gorthi | Comparison of energy minimization methods for 3-D brain tissue classification[END_REF]. Unsurprisingly, results we report are not superior. Specialized methods are expected to, and indeed should, produce better results. Therefore, the present comparison is an indication of how close a generic approach can come to a specialized one on this problem instance. Nevertheless, a practical advantage of our approach is that no manual tissue removal is required to segment the brain volume.

2) Thoracic-abdominal CT image: Lastly, we present a segmentation case study of in vivo acquired data on a thoracicabdominal 3D CT-scan realized during the arterial phase in inhaled position. The image comes from the 3D-IRCADb dataset [START_REF]IRCADb[END_REF] and is accompanied by a manual segmentation of 21 structures created by clinical experts of which we use 16 in evaluation. The volume we use is of size 480×370×167 voxels with 0.961 × 0.961 × 1.8 mm voxel size. We have segmented all visible structures in this image, including those for which no ground-truth estimate was available. The vicinity prior is defined according to the graph model of 34 structures given in Fig. 8b. Seeds were input by a user through a graphical interface by marking few dozens of pixels of any target structure visible in every tenth axial cross-section. Table II gives overall and individual-class DSM measures for Potts and vicinity prior-based Graph Cut segmentation applied to the 5%-clustered image at t = 0.5. The statistical test shows that the improvement brought by vicinity prior on mean DSM is significant at level α = 0.05. Qualitative comparison is given in Fig. 9 on coronal cross-sections of ground-truth, vicinity prior and Potts prior-segmented images. Note that the ground-truth can be relied on for qualitative inspection only for the structures listed in Table II. Fig. 10 gives 3D views of simplified surface meshes generated from different structure volumes via the mesh simplification method in [START_REF] Valette | Approximated centroidal voronoi diagrams for uniform polygonal mesh coarsening[END_REF]. Runtime and memory footprint figures for voxel and 5%clustered image segmentations are given in Table III.

We conclude this section with surface error measurements for Potts and vicinity prior-based segmentations of the 5%clustered thoracic-abdominal CT image. In Fig. 11, we trace the cumulative histograms of error (11) for 6 representative structures from both segmentations. graphs indicate the percentage of segmented surface lying within a given distance from the true surface. We can see that for reasonably low error levels of 0-1 cm, vicinity prior segmentation has retrieved a considerably larger fraction of structure boundaries than its Potts prior counterpart suffering much lower maximum error. However, it has produced some noticeable outliers, notably for the heart and right lung volumes. Such outliers constitute a small percentage of the structure, and can be corrected via "corrective" seeds as mentioned in Section III-B.

V. DISCUSSION AND FUTURE WORK

Unlike the Potts prior model, the vicinity prior model we propose is a piecewise-constant model incurring multiple levels of penalization capturing explicit spatial configuration of objects in multiobject segmentation problems. Compared to shape prior-based segmentation approaches, structural models capturing such spatial configurations are quite robust to shape deformations, because relative positions of objects remain largely stable. Furthermore, such models are relatively easy to define at an arbitrary level of detail out of simple specifications of spatial relationships on pairs of objects, and may also be generated from existing anatomical models for medical applications as we have seen in Section IV.

The shortest-path formulation produces a prior that is submodular which allows us to optimize (5) via the efficient Expansion Moves Graph Cut algorithm. Pairwise priors learned from training sets might not satisfy this property. Therefore some authors truncate nonsubmodular terms in order to optimize with Graph Cut algorithms [START_REF] Winn | The layout consistent random field for recognizing and segmenting partially occluded objects[END_REF], and others resort to arbitrary energy minimizers [START_REF] Lucchi | Are spatial and global constraints really necessary for segmentation?[END_REF].

Qualitative validation and quantitative analysis on synthetic, simulated and real medical images in Section IV confirmed the advantages of our vicinity prior model over the standard Potts model for Graph Cut segmentation, particularly for the segmentation of distinct objects having identical intensities, the accurate placement of structure boundaries in addition to robustness to clustering resolution. However, like most Graph Cut approaches, our segmentation method is prone to a shrinking bias. This is especially true for thin, elongated structures, like vasculature, where the sum of pairwise costs along the boundary of a vessel is higher than the cost of a shortcut through its thin structure, especially when data-cost terms do not incur a sufficiently strong preference for the vessel's class. We can notice the fragmentary segmentation of abdominal and hepatic arteries and portal veins in Fig. 10b. A number of solutions have been proposed to overcome this problem. For example, [START_REF] Vicente | Graph cut based image segmentation with connectivity priors[END_REF] improves an initial singleobject Graph Cut segmentation with connectivity constraints introduced by the user via markings on image regions which need to be connected to the principal object. [START_REF] Pamulapati | Intra-hepatic vessel segmentation and classification in multi-phase CT using optimized graph cuts[END_REF] proposes a Hessian-based "vesselness" likelihood approach for Graph Cut segmentation of hepatic blood vessels and portal veins in images, which could be applied to similar structures such as bronchi and pulmonary arteries.

We are currently investigating the extension of the graphbased prior model to express other types of relations, notably spatial directionality. Since the shortest-path weight formulation on a directed graph would produce asymmetric pairwise prior definitions, such extensions would understandably require revisiting some of the fundamental aspects of the underlying theory of Bayesian labeling.

Another interesting venue for future research is the introduction of uncertainty into the graph-based prior model. As noted in Section III-A, graph models we use to define the vicinity prior are not limited to unit or discrete edge weights. Continuous weighs might be used, for example, to express degrees of adjacency, which might be learned from several model variants and thus capture possible variabilities. His teaching and research interests include image processing, pattern recognition, shape modeling and machine learning.

Fig. 1 .

 1 Fig. 1. A CVT clustering and its dual graph for a circle image.

Fig. 2 .

 2 Fig. 2. Left column: quantitative evaluation of clustering algorithms on the Berkeley segmentation dataset. Right column: qualitative comparison of their output with 300 clusters on a single photo from the dataset. Boundary recall and undersegmentation error are respectively (b) CVT: 0.59, 1.53 (d) TurboPixels: 0.52, 2.01 (f) GC superpixels: 0.59, 2.76.

Fig. 3 .

 3 Fig. 3. (a) Synthetic phantom and (b) the graph model for adjacency relations of its structures.

  R be the set of symmetric binary relations on pairs of distinct labels, R = {r | l a r l b , l a , l b ∈ L, a = b}, representing such relations as adjacency, distance, area of common surface, etc. Assume that each relation r ∈ R has a corresponding characteristic function r which gives a quantitative measure for r; r : L × L \ {(l a , l a )} → R + . R can be equivalently represented by a weighted undirected graph A = L, W on labels L with the set of edges W = {l a , l b } | ∃r ∈ R, l a r l b , a = b where edge weights are defined by w ({l a , l b }) = r (l a , l b ), such that w ({l a , l b }) = ∞, ∃r ∈ R, l a r l b .

Fig. 4 .

 4 Fig. 4. (a) Image generated from the phantom in Fig. 3a with additive Gaussian noise σ = 0.09 along with "seeds" used in all evaluations (b) Potts prior (DSM = 0.937) and (c) vicinity prior (DSM = 0.994) segmentations on pixels of the image (d) Potts prior (DSM = 0.938) and (e) vicinity prior (DSM = 0.986) segmentations on a 10% clustering of the image. Results correspond to best temperature settings. Refer to color code in Fig. 3b

Fig. 5 .Fig. 6 .

 56 Fig. 5. Quantitative evaluation and comparison of Potts, adjacency and vicinity prior-based pixel segmentation on the synthetic phantom dataset.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Quantitative comparison of brain tissue segmentation with [31] on the BrainWeb simulated MRI dataset.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Qualitative comparison of Potts and vicinity-prior based segmentations of the 5%-clustered 3D-IRCADb thoracic-abdominal CT image. Coronal cross-sections are for (a) ground-truth (b) vicinity prior and (c) Potts prior-segmented images. Refer to color code in Fig. 8b. Ellipses on (c) indicate erroneous segmentations of right lung (top left), muscle (top right), air (center left), right and left renal pyramids (center) and bone marrow (bottom). Notice how right and left lungs are correctly segmented as distinct organs with vicinity prior.

  Surface fraction w.r.

Fig. 11 .

 11 Fig. 11. Surface error measurements for (a) Potts and (b) vicinity prior segmentations of the 5%-clustered 3D-IRCADb thoracic-abdominal CT image.

Razmig

  Kéchichian received his BSc and MSc degrees in computer science from the University of Aleppo, Syria and the National Institute of Applied Sciences (INSA) of Lyon, France, respectively in 2007 and 2009. He is currently a PhD candidate in the Images and Models project in the Creatis Laboratory (CNRS UMR 5220, Inserm U1044) at INSA-Lyon. His research interests include biomedical image processing, MRF modeling and combinatorial optimization. Sébastien Valette received his MSc degree from the Electrical Engineering Department at the National Institute for Applied Sciences (INSA) of Lyon, France, in 1998. He obtained the PhD Degree at INSA-Lyon in 2002. He is currently a CNRS researcher in the Images and Models project in the Creatis Laboratory (CNRS UMR 5220, Inserm U1044) at INSA-Lyon. His research interests include geometry and image processing. Michel Desvignes received his MSc degree from the Enginering School of Caen (ISMRA) in 1985 and the PhD degree from the University of Caen in 1990. He is currently a computer science professor at Grenoble Institute of Technology (Grenoble-INP).

TABLE I MEAN

 I AND STANDARD DEVIATIONS FOR OVERALL AND INDIVIDUAL-CLASS DSM MEASURES FOR VOXEL AND 10%-CLUSTERED SEGMENTATIONS OF BRAINWEB SIMULATED MRI DATASET WITH POTTS AND VICINITY PRIORS. HIGHLIGHTED VALUES CORRESPOND TO IMPROVEMENTS SIGNIFICANT AT α = 0.05. ± 0.01 0.98 ± 0.01 0.98 ± 0.01 0.99 ± 0.00 Cerebrospinal fluid 0.25 ± 0.16 0.65 ± 0.06 0.29 ± 0.11 0.67 ± 0.05 Gray matter 0.79 ± 0.09 0.82 ± 0.08 0.80 ± 0.08 0.85 ± 0.06 White matter 0.86 ± 0.08 0.88 ± 0.07 0.87 ± 0.07 0.88 ± 0.06 Fat 0.89 ± 0.05 0.84 ± 0.08 0.87 ± 0.04 0.85 ± 0.06 Muscle 0.73 ± 0.06 0.74 ± 0.12 0.73 ± 0.07 0.79 ± 0.09 Skin 0.70 ± 0.07 0.75 ± 0.09 0.69 ± 0.07 0.75 ± 0.07 Skull 0.76 ± 0.07 0.79 ± 0.05 0.76 ± 0.05 0.80 ± 0.04 Connective tissue 0.71 ± 0.14 0.69 ± 0.21 0.73 ± 0.12 0.75 ± 0.15 Overall 0.81 ± 0.03 0.85 ± 0.04 0.82 ± 0.03 0.87 ± 0.04

		voxels		clusters 10%
	Structures	Potts	vicinity	Potts	vicinity
	Air	0.98

REFERENCES

[1] S. Seifert, A. Barbu, S. K. Zhou, D. Liu, J. Feulner, M. Huber, M. Suehling, A. Cavallaro, and D. Comaniciu, "Hierarchical parsing and semantic navigation of full body CT data," in SPIE Medical Imaging, Lake Buena Vista, FL, USA, Feb. 2009.

TABLE II OVERALL

 II AND INDIVIDUAL-CLASS DSM MEASURES FOR 5%-CLUSTERED SEGMENTATION OF 3D-IRCADB THORACIC-ABDOMINAL CT IMAGE WITH POTTS AND VICINITY PRIORS. HIGHLIGHTED VALUES CORRESPOND TO IMPROVEMENTS SIGNIFICANT AT α = 0.05.TABLE III RUNTIME OF VICINITY PRIOR SEGMENTATION ON VOXELS AND 10% AND 5%-CLUSTERED BRAINWEB MR AND 3D-IRCADB CT IMAGES RESPECTIVELY, RECORDED RESPECTIVELY ON A 4-CORE CPU @ 2.84 GHZ WITH 6 GB OF RAM AND A 16-CORE CPU @ 2.66 GHZ WITH 47 GB OF RAM. LAST 4 COLUMNS GIVE CLUSTERING, GRAPH CONSTRUCTION AND OPTIMIZATION TIMES IN MINUTES AND MEMORY FOOTPRINT IN MB.

	Structures	Potts vicinity Structures		Potts vicinity
	Bone	0.735	0.857 Kidney (R)		0.907	0.967
	Heart	0.918	0.988 Kidney (L)		0.880	0.968
	Lung (R)	0.096	0.976 Adrenal gland (R) 0.765	0.743
	Lung (L)	0.646	0.980 Adrenal gland (L) 0.706	0.707
	Liver	0.990	0.989 Hyperplasia		0.645	0.646
	Spleen	0.965	0.982 Aorta		0.512	0.950
	Gallbladder 0.902	0.923 Vena Cava		0.894	0.923
	Pancreas	0.895	0.959 Portal vein		0.537	0.800
	Overall	0.721	0.962 mean		0.750	0.897
	Image		Segmentation Clust. Constr. Optim.	Mem.
	BrainWeb MRI	voxels clusters 10%	-3.19	0.25 0.18	12.46 2.83	3053.13 376.55
	3D-IRCADb CT	voxels clusters 5%	-15.48	2.10 460.23 23864.38 1.44 14.57 2624.69

have performed some of the evaluations on image pixels (voxels) for comparison purposes with clustered image segmentation. The theoretical framework outlined in Section III still applies by letting variables i ∈ S correspond to image pixels (voxels) v ∈ I and defining a 4-connected (6connected) neighborhood N i for each variable (respectively in

3D).

We have merged the tissue class "glial matter" into the class "gray matter" without affecting anatomical consistency.

Sensitivity and specificity are defined respectively as the true-positive and the true-negative rates for individual label classes.