Apollinaire Barme

Simon Tamayo

Arthur Gaudron

Statistic regression and open data approach for

Keywords: e-commerce, statistical modeling, regression, empirical research

I. INTRODUCTION

-commerce embodies a particularly complex system that challenges logistics systems in today's environment. Online shopping is related to a high fragmentation of receivers and deliveries of smaller quantities. The result is that the last mile distribution is the most expensive link of the supply chain [START_REF] Allen | Understanding the impact of e-commerce on last-mile light goods vehicle activity in urban areas: The case of London[END_REF] [START_REF] Gevaers | Cost Modelling and Simulation of Last-mile Characteristics in an Innovative B2C Supply Chain Environment with Implications on Urban Areas and Cities[END_REF]. As a result, e-commerce pushes the traditional logistic approaches based on massification [START_REF] Yu | E-commerce Logistics in Supply Chain Management: Practice Perspective[END_REF] [4] [START_REF] Jain | e-Commerce and supply chains: Modelling of dynamics through fuzzy enhanced high level petri net[END_REF] in order to remain profitable.

The increasing access to Internet worldwide is allowing consumption through e-commerce in more countries every year. During the past 20 years, although e-commerce in developed countries has been present, it did not grow, as rapidly as expected. It is mostly the past five years, with the adoption of smartphones, tablets and the omnipresence of Internet, that e-commerce has displayed its exponential growth [START_REF] Winter | The rise of cyberinfrastructure and grand challenges for eCommerce[END_REF].

This paper aims at identifying linear relations between national e-commerce sales and major economic and demographic indicators. This problem has not been explicitly addressed in the literature, but previous works identified important elements on which this paper intends to build.

In 2007, Ho et al. [START_REF] Ho | A growth theory perspective on B2C e-commerce growth in Europe: An exploratory study[END_REF] addressed the problem of modeling A. B. Author was with MINES ParisTech at the Center for Robotics (email: apollinairebarme@gmail.com). S. T. Author is with MINES ParisTech at the Center for Robotics (corresponding author, e-mail: simon.tamayo@mines-paristech.fr, phone: +33140519452 )

A. G. Author ia with MINES ParisTech at the Center for Robotics (e-mail: arthur.gaudron@mines-paristech.fr).

the growth of e-commerce revenues at a national level, identifying three underlying mechanisms: endogenous, exogenous and a mixed endogenous-exogenous. The study of Ho et al. highlighted the disparity of influences of internal and external drivers on e-commerce growth across countries.

Moreover in 2011, Ho et al. [START_REF] Ho | Internet-based selling technology and e-commerce growth: A hybrid growth theory approach with cross-model inference[END_REF] formalized a hybrid growth theory with cross-model inference, which allowed to conclude that "both endogenous variables (e.g. internet user penetration, investment in telecommunication) and exogenous variables (international openness) drive the GDP-normalized level of B2C e-commerce revenues" [START_REF] Ho | Internet-based selling technology and e-commerce growth: A hybrid growth theory approach with cross-model inference[END_REF].

Mahmood et al. [START_REF] Mahmood | On-line shopping behavior: Cross-country empirical research[END_REF] studied the online shopping behavior of consumers and concluded that e-commerce sales are influenced by two major indicators: namely: (a) economic condition and (b) trust.

The work of Tan et al. [START_REF] Tan | Business-to-business adoption of eCommerce in China[END_REF] indicated that in China, important factors related to B2B e-commerce adoption are: (a) access to computers, (b) lack of internal trust, (c) lack of enterprise-wide information sharing, (d) intolerance towards failure, and (e) incapability of dealing with rapid change.

Meso et al. [START_REF] Meso | Information infrastructure, governance, and socio-economic development in developing countries[END_REF] studied the relationships between information infrastructure and socio-economic development.

Their study pointed out that e-commerce adoption is highly related to (a) telephone density, (b) Internet density, (c) GDP, (d) purchasing power parity, among others.

Bejarano [START_REF] Bejerano | Challenges of forecasting demand for e_commerce[END_REF] proposed forecasting models for e-commerce demand in the U.S.A. and highlighted the difficulty of finding independent variables to predict an unbiased and accurate ecommerce demand.

The following sections of this paper are organized as follows. Section 2 presents the problem statement and the input data. Section 3 presents the proposed approach in which a theoretical background abound linear regression is given and the main steps of the methodology are explained. Section 4 presents the application of the methodology to a case study of 24 countries. And finally Section 5 concludes and opens future research perspectives

II. PROBLEM STATEMENT AND INPUT DATA

The problem addressed in this paper consists in identifying the major economic and demographic variables that influence e-commerce sales.

A methodology is proposed to relate the data of ecommerce sales of several countries, to a set of openly available indicators Apollinaire Barme, Simon Tamayo, Arthur Gaudron

Statistic regression and open data approach for identifying economic indicators that influence e-commerce E Fig. 1 Overview of the proposed methodology

The presented methodology allows constructing and validating a linear regression model to determine a subset of explanatory variables as shown in Fig 1 . As a result, this paper deals with the following research questions:

• Which major variables influence e-commerce sales?

• How to identify these variables?

• Would it be possible to predict the sales of e-commerce in a country from these variables?

The figures of e-commerce sales of 24 countries of were obtained from the Global B2C E-commerce Report 2016 [START_REF]Global B2C E-commerce Report[END_REF].

Economic and demographic data about these 24 countries were obtained in the "Country Profiles Data" in the World Bank Open Data Portal [START_REF]Country Profiles Data[END_REF]. An initial set of 36 variables was queried. Table I presents a synthesis of these input variables.

III. PROPOSED APPROACH

The proposed approach is based on a general linear model. This section is presented in three parts. First, the hypotheses of the proposed model are listed. Second, an overview of the theoretical background about linear regression is proposed. And third, the proposed methodology is explained.

A. Hypothesis

The proposed methodology is based on a linear regression model. Accordingly, it makes the following hypotheses about the explanatory variables (i.e. indicators) and the dependent variable (i.e. e-commerce sales):

Linearity: the mean of the e-commerce sales result is a linear combination of regression coefficients and the explanatory variables. This assumption is not as restrictive as it seems, because the explanatory variables themselves can be transformed, and in fact multiple copies of the same variable can be added, each one transformed differently [START_REF] Seber | Linear Regression Analysis[END_REF], leading for example to a polynomial regression function.

Homoscedasticity: the errors in the e-commerce sales results have a constant variance, regardless of the values of the explanatory variables.

Independence of errors: the errors of the e-commerce sales result are uncorrelated with each other.

Lack of multicollinearity: the explanatory variables should not be highly correlated [START_REF] Azais | Le modèle linéaire par l'exemple[END_REF]. In order to comply with this

China

E-commerce sales (US$ billions)

hypothesis a correlation test is performed in the data collection step and the input parameters are narrowed down.

B. Theoretical background

The approach proposed in this paper applies some concepts of statistics that are recalled hereafter.

1) General linear model

A general linear model is used to represent the relationships between e-commerce sales denoted Y (modeled as scalar dependent variable) and a set of explanatory variables denoted X (i.e. the economic and demographic information about each country). Let Rn with its scalar product, its norm and its distance be our reference space. The General Linear Model is represented by the following equation: (1) We assume that: X and Y are both known. Y is a column vector with n entries that contain the dependent variable data. In our case there are 24 countries, therefore n=24. X is matrix containing the explanatory variables data with the first column filled ones. In order to ensure the existence of a single solution to the linear problem, the number of explanatory variables (p) needs to be less than n-1. The vector θ contains the coefficients of the linear model. The aim of the regression model is to determineθ .

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + = n p
In order to ensure a uniform variance of the expected errors in the model, the Homoscedasticity criterion should be validated (cf. previous section). This criterion is validated if the average expected value of the noises is equal to zero and their variance is constant.

2) Regression

The regression follows the steps shown in [START_REF] Gevaers | Cost Modelling and Simulation of Last-mile Characteristics in an Innovative B2C Supply Chain Environment with Implications on Urban Areas and Cities[END_REF]:

θ = arginf Y -Xθ ( ) 2 [X] = Xθ,θ ∈ ℜ p+1 { } but [X] is a finite vector space Hence, X θ = argmin d(Y,[X])
To conclude, due to unicity, X θ

= P[X](Y )
Then, the following property is satisfied:

Y-P[X](Y ) ⊥ X Which leads to Y -P[X](Y ), X = 0 We deduce that, X T (Y -X θ ) = 0 Then, X T Y = X T X θ But, X T X ∈ Mp + 1(ℜ) and Rank(X T X) = p +1 Consequently, (X T X) -1 exists To conclude, θ = (X T X) -1 X T Y (2)
Our aim is to find θ . The linear model hypothesis establishes that Y = Xθ +ε with ε unknown. To find θ , we look for Xθ the closest to Y . The set X is known, so finding θ or Xθ is equivalent. We then place [X] as the working space (which is a finite vector space), and we minimize the distance between Y and [X] , which is obtained for X θ , that is, the orthographic projection of Y over [X] , as the straight line is that the shortest path between two points.

3) Validation

Many inputs are involved in a regression, some of which are irrelevant. In order to identify a significant subset of explanatory variables, two types of tests are preformed: Student's and Fisher's.

a) Student's Test The student's test is performed in order to ascertain if a single variable is meaningful to the regression model. This is quantified with the student's p-value of each input. For example, in order to verify if θ 1 is a significant parameter, the next test is performed:

H 0 : C T θ = 0 H1 : C T θ ≠ 0 C T = (0,1,...,0) (3)
We define T the statistic of the test:

T = C T θ σ C T X T X ( ) -1 ( ) -1/2 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ -1 T follows St n -p -1 ( ) (4)
If P(|St(n-p-1)|>T) > 0.05 then we reject H 1 else we reject H 0 . In other words, if H 1 is validated, the parameter is significant. P(|St(n-p-1)|>T) is known as the "p-value" (for more insight on the theoretical aspects of this test see [START_REF] Azais | Le modèle linéaire par l'exemple[END_REF]).

b) Fisher's Test The fisher's test is performed in order to detect if a set of inputs is significant. It is important to note that the student's test can yield a set of parameters that are significant individually, but that are not relevant as a set. This is why the fisher's test often complements the student's test. For example, in order to verify if the set θ 1 ,θ 2 ,θ 3 is significant, the next test is performed:

H 0 : C T θ = 0 H1 : H 0 is false C T = 0 1 -1 0 0 ... 0 0 0 1 -1 0 ... 0 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ (5)
We define T the statistic of the test:

T = θ T C C T X T X ( ) -1 C ( ) -1 C T θ σ -2 T follows F p, n -p -1 ( ) (6
) If P(|St(n-p-1)|>T) > 0.05 then we reject H 1 else we reject H 0 . As result, the fisher's test returns one single p-value result. If this result is less than or equal to 0.05, the set of parameters is considered to be significant.

