[MATLAB]3점에 의한 원의 중심 구하기...

원문 출처 : http://mathworld.wolfram.com/Circle.html 중 일부  입니다.

 

The equation of a circle passing through the three points for , 2, 3 (the circumcircle of the triangle determined by the points) is

(25)

The center and radius of this circle can be identified by assigning coefficients of a quadratic curve

(26)

where and (since there is no cross term). Completing the square gives

(27)

The center can then be identified as

			(28)
			(29)

 

and the radius as

(30)

where

			(31)
			(32)
			(33)
			(34)

Four or more points which lie on a circle are said to be concyclic. Three points are trivially concyclic since three noncollinear points determine a circle.

 

 

아래는 원문중 원의 공식을 이용해 MATLAB과 OpenCV 행렬 함수를 이용해서 구현된 소스입니다.

 

학교 과제할때 참고한 공식으로 유용하게 사용했습니다.^^

 

관심점 3점을 추출하여 중심점을 계산하기 원할때 사용하시면 좋을것 같습니다. 응용은 알아서 하세요^^

 

예1) 원이 찌그러 졌을때 중심위치 구하기

     바운더리 픽셀의 픽셀중심을 이용하는 것보다 곡율이 일정한 부분을 찾아서 관심점 3점을 이용하면 정확한 원의 중심을 계산할 수 있습니다. Blob의 중심 모멘트가 아닌 기하학 성분의 원의 중심을 구하는 거죠^^

 

예2) 정삼각형 무게중심, 삼각형의 외심(외적원)... (정삼각형이 아닌 일반 삼각형의 무게중심은 될수 없습니다.)

 

예3) 기타..

 

[code cpp]
% //MATLAB
% //1. x1, x2, x3, y1, y2, y3 3점의 좌표를 구한다.(관심점 3점은 원의 궤적이 될수 있는 좌표를 입력해야 합니다.)
x1 = 5; y1 = 0;   % 예) 12시 방향 좌표
x2 = 10; y2 = 5 ; % 예) 3시 방향 좌표
x3 = 5; y3  = 10; % 예) 6시 방향 좌표

% //주어짐 관심점 3점을 행렬에 대입..
a = [x1, y1, 1;    x2, y2, 1;    x3, y3, 1];
d = [x1^2+y1^2, y1, 1;    x2^2+y2^2, y2, 1;    x3^2+y3^2, y3, 1];
e = [x1^2+y1^2, x1, 1;    x2^2+y2^2, x2, 1;    x3^2+y3^2, x3, 1];
f = [x1^2+y1^2, x1, y1;    x2^2+y2^2, x2, y2;    x3^2+y3^2, x3, y3];

% //2. 3점에 대한 원의 방정식을 행렬식에 대입하여, centerx, centery,  radius를 구함..
centerx = -(det(-d)/(2 * det(a)))
centery = -(det(e)/(2 * det(a)))
radius = sqrt( ((det(-d)^2 + det(e)^2) / (4*det(a)^2)) - (det(-f)/det(a)))

[/code]

[code cpp]
//OpenCV
//관심점 3점을 원의 방적식에 대입
float x1, x2, x3, y1, y2, y3;
x1 = 5; y1 = 0;
x2 = 10; y2 = 5;
x3 = 5; y3  = 10;

//a= [x1, y1, 1; x2, y2, 1; x3, y3, 1];
CvMat *matA = cvCreateMat(3, 3, CV_32F); // a행렬
cvmSet(matA, 0, 0, x1);
cvmSet(matA, 0, 1, y1);
cvmSet(matA, 0, 2, 1);
cvmSet(matA, 1, 0, x2);
cvmSet(matA, 1, 1, y2);
cvmSet(matA, 1, 2, 1);
cvmSet(matA, 2, 0, x3);
cvmSet(matA, 2, 1, y3);
cvmSet(matA, 2, 2, 1);

//d = [x1^2+y1^2, y1, 1; x2^2+y2^2, y2, 1; x3^2+y3^2, y3, 1];
CvMat *matD = cvCreateMat(3, 3, CV_32F); // d행렬
cvmSet(matD, 0, 0, x1*x1+y1*y1);
cvmSet(matD, 0, 1, y1);
cvmSet(matD, 0, 2, 1);
cvmSet(matD, 1, 0, x2*x2+y2*y2);
cvmSet(matD, 1, 1, y2);
cvmSet(matD, 1, 2, 1);
cvmSet(matD, 2, 0, x3*x3+y3*y3);
cvmSet(matD, 2, 1, y3);
cvmSet(matD, 2, 2, 1);

//e = [x1^2+y1^2, x1, 1; x2^2+y2^2, x2, 1; x3^2+y3^2, x3, 1];
CvMat *matE = cvCreateMat(3, 3, CV_32F); // e행렬
cvmSet(matE, 0, 0, x1*x1+y1*y1);
cvmSet(matE, 0, 1, x1);
cvmSet(matE, 0, 2, 1);
cvmSet(matE, 1, 0, x2*x2+y2*y2);
cvmSet(matE, 1, 1, x2);
cvmSet(matE, 1, 2, 1);
cvmSet(matE, 2, 0, x3*x3+y3*y3);
cvmSet(matE, 2, 1, x3);
cvmSet(matE, 2, 2, 1);

//f = [x1^2+y1^2, x1, y1; x2^2+y2^2, x2, y2; x3^2+y3^2, x3, y3];
CvMat *matF = cvCreateMat(3, 3, CV_32F); // f행렬
cvmSet(matF, 0, 0, x1*x1+y1*y1);
cvmSet(matF, 0, 1, x1);
cvmSet(matF, 0, 2, y1);
cvmSet(matF, 1, 0, x2*x2+y2*y2);
cvmSet(matF, 1, 1, x2);
cvmSet(matF, 1, 2, y2);
cvmSet(matF, 2, 0, x3*x3+y3*y3);
cvmSet(matF, 2, 1, x3);
cvmSet(matF, 2, 2, y3);

//중심점과 반지름 구하기
float centerx, centery, radius;
centerx = -( -cvDet(matD) / (2 * cvDet(matA)) );
centery = -( cvDet(matE) / (2 * cvDet(matA)) );
radius = sqrt( (((-cvDet(matD) * -cvDet(matD)) + (cvDet(matE)*cvDet(matE))) / (4*(cvDet(matA)*cvDet(matA)))) - (-cvDet(matF)/cvDet(matA)) );

//메모리 해제
cvReleaseMat(&matA);
cvReleaseMat(&matD);
cvReleaseMat(&matE);
cvReleaseMat(&matF);

[/code]