[ 최소자승법 ] 원, 타원 측정

* 이글은 다크프로그래머 님의 - [최소자승법의 이해] 를 읽고 필요한 부분을 정리 및 제가 아는 내용을 추가하였습니다.

http://darkpgmr.tistory.com/56

** Least Square Method, fit ellipse, fit circle, python opencv 

가끔 영상을 보면서 원과 타원에 대해서 어떠한 정보를 얻어야 될 때가 있습니다. 

외각선 정보보다, 몇개의 파라미터 값(중심, 사이즈 ...)만 알면 표현이 훨씬 간단할 뿐만 아니라 각 물체(원, 타원, 사각) 간의 관계 표현할 때 이만큼 좋은 정보가 없습니다.

그럼 세부적으로 중심, 지름, 장축, 단축 등 이런 정보들은 어떻게 얻을 수 있을 까요?

이럴 때 많이 사용하는 방식이 피팅(fitting)을 사용해서 원의 방정식 또는 타원의 방정식에 얻어진 외각선( x,y 위치 정보)에 데이터가 얼마나 잘 맞는지 넣어보면서 파라미터(장축,단축, 중심) 를 찾는 방향으로 알 수 있습니다.

영상내 어떤 둥그런 물체(object) 를 필터링 해서 크기가 얼마인지 알고 싶을 때, 

제가 아끼는 마우스패드가 쇼파 위에 이쁘게 놓여있네요, 이럴 때 이 마우스 패드의 중심과 크기를 알려고 하면 간단히 opencv 의 fitEllipse 함수를 사용해도 큰 문제는 없지만, 기본 개념부터 구현해 보았습니다.

□원의 방정식

답을 찾아가는 방법으로 pesudo inverse 방식을 쓰려고 양함수( f(x,y) = z) 모델을 이용하여 수식을 정리하면 아래와 같이 a,b,c 라는 파라미터를 구할 수 있습니다.

$$ \begin{matrix} { (x-x_1)^2 } +  { (y-y_1)^2  } &=& r^2 \\ \\  \Rightarrow x^2+y^2 -2ax - 2by + a^2 + b^2 - r^2 &=& 0 \\ \\ \Rightarrow x^2+y^2 -2ax - 2by + c = 0 &,& c = a^2 + b^2 - r^2 \\ \\ \Rightarrow  -2ax - 2by + c &=& -x^2-y^2 \end{matrix} $$

$$ \begin{matrix} \begin{bmatrix} -2x_1& -2y_1& 1\\ \vdots& \vdots&\vdots  \\ -2x_n& -2y_n& 1 \end{bmatrix}  \begin{bmatrix} a \\ b \\ c \end{bmatrix}  = \begin{bmatrix} -x_1^2-y_1^2 \\ \vdots \\ -x_n^2-y_n^2 \end{bmatrix} \\  \\  J =\begin{bmatrix} -2x_1& -2y_1& 1 \\ \vdots& \vdots&\vdots  \\ -2x_n& -2y_n& 1 \end{bmatrix}, X = \begin{bmatrix} a \\ b \\ c \end{bmatrix}, Y = \begin{bmatrix} -x_1^2-y_1^2 \\ \vdots \\ -x_n^2-y_n^2 \end{bmatrix} \\ \\ J X = Y \\  X = {(J^T J)}^{-1} J^T Y \\ \\ \end{matrix} $$ 

□ 타원의 방정식

중심이 0,0 위치에서의 타원의 방정식

$$ \Large \begin{matrix}{ x^2 \over w^2 } +  { y^2 \over h^2 } = 1  \end{matrix} $$ 

원은 회전된 경우(rotated) 문제 없지만, 타원의 경우는 중요하기 때문엥 이를 방정식에 반영 해주어야 합니다.

접근 방법을 

1. roatated ellipse (회전이 되어 있는 타원) 

2. rotated ellipse 의 회전 각도를 반대로 적용하여 타원의 방정식 폼에 맞게 수정

3. 타원의 중심이 이동 할 수 있음으로 타원의 방정식에 반영

CCW 로 회전이 되어 있다는 가정하에  복원을 위해 CW 방향으로 회전 변환을 적용하겠습니다.

$$ \begin{matrix}   \begin{bmatrix} x' \\ y' \end{bmatrix}  = \begin{bmatrix} cos(\theta) & sin(\theta) \\ -sin(\theta) & cos(\theta) \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} \\ \\ \Rightarrow x' = cos(\theta)x + sin(\theta)y \\ \Rightarrow y' = -sin(\theta)x + cos(\theta)y \end{matrix} $$ 

최종적으로 x1,y1 의 중심에 회전이 반영된 방정식입니다.

