Thomas Leduc
email: thomas.leduc@crenau.archi.fr

Michel Leduc

Minimum-area ellipse bounding an isovist: towards a 2D GIS-based efficient implementation

Keywords: Minimum-area bounding ellipse, isovist, elliptic hull, caliper, diameter

In geographic information science and technology, various methods and studies exist to characterize the linearity, rectangularity, convexity, circularity or compactness, sinuosity or tortuosity of a given spatial shape. Although there is much work on ellipticity in image processing, we do not address, in geomatics, the issue of matching to a reference elliptical shape. Regarding this issue, this article is a contribution to the qualification of urban open spaces. It provides an operating algorithm for determining a minimum-area bounding ellipse for any given polygonal shape. It also proposes an implementation of this algorithm in the context of a Geographic Information System and a Jupyter Notebook. As an application, it focuses on two real urban configurations on fields of about 300 isovists. The results from the application of this approach in two urban areas in France show that the ellipse is a better minimum bounding geometry than are circle or rectangles, at least for the half-dozen descriptors studied. The improvement relatively to the minimum bounding rectangle is particularly significant in terms of correlation concerning the orientation (+20%) and drift (+10%).

Introduction

Shape in GI Science

As [START_REF] Boyce | The concept of shape in geography[END_REF] note in a seminal article, "Although shape has been primarily used in geography as a descriptive device [e.g. the 'boot' shape of Italy], it has also been of value as an analytical tool". According to [START_REF] Galton | Shape. Geographic Information Science & Technology Body of Knowledge[END_REF], the appeal for the morphology is justified "because the shape of a geographical entity can have farreaching effects on significant characteristics of that entity". He illustrates the point by referring to settlement patterns, trading relations, road network layout, etc.

"To characterize the shapes of geographical entities in ways that usefully capture such implications" [START_REF] Galton | Shape. Geographic Information Science & Technology Body of Knowledge[END_REF], three main categories of methods have been developed in the literature such as a) shape descriptors, b) some commonly used shape surrogates, and c) measures of similarity between shapes. In the first category, the objective is, for example, to measure the linearity [START_REF] Stojmenović | Measuring linearity of planar point sets[END_REF], rectangularity (Rosin 1999, Zhang and[START_REF] Zhang | Novel shape indices for vector landscape pattern analysis[END_REF], convexity [START_REF] Rahtu | A new convexity measure based on a probabilistic interpretation of images[END_REF], circularity or compactness [START_REF] Lee | A Method of Measuring Shape[END_REF][START_REF] Maceachren | Compactness of Geographic Shape: Comparison and Evaluation of Measures[END_REF][START_REF] Angel | Ten compactness properties of circles: measuring shape in geography[END_REF][START_REF] Li | An efficient measure of compactness for two-dimensional shapes and its application in regionalization problems[END_REF], sinuosity or tortuosity [START_REF] Adolphe | A simplified model of urban morphology: application to an analysis of the environmental performance of cities[END_REF][START_REF] Grisan | A novel method for the automatic grading of retinal vessel tortuosity[END_REF]) of a given shape. In the second category, regarding shape surrogates, one can mention several minimum bounding geometries, such as the convex hull, the circle, the envelope, and the oriented rectangles. Concerning the latter, [START_REF] Zhang | Novel shape indices for vector landscape pattern analysis[END_REF] introduces several variants, such as the minimal width bounding box or the moment bounding box. According to [START_REF] Galton | Shape. Geographic Information Science & Technology Body of Knowledge[END_REF], the same category also includes skeletons [START_REF] Leborgne | Noise-resistant Digital Euclidean Connected Skeleton for graph-based shape matching[END_REF][START_REF] Sarradin | Comparing sky shape skeletons for the analysis of visual dynamics along routes[END_REF]) and quadtrees [START_REF] Samet | Storing a collection of polygons using quadtrees[END_REF]. Finally, in the third category, regarding geometric shape matching or similarity measurements, various techniques exist such as those derived from signal processing [START_REF] Ai | A shape analysis and template matching of building features by the Fourier transform method[END_REF], or measurements such as Hausdorff distance [START_REF] Sarradin | Comparing sky shape skeletons for the analysis of visual dynamics along routes[END_REF].

Our paper is in line with research on two-dimensional geographic shapes and more precisely on the morphological analysis of urban open spaces. It benefits from a set of works relating to the analysis of urban layout in an immersed perspective, to visualbased studies, to field-oriented approaches, etc. Its purpose is, at every point of the open space in the city, to qualify the two-dimensional shape of the immediately surrounding space, or interstitial space that interweaves between buildings to form the city's connective network. As [START_REF] Marshall | Streets & Patterns[END_REF] points out, streetscape forms "a contiguous network or continuum by which everything is linked to everything else" and this "plenum of the urban space" [START_REF] Couclelis | People Manipulate Objects (but Cultivate Fields): Beyond the Raster-Vector Debate in GIS[END_REF], is an alternative way to consider the study of urban form. It contrasts with the standard atomic approach, which focuses on the study of perfectly identified urban feature configurations (the buildings) with clear and stable boundaries.

