Procedural modeling of plant ecosystems maximizing vegetation cover

Vegetation plays a major role in the realistic display of outdoor scenes. However, manual plant placement can be tedious. For this reason this paper presents a new proposal in the field of procedural modeling of natural scenes. This method creates plant ecosystems that maximizes the covered space by optimizing an objective function subject to a series of constraints defined by a system of inequalities. This system includes the constraints of the environment taking into account characteristics of the terrain and the plant species involved. Once the inequality system has been defined, a solution will be obtained that tries to maximize the radius of the projected area of the trees and therefore the extension of the vegetation cover on the ground. The technique eliminates the trees that do not achieve a minimum growth radius, simulating the typical competitive process of nature. Results show the good performance and the high visual quality of the ecosystems obtained by the proposed technique. The use of this kind of optimization techniques could be used to solve other procedural modeling problems in other fields of application.


Introduction
The representation of outdoor scenes is currently a very popular field of research due to the huge number of applications that require them, such as ecosystem simulation [20,26], video games [27], or geographic information systems [6,39]. Procedural modeling of individual plants and trees has been widely studied in the literature [38]. Nevertheless, distributing them automatically so as to fill the space available is still an open problem. Smelik et al. [41] reviewed some studies that use procedural techniques for the creation of natural environments and virtual worlds. Specifically, they estimate the virtual position of the plants from a procedural calculation of their distribution. They conclude, however, that these techniques achieve unrealistic results, since they lack the ability to translate ecological input data into procedural distributions.
In order to increase the realism of the outdoor scenes, it is necessary to take into account some botanical information in the procedural placement, such as the biotic and abiotic characteristics of the species, as well as the volume that they take up [11,12,36]. Smelik et al. [41] stated that an ecological plant distribution with sufficient density and variety, where each vegetation species is placed in a predetermined position, is a difficult challenge to achieve [2,3,8,14].
According to Lane and Prusinkiewicz [29], ecological modeling techniques can be divided into two groups: local-to-global and global-to-local. On the one hand, the procedural methods classified as local-to-global treat plants as individual entities to which rules are applied in order to simulate factors such as reproduction, growth, competition, interaction, variation or adaptation to the environment [15]. This biological model of plants predicts patterns of vegetation distribution in an ecosystem where plants grow, adapt to the environment, reproduce and compete for space. These methods require a high consumption of resources and the processing time for each simulation cycle is also high. On the other hand, the procedural methods classified as global-to-local calculate the positions of the plants directly from a globally defined environment. Thus, they include aspects related to the terrain such as height and slope, climate characteristics or vegetation cover density.
In order to avoid the problems of local-to-global techniques [29], this work presents a method for the distribution of trees and plants in natural environments through the use of linear programming techniques. In this case, the objective function to be maximized is the projected area of the set of plant species. This means trying to achieve the largest possible vegetation cover with certain constraints. These constraints are defined by a system of inequalities that includes the biological characteristics of the plant species, such as the maximum growth radius or the relationships between the growth radii of the different species.
In this sense, after making an initial irregular distribution of the seed positions where the different plant species can be placed, the user is allowed to indicate some initial growth regions of each of the plant species. Finally, a distribution of the species is made taking into account for each placement position both the characteristics of the plants, as well as the characteristics of inclination and height of the terrain. Once the final distribution of the plant species has been determined, the system of inequalities that includes the radii of each of the plants as variables is created.
The novelty of this methodology is based on providing a global solution to the distribution of plants and trees, avoiding the drawbacks of iterative techniques. Furthermore, the technique can easily be adapted to the definition of other types of objective functions and any modification or inclusion of new constraints. Finally, it should be noted that the implementation of the system is simple and robust, which makes it ideal for solving this type of problem.
The paper is organized as follows. Section 2 considers the state of the art related to the existing techniques in the procedural generation of natural environments. Section 3 briefly describes the method operation scheme. Section 4 analyzes the distribution of seed positions, and Section 5 describes the labeling process for each seed positioned on the terrain. Section 6 describes the space-filling process performed by the inequality systems, which makes it possible to obtain a high degree of realism. Section 7 analyzes the results obtained by testing the method and, lastly, Section 8 presents some conclusions and observations related to future work.

