三维图形的覆盖问题【附matlab代码】

三维图形的覆盖问题

引用本文代码发布网络分享文章，一定要提及本文。

---

前言

碎碎念：这里是记录一下我们小组做的一道校赛数学建模题，文章主要讲述的是我的部分的代码编写过程及思路，代码还没完全修改完，这里会放我当前做到的进度。后期有时间，做完了还会更新。

原题目：用5mm到10mm的球体填充下图中椭球体减去其与球体交集的部分后，剩余的三维图形:
图中下部蓝绿色球体区域是球心在坐标系原点，半径为4cm的球体。上半部分黄色区域，是一个x、y、z三个半轴分别为4cm、3cm、2cm，球心在（0，0，3）的椭球体。

---

提示：我参考的源代码链接也会附上

一、分析题目

在这道题之前还有一小问是二维的椭圆如何被小圆填充，因为网上直接能搜到代码，这里就不说了，直接说上一小问，我们组的搜集原理如何被用到这道题中。

正六边形可以密铺所有平面。圆致密的排列方式一定是要使每一个圆尽可能地靠近，并且是按照一定规律排列。阿克塞尔图早在1910 年就对圆形的六角密堆积进行了证明，对于圆形而言，六角密堆积是最致密结构。
由于问题中是用圆覆盖椭圆，允许圆有一定的重合。修改圆的半径，让原本相切的圆相交。

由图可知，三个圆交点引出的三条弦之间的夹角为120°，弦与半径之间的夹角为60°，由几何关系可得d=√3 r。

用上面的结论进行一层层的堆叠，要达到的效果是密堆积结构：

修改相交的半径，即可得到三维图形堆积所需要的球体排列。
暂时题目只考虑了椭球体内的小球填充

二、代码编写

1.题目图形的复现

clc,clear all
[xx,yy,zz]=sphere(40);
%椭球体半轴
x=xx*40;y=yy*30;z=zz*20;
%椭球体球心坐标
X=0;Y=0;Z=30;
%球体球心坐标
X1=0;Y1=0;Z1=0;
R=40;

%绘制椭球体
surf(x+X,y+Y,z+Z);
hold on 
axis equal
grid on
%绘制球体
%myplotsphere(X1,Y1,Z1,R);
%hold on 
%axis equal
%grid on

结果图：

2.小球排列方式

源代码链接：立方晶胞密堆积的matlab程序
源代码的思路是先在第一层通过错位的方式摆出四个小球，第二层同样用相同的方式摆出四个小球，然后让两层小球错位叠放。不过源代码是通过个数来控制小球的摆放，对后面我们用坐标来控制小球有点不方便，所以进行了修改。在修改的过程中，发现代码的一些数值设定也不太妥，通过计算又对代码进行了大改。最终得到现在的版本，只保留了源代码的基本框架。

首先，对代码摆放小球的球心重新计算。
利用正四面体的顶点来计算球心：

点A(0，0，0)；点B(0，a，0)；
点C(√3a/2，a/2，0)；点D(√3a/6，a/2，√6a/3)
点A和点B为两个球体之间的球心距离，由分析题目的结论，可得，a=√3r。

根据两层球心的摆放图，写出如下代码：

H1=sqrt(2)*R;
H2=sqrt(3)*R/2;
H3=3*R/2;
%第一层
balls(R,0+H2*2*i,    0+2*H3*j,  0+H1*2*k);
balls(R,-H2+H2*2*i,  H3+2*H3*j, 0+H1*2*k);
balls(R,H2+H2*2*i,   H3+2*H3*j, 0+H1*2*k);
balls(R,H2*2+H2*2*i, 0+2*H3*j,  0+H1*2*k);
%第二层
Hj=R/2;   %Hj=d
Hi=sqrt(3)*R/2;
balls(R,0+Hi+H2*2*i,    0+Hj+2*H3*j,  H1+H1*2*k);
balls(R,-H2+Hi+H2*2*i,  H3+Hj+2*H3*j, H1+H1*2*k);
balls(R,H2+Hi+H2*2*i,   H3+Hj+2*H3*j, H1+H1*2*k);
balls(R,H2*2+Hi+H2*2*i, 0+Hj+2*H3*j,  H1+H1*2*k);

这个代码也还是根据个数控制球的摆放

3.球心的限制条件

根据小球球心是否在椭球体内，得到填充球的限制条件。
提示：要考虑每一次绘制出的8个小球的球心，只要有其中一个在椭球体内，都要绘制小球，不过得分好几种情况绘制，后续还没做出来。

---

总结：现有全部代码

调用子函数1：

function balls(r,x,y,z)
[X,Y,Z]=sphere();
X=x+r*X;
Y=y+r*Y;
Z=z+r*Z;
surf(X,Y,Z);
hold on;
colormap([0 0 1]);

调用子函数2：

function  bodycube(i,j,k,R)
H1=sqrt(2)*R;
H2=sqrt(3)*R/2;
H3=3*R/2;
balls(R,0+i,    0+j,  0+30+k);
balls(R,-H2+i,  H3+j, 0+30+k);
balls(R,H2+i,   H3+j, 0+30+k);
balls(R,H2*2+i, 0+j,  0+30+k);

Hj=R/2;
Hi=sqrt(3)*R/2;
balls(R,0+Hi+i,    0+Hj+j,  H1+30+k);
balls(R,-H2+Hi+i,  H3+Hj+j, H1+30+k);
balls(R,H2+Hi+i,   H3+Hj+j, H1+30+k);
balls(R,H2*2+Hi+i, 0+Hj+j,  H1+30+k);

hold on
axis equal
axis on
colormap([.2 0.2 1]);

主函数：

clc,clear all
[xx,yy,zz]=sphere(40);
x=xx*40;
y=yy*30;
z=zz*20;
X=0;Y=0;Z=30;
X1=0;Y1=0;Z1=0;
R=40;

surf(x+X,y+Y,z+Z);
hold on 
axis equal
grid on
%绘制大球体
%myplotsphere(X1,Y1,Z1,R);
%hold on 
%axis equal
%grid on

R=10;
H1=sqrt(2)*R;
H2=sqrt(3)*R/2;
H3=3*R/2;
Hj=R/2;
Hi=sqrt(3)*R/2;
for y=-30-2*H3:2*H3:30+2*H3
    for x=-40-2*H2:2*H2:40+2*H2
        for z=-20-2*H1:2*H1:20+2*H1
            if ((x^2)/1600+(y^2)/900+(z^2)/400<1)||(((-H2+x)^2)/1600+((H3+y)^2)/900+(z^2)/400<1)||(((H2+x)^2)/1600+((H3+y)^2)/900+(z^2)/400<1)||(((H2*2+x)^2)/1600+(y^2)/900+(z^2)/400<1)||(((Hi+x)^2)/1600+((Hj+y)^2)/900+((H1+z)^2)/400<1)||(((-H2+Hi+x)^2)/1600+((H3+Hj+y)^2)/900+((H1+z)^2)/400<1)||(((H2+Hi+x)^2)/1600+((H3+Hj+y)^2)/900+((H1+z)^2)/400<1)||(((H2*2+Hi+x)^2)/1600+((Hj+y)^2)/900+((H1+z)^2)/400<1)
                bodycube(x,y,z,R)
            end
        end
    end
end

hold on
axis equal
grid on

效果图：

提示：会有一部分椭球体没被覆盖，因为小球球心已经超出椭球体了，一些椭球体外的小球也是因为没有分情况绘制导致的结果。