MPFITELLIPSE 함수를 이용한 타원형 궤적 근사

지난 주에는 2차원 공간상의 세 점을 모두 지나가는 원의 궤적을 구하고 그려보는 예제를 소개한 바 있습니다. 오늘은 개념을 좀 더 확장해서, 2차원 공간상에 다수의 점들이 존재할 때 이 점들을 모두 지나가지는 않지만 가장 비슷하게 궤적을 그리는 원(또는 타원)의궤적을 구하고 그려보는 예제를 소개하고자 합니다. 이러한 작업은 이미지상에서 찾고자 하는 형태를 근사하는데 있어서 유용하게 사용될 수 있는 방법이기도 합니다. 원 또는 타원의 형태를 근사하는 작업에 있어서는 IDL 커뮤니티에서 이미 유명한 MPFIT 라이브러리에 있는 MPFITELLIPSE 함수를 사용하고자 합니다. 이 함수를 사용하기 위해서는 MPFIT 라이브러리를 설치하고 경로 설정을 해두어야 하는데요. 이 라이브러리를 받을 수 있는 웹페이지는 다음과 같습니다. 이 웹의 첫 페이지에서 중간쯤 보면 MPFIT 패키지를 다운로드할 수 있는 링크가 있습니다. 여기서 ZIP 파일을 받아서 압축을 푼 후 이 디렉토리를 IDL의 경로에 추가해두면 됩니다.

http://www.physics.wisc.edu/~craigm/idl/fitting.html

우선 예제로 사용할 데이터를 만들고자 합니다. 데이터 자체는 각 포인트들에 대한 X, Y 좌표들로 구성되면 됩니다. 즉 2차원 공간상에 N개의 데이터 포인트들이 있다면, X 좌표값들끼리 모인 배열과 Y 좌표값들끼리 모인 배열만 있으면 됩니다. 물론 각 배열은 N개의 값들로 구성되어 있겠지요. 그래서 전반적으로 원의 궤적을 그릴 것 같은 그럴싸한 가상 데이터를 만들어보기 위하여 다음과 같이 약간의 테크닉을 좀 사용해 봤습니다. 중요한 것은 최종적으로 산출된 x, y 배열이 오늘 사용될 예제 데이터가 된다는 점입니다.

data = HANNING(200, 200)*100 + RANDOMN(seed, 200, 200)*10 + 50

ww = WHERE(data GT 115 AND data LT 120, count)

inds = ARRAY_INDICES(data, ww)

x = inds[0, *]

y = inds[1, *]

이렇게 가상으로 만들어본 x, y 포인트들의 갯수는 대략 800~1000개 정도입니다. 생성에 있어서 RANDOMN 함수를 사용했기 때문에 세부적인 갯수는 매번 조금씩 다릅니다. 분포하는 위치는 X, Y 좌표상으로 대략 50~150의 값 범위내에 있도록 생성해 보았습니다. 이 포인트들의 분포가 실제로 어떤 모습인지를 먼저 보기 위하여 다음과 같이 SCATTERPLOT 함수를 사용하여 그림을 그려 볼 수 있는데, 그 모습은 다음 그림과 같습니다. 눈으로 보기에도 이 점들이 대략적으로 원 또는 타원형의 궤적을 그릴 것 같다는 느낌이 오는 것 같습니다.

win = WINDOW(DIMENSIONS=[500, 500])

pl = SCATTERPLOT(x, y, XRANGE=[50, 150], YRANGE=[50, 150], SYMBOL='circle', $

  /SYM_FILLED, SYM_COLOR='orange', /CURRENT)

이제 이 포인트들을 가장 비슷하게 근사하는 원 또는 타원의 궤적을 구할 차례입니다. 앞서 언급했던 MPFIT 라이브러리가 설치된 상태라면 MPFITELLIPSE 함수를 바로 이용하면 됩니다. 기본적인 이용 방법은 다음과 같이 비교적 간단합니다. 위에서 생성했던 포인트 데이터인 x, y를 인자로 투입하면 됩니다.