Previous work
Designing a realistic outdoor scenario with an ecological distribution of plants is not an easy task to achieve. The biotic and abiotic characteristics of the vegetation species [10,44] are essential in the process of obtaining the exact positions of the plants in a terrain where a high density and variety of species are required [32]. Onrust et al. [34] classified the techniques for plant species placement in two main groups: those based on ecological modeling techniques and those based on the procedural generation of ecosystems.
Studies classified in the first group distribute the plants based on real data related to the terrain or the vegetation species to be represented. These techniques can, in turn, be classified in two sub-groups depending on the source of the data used as input. On the one hand, dynamic ecological modeling techniques simulate ecological processes taking into account data that can be extracted from a raster map, such as height, biomass or vegetation cover [6,20,40,42]. On the other hand, neutral ecological modeling techniques [30] generate a distribution of plants on a terrain following a grid structure that is based on some vegetation species information.
According to Onrust et al. [34], the second group is related to the procedural generation techniques that perform an automatic distribution of plants and trees. Studies in this group estimate their position using different algorithms. Nevertheless, these techniques usually fail to create a distribution of plants with sufficient variety and realism because no ecological input data are considered. Over the years, several studies that address this problem have appeared.
Deussen et al. [15] proposed the framework of the Open L-Systems model [31] to simulate the individual-based ecosystem models. The growth of an individual ecosystem is represented using the dynamic behavior model of Firbank et al. [18]. This model starts from a set of randomly distributed circles with initially random radii that represent plants growing in the ground. In a second phase, control by the user is incorporated through a gray-scale image that represents the density of the plants. Individual plant positions that fit these densities are created using the Floyd-Steinberg error diffusion algorithm [19]. Finally, the points produced by this algorithm are altered slightly to make the distribution more random. Nevertheless, the resulting distribution of the plants is too uniform. Moreover, the cost of the process in terms of time is high, according to Lane et al. [29].
These authors [29] proposed the symbiosis of multiple vegetation species to generate plant distributions. They combine a probability value with a dart-throwing algorithm that allows plants to be placed in their most probable areas. In addition, each plant can produce an inhibition or expansion effect on the remaining plants, by updating the probability field around it. However, the grammatical elaboration method turns out to be quite complex and the computation times are relatively high.
A method called FON (Field Of Neighborhood) was proposed by Bauer et al. [7]. It describes a circular zone of influence around the plant, the radius of which determines the distance from which it can interact with its neighbors. This radius depends on the soil characteristics, size, and nutritional needs of the plant and the space needed for that nutrition. Based on FON, Alsweis et al. [4] proposed a method that simulates plant communities of two different species in evolution. Later, in [5], they calculated the distribution of the plants by combining Poisson Disk Distributions (PDDs) [28] with Wang-tiles [17], although neither the distribution by species nor the interactions between them are considered.
Weier et al. [43] extended the previous method and generated large areas of vegetation using Wang-tiles with reusable distributions that allow a terrain to be covered. Each Wangtile is generated using PDDs with a variable radius. The tiles are organized in a multilevel structure and treated with a nested KD-tree, which allows them to be reused during visualization. In addition, they propose a solution to represent plants whose area of influence is on the edges of the tiles. The main problem with this proposal is that the size of the disks that line the ground tend to be too homogeneous.
Onrust el al. [33] presented a procedural plant placement technique that combines the Poisson Disk Distribution with Wang-tiles, thereby establishing an ecological modeling technique. Results were only tested with isolated plants that had problems with the position of the impostors with respect to the camera. The work proposed by Onrust et al. [34] is based on vegetation location information resulting from the combination of several parameters such as the Normalized Difference Vegetation Index (NDVI), height-map, biomass or coverage maps. Gasch et al. [23] presented a method of automatic vegetation placement that takes as its input the height-map of a terrain where a route is represented [24]. The method extracts the values necessary to distribute the vegetation species from the height-map and combines them with a noise distribution based on Perlin noise [35]. By doing so, it becomes possible to obtain placement patterns with a greater probability of occurrence in favorable points on the map and near the route. However, this method does not take into account the development and growth of plants.
Li et al. [30] proposed a procedural model based on graphs that generates realistic vegetation distributions with different compositions and configurations adjusted to the metric of a specific landscape created from patches. The metrics used include the number of tree species, their proportions in the landscape, the number of patches that have to be populated by each species and their average shape. The seeds are placed randomly and their species are assigned in order to obtain the final distributions of the trees. The method establishes the distribution of each species with a time efficiency that satisfies the requirements to build large-scale virtual forests.
A more generic method to describe growth and simulation of ecosystems is presented in Growth-Grammars [16], which provides the user with access to a tool for developing growth grammars through an interactive modeling platform, GroIMP. This grammar is executed through programs written in XL (a Java-based growth grammar language) and visualized by XL4C4D (an XL-plugin for CINEMA 4D). Although this method allows populations of plants to be generated in small areas, it does not reuse certain components to cover more extensive areas.
Benes et al. [9] described a method to distribute vegetation in urban environments procedurally. They addressed the problem of spatial distribution of vegetation by combining a user-guided simulation system with an interactive procedural system that integrates plants into the urban environment. Starting from the 3D geometry of an urban layout with initial conditions and parameters, procedural rules are inferred. A manageability level is calculated for each area of the urban space, which defines the degree of influence between the simulation of the wild ecosystem, where plants compete freely for resources and seeds, and the ecosystem is managed by the user and the plants grow only under well-defined conditions. The method presented takes some of the weaknesses found in the literature into account and proposes a solution that includes most of the detected shortcomings. How it works is analyzed in the following sections.

