Menu
Is there an algorithm for finding of an enclosing cylinder with the smallest radius for a 3D cloud of dots? I know that 2D case with the smallest enclosing circle is solved (for example this thread Smallest enclosing circle in Python, error in the code), but is there any working approach for 3D?
EDIT1: OBB. Below is an example of an arc-shaped cloud of dots. The smallest enclosing circle was found by this tool https://www.nayuki.io/page/smallest-enclosing-circle
Circle is defined by three dots of which two are lying almost on a diameter, so it is easy to estimate where the central axis is. "Boxing" of dots will yield a center of the box obviously much shifted from the true center.
I conclude, that OBB approach is not general.
EDIT2: PCA. Below is an example of PCA analysis of a tight dot cloud vs. dot cloud with outliers. For the tight dot cloud PCA predicts the cylinder direction satisfactorily. But if there is a small number of outliers, compared to the main cloud, than PCA will basically ignore them, yielding vectors which are very far from the true axis of an enclosing cylinder. In the example below the true geometrical axis of an enclosing cylinder is shown in black.
I conclude that PCA approach is not general.
EDIT3: OBB vs. PCA and OLS. A major difference - OBB relies only on a geometrical shape, while PCA and OLS are dependent from the overall number of points, including those in the middle of the set, which do not affect the shape. In order to make them more efficient, a data preparation step can be included. First, find the convex hull. Second, exclude all internal points. Then, points along the hull can be distributed unevenly. I'd suggest to remove all of them, leaving only the polygonal hull body, and cover it with mesh, where nodes will be new points. Application of PCA or OLS to this new cloud of points should provide much more accurate estimation of the cylinder axis.
All this can be unnecessary, if OBB provides an axis, as much parallel to the enclosing cylinder axis, as possible.
EDIT4: published approaches. @meowgoesthedog: paper by Michel Petitjean ("About the Algebraic Solutions of Smallest Enclosing Cylinders Problems") could help, but I'm insufficiently qualified to convert it to a working program. Author himself did it (module CYL here http://petitjeanmichel.free.fr/itoweb.petitjean.freeware.html). But at the conclusions in the paper he says: "and the present software, named CYL, downloadable for free at http://petitjeanmichel.free.fr/itoweb.petitjean.freeware.html, is neither claimed to offer the best possible implementations of the methods nor is claimed to work better than other cylinder computation softwares." Other phrases from the paper also makes an impression, that it is an experimental approach, which was not thoroughly validated. I'll try to use it, anyway.
@Ripi2: this paper by Timothy M. Chan is also a bit too complicated for me. I'm not an expert of that level in mathematics, to be able to convert to a tool.
@Helium_1s2: probably, it is a good suggestion, however, it is much less detailed compared to two papers above. Also, not validated.
EDIT5: reply for user1717828. Two most distant points vs. cylinder axis. A counter example - 8 points in a shape of cube, fit in a cylinder. The biggest distance between two points - green diagonal. Obviously not parallel to cylinder axis.
"Middle-points" approach by Ripi2: it works only in 2D. In a 3D case the cylinder axis may not intersect a single segment between any two points.
609
Find the two points with the greatest distance h between them. They are on the faces of the cylinder and the line connecting them will be parallel to the axis of the cylinder.
Project all points on the plane perpendicular to that axis.
Find the two points with the greatest distance d between them on that plane. They define the circle whose diameter d is that of the cylinder.
The cylinder with the smallest possible volume* containing all points has
.
* This assumes there is only one pair of points with the greatest distance between them defining the axis of the cylinder. If there's a chance two pairs of points share the greatest value, repeat steps 2-4 for each pair and choose the cylinder with the smallest diameter.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook
from numpy.linalg import norm
from scipy.spatial.distance import pdist, squareform
Generate points if you don't have them already:
np.random.seed(0)
N = 30
M = np.random.randint(-3,3,(N,3))
print(M)
[[ 1 2 -3]
[ 0 0 0]
[-2 0 2]
[-1 1 -3]
[-3 1 -1]
...
[ 1 -3 1]
[ 0 -1 2]]
Calculate the distance between every possible pair of points and select the pair with the greatest distance.
max_dist_pair = list(pd.DataFrame(squareform(pdist(M))).stack().idxmax())
p1 = M[max_dist_pair[0]]
p2 = M[max_dist_pair[1]]
print(f"Points defining cylinder faces: {p1}, {p2}")
print(f"Length of cylinder: {norm(p1-p2)}")
Points defining cylinder faces: [-1 -3 -3], [1 2 2]
Length of cylinder: 7.3484692283495345
Graph the points showing all the points in blue, with the maximally separated points in red.
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(*M.T, c = ['red' if i in max_dist_pair else 'blue' for i in range(N)])
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
plt.show()
Here is the same plot rotated so we're looking along the axis between the two red dots.
The above view is the same thing as the points being projected on the plane perpendicular to the axis of the cylinder. Find the smallest circle that contains points in this plane. We do this by finding the displacement of each point to the axis, then finding the largest distance between two points.
perp_disp = (np.cross(p2-p1, M-p1))/norm(p2-p1) # Get perpendicular displacement vectors.
print(perp_disp)
[[-3.40206909 1.36082763 0. ]
[ 0. -0.13608276 0.13608276]
[ 1.36082763 -2.04124145 1.4969104 ]
[-2.72165527 0. 1.08866211]
[-1.36082763 -1.90515869 2.44948974]
[ 0.68041382 -0.95257934 0.68041382]
[ 2.72165527 0.68041382 -1.76907593]
...
[ 0. 0.27216553 -0.27216553]
[ 0. -0.40824829 0.40824829]
[ 2.72165527 0.27216553 -1.36082763]
[ 2.04124145 -0.68041382 -0.13608276]]
The largest distance is found by doing the same pdist used trick above.
max_perp_disp_pair = list(pd.DataFrame(squareform(pdist(perp_disp))).stack().idxmax())
perp_p1 = M[max_perp_disp_pair[0]]
perp_p2 = M[max_perp_disp_pair[1]]
print(perp_p1, perp_p2)
[ 1 2 -3] [-3 -2 1]
Finally, we get the diameter of the cylinder.
print(norm(perp_p1 - perp_p2))
6.92820323028
The smallest volume of a cylinder that can contain these points is
The maximum distance between points was found usings Numpy's pairwise distance function pdist. Then it was formatted with squareform to throw it into a Pandas DataFrame so the indices of both points could be easily found using idxmax. There is probably a better way to do this without Pandas.
If the np.cross part left you scratching your head, this is simply to find the minimum distance between a point and a line. I can followup with more detail if you're interested, but if you draw a cross product of two lines you get a parallelogram where the two non-parallel sides are given by the lines. This parallelogram has the same area as a rectangle with length equal to one of the lines, and width equal to the distance from the point to the line.
150
compute OBB
so either use PCA or this
to obtain 3D OBB. The code in the link must be ported to 3D but the principle is the same. Here my more advanced 3D OBB approximation using recursive search on cube_map (the code and approach in here is inferior to it).
initial guess
so OBB will give you oriented bounding box. Its biggest side will be parallel to rotation axis of your cylinder. So lets start with cylinder outscribing this OBB. So the central axis will be center of the OBB and parallel to its biggest side. (if you do not have biggest side then you need to check all 3 combinations). The diameter will be the bigger of the remaining sides.
fit cylinder
Now just try "all" combinations of offset and radius (may be also height) enclosing all your points near initial guess and remember the best one (according to your wanted specs). You can use any optimization method for this but my favorite is this:
The validity of the result depends on the fitting process. But do not get too wild with the nested fitting as the complexity goes wild too.
[Edit1] 3D OBB in C++
I was curious and got some time today so I encoded 3D OBB similar to the 2D example linked above. Looks like its working. This is preview:
I used 2 clusters to verify the orientation... Here the code (in form of simple C++ class):
//---------------------------------------------------------------------------
class OBB3D
{
public:
double p0[3],u[3],v[3],w[3]; // origin,3 axises sorted by size asc
double V,l[3]; // volume, and { |u|,|v|,|w| }
double p[8*3]; // corners
OBB3D() {}
OBB3D(OBB3D& a) { *this=a; }
~OBB3D() {}
OBB3D* operator = (const OBB3D *a) { *this=*a; return this; }
//OBB3D* operator = (const OBB3D &a) { ...copy... return this; }
void compute(double *pnt,int num) // pnt[num] holds num/3 points
{
OBB3D o; // temp OBB values
int i,j;
double a,_a,a0,a1,da,ca,sa; int ea; // search angles
double b,_b,b0,b1,db,cb,sb; int eb;
double c,_c,c0,c1,dc,cc,sc; int ec;
double u0[3],v0[3],pmin[3],pmax[3],q,*qq;
const double deg=M_PI/180.0;
p0[0]=0.0; u[0]=1.0; v[0]=0.0; w[0]=0.0; l[0]=0.0; V=-1.0;
p0[1]=0.0; u[1]=0.0; v[1]=1.0; w[1]=0.0; l[1]=0.0;
p0[2]=0.0; u[2]=0.0; v[2]=0.0; w[2]=1.0; l[2]=0.0;
if (num<3) { V=0.0; return; }
a0=0; a1=360.0*deg; da=10.0*deg; _a=a0;
b0=0; b1= 90.0*deg; db=10.0*deg; _b=b0;
c0=0; c1= 90.0*deg; dc=10.0*deg; _c=c0;
// recursively increase precision
for (j=0;j<5;j++)
{
// try all 3D directions with some step da,db
for (ea=1,a=a0;ea;a+=da){ if (a>=a1) { a=a1; ea=0; } ca=cos(a); sa=sin(a);
for (eb=1,b=b0;eb;b+=db){ if (b>=b1) { b=b1; eb=0; } cb=cos(b); sb=sin(b);
// spherical to cartesian direction
o.w[0]=cb*ca;
o.w[1]=cb*sa;
o.w[2]=sb;
// init u,v from cross product
vector_ld(u0,1.0,0.0,0.0);
if (fabs(vector_mul(u0,o.w))>0.75) // |dot(u,w)>0.75| measn near (anti)parallel
vector_ld(u0,0.0,1.0,0.0);
vector_mul(v0,o.w,u0); // v0 = cross(w,u0)
vector_mul(u0,v0,o.w); // u0 = cross(v0,w)
vector_one(u0,u0); // u0/=|u0|
vector_one(v0,v0); // v0/=|v0|
// try all rotations within u0,v0 plane
for (ec=1,c=c0;ec;c+=dc){ if (c>=c1) { c=c1; ec=0; } cc=cos(c); sc=sin(c);
for (i=0;i<3;i++)
{
o.u[i]=(u0[i]*cc)-(v0[i]*sc);
o.v[i]=(u0[i]*sc)+(v0[i]*cc);
}
// now u,v,w holds potential obb axises socompute min,max
pmin[0]=pmax[0]=vector_mul(pnt,o.u); // dot(pnt,u);
pmin[1]=pmax[1]=vector_mul(pnt,o.v); // dot(pnt,v);
pmin[2]=pmax[2]=vector_mul(pnt,o.w); // dot(pnt,w);
for (i=0;i<num;i+=3)
{
q=vector_mul(pnt+i,o.u); if (pmin[0]>q) pmin[0]=q; if (pmax[0]<q) pmax[0]=q;
q=vector_mul(pnt+i,o.v); if (pmin[1]>q) pmin[1]=q; if (pmax[1]<q) pmax[1]=q;
q=vector_mul(pnt+i,o.w); if (pmin[2]>q) pmin[2]=q; if (pmax[2]<q) pmax[2]=q;
}
// compute V,l from min,max
for (o.V=1.0,i=0;i<3;i++) { o.l[i]=pmax[i]-pmin[i]; o.V*=o.l[i]; }
// remember best solution u,v,w,V,l and compute p0
if ((V<0.0)||(V>o.V))
{
*this=o; _a=a; _b=b; _c=c;
for (i=0;i<3;i++) p0[i]=(pmin[0]*u[i])+(pmin[1]*v[i])+(pmin[2]*w[i]);
}
}
}}
a0=(_a-0.5*da); a1=a0+da; da*=0.1;
b0=(_b-0.5*db); b1=b0+db; db*=0.1;
c0=(_c-0.5*dc); c1=c0+dc; dc*=0.1;
}
// sort axises
{ i=0; qq=u; } // w is max
if (l[1]>l[i]){ i=1; qq=v; }
if (l[2]>l[i]){ i=2; qq=w; }
for (j=0;j<3;j++) { q=w[j]; w[j]=qq[j]; qq[j]=q; } q=l[2]; l[2]=l[i]; l[i]=q;
{ i=0; qq=u; } // v is 2nd max
if (l[1]>l[i]){ i=1; qq=v; }
for (j=0;j<3;j++) { q=v[j]; v[j]=qq[j]; qq[j]=q; } q=l[1]; l[1]=l[i]; l[i]=q;
// compute corners from p0,u,v,w,l
for (i=0;i<3;i++)
{
j=i;
p[j]=p0[i] ; j+=3;
p[j]=p0[i]+(l[0]*u[i]) ; j+=3;
p[j]=p0[i]+(l[0]*u[i])+(l[1]*v[i]) ; j+=3;
p[j]=p0[i] +(l[1]*v[i]) ; j+=3;
p[j]=p0[i] +(l[2]*w[i]); j+=3;
p[j]=p0[i]+(l[0]*u[i]) +(l[2]*w[i]); j+=3;
p[j]=p0[i]+(l[0]*u[i])+(l[1]*v[i])+(l[2]*w[i]); j+=3;
p[j]=p0[i] +(l[1]*v[i])+(l[2]*w[i]); j+=3;
}
}
void gl_draw()
{
glBegin(GL_LINES);
glVertex3dv(p+ 0); glVertex3dv(p+ 3);
glVertex3dv(p+ 3); glVertex3dv(p+ 6);
glVertex3dv(p+ 6); glVertex3dv(p+ 9);
glVertex3dv(p+ 9); glVertex3dv(p+ 0);
glVertex3dv(p+12); glVertex3dv(p+15);
glVertex3dv(p+15); glVertex3dv(p+18);
glVertex3dv(p+18); glVertex3dv(p+21);
glVertex3dv(p+21); glVertex3dv(p+12);
glVertex3dv(p+ 0); glVertex3dv(p+12);
glVertex3dv(p+ 3); glVertex3dv(p+15);
glVertex3dv(p+ 6); glVertex3dv(p+18);
glVertex3dv(p+ 9); glVertex3dv(p+21);
glEnd();
}
} obb;
//---------------------------------------------------------------------------
You just call compute with point cloud data where num is 3x number of points. The result is stored as unit basis vectors u,v,w and origin p0 along with sizes l[] per each axis or as 8 corner points of OBB p
The stuff works simply by trying "all" spherical positions with some step for the w axis and then try all u,v polar positions perpendicular to each and w remembering the minimal volume OBB. Then recursively search only positions near found best solution with smaller step to improve accuracy.
I think this should provide fine start point. If you implement the minimal circle instead of u,v rotation (loop for (ec=1,c=c0;ec;c+=dc)) then you might obtain your cylinder directly from this search.
The code is not optimized yet (some parts like w axis check) can be moved to lower layer of nested for loop. But I wanted to keep this simple and understandable as much as I could instead.
[Edit2] 3D OBC in C++
I managed to modify my 3D OBB by replacing U,V search with minimal enclosing circle (hope I implemented it right but it looks like it...) that find the minimal enclosing 2D circle of all the points projected on UV plane which makes it an Oriented Bounding Cylinder parallel to W. I used the first approach from pdf from your link (using bisector). Here the result:
In blue is the 3D OBB and in brown/orange-ish is the found 3D OBC. Here the code:
class OBC3D // 3D Oriented Bounding Cylinder
{
public:
double p0[3],u[3],v[3],w[3]; // basecenter,3 axises
double V,r,h; // volume, radius height
double p1[3]; // other base center
OBC3D() {}
OBC3D(OBC3D& a) { *this=a; }
~OBC3D() {}
OBC3D* operator = (const OBC3D *a) { *this=*a; return this; }
//OBC3D* operator = (const OBC3D &a) { ...copy... return this; }
void compute(double *pnt,int num) // pnt[num] holds num/3 points
{
OBC3D o; // temp OBB values
int i,j,k,kk,n;
double a,_a,a0,a1,da,ca,sa; int ea; // search angles
double b,_b,b0,b1,db,cb,sb; int eb;
double pmin[3],pmax[3],q,qq,*pnt2,p[3],c0,c1,u0,v0,du,dv,dr;
const double deg=M_PI/180.0;
p0[0]=0.0; u[0]=1.0; v[0]=0.0; w[0]=0.0; V=-1.0;
p0[1]=0.0; u[1]=0.0; v[1]=1.0; w[1]=0.0; r=0.0;
p0[2]=0.0; u[2]=0.0; v[2]=0.0; w[2]=1.0; h=0.0;
if (num<3) { V=0.0; return; }
// prepare buffer for projected points
pnt2=new double[num];
a0=0; a1=360.0*deg; da=10.0*deg; _a=a0;
b0=0; b1= 90.0*deg; db=10.0*deg; _b=b0;
// recursively increase precision
for (k=0;k<5;k++)
{
// try all 3D directions with some step da,db
for (ea=1,a=a0;ea;a+=da){ if (a>=a1) { a=a1; ea=0; } ca=cos(a); sa=sin(a);
for (eb=1,b=b0;eb;b+=db){ if (b>=b1) { b=b1; eb=0; } cb=cos(b); sb=sin(b);
// spherical to cartesian direction
o.w[0]=cb*ca;
o.w[1]=cb*sa;
o.w[2]=sb;
// init u,v from cross product
vector_ld(o.u,1.0,0.0,0.0);
if (fabs(vector_mul(o.u,o.w))>0.75) // |dot(u,w)>0.75| measn near (anti)parallel
vector_ld(o.u,0.0,1.0,0.0);
vector_mul(o.v,o.w,o.u); // v0 = cross(w,u0)
vector_mul(o.u,o.v,o.w); // u0 = cross(v0,w)
vector_one(o.u,o.u); // u0/=|u0|
vector_one(o.v,o.v); // v0/=|v0|
// now u,v,w holds potential obb axises so compute min,max and convert to local coordinates
pmin[0]=pmax[0]=vector_mul(pnt,o.u); // dot(pnt,u);
pmin[1]=pmax[1]=vector_mul(pnt,o.v); // dot(pnt,v);
pmin[2]=pmax[2]=vector_mul(pnt,o.w); // dot(pnt,w);
for (i=0;i<num;i+=3)
{
q=vector_mul(pnt+i,o.u); if (pmin[0]>q) pmin[0]=q; if (pmax[0]<q) pmax[0]=q; pnt2[i+0]=q;
q=vector_mul(pnt+i,o.v); if (pmin[1]>q) pmin[1]=q; if (pmax[1]<q) pmax[1]=q; pnt2[i+1]=q;
q=vector_mul(pnt+i,o.w); if (pmin[2]>q) pmin[2]=q; if (pmax[2]<q) pmax[2]=q; pnt2[i+2]=q;
}
// [compute min enclosing circle]
n=0;
// center (u0,v0) = avg( pnt2 )
for (u0=0.0,v0=0.0,i=0;i<num;i+=3)
{
u0+=pnt2[i+0];
v0+=pnt2[i+1];
} q=3.0/double(num); u0*=q; v0*=q;
// r = max(|pnt2 - (u0,v0)|)
for (o.r=0.0,i=0;i<num;i+=3)
{
c0=pnt2[i+0]-u0;
c1=pnt2[i+1]-v0;
q=(c0*c0)+(c1*c1);
if (o.r<q) o.r=q;
} o.r=sqrt(o.r);
for (kk=0;kk<4;kk++)
{
// update edgepoints count n
qq=o.r*o.r;
for (i=n;i<num;i+=3)
{
c0=pnt2[i+0]-u0;
c1=pnt2[i+1]-v0;
q=fabs((c0*c0)+(c1*c1)-qq);
if (q<1e-10)
{
pnt2[n+0]=pnt2[i+0];
pnt2[n+1]=pnt2[i+1];
pnt2[n+2]=pnt2[i+2]; n+=3;
}
}
// compute bisector (du,dv)
for (du=0.0,dv=0.0,i=0;i<n;i+=3)
{
du+=pnt2[i+0]-u0;
dv+=pnt2[i+1]-v0;
} q=1.0/sqrt((du*du)+(dv*dv)); du*=q; dv*=q;
// try to move center towards edge points as much as possible
for (dr=0.1*o.r,j=0;j<5;)
{
u0+=dr*du;
v0+=dr*dv;
// q = max(|pnt2 - (u0,v0)|)
for (qq=0.0,i=0;i<num;i+=3)
{
c0=pnt2[i+0]-u0;
c1=pnt2[i+1]-v0;
q=(c0*c0)+(c1*c1);
if (qq<q) qq=q;
} qq=sqrt(qq);
// recursively increase precision
if (qq>o.r)
{
u0-=dr*du;
v0-=dr*dv;
dr*=0.1;
j++;
}
else o.r=qq;
}
}
// compute h,V
o.h=pmax[2]-pmin[2];
o.V=M_PI*o.r*o.r*o.h;
// remember best solution u,v,w,V,l and compute p0
if ((V<0.0)||(V>o.V))
{
*this=o; _a=a; _b=b;
for (i=0;i<3;i++) p0[i]=(u0*u[i])+(v0*v[i])+(pmin[2]*w[i]);
}
}}
a0=(_a-0.5*da); a1=a0+da; da*=0.1;
b0=(_b-0.5*db); b1=b0+db; db*=0.1;
}
// compute corners from p0,u,v,w,l
for (i=0;i<3;i++) p1[i]=p0[i]+(h*w[i]);
delete[] pnt2;
}
void gl_draw()
{
int i,j,n=36;
double a,da=2.0*M_PI/double(n),p[3],uu,vv;
glBegin(GL_LINES);
glVertex3dv(p0); glVertex3dv(p1);
glEnd();
glBegin(GL_LINE_LOOP);
for (a=0.0,i=0;i<n;i++,a+=da)
{
uu=r*cos(a);
vv=r*sin(a);
for (j=0;j<3;j++) p[j]=p0[j]+(u[j]*uu)+(v[j]*vv);
glVertex3dv(p);
}
glEnd();
glBegin(GL_LINE_LOOP);
for (a=0.0,i=0;i<n;i++,a+=da)
{
uu=r*cos(a);
vv=r*sin(a);
for (j=0;j<3;j++) p[j]=p1[j]+(u[j]*uu)+(v[j]*vv);
glVertex3dv(p);
}
glEnd();
}
};
//---------------------------------------------------------------------------
Usage is the same ... I tested with this:
OBB3D obb;
OBC3D obc;
void compute()
{
int i,n=500;
// random pnt cloud
Randomize();
RandSeed=98123456789;
pnt.allocate(3*n); pnt.num=0;
// random U,V,W basis vectors
double u[3],v[3],w[3],x,y,z,a;
for (i=0;i<3;i++) w[i]=Random()-0.5; // random direction
vector_one(w,w); // w/=|w|
vector_ld(u,1.0,0.0,0.0);
if (fabs(vector_mul(u,w))>0.75) // |dot(u,w)>0.75| measn near (anti)parallel
vector_ld(u,0.0,1.0,0.0);
vector_mul(v,w,u); // v = cross(w,u)
vector_mul(u,v,w); // u = cross(v,w)
vector_one(u,u); // u/=|u|
vector_one(v,v); // v/=|v|
// random cylinder point cloud
for (i=0;i<n;i++)
{
a=2.0*M_PI*Random();
x= 0.5+(0.75*(Random()-0.5))*cos(a);
y=-0.3+(0.50*(Random()-0.5))*sin(a);
z= 0.4+(0.90*(Random()-0.5));
pnt.add((x*u[0])+(y*v[0])+(z*w[0]));
pnt.add((x*u[1])+(y*v[1])+(z*w[1]));
pnt.add((x*u[2])+(y*v[2])+(z*w[2]));
}
obb.compute(pnt.dat,pnt.num);
obc.compute(pnt.dat,pnt.num);
}
Where List<double> pnt is my dynamic array template double pnt[]. Which is not important in here.
Beware that if you chose too big initial step (da,db) for the W direction search you might miss the correct solution by trapping itself inside local minimum.
21
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
