INTERACTIVE LAND-USE OPTIMIZATION USING LAGUERRE VORONOI DIAGRAM WITH DYNAMIC GENERATING POINT ALLOCATION

: In this work, we devise an efﬁcient method for the land-use optimization problem based on Laguerre Voronoi diagram. Previous Voronoi diagram-based methods are more efﬁcient and more suitable for interactive design than discrete optimization-based method, but, in many cases, their outputs do not satisfy area constraints. To cope with the problem, we propose a force-directed graph drawing algorithm, which automatically allocates generating points of Voronoi diagram to appropriate positions. Then, we construct a Laguerre Voronoi diagram based on these generating points, use linear programs to adjust each cell, and reconstruct the diagram based on the adjustment. We adopt the proposed method to the practical case study of Chiang Mai University’s allocated land for a mixed-use complex. For this case study, compared to other Voronoi diagram-based method, we decrease the land allocation error by 62.557%. Although our computation time is larger than the previous Voronoi-diagram-based method, it is still suitable for interactive design.


INTRODUCTION
Land-use optimization problem (LUOP) is one of the most important problems in computational landscape design.In this problem, we want to divide a land into several parts.Each part must have a specific area, and two adjacent parts must have a strong compatibility.There are many methods proposed to solve the problem.Most of them are based on combinatorial optimization techniques such as simulated annealing (Santé-Riveira et al., 2008), tabu search (Qi et al., 2008), genetic algorithms (Karakostas and Economou, 2014), and a hybrid of those methods (Mohammadi et al., 2016).
Those combinatorial optimization algorithms can output a very satisfiable output for the problem in less than 100 hours.When the problem inputs, which include many design issues, are decided beforehand, this computation time is acceptable.However, in many cases, the design process needs to be interactive (Michalek and Papalambros, 2002).Landscape architects usually cannot design the inputs beforehand.They have a rough idea and use that rough idea as an input of an optimization process.After that, they will iteratively improve their inputs based on the outputs they have got from the process.To have a successful interactive design, our optimization process cannot be longer than a few seconds.We have to consider a method that is much faster than those combinatorial optimization techniques.
For this propose, we can use techniques based on computational geometry techniques such as Voronoi diagram (Karimi et al., 2009), multiplicatively weighted Voronoi diagram (Karimi et al., 2009, Dong, 2008), adaptive multiplicatively weighted Voronoi diagram (Reitsma et al., 2007), and LaguerreVoronoi diagram (Song et al., 2015).In this paper, we choose to develop our methods based on the Laguerre-Voronoi diagram because of two reasons.First, the boundaries of each cell in the diagram are straight lines (Imai et al., 1985), which are more desirable in practice than curves obtained from the multiplicatively weighted Voronoi diagrams.Second, we can adjust the weight of cells based on their area requirements.We can give a larger weight to a cell that requires a larger area, and give a small weight to the others.

Our Contribution
All Voronoi diagram-based approaches assume that generating points, which denote centers of all cells, are given as inputs.Indeed, those point positions are not given in many cases.Placing the points in inappropriate positions results in a partition of which cells do not have an expected area or pair of adjacent cells are not compatible with each other.
In this paper, we propose a method to place the generating points to the best positions as possible.
Our inputs include a network of cells, where two cells are adjacent in the network if they are compatible with each other.Then, we use a force-directed graph drawing algorithm (Bannister et al., 2012) to place the cell in a plane.By the algorithm, a pair of cells that is not compatible with each other is pushed to be apart from each other, and a pair of compatible cells is placed at a close position.
Based on the generating point positions obtained from the graph drawing algorithm, we construct a Laguerre Voronoi diagram in a way that the area of each cell is not so different from the desired area, and the cells are adjacent as desired.We formulate the optimization problem to search the Laguerre Voronoi diagram satisfying the properties which is the multiobjective optimization problem.
We adopt our proposed method to the practical case study of Chiang Mai University's allocated land.The land is called as CMU square, and its area is 140,800 sq.m. (Prachachat Online, 2015).In 2015, the university plan to develop the land into a mixed-use complex, which contains areas of arts, culture, nature, and community.Although the land allocation has been finalized, we think that the compatibility between adjacent areas can be improved.
We use the center of each cell in the finalized land allocation as our baseline generating point positions.Compared to the baseline, we decrease the land allocation error, which indicates how much the areas of each cell are different from the desired areas, by 62.557%.In addition to that, the compatibility value is increased by 28.010%.
Our computation time is 839.78 s in Matlab.Because computation times in Matlab is usually more than a hundred times larger than those in C or C++ (Andrews, 2012), we strongly believe that our computation time will reduce to a few seconds when our software is implemented in those languages.Hence, the algorithm is fast enough for the interactive design.

