How to speed up the process of determining if a point is inside a cylinder

I would like to speed up the process of determining if a point is inside a cylinder. Specifically:

A 3D cubic mesh is used to divide a problem. In 3 separate 3D arrays the cell center coordinates are stored (i.e cx, cy, cz). A cylinder is upright either in the x, y, or z direction. Then based on the cap center coordinates, the height of the cylinder, and the radius of the cylinder, the following code is used to determine if a point is inside the cylinder:

% nx, ny, nz is the number of cells in cx, cy or cz.
% range = distance between the cap centers of the cylinder
% For all points in the mesh
for lind = 1:nx
    for jnd = 1:ny
        for knd = 1:nz
            % Get cell center coordinates (x, y, z)
            apoint = [cx(lind,jnd,knd) cy(lind,jnd,knd) cz(lind,jnd,knd)];
            % Project point on center line
            tmin(lind, jnd, knd) = (dot(apoint, range) - dot(cap1, range))/ ...
            dot(range, range);
            linept = cap1 + (tmin(lind, jnd, knd) * range);
            % Distance form center line
            diff = linept - apoint;
            dist(lind, jnd, knd) = sqrt(dot(diff, diff));
        end
    end
end

% Get all the indices of the points that are part of the cylinder
l = find((((tmin>=0) & (tmin<=1)) & (dist<=radius)));

How can I speed up the computation of determining whether a point is part of a cylinder specifically avoid using for-loops?

Any help would be much appreciated!

Answer 1

[Doug Hull]

For loops are not always a problem. Run this through the profiler to find where the slow down is. I would do something like this.

pointXVec, pointYVec, pointZVec are locations TopZ, BotZ, are the locations of top and bottom CenterX, CenterY, are location in X and Y Radius

These can be switched around for other orientations.

flagBelowVec = (pointZVec <= topZ);
flagAboveVec = (pointZVec >= botZ);
radialDistanceSquaredVec = (pointXVec-centerX)^2 + (pointYVec-centerY)^2;
flagInsideVec = (radialDistanceSquaredVec <= radius^2);
flagIsIn = (flagBelowVec & flagAboveVec & flagInsideVec);

[Christos]

Thank you Doug! Absolutely perfect! The profiler showed that my previous code called the dot product function 6524928, for a 96x96x177 mesh, taking up to 263s. Your method ... 0.108s!!!! Thank you again!!!

[sam l]

Hi Doug Hull

I amusing the above code you wrote to get index of of points inside Cylinder.

However, my case is having a cylinder atinclined angle. say 18 degree, with cylinder start point x0,y0,z0

How to change the index syntax for that inclination

[Bruno Luong]

@sam: if you open a new thread, I'll reply to you