[ Lagrange multiplier ] 두 타원 간의 최소거리

* shortest path, Lagrange multiplier, constrained optimization 

* 라그랑주 승수법

이 번주 너무 바쁜 핑계로 후배에게 한가지 일을 시켰습니다. 

" 영상에서 두 blob 간에 최소거리를 측정하는 모듈을 만들어서 줘 " 

후회를 시키지않는 똑똑한 후배는 역시 해당 모듈을 뚝딱 만들어서 보내줬습니다. 

그런데 한가지 단점이 보이기 시작했네요, 속도가 너무 느린 것 같아요. 

두 blob 의 중심을 이용하여 외각선을 추려서 해당 point to point 비교로 최소 거리를 측정하여 

이중 for 문으로 $O(n^2)$ 가 나와 버리는 단점이 생겼네요,,,

이러한  속도 문제점  때문에 numerical method 방식으로 문제를 풀어 볼 수는 없을까요?

제한 조건이 있는 최적화 문제

• 라그랑주 승수법(Lagrange multiplier)

수식으로 표시할 수 있는 제한조건의 최적화 문제는 라그랑주 승수법(Lagrange multiplier)을 사용하여 최적화 할 수 있습니다. 

위키피디아 (https://en.wikipedia.org/wiki/Lagrange_multiplier) 에 나와 있는 내용처럼 제약된 문제를 제약이 없는 문제로 바꿔서 풀 수 있습니다. 

우선 아래와 같이 최적화 함수와 제한 조건을 정의 합니다. 제한 조건의 등호에 유의 하셔야 합니다.

최적화 하려는 함수  $f(x,y)$

제한 조건  $g(x,y) = 0$

f, g 1차 편미분에 대해 연속이라 가정하면, 새로운 변수 λ(Lagrange multiplier)  도입하여 
Lagrange function 정의 할 수 있습니다. 단순히 2개 함수를 하나의 표현식으로 나타낼 수 있습니다. 

$ \mathcal{L}(x,y,\lambda) = f(x,y) -\lambda \cdot g(x,y)$

[ General case ]

$ \mathcal{L}(x_1,\cdots,x_n,\lambda_1,\cdots,\lambda_M) = f(x_1,\cdots,x_n) - \sum_{k=1}^M \lambda_kg_k(x_1,\cdots,x_n) $

위키피디아에서 설명한 것처럼  g=0 의 둘레를 돌다 보면 f 와 평행해 보이는 곳이 나타날 것인데, 이러한 지점들이 minma, maxima 중에 하나가 될 것입니다. 

이 의미를 잘 해석해보면 2가지로 해석이 되는데 

1. f, g 두 함수가 평행해 지는 곳

2. f 함수 자체가 변하지 않아 g=0 의 둘레 어느 곳이든 똑같아 보이는 현상

$ \nabla_{x,y,z} \mathcal{L}(x,y,\lambda) = 0 \iff  \begin{cases} \nabla_{x,y}f(x,y)=\lambda\nabla_{x,y}g(x,y), \\ g(x,y)=0, \end{cases}$

[ General case ]

$ \nabla_{x_1,\cdots,x_n,\lambda_1,\cdots,\lambda_M} \mathcal{L}(x_1,\cdots,x_n,\lambda_1,\cdots,\lambda_M) = 0 \iff  \begin{cases} \nabla f(\mathbf{x}) -\sum_{k=1}^M \lambda_k g_k(\mathbf{x}) , \\ g_1(\mathbf{x})=\cdots=g_M(\mathbf{x}) = 0, \end{cases}$

위 수식에서 변수 n + M  식이 n + M 가 존재하기 때문에 충분히 풀 수가 있습니다. 

다들 아시겠지만 gradient 의 결과는 vector 나오기 때문에 n x 1 의 column vector 로   n 개의 등식이 존재합니다.

(g1= ... = gM = 0  여기서 M 개, gradient f  = lambda g 에서  n 개)

설명을 위해 머나먼 길을 돌아온 것 같네요 ^^

 Lagrange multiplier 를 이용한 최적화 문제

 문제

두 타원의 최소거리를 구하려고 합니다. 물론 각 점들은 타원 1과 타원 2 안에 있는 후보 중 어느 2개의 점입니다.

최적화 함수 

$d^2 = [(x_1 -x_2)^2 + (y_1 - y_2)^2].$

제한 조건

$g(x_1,y_1)=a_1 x_1^2 + b _1y_1^2 + c_1 x _1y_1 + d_1 x_1 + e _1y_1 + f_1 =0.$

$g(x_2,y_2)=a_2 x_2^2 + b _2y_2^2 + c_2 x _2y_2 + d_2 x_2 + e _2y_2 + f_2 =0.$

[ 원본 영상 - 손으로 그려서 조잡합니다...] 

[ 두 타원(blob) 간의 최소 거리 ] 

코드

#!/usr/bin/env python

from __future__ import division
from __future__ import print_function
from __future__ import absolute_import


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

import scipy.optimize as sp

def fit_rotated_ellipse(data):

    """
    Least Square Mean 방식으로 타원 피팅
    model -> ax^2 + by^2 + cxy + dx + ey + f = 0 
    """
    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];

    return np.array([a,b,c,d,e,f])


