MATLAB: Intersection of 3D ray and surface of revolution

3dintersectionMATLABraysurface of revolution

Greetings. I'm looking for an efficient routine to find the 3D intersection point of a ray and a surface of revolution. Assume we have the generic 3D vase of height H, generated by revolving a curve of pre-stored data points. Sit the vase on the x-y plane with the center of the base over the origin. Now shoot a 3D ray from the origin in an az-el direction to an altitude of H, long enough to travel through the vase surface. The shape of the vase is such that the ray will only hit one point on the surface, probably between the predefined data points. Can the community weigh in on the most efficient way to find the intersection point? I have several brute-force solutions that look downright ugly but seem to produce reasonable results. Thanks.

Best Answer

I'll assume, without loss of generality, that the CCW angle is zero which means we don't have to work with the 3D surface, only its pre-revolution 1D profile. Since the geometry is circularly symmetric, you can just rotate the solution by the actual, nonzero CCW angle at the end. Let the revolution profile be defined by (r,z) coordinates where vector R are the distance from the origin of your samples and vector Z are their z-coordinates. Let d be the length of your line projected onto the XY plane. Then, the solution Rmin,Zmin could be obtained by,

       pp=interp1(Z,R,'linear','pp');
      [~,idx]=min(abs(R-Z*d/H));
      Zmin = Z(idx);
      Zmin =  fminsearch(@(z) abs( ppval(pp,z) - z*d/H ),  Zmin);
      Rmin=ppval(pp,Zmin);

Related Solutions

MATLAB: How to find the intersection between a line and a 3D surface

The "surface" exists only in your mind. All you have are points. Isolated points. So until you find away to connect those points into a true surface in some form, you can do nothing.

Simple is to use a convex hull. That is great, if the volume described is truly convex, or even sufficiently close to being convex. Your next choice is to use an alpha shape.

In either of the above cases, the outer boundary is a triangulation. So now the problem reduces to solving for the intersection of a line with a triangulated surface. This part is pretty easy. You just need to solve for the intersection of a line with a triangle, and then repeat for each triangle. Better code will see how to eliminate some of those triangles from consideration, but in any case it is not that hard. Just computational geometry.

MATLAB: Generate a cone with a point and angle

This question initially has absolutely nothing to do with MATLAB, because you have no idea how to solve it. Until you understand the mathematics of a problem, writing code is never possible. But I have a few minutes, so I might as well go through the thinking for you.

When you have a problem like this, break it down. If necessary, think in terms of other coordinate systems. Here, cylindrical coordinates would make sense, since a cone is just a rotationally symmetric surface, rotated around the z-axis.

So you want to create a function, z(x,y) that defines the cone. Since we will start out by working in cylindrical coordinates, I'll create things as z(r,theta). Then the conversion back to cartesian coordinates is trivial.

Assume the vertex of the cone will be at some value z0. And you want the angle at the point to be phi. Since the problem is rotationally symmetric cone, we can ignore theta completely. Just assume we have a line in the plane (r,z) plane, that intersects the z axis at z==z0. You want the angle that the line crosses the z axis to be phi/2, with the line going down for increasing r.

So the result in the plane (r,z) will be two lines , intersecting the z axis. The total angle between the lines will be phi, but the slope of either line will be -cotd(phi/2). (Think about it. I won't do all your thinking for you. Basic trig here. Consider the angle the line makes with the z-axis.) I've used cotd so we can use degrees for the angle.

So we can represent the cone in the (r,z) plane as

z = z0 - abs( r )*cotd(phi/2)

Now, convert to cartesian coordinates.

We should know that r=the equation of a circle of radius r is simply

r^2= x^2 + y^2

So then

abs( r ) = sqrt(x^2 + y^2)

thus the positive branch of the sqrt. This gives you the surface equation as desired.

z(x,y) = z0 - sqrt(x.^2 + y.^2)*cotd(phi/2)

So a rotationally symmetric surface arpund the z-axis. The radius of the circle increases linearly as z decreases. The peak of the cone is at (x,y,z) = (0,0,z0), again as desired. The angle of the cone is clearly phi.

Note that I used .^ in that equation to make sure you understand to use that operator. This will allow vectorized computations to be done.

For example...

conefun = @(x,y,z0,phi) z0 - sqrt(x.^2 + y.^2)*cotd(phi/2);
fsurf(@(x,y) conefun(x,y,10,90))
axis equal

So a conical surface, here passing through z0 = 10, with an angle of 90 degrees. You can see the surface is clearly rotationally symmetric, and we see the conic surface decreases linearly with distance from the axis.

To convince you the surface is rotationally symmetric around the z axis, look at contours of the surface.

fcontour(@(x,y) conefun(x,y,10,90))
axis equal
grid on

So perfectly circular contours, centered around the origin, (x,y) = (0,0). A cone. The conic surface extends to -inf as (x,y) go away from the origin.

Best Answer

Related Solutions

MATLAB: How to find the intersection between a line and a 3D surface

MATLAB: Generate a cone with a point and angle

Related Question