[Insight-users] Binary median filter

[Luis Ibanez]

Hi Bjorn,


Your   itkBinaryMedianImageFilter   has been assimilated !     :-)


If you cvs update, you can find it now under:

   Insight/Code/BasicFilters/itkBinaryMedianImageFilter

A test for the filter can be found under:

    Insight/Testing/Code/BasicFilters/itkBinaryMedianImageFilterTest.cxx


In order to make it more general, we added two ivars for
holding the values associated with the background and
foreground, so the two values in your image don't have
to be 0 and 255, they could be for example:15 and 113.
This is set using

  SetBackgroundValue()
  SetForegroundValue()

Using these values the filter applies the counting strategy
that you proposed.

The implementation was done by modifying the itkMedianImageFilter
which currently use the neighborhood iterators. Only the inner loop was
changed in order to replace the median computation by the counting
approach that you proposed.

As the MedianImageFilter, the new BinaryMedianImageFilter will work
on any dimension and can be templated over any InputimageType and
OutputImageType. Note that a casting will be performed when assigning
the values to the output pixels.

We will be interested in any performace test you can apply to
this filter using your current data.


Thanks lot for your contribution to ITK !


   Luis



=====================================================

Bjorn Hanch Sollie wrote:

> 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 mi n_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());
> sizevolu me->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