MATLAB Answers

Step Function for Checking if Point is inside Spherocylinder

1 view (last 30 days)
Mari Teli
Mari Teli on 15 Apr 2021
Answered: William Rose on 15 Apr 2021
Hello! I'm trying to generate a step function that will check to see if a point is inside the spherocylinder cell. After that is generated then I need to use normrnd to pick a random point within the cell. However, I'm stuck on this, therefore any help would be appreciated.
This is what I have so far:
function [X,Y,Z]=spherocylinder(pos, vec, totalLength, AR);
%% preliminary
points = 20;
vec = vec/norm(vec);
D = (0.65);
radius = D/2;
bottom = -0.675; % where bottom of cylinder is
top = 0.675; % where top of cylinder is
%% Generate cylinder with center at origin
[xcyl,ycyl,zcyl]=cylinder(radius,points);
zcyl(1,:) = bottom;
zcyl(2,:) = top;
%% Generate points for sphere
[x,y,z] = sphere(points);
pdist = points/2+1;
% split sphere into top and bottom hemispheres
xtop = x(pdist:end,:)*radius;
ytop = y(pdist:end,:)*radius;
ztop = z(pdist:end,:)*radius+top;
xbot = x(1:pdist,:)*radius;
ybot = y(1:pdist,:)*radius;
zbot = z(1:pdist,:)*radius+bottom;
%% Combine data for spherocylinder oriented along z-direction
X = [xbot; xcyl; xtop];
Y = [ybot; ycyl; ytop];
Z = [zbot; zcyl; ztop];
%% Plot spherocylinder
surf(X, Y, Z, 'EdgeColor','none',...
'LineStyle','none','FaceColor',[0.85 0.85 0.85]);
axis equal;
% length, width and height of box
L = 0.65; W = 0.65; H = 2;
% pick plenty of random points
x = L*rand(1,1); %x-coordinate of a point
y = W*rand(1,1); %y-coordinate of a point
z = H*rand(1,1); %z-coordinate of a point
scatter3(x, y, z, 'MarkerFaceColor','b','MarkerEdgeColor','b');
box
end
Thank you for any help!
  1 Comment
DGM
DGM on 15 Apr 2021
I would consider the line segment connecting the centers of the two hemispheres. All points within the object should be within distance=radius from that line segment.
As an aside, it's good to include some form of synopsis in your function to describe what the input (and output) arguments are.

Sign in to comment.

Answers (1)

William Rose
William Rose on 15 Apr 2021
@Mari Teli, I do not understand what spherocylinder() is supposed to do, because your question does not seem to match the output from the function. I ran the function. It produced a 3D plot of what appears to be a single point.
That plot was not enlightening, but the function also returned values of arrays x,y,z. I plotted those arrays with "surf(x,y,z)" and got a 3D spherocylinder, as shown on left below. To get the correctly scaled object, as shown on the right, in which the ends look like hemispheres, add the command "axis equal". You could add these to your function, if you want such a plot.
>> clear;
>> [x,y,z]=spherocylinder(1,1,1,1);
>> figure;
>> surf(x,y,z)
>> axis equal
becomes
Three of the arguments to the function are not used at all, an the other argument, vec, is scaled, but it is otherwise not used, so it does not affect the plot made or the output variables. What did you want those four function arguments to do?
You said in your quesiton "I'm trying to generate a step function that will check to see if a point is inside the spherocylinder cell. After that is generated then I need to use normrnd to pick a random point within the cell." That seems backwards. It would make more sense to pick a random point first, then test it to see it it is inside the cell.
The parameters of the cyliner are its radius, top height, and bottom height. Here is a function that retrns 1 if point is inside and 0 if not. I tested it for many combinations and it worked, but you should do more testing, to be sure.
function b = inSpherocylinder(p,radius,th,bh)
% WCRose 20210415
% Return 1 if p is in spherocylinder, else return 0
% Cylinder's axis is assumed to be along z and centered on x,y=0,0.
% th=level where cylinder meets top hemisphere
% bh=level where cylinder meets bottom hemisphere
% p should be a 3-vector: p=[x,y,z].
% th should be greater than or equal to bh.
if length(p)~=3
b=-99;
disp('Error in inCylinder: p should be a 3-element vector.');
return;
elseif bh>th
b=-99;
disp('Error in inCylinder: th should be >= bh.');
end
x=p(1); y=p(2); z=p(3);
if x^2+y^2>radius^2
b=0; %point is outside cylinder or its extension
return;
elseif (z>th) & (x^2+y^2+(z-th)^2>radius^2)
b=0; %point is above the top hemi
return;
elseif (z<bh) & (x^2+y^2+(z-bh)^2>radius^2)
b=0; %point is below the bottom hemi
return;
else
b=1; %point is inside
end
end

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!