将外部支架生成网格以进行 3D 打印
Generating outside supporters into mesh for 3D printing
序言
这是我尝试重新询问已关闭的 Generating supporters for 3D printing，因为它是一个有趣的问题，但缺少重要的细节......这打算作为 Q&A，目前我正在编写答案的代码，但请随时回答（我接受最佳答案）。
问题描述
好的，关于这个问题的一些基本信息：
由于这是一个大问题，我将重点关注通用的 mesh/support-pattern 合并几何问题。
简而言之，如果我们想要打印任何网格，我们只能在它连接到起始平面达到 ~45 度角（+/- 对于不同的打印技术）时才能做到。因此，如果我们得到未连接到该平面的部件，我们需要创建一座桥，将 hold/connect 它连接到它。像这样（图片取自上面链接的页面）：
当然，我们需要添加尽可能少量的 material 并且它仍然足够坚固以将我们的网格固定到位而不会弯曲。最重要的是，我们需要削弱网格附近的支撑，以便在打印后轻松折断。
不要忘记形状和位置取决于许多因素，例如 material 和使用的技术、热流。
问题：
为了将这个庞大的话题缩小为可回答的问题，让我们只关注这个问题：
如何合并 3D 三角网格（边界表示如 STL）与预定义的支撑模式（如 3 边棱镜）从定义的平面垂直连接它？
使用简单的 C++.
好的，让我们从绝对基础开始。
支撑形状
您可以使用任何形状来满足所用印刷技术的具体要求。在 STL 中最容易生成的是 3 边棱柱形，它包含 2 个三角形底面（顶部和底部）和 3 个边，所有边都有 2 个三角形。所以总共有 8 个三角形。
此形状将从某个基本平面 (Z=0) 开始，然后上升直到碰到网格。然而，为了使这项工作有效，支撑必须在网格和自身之间有一个小的 gap，我们将在其中添加带有网格的弱化关节结构。
支持模式
这里有很多选项，所以我选择了最简单的（但不是防犯规的），那就是将支架放置在均匀的网格中，支架之间的距离为 grid。
所以只需从基平面上的每个网格位置向上投射一条射线，并检查与网格的交点。如果找到，请将支撑放置在高度刚好低于交点 gap 的位置。
关节
想法是将非常薄的支撑扇加入锥形连接并覆盖主支撑棱柱上方的支撑表面，角度小于 45 度（因此 gap 应该足够大以覆盖 grid这样的距离）。
这里的主要问题是我们必须细分我们连接的三角形，以便满足 STL 网格属性。为了解决连接问题（避免空洞或破坏 STL 的连接要求），我们可以使用不同的实体作为支撑，并为我们的网格使用不同的实体。这也将使我们能够在没有 re-triangulating 的情况下触摸表面，从而使这项任务变得容易得多。
为了简单起见，我选择了四面体形状，它很容易从三角形构造，并且在 mesh/support 关节处也存在弱点。
所以让我们进行一些测试 STL 网格并将其放置在我们的基准平面上方：
并放置我们的主要支撑：
还有关节：
这里 VCL/C++ 这个 STL3D.h 的代码:
//---------------------------------------------------------------------------
//--- simple STL 3D mesh ----------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _STL3D_h
#define _STL3D_h
//---------------------------------------------------------------------------
#ifdef ComctrlsHPP
TProgressBar *progress=NULL; // loading progress bar for realy big STL files
#endif
void _progress_init(int n);
void _progress (int ix);
void _progress_done();
//---------------------------------------------------------------------------
class STL3D // STL 3D mesh
{
public:
double center[3],size[3],rmax; // bbox center,half sizes, max(size[])
struct _fac
{
float p[3][3]; // triangle vertexes CCW order
float n[3]; // triangle unit normal pointing out
WORD attr;
_fac() {}
_fac(_fac& a) { *this=a; }
~_fac() {}
_fac* operator = (const _fac *a) { *this=*a; return this; }
//_fac* operator = (const _fac &a) { ...copy... return this; }
void compute() // compute normal
{
float a[3],b[3];
vectorf_sub(a,p[1],p[0]);
vectorf_sub(b,p[2],p[1]);
vectorf_mul(n,a,b);
vectorf_one(n,n);
}
double intersect_ray(double *pos,double *dir) // return -1 or distance to triangle and unit ray intersection
{
double p0[3],p1[3],p2[3]; // input triangle vertexes
double e1[3],e2[3],pp[3],qq[3],rr[3]; // dir must be unit vector !!!
double t,u,v,det,idet;
// get points
vector_ld(p0,p[0][0],p[0][1],p[0][2]);
vector_ld(p1,p[1][0],p[1][1],p[1][2]);
vector_ld(p2,p[2][0],p[2][1],p[2][2]);
//compute ray triangle intersection
vector_sub(e1,p1,p0);
vector_sub(e2,p2,p0);
// Calculate planes normal vector
vector_mul(pp,dir,e2);
det=vector_mul(e1,pp);
// Ray is parallel to plane
if (fabs(det)<1e-8) return -1.0;
idet=1.0/det;
vector_sub(rr,pos,p0);
u=vector_mul(rr,pp)*idet;
if ((u<0.0)||(u>1.0)) return -1.0;
vector_mul(qq,rr,e1);
v=vector_mul(dir,qq)*idet;
if ((v<0.0)||(u+v>1.0)) return -1.0;
// distance
t=vector_mul(e2,qq)*idet;
if (t<0.0) t=-1.0;
return t;
}
};
List<_fac> fac; // faces
STL3D() { reset(); }
STL3D(STL3D& a) { *this=a; }
~STL3D() {}
STL3D* operator = (const STL3D *a) { *this=*a; return this; }
//STL3D* operator = (const STL3D &a) { ...copy... return this; }
void reset(){ fac.num=0; compute(); } // clear STL
void draw(); // render STL mesh (OpenGL)
void draw_normals(float size); // render STL normals (OpenGL)
void compute(); // compute bbox
void compute_normals(); // recompute normals from points
void supports(reper &obj); // compute supports with obj placement above base plane z=0
void load(AnsiString name);
void save(AnsiString name);
};
//---------------------------------------------------------------------------
void STL3D::draw()
{
_fac *f; int i,j; BYTE r,g,b;
glBegin(GL_TRIANGLES);
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
glNormal3fv(f->n);
if (f->attr<32768)
{
r= f->attr &31; r<<=3;
g=(f->attr>> 5)&31; g<<=3;
b=(f->attr>>10)&31; b<<=3;
glColor3ub(r,g,b);
}
for (j=0;j<3;j++) glVertex3fv(f->p[j]);
}
glEnd();
}
//---------------------------------------------------------------------------
void STL3D::draw_normals(float size)
{
_fac *f;
int i; float a[3],b[3];
glBegin(GL_LINES);
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
vectorf_add(a,f->p[0],f->p[1]);
vectorf_add(a,a ,f->p[2]);
vectorf_mul(a,a,1.0/3.0);
vectorf_mul(b,f->n,size); glVertex3fv(a);
vectorf_add(b,b,a); glVertex3fv(b);
}
glEnd();
}
//---------------------------------------------------------------------------
void STL3D::compute()
{
_fac *f;
int i,j,k;
double p0[3],p1[3];
vector_ld(center,0.0,0.0,0.0);
vector_ld(size,0.0,0.0,0.0);
rmax=0.0;
if (fac.num==0) return;
// bbox
for (k=0;k<3;k++) p0[k]=fac.dat[0].p[0][k];
for (k=0;k<3;k++) p1[k]=fac.dat[0].p[0][k];
for (f=fac.dat,i=0;i<fac.num;i++,f++)
for (j=0;j<3;j++)
for (k=0;k<3;k++)
{
if (p0[k]>f->p[j][k]) p0[k]=f->p[j][k];
if (p1[k]<f->p[j][k]) p1[k]=f->p[j][k];
}
vector_add(center,p0,p1); vector_mul(center,center,0.5);
vector_sub(size ,p1,p0); vector_mul(size ,size ,0.5);
rmax=size[0];
if (rmax<size[1]) rmax=size[1];
if (rmax<size[2]) rmax=size[2];
// attr repair
for (f=fac.dat,i=0;i<fac.num;i++,f++)
if (f->attr==0) f->attr=32768;
}
//---------------------------------------------------------------------------
void STL3D::compute_normals()
{
_fac *f; int i;
for (f=fac.dat,i=0;i<fac.num;i++,f++) f->compute();
}
//---------------------------------------------------------------------------
void STL3D::supports(reper &obj)
{
_fac *f,ff;
int i,j,k;
double p[3],dp[3],x0,y0,h0,x1,y1,x2,y2,h1,t;
// some config values first
const WORD attr0=31<<10; // support attr should be different than joint
const WORD attr1=31<<5; // joint attr should be different than mesh,support
const double grid0=8.0; // distance between supports
const double grid1=2.0; // distance between joints
const double gap=grid0/tan(45.0*deg);// distance between main support and mesh (joint size)
const double ha=1.0; // main support side size
// do not mess with these
const double hx= ha*cos(60.0*deg); // half size of main support in x
const double hy=0.5*ha*sin(60.0*deg); // half size of main support in y
const double grid2=0.4*hy; // distance between joints bases
const double ga=2.0*grid2*grid1/grid0; // main support side size
const double gx=hx*grid2/grid0; // half size of joint support in x
const double gy=hy*grid2/grid0; // half size of joint support in y
// apply placement obj (may lose some accuracy) not needed if matrices are not used
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
for (j=0;j<3;j++)
{
for (k=0;k<3;k++) p[k]=f->p[j][k]; // float->double
obj.l2g(p,p);
for (k=0;k<3;k++) f->p[j][k]=p[k]; // double->float
}
for (k=0;k<3;k++) p[k]=f->n[k]; // float->double
obj.l2g_dir(p,p);
for (k=0;k<3;k++) f->n[k]=p[k]; // double->float
} compute();
// create supports
for (x0=center[0]-size[0]+(0.5*grid0);x0<=center[0]+size[0]-(0.5*grid0);x0+=grid0)
for (y0=center[1]-size[1]+(0.5*grid0);y0<=center[1]+size[1]-(0.5*grid0);y0+=grid0)
{
// cast ray x0,y0,0 in Z+ direction to check for mesh intersection to compute the support height h0
h0=center[2]+size[2]+1e6;
vector_ld(p,x0,y0,0.0);
vector_ld(dp,0.0,0.0,+1.0);
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
t=f->intersect_ray(p,dp);
if ((t>=0.0)&&(t<h0)) h0=t;
}
if (h0>center[2]+size[2]+1e5) continue; // skip non intersected rays
h0-=gap; if (h0<0.0) h0=0.0;
// main suport prism
ff.attr=attr0;
// sides
ff.attr=attr0;
vectorf_ld(ff.p[0],x0-hx,y0-hy,0.0);
vectorf_ld(ff.p[1],x0+hx,y0-hy,0.0);
vectorf_ld(ff.p[2],x0-hx,y0-hy, h0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0+hx,y0-hy,0.0);
vectorf_ld(ff.p[1],x0+hx,y0-hy, h0);
vectorf_ld(ff.p[2],x0-hx,y0-hy, h0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0-hx,y0-hy, h0);
vectorf_ld(ff.p[1],x0 ,y0+hy,0.0);
vectorf_ld(ff.p[2],x0-hx,y0-hy,0.0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0-hx,y0-hy, h0);
vectorf_ld(ff.p[1],x0 ,y0+hy, h0);
vectorf_ld(ff.p[2],x0 ,y0+hy,0.0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0 ,y0+hy, h0);
vectorf_ld(ff.p[1],x0+hx,y0-hy,0.0);
vectorf_ld(ff.p[2],x0 ,y0+hy,0.0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0 ,y0+hy, h0);
vectorf_ld(ff.p[1],x0+hx,y0-hy, h0);
vectorf_ld(ff.p[2],x0+hx,y0-hy,0.0); ff.compute(); fac.add(ff);
// base triangles
vectorf_ld(ff.p[0],x0 ,y0+hy,0.0);
vectorf_ld(ff.p[1],x0+hx,y0-hy,0.0);
vectorf_ld(ff.p[2],x0-hx,y0-hy,0.0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0-hx,y0-hy, h0);
vectorf_ld(ff.p[1],x0+hx,y0-hy, h0);
vectorf_ld(ff.p[2],x0 ,y0+hy, h0); ff.compute(); fac.add(ff);
// joints
for (x1=x0-(0.5*grid0),x2=x0-(0.5*grid2);x1<=x0+(0.5*grid0);x1+=grid1,x2+=ga)
for (y1=y0-(0.5*grid0),y2=y0-(1.9*grid2);y1<=y0+(0.5*grid0);y1+=grid1,y2+=ga)
{
// cast ray x1,y1,0 in Z+ direction to check for mesh intersection to compute the joint height h1
h1=h0+gap+1e6;
vector_ld(p,x1,y1,0.0);
vector_ld(dp,0.0,0.0,+1.0);
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
t=f->intersect_ray(p,dp);
if ((t>=0.0)&&(t<h1)) h1=t;
}
if (h1>h0+gap+1e5) continue; // skip non intersected rays
// tetrahedron joints
ff.attr=attr1;
// base triangle
vectorf_ld(ff.p[0],x2 ,y2+gy,h0);
vectorf_ld(ff.p[1],x2+gx,y2-gy,h0);
vectorf_ld(ff.p[2],x2-gx,y2-gy,h0); ff.compute(); fac.add(ff);
// sides
vectorf_ld(ff.p[0],x2+gx,y2-gy,h0);
vectorf_ld(ff.p[1],x2 ,y2+gy,h0);
vectorf_ld(ff.p[2],x1 ,y1 ,h1); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x2 ,y2+gy,h0);
vectorf_ld(ff.p[1],x2-gx,y2-gy,h0);
vectorf_ld(ff.p[2],x1 ,y1 ,h1); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x2+gx,y2+gy,h0);
vectorf_ld(ff.p[1],x2-gx,y2-gy,h0);
vectorf_ld(ff.p[2],x1 ,y1 ,h1); ff.compute(); fac.add(ff);
}
}
// reverse placement obj (may lose some accuracy) not needed if matrices are not used
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
for (j=0;j<3;j++)
{
for (k=0;k<3;k++) p[k]=f->p[j][k]; // float->double
obj.g2l(p,p);
for (k=0;k<3;k++) f->p[j][k]=p[k]; // double->float
}
for (k=0;k<3;k++) p[k]=f->n[k]; // float->double
obj.g2l_dir(p,p);
for (k=0;k<3;k++) f->n[k]=p[k]; // double->float
} compute();
}
//---------------------------------------------------------------------------
void STL3D::load(AnsiString name)
{
int adr,siz,hnd;
BYTE *dat;
AnsiString lin,s;
int i,j,l,n;
_fac f;
reset(); f.attr=0;
siz=0;
hnd=FileOpen(name,fmOpenRead);
if (hnd<0) return;
siz=FileSeek(hnd,0,2);
FileSeek(hnd,0,0);
dat=new BYTE[siz];
if (dat==NULL) { FileClose(hnd); return; }
FileRead(hnd,dat,siz);
FileClose(hnd);
adr=0; s=txt_load_str(dat,siz,adr,true);
// ASCII
if (s=="solid")
{
_progress_init(siz); int progress_cnt=0;
for (adr=0;adr<siz;)
{
progress_cnt++; if (progress_cnt>=128) { progress_cnt=0; _progress(adr); }
lin=txt_load_lin(dat,siz,adr,true);
for (i=1,l=lin.Length();i<=l;)
{
s=str_load_str(lin,i,true);
if (s=="solid") { name=str_load_str(lin,i,true); break; }
if (s=="endsolid") break;
if (s=="facet")
{
j=0;
s=str_load_str(lin,i,true);
f.n[0]=str2num(str_load_str(lin,i,true));
f.n[1]=str2num(str_load_str(lin,i,true));
f.n[2]=str2num(str_load_str(lin,i,true));
}
if (s=="vertex")
if (j<3)
{
f.p[j][0]=str2num(str_load_str(lin,i,true));
f.p[j][1]=str2num(str_load_str(lin,i,true));
f.p[j][2]=str2num(str_load_str(lin,i,true));
j++;
if (j==3) fac.add(f);
}
break;
}
}
}
// binary
else{
adr=80;
n=((DWORD*)(dat+adr))[0]; adr+=4;
fac.allocate(n); fac.num=0;
_progress_init(n); int progress_cnt=0;
for (i=0;i<n;i++)
{
if (adr+50>siz) break; // error
progress_cnt++; if (progress_cnt>=128) { progress_cnt=0; _progress(i); }
f.n[0]=((float*)(dat+adr))[0]; adr+=4;
f.n[1]=((float*)(dat+adr))[0]; adr+=4;
f.n[2]=((float*)(dat+adr))[0]; adr+=4;
for (j=0;j<3;j++)
{
f.p[j][0]=((float*)(dat+adr))[0]; adr+=4;
f.p[j][1]=((float*)(dat+adr))[0]; adr+=4;
f.p[j][2]=((float*)(dat+adr))[0]; adr+=4;
}
f.attr=((WORD*)(dat+adr))[0]; adr+=2; // attributes
fac.add(f);
}
}
_progress_done();
delete[] dat;
compute();
}
//---------------------------------------------------------------------------
void STL3D::save(AnsiString name)
{
// ToDo
}
//---------------------------------------------------------------------------
void _progress_init(int n)
{
#ifdef ComctrlsHPP
if (progress==NULL) return;
progress->Position=0;
progress->Max=n;
progress->Visible=true;
#endif
}
//---------------------------------------------------------------------------
void _progress (int ix)
{
#ifdef ComctrlsHPP
if (progress==NULL) return;
progress->Position=ix;
progress->Update();
#endif
}
//---------------------------------------------------------------------------
void _progress_done()
{
#ifdef ComctrlsHPP
if (progress==NULL) return;
progress->Visible=false;
#endif
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
用法很简单：
#include "STL3D.h" // STL mesh (this is the important stuff)
STL3D mesh; // point cloud and tetrahedronal mesh
mesh.load("space_invader_magnet.stl");
mesh.supports(obj); // obj is object holding 4x4 uniform matrix of placement if you STL is already placed than it is not needed
我使用了很多来自我的 OpenGL 引擎的东西，比如动态 List<> 模板：
List<double> xxx; 等同于 double xxx[];
xxx.add(5); 将 5 添加到列表末尾
xxx[7] 访问数组元素（安全）
xxx.dat[7]访问数组元素（不安全但直接访问速度快）
xxx.num是数组实际使用的大小
xxx.reset()清空数组并设置xxx.num=0
xxx.allocate(100) 为 100 项预分配 space
或矢量和矩阵数学（vectorf_ 与 float* 一起工作，vector_ 与 double 一起工作）这不太重要。如果您需要数学，请参阅：
如果 STL 已经放置（无矩阵），则根本不需要放置转换，也不需要 obj。该代码反映了上面的项目符号。我想让它尽可能简单，所以还没有优化。
gap 和 grid 常量在支持函数中 hard-coded，尚未设置为有效值。
[备注]
现在这几乎只涵盖了问题的最基本部分，还有很多边缘情况 un-handled 需要保留这个 "short"。代码本身不检查三角形是否高于 45 度斜率，但可以通过简单的法线角度检查来完成，例如：
if (acos(dot(normal,(0.0,0.0,1.0))<45.0*deg) continue;
还需要在网格的各个部分之间添加支撑，例如，如果您的对象的层数多于基平面仅支持第一层。其余部分必须使用其下方的层......并在支撑两侧使用弱化接头。这类似于放置第一层支撑，你只需要在两个方向上投射光线......的点积）。例如，这是可能需要的网格放置（对于某些技术）：
在设计支撑时请记住，您应该满足打印过程的正确缠绕规则 (CCW) 和法线方向（向外）...
序言
这是我尝试重新询问已关闭的 Generating supporters for 3D printing，因为它是一个有趣的问题，但缺少重要的细节......这打算作为 Q&A，目前我正在编写答案的代码，但请随时回答（我接受最佳答案）。
问题描述
好的，关于这个问题的一些基本信息：
由于这是一个大问题，我将重点关注通用的 mesh/support-pattern 合并几何问题。
简而言之，如果我们想要打印任何网格，我们只能在它连接到起始平面达到 ~45 度角（+/- 对于不同的打印技术）时才能做到。因此，如果我们得到未连接到该平面的部件，我们需要创建一座桥，将 hold/connect 它连接到它。像这样（图片取自上面链接的页面）：
当然，我们需要添加尽可能少量的 material 并且它仍然足够坚固以将我们的网格固定到位而不会弯曲。最重要的是，我们需要削弱网格附近的支撑，以便在打印后轻松折断。
不要忘记形状和位置取决于许多因素，例如 material 和使用的技术、热流。
问题：
为了将这个庞大的话题缩小为可回答的问题，让我们只关注这个问题：
如何合并 3D 三角网格（边界表示如 STL）与预定义的支撑模式（如 3 边棱镜）从定义的平面垂直连接它？
使用简单的 C++.
好的，让我们从绝对基础开始。
支撑形状
您可以使用任何形状来满足所用印刷技术的具体要求。在 STL 中最容易生成的是 3 边棱柱形，它包含 2 个三角形底面（顶部和底部）和 3 个边，所有边都有 2 个三角形。所以总共有 8 个三角形。
此形状将从某个基本平面 (
Z=0) 开始，然后上升直到碰到网格。然而，为了使这项工作有效，支撑必须在网格和自身之间有一个小的gap，我们将在其中添加带有网格的弱化关节结构。支持模式
这里有很多选项，所以我选择了最简单的（但不是防犯规的），那就是将支架放置在均匀的网格中，支架之间的距离为
grid。所以只需从基平面上的每个网格位置向上投射一条射线，并检查与网格的交点。如果找到，请将支撑放置在高度刚好低于交点
gap的位置。关节
想法是将非常薄的支撑扇加入锥形连接并覆盖主支撑棱柱上方的支撑表面，角度小于 45 度（因此
gap应该足够大以覆盖grid这样的距离）。这里的主要问题是我们必须细分我们连接的三角形，以便满足 STL 网格属性。为了解决连接问题（避免空洞或破坏 STL 的连接要求），我们可以使用不同的实体作为支撑，并为我们的网格使用不同的实体。这也将使我们能够在没有 re-triangulating 的情况下触摸表面，从而使这项任务变得容易得多。
为了简单起见，我选择了四面体形状，它很容易从三角形构造，并且在 mesh/support 关节处也存在弱点。
所以让我们进行一些测试 STL 网格并将其放置在我们的基准平面上方：
并放置我们的主要支撑：
还有关节：
这里 VCL/C++ 这个 STL3D.h 的代码:
//---------------------------------------------------------------------------
//--- simple STL 3D mesh ----------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _STL3D_h
#define _STL3D_h
//---------------------------------------------------------------------------
#ifdef ComctrlsHPP
TProgressBar *progress=NULL; // loading progress bar for realy big STL files
#endif
void _progress_init(int n);
void _progress (int ix);
void _progress_done();
//---------------------------------------------------------------------------
class STL3D // STL 3D mesh
{
public:
double center[3],size[3],rmax; // bbox center,half sizes, max(size[])
struct _fac
{
float p[3][3]; // triangle vertexes CCW order
float n[3]; // triangle unit normal pointing out
WORD attr;
_fac() {}
_fac(_fac& a) { *this=a; }
~_fac() {}
_fac* operator = (const _fac *a) { *this=*a; return this; }
//_fac* operator = (const _fac &a) { ...copy... return this; }
void compute() // compute normal
{
float a[3],b[3];
vectorf_sub(a,p[1],p[0]);
vectorf_sub(b,p[2],p[1]);
vectorf_mul(n,a,b);
vectorf_one(n,n);
}
double intersect_ray(double *pos,double *dir) // return -1 or distance to triangle and unit ray intersection
{
double p0[3],p1[3],p2[3]; // input triangle vertexes
double e1[3],e2[3],pp[3],qq[3],rr[3]; // dir must be unit vector !!!
double t,u,v,det,idet;
// get points
vector_ld(p0,p[0][0],p[0][1],p[0][2]);
vector_ld(p1,p[1][0],p[1][1],p[1][2]);
vector_ld(p2,p[2][0],p[2][1],p[2][2]);
//compute ray triangle intersection
vector_sub(e1,p1,p0);
vector_sub(e2,p2,p0);
// Calculate planes normal vector
vector_mul(pp,dir,e2);
det=vector_mul(e1,pp);
// Ray is parallel to plane
if (fabs(det)<1e-8) return -1.0;
idet=1.0/det;
vector_sub(rr,pos,p0);
u=vector_mul(rr,pp)*idet;
if ((u<0.0)||(u>1.0)) return -1.0;
vector_mul(qq,rr,e1);
v=vector_mul(dir,qq)*idet;
if ((v<0.0)||(u+v>1.0)) return -1.0;
// distance
t=vector_mul(e2,qq)*idet;
if (t<0.0) t=-1.0;
return t;
}
};
List<_fac> fac; // faces
STL3D() { reset(); }
STL3D(STL3D& a) { *this=a; }
~STL3D() {}
STL3D* operator = (const STL3D *a) { *this=*a; return this; }
//STL3D* operator = (const STL3D &a) { ...copy... return this; }
void reset(){ fac.num=0; compute(); } // clear STL
void draw(); // render STL mesh (OpenGL)
void draw_normals(float size); // render STL normals (OpenGL)
void compute(); // compute bbox
void compute_normals(); // recompute normals from points
void supports(reper &obj); // compute supports with obj placement above base plane z=0
void load(AnsiString name);
void save(AnsiString name);
};
//---------------------------------------------------------------------------
void STL3D::draw()
{
_fac *f; int i,j; BYTE r,g,b;
glBegin(GL_TRIANGLES);
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
glNormal3fv(f->n);
if (f->attr<32768)
{
r= f->attr &31; r<<=3;
g=(f->attr>> 5)&31; g<<=3;
b=(f->attr>>10)&31; b<<=3;
glColor3ub(r,g,b);
}
for (j=0;j<3;j++) glVertex3fv(f->p[j]);
}
glEnd();
}
//---------------------------------------------------------------------------
void STL3D::draw_normals(float size)
{
_fac *f;
int i; float a[3],b[3];
glBegin(GL_LINES);
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
vectorf_add(a,f->p[0],f->p[1]);
vectorf_add(a,a ,f->p[2]);
vectorf_mul(a,a,1.0/3.0);
vectorf_mul(b,f->n,size); glVertex3fv(a);
vectorf_add(b,b,a); glVertex3fv(b);
}
glEnd();
}
//---------------------------------------------------------------------------
void STL3D::compute()
{
_fac *f;
int i,j,k;
double p0[3],p1[3];
vector_ld(center,0.0,0.0,0.0);
vector_ld(size,0.0,0.0,0.0);
rmax=0.0;
if (fac.num==0) return;
// bbox
for (k=0;k<3;k++) p0[k]=fac.dat[0].p[0][k];
for (k=0;k<3;k++) p1[k]=fac.dat[0].p[0][k];
for (f=fac.dat,i=0;i<fac.num;i++,f++)
for (j=0;j<3;j++)
for (k=0;k<3;k++)
{
if (p0[k]>f->p[j][k]) p0[k]=f->p[j][k];
if (p1[k]<f->p[j][k]) p1[k]=f->p[j][k];
}
vector_add(center,p0,p1); vector_mul(center,center,0.5);
vector_sub(size ,p1,p0); vector_mul(size ,size ,0.5);
rmax=size[0];
if (rmax<size[1]) rmax=size[1];
if (rmax<size[2]) rmax=size[2];
// attr repair
for (f=fac.dat,i=0;i<fac.num;i++,f++)
if (f->attr==0) f->attr=32768;
}
//---------------------------------------------------------------------------
void STL3D::compute_normals()
{
_fac *f; int i;
for (f=fac.dat,i=0;i<fac.num;i++,f++) f->compute();
}
//---------------------------------------------------------------------------
void STL3D::supports(reper &obj)
{
_fac *f,ff;
int i,j,k;
double p[3],dp[3],x0,y0,h0,x1,y1,x2,y2,h1,t;
// some config values first
const WORD attr0=31<<10; // support attr should be different than joint
const WORD attr1=31<<5; // joint attr should be different than mesh,support
const double grid0=8.0; // distance between supports
const double grid1=2.0; // distance between joints
const double gap=grid0/tan(45.0*deg);// distance between main support and mesh (joint size)
const double ha=1.0; // main support side size
// do not mess with these
const double hx= ha*cos(60.0*deg); // half size of main support in x
const double hy=0.5*ha*sin(60.0*deg); // half size of main support in y
const double grid2=0.4*hy; // distance between joints bases
const double ga=2.0*grid2*grid1/grid0; // main support side size
const double gx=hx*grid2/grid0; // half size of joint support in x
const double gy=hy*grid2/grid0; // half size of joint support in y
// apply placement obj (may lose some accuracy) not needed if matrices are not used
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
for (j=0;j<3;j++)
{
for (k=0;k<3;k++) p[k]=f->p[j][k]; // float->double
obj.l2g(p,p);
for (k=0;k<3;k++) f->p[j][k]=p[k]; // double->float
}
for (k=0;k<3;k++) p[k]=f->n[k]; // float->double
obj.l2g_dir(p,p);
for (k=0;k<3;k++) f->n[k]=p[k]; // double->float
} compute();
// create supports
for (x0=center[0]-size[0]+(0.5*grid0);x0<=center[0]+size[0]-(0.5*grid0);x0+=grid0)
for (y0=center[1]-size[1]+(0.5*grid0);y0<=center[1]+size[1]-(0.5*grid0);y0+=grid0)
{
// cast ray x0,y0,0 in Z+ direction to check for mesh intersection to compute the support height h0
h0=center[2]+size[2]+1e6;
vector_ld(p,x0,y0,0.0);
vector_ld(dp,0.0,0.0,+1.0);
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
t=f->intersect_ray(p,dp);
if ((t>=0.0)&&(t<h0)) h0=t;
}
if (h0>center[2]+size[2]+1e5) continue; // skip non intersected rays
h0-=gap; if (h0<0.0) h0=0.0;
// main suport prism
ff.attr=attr0;
// sides
ff.attr=attr0;
vectorf_ld(ff.p[0],x0-hx,y0-hy,0.0);
vectorf_ld(ff.p[1],x0+hx,y0-hy,0.0);
vectorf_ld(ff.p[2],x0-hx,y0-hy, h0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0+hx,y0-hy,0.0);
vectorf_ld(ff.p[1],x0+hx,y0-hy, h0);
vectorf_ld(ff.p[2],x0-hx,y0-hy, h0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0-hx,y0-hy, h0);
vectorf_ld(ff.p[1],x0 ,y0+hy,0.0);
vectorf_ld(ff.p[2],x0-hx,y0-hy,0.0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0-hx,y0-hy, h0);
vectorf_ld(ff.p[1],x0 ,y0+hy, h0);
vectorf_ld(ff.p[2],x0 ,y0+hy,0.0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0 ,y0+hy, h0);
vectorf_ld(ff.p[1],x0+hx,y0-hy,0.0);
vectorf_ld(ff.p[2],x0 ,y0+hy,0.0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0 ,y0+hy, h0);
vectorf_ld(ff.p[1],x0+hx,y0-hy, h0);
vectorf_ld(ff.p[2],x0+hx,y0-hy,0.0); ff.compute(); fac.add(ff);
// base triangles
vectorf_ld(ff.p[0],x0 ,y0+hy,0.0);
vectorf_ld(ff.p[1],x0+hx,y0-hy,0.0);
vectorf_ld(ff.p[2],x0-hx,y0-hy,0.0); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x0-hx,y0-hy, h0);
vectorf_ld(ff.p[1],x0+hx,y0-hy, h0);
vectorf_ld(ff.p[2],x0 ,y0+hy, h0); ff.compute(); fac.add(ff);
// joints
for (x1=x0-(0.5*grid0),x2=x0-(0.5*grid2);x1<=x0+(0.5*grid0);x1+=grid1,x2+=ga)
for (y1=y0-(0.5*grid0),y2=y0-(1.9*grid2);y1<=y0+(0.5*grid0);y1+=grid1,y2+=ga)
{
// cast ray x1,y1,0 in Z+ direction to check for mesh intersection to compute the joint height h1
h1=h0+gap+1e6;
vector_ld(p,x1,y1,0.0);
vector_ld(dp,0.0,0.0,+1.0);
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
t=f->intersect_ray(p,dp);
if ((t>=0.0)&&(t<h1)) h1=t;
}
if (h1>h0+gap+1e5) continue; // skip non intersected rays
// tetrahedron joints
ff.attr=attr1;
// base triangle
vectorf_ld(ff.p[0],x2 ,y2+gy,h0);
vectorf_ld(ff.p[1],x2+gx,y2-gy,h0);
vectorf_ld(ff.p[2],x2-gx,y2-gy,h0); ff.compute(); fac.add(ff);
// sides
vectorf_ld(ff.p[0],x2+gx,y2-gy,h0);
vectorf_ld(ff.p[1],x2 ,y2+gy,h0);
vectorf_ld(ff.p[2],x1 ,y1 ,h1); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x2 ,y2+gy,h0);
vectorf_ld(ff.p[1],x2-gx,y2-gy,h0);
vectorf_ld(ff.p[2],x1 ,y1 ,h1); ff.compute(); fac.add(ff);
vectorf_ld(ff.p[0],x2+gx,y2+gy,h0);
vectorf_ld(ff.p[1],x2-gx,y2-gy,h0);
vectorf_ld(ff.p[2],x1 ,y1 ,h1); ff.compute(); fac.add(ff);
}
}
// reverse placement obj (may lose some accuracy) not needed if matrices are not used
for (f=fac.dat,i=0;i<fac.num;i++,f++)
{
for (j=0;j<3;j++)
{
for (k=0;k<3;k++) p[k]=f->p[j][k]; // float->double
obj.g2l(p,p);
for (k=0;k<3;k++) f->p[j][k]=p[k]; // double->float
}
for (k=0;k<3;k++) p[k]=f->n[k]; // float->double
obj.g2l_dir(p,p);
for (k=0;k<3;k++) f->n[k]=p[k]; // double->float
} compute();
}
//---------------------------------------------------------------------------
void STL3D::load(AnsiString name)
{
int adr,siz,hnd;
BYTE *dat;
AnsiString lin,s;
int i,j,l,n;
_fac f;
reset(); f.attr=0;
siz=0;
hnd=FileOpen(name,fmOpenRead);
if (hnd<0) return;
siz=FileSeek(hnd,0,2);
FileSeek(hnd,0,0);
dat=new BYTE[siz];
if (dat==NULL) { FileClose(hnd); return; }
FileRead(hnd,dat,siz);
FileClose(hnd);
adr=0; s=txt_load_str(dat,siz,adr,true);
// ASCII
if (s=="solid")
{
_progress_init(siz); int progress_cnt=0;
for (adr=0;adr<siz;)
{
progress_cnt++; if (progress_cnt>=128) { progress_cnt=0; _progress(adr); }
lin=txt_load_lin(dat,siz,adr,true);
for (i=1,l=lin.Length();i<=l;)
{
s=str_load_str(lin,i,true);
if (s=="solid") { name=str_load_str(lin,i,true); break; }
if (s=="endsolid") break;
if (s=="facet")
{
j=0;
s=str_load_str(lin,i,true);
f.n[0]=str2num(str_load_str(lin,i,true));
f.n[1]=str2num(str_load_str(lin,i,true));
f.n[2]=str2num(str_load_str(lin,i,true));
}
if (s=="vertex")
if (j<3)
{
f.p[j][0]=str2num(str_load_str(lin,i,true));
f.p[j][1]=str2num(str_load_str(lin,i,true));
f.p[j][2]=str2num(str_load_str(lin,i,true));
j++;
if (j==3) fac.add(f);
}
break;
}
}
}
// binary
else{
adr=80;
n=((DWORD*)(dat+adr))[0]; adr+=4;
fac.allocate(n); fac.num=0;
_progress_init(n); int progress_cnt=0;
for (i=0;i<n;i++)
{
if (adr+50>siz) break; // error
progress_cnt++; if (progress_cnt>=128) { progress_cnt=0; _progress(i); }
f.n[0]=((float*)(dat+adr))[0]; adr+=4;
f.n[1]=((float*)(dat+adr))[0]; adr+=4;
f.n[2]=((float*)(dat+adr))[0]; adr+=4;
for (j=0;j<3;j++)
{
f.p[j][0]=((float*)(dat+adr))[0]; adr+=4;
f.p[j][1]=((float*)(dat+adr))[0]; adr+=4;
f.p[j][2]=((float*)(dat+adr))[0]; adr+=4;
}
f.attr=((WORD*)(dat+adr))[0]; adr+=2; // attributes
fac.add(f);
}
}
_progress_done();
delete[] dat;
compute();
}
//---------------------------------------------------------------------------
void STL3D::save(AnsiString name)
{
// ToDo
}
//---------------------------------------------------------------------------
void _progress_init(int n)
{
#ifdef ComctrlsHPP
if (progress==NULL) return;
progress->Position=0;
progress->Max=n;
progress->Visible=true;
#endif
}
//---------------------------------------------------------------------------
void _progress (int ix)
{
#ifdef ComctrlsHPP
if (progress==NULL) return;
progress->Position=ix;
progress->Update();
#endif
}
//---------------------------------------------------------------------------
void _progress_done()
{
#ifdef ComctrlsHPP
if (progress==NULL) return;
progress->Visible=false;
#endif
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
用法很简单：
#include "STL3D.h" // STL mesh (this is the important stuff)
STL3D mesh; // point cloud and tetrahedronal mesh
mesh.load("space_invader_magnet.stl");
mesh.supports(obj); // obj is object holding 4x4 uniform matrix of placement if you STL is already placed than it is not needed
我使用了很多来自我的 OpenGL 引擎的东西，比如动态 List<> 模板：
List<double> xxx; 等同于 double xxx[];
xxx.add(5); 将 5 添加到列表末尾
xxx[7] 访问数组元素（安全）
xxx.dat[7]访问数组元素（不安全但直接访问速度快）
xxx.num是数组实际使用的大小
xxx.reset()清空数组并设置xxx.num=0
xxx.allocate(100) 为 100 项预分配 space
或矢量和矩阵数学（vectorf_ 与 float* 一起工作，vector_ 与 double 一起工作）这不太重要。如果您需要数学，请参阅：
如果 STL 已经放置（无矩阵），则根本不需要放置转换，也不需要 obj。该代码反映了上面的项目符号。我想让它尽可能简单，所以还没有优化。
gap 和 grid 常量在支持函数中 hard-coded，尚未设置为有效值。
[备注]
现在这几乎只涵盖了问题的最基本部分，还有很多边缘情况 un-handled 需要保留这个 "short"。代码本身不检查三角形是否高于 45 度斜率，但可以通过简单的法线角度检查来完成，例如：
if (acos(dot(normal,(0.0,0.0,1.0))<45.0*deg) continue;
还需要在网格的各个部分之间添加支撑，例如，如果您的对象的层数多于基平面仅支持第一层。其余部分必须使用其下方的层......并在支撑两侧使用弱化接头。这类似于放置第一层支撑，你只需要在两个方向上投射光线......的点积）。例如，这是可能需要的网格放置（对于某些技术）：
在设计支撑时请记住，您应该满足打印过程的正确缠绕规则 (CCW) 和法线方向（向外）...