Abstract

Massive sport, such as unmanned aerial vehicle performance, often needs fast and efficient calculation of formation morphing and individual path planning. This paper introduces a novel fast formation control method of a crowd. First, we get the agents’ location in a 2D polygon with centroidal Voronoi tessellation and L-BFGS techniques. Then, we transform crowd formation shapes with a global shortest motion path pair assignment using earth mover’s distance algorithm. Finally, the repulsing force between agents and obstacles is calculated based on the recursive velocity observer method control agents’ motion. Extensive experimental results show the effectiveness and usefulness of our algorithm in 2D group formation transformation.

1. Introduction

Crowd animation simulation, which is a fundamental problem in computer graphics, augmented, virtual, and mixed reality, has received much attention in recent years [1, 2]. With the rapid development of computer hardware and software technology, especially the continuous improvement of GPU rendering technology, simulating animation of a massive number of virtual agents on a limited computing capability hardware is a problem that has recently become quite popular and crowd formation control has become one of the most important issues in the research field of Multirobot System (MRS). It plays an important role in a wide range of applications, including industry, military affairs, aviation, military drills, urban planning, art and sports, and video games. Especially, the manufacturing industry has an increasing demand for automatic processing of large objects or dangerous goods. In the military, the power and command efficiency can be improved by maintaining the military equipment or the army as reasonable and appropriate formation by the group simulation system. Agents usually interact with the user, and it not only improves operation efficiency but also ensures human safety. With the help of virtual agents, we can effectively complete the task and ensure the safety of humans in dangerous activities, such as exploring unknown environments, demining, patrolling, and reconnaissance. However, the crowd formation transform is still quite limited in terms of both agent number and individual intelligence. Designing and rehearsal of group formation are tedious and time-consuming; the crowd formation transformation simulation system is an effective solution for this problem. Through team formation control, formation design, and crowd simulation, the simulation system can combine with art and sports perfectly.

There are three major classes of algorithms for crowd formation control, namely, formation generation, formation maintenance, and formation transformation. The main difficulties are collision avoidance and path planning. If the motion matching between individuals is not accurate, the collision avoidance between individuals is difficult to be realized. It is complex and unrealistic to use artificial methods to realize the matching. Leader-follower dynamic technique [3] uses a linearization method to design a velocity controller for transforming, but it depends on leader agents. When the leader follower fails, the formation will be out of control. Zheng et al. [4] proposed a group formation control method for formation based on geometric constraints, which can be applied to perform agent pairing, path planning, and collision avoidance for the crowd; the calculating process still needs a tedious iterative step and is time-consuming in real-time simulation.

To form a uniform group formation distribution, efficient calculation of individual motion matching, and smooth formation transformation, we propose a novel group formation transform approach to automatically generate a homogeneously distributed crowd and visually pleasing formation transformation. The group is uniformly sampled by the centroid Voronoi structure centroidal Voronoi tessellation (CVT). Based on the earth mover’s distance (EMD) algorithm, we precompute a global shortest motion path assignment in the offline process. During the real-time transformation simulation, our contributions are summarized as follows:(1)We introduce a novel group formation morphing approach that can automatically generate a global shortest motion path pair assignment based on geometric constraints(2)We introduce a new method to enable a swarm of insects to follow the motion of a moving shape by applying EMD method

2.1. Group Formation Morphing