Plant covering method
The main objective of this method is to maximize the size of the vegetation cover in nature scenes. To obtain a realistic representation, the method considers three main processes: the distribution of positions where seeds can be placed, the assignment or labeling of the seeds to those positions and the resolution of the system of inequalities that will allow to maximize the area of vegetation cover (see Fig. 1).
The initial distribution of the seed positions can be done with any strategy. In this case, a regular distribution has been made whose positions are later slightly modified to obtain a more natural irregular appearance.
Once the positions where the different plant species can grow have been determined, a subset of the positions is chosen and an initial labeling is carried out taking into account some characteristics of the terrain, such as slope and height. In this way, an initial distribution of the plants is already obtained. Subsequently, the final labeling of all the seed placement positions is carried out taking into account the closest neighbors.
The method ends with the creation of a system of inequalities that poses the constraints regarding maximum radius and distance relationships between species. The resolution of the system will determine the radius of each of the trees or plants that maximizes the area of vegetation cover. After obtaining these radii, all the species that have not been able to reach their minimum radius are eliminated and the system of inequalities is redone a second time to obtain in this case a more precise solution.

Seed dissemination process
The method presented here starts by determining the areas without vegetation on the terrain. A criterion to choose these areas can be, for example, zones situated higher than a Fig. 1 Scheme of the proposed model certain height, or areas where some paths, lakes or rivers are located. However, other ways to achieve these non-vegetation areas can be considered in the case of computer-generated terrains.
In this work, the terrain was obtained using the Perlin noise method [35]. The process starts by assigning a value that represents the possibility of being populated by plants (between 0 and 1) to each point in the terrain that has been generated. A threshold is then established to determine the areas that can simulate grasslands, sandy areas, rocky areas, etc.
The initial positioning follows a regular grid distribution where the seeds are separated by a distance of d spacing = 5 meters. For this purpose, a relation of 1 meter = 1 pixel was considered in its representation within the placement map. Then, in order to break up the regular location of the seed positions, the method performs a deviation of the seed positions s p , by modifying their (x, y) coordinates with a random value using a uniform distribution of random number generator. This random value will always be less than half the spacing used in the generation of the grid, in order to prevent two elements from overlapping.
The final result of the redistribution can be seen in Fig. 2a, where areas with no vegetation are represented in black, and a closer view can be observed in Fig. 2b. Both figures show the deviation suffered by the seed positions, thereby eliminating the unrealistic regular appearance and covering all the terrain.

Seed labeling process
The next process is to determine the species to which each seed belongs. Initially, only a sub-set of the distributed seeds are evaluated by taking into account several abiotic features of the terrain and a simulation of the grouping of plants of the same species that happens in nature. Finally, the influence of the neighboring vegetable species determines the label of the rest of the seeds that remain unlabeled on the terrain.
Only the height and the slope of the terrain have been considered in this work. These abiotic characteristics have been included in a linear equation that establishes the appropriate vegetable species depending on the features taken into account. However, more features could easily be added to the method.

