首页 > 解决方案 > 在曲面上绘制螺旋曲线

问题描述

我对在这样的表面上绘制螺旋曲线很感兴趣: 一条裤子

它的等式是:

((1.0-z)*(((x-1.0)*(x-1.0))+y*y-(1.0/3.0))*(((x+1.0)*(x+1.0))+(y*y)-(1.0/3.0)))+(z*((x*x)+(y*y)-(1.0/3.0))) == 0.0

在哪里

x=<-1.0-sqrt(1/3),+1.0+sqrt(1/3)>
y=<-sqrt(1/3),+sqrt(1/3)>
z=< 0.0,1.0>

我在这个网站上找到了等式:

https://www.quora.com/What-is-the-mathematical-expression-which-when-plotted-looks-like-a-pair-of-pants

然后我看到了那个帖子,StackMath但我不明白答案......:

https://math.stackexchange.com/questions/3267636/drawing-a-pair-of-pants-using-python?answertab=active#tab-top

如果有人知道在那种表面上绘制螺旋曲线的过程,那将对我有很大帮助。不要犹豫,问我更多关于这个问题的信息。谢谢你。

标签: 3dgeometrycomputational-geometrygeometry-surface

解决方案


分析方法并不是解决此问题的唯一方法(无论如何都不适合此站点),因此还有很多其他选择。但是您没有指定螺旋应该是什么样子,因为形状有 3 个末端,而螺旋只有 2 个......

为简单起见,我假设螺旋信封(就像你在形状周围缠绕一根绳子)。

那么如何攻击这个。2D <--> 3D通常的方法是在一些矩形 2D 区域(纹理)和您的表面之间创建映射。通常利用形状的拓扑结构,如圆柱/球坐标等。

另一种选择是将您的形状切割成单独的切片并将每个切片处理为 2D ......这大大简化了事情。

要将螺旋线应用到某个曲面上,您要么将其投影,要么在由参数方程定义的螺旋线上的实际点的形状和方向之间找到交点。

这里使用 OpenGL 可视化的 C++ 示例:

//---------------------------------------------------------------------------
List<int> slice;        // start index of each z slice in pnt
List<double> pnt;       // GL_POINTS surface points
List<int> lin;          // LINE_STRIP spiral point indexes
//---------------------------------------------------------------------------
void obj_init()
    {
    const double ry=sqrt(1.0/3.0);
    const double rx=1.0+ry;
    double a,x,y,z,d,N=7.0,da,dx,dy,dz,r,rr,_zero;
    int i,j,i0,i1,ii;
    // get points of analytical surface
    pnt.num=0;
    slice.num=0;
    d=0.02; dz=0.25*d;_zero=1e-2;
    for (z=0.0;z<=1.0;z+=dz,slice.add(pnt.num))
     for (x=-rx;x<=rx;x+=d)
      for (y=-ry;y<=ry;y+=d)
       if (fabs(((1.0-z)*(((x-1.0)*(x-1.0))+y*y-(1.0/3.0))*(((x+1.0)*(x+1.0))+(y*y)-(1.0/3.0)))+(z*((x*x)+(y*y)-(1.0/3.0))))<=_zero)
        {
        pnt.add(x);
        pnt.add(y);
        pnt.add(4.0*z-2.0);
        }
    // spiral line
    lin.num=0;
    da=5.0*M_PI/180.0;
    d=pnt[slice[1]+2]-pnt[slice[0]+2];
    for (a=0.0;a<=2.0*M_PI*N;a+=da)     // go through whole spiral
        {
        z=a/(2.0*M_PI*N);               // z = f(angle,screws)
        z=4.0*z-2.0;

        j=((z+2.0)/d);
        i0=j-1; if (i0<0) i0=0;     if (i0>=slice.num) break;
        i1=j+1; if (i1<0) continue; if (i1>=slice.num) i1=slice.num-1;
        i0=slice.dat[i0];
        i1=slice.dat[i1];

        dx=cos(a);                      // direction vector to the spiral point
        dy=sin(a);
        dz=z;

        for (rr=0.0,i=i0;i<i1;i+=3)         // find point with biggest distance from center and close to a,dx,dy
            {
            x=pnt.dat[i+0];
            y=pnt.dat[i+1];
            r=(x*dx)+(y*dy)+(z*dz);     // dot( (dx,dy,dz) , (x,y,z) )
            if (r>rr){ii=i; rr=r; }     // better point found
            }
        if (ii>=0) lin.add(ii);         // add spiral point to lin[]
        }
    }
