RANSAC(RANdom SAmple Consensus)

1980년도 원전 논문

RANDOM SAMPLE CONSENSUS : a Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography (http://www.dtic.mil/dtic/tr/fulltext/u2/a460585.pdf)

RANSAC(RANdom SAmple Consensus)을 이용한 Line Fitting Example

http://blog.daum.net/pg365/62

RANSAC 은 노이즈가 심한 원본 데이터로부터 모델 파라미터를 예측하는 방법이다.

위의 링크에는 line fitting을 예제로 하였으나, 실제로 cost와 모델링하는 함수만 만든다면

2차 곡선, 3차 곡선, 기타 모델링을 수행할 수 있다.

위의 샘플은 클래스로 구현되어 있지 않아, 클래스로 변경하였다.

기본 알고리즘을 보면,

1. 주어진 전체 샘플에서 상당히 적은 샘플을 수집한다.

2. 주어진 샘플로 모델 파라미터를 구해본다.

3. 전체 모델 집합에서 2에서 구한 모델 파라미터로 피팅해봤을때 일정 범위 안에 만족하는 것만(Inlier)  수집한다.

4. 3의 과정에서 전체 집합이 2에서 구한 파마미터에 얼마나 잘 fitting되는지 cost를 구하게 된다. 만약 구한 cost가 이전에 구한 cost보다 크다면 5번을 수행하고 아니면 다시 1번을 수행한다. ( 여기서 cost라는 것은 약간 혼동이 될 수 있는데, 의미적으로 보면 cost라기 보단, 얼마나 잘 적합되었는지를 나타낸다. 즉 값이 크면 좋은 모델이라는 것, 코딩에 따라서 반대로도 가능하겠지요...)

5. 수집된 Inlier를 가지고 다시 모델 파라미터를 구한다. 구한 파라미터는 저장하고, 너무 많이 시도했다면 끝내고 아니면 1번부터 다시 수행한다.

 

아래 알고리즘을 이용해서 손 영역의 방향을 예측한 결과

 

#include "StdAfx.h"
#include "RANSAC.h"


CRANSAC::CRANSAC(void)
{
}


CRANSAC::~CRANSAC(void)
{

}

double CRANSAC::GetLineParameter(CSvImage& MaskImage, sLine &model, double distance_threshold)
{
	vector<CPoint> vPoint;

	for(int y=0; y<MaskImage.GetHeight(); y++)
	{
		for(int x=0; x< MaskImage.GetWidth(); x++)
		{
			if(MaskImage.GetPixel(x,y) == HANDCOLOR)
			{
				vPoint.push_back(CPoint(x,y));
			}
		}
	}

	return Ransac_Line_Fitting(vPoint, model, distance_threshold);

}
double CRANSAC::Ransac_Line_Fitting(vector<CPoint>& vPoint, sLine &model, double distance_threshold)
{
	const int no_samples = 2;

	if (vPoint.size() < no_samples) {
		return 0.;
	}

	vector<CPoint> vSamplePoint;
	vector<CPoint> vInliersPoint;
	int no_inliers = 0;
	

	sLine 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 (vSamplePoint, no_samples, vPoint);

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

		// 2. Verification

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

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

			compute_model_parameter (vInliersPoint, model);
		}
	}

	return max_cost;
}


void CRANSAC::get_samples (vector<CPoint>& vSamplePoint,int nNumberofSample, vector<CPoint>& vPoint)
{
	// 데이터에서 중복되지 않게 N개의 무작위 셈플을 채취한다.
	for (int i=0; i<nNumberofSample; ) 
	{
		int j = rand()%vPoint.size();

		if (!find_in_samples(vSamplePoint, &vPoint.at(j))) {
			vSamplePoint.push_back(vPoint.at(j));
			++i;
		}
	};
}

bool CRANSAC::find_in_samples (vector<CPoint>& vSamplePoint, CPoint*data)
{
	for (int i=0; i<vSamplePoint.size(); ++i) {
		if (vSamplePoint.at(i).x == data->x && vSamplePoint.at(i).y == data->y) {
			return true;
		}
	}
	return false;
}

int CRANSAC::compute_model_parameter(vector<CPoint>& vInlierPoint, sLine &model)
{
	// PCA 방식으로 직선 모델의 파라메터를 예측한다.

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

	for(int i = 0; i<vInlierPoint.size();++i)
	{
		double x = (double)vInlierPoint.at(i).x;
		double y = (double)vInlierPoint.at(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.mx = cos(theta);
	model.my = sin(theta);

	//center of mass(xc, yc)
	model.sx = sx/sw;
	model.sy = sy/sw;

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

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

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

double CRANSAC::model_verification (vector<CPoint>& vinliers, sLine &estimated_model, vector<CPoint>& vPoint, double distance_threshold)
{
	

	double cost = 0.;

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

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

			vinliers.push_back(vPoint.at(i));
			
			
		}
	}

	return cost;
}