Influence of the abiotic characteristics
The slope and the altitude of the seeds are data available from the height-map. Both characteristics are taken into consideration to analyze the possibility of each vegetable species growing at that location. In their work, Galin et al. [21] stated that the ability of the species to live is limited by multiple conditions, so that although a vegetable species has a perfect soil, it cannot grow if it does not have sufficient temperature or light. In this work, if the position of a seed is not in the range of the altitude or slope values where the vegetable species under consideration can grow, it will not be able to develop and it will die.
Let V be the set of vegetable species to be distributed and v a particular species in that set. Since s p is the position of the seed, the altitude of the seed can be obtained directly in meters from the height-map height (s p ) (Fig. 3a). The slope where a seed is situated, s p , is computed using the following equation: where N is the vector normal to the horizontal plane and N(s p ) is the normal vector at the point where the seed is situated. The slope(s p ) is measured in degrees. Figure 3b shows the slope characteristics of the terrain under consideration.
Analyzing the plant distribution depending on the height, the distribution of plants per species w(s p , v) is given by the following equation: where for each species v ∈ V , a v is the value at the highest point of the peak of the curve of the Gaussian distribution, b v is the average height value, and c v is an approximation of the standard deviation of the heights for that species. Finally, the distribution of plants per species w(s p , v) has to be normalized to obtain how the height w height (s p , v) conditions a vegetable species v to be assigned to a particular seed s p . It is obtained by means of (3): where min(w(v)) is the minimum value of w(s p , v) for species v considering all the seeds on the terrain and max(w(v)) is the maximum. Assuming the same reasoning for the slope, the value w slope (s p , v) would be obtained in a similar way to (3), but in this case considering the function slope(s p ) that measures the slope of the terrain where the seed is positioned, instead of height (s p ) in (2).

Influence of the grouping of the vegetation species
In a real outdoor environment, it is common to find plant groups of the same vegetable species. In order to simulate these groupings, the method allows for the addition of areas where a specific vegetable species has more likelihood of growing. The number of grouping areas where a vegetable species is more profuse is an input parameter of the method. They are represented by certain points on the terrain, called diffusion points, considered as the center of a circular area where the grouping is allocated. These diffusion points can be positioned randomly in the work presented or can be specified by the user. Figure 4a shows an example where some diffusion points for three different species have been distributed. They have been differentiated by three colors: blue, red, and yellow.
For each seed on the map, the method searches to determine whether it is situated in the area of influence of one or more diffusion points. This value will reach its maximum when it is right in one diffusion point and will be reduced until reaching the limit of its area of influence, with a value of 0 just when leaving it. Equation (4) illustrates the calculation performed to measure the influence of each diffusion point diff v belonging to the species v on a seed s p as follows: where d(s p , diff v ) is the distance of the seed under consideration s p from the diffusion point diff v and r(diff v ) is the radius of the area of influence associated to it. Finally, the influence of each diffusion point set by the user for a vegetation species v has to be evaluated on a seed s p following the expression:

Labeling process
The abiotic features and the grouping simulation have to be considered together to perform the labeling process w sp (s p , v). Equation (6) shows how this weight is obtained by performing a linear combination between the influence of the terrain, w height and w slope , and the influence of the diffusion points w diff . The user can modify the degree of importance of each of them by managing several coefficients (c height , c slope , c diff ). By default, the same value is maintained for the terrain features as for the area of influence, i.e., c height + c slope = c diff . If the seed that is being evaluated is outside the area of influence of that species, w diff = 0 and the weight is obtained considering only the height and the slope.
The vegetation species v that will be assigned to that seed s p will be the one that gives w sp (s p , v) its maximum value.
being V the set of the different vegetation species that are going to be present in the outdoor environment.

