There seem to be a number of hits in the literature from other domains that use various methods to segment the image (areas likely to belong to a filament and 'background') followed by filament identification / extraction which you might want to have a look at anyway.
I was intrigued by the choice of PCA as a first indicator of extracting the local gradient and so this is a response with the "low hanging fruit" you can perhaps try (and in fact extremely easily) before jumping to the big guns of pattern recognition.
...the big question here is:
What is the best way to find the filaments and/or their orientation at a given point?
Unless you have prior information (or assumptions) about the direction of the flow / filaments, it is impossible to derive their orientation. In other words, you can get some indication about the local orientation of (something that looks like) an edge, but not the orientation of the flow of the filament. You can assume, for example, that there is a particular intensity profile along the body of the filament which could then be used to determine "which part is the 'head' ".
The first thing to do (anyway) is to obtain the image gradient. The image gradient will give you the local slope of the intensity field (a.k.a image). All code here in Octave but only for demonstration. The operators are easily transferable across to scipy.
%Load the image
Q = imread('someImage.jpg');
%Resize it to 1/4 of its original scale, this is not essential but it helps with visualisation.
Qs = imresize(Q,0.25)
%
% Get the gradient
% [gradMag, gradDir] = imgradient(Qs);
% gradMag is the Magnitude, gradDir is the direction in degrees
% Break down gradDir to its components to process and visualise it
X = cos((gradDir/180.0).*pi);
Y = sin((gradDir/180.0).*pi);
%Let's have a look at the image with its gradient superimposed
imshow(Q2);hold on;quiver(X,Y);
The gradient of the (rescaled) image looks like this:
And if we zoom in (to the plot, not the image) around a bundle, we can see something like this:
So, predictably enough, the arrows are pointing away from the top of a ridge that signifies a filament and this can be used, to an extent, to derive its local orientation.
BUT, the bundle's intensity is not constant and the local gradient is not informed by any prior so that it "knows" which way to point.
Now, once you have moved from a scalar (Qs) to a vector field (X,Y), you might find other vector field operators useful as well, depending on what you are trying to do.
Can we do better?
Yes, we can apply a transformation to the intensity values which dampens the soft values and boosts the strong values. This will 'intensify' the ridges and basically result in a more well defined gradient later on. This is a modification of the image's dynamic range and is also known as 'expanding'.
An exponent is a good choice, it effects the following mapping:
Where we simply raise the pixel value to a power and take care of normalisation so that it does not overshoot the typical image limits.
This transforms the original (scaled) image to something like:
And conversely the local gradient now will be more well defined.
Now, if you are trying to detect broader structures which you would then infer (with some certainty) that they are formed by or belong to filaments, then you can do a couple of things. First of all, you could use Curl to identify "turbulent" (or absolutely no turbulence) areas and from there try to make sense of 'circularity' and then bundle together those areas.
The other thing you can do is something along the lines of line or circle detection through the Hough Transform. The Hough Transform is a "variant" of the Radon Transform. The latter is transforming an image from a $x,y$ plane to a $r, \theta$ plane by obtaining line integrals (the magnitude of which is $r$) at angle $\theta$ around the image. This looks like this (on Qs):
[R,xp] = radon(Qs);
imshow(R./max(max(R)));
The Radon transform converts lines (in Qs) to points. A line will produce a high value integral when $\theta$ is parallel to the bearing of the line. Conversely, arcs (and circles) are converted to arcs in the $r, \theta$ plane and if you follow intense arcs along the $\theta$ direction (corresponding to $x$ in a typical image) you will be able to identify the curled ridges that are formed by the filaments.
This is used by the Hough Transform which essentially operates on a thresholded version of the initial image (e.g. Qs) to detect lines and circles by looking for points with specific spatial relationships in the $r, \theta$ plane.
Hope this helps.
