BTW 复杂度确实是 O(n^3) 降低到你需要的:>

BTW complexity is indeed O(n^3) to lower that you need to:

1. 以某种方式对点进行排序

by x 和或 y 按升序或降序排列.有时使用极坐标也有帮助

by x and or y in ascending or descending order. Also use of polar coordinates can help sometimes

在impera(分而治之)算法中使用除法

通常对于平面几何算法来说，将区域划分为象限和子象限是个好主意，但这些算法很难在矢量图形上编码

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

还有另一种加速可能性

检查所有可能的方向(数量有限，例如仅限360角度)，导致O(n^2).然后计算结果仍然是 O(m^3) 其中 m 是每个测试方向的点的子集.

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.

好的，这是我用 C++ 编写的基本代码，例如:

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[] 是输入 2D 点(我在随机线周围输入随机点 + 随机噪声点)
view.pnt.num 是点数

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[] 是输出行(只用了一行)

view.lin[] are output lines (just one line is used)

准确性

使用 pacc,lmin,lmax 常量来改变行为和计算速度.改变dirs改变方向精度和计算速度

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