The anisotropic diffusion algorithm generates a sequence of images by iteratively smoothing using a finite difference implementation of the anisotropic diffusion equation [1]
Here x denotes spatial indices, and n is the index for the image sequence.
The object is to smooth image intensity variation within the boundaries of
the object in the image and within the background in preference to smoothing
across object/background edges.
Achieving this is attained by choice of the coefficient function
.
A natural choice is for
to be a function of
,
Following [1] 3DMA implements the specific choice
It is expected that within the bondaries of either the object or
background phases,
is comparatively small,
and therefore
.
Equation (1) therefore generates stong diffusive smoothing during the iteration
.
Across edges,
is comparatively large, and
,
generating relatively small diffusive smoothing.
This functional choice privileges high-contrast edges over low-contrast edges.
There are two free parameters in the diffusion algorithm,
N the total number of iterations, and K a noise-level threshold.
With choices made for the pair (N,K) (we return below to a further
discussion of how N and K are chosen), the smoothing algorithm is
applied to the input image
producing a final diffused image
.
The user can choose to diffuse each 2D slice in the image separately
§12.2.12,
or diffuse a complete 3D image §12.2.13.
In the current version of the software package, the 3D image diffusion is
extremely memory intensive.
Depending on RAM and swap space availability as well as image size, one may
have to diffuse the 3D image in several `chunks', or use slicewise segmentation.
Slice-wise segmentation, though faster and much less memory intensive has
the disadvantage of not pickeing up edges that occur between slices.
Of the two free parameters, K is the more critical.
An option in the 3DMA package, §12.2.27, allows the user to
investigate the effects of different choices of K value by smoothing
a single slice with a user chosen set of K values and plotting
each diffused image for visual comparison.
Too low a value of K produces little change in the histogram compared
to
;
too high a value of K produces noticeable degradation/distortion of
the histogram.
An approximately optimum K value is thus revealed which produces the
sharpest resolved histogram.
There are suggestions in the literature that K be chosen, or possible set
automatically from an investigation of the histogram of
.
One direction is the implementation described by Canny [2];
a histogram of
is computed and
set to that value corresponding to the H'th percentile (e.g. 90% )
of the integral of the histogram.
is then used to produce the
next iterate image . To save computing time, the value of K
can be changed only every p iterations. This procedure does not
reduce the number of free variables, the free variable pair now
becomes (N,H).
We have not pursued this possbility in the current version of 3DMA.
With well chosen N and K, a comparison of the intensity histograms of
and
, should show cleaner separation between the two peaks
in
.
Segmenting
based upon choice of a global threshold chosen from
the smoothed histogram should produce an improved segmented image.
An option in 3DMA, §12.2.28, allows the user to investigate the
result of applying a single threshold choice to a slice diffused with
various values of K or to apply several threshold choices to a single slice
diffused with a single value of K.
If requested by the user, the diffused image
is printed out
sliceswise in binary integer format.
Additionally the histograms for each slice of the diffused image is printed.
Automatic Thresholding after Diffusion
Historically, an automatic thresholding procedure was developed for use
with the anisotropic diffusion algorithm.
In the current version of the code, this automatic thresholding is still
activated.
Typically the histogram of any final diffused slice reveals sharpened peaks,
with reduced overlap.
The automatic algorithm picks
corresponding to that
histogram bin of fewest counts in the `valley' between the peaks, and
segments the image based using this threshold.
When diffusing a 3D data set, this automotic mode is reasonable,
though choice of threshold in the valley range (which must be specified
in the input file) is susceptible to minor fluctuations in the histogram
in this range - it just picks the smallest occupancy bin.
This problem is more severe when a 3D image is diffused/segmented slicewise,
the automatic threshold algorithm will generally find different threshold values
for the different slices.
In addition to printing the diffused image files, the automatic thresholding procedure prints the segmented files in the same packed bit format as the global thresholding option.
Experience reveals that for realistic images, choice of valley threshold
is important as each threshold choise will reveal some of the weaker edges
in the image and hide others.
The user may find that he prefers to delete the segmentation files produced
by the automatic threshold procedure, and run the global thresholding option
on the diffused files
which are saved by the diffusion/auto-thresholding
run, using a global threshold choice based upon visual examination of
the histograms also printed for the diffused files.