구글에드센스


Random sample consensus (RANSAC) - 사용중인 SRC - by Skai

사용중인 Lib가 Mil 이기 때문에 Mil 버전으로 업데이트 합니다.
뭐 필요한부분은 많은 분들이 가져다 사용하셔도 무방합니다.
개인적인 바램이 있다면, 덧글 혹은 메일이라도 주셨으면 합니다.

///////////////////////////////////////////////////
///RanSac 
typedef struct {
double mx, my; // x,y 기울기 Radian.
double cx, cy; // 구한 직선의 중심점.
double dTheta;
}RansacLine;


bool find_in_samples(CPoint *samples, int no_samples, CPoint *data);
void get_samples(CPoint *samples, int no_samples, CPoint *data, int no_data);
int compute_model_parameter(CPoint samples[], int no_samples, RansacLine &model);
double compute_distance(RansacLine &line, CPoint &x);
double model_verification(CPoint *inliers, int *no_inliers, RansacLine &estimated_model, CPoint *data, int no_data, double distance_threshold);
double ransac_line_fitting(CPoint *data, int no_data, RansacLine &model, double distance_threshold);
double RanSac_Align(MIL_ID pSrc, RansacLine& line, int distance_threshold, CPoint* cpLT, CPoint* cpRT, int nPos = Align_X);



bool CAOT_Algorithm::find_in_samples(CPoint *samples, int no_samples, CPoint *data)
{
for (int i = 0; i < no_samples; ++i) {
if (samples[i].x == data->x && samples[i].y == data->y) {
return true;
}
}
return false;
}

void CAOT_Algorithm::get_samples(CPoint *samples, int no_samples, CPoint *data, int no_data)
{
// 데이터에서 중복되지 않게 N개의 무작위 셈플을 채취한다.
for (int i = 0; i < no_samples;) {
int j = rand() % no_data;

if (!find_in_samples(samples, i, &data[j])) {
samples[i] = data[j];
++i;
}
};
}

int CAOT_Algorithm::compute_model_parameter(CPoint samples[], int no_samples, RansacLine &model)
{
// PCA 방식으로 직선 모델의 파라메터를 예측한다.

double sx = 0, sy = 0;
double sxx = 0, syy = 0;
double sxy = 0, sw = 0;

for (int i = 0; i < no_samples; ++i)
{
double x = samples[i].x;
double y = samples[i].y;

sx += x;
sy += y;
sxx += x*x;
sxy += x*y;
syy += y*y;
sw += 1;
}

//variance;
double vxx = (sxx - sx*sx / sw) / sw;
double vxy = (sxy - sx*sy / sw) / sw;
double vyy = (syy - sy*sy / sw) / sw;

//principal axis
double theta = atan2(2 * vxy, vxx - vyy) / 2;
model.dTheta = theta;

model.mx = cos(theta);
model.my = sin(theta);

//center of mass(cx, cy)
model.cx = sx / sw;
model.cy = sy / sw;

//직선의 방정식: sin(theta)*(x - sx) = cos(theta)*(y - sy);
return 1;
}

double CAOT_Algorithm::compute_distance(RansacLine &line, CPoint &x)
{
// 한 점(x)로부터 직선(line)에 내린 수선의 길이(distance)를 계산한다.

return fabs((x.x - line.cx)*line.my - (x.y - line.cy)*line.mx) / sqrt(line.mx*line.mx + line.my*line.my);
}

double CAOT_Algorithm::model_verification(CPoint *inliers, int *no_inliers, RansacLine &estimated_model, CPoint *data, int no_data, double distance_threshold)
{
*no_inliers = 0;

double cost = 0.;

for (int i = 0; i < no_data; i++){
// 직선에 내린 수선의 길이를 계산한다.
double distance = compute_distance(estimated_model, data[i]);

// 예측된 모델에서 유효한 데이터인 경우, 유효한 데이터 집합에 더한다.
if (distance < distance_threshold) {
cost += 1.;

inliers[*no_inliers] = data[i];
++(*no_inliers);
}
}

return cost;
}