Paper Organization
This paper is organized as follows: In Section 2, we will explain the techniques which we use in this work, Laguerre Voronoi diagram and force-directed graph drawing.Then, in Section 3, we will formally state our problem definition.In Section 4, we will give our methods, and, in Section 5, we will show the experimental results of our case study.We then conclude the paper in Section 6.

Laguerre Voronoi Diagram
The concepts of Voronoi diagram have been widely studied in many areas as complied in (Okabe et al., 2009).There are many generalizations of Voronoi diagram.One of the generalizations is to put weights to the generators of the Voronoi diagram such as additively and multiplicatively weighted Voronoi diagram.In this study, we mainly focus on the use of Laguerre Voronoi diagram, the special case of the additively weighted Voronoi diagram, which was introduced by (Imai et al., 1985, Aurenhammer, 1987).Let S = {c1, ..., cn} be a set of generating circles in R 2 such that ci = (p i , ri), where p i = (xi, yi) is the position of generator and ri is it weight, which is interpreted as a circle in the plane.The Laguerre distance (or recognized as power distance) is defined by the distance For each adjacent circles ci and cj, Laguerre bisector is defined by which is a straight line in R 2 .The half space of pi dominated by the generator pj is and finally the Laguerre Voronoi region is defined by L(pi) = ∩ j∈In\{i} H(pi, pj).Therefore, the Laguerre Voronoi diagram can be constructed from the set of region L = {L(p1), ..., L(pn)} including their boundaries.The properties were widely studied as shown in (Okabe et al., 2009), and the robust algorithm was presented in (Sugihara, 2000).
In the case of Laguerre Voronoi diagram, remark that cells of the corresponding generators may be lost, or generating points may not be contained in the Voronoi cell as mentioned in (Reitsma et al., 2007).However, (Cheng et al., 2000) gave the necessary condition for keeping the generator positions lay inside the cell by the following theorem.
Theorem 1. (Cheng et al., 2000) For a set S of points, let w(p) ≥ 0 be the weight of point p ∈ S. If, for all q, d(p, q) ≥ 2 • w(p), then no cell of the Laguerre Voronoi diagram is lost, and each point of S lays inside its cell.
Remark that the Laguerre Voronoi diagram has its dual structure called the Laguerre triangulation, which can be considered as a planar graph.If a cell L(pi) is adjacent to cell L(pj), there exists an edge {pi, pj} in the graph whose the vertex set is {p1, ..., pn}.

