모르는 원의 중심

위와 같은 사진이 있다. 저 사진으로부터 녹색 테두리를 갖는 노란색 원의 중심이 어디인지 정확히 결정해야 한다. 좌표값은 픽셀로 주어져 있다.

원본 사진은 보안관계상 공개할 수 없고, 아무튼 내가 얻은 이미지는 저렇게 생겼다.

그래서, 일단 저 원을 이루는 녹색 테두리의 여러 점들을 찍어서 다음과 같은 값을 얻어냈다.

(23, 176)
(156,38)
(33,492)
(8,216)
(43,143)
(42,506)

이 점들은 하나의 원 위에 있는 점이다. (원한다면 수백개의 좌표를 찍어볼 수도 있다.)

일단, 헤론의 공식을 응용한 외접원의 공식을 사용하면 반지름은 얻을 수 있다.

s = (a+b+c)/2
R = abc/(4sqrt(s(s-a)(s-b)(s-c)))

a,b,c는 외접원을 결정하는 삼각형의 세 변의 길이가 된다.

문제는 오차가 생긴다는 점.

최소제곱법을 사용하고 싶은데, 원의 방정식은 다음과 같다.

(x-x0)*(x-x0)+(y-y0)*(y-y0) = r*r

이 원의 x0, y0, r을 얻어낼 수 있다면 좋겠다.

그래서 지금 계산중... (계산 끝나면 이 글은 수정됨.)

2차식을 찾아야 할 때에 최소제곱법을 계산하는 방법을 알아낸 것 같긴 한데, 잘 안된다.

---
결국 정답을 검색했다.
http://www.dtcenter.org/met/users/docs/write_ups/circle_fit.pdf
http://www.ulb.ac.be/assoc/bms/Bulletin/sup962/gander.pdf

난 수학적 재능이 없는가보다.

---
위의 알고리즘을 적용하여 만든 프로그램. 파이썬이다.

import numpy.linalg as la
import numpy

A = numpy.matrix([[23,176],
[156,38],
[33,492],
[8,216],
[43,143],
[42,506],
[4,434],
[180,24],
[41,506],
[5,221],
[253,2]])
avgA = 0
avgB = 0
uu = 0
vv = 0
uv = 0
uuu = 0
vvv = 0
uvv = 0
vuu = 0

for i in A:
    avgA+=i[0,0]
    avgB+=i[0,1]

avgA/=len(A)
avgB/=len(A)

for i in A:
    i[0,0]-=avgA
    i[0,1]-=avgB

for i in A:
    uu += i[0,0]*i[0,0]
    vv += i[0,1]*i[0,1]
    uv += i[0,0]*i[0,1]
    uuu += i[0,0]*i[0,0]*i[0,0]
    vvv += i[0,1]*i[0,1]*i[0,1]
    uvv += i[0,0]*i[0,1]*i[0,1]
    vuu += i[0,1]*i[0,0]*i[0,0]

U = numpy.matrix([[uu,uv],[uv,vv]])
UU = numpy.matrix([[0.5*(uuu+uvv)],[0.5*(vvv+vuu)]])
S = la.inv(U).dot(UU)
print(S[0]+avgA, S[1]+avgB)

그래서 얻은 답은 (323.8, 327.7)이다. 점이 더 많아지면 더 정확해 질수도 있겠지만...