用圆搞定FBB---FBB、Matlab与航天

公众号：理念世界的影子

文不可无观点，观点不可无论据。

转载请注明出处

前两天和朋友聊天，说前几期MATLAB公众号文章，是想了很久的，理解了对怎么写更漂亮的MATLAB程序大有帮助，保证任何书上都不会有，得意之作迫切希望和大家分享、交流、探讨，但阅读量太少更无人探讨。朋友说，你要在Matlab中写点FBB，阅读量就可以上去了。MATLAB和FBB有毛关系，这不难为人嘛，绞尽脑汁，终于编了一个“用圆搞定FBB”的程序，具体如上图。

上图由200个圆(个数可程序自设)，1个圆接1个圆，所有圆都绕前一个圆作匀速转动，最后一个圆的一条固定半径最终勾勒除了F的轨迹。200个圈子真乱，运动起来令人眼花缭乱。圈子少一点行不行？还真不行，如下图，用了9个圈子。

这样写出来的'F'软趴趴的，最上面还勒脖子了，必须更多的圈子才能更直。真是左旋右转不知疲，千匝万周无已时。人间物类无可比，奔车轮缓旋风迟。

这么乱的圈子怎么获得呢？傅立叶级数(助记法：服你爷)，不解释，看文后代码。

完美地画出了'F'，怎么画出'B'呢？将代码第一句中的'F'改为‘B’，发现如下图，绘制的曲线缺少中间的两个洞。

这是因为自动识别轮廓的原因，程序将内轮廓放到元胞数组了，此时需要将程序的第21行改为 edge_pts=[B{1}; B{2}; B{3}]; 即可补充上两个洞，

但由于存在洞，结果是不完美的，如下图，绘制出的图形会有一些跨越和毛刺，而且在外轮廓来回打转。有洞后，边界坐标的自动识别有困难，要画出更为完美的'B'，还得靠网友人工辅助来识别坐标。

用圆搞定FBB演示结束，这和航天（或宇宙）有什么关系呢？有点关系。

古希腊人认为，理念世界是完美的，毕达哥拉斯设想所有的星星固定在天球上，天球是一个标准的圆球，天球带着星星一起转。

但行星的运行极为诡异，不但时快时慢，而且有时会逆行，孤立的圆球模型无法表述。阿波罗尼乌斯提出了“本轮-均轮”模型，托勒密将之发展到极为精巧的程度。如下图，大圆叫均轮，其它所有小圆叫本轮，一个本轮即可解释火星轨道逆行。如果一个本轮形成的轨道不够精确怎么办？增加本轮数量即可。其实这就是傅立叶级数的雏形。

上面用圆搞定FBB的演示也可视为均轮+199个本轮的体系。有点小差别是这里需要让均轮左右往复运动。由于此处仅对Y坐标进行了傅立叶级数展开，展开后就已确定了划过的X轨迹，这个轨迹与原始轨迹不一致，需要人为对齐。是否存在均轮不运动的处理方法和解，笔者仍在思考中，期待大家一起讨论。

扯淡完毕。本篇偷懒了，“本轮-均轮”体系的历史，傅立叶级数，代码的思路都懒得详细写了，后续有闲补上。此处直接上代码。

fbb.m

---

1. str4draw='F'; % 需要写的字符

2. ncircle=200; % 使用圆的数量

3. giffilename='F.gif'; % 输出的gif图像名称

4. clc; close all;

5. %% 直接利用matlab获取字符轮廓坐标，当然，如果有坐标，无需此部分代码

6. % 思路: 1. 用text写出字符并放在坐标系中间

7. %      2. 截图，并提取图像边缘，获得曲线坐标

8. axis off; % 避免坐标线的干扰 

9. h=text(0, 0, str4draw, 'fontsize', 360, 'fontname', 'times new roman', 'VerticalAlignment', 'bottom');

10. range=get(h,'extent'); % 获取文字边缘在坐标系中的位置，但文字大小与坐标系无关，只能将窗口拉大

11. limx=get(gca, 'xlim'); limy=get(gca, 'ylim'); % 坐标系

12. scale=range(3:4)./[diff(limx) diff(limy)]; % 从文字边缘坐标和当前坐标，获取窗口需要拉大的系数，不太准确，但够用了

13. pos=get(gcf, 'pos'); % 获取窗口位置

14. set(gcf, 'pos', [pos(1:2) pos(3:4).*scale]); % 将窗口拉大

15. img=getframe; img=img.cdata; close; % 获得字符截屏

16. bm=edge(rgb2gray(img(:, :, :))); % 获取字符边缘

17. B=bwboundaries(bm); % 将边缘变成曲线

18. edge_pts=B{1}; edge_pts=edge_pts*[0 -1; 1 0]; % 记录曲线坐标，获得的坐标是躺着着，因此顺时针转90°

19. %% 通过傅立叶变换，获取圆的半径、旋转速度和初始相角

20. % 思路: 将y视为函数，通过傅立叶变换获得圆的半径、旋转速度和初始相角

21. len_pts=length(edge_pts); % 数组大小

22. dt=0.1; t_pts=linspace(0, len_pts*dt, len_pts); % t_pts为时间

