I think your problem is in the breath of the search. If you make this a problem of finding sets of three points, then you have trillions of such sets to search over.
So reduce the problem size, not by throwing away points, but by looking for what factor do sets of three collinear points have in common. Not worded well, but you should see the idea once I get going. I want to reduce a doubly or triply nested problem into a single loop.
For example, consider the point cloud
N = 1e5;
xyz = round(rand(N,3)*1000);
xyz = unique(xyz,'rows');
N = size(xyz,1)
I've rounded them onto an integer lattice, because in my example, I am hoping to find at least a few triplets of points that are truly collinear. In this exmple, I ended up with a set of 99996 distinct points.
Consider the first point in that point cloud. Is it collinear with any other two points? If it is, then there exists one line that passses through point pairs (1,i) and (1,j). So compute the slope of the lines that pass through all such pairs, not as a slope, of course, but as an angle. This does away with the issue of a vertical line, which would have an infinite slope. And, of course, this must be done in 3-dimensions. Again, lets just try it...
A = [atan2(xyz(2:end,2) - xyz(1,2),xyz(2:end,1) - xyz(1,1)),atan2(xyz(2:end,3) - xyz(1,3),xyz(2:end,1) - xyz(1,1))];
So A represents the inclination of each line that passes through point 1, and every other point in the cloud. Essentially, I've computed the slopes dy/dx amd dz/dx for each such line, but I've represented that in an angular form using atan2.
Now, we need to find rows of A that are the same, to within a tolerance. uniquetol does some of the work for us.
[Au,Ia,Ic] = uniquetol(A,1e-12,'byrows','true');
Now we are getting someplace. There wer 99942 distinct sets of slopes. But that means there were at least a few sets of points that are collinear, and include point #1. We can find them using the Indicator varisble Ic.
colind = find(accumarray(Ic,1) > 1)
So there were 6 distinct lines passing through point #1, such that the line in question also passes through more than one other point. We can easily now go back and find those sets, saving them away.
Having done this for point #1, we need never look at point #1 again, as we have found every line that includes point # 1 and any other two or more points. This means the array A will get progressively smaller as we run along. But it also means the entire process cn be done using just a single loop over the point which may be in common.
In fact, the above scheme will be pretty efficient, since it is just a single loop. Inside the loop create the array of inclinations, then call uniquetol, and use accuarray to find the sets with common inclinations. This will be eminently doable, even for considerably larger point clouds.