注！这是一篇水文，只有仿真，不讲正事儿~

这个是最近网上挺好玩的一个东西，来源是这篇论文，说白了就是一个很简单的操作，现在空间中随机生成N个点，然后把它们首尾相连，变成一个封闭的多边形，这样就得到N条线段，然后再把N条线段的中点依次相连，得到新的封闭多边形，不断重复上述过程，最后一定会变成一个椭圆；

如果是三个点，那么可以想象最后会变成一个超小的三角，无限接近于一个点，但是我们要知道，一个无限接近于点的点是特殊的圆，圆，则是特殊的椭圆。。。

至于理解嘛，网上已经有大神给出了很好的证明了，我这就不凑热闹了，那篇文章如开篇所言，其实这篇证明很简单；

但是考虑到有想自己研究证明但是还证不出来的高志向的同学，或者完全连一点数学都不想看，但是想了解一下思想以便将来好装逼的少年们，我这里就不涉及任何数学地简单说一下；

其实要想理解，只要把取中点这个操作看成一个矩阵的变换就可以了，被变换的矩阵是一个N×2的矩阵，然后乘以这个取中点的变换矩阵，就得到变换后的结果，那么其实最后的结果就是这个变换矩阵的k次幂再乘以一开始的N×2的矩阵而已，k是迭代次数。

那么我们所需要理解的就是这个变换矩阵的k次幂而已，然后自己写一下，就会发现这个变换矩阵很简单，是一个循环矩阵，啥叫循环矩阵，就是把一个行向量循环右移一位，作为第二行，再循环右移一位，作为第三行，一直做下去，最后就得到一个循环矩阵了。。。

然后假设循环矩阵，也就是变换矩阵为M，然后从线性代数，我们知道一个矩阵如果具有某些特性，就可以写成PΛP^-1的形式，其中Λ是M的特征值构成的对角矩阵；

关于那么M^k就是PΛ^kP^-1，然后就引进一个性质，循环矩阵的性质，很有趣的，大家可以自己验证一下，循环矩阵的特征值就是第一行那个东西的傅里叶变换！！【我是第一次听说，mathematica验证了一下，是真的。。废话。。】

然后，唉，算了，接下来必须涉及数学了，大家还是看原文去吧。。反正接下来就是把傅里叶变换矩阵代进去，发现k趋于无穷大的时候，不为零的就只剩下一项了，然后那一项算一下发现满足椭圆的性质就是了。。

其实我是来贴代码的，不小心又说多了。。

算了~反正我就无聊用mathematica和matlab各写了一个，其实我就是想玩一玩gif生成这个功能而已。。。【然后实验室的师弟用MFC也写了一个。。。】

下面是mathematica的：

mathematica code

Draw[pl_] := Module[{copy = pl}, 
    AppendTo[copy, pl[[1]]]; ListLinePlot[
      copy, Epilog -> {PointSize[Large], 
        Red, Point[copy]}]]; 
GetCenterList[pl_] := Module[{result}, 
    result = (pl + RotateLeft)/2]; 
iter = 400;
NN = 30;
PointList = RandomReal[10, {NN, 2}];
gifl = NestList[GetCenterList, PointList, iter];
plotlist = Table[Draw[gifl[[i]]], {i, 1, iter}];
ListAnimate[plotlist]

list = RandomReal[{-1, 1}, {100, 2}];
Dynamic[Graphics@Polygon[list = (list + RotateLeft@list)/2]]

下面是matlab版本，嗯，比较之下，我比较喜欢mathematica版本的~

function GetElic()
N = 30;
PointList = rand(N,2);
for i = 1 : 250
    plot_wrpa(PointList);
    PointList = GetCenterList(PointList);
    PointList = ZoonOut(PointList);
    drawnow
end

function plot_wrpa(PointList)
PointList = [PointList ; PointList(1,:)];
plot(PointList(:,1),PointList(:,2));
hold on;
plot(PointList(:,1),PointList(:,2),'r.','MarkerSize',20);
hold off;
axis([0,1,0,1]);

function CenterList = GetCenterList(PointList)
CenterList = PointList;
PointList = [PointList ; PointList(1,:)];
for i = 1 : size(CenterList,1)
    CenterList(i,:) = (PointList(i,:)+PointList(i+1,:))/2;
end

function PointList = ZoonOut(PointList)
minx = min(PointList(:,1));
maxx = max(PointList(:,1));
PointList(:,1) = (PointList(:,1)-minx)/(maxx-minx);
miny = min(PointList(:,2));
maxy = max(PointList(:,2));
PointList(:,2) = (PointList(:,2)-miny)/(maxy-miny);

【完】

本文内容遵从CC版权协议，转载请注明出自http://www.kylen314.com