Given n points on a 2D plane, find the maximum number of points that lie on the same straight line

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

1. sort the points somehow

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

2. 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

3. 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