Multilevel thresholding is an extension to Otsu thresholding method, where we search for multiple thresholding intensities. We also mentioned in another post that Otsu’s method searches for these intensities by maximizing between-class variance.

However, because we’re dealing with multiple thresholds , we’ll need to use a different formula for computing between-class variance. The reason for this is, because we need to take into account contribution from all thresholds together.

In case you’re familiar with other segmentation processes we’ve covered on this blog, you might know already that we used an example of hysteresis thresholding. This type of thresholding uses two thresholds. But we’ve yet to use Otsu’s method to find multiple, so this is what we’re going to do here.

We can use Otsu’s method to find an arbitrary number of thresholds, but the process itself starts to lose meaning as we increase this number above 2 or 3. Therefore, I’m going to demonstrate the process for finding only two thresholds.

Furthermore, we use multilevel thresholding only when we can effectively solve a problem with it. Usually, for problems, that require of us to use more than two thresholds, we use other parameters as well. In other words, we use other descriptors such as color.

How does multilevel thresholding work?

There’s a recursive structure to this process and the number of recursions depend on the number of thresholds. I’m going to describe this process for finding two thresholds, which we also call hysteresis thresholding.

Firstly, we need to select the first threshold intensity. Secondly, we iterate through every possible intenisty that is larger than first threshold. As we iterate we compute which combination of the two thresholds yields the largest between-class variance.

When we iterate through all intensities for second threshold, we move on with the first threshold and repeat this process. By iterating both threshold intensites, we’re going to compute between-class variance for all possible combinations for two thresholds.

We use the following formula for each combination. First one is the between-class variance formula for all thresholds. And other two are cumulative sum and mean for each group of pixels.

Code

The following code demonstrates the process I described above.

public static Bitmap HysteresisThreshold(this Bitmap image)
     {
         int w = image.Width;
         int h = image.Height;

         BitmapData image_data = image.LockBits(
             new Rectangle(0, 0, w, h),
             ImageLockMode.ReadOnly,
             PixelFormat.Format24bppRgb);

         int bytes = image_data.Stride * image_data.Height;
         byte[] buffer = new byte[bytes];
         byte[] result = new byte[bytes];

         Marshal.Copy(image_data.Scan0, buffer, 0, bytes);
         image.UnlockBits(image_data);

         //Get normalized histogram
         double[] histogram = new double[256];
         for (int i = 0; i < bytes; i += 3)
         {
             histogram[buffer[i]]++;
         }
         histogram = histogram.Select(x => x / (w * h)).ToArray();

         //Get global mean
         double mg = 0;
         for (int i = 0; i < 256; i++)
         {
             mg += histogram[i] * i;
         }

         int[] thresholds = new int[2];
         double bcv = 0;
         for (int i = 1; i < 254; i++)
         {
             double cs1 = 0;
             double m1 = 0;
             for (int j = 0; j < i; j++)
             {
                 cs1 += histogram[j];
                 m1 += j * histogram[j];
             }
             m1 /= cs1;

             for (int j = i + 1; j < 253; j++)
             {
                 double cs2 = 0;
                 double m2 = 0;
                 for (int k = i; k < j; k++)
                 {
                     cs2 += histogram[k];
                     m2 += k * histogram[k];
                 }
                 m2 /= cs2;

                 double new_bcv = cs1 * Math.Pow(m1 - mg, 2) + cs2 * Math.Pow(m2 - mg, 2);
                 if (new_bcv > bcv)
                 {
                     bcv = new_bcv;
                     thresholds[0] = i;
                     thresholds[1] = j;
                 }
             }
           }

         for (int i = 0; i < bytes; i++)
         {
             if (buffer[i] < thresholds[0])
             {
                 result[i] = 0;
             }

             else if (buffer[i] >= thresholds[0] && buffer[i] < thresholds[1])
             {
                 result[i] = 128;
             }

             else if (buffer[i] > thresholds[1])
             {
                 result[i] = 255;
             }
         }

         Bitmap res_img = new Bitmap(w, h);
         BitmapData res_data = res_img.LockBits(
             new Rectangle(0, 0, w, h),
             ImageLockMode.WriteOnly,
             PixelFormat.Format24bppRgb);
         Marshal.Copy(result, 0, res_data.Scan0, bytes);
         res_img.UnlockBits(res_data);

         return res_img;
     }

Conclusion

I hope this tutorial was helpful in giving you a better understanding on how multilevel thresholding works.

You can also download the demo project and try it out yourself.