How to detect loops in a vector field?

Question

Given a vector field and I'd like to detect if there are streamlines in closed patterns. For example, if the input is a vector field like this:

image = ColorConvert[Rasterize[\[Infinity], ImageSize -> 40],"Grayscale"]
{dx, dy} = Table[ImageAdjust@GaussianFilter[image, {6, 2}, \[Delta]],
  {\[Delta], {{1, 0}, {0, 1}}}];
data = Reverse /@ Transpose[2 ImageData /@ {dx, dy} - 1, {3, 2, 1}];
ListVectorPlot[data, VectorPoints -> 30]

then I'd like to find these three streamline curves that I've drawn in red:

Detecting closed currents in a vector field seems very doable. I'm not sure if I should use autocorrelation.

Related: - Detecting edge/region in contour plot - How to find the center of a circular pattern? - Recovering data points from an image - Output is Input for a Differential with Sign?

Answer 1

Maybe you can transform this into an easier problem: If you were instead looking for lines that are at any point perpendicular to your vector field:

(that's easy, it's just a coordinate transformation)

and if your (transformed) vector field were the gradient of a scalar field (yes, big if), then this would be really easy: These lines are simply the level lines of the scalar field.

So I would try to "guess" a scalar field so that it's gradient is close to your original vector field (properly transformed). Using your code (I've just added padding and made the field a little larger):

image = ColorConvert[
  ImagePad[Rasterize[\[Infinity], ImageSize -> 80], 10, White], 
  "Grayscale"]
{dx, dy} = 
  Table[ImageAdjust@
    GaussianFilter[image, 
     2 {6, 2}, \[Delta]], {\[Delta], {{1, 0}, {0, 1}}}];
data = Reverse /@ Transpose[2 ImageData /@ {dx, dy} - 1, {3, 2, 1}];

I would "rotate" the vector field:

{dx, dy} = Transpose[data, {2, 3, 1}];
{rdx, rdy} = {{0, 1}, {-1, 0}}.{dx, dy};

and then as a first approximation I would simply accumulate the values:

sum = Accumulate[rdx] + Accumulate /@ rdy;
ListContourPlot[sum\[Transpose], PlotRange -> All]

(This is the tricky part: the "rotated" vector field doesn't actually have to be the gradient of this field, as long as the level lines are close to the lines you're looking for. I have no idea if this is possible for your data, and I doubt that the crude Accumulate approach is optimal. Maybe deconvolution can get better results - at least in signal processing theory, it should be the optimal way to "invert" a gradient filter.)

Then I would look for the level lines using image processing:

threshold = 0.5;
img = Image[Rescale@sum];
perimeter = 
 MorphologicalPerimeter[ColorNegate@Binarize[img, threshold], 
  CornerNeighbors -> False]

(where threshold would take values from 0..1)

By finding connected components, I can get (approximations to) the curves you're looking for:

comp = MorphologicalComponents[perimeter];
compPaths = Table[
   With[{pos = Position[comp, i]},
    pos[[#]] & /@ FindCurvePath[pos][[1]]], {i, Max[comp]}];

The idea is, that if there are closed streamlines, these level lines should be close. To check if they actually are closed streamlines, you would have to integrate your original vector field along these level lines:

Function[path,
  {pdx, pdy} = Extract[#, path] & /@ {dx, dy};
  {dirY, dirX} = ListConvolve[{-1, 1}, #, 1] & /@ Transpose[path];
  dirX.pdx + dirY.pdy] /@ compPaths

(which gives a small, but in this case nonzero values - probably because the "rotated vector field is only approximately the gradient of sum)

ListVectorPlot[
 Transpose[{dx, dy}, {3, 1, 2}], VectorPoints -> 30,
 Epilog -> {Red, Line[compPaths]}]

The idea would be to "improve" the paths found this way using e.g. numerical optimization. You now have a starting point, so you can use FindMinimum magic.

Defining some interpolation and starting with a smoothed version of the path:

intDx = ListInterpolation[dx];
intDy = ListInterpolation[dy];

smoothPaths = 
  Transpose[
     Map[GaussianFilter[#, 5, Padding -> "Periodic"] &, 
      Transpose[#]]] & /@ compPaths;
path = smoothPaths[[1, ;; ;; 2]];

I can then define an optimized path with a small offset from this one:

pathTangent = Normalize /@ (path - RotateLeft[path]);
pathNormal = pathTangent.{{0, 1}, {-1, 0}};

optimizePath = 
  Table[path[[i]] + pathNormal[[i]]*ofs[i], {i, Length[path]}];

optimizePathTangent = (optimizePath - RotateLeft[optimizePath]);

...and an optimization objective that minimizes the angle between the path tangent and the vector field:

pathCost = Total[
   Table[
    ({intDx @@ optimizePath[[i]], 
        intDy @@ optimizePath[[i]]}.{{0, 1}, {-1, 
         0}}.optimizePathTangent[[i]])^2
    , {i, Length[path]}]];

and optimize that:

{err, sol} = 
  FindMinimum[pathCost, Table[{ofs[i], 0}, {i, Length[path]}], 
   Method -> "LevenbergMarquardt"];

ListVectorPlot[Transpose[{dx, dy}, {3, 1, 2}], VectorPoints -> 30,
 Epilog -> {Red, Line[optimizePath /. sol]}]

This is not as smooth as I would have expected, but according to the residual error, it follows the streamlines closely.

That's all I've come up with so far. I'm not sure it's a good start, but it might point in the right direction.

Answer 2

The following is far from perfect, just a kickstart:

image = ColorConvert[Rasterize[8, ImageSize -> 40], "Grayscale"]
{dx, dy} =  Table[ImageAdjust@ GaussianFilter[ image, {6, 3}, δ], {δ, {{1, 0}, {0, 1}}}];
data = Reverse /@ Transpose[2 ImageData /@ {dx, dy} - 1, {3, 2, 1}];

p = ListStreamPlot[data, StreamPoints -> {400, .2, 10}, Frame -> False];
arw = Cases[p, {{Arrowheads[__], ___}, ___}, Infinity];
j = Join @@ (Cases[#, Arrow[a : __] :> a, Infinity]) & /@ arw;
lines = Select[j, Equal @@ (FindCurvePath[#][[1, {1, -1}]]) &];
Show[p, Graphics[{Red, Line[{##}] & @@@ lines}]]

Result in 10.0.2 (Thanks @MichaelE2) and V9