Spacing out random walks so they don't overlap

23

5

From an external simulation program, I have lots of particle tracks which follow random walks, all of which start at or near the origin $(0,0)$. This means that when I plot the walks in Mathematica, it looks a bit of a mess because all the tracks overlap with each other.

For example:

randomwalks = Table[
    Accumulate@RandomVariate[StableDistribution[1, 1.7, 0, 0, 2], {512, 2}], 
    {i, 3}];

ListLinePlot[randomwalks, PlotRange -> Full, PlotTheme -> "Minimal", 
    Axes -> None, Frame -> True, AspectRatio -> 1];

This gives the image on the left - the fact that all the walks start near the origin means that they overlap and it's hard to see what is going on.

enter image description here

What I'd like to plot instead is something like the image on the right, where the random walks have been manually translated to stop them overlapping.

Is there a way to achieve this translation automatically with Mathematica?

The direction of each jump in the random walk is important, so I want to keep this in the visualization (hence translations only and no rotations).

Note also the choice of a Lévy-stable distribution for the random walk rather than the normal distribution. This typically gives walks that are more spread-out with longer jumps (due to the fat distribution tails). In turn this leads to a more diverse spread of shapes to try and fit together, and means arranging the shapes in a simple grid (e.g. via Column[Graphics /@ Line /@ randomwalks] is not always as successful as I'd like.


My first thought was to calculate the bounding rectangle for each track, and then use the code in An algorithm to space out overlapping rectangles? to space out those bounding rectangles. However, this might leave large regions of whitespace due to the nature of the random walks, so perhaps a bounding polygon might be better to give a tighter arrangement of the tracks.

Given I have the track data itself, is there perhaps a way to do it without using the bounding rectangles? As suggested in a comment, perhaps using the convex hull of each path and then packing the resulting polygons?

regions = ConvexHullMesh[#] &/@ randomwalks
centroids = RegionCentroid[#] &/@ regions

Once we have the centroids, one approach would be to minimize the Euclidean distance between each centroid whilst ensuring the regions don't overlap - but that would need a function such as RegionOverlapQ[] or similar.

dr.blochwave

Posted 2015-09-10T09:37:38.847

Reputation: 8 483

Maybe it could be done with the help of rasterization and procedures described in 2334 and 88846

– shrx – 2015-09-10T11:23:04.440

5One idea: Take the convex hull of each path, and then your problem is reduced to packing convex polygons. – Joseph O'Rourke – 2015-09-10T11:23:59.637

@JosephO'Rourke "convex hull" - that's the term I was looking for! Nice idea. – dr.blochwave – 2015-09-10T11:30:22.867

5Caveat: Optimal packing of convex polygons into a rectangular box is NP-hard, so you will be forced to use heuristics. (Mma has ConvexHull[ ].) – Joseph O'Rourke – 2015-09-10T11:34:49.580

@JosephO'Rourke sure - it needn't be optimal, just visually pleasing! – dr.blochwave – 2015-09-10T11:39:10.630

2Why not just Column[Graphics /@ Line /@ randomwalks] or something? – Kuba – 2015-09-10T11:41:49.803

2@Kuba because I fancied a challenge? :-) – dr.blochwave – 2015-09-10T11:43:09.887

1why only translations and no rotations? – Nikolay Gromov – 2015-09-10T13:42:04.610

@NikolayGromov because in the output from my simulation program, the direction of each jump is important and I want to keep that in the visualization (good point though, will clarify in the question) – dr.blochwave – 2015-09-10T13:46:23.963

2OMG My code there is mess – Dr. belisarius – 2015-09-10T14:21:56.703

1

One approximative approach to this is related to tiling a region with a set of polynominoes approximating individual random walks. This can be solved, for example, using binary integer linear programming. Time to find instances of these solutions is not necessarily nice. I have a crude solution to this problem, but code is truly awful at this point...

– kirma – 2015-09-10T17:57:13.113

Answers

13

One approach is to rasterize the lines separately and use correlation to measure overlap. The lines are positioned one by one, as close to the centre as possible without overlapping any previously positioned lines. Heike used this in her answer to the word cloud question and I borrowed the idea here.

Here's a quick & dirty implementation:

randomwalks = 
  Table[Accumulate@
    RandomVariate[StableDistribution[1, 1.7, 0, 0, 2], {512, 2}], {i, 15}];

r = 2.05 Max@Abs@randomwalks;

bin[walk_] := Binarize@ColorNegate@Erosion[
    Rasterize[Graphics[Line@walk, PlotRange -> r], ImageSize -> 500], 
    DiskMatrix[3]]

radial = RadialGradientImage[Black, {500, 500}];

iter[im1_, im2_] := Module[{corr, minx, miny},
  corr = ImageAdd[radial, ImageCorrelate[im1, im2]];
  {minx, miny} = 250 - PixelValuePositions[corr, Min@ImageData@corr][[1]];
  Sow[-{minx, miny} r/250];
  ImageAdd[im1, Image@RotateRight[ImageData@im2, {miny, -minx}]]]

offsets = Reap[Fold[iter, bin /@ randomwalks]][[-1, 1]] ~Prepend~ {0, 0};

ListLinePlot[
 Transpose[(Transpose[randomwalks, {1, 3, 2}] + offsets), {1, 3, 2}], 
 PlotRange -> Full, PlotTheme -> "Minimal", Axes -> None, 
 Frame -> True, AspectRatio -> 1]

enter image description here

Simon Woods

Posted 2015-09-10T09:37:38.847

Reputation: 81 905

1Very much different from mine, but I am confident there are many different approaches on this problem (as already outlined on the question). – kirma – 2015-09-12T15:41:14.640

This one has ended up being more appropriate for my needs, although @kirma's answer is a very interesting idea nonetheless! – dr.blochwave – 2015-09-21T12:29:22.370

11

This solution is primarily a proof of concept; its run-time complexity makes it impractical for large amounts of items. Large amount means more than six, in this case.

Nonetheless, the idea is the following: individual plots are converted to polynominoes to be laid out on a rectangular grid. Every possible alignment on a finite grid is encoded to a Boolean equation representing legal layouts (every piece is present on exactly one position, and the pieces don't overlap), and a satisfying layout instance is searched for this equation. Result is baked back to datasets as offsets, which can be plotted now without overlaps.

ClearAll[layOutDatasets];

layOutDatasets[datasets_List, numinstances_Integer, 
   resolution_Integer, maxdim_List] :=
  Module[
   {normalizeDatasets, generateTile, generateAlignments, 
    generateOrientations, atMostOne, exactlyOne},

   normalizeDatasets[
     dss_List] := ((With[{offset = 
            First@CoordinateBoundingBox@#}, # - offset & /@ #]) & /@ 
       dss)/Max[Max@Abs[Subtract @@ CoordinateBoundingBox@#] & /@ dss];

   generateTile[ds_List] :=
    ImageData[
        ImageCrop@
         Binarize[ImageResize[#, resolution + 2], {0.99, 1}]] /. {0 ->
          True, 1 -> False} &@
     ImageReflect@
      Rasterize[
       Style[Graphics[{Thickness[1/(2 resolution)], Line@ds}, 
         PlotRange -> {{-1/resolution, 
            1 + 1/resolution}, {-1/resolution, 1 + 1/resolution}}], 
        Antialiasing -> False], ImageSize -> 8 (resolution + 2)];

   generateAlignments[patt_List] :=
    With[{maxoff = maxdim - Dimensions@patt},
     Flatten[Table[{i, j}, {i, 0, maxoff[[1]]}, {j, 0, maxoff[[2]]}], 
      1]];

   generateOrientations[patt_List] :=
    PadRight[PadLeft[patt, Dimensions[patt] + #, False], maxdim, 
       False] & /@ generateAlignments[patt];

   atMostOne[v_List] := BooleanCountingFunction[{0, 1}, Length@v] @@ v;
   exactlyOne[v_List] := BooleanCountingFunction[{1}, Length@v] @@ v;

   Module[{normdatasets, pieces, placemasks, sols, offsols},
    normdatasets = normalizeDatasets[datasets];
    pieces = generateTile /@ normdatasets;

    placemasks = generateOrientations[#] & /@ pieces;

    With[{vars = 
       Flatten@MapIndexed[c@First@#2 /@ Range@Length@#1 &, placemasks],
      varmaps = 
       Flatten@MapIndexed[
         c[#2[[1]]][#2[[2]]] -> Reverse@#1/resolution &, 
         generateAlignments[#] & /@ pieces, {2}]},

     sols = 
      SatisfiabilityInstances[
       And @@ (atMostOne /@ 
           Flatten[
            Transpose[
             Flatten[
              MapIndexed[#1 && c[#2[[1]]][#2[[2]]] &, 
               placemasks, {4}], 1], {3, 1, 2}], 1]) &&

        And @@ MapIndexed[
          With[{constants = c@First@#2 /@ Range@Length@#1},
            exactlyOne@constants] &, placemasks],
       vars, numinstances];

     offsols = (Pick[vars, #, True] /. varmaps) & /@ sols;

     MapThread[
        Function[{set, off}, # + off & /@ set], {normdatasets, #}] & /@
       offsols]]];

An example of use:

ListLinePlot[#, Axes -> False, AspectRatio -> Automatic] &@
 First@layOutDatasets[
   Table[Accumulate@
     RandomVariate[StableDistribution[1, 1.7, 0, 0, 2], {512, 2}], {i,
      5}], 1, 5, {15, 15}]

enter image description here

Another run:

enter image description here

EDIT: Note that there is no reason this method would be limited to a rectangular grid, regular alignments or even specific rotations. Generating more varied alternatives is always an option (although it increases size of the search space, which increases time complexity exponentially), and generating random alignments instead of regular lattice alignment is also an option. Hexagonal lattice would be quite interesting, but "rasterizing" graphics to such a grid is somewhat cumbersome.

kirma

Posted 2015-09-10T09:37:38.847

Reputation: 13 550

This looks really cool - I'll give it a go later on and report back properly! – dr.blochwave – 2015-09-12T11:29:55.053

@blochwave It's more than a bit unpolished and has hacks all around - but I decided to dump it here anyway, because polishing might have never happened... The core idea is pretty simple, though. – kirma – 2015-09-12T11:33:17.603

1you never know, community action might help polish it! – dr.blochwave – 2015-09-12T11:36:22.513