Geometry & Recognition

기준 좌표계에 대해서 원점을 이동하고 좌표축을 회전시킨 새로운 좌표계에서 점의 좌표는 바뀐다. 원래의 좌표와 바뀐 좌표값 사이의 관계를 주는 변환이 Isometric transformation (isometry)이다. 평면에서 이 변환은 평행이동을 나타내는 파라미터 2개, 그리고 1개의 회전각 파라미터에 의해서 결정이 된다. 회전각이 $\theta$고, 평행이동이 $(t_x,t_y)$인 isometry에 의해서 두 점 $(x,y)$가 $(u,v)$로 연결이 되는 경우에, 아래의 식으로 표현이 된다:

$$u=\cos (\theta )x-\sin (\theta )y+t_x$$

$$v=\sin (\theta )x+\cos (\theta )y+t_y;$$

따라서 isometry로 연결이 되는 두 점의 조합 $\{(x_1,y_1)\to (u_1,v_1),(x_2,y_2)\to (u_2,v_2)\}$ 만 있으면 이들 파라미터를 정확히 결정할 수 있다. 그러나 변환에 필요한 점 정보를 얻는 과정은 필연적으로 노이즈의 영향을 받게 되므로 주어진 모든 점을 정확히 연결하는 변환을 일반적으로 구할 수 없다. 이 경우에는 isometry 파라미터는 일반적으로 최소자승법에 의해서 결정될 수 있다. 

최소자승법을 사용하기 위해서는 회전각 $\theta$보다는 $a=\cos \theta$, $b=\sin \theta$로 정의된 새로운 파라미터로 식을 표현하는 것이 더 편리하다. 그러나 이 경우에 파라미터 $a,b$는 서로 독립적이 아니고 $a^2+b^2=1$의 제한 조건을 만족시켜야 한다.  

평행이동 파라미터는 질량중심의 isometry 관계로 해결이 되므로, 이 전체 계산을 각각의 질량중심을 원점으로 하는 좌표로 옮겨서 적용하면 더 이상 평행이동을 고려할 필요 없이 회전만 계산하면 된다.

최소자승법의 원리에 따라 입력점의 isometry 결과와 대응점 사이의 거리의 제곱 합 $L$을 주어진 제약조건 내에서 최소화시키는 파라미터 $a,b,\lambda$를 찾으면 된다:

$$L=\sum _i[(a{x}_i-b{y}_i-u_i{)}^2+(b{x}_i+a{y}_i-v_i{)}^2]+\lambda (a^2+b^2-1);$$

여기서 $\lambda$는 제한 조건 $a^2+b^2=1$를 넣기 위한 Lagrange multiplier이다. 극값을 찾기 위해서 $L$를 각각 $a,b,\lambda$에 대해서 미분해서 다음 조건을 얻는다:

$$\sum _i(a{x}_i-b{y}_i-u_i)x_i+(b{x}_i+a{y}_i-v_i)y_i+\lambda a=0$$

$$\sum _i(a{x}_i-b{y}_i-u_i)(-y_i)+(b{x}_i+a{y}_i-v_i)x_i+\lambda b=0$$

$$a^2+b^2=1$$

이 식들을  $a,b,\lambda$에 대해서 풀면 다음의 관계식을 얻는다:

$$a=\sum (x_i{u}_i+y_i{v}_i)/\sum (x_i^2+y_i^2+\lambda )$$

$$b=\sum (x_i{v}_i-y_i{u}_i)/\sum (x_i^2+y_i^2+\lambda )$$
또한, Lagrange 멀티플라이어 $\lambda$는

$$A=\sum (x_i{u}_i+y_i{v}_i)$$

$$B=\sum (x_i{v}_i-y_i{u}_i);$$

로 놓으면, $a^2+b^2=1$ 에서

$$\sum (x_i^2+y_i^2+\lambda )=\sqrt{A^2+B^2};$$

임을 쓰면 된다. 따라서 회전각은

$$\cos \theta =a=\frac{A}{\sqrt{A^2+B^2}}$$

