[Insight-users] Binary median filter
Bjorn Hanch Sollie
bhs at pvv . org
Wed, 29 May 2002 00:43:16 +0200 (CEST)
On Sun, 26 May 2002, Luis Ibanez wrote:
I have been working attempting to find an acceptable solution (in
terms of performance and quality) to my problem with smoothing binary
segmented images. I've experimented a bit with the binary filters for
dilation and erosion, and the results I've got have been acceptable so
far, although median filtering still gives superior results in terms
of quality. Thanks again, Luis, for all your help and suggestions.
I decided to solve the performance problems I had with the median
filter by writing some code (see below) to optimize performance for
binary images specifically. What this code does, is basically first
to extract the smallest possible box-shaped volume that contains the
segmented region in the image. It then determines pixel value by
counting rather than sorting or averaging. Lastly, it writes the
filtered subvolume back into a volume the same size as the original
input. Experimenting with a binary volume of 256x256x40 voxels, the
"general" median filter in ITK took 14 minutes to finish, while the
code below did the job in 25 seconds, producing the same results.
Performance is heavily dependent on the nature of the volume in
question, though.
I'm not sure if there's any interest in including a binary median
filter in ITK? In any case I have included my code to share it with
anyone who might be interested and find any use for it. It shouldn't
be too hard to read.
BinaryMedian.cxx:
ImageType::Pointer BinaryMedian(ImageType *input)
{
// Get subvolume for median filtering.
printf("Extracting subvolume for median filtering...\n");
int filter_radius = 2;
int value;
int max_x = 0;
int max_y = 0;
int max_z = 0;
int min_x = input->GetRequestedRegion().GetSize()[0] - 1;
int min_y = input->GetRequestedRegion().GetSize()[1] - 1;
int min_z = input->GetRequestedRegion().GetSize()[2] - 1;
int xdim;
int ydim;
int zdim;
bool in_object = false;
typedef itk::ImageRegionIteratorWithIndex<ImageType> ImageIterator;
ImageIterator minmaxiter(input, input->GetBufferedRegion());
minmaxiter.GoToBegin();
// Find the minimum and maxumum x, y and z values of the region to filter.
while(!minmaxiter.IsAtEnd())
{
if (!in_object && minmaxiter.Get() > 0)
{
in_object = true;
value = minmaxiter.GetIndex()[0];
if (min_x > value)
min_x = value;
value = minmaxiter.GetIndex()[1];
if (min_y > value)
min_y = value;
value = minmaxiter.GetIndex()[2];
if (min_z > value)
min_z = value;
}
if (in_object && minmaxiter.Get() == 0)
{
in_object = false;
value = minmaxiter.GetIndex()[0];
if (max_x < value)
max_x = value;
value = minmaxiter.GetIndex()[1];
if (max_y < value)
max_y = value;
value = minmaxiter.GetIndex()[2];
if (max_z < value)
max_z = value;
}
++minmaxiter;
} // while
// End find maximum and minimum x, y and z values.
// Make room for the median filter iterator.
// Crop the z-direction in our case, though, or the shape will narrow at
// the top and the bottom, which we don't want.
xdim = max_x - min_x + 1 + (2 * filter_radius);
ydim = max_y - min_y + 1 + (2 * filter_radius);
zdim = max_z - min_z + 1; // + (2 * filter_radius);
// Create the smaller image for the subvolume.
ImageType::Pointer subvolume = ImageType::New();
ImageType::RegionType region;
ImageType::IndexType index;
ImageType::SizeType size = {{xdim, ydim, zdim}};
index[0] = index[1] = index[2] = 0; // start index for the region to be processed
region.SetIndex(index);
region.SetSize(size);
subvolume->SetLargestPossibleRegion(region);
subvolume->SetBufferedRegion(region);
subvolume->SetRequestedRegion(region);
subvolume->SetSpacing(input->GetSpacing());
subvolume->Allocate();
ImageIterator subvoliter(subvolume, subvolume->GetBufferedRegion());
subvoliter.GoToBegin();
while(!subvoliter.IsAtEnd())
{
subvoliter.Set(0);
++subvoliter;
}
// End create smaller image for the subvolume.
// Copy the subvolume from the original image into our smaller image.
int p, q, r;
p = filter_radius;
q = filter_radius;
r = 0;
ImageType::IndexType idxinput;
ImageType::IndexType idxoutput;
for(int z = min_z; z <= max_z; ++z)
{
idxinput[2] = z;
idxoutput[2] = r++;
for (int y = min_y; y <= max_y; ++y)
{
idxinput[1] = y;
idxoutput[1] = q++;
for (int x = min_x; x <= max_x; ++x)
{
idxinput[0] = x;
idxoutput[0] = p++;
subvolume->SetPixel(idxoutput, input->GetPixel(idxinput));
}
p = filter_radius;
}
q = filter_radius;
}
r = 0;
// End get subvolume for median filtering.
// Optimized binary median filter
// binary 256x256x40: Conventional: 14 minutes
// binary 256x256x40: Optimized: 25 seconds
printf("Optimized median filtering...\n");
int foreground = 0;
int radval = 4;
int edge = (radval * 2) + 1;
int max_n = edge * edge * edge;
int i = 0;
int n_total;
// Create an image to hold the filtered volume.
ImageType::Pointer smoothed = ImageType::New();
smoothed->SetLargestPossibleRegion(subvolume->GetLargestPossibleRegion());
smoothed->SetBufferedRegion(subvolume->GetBufferedRegion());
smoothed->SetRequestedRegion(subvolume->GetRequestedRegion());
smoothed->SetSpacing(subvolume->GetSpacing());
smoothed->Allocate();
typedef itk::ImageRegionIteratorWithIndex<ImageType> ImageIterator;
ImageIterator smoothiter(smoothed, smoothed->GetBufferedRegion());
smoothiter.GoToBegin();
while(!smoothiter.IsAtEnd())
{
smoothiter.Set(0);
++smoothiter;
}
smoothiter.GoToBegin();
ImageType::SizeType radius = {radval, radval, radval};
itk::SmartNeighborhoodIterator<ImageType> medianiter(radius, subvolume, subvolume->GetRequestedRegion());
medianiter.GoToBegin();
::size_t center = (::size_t) (medianiter.Size() / 2);
n_total = center * 2;
// The ideas here is that since we have a binary image, we don't sum or sort
// the voxel values. We just count them.
while (!medianiter.IsAtEnd())
{
foreground = 0;
i = 0;
for (i = 0; i <= center; i++)
{
if (medianiter.GetPixel(i) > 0)
foreground++;
}
if (foreground == 0)
smoothiter.Set(0);
else if (foreground - 1 == center)
smoothiter.Set(255);
else
{
for (i; i <= n_total; i++)
{
if (medianiter.GetPixel(i) > 0)
foreground++;
}
if (foreground > center + 1)
smoothiter.Set(255);
}
++medianiter;
++smoothiter;
}
// End optimized binary median filter.
// Insert subvolume into volume of original size.
printf("Inserting subvolume into volume of original size...\n");
// Create a volume of same size as the input
ImageType::Pointer sizevolume = ImageType::New();
sizevolume->SetBufferedRegion(input->GetBufferedRegion());
sizevolume->SetRequestedRegion(input->GetRequestedRegion());
sizevolume->SetSpacing(input->GetSpacing());
sizevolume->Allocate();
ImageIterator sizeiter(sizevolume, sizevolume->GetBufferedRegion());
sizeiter.GoToBegin();
while(!sizeiter.IsAtEnd())
{
sizeiter.Set(0);
++sizeiter;
}
// End create volume of same size as the input.
// Copy the subvolume into the bigger volume.
p = q = r = 0;
min_x = min_x - (filter_radius);
max_x = max_x + (filter_radius);
min_y = min_y - (filter_radius);
max_y = max_y + (filter_radius);
min_z = min_z;
max_z = max_z;
ImageType::IndexType idxbigvol;
ImageType::IndexType idxsubvol;
for (int z = min_z; z <= max_z; ++z)
{
idxbigvol[2] = z;
idxsubvol[2] = r++;
for (int y = min_y; y <= max_y; ++y)
{
idxbigvol[1] = y;
idxsubvol[1] = q++;
for (int x = min_x; x <= max_x; ++x)
{
idxbigvol[0] = x;
idxsubvol[0] = p++;
sizevolume->SetPixel(idxbigvol, smoothed->GetPixel(idxsubvol));
}
p = 0;
}
q = 0;
}
r = 0;
// End copy subvolume into the bigger volume.
return sizevolume;
}
-Bjorn
--
The History of the Universe
Chapter 1: Bang! Chapter 2: Sss... Chapter 3: Crunch!
The End