Group formation control is a vital problem of many crowds and multiagent systems. Based on the Velocity Obstacle concept [5], Van den Berg et al. [6] proposed the concept of Reciprocal Velocity Obstacles for real-time multiagent navigation among autonomous decision-making entities, and it is widely used in the crowd simulation system. Anderson et al. [7] introduced a flocks animation simulation method with shape constraints to improve the motion concerning the underlying behavioral model. Ji et al. [8] proposed a method to simulate the movement modeling, path planning, and collision avoidance, which provided an auxiliary tool for the group calisthenics to improve the design of the formation and pattern. Zheng et al. [9] adopt geometric constraint idea for smooth formation animation of regulated crowds. Chen et al. [10] presented a robust multiagent model for large-scale controllable shape constrained simulation of flying insects. Wang et al. [11] proposed a linear formation control method under directed and switching topologies for distributed robot systems. Zheng et al. [4] proposed a geometric constraint group formation control method for group calisthenics, which is combined morphing, Centroidal Voronoi Tessellation, and Lloyd techniques and it is an efficient solution for individual distribution layout, individual matching, path planning, and collision avoidance. Furthermore, Zheng et al. [12] introduced a regulated heterogeneous flock formation control method based on a capacity-constrained centroidal power diagram to model the internal layouts of groups in the constraint shapes. Some researchers convert the formation problem into tracking problem, use the pseudorigid body method to calculate the path of the expected point of the agent, and use the reverse step method to design the formation control of the nonlinear control algorithm [13]. Based on backstepping control and bioheuristic neural network, Ding et al. [14] used a trajectory tracking method for three-dimensional formation control and obstacle avoidance. Henry et al. [15, 16] proposed a crowd control algorithm using deformed grids for single traversal search, and it can be used to control crowd interactive with the environment and collision avoidance. Xu et al. [1] proposed a natural group transformation and control method between source shape and target shape based on Kuhn–Munkres and mutual information method, which can be used for the interaction and collision avoidance of between individuals and environments.

2.2. Voronoi Diagram

In the centroidal Voronoi diagram, the distribution of seeds is very uniform. Due to these properties, centroidal Voronoi graphs can be used in a wide range of applications, including image compression, numerical integration, finite difference, location problem, and meshing [17]. The Voronoi diagram is suitable for not only uniform sampling but also sampling with constraints; it can be used to calculate a distribution as uniform as possible and compute region partition with capacity constraint. Based on the quasi-Newton method, Liu et al. [18] proposed a superlinear convergence algorithm for centroidal Voronoi diagram and successfully applied it to the important problem of remeshing. Yan et al. [19] proposed a method to compute the centroidal Voronoi diagram limited to the curved surface, which improved the efficiency of remeshing. After that, the centroidal Voronoi diagram has been widely studied in digital geometry. Lévy and Liu [20] introduced a centroidal Voronoi diagram method based on p-norm and demonstrated its application in quadrilateral meshing. Pan et al. [21] used the centroidal Voronoi diagram to construct grid surfaces with constant mean curvature. Wang et al. [22] proposed a fast algorithm for calculating the Riemannian center and the center of mass for any geodesic Voronoi diagram. Recently, Liu et al. [23] further improved the calculating speed of the geodesic centroid Voronoi diagram based on the evolutionary algorithm. To compute the Voronoi diagram with capacity constraint, Sorkine et al. [24] and Zheng et al. [25, 26] proposed some new optimization strategies, respectively, and their algorithms have obvious improvement in computing efficiency. The distribution problem with a capacity constraint can be interpreted by the optimal transmission theory, Desbrun et al. [27] provided a first-order convergent optimization algorithm and found that the blue noise generated by the capacity constraint has better spectral characteristics. Xin et al. [28] named this problem as the centroid power diagram problem and proposed an optimization algorithm with superlinear convergence.

3. Our Methods

Figure 1 presents an overview of the pipeline of our approach. Specifically, it can be divided into four phases: (1) sampling number calculation, (2) group distribution sampling, (3) pair assignment precalculation, and (4) movement control phase. The sampling number calculation phase calculates the number of samples for each shape in a multishape formation. The group distribution sampling phase computes the sampling position based on the centroid Voronoi diagram. Pair assignment precalculation phase formulates the globally shortest path matching as a linear assignment optimization problem. The movement control phase combines target-oriented motion and interaction force between nearest agents to ensure collision-free movements for all the individual agents in real time.

3.1. Sampling Number Calculation

For the source and target formations, we first convert the formation shapes to Delaunay triangulation (DT) representations [29] in the interior of the shape. Then, we calculate the area of each triangulated subshape for a formation with multiple shapes. Finally, we use the followed formula to calculate the number of sampled agents in each triangulated subshape:where is the total number of sampled agents and is the number of subshapes.

3.2. Group Distribution Sampling

