Generalizing planar detection to k-plane detection

In this post the planar detection algorithm in $\R^3$ of Bauer and Polthier in Detection of Planar Regions in Volume Data for Topology Optimization is generalized to detect $k$-planes with largest density in $\R^n$. Let $\Omega\subset \R^n$ be the compact support of a piecewise-constant probability density function $\rho:\R^n\to \R_{\geqslant 0}$.

Definition: Let $(G,\rho)$ be a grid, where $G \subset \lambda \Z^n + c \subset \R^n$ is a lattice in $\Omega$.  A cell $x$ of the grid is $B_\infty(x,\lambda/2) = \{y\in \R^n\ :\ ||x-y||_\infty\leqslant\lambda/2\}$, for $x\in G$. Every cell is assigned a value
\[
\int_{B_\infty(x,\lambda/2)}\rho\ dx,
\]
called the mass of the cell, which may be though of as a type of Radon transform of $\rho$.

Assuming that $k$ is a global variable, running $Recursive(G,w,k)$ will give the desired result. This algorithm is the naive generalization of Bauer and Polthier, and suffers from calculating mass along the same $k$-plane several times, whenever $k<n-1$ (as any $k$-plane does not lie in a unique $(k+1)$-plane).

---

Algorithm: $k$PlaneFinder
$Recursive(G,w,k)$:
Input: A grid $(G,\rho)$, a width $w$ of fattened $k$-planes, the current plane dimension $k\leqslant k'<n$
Output: A $k$-planar connected component covering most mass in $G$

discretize the unit $(k'-1)$-hemisphere in an appropriate manner
order the vertices by a Hamiltonian path
for each vertex $\textbf{n}$:
    sort the grid cells in direction $\textbf{n}$
    discretize the range in direction $\textbf{n}$ equidistantly
    for each $k'$-plane $(\textbf{n},d)$:
        collect the cells closer than $w$ to the $k'$-plane into a graph $G'$
        if $k'\neq k$:
            run $Recursive(G',w,k'-1)$
        else:
            compute the connected component having the most mass in $G'$
return the connected $k$-component having most mass (and the corresponding $k$-plane)

---

Measuring along connected components of a $k$-plane works the same way as in the original version, as the gird on $\R^n$ similarly induces a connectivity graph.

Remark: Bauer and Polthier cite Kantaforoush and Shahshahani in evenly sampling points on the unit 2-sphere, but it is not clear how their method (using the inscribed icosahedron) generalizes. Another method would be uniformly sampling random points on $S^{k-1}$ and take all on one hemisphere. A Hamiltonian path could then be taken from an arbitrary point and then using the greedy algorithm (with respect to Euclidean distance) to find consecutive vertices (to keep down the time of consecutive sorting operations).

Recall the Grassmannian $Gr(n,k)$ of all $k$-planes in $\R^n$ through the origin, a compact manifold of dimension $k(n-k)$. Note that any $k$-plane $P\subset \R^n$ is a translation of an element $Q\in Gr(n,k)$ by an element of $Q^\perp$ (we conflate notation for $Q$ and its natural embedding in $\R^n$).

Remark: $Gr(n,k)$ is parametrizable, so by choosing directions in the unit $(n-k)$-hemisphere, the process of choosing $k$-planes in the algorithm may be completely parametrized. The quick sorting of points that was available in Bauer and Polthier's $n=3,k=2$ case may be replaced by an iterated restriction of the original data set through a complete flag $P\subset \cdots \subset \R^n$.

References: Bauer and Polthier (Detection of Planar Regions in Volume Data for Topology Optimization), Katanforoush and Shahshahani (Distributing points on the Sphere 1)