The isovist [START_REF] Benedikt | To take hold of space: isovists and isovist fields[END_REF], as the set of points of the horizontal plane associated with an observer -or a panoptic visual sensor -immersed in the city, is a solution, within the framework of the ontology plenum, for determining the immediate environment associated with the position. It makes it possible to describe the local geometric properties of the surrounding space from a single point of generation (the point of view) by associating it with a single polygonal shape that can itself be qualified.

Unfortunately, the isovist, which is a valuable tool to determine "'how far' and 'how much' one can see" [START_REF] Batty | Exploring Isovist Fields: Space and Shape in Architectural and Urban Morphology[END_REF], does not explicitly provide a dominant orientation of the space or a measure of its flattening or ellipticity. And, to this end, standard solution based on the minimum-area bounding rectangle does not appear to be appropriate. Its sensitivity to the configuration of the surroundings is highlighted in Figure 1, where the orientation of such a rectangle seems unsuitable. In contrast, the main axis of the ellipse of the minimum area enclosing the isovist determines a smoother overall direction, closer to the direction of the isovist diameter. This being said, the more precise purpose of this work is to address the issue of ellipticity and the matching of a polygonal shape to a reference elliptical shape. We will thus pursue the purpose developed by [START_REF] Lasserre | A generalization of Löwner-John's ellipsoid theorem[END_REF]: "'Approximating' data by relatively simple geometrical objects is a fundamental problem with many important applications and the ellipsoid of minimum volume is the most well-known of the associated computational techniques." In order to be consistent with the name "minimum-area bounding rectangle" (resp. minimum-perimeter bounding rectangle), usually abbreviated by the acronym MABR (resp. MPBR), we adopt the naming convention "minimum-area bounding ellipse" and the corresponding acronym MABE, to describe the ellipse of the minimum area enclosing a given geometric shape. From [START_REF] Toussaint | Solving Geometric Problems with the Rotating Calipers[END_REF])'s work, we could also name it the smallest-area enclosing ellipse.

Yet, as Žunić et al. (2017) point out, "a lot of work has been dedicated to solve the appearing, ellipticity1 associated, problems", in particular in "astronomy, astrophysics, nano-particles analysis, and traffic analysis". So, in terms of image processing, "an infinite family of ellipticity measures" ( Žunić et al. 2017) has been recently designed.

About ellipse and ellipticity

Standard geographic information systems or database extender (including ArcGIS Desktop Advanced (ESRI 2014), Grass (GRASS Development Team 2017), QGIS (QGIS 2020), SAGA [START_REF] Conrad | System for Automated Geoscientific Analyses (SAGA) v. 2.1.4[END_REF], and PostGIS (PostGIS 2020)) offer a set of minimum bounding geometry tools capable of extracting, for a given polygonal shape, its convex hull, its oriented minimum bounding rectangle, its minimum bounding circle, and even sometimes its concave hull. Matching to a shape surrogate such as the ellipse or identification of an elliptical minimum bounding geometry is missing in vector processing.

As exhibited in Figure 2, it seems that the ellipse has a stability that justifies tackling the elliptic hull issue in the context of vector processing. Indeed, the given hexagon under study emphasizes that the bounding rectangle with the minimum area ( 7.1 area units, 11.3 length units) is not necessarily the bounding rectangle with the minimum perimeter ( 7.9 area units, 11.2 length units) on the one hand, and that the orientations of these two rectangles differ substantially ( 50°) on the other hand.

Figure 2. Considering the plotted hexagon, its bounding ellipse with the minimum area (solid red contour) passes through the four points P 1 , P 2 , P 3 , and P 5 , whereas the bounding ellipse with the minimum perimeter (dashed red contour) passes through the four points P 1 , P 3 , P 5 , and P 6 . We also represent here the bounding rectangle with the minimum area (solid blue contour) and the bounding rectangle with the minimum perimeter (dashed blue contour).

In contrast, this figure shows that the bounding ellipse with the minimum area ( 7.0 area units, 10.0 length units) is also the "best bidder" compared to these bounding rectangles, whether to minimize the area or minimize the perimeter. Finally, it shows that, in this particular case, the bounding ellipse with the minimum perimeter ( 7.2 area units, 9.8 length units), although different from the minimum-area bounding ellipse, nevertheless remains quite close (there is less than 0.3°of angular deviation between their respective major axes). In this illustration, we have deliberately left the minimum bounding circle aside since it does not allow, by essence, to identify a major orientation.

This raises the question of whether comparison to a reference elliptical shape is an interesting option for characterizing a shape. Aware that there is no single numerical shape descriptor allowing such a characterization, we propose to focus on a set of six measures inspired in particular by [START_REF] Galton | Shape. Geographic Information Science & Technology Body of Knowledge[END_REF]. Area and perimeter are standard and rotation invariant measures, as are the diameter 2 , the stretching -which is defined as the ratio of the minimum Feret diameter3 to the maximum Feret diameter -and the drift4 . These last three descriptors, like the shape orientation5 , which is the sixth descriptor we have chosen, are scale invariant. To measure the convexity, circularity, rectangularity or ellipticity of the shape we proceed by correlating the descriptors of the surrogate shapes. The question raised here is whether the ellipse could be, in a "continuum" from the convex hull to the minimum bounding circle, a good compromise in terms of shape surrogate.

