The Image Processing Corner

The Mean Shift filter is an image smoothing by averaging, which preserves contours.

It has two parameters: a pixel neighborhood and an intensity value which corresponds to the minimum gradient of a contour.

The principle is to apply the following algorithm for each pixel, one after the other :

• Let C be the initial color of the pixel.

• C is replaced by the average neighborhood pixel color, but computed only on pixels which color differs at most of Imin from C.

• The previous step is performed until C does not vary any more (mathematical convergence is ensured, with floating point computations).

• Store the final C color in the result image for the pixel.

The result with a 9x9 pixel window and a contour gradient minimal intensity of 20 per channel is illustrated below. The first row gives the input noisy images, the second row gives the results and the third and fourth rows show a detail with colors which enhance contrast.

Execution time on a AMD 1800+ for the 384x288 image: 0.95 seconds for gaussian and uniform noises, 0.49 seconds for "salt and pepper" and periodic noises.

First remark: this filter preserves contours which gradient exceeds Imin.

Second remark: results on gaussian and uniform noises are excellent, compared to the strong noise added.

Third remark: This filter is not adapted to the "salt and pepper" noise. Indeed, this kind of noise is considered as a contour and is not filtered.

The source code in C of the algorithm core is given below. Undeclared variables and image structure management are without ambiguity and must be added in order to compile the program.

 // Loop on pixels 
 for(int v=0; v<sizeV; ++v) {
  for(int u=0; u<sizeU; ++u) {
   
   // Compute mean box
   int uMin = utMax(u-halfWindowSize, 0);
   int uMax = utMin(u+halfWindowSize, sizeU-1);
   int vMin = utMax(v-halfWindowSize, 0);
   int vMax = utMin(v+halfWindowSize, sizeV-1);

   // Initialize pixel filtering
   int startOffset = u + v*sizeU;
   int currentR = bufferR[startOffset];
   int currentG = bufferG[startOffset];
   int currentB = bufferB[startOffset];

   // Filter the pixel
   int iterationQty = 10;  // Convergence with integers is not ensured
   while(iterationQty--) {
    // Compute mean inside noise sphere
    int qty = 0;
    int sumR = 0;
    int sumG = 0;
    int sumB = 0;
    meanIterationQty += 1;
    for(int v1=vMin; v1<=vMax; ++v1) {
     int offset = uMin + v1*sizeU;
     for(int i=1+uMax-uMin; i; --i, ++offset) {
      if(utAbs(currentR-(int)bufferR[offset])+
         utAbs(currentG-(int)bufferG[offset])+
         utAbs(currentB-(int)bufferB[offset])<noiseIntensity) {
       sumR += bufferR[offset];
       sumG += bufferG[offset];
       sumB += bufferB[offset];
       qty += 1;
      }
     }
    }

    // Update the color
    int oldR = currentR;
    int oldG = currentG;
    int oldB = currentB;
    currentR = sumR/qty;
    currentG = sumG/qty;
    currentB = sumB/qty;

    // Check the convergence
    if(oldR==currentR && oldG==currentG && oldB==currentB) {
     break;
    }
   } // End of iterations on the pixel

   // Store the result for this pixel 
   resultBufR[startOffset] = (unsigned char)currentR;
   resultBufG[startOffset] = (unsigned char)currentG;
   resultBufB[startOffset] = (unsigned char)currentB;

  } // End of loop on columns 
 } // End of loop on rows