How can I find intersection of a cylinder and helical isosurface.

[x3, y3,z3] = meshgrid(linspace(-1, 1),linspace(-1, 1),linspace(0, 4*pi)); f1=x3.^2+y3.^2-1; f2=y3-x3.*tan(z3); [x2, y2] = meshgrid(linspace(-1,1)); z2=atan2(y2,x2);

% Visualize the two surfaces. % patch(isosurface(x3, y3, z3, f1, 0), 'FaceColor', [0.5 1.0 0.5], 'EdgeColor', 'none'); patch(isosurface(x3, y3, z3, f2, 0), 'FaceColor', [1.0 0.5 0.0], 'EdgeColor', 'none'); view(3); camlight; axis vis3d; hold on

Df3=f1-f2;

Df3s = interp3(x3, y3, z3, Df3, x2, y2, z2);

% Find the contour where the difference (on the surface) is zero. C = contours(x2, y2, Df3s, [0 0]);

% Extract the x- and y-locations from the contour matrix C. xL = C(1, 2:end); yL = C(2, 2:end);

% Interpolate on the first surface to find z-locations for the intersection % line. zL = interp2(x2, y2, z2, xL, yL); % Visualize the line. line(xL,yL,zL,'Color','g','LineWidth',100);

Answer 1

[Mike Garrity]

In general, computing the intersection of implicit surfaces is a fairly gnarly problem. I think that the curve tracing method described here is the standard approach.

But if one of your surfaces is always going to be a cylinder, then you should be able to take some shortcuts. I don't know whether you have the Symbolic toolbox, but if you do, you might be able to get an analytic solution.

If not, here's a rough outline of how you might do it using what you already have.

The isosurface function returns a struct with two fields named vertices and faces. The faces array contains 180,146 rows and 3 columns. Each row is 3 indices into the Vertices array. Those are the three vertices at the corners of a triangle. You can take those vertices and turn them into a list of ones which are inside the cylinder and the ones that are outside the cylinder:

[x3, y3,z3] = meshgrid(linspace(-1, 1), ...
                       linspace(-1, 1), ...
                       linspace(0, 4*pi)); 
f2 = y3-x3.*tan(z3); 
hel = isosurface(x3, y3, z3, f2, 0);
r = sqrt(hel.vertices(:,1).^2 + hel.vertices(:,2).^2);
mask = r>1;

Now, we can go through all of the faces and add up the mask elements for the three vertices.

outcount = sum(mask(hel.faces),2);

Whenever this adds up to 3, then we have a triangle which is completely outside the cylinder. When it adds up to 0, we have a triangle which is completely inside the cylinder. We can ignore both of those cases. But if it adds up to either 1 or 2, then we have a triangle which is one the original surface which either has one vertex inside the cylinder and two outside, or it has two inside and one outside. Let's take a look at these triangles:

cross = (outcount == 2) | (outcount == 1);
patch('Vertices',hel.vertices,'Faces',hel.faces(cross,:), ...
      'FaceColor','red','EdgeColor','black')
view(3)

That looks kind of promising. Let's zoom in.

The next step would be to loop over those 5,569 triangles and compute out the line where they cross the cylinder. That's step is pretty straight forward, and I'll leave it for you, if that's OK. The sum of all of those line segments is the curve you want.

Update : When I first hit post, I had a mistake. I was only looking at the tris with a outcount of 2. I edited to add the ones with a outcount of 1.

[Mike Garrity]

I've expanded on this approach over on the MATLAB Graphics blog .