Because there is no implementation we can use6 , this article will first present a method for determining the MABE and its algorithmic implementation details in the context of a geographic information system (namely, QGIS). It will then present an application of this algorithm to the analysis of isovist fields in two real urban cases.

Materials and method

General principle

According to [START_REF] Gruber | John and Loewner Ellipsoids[END_REF], for any compact set in the plane, there exists a unique ellipse of the minimum area containing this set. It is called the Löwner-John ellipsoid of the compact set, or more simply its minimum-area bounding ellipse.

Let us first observe that five vertices of a convex polygon determine at most one ellipse passing through them. The coefficients of its Cartesian equation are easily determined by solving a five-dimensional linear equation (see Subsection 2.10). However, the ellipse thus obtained is not necessarily of the minimum area (see Figure 3). Sometimes three of the five aforementioned vertices are sufficient to identify the MABE.

Figure 3. Let be a set of points of the Cartesian plane given by the following WKT geometry: MULTIPOINT ((-1, -1), (1, -1), (2, 0), (-2, 1), (-2, 0)). The MABE (solid red line) passing through the three points P 1 , P 3 , and P 4 has an area of about 8.46. In contrast, the MABE (dashed black line) passing through the four points P 1 , P 2 , P 3 , and P 4 has an area of about 8.49. Finally, the MABE (dashed blue line) passing through the five points has an area of about 11.22. [START_REF] Rublev | Minimum-Area ellipse containing a finite set of points[END_REF] show that if an ellipse passes through less than three vertices of a convex polygon and contains all the others, then it cannot be of the minimum area because it may be reduced by compression.

So, there are just three cases to examine: either the MABEs passing through 3 vertices, or the MABEs passing through 4 vertices, or the MABEs passing through 5 vertices.

The example in Figure 4 shows that there is a MABE passing through four of the six vertices and three MABEs passing through five of the six vertices. The best solution (minimum area) is generated by one of the 5-combinations.

All ellipsoids being convex, if a set of given points in the Cartesian plane is inside is defined as the greatest chord regardless of the selected pair of contour points.

Figure 4. Let be a set of points of the Cartesian plane given by the following WKT geometry: MULTIPOINT ((0 -1.5), (1 -1), (1 1), (0 1.1), (-1 1), (-1 -1)). The MABE (solid red line) passing through the five points P 1 , P 2 , P 3 , P 5 , and P 6 has an area of about 6.3. In contrast, the MABE (dashed black line) passing through the four points P 1 , P 3 , P 4 , and P 5 has an area of about 10.6. Finally, the MABEs (dashed magenta -resp. blue -line) passing through the five points P 1 , P 2 , P 3 , P 4 , and P 5 (resp. P 1 , P 3 , P 4 , P 5 , and P 6 ) has an area of about 13.5.

an ellipse, then its convex hull is also included in that ellipse. Moreover, if the set is finite, then the vertices of this convex hull are some of the given points.

The general principle of our method, inspired by the approach developed by Rublev and Petunin (1998), can therefore be summarized as follows. The Steiner ellipse of three non-collinear points [START_REF] Weisstein | Steiner Circumellipse. From MathWorld-A Wolfram Web Resource[END_REF] is the MABE containing these three points. This ellipse passes through these three points and has their isobarycentre as its centre. Therefore, if the Steiner ellipse of three vertices contains all the others, it is the MABE containing all the given points. If it is not so, we must examine the two other cases.

Pre-and post-processing

Isovists are, by nature, star-shaped and therefore 1-simply connected. As such, their contours frequently present numerous concavities or collinearities of points that can potentially overload or complicate their processing. To reduce the complexity and limit the number of points to be dealt with, the solution lies in convexity. The first step therefore consists, from the input contour, in calculating its convex hull. To this end, we use the native convexHull() method of the QgsGeometry class of the QGIS Python API. We can note that, in the specific case of isovists, the number of vertices of the convex envelope is generally much lower than that of the contour of the given isovist (see Table 1 in Section 3.1 for a more detailed illustration). Finally, when the isovist is bounded by an artificial horizon, prior pruning of contour vertices may be necessary.

In addition, the data we usually handle are geographical entities based on geographic coordinate reference systems. In these spatial reference systems, and more particularly in metric Cartesian systems such as the EPSG:2154 we use, the coordinates have high magnitudes and large numbers of significant digits.

To avoid possible problems related to computer arithmetic (especially when solving linear systems), the solution implemented is twofold. The first step aims to translate all the points of the convex envelope into a local coordinate system whose origin is an angle-weighted barycentre that will be detailed later. The second step -necessary when the shape diameter to be processed is at least hectometric -consists of a homothetic contraction aimed at uniformly reducing the shape7 .