$$\sin \theta =b=\frac{B}{\sqrt{A^2+B^2}};$$

로 주어진다.

질량중심을 빼기 전 좌표 $(x,y)$의 질량중심과 $(u,v)$의 질량중심은 서로 isometry에 의해서 연결이 되므로, 이 관계에서 평행이동 파라미터 $(t_x,t_y)$는 결정이 된다:
$$(x_c,y_c)\to (u_c,v_c)$$

$$u_c=a{x}_c-b{y}_c+t_x$$

$$v_c=b{x}_c+a{y}_c+t_y;$$

참고:
** affine transformation = similarity transformation + shear;
** similarity transformation = isometry transformation + overall scaling;

/* struct CfPt { double x, y;};
*      u = T[0] * x + T[1] * y +T[4] ;
*      v = T[2] * x + T[3] * y + T[5] ; 
*/
BOOL IsometryTransform(std::vector<CfPt> &A, std::vector<CfPt> &U, double T[6]) {
    // A.size()==U.size();
    double cx = 0, cy = 0;
    double ux = 0, uy = 0;
    for (int i = A.size(); i-->0;) {
        cx += A[i].x ;  cy += A[i].y ;
        ux += U[i].x ;  uy += U[i].y ;
    };
    //center of mass ;
    cx /= A.size(); cy /= A.size();
    ux /= A.size(); uy /= A.size();

    //centering 된 좌표계에서 계산;
    double dot = 0 , cross = 0;
    for (int i = A.size(); i-->0;) {
        double x = A[i].x - cx, y = A[i].y - cy;
        double u = U[i].x - ux, v = U[i].y - uy;
        dot += (x * u + y * v);
        cross += ( x * v - y * u) ;
    };
    double norm = sqrt(dot * dot + cross * cross) ;
    double a = dot / norm ;
    double b = cross / norm ;

    T[0] = a ; T[1] = -b ; T[2] = b; T[3] = a; 
    T[4] = ux - (a * cx - b * cy) ;
    T[5] = uy - (b * cx + b * cy) ;
    return 1;
} ;

2차원 이미지의 기하학적인 변형 중에서 평행이동, 회전 및 전체적인 크기의 변화를 주는 변환이 similarity transformation이다. 이 변환은 두 직선이 이루는 각을 보존하고 길이 비를 유지한다. 따라서 similarity 변환 후 물체의 모양은 변환 전과 같은 형태를 가진다. 이 변환보다도 더 일반적인 2차원의 기하학적인 변환은 affine transformation이다. Affine 변환은 한쪽 방향으로의 밀림(sheer)도 허용한다. 평행한 두 직선은 affine 변환 후에도 여전히 평행하다.

Similarity transformation은 전체적인 크기를 바꾸는 scale parameter($s$) 1개와 회전각($\theta$) 1개, 그리고 $x,y$축으로의 평행이동을 나타내는 parameter ($t_x$, $t_y$) 2 개를 합해서 총 4개가 있어야 한다. 이 parameter에 의해서 원본 이미지의 픽셀 $(x,y)$가 변환된 이미지의 픽셀 $(u,v)$에 대응한다고 하면, 이들 간의 관계는 다음식으로 주어진다.

$$u=s\cos (\theta )x-s\sin (\theta )y+t_x$$

$$v=s\sin (\theta )y+s\cos (\theta )y+t_y$$

따라서 원본 영상의 2점에 대응하는 정보만 주어지면 파라미터 $(s,\theta ,t_x,t_y)$를 유일하게 결정할 수 있다.

$$(x_1,y_1)\to (u_1,v_1)$$

$$(x_2,y_2)\to (u_2,v_2)$$

그러나 많은 경우에는 기준점을 잡는데 에러 등을 고려하여서 일반적으로 원본 영상의 $N(\ge 2)$ 개의 점에 대응하는 정보를 주게 되는데, 이 경우에 변환 관계식은 overdetermined 되어서 해를 구할 수 없는 경우도 있다. 이 경우에는 최소자승법을 써서 변환점과 변환식에 의해서 의해서 주어지는 값의 차이를 최소화시키는 파라미터를 구해서 쓰면 된다.

