Robust orientation of a point cloud

Question

I have 2D point clouds which are 4-way symmetrical (invariant by 90° rotation). The points are usually arranged on the nodes of a square grid, densely populated, but some cases can be more complicated. I do know the point pattern in advance.

I need to find the symmetry axis, which are arbitrarily oriented. I initially thought of the moments of inertia, but in the case of symmetrical clouds the ellipse of inertia is a circle and no orientation information can be obtained.

In addition, the method should be robust because a few points (say <2%) can be missing, displaced or extraneous. There can be between 500 and 2500 points in total.

Any suggestion about how this can done efficiently ? I suspect that some other moments can do.

Typical case:

Answer 1

A point cloud which is 4-way symmetrical will have 4 axes of reflectional symmetry, where each axis is offset from the next by 45 degrees. To find all axes it therefore suffices to find one.

The following home-made algorithm works pretty well at finding that axis. Trial testing of clouds with an average of 500 points gives the following statistics regarding the deviation $D$ between the axis angle found by the algorthm and the true axis angle:

Perfect cloud:

$D \lt 1^\circ$: $100$%

Cloud with 2% missing points:

$D \lt 1^\circ$: $55$%

$D \lt 2^\circ$: $90$%

$D \lt 3^\circ$: $99$%

I'll discuss problem areas at the end of the post, but first to the algorithm.

Algorithm overview

For a 4-way symmetrical cloud, at least one axis of symmetry must exist within any $45^\circ$ sector, as seen from the centroid of the cloud. At the extreme, there will be one axis at the very start of the sector and another axis at the very end of the sector. An axis of symmetry is characterised, in this case, by the fact that all points within a $45^\circ$ sector to the "left" of the axis are mirrored by the points within a $45^\circ$ sector to the "right" of the axis. A very simple (non-robust) way to check if an axis at a given angle is an axis of symmetry, is to see if the number of points in the $45^\circ$ sector to the left of the axis is equal to the number of points in the $45^\circ$ sector to the right of the axis. A better way (which this algorithm uses) is to check if the sum of the distances from the centroid to the points in the $45^\circ$ sector to the left of the axis, is equal to the sum of the distances to the points in the $45^\circ$ sector to the right of the axis.

So, the basic idea of the algorithm is to have an axis move from $91^\circ$ to $44^\circ$ in $1^\circ$ steps, calculating the sum of the distances to the points in the left sector and the sum of the distances to the right sector at each step, and comparing these sums. The angle of the axis where the sums are "most equal" is then the algorithm's axis of symmetry.

The reason for using $91^\circ$ and $44^\circ$ instead of $90^\circ$ and $45^\circ$ is to avoid rounding problems. In addition, I use sectors of $47^\circ$ (which means the left and right sector overlap) instead of sectors of $45^\circ$ to avoid problems with points which are at or close to the dividing line between the left and the right sector.

The above is obviuosly not a robust method in general, but for a cloud where 4-way symmetry is assumed, it works.

Algorithm with details

Assuming all points are in an array with Cartesian coordinates:

1. Find the centroid $(X_c, Y_c)$ of the cloud.

2. Find the polar coordinates of each point, using $(X_c, Y_c)$ as origo, and save the points which have an angle between $91 + 46 = 137$ degrees and $44 - 46 = -2$ degrees in a new array $NewA$.

3. Sort $NewA$ according to angle, from largest to smallest.

4. Run through the angles $A$ from $91^\circ$ down to $44^\circ$ in steps of $1^\circ$. At each angle, calculate the distance sum for points in the left sector (i.e. points with $ A - 1 \lt$ angle $\lt A + 46$) and compare with the distance sum for points in the right sector (i.e. points with $ A + 1 \lt$ angle $\lt A - 46$).

5. The angle $A$ at which this comparison is closest to one, is the algorithm's axis of symmetry angle.

Areas for improvement

1. The algorithm is sensitive to "error" in the coordinates of the centroid. When points are missing or displaced, the centroid is not at the ideal (correct) location.

2. The algorithm is also somewhat sensitive to the 2% errors if they mainly occur in the first and second quadrant (where this algorithm's input comes from).

Not sure what to do about these points.

Answer 2

This is not an answer, but an extended comment that might lead to a workable solution.

I did some experiments, generating histograms of $$\chi_{ij} = -\chi_{ji} = \frac{x_j - x_i}{\lvert x_j - x_i \rvert + \lvert y_j - y_i \rvert}$$ and $$\gamma_{ij} = -\gamma_{ji} = \frac{x_j - x_i}{\lvert x_j - x_i \rvert + \lvert y_j - y_i \rvert}$$ with each sample weighed by $$\omega_{ij} = \omega_{ji} = \frac{1}{(x_j - x_i)^2 + (y_j - y_i)^2}$$ from rotated and translated point clouds with four-way mirror symmetry, with up to 10% additional (non-symmetric) or missing points.

Note that for $$\theta_{ij} = \arctan\left(\frac{x_j - x_i}{y_j - y_i}\right)$$we have $$\theta_{ij} = \begin{cases} \arctan\left(\frac{1}{\chi_{ij}} - 1\right), & 0 \le \theta \lt \frac{\pi}{2} \\ \arctan\left(1 - \frac{1}{\chi_{ij}}\right), & -\frac{\pi}{2} \lt \theta \le 0 \end{cases}$$ and $$\theta_{ij} = \begin{cases} \arctan\left(\frac{\gamma_{ij}}{1 - \gamma_{ij}}\right), & 0 \lt \theta \le \frac{\pi}{2} \\ \arctan\left(\frac{\gamma_{ij}}{1 + \gamma_{ij}}\right), & -\frac{\pi}{2} \le \theta \lt 0 \end{cases}$$

For $N$ points, there are $N(N-1)/2$ unique $i , j$ pairs, and the above is computationally cheap to calculate, and is relatively easily vectorized, too; single-precision or fixed-point integer arithmetic suffices for the purposes of the histogram. Where $\arctan$ or $\operatorname{atan2}$ might be considered too slow, the above should be fast enough.

The reason I think this might lead to a workable solution, is that the histogram shows extremely clear spikes along the symmetry/mirror axes for random point clouds. (Because of the sign differences, you need to look at both to determine the actual symmetry axes.)

For some regular lattices, I saw more spikes (as one would expect), including dual spikes in certain cases, so this definitely needs further math work (either comprehensive numerical research of different point clouds, or analysis on how the four-way mirror symmetry causes the spikes seen in the angular histogram, if properly weighted by the distance) before I'd consider it a solution.