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!

Well, you need a way to tell in which arm lays each star ... It isn't clear enough for me that such a thing could be done for your "galaxy" — Dr. belisarius, Oct 30 '15 at 13:20
How would you decide which arm a point on their intersection belongs to? There's no clear spatial separation between the arms this time, so we can't use spatial clustering. Do you have other attached data that could differentiate them, e.g. velocity or similar? If not, let's think about why it naively seems like it is possible to do this at all. I think that's because the spiral is a simple geometric form that we, as humans, are familiar with and can manipulate. Any manual classification of the points is going to rely on your a priori knowledge of what a spiral is. — Szabolcs, Oct 30 '15 at 13:21
Maybe we should think along these lines and try to fit the shape to some model of a spiral to be able to discriminate between the two arms. Incorporate the a priori knowledge that we are working with spirals into the classification procedure. — Szabolcs, Oct 30 '15 at 13:22
@Szabolcs Currently the data file contains the (x,y) position of the stars at every time. But I could also print the momenta p_x and p_y at every moment. Would this help? — Vaggelis_Z, Oct 30 '15 at 13:24
I don't know. You're the astrophysicist :) Is there a clearer difference in the momenta of the stars belonging to the two arms (when they are near the same position)? — Szabolcs, Oct 30 '15 at 13:25
@Szabolcs I am afraid there is no difference. If somehow we could remove one arm form the data file then we could plot the remaining one with a color and then generate the mirror-symmetrical one with another color. Is this a good idea? — Vaggelis_Z, Oct 30 '15 at 13:29
@Vaggelis_Z It gives me another idea: if two points are each other's mirror image relative to reflection through the origin, then they must be in different clusters. At the same time, try to put nearby points in the same cluster. Not yet sure how to implement this in practice though ... — Szabolcs, Oct 30 '15 at 13:37
sorted=SortBy[Norm][data], then look at ListPlot[{sorted[[1;;-1;;2]], sorted[[2;;-1;;2]]}]. Take each in smallish chunks and find clusters. — LLlAMnYP, Oct 30 '15 at 13:55
Momenta would definitely help. For points where {x,p} are the same, choose color randomly. — Vitaliy Kaurov, Oct 30 '15 at 13:55
@LLlAMnYP I think something is missing in sorted=SortBy[Norm][data]. SortBy needs two arguments. — Vaggelis_Z, Oct 30 '15 at 13:59
Huh? SortBy should work with operator form, but if that's the case, then use SortBy[data, Norm] — LLlAMnYP, Oct 30 '15 at 14:02
@Vaggelis_Z the point is, that then every two successive datapoints in the list will be of the same norm and mirror images of each other. Then it's possible to build up clusters point-by-point. It's a nice problem, although it is quite localized to your very nice dataset where every point has a mirror image. I'm afrai I have to go, otherwise I'd have taken a shot at this. — LLlAMnYP, Oct 30 '15 at 14:04
What I said about the symmetry is not going to work. I know because it turned out to be very easy to implement: g = NearestNeighborGraph[data, 10]; {a, b} = FindGraphPartition[g, 2]; ListPlot[{a, b}] — Szabolcs, Oct 30 '15 at 14:53
@Szabolcs I have v9 and NearestNeighborGraph is not recognized... — Vaggelis_Z, Oct 30 '15 at 15:01
It doesn't matter as it won't give a good result anyway. But if you have v9, you should mention that in the question. Mma came a long way since v9. — Szabolcs, Oct 30 '15 at 15:06
I wonder if some of the answers to this previous question may be useful. — , Oct 31 '15 at 07:11

score 10 · Accepted Answer · answered Oct 31 '15 at 00:50

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.

Recently a gave a talk about two-stream instability in semiconductors. When I prepared this talk I found that this instability can occur in counter-rotating galaxies (e.g. M64) with arm formation (1, 2). It is quite interesting and can be related. — ybeltukov, Oct 31 '15 at 01:04

Coloring spiral structures

1 Answers1

Results:

Phase

Spiral

Separate arms