Multi otsu(multi-thresholding) with openCV

I am trying to carry out multi-thresholding with otsu. The method I am using currently is actually via maximising the between class variance, I have managed to get the same threshold value given as that by the OpenCV library. However, that is just via running otsu method once.

Documentation on how to do multi-level thresholding or rather recursive thresholding is rather limited. Where do I do after obtaining the original otsu's value? Would appreciate some hints, I been playing around with the code, adding one external for loop, but the next value calculated is always 254 for any given image:(

My code if need be:

//compute histogram first
cv::Mat imageh; //image edited to grayscale for histogram purpose
//imageh=image; //to delete and uncomment below;
cv::cvtColor(image, imageh, CV_BGR2GRAY);

int histSize[1] = {256}; // number of bins
float hranges[2] = {0.0, 256.0}; // min andax pixel value
const float* ranges[1] = {hranges};
int channels[1] = {0}; // only 1 channel used

cv::MatND hist;
// Compute histogram
calcHist(&imageh, 1, channels, cv::Mat(), hist, 1, histSize, ranges);

IplImage* im = new IplImage(imageh);//assign the image to an IplImage pointer
IplImage* finalIm = cvCreateImage(cvSize(im->width, im->height), IPL_DEPTH_8U, 1);
double otsuThreshold= cvThreshold(im, finalIm, 0, 255, cv::THRESH_BINARY | cv::THRESH_OTSU );

cout<<"opencv otsu gives "<<otsuThreshold<<endl;

int totalNumberOfPixels= imageh.total();
cout<<"total number of Pixels is " <<totalNumberOfPixels<< endl;


float sum = 0;
for (int t=0 ; t<256 ; t++) 
{
    sum += t * hist.at<float>(t);
}
cout<<"sum is "<<sum<<endl;

float sumB = 0; //sum of background
int wB = 0; // weight of background
int wF = 0; //weight of foreground

float varMax = 0;
int threshold = 0;

//run an iteration to find the maximum value of the between class variance(as between class variance shld be maximise)
for (int t=0 ; t<256 ; t++) 
{
       wB += hist.at<float>(t);               // Weight Background
       if (wB == 0) continue;

       wF = totalNumberOfPixels - wB;                 // Weight Foreground
       if (wF == 0) break;

       sumB += (float) (t * hist.at<float>(t));

       float mB = sumB / wB;            // Mean Background
       float mF = (sum - sumB) / wF;    // Mean Foreground

       // Calculate Between Class Variance
       float varBetween = (float)wB * (float)wF * (mB - mF) * (mB - mF);

       // Check if new maximum found
       if (varBetween > varMax) {
          varMax = varBetween;
          threshold = t;
       }
}

       cout<<"threshold value is: "<<threshold;

Solution 1:

To extend Otsu's thresholding method to multi-level thresholding the between class variance equation becomes:

Please check out Deng-Yuan Huang, Ta-Wei Lin, Wu-Chih Hu, Automatic Multilevel Thresholding Based on Two-Stage Otsu's Method with Cluster Determination by Valley Estimation, Int. Journal of Innovative Computing, 2011, 7:5631-5644 for more information.

http://www.ijicic.org/ijicic-10-05033.pdf

Here is my C# implementation of Otsu Multi for 2 thresholds:

/* Otsu (1979) - multi */

