Finding filaments in an image I am at the moment working on images such as this one:



What you see are filamentous structures / bundles. Other images coming from slightly different experiments could have more sparse / thick / bright filaments, but the order-of-magnitude and qualitative average image is like the one above.

What I would like to do is to extract as many information as possible from this image and the best thing would be to identify each single filament/bundle or a "skeletonization" of the image, but I am happy already with a vector field indicating the local orientation of the filament so that I can compute its divergence etc. etc.

So the question is, do you know of any program / conceptual tool / framework which may help me identyfyinf the filaments and or their local orientation?

I usually work in Python and any example-code is welcomed, but I would be more than satisfied with a purely mathematical / pseudocode / conceptual answer.

What I did so far:


I tried to use PCA (Principal component analysis). I select a ROI, threshold the image, create the dataset $(X, Y)$ of the pixels with high intensity and compute the covariance matrix using the intensity as a weight. Then I diagonalise it finding the "longer" eigenvecotr and assume that that is the local orientation. Works ok-ish, but I am afraid is very paramtere-dependant and may fail when filaments become too thick or too thin. An example is below. It looks good but I fear is too simplistic an approach.
I tried to implement (in Python) the procedure they describe here. Maybe I am tweaking the parameters wrong or is supposed to work only for very thin filaments, but it works very badly in my case.




I know there are pre-made tools online, but I'd rather write everything from scratch.

So summing up, the big question here is:

What is the best way to find the filaments and/or their orientation at a given point?

Thank you all!
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.