//---------------------------------------------------------------------------
void obj_draw()
    {
    int i,j;
    glColor3f(0.3,0.3,0.3);
    glBegin(GL_POINTS);
    for (i=0;i<pnt.num;i+=3) glVertex3dv(pnt.dat+i);
    glEnd();
    glColor3f(0.2,0.9,0.2);
    glBegin(GL_LINE_STRIP);
    for (i=0;i<lin.num;i+=3) glVertex3dv(pnt.dat+lin.dat[i]);
    glEnd();
    }
//---------------------------------------------------------------------------

和预览:

预习

我所做的很简单:

  1. 获得表面点pnt[]

    只需测试(x,y,z)某个BBOX卷中的所有点,如果该点的方程结果接近零,我将其添加到点列表中pnt[]。由于我使用嵌套的 for 循环,其中z轴是外循环,点已经按z坐标排序,并且slice[]列表包含每个z切片的点的起始索引,因此无需在整个列表中缓慢搜索,因为我可以从z到直接切片索引。

  2. 计算参数螺旋

    使用圆柱体参数方程,我可以计算x,y,z任何t=<0,1>如果我知道半径 ( 1.0) 和螺钉数量 ( N)。由于z已经是形状的参数,我可以使用它而不是t......

  3. 计算螺旋和形状之间的交集

    简单地循环整个螺旋。对于它的每个点,找到与螺旋/形状中心轴具有相同方向并且更远的点......这可以通过方向向量之间的简单点积来完成。因此,只需z在形状的螺旋计算切片上使用实际测试点的坐标并测试其所有点,记住点积的最大值。由于中心轴是z不需要叠加的......找到的点索引被存储起来以lin[]备后用......

正如您在预览中看到的那样,螺旋有点像锯齿状。这是由于点的非线性分布。如果您在每个切片上创建点的拓扑并通过它连接所有相邻切片,则可以插入任何表面点来消除此问题。在此之后,您可以使用更少的点来获得更好的结果,并且还可以渲染面部而不仅仅是线框。

另一种选择是使用平均 FIR 滤波器来平滑螺旋点

