Image processing - diameter distribution of fibers
DavideT
Wed, 04/17/2013 - 11:48 pm
I have a bunch of SEM images of a network of fibers. I would like to get the distribution of diameters, but please note that fibers' diameter changes along the fibers' length - so even if there are 10 fibers in the image, the distribution should not be made from 10 points, rather from a large number of step by step detected diameters (or better, width) for each fiber.
I have played around with the image processing demo experiment and not found a suitable way to do that. I was thinking of trying a brute force approach, making the image binary, then scanning pixel by pixel to check if the actual pixel is at the border of a fiber, then take the previous and next 5 pixels on the border and make a line through them, then create the perpendicular line to this one and finally check when this line touches the other side of the fiber.
Maybe difficult to be done and definitely time consuming, scanning millions and millions of pixels in a bunch of images...
Do you think there is a more elegant and time saving way to do that?
Or is it possible to somehow implement this feature as a series of steps in the already existing Image Processing functions?
Or even better, if someone has done it already!
Let me add, this would be an extremely appreciated feature, as my alternative is to buy a very expensive stand-alone software that does just and only that.
Thanks in advance for your answers!
ImageAnalyzeParticles [ flags ] keyword imageMatrixoperation. In particular, the following key word looks interesting:
reported for all particles whose area exceeds the minArea specified by the
/A flag. The results of the measurements are:
V_NumParticles Number of particles that exceed minArea limit.
W_ImageObjArea Area (in pixels) for each particle.
W_ImageObjPerimeter Perimeter (in pixels) of each particle. The perimeter
calculation involves estimates for 45-degree edges resulting in
noninteger values.
This pre-supposes your fibers are isolated and can be treated as "particles". Masking might be needed at the edges of the large image wave. If this works, the mean fiber diameter would be its area divided by the perimeter.
April 18, 2013 at 06:10 am - Permalink
thanks for your answer. I think that does not work, because the fiber must not be considered as a single particle. I need the distribution of diameter along its length. To make an example, if I have one fiber in my image, I don't need 1 diameter - rather a full distribution centered at the most common diameter.
The approach you kindly suggested may work in case it would be possible to change the Wavemetrics code, detecting particles but then not excluding the area covered by this particle from next detection. And the particles, should be forced to be circles. Kind of filling the interior of the fiber with circles. Maybe :)
April 18, 2013 at 07:16 am - Permalink
--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAHuntsville
April 18, 2013 at 08:10 am - Permalink
Jeff Weimer is right- we need an example image.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
April 18, 2013 at 09:13 am - Permalink
Without a picture or more details it would seem to me as if you have an image of a fiber that has a varying diameter as a function of length and you want to map that diameter. If that is the case then a tricky way of getting a "measure" of the diameter is actually using ImageLineProfile. If you set the profile path along the center of the fiber, choose a width w such that w is just a bit larger than the largest diameter and use the /S flag you will get a distribution of the standard deviation of the intensity profile perpendicular to the axis of the fiber. Clearly this method will fail if you have significant background or foreground intensity fluctuations. In which case I'd approach the problem differently:
1. identify the path of the fiber using thresholding followed by morphological operations to get a path corresponding to the axis of the fiber.
2. for points along the axis of the fiber compute local derivative to find the local normal to the axis.
3. use ImageLineProfile for two points on both sides of the local normal.
4. extract the edges of the fiber from the path profile.
I hope this helps,
A.G.
WaveMetrics, Inc.
April 18, 2013 at 09:47 am - Permalink
if the thresholding is good enough ---i.e., if you don't gain much from determining the fibre edge locations from each line profile separately--- a Euclidean distance transform will give you the fiber diameter distribution pretty much directly. The histogram of the EDT is something like the integral of the desired fibre width distribution, because the EDT will have a triangular profile across the width of the fiber, reaching half the fibre width in the middle. There's no need to find the middle of the fibers and determine directions normal to the fiber axis, because the EDT does that. This method would be quite sensitive to outlier background pixels inside the fiber. Some morphological operations to get rid of those first might be wise. Also, getting the distribution from its integral is liable to be numerically unstable (If the radius of curvature of the centerline is smaller than the radius of the fiber, parts of the triangular profile go missing.), so this method will only work well for the distribution of the thickest fibers.
I've posted a function to calculate something close to the EDT in snippets. It's been a while since I used that code, and with those for loops it's probably not nearly as fast as an XOP or built-in function would be, but I seem to remember it worked more than well enough for my needs at the time.
April 18, 2013 at 07:21 pm - Permalink
I am sorry I messed it up with the description: my aim is to measure the diameter of fibers indeed, however in a SEM image these are the "widths of linear elements" as suggested above.
I attached a couple of examples. The diameter ermmm width seems to be the same along the fiber, however this is not a general example.
I will try to have a look at the ImageLineProfile, but my guts feeling is that the task is a little more complicated than what can be done with standard features. As for A.G. comment and his 4 suggested steps... "identify the path of the fiber using thresholding followed by morphological operations to get a path corresponding to the axis of the fiber" seems very difficult, please see the attached examples. I will try somehow, and sorry again if I wasted your time with a not too clear description of the problem. I hope the attachments can help.
Sbossuyt, I didn't really understand what the Euclidean Distance Transform does, so I will just try it and see if it works.
April 21, 2013 at 09:07 am - Permalink
Now that I see the images in question it is clear that what I suggested is not applicable. In fact, I wonder if any simple method is applicable because it appears that the fibers in the image are not necessarily parallel to the plane of the image complicated further by the fact that you have fibers intersecting each other. This makes me wonder if you really need to know the variation of the diameter along the length of individual fibers or if you need to determine the general distribution of diameters.
If you need to find the actual variation along each fiber then you first have to extract the relevant portion of a fiber from the image and only then apply one of the methods mentioned above. If you are only interested in the distribution of fiber diameters within the image I'd explore the possibility of obtaining that information from spatial transforms.
A.G.
WaveMetrics, Inc.
April 22, 2013 at 10:30 am - Permalink
* Draw two perfectly vertical lines at opposite ends of the image frame. They are at x positions of 0 and xI, where xI is the width (in x) of the image.
* Find the intercepts of those lines with the fibers. You will have four points ... LEFT (x = 0) : (yL1, yL2) and RIGHT (x = xI) : (yR1, yR2).
* The equation along the fiber center line is yCL = mCL x + bCL. Find it.
* Use geometry to find the equation for the perpendicular to the center line at the two ends
* Use geometry to get the fiber diameter at the two ends
You could do the same with perfectly horizontal lines instead when you have fibers that are nearly vertical.
This suggests to me that a function could be designed do all of the math when given the (six) coordinates xL, yL1, yL2, xR, yR1, yR2. All that would be left is a GUI to allow you to pick off coordinates on any given fiber from two vertical or two horizontal lines that intercept the fiber at some reasonable separation distance.
An alternative method that comes to mind is to draw straight lines along the two fiber edges, then find the angle between them using vector algebra (dot product). Use geometry to convert the change in angle to a change in diameter per length. With this, and given the diameter at one end of a fiber, you can calculate its diameter at the other end.
--
J. J. Weimer
Chemistry / Chemical & Materials Engineering, UAHuntsville
April 22, 2013 at 10:40 am - Permalink
Given a binary image, i.e., an image where every pixel is labeled as either foreground (1) or background (0), the Euclidean Distance Transform returns a grayscale image where every pixel value is the euclidean distance (i.e., straight-line distance) to the nearest background pixel. For the background pixels, that distance is 0, for pixels directly adjacent to a bacgkground pixel it is 1, for pixels sharing a corner but not an edge with a background pixel it is sqrt(2), etc...
The algorithm actually calculates the distance squared, which is separable according to Pythagoras' theorem, and then takes the square root of that at the end. The algorithm I implemented isn't exact, but it has the benefit that it's easy to follow what's going on: it scans across the image left-to-right and top-to-bottom, and then right-to left and bottom-to-top, keeping a record of the nearest background pixel of each pixel and checking nearest background pixels to the nearest neighbors of the current pixel as possibly nearer background pixels to the current pixel than what it has recorded until then. It should be fairly straightforward to modify the code so that it creates an extra wave flagging where the nearest background pixel is very different to that of a neighbor, for example storing a copy of the distance value from the EDT for those pixels. The binary image indicating the location of those pixels is called the skeletonization of the binary image, and includes the centerlines of the fibers. It also includes side branches where the fiber changes direction. (Interestingly, the combination of the skeletonization and the EDT values for those pixels is enough information to reproduce the original binary image.) The distance values will be larger than the fiber diameter where fibers overlap in the image (because the nearest background pixels will be the corner pixels of where the fibers overlap, rather than the edge pixels of each fiber), but where the fibers do not overlap, the value of the EDT on the centerline is half the diameter of the fiber.
It will likely be non-trivial to deal with the cases where fibers overlap in the image in any sort of automated way. Also, have you thought about the need to avoid sampling bias in the way you measure the fiber width distribution? The number of width values you take from each fiber to build up the distribution should reflect how much you want that fiber to count towards the distribution, rather than how many width values you could get out of the image for that fiber! The fiber crossings will affect how many width values you can get out of the image for a given fiber, for example. Depending on what you want, you may
1) just manually identify the ridges in the EDT image and store a handful of values from that, or
2) you may want to mask out the fiber crossings and side branches from the skeletonization, but otherwise gather statistics on all the EDT values from the pixels picked out by the skeletonization, or
3) mask each fiber individually, and get separate statistics for the width distribution of that fiber along its length from that, or
4) maybe something else still.
On second thought, I just had an idea for what might be a good way to deal with the cases where the fibers cross in the image: if you decide to modify the EDT code I posted, it should be fairly straightforward to identify when the new nearest-background-pixel is (close enough to) diametrically opposed from the old one, and only store a copy of the distance for those pixels instead of doing the skeletonization properly. Rx and Ry in the code give the relative position to the nearest-background pixel from the current pixel. At the very end, of the second pass, you find out what direction that is, take the nearest neighbor pixel to the current pixel in the opposite direction, and if that has Rx and Ry values close to the opposite of the current pixel, you keep a copy of its distance value. If you care about the orientation of the fiber, you might keep a copy of Rx and Ry instead. I don't have time to try this for myself, but I'm curious how well it works. Looking at your example images, I think it should work reasonably well if your method of classifying pixels as foreground or background pixels works well enough. The methods does still introduce sampling bias in the width distribution of the fibers, however.
Kind regards,
-- Sven Bossuyt
April 24, 2013 at 04:42 am - Permalink
thanks for your comments, I have some more ideas now how to approach the thingy - and some interesting steps forward. I will just be stubborn and keep on trying. After all, I realized that the task is difficult and also the alternative software does not really give "the" distribution, rather a good approximation. Hopefully I'll manage to get a rather good aprroximation myself. I could buy tens of Igor licenses (per year!) for all my colleagues, with the money saved from that other software :)
I'll keep you informed of the progress.
April 25, 2013 at 01:15 am - Permalink
here I am again... we switched to a FEG-SEM and used the back scattered electrons detector, so the contrast is better.
Then I wrote a few lines of code that take the image, make it binary with a certain threshold (I used this, so I was able to delete the SEM info without using ROI etc. In the future I will identify the threshold from the histogram. At the moment I pass 120 and it uses it for the threshold), detect the edges of the fibers, and find the centers.
Finding the centers with my own code kind of works, it just "overdetects" at the edges but I am relatively confident I can fix that, adding an if-endif. [corrected editing this post, it doesn't anymore]
I attached a test wave that is one example image, if you run this code you will get three images indeed (binary, edges, edges+centers).
It looks like it's time to use your previous suggestions with line profile or euclidean distance transform.
I thought of posting this, both such that other users can use something similar if needed, and also to possibly get some more ideas from you looking at these examples.
A cool thing would be to get thinner edges.
And... I hope you will be able to load and run this function, working with images and their .tif at the end into the name of the waves gives me headaches... just in case I did something wrong, I also attached the jpeg of how these images look like.
Thanks in advance for any idea!
variable threshold
wave 'area1_001b.tif' //at the moment, it works on this specific image/wave
duplicate/o 'area1_001b.tif', thisone
variable c1,c2,c3,c4,c5 //counters
variable columns=1024
variable rows=1024 //one is wrong but doesn't matter for now
for(c1=0;c1<columns;c1+=1)
for(c2=0;c2<rows;c2+=1)
if(c2<58) // 58 to get the SEM info away
thisone[c1][c2]=0
endif
if(thisone[c1][c2]<=threshold)
thisone[c1][c2]=0
endif
if(thisone[c1][c2]>threshold)
thisone[c1][c2]=254
endif
endfor
endfor
Display/k=1;AppendImage thisone
//MatrixFilter/N=(3)/P=1 median thisone // median filter maybe?
ImageEdgeDetection/M=4/N Kirsch, thisone //edge detection with fuzzy entropy
wave m_imageedges
display/k=1;appendimage m_imageedges
//here on, try detecting fiber center
duplicate/o M_imageedges,edges
duplicate/o edges,centers //to have them superimposed
// slim_it() tried an auxiliary function to slim the edges which quite works, not perfect
for(c1=0;c1<columns;c1+=1)
for(c2=0;c2<rows;c2+=1)
if(edges[c1][c2]<1) //means it is an edge point - edge is 3 pixels thick!?
if(thisone[c1][c2]>0) // and it is not a "void" between fibers
for(c3=1;c3<100;c3+=1) //moves horizontally or vertically, not sure. and scans looking for the next edge
c5=c2+c3
if(thisone[c1][c5]<1) //it didn't touch the other side till this comes true
c4=c2+c3/2 //middle position should be the center
centers[c1][c4]=255
c2=c2+10 //hopefully move away from the identified center
c3=101 //stop the for cycle... loved GOTO in basic
endif
endfor
endif
endif
endfor
endfor
let_edges_be_edges()
display/k=1;appendimage edges
display/k=1;appendimage centers
end
//******************************************************
function let_edges_be_edges()
wave centers,edges
variable c1,c2,c3,c4,c5 //counters
variable columns=1024
variable rows=1024 //one is wrong but doesn't matter for now
for(c1=0;c1<columns;c1+=1)
for(c2=0;c2<rows;c2+=1)
if(centers[c1][c2]>254)
if(edges[c1][c2]<1)
centers[c1][c2]=0
endif
endif
endfor
endfor
end
April 26, 2013 at 09:28 am - Permalink
Do you have an EDS system attached? You could run element maps or better spectral images on those samples and try to map out the fibres based on spectral features, rather than on grey scales. That could make it easier to separate the objects of interest from the background.
Regarding the wave name. Why not pass the wave as parameter to the function, e.g.
wave W_Image //this is your image
variable threshold
//main code
end
then call it:
Ideally get rid of the '.tif' ending during import.
April 26, 2013 at 03:15 am - Permalink
Thanks for the advice about the wave. When the actual algorithm will be done, I will run a loader function that loads and renames all loaded waves as wave0, wave1, wave2 etc so the beginning of the analysis will be to create strings with content "wave0", "wave1", "wave2" etc, duplicate the waves with those names into the wave ***thisone*** and go ahead with the hopefully working code that I pasted above.
I forgot to add, I don't need indeed the distribution of diameters-widths along each fiber, an overall distribution is more than fine. I tried to FFT the image with the edges and create an isotropic Power Spectral Density. That looks pretty much like a SAXS intensity vs q but I wasn't able to actually extract the distribution I am looking for. A spatial transform of sort might work fine indeed, only I don't know what else to try there. I am starting to think that real "measurement" across the centers is the only way. Or perhaps growing circles, point by point in the "centers" image till the circle touches the edge.
April 26, 2013 at 03:55 am - Permalink
Otherwise I think the advice you got above is good, i.e., you should pass the wave as an argument to the function instead of using a hard-wired name.
If you want to get an over-all distribution you might be able to simplify your problem by subdividing the image into sufficiently small rectangles so that each one contains a very small number (ideally one or two) fibers. The other advantage of this approach is that you could assume that there is only a single fiber width in each one of these sub-images and a single characteristic width. To obtain that width I would use thresholding followed by particle analysis and then use the fitting ellipse (or the raw-moments data) to obtain an estimate of the width.
As for the spatial-transform approach: I am not sure how you are analyzing your data. Simple observation of the FFT is not going to be helpful. To use this approach effectively you need to model your fibers and then fit the transform results to the model to obtain an estimate of the distribution. This is theoretically similar to the problem of estimating surface roughness from laser speckle images.
AG
April 26, 2013 at 02:18 pm - Permalink
If you can do this, it might be useful to use Igor's Hough Transform (in ImageTransform) to identify fiber line segments.
April 26, 2013 at 03:29 pm - Permalink
That is equivalent to growing circles from the edge, until they touches a circle from the opposite side. Which is exactly what the Euclidean Distance Transform does, quickly and efficiently, for all edge pixels at once. (It would be much quicker still if it was implemented as a native function or an XOP, but the algorithm I implemented is pretty efficient: it is a two-pass filter using only information from nearest-neighbour pixels, so the run-time per pixel is constant.) The optional parameter bg in the Make_fastEDT function I posted allows to directly pass a grayscale image and set the threshold for the background level.
As I explained earlier, the EDT function internally keeps a record of the relative location of the nearest edge pixel. I've updated the function with a quick test whether the next pixel away from the nearest edge pixel has a different edge pixel as its nearest edge pixel, or one that is nearly diametrically opposed. If you want access to the various waves the function uses internally, just remove or comment out the "KillDatafolder :" line at the end, but with the updated function you can also simply pass another optional parameter skel, which will set all the redundant values of the EDT to NaN. As I understand it, this is what you want. Based on some quick tests on simple geometrical masks and running in on the 'area1_001b.TIF' image you posted, it seems to work quite well. It still does not give an unbiased sample of the fiber diameters, but if you weight the histogram of the skeletonized EDT with the square of the EDT (appropriate if the fibers are round), the histogram looks reasonable.
Attached is a picture of the result of:
•NewImage/K=0 root:M_DistanceTransform
•ModifyImage M_DistanceTransform ctab= {*,*,BlueGreenOrange,0}
and
•SetScale/P x 1,2,"", W_Histogram
•histogram/B=2 M_DistanceTransform, W_Histogram
•W_Histogram*= x^2 // mass-fraction of each diameter, assuming fibers are round
The histogram is based on a bit over 12500 distance measurements from your 1 megapixel image, which are calculated in about 5 seconds on a not-very-new MacBook Pro.
April 27, 2013 at 07:54 am - Permalink
thanks that is great.
These days I have programmed a little and done it in two ways:
1) Real brute force with horizontal imagelinesection of each line, measuring the width, then correcting it for the local angle between the fibre and the horizontal line - takes far too long though
2) Pretty much as described in my previous post, finding the centers and then growing two lines (one horizontal and one vertical) till these two hit the end of the fiber. Then Pitagora's theorem to evaluate the half width, assuming the fiber is locally straight. This takes a reasonable time, ~10 seconds of my ThinkPad per image, detecting 40% to 80% of the diameters. So it asks for a path, loads all 16 images (4 areas at 4 different magnifications) in the folder, recognizes the magnification, runs the analysis, saves the .tif of all produced images (edges, centers, detected diameters) and the distributions. Total 3 minutes per folder.
I will try also your EDT approach now, so the detection will be done in 2 different ways. Just to validate each other and be on the safer side - I would guess this will be the closest thing to the "real" distribution.
Thanks again.
--- Sven, I tried your code... it is marvellous! It works perfectly! ---
April 30, 2013 at 04:16 am - Permalink
A.G., is there any chance of getting the Euclidean Distance Transform added to the Imagemorphology operation? I haven's asked before, because I thought I might be the only one who would use it, but for many image processing problems, using the EDT really is an elegant (and very fast) way to do it. For example, thresholding the EDT is equivalent to binary erosion with a circular structure element, and for large circular structure elements it is much faster. In fact, I wrote this originally because I needed to do morphological operations with circles of radius > 10 pixels and with periodic boundary conditions.
Also, I think this algorithm is just about the worst case scenario for being slower than it needs to be, when written as an Igor user function. The algorithm accesses the matrix values in a systematic way, and all the operations except the square root at the end are simple operations on integers, with integers guaranteed to be smaller than the image's width, height, or square of the diagonal, respectively. It is possible to rewrite it as a combination of MatrixOps and Imagemorphology operations (using structure elements like {1,0,0} to do pixel shifts), but it becomes so complicated that it's almost impossible to read the code, it doesn't preserve how systematically the matrix values are accessed, and it takes so many of the fast operations that writing it that way isn't that much faster in the end.
May 1, 2013 at 09:12 am - Permalink
Some time back I implemented distance transform as part of the ImageTransform operation. This will be available in IP7. There are options for Manhattan, Euclidian and scaled-Euclidian distances. Looking at the code today I see that it would benefit from automatic multithreading which I may implement if I have time before IP7 beta.
AG
May 1, 2013 at 09:41 am - Permalink
I guess that's one more reason to look forward to IP7 then!
May 3, 2013 at 02:42 am - Permalink
Some of us having been wishing, pining, dreaming, praying,p leading,... for it to break cover.
Our hopes are inevitably dashed and we are told to wait.
We are a patient lot.
May 3, 2013 at 02:26 pm - Permalink
Not meant to tease. It was a direct response to:
AG
May 3, 2013 at 05:04 pm - Permalink
Sorry to resurrect an old thread, but I'm trying to accomplish exactly the same task as the original poster.
This is my first day trying out Igor Pro. I've figured out how to generate the M_DistanceTransform image, although I'm not sure how to get that scale that you've got on the side of yours. Again, I've basically been following the steps outlined above verbatim, except of course replacing my own filename in the code.
So my (newb) questions are:
Correlation between histogram and scale bar: It took me a good long while to figure out how to make the histogram properly, even following the example given. But ultimately how do I correlate the values from the histogram into real dimensions (from the scalebar on the original SEM image)?
Giving credit: What would be the preferred way for me to cite this software and the EDT code in any future publications?
Correcting for specimen coating: Also, how is the thickness of the sputtered metal coating taken into account in instances like this (or perhaps it isn't taken into account at all?).
Sorry for the obnoxiously neophyte questions.
I've attached my own image files in case they provide any clues (the cropped grayscale one is the one I'm using for my Wave file, and the original one has the scale bar).
I'm very excited at the possibilities here, just need a little nudging in the right direction. Thank you!!!
July 21, 2014 at 04:49 pm - Permalink
If I'm understanding correctly, the values I'm getting are the number of pixels in the radius of the fibers (x axis) and the frequency count of measurements matching that radius (y axis).
So it seems that all I have to do is convert the number of pixels into a nm dimension, based on the scale bar in my original image.
I'm doing a few spot checks with various images but so far this seems like a great method (assuming I've understood it correctly).
July 21, 2014 at 09:55 pm - Permalink
I guess you could cite the URL for this thread. Using the EDT for this problem seems too small a contribution to write up as a more traditional publication.
Note that the algorithm used gives only an approximation of the Euclidean Distance Transform. The algorithm is described in the paper cited at the EDT code snippet, but frankly that paper is not very good: it offers a proof that the algorithm is exact, but that is not in fact the case. The error in the proof arises because a pixel's nearest background point is not necessarily the nearest background point of any of the pixel's nearest-neighbor pixels. For example, the algorithm fails for a pixel whose nearest point is at offset (-4,1) if there are also points at (-3,3) and (-5,0), erroneously returning sqrt(18) instead of sqrt(17) for that pixel.
August 6, 2014 at 03:14 am - Permalink
That would be correct if the distance to the nearest background pixel could be relied upon to equal the radius of the fiber. For a straight fiber segment without spurious background pixels, that is the case. Where fibers cross each other or the edge of the image, it is not, but maybe the errors are small enough to be acceptable.
Issues at the edge of the image you can deal with by cropping the EDT image, discarding a border equal to the maximum fiber diameter, before taking the histogram.
Where fibers cross, you are systematically under-counting the smaller of the two fibers in the histogram and inflating the diameter of the larger fiber. For two equal-diameter fibers crossing perpendicularly, the nearest background pixel to the center of the crossing is at a distance of sqrt(2) times the fiber diameter. Unless you get multiple fibers crossing, that is the worst case for inflated values of the fiber diameter. You can estimate for yourself, from looking at the images, how badly those two caveats are likely to affect the histogram.
That is correct.
August 6, 2014 at 03:46 am - Permalink