The purpose of post-processing is, by reverse dilation and translation, to move and transform the centre of the ellipse, the two semi-axis lengths and the azimuth of its major axis in the input coordinate system and give them the required magnitude.

Main workflow

The purpose of our processing is to determine an ellipse from its centre, its semi-axis lengths, and the azimuth of its major axis. These parameters characterize the ellipse and allow to plot it easily. To assess these coefficients, we use the general principle set out in Subsection 2.1, which aims, from the set of vertices of the convex hull, to determine its unique MABE as the minimum-area ellipses passing through 3 or 4 and 5 vertices.

As shown in the workflow in Figure 5, regardless of the value of k and the kcombination of vertices, the objective is first to identify the centre of the ellipse in order to remove the linear coefficients from Equation A2 (see Appendix A8 ) and transform it into Equation A5. The next challenge consists in identifying the three quadratic coefficients a, b, and c of Equation A5 and the corresponding area. It should be noted that, for efficiency reasons, the minimum-area condition is checked before the enclosing condition. For the same reasons, semi-axis lengths and azimuth of the major axis are calculated only once, at the very end of the process, when the solution -the bounding ellipse of minimum area -has been identified. It is important to note that the initial centring on the entire convex hull is by no means universal. Moving into a coordinate system centred at the isobarycentre of the k vertices allows writing Equation A2. This naturally requires the reverse translation at the end of the entire process.

The workflow presented in Figure 5 summarizes the general principle we have implemented in Python in the context of QGIS. With the exception of the Pentagon case, the underlying principle is that any surface is modified proportionally9 by applying an affine isomorphism.

Shared processing: determination of the azimuth and semi-axes values

In this section, we assume that the ellipse is centred at the origin of the coordinate system. In such a coordinate system, the general formula of an ellipse is that of Equation A5. As mentioned in Appendix A, if a, b, and c are known parameters, we can assess the azimuth of the ellipse axes from Equation A9. Then, this azimuth having been determined, we can calculate the two semi-axis lengths using Equations A10 and A11.

Shared processing: test of vertices inclusion in a shifted ellipse

Any point with (x, y) as coordinates is enclosed in the ellipse defined by Equation A5 if and only if ax 2 + bxy + cy 2 ≤ 1. The "enclosing" test in the corresponding ellipse is based on the values of the coordinates of the centre and that of the three quadratic variables a, b, and c of Equation A5.

Shared processing: initial sorting and shifting procedures

The sorting procedure we have implemented aims to list vertices according to their increasing inner angle values. The heuristic we have implemented assumes that the minimum inner angle vertices are, in probability, the structuring vertices of the minimumarea ellipse.

In the same procedure, we are looking for a centroid. That is to say, a central point representative of all the points given. To prevent a possible shift of the centre due to an over-representation of points concentrated in very close locations, our barycentre calculation is weighted by the value of the angle at each vertex.

The triangle case

The case of a 3-combination of vertices or triangle case first requires identifying the centre of the corresponding Steiner (circum)ellipse. This centre is calculated as the isobarycentre of the three vertices of the triangle. Let us define P 1 = (x 1 , y 1 ), P 2 = (x 2 , y 2 ), and P 3 the points resulting from the translation which aims to move their isobarycentre at the origin of the Cartesian coordinate system.

To determine the quadratic coefficients, we have developed a purely geometric method (to bypass the resolution of a 3x3 linear system and therefore avoid possible computer arithmetic instabilities) that we will reuse further with the parallelograms. It consists first of all by applying an appropriate linear operator M, to transform the given triangle (P 1 , P 2 , P 3 ), into the equilateral triangle of vertices P 1 = (1, 0),

P 2 = -1 2 , √ 3
2 , and

P 3 = -1 2 , - √ 3 2
. This equilateral triangle has trivially the unit disc for MABE. Any point P of this disk is, by linear transformation of operator M, transformed into a point P = (x, y) and then, is written as follows:

P = M × P = m 11 m 12 m 21 m 22 × x y (1)
Providing the definition of the operator M:

M = m 11 m 12 m 21 m 22 = 1 -1 2 0 √ 3 2 × x 1 x 2 y 1 y 2 -1
If P satisfies the unit disc equation then P satisfies the following equation:

(m 11 x + m 12 y) 2 + (m 21 x + m 22 y) 2 = 1 (2)
By developing and term-to-term matching to Equation A5, we then obtain the following formulas for the three quadratic coefficients:

a = m 2 11 + m 2 21 b = 2 (m 11 m 12 + m 21 m 22 ) c = m 2 12 + m 2 22 (3)
Therefore, from the three quadratic parameters of Equation A5, we can analytically determine the semi-axis lengths, see Equations A10 and A11, and the azimuth of the ellipse axes, see Equation A9.

Non-parallelogram quadrangular case

Let P 1 = (x 1 , y 1 ), P 2 = (x 2 , y 2 ), P 3 = (x 3 , y 3 ), and P 4 = (x 4 , y 4 ) be the quadruplet of vertices resulting from the translation which aims to move their isobarycentre at the origin of the Cartesian coordinate system. Let M 13 and M 24 be the following two matrices from the coordinates of these vertices:

