CVT 1D Milestone

http://people.sc.fsu.edu/~jburkardt/classes/urop_2016/cvt_1d_milestone.html

So far, we have only used the information from a Voronoi diagram of a set of points G as a way to illustrate how the space in a region R is divided up among the points. Now, we want to use that information to try to adjust the G points so that the division is more even.

It turns out that a mathematically reasonable definition of an even subdivision of a region requires that the points G each be the centroid of the subregion they induce. This is an easy condition to check, but not so easy to enforce. However, if we are willing to carry out an interative process, we can gradually modify our point set G to approximate this special arrangement, which is known as a centroidal Voronoi Tessellation, or, because no one wants to say that twice, a "CVT".

As we have already seen in the Voronoi notebook, if we are working in a 1-dimensional region, then it is possible to compute the exact Voronoi diagram information. Therefore, it will be convenient to start our discussion of CVT's by working on a 1D problem.

Let's assume that we have constructed a Voronoi diagram for a set of generator points G that are inside a region R which is actually the unit interval [0,1]. This will define, for each point G(i), a Voronoi polygon, or in this case, simply a subinterval, P(i) of points, extending from A(i) to B(i), that are closer to G(i) than to any other generator.

We can imagine that every point in the region wants to be as close as possible to a generator. A reasonable measure of closeness is the square of the distance. In order to consider all the points in subinterval P(i), we have to compute the integral of this squared distance, which we label E(i), called the "energy" of subinterval P(i). If we sum up the energies over all the subintervals, we have the total energy E for this arrangement.

Clearly, except in a few very peculiar cases, E will not be zero. Moreover, the value of E clearly depends on the locations of the generators G. And it's interesting to wonder whether we can reduce E by using a different set of generators. But in fact, this is easy to do. For each subinterval P(i), if we simply replace G(i) by the centroid of P(i), you can show that the value of E(i) goes down, and so the total energy E goes down. Hence, this arrangement reduces the clustering energy.

But we can reduce the energy even more; because we moved the points, we should recompute the Voronoi diagram, and this effectively guarantees that some points reduce their energy simply because the new generator locations are closer to them.

It turns out that by repeatedly using points G to compute Voronoi subintervals, and using the centroids of the Voronoi subintervals to replace G, the sequence of points G approaches a centroidal Voronoi Tessellation (CVT). Such a set of points G has the property that when you use the points G to compute the Voronoi diagram, the centroids of the subintervals are identical to G.

Of course, having made this sound like a hard problem, it turns out that we can see the answer right away. If all the subintervals should be the same size, and the points G should be in the centers, then we just divide the unit interval into equal subintervals and place a G at the center of each. If we use 5 points, for instance, then G = (1/10, 3/10, 5/10, 7/10, 9/10) and each subinterval is 1/5 of a unit in width.

OK, we know this, but the program we write doesn't know that. Let's start with a random set of G's, and see how fast the iteration is able to approximate the answer we already know!

Topics:

• A single step of the CVT iteration accepts the current generators G, computes the subintervals [A,B], and returns the centroids C. Write a Python function cvt_1d_step() to do this, assuming that the region R is [0,1], and that the points G are sorted, and in ascending order.
• Write a Python function cvt_1d() which accepts a starting set of generator points G, sorts G, applies a number of CVT steps, replacing G the values of C, and when done, returns in G the improved values. This function should accept as input the points G, and the number of steps to take.
• Run a sample calculation by picking 10 random G values to start, and running 20 CVT steps. Compare your answers to the exact solution we already know.
• Modify cvt_1d by computing the total energy at the start, and after each step. Return this as an additional output argument E. Make a plot of the total energy as it changes from steps 0, 1, 2, ..., 20.
• Modify cvt_1d by storing the location of all the generators at steps 0 through 20. This means we have an array G_SAVE of dimension 21xN, where row I is a copy of the values of G on step I. Have cvt_1d() return this array, and then try to plot the individual paths of each generator, all together on one plot. If you do this correctly, you will see 10 lines that start out at random points and wiggle upwards toward their proper locations, suggesting how the points G move towards their proper locations in the CVT.