Dennis Van Der Meer
email: dennis.vandermeer@minesparis.psl.eu

Simon Camal
email: simon.camal@minesparis.psl.eu

Georges Kariniotakis
email: georges.kariniotakis@minesparis.psl.eu

Generalizing Renewable Energy Forecasting Using Automatic Feature Selection and Combination

important to select important features and discard redundant ones in order to retain accuracy and reduce time complexity. Often, feature selection is classified into three categories, namely: (i) filtering methods; (ii) wrapper methods; and (iii) embedded methods [START_REF] Guyon | An introduction to variable and feature selection[END_REF]. Filtering methods rank features based on a score, e.g., Pearson correlation, to subset the features, and are therefore model agnostic. Recently, an ultra-fast algorithm for similarity search was proposed to preselect relevant features [START_REF] Yang | Ultra-fast preselection in lasso-type spatio-temporal solar forecasting problems[END_REF]. Wrapper methods aim to find the optimal feature subset by subsetting the feature set, learning a predictor on this subset and comparing its performance to the same predictor on other feature subsets. Wrapper methods are not common due to their computational cost although a brute-force attempt to find a relevant subset of endogenous variables was made in [START_REF] Van Der Meer | Probabilistic forecasting of electricity consumption, photovoltaic power generation and net demand of an individual building using gaussian processes[END_REF]. Finally, the embedded methods embed the feature selection method, such as Lasso regression [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. A multitude of studies in solar forecasting have applied embedded methods, e.g, Lasso combined with quantile regression [START_REF] Agoua | Probabilistic Model for Spatio-Temporal Photovoltaic Power Forecasting[END_REF].

As the penetration of RESs increases, it becomes more pertinent to aggregate these assets as VPPs in order to participatewith high reliability-in market bidding and power system regulation. Large portfolio end-users often acquire multiple forecasts originating from different service providers or different NWPs etc. In this paper we combine feature selection methods with forecast combination techniques to derive a model that outperforms the state of the art. It is well-known that forecast combination can reduce modeling uncertainty although it is not commonly researched in solar forecasting [START_REF] Yang | Post-processing in solar forecasting: Ten overarching thinking tools[END_REF]. Probabilistic forecasts can be combined using heuristic, linear and nonlinear combination methods; see [START_REF] Yang | Post-processing in solar forecasting: Ten overarching thinking tools[END_REF] for an overview. For instance, the opinion linear pool (OLP) is a heuristic method that linearly combines predictive distributions by assigning equal weights to each expert. In solar forecasting, minimization of the continuous ranked probability score (CRPS) has been used to learn the weights to linearly combine predictive distributions in a traditional linear pool (TLP) [START_REF] Bracale | A Probabilistic Competitive Ensemble Method for Short-Term Photovoltaic Power Forecasting[END_REF]. Thorey et al. [START_REF] Thorey | Ensemble forecast of photovoltaic power with online CRPS learning[END_REF] developed an online optimization of the CRPS to linearly combine forecasts, which improved the calibration although underdispersion remained. Besides TLP, nonlinear combination methods have been proposed, i.e., the spread-adjust linear pool (SLP) and the beta-transformed linear pool (BLP) [START_REF] Gneiting | Combining predictive distributions[END_REF]. Möller and Groß used SLP to calibrate the ensemble prediction system of the European Center for Medium-range Weather Forecasts (ECMWF) [START_REF] Möller | Probabilistic temperature forecasting with a heteroscedastic autoregressive ensemble postprocessing model[END_REF], while BLP Abstract-Spatially aggregating renewable power plants is beneficial w hen p articipating i n e lectricity m arkets. I n t his context, a substantial number of features is available from various data sources. In machine learning, feature selection is common so as to relieve the curse of dimensionality and avoid overfitting. However, there is no guarantee that the selected features result in reliable forecasts and post-processing can therefore be valuable. In this study, we combine model agnostic feature selection with linear and nonlinear probabilistic forecast combination techniques. Moreover, the filters a utomatically c ompute t he w eights f or our analog ensemble (AnEn) forecast model. We verify our model chain by generating intra-day forecasts of the aggregated output of 20 photovoltaic power plants using 831 input features in total. We show that the collection of filters s elects a heterogeneous feature set but that each individual AnEn-filter combination results in underdispersed forecasts, which is efficiently remedied by the forecast combination techniques.

Index Terms-Filtering, forecast combination, highdimensional data, virtual power plant, probabilistic forecasts

I. INTRODUCTION

The power generated by renewable energy sources (RESs) such as photovoltaic (PV) systems or wind turbines is highly variable and it is well-established that probabilistic forecasts contribute to optimising power system management when used in functions like unit commitment [START_REF] Wang | Wind power forecasting uncertainty and unit commitment[END_REF]. Unit commitment is usually performed on the day-ahead market (DAM) and real-time market (RTM), where penalties are incurred based on the forecast errors. Consequently, the accuracy of the probabilistic forecasts is important, which is challenged by the aforementioned variability. The variability can be relieved by aggregating spatially distributed power plants. The resulting virtual power plant (VPP) outputs a smoother power profile caused by the reduced spatial correlation. However, forecasting the VPP output requires more input features as the VPP likely covers large geographical areas reflected b y multiple grid points of Numerical Weather Prediction (NWP) models or satellite images.

There are at least two challenges that arise when including vastly more input features: (i) features may be multicollinear; and (ii) the curse of dimensionality. In other words, it is

The present research was carried as part of the Smart4RES Project (European Union's Horizon 2020, No. 864337). was used to combine two parametric predictive densities from regression models based solely on lagged observations [START_REF] Fatemi | Parametric methods for probabilistic forecasting of solar irradiance[END_REF]. Consequently, there is need for additional study into nonlinear combination methods applied to solar forecasting.

In this study, we focus on operational forecasts of the intra-day aggregated output of a VPP that consists of 20 PV plants with a combined nominal capacity of 4005 kW, located in mid-west France. The total number of features is 831: 289 satellite pixels and 108 NWP grid points with 5 variables and 2 astronomical variables (cf Fig. 1). We generate the probabilistic forecasts with an Analog Ensemble (AnEn), which is an established forecasting method that compares the current weather state to historical weather states and outputs the observed power generation corresponding to the most similar historical states as an ensemble forecast [START_REF] Alessandrini | An analog ensemble for short-term probabilistic solar power forecast[END_REF]. We use 6 filtering methods based on mutual information to select the input features for our AnEn model and automatically compute their weights. Since a recent benchmark study of 22 filter methods concludes that "there is no subset of filter methods that outperforms all other filter methods" [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF], we combine the predictive distributions using 1 linear and 2 nonlinear combination methods. The contributions of this study can therefore be summarized as follows: (i) we employ 6 filtering methods based on mutual information to account for the fact that no filter works best in all circumstances [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF] and to account for the nonlinearity between PV power generation and the features; (ii) we apply linear and nonlinear probabilistic combination methods to minimize model uncertainty, which is an underserved area of research in solar forecasting [START_REF] Yang | Post-processing in solar forecasting: Ten overarching thinking tools[END_REF]; and (iii) these methods combined generalize existing probabilistic forecasting approaches by automating filtering and model combination.

The remainder of this paper is organized as follows. Section II details the methods used in this study. Then, Section III presents the results and Section IV concludes this paper.

II. METHODOLOGY

This section describes the methodology of the paper and starts with describing the filtering methods and the forecast combination methods. Then, the forecast evaluation tools used in this study are described. Figure 2 presents an overview of the methodology.

A. Filtering

Filtering is a category of feature selection in which a score is computed to rank the importance of features. The goal is to find a feature subset S from the complete feature set X that represents the target variable Y with high accuracy and minimal residual uncertainty [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF]. In this study, the filtering methods are based on mutual information, which originates from information theory and describes the amount of information that can be known about random variable Y when knowing random variable X [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF].

For the remainder, it is helpful to introduce the entropy of a discretized random variable Y with marginal probability mass function p [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF]:

H(Y ) = - y p(y)log 2 (p(y)) , (1)
which measures the uncertainty of Y . In addition, the conditional entropy of Y given X is defined as [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF]:

H(Y |X) = x p(x) - y p(y|x)log 2 (p(y|x)) . (2)
The mutual information between Y and X is then defined as I(Y ; X) = H(Y ) -H(Y |X) and describes how the uncertainty of Y is lowered by knowing X [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF]. In this study, we use the praznik package [START_REF] Kursa | praznik: Tools for Information-Based Feature Selection[END_REF] in the statistical software R. The praznik toolbox discretizes the range of the continuous features into max{min{ n 3 , 10}, 2} intervals that are equally spaced, where n is the number of training samples [START_REF] Kursa | praznik: Tools for Information-Based Feature Selection[END_REF].

As mentioned above, we use 6 different filtering methods based on mutual information. The first filter, mutual information maximization (MIM), computes the mutual information between feature k and target variable Y as:

J MIM (X (k) ) = I(Y ; X (k) ), (3)
and returns a predetermined number of features that maximizes J [START_REF] Kursa | praznik: Tools for Information-Based Feature Selection[END_REF]. The remaining filters greedily add a feature to S that at each iteration maximizes the score J(X (k) ). The first of those is minimal conditional mutual information (CMIM), which uses the score

J CMIM (X (k) ) = min I(Y ; X (k) ), min X (j) ∈S I(Y ; X (k) |X (j) ) , (4)
where [START_REF] Kursa | praznik: Tools for Information-Based Feature Selection[END_REF]. Similar to MIM, CMIM describes how knowing X (k) lowers the uncertainty of Y given that we already know X (j) . The third filter is conditional mutual information (CMI), which ranks features according to the score:

I(Y ; X (k) |X (j) ) = H(Y |X (j) ) -H(Y |X (k) , X (j) )
J CMI (X (k) ) = I(Y ; X (k) |S), (5)
and essentially evaluates the added value of feature X (k) given the already selected features. The fourth filter is double input symmetrical relevance (DISR), which uses the score [START_REF] Kursa | praznik: Tools for Information-Based Feature Selection[END_REF]:

J DISR (X (k) ) = X (j) ∈S I(Y ; X (k) , X (j) ) H(Y, X (k) , X (j) ) , (6)
where I(Y ; X (k) , X (j) ) evaluates the complementary information that X (k) and X (j) provide for Y , whereas H(Y, X (k) , X (j) ) reduces the possibility to choose highly variable features [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF]. Maximum relevance and minimum redundancy (MRMR) is the fifth filter and uses the following score to rank the features:

J MRMR (X (k) ) = I(Y ; X (k) ) - 1 |S| X (j) ∈S I(X (k) ; X (j) ), (7)
which is designed to ensure that there is as little redundancy as possible between X (k) and S (second term) and providing maximal information about Y [START_REF] Kursa | praznik: Tools for Information-Based Feature Selection[END_REF]. The sixth and final filter is minimal normalized joint mutual information (NJMIM), which scores the features according to [START_REF] Kursa | praznik: Tools for Information-Based Feature Selection[END_REF]:

J NJMIM (X (k) ) = min X (j) ∈S I(Y ; X (k) , X (j) ) H(Y, X (k) , X (j) ) , (8)
which is a modification of filter DISR and evaluates the minimal relative information between Y , X (k) and already selected features instead of the sum [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF].

Our algorithm runs over testing time instances t = n + 1, • • • , n ′ . The data that the filters are applied to comprise the d previous days recorded during the same time stamp as t (expressed as "HH:MM"). This setup allows us to include the most recent data in our filtering method. In order to account for temporal dependencies, we also include one time instant prior to "HH:MM" and one time instant after "HH:MM" (except for t + 1). Consequently, the filters have a total of 3 • d -1 time instances at their disposal.

These filters have been selected because (i) they can uncover nonlinear relationships; (ii) do not require the random variables to be on the same scale; (iii) performed decently as [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF] concludes; (iv) tend to select a different order of features [START_REF] Bommert | Benchmark for filter methods for feature selection in high-dimensional classification data[END_REF], which could benefit forecast diversity; and (v) except for MIM, take into account interactions between features.

B. Forecast Combination

As mentioned above, OLP is a heuristic method in which m predictive distributions F j are linearly combined as G = m j=1 w j F j where w j ≥ 0 and m j=1 w j = 1 with equal weights w j = 1/m. In this work, however, we determine the weights w j based on a scoring rule; specifically, the logarithmic score (IGN) as proposed by Gneiting and Ranjan [START_REF] Gneiting | Combining predictive distributions[END_REF]. The following subsections describe the three combination methods in detail. In addition, OLP is used as a reference method to compare the combination methods below with.

1) Linear Combination: The traditional linear pool (TLP) is similar to the opinion pool, except that we determine the weights by optimizing the logarithmic score. The weights can be determined by evaluating g i (y) = m j=1 w j f ij (y) and minimizing over n fitting samples:

IGN = - 1 n n i=1 log   m j=1 w j f ij (y i )   , (9)
which is equivalent to maximizing the log-likelihood [START_REF] Yang | Post-processing in solar forecasting: Ten overarching thinking tools[END_REF]. The main challenge of TLP is that it increases dispersion (cf. Theorem 3.1 in [START_REF] Gneiting | Combining predictive distributions[END_REF]) and is therefore not suited when the forecast members are probabilistically calibrated or overdispersed.

2) Nonlinear combination: In this study, we consider two nonlinear combination methods. The first is the spreadadjusted linear pool (SLP) that includes a strictly positive parameter c which adjusts the spread of the predictive distributions. The combined predictive distribution is defined as follows [START_REF] Gneiting | Combining predictive distributions[END_REF]:

