BTW complexity is indeed O(n^3)
to lower that you need to:
sort the points somehow
by x
and or y
in ascending or descending order. Also use of polar coordinates can help sometimes
use divide at impera (divide and conquer) algorithms
usually for planar geometry algorithms is good idea to divide area to quadrants and sub-quadrants but these algorithms are hard to code on vector graphics
Also there is one other speedup possibility
check against all possible directions (limited number of them for example to 360
angles only) which leads to O(n^2)
. Then compute results which is still O(m^3)
where m
is the subset of points per the tested direction.
Ok here is something basic I coded in C++ for example:
void points_on_line()
{
const int dirs =360; // num of directions (accuracy)
double mdir=double(dirs)/M_PI; // conversion from angle to code
double pacc=0.01; // position acc <0,1>
double lmin=0.05; // min line size acc <0,1>
double lmax=0.25; // max line size acc <0,1>
double pacc2,lmin2,lmax2;
int n,ia,ib;
double x0,x1,y0,y1;
struct _lin
{
int dir; // dir code <0,dirs>
double ang; // dir [rad] <0,M_PI>
double dx,dy; // dir unit vector
int i0,i1; // index of points
} *lin;
glview2D::_pnt *a,*b;
glview2D::_lin q;
_lin l;
// prepare buffers
n=view.pnt.num; // n=number of points
n=((n*n)-n)>>1; // n=max number of lines
lin=new _lin[n]; n=0;
if (lin==NULL) return;
// precompute size of area and update accuracy constants ~O(N)
for (a=view.pnt.dat,ia=0;ia<view.pnt.num;ia++,a++)
{
if (!ia)
{
x0=a->p[0]; y0=a->p[1];
x1=a->p[0]; y1=a->p[1];
}
if (x0>a->p[0]) x0=a->p[0];
if (x1<a->p[0]) x1=a->p[0];
if (y0>a->p[1]) y0=a->p[1];
if (y1<a->p[1]) y1=a->p[1];
}
x1-=x0; y1-=y0; if (x1>y1) x1=y1;
pacc*=x1; pacc2=pacc*pacc;
lmin*=x1; lmin2=lmin*lmin;
lmax*=x1; lmax2=lmax*lmax;
// precompute lines ~O((N^2)/2)
for (a=view.pnt.dat,ia=0;ia<view.pnt.num;ia++,a++)
for (b=a+1,ib=ia+1;ib<view.pnt.num;ib++,b++)
{
l.i0=ia;
l.i1=ib;
x0=b->p[0]-a->p[0];
y0=b->p[1]-a->p[1];
x1=(x0*x0)+(y0*y0);
if (x1<=lmin2) continue; // ignore too small lines
if (x1>=lmax2) continue; // ignore too big lines
l.ang=atanxy(x0,y0);
if (l.ang>M_PI) l.ang-=M_PI; // 180 deg is enough lines goes both ways ...
l.dx=cos(l.ang);
l.dy=sin(l.ang);
l.dir=double(l.ang*mdir);
lin[n]=l; n++;
// q.p0=*a; q.p1=*b; view.lin.add(q); // just visualise used lines for testing
}
// test directions
int cnt,cntmax=0;
double t;
for (ia=0;ia<n;ia++)
{
cnt=1;
for (ib=ia+1;ib<n;ib++)
if (lin[ia].dir==lin[ib].dir)
{
a=&view.pnt[lin[ia].i0];
if (lin[ia].i0!=lin[ib].i0)
b=&view.pnt[lin[ib].i0];
else b=&view.pnt[lin[ib].i1];
x0=b->p[0]-a->p[0]; x0*=x0;
y0=b->p[1]-a->p[1]; y0*=y0;
t=sqrt(x0+y0);
x0=a->p[0]+(t*lin[ia].dx)-b->p[0]; x0*=x0;
y0=a->p[1]+(t*lin[ia].dy)-b->p[1]; y0*=y0;
t=x0+y0;
if (fabs(t)<=pacc2) cnt++;
}
if (cntmax<cnt) // if more points on single line found
{
cntmax=cnt; // update point count
q.p0=view.pnt[lin[ia].i0]; // copy start/end point
q.p1=q.p0;
q.p0.p[0]-=500.0*lin[ia].dx; // and set result line as very big (infinite) line
q.p0.p[1]-=500.0*lin[ia].dy;
q.p1.p[0]+=500.0*lin[ia].dx;
q.p1.p[1]+=500.0*lin[ia].dy;
}
}
if (cntmax) view.lin.add(q);
view.redraw=true;
delete lin;
// Caption=n; // just to see how many lines actualy survive the filtering
}
It is from my geometry engine so here is some stuff explained:
glview2D::_pnt
view.pnt[]
are input 2D points (I feed random points around random line + random noise points)
view.pnt.num
is number of points
glview2D::_lin
view.lin[]
are output lines (just one line is used)
accuracy
Play with pacc,lmin,lmax
constants to change the behavior and computation speed. Change dirs
to change direction accuracy and computation speed
Complexity estimation is not possible due to big dependency on input data
But for my random test points are the runtimes like this:
[ 0.056 ms]Genere 100 random 2D points
[ 151.839 ms]Compute 100 points on line1 (unoptimized brute force O(N^3))
[ 4.385 ms]Compute 100 points on line2 (optimized direction check)
[ 0.096 ms] Genere 200 random 2D points
[1041.676 ms] Compute 200 points on line1
[ 39.561 ms] Compute 200 points on line2
[ 0.440 ms] Genere 1000 random 2D points
[29061.54 ms] Compute 1000 points on line2
与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…