위와 같은 사진이 있다. 저 사진으로부터 녹색 테두리를 갖는 노란색 원의 중심이 어디인지 정확히 결정해야 한다. 좌표값은 픽셀로 주어져 있다.
원본 사진은 보안관계상 공개할 수 없고, 아무튼 내가 얻은 이미지는 저렇게 생겼다.
그래서, 일단 저 원을 이루는 녹색 테두리의 여러 점들을 찍어서 다음과 같은 값을 얻어냈다.
(156,38)
(33,492)
(8,216)
(43,143)
(42,506)
이 점들은 하나의 원 위에 있는 점이다. (원한다면 수백개의 좌표를 찍어볼 수도 있다.)
일단, 헤론의 공식을 응용한 외접원의 공식을 사용하면 반지름은 얻을 수 있다.
R = abc/(4sqrt(s(s-a)(s-b)(s-c)))
a,b,c는 외접원을 결정하는 삼각형의 세 변의 길이가 된다.
문제는 오차가 생긴다는 점.
최소제곱법을 사용하고 싶은데, 원의 방정식은 다음과 같다.
이 원의 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
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)이다. 점이 더 많아지면 더 정확해 질수도 있겠지만…
완전한 원 위의 점이 아니라 확률분포로 다소 흩어져있군요..
뭐 그렇다면… 최소자승법까지 가지 않더라도
임의의 세 점을 계속 뽑아가며 그 세 점을 지나는 원의 중심을 각각 구하는 것을 반복해나가면서
구해진 각각의 중심점들의 평균검을 구하는 것도 좋은 방법이 되지 않을까 싶네요.
그렇게 풀 수는 있는데, 저 점들이 정확히 정의된 원 위에 있는게 아니라 어느정도 오차를 갖고 있기 때문에, 최적값을 구하는 게 좋습니다.
r도 필요하신거면… a, b값을 식 ⓐ에 다시 대입하면
r=328.2990257
그냥 이렇게 풀면 안될까요..?
원의 공식 (x-a)^2 + (y-b)^2 = r^2
임의의 세 점을 지나는 원은 단 하나만 존재하므로,
점 세개만을 택하고 식에 대입하면,
ⓐ: (23, 176) -> (23-a)^2 + (176-b)^2 = r^2
ⓑ: (156,38) -> (156-a)^2 + (38-b)^2 = r^2
ⓒ: (33,492) -> (33-a)^2 + (492-b)^2 = r^2
ⓐ를 풀면 a^2+b^2-r^2-46a-352b+31505=0
ⓑ를 풀면 a^2+b^2-r^2-312a-76b+25780=0
ⓒ를 풀면 a^2+b^2-r^2-66a-984b+243153=0
ⓐ-ⓑ하면 266a-276b+5725=0
그러므로 b=(266a+5725)/276 —-①
ⓑ-ⓒ하면 -246a+908b-217373=0
그러므로 a=(908b-217373)/246
여기에 ①을 대입하면
a={908*(266a+5725)/276)-217373}/246
a=3.557322965a-807.0673972
a=315.5907206 ———②
다시 ①에 ②를 대입하면
b=324.8990278
그러므로 원의 중심은 (315.5907206, 324.8990278)
(r도 대입하여 풀 수 있지만 필요 없으므로 생략)