G c i (y) = m j=1 w j F 0 ij y -µ ij c , (10)
where

F ij = F 0 ij (y -µ ij
) and µ ij is the median of predictive distribution of F ij . Equation 10shows that, in case c = 1, SLP reduces to TLP. Setting c < 1 tends to improve calibration when the individual predictive distributions are overdispersed or neutrally dispersed and setting c ≥ 1 may benefit underdispersed predictive distributions [START_REF] Gneiting | Combining predictive distributions[END_REF]. The flexibility of SLP is limited but it can be effective when the component forecasts are underdispersed or neutrally dispersed [START_REF] Gneiting | Combining predictive distributions[END_REF].

Therefore, we additionally consider the beta-transformed linear pool. Generalized by Gneiting and Ranjan [START_REF] Gneiting | Combining predictive distributions[END_REF], the predictive cumulative distribution function (CDF) of the betatransformed linear pool is defined as

G α,β i (y) = B α,β   m j=1 w j F ij (y)   , (11)
where B α,β represents the beta CDF with parameters α > 0 and β > 0, and BLP reduces to TLP when α = β = 1. Gneiting and Ranjan [START_REF] Gneiting | Combining predictive distributions[END_REF] prove that when the weights w 1 , . . . , w m are fixed, the variance of the probability integral transform random variable can attain any value in the open interval 0, 1 4 , with 1 12 indicating neutrally dispersed probabilistic forecasts.

