【OpenCV】根据点集（轮廓）拟合圆

最小二乘法拟合圆公式推导

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/opencv.hpp>

#include <vector>
#include <iostream>

using namespace cv;
using namespace std;


typedef struct Box
{
	double x;
	double y;
	double r;
}Box;

Box circleLeastFit(const vector<CvPoint *> &points)
{

	Box box;
	box.x = 0.0f;
	box.y = 0.0f;
	box.r = 0.0f;

	if (points.size() < 3)
	{
		return box;
	}

	int i=0;

	double X1=0;
	double Y1=0;
	double X2=0;
	double Y2=0;
	double X3=0;
	double Y3=0;
	double X1Y1=0;
	double X1Y2=0;
	double X2Y1=0;

	int Sum = points.size();
	for (i=0;i<Sum;i++)
	{
		X1 = X1 + points[i]->x;
		Y1 = Y1 + points[i]->y;
		X2 = X2 + points[i]->x*points[i]->x;
		Y2 = Y2 + points[i]->y*points[i]->y;
		X3 = X3 + points[i]->x*points[i]->x*points[i]->x;
		Y3 = Y3 + points[i]->y*points[i]->y*points[i]->y;
		X1Y1 = X1Y1 + points[i]->x*points[i]->y;
		X1Y2 = X1Y2 + points[i]->x*points[i]->y*points[i]->y;
		X2Y1 = X2Y1 + points[i]->x*points[i]->x*points[i]->y;
	}

	double C,D,E,G,H,N;
	double a,b,c;
	N = points.size();
	C = N*X2 - X1*X1;
	D = N*X1Y1 - X1*Y1;
	E = N*X3 + N*X1Y2 - (X2+Y2)*X1;
	G = N*Y2 - Y1*Y1;
	H = N*X2Y1 + N*Y3 - (X2+Y2)*Y1;
	a = (H*D-E*G)/(C*G-D*D);
	b = (H*C-E*D)/(D*D-G*C);
	c = -(a*X1 + b*Y1 + X2 + Y2)/N;

	double A,B,R;
	A = a/(-2);
	B = b/(-2);
	R = sqrt(a*a+b*b-4*c)/2;

	box.x= A;
	box.y = B;
	box.r = R;
	
	return box;
}

//获取自定义核
void f_GetElement(int **p,IplImage *src)//为什么p一定要是二级指针？？？
{
	int length = (src->width)*(src->height);//图片大小
	*p = (int*)malloc(sizeof(int)*length);//动态分配内存，大小为(sizeof(int)*length)
	memset(*p,0,sizeof(int)*length);//初始化内存区间 p为0

	for(int y=0; y<src->height; y++)
	{
		uchar *ptr = (uchar*)(src->imageData + y*src->widthStep);

		int b = y*src->widthStep;//???
		for(int x=0; x<src->width; x++)
		{
			if(ptr[x]==255)
			{
				(*p)[b+x] = 1;
			}
			else 
			{
				(*p)[b+x] = 0;
			}
		}
	}
}


int main()
{ 
	
	IplImage* src = cvLoadImage("5.bmp",1);
	
	IplImage *gray_image = cvCreateImage(cvGetSize(src),8,1);
	IplImage *rgb_image = cvCreateImage(cvGetSize(src),8,3);
	cvZero(gray_image);
	cvZero(rgb_image);

	cvCvtColor(src,gray_image,CV_RGB2GRAY);
		
	cvCvtColor(gray_image,rgb_image,CV_GRAY2BGR);
	

	CvMemStorage *storage = cvCreateMemStorage(0);
	CvSeq *first_contours = NULL,*temp_contours = NULL;

	int n = cvFindContours(gray_image,storage,&first_contours,sizeof(CvContour), CV_RETR_LIST, CV_CHAIN_APPROX_SIMPLE);
	printf("n = %d\n",n);
	
	temp_contours = first_contours;
	
	for (;temp_contours != NULL; temp_contours = temp_contours->h_next)
	{
		vector <CvPoint*> pt_vec;
		CvPoint *pt;
		//cvDrawContours(rgb_image, temp_contours, CV_RGB(0, 0, 255), CV_RGB(0, 0, 255), 0,2,CV_FILLED, cvPoint(0, 0));	
 		CvRect rc =cvBoundingRect(temp_contours,0);    
		
		cvRectangle(rgb_image, cvPoint(rc.x, rc.y), cvPoint(rc.x + rc.width, rc.y + rc.height),CV_RGB(255, 0, 0), 1, 8, 0);
		
		for (int i = 0; i < temp_contours->total;i++)
		{
			pt = (CvPoint *)cvGetSeqElem(temp_contours,i);
			pt_vec.push_back(pt);
		}

		double x = 0.0;
		double y = 0.0;
		double r = 0.0;
		vector<Box *>  box;
		Box box1;
		box1 = circleLeastFit(pt_vec);
		printf("%lf %lf %lf \n",box1.x,box1.y,box1.r);
		cvCircle(rgb_image,cvPoint(box1.x,box1.y),box1.r,cvScalar(155,255,0),5,8,0);
		cvCircle(rgb_image,cvPoint(box1.x,box1.y),5,cvScalar(155,255,155),5,8,0);

	}

 	cvNamedWindow("32");
 	cvShowImage("32",rgb_image);
 	cvWaitKey(0);

	return 0;
}