C. Proposed methodology

The proposed methodology is divided in three main steps:

(1) Data collection; (2) Regression and selection; and (3) Validation, presented as shown in Fig. 3.

4) Data collection

The presented analysis is based on a sample of 24 countries (n=24). In order to be exhaustive, we queried as many parameters as we could. Nonetheless, as indicated previously, the number of explanatory variables needs to be less than n-1. As a result only 22 indicators were selected. We measured the correlation between the data and removed the highly correlated indicators. This was done by calculating the correlations matrix and eliminating the correlations higher than 85%. This procedure was repeated iteratively until 22 inputs were obtained.

5) Regression

Once the data is selected, we perform the linear regression following the steps in (2) and we compute the Student's pvalues for the different parameters, the Fisher's for the set and the R², this latter, corresponds to the percentage of variances explained, which is a measure of how the model is fitting the observations. The Fisher's p-value is calculated for the entire set of input variables (22). If the p-value is above 0.05 the set is globally irrelevant. In this case, we use the Bayesian Information Criterion (BIC), which allows selecting the best-fitted set of variables. The BIC returns a more relevant subset of variables that are fed again to the regression model.

6) Validation

Once the subset of parameters is validated in terms of pvalues, three graphs associated with the regression are analyzed:

1. "Residuals versus fitted" shows the positions of the residuals; 2. "Scale location" shows the variances of the residuals; 3. "Residuals vs. leverage" allows detecting aberrant points (if any). If the variances of the residuals are not uniformly distributed along the axis, the residuals are considered to be heteroscedastic. In this case, the linear regression cannot be applied (cf. hypotheses of the model in section III.A). If such is the case, we perform a Box-Cox transformation on the dependent variable in order to obtain homoscedastic residuals.

The Box-Cox transformation is a power transformation, which determinates the proper power on the variables to ameliorate the regression. After the Box-Cox succeeds, regression is re-applied with the transformed variables.

At the end of the methodology, if all p-values are below 0.05 (Fisher's and Student's), the R² is close to 1, the epsilons are homoscedastic and normally distributed, the regression is validated and the final subset of explanatory indicators is obtained.

IV. APPLICATION TO E-COMMERCE SALES

The application presented in this section was implemented within the "R" software using the libraries "MASS" and "CAR".

A. Data collection

The proposed methodology was applied to a sample of 24 countries. An initial set of 36 explanatory variables was considered (cf. section II). After filtering the input data, significant correlations were identified, and the initial set was reduced to 22 variables on the grounds of the correlation coefficients.

B. Regression

A linear regression was performed following the steps shown in (2) and we obtained an adjusted R2 equal to 0.4654 (which is not yet satisfactory) and a Fisher's p-value of 0.5231. Table II presents a summary of the regression results for the 22 variables. The Fisher's p-value result (higher than 0.05) indicates that this set of variables is not significant (cf. section III.B.2). As result, we use the BIC algorithm in order to filter the input variables. After running the BIC procedure, we obtained a subset of 14 variables, presented in table III.

The Bayesian Criterion indicates that these variables should be explanatory inputs for the sales observation. The BIC column in Table III indicates that the 14 variables are relevant (i.e. the BIC result of each variable is greater than that of "none").

Evaluate correlation

Filter linearly correlated variables

Variables correlated?

P-values ok?

Compute Bayesian information criterion

Validate normality and homoscedasticity

Results

ok?

End

Filter and transform variables

Yes

No

Run linear regression

No

Yes

No

Yes

Input variables Once this new set is established we run once more the linear regression and we obtained an adjusted R 2 equal to 0.9139 and a Fisher's p-value of 6.306e-05. The results of this second regression are shown in Table IV.

E-commerce sales

Data collection

Regression and selection

Validation

The Fisher's p-value is below 0.05, as a result, the sample is significant. The R² is close to 1, nonetheless some variables have p-values bigger than 0.05 ("Population", "Gross capital formation" and "mobile cellular subscriptions"). We removed those variables one by one verifying that the resulting set had a valid fisher's p-value, and we obtained a new subset of 11 indicators. All the variables have their p-values below 0.05, the Fisher's p-value is adequate and the R² is more than sufficient. The summary suggests that the regression is meaningful. We now analyze the graphs to validate these results.

C. Validation

The scale location graph (Fig. 4) shows the homogeneity between the variances of the residuals. That is, the epsilons are uniformly spread along the axis, which indicates that they are homoscedastic. Moreover, the Breush-Pagan test [START_REF] Society | A Simple Test for Heteroscedasticity and Random Coefficient Variation Author ( s )[END_REF] reinforces the hypothesis of homoscedasticity of the residuals. From a practical point of view, this means that in these two countries e-commerce sales follow a different logic than that of the rest of the sample.

