3점을 지나는 원은 하나만 만들어집니다.
물론 3점이 직선이 아닐때의 얘기지요
옛날 프로젝트 작업할때 필요해서 만들었던 겁니다.
버그가 없다고 장담할 수는 없지만 그냥 갖다 써도
큰 문제는 없을듯 하네요
버그 찾아주시면 더 좋고..
소스에 주석이 들어있으니 사용하는 방법은 한번 보면 알 수 있습니다.
이런 거 필요한 사람이 없을줄 알았는데..
델파이 사이트에서 누가 필요하다고 올린 글을 보고 생각나서 여기에 올려둡니다.
델마당에는 가입을 안해서 글쓰기가 안되네요.. 가입하기도 귀찮구.. ^^;
#include <stdio.h>
#include "XMathUtil.h"
#ifndef M_PI
#define M_PI 3.1415926535
#endif
// 행렬 치환 명령 함수
// 행을 d1 * f1 + d2 * f2 로 계산해서 d1행렬에 입력
bool ChangeMat(double d1[], double f1 , double d2[], double f2)
{
for(int i=0; i<MAT_COL; i++)
{
d1[i] = f1 * d1[i] + f2 * d2[i];
}
return true;
}
//--------------------------------------------------------------------
// 3X4 행렬의 가우스 소거법으로 근을 구하는 함수
//--------------------------------------------------------------------
#define MIN_DIVIDED_VALUE 0.000001
#define MIN_SLOPE 0.0001
bool CalcGauss(double data[][MAT_COL])
{
int i, j, k;
double mat[MAT_ROW][MAT_COL];
// 세점의 경사도 측정 하여 한 직선위에 있을 경우 에러 리턴
// ( y2-y1) / ( x2-x1) == ( y3 - y2 ) / (x3 - x2)
// --> ( y2 - y1 ) * ( x3 - x2 ) == ( y3 - y2 ) * ( x2 - x1 )
double y2_y1 = data[1][1]-data[0][1];
double x3_x2 = data[2][0]-data[1][0];
double y3_y2 = data[2][1]-data[1][1];
double x2_x1 = data[1][0]-data[0][0];
if( fabs( y2_y1 * x3_x2 - y3_y2 * x2_x1 ) < MIN_SLOPE )
{
return false;
}
//----------------------------------------------------------------------------
// 대각선행렬값(1,1),(2,2),(3,3)에 0이 아닌 값이 배열되도록 한다.
// Gauss 계산 가능한 mat 행렬 만들기
for(k=0; k<MAT_COL-1; k++)
{
for(i=0; i<MAT_ROW; i++)
{
if(data[i][k] != 0.0)
{
// M(i, k) 행열의 값이 0.0이 아닌 열의 데이터 복사
// mat에 데이터 입력후 data원소 0.0으로 초기화하여 중복를 방지함
for(j=0; j<MAT_COL; j++)
{
mat[k][j] = data[i][j];
data[i][j] = 0.0;
}
break;
}
}
// 에러처리
// k열에 0이 아닌 값이 없으면 에러 처리함
// k가 0일때는 x좌표가 모두 0.0 이므로 직선이고
// k가 1일때는 y좌표가 모두 0.0인 것이므로 세점은 윈위의 좌표가 아님
if(i == MAT_ROW) // 에러
return false;
}
//----------------------------------------------------------------------------
// 0.0 으로 나누는 문제 에러 처리
if( fabs(mat[0][0]) < MIN_DIVIDED_VALUE )
{
return false;
}
// 첫번째 열을
// ( 1 , ? , ? )
// ( 0 , ? , ? )
// ( 0 , ? , ? )
// 으로 만듦
ChangeMat(mat[0], 1.0 / mat[0][0], mat[0] , 0.0);
ChangeMat(mat[1], 1.0, mat[0] , -mat[1][0]);
ChangeMat(mat[2], 1.0, mat[0] , -mat[2][0]);
//----------------------------------------------------------------------------
// 0.0 으로 나누는 문제 에러 처리
if( fabs(mat[1][1]) < MIN_DIVIDED_VALUE )
{
return false;
}
// 두번째 열을
// ( 1 , 0 , ? )
// ( 0 , 1 , ? )
// ( 0 , 0 , ? )
// 으로 만듦
ChangeMat(mat[1], 1.0 / mat[1][1], mat[1] , 0.0);
ChangeMat(mat[0], 1.0, mat[1] , -mat[0][1]);
ChangeMat(mat[2], 1.0, mat[1] , -mat[2][1]);
//----------------------------------------------------------------------------
// 0.0 으로 나누는 문제 에러 처리
if( fabs(mat[2][2]) < MIN_DIVIDED_VALUE )
{
return false;
}
// 세번째 열을
// ( 1 , 0 , 0 )
// ( 0 , 1 , 0 )
// ( 0 , 0 , 1 )
// 으로 만듦
ChangeMat(mat[2], 1.0 / mat[2][2], mat[2] , 0.0);
ChangeMat(mat[0], 1.0, mat[2] , -mat[0][2]);
ChangeMat(mat[1], 1.0, mat[2] , -mat[1][2]);
for(i=0; i<MAT_ROW; i++)
for(j=0; j<MAT_COL; j++)
data[i][j] = mat[i][j];
return true;
}
//--------------------------------------------------------------------
// 3점을 이용해서 원의 방정식을 구하는 함수
//--------------------------------------------------------------------
// x^2 + y^2 + Ax + By + C = 0 원의 방정식 일반형
// 여기서 A, B , C의 계수를 구해서 가우스 소거법 적용
// A의 계수 : x
// B의 계수 : y
// C의 계수 : 1
// D의 계수 : -x^2 - y^2
// A,B,C값의 근은 행렬의 네번째 열에 저장됨
//--------------------------------------------------------------------
// A,B,C로 원의 방정식 구하기
//--------------------------------------------------------------------
// 원의 중심 ( cx, cy )
// cx = - A / 2
// cy = - B / 2
// 반지름 R = Sqrt( -C - (A^2) / 4 - (B^2) / 4 )
// angle 포인터에 시점, 중점, 종점의 각도 입력
bool Get3PointCircle(/*IN*/double *x, /*IN*/double *y
, /*OUT*/double &cx, /*OUT*/double &cy, /*OUT*/double &r, double *angle)
{
double data[MAT_ROW][MAT_COL];
for (int i=0; i<MAT_ROW; i++)
{
data[i][0] = x[i];
data[i][1] = y[i];
data[i][2] = 1.0;
data[i][3] = - x[i]*x[i] - y[i]*y[i];
}
if( CalcGauss(data) == false)
return false;
double A, B, C;
A = data[0][3];
B = data[1][3];
C = data[2][3];
cx = - A / 2.0;
cy = - B / 2.0;
r = sqrt( -C + (A*A /4.0) + (B*B/4.0) ) ;
if(angle!=NULL)
{
double dx, dy;
for (int i=0; i<MAT_ROW; i++)
{
dx = x[i] - cx;
dy = y[i] - cy;
angle[i] = atan2(dy, dx) * 180.0 / M_PI;
}
}
return true;
}
//---------------------------------------------------------------------------
/*
bool Get3PointCircle( POINT &p1, POINT &p2, POINT &p3,
double &cx, double &cy, double &r, double &angle)
{
static const size_t lDimension = 2;
double data[lDimension][MAT_ROW];
memset( data, 0, sizeof(data) );
data[0][0] = p1.x;
data[1][0] = p1.y;
data[0][1] = p2.x;
data[1][1] = p2.y;
data[0][2] = p3.x;
data[1][2] = p3.y;
return Get3PointCircle( data[0], data[2], cx, cy, r, &angle );
}
*/
//--------------------------------------------------------------------
// 중심좌표(cx, cy) 시작점 (sx, sy)와 각도(angle)를 입력 하면
// ARC의 마지막 좌표 (ex, ey)를 계산해 주는 함수
//--------------------------------------------------------------------
bool GetArcEndPos(/*IN*/double cx, /*IN*/double cy, /*IN*/double sx, /*IN*/double sy
, /*IN*/double angle
, /*OUT*/double &ex, /*OUT*/double &ey )
{
if(cx == sx && cy==sy)
return false;
double dx = sx - cx;
double dy = sy - cy;
// 반지름 구하기
double r = sqrt( dx*dx + dy*dy);
//str.sprintf("반지름 : %2.4lf", r);
//Form1->Memo1->Lines->Add(str);
// 시작점의 각도
double sangle = atan2(dy, dx) * 180.0 / M_PI;
//str.sprintf("시작점 각도 : %2.4lf", sangle);
//Form1->Memo1->Lines->Add(str);
// 마지막 점의 각도, 좌표 계산
double eAngle = sangle + angle;
ex = r * cos( eAngle * M_PI / 180.0 ) + cx;
ey = r * sin( eAngle * M_PI / 180.0 ) + cy;
//str.sprintf("종점 좌표 (%2.4lf, %2.4lf) , 종점 각도 : %2.4lf", ex, ey, eAngle);
//Form1->Memo1->Lines->Add(str);
char buf[100];
sprintf(buf, "kkol r=%2.3lf, sangle=%2.3lf, eangle=%2.3lf, end(%2.3lf , %2.3lf)",
r, sangle, eAngle, ex, ey);
return true;
}
//--------------------------------------------------------------------
// 중심좌표(cx, cy) 시작점 (sx, sy)와 각도(angle)를 입력 하면
// ARC의 쭝간 좌표 (ex, ey)를 추출해주는 주는 함수
//--------------------------------------------------------------------
bool GetArcMidPos(/*IN*/double cx, /*IN*/double cy, /*IN*/double sx, /*IN*/double sy
, /*IN*/double angle
, /*OUT*/double &ex, /*OUT*/double &ey )
{
return GetArcEndPos( cx, cy, sx, sy, angle/2, ex, ey );
}
//----------------------------------------------------------------------------
// cx, cy를 중심으로 angle 만큼 회전 변환 처리
void CalcRotation(double cx, double cy, double &x, double &y, double angle)
{
double ox = x-cx;
double oy = y-cy;
// 라디안으로 변환
double r = angle * M_PI / 180.0 ;
x = cos(r) * ox - sin(r) * oy;
y = sin(r) * ox + cos(r) * oy;
x += cx;
y += cy;
}
//----------------------------------------------------------------------------
void swap(double &a, double &b)
{
double x;
x = a;
a = b;
b = x;
}
//----------------------------------------------------------------------------
// 두점이 만드는 직선의 각도 리턴
#define MOTION_TOLERANCE (0.01)
double CalcAngle(double x1, double y1, double x2, double y2)
{
double angle;
// 데이터 정열 (x1, y1)이 더 작은 값을 가지도록 정열
if(x1 > x2 + MOTION_TOLERANCE)
{
swap(x1, x2);
swap(y1, y2);
}
else if(fabs(x1 - x2) < MOTION_TOLERANCE) // x1, x2가 비슷하면...
{
if(y1 > y2 + MOTION_TOLERANCE)
{
swap(x1, x2);
swap(y1, y2);
}
}
double dx = x2 - x1;
double dy = y2 - y1;
// 두점이 같다면..에러
if(fabs(dx)< MOTION_TOLERANCE && fabs(dy)<MOTION_TOLERANCE )
return 0.0;
// 각도 계산
angle = atan2(dy , dx) * 180 / M_PI;
return angle;
}
//----------------------------------------------------------------------------