The Voronoi diagram is a partitioning of a given domain into regions based on distance to a set of points in a specific subset of the domain. That set of seeds is specified beforehand, and the region containing point is called the Voronoi cells; for each seed, there is a Voronoi cell consisting of all points closer to that seed than to any other. represents the discrete point set on the domain ; the Voronoi cell of the site can be denoted as follows:

Let the domain be endowed with a density function which is bounded; then, the centroid of each Voronoi cell is given bythat is, the points that serve as generators for the Voronoi regions are themselves the mass centroids of those regions. We call such a region a centroidal Voronoi diagram. Then, the CVT energy function of the Voronoi tessellation is defined as [17]

In this case, minimizing the energy function , we can get all the centroids of the tessellation which directly ensures that each position is a representative of . Liu et al. [18] and Xin et al. [28] use gradient descent method to calculate the center of CVT rapidly; the gradient of the energy with regard to is defined aswhere is the centroid and is the mass of Voronoi regions .

It is worth noting that all point sampling is located inside of shape using the above CVT method. However, groups who are strictly distributed according to user-designed shapes need to evenly distribute a part of the sampling on the shape boundary to achieve a more orderly crowd distribution in crowd formation transform simulation. Therefore, the sampling points of the group formation need to be calculated in two parts. In the first part, you need to sample points on the boundary; namely, it needs to generate points on the boundary of the shape and the Voronoi region of denoted as . The remaining of the sampling points need to be distributed in the internal region of shape and represents the Voronoi region of , as Figure 2 shows. The whole energy function of CVT is defined aswhose gradient with regard to iswhere , , is the surface’s area of shape, and is the boundary’s circumference of shape. We are able to adopt the L-BFGS solver [30] to minimize the objective function, leading to empirically superlinear convergence rate. The pseudocode of uniform sampling in shape with CVT is shown in Algorithm 1

Input: collective crowd shape; an error tolerance .
Output: group sampling points .
Generate the initial points as the seed;
Compute the Voronoi region for each seed
Compute the centroid of each Voronoi region ;
While do
L-BFGS
 Compute the energy function in equation (6)
 Compute the gradient w.r.t. x; see equation (7)
 Update x by step search
end
return the final centroid points
3.3. Position Assignment

The goal of group formation transformation is to change all the agent from source group formation to target group formation. The formation transformation not only realizes collision avoidance during crowd movement but also achieves the global shortest moving distance. In group formation transformation simulation, agents need to observe the movement and position of other agents, maintain the coordination of the team, and avoid the free and scattered individual behavior. Therefore, it is necessary to realize the optimal movement control of the group and keep the number of the crowds unchanged. We apply the classical earth mover’s distance (EMD) algorithm [31] to solve the global position assignment problem. The group formation transformation objective function is defined aswhere is the -th position of from source group, is the -th position of target group, is the distance from to its nearest neighboring point, and is the distance from to its nearest neighboring point. The limiting condition of objective optimization function is subject towhere denotes that each agent in the source group can only move to one target position and means that a position of the target group can only be occupied by one agent, as shown in Figure 3.

3.4. Movement Control

In movement control phase, we can also mix our method with other methods for agent-based planning [6, 32]. Recursive Velocity Observer (RVO) method is used to calculate the actual speed of the agent. In addition to collision avoidance, the agent may change the target exit according to the updated environmental information or change the path through the area and then need to replan the path. For each time step, we first advance the individually controlled agents and then treat them as obstacles in the continuum model so that the crowd motion will flow around them.

4. Results and Performance

We implemented and experimented with our algorithm on a computer with a 64-bit version of Win7 system, a 3.07 GHz Intel(R) Core(TM) i7 CPU, and 6-GB memory. The coding language is C++. We shall use extensive experimental results to exhibit and demonstrate the algorithm’s efficiency, as well as its insensitiveness to scene resolution. To demonstrate the effectiveness and practicability of our method, we carry out the evacuation simulation in different scenarios.

In our method, the most computationally time-consuming step is the position assignment, which has the time complexity of ; we precalculate the result of the position assignment in the offline process. We recorded the time cost for the different simulation scenarios (Table 1); the numbers of agents are , and 1000, respectively. As shown in Table 1, the cost time difference is very small on a different shape which contains the same number of agents. For example, the solving time of the objective optimization function is about 3 seconds for pair assignment of a group with 200 agents, while the cost time is 5 seconds and 45 seconds for groups which contain 300 and 1000 agents, respectively. It is suitable for offline formation preordering operations.

