## Intensity Variation Measure

This method is based on the 3D version of the spatial gray-level dependence histogram (SGLDH). For this purpose, Chetverikov's 2D approach [3] is generalized to three dimensions for arbitrary relative positions of the compared density values.

First we define the five-dimensional cooccurrence histogram with elements that count the number of pairs of voxels that appear at a certain relative position with respect to each other and have certain gray values. For example, element h(i, j, z, d) of this histogram indicates the number of pairs of voxels that were found to be distance d apart, in the z) direction, and have density values i and j, respectively. The value of the image at noninteger positions is obtained by trilinear interpolation: The gray value at a general position (x, y, z) is assumed to be a trilinear function of the position itself, i.e.,

FIGURE 3 Two versions of the gradient-based indicatrix: (a) An original CT vertebral column image of size 159 x 159 x 289. (b) Orientation histogram where the absolute values of the gradient vectors are added up. (c) Orientation histogram where simply vectors in each orientation cone are counted.

where a1,..., a8 are some parameters. The values of these parameters are determined from the known values at the eight integer positions that surround the general position (x, y, z): ([x], [y], [z]), ([x], [y], [z] +1), ([x], [y] +1, [z]), ([x] + l, [y], [z]), ([x], [y] + 1, [z] + l), ([x] + 1, [y] + 1, [z]), ([x] + 1, [y], [z] + 1), ([x] + 1, [y] + 1, [z] + 1), where [w] means integer part of w.

Following Conners and Harlow [5], we define the inertia of the image with the help of this histogram as

i=o j = o where Ng is the number of gray values in the image, and we have used a semicolon to separate d from the rest of the variables in order to indicate that in our calculations we shall keep d fixed. We use Iz; d) to characterize the texture of the image and visualize it as a 3D structure showing the magnitude of Iz; d) in the direction z).

In practice, of course, we never construct the five-dimensional array h(i, j, z; d). Instead, we fix the value of d and the values of the directions z) at which function Iz; d) will be sampled. Then, each voxel is visited in turn and its pair location is found at the given distance and given direction; the gray value at that location is calculated, its difference from the gray value of the voxel is found, and it is squared and accumulated. This is repeated for all chosen sampling directions z). These directions could, for example, be the centers of the bins used for the GD method. However, there is a fundamental difference in the nature of the two approaches: In the GD approach, the more directions we chose, the larger the number of bins we have to populate by the same, fixed number of gradient vectors. As a result, the more bins we have, the more rough and noisy the indicatrix looks because the number of vectors we have to populate them becomes less and less sufficient. Here, the number of directions we choose are simply points at which we sample a continuous function of z). The more sampling points we choose, the better the function is sampled and the smoother it looks.