analytical approach is not the only way how to attack this (an not suited for this site anyway) so there are quite a few other options. But you did not specify how the spiral should look like as there are 3 ends of the shape and spiral have just 2...
For simplicity I assuming spiral envelope (like you wrap a rope around the shape).
So how to attack this. The usual approach is to create mapping 2D <--> 3D between some rectangular 2D area (texture) and your surface. Usually exploiting topology of the shape like cylindircal/spherical coordinates etc.
Another option is to cut your shape into individual slices and handle each slice as 2D ... Which simplifies things a lot.
To apply a spiral onto some curved surface you either project it or find intersection between shape and direction to actual point on spiral defined by parametric equation.
Here C++ example using OpenGL visualization:
//---------------------------------------------------------------------------
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();
}
//---------------------------------------------------------------------------
and preview:
What I did is simple:
obtain surface points pnt[]
simply test all (x,y,z) point in some BBOX volume and if your equation result for the point is close to zero I add it to the point list pnt[]. Thanks to the fact I use nested for loops with z axis being the outer loop the points are already sorted by z coordinate and the slice[] list holds the start indexes of points for each z slice so no need to slow searching in the whole list as I can go from z to slice index directly.
compute parametric spiral
using cylinder parametric equation I can compute x,y,z for any t=<0,1> If I know radius (1.0) and number of screws (N). As z is already the parameter of the shape I can use it instead of t ...
computing intersection between spiral and shape
simply loop through whole spiral. For each of its points find point that has the same direction from spiral/shape central axis and is further away... This can be done by simple dot product between the direction vectors. So simply using z coordinate of the actual tested point on spiral compute slice of the shape and test all its points remembering the max of the dot product. As the central axis is z no superpositions are needed ... The found point indexes are stored in lin[] for later use...
As you can see on the preview the spiral is a bit jagged like. This is due to non linear distribution of points. If you create a topology of the points on each slice and connect by it all the neighbor slices you can then interpolate any surface point eliminating this problem. After this you can use much less points with still better result and also render the faces instead of just wire-frame.
Another option is just to smooth the points of spiral using averaging FIR filter
Beware in the code I rescaled the z axis (z' = 4*z - 2) so it matches the plot in your link ...
[Edit1] path aligned spiral
You can create a curve/polyline or whatever describing the center of the spiral/helix. I was too lazy to construct such control points so I used the shape center points instead (computed as 0.5*(min+max) of x,y coordinate per each slice but just for a half of shape (x<0.0) for parts where the "pants" are 2 tubes... per slice. The rest I just interpolated by quadratic curve ...
From this just process all the centers of the spiral/helix, compute local TBN matrix (tangent,normal,binormal) where tangent is also tangent of the centers path and one of the axis is aligned to Y axis as there is no change (y=0) so the spiral angle will be aligned to the same axis and not twisting. From that just compute spiral angle (the center arclength from path start is the parameter) and apply the cylindrical coordinates on the t,b,n basis vectors instead of directly on x,y,z coordinates.
After that just cast ray from center into the helix direction and when hit intersection with the surface render/store the point ...
Here preview of the result:
And updated C++/VCL/GL code:
//---------------------------------------------------------------------------
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);
}
//---------------------------------------------------------------------------
I also smoothed the path and compute just half of it ... as the rest is just copy/mirror ...
This approach should work for any analytical surface and center path ...
I also use mine dynamic list template so:
List<double> xxx; is the same as double xxx[];
xxx.add(5); adds 5 to end of the list
xxx[7] access array element (safe)
xxx.dat[7] access array element (unsafe but fast direct access)
xxx.num is the actual used size of the array
xxx.reset() clears the array and set xxx.num=0
xxx.allocate(100) preallocate space for 100 items
Some related QA readings:
