【matlab】震撼！只能用这两个字了，真是奇女子

一千个圆零一个姑娘http://www.withparadox2.com/archives/90#comments

我参考Hanson在1960年发表的一篇论文（点击这里查看）做了这个动画。如果你看过我之前的一篇名为《图片迷宫》的文章，便不会对下面这位女同学感到陌生。与上次不一样，这一次我先征得了她的同意方才对照片进行了大刀阔斧的处理，由于手艺不好，勾勒出的人物效果也不尽人意，在此对她表达深深的歉意。

这个动画实际是一套行星圆系统，设想太阳不动，地球绕着太阳转，月球绕着地球转，XX绕着月球转……圆越多，得出的形状也就越复杂。Hanson在论文中论证了，如果给出的圆足够多，理论上是可以得出任何形状的轨迹线。下面简要介绍一下基本的制作步骤。

1，处理图片
首先需要勾勒出线条图，最好能够一笔从头画到尾。但在这个图里，由于人物的眉毛、嘴巴的存在，没办法一笔画出来，我采取的办法是加上辅助线，最后作图时将辅助线抹掉。

然后通过下面这段代码将上面的位图转化为坐标，但这里线段的顺序是乱的，需要很大的耐心去手动调整。目前没找到好的方法。

pic = Import["location of your picture"];
lcp = ListCurvePathPlot@
   Position[ImageData[Thinning@Binarize[ColorNegate@pic]], 1];
lines = Cases[lcp[[1, 1, 3]], _Line] /. {x_Real, y_Real} :> {y, -x};

2，复数及傅里叶变换
下面仅简单介绍一下与本文有关的知识，更具体的内容参看其他资料。
对于平面上的一点，可以用复数 z=x+iy 来表示，极坐标形式为 z=reiθ 。对于一个以坐标 c=x+iy 为圆心，以 ρ 为半径，以 T 为周期，以 α 为初始角的圆周运动，可用下式表示（其中t代表时间）：

z = c + ρ e i (2 π t / T + α)

例如当

t=t0 时，点的坐标为：

z = c + ρ e i (2 π t 0 / T + α) = x + i y + ρ \cdot c o s (2 π t / T + α) + i ρ \cdot s i n (2 π t / T + α) = (x + ρ \cdot c o s (2 π t / T + α) + i (y + ρ \cdot s i n (2 π t / T + α)))

若

f(t) 表示运动的点的轨迹，那么以此点为圆心做圆周运动的点可以表示为：

z = f (t) + ρ 1 e i (2 π t / T 1 + α 1)

对于N个圆，除了第一个圆心是固定点（

c=x+iy ）外，还剩下N个运动的点（除了N-1个运动的圆心，第N个圆上还有一个运动的点），由上可得出第N个圆上点的轨迹：

z = c + \sum k = 1 N - 1 ρ k e i (2 π t / T k + α k)

事实上，这已经和复数的离散傅里叶变换有了绝大部分的相似，下面就越来越简单了。首先，将之前得到的图片坐标转化为复数坐标，即：

(x,y)→x+iy ，然后对这N个复数进行傅里叶变换，得到

X[k] ，然后再用逆变换通过另一种形式得到原来的坐标：

x [n] = \sum k = 0 N - 1 X [k] e j 2 π k n / N

令

n/N=t ，

k=1/Tk 。上面

ρk 为实数，这里

X[k] 为复数，好像不一样。那是因为

X[k] 中既包括了半径，也包括了初始角，进一步变换可知：

x [t] = \sum k = 0 N - 1 X [k] e j 2 π t / T k = X [0] + \sum k = 1 N - 1 | X [k] | X [ k ] | X [ k ] | e j 2 π t / T k = X [0] + \sum k = 1 N - 1 | X [k] | e j α k e j 2 π t / T k = X [0] + \sum k = 1 N - 1 | X [k] | e j (2 π t / T k + α k) = X [0] + \sum k = 1 N - 1 | X [k] | e j (2 π t / T k + α k)

这里的

|X[k]| 便是圆的半径，通过求和得到圆心坐标，当改变

t(0⩽t⩽1) ，圆心坐标就会得到更新，从而生成连续动画。

在作图时，可以根据半径大小排序，将半径小的作为“笔尖”，并且可以去掉一些半径太小的圆。

3，一千个圆零一个姑娘【看原址，辣鸡CSDN传不上】