$$\begin{array}{rl}L=\sum _i& [{|u_i-(s\cos (\theta )x_i-s\sin (\theta )y_i+t_x)|}^2\\ & +{|v_i-(s\sin (\theta )x_i+s\cos (\theta )y_i+t_y)|}^2]\end{array}$$

$$(s,\theta ,t_x,t_y)=\text{argmin}(L)$$

이 식을 최소화시키는 파라미터는 $(a=s\cos (\theta ),b=s\sin (\theta )$로 놓으면) $a,b,t_x,t_y$에 대해서 극값을 가질 조건에서 얻을 수 있다. $$\frac{\partial L}{\partial a}=0:\sum _i(u_i-(a{x}_i-b{y}_i+t_x))(-x_i)+(v_i-(b{x}_i+a{y}_i+t_y))(-y_i)=0$$

$$\frac{\partial L}{\partial b}=0:\sum _i(u_i-(a{x}_i-b{y}_i+t_x))(y_i)+(v_i-(b{x}_i+a{y}_i+t_y))(-x_i)=0$$

$$\frac{\partial L}{\partial t_x}=0:\sum _i(u_i-(a{x}_i-b{y}_i+t_x))=0$$

$$\frac{\partial L}{\partial t_y}=0:\sum _i(vi-(b{x}_i+a{y}_i+t_y))=0.$$

따라서, $S_u=\sum _i{u}_i$, $S_v=\sum _i{v}_i$, $S_{ux}=\sum _i{u}_i{x}_i$, $S_{uy}=\sum _i{u}_i{y}_i$, $S_{vx}=\sum _i{v}_i{x}_i$, $S_{vy}=\sum _i{v}_i{y}_i$, $S_x=\sum x_i$, $S_y=\sum _i{y}_i$, $S_{xx}=\sum _i{x}_i^2$, $S_{xy}=\sum _i{x}_i{y}_i$, $S_{yy}=\sum _i{y}_i^2$ 라고 하면,

$$-S_{ux}+a{S}_{xx}+t_x{S}_x-S_{vy}+a{S}_{yy}+t_y{S}_y=0$$

$$S_{uy}+b{S}_{yy}-t_x{S}_y-S_{vx}+b{S}_{xx}+t_y{S}_x=0$$

$$S_u-a{S}_x+b{S}_y-t_{x}N=0$$

$$S_v-b{S}_x-a{S}_y-t_{y}N=0$$

의 4개의 식을 얻으므로 $(a,b,t_x,t_y)$에 대한 1차 연립방정식을 풀면 된다.

$$\left[\begin{array}{cccc}{S}_x& -S_y& N& 0\\ S_y& S_x& 0& N\\ S_{xx}+S_{yy}& 0& S_x& S_y\\ 0& S_{xx}+S_{yy}& -S_y& S_x\end{array}\right]\left[\begin{array}{c}a\\ b\\ t_x\\ t_y\end{array}\right]=\left[\begin{array}{c}{S}_u\\ S_v\\ S_{ux}+S_{vy}\\ S_{vx}-S_{uy}\end{array}\right]$$

$$\text{or}\mathbf{A}\cdot \mathbf{x}=\mathbf{b}$$

$4\times 4$ 행렬 $\mathbf{A}$의 역행렬은 다음과 같이 쉽게 구해진다. 

$${\mathbf{A}}^{-1}=\frac{1}{S_x^2+S_y^2-N(S_{xx}+S_{yy})}\left[\begin{array}{cccc}{S}_x& S_y& -N& 0\\ -S_y& S_x& 0& -N\\ -(S_{xx}+S_{yy})& 0& S_x& -S_y\\ 0& -(S_{xx}+S_{yy})& S_y& S_x\end{array}\right]$$

아래의 코드는 이것을 구현한 것이다. 물론, $N=2$개인 경우에는 파라미터는 유일하게 정해지고 이보다도 더 간단한 식으로 주어진다.