C. Analog Ensemble

The work in [START_REF] Alessandrini | An analog ensemble for short-term probabilistic solar power forecast[END_REF] is one of the first studies that applied AnEn to solar forecasting. Therein, the authors compare the most recent NWP forecast to historical NWP forecasts and select the most similar NWP historical forecasts (the analogs), of which the corresponding PV power measurements form the ensemble forecast. In this study, we modify the original similarity metric to include additional data sources, similar to [START_REF] Carriere | A novel approach for seamless probabilistic photovoltaic power forecasting covering multiple time frames[END_REF]. Furthermore, we leverage the algorithmic efficiency of k-dimensional tree (kd-tree), as suggested by [START_REF] Yang | Ultra-fast analog ensemble using kd-tree[END_REF]. For this, we further modify the original similarity metric such that the lags of the features are distinct features [START_REF] Yang | Post-processing in solar forecasting: Ten overarching thinking tools[END_REF]:

d(X h t , A h i ) = J j=1 w h j x (j) t -x (j) i 2 , (12)
where J = N v • (2 t + 1) is the total dimension, w h j is the j th weight for the h th forecast horizon and x (j) is the j th feature. Vector X h t contains the most recent filtered features, which comprises the NWP forecasts and clear-sky data (CLS) valid at time t, as well as the satellite derived clear-sky index and PV power measurement observed h time steps before. Vector A h i contains the same filtered features as X h t except historical ones. Note that the features in A h i are scaled and centered and these scaling factors are used to scale and center X h t . Herein, the weights w h j in eq. 12 are taken directly from the scores that the filtering methods assign to the features. The weights are subsequently normalized by the sum of all feature scores so that the weights sum to 1. This allows us to dynamically update the feature weights instead of employing a costly wrapper.

D. Complete History Persistence Ensemble

We employ the complete history persistence ensemble (CH-PeEn) to evaluate the forecast methods described above [START_REF] Yang | A universal benchmarking method for probabilistic solar irradiance forecasting[END_REF]. CH-PeEn is an external multiple-valued climatology, is independent of forecast horizon and its output is calibrated by design [START_REF] Yang | A universal benchmarking method for probabilistic solar irradiance forecasting[END_REF]. Note that the CH-PeEn is constructed using detrended PV power measurements, as Section II-F details.

E. Forecast Verification

An important aspect of forecasts is their quality, which in the case of probabilistic forecasts consists of calibration, resolution and sharpness. To evaluate the probabilistic calibration of our forecasts, we compute the probability integral transform (PIT) z t = F t (y t ) and plot these as a histogram for visual verification [START_REF] Gneiting | Combining predictive distributions[END_REF]. In case of probabilistic calibration, the resulting histogram is close to flat and has a variance close to 1/12.

In addition, we evaluate the forecasts numerically and normalize the results using the nominal capacity of the VPP. The first score is the CRPS, which is a strictly proper and a negatively-oriented score. Interestingly, it is possible to decompose the average CRPS into reliability, uncertainty and resolution terms as CRPS = Reliability + Uncertainty -Resolution [START_REF] Hersbach | Decomposition of the continuous ranked probability score for ensemble prediction systems[END_REF]. Additionally, we use the CRPS to compute a skill score, i.e., the relative improvement over CH-PeEn, that has a maximum value of 100%. Finally, we use the square root of the average of the variance (RMV) of the predictive distributions as a measure of the sharpness [START_REF] Gneiting | Combining predictive distributions[END_REF].

Although this study focuses on probabilistic forecasts, the root mean squared error (RMSE) is arguably one of the most common verification metrics and is therefore included in our analysis.

F. Data

The available PV power measurements extend from 2019-01-01 until 2020-09-30 and measurements recorded at a zenith angle larger than 85°are removed. The period from 2019-01-01 until 2019-09-30 is used for training, while 2019-10-01 until 2019-12-31 is used to optimize the combination weights. The remainder of the data is used for testing. The PV plants have the same orientation and tilt, therefore we compute the clear-sky global tilted irradiance (GTI cs ) at the averaged coordinates. We detrend the PV power measurements by dividing them with GTI cs and note that more elegant solutions exist but they require additional system information that is unavailable to us. We include the last PV power observation as a feature (referred to as "VPP"), as well as the zenith angle and solar azimuth as features (referred to as "CLS").

