原文地址:http://www.cnblogs.com/skyseraph/archive/2010/12/21/1913058.html
代码
/*===============================图像分割=====================================*/
/*---------------------------------------------------------------------------*/
/*手动设置阀值*/
IplImage* binaryImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1);
cvThreshold(smoothImgGauss,binaryImg,71,255,CV_THRESH_BINARY);
cvNamedWindow("cvThreshold", CV_WINDOW_AUTOSIZE );
cvShowImage( "cvThreshold", binaryImg );
//cvReleaseImage(&binaryImg);
/*---------------------------------------------------------------------------*/
/*自适应阀值 //计算像域邻域的平均灰度,来决定二值化的值*/
IplImage* adThresImg = cvCreateImage(cvSize(w, h),IPL_DEPTH_8U, 1);
double max_value=255;
int adpative_method=CV_ADAPTIVE_THRESH_GAUSSIAN_C;//CV_ADAPTIVE_THRESH_MEAN_C
int threshold_type=CV_THRESH_BINARY;
int block_size=3;//阈值的象素邻域大小
int offset=5;//窗口尺寸
cvAdaptiveThreshold(smoothImgGauss,adThresImg,max_value,adpative_method,threshold_type,block_size,offset);
cvNamedWindow("cvAdaptiveThreshold", CV_WINDOW_AUTOSIZE );
cvShowImage( "cvAdaptiveThreshold", adThresImg );
cvReleaseImage(&adThresImg);
/*---------------------------------------------------------------------------*/
/*最大熵阀值分割法*/
IplImage* imgMaxEntropy = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
MaxEntropy(smoothImgGauss,imgMaxEntropy);
cvNamedWindow("MaxEntroyThreshold", CV_WINDOW_AUTOSIZE );
cvShowImage( "MaxEntroyThreshold", imgMaxEntropy );//显示图像
cvReleaseImage(&imgMaxEntropy );
/*---------------------------------------------------------------------------*/
/*基本全局阀值法*/
IplImage* imgBasicGlobalThreshold = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
cvCopyImage(srcImgGrey,imgBasicGlobalThreshold);
int pg[256],i,thre;
for (i=0;i<256;i++) pg[i]=0;
for (i=0;i<imgBasicGlobalThreshold->imageSize;i++) // 直方图统计
pg[(BYTE)imgBasicGlobalThreshold->imageData[i]]++;
thre = BasicGlobalThreshold(pg,0,256); // 确定阈值
cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<thre<<endl;//输出显示阀值
cvThreshold(imgBasicGlobalThreshold,imgBasicGlobalThreshold,thre,255,CV_THRESH_BINARY); // 二值化
cvNamedWindow("BasicGlobalThreshold", CV_WINDOW_AUTOSIZE );
cvShowImage( "BasicGlobalThreshold", imgBasicGlobalThreshold);//显示图像
cvReleaseImage(&imgBasicGlobalThreshold);
/*---------------------------------------------------------------------------*/
/*OTSU*/
IplImage* imgOtsu = cvCreateImage(cvGetSize(imgGrey),IPL_DEPTH_8U,1);
cvCopyImage(srcImgGrey,imgOtsu);
int thre2;
thre2 = otsu2(imgOtsu);
cout<<"The Threshold of this Image in Otsu is:"<<thre2<<endl;//输出显示阀值
cvThreshold(imgOtsu,imgOtsu,thre2,255,CV_THRESH_BINARY); // 二值化
cvNamedWindow("imgOtsu", CV_WINDOW_AUTOSIZE );
cvShowImage( "imgOtsu", imgOtsu);//显示图像
cvReleaseImage(&imgOtsu);
/*---------------------------------------------------------------------------*/
/*上下阀值法:利用正态分布求可信区间*/
IplImage* imgTopDown = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );
cvCopyImage(srcImgGrey,imgTopDown);
CvScalar mean ,std_dev;//平均值、 标准差
double u_threshold,d_threshold;
cvAvgSdv(imgTopDown,&mean,&std_dev,NULL);
u_threshold = mean.val[0] +2.5* std_dev.val[0];//上阀值
d_threshold = mean.val[0] -2.5* std_dev.val[0];//下阀值
//u_threshold = mean + 2.5 * std_dev; //错误
//d_threshold = mean - 2.5 * std_dev;
cout<<"The TopThreshold of this Image in TopDown is:"<<d_threshold<<endl;//输出显示阀值
cout<<"The DownThreshold of this Image in TopDown is:"<<u_threshold<<endl;
cvThreshold(imgTopDown,imgTopDown,d_threshold,u_threshold,CV_THRESH_BINARY_INV);//上下阀值
cvNamedWindow("imgTopDown", CV_WINDOW_AUTOSIZE );
cvShowImage( "imgTopDown", imgTopDown);//显示图像
cvReleaseImage(&imgTopDown);
/*---------------------------------------------------------------------------*/
/*迭代法*/
IplImage* imgIteration = cvCreateImage( cvGetSize(imgGrey), IPL_DEPTH_8U, 1 );
cvCopyImage(srcImgGrey,imgIteration);
int thre3,nDiffRec;
thre3 =DetectThreshold(imgIteration, 100, nDiffRec);
cout<<"The Threshold of this Image in imgIteration is:"<<thre3<<endl;//输出显示阀值
cvThreshold(imgIteration,imgIteration,thre3,255,CV_THRESH_BINARY_INV);//上下阀值
cvNamedWindow("imgIteration", CV_WINDOW_AUTOSIZE );
cvShowImage( "imgIteration", imgIteration);
cvReleaseImage(&imgIteration);
Otsu代码一
代码
/*======================================================================*/
/* OTSU global thresholding routine */
/* takes a 2D unsigned char array pointer, number of rows, and */
/* number of cols in the array. returns the value of the threshold */
/*parameter:
*image --- buffer for image
rows, cols --- size of image
x0, y0, dx, dy --- region of vector used for computing threshold
vvv --- debug option, is 0, no debug information outputed
*/
/*
OTSU 算法可以说是自适应计算单阈值(用来转换灰度图像为二值图像)的简单高效方法。
下面的代码最早由 Ryan Dibble提供,此后经过多人Joerg.Schulenburg, R.Z.Liu 等修改,补正。
算法对输入的灰度图像的直方图进行分析,将直方图分成两个部分,使得两部分之间的距离最大。
划分点就是求得的阈值。
*/
/*======================================================================*/
int otsu (unsigned char*image, int rows, int cols, int x0, int y0, int dx, int dy, int vvv)
{
unsigned char*np; // 图像指针
int thresholdValue=1; // 阈值
int ihist[256]; // 图像直方图,256个点
int i, j, k; // various counters
int n, n1, n2, gmin, gmax;
double m1, m2, sum, csum, fmax, sb;
// 对直方图置零
memset(ihist, 0, sizeof(ihist));
gmin=255; gmax=0;
// 生成直方图
for (i = y0 +1; i < y0 + dy -1; i++)
{
np = (unsigned char*)image[i*cols+x0+1];
for (j = x0 +1; j < x0 + dx -1; j++)
{
ihist[*np]++;
if(*np > gmax) gmax=*np;
if(*np < gmin) gmin=*np;
np++; /* next pixel */
}
}
// set up everything
sum = csum =0.0;
n =0;
for (k =0; k <=255; k++)
{
sum += (double) k * (double) ihist[k]; /* x*f(x) 质量矩*/
n += ihist[k]; /* f(x) 质量 */
}
if (!n)
{
// if n has no value, there is problems...
fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
return (160);
}
// do the otsu global thresholding method
fmax =-1.0;
n1 =0;
for (k =0; k <255; k++)
{
n1 += ihist[k];
if (!n1)
{
continue;
}
n2 = n - n1;
if (n2 ==0)
{
break;
}
csum += (double) k *ihist[k];
m1 = csum / n1;
m2 = (sum - csum) / n2;
sb = (double) n1 *(double) n2 *(m1 - m2) * (m1 - m2);
/* bbg: note: can be optimized. */
if (sb > fmax)
{
fmax = sb;
thresholdValue = k;
}
}
// at this point we have our thresholding value
// debug code to display thresholding values
if ( vvv &1 )
fprintf(stderr,"# OTSU: thresholdValue = %d gmin=%d gmax=%d\n",
thresholdValue, gmin, gmax);
return(thresholdValue);
}
Otsu代码二
代码
/*======================================================================*/
/* OTSU global thresholding routine */
/*======================================================================*/
int otsu2 (IplImage *image)
{
int w = image->width;
int h = image->height;
unsigned char*np; // 图像指针
unsigned char pixel;
int thresholdValue=1; // 阈值
int ihist[256]; // 图像直方图,256个点
int i, j, k; // various counters
int n, n1, n2, gmin, gmax;
double m1, m2, sum, csum, fmax, sb;
// 对直方图置零...
memset(ihist, 0, sizeof(ihist));
gmin=255; gmax=0;
// 生成直方图
for (i =0; i < h; i++)
{
np = (unsigned char*)(image->imageData + image->widthStep*i);
for (j =0; j < w; j++)
{
pixel = np[j];
ihist[ pixel]++;
if(pixel > gmax) gmax= pixel;
if(pixel < gmin) gmin= pixel;
}
}
// set up everything
sum = csum =0.0;
n =0;
for (k =0; k <=255; k++)
{
sum += k * ihist[k]; /* x*f(x) 质量矩*/
n += ihist[k]; /* f(x) 质量 */
}
if (!n)
{
// if n has no value, there is problems...
//fprintf (stderr, "NOT NORMAL thresholdValue = 160\n");
thresholdValue =160;
goto L;
}
// do the otsu global thresholding method
fmax =-1.0;
n1 =0;
for (k =0; k <255; k++)
{
n1 += ihist[k];
if (!n1) { continue; }
n2 = n - n1;
if (n2 ==0) { break; }
csum += k *ihist[k];
m1 = csum / n1;
m2 = (sum - csum) / n2;
sb = n1 * n2 *(m1 - m2) * (m1 - m2);
/* bbg: note: can be optimized. */
if (sb > fmax)
{
fmax = sb;
thresholdValue = k;
}
}
L:
for (i =0; i < h; i++)
{
np = (unsigned char*)(image->imageData + image->widthStep*i);
for (j =0; j < w; j++)
{
if(np[j] >= thresholdValue)
np[j] =255;
else np[j] =0;
}
}
//cout<<"The Threshold of this Image in Otsu is:"<<thresholdValue<<endl;
return(thresholdValue);
}
最大熵阀值
代码
/*============================================================================
= 代码内容:最大熵阈值分割
= 修改日期:2009-3-3
= 作者:crond123
= 博客:http://blog.csdn.net/crond123/
= E_Mail:crond123@163.com
===============================================================================*/
// 计算当前位置的能量熵
double caculateCurrentEntropy(CvHistogram * Histogram1,int cur_threshold,entropy_state state)
{
int start,end;
int total =0;
double cur_entropy =0.0;
if(state == back)
{
start =0;
end = cur_threshold;
}
else
{
start = cur_threshold;
end =256;
}
for(int i=start;i<end;i++)
{
total += (int)cvQueryHistValue_1D(Histogram1,i);//查询直方块的值 P304
}
for(int j=start;j<end;j++)
{
if((int)cvQueryHistValue_1D(Histogram1,j)==0)
continue;
double percentage = cvQueryHistValue_1D(Histogram1,j)/total;
/*熵的定义公式*/
cur_entropy +=-percentage*logf(percentage);
/*根据泰勒展式去掉高次项得到的熵的近似计算公式
cur_entropy += percentage*percentage;*/
}
return cur_entropy;
// return (1-cur_entropy);
}
//寻找最大熵阈值并分割
void MaxEntropy(IplImage *src,IplImage *dst)
{
assert(src != NULL);
assert(src->depth ==8&& dst->depth ==8);
assert(src->nChannels ==1);
CvHistogram * hist = cvCreateHist(1,&HistogramBins,CV_HIST_ARRAY,HistogramRange);//创建一个指定尺寸的直方图
//参数含义:直方图包含的维数、直方图维数尺寸的数组、直方图的表示格式、方块范围数组、归一化标志
cvCalcHist(&src,hist);//计算直方图
double maxentropy =-1.0;
int max_index =-1;
// 循环测试每个分割点,寻找到最大的阈值分割点
for(int i=0;i<HistogramBins;i++)
{
double cur_entropy = caculateCurrentEntropy(hist,i,object)+caculateCurrentEntropy(hist,i,back);
if(cur_entropy>maxentropy)
{
maxentropy = cur_entropy;
max_index = i;
}
}
cout<<"The Threshold of this Image in MaxEntropy is:"<<max_index<<endl;
cvThreshold(src, dst, (double)max_index,255, CV_THRESH_BINARY);
cvReleaseHist(&hist);
}
基本全局阀值法
代码
/*============================================================================
= 代码内容:基本全局阈值法
==============================================================================*/
代码
/*============================================================================
= 代码内容:基本全局阈值法
==============================================================================*/
int BasicGlobalThreshold(int*pg,int start,int end)
{ // 基本全局阈值法
int i,t,t1,t2,k1,k2;
double u,u1,u2;
t=0;
u=0;
for (i=start;i<end;i++)
{
t+=pg[i];
u+=i*pg[i];
}
k2=(int) (u/t); // 计算此范围灰度的平均值
do
{
k1=k2;
t1=0;
u1=0;
for (i=start;i<=k1;i++)
{ // 计算低灰度组的累加和
t1+=pg[i];
u1+=i*pg[i];
}
t2=t-t1;
u2=u-u1;
if (t1)
u1=u1/t1; // 计算低灰度组的平均值
else
u1=0;
if (t2)
u2=u2/t2; // 计算高灰度组的平均值
else
u2=0;
k2=(int) ((u1+u2)/2); // 得到新的阈值估计值
}
while(k1!=k2); // 数据未稳定,继续
//cout<<"The Threshold of this Image in BasicGlobalThreshold is:"<<k1<<endl;
return(k1); // 返回阈值
}