Force-Directed Graph Drawing
Graph drawing is a technique to visualize a given network (V = {v1, . . ., vn}, E) on a 2-dimensional plane by a map function f : V → R 2 .To find the best visualization for each network, one can define energy for all map functions, and search for a function that minimizes the defined energy.To have generating point positions that satisfies Theorem 1, in this paper, we choose the energy function defined in (Kamps et al., 1995).
The energy of each function is defined based on two types of energy.Given a parameter α, we wants the Euclidean distance between vi and vj to be α • d(vi, vj) when d(vi, vj) is a shortest path length between vi and vj in (V, E).The energy between vi and vj, which is defined by the expected Euclidean distance, can be defined as follows: ρi,j(f where and the energy becomes large when the difference between two distances becomes large. The first type of energy for mapping function f , defined as ρ(f ), can be defined as follows: The second type of energy comes from the fact that we draw each node vi as a rectangle on the plane, and we want to have an overlap area between two rectangles to be as small as possible.We assume that the size of a rectangle for vi is given as si, and Rec(f (vi), si) be a rectangle centered at f (vi) with size si.Given a parameter ω, the energy between vi and vj, which is defined by the overlap area, can be defined as follows: The second type of energy, γ(f ), is defined as follows: We are now ready to define the energy function, denoted by E(f ), in the following equation: In most graph drawing settings, there are some points with given fixed positions.Let V ⊆ V be the set, and f : V → R 2 be a function that give those fixed positions.The set of possible function, denoted by F, is Because the function E is differentiable at most points, to find f ∈ F such that E(f ) is minimized, it is shown in the paper that the gradient descent method can be used.

Nelder-Mead Method for Optimization
The Nelder-Mead method is a technique for solving the optimization problem proposed by Nelder and Mead (Nelder and Mead, 1965), which is known as a derivative-free optimization method.
In brief, for a function f : R d → R which is supposed to find the minimum of f (X) where X ∈ R d , a simplex consisting of d + 1 vertices of dimension d is chosen.In the optimization processes, the simplex is expected to move and converge to the local minimum.
From a set of simplex vertices {x1, ..., x d+1 }, the set {f (xi)} d+1 i=1 is ordered and relabelled with respect to the sequence f (x1) ≤ ... ≤ f (x d+1 ).After ordering and relabelling, x1 is said to be the best vertex and x d+1 the worst vertex.From the worst vertex, we try to replace the worst vertex with a new vertex which is acquired from the reflection, expansion, or contraction through the direction of the centroid among the remaining points, or shrinkage all vertices (except for the best vertex) to the direction of the best vertex.The parameters for reflection, expansion, contraction, and shrinkage are written as {α, β, γ, δ}, and normally set as {α, β, γ, δ} = {1, 2, 1/2, 1/2}.
It is remarkable that the Nelder-Mead method is easy to implement.The convergent properties of the Nelder-Mead in low dimension was proved in (Lagarias et al., 1998).

PROBLEM DEFINITION
In this section, we will give a formal definition for our LUOP.We assume that a designer gives us a landscape object of region (V, E).Each vi ∈ V denotes an object, and, for any v1, v2 ∈ V , {v1, v2} ∈ E if v1 and v2 are compatible with each other.For each vi, the designer also gives us an expected area of the region vi, denoted by A * i , and a land before an allocation R.
Our input defined in the previous paragraph are a particular case of the input of the LUOP in (Mohammadi et al., 2016).In that paper, instead of defining if there is a compatibility between two areas, they give a compatibility value between all pairs of areas.Although their input provides more flexibilities to designers, putting values for all pairs of areas is a tedious task, and is not possible for the interactive design.Because of that, we limit the values to only 0 and 1, and make the inputs be a network, which can be entered interactively.
We aim to find a partition of regions R = {R1, . . ., Rn}, where Ri is a region corresponding to object vi.The regions must cover a whole land area, i.e.R = n i=1 Ri.Two regions are adjacent to each other if they share some boundary lines.
Definition 1 (Compatibility of Partition).Let Ri ⊆ R be a set of regions that is adjacent to Ri.A compatibility of R, denoted by C(R), can be defined as follows: The previous compatibility definition is slightly different from the definition given in (Mohammadi et al., 2016).In the previous definition, the geographic properties of the regions, such as lengths of lines that two regions share each other, are also considered.We do not consider those geographic properties because calculating those properties is time-consuming.If we have to calculate their compatibility value, our method might not be efficient enough for interactive design.Also, we strongly believe that, when two landscape objects are allocated adjacent to each other, users of the objects can easily walk between them.Having them sharing a long line does not make the walk between two regions more convenient.
Next, we will define the land allocation error of partition R := {R1, . . ., Rn} in the following definition.The definition is same as that defined in (Karimi et al., 2009).
Definition 2 (Land Allocation Error).Recall that the area of Ri is |Ri|.The land allocation error of R, denoted by err(R), as follows: In this paper, we aim to find a partition of regions that maximize the compatibility and minimize the land allocation error.

Graph Drawing-Based Point Placement
Our point placement method is based on the graph drawing techniques explained in Section 2.2.In the following subsection, we construct Laguerre Voronoi diagram from f (v1), . . ., f (vn) obtained from the techniques.Let wi be the weight corresponding to vi in the construction, and Ri be the region obtained.Because we do not want Ri to be an empty region, by Theorem 1, we require that the distance from f (vi) to any other points have to be at least 2 • wi.We can have that requirement, when no rectangle overlap with Rec(f (vi), 4 • wi).Because of that, we choose si in Section 2.2 to 4 • wi, and, to guarantee that the overlapping does not happen, we set the parameter ω to a very large real number.
On the other hand, during the interactive design, we ask designers to choose the appropriate α for their design.A larger α will make the generating points spread more on the plane.
While there is not a drawing boundary in the graph drawing technique, we have to place all points in a given boundary.We update the energy function E(f ) to include the constraint in our point placement.Let 1, . . ., m be line segments that define the boundary, and let D(f (vi), j ) be a shortest distance from f (vi) to j .Given a parameter θ, we define another type of energy for vi, denoted by φi(f ), as follows: The energy of φi(f ) is very large when the point f (vi) is very close to some boundary.The third type of energy φ(f ) and the updated energy function E is defined as in the following equation: Similar to function E, function E is also derivable, and we can use the gradient descent method to search for f ∈ F that minimize E .

Laguerre Approximation
Let S = {c1, ..., cn} be a set of circles such that ci = (f (vi), wi), where f (vi) is a of point which is determined by the graph drawing in section 4.1, and wi is a weight corresponding to the position f (vi).
In this process, we would like to find the Laguerre Voronoi diagram R = {R1, ..., Rn} satisfying the following properties.
P1 In the set E denoted from a designer, the Laguerre triangulation of R contains edges of E as much as possible.
P2 Area of Ri is not much different to the expected area A * i ; Therefore, the search of the Laguerre Voronoi diagram which satisfies the mentioned properties is necessary.Therefore, we formulate this issue to the optimization problem.
For satisfying the property [P1], we use the compatibility as defined in Definition 1 for determining the connection between any two cells as the desired structure in the graph {V, E} from the designer.Also, we consider the property [P2] when the i-th area of Laguerre Voronoi diagram is computed as |Ri| which is the area of the convex polygon.Then we use the criterion of land allocation error as defined in Definition 2 to the Laguerre Voronoi diagram R.
Since the Laguerre Voronoi diagram is constructed with respect to the set of of generators S, we can consider the compatablity and land allocation error as a function of generator set, i.e.C(S) and err(S) is a function from R 3d to R defined by. and We would like to find the set S * which maximizes C(S) which is equivalent to the minimization problem min −C(S), and minimizes err(S).This problem leads to the multi-objective optimization, that is, to find In this study, we solve the multi-optimization problem (15) using the weighting method.Therefore, we define a function and the minimization is done by solving min F (S).
We treat the function F (S) defined in ( 16) as the black-box function, i.e. we can implicitly compute the value of the function without knowing the explicit formula.Therefore, we solve the optimization problem using the derivative-free method.In this study, we will use the Nelder-Mead method for finding the minimum of F (S).

CASE STUDY
Our application of the experiment will be conducted on the plot of land in Chiang Mai, Thailand, owned by Chiang Mai University.
Originally, the land is utilized by the Faculty of Agriculture for education and experiments, however, will later be allocated as an area to promote education and Lanna's historical heritages.This project is called CMU Square as announced in (Prachachat Online, 2015).
Along with the original programmatic scheme design, the proposed master plan seems like the design went through without any regards to the surrounding context as shown in Figure 1.Hence, this is where the external anchor points on the graph come to play.To create a contextually related master plan, one must be mindful of what the surrounding contents are.The selected site is located between various important landmarks of Chiang Mai, whether it be the Chiang Mai University to the west, the university's hall to the east, and even the famous Nimmanhaemin road towards the north.With multiple characters surrounding the site, it is safe to assume the local demographic of the users, as well as where and how will they access the area.

Our Assumptions
Our method's objective is to maximize the usable spaces while minimizing hard to reach areas resulting from a long clean-cut zoning plan.While the original design portrays a solid and organized, there are certain portions of the zoning that could be underutilized.On the other hand, our experiment method will allow us to optimize these space into polygons surrounding the center point of the zoning which would be the structure, meaning that the circulation distance within the zoning would be equally distributed rather than taking a long distance walk from a corner of a rectangle to another.
Using our proposed method, we define the graph G = (V, E) consisting of the vertices represented as the facilities and edges The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W7, 2017 ISPRS Geospatial Week 2017, 18-22 September 2017, Wuhan, China  2 where the shaded area is the area for allocation.We assume that the area for allocation is a convex polygon written as P .
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W7, 2017 ISPRS Geospatial Week 2017, 18-22 September 2017, Wuhan, China In this experiment, there were 21 vertices in the graph whose seven vertices represented as the existing buildings.Five vertices among all fixed vertices are fixed points laid outside the shaded area.The weight wi for point i is computed by wi = A * i /π.In the case of outer fixed points, we assume that weights are equal to zero.
We follow the process discussed in Section 2.2 and 4.1 to find an appropriate generator positions.After that, we used the resulting positions in the optimization processes using the Nelder-Mead method.In this case, the optimizations were done among all non-fixed points with parameters of generator positions and their weights.In the case of fixed points, we only optimized radii of the inner fixed-points.In the step of the Laguerre Voronoi diagram construction, since the Laguerre Voronoi regions of outer fixed points may affect the Laguerre Voronoi diagram cell inside the polygon P , the construction of Laguerre Voronoi diagram were done only the vertices of the polygon P .
Hence, this experiment considered the function F (S) as a function from R 44 to R when the parameter ω was chosen.In this sense, ω was determined to show the priority between the compatibility and land allocation error.
We optimized F (S) using Nelder-Mead method with a chosen parameter ω.Remark that the simplex contained 45 vertices.For the initial simplex, the first vertex S0 ∈ R 44 was chosen from the resulting parameters of graph drawing.The remaining vertices were chosen from by Si = S0 + kSR, where each coordinate sj of SR was a random number from 0 to 1, and k was chosen for scaling the size of the simplex.In this study, the Nelder-Mead method was terminated when the standard deviation of computed function values of all simplex vertices was smaller than a fixed threshold .

Experiments using the proposed framework
In the experiment, coordinate of points in the initial step were located using GeoGebra 5.0.The graph drawing was implemented using Java.The construction of the Laguerre Voronoi diagram and Nelder-Mead optimization were done using the software MATLAB R R2017a with libraries 'Bounded Power Diagram' and 'Nelder-Mead Simplex Optimization' published in file exchange system of MATLAB R .
During the graph drawing process, we set α = 0.5, ω = 20.0, and θ = 0.5.For the Nelder-Mead optimization to minimize F (S), we set ω = 1/8 with the simplex size k = 0.1.The initial function value from the initial point S0 is 15.606.The land allocation error value was 18.170, and the compatibility was 2.3417.
After 4,915 iterations with threshold = 10 −5 , the function value of F (S) converges to 5.5773, as shown in Figure 4.The land allocation error was 6.8033, which decreased by 62.557%, and the compatibility was 2.9976 which increased by 28.010%.The optimization with Nelder-Mead method took 839.78 seconds.The resulting diagram plotted from the converged value of the function is shown in Figure 3.
The experimental result shows that we can automatically locate the facilities using graph drawing schemes.Also, we can allocate the area to facilities using the Laguerre Voronoi diagram in a way that land allocation error is minimized, and the compatibility is maximized.
Remark that the computational time can be shortened when we use other compilers such as C or C++, as discussed in (Andrews, 2012).In many cases, the computations in C or C++ is more than a hundred time faster than those in Matlab.Hence, we strongly believe that our computation time will decrease to a few seconds if we implement our software in those languages.

Comparison with the previous studies
We validate how useful our point allocation method by comparing our result with point positions extracted from the master plan of CMU square (Prachachat Online, 2015), which is designed by a professional architect.The point positions are shown in Figure 5.Then, we generate Laguerre Voronoi diagram and multiplicative Voronoi diagram from the point positions in the master plan, and compare their land allocation error and compatibility with our results.
When we construct Laguerre Voronoi diagram from the point positions in the master plan, the land allocation error was 17.1373, which was worse than the result from the framework by 86.3037%.On the other hand, the compatibility was 3.5881, which was better than the optimal result 17.9328% difference.Although the compatibility obtained from the professional architect is slightly larger than our automatically generated point positions, their land allocation error is much larger than ours.Because the whole area is about 42.21, we strongly believe that the land allocation error of 17.14, more than 30% of the whole area, is not acceptable.
The results are even worse in the multiplicative Voronoi diagram.The compatibility obtained in the setting is 4.0, 28.6498% better than ours and slight better than Laguerre Voronoi diagram's.However, the land allocation error is as unacceptable as 33.0.By the multiplicative Voronoi diagram, many areas are significantly larger than the others.For example, we found that the region for Mall1 occupies almost all areas in the upper half of the space.We obtain a large compatibility.because those large regions connect to almost every regions, not because the large value do not indicate that the point positions follows the expected topology.
Although the experimental results cannot be used to judge whether or not our proposed method is better than the previous methods, the tendency of the results shows that the dynamic generating point allocation method together with the optimization can provide the smaller land allocation error as well.

CONCLUSION
Ultimately, the experiment is conducted to acquire the balance between maximizing the functionality of the land and the aesthetics of the design as a whole.Tackling the design project with computer-aided design software would reduce the time required for creating a complex design while still maintaining human errors to a minimum.
The allocation of functions, initially done by hand, is gradually optimized through the process by a design tool to maximize the compatibility for the adjacent land use.
For a further experiment, our focus is the development of simple but useful software to instantly allocate land uses according to their contributing factors such as the external context, internal relationship, and size.In parallel, a survey shall be conducted with a selected target group to acquire inputs for further development.

Figure 1 .Figure 2 .
Figure 1.The master plan of CMU square from a newspaper.The figure is adapted from Google R map.

Figure 3
Figure 3. (Left) The Laguerre Voronoi diagram of the generator set generated from the graph drawing-based point placement; (Right) The resulting Laguerre Voronoi diagram after employing the optimization with parameter ω = 7/8

Figure 4 .
Figure 4.The graph shows the convergence of the function F (S) using the Nelder-Mead method.