opencv——最大阈值分割

#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;
}

#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);