otsu threshold and openCV2--00 : implementation of otsu threshold

    First of all, otsu threshold was proposed by Nobuyuki otsu in their paper.This post will talk about the implementation of the otsu threshold of openCV2.

     Following post will introduce the basic idea of otsu threshold based on the text book Digital image processing 3rd, this could help us understand what are the codes trying to do, I would not try to explain the details of otsu threshold becaus there are many good resouces to study, you should read the paper or find a good text book like Digital image processing 3rd if you want to know the whole story.

     Objective of otsu threshold is maximize the between-class variance(9).The basic idea is simple, if the threshold is good enough, the intensity of different class should be pretty different;In other words, if we could make the intensities between classes have biggest different, the threshold should be the best one.

    Class is the component of pixels within specific range.In otsu threshold, we classify the pixels into two classes.If C1 range is within [0, k], then C2 is within [k + 1, L - 1](assume the pixels intensity is within [0,  L - 1]).
   

   I want to give some explanation of the equations which will be used by the implementation of openCV2, without some investigation to these equations, it is almost impossible to understand what are the codes doing.

(1) is the total probability of C1

(2) is the total probability of C2.

(3) and (4) is the mean intensity of C1 and C2.

(5) is the global mean of the whole image.

(7) is deduced by (6) and will be used to compute (4).

(8) is the total variance combine by two classes

(9) is the  between-class variance.

  

    Now we should be able to study the codes.The implementation of the otsu threshold of openCV2 could be found on otsu threshold.Below is the code from openCV2(I do some small modification to it ).

static double
getThreshVal_Otsu_8u( const Mat& _src )
{
    Size size = _src.size();
    if( _src.isContinuous() )
    {
        size.width *= size.height;
        size.height = 1;
    }
    int const N = 256;
    int h[N] = {0}; //histogram
    for(int i = 0; i < size.height; ++i ) //#1
    {
        uchar const *src = _src.data + _src.step*i;
        
        for(int j = 0; j < size.width; ++j )
            ++h[src[j]];
    }

    double mu = 0, scale = 1./(size.width*size.height);
    for(int i = 0; i < N; ++i )
        mu += i * static_cast<double>(h[i]);

    mu *= scale;  //#2
    double mu1 = 0, q1 = 0;
    double max_sigma = 0, max_val = 0;

    for(int i = 0; i < N; ++i )
    {       
        double const p_i = h[i] * scale; //#3   
        q1 += p_i; //#4
        double const q2 = 1. - q1; //#5

        
        if( std::min(q1,q2) < FLT_EPSILON || 
            std::max(q1,q2) > 1. - FLT_EPSILON )
            continue;  //#6

        mu1 = (mu1 * q1 + i * p_i) / q1; //#7
        double const mu2 = (mu - q1*mu1)/q2; //#8
        double const sigma = q1 * q2 *
             (mu1 - mu2)*(mu1 - mu2); //#9
        if( sigma > max_sigma )
        {
            max_sigma = sigma;
            max_val = i;
        }
    }

    return max_val;
}

 

The algorithm of openCV2 adopt is pretty straigtforward.  

Step 1 : build the histogram of the 8bit input image  

Step 2 : find out global mean of equation(4)  

Step 3 : iterate i from 0~255, classify the classes from [0, i], [i, L-1]. Find out the i which make the equation(9) become biggest.  
 

Anatomy of codes
 

  The loop of #1 build the histogram.  

  mu(#2) equal to equation(5).

  #3 is the probability of pixel i.

  #4 is equal to equation(1).

  #5 use equation(8) to find out the result of equation(2).

  #6 pass the loop if one of the class probability is too small, this make sense since m1-m2
   would be very small in this case. If you wonder why do the codes check the min and max but
   not just checking min or max, please give this site(otsu problem) a visit, I have the same
   question too.

   #7 find out the result of equation(3)

   #8 equal to equation(7)

   #9 equal to equation(9)

Next post of otsu
 

  There are an extension of otsu--multi otsu method, I would like to implement it by openCV2 next time(openCV2 haven't implemented this algorithm yet when I write this post).