double CAOT_Algorithm::ransac_line_fitting(CPoint *data, int no_data, RansacLine &model, double distance_threshold)
{
const int no_samples = 2;

if (no_data < no_samples) {
return 0.;
}

CPoint *samples = new CPoint[no_samples];

int no_inliers = 0;
CPoint *inliers = new CPoint[no_data];

RansacLine estimated_model;
double max_cost = 0.;

int max_iteration = (int)(1 + log(1. - 0.99) / log(1. - pow(0.5, no_samples)));

for (int i = 0; i < max_iteration; i++) {
// 1. hypothesis

// 원본 데이터에서 임의로 N개의 셈플 데이터를 고른다.
get_samples(samples, no_samples, data, no_data);

// 이 데이터를 정상적인 데이터로 보고 모델 파라메터를 예측한다.
compute_model_parameter(samples, no_samples, estimated_model);

// 2. Verification

// 원본 데이터가 예측된 모델에 잘 맞는지 검사한다.
double cost = model_verification(inliers, &no_inliers, estimated_model, data, no_data, distance_threshold);

// 만일 예측된 모델이 잘 맞는다면, 이 모델에 대한 유효한 데이터로 새로운 모델을 구한다.
if (max_cost < cost) {
max_cost = cost;

compute_model_parameter(inliers, no_inliers, model);
}
}

delete[] samples;
delete[] inliers;

return max_cost;
}

double CAOT_Algorithm::RanSac_Align(MIL_ID pSrc, RansacLine& line, int distance_threshold, CPoint* cpLT, CPoint* cpRT, int nPos)
{
/////////////////////////////////////////////////////////////////////////
/// 1. 입력영상은 무조건 이진 영상만 입력해야합니다.
/// 2. 입력영상은 전체 영상이 아닌 상단만 잘라서 입력 해주면 좋습니다.

double ReturnRadian = -1.;
register int y = 0, x = 0, i = 0, j = 0;
int nWidth = MbufInquire(pSrc, M_SIZE_X, M_NULL);
int nHeight = MbufInquire(pSrc, M_SIZE_Y, M_NULL);
CPoint *data = NULL;
//RansacLine line;

///////////////////////////////////////////////////////////////
/////////////
// 영상이 까맣게 들어오면 무한루프 빠짐. 예외처리 합니다.
MIL_INT AvgValue = 0;
MIL_ID ImageResult = MimAllocResult(clsMil.m_MilSystem_Penetrate, 1, M_STAT_LIST, M_NULL);
MimStat(pSrc, ImageResult, M_MEAN, M_NULL, M_NULL, M_NULL);
MimGetResult(ImageResult, M_MEAN + M_TYPE_MIL_INT, &AvgValue);
MimFree(ImageResult);

if (AvgValue == 0)
return ReturnRadian;
/////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////
BYTE* nBuf = new BYTE[nWidth * nHeight];
memset(nBuf, 0, nWidth * nHeight * sizeof(BYTE));
MbufGet2d(pSrc, 0, 0, nWidth, nHeight, nBuf);
//////////////////////////////////////////////////

double cost = 0.0;
if (nPos == Align_X)
{
data = new CPoint[nWidth];
for (x = 0; x < nWidth; x++)
{
for (y = 0; y < nHeight; y++)
{
if (nBuf[IMGARR1(y, x, nWidth)] == 255)
{
data[x].x = x;
data[x].y = y;
break;
}
else
{
data[x].x = 0;
data[x].y = 0;
}
}
}
cost = ransac_line_fitting(data, nWidth, line, distance_threshold);
}
else
{
data = new CPoint[nHeight];
for (y = 0; y < nHeight; y++)
{
for (x = 0; x < nWidth; x++)
{
if (nBuf[IMGARR1(y, x, nWidth)] == 255)
{
data[y].x = x;
data[y].y = y;
break;
}
else
{
data[y].x = 0;
data[y].y = 0;
}

}
}

cost = ransac_line_fitting(data, nHeight, line, distance_threshold);
}
cpLT->x = data[0].x;
cpLT->y = data[0].y;
if (nPos == Align_X)
{
ReturnRadian = (RADIAN_TO_DEGREE(line.my));
}
else
{
ReturnRadian = (RADIAN_TO_DEGREE(line.mx));
}
return ReturnRadian;
}



덧글

댓글 입력 영역