How to perform batch image thresholding with variable threshold for each image?

I have a 3d matrix of grayscale voxels that I am attempting to 'slice' into a stack of 2D images, and then perform image segmentation on each image, then recombine the result.

The problem I am having is that in the histogram, the intensity of the desired section (bone) appears to vary from image to image. Global threshold doesn't seem to work -- it includes the skin. Manual threshold works for images similar to the image in I use as the basis for determining the threshold, but fails on others. Adaptive threshold includes many portions outside the region I want, and not all the ones inside. I've also tried triangle thresholding, and hysteresis thresholding (I've had some degree of success with this but it's less than ideal.)

There are two peaks in the histogram, the dark area of empty space surrounding the subject, and the subject, respectively. The portion of the image where the intensity is beyond the second peak is the information I want.

Is anyone aware of an image thresholding technique that would enable me to output only these areas of the image?

Example (in the first image I'd want the portion where intensity > 0.425, the second > 0.5):

Image 59

Image 106

Answer 1

[Erik Johnson]

I figured out how to do this based on the second peak in the histogram. It feels kind of hackish and unwieldy, but it works.

Assuming MZ is your two dimensional image:

     hist = imhist(MZ);

     % Get counts and bins for histogram
     [counts, binlocs] = imhist(MZ);

     % Find peaks of histogram
     [peaks,locs] = findpeaks(hist);

     % Get 2nd peak by sorting in descending order and choosing second index
     peaks_sorted = sort(peaks,'descend');
     second_peak = peaks_sorted(2);

     %Find index of 2nd max
     second_peak_index = find(counts==second_peak);
     second_peak_index = max(second_peak_index); 
     % If multiple instances of intensity of second_peak exist, we want the 
     % highest index.  There is probably a better way to do this.

     %Determine intensity value of second peak
     second_peak_intensity = binlocs(second_peak_index);

     %Choose threshold value based on second peak
     thresh_value = second_peak_intensity + 0.01;

     %Create image segmentation based on dynamic threshold value
     BW = MZ > thresh_value;         
     seg = MZ;
     seg(~BW) = 0;

Answer 2

[Image Analyst]

The standard Otsu method in graythresh works well for these kinds of high contrast, bimodal images to find thresholds between the humps. However you don't want a threshold between the humps. Triangle thresholding would be perfect for getting the thresholds you want. Maybe you're not using a good implementation of it. I wrote my own. Please upload an image so I can use it.

[Erik Johnson]

I don't want it between the peaks, I want the data beyond the second peak. Since the baseline for that portion varies across images, that's where I'm having the difficulty.

[Erik Johnson]

I think what might work is a manual threshold algorithm that dynamically adjusts itself based on the the intensity of the second peak.

Something like: 1. Find 2nd peak of histogram 2. Add 0.5 to intensity level of 2nd peak 3. Use this as manual threshold minimum

I'm not well-versed in MATLAB's statistics functionality, though. Any thoughts on how I would find the intensity value that of the second peak of the histogram?

[Erik Johnson]

Here's an image:

I attempted to use: http://www.mathworks.com/matlabcentral/fileexchange/28047-gray-image-thresholding-using-the-triangle-method

But I either used it incorrectly, or it didn't work properly.

[Erik Johnson]

I fixed an error in my code that was adversely affecting the intensity values, and the histogram differences are vastly more regular now. The reconstructed 3d matrix looks much better now, but not quite there yet. I'd still like to know if you have any thoughts as to how I can dynamically determine the 'center' of the 2nd hump and use it as a baseline.

[Image Analyst]

Can you give the original image slice without all that huge white surrounding area?

[Erik Johnson]

The large surrounding area is a piece I am trying to segment out. It is present in the original slice.

[Image Analyst]

Well you can just crop that out since I'm sure it's a constant. I was going to show you my triangle threshold method but it looks like you've accepted an answer that works for you. Anyway good luck with the rest of your project.

[Erik Johnson]

Sorry, I didn't see a response from you for a few days so I went ahead and accepted my own answer. I would still like to see the triangle threshold method, if you can show me. The variable threshold method works to some degree, but the 'fudge factor' is still static and requires manual intervention, which is less than ideal for my use case.