Computing neighboring information
Once the initial subset has been labeled, the rest of them have to be processed. In nature, it has been demonstrated that plants can condition vegetable species that grow near them, because of certain facts such as the dissemination of seeds from them. Some studies in the literature have simulated these clusters by means of different techniques. For instance, [29] used kernel deformation to create an algorithm to label species by neighborhood. In the method presented here, it has been considered that each seed conditions the label of its neighbors to create this influence. Consequently, this has been taken into account in the labeling process: the method uses a k-nearest neighbor classifier [13], in our case with k = 3 neighbors.
Once the sub-set has been labeled, the process assesses the rest of the seeds by analyzing the label of the three nearest seeds in the initial set. The label is assigned according to the vegetation species with the majority voting among these three neighboring seeds. This makes it possible to simulate the neighborhood that occurs in a real environment, allowing some plants to condition the species of the nearby seeds. Figure 5a shows the labels applied to the seeds in the initial set, a different color being assigned to each vegetation species. If the position of a seed is not appropriate for any of the vegetation species to be distributed, it is marked in black and an empty label is assigned to it. Figure 5b illustrates the labels of all the seeds on the terrain after performing the labeling process that takes the neighboring vegetable species into account. This figure shows some areas where a vegetation species has a strong influence, and some others where different vegetation species appear in the middle. This fact simulates the randomness of nature.

Inequality systems to maximize plant representations
The main contribution of this work is the design of inequality systems adapted to each seed in order to compute the size of the plant representation. Once a vegetable species has been assigned to every seed, the size of each plant representation has to be obtained in order to cover most of the terrain. This goal was accomplished by designing a set of inequality systems that consider different (biotic and abiotic) features to estimate the maximum size of each plant without overlapping between neighboring plants. This size is represented by a circle around the seed, similar to the one used in the FON model [4,5], which determines the maximum area that can be reached by the plant representation. Its radius can be estimated tacking into account that the projected plant canopy on the ground is maximized. This problem can be treated using an optimization technique that establishes relationships between distances with constraints as presented for foams in [25]. In this paper, inequality systems consider some biological data about the vegetable species assigned to the seed in addition to some general constraints: -The circle of each plant cannot intersect with any other plant. This is not true in nature, where a certain overlapping between the tree canopies does exist, but this overlap will be performed during the final visualization process. -The radius of the circle of each plant cannot exceed the maximum radius that the vegetation species is biologically capable of reaching. -If the circle associated with the plant does not reach the minimum radius for the species, that plant cannot survive.
Given a distribution of points that represent the seeds, a distance matrix is constructed. Based on this information, we can estimate the area occupied by a set of plants on a terrain or domain Ω. Let C r (s p ) be any circle centered on point s p with radius r and whose area A is bounded by: Let P = {s p1 , · · · , s pn } be a set of n isolated seed positions representing the coordinates where the plants have been planted. Considering Ω as the possible locations on the terrain, i.e., Ω ⊂ R 2 , this set P ⊂ Ω is associated with the subset of R 2 , given by (9).
D Ω (P ) = (r 1 , r 2 , · · · , r n ) ⊂ R n | C r i (s pi ) ⊂ Ω for i = 1, · · · , n and C r i (s pi ) ∩ C r j (s pj ) = ∅ for all i = j (9) where D Ω (P ) corresponds to the solution space to which all the radii found by themethod belong. The solution D Ω (P ) found from the above equation guarantees that the circular areas of the plants do not intersect each other. These regions called circular areas without intersection (cawi) in set P are determined according to The following example clarifies the procedure.  Figure 6 illustrates these data by means of a graph that relates all these distances between a plant of each vegetation species and the vector M, which represents the maximum radius that each plant can reach. From the graph, the following system of inequalities can be extracted: This system of inequalities can be identified by S ·r ≤ b, where S is a matrix of zeros and ones that relates the different plants to each other and to the vector M of maximum radii; r is the radii vector to be solved, in this case (r A , r B , r C ); and b contains the independent terms of the inequalities.
In this example, the solution that maximizes the cawi value for the three plants is given by (12).
On observing the system of inequalities, it can be seen that the constraints do not include the minimum radius that each vegetation species requires in order to live. The reason for this is to favor the growth of large trees rather than small trees, since these provide the cawi value with a greater numerical result. This decision allows plants with a radius lower than their minimum radius to be eliminated in the same way as was performed with a self-thinning model in [29]. Finally, the whole process is repeated a second time so that the space left by the plants that are removed can be covered by the rest of the plants. Figure 7 shows the result of applying this growth method in the terrain used as an example by solving the system of inequalities. As can be observed, the previous points with the same size shown in Fig. 5 now appear re-sized with different magnitudes.
The proposed method is flexible when it comes to generating different plant distributions depending on the user requirements. If required, a more regular distribution could be generated on the terrain, for example by homogenizing the size of the plant population by not allowing the existence of small plants, or simulating the possibility that if a plant reaches its maximum size, it can die of old age and then its space can be occupied by the plants around it.