def find_shorted_ellipses(contours,img_path):

    """
    opencv 에서 제공하는 fitEllipse 함수와
    실제 어느정도 차이가 나는지 확인을 할 수 있도록
    fitEllipse 함수를 통해서 타원을 그려보겠습니다.
    """

    ellipse1 = cv2.fitEllipse(contours[0])
    ellipse2 = cv2.fitEllipse(contours[1])


    e1 = fit_rotated_ellipse(contours[0].reshape(-1,2))
    e2 = fit_rotated_ellipse(contours[1].reshape(-1,2))


    # subject
    distance2e = lambda x : (x[0]-x[2])**2 + (x[1]-x[3])**2


    # constraint equation 1
    e1const = lambda x : \
            e1[0]*(x[0]**2) + e1[1]*(x[0]*x[1]) + \
            e1[2]*(x[1]**2) + e1[3]*x[0] + e1[4]*x[1] + e1[5]


    # constraint equation 2
    e2const = lambda x : \
            e2[0]*(x[2]**2) + e2[1]*(x[2]*x[3]) + \
            e2[2]*(x[3]**2) + e2[3]*x[2] + e2[4]*x[3] + e2[5]

    rst = sp.fmin_slsqp(distance2e,np.array([450,600,550,200]),eqcons=[e1const,e2const])
    print(rst)

    src = cv2.imread(img_path,cv2.IMREAD_COLOR)
    cv2.ellipse(src,ellipse1,(200,30,30),2)
    cv2.ellipse(src,ellipse2,(200,30,30),2)
    cv2.line(src,(int(rst[0]),int(rst[1])),(int(rst[2]),int(rst[3])),(20,250,20),3)
    cv2.imwrite('shorted.jpg',src)



    return

def find_contours(img_path):
    """
    img_path : 입력 이미지 path
    return : 2 개의 타원의 외각선을 찾아 retrun [외각선1,외각선2]
    """

    src = cv2.imread(img_path,cv2.IMREAD_GRAYSCALE)

    h, w = src.shape[:2]
    mask = np.zeros((h+2, w+2), np.uint8)

    _,thv = cv2.threshold(src,0,255,cv2.THRESH_BINARY_INV+cv2.THRESH_OTSU)
    cv2.floodFill(thv,mask,(w-1,h-1),0)

    ret_contours= []

    _,cts,_= cv2.findContours(thv.copy(),cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_NONE)
    for con in cts:
        p = len(con)
        a = cv2.contourArea(con)
        if(a > 2000 and (1.0-p**2)/(4.0*np.pi*a) < 0.1):
            ret_contours.append(con)

    return ret_contours


if __name__ == '__main__':

    """
    usage> python  two_ellipse.py --img_path=...
    """
    parser = argparse.ArgumentParser()
    parser.add_argument("--img_path",help="the image path")
    if(len(sys.argv) != 2):
        parser.print_help()
        parser.exit()

    args = parser.parse_args()
    cons = find_contours(args.img_path)
    find_shorted_ellipses(cons,args.img_path)

[ 참고 ]

Fabian, how to calculate minimum distance between two arbitrary ellipses in 2d, 

MATHEMATICS,LINK [2012.9.10]

데이터 사이언스 스쿨, 제한조건이 있는 최적화 문제. 2017-06-09, LINK