一、工具:VC+OpenCV
二、语言:C++
三、原理
otsu法(最大类间方差法,有时也称之为大津算法)使用的是聚类的思想,把图像的灰度数按灰度级分成2个部分,使得两个部分之间的灰度值差异最大,每个部分之间的灰度差异最小,通过方差的计算来寻找一个合适的灰度级别 来划分。 所以 可以在二值化的时候 采用otsu算法来自动选取阈值进行二值化。otsu算法被认为是图像分割中阈值选取的最佳算法,计算简单,不受图像亮度和对比度的影响。因此,使类间方差最大的分割意味着错分概率最小。
设t为设定的阈值。
wo: 分开后 前景像素点数占图像的比例
uo: 分开后 前景像素点的平均灰度
w1:分开后 被景像素点数占图像的比例
u1: 分开后 被景像素点的平均灰度
u=w0*u0 + w1*u1 :图像总平均灰度
从L个灰度级遍历t,使得t为某个值的时候,前景和背景的方差最大, 则 这个 t 值便是我们要求得的阈值。
其中,方差的计算公式如下:
g=wo * (uo - u) * (uo - u) + w1 * (u1 - u) * (u1 - u)
[ 此公式计算量较大,可以采用: g = wo * w1 * (uo - u1) * (uo - u1) ]
由于otsu算法是对图像的灰度级进行聚类,so 在执行otsu算法之前,需要计算该图像的灰度直方图。
迭代法原理:迭代选择法是首先猜测一个初始阈值,然后再通过对图像的多趟计算对阈值进行改进的过程。重复地对图像进行阈值操作,将图像分割为对象类和背景类,然后来利用每一个类中的灰阶级别对阈值进行改进。
图像阈值分割---迭代算法
1 .处理流程:
1.为全局阈值选择一个初始估计值T(图像的平均灰度)。
2.用T分割图像。产生两组像素:G1有灰度值大于T的像素组成,G2有小于等于T像素组成。
3.计算G1和G2像素的平均灰度值m1和m2;
4.计算一个新的阈值:T = (m1 + m2) / 2;
5.重复步骤2和4,直到连续迭代中的T值间的差小于一个预定义参数为止。
适合图像直方图有明显波谷
四、程序
主程序(核心部分)
阈值分割 1 /*===============================图像分割=====================================*/ 2 /*---------------------------------------------------------------------------*/ 3 /*手动设置阀值*/ 4 IplImage* binaryImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1); 5 cvThreshold(smoothImgGauss,binaryImg,71,255,CV_THRESH_BINARY); 6 cvNamedWindow("cvThreshold", CV_WINDOW_AUTOSIZE ); 7 cvShowImage( "cvThreshold", binaryImg ); 8 //cvReleaseImage(&binaryImg); 9 /*---------------------------------------------------------------------------*/ 10 /*自适应阀值 //计算像域邻域的平均灰度,来决定二值化的值*/ 11 IplImage* adThresImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1); 12 double max_value=255; 13 int adpative_method=CV_ADAPTIVE_THRESH_GAUSSIAN_C;//CV_ADAPTIVE_THRESH_MEAN_C 14 int threshold_type=CV_THRESH_BINARY; 15 int block_size=3;//阈值的象素邻域大小 16 int offset=5;//窗口尺寸 17 cvAdaptiveThreshold(smoothImgGauss,adThresImg,max_value,adpative_method,threshold_type,block_size,offset); 18 cvNamedWindow("cvAdaptiveThreshold", CV_WINDOW_AUTOSIZE ); 19 cvShowImage( "cvAdaptiveThreshold", adThresImg ); 20 cvReleaseImage(&adThresImg); 21 /*---------------------------------------------------------------------------*/ 22 /*最大熵阀值分割法*/ 23 IplImage* imgMaxEntropy = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1); 24 MaxEntropy(smoothImgGauss,imgMaxEntropy); 25 cvNamedWindow("MaxEntroyThreshold", CV_WINDOW_AUTOSIZE ); 26 cvShowImage( "MaxEntroyThreshold", imgMaxEntropy );//显示图像 27 cvReleaseImage(&imgMaxEntropy ); 28 /*---------------------------------------------------------------------------*/ 29 /*基本全局阀值法*/ 30 IplImage* imgBasicGlobalThreshold = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1); 31 cvCopyImage(srcImgGrey,imgBasicGlobalThreshold); 32 int pg[256],i,thre; 33 for (i=0;i<256;i++) pg[i]=0; 34 for (i=0;i<imgBasicGlobalThreshold->imageSize;i++) // 直方图统计 35 pg[(BYTE)imgBasicGlobalThreshold->imageData[i]]++; 36 thre = BasicGlobalThreshold(pg,0,256); // 确定阈值 37 cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<thre<<endl;//输出显示阀值 38 cvThreshold(imgBasicGlobalThreshold,imgBasicGlobalThreshold,thre,255,CV_THRESH_BINARY); // 二值化 39 cvNamedWindow("BasicGlobalThreshold", CV_WINDOW_AUTOSIZE ); 40 cvShowImage( "BasicGlobalThreshold", imgBasicGlobalThreshold);//显示图像 41 cvReleaseImage(&imgBasicGlobalThreshold); 42 /*---------------------------------------------------------------------------*/ 43 /*OTSU*/ 44 IplImage* imgOtsu = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1); 45 cvCopyImage(srcImgGrey,imgOtsu); 46 int thre2; 47 thre2 = otsu2(imgOtsu); 48 cout<<"The Threshold of this Image in Otsu is:"<<thre2<<endl;//输出显示阀值 49 cvThreshold(imgOtsu,imgOtsu,thre2,255,CV_THRESH_BINARY); // 二值化 50 cvNamedWindow("imgOtsu", CV_WINDOW_AUTOSIZE ); 51 cvShowImage( "imgOtsu", imgOtsu);//显示图像 52 cvReleaseImage(&imgOtsu); 53 /*---------------------------------------------------------------------------*/ 54 /*上下阀值法:利用正态分布求可信区间*/ 55 IplImage* imgTopDown = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 ); 56 cvCopyImage(srcImgGrey,imgTopDown); 57 CvScalar mean ,std_dev;//平均值、 标准差 58 double u_threshold,d_threshold; 59 cvAvgSdv(imgTopDown,&mean,&std_dev,NULL); 60 u_threshold = mean.val[0] +2.5* std_dev.val[0];//上阀值 61 d_threshold = mean.val[0] -2.5* std_dev.val[0];//下阀值 62 //u_threshold = mean + 2.5 * std_dev; //错误 63 //d_threshold = mean - 2.5 * std_dev; 64 cout<<"The TopThreshold of this Image in TopDown is:"<<d_threshold<<endl;//输出显示阀值 65 cout<<"The DownThreshold of this Image in TopDown is:"<<u_threshold<<endl; 66 cvThreshold(imgTopDown,imgTopDown,d_threshold,u_threshold,CV_THRESH_BINARY_INV);//上下阀值 67 cvNamedWindow("imgTopDown", CV_WINDOW_AUTOSIZE ); 68 cvShowImage( "imgTopDown", imgTopDown);//显示图像 69 cvReleaseImage(&imgTopDown); 70 /*---------------------------------------------------------------------------*/ 71 /*迭代法*/ 72 IplImage* imgIteration = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 ); 73 cvCopyImage(srcImgGrey,imgIteration); 74 int thre3,nDiffRec; 75 thre3 =DetectThreshold(imgIteration, 100, nDiffRec); 76 cout<<"The Threshold of this Image in imgIteration is:"<<thre3<<endl;//输出显示阀值 77 cvThreshold(imgIteration,imgIteration,thre3,255,CV_THRESH_BINARY_INV);//上下阀值 78 cvNamedWindow("imgIteration", CV_WINDOW_AUTOSIZE ); 79 cvShowImage( "imgIteration", imgIteration); 80 cvReleaseImage(&imgIteration); 迭代 1 /*======================================================================*/ 2 /* 迭代法*/ 3 /*======================================================================*/ 4 // nMaxIter:最大迭代次数;nDiffRec:使用给定阀值确定的亮区与暗区平均灰度差异值 5 int DetectThreshold(IplImage*img, int nMaxIter, int& iDiffRec) //阀值分割:迭代法 6 { 7 //图像信息 8 int height = img->height; 9 int width = img->width; 10 int step = img->widthStep/sizeof(uchar); 11 uchar *data = (uchar*)img->imageData; 12 13 iDiffRec =0; 14 int F[256]={ 0 }; //直方图数组 15 int iTotalGray=0;//灰度值和 16 int iTotalPixel =0;//像素数和 17 byte bt;//某点的像素值 18 19 uchar iThrehold,iNewThrehold;//阀值、新阀值 20 uchar iMaxGrayValue=0,iMinGrayValue=255;//原图像中的最大灰度值和最小灰度值 21 uchar iMeanGrayValue1,iMeanGrayValue2; 22 23 //获取(i,j)的值,存于直方图数组F 24 for(int i=0;i<width;i++) 25 { 26 for(int j=0;j<height;j++) 27 { 28 bt = data[i*step+j]; 29 if(bt<iMinGrayValue) 30 iMinGrayValue = bt; 31 if(bt>iMaxGrayValue) 32 iMaxGrayValue = bt; 33 F[bt]++; 34 } 35 } 36 37 iThrehold =0;// 38 iNewThrehold = (iMinGrayValue+iMaxGrayValue)/2;//初始阀值 39 iDiffRec = iMaxGrayValue - iMinGrayValue; 40 41 for(int a=0;(abs(iThrehold-iNewThrehold)>0.5)&&a<nMaxIter;a++)//迭代中止条件 42 { 43 iThrehold = iNewThrehold; 44 //小于当前阀值部分的平均灰度值 45 for(int i=iMinGrayValue;i<iThrehold;i++) 46 { 47 iTotalGray += F[i]*i;//F[]存储图像信息 48 iTotalPixel += F[i]; 49 } 50 iMeanGrayValue1 = (uchar)(iTotalGray/iTotalPixel); 51 //大于当前阀值部分的平均灰度值 52 iTotalPixel =0; 53 iTotalGray =0; 54 for(int j=iThrehold+1;j<iMaxGrayValue;j++) 55 { 56 iTotalGray += F[j]*j;//F[]存储图像信息 57 iTotalPixel += F[j]; 58 } 59 iMeanGrayValue2 = (uchar)(iTotalGray/iTotalPixel); 60 61 iNewThrehold = (iMeanGrayValue2+iMeanGrayValue1)/2; //新阀值 62 iDiffRec = abs(iMeanGrayValue2 - iMeanGrayValue1); 63 } 64 65 //cout<<"The Threshold of this Image in imgIteration is:"<<iThrehold<<endl; 66 return iThrehold; 67 } 68 Otsu代码一 1 /*======================================================================*/ 2 /* OTSU global thresholding routine */ 3 /* takes a 2D unsigned char array pointer, number of rows, and */ 4 /* number of cols in the array. returns the value of the threshold */ 5 /*parameter: 6 *image --- buffer for image 7 rows, cols --- size of image 8 x0, y0, dx, dy --- region of vector used for computing threshold 9 vvv --- debug option, is 0, no debug information outputed 10 */ 11 /* 12 OTSU 算法可以说是自适应计算单阈值(用来转换灰度图像为二值图像)的简单高效方法。 13 下面的代码最早由 Ryan Dibble提供,此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改,补正。 14 算法对输入的灰度图像的直方图进行分析,将直方图分成两个部分,使得两部分之间的距离最大。 15 划分点就是求得的阈值。 16 */ 17 /*======================================================================*/ 18 int otsu (unsigned char*image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv) 19 { 20 21 unsigned char*np; // 图像指针 22 int thresholdValue=1; // 阈值 23 int ihist[256]; // 图像直方图,256个点 24 25 int i, j, k; // various counters 26 int n, n1, n2, gmin, gmax; 27 double m1, m2, sum, csum, fmax, sb; 28 29 // 对直方图置零 30 memset(ihist, 0, sizeof(ihist)); 31 32 gmin=255; gmax=0; 33 // 生成直方图 34 for (i = y0 +1; i < y0 + dy -1; i++) 35 { 36 np = (unsigned char*)image[i*cols+x0+1]; 37 for (j = x0 +1; j < x0 + dx -1; j++) 38 { 39 ihist[*np]++; 40 if(*np > gmax) gmax=*np; 41 if(*np < gmin) gmin=*np; 42 np++; /* next pixel */ 43 } 44 } 45 46 // set up everything 47 sum = csum =0.0; 48 n =0; 49 50 for (k =0; k <=255; k++) 51 { 52 sum += (double) k * (double) ihist[k]; /* x*f(x) 质量矩*/ 53 n += ihist[k]; /* f(x) 质量 */ 54 } 55 56 if (!n) 57 { 58 // if n has no value, there is problems... 59 fprintf (stderr, "NOT NORMAL thresholdValue = 160\n"); 60 return (160); 61 } 62 63 // do the otsu global thresholding method 64 fmax =-1.0; 65 n1 =0; 66 for (k =0; k <255; k++) 67 { 68 n1 += ihist[k]; 69 if (!n1) 70 { 71 continue; 72 } 73 n2 = n - n1; 74 if (n2 ==0) 75 { 76 break; 77 } 78 csum += (double) k *ihist[k]; 79 m1 = csum / n1; 80 m2 = (sum - csum) / n2; 81 sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2); 82 /* bbg: note: can be optimized. */ 83 if (sb > fmax) 84 { 85 fmax = sb; 86 thresholdValue = k; 87 } 88 } 89 90 // at this point we have our thresholding value 91 92 // debug code to display thresholding values 93 if ( vvv &1 ) 94 fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n", 95 thresholdValue, gmin, gmax); 96 97 return(thresholdValue); 98 } Otsu代码二 1 /*======================================================================*/ 2 /* OTSU global thresholding routine */ 3 /*======================================================================*/ 4 int otsu2 (IplImage *image) 5 { 6 int w = image->width; 7 int h = image->height; 8 9 unsigned char*np; // 图像指针 10 unsigned char pixel; 11 int thresholdValue=1; // 阈值 12 int ihist[256]; // 图像直方图,256个点 13 14 int i, j, k; // various counters 15 int n, n1, n2, gmin, gmax; 16 double m1, m2, sum, csum, fmax, sb; 17 18 // 对直方图置零... 19 memset(ihist, 0, sizeof(ihist)); 20 21 gmin=255; gmax=0; 22 // 生成直方图 23 for (i =0; i < h; i++) 24 { 25 np = (unsigned char*)(image->imageData + image->widthStep*i); 26 for (j =0; j < w; j++) 27 { 28 pixel = np[j]; 29 ihist[ pixel]++; 30 if(pixel > gmax) gmax= pixel; 31 if(pixel < gmin) gmin= pixel; 32 } 33 } 34 35 // set up everything 36 sum = csum =0.0; 37 n =0; 38 39 for (k =0; k <=255; k++) 40 { 41 sum += k * ihist[k]; /* x*f(x) 质量矩*/ 42 n += ihist[k]; /* f(x) 质量 */ 43 } 44 45 if (!n) 46 { 47 // if n has no value, there is problems... 48 //fprintf (stderr, "NOT NORMAL thresholdValue = 160\n"); 49 thresholdValue =160; 50 goto L; 51 } 52 53 // do the otsu global thresholding method 54 fmax =-1.0; 55 n1 =0; 56 for (k =0; k <255; k++) 57 { 58 n1 += ihist[k]; 59 if (!n1) { continue; } 60 n2 = n - n1; 61 if (n2 ==0) { break; } 62 csum += k *ihist[k]; 63 m1 = csum / n1; 64 m2 = (sum - csum) / n2; 65 sb = n1 * n2 *(m1 - m2) * (m1 - m2); 66 /* bbg: note: can be optimized. */ 67 if (sb > fmax) 68 { 69 fmax = sb; 70 thresholdValue = k; 71 } 72 } 73 74 L: 75 for (i =0; i < h; i++) 76 { 77 np = (unsigned char*)(image->imageData + image->widthStep*i); 78 for (j =0; j < w; j++) 79 { 80 if(np[j] >= thresholdValue) 81 np[j] =255; 82 else np[j] =0; 83 } 84 } 85 86 //cout<<"The Threshold of this Image in Otsu is:"<<thresholdValue<<endl; 87 return(thresholdValue); 88 } 最大熵阀值 1 /*============================================================================ 2 = 代码内容:最大熵阈值分割 3 = 修改日期:2009-3-3 4 = 作者:crond123 5 = 博客:http://blog.csdn.net/crond123/ 6 = E_Mail:crond123@163.com 7 ===============================================================================*/ 8 // 计算当前位置的能量熵 9 double caculateCurrentEntropy(CvHistogram * Histogram1,int cur_threshold,entropy_state state) 10 { 11 int start,end; 12 int total =0; 13 double cur_entropy =0.0; 14 if(state == back) 15 { 16 start =0; 17 end = cur_threshold; 18 } 19 else 20 { 21 start = cur_threshold; 22 end =256; 23 } 24 for(int i=start;i<end;i++) 25 { 26 total += (int)cvQueryHistValue_1D(Histogram1,i);//查询直方块的值 P304 27 } 28 for(int j=start;j<end;j++) 29 { 30 if((int)cvQueryHistValue_1D(Histogram1,j)==0) 31 continue; 32 double percentage = cvQueryHistValue_1D(Histogram1,j)/total; 33 /*熵的定义公式*/ 34 cur_entropy +=-percentage*logf(percentage); 35 /*根据泰勒展式去掉高次项得到的熵的近似计算公式 36 cur_entropy += percentage*percentage;*/ 37 } 38 return cur_entropy; 39 // return (1-cur_entropy); 40 } 41 42 //寻找最大熵阈值并分割 43 void MaxEntropy(IplImage *src,IplImage *dst) 44 { 45 assert(src != NULL); 46 assert(src->depth ==8&& dst->depth ==8); 47 assert(src->nChannels ==1); 48 CvHistogram * hist = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);//创建一个指定尺寸的直方图 49 //参数含义:直方图包含的维数、直方图维数尺寸的数组、直方图的表示格式、方块范围数组、归一化标志 50 cvCalcHist(&src,hist);//计算直方图 51 double maxentropy =-1.0; 52 int max_index =-1; 53 // 循环测试每个分割点,寻找到最大的阈值分割点 54 for(int i=0;i<HistogramBins;i++) 55 { 56 double cur_entropy = caculateCurrentEntropy(hist,i,object)+caculateCurrentEntropy(hist,i,back); 57 if(cur_entropy>maxentropy) 58 { 59 maxentropy = cur_entropy; 60 max_index = i; 61 } 62 } 63 cout<<"The Threshold of this Image in MaxEntropy is:"<<max_index<<endl; 64 cvThreshold(src, dst, (double)max_index,255, CV_THRESH_BINARY); 65 cvReleaseHist(&hist); 66 } 基本全局阀值法 1 /*============================================================================ 2 = 代码内容:基本全局阈值法 3 ==============================================================================*/ 4 int BasicGlobalThreshold(int*pg,int start,int end) 5 { // 基本全局阈值法 6 int i,t,t1,t2,k1,k2; 7 double u,u1,u2; 8 t=0; 9 u=0; 10 for (i=start;i<end;i++) 11 { 12 t+=pg[i]; 13 u+=i*pg[i]; 14 } 15 k2=(int) (u/t); // 计算此范围灰度的平均值 16 do 17 { 18 k1=k2; 19 t1=0; 20 u1=0; 21 for (i=start;i<=k1;i++) 22 { // 计算低灰度组的累加和 23 t1+=pg[i]; 24 u1+=i*pg[i]; 25 } 26 t2=t-t1; 27 u2=u-u1; 28 if (t1) 29 u1=u1/t1; // 计算低灰度组的平均值 30 else 31 u1=0; 32 if (t2) 33 u2=u2/t2; // 计算高灰度组的平均值 34 else 35 u2=0; 36 k2=(int) ((u1+u2)/2); // 得到新的阈值估计值 37 } 38 while(k1!=k2); // 数据未稳定,继续 39 //cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<k1<<endl; 40 return(k1); // 返回阈值 41 }