As mentioned in the theoretical background, the "residuals versus leverage" (Fig. 6) graph indicates atypical points (those that are behind the Cook's distance boundary [15] [16]). Fig. 6 points out that the regression would be improved if the first and the fourth observations (i.e. China and U.S.A.) were removed. Nonetheless, because of the size of the sample (only 24 observations) it is not possible to remove these points. It is important to note that the R 2 and the p-values indicate that the regression is relevant event with these two abnormal observations. Moreover, in Figures 4 and5, U.K. appears to be a remote point. Yet, in Fig. 6 the cook's distance shows that the U.K. is within the regular boundaries. As result, U.K. is considered to be an extreme point that remains relevant to the regression.

D. Discussion

36 variables were included in our initial data set. Due to this large number of parameters, we use the correlation between them to erase the correlated entries. Once the set was reduced to 22 variables, we proceeded to the regression, the result was unequivocal: the set was still not reliable. The BIC algorithm allowed to the narrow down the set to 14 variables. The regression was run again and three variables with high pvalues, were also removed. Finally, we obtained a final set of 11. A last regression indicated that this was a valid set (adjusted R 2 = 0.8915 and Fisher's p-value = 8.296e-06).

The regression indicates that the indicators with the higher influence on e-commerce sales are: (1) GDP, (2) Individuals using Internet, (3) High technology exports and (4) Population density. This is observed from the student's p-values results in Table V. It is important to note, that these indicators were also identified as influential parameters when assessing B2B ecommerce adoption in the works of Ho et al. [START_REF] Ho | Internet-based selling technology and e-commerce growth: A hybrid growth theory approach with cross-model inference[END_REF] and Meso et al. [START_REF] Meso | Information infrastructure, governance, and socio-economic development in developing countries[END_REF]. GDP is the most significant indicator the lowest pvalue of 3.35e-05. Fig. 7 overlaps the graphs of E-commerce sales and GDP for the studied countries, in which the linear correspondence is remarkable.

Other interesting variables, such as "Adolescent fertility", "Immunization measles" and "Terrestrial and marine protected" are present in the regression, these are surprising results, but that may quantify in a certain way the development of the countries.

It is to note that regression identifies China and U.S.A. as aberrant points (cf. Fig. 4 and Fig. 5). However it is not surprising that the homelands of Amazon and Alibaba have a behavior different than the average.

GDP (US$ billions) E-commerce sales (US$ millions)

SALES GDP

commerce sales to economic and demographic indicators. This paper addressed three research questions that are reviewed hereafter.

Which major variables influence e-commerce sales? A set of 11 explanatory variables has been established. It includes the following indicators in order of significance: (1) GDP, (2) Individuals using Internet, (3) High technology exports, (4) Adolescent fertility rate, (5) Population density, (6) Population growth, (7) Immunization, measles, (8) Fertility rate, (9) Surface area, [START_REF] Tan | Business-to-business adoption of eCommerce in China[END_REF] Terrestrial and marine protected areas and (11) Imports of goods and services.

How to identify these variables? A methodology was proposed in order to validate the significance of a set of input indicators. First the indicators are filtered on the grounds of linear correlation between them. Then a subset of independent indicators is fed into the linear regression model, which in turn, iteratively narrows these indicators in order to validate the hypotheses of linearity, independence and homoscedasticity. The methodology is presented in thee main steps: Data collection, Regression and selection and Validation.

Would it be possible to predict the aggregate sales of ecommerce in a certain country from these variables? A low pvalue and high R 2 results indicated that the resulting model is adequate. However, the input set of observations is still very limited: only 24 countries form the set of dependent variables. A sample size this small cannot yield relevant prediction results. In order to assess the predictive capabilities of the methodology, the model was build with the data of 23 countries, and the sales the remaining one were predicted. This was done with different combinations of countries, and the average error was about 40%.

A future version of this approach should take into account more sales observations and more input indicators. Also it would be interesting to perform similar regressions data of consecutive years in order to deduce parameters linked to ecommerce growth, such as the arrival of a big player, a significant improvement in the access to the Internet, etc.

The proposed approach has two important limitations. First, the linear model hypothesis supposes that the influence of the explanatory variables is uniform regardless of the country. For example, it admits that the GDP has as much influence on ecommerce sales in China as in Jordan, yet this assumption can be challenged. Second, the input data set is far from ideal: on the one hand, the sample size is small, which limits the number of indicators that can be considered (cf. section III) and on the other hand, there is no measure of the quality of the data.

Fig. 2 .

 2 Fig. 2. E-commerce sales (billions of US dollars) in 2015 per country.Source: "Global B2C E-commerce Report 2016"[START_REF]Global B2C E-commerce Report[END_REF]

Fig. 3

 3 Fig. 3 Detail of the proposed methodology

.

Fig. 4

 4 Fig. 4 Scale location

Fig. 7

 7 Fig.7 Overlap of the e-commerce sales and the GDP

TABLE I

 I

		.
	INPUT VARIABLES (SOURCE: WORLD BANK OPEN DATA)
	Variable	Unit
	1. Population	Millions
	2. Population growth	Annual %
	3. Surface area	Thousands of sq. km
	4. Population density	People per sq. km of land area
	5. GNI, Atlas method	Current US$ billions
	6. GNI per capita, Atlas method	Current US$
	7. GNI, PPP	Current international $ billions
	8. GNI per capita, PPP	Current international $
	9. Life expectancy at birth	Years
	10. Fertility rate	Births per woman
	11. Adolescent fertility rate	Births per 1,000 women ages 15-19
	12. Mortality rate, under-5	Deaths per 1,000 live births
	13. Immunization, measles	% of children ages 12-23 months
	14. School enrollment, primary	% gross
	15. School enrollment, secondary	% gross
	16. Forest area	Thousands of sq. km
	17. Terrestrial and marine protected areas % of total territorial area
	18. Improved sanitation facilities	% of population with access
	19. Urban population growth	Annual %
	20. Energy use	Kg of oil equivalent per capita
	21. CO2 emissions	Metric tons per capita
	22. Electric power consumption	KWh per capita
	23. GDP	Current US$ billions
	24. GDP growth	Annual %
	25. Inflation, GDP deflator	Annual %
	26. Exports of goods and services	% of GDP
	27. Imports of goods and services	% of GDP
	28. Gross capital formation	% of GDP
	29. Time required to start a business	Days
	30. Mobile cellular subscriptions	Subscriptions per 100 people
	31. Individuals using the Internet	% of population
	32. High-technology exports	% of manufactured exports
	33. Merchandise trade	% of GDP
	34. Net barter terms of trade index	Index 2000 = 100
	35. Personal remittances, received	Current US$ millions
	36. Foreign direct investment net inflows Current US$ millions

TABLE IV

 IV final regression was run with this final subset and we obtained an adjusted R 2 equal to 0.8915 and a Fisher's p-value of 8.296e-06, which are both satisfactory. The results of this third regression are shown in TableV.

		.
	RESULTS OF THE SECOND REGRESSION
		Estimate	Std. Error	Pr(>|t|)
	(Intercept)	2.748e+05	3.516e+05	0.454586
	Population	1.238e+02	7.877e+01	0.150571
	Population growth	1.206e+05	5.276e+04	0.048152
	Surface area	-1.743e+01	5.298e+00	0.009384
	Population density	-8.596e+02	2.186e+02	0.003446
	Fertility rate	-1.518e+05	6.972e+04	0.057413
	Adolescent fertility	-4.945e+03	1.811e+03	0.023218
	Immunization measles	6.920e+03	2.782e+03	0.034541
	Terrestrial and marine protected -5.989e+03	2.119e+03	0.019865
	GDP	2.435e-02	4.734e-03	0.000608
	Imports of goods	-6.371e+03	2.786e+03	0.048002
	Gross capital formation	-4.446e+03	3.435e+03	0.227746
	Mobile cellphone subscriptions	9.998e+02	4.914e+02	0.072391
	Individuals using the internet	-6.251e+03	2.643e+03	0.042264
	High technology exports	1.716e+04	4.378e+03	0.003515

A

TABLE V

 V

		.
	RESULTS OF THE FINAL REGRESSION
		Estimate	Std. Error	Pr(>|t|)
	(Intercept)	1.138e+05	2.329e+05	0.63382
	Population growth	8.968e+04	2.904e+04	0.00940
	Surface area	-1.229e+01	4.957e+00	0.02897
	Population density	-6.228e+02	1.825e+02	0.00516
	Fertility rate	-1.087e+05	4.149e+04	0.02234
	Adolescent fertility rate	-4.062e+03	1.158e+03	0.00432
	Immunization measles	7.285e+03	2.577e+03	0.01524
	Terrestrial and marine protected -4.148e+03	1.744e+03	0.03489
	GDP	2.889e-02	4.506e-03	3.35e-05
	Imports of goods	-5.082e+03	2.193e+03	0.03895
	Individuals using the internet	-5.770e+03	1.308e+03	0.00085
	High technology exports	1.266e+04	3.060e+03	0.00138

ACKNOWLEDGMENT

This work is supported by the Urban Logistics Chair at MINES ParisTech, sponsored by ADEME (Agency for the Environment and Energy Management), La Poste Group, Mairie de Paris (Paris City Hall), Pomona Group and RENAULT.


Jonàs Martínez

Samuel Hornus

Haichuan Song

Sylvain Lefebvre

Polyhedral Voronoi diagrams for additive manufacturing

Keywords: CCS Concepts:, Computing methodologies → Shape modeling; Voronoi diagram, 3D printing, additive manufacturing

Density

Anisotropy & Angle D e n si ty Euclidean Our: polygonal Fig. 1. Left: We propose a novel type of microstructures well suited for fabrication on continuous deposition processes such as fused filament fabrication. Our microstructures afford for a wide range of elastic behaviors, from isotropic to orthotropic. As illustrated their geometry can be spatially graded. Right: We model the foams as Voronoi diagrams of a special kind, which contrary to Euclidean Voronoi diagrams produce geometry enforcing all fabrication constraints.

A critical advantage of additive manufacturing is its ability to fabricate complex small-scale structures. These microstructures can be understood as a metamaterial: they exist at a much smaller scale than the volume they fill, and are collectively responsible for an average elastic behavior different from that of the base printing material making the fabricated object lighter and/or flexible along specific directions. In addition, the average behavior can be graded spatially by progressively modifying the microstructure geometry.

The definition of a microstructure is a careful trade-off between the geometric requirements of manufacturing and the properties one seeks to obtain within a shape: in our case a wide range of elastic behaviors. Most existing microstructures are designed for stereolithography (SLA) and laser sintering (SLS) processes. The requirements are however different than those of continuous deposition systems such as fused filament fabrication (FFF), for which there is currently a lack of microstructures enabling graded elastic behaviors.

In this work we introduce a novel type of microstructures that strictly enforce all the requirements of FFF-like processes: continuity, self-support and overhang angles. They offer a range of orthotropic elastic responses that can be graded spatially. This allows to fabricate parts usually reserved to the most advanced technologies on widely available inexpensive printers that also benefit from a continuously expanding range of materials.

INTRODUCTION

Novel capabilities of additive manufacturing processes are challenging the traditional separation between shape and material. Tool accessibility and mold extraction constraints no longer apply, and shapes with intricate internal details can be fabricated. These internal microstructures modify the large-scale material properties, making parts lighter, porous or flexible and resilient. In addition, these microstructures can be spatially graded: the same object may be rigid where external stresses apply, and lighter in other regions.

This has spawned an intense research effort towards the geometric design of microstructures that can be embedded into volumes while triggering varying properties. The microstructures are typically considered as a metamaterial: at large-scale they behave equivalently to a uniform material characterized by a homogeneous elasticity tensor. This average behavior can be analyzed, abstracting away the fine-scale details, in a process called numerical homogenization. In addition, the microstructures are often parametric: the geometry of their fine scale details is controlled through high-level parameters, such as thickness or orientation. These parameters impact the large-scale elastic behavior, and thus the observed elasticity tensor can be controlled by manipulating these parameters directly. The main interest in using such parametric metamaterials is that their elasticity can be spatially controlled -graded -by a control field, varying mechanical properties within the shape volume. The control field can be either manually painted through a dedicated user interface [START_REF] Ion | Metamaterial Mechanisms[END_REF] or automatically computed through, e.g. topology optimization [START_REF] Zhu | Two-Scale Topology Optimization with Microstructures[END_REF].

In this paper we focus on defining a novel class of parametric metamaterials. A key challenge -besides triggering the desired elastic behavior -is to enforce constraints due to the additive manufacturing processes. Most existing works on the topic have focused on processes such as selective laser sintering (SLS) and stereolithography (SLA). These technologies have different constraints than the widely available and inexpensive fused filament fabrication (FFF) 3D printers. In particular, most recent works produce geometries exhibiting truss structures (detailed in Section 2). After slicing, these structures produce many isolated ellipses in each slice. These are extremely challenging to print on FFF systems, where continuity of deposition, low overhang angles, and support from below are strict requirements (for examples of actual trials we recommend the online article by [START_REF] Zhou | A study in fabricating microstructures[END_REF]). This is especially unfortunate, as FFF is not only widely available but allows to fabricate in a wide spectrum of materials, from inexpensive plastics, food-safe materials, nylons, to metal-polymer compounds and carbon fiber reinforced materials (see filament by e.g. Proto-Plasta). In addition, other technologies such as wire arc metal printing1 and contour crafting2 share the same constraints. These constraints are very restrictive and have so far prevented the development of metamaterials for these processes.

Contributions.

Our paper aims at filling this gap, defining a novel parametric metamaterial with the following properties:

• closed-cell, stochastic foam-like geometry readily printable on FFF-like systems, with strictly enforced overhang angles and that is everywhere supported from below, • density and anisometry controllable through user provided fields that can vary arbitrarily, • defines orthotropic materials that can be freely oriented in the plane orthogonal to the fabrication direction, • a purely procedural computation that scales to arbitrarily large shapes. This is made possible by defining procedural Voronoi diagrams based on polyhedral distances. We identify a class of parameterized polyhedral distances providing guarantees on the maximum overhang angles, and ensuring that all deposited material is supported from below. These constraints remain enforced when varying density and orientation parameters. We detail the procedural generation of our foams, analyze their elastic behavior and the link between the underlying parameters and the obtained elasticity tensors, and produce a variety of results demonstrating applications.

RELATED WORK

Researchers in additive manufacturing have recognized the need for internal fill patterns early on [START_REF] Mcmains | Layered Manufacturing of Thin-Walled Parts[END_REF]]. Indeed, by avoiding densely filling an interior, the print time and material use are greatly reduced. Thus, early infill patterns were essentially targeted at reducing material density, while being fast to generate and efficient to fabricate. Typical patterns for continuous deposition encompass parallel lines with varying spacing, but more elaborate patterns have been devised in an effort to achieve better strength to weight ratios [START_REF] Livesu | From 3D Models to 3D Prints: An Overview of the Processing Pipeline[END_REF]].

In recent years, the focus has shifted to producing infill patterns with controllable mechanical properties. Most approaches cast the problem as a two-scale optimization [START_REF] Zhu | Two-Scale Topology Optimization with Microstructures[END_REF]. A parametric metamaterial is defined either as a set of (parametric) periodic microstructures [START_REF] Panetta | Worst-case Stress Relief for Microstructures[END_REF][START_REF] Panetta | Elastic Textures for Additive Fabrication[END_REF][START_REF] Schumacher | Microstructures to Control Elasticity in 3D Printing[END_REF] or as a random process producing geometry [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF][START_REF] Martínez | Orthotropic k-nearest foams for additive manufacturing[END_REF]]. The link between microstructure parameters and their average elastic behavior is analyzed, defining a material space. Each entry of the material space relates a choice of structure parameters to an elasticity tensor. A user or an algorithm can directly work at a coarser scale and specify the desired elastic behavior in different regions of space. This is translated into microstructure parameters used to produce a final geometry. The definition of the geometry is often procedural, such that it is instantiated only when sent to the printer, slice by slice [START_REF] Pasko | Procedural function-based modelling of volumetric microstructures[END_REF][START_REF] Vidimče | OpenFab: A Programmable Pipeline for Multi-material Fabrication[END_REF]. Our work is inscribed in this direction of research.

Most results in this area produce microstructures that are fabricable only with SLS/SLA, which have constraints different from FFF. The challenges in designing infill patterns for continuous deposition are threefold. First, material can only be deposited on top of already solidified layers: it is not possible to start a new feature in 'mid-air'. Similarly, while slanted walls can be fabricated, there is a limit to the maximum overhang angle (typically around 45 degrees, varies with layer and deposition thicknesses). Second, the extrusion process should be interrupted as little as possible, as each start/stop produces small deposition defects. Third, the deposition head is typically a relatively heavy device, and acceleration rates are limited by mechanical capabilities.