// dst = (S|T)(src)
BOOL SimilarTransParams(std::vector<CPoint>& src, std::vector<CPoint>& dst, double ST[4]) {
    double Sx = 0, Sy = 0, Sxx = 0, Syy = 0;
    double Su = 0, Sv = 0, Sxu = 0, Sxv = 0, Syu = 0, Syv = 0;
    for (int i = srcPts.size(); i-->0;) {
        double x = src[i].x, y = src[i].y;
        double u = dst[i].x, v = dst[i].y;
        Sx  += x;        Sy  += y;
        Sxx += (x * x);  Syy += (y * y);
        Su  += u;        Sv  += v;
        Sxu += (x * u);  Syv += (y * v);
    }
    double Z = Sxx + Syy;
    double denorm = Sx * Sx + Sy * Sy - src.size() * Z;
    // det = (denorm)^2;
    if (denorm == 0) return FALSE;
    invA[16] = { Sx, Sy, -src.size(),           0,
                -Sy, Sx,           0, -src.size(),
                -Z,   0,          Sx,         -Sy,
                 0,  -Z,          Sy,          Sx};
    for (int i = 0; i < 16; i++) invA[i] /= denorm;
    //
    double b[4] = {Su, Sv, Sxu + Syv, Sxv - Syu};
    for (int i = 0; i < 4; i++) {
    	double s = 0;
        for (int j = 0; j < 4; j++) 
            s += invA[i * 4 + j] * b[j];
        ST[i] = s;
    }
    return TRUE ;
};

InvertMatrix4x4()는 4x4행렬의 역행렬을 구한다(OpenCV에서)

BOOL InvertMatrix4x4_d(double* srcMatr, double* dstMatr) {
    double di = srcMatr[0];
    double d = 1.0 / di;

    dstMatr[0] = d;
    dstMatr[4] = srcMatr[4] * -d;
    dstMatr[8] = srcMatr[8] * -d;
    dstMatr[12] = srcMatr[12] * -d;
    dstMatr[1] = srcMatr[1] * d;
    dstMatr[2] = srcMatr[2] * d;
    dstMatr[3] = srcMatr[3] * d;
    dstMatr[5] = srcMatr[5] + dstMatr[4] * dstMatr[1] * di;
    dstMatr[6] = srcMatr[6] + dstMatr[4] * dstMatr[2] * di;
    dstMatr[7] = srcMatr[7] + dstMatr[4] * dstMatr[3] * di;
    dstMatr[9] = srcMatr[9] + dstMatr[8] * dstMatr[1] * di;
    dstMatr[10] = srcMatr[10] + dstMatr[8] * dstMatr[2] * di;
    dstMatr[11] = srcMatr[11] + dstMatr[8] * dstMatr[3] * di;
    dstMatr[13] = srcMatr[13] + dstMatr[12] * dstMatr[1] * di;
    dstMatr[14] = srcMatr[14] + dstMatr[12] * dstMatr[2] * di;
    dstMatr[15] = srcMatr[15] + dstMatr[12] * dstMatr[3] * di;
    di = dstMatr[5];
    dstMatr[5] = d = 1.0 / di;
    dstMatr[1] *= -d;
    dstMatr[9] *= -d;
    dstMatr[13] *= -d;
    dstMatr[4] *= d;
    dstMatr[6] *= d;
    dstMatr[7] *= d;
    dstMatr[0] += dstMatr[1] * dstMatr[4] * di;
    dstMatr[2] += dstMatr[1] * dstMatr[6] * di;
    dstMatr[3] += dstMatr[1] * dstMatr[7] * di;
    dstMatr[8] += dstMatr[9] * dstMatr[4] * di;
    dstMatr[10] += dstMatr[9] * dstMatr[6] * di;
    dstMatr[11] += dstMatr[9] * dstMatr[7] * di;
    dstMatr[12] += dstMatr[13] * dstMatr[4] * di;
    dstMatr[14] += dstMatr[13] * dstMatr[6] * di;
    dstMatr[15] += dstMatr[13] * dstMatr[7] * di;
    di = dstMatr[10];
    dstMatr[10] = d = 1.0 / di;
    dstMatr[2] *= -d;
    dstMatr[6] *= -d;
    dstMatr[14] *= -d;
    dstMatr[8] *= d;
    dstMatr[9] *= d;
    dstMatr[11] *= d;
    dstMatr[0] += dstMatr[2] * dstMatr[8] * di;
    dstMatr[1] += dstMatr[2] * dstMatr[9] * di;
    dstMatr[3] += dstMatr[2] * dstMatr[11] * di;
    dstMatr[4] += dstMatr[6] * dstMatr[8] * di;
    dstMatr[5] += dstMatr[6] * dstMatr[9] * di;
    dstMatr[7] += dstMatr[6] * dstMatr[11] * di;
    dstMatr[12] += dstMatr[14] * dstMatr[8] * di;
    dstMatr[13] += dstMatr[14] * dstMatr[9] * di;
    dstMatr[15] += dstMatr[14] * dstMatr[11] * di;
    di = dstMatr[15];
    dstMatr[15] = d = 1.0 / di;
    dstMatr[3] *= -d;
    dstMatr[7] *= -d;
    dstMatr[11] *= -d;
    dstMatr[12] *= d;
    dstMatr[13] *= d;
    dstMatr[14] *= d;
    dstMatr[0] += dstMatr[3] * dstMatr[12] * di;
    dstMatr[1] += dstMatr[3] * dstMatr[13] * di;
    dstMatr[2] += dstMatr[3] * dstMatr[14] * di;
    dstMatr[4] += dstMatr[7] * dstMatr[12] * di;
    dstMatr[5] += dstMatr[7] * dstMatr[13] * di;
    dstMatr[6] += dstMatr[7] * dstMatr[14] * di;
    dstMatr[8] += dstMatr[11] * dstMatr[12] * di;
    dstMatr[9] += dstMatr[11] * dstMatr[13] * di;
    dstMatr[10] += dstMatr[11] * dstMatr[14] * di;
    return TRUE;
}