M 13 = x 1 x 3 y 1 y 3 M 24 = x 2 x 4 y 2 y 4
If det(M 13 ) = det(M 24 ) = 0 then the given quadrangle is a parallelogram. This particular case is handled in Section 2.9. Since we consider here that the quadrangle studied is not a parallelogram, at least one of these two determinants is non-zero. Let us consider that det(M 13 ) = 0 (otherwise the following reasoning is to be applied to M 24 ). Matrix M 13 is thus invertible.

In the manner of the triangle case, we can define an appropriate linear operator M aiming at transforming respectively the points P 1 and P 3 into P 1 = (1, 0) and P 3 = (0, 1). This operator is written as follows:

M = x 1 x 3 y 1 y 3 -1
The images P 2 and P 4 of P 2 and P 4 verify P 2 + P 4 = -(P 1 + P 2 ) = -(1, 1). These two vertices are therefore symmetrical relatively to (-0.5, -0.5). The convexity of the image quadrangle implies that they are each in one of the two grey sectors of Figure 6. Let us denote (u,v) the coordinates of the one in the northeast sector. The other is then -(1 + u, 1 + v). If the ellipse of Equation A2 passes through P 1 , P 2 , P 3 , and P 4 then the following system of equations is verified:

           a 2 + d 2 = 1 a 2 u 2 + b 2 uv + c 2 v 2 + d 2 u + e 2 v = 1 c 2 + e 2 = 1 a 2 (1 + u) 2 + b 2 (1 + u)(1 + v) + c 2 (1 + v) 2 -d 2 (1 + u) -e 2 (1 + v) = 1
Appendix B explains how to obtain the minimum-area ellipse passing through these four vertices.

Figure 6. In order for the four vertices P 1 , P 2 , P 3 , and P 4 to form a convex quadrangle described in the counter clockwise direction (P 1 and P 3 being fixed), P 2 must be located in the northeast grey sector resulting from the intersection of the three given half-planes.

The coefficients a 2 , b 2 , c 2 , d 2 , and e 2 having been determined, we still have to calculate the coordinates of the centre of the ellipse (see Formulas A3 and A4) as well as the quadratic coefficients of the centred ellipse (see Formulas A7).

As an additional remark, to prevent possible numerical instability problems when searching for the roots of Equation B6, it may be necessary to perform permutations of the image vertices described in counter clockwise order.

Specific quadrangular case: the parallelogram case

As already mentioned, the parallelogram case has symmetries and particularities that justify a specific treatment. Let P 1 = (x 1 , y 1 ), P 2 = (x 2 , y 2 ), P 3 , and P 4 be the vertices of the studied parallelogram, resulting from the translation which aims to move their isobarycentre at the origin of the Cartesian coordinate system. Let M be the matrix of the linear application that transforms the point P 1 into the point P 1 = (1, 0), the point P 2 into the point P 2 = (0, 1), and consequently the point P 3 into the point P 3 = (-1, 0), and the point P 4 into the point P 4 = (0, -1). One can notice that the MABE of the points P i resulting from this linear transformation is trivially the unit circle.

From this observation and the reasoning implemented in Equations 1 and 2, we quickly find the analytical formulations of the three quadratic coefficients (see Equations 3). It should be noted here that the operator M is simpler to write than in the case of the triangle (as it is a question of multiplying by the identity matrix):

M = m 11 m 12 m 21 m 22 = x 1 x 2 y 1 y 2 -1

The pentagon case

In the pentagon case, we have as input data five vertices of the convex hull. Provided that the origin of the Cartesian plane is not a point on the contour of the ellipse (which is guaranteed by the second translation specific to this k-combination of vertices), each coordinate of these five vertices satisfies the implicit Equation A2. From five vertices of the ellipse, we can then construct and solve 10 a linear system with five equations to determine the values of the five variables a 2 , b 2 , c 2 , d 2 , and e 2 . The values of these five variables having been determined, we then return to the purely quadratic Equation A5.

As mentioned in Appendix A, if a 2 , b 2 , c 2 , d 2 , and e 2 are known parameters, we can derive the coordinates of the centre of the ellipse and the quadratic parameters a, b, and c from Equations A3, A4, and A7. This is the processing we implemented.

Experimental validation

Case study: sites location

To illustrate the relevance of matching a given shape (and, in this case, an isovist as part of the open space) to the minimum-area bounding ellipse, we chose two French cities with various spatial configurations.

Regarding Le Havre city centre (see Figure 7a), we chose an urban area of 2 km 2 with 6,435 building footprints. This area has two distinct morphologies. It includes, in the south, a sub-area entirely rebuilt by Auguste Perret after the Second World War. In a totally planned layout, it is composed of a regular grid of relatively highrise buildings separated by large open spaces, right-angled streets, and quadrangular plazas of various sizes. On the other hand, the sub-area to the Northeast (Danton borough) has a more traditional morphology. In this area, the built density is higher and the street layout is more organic, with contrasting street widths, unconventional intersections, and unaligned facades. The second case is Paris 9 th borough (see Figure 7b). This area, of about 2.2 km 2 with 3,923 building footprints, has a complex and sometimes-ambiguous layout. Indeed, Haussmann's urban renewal in the late 19 th century has intensively influenced 10 Using the routine solve of the numpy.linalg package. its spatial configuration by superimposing broad rectilinear boulevards with narrow and more organic streets.