Tuple < int, int > otsuMulti(object sender, EventArgs e) {
    //image histogram
    int[] histogram = new int[256];

    //total number of pixels
    int N = 0;

    //accumulate image histogram and total number of pixels
    foreach(int intensity in image.Data) {
        if (intensity != 0) {
            histogram[intensity] += 1;
            N++;
        }
    }

    double W0K, W1K, W2K, M0, M1, M2, currVarB, optimalThresh1, optimalThresh2, maxBetweenVar, M0K, M1K, M2K, MT;

    optimalThresh1 = 0;
    optimalThresh2 = 0;

    W0K = 0;
    W1K = 0;

    M0K = 0;
    M1K = 0;

    MT = 0;
    maxBetweenVar = 0;
    for (int k = 0; k <= 255; k++) {
        MT += k * (histogram[k] / (double) N);
    }


    for (int t1 = 0; t1 <= 255; t1++) {
        W0K += histogram[t1] / (double) N; //Pi
        M0K += t1 * (histogram[t1] / (double) N); //i * Pi
        M0 = M0K / W0K; //(i * Pi)/Pi

        W1K = 0;
        M1K = 0;

        for (int t2 = t1 + 1; t2 <= 255; t2++) {
            W1K += histogram[t2] / (double) N; //Pi
            M1K += t2 * (histogram[t2] / (double) N); //i * Pi
            M1 = M1K / W1K; //(i * Pi)/Pi

            W2K = 1 - (W0K + W1K);
            M2K = MT - (M0K + M1K);

            if (W2K <= 0) break;

            M2 = M2K / W2K;

            currVarB = W0K * (M0 - MT) * (M0 - MT) + W1K * (M1 - MT) * (M1 - MT) + W2K * (M2 - MT) * (M2 - MT);

            if (maxBetweenVar < currVarB) {
                maxBetweenVar = currVarB;
                optimalThresh1 = t1;
                optimalThresh2 = t2;
            }
        }
    }

    return new Tuple(optimalThresh1, optimalThresh2);
}

And this is the result I got by thresholding an image scan of soil with the above code:

(T1 = 110, T2 = 147).

Otsu's original paper: "Nobuyuki Otsu, A Threshold Selection Method from Gray-Level Histogram, IEEE Transactions on Systems, Man, and Cybernetics, 1979, 9:62-66" also briefly mentions the extension to Multithresholding.

https://engineering.purdue.edu/kak/computervision/ECE661.08/OTSU_paper.pdf

Hope this helps.

Solution 2:

Here is a simple general approach for 'n' thresholds in python (>3.0) :

# developed by- SUJOY KUMAR GOSWAMI
# source paper- https://people.ece.cornell.edu/acharya/papers/mlt_thr_img.pdf

import cv2
import numpy as np
import math

img = cv2.imread('path-to-image')
img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
a = 0
b = 255
n = 6 # number of thresholds (better choose even value)
k = 0.7 # free variable to take any positive value
T = [] # list which will contain 'n' thresholds

def sujoy(img, a, b):
    if a>b:
        s=-1
        m=-1
        return m,s

    img = np.array(img)
    t1 = (img>=a)
    t2 = (img<=b)
    X = np.multiply(t1,t2)
    Y = np.multiply(img,X)
    s = np.sum(X)
    m = np.sum(Y)/s
    return m,s

for i in range(int(n/2-1)):
    img = np.array(img)
    t1 = (img>=a)
    t2 = (img<=b)
    X = np.multiply(t1,t2)
    Y = np.multiply(img,X)
    mu = np.sum(Y)/np.sum(X)

    Z = Y - mu
    Z = np.multiply(Z,X)
    W = np.multiply(Z,Z)
    sigma = math.sqrt(np.sum(W)/np.sum(X))

    T1 = mu - k*sigma
    T2 = mu + k*sigma

    x, y = sujoy(img, a, T1)
    w, z = sujoy(img, T2, b)

    T.append(x)
    T.append(w)

    a = T1+1
    b = T2-1
    k = k*(i+1)

T1 = mu
T2 = mu+1
x, y = sujoy(img, a, T1)
w, z = sujoy(img, T2, b)    
T.append(x)
T.append(w)
T.sort()
print(T)

For full paper and more informations visit this link.

Solution 3:

I've written an example on how otsu thresholding work in python before. You can see the source code here: https://github.com/subokita/Sandbox/blob/master/otsu.py

In the example there's 2 variants, otsu2() which is the optimised version, as seen on Wikipedia page, and otsu() which is more naive implementation based on the algorithm description itself.

If you are okay in reading python codes (in this case, they are pretty simple, almost pseudo code like), you might want to look at otsu() in the example and modify it. Porting it to C++ code is not hard either.