在曲面上绘制螺旋曲线

Question

我对在表面上绘制螺旋曲线很感兴趣，它看起来像这样：

它的等式是：

((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中看到post但是我不明白答案...：

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

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

Answer 1

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

为简单起见，我假设为螺旋形包络线（就像用绳子缠绕形状一样）。

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

另一种选择是将您的形状切割成单独的切片并将每个切片处理为 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();
    }
//---------------------------------------------------------------------------

并预览：

我做的很简单：

获取曲面点数pnt[]

只需测试一些 BBOX 体积中的所有 (x,y,z) 点，如果该点的方程结果接近于零，我将其添加到点列表 pnt[].由于我使用 z 轴作为外循环的嵌套 for 循环，点已经按 z 坐标排序，并且 slice[] 列表包含每个点的起始索引 z slice 所以不需要在整个列表中缓慢搜索，因为我可以直接从 z 到切片索引。
计算参数螺旋

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

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

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

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

请注意我重新调整了 z 轴 (z' = 4*z - 2) 的代码，使其与您 link ...

中的绘图相匹配

[Edit1] 路径对齐螺旋

您可以创建一个 curve/polyline 或任何描述 spiral/helix 中心的内容。我懒得构造这样的控制点，所以我使用了形状中心点（计算为每个切片的 x,y 坐标的 0.5*(min+max)，但只是形状的一半（x<0.0） “裤子”的部分是 2 管......每片。其余部分我只是用二次曲线插值......

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

之后，只需将光线从中心投射到螺旋方向，当与表面相交时 render/store 点 ...

此处预览结果：

并更新了 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);
    }
//---------------------------------------------------------------------------

我还平滑了路径并只计算了一半...因为其余部分只是 copy/mirror ...

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

我也使用我的动态列表模板：

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

一些相关的 QA 读物：

Smoothly connecting circle centers 用于类似目的的基向量
使用矢量数学

在曲面上绘制螺旋曲线

Draw spiral curves on a surface

3d

geometry

computational-geometry

geometry-surface