The first constraint (overhangs, support) relates to the geometry of the patterns and requires defining self-supporting geometry (e.g. not requiring auxiliary support for fabrication). The two other constraints (continuity, acceleration) relate to deposition efficiency and quality. In particular, truss structures which print well on SLA/SLS are ill-suited for continuous deposition due to the last two constraints: they require printing many small discs spread throughout each layer. This triggers many start/stop and abrupt speed changes -the low acceleration limiting the maximum reachable speed. As a consequence, the beams have lower quality and print slower than their volumes would indicate. Unsurprisingly, most patterns for FFF are comprised of as long as possible continuous lines.

Nevertheless, recent research has considered improved infill patterns for continuous deposition. [START_REF] Leary | Feasible Build Orientations for Self-Supporting Fused Deposition Manufacture: A Novel Approach to Space-Filling Tesselated Geometries[END_REF] advocate for the use of self-supported spatial tessellations, considering maximum overhang angles. [START_REF] Wu | Selfsupporting rhombic infill structures for additive manufacturing[END_REF] explore how sub-dividable rhombic infill patterns [START_REF] Lefebvre | 3D infilling: faster, stronger, simpler[END_REF]] can be optimized for rigidity. [START_REF] Lee | Block-based inner support structure generation algorithm for 3D printing using fused deposition modeling[END_REF] optimize similar patterns to produce parts that are as empty as possible -while producing self-supported structures. [START_REF] Hornus | Tight Printable Enclosures and Support Structures for Additive Manufacturing[END_REF], [START_REF] Xie | Support-free interior carving for 3D printing[END_REF] and [START_REF] Wang | Support-Free Hollowing[END_REF] consider how to maximally carve parts while producing self-supported inner cavities. [START_REF] Lu | Build-to-last: Strength to Weight 3D Printed Objects[END_REF] optimize the positions of the point sites of a Voronoi diagram to strengthen an object. Interestingly, a soluble support material had to be used within the cells (page 6, paragraph Physical test of [START_REF] Lu | Build-to-last: Strength to Weight 3D Printed Objects[END_REF]) -that is because usual Voronoi diagrams cannot form self-supported structures (see Section 3). [START_REF] Lee | Support-Free Hollowing for 3D Printing via Voronoi Diagram of Ellipses[END_REF] pack many ellipses within a volume interior to obtain a self-supported interior fill.

However, none of these techniques afford for the definition of a material space, in the same sense as the microstructures employed with SLS/SLA [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF][START_REF] Panetta | Elastic Textures for Additive Fabrication[END_REF][START_REF] Schumacher | Microstructures to Control Elasticity in 3D Printing[END_REF]]. One exception is the work of Martinez et al. [2017] on orthotropic foams, where orthotropic 2D patterns are fabricated on FFF printers. However, this only applies to 2D patterns extruded vertically, not 3D structures. Subdivisible infills [START_REF] Medeiros E Sá | Adaptive voids[END_REF][START_REF] Wu | Selfsupporting rhombic infill structures for additive manufacturing[END_REF]] come close, but these cannot produce smooth spatial variations in orientation and density due to their underlying regular nature. We discuss them in more details and provide elements of comparison in Section 5.5.

Our work aims at introducing infill patterns which are well suited for fabrication by continuous deposition, can be graded spatially, and provide a wide material space for design.

OVERVIEW

The microstructures we explore belong to the family of closed-cell foams. They are defined as a stochastic process, which parameters drive the small-scale material arrangement, and ultimately the average elastic behavior.

Our foams are defined as the thickening of the facets of Voronoi diagrams of a special kind. We give below some more specific background on the use of Voronoi diagrams to represent foams and microstructures, explain the challenges and give some intuition on our technique. The precise definition and analysis of our microstructures are given in Section 4.

Voronoi foams in additive manufacturing

In material science, naturally occurring cellular solids are often idealized as sets of Voronoi edges (open-cell) or facets (closed-cell) [START_REF] Lorna | Cellular solids: structure and properties[END_REF]. These have been well studied as they exhibit desirable properties such as an isotropic elastic behavior [START_REF] Mathias | Numerical simulations of 3D open cell structures -influence of structural irregularities on elastoplasticity and deformation localization[END_REF]] with a quasi-linear link between density and Young's modulus [START_REF] Lorna | Cellular solids: structure and properties[END_REF][START_REF] Roberts | Elastic properties of model random threedimensional open-cell solids[END_REF].

Voronoi open-cell foams are used for fabrication on SLA/SLS processes, as they form fully connected truss networks well suited to resin and powder solidification. They are amenable to efficient procedural synthesis [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF]]. Open-cell foams are however not very suitable for continuous deposition, due to the constraints mentioned in Section 2. For continuous deposition, closedcell foams are more promising: they form walls that would afford for continuous deposition and are also known to exhibit good mechanical behaviors -in fact, for the same volume closed-cell foams are more rigid than open-cell foams (see [Gibson and Ashby 1999, Formulas 5.6, 5.13] and [START_REF] Sigmund | On the (non-)optimality of Michell structures[END_REF]). The reason closedcell foams are not used with SLA/SLS is that they would trap non solidified material in closed cells. This, however, is not an issue on FFF-like processes.

Unfortunately, as we discuss next, the usual Voronoi diagrams do not define self-supporting closed-cell structures. Besides, their elastic behavior would remain limited to isotropic elasticity, while we seek a wider material space.

Challenges of usual Voronoi diagrams

We now discuss why usual (Euclidean) Voronoi diagrams cannot be readily used. Let us consider such a diagram defined from a set of point sites s i ∈ S. Each Euclidean Voronoi cell is defined by:

V P (s i ) = p ∈ R 3 : ||s i -p|| 2 < ||s j -p|| 2 , ∀s i s j ∈ S (1)
The closed-cell foam is formed by thickening the facets of the diagram. Figure 2 shows 2D examples where facets are edges. In these side views, the build direction goes upward. As can be seen in the Euclidean diagram (leftmost), the facets/edges can take an arbitrary orientation and many of them violate the overhang constraint: they cannot be printed on continuous deposition systems. In addition, when the point sites are uniformly distributed, the cells tend to be isotropic, which prevents the emergence of anisotropic elastic behaviors (see [START_REF] Martínez | Orthotropic k-nearest foams for additive manufacturing[END_REF]). 3.2.1 Changing the distance. Voronoi diagrams can be defined using different distances. Figure 2 shows the diagram obtained using the L 1 and L ∞ distances on the same point sites.

Euclidean ⃝ L 1 L ∞ □ Cone △
While these diagrams are not printable -due to the horizontal facets -an interesting emerging property is that the facets take only a fixed number of different angles. This, in fact, is a property of any Voronoi diagram under a polyhedral distance: a distance induced by a convex polyhedron [START_REF] Icking | A Tight Bound for the Complexity of Voronoi Diagrams Under Polyhedral Convex Distance Functions in 3D[END_REF][START_REF] Ma | Bisectors and Voronoi diagrams for convex distance functions[END_REF].

Our intuition is that some specific convex polyhedron must exist that induces polyhedral Voronoi diagrams whose facet set is fully fabricable. Figure 2 (rightmost) shows an initial test that confirmed this intuition: using a cone distance produces a fully printable diagram. This idea forms the basis of our approach.

We detail how we build upon this initial intuition to define parameterized 3D Voronoi diagrams in Section 4. We analyze the mechanical behavior of the structures in Section 5 and demonstrate applications to additive manufacturing in Section 6.

METHOD

We now describe the use of polyhedral Voronoi diagrams for producing microstructures with spatially varying elasticity, well suited for continuous deposition (FFF-like processes).

Let B d = {p : p z = d} be a plane orthogonal to the printing direction z, and B - d = {p : p z < d} be a half-space below B d . We call slope of a facet or slope of a plane the dihedral angle between the plane supporting the facet and any horizontal (xy) plane. For example a vertical wall has slope π /2, the printing plane slope zero.

Our microstructure is defined by a mesh K that must satisfy the following two fabrication constraints:

(1) Angle constraint: Each facet f of K has to be "sufficiently close to vertical" so that it is printable with FFF. We model this constraint by imposing a lower-bound θ ⋆ ∈ [0, π /2] on the slope of facet f . If n is the normal vector of facet f , it must satisfy |n z | ≤ cos θ ⋆ . Reasonable values of θ ⋆ may range, say, from π /6 to π /3 and depends on the fabrication hardware at hand. (2) No-local-minimum constraint: Each point p ∈ K has to be "supported from below" so that it does not appear mid-air during fabrication; i.e. we have that K ∩ N (p) ∩ B - These two constraints are well-suited for FFF-like processes, as they capture the requirement that deposited material has to be well supported from below. For our purpose in this paper, we model the mesh K as the union of the facets of a polyhedral Voronoi diagram, namely, the Voronoi diagram of a set of point sites that arises when using a specific polyhedral distance instead of the usual Euclidean one (Section 4.1.1). We show that some choices of distance polytopes guarantee that K satisfies both fabrication constraints. In particular, we show how one can easily check that a given polyhedral distance entails a fabricable mesh K (Section 4.1.4) and we exhibit a simple family of such polytopes (cones, Section 4.2). We discuss which parameters can be spatially graded and how (Section 4.3). Finally, in Section 4.4 we detail a simple algorithm to extract the deposition paths delineating to the Voronoi facets.

Polyhedral Voronoi diagrams

We now recall the definition of polyhedral distances, polyhedral bisectors, and their use in Voronoi diagrams. For simplicity, some of the following illustrations consider a planar polygonal distance.

4.1.1 Polyhedral distances. Let P ⊂ R 3 be a compact convex polyhedron of combinatorial size k containing the origin in its interior. The polyhedral distance induced by P from a point p to q is (e.g., [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF]):

d P (p, q) = min {t ≥ 0 : q ∈ p + tP } .
(2)

The value of d P (p, q) is the minimal scaling factor t applied to P so that p +tP touches q (Figure 3). Note that the L 1 and L ∞ distances are particular cases of polyhedral distances. Since P is convex, the triangle inequality holds: 0 ≤ d P (p, r ) ≤ d P (p, q) + d P (q, r ) (proof in supplemental material). However, d P is symmetric if and only if P is centrally symmetric with respect to the origin.

The distance d P (p, q) can be computed as follows [START_REF] Icking | Convex distance functions in 3-space are different[END_REF]. Let q ′ be the unique point of p + P intersected by the ray from p to q (see Figure 3). Then d P (p, q) = ∥q-p ∥ ∥q ′ -p ∥ . In our implementation, we compute the closest intersection of the ray with the planes supporting the facets of P. 4.1.2 Polyhedral bisectors. The bisector between two point sites p and q is the set of points equidistant to p and q under the distance d P , that is:

B(p, q) = r ∈ R 3 : d P (p, r ) = d P (q, r ) (3)
Note that the polyhedral distance is taken from p and q to r because d P is not symmetric.