请注意在代码中我重新调整了z轴(z' = 4*z - 2),使其与链接中的图相匹配......

[Edit1] 路径对齐螺旋

您可以创建曲线/多段线或描述螺旋/螺旋中心的任何内容。我懒得构造这样的控制点,所以我改用形状中心点(按每个切片0.5*(min+max)x,y坐标计算,但对于“裤子”是 2 个管子的部分……每切片只计算一半的形状(x<0.0)。其余的我只是用二次曲线插值......

从此只需处理螺旋/螺旋的所有中心,计算局部 TBN 矩阵(切线、法线、副法线),其中切线也是中心路径的切线,并且其中一个轴与轴对齐,Y因为没有变化(y=0)所以螺旋角将与同一轴对齐并且不会扭曲。从中只需计算螺旋角(从路径开始的中心弧长是参数)并将圆柱坐标应用于t,b,n基向量而不是直接应用于x,y,z坐标。

之后,只需将光线从中心投射到螺旋方向,并在与表面相交时渲染/存储点......

这里预览结果:

路径对齐螺旋

并更新了 C++/VCL/GL 代码:

//---------------------------------------------------------------------------
List<int> slice;        // start index of each z slice in pnt
List<double> pnt;       // GL_POINTS surface points
List<double> path;      // LINE_STRIP curved helix center points
List<double> spir;      // LINE_STRIP spiral points
//---------------------------------------------------------------------------
void obj_init()
    {
    const double ry=sqrt(1.0/3.0);
    const double rx=1.0+ry;
    double a,x,y,z,zz,d,N=12.0,da,dx,dy,dz,ex,ey,ez,r,rr,_zero;
    int i,j,i0,i1,ii;

    // [get points of analytical surface]
    pnt.num=0;
    slice.num=0;
    d=0.02; dz=0.25*d;_zero=1e-2;
    for (z=0.0;z<=1.0;z+=dz,slice.add(pnt.num))
     for (x=-rx;x<=rx;x+=d)
      for (y=-ry;y<=ry;y+=d)
        if (fabs(((1.0-z)*(((x-1.0)*(x-1.0))+y*y-(1.0/3.0))*(((x+1.0)*(x+1.0))+(y*y)-(1.0/3.0)))+(z*((x*x)+(y*y)-(1.0/3.0))))<=_zero)
            {
            pnt.add(x);
            pnt.add(y);
            pnt.add(4.0*z-2.0);
            }

    // [helix center path] as center point of half-slice
    path.num=0;
    for (i1=0,j=0;j<slice.num;j++)      // all slices
        {
        i0=i1; i1=slice.dat[j];
        dx=+9; ex=-9;
        dy=+9; ey=-9;
        dz=+9; ez=-9;
        for (i=i0;i<i1;i+=3)            // single slice
            {
            x=pnt.dat[i+0];
            y=pnt.dat[i+1];
            z=pnt.dat[i+2];
            if (x<=0.0)                 // just left side of pants
                {
                if (dx>x) dx=x; if (ex<x) ex=x; // min,max
                if (dy>y) dy=y; if (ey<y) ey=y;
                if (dz>z) dz=z; if (ez<z) ez=z;
                }
            }
        if (dz>0.25) break;             // stop before pants join
        path.add(0.5*(dx+ex));
        path.add(0.5*(dy+ey));
        path.add(0.5*(dz+ez));
        }
    // smooth by averaging
    for (j=0;j<20;j++)
     for (i=3;i<path.num;i+=3)
        {
        path.dat[i+0]=0.75*path.dat[i+0]+0.25*path.dat[i-3+0];
        path.dat[i+1]=0.75*path.dat[i+1]+0.25*path.dat[i-3+1];
        path.dat[i+2]=0.75*path.dat[i+2]+0.25*path.dat[i-3+2];
        }
    // interpolate bridge between pants from last path to (0.0,0.0,0.75)
    i=path.num-3;
    dx=path.dat[i+0];
    dy=path.dat[i+1];
    dz=path.dat[i+2];
    for (a=0.0;a<1.0;a+=d)
        {
        x=dx*(1.0-a*a);
        y=dy;
        z=dz+0.5*a;
        path.add(x);
        path.add(y);
        path.add(z);
        }
    // mirror/reverse other half
    for (i=path.num-3;i>=0;i-=3)
        {
        path.add(-path.dat[i+0]);
        path.add( path.dat[i+1]);
        path.add( path.dat[i+2]);
        }

    // [path aligned spiral line envelope]
    spir.num=0; _zero=1e-2; d=0.01;
    double *p1,*p0,t[3],b[3],n[3];      // spiral center (actual,previous), TBN (tangent,binormal,normal,tangent)
    double u[3],v[3],p[3],dp[3];
    vector_sub(p,path.dat+3,path.dat);  // mirro first path point
    vector_sub(p,path.dat,p);
    p1=p;
    for (j=0;j<path.num;j+=3)           // go through whole path
        {
        // path center previous,actual
        p0=p1;
        p1=path.dat+j;
        // TBN basis vectors of the spiral slice
        vector_sub(t,p1,p0);    vector_one(t,t);    // tangent direction to next center o npath
        vector_ld(n,0.0,1.0,0.0);
        vector_mul(b,n,t);      vector_one(b,b);    // binormal perpendicular to Y axis (path does not change in Y) and tangent
        vector_mul(n,t,b);      vector_one(n,n);    // normal perpendiculer to tangent and binormal
        // angle from j as parameter and screws N
        a=N*2.0*M_PI*double(j)/double(path.num-3);
        // dp = direction to spiral point
        vector_mul(u,n,sin(a));
        vector_mul(v,b,cos(a));
        vector_add(dp,u,v);
        vector_mul(dp,dp,d);
        // center, step
        x=p1[0]; dx=dp[0];
        y=p1[1]; dy=dp[1];
        z=p1[2]; dz=dp[2];
        // find intersection between p+t*dp and surface (ray casting)
        for (r=0.0;r<2.0;r+=d,x+=dx,y+=dy,z+=dz)
            {
            zz=(z+2.0)*0.25;
            if (fabs(((1.0-zz)*(((x-1.0)*(x-1.0))+y*y-(1.0/3.0))*(((x+1.0)*(x+1.0))+(y*y)-(1.0/3.0)))+(zz*((x*x)+(y*y)-(1.0/3.0))))<=_zero)
                {
                spir.add(x);
                spir.add(y);
                spir.add(z);
                break;
                }
            }
        }
    }
//---------------------------------------------------------------------------
void obj_draw()
    {
    int i,j;
    glColor3f(0.3,0.3,0.3);
    glBegin(GL_POINTS);
    for (i=0;i<pnt.num;i+=3) glVertex3dv(pnt.dat+i);
    glEnd();
    glLineWidth(5.0);
    glColor3f(0.0,0.0,0.9);
    glBegin(GL_LINE_STRIP);
    for (i=0;i<path.num;i+=3) glVertex3dv(path.dat+i);
    glEnd();
    glColor3f(0.9,0.0,0.0);
    glBegin(GL_LINE_STRIP);
    for (i=0;i<spir.num;i+=3) glVertex3dv(spir.dat+i);
    glEnd();
    glLineWidth(1.0);
    }
//---------------------------------------------------------------------------

我还平滑了路径并只计算了一半......因为其余的只是复制/镜像......

这种方法应该适用于任何分析表面和中心路径......

我也使用我的动态列表模板,所以:


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

一些相关的质量保证读数:


推荐阅读