Interpolate Helix in 3d

I have three helices with 39 data points for each helix:

t=linspace(0,8.6*pi,39);

plot3(35*cos(t)+35,35*sin(t)+35,t,'LineWidth',3);

hold on;

plot3(25*cos(t)+35,25*sin(t)+35,t,'LineWidth',3);

plot3(15*cos(t)+35,15*sin(t)+35,t,'LineWidth',3);

hold off;

I am wondering how I could interpolate the area between the helices in the x,y,z direction? Thanks!

Answer 1

[John D'Errico]

If you want to interpolate the helices smoothly, download interparc from the file exchange.

t=linspace(0,8.6*pi,39);
x = 35*cos(t)+35;
y = 35*sin(t)+35;

xyzi = interparc(1000,x,y,t,'spline');
plot3(x,y,t,'ro')
hold on
plot3(xyzi(:,1),xyzi(:,2),xyzi(:,3),'-b')
box on
grid on

But re-reading your question, you asked to interpolate the area BETWEEN the helices.

Are you asking how to make the three curves into a plottable surface, connecting them together? This is quite doable, but I won't do so without affirmation from you as to my guess.

[Robert Wu]

Sorry for the confusion. I would like to interpolate vertically between data points in a helix as a start - so the interpolation would pretty much look like a cylinder. I would then like to interpolate horizontally/diagonally between the data points in all 3 helices.

[Robert Wu]

I would like something like this, where the the coloured area will be the space interpolated from the 39 red data points:

Please note the cylinder was generated by means that will not help me in this problem.

[John D'Errico]

Good, in one sense. At least I did not waste time on a mis-guess. :)

So you have some general helix, not a nice, easily defined, circular, cylindrical helix. It is defined only by a set of ordered triplets in 3-d. And now you want to interpolate up along the axis of the helix. I suppose, just to make things even worse, you won't even be willing to assume the axis is parallel to the z-axis?

The point is, if you were willing to parameterize the domain in cylindrical coordinates, then you could write the helix as an ordered set of triplets (r,theta,z). Now, your goal would be to create the surface r(theta,z). This would be doable. If it turned out that the cylinder axis were not parallel to the z-axis, then I would find the axis, and do a rotation of the coordinate space so that is now true. At least that part is totally simple.

So, I might start like this:

z = t;
xbar = mean(x);
ybar = mean(y);
theta = atan2(y-ybar,x-xbar);
r = sqrt((x-xbar).^2 + (y-ybar).^2);
plot3(theta,z,r,'o')
grid on
box on
view(-7,18)

I'm not sure you can see that as obviously a surface, if you rotate the figure around a bit, you will see it as such.

[rgrid,thetagrid,zgrid] = gridfit(theta,z,r,linspace(-pi,pi,100),linspace(min(z),max(z),100));

figure
plot3(theta,z,r,'ro')
hold on
surf(thetagrid,zgrid,rgrid)
view(-7,18)
box on
grid on

xlabel theta
ylabel z
zlabel r

Ok, this is better. We have a surface, that in theory could be converted back into the original coordinate system. But for one small problem, that when I built the gridded surface, I was unable to explicitly put in the information that it was periodic in theta.

So you can see there is some crap happening at the boundary where theta transitions from -pi to pi.

Gridfit is unable to handle a periodic boundary condition like that, and as you can see, the edge in that region got a bit ragged, since gridfit thinks it is forced to extrapolate.

Lets recover the original coordinate system.

xgrid = xbar + rgrid.*cos(thetagrid);
ygrid = ybar + rgrid.*sin(thetagrid);
surf(xgrid,ygrid,zgrid,'edgecolor','none')
hold on
plot3(x,y,z,'ro')
view(-171,78)
box on
grid on

It is pretty good, except that the periodic boundary fails periodicity. I've rotated the view to show that problem.

If that were not a problem, then I would tell you to simply download gridfit from the File Exchange and use it as is.

It is a start at least. :) That gives you a surface.

Then your next comment was you wanted to next interpolate between these surfaces. That would create what is essentially a 3-d volume, a cylinder with a thick wall. I'm not sure what you mean by interpolating there.

[Robert Wu]

Thanks John!

The z axis is parallel to the axis of the helix. The space at the pi boundary may become a problem in the future.

I downloaded gridfit and tried to use my 39 data points

48 48 48 49 49 48 48 48 49 50 49 50 49 45 36 10 11 19 20 19 19 19 19 21 22 22 22 21 21 22 22 23 23 23 23 23 23 24 24

with the code you have provided above. It did not not work out too well; I think there was complications when using gridfit.

[John D'Errico]

I don't know what you have shown me there. You need to be careful in what variables you treat as the dependent and independent variables. Perhaps you might attach the actual data that you have.

The -pi/pi boundary will be an issue here, as gridfit would need to be modified to solve your problem. It would almost be easier to write a hacked (and very simple version) that did have periodic boundary conditions in one variable.

Anyway, there are other solutions one might look at. For example, given all three curves, generate a 3-d delaunay tessellation. Of course, that will include the volume internal to the central helix. Now, delete any tetrahedra that cross that internal region. The result will be a tessellation of the annular volume in question. The only issue would be how to delete the stuff in the middle. That would be doable though, and it completely eliminates the periodic boundary problem.

[Robert Wu]

Thanks for the info! That row of numbers in the previous comment is the data points that winds down a cylinder in a helical fashion. I want to make something like this 2D plot, but on 3D cylinder:

This is a cross-section view of the cylinder. The 3 sine waves are the helices in 2D. The dots on the wave is where the data points are. The area is interpolated by what the value of the dots are.

What I am having trouble with is representing this in 3D.