Coloring spiral structures

Question

In a previous question Original post @Julian S. proposed a method for coloring spiral structures. However this method works well only if the two spirals do not merge together.

For example, for these data we have

data = Import["spirals.dat", "Table"];
radius = 10.64;

data1 = FindClusters[Select[data, Norm[#] > radius &], 2, 
  Method -> "Agglomerate"];
L0 = Show[ListPlot[data1, PlotRange -> Full, AspectRatio -> 1, 
     PlotStyle -> {Directive[Darker[Green], PointSize[0.003]], 
     Directive[Red, PointSize[0.003]]}]]

As we can see the method fails here. There are two spiral arms; one starting from x_0 = -10.64 and the other one from x_0 = 10.64. I want the first one to be in red color and the second one in green like the following:

Surely you want to know how I obtained it. Well, I cheated! Since the two spirals are symmetrical I integrated the initial conditions of only the one spiral. Then I generated the symmetrical initial conditions for the other spiral arm and I merged the two plots. However, this is NOT a solution. In the case that the two arms are not symmetrical the cheat does not work, so the question remains.

Any suggestions?

Many thanks in advance!

Answer 1

One can use the following simple approach. For each radius $r$ one can decompose the density $\rho(r,\theta)$ to waves $e^{im\theta}$. If $m$ is equal to the number of arms we obtain the phase of the density wave, i.e. the spiral function $\varphi(r)$. Then we can use this phase to separate all stars to arms with a certain (I hope, acceptable) accuracy.

arms = 2;
data = Developer`ToPackedArray@Import["spirals.dat"]; 
r = Sqrt@Total[data^2, {2}];
(data = data[[#]]; r = r[[#]] ) &@Ordering@r; 
z = (data.{1, I})^arms;
φ = Arg@GaussianFilter[z, 200]; 
φ[[2 ;;]] += Accumulate@Round[Most@φ - Rest@φ, 2 π];
smoothR = GaussianFilter[
    Join[2 First@r - Most@Reverse@r, r, 2 Last@r - Rest@Reverse@r], 
    10][[Length@r ;; 2 Length@r - 1]];
f = Interpolation@Transpose@{smoothR, φ/arms};
arm = 1 + Floor[Mod[ArcTan @@ Transpose@data - f@r + π (2 n + 1)/arms, 
      2 π] arms/2/π];
split = Table[Pick[data, arm, n], {n, arms}];

Results:

Phase

Plot[f[r1], {r1, Min@r, Max@r}, AxesLabel -> {"r", "φ"}]

Spiral

Show[ListPlot[data, AspectRatio -> Automatic], 
 ParametricPlot[Evaluate@Table[{r1 Cos[f@r1 + 2 π n/arms], 
     r1 Sin[f@r1 + 2 π n/arms]}, {n, arms}], {r1, Min@r, Max@r}, 
  PlotStyle -> Red]]

Separate arms

ListPlot[split, AspectRatio -> Automatic]

There is a small difference with the expected result. However there is no straightforward method in a general case without additional information. Also one can use higher angular harmonics to enhance the result.