4.1. Performance

Figure 4 illustrates the process of motion transformation of different groups in a single shape. Figure 4(a) shows that the “A” word shape formation transformed into “H” word formation and then into “H” word formation with 200 agents. Figure 4(b) shows a 300-person square formation transformed into a five-star formation and then into a circular formation. Figures 4(c) and 4(d) illustrate the group with 400-person formation from pig shape to rabbit shape and rabbit shape to skiing logo shape, respectively.

Figure 5 shows the process of motion transformation of crowds in multiple shapes. Figure 5(a) shows that the group has 200-agent formation from “2018” word with 4 subshapes to the Olympic rings logo which contains 5 subshapes. Figure 5(b) shows the transform process from the shapes of English word “welcome” to the shapes of Chinese word “welcome.” From Figure 5, we can observe that our method has good adaptability in nonconvex and irregular formations with multiple shapes, and the transformation process is smooth and without any confusion.

Figure 6 shows the process of formation transformation with 400 agents in diamond shape under nonuniform density. We set the shape with an e density function which is defined as [4] . Then, the diamond can also be decomposed into four symmetrical shapes. It can be seen that our algorithm can get different formation transformation effects.

4.2. Evaluations and Comparisons

To quantitatively analyze the effect of a crowd formation transformation, we introduce two different objective measures, as described as follows.

4.2.1. The Rangeability of Sampling Points

For agent , the neighboring points set of within its the nuclear radius is denoted as . Let be its the minimum distance from position to its neighbors , and is the current position of in the target group . The source position will move to its target position . The standard deviation and the average value can be used to quantify the stability rangeability of source group sampling during the transformation:

Clearly, the lower value of indicates that the agents have more similar distances from their neighbors.

4.2.2. Transform Distance from Source to Target Group

We also employ the standard deviation and the average value of transform distance to quantify the motion effect of agents. For an agent , its transform distance is defined as follows:

Here, and denote its source position and its current position, respectively. The average value of all the agents’ transform distance is defined as

The standard deviation of transform distance is calculated as follows:

Clearly, the lower value of and indicates that the agents have less time consumption during the crowd transformation process.

We compare our method with the other three crowd formation transform simulation algorithms including RVO [6], geometry-constrained mechanism method [4], and mutual information-based method [1]. As shown in Figure 7, the four different methods are applied to the same source and target formations. The visual results show that our approach can achieve the best visual effect among the four approaches. Table 2 shows the comparison of results from our method on crowd simulation and other algorithms using simulation time, , , , and . From this table, we can observe that our simulation method on crowd formation transform outperforms other approaches no matter what kind of evaluation metric is used, except the Zheng method. When comparing with the Zheng method, our algorithm takes less simulation time.

5. Conclusions

In this paper, we introduce a novel crowd formation transform method based on Earth mover’s distance. The main precompute step of our approach is to sample shapes with CVT technique and use EMD to calculate the global shortest position assignment. In the real-time simulation, the interactions between individuals and obstacles are calculated based on RVO which realize the collision avoidance effectively. Large experimental results show that our algorithm has an obvious speed advantage in the real-time simulation process, and our method can also be used to realize the formation transformation in three dimensions.

For large-scale group, formation transformation simulation, the precomputation process, and the collision avoidance calculation process are still time-consuming. In the future, we will keep improving the performance by advanced parallelization techniques. We will also explore more applications in visual computing.

Data Availability

The simulation data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

All authors contributed equally to this work.

Acknowledgments

This work was supported by NSF of Zhejiang (LY19F020001 and LQ18F010008), NSF of Ningbo city (2017A610115, 2018A610092, and 2019C50008), Opening Foundation of Zhejiang Mobile Network Application Technology Key Laboratory (F2018007), NSF of Huzhou city (2016YZ02), Science and Technology Plan Project of Zhejiang (LGG18F020001), and Social Sciences Project of Zhejiang (19NDJC0103YB).