Draw spiral curves on a surface
I am interesting in drawing spiral curves on a surface which look like that:
The equation of it is:
((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
where
x=<-1.0-sqrt(1/3),+1.0+sqrt(1/3)>
y=<-sqrt(1/3),+sqrt(1/3)>
z=< 0.0,1.0>
I found the equation on this website:
Then I saw that post in StackMath but I don't understand the answer...:
If someone have any idea of a process to draw spiral curves on that kind of surface, it would greatly help me. Do not hesitate to ask me more about the problem. Thank you.
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 listpnt[]. Thanks to the fact I use nested for loops withzaxis being the outer loop the points are already sorted byzcoordinate and theslice[]list holds the start indexes of points for eachzslice so no need to slow searching in the whole list as I can go fromzto slice index directly.compute parametric spiral
using cylinder parametric equation I can compute
x,y,zfor anyt=<0,1>If I know radius (1.0) and number of screws (N). Aszis already the parameter of the shape I can use it instead oft...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
zcoordinate 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 iszno superpositions are needed ... The found point indexes are stored inlin[]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:
- Smoothly connecting circle centers basis vectors used for similar purpose
- Understanding 4x4 homogenous transform matrices vector math used