$$ { (x − x_1 )^2 \over a^2 } +  { (y − y_1 )^2 \over b^2} = 1 \\  { [cos(\theta)x  + sin(\theta)y − x_1 ]^2 \over a^2 } +  { [−sin(\theta)x + cos(\theta)y −y_1]^2 \over b^2} = 1 $$ 

위 수식을 풀어서 계수를 정리하면 다음과 같습니다.  여기에 간단하게 회전각도에 대해서만 유도를 해보면

( 추가된 내용)

$$  x^2[b^2cos^2(\theta) + a^2sin^2(\theta)] + xy[2b^2sin(\theta)cos(\theta) − 2a^2sin(\theta)cos(\theta)] +y^2[b^2sin^2(\theta)+a^2cos^2(\theta)] +\cdots = a^2b^2 \\ \\ c_1 x^2 + c_2 xy +c_3 y^2 +c_4 x + c_5 y + c6 = 0 $$

$$ \begin{array}{l} { c2 \over (c1 − c3) } = { 2b^2sin\theta cos\theta − 2a^2sin\theta cos\theta \over b^2cos^2\theta + a^2sin^2\theta −  b^2cos^2\theta −  a^2sin^2\theta} = {2sin\theta cos\theta(b^2 − a^2) \over(cos^2\theta − sin^2\theta)(b^2 − a^2)}  \\ = { 2sin\theta cos\theta \over cos^2\theta \; sin^2\theta} = { sin2\theta \over cos 2\theta } = tan2\theta \\ \therefore \theta = {1 \over 2} tan^{−1}({c2 \over {c1 − c3}}) \end{array}$$

아래 코드에 시간이 되시는 분들은 한번 생각해서 유도 해보시면 큰 도움이 되실 것 같아요.

$$ \begin{matrix}  ax^2 + bxy + cy^2 + dx + ey + f &=& 0 \\ \\ \Rightarrow x^2 + b'xy + c'y^2 + d'x + e'y + f' &=& 0  \end{matrix} \tag{ Dividing by a} $$

$$\begin{matrix}  \begin{bmatrix} x_1y_1& {y_1}^2 & x_1 & y_1 & 1 \\ \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ x_ny_n& {y_n}^2 & x_n & y_n & 1 \end{bmatrix} \begin{bmatrix} b'\\ c'\\ d'\\ e'\\  f' \end{bmatrix} &=& \begin{bmatrix} - {x_1}^2\\ \vdots \\ -{x_n}^2 \end{bmatrix} \\ \end{matrix}  \tag{ explicit function form, y=f(x,y)} $$

$$ \begin{matrix} J X &=& Y \\ X &=& {(J^T J)}^{-1} J^T Y \end{matrix} \tag{Pesudo inver }$$

$$ \Large \begin{matrix}   \theta = { 1 \over 2 } \tan^{-1}( { b \over a-c })  \\ c_x = {  2cd - be \over b^2 - 4ac } \\ c_y = {  2ae- bd \over b^2 - 4ac } \\  w = \sqrt{ {ac_x^2 + bc_xc_y+ cc_y^2 -f \over a \cos^2\theta + b\cos\theta\sin\theta + c\sin^2\theta } }  \\ h = \sqrt{ {ac_x^2 + bc_xc_y+ cc_y^2 -f \over a \sin^2\theta - b\cos\theta\sin\theta + c\cos^2\theta } }  \end{matrix} $$

노랑색의 경우 원으로 피팅한 결과 입니다. 당연히 LSM(Least Square Method) 방식으로 맞추다 보니 전체의 모양이 반영이 안되었네요, 파랑색의 경우는 타원의 피팅한 결과 입니다. fitEllipse 함수와 거의 유사하게 나오는 것을 확인 하 실 수 있습니다.

#dependency opencv-python 3.2 
#python version 3.6

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


def fit_circle(data):
	#data [[x1,y],[x2,y2]...]
	xs = data[:,0].reshape(-1,1) 
	ys = data[:,1].reshape(-1,1)

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

	cx = X[0,0]
	cy = X[1,0]
	c = X[2,0]
	r = np.sqrt(cx ** 2 + cy ** 2 - c)
	return (cx,cy,r)



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];
        # To do implementation
        #a,b,c,d,e,f 를 통해 theta, 중심(cx,cy) , 장축(major), 단축(minor) 등을 뽑아 낼 수 있어요

	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))

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

def main(img_path):

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

	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) 
	_, contours , _ = cv2.findContours(th,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)

	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):

			a,b,r = fit_circle(con.reshape(-1,2))
			#cv2.drawContours(src,[con],0,color_list[2],2)
			cv2.circle(src,(int(a),int(b)),int(r),color_list[1])
			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[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])