Goal:
Search for a method to extract extreme points (local min/max) from a 3D set of volumetric data
Summary:
We performed a lengthy literature search where many methods were discovered to find various features but none for finding all local min/max points on the surface of a volumetric data set. We developed an algorithm for taking an input volumetric data set and extracting the local min/max points.
Our literatures search is summarized elsewhere:
http://www.coe.uncc.edu/~dswarren/readingFall2001.htm
Discussion:
Most feature extraction algorithms involve computing derivative approximations on the original data. These methods are sufficient for finding min/max points if only the appropriate points are considered. For example, when given a function and desiring to find local min/max points many algorithms are applicable but all only consider points that actually belong to that function. When looking at volumetric data there are many points that do not belong to the “surface” of the object of interest so are not useful when comparing properties for min/max searching.
Our first task was to ensure we only considered surface points. To do this we filtered the data by setting anything with a value < 100 to 0 and anything greater to a value of 255 (this was grayscale image data). This nicely highlights the boundary but any edge detection technique could be used since in essence this is just an attempt to consider only edge points.
Second we computed gradient vectors for the set of 0/255 data points. Gradients will only be non-zero along the edge. Next compare the gradient vectors of nearest neighbors and find the maximum change (this will be an area of highest curvature). The critical notion is to ignore points with zero gradient. This ensures that only gradient vectors belonging to the surface points are considered. A simple distance between gradient vectors is sufficient. This set of maximal gradient changes is used for extracting extrema (note that this set has all edge points as non-zero and any other point is zero. The local “brightness” is proportional to the extremity of a point. In other words, extreme points in this set are local maxima on the scale 0 to 255.)
Algorithm in Matlab (% denotes comment for rest of line):
% get gradients and perform local max gradient analysis
[FX, FY, FZ] = gradient(volDataBW); % volDataBW is 0,255 data
toothMaxPts = zeros(52,50,81); % size of original volume
% get maximum gradient differences in a point's immediate neighborhood
for i=2:51 % need to use image dimensions - hardcoded here for this image
for j=2:49
for k=2:80
% search local neighborhood for maximum difference - non-zero only
maxDiff = 0;
if (FX(i,j,k)^2 + FY(i,j,k)^2 + FZ(i,j,k)^2) > 0 % don't calculate if this point is zero
for nX = -1:1
for nY = -1:1
for nZ = -1:1
if nX ~= 0 | nY ~= 0 | nZ ~= 0
r = i + nX;
s = j + nY;
t = k + nZ;
if ((FX(r,s,t)^2 + FY(r,s,t)^2 + FZ(r,s,t)^2) > 0)
% calculate difference
diff = sqrt((FX(r,s,t) - FX(i,j,k))^2 + (FY(r,s,t) - FY(i,j,k))^2 + (FZ(r,s,t) - FZ(i,j,k))^2);
if diff > maxDiff
maxDiff = diff;
end
end
end
end
end
end
end
toothMaxPts(i,j,k) = maxDiff; % only keep the maximum difference
end
end
end
The above set is too big to be useful. Simply searching using threshold values are not sufficient since many local maxima are not globally maximal. We use the following algorithm to extract the local maxima from the above data points. This algorithm merely searches a desired block size for a local maximum. A block contains a local maxima if all other points within that block are lower than the local maxima. Unfortunately, this criteria was too harsh (resulting in no points) so it was relaxed to local maxima if within the block all other points except 1 other point was less than the original. This new set of points contains only local min/max points for the original volume.
for i=2:51 % need to use image dimensions - hardcoded here for this image
for j=2:49
for k=2:80
% search local neighborhood for maximum difference - non-zero only
peakPt = toothMaxPts(i,j,k);
foundPeak = 1.0;
count = 0.0; % check for equal height points
if peakPt > 0 % don't calculate if this point is zero
for nX = -1:1
for nY = -1:1
for nZ = -1:1
if nX ~= 0 | nY ~= 0 | nZ ~= 0
r = i + nX;
s = j + nY;
t = k + nZ;
if ( toothMaxPts(r,s,t) > peakPt)
foundPeak = 0.0;
end
if (toothMaxPts(r,s,t) == peakPt)
count = count + 1;
end
end
end
end
end
if foundPeak > 0.0 & count < 3.0 % often have 2 local max points in a region
toothPeakPts(i,j,k) = 1.0;
end
end
end
end
end
For a picture of the extreme points calculated on a tooth:
http://www.coe.uncc.edu/~dswarren/toothDataPicture.tif
Conclusion:
Our new method is a good first start for finding all local min/max points on a volumetric dataset. It does depend on the type of filtering used for edge detection but after a desired edge set is extracted, this method works well. This test set had many double pixeled maxima (where peaks consisted of 2 points not just 1) which might not be the case when using other types of volumetric data.