To analyse these two examples, we use the standard topographic data sets provided by the IGN BD Topo® database (March 2017 edition). This is a 3D metric and vectorial description (structured in objects) of the elements of the territory and its infrastructures, usable on scales ranging from 1:5000 to 1:50,000, compliant with the European INSPIRE directive.

Using a technique known as open space skeletonization (which was implemented, in particular, in [START_REF] Rodler | Local climate zone approach on local and micro scales: Dividing the urban open space[END_REF])), we positioned a set of about 300 virtual sensors on each of these two urban layouts. These sensors immersed in the city correspond to as many viewpoints. In each of them an isovist is generated, which presents a more or less great complexity of contour according to the urban environment. Table 1 presents a statistical study of the number of contour points of each of the 291 (resp. 300) isovists, and the number of vertices of their respective convex hulls, for Le Havre city centre (resp. Paris 9 th borough).

Table 1. Descriptive statistics of the number of contour points for the 291 (resp. 300) isovists and the number of vertices of their respective convex hulls, without the prior pruning of contour vertices mentioned in 2.2.

Site

Descriptive analysis

The boxplots in Figure 8 present a summary of the distribution of the six shape descriptors selected, after normalization of their measures, for various shape surrogates constructed from the 591 isovists under study. We have deliberately chosen not to represent the distribution of areas, perimeters and drift of the isovists -descriptors easy to calculate -in these diagrams, considering that the minimum boundary geometries are all derived from the convex hull, but also for readability reasons. Although, at first glance, the set of measures share broadly identical orders of magnitude, there are some disparities that we would like to highlight. The first point concerns area measurement. By observing the minimum, the 1st quartile, the median, the 3rd quartile and the maximum, we see that the minimum bounding circle is obviously less correlated to the convex hull than the rectangles and the ellipse. This lower correlation is confirmed (less obviously) with the perimeter and drift boxplots. As far as the drift is concerned, the minimum and median values show that the distance from the centre of the ellipse to the viewpoint is a slightly better estimator of the distance from the centroid of the convex hull to the viewpoint than are the rectangles.

In terms of stretching, the ratio of the two axis lengths of the MABE is a better estimator than the stretching of the rectangles. The median and 3rd quartile of the major axis of the ellipse correlates better with the diameter of the bounding circle than do the longest sides of the bounding rectangles. Finally, in terms of orientation, it seems that the ellipse is even better (in the sense that it correlates better with the direction indicated by the diameter of the isovist) than the rectangles.

An urban fabric dependent statistical relationship

This better correlation (more or less obvious depending on the considered descriptor) of the MABE is confirmed by the scatterplots and the associated correlation coefficients in Figures 9 and10. In all cases, the MABE is the surrogate shape that best correlates with the so-called reference shape (i.e., depending on the selected measure, the convex hull for the area, the perimeter and the drift; the caliper for the stretching; the MBC for the length; the diameter for the orientation). This trend is significant for orientation (21.3% in Le Havre city centre and 18.6% in Paris 9 th ) or drift ( 10.7%), of the order of one percent for length, and marginal for perimeter, stretching or area. Finally, it should be noted that the correlation gain is systematically stronger in Le Havre than in Paris for the three measures in Figure 10, whereas this trend is reversed in the case of Paris, where the ellipse correlates slightly better with the convex hull of the isovist.

Zoom on the orientation and drift descriptors

As already noted, orientation ( 20%) and drift ( 10%) are the two descriptors for which, in terms of minimum bounding geometry, the ellipse stands out clearly from the rectangles for correlation to the reference shape. The purpose of Figures 11 and12 is to highlight the fact that this better correlation is not only marked but also fairly localized. Figure 11 shows a distribution of the deviations between the orientation of the diameter and the orientation of the MABR for Le Havre city centre (a) and Paris 9 th (c) on the one hand and a distribution of the deviations between the orientation of the diameter and the orientation of the ellipse for Le Havre city centre (b) and Paris 9 th (d) on the other hand. It is clear from these maps that the strongest orientation deviations are uniformly distributed over the two studied areas. In Paris 9 th , deviations of more than 30°(i.e. 20% of the cases) are almost all located at street intersections or on squares or plazas. This clear trend is confirmed in Le Havre city centre (more than a third of the cases have an orientation deviation of more than 30°) where the space is much more open in the southern part of the city.

Figure 12 shows a distribution of the deviations between the drift of the convex hull and the drift of the MABR for Le Havre city centre (a) and Paris 9 th (c) on the one hand and a distribution of the deviations between the drift of the convex hull and the We note that the distance from the centroid of the rectangle to the corresponding viewpoint is significant at T-or Y-intersections in the case of Paris.

Focus on a few viewpoints

Discussion

Towards a better understanding of the city layout