4.1.3 Polyhedral Voronoi diagrams. Equipped with the polyhedral distance, we now describe their use in Voronoi diagrams. Let S be a set of n point sites in R 3 . The polyhedral Voronoi diagram Vor P (S) is the decomposition of R 3 into Voronoi cells. A cell V P (s i ) for a site s i ∈ S is defined as:

V P (s i ) = p ∈ R 3 : d P (s i , p) ≺ d P (s j , p), ∀s i s j ∈ S (4)
where the lexicographical distance relation ≺ is necessary to deal with cases where a full region space belongs to the bisector (nongeneral position cases [START_REF] Klein | Voronoi diagrams based on general metrics in the plane[END_REF]), see Figure 4:

d P (s i , p) ≺ d P (s j , p) ⇐⇒ d P (s i , p) < d P (s j , p) ∨ (d P (s i , p) = d P (s j , p) ∧ (i < j)) (5)
Fig. 4. A non-general position case: the line passing through s 1 s 2 is parallel to the line containing the facet H 1 of P . The gray region are points of the polyhedral bisector B(s 1 , s 2 ). The lexicographical distance of Equation ( 5) selects one of the boundaries of the gray region.

The Voronoi facets form the boundary set of Voronoi cells. Due to the triangle inequality, the Voronoi cells are star-shaped polyhedrons [START_REF] Chew | Voronoi Diagrams Based on Convex Distance Functions[END_REF]: for all q ∈ V P (s) the line segment from s to q is in V P (s). Regarding the combinatorial complexity, assuming the complexity k of P is constant, the worst case complexity of the Voronoi diagram is Θ(n 2 ) [START_REF] Icking | A Tight Bound for the Complexity of Voronoi Diagrams Under Polyhedral Convex Distance Functions in 3D[END_REF]

(n = |S |).
An important property for our purpose is that the facets of Vor P (S) form a subset of the arrangement of the bisectors of each pair of sites in S × S (see, e.g., [START_REF] Ma | Bisectors and Voronoi diagrams for convex distance functions[END_REF]). Thus, if all the bisectors satisfy both fabricability constraints then the facets of the Voronoi diagram form a fabricable mesh as well. We now analyse polyhedral bisectors in more detail.

Properties of polyhedral bisectors.

In this section, we recall known facts about polyhedral bisectors and explain how these properties can be used to check whether the geometry of P induces fabricable bisectors. Various works have already studied polyhedral bisectors [START_REF] Corbalan | Geometry of bisectors for strictly convex distances[END_REF]Icking et al. 2001;[START_REF] Martini | The geometry of Minkowski spaces -A survey[END_REF]. For a more detailed exposition, we refer the reader to the thesis of [START_REF] Ma | Bisectors and Voronoi diagrams for convex distance functions[END_REF] and the references therein.

For simplicity, we assume that two points p and q are in general position with respect to P. That is, the line through p and q is not parallel to any line segment contained in the boundary of P [START_REF] Boissonnat | Voronoi Diagrams in Higher Dimensions Under Certain Polyhedral Distance Functions[END_REF]. Then, B(p, q) is piecewise linear and homeomorphic to a plane [START_REF] Ma | Bisectors and Voronoi diagrams for convex distance functions[END_REF]]. We give a method to construct the bisector of two point sites in the supplemental material.

We define a feature as either a vertex, a (straight) edge or a (planar) facet. We say that two features X ⊂ P and Y ⊂ P on the boundary of the distance polytope P generate a feature B of the bisector B(p, q) when B = λ>0 ((p + λX ) ∩ (q + λY )) . B is either empty or does not depend on the relative position of the point sites p and q. It follows from the constructions of polyhedral bisector (detailed in the supplemental material) and Proposition 4.1 below, that each feature of a bisector is generated by two features on P.

To determine fabricability, we have to analyse the slope of the bisector facets. This builds upon the following property. Let B i j be a facet of the bisector B(p, q) generated by the facets f i and f j of P. Let n i be the outward normal vector to facet f i . Let α i > 0 be the distance from the origin to the plane

H i supporting f i . (In particular, if r is a point of f i then α i = r • n i .)
Proposition 4.1. Write H i j for the plane that spans the origin and the line H i ∩ H j . Then, the bisector facet B i j (generated by f i and f j ) is parallel to the plane H i j (Figure 5). The normal vector of B i j is proportional to α j n i -α i n j .

The proof is in [START_REF] Icking | A Tight Bound for the Complexity of Voronoi Diagrams Under Polyhedral Convex Distance Functions in 3D[END_REF]. The derivation of the formula for the normal vector can be found in the supplemental material. Given Proposition 4.1, we now explain how to check if both fabrication constraints are respected.

Verifying the angle constraint. Each planar facet B i j of the bisector B(p, q) is defined through the interaction between two facets f i and f j of P (see Figure 5). Proposition 4.1 tells us that the normal vector of B i j depends only on f i and f j and is independent of the relative position of p and q. Thus, there are at most k 2 different orientations for the bisector facets. We can check that they satisfy the angle constraint by enumerating all pairs of facets of P.

Verifying the no-local-minimum constraint. Similarly, the vertices of all bisectors exhibit a finite number of configurations that can be enumerated in time O(k 2 ), to verify that no vertex of any bisector is a local minimum with respect to the build direction. Each bisector vertex is generated by either an edge-edge pair or a vertex-facet pair (see Figure 6). Details are given in the supplemental material. In summary, for a given polyhedral distance d P induced by a polytope P of size k, we can check in time O(k 2 ) if the polyhedral bisectors under d P satisfy both fabrication constraints. Fig. 5. Let f i , f j be facet of P . Let H i and H j be the supporting planes of f i and f j , respectively. The feature of B(p, q) generated by f i and f j is drawn bold and green. A construction for its supporting plane is shown. Lines with the same color are parallels. The bold green bisector feature is indeed parallel to the top green line, as constructed in Proposition 4.1. See Lemma 3.2.1.2 in [START_REF] Ma | Bisectors and Voronoi diagrams for convex distance functions[END_REF]]. Fig. 6. Polygonal bisector between two point sites p and q. Each facet B i j of the bisector is generated by facets f i on p + P and f j on q + P . The dark-red bisector vertex is generated by the two light-red features of P . The dark-green bisector edge is generated by the two light-green features of P . A section in the supplemental material details this geometric construction.

Polyhedral cone distances

We now describe our choice of a particular family of polyhedral distances that induce fabricable Voronoi diagrams. The ideal family has few parameters and the influence of each parameter should be well understood, and easy to relate to the final elastic behavior. We prefer a distance that is symmetric around the z-axis so that it is easier to reason about. The most difficult property to achieve is the absence of local minimum.

We determined that a polyhedral cone achieves these objectives, and we now describe our family of cone distances. Given the parameterized cone geometry (Section 4.2.1), we explain how the fabrication constraints are satisfied (Section 4.2.2) and extend the family with anisotropic cones (Section 4.2.3). We highlight the influence of the parameters on the induced Voronoi diagram in Section 4.3.

4.2.1

Geometry of the cone. The base facet of the cone is a horizontal regular k-gon (Figure 7-bottom-left). The cone has k side facets (typically k = 8). The scale of the distance polytope does not affect the Voronoi diagram that it induces, so we assume that the base facet lies at z = -1. The cone has apex A = (0, 0, A z ) with A z > 0. (Figure 7-right). The vertices of the base facet have y, -x-coordinates: r cos a i and r sin a i where a i = (2i + 1) π k for i = 0..k -1, and r is the radius.

Size: k Slopes: θ , µ Anisotropy: σ , ζ 3 ≤ k θ ⋆ ≤ θ < π /2 0 < µ < 1 0 < σ < 1 ζ ∈ R O A θ α 1 A z k = 8 y x r = tan(π/2 -θ) = r cos(π/k) φ α = cos φ Az 1+Az = µ z x
We parameterize the cone with parameters k, θ , µ (Figure 7-top) and add parameters σ , ζ for controlling anisotropy (Section 4.2.3). All other quantities are determined from these parameters to enforce fabrication constraints (Section 4.2.2). Examples of cones within our family are shown in Figure 8.

4.2.2

Cones for fabricable Voronoi diagrams. We now consider constraints on the slopes of the bisector features. All the side facets contain the apex A of the cone, thus if H i and H j are planes supporting two side facets f i and f j , their intersection line contains A as well. By Proposition 4.1, H i j contains both A and the origin O, so that it is a vertical plane. We conclude that any bisector feature generated by two side facets is vertical and therefore satisfies the angle constraint.

Let us now consider a feature generated by a side facet and the base facet. By Proposition 4.1, this feature is parallel to a plane spanning an edge of the base facet and the origin. Its slope θ (Figure 7-bottom) is the same for all side facets. In order to satisfy the angle constraint, we choose the value of θ in the range [θ ⋆ , π /2), where θ ⋆ is the angle constraint bound. This fixes the value of ℓ = tan(π /2 -θ ) and the radius r = ℓ/cos(π /k) of the base facet.

The last value to determine is A z . It comes from our parameter µ, which controls the distance α from the origin to a side facet. By expressing the sine of angle ϕ, we compute α = ℓ cos ϕ A z 1+A z . In order to decouple parameters θ and µ, we interpret our parameter µ as µ = α/ℓ. Then, by expressing the squared tangent of angle ϕ

as α 2 A 2 z -α 2 = ℓ 2 (1+A z ) 2 , we obtain A z = µ 1-µ 2 µ + 1 + ℓ 2 (1 -µ 2 ) .
Figure 9 shows a 2D example cone (a triangle). The table of possible bisector vertices shows the impossibility of a local minimum. We provide in the supplementary material an interactive applet illustrating the induced bisectors in 3D.

Anisotropic cones.

We vary the anisotropy of the polyhedral Voronoi diagram by shrinking the distance cone along the x axis followed by rotating it around the z axis. We do not allow other axes of rotation so as to guarantee that the Voronoi mesh stays fabricable: the base of the cone becomes a polygonal approximation of an ellipse and stays horizontal. The shrink is obtained by applying to the cone the linear transformation diag(σ , 1, 1) where 0 < σ < 1. 3 The rotation angle is specified with parameter ζ . The analysis in Section 4.2.2 still applies. A pair of side facets generate a vertical bisector feature and the base facet generates features with slopes that vary for each side facet, but are always larger than or equal to θ . Figure 8 illustrates how the parameterization impacts the corresponding polyhedral Voronoi facets.

Spatial variations

The capacity to obtain gradients of properties within fabricated shapes is a major advantage offered by additive manufacturing. In this section we discuss how the geometric properties of our microstructures can vary spatially.

4.3.1

Varying the point sites density. The microstructures that we fabricate are materializations of the Voronoi mesh K. We control their density by changing the set of point sites S on which the Voronoi diagram is computed, following an input density field R 3 → D where D ⊂ R is a range of fabricable densities. In practice we follow the implementation of [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF] to produce point sites in a multi-resolution jittered grid.