Results
Several experiments were performed using a computer with an Intel Xeon E5 v2 3.70GHz processor. The trees and textures used in the examples belong to the Landscape Automaterial plugin [22]. Unity 3D version 2018.2.7 was used to view the resulting outdoor environments with a realistic visualization. Fig. 7 (a) The cawi obtained is shown for each seed on the terrain. (b) The detail on the right shows the different sizes of the circles Three different vegetation species were considered in the experiments: Pinus Sylvestris, Abies Alba, and Sorbus Aucuparia (Fig. 8). They have similar characteristics in terms of the type of ecosystem they need to live in, but conversely they differ in size and in the characteristics of the terrain (height and slope) on which they can grow. The data requirements for each of them were acquired from real data [1,36,37] and are summarized in Table 1. This table provides information about the minimum and the maximum radii that a plant of that species can reach, the minimum and the maximum height where it can live, and the maximum slope where that species can appear.
In order to evaluate the influence of the terrain features on the labeling of a seed, the values a v , b v , c v are required for each species so as to be able to compute (2). These values were obtained from Growth-Grammars [16]. They are shown in Table 2, classified according to the height and the slope of the terrain.
Regarding the diffusion points, the ones that were considered are those shown in Fig. 4. The labeling process was performed following a 3 nearest neighbor classifier, as was previously reported. Finally, the system of inequalities was solved by applying the Simplex optimization technique, since this technique makes it possible to obtain the best solution without constraint in the number of variables, thus allowing these variables to be expanded should it be necessary to do so in the future.

