``3D_Adaptive_thresholding'' - automatically threshold a 3-D volume

VFX name:

iAdaptiveThresh

IMPROMPTU Equivalent:

3D_Adaptive_thresholding

purpose:

To automatically threshold a 3-D volume. The brightest region is thresholded out.

The main idea of the function is the following: (1) start with a ``good'' cue slice and find a threshold using a histogram-search procedure [2]; (2) for other slices around the cue slice, use an adaptive and iteratively relaxed histogram-search procedure to calculate thresholds. See [12] for more.

input:

An 8-bit grayscale volume.

output:

An 8-bit binary-valued volume and a file called ``thXXXXXX.stats''; the string ``XXXXXX'' is some computer-generated random string like c67109. The file
``thXXXXXX.stats'' can be printed out or plotted with the ANALYZE PLOT program. Its first column gives the slice-group number and the second gives the threshold computed by the autothresholder for the slice group.

parameters:

Default parameters menu is

mode = restrict_rise
low_threshold = 120
cue_slice = 1
slices_per_group = 2
max_intergroup_change = 20
max_thresh_rise = 20
area_of_interest = no
expansion_factor = 100
levels_per_bin = 2
min_peak_sep = 30
min_peak_area = 2.00
hist_valley_dev = 50
abs_min_peak_sep = 16
abs_min_peak_area = 1
smooth_histogram = yes

The cue slice, specified by the parameter ``cue_slice'', should be a slice where the region of interest is large and has high contrast. (Typically, for a DSR LV volume, this is around slice 35.) Thresholds generally are computed for groups of slices, called slice groups. The number of slices per group is set by the parameter ``slices_per_group''. Thus, if the cue slice is 35 and we pick 3 slices per group, then one threshold is calculated for slices 35-37, the second threshold is computed for slices 38-40, etc. (Similarly, thresholds are computed for slices 32-34, 29-31, etc.)

Thresholds can go up or down from one slice group to the next. The parameters ``mode'', ``max_intergroup_change'', and ``max_thresh_rise'' control the amount that a threshold can change.

For ``mode = usual'', the parameter ``max_intergroup_change'' determines how many graylevels a new threshold can go up or down.

For ``mode = restrict rise'', the parameter ``max_intergroup_change'' determines how many graylevels a new threshold can go down, while the parameter ``max_thresh_rise'' determines how many graylevels a threshold can go up. This mode is useful for situations where the region of interest is weakly separated from the background on a few slices.

The majority of the parameters regulate the histogram search procedure. For a particular group of slices, one grayscale histogram is computed. The parameter ``levels_per_bin'' sets the x-axis resolution. (I recommend setting this to > 1 to give a fuller histogram.) The parameter ``smooth_histogram'' gives one the option of smoothing the histogram - I recommend this.

Figure depicts a grayscale histogram and its associated parameters. To find a threshold, the histogram search procedure first locates the two highest (in grayscale) peaks and then finds the middle of the deepest valid valley between these peaks. The highest peak must cover a minimum area of the histogram - this keeps spurious very high peaks (around, say, 250) from being considered as peaks and throwing off the histogram search. The parameter ``min_peak_area'' determines this. This restriction only applies to the highest peak.

The second highest peak must be a minimum # of graylevels from the highest peak. The parameter ``min_peak_sep'' determines this.

Once the peaks are found, the deepest valid valley between these peaks is found. The deepest valid valley is the longest contiguous set of histogram bins that fall within the grayscale values $\rm G_1$ and $\rm G_2$ (see Figure ). $\rm G_1$ equals the minimum histogram value between the two peaks and

$$\rm G_2 = \left( P_1 - G_1 \right) \times X / 100$$

where X is set with the parameter ``hist_valley_dev''. The threshold equals the grayscale value corresponding to the middle of this valley.

Only the threshold picked for the cue-slice group must satisfy the conditions set by

$$\rm \lq\lq min\_peak\_sep''$$ (1)

$$\rm \lq\lq min\_peak\_area''$$ (2)

Other adjacent groups can have these conditions relaxed. If a threshold calculated for a new group is too far from the threshold from the previous group (as determined by ``mode'', ``max_intergroup_change'', and ``max_thresh_rise''), then the conditions of equations (1-2) are loosened: the peak separation is lowered by 2 graylevels and the % of total area is lowered by 0.5%. Then, a new threshold is calculated. If the new threshold still doesn't meet the intergroup conditions, the conditions of equations (1-2) are again relaxed. This process continues until a good threshold is found or until the lower limits of the histogram search are reached, per the settings of the parameters

$$\begin{eqnarray*}\rm \lq\lq abs\_min\_peak\_sep\''\\ \rm \lq\lq abs\_min\_peak\_area'' \end{eqnarray*}$$

If the limits are reached and a good threshold is not found, then by default the threshold from the previous group is used.

No threshold can drop below the value set by the parameter ``low_threshold''. This is a good safeguard against anomolies.

In general, all of the voxels in a slice group are used in the histogram analysis. If you elect mbc (minimum-bounding cuboid) analysis, per the parameter ``area_of_interest'', then the number of voxels used can be restricted to a region about the region of interest.

This works as follows. For the cue-slice group, calculate a threshold using all of the voxels in the slices. Threshold the slices and find the minimum bounding rectangle that encompasses all of the large connected clumps of voxels (small clumps are deleted by the autothresholder to reduce spurious points). We only care about the x and y dimensions of this rectangle. Expand each side of this rectangle by Y% per the setting of the parameter ``expansion_factor''. The next group of slices will only use the voxels in this x-y area on each slice of the group. The process is then repeated (get threshold, threshold the slices , get minimum bounding rectangle, etc.) for other slice groups. This method helps focus the histogram search procedure on ``important'' areas - i.e., the objects of interest, nominally assumed to be connected structures over a number of slices.

comments:

See below.

1.

A medium-speed function, taking roughly 2 minutes to run on a 100x100x100 volume.

2.

I recommend picking at least 2 slices per slice group and more than one graylevel per histogram bin; this helps make fuller histograms for the autothresholder's histogram analyzer. I also recommend always using histogram smoothing.

3.

Objects of different grayscales, provided they are reasonably large, can be separated out of a volume using the autothresholder and a few of the other ROI functions. This can be done as follows:

(a)

Copy volume 0 to volume 1.

(b)

Autothreshold volume 0; this gives the brightest object.

(c)

Compute the gray_and of volume 0 and volume 1, putting the result in volume 0 (use ``two-volume ops'' in the ``other'' category); this gives the grayscale values for the brightest object.

(d)

Subtract volume 0 from volume 1 and put the result in volume 1; this gives a volume with the brightest object removed.

(e)

Copy volume 1 to volume 2, authothreshold volume 1, etc.

4.

Don't try to threshold out things that are too small (i.e., < 0.5 % of the total volume). The histogram-search may not work too well.

5.

The mbc analysis is tricky to use. I would recommend using an expansion factor of at least 50%. MBC analysis makes the autothresholder about 50% slower.