4.3.2

Varying the polyhedral distance. We now consider that each of the parameters θ , µ, σ and ζ are fields in R 3 . The value of the parameter is picked at the location of a point site s ∈ S, and used to define the distance polyhedron at s. Thus, we are effectively using a different distance for each site in S. We discuss two subsets of parameters that can be spatially varied together while still guaranteeing fabricability in Section 4.3.3 and Section 4.3.4. We always consider k ≥ 3 and θ ≥ θ ⋆ constant.

Remark. When changing the distance at each site, the fact that the bisector between two sites p and q is homeomorphic to a plane does not hold anymore. In fact, in such a situation, the bisector can have more that one connected component [Icking et al. 2001]. We never observed this behavior in our experiments. This is due to the fact that point sites are distributed all over space so that pieces of bisectors that may appear far from their two parent sites are "taken over" by the Voronoi cells of other sites.

Varying µ.

The apex of the cone moves vertically from just above the origin when µ ≈ 0 to infinitely high as µ approaches 1. Consider two side facets of two identical cones and start raising the apex of one cone (but not the other). Then, the feature that those facets generate becomes increasingly horizontal. The slope decreases but fortunately never reaches below θ .

While we do not have a formal proof, we can verify it numerically. Consider a side facet f i of a cone with a value µ < 1 together with each side facet f j of a cone with µ = 1, i.e., f j is vertical. Each pair generates the lowest possible angle for this "side j" of the cone.

k = 4 θ = 25 • σ = 0.4 µ = 0.2 µ = 0.8 σ = 0.7 θ = 65 • k = 16
For each pair, the construction of H i j from Proposition 4.1 lets us compute this lowest slope. We did check numerically that it is indeed always greater than θ .

Changes in the Voronoi mesh. A 2D illustration of how the Voronoi diagram changes is shown in Figure 10. Varying the parameter µ modifies the distribution of the slopes of the facets of the Voronoi diagram, from "all facets have slope θ " to "most facets are vertical. " Indeed, when µ = 0 every bisector features are generated by a base facet and have slope θ . When µ is close to 1 the situation is exactly the opposite, most bisector features are generated by side facets resulting in vertical walls.

4.3.4

Varying σ and ζ . The parameters σ and ζ , that control anisotropy and rotation, can vary spatially together; the fabricability is maintained. Indeed, other parameters being fixed, the cones in this 2-parameters family all share the same apex point. As we mentioned before, this implies that a bisector feature generated by side facets of two different cones remains vertical. The features generated by a base facet are the same as in Section 4.2.3. Changes in the Voronoi mesh. Geometrically, as σ decreases, the Voronoi cells become flatter in a direction specified by ζ . This flattening also increases the slope of some bisectors (as we have seen in Section 4.2.3), especially those parallel to the anisotropy direction. Compare, for example, Figure 8 center and top-right. 4.3.5 Varying σ , ζ , µ and θ . It is possible to vary all four parameters simultaneously. In that case however, µ and θ are interdependent, and function as a single parameter. The supplemental material gives more details. We have not experimented with this 3-parameters family yet, and set it as future work.

Algorithm and implementation

For fabrication we seek to efficiently generate single slices of the microstructures. In order to print with the smallest available thickness, we produce paths delineating the Voronoi facets -these paths have no thickness, since the thickening comes from the physical material deposition. To guarantee scalability, we devise a procedural approach relying on discretization, as illustrated in Figure 11. We overlay a 2D regular grid over the slice, called the labels grid and label each grid square with the site of S closest to the square center, as measured by the polyhedral distance. Then, we extract all the boundary edges between adjacent grid squares with different labels and simplify the extracted jagged paths in order to produce the final deposition paths for fabricating the slice of microstructure. This affords for a simple space-tiling approach, limiting required memory and enabling parallelism.

There are many interesting details on how to implement this efficiently, and we refer the interested reader to the supplemental material. Our implementation uses OpenCL for parallel grid labeling, and C++ for path extraction. On a GeForce GTX 1080 and an Intel i7-6800K, the V oronoi labeling for 1000 2 labels takes around 20 ms and the path extraction 1600 ms (the cone parameters are θ = 45 • , k = 8, µ = 0.6, σ = 1). The current bottleneck lies in the path extraction since we use an unoptimized sequential method.

ELASTIC BEHAVIOR

We analyze the elasticity of polyhedral Voronoi diagrams with numerical periodic homogenization akin to [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF][START_REF] Panetta | Elastic Textures for Additive Fabrication[END_REF][START_REF] Schumacher | Microstructures to Control Elasticity in 3D Printing[END_REF]]. Periodic homogenization seeks to find the elasticity tensor characterizing a periodic composite material defined from a unit periodic cell. We do not reproduce here all the details and refer the interested reader to the aforementioned publications and our supplemental material. One difference with prior works is that we use the CrAFT software [START_REF] Boittin | CraFT[END_REF], based on the fast Fourier transform (FFT) [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]]. In comparison to finite element methods, it allows us to homogenize high resolution volumes, achieving higher precision. CrAFT has been successfully used to homogenize porous microstructures such as Voronoi foams [START_REF] Boittin | Numerical simulations and modeling of the effective plastic flow surface of a biporous material with pressurized intergranular voids[END_REF]].

After analysis, we consider the parameters of the fitted orthotropic material, primarily the Young's moduli E x , E y , E z along respectively the x,y and z axis (measure of stiffness). We provide analysis of the shear moduli G and Poisson's ratio v in supplemental material, as well as the homogenized tensors obtained by this procedure.

We are also interested in evaluating the overall compressive strength. To do so, we consider the bulk modulus K, that characterizes the tendency of a material to deform in all directions when uniformly loaded in all directions. A high bulk modulus implies higher incompressibility. We consider the Voigt-Reuss average K VR of the bulk modulus bounds [START_REF] Hill | The elastic behaviour of a crystalline aggregate[END_REF]] (details in supplemental material).

All analysis is performed with a base isotropic material having Young's modulus E = 1 (normalized) and Poisson's ratio v = 0.3. Thus, it is expected that E i ∈ [0, 1]. The Poisson's ratio of an orthotropic material is unbounded [START_REF] Ting | Poisson's ratio for anisotropic elastic materials can have no bounds[END_REF]. The CrAFT convergence threshold is set such that the modulus of the divergence of the stress field is lower than 0.005, providing convergence and a reasonable trade-off between accuracy and performance.

Material space exploration

We explore the material space spanned by the parameterized cone distance (Section 4.2). For all tests the unit periodic cell is discretized into 200 3 voxels and the thickness of Voronoi facets is 0.03.

We first consider a point density ρ of 125 point sites per unit cell (given by a 5 × 5 × 5 jittered grid point distribution). We evaluate the Cartesian product of the following cone parameters, for 5 different random realizations of point sites, performing in total 2205 tests:

• θ ∈ {25 • , 45 • , 65 • } • k ∈ {4, 8, 16}
• µ ∈ {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8} • σ ∈ {1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4}

The material space we explored spans a range of relative volumes V ∈ [0.413, 0.634] for unit cells. The results for the Young's moduli are shown in Figure 12, offering a wide range of possibilities along each orthotropy axis. The results for Poisson's ratio and shear moduli are shown in the supplemental material, spanning a variable range of elastic behaviors, with an always positive Poisson's ratio.

We observe a good correlation between the cone parameters and the orthotropic elastic constants, as detailed in Section 5.2. The variance of the logarithmic Euclidean distance [START_REF] Moakher | The Closest Elastic Tensor of Arbitrary Symmetry to an Elasticity Tensor of Lower Symmetry[END_REF] between each 5 random realizations of points sites having the same parameters remains low at 0.0006, indicating that randomness has a limited impact. In addition, we verified that zero entries of the ideal elasticity orthotropic tensor differ by at most 0.005 from the homogenized one, indicating that an orthotropic material constitutes a plausible approximation.

Influence of parameters

We now consider how each parameter of the polyhedral distance impacts the mechanical response.

Influence of θ . Increasing θ increases further E z , since the Voronoi facets become closer to being vertical (see Figure 8).

Influence of µ. Increasing µ mainly increases E z , since the area of vertical Voronoi facets increases (see Figure 8). In addition, increasing µ mainly decreases v yz , increases v zx , and decreases G yz , G zx in a different degree depending on θ (see supplemental material).

Influence of σ . Decreasing σ decreases E x and increases E y , E z , since the Voronoi facets become comparatively more elongated in the y and z direction than in the x direction (see Figure 8). In addition, decreasing σ mainly increases v zx , decreases v xy , and increases G xy (see supplemental material).

Influence of k. We observe that the value of k impacts on the following two properties:

• Increasing k decreases the variance of the homogenized elasticity under different random distributions of the point sites.

For each 5 random realizations with same parameters, the average logarithmic Euclidean distance between any pair of tensors is 11% lower for k = 8 (similar for k = 16) with respect to k = 4. • Increasing k increases the Voigt-Reuss average bulk modulus.

In particular, for k = 4 the average unit cell volume is V = 0.481 and K VR = 0.107, for k = 8 we have V = 0.469 and K VR = 0.11. Thus, k = 8 exhibits a higher average K VR (similar for k = 16) while having an even lower V.

Thus, a value of k ≥ 8 is encouraged to achieve better strength to weight ratio, and being closer to the homogenized analysis.

Influence of point density. Increasing the point density ρ increases the overall elastic moduli, since the unit cell volume V increases [START_REF] Lorna | Cellular solids: structure and properties[END_REF]. This is illustrated in Figure 13, where we evaluate additional point densities ρ.

Material space control

An orthotropic material is defined by nine independent constants (see supplemental material). For a fixed θ , k, and ζ , we have three independent parameters (ρ, µ, σ ) controlling the microstructure geometry. Thus, our microstructure only spans a subset of the space of orthotropic materials.

In the context of material design, one usually seeks to independently control the three orthogonal Young's moduli E x , E y , E z . In order to unveil the correlation of the different microstructure parameters and the other elastic constants we propose to visualize each as a function of E x , E y , E z . Results are reported in Figure 14.

To provide a dense visualization from the discrete sampling of the material space, we interpolate the space of variables (see Section 5.1). We seek to approximate a function f (E x , E y , E z ) = (ρ, µ, σ ). We perform a multivariate interpolation of the function f implemented by the Python package scipy.interpolate. Other triplets of physical parameters could be chosen to perform this analysis. However, the space of Young's moduli provides the largest design freedom.

Design freedom. From the designer perspective the material space offers the following controls. The Young's moduli may be freely changed along each axis of the orthotropic frame. The z axis of the orthotropic frame has to remain aligned with the build direction. The other axes may otherwise be freely rotated (parameter ζ ). These variations may be spatially graded. Currently we limit gradations to either µ (Section 4.3.3) or σ and ζ (Section 4.3.4), but full variations are possible (Section 4.3.5). We show in Figure 15 the effect of fixing either µ or σ and varying the other. Fixing one spans a subset curve within E x , E y , E z .

