14
$\begingroup$

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]

enter image description here

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

enter image description here

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?

$\endgroup$
5
  • 1
    $\begingroup$ Do you have any particular application for this? Usually, whether flow lines are closed or not isn't really of practical importance, but one can sometimes make up special interpretations for them, see here. $\endgroup$
    – Jens
    Aug 27, 2015 at 4:10
  • $\begingroup$ The application here is that detection of loops in the phase space means something about the state of the physical system. $\endgroup$
    – M.R.
    Aug 27, 2015 at 4:36
  • 2
    $\begingroup$ But the example you're using is not a physical system, and is in fact just not even a vector field because it is a discrete set of data points. So We first have to construct an interpolation. Are there any physical properties of the vector field that have to be satisfied in the interpolation? Whether lines will be closed for a given starting point will in general depend sensitively on the interpolation. $\endgroup$
    – Jens
    Aug 27, 2015 at 5:04
  • $\begingroup$ @Jens: Sure, "contains loops" would have to be interpreted as e.g. "if you follow the vector field, you come back to your point of origin again, up to some numerical error epsilon". But isn't that true of practically any boolean property when you're doing numerical calculations? $\endgroup$ Aug 27, 2015 at 15:12
  • $\begingroup$ @nikie I was just trying to get as much info as possible before going in the wrong direction. I see you've already answered based on an assumption such as what I was asking for... There could be different assumptions - it could bee the curl of a vector, or you could be thinking of limit cycles in a nonlinear system, etc. $\endgroup$
    – Jens
    Aug 27, 2015 at 15:32

2 Answers 2

20
$\begingroup$

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:

enter image description here

(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]

enter image description here

(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)

enter image description here

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]}]

enter image description here

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]}]

enter image description here

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.

$\endgroup$
2
  • $\begingroup$ Amazing effort! $\endgroup$ Aug 27, 2015 at 11:02
  • $\begingroup$ Very nice indeed. $\endgroup$ Aug 27, 2015 at 14:42
10
$\begingroup$

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

Mathematica graphics

$\endgroup$
2
  • $\begingroup$ The result is slightly different (better, imo) in V10.2. $\endgroup$
    – Michael E2
    Aug 27, 2015 at 11:24
  • $\begingroup$ @MichaelE2 Thanks a lot, added $\endgroup$ Aug 27, 2015 at 14:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.