opencv——最大阈值分割

opencv1.0版本

#include "stdio.h"  
#include "cv.h"  
#include "highgui.h"  
#include "Math.h"
int Otsu(IplImage* src);


int main(int argc, char*argv[])
{
    IplImage* img = cvLoadImage("D:\\chengxu\\opencv\\Desert.jpg", 0);
    IplImage* dst = cvCreateImage(cvGetSize(img), 8, 1);
    int threshold = Otsu(img);
    printf("threshold = %d\n", threshold);
    cvThreshold(img, dst, threshold, 255, CV_THRESH_BINARY);

    cvNamedWindow("img", 1);
    cvShowImage("img", dst);


    cvWaitKey(-1);

    cvReleaseImage(&img);
    cvReleaseImage(&dst);

    cvDestroyWindow("dst");
    return 0;
}



int Otsu(IplImage* src)
{
    int height = src->height;
    int width = src->width;

    //histogram    
    float histogram[256] = { 0 };
    for (int i = 0; i < height; i++)
    {
        unsigned char* p = (unsigned char*)src->imageData + src->widthStep * i;
        for (int j = 0; j < width; j++)
        {
            histogram[*p++]++;
        }
    }
    //normalize histogram    
    int size = height * width;
    for (int i = 0; i < 256; i++)
    {
        histogram[i] = histogram[i] / size;
    }

    //average pixel value    
    float avgValue = 0;
    for (int i = 0; i < 256; i++)
    {
        avgValue += i * histogram[i];  //整幅图像的平均灰度  
    }

    int threshold;
    float maxVariance = 0;
    float w = 0, u = 0;
    for (int i = 0; i < 256; i++)
    {
        w += histogram[i];  //假设当前灰度i为阈值, 0~i 灰度的像素(假设像素值在此范围的像素叫做前景像素) 所占整幅图像的比例  
        u += i * histogram[i];  // 灰度i 之前的像素(0~i)的平均灰度值： 前景像素的平均灰度值  

        float t = avgValue * w - u;
        float variance = t * t / (w * (1 - w));
        if (variance > maxVariance)
        {
            maxVariance = variance;
            threshold = i;
        }
    }

    return threshold;
}

opencv2.0版本

#include<opencv2/opencv.hpp>
#include"cv.h"
#include"highgui.h"
#include <stdlib.h>
#include <stdio.h>

using namespace cv;
using namespace std;

int otsu(const Mat &img)//最大类间方差阈值分割
{
    float histogram[256] = { 0 };
    for (int i = 0; i<img.rows; i++)
    {
        const unsigned char* p = (const unsigned char*)img.ptr(i);
        for (int j = 0; j<img.cols; j++)
        {
            histogram[p[j]]++;
        }
    }

    float avgValue = 0;
    int numPixel = img.cols*img.rows;
    for (int i = 0; i<256; i++)
    {
        histogram[i] = histogram[i] / numPixel;
        avgValue += i*histogram[i];
    }

    int threshold = 0;
    float gmax = 0;
    float wk = 0, uk = 0;
    for (int i = 0; i<256; i++) {

        wk += histogram[i];
        uk += i*histogram[i];

        float ut = avgValue*wk - uk;
        float g = ut*ut / (wk*(1 - wk));

        if (g > gmax)
        {
            gmax = g;
            threshold = i;
        }
    }
    return threshold;
}
int main( int, char** argv )
{
    Mat src,gray;
    int width = 0;
    int height = 0;
    int channels = 0;
    const char* window_name = "erzhitu";
    src = imread("C:\\Program Files (x86)\\MesaImaging\\Swissranger\\matlab\\swissranger\\xx.jpg");
    if (src.empty())
    {
        return -1;
    }
    width = src.cols;//求图像的行列及通道
    height = src.rows;
    channels = src.channels();
    cvtColor(new_src, gray, CV_BGR2GRAY);//灰度化

    int gray_threshold = 0;//阈值分割
    gray_threshold = otsu(gray);

    cout << gray_threshold << endl;
    Mat erzhitu;
    threshold(gray, erzhitu, gray_threshold, 255,1);
    namedWindow(window_name, WINDOW_AUTOSIZE);
    imshow(window_name, erzhitu);