利用函数来画任意图形

dft

最近我们经常看到像鸟叔,初音之类的通过函数图像来画出来,看上去十分神奇的样子,wolframalpha这里有大量的通过函数图像来画人物的例子,大家可以去围观,而且最上面我这几个字也是我用函数图像画出来的,今天我们就说说这是怎么做到的。

首先我画的图形的函数是这个样子的

x(t)= 3.69696969697 *cos( 0.0 *t)- 1.78787878788 *sin( 0.0 *t) + -0.608631557183 *cos( 0.190399554763 *t)- 0.637886769069 *sin( 0.190399554763 *t) + -1.00348014017 *cos( 0.380799109526 *t)- -0.251077539499 *sin( 0.380799109526 *t) + 0.0858300019403 *cos( 0.571198664289 *t)- 0.0312569171374 *sin( 0.571198664289 *t) + 0.0403938878342 *cos( 0.761598219052 *t)- -0.424077642684 *sin( 0.761598219052 *t) + -0.029533908125 *cos( 0.951997773815 *t)- -0.612124881382 *sin( 0.951997773815 *t) + -0.196190215552 *cos( 1.14239732858 *t)- 0.0779263864101 *sin( 1.14239732858 *t) + -0.317379458544 *cos( 1.33279688334 *t)- -0.107897599954 *sin( 1.33279688334 *t) + 0.195824060124 *cos( 1.5231964381 *t)- -0.0603270845874 *sin( 1.5231964381 *t) + 0.0206709920939 *cos( 1.71359599287 *t)- -0.192643145398 *sin( 1.71359599287 *t) + -0.0390613766898 *cos( 1.90399554763 *t)- -0.0233815418361 *sin( 1.90399554763 *t) + 0.120124291368 *cos( 2.09439510239 *t)- 0.165578836822 *sin( 2.09439510239 *t) + 0.118680071019 *cos( 2.28479465716 *t)- 0.192116690811 *sin( 2.28479465716 *t) + 0.000254592255497 *cos( 2.47519421192 *t)- 0.0310134506924 *sin( 2.47519421192 *t) + 0.082560849545 *cos( 2.66559376668 *t)- -0.0576194138539 *sin( 2.66559376668 *t) + -0.0926779308527 *cos( 2.85599332145 *t)- -0.0839813133077 *sin( 2.85599332145 *t) + -0.192145472441 *cos( 3.04639287621 *t)- -0.00294302065064 *sin( 3.04639287621 *t) + -0.0467279999084 *cos( 3.23679243097 *t)- 0.0312620057434 *sin( 3.23679243097 *t) + -0.131380808582 *cos( 3.42719198573 *t)- 0.0240572405524 *sin( 3.42719198573 *t) + -0.146394824366 *cos( 3.6175915405 *t)- 0.0428285063991 *sin( 3.6175915405 *t) + -0.0688190483421 *cos( 3.80799109526 *t)- -0.02541184382 *sin( 3.80799109526 *t) + -0.0901484634214 *cos( 3.99839065002 *t)- 0.181258333651 *sin( 3.99839065002 *t) + -0.0898212610648 *cos( 4.18879020479 *t)- -0.0443667156102 *sin( 4.18879020479 *t) + -0.0750426044852 *cos( 4.37918975955 *t)- 0.0736544316679 *sin( 4.37918975955 *t) + -0.119525232263 *cos( 4.56958931431 *t)- -0.0355835823163 *sin( 4.56958931431 *t) + -0.193150001399 *cos( 4.75998886908 *t)- 0.0319541260203 *sin( 4.75998886908 *t) + 0.25777892582 *cos( 4.95038842384 *t)- -0.000416186591487 *sin( 4.95038842384 *t) + 0.171027143418 *cos( 5.1407879786 *t)- 0.218778298585 *sin( 5.1407879786 *t) + -0.16745165654 *cos( 5.33118753336 *t)- -0.127545487283 *sin( 5.33118753336 *t) + 0.242332056731 *cos( 5.52158708813 *t)- -0.263812652924 *sin( 5.52158708813 *t) + 0.0367447452299 *cos( 5.71198664289 *t)- -0.534397947337 *sin( 5.71198664289 *t) + -0.325288006456 *cos( 5.90238619765 *t)- 0.163558793585 *sin( 5.90238619765 *t) + -1.13634134796 *cos( 6.09278575242 *t)- -0.84340197599 *sin( 6.09278575242 *t)
y(t)= 1.78787878788 *cos( 0.0 *t)+ 3.69696969697 *sin( 0.0 *t) + 0.637886769069 *cos( 0.190399554763 *t)+ -0.608631557183 *sin( 0.190399554763 *t) + -0.251077539499 *cos( 0.380799109526 *t)+ -1.00348014017 *sin( 0.380799109526 *t) + 0.0312569171374 *cos( 0.571198664289 *t)+ 0.0858300019403 *sin( 0.571198664289 *t) + -0.424077642684 *cos( 0.761598219052 *t)+ 0.0403938878342 *sin( 0.761598219052 *t) + -0.612124881382 *cos( 0.951997773815 *t)+ -0.029533908125 *sin( 0.951997773815 *t) + 0.0779263864101 *cos( 1.14239732858 *t)+ -0.196190215552 *sin( 1.14239732858 *t) + -0.107897599954 *cos( 1.33279688334 *t)+ -0.317379458544 *sin( 1.33279688334 *t) + -0.0603270845874 *cos( 1.5231964381 *t)+ 0.195824060124 *sin( 1.5231964381 *t) + -0.192643145398 *cos( 1.71359599287 *t)+ 0.0206709920939 *sin( 1.71359599287 *t) + -0.0233815418361 *cos( 1.90399554763 *t)+ -0.0390613766898 *sin( 1.90399554763 *t) + 0.165578836822 *cos( 2.09439510239 *t)+ 0.120124291368 *sin( 2.09439510239 *t) + 0.192116690811 *cos( 2.28479465716 *t)+ 0.118680071019 *sin( 2.28479465716 *t) + 0.0310134506924 *cos( 2.47519421192 *t)+ 0.000254592255497 *sin( 2.47519421192 *t) + -0.0576194138539 *cos( 2.66559376668 *t)+ 0.082560849545 *sin( 2.66559376668 *t) + -0.0839813133077 *cos( 2.85599332145 *t)+ -0.0926779308527 *sin( 2.85599332145 *t) + -0.00294302065064 *cos( 3.04639287621 *t)+ -0.192145472441 *sin( 3.04639287621 *t) + 0.0312620057434 *cos( 3.23679243097 *t)+ -0.0467279999084 *sin( 3.23679243097 *t) + 0.0240572405524 *cos( 3.42719198573 *t)+ -0.131380808582 *sin( 3.42719198573 *t) + 0.0428285063991 *cos( 3.6175915405 *t)+ -0.146394824366 *sin( 3.6175915405 *t) + -0.02541184382 *cos( 3.80799109526 *t)+ -0.0688190483421 *sin( 3.80799109526 *t) + 0.181258333651 *cos( 3.99839065002 *t)+ -0.0901484634214 *sin( 3.99839065002 *t) + -0.0443667156102 *cos( 4.18879020479 *t)+ -0.0898212610648 *sin( 4.18879020479 *t) + 0.0736544316679 *cos( 4.37918975955 *t)+ -0.0750426044852 *sin( 4.37918975955 *t) + -0.0355835823163 *cos( 4.56958931431 *t)+ -0.119525232263 *sin( 4.56958931431 *t) + 0.0319541260203 *cos( 4.75998886908 *t)+ -0.193150001399 *sin( 4.75998886908 *t) + -0.000416186591487 *cos( 4.95038842384 *t)+ 0.25777892582 *sin( 4.95038842384 *t) + 0.218778298585 *cos( 5.1407879786 *t)+ 0.171027143418 *sin( 5.1407879786 *t) + -0.127545487283 *cos( 5.33118753336 *t)+ -0.16745165654 *sin( 5.33118753336 *t) + -0.263812652924 *cos( 5.52158708813 *t)+ 0.242332056731 *sin( 5.52158708813 *t) + -0.534397947337 *cos( 5.71198664289 *t)+ 0.0367447452299 *sin( 5.71198664289 *t) + 0.163558793585 *cos( 5.90238619765 *t)+ -0.325288006456 *sin( 5.90238619765 *t) + -0.84340197599 *cos( 6.09278575242 *t)+ -1.13634134796 *sin( 6.09278575242 *t)

