(Optimization) How to quickly get coordinates (coordinates are integer) of all possible points inside a specific cylinder?

Let's say I know coordinates of the center points of bottom surface and top surface; radius (R) and height (H) of cylinder. I need to write down all points which are inside the cylinder and also have integer coordinates.

The calculation speed could not be too slow since I have lots of calculation to run and need my simulation work down in rational time!

_____________________________________________________________________________

Roger Stafford has given a right method. Is any one can optimize the method for improving calculation speed? Since I need to do lots of that calculation. _____________________________________________________________________________

Answer 1

[Roger Stafford]

The hard part of your question is formulating the equations of the cylinder's bounding surfaces.

If P1 = [x1,y1,z1] and P2 = [x2,y2,z2] are the center points you referred to, and R is the cylinder's radius, then the equation of any point P = [x,y,z] on the infinite cylinder's surface is:

   dot(P-P1,P-P1)-dot(P-P1,P2-P1)^2/dot(P2-P1,P2-P1) = R^2

The equation of the infinite plane of the end containing P1 is:

   dot(P-P1,P2-P1) = 0

and for the other end

   dot(P-P2,P2-P1) = 0.

Nested for-loops can find all the points within the cylinder with integer-valued coordinates:

   in = zeros(ceil(pi*(R+1)^2*(norm(P2-P1)+2)),3);
   k = 0;
   H2 = dot(P2-P1,P2-P1);
   for x = floor(-R+min(x1,x2)):ceil(R+max(x1,x2))
    for y = floor(-R+min(y1,y2)):ceil(R+max(y1,y2))
     for z = floor(-R+min(z1,z2)):ceil(R+max(z1,z2))
      P = [x,y,z];
      if dot(P-P1,P-P1)-dot(P-P1,P2-P1)^2/H2<=R^2 & ...
         dot(P-P1,P2-P1)>=0 & dot(P-P2,P2-P1)<=0
        k = k + 1;
        in(k,:) = P;
      end
     end
    end
   end
   in(k+1:end,:) = [];

[Corrected]

[Yi]

Thank you very much. I am very appreciate about your answer. But I still have some concerns about the calculation speed.

Another difficulty of my problem is the optimization of calculation speed. Since I need to do lots of such calculation in my work and the distance between x1 and x2, y1 and y2, z1 and z2 goes to more than 1000 sometimes, the loop of 10^9 calculation require lots of time. I was wondering if there are any way for optimization ?

I can use meshgrid() and reshape() to eliminate loops (as following), but is there any way to minimize the number of points need to be calculated ? Now, the number of points need to be calculated is LX*LY*LZ; Where

LX=length(floor(-R+min(x1,x2)):ceil(R+max(x1,x2))); 
LY=length(floor(-R+min(y1,y2)):ceil(R+max(y1,y2))); 
LZ=length(floor(-R+min(z1,z2)):ceil(R+max(z1,z2)))

Can we minimize this number? If there is any way to do that, the calculation time would be decreased significantly.

BTW, in my data, distance between x1 and x2 might be big, but R is usually very small, usually R=0.5 ;

%%%%%%%%%%%%%%%%%%%%%%%%%

[X,Y,Z]=meshgrid(floor(-R+min(x1,x2)):ceil(R+max(x1,x2)), floor(-R+min(y1,y2)):ceil(R+max(y1,y2)), floor(-R+min(z1,z2)):ceil(R+max(z1,z2)));

LX=length(floor(-R+min(x1,x2)):ceil(R+max(x1,x2)));

LY=length(floor(-R+min(y1,y2)):ceil(R+max(y1,y2)));

LZ=length(floor(-R+min(z1,z2)):ceil(R+max(z1,z2)));

NUM=LX*LY*LZ;

P=[reshape(X,NUM,1),reshape(Y,NUM,1),reshape(Z,NUM,1)];

P1=[ones(NUM,1)*P1(1,1) , ones(NUM,1)*P1(1,2) , ones(NUM,1)*P1(1,3)];

P2=[ones(NUM,1)*P2(1,1) , ones(NUM,1)*P2(1,2) , ones(NUM,1)*P2(1,3)];

P_extract=(dot(P-P1,P-P1,2)-dot(P-P1,P2-P1,2).^2/H2<=R^2) .* (dot(P-P1,P2-P1,2)>=0) .* (dot(P-P2,P2-P1,2)<=0) .* (sum(P'))';

index=find(P_extract>0); %%this might miss the origin point...

P=P(index,:);

%%%%%%%%%%%%%%%%%%%%%%%%%

My first version code actually failed because of too much calculation time. (Matlab code is much easy but takes too much calculation time than C++ or C)

Your code is much better than mine. I am gonna accepted this if no one can provide optimization~

=_= ^_^ ^v^ I really need to optimize... >_<

Again, thank you so much. Very appreciated!

(Corrected)

[yogesh jain]

Thank you very much Yi , this code helped me very much for different purpose . :)

[yogesh jain]

If I want to modify this for cone then what parameters should I change ?

Answer 2

[Image Analyst]

Why are the (x,y) coordinates at some z level any different than those at the end caps? Why can't you just use repmat() to copy them?

[Yi]

Because the cylinder could be oblique and end cap of the cylinder is not always parallel to x-y plane, you can not just copy it. And calculation for every z level is needed.

[Image Analyst]

For corresponding points on the upper cap and lower cap, you can simply use the point slope formula of a line to get points at any location along a line connecting them.

[Yi]

yeah, that's a great idea. But do I know How many lines I need for get all the points with integer points. More specifically, which line has points with integer coordinates on it ?

[Image Analyst]

Possibly the easiest way is to just rotate the cap so that it is parallel to the x-y plane (it's all in 1 z level), then use repmat(). Then rotate back.

[Yi]

Yeah, that would definitely be much easier for calculation. But how could we know which point is gonna have integer coordinates after rotating back? since in the rotation, a point with integer coordinates might no longer have integer coordinates. My point is that if we are not calculating points with integer coordinates any more, how do we know all possible points are calculated and which point should be calculate?