Ellipticity is a shape characteristic that is intermediate between compactness (which is a measure of the complexity of the interface between the study area and its surroundings) and rectangularity (which is correlated to the notion of flattening). Therefore, the properties of the ellipse and its major axis make it possible to determine in particular an orientation and drift for the convex hull of the open space surrounding any point of view of the urban space. These orientation and drift, which draw major trends in the spatial configuration of open urban space, can also be studied in relation to considerations of exposure to direct solar beams or to wind (speed and kinetic energy).

A highly simplifying indicator

It must be noted that, by focusing on the main orientation of the MABE, we are achieving a drastic dimension reduction (from a polygonal shape to a single floating-point value). This indicator is all the more simplistic because it is based, for its determination, not on the shape itself but on its convex envelope. Two distinct isovists with the same convex hull will thus produce identical MABE. This indicator is therefore not a tool for measuring the concavities or convexity defect of the studied isovist. However, one can notice that the MABR and the MBC face the same limits.

About complexity

Compared to the complexity of the algorithm for calculating the MABR (which is known to be a linear-time algorithm), the complexity of our algorithm for determining the MABE is of a complexity at worst in O(n 5 ), where n is the number of vertices of the convex hull. Indeed, if none of the 3-combinations of vertices allows to identify an ellipse enclosing the given shape (complexity at worst in O(n 3 )), then it is necessary to review all the ellipses resulting from the 4-and 5-combinations of vertices.

As shown in Table 2, in 97 cases (33.3%) in Le Havre city centre, we have a complexity at worst in O(n 3 ). The remaining 66.6% have a complexity in O(n 5 ). Similarly, in 96 cases (32%) in Paris 9 th , we have a complexity at worst in O(n 3 ). The remaining 68% have a complexity in O(n 5 ). The number of vertices of the convex hull is decisive for the determination of the MABE and its complexity. This number and its dispersion (see Table 1) are clearly smaller in Paris 9 th than in Le Havre city centre. A more precise study of the Le Havre case shows that, in terms of the MABE calculation, the complexity of the Danton borough (Northeast) is close to that of Paris 9 th , while that of the Auguste Perret borough (South) is much higher due to artificial horizon generated by long perspectives.

Conclusion

This article is a contribution to the qualification of urban open spaces. It provides an operating algorithm for determining a minimum-area bounding ellipse for any given polygonal shape. More precisely, using the theorem of existence and uniqueness of Löwner-John's ellipsoid, we develop a method based on three, four, and five relevant combinations of vertices of the convex hull of the given isovist field. The application to two real urban configurations shows that the ellipse is a minimum bounding geometry that performs better than the circle or rectangle, at least for the half-dozen studied descriptors. This improvement, substantial as far as orientation or drift is concerned, has however a drawback in terms of computing time since the implementation we propose has an algorithmic complexity at worst in O(n 5 ).

In the context of a systematic morphological analysis of cities, the question is raised as to whether this combined approach would make it possible to massively produce city typologies through open spaces' classification (in particular with regard to the issue of the universality of the Zipf's law raised by [START_REF] Jiang | Zipf's law for all the natural cities around the world[END_REF]). We may also wonder whether the method of automatic vehicle extraction from point cloud based on three descriptors, namely surface area, rectangularity and elongation, developed by [START_REF] Zhang | Automatic Vehicle Extraction from Airborne LiDAR Data Using an Object-Based Point Cloud Analysis Method[END_REF], could also be extended to the detection of other features by taking advantage of the contributions of our work on ellipticity.

Data and codes availability statement

The data and codes that support the findings of this study are available with a DOI at http://doi.org/10.6084/m9.figshare.c.4964291. In particular, we provide a Jupyter notebook that explains step by step how to use the Python plugin (released under the GPL v3 license) that implements the method presented in this article.

Figure 1 .

 1 Figure 1. Comparisons of three dominant directions (main axis of the minimum-area bounding ellipse -black dashed segment, the minimum-area bounding rectangle -blue dashed polygon, and the isovist diameter -red dashed segment) for a given spatial configuration in Paris 9 th borough, France. The isovist corresponding to the viewpoint (black point labelled 384) is represented by a yellow polygon.

Figure 5 .

 5 Figure5. Workflow (for instance, the abbreviation "Sec. 2.10" refers to the Section 2.10 of this article).

Figure 7 .

 7 Figure 7. a) Map of Le Havre city centre; b) Map of Paris 9 th borough.

Figure 8 .

 8 Figure 8. Boxplots of the six normalized descriptors for various shape surrogates derived from the set of 591 isovists examined in Le Havre city centre and Paris 9 th .

Figure 9 .

 9 Figure 9. Correlation matrices for area, perimeter and drift descriptors in Le Havre city centre (left) and in Paris 9 th (right).

Figure 10 .

 10 Figure 10. Correlation matrices for stretching, length and orientation descriptors in Le Havre city centre (left) and in Paris 9 th (right).

