首页 > 数学 > 任意点变换成椭圆

任意点变换成椭圆

2013年10月11日 发表评论 阅读评论

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

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

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

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

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

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

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

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

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

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

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

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

下面是mathematica的:

mathematica code


另外,我看到网上也有流传这个版本:

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);

b


【完】

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

分类: 数学 标签: ,
  • 博主好厉害~问下那段MMAcode~我粘到我的MMA9.0上运行的时候,卡在最后一句跑了五分钟没跑出来还崩溃了大概是哪儿有问题啊- -

  • Hi~ o(* ̄▽ ̄*)ブ

    start(timer('timer',{@(~,~,p)set(p,...'verti',(p.Vertices+circshift(p.Vertices,-1))/2),... fill(rand(100,1),rand(100,1),'r')},'P',.01,'Ex','FixedR'))