Can we automatically find the threshold from CT images having 3 different materials?

I am having 551 CT images (Fractal characteristics of unsaturated sands (figshare.com)) and i have to reconstruct the into 3d. I need to find the thresohold automatically, as each slice of ct image have different thresold values. So in order to get better volume construction using volumer viewer matlab app. I am attaching the code but i dont want to put the treshold values manually (135 and 41), as for more accurate volume construction the treshold values should be calculated differently for each slice based on histogram of each slice.

Note : Please download the images from the link provided.

Link for CT images:https://figshare.com/articles/figure/Fractal_characteristics_of_unsaturated_sands/6938261.

clear;clc;

% Download the CT Images (https://figshare.com/articles/figure/Fractal_characteristics_of_unsaturated_sands/6938261)

rfiles = dir('*.tif');

for i = 1:length(rfiles)

im = imread(rfiles(i).name);

imc(:,:,i) = im;

end

% h= histogram(imc);

%%% I want threshold no. to be automated based on each CT images inside the

%%% loop. NOT manually inserting it

imbin_w = imc <= 135 & imc >= 41; % these threshold should be different for each slice based on histogram for each slice

imbin_particles = imc > 135;

imbin_gas = imc <41;

%% Use this once you export from volume viewer

volshow(imc,objectConfig)

Answer 1

[Image Analyst]

Did you try multithresh?

help multithresh

 MULTITHRESH Multi-level image thresholds using Otsu's method.
    THRESH = MULTITHRESH(A) computes a single threshold for image A using
    Otsu's method and returns it in THRESH. THRESH can be used to convert A
    into a two-level image using IMQUANTIZE.
 
    THRESH = MULTITHRESH(A, N) computes N thresholds for image A using the
    Otsu's method and returns them in THRESH. THRESH is a 1xN vector which
    can be used to convert A into an image with (N+1) discrete levels using
    IMQUANTIZE. The maximum value allowed for N is 20.
 
    [THRESH, METRIC] = MULTITHRESH(A,...) returns the effectiveness metric
    as the second output argument. METRIC is in the range [0 1] and a
    higher value indicates greater effectiveness of the thresholds in
    separating the input image into N+1 regions based on Otsu's objective
    criterion.
 
    Class Support
    -------------
    The input image A is an array of one of the following classes: uint8,
    uint16, int16, single, or double. It must be nonsparse. N is a positive
    integer between 1 and 20, and can be any numeric class. THRESH is a 1xN
    numeric vector of the same class as A. METRIC is a scalar of class
    double.
 
    Notes
    -----
    1. A can be an array of any dimension. MULTITHRESH finds the thresholds
       based on the aggregate histogram of the entire array. An RGB image
       is considered a 3D numeric array and the thresholds are computed
       for the combined data from all the three color planes.
 
    2. MULTITHRESH uses the range of the input A [min(A(:)) max(A(:))] as
       the limits for computing image histogram which is used in subsequent
       computations. Any NaNs in A are ignored in computation. Any Infs and
       -Infs in A are counted in the first and last bin of the histogram,
       respectively.
 
    3. For N > 2, MULTITHRESH uses search-based optimization of Otsu's
       criterion to find the thresholds. The search-based optimization
       guarantees only locally optimal results. Since the chance of
       converging to local optimum increases with N, it is preferable to
       use smaller values of N, typically N < 10.
 
    4. For degenerate inputs where the number of unique values in A is less
       than or equal to N, there is no viable solution using Otsu's method.
       For such inputs, MULTITHRESH returns THRESH such that it includes
       all the unique values from A and possibly some extra values that are
       chosen arbitrarily.
 
    5. The thresholds (THRESH) returned by MULTITHRESH are in the same
       range as the input image A. This is unlike GRAYTHRESH which returns
       a normalized threshold in the range [0, 1].
 
    Reference
    ---------
    N. Otsu, "A Threshold Selection Method from Gray-Level Histograms," IEEE
    Transactions on Systems, Man, and Cybernetics, Vol. 9, No. 1, pp. 62-66,
    1979.
 
    Example
    ---------
      % This example computes multiple thresholds for an image using
      % MULTITHRESH and applies those thresholds to the image using IMQUANTIZE
      % to get segment labels.
 
      I = imread('circlesBrightDark.png');
      imshow(I, [])
      title('Original Image');
 
      % Compute the thresholds
      thresh = multithresh(I,2);
 
      % Apply the thresholds to obtain segmented image
      seg_I = imquantize(I,thresh);
 
      % Show the various segments in the segmented image in color
      RGB = label2rgb(seg_I);
      figure, imshow(RGB)
      title('Segmented Image');
 
    See also GRAYTHRESH, IMBINARIZE, IMQUANTIZE, RGB2IND.

    Documentation for multithresh
       doc multithresh

[Sahil Wani]

Thanks. It worked.

[Image Analyst]

Since this Answer solves your question, then could you please click the "Accept this answer" link to award me with "reputation points" for their efforts in helping you? They'd appreciate it. Thanks in advance. 🙂 Note: you can only accept one answer (so pick the best one) but you can click the "Vote" icon for as many Answers as you want. Voting for an answer will also award reputation points.

Answer 2

[Shreeya]

Hello

According to my understanding, you want to determine the thresholds in the images without any manual intervention. If you can determine an objective function which needs to be minimized for obtaining the required thresholding parameters, maybe bayesian optimization can help you. I have linked the resouces for helping with the implementation:

https://www.mathworks.com/help/stats/bayesian-optimization-algorithm.html