Figure 11 .

 11 Figure 11. Comparisons of the deviation of the orientation of the rectangle to the orientation of the diameter for a) Le Havre city centre and c) Paris 9 th with the deviation of the orientation of the ellipse to the orientation of the diameter for b) Le Havre city centre and d) Paris 9 th .

Figure 12 .

 12 Figure 12. Comparisons of the deviation of the drift of the rectangle to the drift of the convex envelope for a) Le Havre city centre and c) Paris 9 th with the deviation of the drift of the ellipse to the drift of the convex envelope for b) Le Havre city centre and d) Paris 9 th .

Figure 13

 13 Figure 13 maps four urban configurations in which deviations from reference, both orientation and drift, are particularly strong.

Figure 13 .

 13 Figure 13. Zoom on some spatial configurations with strong deviations of direction in a) Le Havre city centre or b) Paris 9 th and strong deviations of drift in c) Le Havre city centre or d) Paris 9 th .

  Feature Min. Median Mean Max. Std. dev.

	Le Havre	Isovist	10	146	190	812	165
		Conv. hull	4	12	20	95	18
	Paris 9 th	Isovist	11	83	94	256	48
		Conv. hull	4	9	10	24	4

Table 2 .

 2 Number of MABE generated by 3, 4, or 5 vertices for the isovist fields in Le Havre city centre (resp. Paris 9 th ), France.

	Site	Tri. Quadri. Penta.
	Le Havre	97	142	52
	Paris 9 th	96	172	32

We adoptŽunić et al. (2017)'s definition: "Notice that by the shape ellipticity we mean the similarity of a given shape to the planar region bounded by an ellipse".

Let's call chord length the distance in a straight line between two points of the shape contour. The diameter

The Feret diameter, also known as caliper diameter, is defined as the distance between two parallel tangents of the shape contour. It depends on the chosen direction[START_REF] Toussaint | Solving Geometric Problems with the Rotating Calipers[END_REF].

The shapes under study in this paper all derive from the isovist and as such are generated from a viewpoint. By extending the concept defined by[START_REF] Conroy | Spatial Navigation In Immersive Virtual Environments[END_REF]), we call "drift" the Euclidean distance between the generating viewpoint and the centroid of the corresponding shape.

For any non-circular polygon shape, orientation is defined as the measure of the azimuth of its diameter in the range of 0 to 180°.

Indeed, the QgsEllipse class as defined in QGIS 3.0 core library is a container for elliptical geometries defined by a centre point, a semi-major axis, a semi-minor axis and an azimuth. This class does not in any way provide a solution for determining a minimum-area bounding ellipse.

The underlying idea is to solve linear systems as well-conditioned as possible.

In Appendix A, we have grouped together the few general equations related to ellipses that we use in this paper. All equations with labels beginning with the letter A refer to this appendix.

The proportionality factor being equal to the absolute value of its determinant.

Appendix A. Some general equations related to ellipses

In analytical geometry, if (x, y) are the coordinates of a point M in the Cartesian plane, the ellipse described using M is defined by the general equation of a conic as:

If the origin of the Cartesian plane is not a point on the contour of the ellipse, Equation A1 can be simplified as follows:

By a vector translation (x 0 , y 0 ), with:

the Equation A2 is simplified as follows (ellipse centered at the origin):

Its area is then:

Minimizing the area of the ellipse is equivalent to maximizing the square of the denominator of the Equation A6. Let us define:

the three parameters of Equation A5 become:

and the square of the area as a function of the non-centered coefficients a 2 , b 2 , c 2 , d 2 , and e 2 of the Equation A2 becomes:

Moreover, by a proper rotation of θ angle, where:

Equation A5 is simplified as follows: αx 2 + βy 2 = 1. The axes of the ellipse then overlap the axes of the Cartesian coordinate system. We can deduce the semi-axis lengths using the following formulas:

Appendix B. Some specific equations related to the general quadrangular case

The ellipses passing through the four points (1, 0), (u, v), (0, 1), and (-1 -u, -1 -v) form a one-parameter bundle. We choose e 2 as parameter which we denote t. The equations at the end of Subsection 2.8 allow us to write the following system of polynomial linear equations:

Let us denote respectively:

The previous system provides the coefficients a 2 , b 2 , c 2 , d 2 , and e 2 as the following 1st degree polynomials:

An area being positive, minimizing an area is equivalent to minimizing its square. The function to be minimized is therefore the following rational fraction:

Let ∆ and Γ be the two following polynomials:

They are respectively of degrees at most 2 and 5. And the rational fraction to minimize becomes:

Its derivative is the following rational fraction:

To minimize the area of the ellipse, it is necessary to find the roots of the derivative of the rational fraction of Formula B6, i.e. find the roots of the numerator (2Γ ∆ -5Γ∆ ) which is a polynomial of degree at most 5 (the nullity of its coefficient of degree 6 is trivial). Out of the five possible corresponding roots, the solution t 0 is the real one for which ∆ is positive on the one hand and which minimizes the Γ 2 /∆ 5 ratio on the hand. The value t 0 having been identified, we obtain by Equations B1-B5 the values of the coefficients a 2 , b 2 , c 2 , d 2 , and e 2 .