Experimental verification

We performed experimental verification of the predicted elastic behavior. We print three samples with varying orthotropy and measure their linear elastic response (Young's modulus) using a compression test along each axis. The results are reported in Figure 17. To mitigate the impact of the choice of printer, we use three different machines (see Figure for details).

The experimental results are consistent with numerical results, even though there are differences. The main one is a lower stiffness along the build direction, on average. This is to be expected since layered materials are anisotropic. As future work, we would like to investigate homogenization with an anisotropic material for AM [START_REF] Liu | Homogenization of material properties in additively manufactured structures[END_REF]. Nevertheless, the discrepancies are limited and we observe a direct correlation between control parameters (µ, σ ) and the measured Young's moduli. We provide curves for one sample in Figure 18, revealing the non-linear behavior under large compressions. The curve reveals a smooth response with the expected flattening. We also performed a repeatability test, taking five measurements along each axis of a same sample. We observe variations of ± 3% without any noticeable trend.

Comparison to tessellations by planes

Self-supported structures made of interleaved sets of parallel planes have become a popular choice for filling volumes in the context of FFF [START_REF] Lee | Block-based inner support structure generation algorithm for 3D printing using fused deposition modeling[END_REF][START_REF] Lefebvre | 3D infilling: faster, stronger, simpler[END_REF][START_REF] Wu | Selfsupporting rhombic infill structures for additive manufacturing[END_REF]]. 4 To the best of our knowledge, the material space spawned by these infills has ) and experimental compression tests performed on an Instron 3345 testing machine. Experimental Young's moduli are given in megapascals (MPa). For all three tests we set k = 8 and θ = 45 • , and we vary µ and σ . We have printed each sample (30 mm cubes) on three different printers (Ultimaker 3, and Ultimaker 2, and Prusa i3) using different rolls of filaments (Ninjatek Semiflex) in diameters of either 1.75 mm (Prusai3) or 2.85 mm (Ultimakers), with the same layer thickness (0.2 mm) and print speed (15 mm/sec). We apply a compression test to each direction, to obtain E x , E y , and E z . For comparison we compute ratios to the average value across all nine tests of the same printer (i.e. how much more rigid/flexible each is with respect to the average). As can be seen the ratios agree well overall. We provide in the supplemental material the full expression of the homogenized elasticity tensor of each compression test.

V = 0.22 V = 0.27 V = 0.25 E x E y E z E x E y E z E x E y E z
not been studied in the context of microstructures. Nevertheless, we provide here some elements of comparison.

A first important limitation of rhombic structures is that they cannot be graded spatially in orientation or anisotropy, due to their regular structure. Therefore, they offer significantly less control than our approach. Their density can be varied by subdivision, for x y z Fig. 18. Compression curve for the middle sample from Figure 17, printed on Prusai3. The material enters into the non-linear regime for the y and z axes when the load approaches 300N , it remains linear for the x axis (low stiffness). We stopped the test for y as the response flattened.

instance maximizing rigidity [START_REF] Wu | Selfsupporting rhombic infill structures for additive manufacturing[END_REF]]. However, as illustrated in Figure 25, grading by subdivision introduces sudden jumps that impair the smoothness of the final elastic behavior.

Figure 16 compares the bulk modulus of the structures, for the two typical cases of using three or fours sets of parallel planes. Rhombic structures are 20% stiffer for a same volume, but this advantage vanishes if orthotropy is desired: our structures become significantly stiffer in the preferred directions for a same volume.

The rhombic structures, thanks to the use of straight lines, print slightly faster than our structures (14% difference on a Prusai3 on default settings, may vary with acceleration settings).

APPLICATIONS

In this section we present several designs modeled with our technique. The control fields are created either through a dedicated painting interface, or procedurally (we use a domain specific language for this purpose).

The designs we present are illustrations, demonstrating the versatility of our technique. To reveal the foam we print them without an external perimeter -this means that some small unsupported regions on the sides exist. In practice we observe no detrimental impact. Real designs could add an external thin skin where needed.

A striking difference between our designs and those produced with SLA/SLS is the ease of production -thanks to the use of FFF printers -as well as the robustness of the final part. In particular, parts printed with filaments having varying degrees of elasticity (Ninjatek SemiFlex and Cheetah) can flex very significantly without breaking. In our experience, foams are much more difficult to Pictures, bottom row: The left side bends easily while the right side remains rigid. On the contrary, the right side can be compressed sideways (bottom image), while the left side is much more rigid (top image). Fig. 20. Cylinder with a varying direction of orthotropy (from 0 to 90 degrees between the handles), vertically rigid but allowing rotation in one direction.

produce on SLA printers, require a careful choice of printer and parameters [START_REF] Zhou | A study in fabricating microstructures[END_REF]], and are more fragile after curing.

Our first example is a bar with varying orthotropy, shown in Figure 19. An abrupt change in orthotropy direction in the middle of the bar triggers very different behavior on both sides.

Our second example is a cylinder with a direction of orthotropy that rotates 90 degrees along the main axis, shown Figure 20. This creates an interesting case where, under rotation in the proper orientation, the cylinder twists and flatten across its section. It is otherwise rigid vertically.

Our third example is an illustration of a prosthetic finger, reproducing a result from [START_REF] Martínez | Procedural Voronoi Foams for Additive Manufacturing[END_REF]]. The two joints are obtained by a strong orthotropy aligned with the joint "axis. " Density also varies between the outside (denser) and the inside of the bend. The result offers a preferred direction of flex that corresponds better to a natural joint, in particular it is harder to bend sideways (by orthotropy) or reverse (by density). See Figure 21.

Our fourth example is a wheel that can be mounted on an actual RC car (a Monster Beetle by Tamiya), see Figure 22. The wheel prints as a single part, is lighter than the original (109 g versus 130 g), and deforms radially (shock absorbent). Of course, the design would have to be refined to act as an efficient tire. Our technique makes such applications accessible to hobbyists.

Our fifth example, Figure 23, illustrates how orthotropy can be arranged radially to achieve different properties in a pipe design. Our final example is a shoe sole with controlled elasticity, shown in Figure 24. By modifying density, angle and anisotropy we are able to control both the "softness" and "bending" of the sole. While this is a purely illustrative example, we hope our method will help produce orthopaedic footwear at much lower price points than those produced on high-end SLS/SLA printers.

LIMITATIONS, FUTURE WORK

To the best of our knowledge, the presented microstructures form the first metamaterial that is well suited for FFF and offers a wide range of controllable elastic behaviors. Their properties can be graded following custom control fields, while still satisfying strict guarantees on manufacturability.

There are however a number of limitations and areas of improvement. In terms of overall rigidity (bulk modulus), our structures are less effective than plane tessellations (see Figure 16). It remains an open question whether there exists other polyhedral distances improving bulk modulus, while enforcing fabrication constraints.

The geometry of the cone distance that we use entails a material with a orthotropy axis that can not be moved aways from the vertical axis. This lost degree of freedom may be desirable in some application. We believe further exploration of the space of polyhedral distances could let us recover this degree of freedom even while respecting the fabricability constraints of FFF-like processes.

The stochastic geometry has an impact on printing speed (about 15% slower than straight line infill on our printers). This possibly could be mitigated by optimizing the traversal of the paths; our slicer currently uses a simple "go to next closest path" heuristic.

A minor issue related to design is that the point sites are not always close to the centroid of the Voronoi cell. Therefore, at low densities and high anisotropy, there can be an offset between the control field and its actual effect of the structure. It would be interesting to compensate for this.

Finally, while we efficiently extract slices, we did not investigate interactive visualization of the microstructures -one possible approach is ray-marching. Efficient ray-traversal of the structures is an interesting topic for future work.

CONCLUSIONS

By carefully choosing the distance functions used to model Voronoi diagrams, we are able to produce stochastic closed-cell foams that enforce all manufacturability constraints on continuous deposition processes. We proposed a parameterized polyhedral distance based on cones, that provides good correlation between well understood geometric properties and observed elastic behaviors.

We envision that our method could also be helpful in optimizing periodic microstructures [START_REF] Panetta | Elastic Textures for Additive Fabrication[END_REF][START_REF] Schumacher | Microstructures to Control Elasticity in 3D Printing[END_REF] since it provides a space of geometric structures that are certain to be fabricable. It may also be interesting to use it as a replacement in approaches that optimize Voronoi diagrams for additive manufacturing, such as the work of [START_REF] Lu | Build-to-last: Strength to Weight 3D Printed Objects[END_REF].

There is no direct way to interpret our structures as open-cell foams. The Voronoi edges could be considered but they do not enforce fabrication constraints. Therefore, our technique is complementary to existing works targeting powder and resin based systems, as it answers the need for parameterized metamaterials applicable to continuous deposition.

Fig. 2 .

 2 Fig. 2. The Voronoi diagram of a set of points under various distances, with their unit ball drawn next to the name.

∅

  where N (p) is any open neighborhood of p.

Fig. 3 .

 3 Fig. 3. A polygonal distance function.

Fig. 7 .

 7 Fig. 7. Top. The parameters describing our family of distances. Bottom. The geometry of the cone for k = 8. The right triangle is the half of the cone cutaway in the plane (O, x, z). We indicate remarkable quantities used in computations.

Fig. 8 .

 8 Fig. 8. Illustration of cone parameters and effect in the corresponding polyhedral Voronoi diagram. For better understanding, the boundary of the cube is highlighted in orange color. The central sample has parameters k = 8, θ = 45 • , σ = 1.0 , and µ = 0.5. Each neighbouring sample vary a single parameter. We provide in the supplementary material illustrations for a wider variation of parameters.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Polygonal distance defined by a triangle. The three different bisector vertex configurations are shown in the table, and in context below. The polygonal bisectors do not exhibit any local minimum, and form an angle of at least 45 • with respect to the x axis.