如果像我们这里只用cos和sin的话,我们可以画出任意我们想画的闭合曲线,至于其他图像那样包含很多闭合曲线的是用了step function的技巧,这里我就暂时不说了,也就说我这里是说明如何画出任意的闭合曲线的。当然了也就是说只要你能一笔画的东西都可以,线当然是可以交叉或者重合的。

因为我们是要画闭合曲线,并且注意到我们的函数最后是x(t)和y(t)这种形式,所以也就是说,我们可以把x和y分别当作是周期函数来对待,说道周期函数当然就是想到傅里叶级数了,因为傅里叶级数可以逼近任意的周期函数。所以这里我们就要用到DFT离散傅立叶变换

我们这里就要用到离散傅里叶变换和逆变换的公式

\(\hat{x}[k]=\sum_{n=0}^{N-1}e^{-i\frac{2\pi}{N}nk}x[n]\)
\(x[n]=\frac{1}{N}\sum_{k=0}^{N-1}e^{-i\frac{2\pi}{N}nk}\hat{x}[k]\)

也就是说我们可以把我们想要画的图形画出来,然后依次找出我们要连接的这些点,这些点就是x[n],带入到上面的第一个公式,我们可以得到x[k],于是再调用第二个公式也就是离散傅里叶的逆变换,最后利用欧拉公式把e^ix=cosx+isinx展开,就能得到我们最后的函数形式了。