Validation of the method
The method has been tested using a height-map with a resolution of 512 × 512 meters. As has been said in previous sections, the relation established between pixels and meters has been 1 meter = 1 pixel. Regarding the height of the terrain, a height level difference of 500 meters has been represented between the lowest and highest areas, with a range of values that vary between 800 and 1,300 meters.
Some tests were performed to validate the procedurally created outdoor environments, specific modifications being applied in each execution. Initially, the areas with no vegetation are established by limiting the height values: areas higher than 3,000 meters are considered dry zones. The same map was used in all the experiments to make it easier to notice the visual differences when some parameter values change.
Three different vegetation species were considered and differentiated by colors: Abies Alba is represented in red, Pinus Sylvestris is shown in blue, and Sorbus Aucuparia is shown in yellow. Different parameter variations were tested in order to see how the distribution  The minimum and maximum radii of the bounding box that the plant can reach and, regarding the terrain, the minimum and maximum height and the maximum slope where each plant can grow. All these data are measured in meters and the slope in degrees changes according to the user's preferences. All the data involved in the tests are shown in Table 3 and the results obtained can be seen in Fig. 9.
The first test has been called Basic map (Fig. 9a). Initially, the values that weight the parameters involved in the labeling process of the training test were set to (c 1 + c 2 = c 3 ). They give the same relevance to the terrain characteristics as the areas of influence established by the user. In this case, the values considered were c 1 = 0.25, c 2 = 0.25 and c 3 = 0.5. The visual result shows that the different plants appear grouped by species because the area of influence has played the leading role. As indicated in Table 3, 7,581 seeds were initially distributed, with a spacing between the plants of 5 pixels, but after applying the requirements for growth, 2,145 plants were eliminated. This test required the resolution of a system of inequalities made up of 7,581 inequalities.
The second test, called Height relevance in the table, evaluates the visual results when the influence of the height is increased. The values given to the coefficients in this case were c 1 = 0.6, c 2 = 0.2, c 3 = 0.2. The same distribution and number of seeds were set. These data also produced the same number of inequalities as in the previous test. However, finally, more plants were eliminated: 2,476. Figure 9b illustrates the visual result. Plants of different species appear to be more mixed than in the previous test. This is because all of them have a similar range of height where they can grow, as can be seen in Table 1.
The third test was performed by modifying the range of heights of the terrain. Maintaining the variation of 500 meters, the scale of heights was modified from 1,900 to 2,400 meters. Regarding the data shown in Table 3, it can be observed that on establishing the same number of initial seeds, the number of plants eliminated was quite similar to that of previous tests (2,434). Figure 9c shows the visual results.
On the one hand, an area without any plants can be observed. This area is the highest one on the map and so, by increasing the height range above the maximum height where the three species can grow, an empty area appears. On the other hand, the plant represented in yellow, Sorbus, has a considerably increased number of instances because the maximum height where it can grow is higher than the other two vegetable species. In the low-lying areas, however, the dominant vegetable species becomes Pinus Sylvestris (blue), given that its minimum height is lower than that of the red element, Albies Alba, and it has a higher slope value than the other two species.
Finally, it can be seen in this figure that despite the fact that there is less area covered by vegetation, the number of eliminated seeds is similar to the rest of the tests. This is because the dominant vegetable species has a smaller maximum size than the other two and, therefore, more elements fit in the same area than in the other tests with other dominant species.
Regarding the fourth test, the initial spacing between seeds increased twofold, resulting in 10. This has been called "Increasing spacing". In this case, as more space is required, a lower number of seeds were distributed: 1,933. This produced only 9,665 inequalities, which, after solving them, meant that 242 plants eventually died. Figure 9d illustrates the visual result of this test. It shows that if the plants have more space available around them, most of them can grow to their maximum size, thus reducing the number of those eliminated.
The fifth test doubled the number of vegetation species involved in the experiment, resulting in 6 vegetation species. The denomination "Doubling no. species" was applied in the table. Figure 9e shows the visual result of the experiment. In this case, the new vegetation species considered are: Pinus pinea (represented in magenta), Quercus robur (cyan), and Betula pendula (green). The characteristics regarding the radius of the crown are obtained from [36] and the range of the terrain features (height and slope) were considered to be similar to the three vegetation species shown in Table 1. Figure 9e shows how the competition increases and the size of the groups of vegetation species are smaller.
Finally, the last test was performed, doubling the resolution of the terrain. In this case, a resolution of 1024 × 1024 meters was established. On considering a double-sized terrain and maintaining the same spacing between seeds, 5 meters, the number of seeds distributed increased to 35,969, and so the number of inequalities reached 179,845. In the same way, the number of plants that were finally deleted increased to 5,428. Figure 9f shows the visual results obtained by applying the method on this large map. This image shows that there is a greater number of elements, but the competition between them is not very different from the results obtained using the initial map size.

Computational time analysis
The computation time of the method is conditioned by the number of inequalities that have to be solved in order to perform the space-filling process. As stated in Section 6, they are stored as S · r ≤ b, where S is the matrix storing 1 or 0, r is the radii vector to be solved, and b represents the independent terms of the inequalities. Some of the inequalities considered control the space each plant occupies according to its maximum radius, while others control the distance from that plant to the rest of the plants distributed over the terrain. On analyzing the size of this matrix S, it is observed that given a set of n plants, the number of inequalities that should be incorporated if all possible distances are considered is determined by 2 · n + n · (n − 1) 2 (13) However, the number of inequalities can be reduced by considering that the resulting value of the plant radius is directly related to the distance to its neighboring plants instead of to all of them. This approach is a critical aspect in order to reduce the number of inequalities and make the method effective at the computational level.
Several tests were carried out to choose the appropriate number of neighbors without producing overlapping of crowns between plants. According to their results, for each vegetation species only the distances related to 5-nearest neighboring plants are required. This approximation reduces the number of inequalities to the value shown in (14). Nevertheless, it has been demonstrated experimentally that this number can be reduced by up to 40% if the distances that are repeated between neighbors are eliminated.
2 · n + 5 · n (14) Figure 10 shows an example of the distribution of plants of the vegetation species Pinus Sylvestris. The terrain was considered to be a 40 × 60 meter rectangle where the seed placement followed a uniform distribution. Figure 10a illustrates the results of this distribution, where all distances are considered, while in Fig. 10b only the 5-nearest neighbors are taken into account. The Fig. 10b shows that with only the 5-nearest plants, the method associates to each seed a circular area that does not intersect with any other and produces a good growth proposal.
Finally, in order to analyze the computational cost regarding the number of plants, some populations of Pinus Sylvestris were placed over an area of 512 × 512 meters using this assumption of considering the 5-nearest neighbors. Figure 11 analyzes the results of this experiment. It can be observed how the time in seconds tends to grow exponentially with the number of seeds distributed over the terrain.

