How to determine orientation distribution in a fiber network ?

I am processing SEM images of cellulose network. I am not sure as to how to go about determining the orientation distribution of the fiber network. The algorithm that I thought of is:

1. Generate random pixels
2. If the value at that pixel is '1' I have landed inside the fiber network. Draw arbitrary lines through those pixels (at angle step of ,say, 5 or 10 degrees)
3. The shortest of all those distances where value changes from '1' to '0',in both the directions, would give radial direction at that point
4. At any pixel the orientation of the fiber would be perpendicular to the direction that gives me the shortest distance

An algorithm for diameter distribution was suggested previously but that approach did not take care of the spatial orientation of the fiber network.

https://www.mathworks.com/matlabcentral/answers/108856-how-to-draw-line-with-arbitrary-orientation-through-a-random-pixel-and-calculate-the-distance-along

A sample image is given below:

How should I go for the implementation my algorithm? Any alternate algorithm suggestion would be appreciated.

Regards,

Rajesh

[Walter Roberson]

Have you considered skeletonization?

http://www.mathworks.com/matlabcentral/answers/88284

http://www.mathworks.com/matlabcentral/fileexchange/25865-euclidean-skeleton

Contour-Pruned Skeletonization at http://cs.smith.edu/~nhowe/research/code/

[Rajesh]

Yes I did skeletonize the image and then used an ImageJ plugin that allowed to analyze the skeleton image. I obtained the avg. fiber length, the start and end point coordinates of each fiber, the triple points, quadruple points, end points etc. But I need to construct a computer model that can be read by COMSOL through the COMSOL/MATLAB interface Livelink. So I needed to implement this in matlab. How could that be acheived?

[Walter Roberson]

Once you have skeletonized in MATLAB, you can find the branch points, and the connectivity information. Then from two (x,y) pairs, you can compute the angle -- it is atan2( (y2 - y1), (x2 - x1) ).

[Rajesh]

@Walter: Thanks!! I get your point. But could you please clarify this:

1. How do I select the branch points of interest (meaning to say that I cannot pick up any two branch points, they need to be connected)? How do I ensure this? More over this process has to be automated over the entire skeleton.

2. Further more the fiber may not be a straight one. So while calculating the length, I need to traverse along the curve joining two branch points in the fiber network and not just take the euclidean distance. How do take care of this?

Answer 1

[Image Analyst]

Use imgradient() to get the gradient direction. Multiply by the skeleton to get the gradient direction only at fiber centerlines. Then take the histogram to get the distribution of directions.

[Rajesh]

By the definition of gradient direction it should be radially outward at the edges for a binary image. I am a bit confused as to what gradient direction would signify in the fiber region towards the center (for a binary image), given that the neighboring pixels have the same intensity. Could you please clarify this?

[Image Analyst]

imgradient() will give you the direction of fastest slope. So I was thinking that if you take the gradient of the euclidean distance transform , at the centerlines of the fibers, the slope would point towards the edge, right? You want the distribution of just the gradient along the centerlines so we need to extract out just the gradients there. So if we multiply by the skeleton we'll get zero everywhere except at the centerlines where we'll have the directions of the centerline. Follow me so far? So then if we take the histogram of that then we'll have the distribution I think you're looking for. If you want the direction along the fiber instead of perpendicular to it, then that's just 90 degrees added to the histogram. Does it make more sense now. I haven't done this so if I'm totally wrong, let me know.

[Rajesh]

I took a slightly different approach. I detected the edge of the fiber network. Took the inverted image of the one cited in the thread. Found its imgradient(). This gives me the gradient direction all along the edge (radially outward). Now I multiplied the two. This gives me the direction distribution of the network along the edges. Adding up 90 degrees would give the orientation along the length of the fiber. This approach would take care of the fact if the fiber has a different orientation/behavior on either edge. Correct me if I am wrong!!