一种基于Voronoi图的曲边化随机颗粒模型生成方法

简介 

 在数值模拟中，除了对物体的整体模拟，还要对物体的细观层面进行模拟，以探讨裂纹的生产、扩展；或者细观层面的变化情况。特别是对于某些颗粒增强型复合材料。

楼主查询到目前有关颗粒增强型复合材料的细观模型主要有两类：

第一类是基于蒙特卡洛法随机生成圆形颗粒模型。这种方法生成的细观模型圆形颗粒的粒径可控，颗粒的大小也可控。但是由于生成的颗粒粒径太过规范，与实际不相符。对于某些大颗粒的材料往往外形表面不是如此规整的圆。第二类方法基于Voronoi图，构建多边形颗粒的模型，颗粒的形状具有随机性，颗粒的分布也具有随机性，但是颗粒的边缘过于尖锐，这与颗粒的形状也不相符。

（注：电镜图来源于网络，侵删）

一、方法简介

在泰森多边形的基础上，对多边形进行曲边化处理，使其变为一个曲边化图形。结果如下图所示：

 

使用工具：Matlab；Sketchup（包含rubby语言编译器，圆锥曲线插件） 

二、算法实现

2.1使用matlab生成泰森多边形并获取多边形的每个顶点信息（matlab源码如下：具体不在多述，详情见连接

https://blog.csdn.net/weixin_42943114/article/details/82319332?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2~all~first_rank_v2~rank_v25-2-82319332.nonecase&utm_term=%E7%94%A8matlab%E5%AE%9E%E7%8E%B0%E6%B3%B0%E6%A3%AE%E5%A4%9A%E8%BE%B9%E5%BD%A2&spm=1000.2123.3001.4430）

①matlab生成泰森多边形源码 

%实际使用矩阵 x = rand(20,1); y = rand(20,1); xy=[x,y]; %测试所用矩阵 % x=gallery('uniformdata',[5 1],3); % y=gallery('uniformdata',[5 1],2); % xy=[x,y]; figure(1) TRI = delaunay(x,y); triplot(TRI,x,y,'k'); axis([0 1 0 1]); title('Delaulay三角形'); hold on figure(2) voronoi(x,y); axis([0 1 0 1]); title('泰森多边形')

②需要得到泰森多边形的顶点坐标信息和每个点的顺序，代码如下，需要提取泰森多边形的顶点坐标[vx,vy]；以及连接顺序[V,R],详细介绍见https://ww2.mathworks.cn/help/matlab/ref/delaunaytriangulation.voronoidiagram.html。在matlab的工作区获取[vx,vy]，以及V，R数组。

[vx,vy] = voronoi(x,y); %vx和vy是我们所需要的voronoi图所有点（即包含delaunay三角形外心，也包含inf点） dt = delaunayTriangulation(x,y); CC = circumcenter(dt); [V,R] = voronoiDiagram(dt); % R1=R(10,1); % R2=R{10,1}; % R3=R{10,1}(1,2); %V是所有delaunay三角形的外心，但不包含inf点（区别vx和vy） %R代表点的组成方式，是一个元胞数组；元胞数组运算如下： %①R{m,n}返回m行，n列的具体内容，如果内部是一个数组，那R{m,n}的类型也为数组， %此时R{m，n}符合数组的类型规范，如R{m,n}(i,j)，调用胞元数组m行，n列的数组的i行j列 %②R(m,n)返回m行，n列元胞数组内容，返回值仍然是元胞数组，而不是具体内容。

③为了直接生成Ruby原码，在matlab对数据[V,R]中R进行处理，代码如下： 

sz1=size(R,1); % size（数组，1）取得某个数组的行数 % size（数组，2）取得某个数组的列数 % size（数组）取得某个数组的行数和列数，返回值为一个[m,n]一维矩阵 %以下为sketchup 中RUBY语句 for r1=1:1:sz1 sz2=size(R{r1},2); if R{r1}(1)==1 fprintf('l%d= entities.add_line(',r1) else fprintf('F%d= entities.add_face(',r1) end %以上if语句判断是要连成面还是连线 if R{r1}(1)==1 k=2; else k=1; end %以上if语句判断是消除胞元数组R中第一组数据inf数据 for r2=k:1:sz2 if r2<sz2 fprintf('p%d,',R{r1}(r2)) else fprintf('p%d',R{r1}(r2)) end %以上if语句为了消除（p1,p2,p3,）这样的形式出现保证（p1,p2,p3）正确格式 end fprintf(')\n') %最后ruby所得的语句直接从命令窗口复制过去粘贴，再结合excel所处理的就可以了 end

 上述代码一起复制到matlab中，可直接使用。

