Justifying the selection of threshold value to exclude region not to be analyzed

Hi everyone,

I'm trying to calculate porosity (pore area/(pore + matrix area)). It appears that the image below has multiple depths, and so a single threshold value can not separate all the pores from the matrix. I tried excluding the regions with intensity > a certain threshold (called upper threshold), assuming that these regions are more foreground and tend to have more sample artifacts. For the threshold that separates the pores out (lower threshold), I kept trying different values until the maximum number of pores is reached. My questions:

1. How do I modify my code so that the process of finding the lower threshold can be automatic (find lower threshold that gives maximum number of objects)?

2. How to justify my selection of the upper threshold? The intensity histograms of many of my images have a "shoulder" to the right side of the distribution. If I chose to exclude this shoulder, the ratios between the lower and the upper thresholds in many of my images were between 0.78-0.82. Thus, I now constantly use 0.8 as the ratio. In this way, the process is quite automatic. I firstly find the lower threshold which gives maximum number of pores, and the upper threshold = lower threshold/0.8. This so far gives me pretty good reproducibility among different images. But I'm wondering if this is a reasonable method, and if there's any way to justify it? Below is my code. I also applied tophat, wiener2, imclose and watershed segmentation to help detecting the pores. I appreciate any comment on my codes. Thanks very much!

imdir = 'C:\Users\wkuo7\Documents\MATLAB\20130724\'; imfile1 = '401.tif'; Iori = imread([imdir, imfile1]); imshow(Iori,'border','tight', 'initialMagnification','fit');title('original image') Icrp = imcrop(Iori,[1 1 1424 845]) imshow(Icrp,'border','tight', 'initialMagnification','fit');title('cropped image') %manually remove undesired region Iply = impoly(); Imk3 = Iply.createMask(); imshow(Imk3,'border','tight', 'initialMagnification','fit');title('manually selected mask') % apply tophat filter to even the background se = strel('disk',100); Itop = imtophat(Icrp,se); imshow(Itop,'border','tight', 'initialMagnification','fit');title('tophat filtered to even the background intensity') % average to smooth the edge of pores Iavg = wiener2(Itop,[2 2]); % normalize the grayscale intensity to between 0 & 1 Inml = mat2gray(Iavg) %normalizing intensity to 0~1 imshow(Inml,'border','tight', 'initialMagnification','fit');title('averaged by 2x2') imhist(Inml);axis auto % apply lower threshold to Inml with Imk3 as mask f = @(x) im2bw(x,0.15); Imbw = roifilt2(Inml,~Imk3,f); % use imclose to erode trivial objects NHOOD=[0,1,0;1,0,1;0,1,0]; se = strel('arbitrary', NHOOD) Icls = imclose(Imbw,se); % apply watershed segmentation to separated conjoint pores Idst = -bwdist(Icls, 'chessboard'); Iwts=watershed(Idst,8); Icls(Iwts == 0) = 1; Ilbl = bwlabel(~Icls, 4); % calculate number of labeled pores [Ilbl, num] = bwlabel(~Icls, 4) % apply upper threshold to exclude foreground regions Imk1 = im2bw(Inml, 0.2188); % overlay mask Iov1=showMaskAsOverlay(0.3,Imk1, 'y',Icrp,0); I7 = Inml; I7(~Imk3) = 0; I8=I7 I8(I7>0.2188)=1 I8(I7<0.2188)=0 imshow(I8,'border','tight', 'initialMagnification','fit');title('regions above upper threshold in manual mask') % calculate porosity porosity = bwarea(Ilbl)/(bwarea(~Imk1)-(bwarea(Imk3)-bwarea(I8)))*100 % overaly excluded regions and analyzed pores, surfaces Iov1=showMaskAsOverlay(0.3,Imk1, 'y',Icrp,0); showMaskAsOverlay(1,~Imk1,'w',[],0); showMaskAsOverlay(1,Ilbl,'k',[],0); showMaskAsOverlay(0.3,Imk3,'y',[],0); title('matrix in white & pores in black; region not analyzed in yellow')

[Image Analyst]

Please upload the original tif file, not the jpg file.

[himey]

Sorry, I don't know how to upload a tif file to here. Could you please download it from the link below? Thanks! https://www.dropbox.com/sh/ceowu2n4xw0oaci/bXzP7U7u5a

[Image Analyst]

Did the paper clip icon above the text box when you typed your comment in not work?

[himey]

It said the following words when I tried to attach the tif file

"Cannot attach this file because:

    File format is unsupported. Use any of these formats: .bmp, .csv, .fig, .gif, .jpg, .jpeg, .m, .mat, .mdl, .pdf, .png, .txt, .xls, .zip"

[Image Analyst]

OK, strange. So when you converted it or zipped it up, did it still complain? Or did you just give up? I can't get to dropbox due to company policy - it's blocked.

[Jeff E]

Below is an uncompressed PNG from the Dropbox file.

I can't really follow the code you posted all that well, or what criteria you want to use to reject thick sections. But in any case, the application of your thresholds will be more robust if they are not performed on an image that has been filtered with imtophat.

[himey]

Thank you, Jeff! Do you mean that I should take the imtophat out from my entire code, if my binary image comes from thresholding?

[himey]

Here's the zip file. I hope it will work!

[Jeff E]

imtophat will definitely change the histogram properties of your image. Indeed, using it as you say to "even the background" runs counter to your stated goal of finding, and then excluding, bright areas.

That being said, I wonder if you have looked at using a much smaller kernel size with imtophat. Say, just slightly larger than your largest pore...maybe with a radius of 15 pixels, and then using a global threshold. I got what looked to be decent results just defining pore boundaries in this manner, but I'm still unclear how you want to go about counting each one individually.

Answer 1

[Image Analyst]

For skewed histograms, I like the triangle method.

See my interactive thresholding app in my File Exchange: http://www.mathworks.com/matlabcentral/fileexchange/?term=authorid%3A31862

[Image Analyst]

There is a triangle thresholding method in the File Exchange: http://www.mathworks.com/matlabcentral/fileexchange/28047-gray-image-thresholding-using-the-triangle-method

Answer 2

[Image Analyst]

I suggest you take an entirely different approach. The diameter of a cell is not well defined. Even a person can not be sure what is a window to the cell. What you really need to do is to find some metric that distinguishes the performance of different kinds of materials. Maybe one kind of material looks like it have larger pores and is stiffer, more porous, more absorbent, more abrasive, or whatever it is that that material is supposed to be and do. Maybe you can distinguish them by measuring texture or something. So find something you can measure, and go with that. I suggest you look up how measurements of open cell and closed cell foams are done.

[himey]

Thank you, Image Analyst! I also measure the texture and some other physicochemical properties of my sample. This is hydrogel. So it's probably hard to measure pore properties using those methods for those dry samples such as foams or porous silica. But I still very appreciate your suggestions and help!