2개의 대응점만 주어진 경우 $(x_1,y_1),(x_2,y_2)\to (u_1,v_1),(u_2,v_2)$;

bool SimilarTransParams(double x1, double y1, double x2, double y2, 
                        double u1, double v1, double u2, double v2,
                        double ST[4]) {
    double x21 = x2 - x1, y21 = y2 - y1;
    double u21 = u2 - u1, v21 = v2 - v1;
    double det = x21 * x21 + y21 * y21;
    if (det == 0.) return false;
    double a = (x21 * u21 + y21 * v21) / det ;
    double b = (x21 * v21 - y21 * u21) / det ;
    double tx = u1 - a * x1 + b * y1;
    double ty = v1 - b * x1 - a * y1;
    ST[0] = a; ST[1] = b; ST[2] = tx; ST[3] = ty;
    return true;
};

얼굴인식용 training data set을 만들기 위해서 얼굴을 정렬시키는 데 사용한 예:
- 양 눈의 위치 변환: (70,93), (114, 84) --> (30,45), (100,45)로 변환( linear interpolation사용)
- 실제로 사용되는 변환은 정해진 dst영역으로 매핑하는 src영역을 찾아야 하므로, 역변환이 필요하다.
- 필요한 역변환은 src와 dst의 역할만 바꾸면 쉽게 구할 수 있다.

RANSAC 알고리즘을 써서 주어진 2차원 점집합에서 원을 추정한다. 원을 만들기 위해서는 최소한 3점이 필요하고, 또 일직선에 있지 않아야 한다. 이렇게 만들어진 원은 세 점을 꼭짓점으로 하는 삼각형의 외접원(circumcircle)이다. 주어진 외접원에서 크게 벗어나지 않는 inliers를 찾는다(추가로 이 inliers에 대해 최소자승법으로 원의 중심과 반지름을 다시 구해서 보다 정밀하게 추정하는 과정을 넣을 수도 있다). 무작위로 선택된 세 점에 대해 위의 과정을 반복 시행해서 구한 원 중에서 가장 많은 inliers를 포함하는 원을 결과로 사용한다.