p = MPFITELLIPSE(x, y, /QUIET, /TILT)

여기서 /QUIET 키워드는 타원형으로 근사해나가는 과정을 출력하는냐 아니면 조용히 작업하느냐를 결정하는 역할을 합니다. 내부적으로 여러 회차의 iteration을 거치게 되는데 이 과정이 궁금하다면 이 키워드를 사용하지 않으면 됩니다. 그리고 /TILT 키워드는 근사되는 타원의 형태가 기울기를 가질 수 있도록 하는 역할을 합니다. 만약 기울어지지 않은 타원의 형태로 근사하고자 한다면 이 키워드를 사용하지 않으면 됩니다. 물론 이외에도 여러가지 키워드들이 지원이 되는데, 자세한 내용은 소스코드에 있는 주석문의 내용을 참조하면 됩니다. 이 MPFITELLIPSE 함수가 사용이 되면 그 결과는 기본적으로 5개의 값을 담은 배열의 형태로 전달됩니다. 이 값들의 세부내용은 소스코드의 주석문을 인용하면 다음과 같이 근사된 타원의 기하학적 인자값들에 해당됨을 알 수 있습니다. 다만 유의할 점은, /TILT 키워드가 사용될 경우 5개의 인자들이 산출되지만, /TILT 키워드가 사용되지 않을 경우에는 4개, /CIRCLE 키워드를 사용하여 아예 원으로 근사한 경우에는 3개의 인자들만 산출된다는 점입니다. 매 경우마다 기하학적으로 산출 가능한 인자들의 갯수가 약간씩 달라지기 때문입니다.

;      P[0]   Ellipse semi axis 1

;      P[1]   Ellipse semi axis 2   ( = P[0] if CIRCLE keyword set)

;      P[2]   Ellipse center - x value

;      P[3]   Ellipse center - y value

;      P[4]   Ellipse rotation angle (radians) if TILT keyword set

이와 같은 방법으로 타원의 궤적을대변하는 기하학적 인자들인 장축/단축 길이, 중심좌표, 기울기 등을 구했다면 다음 단계는 무엇이 되어야 할까요? 당연히 이 궤적을 그림을 그려보는 과정이 필요합니다. 이를 위해서 타원의 궤적을 구성하는 점들의 좌표를 구하고이를 POLYLINE 함수에 투입하는 다음과 같은 방법을 사용하였습니다. 약간 복잡하게 보일 여지가 있긴 하지만, 그냥 정형화된 방법이라고 간주하고 그대로 사용하시면 되지 않을까 싶습니다. 어쨌든 이와 같은 방법으로 궤적을 선으로 그려보면 그 결과는 다음 그림과 같습니다.

phi = FINDGEN(101)*2*!pi/100

xm = p[2] + p[0]*COS(phi)*COS(p[4]) + p[1]*SIN(phi)*SIN(p[4])

ym = p[3] - p[0]*COS(phi)*SIN(p[4]) + p[1]*SIN(phi)*COS(p[4])

line2 = POLYLINE(xm, ym, COLOR='crimson', THICK=2, /DATA)

예제로 만든 가상 데이터 포인트들의 분포 자체가 거의 원형에 가깝다보니, /TILT 키워드를 사용했음에도 불구하고 기울기가 아주 두드러지게 눈에 띄지는 않습니다. 만약 /TILT 키워드를 사용하지 않고 결과를 얻은 상황이라면, 위의 식에서 p[4]를 0으로 대체해버리면 됩니다. 어쨌든 포인트들의 분포에 가장 근접하는 궤적의 타원을 얻었습니다. 이번에는 또 다른 형태의 가상 데이터 포인트를 생성해서 유사한 작업을 진행해보고자 합니다. 가상 데이터 생성을 위하여 맨 처음에 소개했던 예제 코드의 내용에서 WHERE 함수가 사용된 줄의 내용만 다음과 같이 변경해 봅시다. 나머진 그대로입니다.

