Melotopia

모르는 원의 중심

작성자

snowall

in"의 한국어 번역은 "안"입니다.

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

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

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

(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차식을 찾아야 할 때에 최소제곱법을 계산하는 방법을 알아낸 것 같긴 한데, 잘 안된다.

—

결국 정답을 검색했다.

circle_fit.pdf에 액세스하려면 클릭하세요.



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)이다. 점이 더 많아지면 더 정확해 질수도 있겠지만…

코멘트

“모르는 원의 중심”에 대한 4개 응답

  1. 푸른달팽이 아바타
    푸른달팽이

    완전한 원 위의 점이 아니라 확률분포로 다소 흩어져있군요..

    뭐 그렇다면… 최소자승법까지 가지 않더라도

    임의의 세 점을 계속 뽑아가며 그 세 점을 지나는 원의 중심을 각각 구하는 것을 반복해나가면서

    구해진 각각의 중심점들의 평균검을 구하는 것도 좋은 방법이 되지 않을까 싶네요.

    응답
  2. snowall 아바타
    snowall

    그렇게 풀 수는 있는데, 저 점들이 정확히 정의된 원 위에 있는게 아니라 어느정도 오차를 갖고 있기 때문에, 최적값을 구하는 게 좋습니다.

    응답
  3. 푸른달팽이 아바타
    푸른달팽이

    r도 필요하신거면… a, b값을 식 ⓐ에 다시 대입하면

    r=328.2990257

    응답
  4. 푸른달팽이 아바타
    푸른달팽이

    그냥 이렇게 풀면 안될까요..?

    원의 공식 (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도 대입하여 풀 수 있지만 필요 없으므로 생략)

    응답

댓글 남기기

이 사이트는 Akismet을 사용하여 스팸을 줄입니다. 댓글 데이터가 어떻게 처리되는지 알아보세요.

더 많은 게시물