Fig. 11 .

 11 Fig. 11. Overview of the path extraction. The input solid C is a bishop chess piece (https://www.thingiverse.com/thing:378322/). (a) The bishop cut by the plane B d . (b) We label each square in the labels grid with its closest site in S . Each Voronoi cell is shown with a different color. (c) We extract the axis-aligned edges between Voronoi cells. (d) We generate simplified deposition paths inside C.

Fig. 12 .

 12 Fig. 12. Material space for k = 8 and ρ = 125. Each column corresponds to a different minimal bisector slope θ . Each plot dot denotes a single test. The plot axes correspond the three orthogonal Young's moduli. The color of dots either depicts µ (first row) or σ (second row). For better interpretation, each dot is projected in gray color to the E z and E y planes.

Fig. 13 .

 13 Fig. 13. Material space for k = 8 and θ = 45 • . We consider same range of parameters σ and µ used in Figure 12, while varying in addition the point density for {64, 125, 512, 1728, 4096} point sites per unit cell. The data for 125 point sites per unit cell (green points) corresponds exactly to θ = 45 • in Figure 12 (middle row). For better interpretation, each dot is projected in gray color to the E z and E y planes.

Fig. 14 .σ

 14 Fig. 14. Interpolation given by the material space sampling shown in Figure 13 (θ = 45 • , k = 8). The plot axes correspond the three orthogonal Young's moduli. Both the Poisson's ratio v x y , v zx , v yz (always positive), and the shear moduli G yz , G zx , G x y (G max = E 2(1+v ) ≈ 0.38 for base material E = 1, v = 0.3) remain fixed for a desired triplet of Young's moduli. The plots unveil the good correlation between the parameters ρ, µ, σ and the Young's moduli. σ ∈ [0.4, 0.7, 1.0] µ ∈ [0.2, 0.55, 0.8]

Fig. 16 .

 16 Fig.16. We compare infills made of tessellations made of planes with polyhedral Voronoi diagrams. We set a minimum slope of θ = 45 • for all tests. We give in the supplemental material the full expression of the homogenized elasticity tensors of each test. Left table:We determine the volume V required to achieve similar K VR values (±0.001). Since polyhedral Voronoi diagrams span a wide space of material, we select the closest sample (with respect to K VR ) in the database of Figure12. However, note that other samples having similar K VR are possible. (a) and (c) correspond, respectively, to the closest match for (b) and (d). As seen in the third and sixth row, our infills require 20% more material for the same K VR . Right table: When orthotropy is desired, for instance to mimic a laminate much stiffer in y, z than in x , polyhedral Voronoi diagrams are better suited. As seen in the third row, for a lower volume (e) is significantly more stiff in y than (b), while in the fifth row (e) is significantly stiffer in both y, z than (d) for a similar volume.

Fig. 17 .

 17 Fig. 17. Comparison between numerical homogenization (Young's moduli in the normalized range [0, 1]) and experimental compression tests performed on an Instron 3345 testing machine. Experimental Young's moduli are given in megapascals (MPa). For all three tests we set k = 8 and θ = 45 • , and we vary µ and σ . We have printed each sample (30 mm cubes) on three different printers (Ultimaker 3, and Ultimaker 2, and Prusa i3) using different rolls of filaments (Ninjatek Semiflex) in diameters of either 1.75 mm (Prusai3) or 2.85 mm (Ultimakers), with the same layer thickness (0.2 mm) and print speed (15 mm/sec). We apply a compression test to each direction, to obtain E x , E y , and E z . For comparison we compute ratios to the average value across all nine tests of the same printer (i.e. how much more rigid/flexible each is with respect to the average). As can be seen the ratios agree well overall. We provide in the supplemental material the full expression of the homogenized elasticity tensor of each compression test.

Fig. 19 .

 19 Fig. 19. Top row: Orientation field ζ (left) and anisotropy field σ (right).Pictures, bottom row: The left side bends easily while the right side remains rigid. On the contrary, the right side can be compressed sideways (bottom image), while the left side is much more rigid (top image).

Fig. 21 .

 21 Fig. 21. The articulations of this finger prosthetic are obtained by spatially varying density and anisotropy. The design is very robust and easy to print: it uses a rubber-like filament (Cheetah by NinjaTek) printed on a standard FFF printer (Prusai3).

Fig. 22 .

 22 Fig. 22. Top left: Orientation field ζ (top) and density field (below). The anisotropy is constant at σ = 0.3. Bottom left: Wheel being printed on an Ultimaker 2 with Semiflex filament at 0.3 mm thickness. The print took 19 hours. Middle column: The finished wheel and its deformation under load. Right: Ready to run!

Fig. 24

 24 Fig. 24. A 3D printed shoe sole. Left: Control fields used on the model, density (top), orthotropy strength (middle) and angle (bottom). Right: Printed shoe, top, side and bending. The shoe is printed without any skin to reveal the foam structure.

Fig. 25 .

 25 Fig. 25. Compared to subdividing rhombic infill (left), Voronoi polyhedral diagrams produce a smoother grading (right) without marked transitions, triggering a smoother elastic behavior. Both structures use the same amount of material. Ours printed 14% slower due to a more complex geometry.

ACM Trans. Graph., Vol. 37, No. 4, Article 129. Publication date: August 2018.

See e.g. https://waammat.com/about/waam

See https://en.wikipedia.org/wiki/Contour c rafting

Let η be the minimum slope of the bisector facets. In order to have η stay constant when we shrink the cone along the x axis, we generate the circle in the y, -x frame so that the cone has at least one side facet perfectly parallel to the x axis.ACM Trans. Graph., Vol. 37, No.

4, Article 129. Publication date: August 2018.

See also cubic/tetrahedral infill in software Cura. ACM Trans. Graph., Vol. 37, No. 4, Article 129. Publication date: August 2018.

European Research Council

Agence nationale de la recherche

Fond européen de développement régional ACM Trans. Graph., Vol. 37, No. 4, Article 129. Publication date: August 2018.

ACKNOWLEDGMENTS

We thank Noémie Vennin and Jimmy Etienne for their help in designing and printing results. This work was supported by ERC 5 grant ShapeForge (StG-2012-307877), ANR 6 MuFFin (ANR-17-CE10-0002), Région Lorraine and FEDER. 7


Michel Calmet

Emanuela Pierantozzi

Stanislaw Sterkowicz

Monica Yuri Takito

Emerson Franchini

come   L'archive ouverte pluridisciplinaire

C h a n g e s r u l e s a n d e v o l u t i o n o f r e s u l t s i n j u d o , a n a n a l y s i s :

of the 2012 and 2016 Olympic Games and 2015 and 2017 World Championships. ) before the rule changes and after the rule changes, considering the different weight categories and competition phases. Professional magazines are interested by these results, and our objective, within the framework of the popularization of the scientific works, is to present an abstract of these two papers, providing a visual point of view of these two papers. A timeline was built to resume the different changes of the rules within the period 2009-2017, and results are presented with two tables and arrows to show evolutions. The effects between the two Olympics were negative, they became better between the two world championships. At that time, rule changes did not achieve the goal to increase scores and decrease penalties, but that this seems to be changing in the last modifications made for the 2017 World Championship. Beyond rule changes, these effects, (while the periods of these 4 championships overlap) could be explained by the frequency of the events and their own modes of selection.

Bibliography: Balafoutas, L. Lindner, F. Sutter, M. ( 2013 Rules are available immediately and can "overlap" several events (World Championships or Olympic Games) which makes the comparison difficult. Changed scores or penalties (Yuko counted in Astana in 2015 and Rio in 2016 but did not exist in Budapest in 2017), frequency of events (Olympic Games every 4 years, World Championships every year or two years) and the judokas selection mode (in or out the ranking list), brought to compare Olympic Games between them and World Championships between them.

The analyzes carried out show that rules change can:

1. Improve the situation (green inclined arrows):

• The number of games without penalties increases (green arrow rising)

• The number of penalties decreases (green arrow down) 2. Do not make any changes (orange horizontal arrows) 3. Deteriorate the situation (red arrows):

• The number of yuko decreases (red arrow down)

• The number of hansoku make increases (HSK) (rising red arrow)

Rule change and results in competition

  Changes rules and evolution of results in judo, an analysis of the 2012 and 2016 Olympic Games and 2015 and 2017 World Championships. In judo, rules have changed in recent years to improve positive judo. Results in competition have been analyzed to see if this goal is achieved. A chronological frieze facilitates how to situate the changes compared to the most important competitions.

  M. Calmet 1* , E. Pierantozzi 2 , S. Sterkowicz 3 , M. Yuri Takito 4 , E. Franchini 5 1 -Université de Montpellier, Université d'Aix-Marseille, michel.calmet@umontpellier.fr 2 -University of Genoa, Italy. emanuela.pierantozzi@unige.it 3 -Combat Sports Unit, Department of Theory of Sport and Kinesiology, University of Physical Education in Cracow, Poland stanislaw.sterkowicz@awf.krakow.pl 4 -Human Movement Pedagogy Department, Group, School of Physical Education and Sport, University of São Paulo, Brazil. mytakito@gmail.com 5 -Martial Arts and Combat Sports Research Group, School of Physical Education and Sport, University of São Paulo, Brazil. emersonfranchini@hotmail.com

	Abstract:
	In the 2009-2017 period, many changes were added in judo rules, especially after London 2012
	Olympics to promote "positive" actions (constant search for maximum score) instead of "negative"
	actions (avoiding opponents' attacks). International Judo Federation publishes data and general
	statistics (www.judobase.org). Two published papers (*) used these data to analyze the total match
	duration, frequency of scores (ippon, waza-ari, yuko) and penalties (shido) between two Olympics
	Games (London 2012 and Rio 2016) and between two World Championships (Astana 2015 and
	Budapest 2017

Olympic Games analyses: London (2012) and Rio de Janeiro (2016)

	Items	Number of ippon	Number of waza-ari	Number of yuko	Number of penalties Female	Number of penalties Male	Number of hansoku make (F and M)	Number of match without penalty	Number of golden score (F and M)	Duration of match F (5min -> M (5min -> 4min) 5min)	Number of penalties Extralight -Heavy
	Evolution
							Not studied.		Not studied.	☺ If
	Comments							☺

more ippon If more Hansoku Make No difference in categories ☺ Female receives fewer penalties than Male Learn

  more: Rule change and Olympic judo scores, penalties and match duration, Calmet M, Pierantozzi E, Sterkowicz S, Challis B, Franchini E, International Journal of Performance Analysis inSport, 2017, 17 -4, 458-465.

	World Championships analyses:
					Astana (2015) and Budapest (2017)
	Items	Number of ippon	Number of waza-ari	Number of penalties Female	Number of penalties Male	Number of hansoku make (F and M)	Number of match without penalty	Number of golden score (F and M)	Duration of match M (5min -> F (4min) 4min)	Number of penalties Extralight -Heavy
	Evolution
	Comments		☺	☺ ☺ Female receives fewer ☺ penalties than Male	Possible explanation: 3 shidos = Hansoku Make		Reduced time for M		The heavier the category, more the judokas received penalties

Learn more: Judo rules: searching for a wind of changes; Calmet M, Pierantozzi E, Sterkowicz S, Takito MY, Franchini E; International Journal of Performance Analysis in Sport 2017, https://doi.org/10.

1080/24748668.2017.1405612

Conclusions

The effects between the two Olympics were negative, they became better between the two world championships: at that time, rule changes did not achieve the goal to increase scores and decrease penalties, but that this seems to be changing in the last modifications made for the 2017 World Championship. Beyond rule changes, these effects, (while the periods of these 4 championships overlap) could be explained by the frequency of the events and their own modes of selection.