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]