23. Y = fft(edge_pts(:,2))/len_pts; Y=2*Y(1:len_pts/2+1); % 傅里叶变换获取幅值

24. f=1/dt/2*linspace(0,1,len_pts/2+1)'; % 对应频率

25. [ans index]=max(Y); if(index==1) y0=Y(1)/2; Y(1)=[]; f(1)=[]; end % 如果a0不等于0，将之单独列出

26. [ans index]=sort(Y, 1, 'descend'); index=index(1:ncircle)'; % 将所有圆按直径大小排序

27. r_eachcircle=abs(Y(index)); % 圆系列尺寸

28. speed_eachcircle=f(index)*360; % 圆系列旋转速度

29. initialtheta_eachcircle=angle(Y(index))*180/pi+90; % 圆系列初始相角

30. % [r_eachcircle speed_eachcircle initialtheta_eachcircle]

31. buf=[0:2:360]'; unit_circle=[[cosd(buf) sind(buf)]; 0 0; 1 0]; % 一个标准圆即一根半径

32. t_rotate=linspace(0, 360/(min(speed_eachcircle)), 100); % 旋转时用的时间

33. color={'b', 'c', 'm'}; color=repmat(color, 1, length(r_eachcircle)); % 每个圆使用不同颜色

34. %% 绘制各圆旋转图和划过的轨迹图

35. % set(gcf,'outerposition',get(0,'screensize')); % 全屏

36. set(gcf, 'unit', 'centi', 'pos', [5 5 20 20]); % 输出为适合贴公众号的尺寸

37. gif=[]; xy_writestr=[]; allcircles=[];

38.     axis equal; axis off; set(gcf, 'color', 'w'); % 显示图形等比例，去掉坐标轴等

39.     limx=[min(edge_pts(:,1)) max(edge_pts(:,1))]; limy=[min(edge_pts(:,2)) max(edge_pts(:,2))];

40.     limx=3.2*(limx-mean(limx))+mean(limx); limy=1.5*(limy-mean(limy))+mean(limy);

41.     xlim(limx); ylim(limy); % 设置图形范围，根据不同字符可能需要调整，如果先记录所有时刻的allcircles，则可直接从allcircles中提取范围

42. for t=t_rotate

43.     theta=speed_eachcircle*t+initialtheta_eachcircle; % 旋转角度

44.     xy0=[0 y0]; % 第一个圆圆心位置

45.     %% 计算所有圆的新位置

46.     for i=1:length(r_eachcircle)

47.         % 下句为旋转矩阵，以前一个圆旋转后位置为中心，进行下一个圆的旋转

48.         allcircles(i, :, :)=repmat(xy0, length(unit_circle), 1)+unit_circle*r_eachcircle(i)*[cosd(theta(i)) sind(theta(i)); -sind(theta(i)) cosd(theta(i))];

49.         xy0=allcircles(i, end, :); xy0=xy0(:)';

50.     end

51.    %% 考虑到y坐标满足了，x坐标被带走，因此用原来坐标重新标定，表现的效果就是第一个圆的平动

52.     xy0=[interp1(t_pts, edge_pts(:,1)', t)-xy0(1) y0];

53.     %% 下面又计算了一次，也许会有更好的办法

54.     for i=1:length(r_eachcircle)

55.         allcircles(i, :, :)=repmat(xy0, length(unit_circle), 1)+unit_circle*r_eachcircle(i)*[cosd(theta(i)) sind(theta(i)); -sind(theta(i)) cosd(theta(i))];

56.         xy0=allcircles(i, end, :); xy0=xy0(:)';

57.     end

58.     % 最末点

59.     xy_writestr=[xy_writestr; xy0]; % 最后一个圆划过的轨迹

60.     cla; hold on;

61.     plot(limx, -r_eachcircle(1)*[1 1]+y0, 'b', edge_pts(:,1), edge_pts(:,2), 'y--'); % 绘制轨道和原曲线

62.     for i=1:length(r_eachcircle)

63.         plot(allcircles(i, :,1), allcircles(i, :,2), color{i}, 'linewidth', 1); % 绘制各个圆

64.     end

65.     plot(xy0(:,1), xy0(:, 2), 'o', xy_writestr(:,1), xy_writestr(:, 2), 'r', 'linewidth', 4); % 绘制最后一个圆划出的曲线

66.     hold off;

67.     buf=getframe; gif=[gif buf]; % 录制屏幕

68.     %% 输出gif

69.     I=frame2im(buf);

70.     [I map]=rgb2ind(I, 256);

71.     if(t==0)

72.         imwrite(I, map, giffilename, 'gif', 'WriteMode', 'overwrite', 'DelayTime', 0.2, 'Loopcount', inf);

73.     else

74.         imwrite(I, map, giffilename, 'gif', 'WriteMode', 'append', 'DelayTime', 0.2);

75.     end

76. end

77. movie(gif, 100, 100); % 播放动画

---

往期文章：

MATLAB

微信扫一扫

关注“理念世界的影子”

版权声明：本文是"洞穴之外"作者原创文章，欢迎转载，须署名并注明来自“理念世界的影子”公众号。