[ RANSAC ] 타원 피팅 (fit ellipse using ransac)

또 등장하는 제 마우스 패드와 선물 받은 큐브를 통해서 

가려진 타원을 어떻게 복원할지에 대해서 설명을 드리고자 합니다.

일반적으로 circle 과 기존의 ellipse fitting 으로는 아래 마우스 패드의 외각선을 추출해 진행해보면 이상한 결과를 얻을 수 있습니다. 물론,, 아래 마우스 패드의 외각선 정보에 이상한 그림이 겹쳐져 잘못된 정보를 포함해 버렸기 때문이죠 ^^

* fit ellipse using RANSAC (random sample consensu)

RANSAC(random sample consensu) 이란 ?

그림의 마우스 패드를 검출 하려고 보면 외각선 정보에 큐브 외각선이 겹쳐 나타나게 됩니다.  어떻게 하면 이 잘못된 정보 (outlier) 를 제외하고 피팅을 진행 할 수 있을까요? RANSAC은 데이터 집합 중에 일정확률을 가지고 샘플을 추려서 모델에 맞춰보고 해당 샘플 중에 가장 데이터 집합을 많이 대변하는 모델을 추려내는 작업이라고 생각할 수 있습니다.

* 말로 설명하려니 어렵네요. 아래 한번 RANSAC 에 대한 접근 방법을 보시면 바로 이해하실 수 있습니다.

모델 :  타원(rotated ellipse)

데이터 :  마우스 패드 외각선(Inlier) 와 큐브 외각선(outlier)이 포함된 외각선 집합

1. 일정 확률을 가지고 DATA 집합에서 임의 추출 진행

2. 추출진 sample 을 통해서 모델을 예측

3. 예측 모델에서 어느정도의 offset 범위 안에 들어오는 기존보다 클 경우, 해당 (x,y) 값들을 저장

4. 1 - 3 을 N 번 반복 수행

5. 최종 저장된 (X,Y) 값을 통해서 fitting 진행

import cv2
import os,re,sys
import numpy as np


def fit_rotated_ellipse_ransac(data,iter=30,sample_num=10,offset=80.0):

	count_max = 0
	effective_sample = None

	for i in range(iter):
		sample = np.random.choice(len(data), sample_num, replace=False)

		xs = data[sample][:,0].reshape(-1,1)
		ys = data[sample][:,1].reshape(-1,1)

		J = np.mat( np.hstack((xs*ys,ys**2,xs, ys, np.ones_like(xs,dtype=np.float))) )
		Y = np.mat(-1*xs**2)
		P= (J.T * J).I * J.T * Y

		# fitter a*x**2 + b*x*y + c*y**2 + d*x + e*y + f = 0
		a = 1.0; b= P[0,0]; c= P[1,0]; d = P[2,0]; e= P[3,0]; f=P[4,0];
		ellipse_model = lambda x,y : a*x**2 + b*x*y + c*y**2 + d*x + e*y + f

		# threshold 
		ran_sample = np.array([[x,y] for (x,y) in data if np.abs(ellipse_model(x,y)) < offset ])

		if(len(ran_sample) > count_max):
			count_max = len(ran_sample) 
			effective_sample = ran_sample

	return fit_rotated_ellipse(effective_sample)


def fit_rotated_ellipse(data):

	xs = data[:,0].reshape(-1,1) 
	ys = data[:,1].reshape(-1,1)

	J = np.mat( np.hstack((xs*ys,ys**2,xs, ys, np.ones_like(xs,dtype=np.float))) )
	Y = np.mat(-1*xs**2)
	P= (J.T * J).I * J.T * Y

	a = 1.0; b= P[0,0]; c= P[1,0]; d = P[2,0]; e= P[3,0]; f=P[4,0];
	theta = 0.5* np.arctan(b/(a-c))  
	
	cx = (2*c*d - b*e)/(b**2-4*a*c)
	cy = (2*a*e - b*d)/(b**2-4*a*c)

	cu = a*cx**2 + b*cx*cy + c*cy**2 -f
	w= np.sqrt(cu/(a*np.cos(theta)**2 + b* np.cos(theta)*np.sin(theta) + c*np.sin(theta)**2))
	h= np.sqrt(cu/(a*np.sin(theta)**2 - b* np.cos(theta)*np.sin(theta) + c*np.cos(theta)**2))

	ellipse_model = lambda x,y : a*x**2 + b*x*y + c*y**2 + d*x + e*y + f

	error_sum = np.sum([ellipse_model(x,y) for x,y in data])
	print('fitting error = %.3f' % (error_sum))

	return (cx,cy,w,h,theta)

def main(img_path):

	color_list = [(238,0,0),(0,252,124),(142,56,142),(10,20,0),(245,245,245)]

	src = cv2.imread(img_path, cv2.IMREAD_COLOR)
	gray = cv2.cvtColor(src,cv2.COLOR_RGB2GRAY)
	retv, th = cv2.threshold(gray,0,255,cv2.THRESH_BINARY + cv2.THRESH_OTSU) 

	color_th = cv2.cvtColor(th,cv2.COLOR_GRAY2RGB)

	_, contours , _ = cv2.findContours(th,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)
	zeros = np.zeros(np.shape(th),dtype=np.uint8)

	for con in contours:
		approx = cv2.approxPolyDP(con, 0.01 * cv2.arcLength(con,True),True)
		area = cv2.contourArea(con)
		if(len(approx) > 10 and area > 30000):

			# fit ellipse
			cx,cy,w,h,theta = fit_rotated_ellipse(con.reshape(-1,2))
			cv2.ellipse(src,(int(cx),int(cy)),(int(w),int(h)),theta*180.0/np.pi,0.0,360.0,color_list[2],2)

			# fit ellipse using ransac
			cx,cy,w,h,theta = fit_rotated_ellipse_ransac(con.reshape(-1,2))
			cv2.ellipse(src,(int(cx),int(cy)),(int(w),int(h)),theta*180.0/np.pi,0.0,360.0,color_list[0],2)

	cv2.imwrite('out.jpg',src)

if __name__ == '__main__':
	if(len(sys.argv) != 2):
		print('usage : {0} <image_abs_path>'.format(sys.argv[0]))
		exit(0)
	main(sys.argv[1])