Find a line with maximum points from N points You are given N points and you want to draw a line such that maximum points lie on the line. What is the efficient way to find the number of maximum points?
 A: The trick is called duality. You map a point $(m,b)$ to a line that has that slope and intercept, and the inverse transformation is obvious to map a line to a point. So you compute the arrangement of lines, putting a bounding box big enough around the outside so you have a finite planar graph. Then you find the intersection point that has the most lines passing through it, and you are done. Note the number of intersection points is linear because the graph is planar.
A: Follows a script in MATHEMATICA which implements a procedure to obtain the line passing through more points.
NOTE
This script admits many improvements regarding computational efficiency. It was released in this basic form with the proposal to show the geometric ideas involved in the algorithm. Calling $p_i, i=1,\cdots n$ the set of points the first step is to form a matrix $V$ such that $V_{i,j} = p_i-p_j$. After, we construct a list $\mathcal{L}$ (list0) with the element indexes in $V$ such that for $i\ne j\ne k$ we have $V_{i,k}\times V_{j,k}=0$. The last operation consists in joining elements on $\mathcal{L}=\{\mathcal{l}_1,\cdots,\mathcal{l}_p,\cdots,\mathcal{l}_q,\cdots,\mathcal{l}_m\}$ such that $\mathcal{l}_p,\mathcal{l}_q$ have at least two indexes (points) in common. The biggest resulting elements in this union are the solutions.
n = 30;
SeedRandom[5]
list = {};
points = Union[RandomInteger[{0, 10}, {n, 2}]];
n = Length[points]
V = Table[Join[points[[i]] - points[[j]], {0}], {i, 1, n}, {j, 1, n}];
For[i = 1, i <= n, i++,
 For[j = 1, j <= n, j++,
  For[k = 1, k <= n, k++,
   If[i != k && j != k && i != j, 
    cross = Cross[V[[i, k]], V[[j, k]]];
    If[cross == {0, 0, 0}, AppendTo[list, Sort[{i, j, k}]]]]
  ]
 ]
]
list0 = Union[list];

lines = {};
npts = 0;
For[i = 1, i <= Length[list0], i++,
 lini = list0[[i]];
 line = lini;
 For[j = 1, j < i, j++,
  linj = list0[[j]];
  If[Length[Intersection[lini, linj]] == 2, line = Union[line,linj]]
 ];
 nl = Length[line];
 If[nl >= npts,
  AppendTo[lines, line];
  npts = nl
 ]
]
nlin = Length[lines]
LINES = {};
For[k = 1, k <= nlin, k++, 
 If[Length[lines[[k]]] == npts, AppendTo[LINES, lines[[k]]]]]
LINES = Union[LINES]
NL = Length[LINES];
gr1 = ListPlot[points];
gr2 = Flatten[Table[Graphics[{Red, PointSize[0.02], Point[points[[list0[[j, k]]]]]}], {k, 1, 3}, {j, 1, Length[list0]}]];
gr2a = Flatten[Table[Table[ParametricPlot[lambda points[[LINES[[j, k]]]] + (1 - lambda) points[[LINES[[j, k + 1]]]], {lambda, 0, 1}], {k, 1, Length[LINES[[j]]] - 1}], {j, 1, NL}]];
Show[gr1, gr2a, gr2, PlotRange -> All, Axes -> False, AspectRatio -> 1]


with SeedRandom[51]

and with SeedRandom[52]

