MATLAB: How to perform batch image thresholding with variable threshold for each image

image processingimage segmentationMATLAB

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

Best Answer

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;
Related Question