Clustering with points lying along different 3D planes

I would try using a Hough Transform, easy to find on Wikipedia. There's a thorough discussion there about using the HT to classify lines from noisy images (a 2-parameter problem). Your task of classifying planes from noisy point clouds is similar, but it's a 3-parameter problem .. two spherical angles can describe the direction of the plane normal, then one more scalar to represent the distance from the origin.

This approach is probably a little old-school, a computer vision expert might have better ideas.

EDIT: If the plane always contains the origin you only have the angles as (two) free parameters. This is fortunate, because the HT scales terribly with dimension anyway.

Observe that any data-point $$x_i = begin{bmatrix} x_{i1}x_{i2}x_{i3} end{bmatrix}$$ can be interpreted as the three-vector pointing from the origin to the point itself (position vector).

1. For any two data points $x_i$ and $x_j$, form the cross-product vector $${u}_{ij} = x_i times x_j$$ Observe that the unique plane that passes through the origin and the two points $x_i$ and $x_j$ is also the unique plane that passes through the origin and is orthogonal to $u_{ij}$. If the latter turns out to be zero, then get rid of one of the points.

2. For each vector $u_{ij}$, create a pair of points ${- p_{ij}, , p_{ij}}$ where $$p_{ij} = frac{u_{ij}}{|u_{ij}|}$$ The points $p_{ij}$ and $-p_{ij}$ lie on the unit sphere $mathbb{S}^2 = {x in mathbb{R}^3 , : , |x| = 1 }$ and are antipodal. You obtain a set of points ${p_{ij}}$ on the unit sphere, which are symmetric with respect to the antipodal symmetry.

3. This point it the difficult one. You need to choose a clustering algorithm for the points ${p_{ij}}$ on the sphere, so that you can detect high accumulation ofpoints. To that end, you could use a convex hull algorithm for ${p_{ij}}$, which will give you the Delaunay decomposition of ${p_{ij}}$ on the unit sphere. It's dual is the Voronoi (nearest neughbour) decomposition of ${p_{ij}}$. Then you can use this information in your clustering methodology.

4. As a result of the (non-trivial) step 3, you end up with a bunch of point-clusters $N_1, , N_2, , ..., , N_{2k}$ (you should expect them to have the antipodal symmetry). Take only half of the clusters, say $N_1, , N_2, , ..., , N_{2k}$ so that no two clusters are antipodally symmetric, and ignore their antipodal images. Take the average point (unit) vector for each cluster $N_s$: $$n_s = frac{1}{|N_s|} , sum_{p_{ij} in N_s} , p_{ij}$$ The normal unit vectors $n_s$ should be the normal vectors of the planes that you seek.

5. You may and up with some redundant planes, so you may want to get rid of some of them. For each plane, represented by a unit normal vector $n_s$, calculate the distances between the points and the planes, i.e. for each $x_i$ calculate $h_{is} = |(x_i circ n_s)|$ and, say using some reasonable treshold, select only the points which are close enough to the plane of $n_s$ (so those with small h_{is}), and project them onto the plane: $x_{is} = x_i - (x_icirc n_s), n_s$. If the points $x_{is}$ form a fairly well spread-out, diffuse, pattern on the plane $n_s$, you keep the plane. If however, the points $x_{is}$ form something that looks very much like a smeared (cloudy), set of lines passing through the point of origin, which is on the plane, then you get rid of the plane.

You could also work in the projective plane instead of the unit sphere, it is just that the notion of distance might be more cumbersome, but I am not sure.

Looking at the other answers, the recommendation is something called "Hough transform". Its philosophy however is not that much different from what I sketched above. I recommend the use of projective duality (like what I am using in the methodology outlined above) from projective geometry point of view. In two dimensions, it utilizes the construction that the family of all lines passing through a given point is equivalent to a (dual) line. Consequently, given two points $p_1$ and $p_2$ in the plane, the family of all lines passing through the point $p_1$ forms a dual line $l_1$, the family of all lines passing through the point $p_2$ forms a dual line $l_2$, so the intersection of the two (dual) lines $l_1$ and $l_2$ is a (dual) point, which represents the common line between two two original points $p_1$ and $p_2$. This has the same philosophy as this "Hough transform" but is linear in nature, while the "Hough transform" uses curves, which I think is more cumbersome. Even with the Hough transform, you would need to find a way to detect accumulations/clusters of points.

Like rchilton1980, I would suggest the Hough transform. In your case, d does not to be quantized because it is known to be zero. When quantizing the two angles, beware that the naive quantization ($phi + ncdotDeltaphi$, $theta +ncdotDeltatheta$) will result in an uneven angle distirbution. I would recommend to use the regularized quantiziation by Jeltsch (the title refers to line detection, but the quantization method applies to any 3D angle quantization):

Jeltsch et al.: "Hough Parameter Space Regularisation for Line Detection in 3D." VISAPP 2016, pp. 345-352 (2016)

For improving the result, an orthogonal least squares fit (which is mathematically equivalent to take the two main components of a PCA) through the points of each plane is a good idea. The plane normal is then simply the eigenvector of teh covariance matrix corresponding to the smallest eigenvalue:

Späth: "Orthogonal least squares fitting with linear manifolds." Numerische Mathematik 48, pp. 441–445 (1986)