// 참고: http://en.wikipedia.org/wiki/RANSAC

// 2024.5.31 재작성;
double dist_deviate(const CfPt& pts, double cparam[3]) {
    double dx = pts.x - cparam[0];
    double dy = pts.y - cparam[1];
    return fabs(hypot(dx, dy) - cparam[2]);
}
int circumcircle(CfPt pts[3], double cparam[3]) {
    double x1 = pts[0].x, x2 = pts[1].x, x3 = pts[2].x;
    double y1 = pts[0].y, y2 = pts[1].y, y3 = pts[2].y;
    double bax = x2 - x1, bay = y2 - y1;
    double cax = x3 - x1, cay = y3 - y1;
    double E = bax * (x1 + x2) + bay * (y1 + y2);
    double F = cax * (x1 + x3) + cay * (y1 + y3);
    double G = 2. * (bax * (y3 - y2) - bay * (x3 - x2));
    if (G == 0.) return 0;    //error;
    //assert(fabs(G)>small_epsilon); //to prevent collinear or degenerate case;
    cparam[0] = (cay * E - bay * F) / G; //cx;
    cparam[1] = (bax * F - cax * E) / G; //cy;
    cparam[2] = hypot(cparam[0]-x1, cparam[1]-y1); //rad;
    return 1;
};
int num_sampling3(double prob_fail, double inlier_ratio) {
    return int(log(prob_fail)/log(1-pow(inlier_ratio, 3))); 
}
std::vector<int> Ransac_CircleFit(std::vector<CfPt>& points, double circle_param[3]) {
    if (points.size() < 3)
        return std::vector<int> (); //return null_vector;

    CfPt center; double inv_scale;
    // normalize input points for the sake of numerical stability;
    std::vector<CfPt> nor_pts = normalize(points, inv_scale, center);
    // distance threshold;
    double distance_thresh = sqrt(double(points.size())) * inv_scale;
    //ransac
    int sample_num = 1000;	//number of sample
    int ransac_count = 0;
    const double prob_fail = 0.01;
    double best_cparam[3] = {0};
    std::vector<int> best_inliers;
    while (sample_num > ransac_count) {
        // pick random 3 indices:[0,points.size()-1];
        int triple[3];
        random_triple(points.size()-1, triple);
        CfPt selected[3];
        for (int i = 0; i < 3; i++) 
            selected[i] = nor_pts[triple[i]];
        // circumcircle of 3 points;
        if (circumcircle(selected, circle_param)) {
            // find inliers;
            std::vector<int> inliers;
            inliers.reserve(points.size());
            for (int i = nor_pts.size(); i-->0;) {
                // error measure = algebric distance;
                double deviation = dist_deviate(nor_pts[i], circle_param);
                if (fabs(deviation) < distance_thresh)
                    inliers.push_back(i);
            }
            if (inliers.size() > best_inliers.size()) {			
                // update sampling_num;
                sample_num = num_sampling3(prob_fail, double(inliers.size())/points.size());
                // update best_inliers;
                best_inliers.swap(inliers);
                // update best circle param;
                for (int i = 0; i < 3; i++) 
                    best_cparam[i] = circle_param[i];
            }
        }
        if (++ransac_count > 1500) {
            TRACE("error! ransac_count exceed!\n");
            break;
        }
    }
    // recover original coordinate and scale;
    denormalize(best_cparam, best_cparam, inv_scale, center);
    if (best_cparam[0] > 0 && best_cparam[1] > 0) {
        for (int i = 0; i < 3; i++)
            circle_param[i] = best_cparam[i];
        TRACE("circle_found(%d, %d)\n", sample_num, ransac_count);
        // more accurate estimation needed at this stage;
    } else 
        best_inliers.clear();
    return best_inliers;
}

https://kipl.tistory.com/207

https://kipl.tistory.com/357