预处理中,使用了OTSU阈值分割,opencv中cvThreshold函数可以一步实现,但发现效果不理想,改用了先前自己实现的的办法。轮廓只保留面积最大的,然后找出圆心并在图上标注。
代码:
#include "cv.h"
#include "highgui.h"
#include <iostream>
using namespace std;
int otsu(IplImage* src)
{
CvScalar cs = cvSum(src);
int hist_size = 256; //直方图尺寸
int hist_height = 256;
float range[] = {0,255}; //灰度级的范围
float* ranges[]={range};
//创建一维直方图,统计图像在[0 255]像素的均匀分布
CvHistogram* gray_hist = cvCreateHist(1,&hist_size,CV_HIST_ARRAY,ranges,1);
//计算灰度图像的一维直方图
cvCalcHist(&src,gray_hist,0,0);
//归一化直方图
cvNormalizeHist(gray_hist,1.0);
int Width = src->width;
int Height = src->height;
double U_t = cs.val[0]/Width*Height; //图像的平均灰度
int threshold = 0;
double delta = 0;
for(int k=0; k<256; k++)
{
double u0 = 0;
double u1 = 0;
double w0 = 0;
double w1 = 0;
double delta_tmp = 0;
for(int i=0; i<=k; i++)
{
u0 += cvQueryHistValue_1D(gray_hist,i)*i; //灰度小于阈值k的像素的平均灰度值
w0 += cvQueryHistValue_1D(gray_hist,i); //灰度小于阈值k的像素的概率
}
for(int j=k+1; j<256; j++)
{
u1 += cvQueryHistValue_1D(gray_hist,j)*j; //灰度大于阈值k的像素的平均灰度值
w1 += cvQueryHistValue_1D(gray_hist,j); //灰度小于阈值k的像素的概率
}
delta_tmp = w0*pow(u0-U_t,2) + w1*pow(u1-U_t,2); //类间方差
if(delta_tmp > delta)
{
delta = delta_tmp;
threshold = k;
}
}
return threshold;
}
CvPoint process_image(IplImage* srcImage)
{
CvMemStorage* stor;
CvSeq* cont;
CvPoint center;
CvBox2D32f* box;
CvPoint* PointArray;
CvPoint2D32f* PointArray2D32f;
IplImage* copy1 = cvCloneImage(srcImage);
// 创建动态结构序列
stor = cvCreateMemStorage(0);
cont = cvCreateSeq(CV_SEQ_ELTYPE_POINT, sizeof(CvSeq), sizeof(CvPoint) , stor);
// 二值话图像.
cvThreshold( srcImage, copy1, otsu(srcImage), 255, CV_THRESH_BINARY );
// 寻找所有轮廓.
cvFindContours( copy1, stor, &cont, sizeof(CvContour),
CV_RETR_LIST, CV_CHAIN_APPROX_NONE, cvPoint(0,0));
// Clear images. IPL use.
cvZero(copy1);
double maxarea = 0;
CvSeq* pMax;
//寻找面积最大轮廓
for(;cont;cont=cont->h_next)
{
double tmparea=fabs(cvContourArea(cont));
if(tmparea > maxarea)
{
pMax = cont;
maxarea = tmparea;
}
}
int i; // Indicator of cycle
int count = pMax->total;
CvSize size;
// Number point must be more than or equal to 6 (for cvFitEllipse_32f).
if( count >= 6 )
{
// Alloc memory for contour point set.
PointArray = (CvPoint*)malloc( count*sizeof(CvPoint) );
PointArray2D32f= (CvPoint2D32f*)malloc( count*sizeof(CvPoint2D32f) );
// Alloc memory for ellipse data.
box = (CvBox2D32f*)malloc(sizeof(CvBox2D32f));
// Get contour point set.
cvCvtSeqToArray(pMax, PointArray, CV_WHOLE_SEQ);
// Convert CvPoint set to CvBox2D32f set.
for(i=0; i<count; i++)
{
PointArray2D32f[i].x = (float)PointArray[i].x;
PointArray2D32f[i].y = (float)PointArray[i].y;
}
//拟合当前轮廓.
cvFitEllipse(PointArray2D32f, count, box);
// 绘制当前轮廓.
cvDrawContours(srcImage,pMax,CV_RGB(255,255,255),
CV_RGB(255,255,255),0,1,8,cvPoint(0,0));
// Convert ellipse data from float to integer representation.
center.x = cvRound(box->center.x);
center.y = cvRound(box->center.y);
size.width = cvRound(box->size.width*0.5);
size.height = cvRound(box->size.height*0.5);
box->angle = -box->angle;
// Draw ellipse.
cvEllipse(srcImage, center, size, box->angle, 0, 360, CV_RGB(0,0,255), 1, CV_AA, 0);
}
// Free memory.
free(PointArray);
free(PointArray2D32f);
free(box);
cvReleaseImage(©1);
return center;
}
int main()
{
CvPoint pt1,pt2,center;
IplImage* src = cvLoadImage("E:\\激光\\228.bmp", 1);
IplImage* srcCopy = cvCreateImage(cvGetSize(src), IPL_DEPTH_8U, 1);
/* 转换成单通道 */
cvCvtColor(src, srcCopy, CV_RGB2GRAY);
IplImage* thresholdImage = cvCreateImage(cvGetSize(srcCopy),srcCopy->depth,srcCopy->nChannels);
cvThreshold(srcCopy,thresholdImage,otsu(srcCopy),255,CV_THRESH_TOZERO);
//cvThreshold(srcCopy, thresholdImage, 0, 255, CV_THRESH_BINARY | CV_THRESH_OTSU);
// Create windows.
cvNamedWindow("Threshold", 1);
cvNamedWindow("Source", 1);
cvNamedWindow("Result", 1);
// Show the image.
cvShowImage("Source", srcCopy);
cvShowImage("Threshold", thresholdImage);
//最小二乘法椭圆拟合
center = process_image(thresholdImage);
cout<<center.x<<endl;
cout<<center.y<<endl;
//画彩色点需要转化为RGB
IplImage* thresholdImageCopy = cvCreateImage(cvGetSize(thresholdImage), IPL_DEPTH_8U, 3);
cvCvtColor(thresholdImage, thresholdImageCopy, CV_GRAY2BGR);
pt1.x=center.x-1;pt1.y=center.y;
pt2.x=center.x+1;pt2.y=center.y;
cvLine(thresholdImageCopy, pt1, pt2, CV_RGB(255,0,0), 2, CV_AA, 0 );
pt1.x=center.x;pt1.y=center.y-1;
pt2.x=center.x;pt2.y=center.y+1;
cvLine(thresholdImageCopy, pt1, pt2, CV_RGB(255,0,0), 2, CV_AA, 0 );
cvShowImage("Result", thresholdImageCopy);
// Wait for a key stroke; the same function arranges events processing
cvWaitKey(0);
cvReleaseImage(&src);
cvReleaseImage(&thresholdImage);
cvReleaseImage(&thresholdImageCopy);
cvReleaseImage(&srcCopy);
cvDestroyWindow("Source");
cvDestroyWindow("Result");
cvDestroyWindow("Threshold");
return 0;
}