The NWP forecasts come from the ECMWF, issued daily at 00:00 UTC with a spatial resolution of 0.1°× 0.1°. The 5 features that we extract are surface solar radiation downwards (SSRD), 100 m U-and V-wind speed (U100 and V100), 2 m temperature (2T) and total cloud cover (TCC). The SSRD is converted to global horizontal irradiance (GHI) and detrended to the clear-sky index using the McClear clear-sky model [START_REF] Lefèvre | McClear: a new model estimating downwelling solar radiation at ground level in clear-sky conditions[END_REF].

Finally, we include satellite derived GHI from Meteosat Second Generation (MSG) satellite images (referred to as "SAT"). The GHI is derived using an improvement over the Heliosat-2 method [START_REF] Zarzalejo | A new statistical approach for deriving global solar radiation from satellite images[END_REF] and stored as a GHI time series for each pixel. Similar to NWP GHI forecasts, the satellite derived GHI is converted to the clear-sky index. Figure 1 indicates where the satellite and NWP grid points and PV systems are located.

III. RESULTS

The first contribution of this study pertains to feature selection by means of the filters described in Section II-A. Figure 3 presents the normalized weights assigned by each filter to the feature groups as a function of the forecast horizon and for all forecast issue times. As can be seen, there are some notable differences between the filters. For instance, filter CMI includes feature group CLS even though the input and output features have been detrended. Note also that filters DISR and MIM to a large extent ignore the most recent power measurement. This can be expected of MIM since it does not consider interaction between input features and is therefore likely to select many similar features if they are relevant to Y . Since DISR evaluates the complementary information between input features, it does not specifically promote feature diversity. As mentioned above, filter NJMIM is a modification of DISR by evaluating the minimal relative information between features, which likely helps it to select more diverse features. As a final note, CMI and MRMR select diverse features as they are designed to reduce redundancy.

The second contribution concerning forecast combination aims at improving the calibration of the component forecasts. Table I presents the combination weights and SLP and BLP parameters, expressed as "mean ± standard deviation" over all forecast horizons. The differences between the combination methods are not drastic although it is interesting to note that BLP favors MIM and DISR, whereas TLP and SLP assign less importance to these filters. Interestingly, SLP parameter c is lower than 1, which may suggest overdispersed component forecasts [START_REF] Gneiting | Combining predictive distributions[END_REF]. In this case, we reason that the combination introduces overdispersion and parameter c tries to compensate.

Figure 4 presents the PIT histograms of the AnEn forecasts generated with the inputs provided by the filters and their combinations. The value in each subfigure corresponds to the variance, which should be close to 1/12 (≈ 0.083) for neutral dispersion [START_REF] Gneiting | Combining predictive distributions[END_REF]. The figure shows that the predictive distributions are calibrated in the main distribution but deviate at the extreme quantiles, which could be attributed to the selection of too few or redundant features. Linear forecast combination increases the dispersion [START_REF] Gneiting | Combining predictive distributions[END_REF], and Fig. 4 confirms this; OLP and TLP are both much closer to uniformity and their variance is closer to 1/12. In our setting, SLP combines the component forecasts most effectively, i.e., results in the most desirable calibration. Interestingly, BLP introduces negative bias, which could be due to the fact that the component forecasts are underdispersed, due to limited training instances or due to a discrepancy between the validation and verification sets.

The numerical results are presented in Fig. 5. It can be seen that the combination methods, especially OLP, TLP and SLP, improve upon the component forecasts across all scores and forecast horizons. Interesting to note is that OLP tends to perform only slightly worse when compared to TLP and SLP, which indicates that linear forecast combination can be highly effective with minimal computational effort as long as the component forecasts are underdispersed. Furthermore, it can be seen that the forecasts resulting from the filters are only slightly less sharp than the combination methods, which indicates that the PIT histogram shape is caused by the lowest (highest) quantile frequently being above (below) the target.

IV. CONCLUDING DISCUSSION

In this study, we generated intra-day forecasts of the aggregated output of a virual power plant with a nominal capacity of 4005 kW. We applied the analog ensemble (AnEn) as a base model as it requires no training. Given the large number of inputs comprising numerical weather prediction, satellite, endogenous and astronomical data, we proposed to use a variety of filters to (i) reduce the problem dimensionality and (ii) automatically and dynamically learn the AnEn feature weights. The results showed that the collection of filters selects a rather heterogeneous feature set but that each individual AnEn-filter combination suffers from calibration issues.

We employed linear and nonlinear forecast combination to improve the component forecasts, since that is an underserved area of research in solar forecasting. We showed that even a naive linear combination significantly improves the calibration, which is supported by theoretical results [START_REF] Gneiting | Combining predictive distributions[END_REF]. In our case, the spread-adjusted linear pool proved to be most effective at calibrating component forecasts as well as maximizing sharpness. Regardless, we recommend forecast practitioners to at least include naive linear forecast combination. An interesting direction for future research on feature selection would be to augment the filters with a time-dependent factor to weigh recent data more. Regarding forecast combination, interesting directions for future research include (i) combining a wider variety of forecast models and (ii) developing methods to compute the combination weights online.

Fig. 1 .

 1 Fig. 1. Overview of the satellite and NWP grid points and the PV systems.

Fig. 2 .

 2 Fig. 2. A flowchart of the pre-process, forecast, and post-process steps.

Fig. 3 .

 3 Fig.3. Normalized scores of the filters as a function of forecast horizon, used as dynamic weights in the AnEn model. Recall that VPP represents PV power at t -h, SAT satellite derived clear-sky index at t -h and CLS the zenith angle and solar azimuth at t.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Histograms of the PIT variables. The value in the subfigure corresponds to the variance, which should be close to 1/12 for neutral dispersion.

TABLE I THE

 I COMBINATION WEIGHTS AND SLP AND BLP PARAMETERS EXPRESSED AS "MEAN ± STANDARD DEVIATION" OVER ALL FORECAST HORIZONS. ± 0.065 0.206 ± 0.099 0.087 ± 0.065 0.068 ± 0.078 0.253 ± 0.057 0.109 ± 0.046 0.931 ± 0.028

				CMI				CMIM					DISR					MIM				MRMR
				0.098				0.099					0.102					0.097					0.101
			NJMIM				OLP					TLP					SLP					BLP
				0.102				0.074					0.076					0.079					0.075
		0.00	0.25	0.50	0.75	1.00	0.00	0.25	0.50	0.75	1.00	0.00	0.25	0.50	0.75	1.00	0.00	0.25	0.50	0.75	1.00	0.00	0.25	0.50	0.75	1.00
	Model		CMI		CMIM			DISR			MIM		MRMR			NJMIM		c			α			β
	TLP	0.267 ± 0.07	0.214 ± 0.097 0.076 ± 0.056	0.06 ± 0.084		0.274 ± 0.06	0.108 ± 0.033		-			-			-
	SLP	0.277 -			-
	BLP	0.227 ± 0.032 0.206 ± 0.066		0.11 ± 0.063	0.131 ± 0.069		0.224 ± 0.05	0.102 ± 0.026		-			0.855 ± 0.039 1.388 ± 0.069


