Lets start with axis aligned cylinder. For such I see it like this:
definition
let XY plane be the base and cylinder starts at (0,0,0) and grows into +Z direction up to distance l and have radius r. Also let me define l0,l1 be the min and max distance between nodes.
create main path
simply put a chain of connected nodes comming from start of cylinder ot its end. These will be used later to grow the clusters. This also ensures that path from start to end exist. So just add some randomized increments to z in range <l0,~l1> and use x,y as a randomized point inside some smaller circle than r to stay on path (I use 10% of r).
grow clusters
simply take randomly any already placed point add a random displacement to it of size <l0,l1> and if still inside cylinder and not too close to any other point add it to your data and link it to the picked point. This can be speed up if you use points sorted by z so you can get rid of the O(n) search and use O(log(n)) instead.
After this you just convert the axis aligned data into the final position and orientation you want. For example if you define your cylinder as 2 endpoints and radius then you can compute l as their distance, l0,l1 as its fraction. Also you can compute 3 perpendicular basis vectors from it using simple vector math (2 representing XY plane and one representing cylinder axis Z) lets call them u,v,w. From that its just matter of vector math to convert... You can also construct 4x4 transform matrix from those and use that.
Here small C++ example for this:
//---------------------------------------------------------------------------
const int N=200; // points to generate
double pnt[N][3]; // random point
int lnk[N]; // -1 or pnt[i] is linked to pnt[lnk[i]]
//---------------------------------------------------------------------------
void vector_mul(double *c,double *a,double *b) // c[3] = cross(a[3],b[3])
{
double q[3];
q[0]=(a[1]*b[2])-(a[2]*b[1]);
q[1]=(a[2]*b[0])-(a[0]*b[2]);
q[2]=(a[0]*b[1])-(a[1]*b[0]);
for(int i=0;i<3;i++) c[i]=q[i];
}
//---------------------------------------------------------------------------
void generate(double *p0,double *p1,double r)
{
int i,j,k,ok;
double u[3],v[3],w[3]; // basis vectors
double a,dx,dy,dz,x,y,z,z0;
double l; // cylinder size |p1-p0|
double l0=0.03; // min distance between major nodes <0,1>
double l1=0.06; // max distance between major nodes <0,1>
double ll0=l0*l0,ll1=l1*l1,rr=r*r;
Randomize();
// basis vectors from endpoints
for (l=0.0,i=0;i<3;i++){ w[i]=p1[i]-p0[i]; l+=w[i]*w[i]; } // w = (p1-p0)
l=sqrt(l); l0*=l; l1*=l; // l=|w| , convert l0,l1 to units
for (i=0;i<3;i++) w[i]/=l; // w/=|w|
if (fabs(w[0])<0.75){ u[0]=1.0; u[1]=0.0; u[2]=0.0; } // u=(1,0,0) or (0,1,0) so it is not paralel to w
else { u[0]=0.0; u[1]=1.0; u[2]=0.0; }
vector_mul(v,u,w); // v = cross(u,w)
// [axis aligne d cylindric data]
// random major path
for (z0=0,i=0;i<N;)
{
x=2.0*r*Random()-r; x*=0.1; // use only 10% of x,y deviation to not sray too much
y=2.0*r*Random()-r; y*=0.1;
z=z0+l0+(0.75*(l1-l0)*Random()); if (z>l) break;
// inside cylinder ?
if ((z<0)||(z>l)) continue;
if ((x*x)+(y*y)>rr) continue;
// no point closer than l0 ?
for (ok=1,j=0;j<i;j++)
{
dx=pnt[j][0]-x;
dy=pnt[j][1]-y;
dz=pnt[j][2]-z;
if ((dx*dx)+(dy*dy)+(dz*dz)<ll0){ ok=0; break; }
}
if (!ok) continue;
// add if valid point
pnt[i][0]=x;
pnt[i][1]=y;
pnt[i][2]=z; lnk[i]=i-1; i++; z0=z;
}
// grow clusters
for (;i<N;)
{
// random 3D displacement <l0,l1>
for (;;)
{
dx=Random()-0.5;
dy=Random()-0.5;
dz=Random()-0.5;
a=(dx*dx)+(dy*dy)+(dz*dz);
if (a>1e-3) break;
}
a=(l0+((l1-l0)*Random()))/sqrt(a); dx*=a; dy*=a; dz*=a;
// convert to position
for (k=0;k<10;k++)
{
// add it to random point already placed
j=Random(i); lnk[i]=j; ok=1;
x=pnt[j][0]+dx;
y=pnt[j][1]+dy;
z=pnt[j][2]+dz;
// inside cylinder ?
if ((z<0)||(z>l)){ ok=0; break; }
if ((x*x)+(y*y)>rr){ ok=0; break; }
// no point closer than l0 ?
for (j=0;j<i;j++)
{
dx=pnt[j][0]-x;
dy=pnt[j][1]-y;
dz=pnt[j][2]-z;
if ((dx*dx)+(dy*dy)+(dz*dz)<ll0){ ok=0; break; }
}
if (ok) break; // valid point
}
if (!ok) continue;
// add if valid point
pnt[i][0]=x;
pnt[i][1]=y;
pnt[i][2]=z; i++;
}
// [convert to final position and orientation]
for (i=0;i<N;i++)
{
x=pnt[i][0];
y=pnt[i][1];
z=pnt[i][2];
for (j=0;j<3;j++) pnt[i][j]=p0[j]+(x*u[j])+(y*v[j])+(z*w[j]);
}
}
//---------------------------------------------------------------------------
and usage:
double p0[3]={-1.7,-0.5,-0.2};
double p1[3]={+1.7,+0.5,+0.4};
generate(p0,p1,0.5);
preview:
beware if you set too big N the second main loop can loop forever. So you might want to add some ending condition like if continue was hit more than 2*i times without i change stop. That is because l0 constraint limits the max density of points and if N is bigger than that you can not add more points ...
Now if you want random radius spheres instead of the points than just add some random radius but do not forget to adjust the inside cylinder and closest distance tests by the radius ...