2.2将Voronoi图绘制在sketchup中

方法一、由于sketchup无法直接根据点的坐标直接输入点，只能通过的程序控制台Ruby控制台输入参考点，之前获取的[vx,vy]坐标信息排上用场了，如下图，X,Y分别代表vx坐标（横坐标）和vy坐标（纵坐标）,按顺序（先复制粘贴ruby指令1，在复制粘贴ruby指令2）复制Ruby指令到ruby编译器，运行后就可以在sketchup中创建泰森多边形的所有顶点。

当在sketchup中得到所有的泰森多边形顶点后 ，使用自带的连线工具将所有的点按照matlab生成的voronoi图的形状连接起来，然后再创建一个1×1的正方形区域，将不必要的部分剔除，就可以得到一个和matlab一样的泰森多边形。

方法二、

 方法一生成的泰森多边形在后续的处理中可能会造成无法圆锥化问题，原因是有相同的点重叠问题。方法二的思路是直接通过ruby程序指令生成泰森多边形，然后再进行二次加工

①将matlab所得的V数组变量（不是[vx,vy]。该变量也为坐标变量，和[vx,vy]的区别后续再提及）复制到excel中的X Y中，自动生成Ruby指令，然后将excel中的rubby指令1从p2开始复制到rubby编译器中（p2~pn）。（注意使用V数组变量中第一组数据为inf，把inf也复制到excel中）

 ②在2.1章节③部分中，利用matlab直接生成rubby源码并输出在命令行窗口。如下

 

Ruby代码如下，表示生成相应的面或者线 （这些点的连接顺序会随着matlab生成的图形而随机变化）

F1= entities.add_face(p19,p25,p20,p22,p32) l2= entities.add_line(p31,p33,p29,p27) F3= entities.add_face(p22,p24,p33,p31,p32) l4= entities.add_line(p10,p4,p14) F5= entities.add_face(p2,p12,p23,p18,p13,p3,p8) l6= entities.add_line(p27,p7,p11,p10) F7= entities.add_face(p12,p28,p30,p24,p22,p20,p23) F8= entities.add_face(p7,p27,p29,p30,p28) F9= entities.add_face(p24,p30,p29,p33) F10= entities.add_face(p5,p16,p17,p26,p15) F11= entities.add_face(p2,p8,p9,p4,p10,p11) l12= entities.add_line(p14,p6,p3,p13,p5,p15) l13= entities.add_line(p26,p17,p19,p32,p31) F14= entities.add_face(p4,p9,p6,p14) F15= entities.add_face(p5,p13,p18,p21,p16) F16= entities.add_face(p16,p21,p25,p19,p17) F17= entities.add_face(p3,p6,p9,p8) F18= entities.add_face(p18,p23,p20,p25,p21) F19= entities.add_face(p2,p11,p7,p28,p12) l20= entities.add_line(p15,p26)

③即使通过以上步骤能够直接生成泰森多边形，但是还是有一部分的图形无法完整地描述 。

这里要强调一下两类坐标顶点V和[vx,vy]的区别。

根据泰森多边形建立的原则，每个多边形的顶点都是delaunay三角形外接圆圆心，由这些圆心相连即可得到相应的泰森多边形，但是结果如下：这种图形看着很怪，而且不在一个正方形封闭区域内。怎么办？matlab自动生成相关点构建泰森多边形的完整性，这些多余的点会自动生成。

改图为使用Ruby代码直接生成的图形

V代表的坐标数组就是所有denaunay三角形的外接圆圆心，而那些自动生成的点统统用inf表示。[vx,vy]代表的坐标则是所有的有关点的坐标，包含matlab自动生成的点。 

因此我们还要进一步创建V中的inf点，还需要进一步使用Rubby代码生成。

继续按照方法一将[vx,vy]坐标绘制在sketchup中，并将其保存在另一个图层中（方便修改），然后将与V数组重复的点删除，按照matlab生成的voronoi图连接。然后加上方框。如下：

这样就把泰森图绘制出来了

2.3基于sketchup圆锥曲线插件，将泰森图曲边化

 首先声明，不是所有的图形都能够曲边化，有些图形需要进行一定修改，确保不会出现过大或者过小的角，若果过大或者过小都无法自动生成响应的曲线。

首先安装圆锥曲线插件，然后按下Tab键确定相关参数，最后选择需要曲边化的图形区域点击鼠标左键即可。结果如下：

 然后选择模型→右键选择炸开模型→删除填充部分→得到颗粒填充材料细观模型，如下：