David Coeurjolly
email: david.coeurjolly@liris.cnrs.fr

Jacques-Olivier Lachaud
email: jacques-olivier.lachaud@univ-smb.fr

A Simple Discrete Calculus for Digital Surfaces ⋆

Keywords: Discrete calculus, differential operators, digital surface

Computing differential quantities or solving partial derivative equations on discrete surfaces is at the core of many geometry processing and simulation tasks. For digital surfaces in Z 3 (boundary of voxels), several challenges arise when trying to define a discrete calculus framework on such surfaces mimicking the continuous one: the vertex positions and the geometry of faces do not capture well the geometry of the underlying smooth Euclidean object, even when refined asymptotically. Furthermore, the surface may not be a combinatorial 2-manifold even for discretizations of smooth regular shape. In this paper, we adapt a discrete differential calculus defined on polygonal meshes to the specific case of digital surfaces. We show that this discrete differential calculus accurately mimics the continuous calculus operating on the underlying smooth object, through several experiments: convergence of gradient and weak Laplace operators, spectral analysis and geodesic computations, mean curvature approximation and tolerance to non-manifold locii.

Introduction

In many geometry processing and simulation applications, solving variational problems or simulating partial differential equations on the object boundary can be a critical step. Involved shapes are generally represented as discrete meshes. It is thus necessary to have a calculus that operates consistently and accurately on these discrete objects. For embedded graphs or triangular manifold meshes, finite element methods [START_REF] Arnold | Differential complexes and stability of finite element methods i. the de rham complex[END_REF][START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF] or discrete exterior calculus (DEC) [START_REF] Desbrun | Discrete exterior calculus[END_REF][START_REF] Grady | Discrete calculus: Applied analysis on graphs for computational science[END_REF] have played a crucial role in many applied mathematics and geometry processing applications. From a broad perspective, such models provide consistent differential operators to process scalar, vector or tensor functions on meshes or embedded graphs (consistency given by satisfying Stokes' theorem for all discrete elements). These frameworks induce several convergence results for PDE solutions, but with strong hypotheses on the discrete-continuous mapping [START_REF] Hildebrandt | On the convergence of metric and geometric properties of polyhedral surfaces[END_REF]. On generic non-triangular meshes, Virtual Element Methods [START_REF] Beirão Da Veiga | The hitchhiker's guide to the virtual element method[END_REF], or ad hoc operators exist [START_REF] De Goes | Discrete differential operators on polygonal meshes[END_REF]. However, these models assume a Euclidean embedding interpolating the smooth manifold.

Digital objects and surfaces corresponds to discrete approximation of continuous objects through a discretization process [START_REF] Klette | Digital geometry: Geometric methods for digital picture analysis[END_REF], or to partitions in images. In terms of perturbation and stability, digital surfaces, made of isothetic unit squares in Z 3 , are very specific: vertices do not interpolate the continuous object, and geometric normals poorly approximate the continuous normal bundle. Even worse, the primal quad surface may not be a combinatorial 2-manifold. To design stable geometric estimators, a key ingredient is the digitization grid step h used to represent objects. Hence, a multigrid process, as a function of h, can be designed to relate a digital surface to the underlying smooth surface [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF]. Stable geometric estimators with convergence properties can then be obtained (i.e. the estimation converges to the expected one on the smooth manifold as h tends to zero) for various quantities: surface area [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF], curvature tensor [START_REF] Lachaud | Robust and Convergent Curvature and Normal Estimators with Digital Integral Invariants[END_REF], or even higher order functional such as the Laplace-Beltrami operator [START_REF] Caissard | Laplace-beltrami operator on digital surfaces[END_REF].

A combinatorial DEC can be constructed on digital surfaces, but their specific geometry makes it poorly reflects the continuous calculus. This article presents a new discrete calculus framework dedicated to digital surfaces, which relies on two ingredients: (i) a convergent normal vector field u (e.g. the integral invariant normal estimator [START_REF] Lachaud | Robust and Convergent Curvature and Normal Estimators with Digital Integral Invariants[END_REF]), which is used to correct the embedding of (ii) a polygonal differential calculus model of de Goes et al. [START_REF] De Goes | Discrete differential operators on polygonal meshes[END_REF]. Several methods exploits this idea of correcting the embedding with a normal vector field [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF][START_REF] Caissard | Laplace-beltrami operator on digital surfaces[END_REF][START_REF] Lachaud | Interpolated corrected curvature measures for polygonal surfaces[END_REF]. Mercat [START_REF] Mercat | Discrete complex structure on surfel surfaces[END_REF] follows this idea with a theory of conformal parametrization and differential operators for digital surfaces restricted to combinatorial 2-manifolds. Our proposal shares some ideas with these works and defines a calculus on generic digital surfaces with a simple per-face construction of the operators.

The paper is organized as follows. We first describe the operator construction from [START_REF] De Goes | Discrete differential operators on polygonal meshes[END_REF] ( §2). §3 describes how to correct the geometrical embedding of the surface elements. Finally, we evaluate the performance of the framework on various variational problems ( §4).

Polygonal differential calculus

We focus here on the formulation proposed by de Goes et al. for polygonal surfaces embedded in R 3 [START_REF] De Goes | Discrete differential operators on polygonal meshes[END_REF]. It defines differential operators per face without assumption on the face geometry, which could be non-convex or even non-planar.

Let M be a mesh with vertices X and faces F. For a given face f with n f vertices, we denote by X f , the vertex positions encoded as a n f × 3 matrix (the i-th row corresponds to the position of the vertex i of the face f ). In this calculus, we consider scalar functions ϕ defined on vertices of M (see Figure 1-(a-b)). For the face f , we denote by ϕ f the restriction of ϕ to the face vertices, represented as a vector of size n f . As all discrete differential estimators will be linear in the vertex positions and the scalar function values defined at vertices, the matrix representation will easily combine operators with matrix products. For instance, the centroid of a face is given by c

f := 1 n f X t f 1 f , with 1 f = (1 . . . 1) t of size n f .
First per face local operators are constructed, starting from the gradient. To cope with (non-planar) polygonal faces, the weak form of the gradient is sought (i.e. f ∇ϕ(x)dx in the continuous setting) leading to a constant gradient integrated per face. Solving the integral on the discrete structure would require an interpolation scheme on the non-planar face. The key ingredient is to focus on the co-gradient ∇ϕ ⊥ that leads to a simpler expression of its weak form using Stoke's theorem:

f ∇ϕ ⊥ (x)dx = ∂f ϕ(x)t(x)dx ,
t(x) being the unit tangent vector at x ∈ ∂f . When discretizing this form for a polygonal face f , t(x) corresponds to the (normalized) edge vectors. Now, let us consider a linear function ϕ f = X f s + 1 f r with s ∈ R s and r ∈ R, integrating ϕ f on an edge correspond to averaging the function value at the endpoints. Hence, the discrete operator for the (integrated) co-gradient is the 3 × n f matrix

G ⊥ f := E t f A f ,
E f being the matrix encoding the face edges, and A f the operator that averages two consecutive vertex values. To be explicit, A f is the n f × n f matrix with A i,i = 1/2, A i,i+1 = 1/2 (zero otherwise), and D f denotes the n f ×n f derivative matrix with

D i,i = -1, D i,i+1 = 1 (zero otherwise). So E f = D f X f .
Given a vector ϕ f , E t f A f ϕ f is a vector of size 3 corresponding to the Euclidean embedding of the integrated co-gradient of ϕ f on f . For a polygonal face f , the vector area a f can be defined as [START_REF] Sullivan | Curvatures of smooth and discrete surfaces[END_REF]:

a f := 1 2 vi∈f x vi × x vi+1 .
The normal vector to f is simply n f := a f /a f , with a f := ∥a f ∥ 2 . Now, we can finally get the expression of the gradient as ∇ϕ ⊥ = [n(x)]∇ϕ, where [n] is the π/2 rotation matrix such that [n]q = n × q for any 3d vector q:

G f := 1 a f [n f ]E t f A f .
The gradient operator is a 3 × n f matrix outputting a 3d vector in 3d when applied to a scalar vector ϕ f (see Figure 1-(c)). Following such per-face construction, de Goes et al. define several differential operators listed in Table 1 and represented as matrices acting on scalar functions, vectors, or discrete differential forms (see [START_REF] Desbrun | Discrete exterior calculus[END_REF][START_REF] Grady | Discrete calculus: Applied analysis on graphs for computational science[END_REF] for an introduction). We do not discuss the rationale and the construction of each operator, please refer to [START_REF] De Goes | Discrete differential operators on polygonal meshes[END_REF] for details. We conclude with the presentation of the Laplace-Beltrami operator ∆, a fundamental tool in many geometry processing applications [START_REF] Lévy | Spectral mesh processing[END_REF]. To solve global PDE involving this operator (e.g. for instance solving a Laplace or Poisson problem ∆u = f ), it is interesting to aggregate the per-face operators into a global operator L. We simply sum up the L f matrices with a global indexing of the vertices. As every L f matrix is negative semi-definite, so is the global Laplace-Beltrami operator. Furthermore, the global operator is very sparse. In opposition with [START_REF] De Goes | Discrete differential operators on polygonal meshes[END_REF], we define it with a negative sign, since usually both the Laplacian and the Laplace-Beltrami are negative operators (with negative eigenvalues). We do not go further into the details of this approach, several results are given in [START_REF] De Goes | Discrete differential operators on polygonal meshes[END_REF] to justify the construction of these operators. For short, they match with classical discrete exterior calculus or finite element ones for triangular meshes [START_REF] Desbrun | Discrete exterior calculus[END_REF], and for the polygonal case, they resemble to Virtuel Element Methods operators for 2d structures embedded in R 3 [START_REF] Beirão Da Veiga | Basic principles of virtual element methods[END_REF]. However, all these relationships between discrete calculus operators and their continuous counterpart on the smooth object only make sense when the discrete structure interpolates the smooth object, with a close normal vector field. Applying this calculus as is on digital surfaces provide poor results. We need to adapt the surface embedding to correct it.

Operator

Size Description V f := E f (I3×3 -n f n t f ) n f × 3 flat operator that maps a vector to a 1-form.

U f := 1 a f [n f ](B t f -c f 1 t f ) 3
× n f sharp operator that maps a 1-form to a 3d vector.

P f := In f ×n f -V f U f n f × n f projection operator acting on a 1-form. M 0 f := a f n f In f ×n f n f × n f inner product for discrete 0-forms (functions). M 1 f := a f U t f U f + λP t f P f n f × n f inner product for discrete 1-forms (for some λ > 0). Div f := -D t f M 1 f n f × n f integrated divergence operator from a 1-form. L f := -D t f M 1 f D f n f × n f weak Laplace-Beltrami operator.
Table 1. Summary of local per-face operators of polygonal calculus [START_REF] De Goes | Discrete differential operators on polygonal meshes[END_REF].

Projected digital surface embedding

As discussed in the introduction, the classical approach is to relate a digital object to its continuous counterpart through the Gauss digitization process [START_REF] Klette | Digital geometry: Geometric methods for digital picture analysis[END_REF]. In this setting, many multigrid convergence results have been obtained for various integral and differential quantities such as the length in 2d [START_REF] Coeurjolly | A comparative evaluation of length estimators of digital curves[END_REF], the surface area in 3d [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF], the curvature tensor [START_REF] Lachaud | Robust and Convergent Curvature and Normal Estimators with Digital Integral Invariants[END_REF], or the Laplace-Beltrami operator [START_REF] Caissard | Laplace-beltrami operator on digital surfaces[END_REF]. Among these techniques, the convergent estimation of the normal vector bundle is the cornerstone of many followup results (e.g. [START_REF] Lachaud | Properties of gauss digitized shapes and digital surface integration[END_REF][START_REF] Caissard | Laplace-beltrami operator on digital surfaces[END_REF]).

Fig. 2. Input digital surface (left), its estimated normal vectors using [START_REF] Lachaud | Robust and Convergent Curvature and Normal Estimators with Digital Integral Invariants[END_REF] (middle), projections of the face vertices onto their respective estimated tangent plane (right).

In the polygonal calculus described above, a key ingredient is that operators are first defined per-face, and later summed up globally with the vertex indexing. This suggests that correcting the embedding of the face vertices during the local construction of the operators would lead to globally corrected operators. As a first-order approximation, we propose to define our projected polygonal calculus by relying on an external estimation of the face tangent space, and by (implicitly) projecting each face vertex onto this tangent plane when constructing the operators. The idea of incorporating external normal vectors to correct the embedding of a discrete structure follows the idea of Lachaud et al [START_REF] Lachaud | Interpolated corrected curvature measures for polygonal surfaces[END_REF]. Let {u f } be a normal vector field estimated on the digital surface, for instance using the multigrid convergent approach of [START_REF] Lachaud | Robust and Convergent Curvature and Normal Estimators with Digital Integral Invariants[END_REF].

Definition 1 (Embedding operator). The per-face f projection operator onto the tangent plane orthogonal to u f passing through the origin is the 3 × 3 matrix

Π f := (I 3×3 -u f u t f ).
Recalling that X f is the n f × 3 matrix of the vertices position of the face f , the new positions are given by X * f := X f Π f . Note that the intercept does not need to be specified as all differential operators are invariant by translation. For illustration purpose, projected faces have been translated to keep their centroid invariant in Figure 2, but this is meaningless in the calculus. Another observation is that for a given vertex v in M, its embedding generally differs for its adjacent faces. This follows from the per-face construction of the operators. The global continuity of scalar functions, or their processing (i.e their Laplacian) comes from the global indexing of vertices, and global operators as defined above.

Experiments

We demonstrate the interest of the projected calculus model for the processing of scalar functions on digital surfaces. More precisely, we first evaluate the relevance of the new embedding for the gradient computation of a scalar function. Next, we study the Laplace-Beltrami operator and more advanced processing. All operators are implemented in C++ and available in DGtal [START_REF]The DGtal Project[END_REF].

Gradient of scalar functions. To evaluate the relevance of including external tangent information in the calculus, let us consider the simple case of the gradient evaluation on a quadratic scalar field defined on a tilted digital plane (with normal vector 1/3

• ( √ 3, √ 3, √ 3
) t , Fig. 3). In a multigrid setting, we compute the per-face gradient of the scalar function (Euclidean distance in ambient space with respect to the plane center, Fig. 3-(a)), and we compare the expected gradient vector direction (Fig. 3-(b)), to the estimated gradient using the classical polygonal calculus (Fig. 3-(c)), and the projected one (Fig. 3-(d)). Gradient vectors are illustrated from their projection onto the Euclidean plane. The error metric used here is the l 2 norm of the difference between estimated and expected normalized gradient vectors. Figures 3-(e) and (f ) detail relative errors. In Figure 3-(g), we have considered a decreasing grid step h from 1 to 1.5 × 10 -2 (3 748 021 faces), demonstrating multigrid convergence of the gradient direction for the projected calculus. Laplacian of a function. First, let us consider an evaluation of the output of the Laplace-Beltrami function on analytical functions on a domain for which we already know the expected Laplacian of that function. Following the setting in [START_REF] Caissard | Laplace-beltrami operator on digital surfaces[END_REF], we evaluate L on a digital sphere with the simple quadratic function ϕ(x) := x 2

x on its boundary, and expect ∆ϕ(x) = 2 -6 cos 2 (ψ) 2 (1 -cos(θ) 2 ) in spherical coordinates. Note that Caissard et al. [START_REF] Caissard | Laplace-beltrami operator on digital surfaces[END_REF] achieve strong consistency of the operator (i.e. pointwise convergence of the Laplacian values) at the price of a convolution process on the surface elements, leading to a quadratic algorithm for the operator construction, and a non-sparse operator. Clearly it is unlikely that the purely local construction of our projected polygonal calculus provides pointwise convergence. We observe nevertheless that the projection process offers reasonable estimations while being purely local (see Figure 4). [START_REF] Caissard | Laplace-beltrami operator on digital surfaces[END_REF] in (c). In (d) we have the results of the original local polygonal calculus operator, and in (e) the result of our projected one. In (f ) we have performed a multigrid convergence test with decreasing gridstep h. Dashed lines correspond to the maximum absolute error when estimating ∆f , the solid ones correspond to mean absolute errors. Spectral analysis. We evaluate the stability of the spectrum of the pointwise Laplace-Beltrami operator (M 0 ) -1 L. Eigenvalues and eigenvectors have a fundamental role in many geometry processing applications [START_REF] Lévy | Spectral mesh processing[END_REF] as they define a (spectral) basis to represent scalar functions on surfaces. Figure 5 compares the first eigenvectors (corresponding to the largest eigenvalues) of the original polygonal calculus Laplace-Beltrami operator, to the ones we propose. Our embedding operator has a clear positive impact on the smoothness of basis functions. Fig. 5. First ten eigenvectors of the Laplace-Betrami operator on a digital sphere (and their projection onto the Euclidean sphere): top two rows correspond to the original polygonal calculus, bottom two rows correspond to the projected polygonal calculus we propose (note that due to proximity of eigenvalues and numerical issues, eigenvectors may not appear in the exact same order).

Mean curvature estimation. A well-known result related to the Laplace-Beltrami operator is that, if x is the function associating to any point of M its coordinates, then ∆x = 2H(x)n(x) with H(x) the mean curvature at x and n(x) the normal to M. We check this method for computing the mean curvature with our calculus. We compute the pointwise Laplace-Beltrami of coordinates functions (X, Y, Z) = x and deduce the mean curvature H by scalar product with the estimated normal u. Otherwise said, the estimated curvature Ĥ is

Ĥ = 1 2 (M 0 ) -1 LX LY LZ u t ,
where u is a n f × 3 matrix storing per row the corrected normal to the vertex.

To account for the fact that the Laplace-Beltrami is only weakly convergent, we can diffuse for a short time the result to simulate a local integration. Figure 6 illustrates this method for approximating the mean curvature, and also confirms that the calculus must be corrected to get meaningful result.

Geodesic distance estimation. To evaluate the proposed calculus on a more complex example, let us consider the geodesic distance estimation problem on digital surfaces. We consider the PDE approach of Crane et al. [START_REF] Crane | The heat method for distance computation[END_REF] that uses heat diffusion to estimate the geodesic distance function from sources. Sources are defined via a characteristic function u 0 in the domain. The heat method consists in three steps: (i) integrate the heat flow u = ∆u for some fixed time t with initial condition u(x, 0) = u 0 (x), (ii) normalize the gradient vector field d := -∇u t /∥∇u t ∥ 2 , and (iii) solve the Poisson equation ∆ϕ = ∇ • d (∇• denotes the divergence operator on vector fields). The scalar function ϕ -min x (ϕ) is a good approximation of the geodesic distance from the sources to any point [START_REF] Crane | The heat method for distance computation[END_REF]. Fig. 6. Approximation of mean curvature H with Laplace-Beltrami operator. For the computation of the diffused approximation Ĥ, the diffusion time t is set to 0.04 < 1 6 h, so very short. The projected calculus uses integral invariant normal estimator (default parameter r = 3h 1/2 ). Projected calculus results are more accurate and stable.

In our discrete calculus framework, the heat diffusion step (i) is obtained with a single backward Euler step for timestep t solving

(u t -u 0 )/t = (M 0 ) -1 Lu t
or more simply (M 0 -tL)u t = M 0 u 0 , which involves the global mass matrix M 0 and the weak Laplace-Beltrami L.

Then, denoting w the solution u t , step (ii) computes per-face the normalized opposite gradient d f of w f for face f as d f := -G f w f /∥G f w f ∥ 2 . Bringing back this vector to a one-form with flat operator V f , the per-face weak divergence is then δ f := Div f V f d f . Vector δ f is a n f × 1 vector associated to face f , assigning values to each vertex of f . The global divergence map δ simply consists in summing up per vertex the divergence contribution of each face.

Finally, step (iii) computes the geodesic distance function ϕ by solving the Poisson linear system Lϕ = δ, which must be shifted by -min x (ϕ) so that the result starts from distance 0.

Figure 7 illustrates geodesic distance computations on a digitized sphere, both with naive calculus and projected calculus (we use integral invariant normal estimations with default parameters). As shown by projecting back distances on the true smooth surface, the projected calculus builds more isotropic isodistance lines while being much more numerically accurate.

Distances on surfaces with boundaries. If the digital surface has boundaries, then the heat diffusion as it is written follows Neumann boundary conditions: the Laplace-Beltrami operator L is built per-face and ignores boundaries. Isodistance lines will tend to be orthogonal to boundaries (see Figure 8 for an example on a half-sphere). As suggested in [START_REF] Crane | The heat method for distance computation[END_REF], we compute a second diffusion in step (i) that assumes Dirichlet zero boundary conditions on the digital surface boundaries. This is done by restricting the linear system (M 0 -tL)u t = M 0 u 0 to non-boundary vertices, simply by assuming value 0 for u on boundary vertices. We then average the two solutions to define the vector w that is given to step (ii) of the heat method. Figure 8 confirms the soundness of this approach.

Processing on non-manifold surfaces. Since every calculus operator is defined per-face, they are quite oblivious to non-manifold edges and vertices of the digital surface. We illustrate this fact on Figure 9 by computing the geodesic distance to a source point on a digital surface approximating the triple junction of three planes. Not a single line of code is changed.

Conclusion

Our contribution is a simple and easily implementable discrete differential calculus for the processing of scalar functions on digital surfaces. It relies on changing the natural embedding of the digital surface when constructing per face differential operators, with the use of an external, multigrid convergent, normal vector estimation. Many experiments back up the effectiveness of this new calculus.

Although we demonstrate its interest to solve integrated problems such as Poisson problem (and the geodesic distance estimation is a perfect example of this class), the proposed local construction does not achieve pointwise convergence for second order operators like the Laplace-Beltrami one. Interesting challenges include the design of nonlocal operators, similarly to [START_REF] Caissard | Laplace-beltrami operator on digital surfaces[END_REF] or [START_REF] Hildebrandt | On approximation of the laplace-beltrami operator and the willmore energy of surfaces[END_REF], or the use of higher order embeddings for the digital surface.

Fig. 1 .

 1 Fig. 1. Illustration of a non-planar face f (a) equipped with a vertex valued scalar function ϕ f (b), and the co-gradient and gradient vectors of ϕ f (c)

Fig. 3 .

 3 Fig. 3. Convergence of the gradient operator: Input scalar function on a tilted plane (a), expected normalized gradient projected onto the Euclidean plane (b), gradient vectors using the original polygonal calculus with many aliasing artifacts (c), gradient vectors using the projected calculus (d). Figures (e) and (f ) details relative errors of gradient estimations (values respectively in [0, 0.28] and [0, 0.62]), and (g) illustrates l2 norm of the error when decreasing the grid step h (from right to left).

Fig. 4 .

 4 Fig.[START_REF] Coeurjolly | A comparative evaluation of length estimators of digital curves[END_REF]. Given a function ϕ(x) in (a) and its analytical Laplacian (b), we illustrate the result of the pointwise convergence heat kernel based Laplace-Beltrami operator[START_REF] Caissard | Laplace-beltrami operator on digital surfaces[END_REF] in (c). In (d) we have the results of the original local polygonal calculus operator, and in (e) the result of our projected one. In (f ) we have performed a multigrid convergence test with decreasing gridstep h. Dashed lines correspond to the maximum absolute error when estimating ∆f , the solid ones correspond to mean absolute errors.

Sphere of radius 1 ,Fig. 7 .

 17 Fig.7. Distance computation to a source point with geodesic in heat method on a digital sphere. Scale is deep blue (0) to red (max distance) with black isolines every 0.1 and red isolines every 1. The expected max distance is π. Top row displays distances computed with naive polygonal calculus. Bottom row displays distances computed with projected calculus with normals u estimated with integral invariant method. All computations made with parameter λ = 0.25 and heat diffusion time t = h 2 .

Fig. 10 .

 10 Fig. 10. Various computations on a more complex voxel object: (from left to right) the gradients of a scalar function, mean curvature and geodesic distances.

  Goursat polynomial surface, digitized at gridstep h = 0.25

	Smooth surface	Naive calculus	Projected calculus
	True H	Pointwise Ĥ	Diffused Ĥ	Pointwise Ĥ	Diffused	Ĥ
	∥H -Ĥ∥2	0.059	0.051	0.025	0.003
	∥H -Ĥ∥∞	0.197	0.146	0.155	0.022

⋆ This work has been partly funded by CoMeDiC ANR-15-CE40-0006 research grant.


Amalia Todirascu

Rodrigo Wilkens

Eva Rolin

Thomas François

Delphine Bernhard

Nuria Gala

A Hybrid

Text Simplification

HECTOR: A Hybrid TExt SimplifiCation TOol for Raw Texts in French

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