Visual results
The process of obtaining the final terrain is illustrated in Fig. 12, in accordance with the experiment called Basic test. The terrain has been textured to show the placement of the different vegetation species and the difference in growth between the plants. Colored circles have been used to offer a visual representation of the relationship between the coverage    Figure 12a displays a top view of this coverage map and Fig. 12b shows a perspective view, which allows the characteristics of the terrain to be observed visually.
Lastly, Fig. 12c illustrates the final visual result, by changing the circles for the corresponding trees. The realistic model where each plant is represented allows us to value how the characteristics of the terrain have conditioned the final results. Apart from the initially established areas without vegetation, the feature that has determined the lack of vegetation was the slope. Height was not relevant because the maximum height considered in the experiments (1,300 meters) has allowed growth of all the species in the test.
In order to verify the good performing of the method with real terrain, a similar simulation has been carried out using three plant species and a height map obtained from a GIS (see Fig. 13). Figure 14 illustrates the results in more detail. The Landscape Auto-material plugin has been used in the representation, so the ecosystem generated can be observed better. This Figure shows three captures of the visual results produced at different distances. Figure 14a shows the different vegetation species in the outdoor environment at a long distance. Moving closer to the camera, Fig. 14b shows the scenario at a medium distance. Here, the lack of vegetation can be seen in areas with high slopes. Finally, Fig. 14c illustrates a shot from the ground, where the good quality of the visual results can be seen.

Conclusions
This article presents a procedural modeling method for natural environments that obtain the size of each plant representation using inequality systems. They consider biological restrictions, so if the size obtained does not reach the minimum to live, this plant is eliminated. Following this process, the proposed method obtains realistic visual results. The inequalities that control the size of each plant can be customized, adapting the final representation to the requirements of the final application. For instance, overlapping has been avoided in the work presented, although this restriction can be changed to allow canopies to mix, as occurs in nature.
The ecosystem modeling method proposed here starts by performing an automatic positioning of seeds, which are later labeled as a particular vegetation species. Different variables condition these labels, mainly based on the characteristic of the position of the seeds on the ground and on the biological requirements of each vegetation species. The labels of the k-neighboring seeds are also considered to represent the typical grouping of plants present in nature. Users are allowed to partially condition the final labeling process by establishing certain areas where these groups of plants have to be allocated.
The final realistic representation is obtained by performing the calculation of the final plant representation that grows from the seeds solving the inequality systems. Maximum coverage of the terrain is achieved by allowing the plants to reach the maximum size, and eliminating them if they do not reach the minimum size. Intersection between the volumes assigned to each plant is not allowed in this work, these restrictions being included as inequalities, although some branches can reach other plants in the final realistic representation.
Results show that the outdoor environment generated meets the requirements given for each vegetation species. Furthermore, the visual results are realistic and natural due to the heterogeneity in the size and position of the plants, the grouping achieved between species, and the fact that some elements of another species grow among them.
As future work, this methodology could be extended with the aim of covering large areas of vegetation. The use of Wang-tiles with KD-trees [43] has been considered as a way to solve this problem. These multilevel data structures partition the space into instantiated representations where each tile is filled with PDD. It could easily be extended to use Wangtiles in order to generate patterns of plant distribution.
Another line of research to work on is to improve the method by adding other biological characteristics of the species, such as soil moisture, temperature, and the position of the sun. In addition, the possibility of adapting the plant representation according to the age that is appropriate to each size is being analyzed. Only one representation has been used in this work, but more biological precision would be obtained if different plant models representing different levels of growth were used for each vegetation species.