• Open CV 七种常见阈值分割


    整理摘录 skyseraph 代码所得

    出处:http://www.cnblogs.com/skyseraph/

    复制The All 可直接运行

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

    运行结果:

                      

           Orignal Img                                                      Threshold Img                                    AdaptiveThreshold Img  

                      

             BasicGlobalThreshold                                      MaxEntropy                                                  Otsu

             

                   TopDown                                                 Iteration                           

                                                                                                                                                                                                                       Hai

                                                                                                                                                                                                                      0730

  • 相关阅读:
    算法学习之基础(背包 列队 栈) 习题1.3.9 补全括号
    LVS负载均衡DR模式部署
    SQL更改表架构
    BCP导入导出MsSql
    E. Holes(分块)
    hdu6230 Palindrome(manacher+树状数组)
    Suffix(hash+lcp+二分)
    k近邻法( k-nearnest neighbor)
    《机器学习》第三章——LDA
    《机器学习》第三章——对率回归
  • 原文地址:https://www.cnblogs.com/gaohai/p/5715536.html
Copyright © 2020-2023  润新知