This is a description of the ring-filter order statistic algorithm published in:
Ring-filter(radius,k,image) applies the kth-order ring filter of the specified radius to the input and returns a new image. (For efficiency, you may wish to pass in the array for the output image, so that the top-level code described below can re-use temporary arrays rather than repeatedly allocating new ones.)
Ring-filter(radius,k,input-image)
diagonal-radius = round(radius/sqrt(2))
integer-radius = round(radius)
output-image = allocate new image of the same size as input-image
at each location P = (Px,Py) in input-image
extract values from the following 9 image locations
- P
- the 4 locations at distance integer-radius from P
in the horizontal, vertical, left, and right directions
- the 4 locations (Px +/- diagonal-radius, Py +/- diagonal-radius)
compute the kth order statistic of these 9 numbers and save it
in location (Px,Py) of output-image
You will need to use your favorite conventions for handling locations near the edges of the input image, and for handling locations (if any) which contain no value (e.g. the corners of a rotated image).
Calling a general-purpose sorting routine to sort 9 numbers is clearly inefficient. The following code finds the kth of 9 numbers, for k in the range [2,8]. (Readers should have no difficulty supplying the much simpler special-case code for k=0 and k=9.)
int nineorder (order, a, b, c, d, e, f, g, h, i)
int order, a, b, c, d, e, f, g, h, i;
{
int up, down;
/* sort the first 8 values
uses 25 comparisons: not optimal, but ok and easy to code */
sortfour(&a, &b, &c, &d);
sortfour(&e, &f, &g, &h);
sortfour(&a, &b, &e, &f);
sortfour(&c, &d, &g, &h);
sortfour(&c, &d, &e, &f);
/* move the two values flanking the desired position into up and down */
if (order >= 5){
if (order >= 7){
if (order == 8) {down = g; up = h;}
else {down = f; up = g;}} /* order = 7 */
else{
if (order == 5) {down = d; up = e;}
else {down = e; up = f;}}} /* order = 6 */
else
if (order == 2) {down = a; up = b;}
else if (order == 3) {down = b; up = c;}
else {down = c; up = d;}
/* compare the two flanking values to i and return the correct one */
if (i <= down) return(down);
else if (i <= up) return(i);
else return(up);
}
This uses a special-purpose utility function which destructively sorts its four inputs. (That is, when it is finished, *a <= *b <= *c <= *d.)
int sortfour(a, b, c, d)
int *a, *b, *c, *d;
{
int t1,t2,t3;
if (*a > *b){t1 = *a; *a = *b;} else t1 = *b;
if (*c > *d){t2 = *d; *d = *c;} else t2 = *c;
/* now have: a <= t1 and t2 <= d */
if (*a > t2){t3 = *a; *a = t2;} else t3 = t2;
if (*d < t1){t2 = *d; *d = t1;} else t2 = t1;
/* now have a: <= {t3,t2} <= d */
if (t3 <= t2) {*b = t3; *c = t2;} else {*b = t2; *c = t3;}
return(1);
}
All of the order statistic filters are computed by applying ring order statistic filters with a descending sequence of radii. For the max and min filters, ascending sequences of radii also work. However, for the other order statistics, the two methods are not equivalent and the sequence must descend.
The factor (downstep) controlling the rate of descent has been determined partly by experimentation. Theoretically, it is easy to show that it must be no larger than 3. However, in practice, good performance requires numbers significantly lower than 3.
The max and min filters are computed as follows. We have had good success with downstep=2.5. (Downstep can be set as high as 3 in the ideal continuous case, but we set it lower to allow for digitization effects in the real implementation.)
Max-filter (radius,downstep,image)
if (radius < 1) return ring-filter(1,1,image)
else
let short = ceiling(radius/downstep)
return max-filter(short,downstep,
ring-filter(radius - short,1,image))
The median filter is computed as follows. We have had good success setting downstep to be sqrt(2).
Median-filter (radius,downstep,image)
if (radius < 1) return ring-filter(1,5,image)
else
let short = min(ceiling(radius/downstep),radius-1)
return median-filter(short,downstep,
ring-filter(radius - short,5,image))
The other order statistics are computed exactly like the median, except that the ring-filter is run with a different order in the first iteration. That is
Order-filter (radius,downstep,image,order)
if (radius < 1) return ring-filter(1,order,image)
else
let short = min(ceiling(radius/downstep),radius-1)
return median-filter(short,downstep,
ring-filter(radius - short,order,image))
The above algorithm works well on images containing normal image noise. If you are processing images containing noise-free or low-noise textures with few (e.g. two) distinct, well-separated intensities, artifacts can happen with order statistics other than max and min. This is not a property of this algorithm, but a general problem with order-statistic filters: a sensible output value would be an intermediate intensity value which does not exist in the input image.
There are (at least) three ways to avoid artifacts. First, when extracting values for diagonal locations in the ring-filter, you can interpolate the values rather than rounding to the nearest integer location. Alternatively, you can add noise to the image before running the order statistic filter. Finally, you can resample the image with randomly perturbed locations prior to running the filter.
This page is maintained by Margaret Fleck.