按照这个思路实现的python代码如下

  1. import math
  2.  
  3. N = int(raw_input())
  4.  
  5. f = []
  6. for i in range(N):
  7.     (x, y) = map(float, raw_input().split())
  8.     f.append(complex(x, y))
  9.  
  10. F = []
  11. for i in range(N):
  12.     ang = -2 * 1j * math.pi * i / N
  13.     r = 0
  14.     for j in range(N):
  15.         r += (math.e ** (ang * j)) * f[j]
  16.     F.append(r)
  17.  
  18. print "set parametric"
  19. print "set samples", N + 1
  20.  
  21. print "x(t)=",
  22. for i in range(N):
  23.     ang = 2 * math.pi * i / N
  24.     if i > 0:
  25.         print "+",
  26.     print F[i].real / N, "*cos(", ang, "*t)-",
  27.     print F[i].imag / N, "*sin(", ang, "*t)",
  28.  
  29. print
  30. print "y(t)=",
  31. for i in range(N):
  32.     ang = 2 * math.pi * i / N
  33.     if i > 0:
  34.         print "+",
  35.     print F[i].imag / N, "*cos(", ang, "*t)+",
  36.     print F[i].real / N, "*sin(", ang, "*t)",
  37.  
  38. print
  39. print "plot [t=0:", N, "] x(t), y(t)"
  40. print "pause 60"

也可以在这里看代码,把代码存成dft.py然后运行python dft.py输入点的个数以及点的位置,运行程序就可以看到生成的函数了,把结果输入到gnuplotli就可以看到函数图像了。

比如我生成的那个函数图像的输入是这个样子的

33
0 0
2 0
1 0
1 3
0 3
2 3
2 2
3.5 0
5 2
5 3
4 3
3.5 2
3 3
2 3
3 3
3.5 2
4 3
5 3
5 1
5.5 0
6.5 0
7 1
7 3
7 1
6.5 0
5.5 0
5 1
5 3
4 3
3.5 2
3 3
1 3
1 0

可以把输入文件保存到input.txt,然后运行cat input.txt | python dft.py | gnuplot就可以看到绘制好的函数了。

也就是说只要把你要画的东西的点描绘出来,输入到程序中就可以生成不可思议的函数了!

参考资料:
http://mathematica.stackexchange.com/questions/17704/how-to-create-new-person-curve
http://tieba.baidu.com/p/2156093774
http://www.quora.com/Mathematics/How-is-the-Gangnam-Style-mathematical-plot-made

您可能喜欢:
我猜您可能还喜欢:
,

《利用函数来画任意图形》有 11 条评论

  1. 第一次用 baidu in put for mac | #1

    lz,你的back to top是用哪个插件啊?

  2. jeff | #2

    我有个问题哈,傅里叶那个公式里面的X[n],在代码中怎么是有图像中的坐标(x,y)构成的x+iy构成的呢?不是说分开考虑X[n]和Y[n]吗?这里我没有搞懂,求解释一下

    • isnowfy

      @jeff: 确实,x和y是不相关的做两次虚部都是0的构造就可以,正是因为不相关,所以我可以把y放到虚部,这样只做一次,就能同时得到x和y的公式,而不用做两次

  3. DreamxForest | #3

    没看懂代码,但是X[n]求出来在用欧拉变换应该是cos+i*sin的,有实和虚,怎么弄

  4. DreamxForest | #4

    怎么我用你给的x和y画出来的是一堆圈,t我取了0到60

  5. ICTCLAS用过没,记得做毕设的使用用到过,就是不开源,当时是自己用java写的,最大正向匹配结合最大逆向匹配正确率会更高写

发表评论