ww = WHERE(data GT 150, count)

이렇게 변경된 가상 데이터를 아까와 마찬가지로 SCATTERPLOT 함수로 표출해보면 그 모습은 다음 그림과 같습니다. 이번에 생성한 가상 데이터 포인트들의 분포는 아까와는 달리 속이 비어있지 않고 거의 원형에 가깝습니다.

이후의 과정들은 아까와 동일한 방식으로 수행하면 됩니다. 즉, MPFITELLIPSE 함수를 사용한 후 그 결과를 바탕으로 타원의 궤적을 POLYLINE 함수를 사용하여 그려보면 그 모습은 다음 그림과 같습니다. 대략적으로 원형의 외곽 경계 부분에 해당되는 궤적이 그려지는 것을 확인할 수 있습니다.

이번에는 MPFITELLIPSE 함수에서 지원되는 WEIGHTS라는 키워드를 주목해보고자 합니다. 단어 자체에서 짐작할 수 있듯이, 대상이 되는 모든 포인트들을 동등하게 취급하는 것이 아니라 weight를 줌으로써 차별적으로 비중을 두고자 할 경우 사용 가능한 키워드입니다. 물론 이 경우에는 weight를 대변하는 XY 공간상의 함수가 준비되어야 합니다. 어떤 형태의 함수를 사용할 것인가는 전적으로 사용자의 판단에 따른 문제입니다. 그래서 여기서는 중심부에 위치한 포인트들의 비중을 높이고 가장자리에 위치한 포인트들의 비중을 낮추는 방식으로 작업을 해보겠습니다. 이를 위하여 다음과 같은 형태의 weights라는 함수를 만들어 보았습니다. 분포의 중심이 대략 (100, 100) 부근이기 때문에 이 부분의 비중을 높이고 주변의 비중을 낮추는 형태로 가상으로 만들어본 함수입니다. 그리고 이 함수의 형태를 보기 위하여 다음과 같이 SURFACE 함수를 사용해 보았습니다. 이렇게 표출해본 weight 함수의 형태는 다음 그림과 같습니다.

weights = 1.0/EXP(((x-100)^2+(y-100)^2)/5000)

winw = WINDOW(DIMENSIONS=[500, 500])

sf = SURFACE('1.0/EXP(((x-100)^2+(y-100)^2)/5000)', $

  XRANGE=[50, 150], YRANGE=[50, 150], /CURRENT)

이제 이 weight 함수를 MPFITELLIPSE 함수에 투입해야 하는데, 다음과 같이 WEIGHTS 키워드를 사용하면 됩니다.  이렇게 하여 얻은 결과를 표출하는 방법도 앞서 했던 과정들과 동일합니다. 이렇게 결과를 표출해보면 그 모습은 다음 그림과 같습니다.

p = MPFITELLIPSE(x, y, /QUIET, /TILT, WEIGHTS=weights)

이전의 결과와 비교해서 보면, 우리가 가상으로 만들어서 투입한 weight 함수의 형태가 반영되어 있음을 충분히 알 수 있습니다. 즉, 중심부에 위치한 포인트들의 비중을 상대적으로 더욱 높였기 때문에, 근사되는 타원의 궤적이 중심부에 좀 더 가까운 형태로 얻어진 것입니다. 따라서 대상이 되는 포인트들 중 중심 부근의 것들이 비중을 둘 것인지 아니면 반대로 가장자리 부근의 것들에 비중을 둘 것인지에 따라 다양한 형태의 weight 함수를 적용하여 근사하는 것도 가능합니다.

여기서 사용한 MPFITELLIPSE 함수에 관한 설명을 보면, 이 함수가 항상 잘 작동하지는 않을 수 있으며 특히 이심률이 큰 길쭉한 형태의 타원으로 근사하는 것은 어려울 수 있다는 언급이 있습니다. 따라서 극단적인 경우가 아니라면 이 함수를 사용하여 원 또는 타원의 형태를 얻는 것이 충분히 가능하리라 생각이 됩니다.