函数共 268 页 11783 行，前后折腾了近半个月。可惜还是有些小瑕疵，不然文章标题就该叫《震惊！死宅竟用三角函数画出可爱的老婆》了......
前言
在推上看到这样一张图：
觉着用函数画妹子什么的好有趣，决定自己试着实现一下。
思路
从图中可以看到，绘制妹子的函数为三角函数项的和。不难发现，函数以 t 为参数，两个一组分别为正弦级数和余弦级数，首项(a0/b0)为常数，很容易看出来，函数应该和傅里叶级数有关系。
大体思路如下：
- 预处理图片。用 MATLAB 将想要用函数绘制的图片读入
- 提取轮廓
- 获取轮廓的坐标
- 处理坐标。获取到的坐标为 N 点长的离散序列，将其进行 N 点离散傅里叶级数展开（DFS）
- 验证。将输出结果的函数绘制成图像，观察效果
实现
先画个简单点的。桌面上有张可爱的 Pusheen 图片，就用它啦。
1. 预处理图片
将 Pusheen 图片读入MATLAB。
cd C:\Users\Desktop
pic = imread('C:\Users\Ascend\Desktop\pusheen.jpg');
转换为二值图像。
pic = im2bw(pic);
show(pic)
转换为二值图像的原因是便于提取轮廓。看下效果。
2. 提取轮廓
这里直接调用了 MATLAB 自带提取轮廓的算法。
pic = edge(pic,'sobel');
imshow(pic);
效果如下图，看起来不错，轮廓分明线条干净。
3. 获取轮廓坐标
获取轮廓坐标，并将能够围成封闭曲线的坐标分组，则每组为一个有限长的离散序列，再将每个序列进行 DFS 展开。
获取轮廓坐标并不难，难在拟合封闭曲线这里。想了想毫无头绪，上网搜索找到了 MATLAB 自带的 imcontour 函数。精度一般，使用条件苛刻，后文再叙，先凑合用。
point = imcontour(pic);
得到一个 2x126494 的矩阵 point，里面有我们需要的坐标。
4. 处理坐标
处理矩阵
坐标已经得到，先试着用坐标画一下图像。
width=point(2,:);
height=point(1,:);
plot(width,height)
却得到了一副非常奇怪的图像。
按道理，坐标描点得到的图像就应该为 Pusheen 猫的轮廓。注意到图像有几个异常点，横坐标特别大。查看矩阵 point 的结构。
第一行为纵坐标 height ，第二行为横坐标 width，每一列为一个点 (height,width)。
第一个点的数值很奇怪。横坐标为 2435 ，纵坐标为 0.1，怎么看都不是轮廓上的点。也许是某种标记？于是向后翻了2435个点，看到第 2437 列：
又是一个 (0.1,2427)。再向后翻 2427 个点：
第 4865 列，又是一个 (0.1,2531)。于是这个点的意义就清楚了：横坐标为 0.1 的点为特殊标记，其值为该组点的个数。依此可以将这个 2x126494 的矩阵拆成多个小矩阵，每个矩阵为一组可以拟合成封闭曲线的坐标(imcontour函数认为的)。
拆分矩阵
%判断point(1,:)第一行 0.1 的个数，从而判断length的长度
length=0;
for i=1:size(point,2)
if point(1,i)==0.1
length=length+1;
end
end
%开始分段，拆成a1，a2，a3，...，一共length个二维数组
for i=1:length
temp=point(2,1);
point(:,1)=[];
eval(sprintf('a%d=point(1:2,1:temp);', i)); %将第i段数据送入第i个数组
point(:,1:temp)=[];
end
%纵坐标取负数，将离散数据转为坐标
for i=1:length
eval(sprintf('a%d(2,:)=-a%d(2,:);',i,i));
end
离散傅里叶级数展开
离散傅里叶级数的公式如下，根据公式，写出 MATLAB 代码。
file = fopen('outer.txt', 'w'); %I/O方式输出结果，结果输出到桌面的outer.txt中
for count=1:length
eval(sprintf('a=a%d;',count));
N = size(a,2);
%a->ax(虚数组)
are=a(1,:);
aim=a(2,:);
ax=are+aim*j;
clear are aim;
%A:傅里叶正变换数组
A=[];
for i = 0:N-1
ang = -2 * 1j * pi * i / N;
r = 0;
for j=0:N-1
r = r + exp(ang * j) * ax(j+1);
end
A(i+1)=r;
end
clear ang i j r;
fprintf(file,['t=1:',num2str(N),';\n']);
fprintf(file,'x(t)=');
for i = 0:N-1
ang = 2 * pi * i / N;
fprintf(file,'+');
fprintf(file,[num2str(real(A(i+1) / N)),'*cos(', num2str(ang), '*t)-']);
fprintf(file,[num2str(imag(A(i+1) / N)),'*sin(', num2str(ang), '*t)']);
end
fprintf(file,';\n');
fprintf(file,'y(t)=');
for i = 0:N-1
ang = 2 * pi * i / N;
fprintf(file,[num2str(imag(A(i+1) / N)),'*cos(', num2str(ang), '*t)+']);
fprintf(file,[num2str(real(A(i+1) / N)),'*sin(', num2str(ang), '*t)+']);
end
fprintf(file,'0;\n');
fprintf(file,'plot(x,y);hold on;\n');
end
fclose(file);
5. 验证
输出的函数在桌面的 outer.txt 文件中。
将函数输入到 MATLAB 中画图像，结果如下。
验证无误。输出那一堆就是能画出可爱的 Pusheen 的函数啦 ~
其他
卡通小猫都画出来了，打算也画一个卡通妹子。在网上找了张艾米莉亚线稿
使用了用同样的方法，结果如下。
一脸茫然，大概轮廓倒能看出来，可怎么会乱成这样子[笑cry]
仔细想了想，原因可能出在闭合曲线分组那儿。卡通猫轮廓上的蜜汁线条也可能因为此。
对卡通猫的函数进行 Debug：
这些就差不多能看出原因了。imcontour 函数将鼻子、眼睛也认为是外轮廓的一部分，错误地将其分为一组数据。对于像 Pusheen 这种线条简单的图形来说，精度一般，分错组数较少。对于艾米莉亚这种复杂的人像来说，精度远远不够，错误就太多了。
回到从推上看到的图片。人家是如何画出那么完美的人像的呢？翻了翻博主的其他推文，发现人家利用游戏引擎 Hot Soup Processor 设计了一个画图板程序，在程序上手绘一组组封闭曲线，游戏引擎会用傅里叶级数将这组封闭曲线展开为函数。这个过程不需要用程序将离散序列分组，精度极高。
在没找到更好的分组方法前，就只能画画轮廓简单的图片啦[笑cry]
一些其他的函数曲线：
Twitter Curve：
1:
参考资料
https://mathematica.stackexchange.com/questions/17704/how-to-create-a-new-person-curve
http://blog.wolfram.com/2013/05/17/making-formulas-for-everything-from-pi-to-the-pink-panther-to-sir-isaac-newton/
http://blog.wolframalpha.com/2013/06/10/using-formulas-for-everything-from-a-complex-analysis-class-to-political-cartoons-to-music-album-covers/
http://blog.wolfram.com/2013/08/15/even-more-formulas-for-everything-from-filled-algebraic-curves-to-the-twitter-bird-the-american-flag-chocolate-easter-bunnies-and-the-superman-solid/
http://www.isnowfy.com/generate-any-function-curve/
https://yvt.jp/contours/
