首页 > 代码库 > 七种常见阈值分割代码(Otsu、最大熵、迭代法、自适应阀值、手动、迭代法、基本全局阈值法)

七种常见阈值分割代码(Otsu、最大熵、迭代法、自适应阀值、手动、迭代法、基本全局阈值法)

转自:http://blog.csdn.net/xw20084898/article/details/17564957

一、工具: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值间的差小于一个预定义参数为止。

       适合图像直方图有明显波谷

 

 

四、程序

主程序(核心部分) 

 

[cpp] view plain copy
 
    1. 阈值分割  
    2. /*===============================图像分割=====================================*/  
    3. /*---------------------------------------------------------------------------*/  
    4. /*手动设置阀值*/  
    5. 4 IplImage* binaryImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1);  
    6. 5 cvThreshold(smoothImgGauss,binaryImg,71,255,CV_THRESH_BINARY);   
    7. 6 cvNamedWindow("cvThreshold", CV_WINDOW_AUTOSIZE );  
    8. 7 cvShowImage( "cvThreshold", binaryImg );  
    9. //cvReleaseImage(&binaryImg);   
    10. 9  /*---------------------------------------------------------------------------*/  
    11. 10 /*自适应阀值 //计算像域邻域的平均灰度,来决定二值化的值*/  
    12. 11 IplImage* adThresImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1);  
    13. 12 double max_value=http://www.mamicode.com/255;
    14. 13 int adpative_method=CV_ADAPTIVE_THRESH_GAUSSIAN_C;//CV_ADAPTIVE_THRESH_MEAN_C  
    15. 14  int threshold_type=CV_THRESH_BINARY;  
    16. 15 int block_size=3;//阈值的象素邻域大小  
    17. 16  int offset=5;//窗口尺寸  
    18. 17   cvAdaptiveThreshold(smoothImgGauss,adThresImg,max_value,adpative_method,threshold_type,block_size,offset);  
    19. 18 cvNamedWindow("cvAdaptiveThreshold", CV_WINDOW_AUTOSIZE );  
    20. 19 cvShowImage( "cvAdaptiveThreshold", adThresImg );  
    21. 20 cvReleaseImage(&adThresImg);  
    22. 21 /*---------------------------------------------------------------------------*/  
    23. 22 /*最大熵阀值分割法*/   
    24. 23 IplImage* imgMaxEntropy = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);  
    25. 24 MaxEntropy(smoothImgGauss,imgMaxEntropy);  
    26. 25 cvNamedWindow("MaxEntroyThreshold", CV_WINDOW_AUTOSIZE );  
    27. 26 cvShowImage( "MaxEntroyThreshold", imgMaxEntropy );//显示图像  
    28. 27   cvReleaseImage(&imgMaxEntropy );   
    29. 28 /*---------------------------------------------------------------------------*/  
    30. 29 /*基本全局阀值法*/  
    31. 30 IplImage* imgBasicGlobalThreshold = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);  
    32. 31 cvCopyImage(srcImgGrey,imgBasicGlobalThreshold);  
    33. 32 int pg[256],i,thre;   
    34. 33 for (i=0;i<256;i++) pg[i]=0;  
    35. 34 for (i=0;i<imgBasicGlobalThreshold->imageSize;i++) // 直方图统计  
    36. 35   pg[(BYTE)imgBasicGlobalThreshold->imageData[i]]++;   
    37. 36 thre = BasicGlobalThreshold(pg,0,256); // 确定阈值  
    38. 37   cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<thre<<endl;//输出显示阀值  
    39. 38   cvThreshold(imgBasicGlobalThreshold,imgBasicGlobalThreshold,thre,255,CV_THRESH_BINARY); // 二值化   
    40. 39   cvNamedWindow("BasicGlobalThreshold", CV_WINDOW_AUTOSIZE );  
    41. 40 cvShowImage( "BasicGlobalThreshold", imgBasicGlobalThreshold);//显示图像  
    42. 41   cvReleaseImage(&imgBasicGlobalThreshold);  
    43. 42 /*---------------------------------------------------------------------------*/  
    44. 43 /*OTSU*/  
    45. 44 IplImage* imgOtsu = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);  
    46. 45 cvCopyImage(srcImgGrey,imgOtsu);  
    47. 46 int thre2;  
    48. 47 thre2 = otsu2(imgOtsu);  
    49. 48 cout<<"The Threshold of this Image in Otsu is:"<<thre2<<endl;//输出显示阀值  
    50. 49 cvThreshold(imgOtsu,imgOtsu,thre2,255,CV_THRESH_BINARY); // 二值化   
    51. 50 cvNamedWindow("imgOtsu", CV_WINDOW_AUTOSIZE );  
    52. 51 cvShowImage( "imgOtsu", imgOtsu);//显示图像   
    53. 52 cvReleaseImage(&imgOtsu);  
    54. 53 /*---------------------------------------------------------------------------*/  
    55. 54 /*上下阀值法:利用正态分布求可信区间*/  
    56. 55 IplImage* imgTopDown = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );  
    57. 56 cvCopyImage(srcImgGrey,imgTopDown);  
    58. 57 CvScalar mean ,std_dev;//平均值、 标准差  
    59. 58 double u_threshold,d_threshold;  
    60. 59 cvAvgSdv(imgTopDown,&mean,&std_dev,NULL);   
    61. 60 u_threshold = mean.val[0] +2.5* std_dev.val[0];//上阀值  
    62. 61 d_threshold = mean.val[0] -2.5* std_dev.val[0];//下阀值  
    63. 62 //u_threshold = mean + 2.5 * std_dev; //错误  
    64. 63 //d_threshold = mean - 2.5 * std_dev;  
    65. 64 cout<<"The TopThreshold of this Image in TopDown is:"<<d_threshold<<endl;//输出显示阀值  
    66. 65 cout<<"The DownThreshold of this Image in TopDown is:"<<u_threshold<<endl;  
    67. 66 cvThreshold(imgTopDown,imgTopDown,d_threshold,u_threshold,CV_THRESH_BINARY_INV);//上下阀值  
    68. 67 cvNamedWindow("imgTopDown", CV_WINDOW_AUTOSIZE );  
    69. 68 cvShowImage( "imgTopDown", imgTopDown);//显示图像   
    70. 69 cvReleaseImage(&imgTopDown);  
    71. 70 /*---------------------------------------------------------------------------*/  
    72. 71 /*迭代法*/  
    73. 72 IplImage* imgIteration = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );  
    74. 73 cvCopyImage(srcImgGrey,imgIteration);  
    75. 74 int thre3,nDiffRec;  
    76. 75 thre3 =DetectThreshold(imgIteration, 100, nDiffRec);  
    77. 76 cout<<"The Threshold of this Image in imgIteration is:"<<thre3<<endl;//输出显示阀值  
    78. 77 cvThreshold(imgIteration,imgIteration,thre3,255,CV_THRESH_BINARY_INV);//上下阀值  
    79. 78 cvNamedWindow("imgIteration", CV_WINDOW_AUTOSIZE );  
    80. 79 cvShowImage( "imgIteration", imgIteration);  
    81. 80 cvReleaseImage(&imgIteration);  
    82. 迭代  
    83. /*======================================================================*/  
    84. /* 迭代法*/  
    85. /*======================================================================*/  
    86. // nMaxIter:最大迭代次数;nDiffRec:使用给定阀值确定的亮区与暗区平均灰度差异值  
    87. int DetectThreshold(IplImage*img, int nMaxIter, int& iDiffRec) //阀值分割:迭代法  
    88. 6 {  
    89. //图像信息  
    90. int height = img->height;  
    91. int width = img->width;  
    92. 10 int step = img->widthStep/sizeof(uchar);  
    93. 11 uchar *data = (uchar*)img->imageData;  
    94. 12   
    95. 13 iDiffRec =0;  
    96. 14 int F[256]={ 0 }; //直方图数组  
    97. 15 int iTotalGray=0;//灰度值和  
    98. 16 int iTotalPixel =0;//像素数和  
    99. 17 byte bt;//某点的像素值  
    100. 18   
    101. 19 uchar iThrehold,iNewThrehold;//阀值、新阀值  
    102. 20 uchar iMaxGrayValue=http://www.mamicode.com/0,iMinGrayValue=255;//原图像中的最大灰度值和最小灰度值  
    103. 21 uchar iMeanGrayValue1,iMeanGrayValue2;  
    104. 22   
    105. 23 //获取(i,j)的值,存于直方图数组F  
    106. 24 for(int i=0;i<width;i++)  
    107. 25 {  
    108. 26 for(int j=0;j<height;j++)  
    109. 27 {  
    110. 28 bt = data[i*step+j];  
    111. 29 if(bt<iMinGrayValue)  
    112. 30 iMinGrayValue = bt;  
    113. 31 if(bt>iMaxGrayValue)  
    114. 32 iMaxGrayValue = bt;  
    115. 33 F[bt]++;  
    116. 34 }  
    117. 35 }  
    118. 36   
    119. 37 iThrehold =0;//  
    120. 38 iNewThrehold = (iMinGrayValue+iMaxGrayValue)/2;//初始阀值  
    121. 39 iDiffRec = iMaxGrayValue - iMinGrayValue;  
    122. 40   
    123. 41 for(int a=0;(abs(iThrehold-iNewThrehold)>0.5)&&a<nMaxIter;a++)//迭代中止条件  
    124. 42 {  
    125. 43 iThrehold = iNewThrehold;  
    126. 44 //小于当前阀值部分的平均灰度值  
    127. 45 for(int i=iMinGrayValue;i<iThrehold;i++)  
    128. 46 {  
    129. 47 iTotalGray += F[i]*i;//F[]存储图像信息  
    130. 48 iTotalPixel += F[i];  
    131. 49 }  
    132. 50 iMeanGrayValue1 = (uchar)(iTotalGray/iTotalPixel);  
    133. 51 //大于当前阀值部分的平均灰度值  
    134. 52 iTotalPixel =0;  
    135. 53 iTotalGray =0;  
    136. 54 for(int j=iThrehold+1;j<iMaxGrayValue;j++)  
    137. 55 {  
    138. 56 iTotalGray += F[j]*j;//F[]存储图像信息  
    139. 57 iTotalPixel += F[j];   
    140. 58 }  
    141. 59 iMeanGrayValue2 = (uchar)(iTotalGray/iTotalPixel);  
    142. 60   
    143. 61 iNewThrehold = (iMeanGrayValue2+iMeanGrayValue1)/2; //新阀值  
    144. 62 iDiffRec = abs(iMeanGrayValue2 - iMeanGrayValue1);  
    145. 63 }  
    146. 64   
    147. 65 //cout<<"The Threshold of this Image in imgIteration is:"<<iThrehold<<endl;  
    148. 66 return iThrehold;  
    149. 67 }  
    150. 68  
    151. Otsu代码一  
    152. /*======================================================================*/  
    153. /* OTSU global thresholding routine */  
    154. /* takes a 2D unsigned char array pointer, number of rows, and */  
    155. /* number of cols in the array. returns the value of the threshold */  
    156. /*parameter:  
    157. 6 *image --- buffer for image 
    158. 7 rows, cols --- size of image 
    159. 8 x0, y0, dx, dy --- region of vector used for computing threshold 
    160. 9 vvv --- debug option, is 0, no debug information outputed 
    161. 10 */  
    162. 11 /* 
    163. 12 OTSU 算法可以说是自适应计算单阈值(用来转换灰度图像为二值图像)的简单高效方法。 
    164. 13 下面的代码最早由 Ryan Dibble提供,此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改,补正。 
    165. 14 算法对输入的灰度图像的直方图进行分析,将直方图分成两个部分,使得两部分之间的距离最大。 
    166. 15 划分点就是求得的阈值。 
    167. 16 */  
    168. 17 /*======================================================================*/  
    169. 18 int otsu (unsigned char*image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv)  
    170. 19 {  
    171. 20   
    172. 21 unsigned char*np; // 图像指针  
    173. 22 int thresholdValue=http://www.mamicode.com/1; // 阈值  
    174. 23 int ihist[256]; // 图像直方图,256个点  
    175. 24   
    176. 25 int i, j, k; // various counters  
    177. 26 int n, n1, n2, gmin, gmax;  
    178. 27 double m1, m2, sum, csum, fmax, sb;  
    179. 28   
    180. 29 // 对直方图置零  
    181. 30 memset(ihist, 0, sizeof(ihist));  
    182. 31   
    183. 32 gmin=255; gmax=0;  
    184. 33 // 生成直方图  
    185. 34 for (i = y0 +1; i < y0 + dy -1; i++)   
    186. 35 {  
    187. 36 np = (unsigned char*)image[i*cols+x0+1];  
    188. 37 for (j = x0 +1; j < x0 + dx -1; j++)  
    189. 38 {  
    190. 39 ihist[*np]++;  
    191. 40 if(*np > gmax) gmax=*np;  
    192. 41 if(*np < gmin) gmin=*np;  
    193. 42 np++; /* next pixel */  
    194. 43 }  
    195. 44 }  
    196. 45   
    197. 46 // set up everything  
    198. 47 sum = csum =0.0;  
    199. 48 n =0;  
    200. 49   
    201. 50 for (k =0; k <=255; k++)   
    202. 51 {  
    203. 52 sum += (double) k * (double) ihist[k]; /* x*f(x) 质量矩*/  
    204. 53 n += ihist[k]; /* f(x) 质量 */  
    205. 54 }  
    206. 55   
    207. 56 if (!n)   
    208. 57 {  
    209. 58 // if n has no value, there is problems...  
    210. 59 fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");  
    211. 60 return (160);  
    212. 61 }  
    213. 62   
    214. 63 // do the otsu global thresholding method  
    215. 64 fmax =-1.0;  
    216. 65 n1 =0;  
    217. 66 for (k =0; k <255; k++)  
    218. 67 {  
    219. 68 n1 += ihist[k];  
    220. 69 if (!n1)   
    221. 70 {   
    222. 71 continue;   
    223. 72 }  
    224. 73 n2 = n - n1;  
    225. 74 if (n2 ==0)  
    226. 75 {   
    227. 76 break;   
    228. 77 }  
    229. 78 csum += (double) k *ihist[k];  
    230. 79 m1 = csum / n1;  
    231. 80 m2 = (sum - csum) / n2;  
    232. 81 sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);  
    233. 82 /* bbg: note: can be optimized. */  
    234. 83 if (sb > fmax)   
    235. 84 {  
    236. 85 fmax = sb;  
    237. 86 thresholdValue = k;  
    238. 87 }  
    239. 88 }  
    240. 89   
    241. 90 // at this point we have our thresholding value  
    242. 91   
    243. 92 // debug code to display thresholding values  
    244. 93 if ( vvv &1 )  
    245. 94 fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n",  
    246. 95 thresholdValue, gmin, gmax);  
    247. 96   
    248. 97 return(thresholdValue);  
    249. 98 }   
    250. Otsu代码二  
    251. /*======================================================================*/  
    252. /* OTSU global thresholding routine */  
    253. /*======================================================================*/  
    254. int otsu2 (IplImage *image)  
    255. 5 {  
    256. int w = image->width;  
    257. int h = image->height;  
    258. 8   
    259. 9 unsigned char*np; // 图像指针  
    260. 10 unsigned char pixel;  
    261. 11 int thresholdValue=http://www.mamicode.com/1; // 阈值  
    262. 12 int ihist[256]; // 图像直方图,256个点  
    263. 13   
    264. 14 int i, j, k; // various counters  
    265. 15 int n, n1, n2, gmin, gmax;  
    266. 16 double m1, m2, sum, csum, fmax, sb;  
    267. 17   
    268. 18 // 对直方图置零...  
    269. 19 memset(ihist, 0, sizeof(ihist));  
    270. 20   
    271. 21 gmin=255; gmax=0;  
    272. 22 // 生成直方图  
    273. 23 for (i =0; i < h; i++)   
    274. 24 {  
    275. 25 np = (unsigned char*)(image->imageData + image->widthStep*i);  
    276. 26 for (j =0; j < w; j++)   
    277. 27 {  
    278. 28 pixel = np[j];  
    279. 29 ihist[ pixel]++;  
    280. 30 if(pixel > gmax) gmax= pixel;  
    281. 31 if(pixel < gmin) gmin= pixel;  
    282. 32 }  
    283. 33 }  
    284. 34   
    285. 35 // set up everything  
    286. 36 sum = csum =0.0;  
    287. 37 n =0;  
    288. 38   
    289. 39 for (k =0; k <=255; k++)   
    290. 40 {  
    291. 41 sum += k * ihist[k]; /* x*f(x) 质量矩*/  
    292. 42 n += ihist[k]; /* f(x) 质量 */  
    293. 43 }  
    294. 44   
    295. 45 if (!n)   
    296. 46 {  
    297. 47 // if n has no value, there is problems...  
    298. 48 //fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");  
    299. 49 thresholdValue =160;  
    300. 50 goto L;  
    301. 51 }  
    302. 52   
    303. 53 // do the otsu global thresholding method  
    304. 54 fmax =-1.0;  
    305. 55 n1 =0;  
    306. 56 for (k =0; k <255; k++)   
    307. 57 {  
    308. 58 n1 += ihist[k];  
    309. 59 if (!n1) { continue; }  
    310. 60 n2 = n - n1;  
    311. 61 if (n2 ==0) { break; }  
    312. 62 csum += k *ihist[k];  
    313. 63 m1 = csum / n1;  
    314. 64 m2 = (sum - csum) / n2;  
    315. 65 sb = n1 * n2 *(m1 - m2) * (m1 - m2);  
    316. 66 /* bbg: note: can be optimized. */  
    317. 67 if (sb > fmax)  
    318. 68 {  
    319. 69 fmax = sb;  
    320. 70 thresholdValue = k;  
    321. 71 }  
    322. 72 }  
    323. 73   
    324. 74 L:  
    325. 75 for (i =0; i < h; i++)   
    326. 76 {  
    327. 77 np = (unsigned char*)(image->imageData + image->widthStep*i);  
    328. 78 for (j =0; j < w; j++)   
    329. 79 {  
    330. 80 if(np[j] >= thresholdValue)  
    331. 81 np[j] =255;  
    332. 82 else np[j] =0;  
    333. 83 }  
    334. 84 }  
    335. 85   
    336. 86 //cout<<"The Threshold of this Image in Otsu is:"<<thresholdValue<<endl;  
    337. 87 return(thresholdValue);  
    338. 88 }  
    339. 最大熵阀值  
    340. /*============================================================================ 
    341. 2 = 代码内容:最大熵阈值分割  
    342. 3 = 修改日期:2009-3-3  
    343. 4 = 作者:crond123  
    344. 5 = 博客:http://blog.csdn.net/crond123/ 
    345. 6 = E_Mail:crond123@163.com  
    346. 7 ===============================================================================*/  
    347. // 计算当前位置的能量熵  
    348. double caculateCurrentEntropy(CvHistogram * Histogram1,int cur_threshold,entropy_state state)  
    349. 10 {  
    350. 11 int start,end;  
    351. 12 int total =0;  
    352. 13 double cur_entropy =0.0;  
    353. 14 if(state == back)   
    354. 15 {  
    355. 16 start =0;  
    356. 17 end = cur_threshold;   
    357. 18 }  
    358. 19 else   
    359. 20 {  
    360. 21 start = cur_threshold;  
    361. 22 end =256;   
    362. 23 }   
    363. 24 for(int i=start;i<end;i++)   
    364. 25 {  
    365. 26 total += (int)cvQueryHistValue_1D(Histogram1,i);//查询直方块的值 P304  
    366. 27 }  
    367. 28 for(int j=start;j<end;j++)  
    368. 29 {  
    369. 30 if((int)cvQueryHistValue_1D(Histogram1,j)==0)  
    370. 31 continue;  
    371. 32 double percentage = cvQueryHistValue_1D(Histogram1,j)/total;  
    372. 33 /*熵的定义公式*/  
    373. 34 cur_entropy +=-percentage*logf(percentage);  
    374. 35 /*根据泰勒展式去掉高次项得到的熵的近似计算公式 
    375. 36 cur_entropy += percentage*percentage;*/   
    376. 37 }  
    377. 38 return cur_entropy;  
    378. 39 // return (1-cur_entropy);  
    379. 40 }  
    380. 41   
    381. 42 //寻找最大熵阈值并分割  
    382. 43 void MaxEntropy(IplImage *src,IplImage *dst)  
    383. 44 {  
    384. 45 assert(src != NULL);  
    385. 46 assert(src->depth ==8&& dst->depth ==8);  
    386. 47 assert(src->nChannels ==1);  
    387. 48 CvHistogram * hist = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);//创建一个指定尺寸的直方图  
    388. 49 //参数含义:直方图包含的维数、直方图维数尺寸的数组、直方图的表示格式、方块范围数组、归一化标志  
    389. 50 cvCalcHist(&src,hist);//计算直方图  
    390. 51 double maxentropy =-1.0;  
    391. 52 int max_index =-1;  
    392. 53 // 循环测试每个分割点,寻找到最大的阈值分割点  
    393. 54 for(int i=0;i<HistogramBins;i++)   
    394. 55 {  
    395. 56 double cur_entropy = caculateCurrentEntropy(hist,i,object)+caculateCurrentEntropy(hist,i,back);  
    396. 57 if(cur_entropy>maxentropy)  
    397. 58 {  
    398. 59 maxentropy = cur_entropy;  
    399. 60 max_index = i;  
    400. 61 }  
    401. 62 }  
    402. 63 cout<<"The Threshold of this Image in MaxEntropy is:"<<max_index<<endl;  
    403. 64 cvThreshold(src, dst, (double)max_index,255, CV_THRESH_BINARY);  
    404. 65 cvReleaseHist(&hist);  
    405. 66 }  
    406. 基本全局阀值法  
    407. /*============================================================================ 
    408. 2 = 代码内容:基本全局阈值法  
    409. 3 ==============================================================================*/  
    410. int BasicGlobalThreshold(int*pg,int start,int end)  
    411. 5 { // 基本全局阈值法  
    412. int i,t,t1,t2,k1,k2;  
    413. double u,u1,u2;   
    414. 8 t=0;   
    415. 9 u=0;  
    416. 10 for (i=start;i<end;i++)   
    417. 11 {  
    418. 12 t+=pg[i];   
    419. 13 u+=i*pg[i];  
    420. 14 }  
    421. 15 k2=(int) (u/t); // 计算此范围灰度的平均值   
    422. 16 do   
    423. 17 {  
    424. 18 k1=k2;  
    425. 19 t1=0;   
    426. 20 u1=0;  
    427. 21 for (i=start;i<=k1;i++)   
    428. 22 { // 计算低灰度组的累加和  
    429. 23 t1+=pg[i];   
    430. 24 u1+=i*pg[i];  
    431. 25 }  
    432. 26 t2=t-t1;  
    433. 27 u2=u-u1;  
    434. 28 if (t1)   
    435. 29 u1=u1/t1; // 计算低灰度组的平均值  
    436. 30 else   
    437. 31 u1=0;  
    438. 32 if (t2)   
    439. 33 u2=u2/t2; // 计算高灰度组的平均值  
    440. 34 else   
    441. 35 u2=0;  
    442. 36 k2=(int) ((u1+u2)/2); // 得到新的阈值估计值  
    443. 37 }  
    444. 38 while(k1!=k2); // 数据未稳定,继续  
    445. 39 //cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<k1<<endl;  
    446. 40 return(k1); // 返回阈值  
    447. 41 }  

七种常见阈值分割代码(Otsu、最大熵、迭代法、自适应阀值、手动、迭代法、基本全局阈值法)